BOSS 7.0.7
BESIII Offline Software System
Loading...
Searching...
No Matches
ElectronWidget/ElectronCalibration.cc
Go to the documentation of this file.
1#include "ElectronCalibration.h"
2
3ElectronCalibration::ElectronCalibration() :
4 m_filename(""),
5 m_constfilename(""),
6 m_mcFlag(1),
7 m_type(0),
8 m_fits(0),
9 m_docabins(20),
10 m_lowerdoca(-1.0),
11 m_upperdoca(1.0),
12 m_entabins(20),
13 m_lowerenta(-1.0),
14 m_upperenta(1.0),
15 m_costhbins(20)
16{}
17
18ElectronCalibration::ElectronCalibration( TString infile, TString constfile, int mcFlag, int type, int fits, int docabins, double upperdoca, double lowerdoca, int entabins, double upperenta, double lowerenta, int costhbins ) :
19 m_filename(infile),
20 m_constfilename(constfile),
21 m_mcFlag(mcFlag),
22 m_type(type),
23 m_fits(fits),
24 m_docabins(docabins),
25 m_upperdoca(upperdoca),
26 m_lowerdoca(lowerdoca),
27 m_entabins(entabins),
28 m_upperenta(upperenta),
29 m_lowerenta(lowerenta),
30 m_costhbins(costhbins)
31{}
32
33void
34ElectronCalibration::fitRunGains( TFile* outfile ){
35
36 // supress TCanvas::Print messages
37 // gErrorIgnoreLevel = kWarning;
38
39 TFile* infile = new TFile(m_filename);
40 TTree* electron = (TTree*)infile->Get("track");
41
42 std::cout << "ElectronCalibration: getting run gains for file " << m_filename << std::endl;
43
44 // --------------------------------------------------
45 // INITIALIZE CONTAINERS
46 // --------------------------------------------------
47
48 int run; // run number per event
49 int nhits; // number of (all) hits on this track
50 double dedx; // dE/dx with electron saturation correction
51 double costh; // cosine of track polar angle
52
53 // Belle II variables
54 int b2nhit; // number of hits on this track
55
56 // BES III variables
57 double b3nhit; // number of hits on this track
58
59 electron->SetBranchAddress("run", &run);
60 electron->SetBranchAddress("dedxsat", &dedx);
61 electron->SetBranchAddress("costh", &costh);
62
63 if( m_type == 0 )
64 electron->SetBranchAddress("numGoodHits", &b3nhit);
65 else if( m_type == 1 )
66 electron->SetBranchAddress("numGoodLayerHits", &b2nhit);
67
68 // --------------------------------------------------
69 // LOOP OVER EVENTS AND FILL CONTAINERS
70 // --------------------------------------------------
71
72 // Fill histograms and fit with Gaussian functions
73 // to extract the means and errors
74 TTree* satTree = new TTree("electron","dE/dx means and errors");
75
76 int fitrun; // run number for this bin
77 double fitmean; // mean dE/dx value for this bin
78 double fitmeanerr; // error on mean dE/dx value for this bin
79 double fitwidth; // dE/dx width for this bin
80 double fitwidtherr; // error on dE/dx width for this bin
81
82 satTree->Branch("run",&fitrun,"run/I");
83 satTree->Branch("dedx",&fitmean,"dedx/D");
84 satTree->Branch("dedxerr",&fitmeanerr,"dedxerr/D");
85 satTree->Branch("width",&fitwidth,"width/D");
86 satTree->Branch("widtherr",&fitwidtherr,"widtherr/D");
87
88 // Create the histograms to be fit (before correction)
89 TH1F* hdedx_run = new TH1F("dedx_run","dedx_run",200,0,2);
90
91 // Container to store the constants
92 std::map<int, double> m_rungains;
93
94 // Fill the histograms to be fitted
95 double lastrun = -1;
96 int nentries = electron->GetEntries();
97 for( unsigned int index = 0; index < nentries; ++index ){
98 electron->GetEvent(index);
99
100 int nhit = 0;
101 if( m_type == 0 )
102 nhit = std::floor(b3nhit);
103 else if( m_type == 1 )
104 nhit = b2nhit;
105
106 // clean up bad events and restrict the momentum range
107 if( nhit < 0 || nhit > 100 || dedx <= 0 || costh != costh )
108 continue;
109
110 if( lastrun == -1 ) lastrun = run;
111
112 if( run == lastrun && index < (nentries-1) ) hdedx_run->Fill(dedx);
113 else{
114 if( hdedx_run->Integral() < 50 ){
115 lastrun = run;
116 continue;
117 }
118
119 // fill some details for this bin
120 fitrun = lastrun;
121 // fit the dE/dx distribution in bins of run number
122 hdedx_run->Fit("gaus","ql");
123 fitmean = hdedx_run->GetFunction("gaus")->GetParameter(1);
124 fitmeanerr = hdedx_run->GetFunction("gaus")->GetParError(1);
125 fitwidth = hdedx_run->GetFunction("gaus")->GetParameter(2);
126 fitwidtherr = hdedx_run->GetFunction("gaus")->GetParError(2);
127 satTree->Fill();
128
129 m_rungains[fitrun] = fitmean;
130
131 hdedx_run->Reset();
132 hdedx_run->Fill(dedx);
133
134 lastrun = run;
135 }
136 }
137
138 satTree->Print();
139 /*
140 outfile->cd();
141 TH2F* rungains = new TH2F("rungains","",m_rungains.size(),0,m_rungains.size());
142 rungains->SetStats(0);
143 satTree->Project("runsigma","width:run");
144 TH1F* runsigma = (TH1F*)outfile->Get("runsigma");
145 runsigma->SetStats(0);
146
147 TCanvas* can = new TCanvas("electrons","",1200,600);
148 can->Divide(2,1);
149 can->cd(1);
150 rungains->SetTitle(";Run number;dE/dx");
151 rungains->Draw("P");
152
153 can->cd(2);
154 runsigma->SetTitle(";Run number;dE/dx #sigma");
155 runsigma->Draw("P");
156 can->SaveAs("plots/rungains.eps");
157
158 delete hdedx_run;
159 delete can;
160 */
161
162 // write out the data to file
163 outfile->cd();
164 // rungains->Write();
165 // runsigma->Write();
166 satTree->Write();
167 outfile->Close();
168 infile->Close();
169}
170
171void
172ElectronCalibration::plotRunGains( TString filename ){
173
174 std::cout << "ElectronCalibration: plotting run gains for file " << m_filename << std::endl;
175
176 TFile* infile = new TFile(filename);
177 TTree* intree = (TTree*)infile->Get("electron");
178
179 intree->Project("rungains","dedx:run");
180 TH1F* rungains = (TH1F*)infile->Get("rungains");
181 intree->Project("runsigma","dedx:width");
182 TH1F* runsigma = (TH1F*)infile->Get("runsigma");
183
184 TCanvas* can = new TCanvas("electrons","",1200,600);
185 can->Divide(2,1);
186 can->cd(1);
187 rungains->SetTitle(";Run number;dE/dx");
188 rungains->Draw("P");
189
190 can->cd(2);
191 runsigma->SetTitle(";Run number;dE/dx #sigma");
192 runsigma->Draw("P");
193 can->SaveAs("plots/rungains.eps");
194
195 delete can;
196
197 infile->Close();
198}
199
200void
201ElectronCalibration::TwoDCorrection( TFile* outfile ){
202
203 std::cout << "ElectronCalibration: reading in file " << m_filename << std::endl;
204
205 TFile* infile = new TFile(m_filename);
206 TTree* track = (TTree*)infile->Get("track");
207
208 // --------------------------------------------------
209 // INITIALIZE CONTAINERS
210 // --------------------------------------------------
211
212 int kMaxEntries = 100;
213 double doca[kMaxEntries]; // distance of closest approach for each hit
214 double enta[kMaxEntries]; // CDC cell entrance angle for each hit
215 double dedxhit[kMaxEntries]; // dE/dx for each hit
216 int nhits; // number of hits on this track
217
218 // Belle II variables
219 int b2nhit; // number of hits on this track
220
221 // BES III variables
222 double b3nhit; // number of hits on this track
223
224 track->SetBranchAddress("doca", doca);
225 track->SetBranchAddress("enta", enta);
226 track->SetBranchAddress("dedxhit", dedxhit);
227
228 if( m_type == 0 )
229 track->SetBranchAddress("numGoodHits", &b3nhit);
230 else if( m_type == 1 )
231 track->SetBranchAddress("numGoodLayerHits", &b2nhit);
232
233 double docastep = (m_upperdoca-m_lowerdoca)/m_docabins;
234 double entastep[m_entabins];
235
236 // --------------------------------------------------
237 // SET THE ENTA BIN SIZES
238 // --------------------------------------------------
239
240 // Fill the histograms to be fitted
241 TH1F* totalenta = new TH1F("totalenta","",314,0,0.785);
242 for( unsigned int index = 0; index < track->GetEntries(); ++index ){
243 track->GetEvent(index);
244
245 if( m_type == 0 )
246 nhits = std::floor(b3nhit);
247 else if( m_type == 1 )
248 nhits = b2nhit;
249 for( int hit = 0; hit < nhits; ++hit ){
250
251 // use approx. rotational symmetry to map to [-pi/4,pi/4]
252 if( enta[hit] > PI/4.0 ) enta[hit] -= PI/2.0;
253 else if( enta[hit] < -1.0*PI/4.0 ) enta[hit] += PI/2.0;
254 totalenta->Fill(std::abs(enta[hit]));
255 }
256 }
257 entastep[0] = -0.785;
258 entastep[39] = 0.785;
259 int bin = 20; int binmin = 1;
260 double total = totalenta->Integral()*2.0/(m_entabins+1);
261 for( int i = 2; i <= 314; ++i ){
262 if( totalenta->Integral(binmin,i) >= total ){
263 entastep[bin] = 0.0025*i;
264 bin++;
265 binmin = i+1;
266 }
267 }
268 for( int i = 0; i < m_entabins; ++i ){
269 if( i > 0 && i < 20 ) entastep[i] = -1.0*entastep[39-i];
270 }
271
272 // --------------------------------------------------
273 // OUTPUT TTREE
274 // --------------------------------------------------
275
276 TTree* outTree = new TTree("electrons","2D means and errors");
277
278 double satdoca; // DOCA for this bin
279 double satenta; // entrance angle for this bin
280 double satdedx; // mean dE/dx value for this bin
281 double saterror; // error on ^
282
283 outTree->Branch("doca",&satdoca,"doca/D");
284 outTree->Branch("enta",&satenta,"enta/D");
285 outTree->Branch("dedx",&satdedx,"dedx/D");
286 outTree->Branch("error",&saterror,"error/D");
287
288 TH2F* twod = new TH2F("twod","",m_docabins,m_lowerdoca,m_upperdoca,m_entabins,m_lowerenta,m_upperenta);
289
290 // --------------------------------------------------
291 // IF DOING TRUNCATION, FILL CONTAINERS AND TRUNCATE
292 // --------------------------------------------------
293
294 if( m_fits == 0 ){
295
296 std::vector<double> docaent[m_docabins][m_entabins];
297 int docaenthits[m_docabins][m_entabins];
298 for( int i = 0; i < m_docabins; ++i ){
299 for( int j = 0; j < m_docabins; ++j ){
300 docaent[i][j].clear();
301 docaenthits[i][j] = 0;
302 }
303 }
304
305 for( unsigned int index = 0; index < track->GetEntries(); ++index ){
306 track->GetEvent(index);
307
308 if( m_type == 0 )
309 nhits = std::floor(b3nhit);
310 else if( m_type == 1 )
311 nhits = b2nhit;
312 for( int hit = 0; hit < nhits; ++hit ){
313
314 // use approx. rotational symmetry to map to [-pi/4,pi/4]
315 if( enta[hit] > PI/4.0 ) enta[hit] -= PI/2.0;
316 else if( enta[hit] < -1.0*PI/4.0 ) enta[hit] += PI/2.0;
317
318 if( doca[hit] > m_upperdoca || doca[hit] < m_lowerdoca ) continue;
319 if( enta[hit] > m_upperenta || enta[hit] < m_lowerenta ) continue;
320
321 int docaBin = (int)((doca[hit]-m_lowerdoca)/(m_upperdoca-m_lowerdoca) * m_docabins);
322 int entaBin = (int)((enta[hit]-m_lowerenta)/(m_upperenta-m_lowerenta) * m_entabins);
323
324 docaent[docaBin][entaBin].push_back(dedxhit[hit]);
325 docaenthits[docaBin][entaBin]++;
326 } // end of hit loop
327 } // end of event loop
328
329 double tmean;
330 ElectronCorrection ecor;
331 for( int i = 0; i < m_docabins; ++i ){
332 for( int j = 0; j < m_entabins; ++j ){
333
334 // fill some details for this bin
335 satdoca = m_lowerdoca+0.5*docastep+i*docastep;
336 satenta = m_lowerenta+0.5*entastep[j]+j*entastep[j];
337
338 // calculate the truncated mean and error
339 ecor.calculateMeans(tmean,satdedx,saterror,docaent[i][j],docaenthits[i][j]);
340 twod->SetBinContent(i,j,satdedx);
341
342 std::cout << i << "\t" << j << "\t" << satdedx << "\t" << saterror << std::endl;
343
344 // fill the output TTree
345 outTree->Fill();
346 }
347 }
348 }
349
350 // --------------------------------------------------
351 // IF DOING FITS, LOOP OVER EVENTS AND FILL CONTAINERS
352 // --------------------------------------------------
353
354 if( m_fits == 1 ){
355
356 // Create the histograms to be fit (before correction)
357 TH1F* hdoca_enta[m_docabins][m_entabins];
358
359 // initialize the histograms
360 for( int i = 0; i < m_docabins; ++i ){
361 for( int j = 0; j < m_entabins; ++j ){
362 char histname[100];
363 sprintf(histname,"doca_%d_enta_%d",i,j);
364
365 hdoca_enta[i][j] = new TH1F(histname,histname,200,0,200);
366 }
367 }
368
369 // Fill the histograms to be fitted
370 for( unsigned int index = 0; index < track->GetEntries(); ++index ){
371 track->GetEvent(index);
372
373 if( m_type == 0 )
374 nhits = std::floor(b3nhit);
375 else if( m_type == 1 )
376 nhits = b2nhit;
377 for( int hit = 0; hit < nhits; ++hit ){
378
379 // use approx. rotational symmetry to map to [-pi/4,pi/4]
380 if( enta[hit] > PI/4.0 ) enta[hit] -= PI/2.0;
381 else if( enta[hit] < -1.0*PI/4.0 ) enta[hit] += PI/2.0;
382
383 if( doca[hit] > m_upperdoca || doca[hit] < m_lowerdoca ) continue;
384 if( enta[hit] > m_upperenta || enta[hit] < m_lowerenta ) continue;
385
386 int docaBin = (int)((doca[hit]-m_lowerdoca)/(m_upperdoca-m_lowerdoca) * m_docabins);
387 int entaBin = (int)((enta[hit]-m_lowerenta)/(m_upperenta-m_lowerenta) * m_entabins);
388
389 hdoca_enta[docaBin][entaBin]->Fill(dedxhit[hit]);
390 } // end of hit loop
391 } // end of event loop
392
393
394 // --------------------------------------------------
395 // FIT IN BINS OF DOCA AND ENTRANCE ANGLE
396 // --------------------------------------------------
397
398 // fit the histograms
399 int fs1; // make sure none of the fits fail...
400 for( int i = 0; i < m_docabins; ++i ){
401 for( int j = 0; j < m_entabins; ++j ){
402
403 // fill some details for this bin
404 satdoca = m_lowerdoca+0.5*docastep+i*docastep;
405 satenta = m_lowerenta+0.5*entastep[j]+j*entastep[j];
406
407 int fs; // check if the fit succeeded
408 fs = hdoca_enta[i][j]->Fit("landau","ql");
409 double mean = hdoca_enta[i][j]->GetFunction("landau")->GetParameter(1);
410 double width = hdoca_enta[i][j]->GetFunction("landau")->GetParameter(2);
411 // fs = hdoca_enta[i][j]->Fit("landau","ql","",mean-2.5*width,mean+2.5*width);
412 fs = hdoca_enta[i][j]->Fit("landau","ql");
413 if( fs != 0 ) fs1++;
414
415 satdedx = hdoca_enta[i][j]->GetFunction("landau")->GetParameter(1);
416 saterror = hdoca_enta[i][j]->GetFunction("landau")->GetParError(1);
417
418 if( fs == 0 ) twod->SetBinContent(i,j,satdedx);
419
420 // fill the tree for this bin
421 outTree->Fill();
422 }
423 }
424 if( fs1 != 0 ) std::cout << "\t\t" <<"MEAN FITS FAILED: " << fs1 << std::endl;
425
426
427 // Print the histograms for quality control
428 TCanvas* ctmp = new TCanvas("tmp","tmp",900,900);
429 ctmp->Divide(3,3);
430 std::stringstream psname; psname << "plots/twoDFits.ps[";
431 ctmp->Print(psname.str().c_str());
432 psname.str(""); psname << "plots/twoDFits.ps";
433 for( int i = 0 ; i < m_docabins; ++i ){
434 for( int j = 0; j < m_entabins; ++j ){
435 ctmp->cd(j%9+1);
436 hdoca_enta[i][j]->Draw();
437 if((j+1)%9==0)
438 ctmp->Print(psname.str().c_str());
439 }
440 }
441 psname.str(""); psname << "plots/twoDFits.ps]";
442 ctmp->Print(psname.str().c_str());
443 delete ctmp;
444
445 for( int i = 0; i < m_docabins; ++i ){
446 for( int j = 0; j < m_entabins; ++j ){
447 delete hdoca_enta[i][j];
448 }
449 }
450 delete totalenta;
451 }
452
453 // write out the data to file
454 outfile->cd();
455 outTree->Write();
456 twod->Write();
457 infile->Close();
458}
459
460void
461ElectronCalibration::SaturationCorrection( TFile* outfile ){
462
463 TFile* infile = new TFile(m_filename);
464 TTree* electron = (TTree*)infile->Get("track");
465
466 std::cout << "ElectronCalibration: electron saturation correction " << m_filename << std::endl;
467
468 // --------------------------------------------------
469 // INITIALIZE CONTAINERS
470 // --------------------------------------------------
471
472 int run; // run number per event
473 double dedxpub; // dE/dx without electron saturation correction
474 double dedx; // dE/dx with electron saturation correction
475 double costh; // cosine of track polar angle
476
477 // Belle II variables
478 int b2nhit; // number of hits on this track
479
480 // BES III variables
481 double b3nhit; // number of hits on this track
482
483 electron->SetBranchAddress("run", &run);
484 electron->SetBranchAddress("dedx", &dedxpub);
485 electron->SetBranchAddress("dedxsat", &dedx);
486 electron->SetBranchAddress("costh", &costh);
487
488 if( m_type == 0 )
489 electron->SetBranchAddress("numGoodHits", &b3nhit);
490 else if( m_type == 1 )
491 electron->SetBranchAddress("numGoodLayerHits", &b2nhit);
492
493 // --------------------------------------------------
494 // LOOP OVER EVENTS AND FILL CONTAINERS
495 // --------------------------------------------------
496
497 // Fill histograms and fit with Gaussian functions
498 // to extract the means and errors
499 TTree* satTree = new TTree("electron","saturation correction");
500
501 double costhbin; // center of cos(theta) bin
502 double fitmean; // mean dE/dx value for this bin
503 double fitmeanerr; // error on mean dE/dx value for this bin
504 double fitwidth; // dE/dx width for this bin
505 double fitwidtherr; // error on dE/dx width for this bin
506
507 satTree->Branch("costh",&costhbin,"costh/D");
508 satTree->Branch("dedx",&fitmean,"dedx/D");
509 satTree->Branch("dedxerr",&fitmeanerr,"dedxerr/D");
510 satTree->Branch("width",&fitwidth,"width/D");
511 satTree->Branch("widtherr",&fitwidtherr,"widtherr/D");
512
513 // Create the histograms to be fit (before correction)
514 TH1F* hdedx_costh[m_costhbins];
515 for( int i = 0; i < m_costhbins; ++i ){
516 hdedx_costh[i] = new TH1F(TString::Format("dedx_costh%i",i),";dE/dx;cos(#theta)",200,0,2);
517 }
518
519 // Fill the histograms to be fitted
520 for( unsigned int index = 0; index < electron->GetEntries(); ++index ){
521 electron->GetEvent(index);
522
523 int nhit = 0;
524 if( m_type == 0 )
525 nhit = std::floor(b3nhit);
526 else if( m_type == 1 )
527 nhit = b2nhit;
528
529 // clean up bad events and restrict the momentum range
530 if( nhit < 0 || nhit > 100 || dedx <= 0 || costh != costh )
531 continue;
532
533 int bin = floor((costh+1.0)/2*m_costhbins);
534 hdedx_costh[bin]->Fill(dedxpub);
535 }
536
537 // Fit the histograms
538 for( unsigned int i = 0; i < m_costhbins; ++i ){
539 costhbin = (2.0*i+1)/m_costhbins-1;
540 hdedx_costh[i]->Fit("gaus","ql");
541 fitmean = hdedx_costh[i]->GetFunction("gaus")->GetParameter(1);
542 fitmeanerr = hdedx_costh[i]->GetFunction("gaus")->GetParError(1);
543 fitwidth = hdedx_costh[i]->GetFunction("gaus")->GetParameter(2);
544 fitwidtherr = hdedx_costh[i]->GetFunction("gaus")->GetParError(2);
545 satTree->Fill();
546 }
547
548 outfile->cd();
549 satTree->Project("satmean","dedx:costh");
550 TH1F* satmean = (TH1F*)outfile->Get("satmean");
551 satmean->SetStats(0);
552 satTree->Project("satsigma","width:costh");
553 TH1F* satsigma = (TH1F*)outfile->Get("satsigma");
554 satsigma->SetStats(0);
555
556 TCanvas* can = new TCanvas("electrons","",1200,600);
557 can->Divide(2,1);
558 can->cd(1);
559 satmean->SetTitle(";cos(#theta);dE/dx");
560 satmean->Draw("P");
561
562 can->cd(2);
563 satsigma->SetTitle(";cos(#theta);dE/dx sigma");
564 satsigma->Draw("P");
565 can->SaveAs("plots/satmean.eps");
566
567 for( int i = 0; i < m_costhbins; ++i ){
568 delete hdedx_costh[i];
569 }
570 delete can;
571
572 // write out the data to file
573 outfile->cd();
574 satmean->Write();
575 satsigma->Write();
576 satTree->Write();
577 outfile->Close();
578
579 infile->Close();
580}
electron
@ electron
Definition: DstMdcDedx.h:9
bin
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per bin
Definition: FoamA.h:85
PI
const double PI
Definition: PipiJpsi.cxx:55
ElectronCorrection::calculateMeans
void calculateMeans(double &mean, double &truncatedMean, double &truncatedMeanErr, double dedx[], int size) const
Definition: ElectronWidget/ElectronCorrection.cc:168
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