1#include "GaudiKernel/MsgStream.h"
12#include "DedxCalibAlg/DedxCalibEAng.h"
19 declareProperty(
"Truncate", m_truncate = 0.7);
20 declareProperty(
"PMin", pMin=0.2);
21 declareProperty(
"PMax", pMax=2.3);
26 MsgStream log(
msgSvc(), name());
27 log<<MSG::INFO<<
"DedxCalibEAng::BookHists()"<<endreq;
39 hlabel<<
"dEdx"<<
"_eang"<<i;
42 hlabel<<
"trunc_dEdx"<<
"_eang"<<i;
45 hlabel<<
"pos_dEdx"<<
"_eang"<<i;
48 hlabel<<
"neg_dEdx"<<
"_eang"<<i;
55 hlabel<<
"trunc_dEdxVsEAng";
58 hlabel<<
"pos_dEdxVsEAng";
61 hlabel<<
"neg_dEdxVsEAng";
68 MsgStream log(
msgSvc(), name());
69 log<<MSG::INFO<<
"DedxCalibEAng::FillHists()"<<endreq;
76 float runNO=0,evtNO=0,runFlag=0,pathlength=0,wid=0,layid=0,dd_in=0,driftdist=0,eangle=0,zhit=0,
costheta=0,tes=0,
ptrk=0,
charge=0;
79 vector<double> m_phlist(0);
80 vector<double> m_phlist_test(0);
81 vector<int> m_iEAng_list(0);
82 double old_costheta(0), old_ptrk(0);
87 cout<<
"runlist: "<<runlist.c_str()<<endl;
88 f =
new TFile(runlist.c_str());
89 n102 = (TTree*)
f->Get(
"n102");
103 n102->SetBranchAddress(
"ptrk1",&
ptrk);
104 n102->SetBranchAddress(
"charge",&
charge);
105 for(
int j=0;j<n102->GetEntries();j++)
109 if((fabs(
costheta-old_costheta)>0.001 || fabs(
ptrk-old_ptrk>0.001)) && m_phlist.size()>=2){
110 sort(m_phlist_test.begin(), m_phlist_test.end());
111 unsigned int nsample = (
unsigned int)(m_phlist_test.size()*m_truncate);
112 double m_trunc_value = m_phlist_test[nsample];
113 for(
unsigned int i=0; i<m_phlist.size(); i++){
114 if(m_phlist[i]>m_trunc_value)
continue;
115 m_trunc_eangle[m_iEAng_list[i]]->Fill(m_phlist[i]);
118 m_phlist_test.clear();
119 m_iEAng_list.clear();
124 if(eangle >0.25) eangle = eangle -0.5;
125 else if(eangle <-0.25) eangle = eangle +0.5;
128 if(
ptrk>pMax ||
ptrk<pMin)
continue;
129 if(tes>1400)
continue;
131 if(layid <4)
continue;
138 dedx =
exsvc->
StandardHitCorrec(0,runFlag,2,runNO,evtNO,pathlength,wid,layid,dedx,dd_in,eangle,
costheta);
140 m_eangle[iEAng]->Fill(dedx);
141 m_phlist.push_back(dedx);
142 m_phlist_test.push_back(dedx);
143 m_iEAng_list.push_back(iEAng);
144 if(
charge>0) m_pos_eangle[iEAng]->Fill(dedx);
145 if(
charge<0) m_neg_eangle[iEAng]->Fill(dedx);
153 MsgStream log(
msgSvc(), name());
154 log<<MSG::INFO<<
"DedxCalibEAng::AnalyseHists()"<<endreq;
157 Double_t entry=0,mean=0,rms=0;
160 gStyle->SetOptFit(1111);
163 entry = m_eangle[ip]->Integral();
164 if(entry<10)
continue;
165 mean = m_eangle[ip]->GetMean();
166 rms = m_eangle[ip]->GetRMS();
167 binmax = m_eangle[ip]->GetMaximumBin();
168 lp[1] = m_eangle[ip]->GetBinCenter(binmax);
170 lp[0] = (m_eangle[ip]->GetBinContent(binmax)+m_eangle[ip]->GetBinContent(binmax-1)+m_eangle[ip]->GetBinContent(binmax+1))/3;
174 func =
new TF1(
"mylan",
mylan,10,3000,4);
175 func->SetParLimits(0, 0, entry);
176 func->SetParLimits(1, 0, mean);
177 func->SetParLimits(2, 10, rms);
178 func->SetParLimits(3, -1, 0);
179 func->SetParameters(entry/36., mean-200, rms/2.3, -0.5);
183 func =
new TF1(
"landaun",
landaun,10,3000,3);
184 func->SetParameters(lp[0], lp[1], lp[2] );
188 func =
new TF1(
"Landau",
Landau,10,3000,2);
189 func->SetParameters(0.7*mean, rms );
193 func =
new TF1(
"Vavilov",
Vavilov,10,3000,2);
194 func->SetParameters(0.05, 0.6);
198 func =
new TF1(
"AsymGauss",
AsymGauss,10,3000,4);
199 func->SetParameters(entry, mean, rms, 1.0 );
201 func->SetLineWidth(2.1);
202 func->SetLineColor(2);
204 m_eangle[ip]->Fit(func,
"rq" );
205 m_pos_eangle[ip]->Fit(func,
"rq");
206 m_neg_eangle[ip]->Fit(func,
"rq");
212 MsgStream log(
msgSvc(), name());
213 log<<MSG::INFO<<
"DedxCalibEAng::WriteHists()"<<endreq;
221 double Norm(0), trunc_Norm(0);
228 entryNo[ip]= m_eangle[ip]->Integral();
229 trunc_mean[ip] = m_trunc_eangle[ip]->GetMean();
230 trunc_sigma[ip] = m_trunc_eangle[ip]->GetRMS();
231 pos_entryNo[ip]= m_pos_eangle[ip]->Integral();
232 pos_mean[ip] = m_pos_eangle[ip]->GetMean();
233 pos_sigma[ip] = m_pos_eangle[ip]->GetRMS();
234 neg_entryNo[ip]= m_neg_eangle[ip]->Integral();
235 neg_mean[ip] = m_neg_eangle[ip]->GetMean();
236 neg_sigma[ip] = m_neg_eangle[ip]->GetRMS();
237 if(entryNo[ip]<10 || pos_entryNo[ip]<10 || neg_entryNo[ip]<10)
continue;
241 fitmean[ip] = m_eangle[ip]->GetFunction(
"mylan")->GetParameter(1);
242 fitmeanerr[ip] = m_eangle[ip]->GetFunction(
"mylan")->GetParError(1);
243 fitsigma[ip] = m_eangle[ip]->GetFunction(
"mylan")->GetParameter(2);
244 chisq[ip] = (m_eangle[ip]->GetFunction(
"mylan")-> GetChisquare())/(m_eangle[ip]->GetFunction(
"mylan")-> GetNDF());
246 pos_fitmean[ip] = m_pos_eangle[ip]->GetFunction(
"mylan")->GetParameter(1);
247 pos_fitmeanerr[ip] = m_pos_eangle[ip]->GetFunction(
"mylan")->GetParError(1);
248 pos_fitsigma[ip] = m_pos_eangle[ip]->GetFunction(
"mylan")->GetParameter(2);
249 pos_chisq[ip] = (m_pos_eangle[ip]->GetFunction(
"mylan")-> GetChisquare())/(m_pos_eangle[ip]->GetFunction(
"mylan")-> GetNDF());
251 neg_fitmean[ip] = m_neg_eangle[ip]->GetFunction(
"mylan")->GetParameter(1);
252 neg_fitmeanerr[ip] = m_neg_eangle[ip]->GetFunction(
"mylan")->GetParError(1);
253 neg_fitsigma[ip] = m_neg_eangle[ip]->GetFunction(
"mylan")->GetParameter(2);
254 neg_chisq[ip] = (m_neg_eangle[ip]->GetFunction(
"mylan")-> GetChisquare())/(m_neg_eangle[ip]->GetFunction(
"mylan")-> GetNDF());
258 fitmean[ip] = m_eangle[ip]->GetFunction(
"landaun")->GetParameter(1);
259 fitmeanerr[ip] = m_eangle[ip]->GetFunction(
"landaun")->GetParError(1);
260 fitsigma[ip] = m_eangle[ip]->GetFunction(
"landaun")->GetParameter(2);
261 chisq[ip] = (m_eangle[ip]->GetFunction(
"landaun")-> GetChisquare())/(m_eangle[ip]->GetFunction(
"mylan")-> GetNDF());
263 pos_fitmean[ip] = m_pos_eangle[ip]->GetFunction(
"landaun")->GetParameter(1);
264 pos_fitmeanerr[ip] = m_pos_eangle[ip]->GetFunction(
"landaun")->GetParError(1);
265 pos_fitsigma[ip] = m_pos_eangle[ip]->GetFunction(
"landaun")->GetParameter(2);
266 pos_chisq[ip] = (m_pos_eangle[ip]->GetFunction(
"landaun")-> GetChisquare())/(m_pos_eangle[ip]->GetFunction(
"landaun")-> GetNDF());
268 neg_fitmean[ip] = m_neg_eangle[ip]->GetFunction(
"landaun")->GetParameter(1);
269 neg_fitmeanerr[ip] = m_neg_eangle[ip]->GetFunction(
"landaun")->GetParError(1);
270 neg_fitsigma[ip] = m_neg_eangle[ip]->GetFunction(
"landaun")->GetParameter(2);
271 neg_chisq[ip] = (m_neg_eangle[ip]->GetFunction(
"landaun")-> GetChisquare())/(m_neg_eangle[ip]->GetFunction(
"landaun")-> GetNDF());
275 fitmean[ip] = m_eangle[ip]->GetFunction(
"Landau")->GetParameter(1);
276 fitmeanerr[ip] = m_eangle[ip]->GetFunction(
"Landau")->GetParError(1);
277 fitsigma[ip] = m_eangle[ip]->GetFunction(
"Landau")->GetParameter(2);
278 chisq[ip] = (m_eangle[ip]->GetFunction(
"Landau")-> GetChisquare())/(m_eangle[ip]->GetFunction(
"Landau")-> GetNDF());
280 pos_fitmean[ip] = m_pos_eangle[ip]->GetFunction(
"Landau")->GetParameter(1);
281 pos_fitmeanerr[ip] = m_pos_eangle[ip]->GetFunction(
"Landau")->GetParError(1);
282 pos_fitsigma[ip] = m_pos_eangle[ip]->GetFunction(
"Landau")->GetParameter(2);
283 pos_chisq[ip] = (m_pos_eangle[ip]->GetFunction(
"Landau")-> GetChisquare())/(m_pos_eangle[ip]->GetFunction(
"Landau")-> GetNDF());
285 neg_fitmean[ip] = m_neg_eangle[ip]->GetFunction(
"Landau")->GetParameter(1);
286 neg_fitmeanerr[ip] = m_neg_eangle[ip]->GetFunction(
"Landau")->GetParError(1);
287 neg_fitsigma[ip] = m_neg_eangle[ip]->GetFunction(
"Landau")->GetParameter(2);
288 neg_chisq[ip] = (m_neg_eangle[ip]->GetFunction(
"Landau")-> GetChisquare())/(m_neg_eangle[ip]->GetFunction(
"Landau")-> GetNDF());
292 fitmean[ip] = m_eangle[ip]->GetFunction(
"Vavilov")->GetParameter(1);
293 fitmeanerr[ip] = m_eangle[ip]->GetFunction(
"Vavilov")->GetParError(1);
294 fitsigma[ip] = m_eangle[ip]->GetFunction(
"Vavilov")->GetParameter(2);
295 chisq[ip] = (m_eangle[ip]->GetFunction(
"Vavilov")-> GetChisquare())/(m_eangle[ip]->GetFunction(
"Vavilov")-> GetNDF());
297 pos_fitmean[ip] = m_pos_eangle[ip]->GetFunction(
"Vavilov")->GetParameter(1);
298 pos_fitmeanerr[ip] = m_pos_eangle[ip]->GetFunction(
"Vavilov")->GetParError(1);
299 pos_fitsigma[ip] = m_pos_eangle[ip]->GetFunction(
"Vavilov")->GetParameter(2);
300 pos_chisq[ip] = (m_pos_eangle[ip]->GetFunction(
"Vavilov")-> GetChisquare())/(m_pos_eangle[ip]->GetFunction(
"Vavilov")-> GetNDF());
302 neg_fitmean[ip] = m_neg_eangle[ip]->GetFunction(
"Vavilov")->GetParameter(1);
303 neg_fitmeanerr[ip] = m_neg_eangle[ip]->GetFunction(
"Vavilov")->GetParError(1);
304 neg_fitsigma[ip] = m_neg_eangle[ip]->GetFunction(
"Vavilov")->GetParameter(2);
305 neg_chisq[ip] = (m_neg_eangle[ip]->GetFunction(
"Vavilov")-> GetChisquare())/(m_neg_eangle[ip]->GetFunction(
"Vavilov")-> GetNDF());
309 fitmean[ip] = m_eangle[ip]->GetFunction(
"AsymGauss")->GetParameter(1);
310 fitmeanerr[ip] = m_eangle[ip]->GetFunction(
"AsymGauss")->GetParError(1);
311 fitsigma[ip] = m_eangle[ip]->GetFunction(
"AsymGauss")->GetParameter(2);
312 chisq[ip] = (m_eangle[ip]->GetFunction(
"AsymGauss")-> GetChisquare())/(m_eangle[ip]->GetFunction(
"AsymGauss")-> GetNDF());
314 pos_fitmean[ip] = m_pos_eangle[ip]->GetFunction(
"AsymGauss")->GetParameter(1);
315 pos_fitmeanerr[ip] = m_pos_eangle[ip]->GetFunction(
"AsymGauss")->GetParError(1);
316 pos_fitsigma[ip] = m_pos_eangle[ip]->GetFunction(
"AsymGauss")->GetParameter(2);
317 pos_chisq[ip] = (m_pos_eangle[ip]->GetFunction(
"AsymGauss")-> GetChisquare())/(m_pos_eangle[ip]->GetFunction(
"AsymGauss")-> GetNDF());
319 neg_fitmean[ip] = m_neg_eangle[ip]->GetFunction(
"AsymGauss")->GetParameter(1);
320 neg_fitmeanerr[ip] = m_neg_eangle[ip]->GetFunction(
"AsymGauss")->GetParError(1);
321 neg_fitsigma[ip] = m_neg_eangle[ip]->GetFunction(
"AsymGauss")->GetParameter(2);
322 neg_chisq[ip] = (m_neg_eangle[ip]->GetFunction(
"AsymGauss")-> GetChisquare())/(m_neg_eangle[ip]->GetFunction(
"AsymGauss")-> GetNDF());
324 if(entryNo[ip]>1500 && fitmean[ip]>0 && fitmeanerr[ip]>0 && fitmeanerr[ip]<10 && fitsigma[ip]>0 && fitsigma[ip]<200){
326 trunc_Norm += trunc_mean[ip];
331 trunc_Norm = trunc_Norm/
count;
333 cout<<
"count= "<<
count<<
" average value: "<< Norm <<
" truncted: " << trunc_Norm << endl;
336 if(entryNo[i]>10 && fitmean[i]>0 && fitsigma[i]>0) gain[i]=fitmean[i]/Norm;
338 if(entryNo[i]>10) trunc_gain[i]=trunc_mean[i]/trunc_Norm;
339 else trunc_mean[i] = 1;
340 if(pos_entryNo[i]>10 && pos_fitmean[i]>0 && pos_fitsigma[i]>0) pos_gain[i]=pos_fitmean[i]/Norm;
341 else pos_gain[i] = 1;
342 if(neg_entryNo[i]>10 && neg_fitmean[i]>0 && neg_fitsigma[i]>0) neg_gain[i]=neg_fitmean[i]/Norm;
343 else neg_gain[i] = 1;
345 m_EAngAverdE -> SetBinContent(i+1, gain[i]);
346 m_EAngAverdE->SetBinError(i+1, fitmeanerr[i]/Norm);
348 if(trunc_gain[i] != 1){
349 m_trunc_EAngAverdE -> SetBinContent(i+1, trunc_gain[i]);
350 m_trunc_EAngAverdE->SetBinError(i+1, trunc_sigma[i]/trunc_Norm);
352 if(pos_gain[i] != 1){
353 m_pos_EAngAverdE -> SetBinContent(i+1, pos_gain[i]);
354 m_pos_EAngAverdE->SetBinError(i+1, pos_fitmeanerr[i]/Norm);
356 if(neg_gain[i] != 1){
357 m_neg_EAngAverdE -> SetBinContent(i+1, neg_gain[i]);
358 m_neg_EAngAverdE->SetBinError(i+1, neg_fitmeanerr[i]/Norm);
362 log<<MSG::INFO<<
"begin getting calibration constant!!! "<<endreq;
364 int denangle_entry[1] = {100};
365 m_EAngAverdE->Fit(
"pol1",
"r",
"",-0.25,0.25);
366 double par0 = m_EAngAverdE->GetFunction(
"pol1")->GetParameter(0);
367 double par1 = m_EAngAverdE->GetFunction(
"pol1")->GetParameter(1);
371 for(
int i=0;i<
NumSlices;i++) denangle[i] = gain[i];
373 log<<MSG::INFO<<
"begin generating root file!!! "<<endreq;
374 TFile*
f =
new TFile(
m_rootfile.c_str(),
"recreate");
376 m_eangle[ip]->Write();
377 m_trunc_eangle[ip]->Write();
378 m_pos_eangle[ip]->Write();
379 m_neg_eangle[ip]->Write();
381 m_EAngAverdE->Write();
382 m_trunc_EAngAverdE->Write();
383 m_pos_EAngAverdE->Write();
384 m_neg_EAngAverdE->Write();
386 TTree* entracalib =
new TTree(
"entracalib",
"entracalib");
387 entracalib->Branch(
"1denangle", denangle,
"denangle[100]/D");
388 entracalib->Branch(
"1denangle_entry", denangle_entry,
"denangle_entry[1]/I");
390 entracalib->
Write();
392 TTree* ddgcalib =
new TTree(
"ddgcalib",
"ddgcalib");
393 ddgcalib ->
Branch(
"hits", entryNo,
"entryNo[100]/D");
394 ddgcalib ->
Branch(
"truncmean", trunc_mean,
"trunc_mean[100]/D");
395 ddgcalib ->
Branch(
"truncsigma", trunc_sigma,
"trunc_sigma[100]/D");
396 ddgcalib ->
Branch(
"fitmean", fitmean,
"fitmean[100]/D");
397 ddgcalib ->
Branch(
"fitmeanerr", fitmeanerr,
"fitmeanerr[100]/D");
398 ddgcalib ->
Branch(
"chi", fitsigma,
"fitsigma[100]/D");
399 ddgcalib ->
Branch(
"gain", gain,
"gain[100]/D");
400 ddgcalib ->
Branch(
"truncgain", trunc_gain,
"trunc_gain[100]/D");
401 ddgcalib ->
Branch(
"chisq", chisq,
"chisq[100]/D");
402 ddgcalib ->
Branch(
"pos_hits", pos_entryNo,
"pos_entryNo[100]/D");
403 ddgcalib ->
Branch(
"pos_mean", pos_mean,
"pos_mean[100]/D");
404 ddgcalib ->
Branch(
"pos_sigma", pos_sigma,
"pos_sigma[100]/D");
405 ddgcalib ->
Branch(
"pos_fitmean", pos_fitmean,
"pos_fitmean[100]/D");
406 ddgcalib ->
Branch(
"pos_fitmeanerr", pos_fitmeanerr,
"pos_fitmeanerr[100]/D");
407 ddgcalib ->
Branch(
"pos_chi", pos_fitsigma,
"pos_fitsigma[100]/D");
408 ddgcalib ->
Branch(
"pos_gain", pos_gain,
"pos_gain[100]/D");
409 ddgcalib ->
Branch(
"pos_chisq", pos_chisq,
"pos_chisq[100]/D");
410 ddgcalib ->
Branch(
"neg_hits", neg_entryNo,
"neg_entryNo[100]/D");
411 ddgcalib ->
Branch(
"neg_mean", neg_mean,
"neg_mean[100]/D");
412 ddgcalib ->
Branch(
"neg_sigma", neg_sigma,
"neg_sigma[100]/D");
413 ddgcalib ->
Branch(
"neg_fitmean", neg_fitmean,
"neg_fitmean[100]/D");
414 ddgcalib ->
Branch(
"neg_fitmeanerr", neg_fitmeanerr,
"neg_fitmeanerr[100]/D");
415 ddgcalib ->
Branch(
"neg_chi", neg_fitsigma,
"neg_fitsigma[100]/D");
416 ddgcalib ->
Branch(
"neg_gain", neg_gain,
"neg_gain[100]/D");
417 ddgcalib ->
Branch(
"neg_chisq", neg_chisq,
"neg_chisq[100]/D");
418 ddgcalib ->
Branch(
"Ip_eangle", Ip_eangle,
"Ip_eangle[100]/D");
419 ddgcalib ->
Branch(
"eangle", eangle,
"eangle[100]/D");
424 TCanvas c1(
"c1",
"canvas", 500, 400);
426 m_EAngAverdE->Draw();
428 m_trunc_EAngAverdE->Draw();
430 m_pos_EAngAverdE->Draw();
432 m_neg_EAngAverdE->Draw();
436 m_eangle[ip]->Draw();
441 m_trunc_eangle[ip]->Draw();
446 m_pos_eangle[ip]->Draw();
451 m_neg_eangle[ip]->Draw();
459 delete m_trunc_eangle[ip];
460 delete m_pos_eangle[ip];
461 delete m_neg_eangle[ip];
464 delete m_trunc_EAngAverdE;
465 delete m_pos_EAngAverdE;
466 delete m_neg_EAngAverdE;
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)
#define Iner_DriftDistCut
double landaun(double *x, double *par)
double AsymGauss(double *x, double *par)
double mylan(double *x, double *par)
DedxCalibEAng(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")