BOSS 7.0.5
BESIII Offline Software System
Loading...
Searching...
No Matches
HadronWidget/HadronCalibration.cc
Go to the documentation of this file.
1#include "HadronCalibration.h"
2
3HadronCalibration::HadronCalibration() {}
4
5void
6HadronCalibration::fitBGCurve( std::vector< TString > particles, TString filename, std::string paramfile ){
7
8 // read in a file that contains fit results for bg bins
9 TFile* infile = new TFile(filename);
10
11 // multigraphs to hold the curve and residual results
12 TMultiGraph* gr = new TMultiGraph("gr",";#beta#gamma;dE/dx");
13 TMultiGraph* gr_res = new TMultiGraph("gr_res",";dE/dx;#sigma");
14 TMultiGraph* grp = new TMultiGraph("grp",";Momentum;dE/dx");
15 TMultiGraph* grp_res = new TMultiGraph("grp_res",";Momentum;#sigma");
16
17 // multigraphs to hold the chi and sigma distributions
18 TMultiGraph* grchi = new TMultiGraph("grchi",";#beta#gamma;#chi mean");
19 TMultiGraph* grsigma = new TMultiGraph("grsigma",";#beta#gamma;#sigma");
20 TMultiGraph* grchip = new TMultiGraph("grchip",";Momentum;#chi mean");
21 TMultiGraph* grsigmap = new TMultiGraph("grsigmap",";Momentum;#sigma");
22
23 const int npart = 5;
24 TGraphErrors particle_dedx[npart];
25 TGraph particle_res[npart];
26 TGraphErrors particle_dedxp[npart];
27 TGraph particle_resp[npart];
28 TGraph newchi[npart];
29 TGraph newsigma[npart];
30 TGraph newchip[npart];
31 TGraph newsigmap[npart];
32
33 // --------------------------------------------------
34 // FILL BG CURVE VALUES
35 // --------------------------------------------------
36
37 double mass = 0.0;
38 for( int i = 0; i < particles.size(); ++i ){
39 TString particle = particles[i];
40
41 if( particle == "pion" ) mass = Widget::mpion;
42 else if( particle == "kaon" ) mass = Widget::mkaon;
43 else if( particle == "proton" ) mass = Widget::mproton;
44 else if( particle == "muon" ) mass = Widget::mmuon;
45 else if( particle == "electron" ) mass = Widget::melectron;
46 if( mass == 0.0 ) exit(1);
47
48 TTree* hadron = (TTree*)infile->Get(particle);
49
50 std::cout << "HadronCalibration: reading " << particle << " in file " << filename << std::endl;
51
52 double dedx_pub; // dE/dx without electron saturation correction
53 double dedx; // dE/dx with electron saturation correction
54 double bg; // track beta-gamma
55
56 double dedxerr; //
57 double dedx_res; //
58 // double nhit_avg; //
59 double sin_avg; //
60 double chi; //
61 double chi_pub; //
62 double sigma; //
63 double sigma_pub; //
64
65 hadron->SetBranchAddress("dedx",&dedx);
66 hadron->SetBranchAddress("dedx_pub",&dedx_pub);
67 hadron->SetBranchAddress("bg_avg",&bg);
68 hadron->SetBranchAddress("error",&dedxerr);
69 hadron->SetBranchAddress("dedxres_avg",&dedx_res);
70 // hadron->SetBranchAddress("nhit_avg",&nhit_avg);
71 hadron->SetBranchAddress("sinth_avg",&sin_avg);
72 hadron->SetBranchAddress("chi",&chi);
73 hadron->SetBranchAddress("chi_pub",&chi_pub);
74 hadron->SetBranchAddress("sigma",&sigma);
75 hadron->SetBranchAddress("sigma_pub",&sigma_pub);
76
77 for( int j = 0; j < hadron->GetEntries(); ++j ){
78 hadron->GetEvent(j);
79
80 particle_dedx[i].SetPoint(j,bg,dedx);
81 particle_dedx[i].SetPointError(j,0,dedxerr);
82 particle_res[i].SetPoint(j,dedx,dedx_res);
83
84 particle_dedxp[i].SetPoint(j,bg*mass,dedx);
85 particle_dedxp[i].SetPointError(j,0,dedxerr);
86 particle_resp[i].SetPoint(j,bg*mass,dedx_res);
87
88 newchi[i].SetPoint(j,bg,chi);
89 newsigma[i].SetPoint(j,bg,sigma);
90 newchip[i].SetPoint(j,bg*mass,chi);
91 newsigmap[i].SetPoint(j,bg*mass,sigma);
92 }
93
94 newchi[i].SetMarkerStyle(4);
95 newchi[i].SetMarkerColor(i+1);
96 if( i == 4 ) newchi[i].SetMarkerColor(i+2);
97 newchip[i].SetMarkerStyle(4);
98 newchip[i].SetMarkerColor(i+1);
99 if( i == 4 ) newchip[i].SetMarkerColor(i+2);
100
101 newsigma[i].SetMarkerStyle(4);
102 newsigma[i].SetMarkerColor(i+1);
103 if( i == 4 ) newsigma[i].SetMarkerColor(i+2);
104 newsigmap[i].SetMarkerStyle(4);
105 newsigmap[i].SetMarkerColor(i+1);
106 if( i == 4 ) newsigmap[i].SetMarkerColor(i+2);
107
108 particle_dedx[i].SetMarkerStyle(4);
109 particle_dedx[i].SetMarkerColor(i+1);
110 if( i == 4 ) particle_dedx[i].SetMarkerColor(i+2);
111 particle_dedxp[i].SetMarkerStyle(4);
112 particle_dedxp[i].SetMarkerColor(i+1);
113 if( i == 4 ) particle_dedxp[i].SetMarkerColor(i+2);
114
115 particle_res[i].SetMarkerStyle(4);
116 particle_res[i].SetMarkerColor(i+1);
117 if( i == 4 ) particle_res[i].SetMarkerColor(i+2);
118 particle_resp[i].SetMarkerStyle(4);
119 particle_resp[i].SetMarkerColor(i+1);
120 if( i == 4 ) particle_resp[i].SetMarkerColor(i+2);
121
122 gr->Add(&particle_dedx[i]);
123 gr_res->Add(&particle_res[i]);
124 grp->Add(&particle_dedxp[i]);
125 grp_res->Add(&particle_resp[i]);
126 grchi->Add(&newchi[i]);
127 grsigma->Add(&newsigma[i]);
128 grchip->Add(&newchip[i]);
129 grsigmap->Add(&newsigmap[i]);
130 }
131
132
133 // --------------------------------------------------
134 // FIT BG CURVE
135 // --------------------------------------------------
136
137 WidgetParameterization gpar(paramfile);
138 WidgetCurve* gc = new WidgetCurve();
139
140 double bgmin1 = 0.1, bgmax1 = 5;
141 double bgmin2 = 4, bgmax2 = 12;
142 double bgmin3 = 8, bgmax3 = 2000;
143
144 TF1* fdedx1 = new TF1("fdedx1",gc,bgmin1,bgmax1,9,"WidgetCurve");
145 fdedx1->SetParameter(0,1);
146 for( int i = 1; i < 9; ++i ){
147 fdedx1->SetParameter(i,gpar.getCurvePars(i-1));
148 }
149 // fix one of the redundant parameters
150 fdedx1->FixParameter(4,1);
151
152 TF1* fdedx2 = new TF1("fdedx2",gc,bgmin2,bgmax2,5,"WidgetCurve");
153 fdedx2->SetParameter(0,2);
154 for( int i = 1; i < 5; ++i ){
155 fdedx2->SetParameter(i,gpar.getCurvePars(7+i));
156 }
157
158 TF1* fdedx3 = new TF1("fdedx3",gc,bgmin3,bgmax3,5,"WidgetCurve");
159 fdedx3->SetParameter(0,3);
160 for( int i = 1; i < 5; ++i ){
161 fdedx3->SetParameter(i,gpar.getCurvePars(11+i));
162 }
163
164 TF1* fdedx4 = new TF1("fdedx4","1",100,100000);
165
166 TCanvas* bandcan = new TCanvas("bandcan","dE/dx",820,750);
167 grp->Draw("APE");
168 bandcan->SaveAs("plots/dedxbands.eps");
169 delete bandcan;
170
171 TCanvas* logbandcan = new TCanvas("logbandcan","dE/dx",820,750);
172 logbandcan->cd()->SetLogy();
173 grp->Draw("APE");
174 logbandcan->SaveAs("plots/dedxbands_log.eps");
175 delete logbandcan;
176
177 TCanvas* bgcurvecan = new TCanvas("bgcurvecan","dE/dx",820,750);
178 bgcurvecan->cd()->SetLogy();
179 bgcurvecan->cd()->SetLogx();
180 gr->Draw("APE");
181
182 int stat1 = gr->Fit("fdedx1","","",bgmin1,bgmax1);
183 for( int i = 0; i < 50; ++i ){
184 if( stat1 == 0 ) break;
185 stat1 = gr->Fit("fdedx1","","",bgmin1,bgmax1);
186 }
187 int stat2 = gr->Fit("fdedx2","","",bgmin2,bgmax2);
188 for( int i = 0; i < 50; ++i ){
189 if( stat2 == 0 ) break;
190 stat2 = gr->Fit("fdedx2","","",bgmin2,bgmax2);
191 }
192 int stat3 = gr->Fit("fdedx3","","",bgmin3,bgmax3);
193 for( int i = 0; i < 50; ++i ){
194 if( stat3 == 0 ) break;
195 stat3 = gr->Fit("fdedx3","","",bgmin3,bgmax3);
196 }
197
198 // if the fit was successful, write out the updated parameters
199 if( stat1 != 0 ) std::cout << "WARNING: BG FIT 1 FAILED..." << std::endl;
200 for( int i = 1; i < 9; ++i ){
201 gpar.setCurvePars(i-1,fdedx1->GetParameter(i));
202 }
203 if( stat2 != 0 ) std::cout << "WARNING: BG FIT 2 FAILED..." << std::endl;
204 for( int i = 1; i < 5; ++i ){
205 gpar.setCurvePars(7+i,fdedx2->GetParameter(i));
206 }
207 if( stat3 != 0 ) std::cout << "WARNING: BG FIT 3 FAILED..." << std::endl;
208 for( int i = 1; i < 5; ++i ){
209 gpar.setCurvePars(11+i,fdedx3->GetParameter(i));
210 }
211
212 fdedx1->SetLineColor(kBlack);
213 fdedx1->Draw("same");
214 fdedx2->SetLineColor(kBlue);
215 fdedx2->Draw("same");
216 fdedx3->SetLineColor(kRed);
217 fdedx3->Draw("same");
218 fdedx4->SetLineColor(kGray);
219 fdedx4->Draw("same");
220
221 TLegend* tleg = new TLegend(0.4,0.6,0.6,0.8);
222 for( int i = 0; i < 5; ++i ){
223 tleg->AddEntry(&particle_dedx[i],particles[i],"p");
224 }
225 tleg->Draw("same");
226
227 bgcurvecan->SaveAs("plots/bgcurve.eps");
228 delete bgcurvecan;
229
230
231 // --------------------------------------------------
232 // GET RESIDUALS AND CHIS
233 // --------------------------------------------------
234
235 double A = 4.5, B = 10;
236
237 TMultiGraph* fit_res = new TMultiGraph("fit_res",";#beta#gamma;Residual");
238 TGraph particle_residual[npart];
239 int respoint = 1;
240 double rmin = 1.0, rmax = 1.0;
241 for( int i = 0; i < npart; ++i ){
242 for( int j = 0; j < particle_dedx[i].GetN(); ++j ){
243 double x, y, fit;
244 particle_dedx[i].GetPoint(j,x,y);
245 if( y == 0 ) continue;
246 if( x < A )
247 fit = fdedx1->Eval(x);
248 else if( x < B )
249 fit = fdedx2->Eval(x);
250 else
251 fit = fdedx3->Eval(x);
252
253 // the curve is just 1 for electrons...
254 if( npart == 4 ) fit = 1.0;
255
256 particle_residual[i].SetPoint(respoint++,x,fit/y);
257 if( fit/y < rmin ) rmin = fit/y;
258 else if( fit/y > rmax ) rmax = fit/y;
259 }
260
261 particle_residual[i].SetMarkerStyle(4);
262 particle_residual[i].SetMarkerColor(i+1);
263 if( i == 4 ) particle_residual[i].SetMarkerColor(i+2);
264 fit_res->Add(&particle_residual[i]);
265 }
266 fit_res->SetMinimum(rmin*0.98);
267 fit_res->SetMaximum(rmax*1.02);
268
269 TCanvas* bgrescan = new TCanvas("bgrescan","dE/dx",820,750);
270 bgrescan->cd()->SetLogx();
271 fit_res->GetXaxis()->SetRange(0.05,5000);
272 fit_res->Draw("AP");
273 tleg->Draw("same");
274
275 bgrescan->SaveAs("plots/bgresidual.eps");
276 delete bgrescan;
277
278
279 // --------------------------------------------------
280 // FIT SIGMA VS DEDX
281 // --------------------------------------------------
282
283 TF1* sigvsdedx = new TF1("sigvsdedx","[0]+[1]*x",0.0,10.0);
284 sigvsdedx->SetParameter(0,gpar.getDedxPars(0));
285 sigvsdedx->SetParameter(1,gpar.getDedxPars(1));
286
287 TCanvas* sigcan = new TCanvas("sigcan","dE/dx",820,750);
288 sigcan->cd();
289 gr_res->Draw("APE");
290 tleg->Draw("same");
291
292 int status = gr_res->Fit("sigvsdedx","qm","",0.0,4.0);
293 for( int i = 0; i < 10; ++i ){
294 if( status == 0 ) break;
295 status = gr_res->Fit("sigvsdedx","q","",0.0,4.0);
296 }
297 sigcan->SaveAs("plots/sigmavsdedx.eps");
298 delete sigcan;
299
300 if( status != 0 ) std::cout << "WARNING: SIGMA VS DEDX FIT FAILED..." << std::endl;
301 gpar.setDedxPars(0,sigvsdedx->GetParameter(0));
302 gpar.setDedxPars(1,sigvsdedx->GetParameter(1));
303
304 // write out the (possibly) updated parameters to file
305 gpar.printParameters("parameters.bgcurve.fit");
306
307 sigvsdedx->Draw("same");
308 tleg->Draw("same");
309
310
311 // --------------------------------------------------
312 // PLOT CHI MEAN VS BETA-GAMMA
313 // --------------------------------------------------
314
315 TLine* line0 = new TLine(0,0,10000,0);
316 line0->SetLineStyle(kDashed);
317 line0->SetLineColor(kRed);
318 TLine* line1 = new TLine(0,1,10000,1);
319 line1->SetLineStyle(kDashed);
320 line1->SetLineColor(kRed);
321
322 TCanvas* cbgcan = new TCanvas("cbgcan","Mean #chi vs. #beta#gamma");
323 grchi->SetMinimum(-0.5);
324 grchi->SetMaximum(+0.5);
325 grchi->Draw("AP");
326 line0->Draw("same");
327 tleg->Draw("same");
328 cbgcan->SaveAs("plots/chimeanVSbetagamma.eps");
329
330 grchi->GetXaxis()->SetLimits(0,25);
331 grchi->Draw("AP");
332 line0->Draw("same");
333 tleg->Draw("same");
334 cbgcan->SaveAs("plots/chimeanVSbetagamma_zoom.eps");
335
336 // grsigma->GetXaxis()->SetTitle("#beta#gamma");
337 // grsigma->GetYaxis()->SetTitle("#chi mean");
338 grsigma->SetMinimum(0.5);
339 grsigma->SetMaximum(1.5);
340 grsigma->Draw("AP");
341 line1->Draw("same");
342 tleg->Draw("same");
343 cbgcan->SaveAs("plots/sigmaVSbetagamma.eps");
344
345 grsigma->GetXaxis()->SetLimits(0,25);
346 grsigma->Draw("AP");
347 line1->Draw("same");
348 tleg->Draw("same");
349 cbgcan->SaveAs("plots/sigmaVSbetagamma_zoom.eps");
350
351 TCanvas* cpcan = new TCanvas("cpcan","Mean #chi vs. momentum");
352 grchip->SetMinimum(-0.5);
353 grchip->SetMaximum(+0.5);
354 grchip->Draw("AP");
355 line0->Draw("same");
356 tleg->Draw("same");
357 cpcan->SaveAs("plots/chimeanVSmomentum.eps");
358
359 grsigmap->SetMinimum(0.5);
360 grsigmap->SetMaximum(1.5);
361 grsigmap->Draw("AP");
362 line1->Draw("same");
363 tleg->Draw("same");
364 cpcan->SaveAs("plots/sigmaVSmomentum.eps");
365
366 delete cbgcan;
367 delete cpcan;
368 delete line0;
369 delete line1;
370 delete gr;
371 delete gr_res;
372 delete grp;
373 delete grp_res;
374 delete grchi;
375 delete grsigma;
376 delete grchip;
377 delete grsigmap;
378}
379
380void
381HadronCalibration::plotBGCurve( std::vector< TString > particles, TString filename, std::string paramfile ){
382
383 // read in a file that contains fit results for bg bins
384 TFile* infile = new TFile(filename);
385
386 // multigraphs to hold the curve and residual results
387 TMultiGraph* gr = new TMultiGraph("gr",";#beta#gamma;dE/dx");
388 TMultiGraph* gr_res = new TMultiGraph("gr_res",";dE/dx;#sigma");
389 TMultiGraph* grp = new TMultiGraph("grp",";Momentum;dE/dx");
390 TMultiGraph* grp_res = new TMultiGraph("grp_res",";Momentum;#sigma");
391
392 // multigraphs to hold the chi and sigma distributions
393 TMultiGraph* grchi = new TMultiGraph("grchi",";#beta#gamma;#chi mean");
394 TMultiGraph* grsigma = new TMultiGraph("grsigma",";#beta#gamma;#sigma");
395 TMultiGraph* grchip = new TMultiGraph("grchip",";Momentum;#chi mean");
396 TMultiGraph* grsigmap = new TMultiGraph("grsigmap",";Momentum;#sigma");
397
398 const int npart = 5;
399 TGraphErrors particle_dedx[npart];
400 TGraph particle_res[npart];
401 TGraphErrors particle_dedxp[npart];
402 TGraph particle_resp[npart];
403 TGraph newchi[npart];
404 TGraph newsigma[npart];
405 TGraph newchip[npart];
406 TGraph newsigmap[npart];
407
408 // --------------------------------------------------
409 // FILL BG CURVE VALUES
410 // --------------------------------------------------
411
412 double mass = 0.0;
413 for( int i = 0; i < particles.size(); ++i ){
414 TString particle = particles[i];
415
416 if( particle == "pion" ) mass = Widget::mpion;
417 else if( particle == "kaon" ) mass = Widget::mkaon;
418 else if( particle == "proton" ) mass = Widget::mproton;
419 else if( particle == "muon" ) mass = Widget::mmuon;
420 else if( particle == "electron" ) mass = Widget::melectron;
421 if( mass == 0.0 ) exit(1);
422
423 TTree* hadron = (TTree*)infile->Get(particle);
424
425 std::cout << "HadronCalibration: reading " << particle << " in file " << filename << std::endl;
426
427 double dedx_pub; // dE/dx without electron saturation correction
428 double dedx; // dE/dx with electron saturation correction
429 double bg; // track beta-gamma
430
431 double dedxerr; //
432 double dedx_res; //
433 // double nhit_avg; //
434 double sin_avg; //
435 double chi; //
436 double chi_pub; //
437 double sigma; //
438 double sigma_pub; //
439
440 hadron->SetBranchAddress("dedx",&dedx);
441 hadron->SetBranchAddress("dedx_pub",&dedx_pub);
442 hadron->SetBranchAddress("bg_avg",&bg);
443 hadron->SetBranchAddress("error",&dedxerr);
444 hadron->SetBranchAddress("dedxres_avg",&dedx_res);
445 // hadron->SetBranchAddress("nhit_avg",&nhit_avg);
446 hadron->SetBranchAddress("sinth_avg",&sin_avg);
447 hadron->SetBranchAddress("chi",&chi);
448 hadron->SetBranchAddress("chi_pub",&chi_pub);
449 hadron->SetBranchAddress("sigma",&sigma);
450 hadron->SetBranchAddress("sigma_pub",&sigma_pub);
451
452 for( int j = 0; j < hadron->GetEntries(); ++j ){
453 hadron->GetEvent(j);
454
455 particle_dedx[i].SetPoint(j,bg,dedx);
456 particle_dedx[i].SetPointError(j,0,dedxerr);
457 particle_res[i].SetPoint(j,dedx,dedx_res);
458
459 particle_dedxp[i].SetPoint(j,bg*mass,dedx);
460 particle_dedxp[i].SetPointError(j,0,dedxerr);
461 particle_resp[i].SetPoint(j,bg*mass,dedx_res);
462
463 newchi[i].SetPoint(j,bg,chi);
464 newsigma[i].SetPoint(j,bg,sigma);
465 newchip[i].SetPoint(j,bg*mass,chi);
466 newsigmap[i].SetPoint(j,bg*mass,sigma);
467 }
468
469 newchi[i].SetMarkerStyle(4);
470 newchi[i].SetMarkerColor(i+1);
471 if( i == 4 ) newchi[i].SetMarkerColor(i+2);
472 newchip[i].SetMarkerStyle(4);
473 newchip[i].SetMarkerColor(i+1);
474 if( i == 4 ) newchip[i].SetMarkerColor(i+2);
475
476 newsigma[i].SetMarkerStyle(4);
477 newsigma[i].SetMarkerColor(i+1);
478 if( i == 4 ) newsigma[i].SetMarkerColor(i+2);
479 newsigmap[i].SetMarkerStyle(4);
480 newsigmap[i].SetMarkerColor(i+1);
481 if( i == 4 ) newsigmap[i].SetMarkerColor(i+2);
482
483 particle_dedx[i].SetMarkerStyle(4);
484 particle_dedx[i].SetMarkerColor(i+1);
485 if( i == 4 ) particle_dedx[i].SetMarkerColor(i+2);
486 particle_dedxp[i].SetMarkerStyle(4);
487 particle_dedxp[i].SetMarkerColor(i+1);
488 if( i == 4 ) particle_dedxp[i].SetMarkerColor(i+2);
489
490 particle_res[i].SetMarkerStyle(4);
491 particle_res[i].SetMarkerColor(i+1);
492 if( i == 4 ) particle_res[i].SetMarkerColor(i+2);
493 particle_resp[i].SetMarkerStyle(4);
494 particle_resp[i].SetMarkerColor(i+1);
495 if( i == 4 ) particle_resp[i].SetMarkerColor(i+2);
496
497 gr->Add(&particle_dedx[i]);
498 gr_res->Add(&particle_res[i]);
499 grp->Add(&particle_dedxp[i]);
500 grp_res->Add(&particle_resp[i]);
501 grchi->Add(&newchi[i]);
502 grsigma->Add(&newsigma[i]);
503 grchip->Add(&newchip[i]);
504 grsigmap->Add(&newsigmap[i]);
505 }
506
507 TCanvas* bandcan = new TCanvas("bandcan","dE/dx",820,750);
508 grp->Draw("APE");
509 bandcan->SaveAs("plots/dedxbands.eps");
510 delete bandcan;
511
512 TCanvas* logbandcan = new TCanvas("logbandcan","dE/dx",820,750);
513 logbandcan->cd()->SetLogy();
514 grp->Draw("APE");
515 logbandcan->SaveAs("plots/dedxbands_log.eps");
516 delete logbandcan;
517
518 TCanvas* bgcurvecan = new TCanvas("bgcurvecan","dE/dx",820,750);
519 bgcurvecan->cd()->SetLogy();
520 bgcurvecan->cd()->SetLogx();
521 gr->Draw("APE");
522
523 // --------------------------------------------------
524 // PLOT CHI MEAN VS BETA-GAMMA
525 // --------------------------------------------------
526
527 TLine* line0 = new TLine(0,0,10000,0);
528 line0->SetLineStyle(kDashed);
529 line0->SetLineColor(kRed);
530 TLine* line1 = new TLine(0,1,10000,1);
531 line1->SetLineStyle(kDashed);
532 line1->SetLineColor(kRed);
533
534 TCanvas* cbgcan = new TCanvas("cbgcan","Mean #chi vs. #beta#gamma");
535 grchi->SetMinimum(-0.5);
536 grchi->SetMaximum(+0.5);
537 grchi->Draw("AP");
538 line0->Draw("same");
539 TLegend* tleg = new TLegend(0.4,0.6,0.6,0.8);
540 for( int i = 0; i < 5; ++i ){
541 tleg->AddEntry(&particle_dedx[i],particles[i],"p");
542 }
543 tleg->Draw("same");
544 cbgcan->SaveAs("plots/chimeanVSbetagamma.eps");
545
546 grchi->GetXaxis()->SetLimits(0,25);
547 grchi->Draw("AP");
548 line0->Draw("same");
549 tleg->Draw("same");
550 cbgcan->SaveAs("plots/chimeanVSbetagamma_zoom.eps");
551
552 // grsigma->GetXaxis()->SetTitle("#beta#gamma");
553 // grsigma->GetYaxis()->SetTitle("#chi mean");
554 grsigma->SetMinimum(0.5);
555 grsigma->SetMaximum(1.5);
556 grsigma->Draw("AP");
557 line1->Draw("same");
558 tleg->Draw("same");
559 cbgcan->SaveAs("plots/sigmaVSbetagamma.eps");
560
561 grsigma->GetXaxis()->SetLimits(0,25);
562 grsigma->Draw("AP");
563 line1->Draw("same");
564 tleg->Draw("same");
565 cbgcan->SaveAs("plots/sigmaVSbetagamma_zoom.eps");
566
567 TCanvas* cpcan = new TCanvas("cpcan","Mean #chi vs. momentum");
568 grchip->SetMinimum(-0.5);
569 grchip->SetMaximum(+0.5);
570 grchip->Draw("AP");
571 line0->Draw("same");
572 tleg->Draw("same");
573 cpcan->SaveAs("plots/chimeanVSmomentum.eps");
574
575 grsigmap->SetMinimum(0.5);
576 grsigmap->SetMaximum(1.5);
577 grsigmap->Draw("AP");
578 line1->Draw("same");
579 tleg->Draw("same");
580 cpcan->SaveAs("plots/sigmaVSmomentum.eps");
581
582 TFile* outfile = new TFile("dedxbands.root","RECREATE");
583 outfile->cd();
584 gr->Write();
585 gr_res->Write();
586 grp->Write();
587 grp_res->Write();
588 grchi->Write();
589 grsigma->Write();
590 grchip->Write();
591 grsigmap->Write();
592 outfile->Close();
593
594 delete cbgcan;
595 delete cpcan;
596 delete line0;
597 delete line1;
598 delete gr;
599 delete gr_res;
600 delete grp;
601 delete grp_res;
602 delete grchi;
603 delete grsigma;
604 delete grchip;
605 delete grsigmap;
606}
607
608void
609HadronCalibration::fitSigmaVsNHit( TString filename, std::string paramfile, int mcflag = 0, int type = 0 ){
610
611 TFile* infile = new TFile(filename);
612 TTree* hadron = (TTree*)infile->Get("track");
613
614 std::cout << "HadronCalibration: reading in file " << filename << std::endl;
615
616 // --------------------------------------------------
617 // INITIALIZE CONTAINERS
618 // --------------------------------------------------
619
620 double dedxpub; // dE/dx without electron saturation correction
621 double dedx; // dE/dx with electron saturation correction
622 double p; // track momentum
623 double bg; // track beta-gamma
624 double costh; // cosine of track polar angle
625 double db; // the nearest distance to the IP in the xy plane
626 double dz; // the nearest distance to the IP in the z plane
627 double chiPi; // PID chi value
628 double eop; // energy over momentum in the calorimeter
629
630 // Belle II variables
631 int b2nhit; // number of hits on this track
632
633 // BES III variables
634 double b3nhit; // number of hits on this track
635
636 hadron->SetBranchAddress("dedx", &dedxpub);
637 hadron->SetBranchAddress("dedxsat", &dedx);
638 hadron->SetBranchAddress("pF", &p);
639 hadron->SetBranchAddress("costh", &costh);
640 hadron->SetBranchAddress("db", &db);
641 hadron->SetBranchAddress("dz", &dz);
642 hadron->SetBranchAddress("chiPi", &chiPi);
643 hadron->SetBranchAddress("eopst", &eop);
644
645 if( type == 0 )
646 hadron->SetBranchAddress("numGoodHits", &b3nhit);
647 else if( type == 1 )
648 hadron->SetBranchAddress("lNHitsUsed", &b2nhit);
649 // hadron->SetBranchAddress("numGoodLayerHits", &b2nhit);
650
651 const int m_nhitbins = 20;
652 const double m_uppernhit = 32.5, m_lowernhit = 5.5;
653
654 double nhitstep = (m_uppernhit-m_lowernhit)/m_nhitbins;
655
656 // Create the histograms to be fit
657 TH1F* hdedx_nhit[m_nhitbins];
658 TH1F* hdedx_nhit_pub[m_nhitbins];
659
660 // initialize the histograms
661 for( int i = 0; i < m_nhitbins; ++i ){
662 char histname[100],histname2[100],histname3[100],histname4[100],histname5[100],histname6[100];
663 sprintf(histname, "dedx_nhit_%d",i);
664 sprintf(histname4, "dedxpub_nhit_%d",i);
665
666 hdedx_nhit[i] = new TH1F(histname,histname,200,-1,1);
667 hdedx_nhit_pub[i] = new TH1F(histname4,histname4,200,-1,1);
668 }
669
670 // Create some containers to calculate averages
671 double sumnhit[m_nhitbins];
672 double sumbg[m_nhitbins];
673 double sumsin[m_nhitbins];
674 int sumsize[m_nhitbins];
675 for( int i = 0; i < m_nhitbins; ++i ){
676 sumnhit[i] = 0;
677 sumbg[i] = 0;
678 sumsin[i] = 0;
679 sumsize[i] = 0;
680 }
681
682
683 // get the hadron saturation parameters
684 // if the parameters do not exist, use the values in the default constructor
685 double hadpar[5];
686 std::ifstream parfile("sat-pars.txt");
687 if( !parfile.fail() ){
688 for( int i = 0; i < 5; ++i ){
689 parfile >> hadpar[i];
690 }
691 parfile.close();
692 m_hadsat.setParameters(hadpar);
693 }
694 WidgetParameterization gpar(paramfile);
695
696 // --------------------------------------------------
697 // LOOP OVER EVENTS AND FILL CONTAINERS
698 // --------------------------------------------------
699
700 // Fill the histograms to be fitted
701 for( unsigned int index = 0; index < hadron->GetEntries(); ++index ){
702 hadron->GetEvent(index);
703
704 bg = fabs(p)/Widget::mpion;
705
706 int nhit;
707 if( type == 0 )
708 nhit = std::floor(b3nhit);
709 else if( type == 1 )
710 nhit = b2nhit;
711
712 if( dedx <= 0 || nhit <= m_lowernhit || nhit >= m_uppernhit )
713 continue;
714
715 int nhitBin = (int)((nhit-m_lowernhit)/(m_uppernhit-m_lowernhit) * m_nhitbins);
716
717 double dedx_new = dedx;
718 if( !mcflag ) dedx_new = m_hadsat.D2I(costh,m_hadsat.I2D(costh,1.0)*dedx);
719 double dedx_cur = gpar.dedxPrediction(bg);
720 double dedx_res = gpar.sigmaPrediction(dedx,nhit,sqrt(1-costh*costh));
721 double chi_new = (dedx_new-dedx_cur)/dedx_res;
722
723 double res_cor = gpar.sinPrediction(sqrt(1-costh*costh));
724 int i = (int) ( (fabs(nhit)-m_lowernhit)/nhitstep );
725 hdedx_nhit[i]->Fill((dedx_new-dedx_cur)/res_cor);
726 hdedx_nhit_pub[i]->Fill((dedx_new-dedx_cur)/res_cor);
727
728 if( fabs(nhit-(m_lowernhit+0.5*nhitstep+nhitBin*nhitstep)) < 0.5*nhitstep ){
729 sumnhit[nhitBin] += nhit;
730 sumbg[nhitBin] += bg;
731 sumsin[nhitBin] += sqrt(1-costh*costh);
732 sumsize[nhitBin] += 1;
733 }
734 }// end of event loop
735
736
737 // --------------------------------------------------
738 // FIT IN BINS OF NHIT
739 // --------------------------------------------------
740
741 // fit the histograms with Gaussian functions
742 // and extract the means and errors
743 TTree* nhitTree = new TTree("sigmavsnhit","dE/dx means and errors");
744
745 double nhitdedx; // dE/dx mean value for this bin
746 double nhitsigma; // dE/dx width value for this bin
747 double nhitdedxpub; // dE/dx mean value for this bin
748 double nhitsigmapub; // dE/dx width value for this bin
749 double nhitavg; // average nhit value for this bin
750 double nhitbgavg; // average bg value for this bin
751 double nhitsinavg; // average sin(theta) value for this bin
752
753 nhitTree->Branch("dedx",&nhitdedx,"dedx/D");
754 nhitTree->Branch("width",&nhitsigma,"width/D");
755 nhitTree->Branch("dedx_pub",&nhitdedxpub,"dedx_pub/D");
756 nhitTree->Branch("sigma_pub",&nhitsigmapub,"sigma_pub/D");
757 nhitTree->Branch("nhit_avg",&nhitavg,"nhit_avg/D");
758 nhitTree->Branch("bg_avg",&nhitbgavg,"bg_avg/D");
759 nhitTree->Branch("sin_avg",&nhitsinavg,"sin_avg/D");
760
761 double nhits[m_nhitbins], nhitserr[m_nhitbins], nhitres[m_nhitbins], nhitreserr[m_nhitbins];
762
763 // Fit the histograms
764 double avg_sigma = 0.0;
765 for( int i = 0; i < m_nhitbins; ++i ){
766
767 // fill some details for this bin
768 nhitavg = sumnhit[i]/sumsize[i];
769 nhitbgavg = sumbg[i]/sumsize[i];
770 nhitsinavg = sumsin[i]/sumsize[i];
771
772 // fit the dE/dx distribution in bins of nhit
773 int fs = hdedx_nhit[i]->Fit("gaus","ql");
774 if( fs != 0 ) std::cout << "\t\t" <<"MEAN FIT STATUS " << fs << std::endl;
775 double mean = hdedx_nhit[i]->GetFunction("gaus")->GetParameter(1);
776 double width = hdedx_nhit[i]->GetFunction("gaus")->GetParameter(2);
777 fs = hdedx_nhit[i]->Fit("gaus","ql","",mean-2.5*width,mean+2.5*width);
778 if( fs != 0 ) std::cout << "\t\t" <<"MEAN FIT STATUS " << fs << std::endl;
779
780 nhitdedx = hdedx_nhit[i]->GetFunction("gaus")->GetParameter(1);
781 nhitsigma = hdedx_nhit[i]->GetFunction("gaus")->GetParameter(2);
782
783
784 // fit the dE/dx distribution in bins of nhit
785 fs = hdedx_nhit_pub[i]->Fit("gaus","ql");
786 if( fs != 0 ) std::cout << "\t\t" <<"MEAN FIT STATUS " << fs << std::endl;
787 mean = hdedx_nhit_pub[i]->GetFunction("gaus")->GetParameter(1);
788 width = hdedx_nhit_pub[i]->GetFunction("gaus")->GetParameter(2);
789 fs = hdedx_nhit_pub[i]->Fit("gaus","ql","",mean-2.5*width,mean+2.5*width);
790 if( fs != 0 ) std::cout << "\t\t" <<"MEAN FIT STATUS " << fs << std::endl;
791
792 nhitdedxpub = hdedx_nhit_pub[i]->GetFunction("gaus")->GetParameter(1);
793 nhitsigmapub = hdedx_nhit_pub[i]->GetFunction("gaus")->GetParameter(2);
794
795
796 // same some information for fitting sigma vs nhit
797 nhits[i] = nhitavg;
798 nhitserr[i] = 0.0;
799 nhitres[i] = nhitsigma;
800 nhitreserr[i] = hdedx_nhit[i]->GetFunction("gaus")->GetParError(2);
801
802 avg_sigma += nhitres[i];
803
804 // fill the tree for this bin
805 nhitTree->Fill();
806 }
807
808 avg_sigma = avg_sigma/m_nhitbins;
809 for( int i = 0; i < m_nhitbins; ++i ){
810 nhitres[i] = nhitres[i]/avg_sigma;
811 nhitreserr[i] = nhitreserr[i]/avg_sigma;
812 }
813
814 // Print the histograms for quality control
815 TCanvas* ctmp = new TCanvas("tmp","tmp",900,900);
816 ctmp->Divide(3,3);
817 std::stringstream psname; psname << "plots/sigmavsnhit_fits.ps[";
818 ctmp->Print(psname.str().c_str());
819 psname.str(""); psname << "plots/sigmavsnhit_fits.ps";
820 for( int i = 0 ; i < m_nhitbins; ++i ){
821 ctmp->cd(i%9+1);
822 hdedx_nhit[i]->Draw();
823 if((i+1)%9==0)
824 ctmp->Print(psname.str().c_str());
825 }
826 psname.str(""); psname << "plots/sigmavsnhit_fits.ps]";
827 ctmp->Print(psname.str().c_str());
828 delete ctmp;
829
830
831 // --------------------------------------------------
832 // FIT SIGMA VS NHIT CURVE
833 // --------------------------------------------------
834
835
836 TCanvas* sigvsnhitcan = new TCanvas("sigvsnhitcan","#sigma vs. nHit",600,600);
837
838 TGraphErrors* gr = new TGraphErrors(m_nhitbins,nhits,nhitres,nhitserr,nhitreserr);
839 gr->SetMarkerStyle(8);
840 gr->SetMarkerSize(0.3);
841 gr->SetTitle(";nHit;#sigma");
842 /*
843 // Set stat options
844 gStyle->SetOptStat(1111111);
845 // Set y-position (fraction of pad size)
846 gStyle->SetStatY(0.65);
847 // Set x-position (fraction of pad size)
848 gStyle->SetStatX(0.9);
849 // Set width of stat-box (fraction of pad size)
850 gStyle->SetStatW(0.15);
851 // Set height of stat-box (fraction of pad size)
852 gStyle->SetStatH(0.15);
853 */
854 gr->Draw("AP");
855
856 WidgetSigma* gs = new WidgetSigma();
857
858 double nhitmin = 2, nhitmax = 50;
859 TF1* fsigma = new TF1("fsigma",gs,nhitmin,nhitmax,6,"WidgetSigma");
860 fsigma->SetParameter(0,2);
861 for( int i = 1; i < 6; ++i ){
862 fsigma->SetParameter(i,gpar.getNHitPars(i-1));
863 }
864
865 // if the fit succeeds, write out the new parameters
866 int status = gr->Fit("fsigma","qm","",nhitmin,nhitmax);
867 for( int i = 0; i < 10; ++i ){
868 if( status == 0 ) break;
869 status = gr->Fit("fsigma","q","",nhitmin,nhitmax);
870 }
871
872 if( status != 0 ) std::cout << "WARNING: SIGMA VS NHIT FIT FAILED..." << std::endl;
873 for( int j = 1; j < 6; ++j ){
874 gpar.setNHitPars(j-1,fsigma->GetParameter(j));
875 }
876
877 fsigma->Draw("same");
878 sigvsnhitcan->SaveAs("plots/sigmavsnhit.eps");
879 delete sigvsnhitcan;
880
881 // write out the (possibly) updated parameters to file
882 gpar.printParameters("parameters.sigmanhit.fit");
883}
884
885void
886HadronCalibration::fitSigmaVsSin( TString filename, std::string paramfile, int mcflag = 0, int type = 0 ){
887
888 TFile* infile = new TFile(filename);
889 TTree* hadron = (TTree*)infile->Get("track");
890
891 std::cout << "HadronCalibration: reading in file " << filename << std::endl;
892
893 // --------------------------------------------------
894 // INITIALIZE CONTAINERS
895 // --------------------------------------------------
896
897 double dedxpub; // dE/dx without electron saturation correction
898 double dedx; // dE/dx with electron saturation correction
899 double p; // track momentum
900 double bg; // track beta-gamma
901 double costh; // cosine of track polar angle
902 double db; // the nearest distance to the IP in the xy plane
903 double dz; // the nearest distance to the IP in the z plane
904 double chiPi; // PID chi value
905 double eop; // energy over momentum in the calorimeter
906
907 // Belle II variables
908 int b2nhit; // number of hits on this track
909
910 // BES III variables
911 double b3nhit; // number of hits on this track
912
913 hadron->SetBranchAddress("dedx", &dedxpub);
914 hadron->SetBranchAddress("dedxsat", &dedx);
915 hadron->SetBranchAddress("pF", &p);
916 hadron->SetBranchAddress("costh", &costh);
917 hadron->SetBranchAddress("db", &db);
918 hadron->SetBranchAddress("dz", &dz);
919 hadron->SetBranchAddress("chiPi", &chiPi);
920 hadron->SetBranchAddress("eopst", &eop);
921
922 if( type == 0 )
923 hadron->SetBranchAddress("numGoodHits", &b3nhit);
924 else if( type == 1 )
925 hadron->SetBranchAddress("lNHitsUsed", &b2nhit);
926 // hadron->SetBranchAddress("numGoodLayerHits", &b2nhit);
927
928 const int m_sinbins = 20;
929 const double m_uppersin = 1.0, m_lowersin = 0.36;
930
931 double sinstep = (m_uppersin-m_lowersin)/m_sinbins;
932
933 // Create the histograms to be fit
934 TH1F* hdedx_sin[m_sinbins];
935 TH1F* hdedx_sin_pub[m_sinbins];
936
937 // initialize the histograms
938 for( int i = 0; i < m_sinbins; ++i ){
939 char histname[100],histname2[100],histname3[100],histname4[100],histname5[100],histname6[100];
940 sprintf(histname, "dedx_sin_%d",i);
941 sprintf(histname4, "dedxpub_sin_%d",i);
942
943 hdedx_sin[i] = new TH1F(histname,histname,200,-1,1);
944 hdedx_sin_pub[i] = new TH1F(histname4,histname4,200,-1,1);
945 }
946
947 // Create some containers to calculate averages
948 double sumnhit[m_sinbins];
949 double sumbg[m_sinbins];
950 double sumsin[m_sinbins];
951 int sumsize[m_sinbins];
952 for( int i = 0; i < m_sinbins; ++i ){
953 sumnhit[i] = 0;
954 sumbg[i] = 0;
955 sumsin[i] = 0;
956 sumsize[i] = 0;
957 }
958
959
960 // get the hadron saturation parameters
961 // if the parameters do not exist, use the values in the default constructor
962 double hadpar[5];
963 std::ifstream parfile("sat-pars.txt");
964 if( !parfile.fail() ){
965 for( int i = 0; i < 5; ++i ){
966 parfile >> hadpar[i];
967 }
968 parfile.close();
969 m_hadsat.setParameters(hadpar);
970 }
971 WidgetParameterization gpar(paramfile);
972
973 // --------------------------------------------------
974 // LOOP OVER EVENTS AND FILL CONTAINERS
975 // --------------------------------------------------
976
977 // Fill the histograms to be fitted
978 for( unsigned int index = 0; index < hadron->GetEntries(); ++index ){
979 hadron->GetEvent(index);
980
981 if( costh != costh ) continue;
982
983 bg = fabs(p)/Widget::mpion;
984 double sin = sqrt(1-costh*costh);
985
986 int nhit;
987 if( type == 0 )
988 nhit = std::floor(b3nhit);
989 else if( type == 1 )
990 nhit = b2nhit;
991
992 if( dedx <= 0 || nhit <= 0 || nhit >= 100 || sin > m_uppersin || sin < m_lowersin )
993 continue;
994
995 int sinBin = (int)((sin-m_lowersin)/(m_uppersin-m_lowersin) * m_sinbins);
996
997 double dedx_new = dedx;
998 if( !mcflag ) dedx_new = m_hadsat.D2I(costh,m_hadsat.I2D(costh,1.0)*dedx);
999 double dedx_cur = gpar.dedxPrediction(bg);
1000 double dedx_res = gpar.sigmaPrediction(dedx,nhit,sqrt(1-costh*costh));
1001 double chi_new = (dedx_new-dedx_cur)/dedx_res;
1002
1003 double res_cor = gpar.nhitPrediction(nhit);
1004 int i = (int) ( (sin-m_lowersin)/sinstep );
1005 hdedx_sin[i]->Fill((dedx_new-dedx_cur)/res_cor);
1006 hdedx_sin_pub[i]->Fill((dedx_new-dedx_cur)/res_cor);
1007
1008 if( fabs(sin-(m_lowersin+0.5*sinstep+sinBin*sinstep)) < 0.5*sinstep ){
1009 sumnhit[sinBin] += nhit;
1010 sumbg[sinBin] += bg;
1011 sumsin[sinBin] += sqrt(1-costh*costh);
1012 sumsize[sinBin] += 1;
1013 }
1014 }// end of event loop
1015
1016
1017 // --------------------------------------------------
1018 // FIT IN BINS OF SIN
1019 // --------------------------------------------------
1020
1021 // fit the histograms with Gaussian functions
1022 // and extract the means and errors
1023 TTree* sinTree = new TTree("sigmavssin","dE/dx means and errors");
1024
1025 double sindedx; // dE/dx mean value for this bin
1026 double sinsigma; // dE/dx width value for this bin
1027 double sindedxpub; // dE/dx mean value for this bin
1028 double sinsigmapub; // dE/dx width value for this bin
1029 double sinnhitavg; // average nhit value for this bin
1030 double sinbgavg; // average bg value for this bin
1031 double sinavg; // average sin(theta) value for this bin
1032
1033 sinTree->Branch("dedx",&sindedx,"dedx/D");
1034 sinTree->Branch("width",&sinsigma,"width/D");
1035 sinTree->Branch("dedx_pub",&sindedxpub,"dedx_pub/D");
1036 sinTree->Branch("sigma_pub",&sinsigmapub,"sigma_pub/D");
1037 sinTree->Branch("nhit_avg",&sinnhitavg,"nhit_avg/D");
1038 sinTree->Branch("bg_avg",&sinbgavg,"bg_avg/D");
1039 sinTree->Branch("sin_avg",&sinavg,"sin_avg/D");
1040
1041 double sins[m_sinbins], sinserr[m_sinbins], sinres[m_sinbins], sinreserr[m_sinbins];
1042
1043 // Fit the histograms
1044 double avg_sigma = 0.0;
1045 for( int i = 0; i < m_sinbins; ++i ){
1046
1047 // fill some details for this bin
1048 sinnhitavg = sumnhit[i]/sumsize[i];
1049 sinbgavg = sumbg[i]/sumsize[i];
1050 sinavg = sumsin[i]/sumsize[i];
1051
1052
1053 // fit the dE/dx distribution in bins of sin
1054 int fs = hdedx_sin[i]->Fit("gaus","ql");
1055 if( fs != 0 ) std::cout << "\t\t" <<"MEAN FIT STATUS " << fs << std::endl;
1056 double mean = hdedx_sin[i]->GetFunction("gaus")->GetParameter(1);
1057 double width = hdedx_sin[i]->GetFunction("gaus")->GetParameter(2);
1058 fs = hdedx_sin[i]->Fit("gaus","ql","",mean-2.5*width,mean+2.5*width);
1059 if( fs != 0 ) std::cout << "\t\t" <<"MEAN FIT STATUS " << fs << std::endl;
1060
1061 sindedx = hdedx_sin[i]->GetFunction("gaus")->GetParameter(1);
1062 sinsigma = hdedx_sin[i]->GetFunction("gaus")->GetParameter(2);
1063
1064
1065 // fit the dE/dx distribution in bins of sin
1066 fs = hdedx_sin_pub[i]->Fit("gaus","ql");
1067 if( fs != 0 ) std::cout << "\t\t" <<"MEAN FIT STATUS " << fs << std::endl;
1068 mean = hdedx_sin_pub[i]->GetFunction("gaus")->GetParameter(1);
1069 width = hdedx_sin_pub[i]->GetFunction("gaus")->GetParameter(2);
1070 fs = hdedx_sin_pub[i]->Fit("gaus","ql","",mean-2.5*width,mean+2.5*width);
1071 if( fs != 0 ) std::cout << "\t\t" <<"MEAN FIT STATUS " << fs << std::endl;
1072
1073 sindedxpub = hdedx_sin_pub[i]->GetFunction("gaus")->GetParameter(1);
1074 sinsigmapub = hdedx_sin_pub[i]->GetFunction("gaus")->GetParameter(2);
1075
1076
1077 // same some information for fitting sigma vs sin
1078 sins[i] = sinavg;
1079 sinserr[i] = 0.0;
1080 sinres[i] = sinsigma;
1081 sinreserr[i] = hdedx_sin[i]->GetFunction("gaus")->GetParError(2);
1082
1083 avg_sigma += sinres[i];
1084
1085 // fill the tree for this bin
1086 sinTree->Fill();
1087 }
1088
1089 avg_sigma = avg_sigma/m_sinbins;
1090 for( int i = 0; i < m_sinbins; ++i ){
1091 sinres[i] = sinres[i]/avg_sigma;
1092 sinreserr[i] = sinreserr[i]/avg_sigma;
1093 }
1094
1095 // Print the histograms for quality control
1096 TCanvas* ctmp = new TCanvas("tmp","tmp",900,900);
1097 ctmp->Divide(3,3);
1098 std::stringstream psname; psname << "plots/sigmavssin_fits.ps[";
1099 ctmp->Print(psname.str().c_str());
1100 psname.str(""); psname << "plots/sigmavssin_fits.ps";
1101 for( int i = 0 ; i < m_sinbins; ++i ){
1102 ctmp->cd(i%9+1);
1103 hdedx_sin[i]->Draw();
1104 if((i+1)%9==0)
1105 ctmp->Print(psname.str().c_str());
1106 }
1107 psname.str(""); psname << "plots/sigmavssin_fits.ps]";
1108 ctmp->Print(psname.str().c_str());
1109 delete ctmp;
1110
1111
1112 // --------------------------------------------------
1113 // FIT SIGMA VS SIN CURVE
1114 // --------------------------------------------------
1115
1116
1117 TCanvas* sigvssincan = new TCanvas("sigvssincan","#sigma vs. sin(#theta)",600,600);
1118
1119 TGraphErrors* gr = new TGraphErrors(m_sinbins,sins,sinres,sinserr,sinreserr);
1120 gr->SetMarkerStyle(8);
1121 gr->SetMarkerSize(0.3);
1122 gr->SetTitle(";sin(#theta);#sigma");
1123 gr->Draw("AP");
1124
1125 WidgetSigma* gs = new WidgetSigma();
1126
1127 TF1* fsigma = new TF1("fsigma",gs,m_lowersin,m_uppersin,6,"WidgetSigma");
1128 fsigma->SetParameter(0,2);
1129 for( int i = 1; i < 6; ++i ){
1130 fsigma->SetParameter(i,gpar.getSinPars(i-1));
1131 }
1132
1133 // if the fit succeeds, write out the new parameters
1134 int status = gr->Fit("fsigma","q","",m_lowersin,m_uppersin);
1135 for( int i = 0; i < 10; ++i ){
1136 if( status == 0 ) break;
1137 status = gr->Fit("fsigma","q","",m_lowersin,m_uppersin);
1138 }
1139
1140 if( status != 0 ) std::cout << "WARNING: SIGMA VS SIN FIT FAILED..." << std::endl;
1141 for( int j = 1; j < 6; ++j ){
1142 gpar.setSinPars(j-1,fsigma->GetParameter(j));
1143 }
1144
1145 fsigma->Draw("same");
1146 sigvssincan->SaveAs("plots/sigmavssin.eps");
1147 delete sigvssincan;
1148
1149 // write out the (possibly) updated parameters to file
1150 gpar.printParameters("parameters.sigmasin.fit");
1151}
1152
1153
1154
1155
1156void
1157HadronCalibration::plotEfficiency( TString filenames[5], TString saveas, std::string paramfile, int mcflag = 0, int type = 0 ){
1158
1159 TFile* pifile = new TFile(filenames[0]);
1160 TFile* kfile = new TFile(filenames[1]);
1161 TTree* pion = (TTree*)pifile->Get("track");
1162 TTree* kaon = (TTree*)kfile->Get("track");
1163
1164 std::cout << "HadronCalibration: reading in file " << filenames[0] << std::endl;
1165
1166 // --------------------------------------------------
1167 // INITIALIZE CONTAINERS
1168 // --------------------------------------------------
1169
1170 double dedxpub; // dE/dx without electron saturation correction
1171 double dedx; // dE/dx with electron saturation correction
1172 double p; // track momentum
1173 double bg; // track beta-gamma
1174 double costh; // cosine of track polar angle
1175 double db; // the nearest distance to the IP in the xy plane
1176 double dz; // the nearest distance to the IP in the z plane
1177 double chiPi; // PID chi value
1178 double eop; // energy over momentum in the calorimeter
1179
1180 // Belle II variables
1181 int b2nhit; // number of hits on this track
1182
1183 // BES III variables
1184 double b3nhit; // number of hits on this track
1185
1186 pion->SetBranchAddress("dedx", &dedxpub);
1187 pion->SetBranchAddress("dedxsat", &dedx);
1188 pion->SetBranchAddress("pF", &p);
1189 pion->SetBranchAddress("costh", &costh);
1190 pion->SetBranchAddress("db", &db);
1191 pion->SetBranchAddress("dz", &dz);
1192 pion->SetBranchAddress("chiPi", &chiPi);
1193 pion->SetBranchAddress("eopst", &eop);
1194
1195 if( type == 0 )
1196 pion->SetBranchAddress("numGoodHits", &b3nhit);
1197 else if( type == 1 )
1198 pion->SetBranchAddress("lNHitsUsed", &b2nhit);
1199 // pion->SetBranchAddress("numGoodLayerHits", &b2nhit);
1200
1201
1202 kaon->SetBranchAddress("dedx", &dedxpub);
1203 kaon->SetBranchAddress("dedxsat", &dedx);
1204 kaon->SetBranchAddress("pF", &p);
1205 kaon->SetBranchAddress("costh", &costh);
1206 kaon->SetBranchAddress("db", &db);
1207 kaon->SetBranchAddress("dz", &dz);
1208 kaon->SetBranchAddress("chiPi", &chiPi);
1209 kaon->SetBranchAddress("eopst", &eop);
1210
1211 if( type == 0 )
1212 kaon->SetBranchAddress("numGoodHits", &b3nhit);
1213 else if( type == 1 )
1214 kaon->SetBranchAddress("lNHitsUsed", &b2nhit);
1215 // kaon->SetBranchAddress("numGoodLayerHits", &b2nhit);
1216
1217 const int m_nhitbins = 20;
1218 const double m_uppernhit = 32.5, m_lowernhit = 5.5;
1219
1220 double nhitstep = (m_uppernhit-m_lowernhit)/m_nhitbins;
1221
1222
1223 // get the hadron saturation parameters
1224 // if the parameters do not exist, use the values in the default constructor
1225 double hadpar[5];
1226 std::ifstream parfile("sat-pars.txt");
1227 if( !parfile.fail() ){
1228 for( int i = 0; i < 5; ++i ){
1229 parfile >> hadpar[i];
1230 }
1231 parfile.close();
1232 m_hadsat.setParameters(hadpar);
1233 }
1234 else
1235 std::cout << "WARNING: NO SATUARTION PARAMETERS!!!" << std::endl;
1236
1237 WidgetParameterization gpar(paramfile);
1238
1239 const int nbins = 20;
1240 TFile* outfile = new TFile(saveas,"RECREATE");
1241
1242 TH1F* pinumer = new TH1F("pinumer","",nbins,0,2);
1243 TH1F* pidenom = new TH1F("pidenom","",nbins,0,2);
1244
1245 TH2F* pihdedx = new TH2F("pihdedx","",100,0,2,100,0,4);
1246 TH2F* pipredpi = new TH2F("pipredpi","",100,0,2,100,0,4);
1247 TH2F* pipredk = new TH2F("pipredk","",100,0,2,100,0,4);
1248
1249 TH2F* pichivp = new TH2F("pichivp","",100,0,2,100,-10,10);
1250
1251 TH1F* knumer = new TH1F("knumer","",nbins,0,2);
1252 TH1F* kdenom = new TH1F("kdenom","",nbins,0,2);
1253
1254 TH2F* khdedx = new TH2F("khdedx","",100,0,2,100,0,4);
1255 TH2F* kpredpi = new TH2F("kpredpi","",100,0,2,100,0,4);
1256 TH2F* kpredk = new TH2F("kpredk","",100,0,2,100,0,4);
1257
1258 TH2F* kchivp = new TH2F("kchivp","",100,0,2,100,-10,10);
1259
1260 // --------------------------------------------------
1261 // LOOP OVER EVENTS AND FILL CONTAINERS
1262 // --------------------------------------------------
1263
1264 // Fill the pions
1265 for( unsigned int index = 0; index < pion->GetEntries(); ++index ){
1266 pion->GetEvent(index);
1267
1268 int nhit;
1269 if( type == 0 )
1270 nhit = std::floor(b3nhit);
1271 else if( type == 1 )
1272 nhit = b2nhit;
1273
1274 if( dedx <= 0 || nhit <= m_lowernhit || nhit >= m_uppernhit )
1275 continue;
1276
1277 double dedx_new = dedx;
1278 if( !mcflag ) dedx_new = m_hadsat.D2I(costh,m_hadsat.I2D(costh,1.0)*dedx);
1279 double dedx_res = gpar.sigmaPrediction(dedx,nhit,sqrt(1-costh*costh));
1280
1281 pihdedx->Fill(fabs(p),dedx_new);
1282
1283 double pibg = fabs(p)/Widget::mpion;
1284 double dedx_cur_pi = gpar.dedxPrediction(pibg);
1285 double chi_new_pi = (dedx_new-dedx_cur_pi)/dedx_res;
1286
1287 pipredpi->Fill(fabs(p),dedx_cur_pi);
1288
1289 double kbg = fabs(p)/Widget::mkaon;
1290 double dedx_cur_k = gpar.dedxPrediction(kbg);
1291 double chi_new_k = (dedx_new-dedx_cur_k)/dedx_res;
1292
1293 pipredk->Fill(fabs(p),dedx_cur_k);
1294
1295 pichivp->Fill(fabs(p),(chi_new_k*chi_new_k-chi_new_pi*chi_new_pi));
1296
1297 pidenom->Fill(fabs(p));
1298 if( chi_new_k*chi_new_k > chi_new_pi*chi_new_pi ) pinumer->Fill(fabs(p));
1299 }
1300
1301
1302 // Fill the kaons
1303 for( unsigned int index = 0; index < kaon->GetEntries(); ++index ){
1304 kaon->GetEvent(index);
1305
1306 int nhit;
1307 if( type == 0 )
1308 nhit = std::floor(b3nhit);
1309 else if( type == 1 )
1310 nhit = b2nhit;
1311
1312 if( dedx <= 0 || nhit <= m_lowernhit || nhit >= m_uppernhit )
1313 continue;
1314 if( fabs(p) < 0.5 and dedx < 1 ) continue;
1315
1316 double dedx_new = dedx;
1317 if( !mcflag ) dedx_new = m_hadsat.D2I(costh,m_hadsat.I2D(costh,1.0)*dedx);
1318 double dedx_res = gpar.sigmaPrediction(dedx,nhit,sqrt(1-costh*costh));
1319
1320 khdedx->Fill(fabs(p),dedx_new);
1321
1322 double pibg = fabs(p)/Widget::mpion;
1323 double dedx_cur_pi = gpar.dedxPrediction(pibg);
1324 double chi_new_pi = (dedx_new-dedx_cur_pi)/dedx_res;
1325
1326 kpredpi->Fill(fabs(p),dedx_cur_pi);
1327
1328 double kbg = fabs(p)/Widget::mkaon;
1329 double dedx_cur_k = gpar.dedxPrediction(kbg);
1330 double chi_new_k = (dedx_new-dedx_cur_k)/dedx_res;
1331
1332 kpredk->Fill(fabs(p),dedx_cur_k);
1333
1334 kchivp->Fill(fabs(p),(chi_new_k*chi_new_k-chi_new_pi*chi_new_pi));
1335
1336 kdenom->Fill(fabs(p));
1337 if( chi_new_k*chi_new_k < chi_new_pi*chi_new_pi ) knumer->Fill(fabs(p));
1338 }
1339
1340 TEfficiency* piteff = new TEfficiency(*pinumer,*pidenom);
1341 TEfficiency* kteff = new TEfficiency(*knumer,*kdenom);
1342
1343 outfile->cd();
1344 pinumer->Write();
1345 pidenom->Write();
1346 piteff->Write();
1347 pihdedx->Write();
1348 pipredpi->Write();
1349 pipredk->Write();
1350 pichivp->Write();
1351
1352 knumer->Write();
1353 kdenom->Write();
1354 kteff->Write();
1355 khdedx->Write();
1356 kpredpi->Write();
1357 kpredk->Write();
1358 kchivp->Write();
1359 outfile->Close();
1360}
nbins
const int nbins
Definition: CheckRung/check/dEdxVsCos.cxx:18
mass
double mass
Definition: CosmicGenerator.cxx:138
sigma
TTree * sigma
Definition: CurConst_dEdXCal.cxx:6
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
gr
TGraph * gr
Definition: DataBase/tau_mode.c:69
sin
double sin(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:210
HadronCalibration::fitSigmaVsSin
void fitSigmaVsSin(TString filename, std::string paramfile, int mcFlag, int type)
Definition: HadronWidget/HadronCalibration.cc:886
HadronCalibration::plotEfficiency
void plotEfficiency(TString filenames[5], TString saveas, std::string paramfile, int mcFlag, int type)
Definition: HadronWidget/HadronCalibration.cc:1157
HadronCalibration::plotBGCurve
void plotBGCurve(std::vector< TString > particles, TString filename, std::string paramfile)
Definition: HadronWidget/HadronCalibration.cc:381
HadronCalibration::fitBGCurve
void fitBGCurve(std::vector< TString > particles, TString filename, std::string paramfile)
Definition: HadronWidget/HadronCalibration.cc:6
HadronCalibration::fitSigmaVsNHit
void fitSigmaVsNHit(TString filename, std::string paramfile, int mcFlag, int type)
Definition: HadronWidget/HadronCalibration.cc:609
HadronSaturation::I2D
double I2D(double cosTheta, double I, double alpha, double gamma, double delta, double power, double ratio) const
Definition: HadronSaturation.cc:133
HadronSaturation::D2I
double D2I(double cosTheta, double D, double alpha, double gamma, double delta, double power, double ratio) const
Definition: HadronSaturation.cc:119
HadronSaturation::setParameters
void setParameters(double par[])
Definition: HadronSaturation.cc:37
WidgetParameterization::setCurvePars
void setCurvePars(int i, double val)
Definition: WidgetParameterization.h:62
WidgetParameterization::setNHitPars
void setNHitPars(int i, double val)
Definition: WidgetParameterization.h:74
WidgetParameterization::sinPrediction
double sinPrediction(double sin)
Definition: WidgetParameterization.cc:194
WidgetParameterization::printParameters
void printParameters(std::string infile)
Definition: WidgetParameterization.cc:44
WidgetParameterization::sigmaPrediction
double sigmaPrediction(double dedx, double nhit, double sin)
Definition: WidgetParameterization.cc:130
WidgetParameterization::nhitPrediction
double nhitPrediction(double nhit)
Definition: WidgetParameterization.cc:216
WidgetParameterization::setDedxPars
void setDedxPars(int i, double val)
Definition: WidgetParameterization.h:68
WidgetParameterization::setSinPars
void setSinPars(int i, double val)
Definition: WidgetParameterization.h:80
y
double y[1000]
Definition: draw_charge_space_coarse_II.cxx:8
sprintf
sprintf(cut,"kal_costheta0_em>-0.93&&kal_costheta0_em<0.93&&kal_pxy0_em>=0.05+%d*0.1&&kal_pxy0_em<0.15+%d*0.1&&NGch>=2", j, j)
nhits
int nhits
Definition: prepare_data.cxx:29
bg
float bg
Definition: prepare_data.cxx:28