1#include "GaudiKernel/MsgStream.h"
13#include "DedxCalibAlg/DedxCalibDocaEAng.h"
19 declareProperty(
"Debug", m_debug=
false);
20 declareProperty(
"Idoca", m_idoca_test=1);
21 declareProperty(
"Ieang", m_ieang_test=1);
22 declareProperty(
"PMin", pMin=0.2);
23 declareProperty(
"PMax", pMax=2.3);
24 declareProperty(
"Truncate", m_truncate=0.7);
29 MsgStream log(
msgSvc(), name());
30 log<<MSG::INFO<<
"DedxCalibDocaEAng::BookHists()"<<endreq;
42 hlabel<<
"dEdx_doca"<<i<<
"_eang"<<j;
47 hlabel<<
"dEdxVsDocaEAng";
53 MsgStream log(
msgSvc(), name());
54 log<<MSG::INFO<<
"DedxCalibDocaEAng::FillHists()"<<endreq;
57 TTree* n103, *t_test(0);
60 double doca_test(0), eang_test(0);
62 double dEdx[100]={0},pathlength[100]={0},wid[100]={0},layid[100]={0},dd_in[100]={0},eangle[100]={0},zhit[100]={0};
65 if(m_debug){ f_test =
new TFile(
"doca_eangle_test.root",
"recreate");
66 t_test =
new TTree(
"inftest",
"information for test doca eangle calibration");
67 t_test->Branch(
"idoca", &m_idoca_test,
"idoca/I");
68 t_test->Branch(
"ieang", &m_ieang_test,
"ieang/I");
69 t_test->Branch(
"doca", &doca_test,
"doca/D");
70 t_test->Branch(
"eang", &eang_test,
"eang/D");
75 cout<<
"runlist: "<<runlist.c_str()<<endl;
76 f =
new TFile(runlist.c_str());
77 n103 = (TTree*)
f->Get(
"n103");
78 if(
f==0){ cout <<
"the file is empty" << endl;
continue;}
79 if(n103==0){ cout <<
"the tree is empty" << endl;
continue;}
98 vector<double> m_phlist_test(0);
99 for(
int j=0;j<n103->GetEntries();j++){
101 if(
ptrk>pMax ||
ptrk<pMin)
continue;
102 if(tes>1400)
continue;
117 for(
int k=0;k<ndedxhit;k++){
119 dd_in[k] = dd_in[k]/doca_norm[int(layid[k])];
120 if(eangle[k] >0.25) eangle[k] = eangle[k] -0.5;
121 else if(eangle[k] <-0.25) eangle[k] = eangle[k] +0.5;
124 if(layid[k] <4)
continue;
131 for(
int i =0; i<14; i++)
133 if(eangle[k] < Eangle_value_cut[i] || eangle[k] >= Eangle_value_cut[i+1])
continue;
134 else if(eangle[k] >= Eangle_value_cut[i] && eangle[k] < Eangle_value_cut[i+1])
136 for(
int m =0; m<i; m++)
138 iEAng += Eangle_cut_bin[m];
140 double eangle_bin_cut_step = (Eangle_value_cut[i+1] - Eangle_value_cut[i])/Eangle_cut_bin[i];
141 int temp_bin = int((eangle[k]-Eangle_value_cut[i])/eangle_bin_cut_step);
145 if(m_debug && (iDoca==m_idoca_test) && (iEAng==m_ieang_test)){
146 doca_test = dd_in[k];
147 eang_test = eangle[k];
149 cout <<
"dedx before cor: " << dEdx[k]*1.5/pathlength[k] << endl;
151 dEdx[k] =
exsvc->
StandardHitCorrec(0,runFlag,2,runNO,evtNO,pathlength[k],wid[k],layid[k],dEdx[k],(dd_in[k])*(doca_norm[
int(layid[k])]),eangle[k],
costheta);
153 m_docaeangle[iDoca][iEAng]->Fill(dEdx[k]);
154 if(m_debug && (iDoca==m_idoca_test) && (iEAng==m_ieang_test)){
155 cout <<
"dedx after cor: " << dEdx[k] << endl << endl;
160 f->Close();
f = 0; n103 = 0;
162 if(m_debug){ f_test->cd(); t_test->Write(); f_test->Close(); }
167 MsgStream log(
msgSvc(), name());
168 log<<MSG::INFO<<
"DedxCalibDocaEAng::AnalyseHists()"<<endreq;
171 Double_t entry=0,mean=0,rms=0;
174 gStyle->SetOptFit(1111);
179 entry = m_docaeangle[id][ip]->GetEntries();
180 if(entry<10)
continue;
181 mean = m_docaeangle[id][ip]->GetMean();
182 rms = m_docaeangle[id][ip]->GetRMS();
183 binmax = m_docaeangle[id][ip]->GetMaximumBin();
184 lp[1] = m_docaeangle[id][ip]->GetBinCenter(binmax);
186 lp[0] = (m_docaeangle[id][ip]->GetBinContent(binmax)+m_docaeangle[id][ip]->GetBinContent(binmax-1)+m_docaeangle[id][ip]->GetBinContent(binmax+1))/3;
190 func =
new TF1(
"mylan",
mylan,10,3000,4);
191 func->SetParameters(entry, mean, rms, -0.8);
192 func->SetParLimits(0, 0, entry);
193 func->SetParLimits(1, 0, mean);
194 func->SetParLimits(2, 10, rms);
195 func->SetParLimits(3, -1, 0);
196 func->SetParameters(entry/36., mean-200, rms/2.3, -0.5);
200 func =
new TF1(
"landaun",
landaun,10,3000,3);
201 func->SetParameters(lp[0], lp[1], lp[2] );
205 func =
new TF1(
"Landau",
Landau,10,3000,2);
206 func->SetParameters(0.7*mean, rms );
210 func =
new TF1(
"Vavilov",
Vavilov,10,3000,2);
211 func->SetParameters(0.05, 0.6);
215 func =
new TF1(
"AsymGauss",
AsymGauss,10,3000,4);
216 func->SetParameters(entry, mean, rms, 1.0 );
218 func->SetLineWidth(2.1);
219 func->SetLineColor(2);
221 m_docaeangle[id][ip]->Fit(func,
"rq");
228 MsgStream log(
msgSvc(), name());
229 log<<MSG::INFO<<
"DedxCalibDocaEAng::WriteHists()"<<endreq;
232 double Out_entryNo[totalNo]={0},Out_mean[totalNo]={0},Out_sigma[totalNo]={0},Out_fitmean[totalNo]={0},Out_fitmeanerr[totalNo]={0},Out_fitsigma[totalNo]={0},Out_gain[totalNo]={0};
233 double Iner_entryNo[totalNo]={0},Iner_mean[totalNo]={0},Iner_sigma[totalNo]={0},Iner_fitmean[totalNo]={0},Iner_fitmeanerr[totalNo]={0},Iner_fitsigma[totalNo]={0},Iner_gain[totalNo]={0};
234 double Id_doca[totalNo]={0},Ip_eangle[totalNo]={0},doca[totalNo]={0},eangle[totalNo]={0},chisq[totalNo]={0};
236 double Norm(0),
count(0);
247 if(ip>=iE && ip<iE+Eangle_cut_bin[i])
break;
248 iE=iE+Eangle_cut_bin[i];
251 eangle[
id*
NumSlicesEAng+ip] = Eangle_value_cut[i] + (ip-iE+0.5)*(Eangle_value_cut[i+1]-Eangle_value_cut[i])/Eangle_cut_bin[i];
254 Out_entryNo[
id*
NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetEntries();
255 Out_mean[
id*
NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetMean();
256 Out_sigma[
id*
NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetRMS();
262 Out_fitmean[
id*
NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction(
"mylan")->GetParameter(1);
263 Out_fitmeanerr[
id*
NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction(
"mylan")->GetParError(1);
264 Out_fitsigma[
id*
NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction(
"mylan")->GetParameter(2);
266 chisq[
id*
NumSlicesEAng+ip] = (m_docaeangle[id][ip]->GetFunction(
"mylan")-> GetChisquare())/(m_docaeangle[
id][ip]->GetFunction(
"mylan")-> GetNDF());
270 Out_fitmean[
id*
NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction(
"landaun")->GetParameter(1);
271 Out_fitmeanerr[
id*
NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction(
"landaun")->GetParError(1);
272 Out_fitsigma[
id*
NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction(
"landaun")->GetParameter(2);
273 chisq[
id*
NumSlicesEAng+ip] = (m_docaeangle[id][ip]->GetFunction(
"landaun")-> GetChisquare())/(m_docaeangle[
id][ip]->GetFunction(
"landaun")-> GetNDF());
277 Out_fitmean[
id*
NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction(
"Landau")->GetParameter(1);
278 Out_fitmeanerr[
id*
NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction(
"Landau")->GetParError(1);
279 Out_fitsigma[
id*
NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction(
"Landau")->GetParameter(2);
280 chisq[
id*
NumSlicesEAng+ip] = (m_docaeangle[id][ip]->GetFunction(
"Landau")-> GetChisquare())/(m_docaeangle[
id][ip]->GetFunction(
"Landau")-> GetNDF());
284 Out_fitmean[
id*
NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction(
"Vavilov")->GetParameter(1);
285 Out_fitmeanerr[
id*
NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction(
"Vavilov")->GetParError(1);
286 Out_fitsigma[
id*
NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction(
"Vavilov")->GetParameter(2);
287 chisq[
id*
NumSlicesEAng+ip] = (m_docaeangle[id][ip]->GetFunction(
"Vavilov")-> GetChisquare())/(m_docaeangle[
id][ip]->GetFunction(
"Vavilov")-> GetNDF());
291 Out_fitmean[
id*
NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction(
"AsymGauss")->GetParameter(1);
292 Out_fitmeanerr[
id*
NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction(
"AsymGauss")->GetParError(1);
293 Out_fitsigma[
id*
NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction(
"AsymGauss")->GetParameter(2);
294 chisq[
id*
NumSlicesEAng+ip] = (m_docaeangle[id][ip]->GetFunction(
"AsymGauss")-> GetChisquare())/(m_docaeangle[
id][ip]->GetFunction(
"AsymGauss")-> GetNDF());
297 m_DocaEAngAverdE -> SetBinContent(
id+1, ip+1, Out_fitmean[
id*
NumSlicesEAng+ip]);
305 cout<<
"count= "<<
count<<
" average value of doca eangle bins: "<<Norm<<endl;
306 for(Int_t i=0;i<totalNo;i++)
308 if(Out_entryNo[i]>10 && Out_fitmean[i]>0 && Out_fitsigma[i]>0) Out_gain[i]=Out_fitmean[i]/Norm;
309 else Out_gain[i] = 0;
312 log<<MSG::INFO<<
"begin generating root file!!! "<<endreq;
313 TFile*
f =
new TFile(
m_rootfile.c_str(),
"recreate");
318 m_docaeangle[id][ip]->Write();
321 m_DocaEAngAverdE->Write();
323 TTree* ddgcalib =
new TTree(
"ddgcalib",
"ddgcalib");
324 ddgcalib ->
Branch(
"Out_hits", Out_entryNo,
"Out_entryNo[1600]/D");
325 ddgcalib ->
Branch(
"Out_mean", Out_mean,
"Out_mean[1600]/D");
326 ddgcalib ->
Branch(
"Out_sigma", Out_sigma,
"Out_sigma[1600]/D");
327 ddgcalib ->
Branch(
"Out_fitmean", Out_fitmean,
"Out_fitmean[1600]/D");
328 ddgcalib ->
Branch(
"Out_fitmeanerr", Out_fitmeanerr,
"Out_fitmeanerr[1600]/D");
329 ddgcalib ->
Branch(
"Out_chi", Out_fitsigma,
"Out_fitsigma[1600]/D");
330 ddgcalib ->
Branch(
"Out_gain", Out_gain,
"Out_gain[1600]/D");
331 ddgcalib ->
Branch(
"Iner_hits", Iner_entryNo,
"Iner_entryNo[1600]/D");
332 ddgcalib ->
Branch(
"Iner_mean", Iner_mean,
"Iner_mean[1600]/D");
333 ddgcalib ->
Branch(
"Iner_sigma", Iner_sigma,
"Iner_sigma[1600]/D");
334 ddgcalib ->
Branch(
"Iner_fitmean", Iner_fitmean,
"Iner_fitmean[1600]/D");
335 ddgcalib ->
Branch(
"Iner_fitmeanerr", Iner_fitmeanerr,
"Iner_fitmeanerr[1600]/D");
336 ddgcalib ->
Branch(
"Iner_chi", Iner_fitsigma,
"Iner_fitsigma[1600]/D");
337 ddgcalib ->
Branch(
"Iner_gain", Iner_gain,
"Iner_gain[1600]/D");
338 ddgcalib ->
Branch(
"Id_doca", Id_doca,
"Id_doca[1600]/D");
339 ddgcalib ->
Branch(
"Ip_eangle", Ip_eangle,
"Ip_eangle[1600]/D");
340 ddgcalib ->
Branch(
"chisq", chisq,
"chisq[1600]/D");
341 ddgcalib ->
Branch(
"doca", doca,
"doca[1600]/D");
342 ddgcalib ->
Branch(
"eangle", eangle,
"eangle[1600]/D");
347 TCanvas c1(
"c1",
"canvas", 500, 400);
353 m_docaeangle[id][ip]->Draw();
363 delete m_docaeangle[id][ip];
366 delete m_DocaEAngAverdE;
curve Branch("CurveSize",&CurveSize,"CurveSize/I")
data SetBranchAddress("time",&time)
DOUBLE_PRECISION count[3]
double Landau(double *x, double *par)
double Vavilov(double *x, double *par)
double landaun(double *x, double *par)
double AsymGauss(double *x, double *par)
double mylan(double *x, double *par)
DedxCalibDocaEAng(const std::string &name, ISvcLocator *pSvcLocator)
vector< string > m_recFileVector
virtual double StandardHitCorrec(int calib_rec_Flag, int runFlag, int ntpFlag, int runNO, int evtNO, double pathl, int wid, int layid, double adc, double dd, double eangle, double costheta) const =0
virtual double StandardTrackCorrec(int calib_rec_Flag, int typFlag, int ntpFlag, int runNO, int evtNO, double ex, double costheta, double t0) const =0
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")