BOSS 7.0.9
BESIII Offline Software System
Loading...
Searching...
No Matches
HadronPrep.cc
Go to the documentation of this file.
1#include "HadronPrep.h"
2
3HadronPrep::HadronPrep() :
4 m_filename("widget.gen.root"),
5 m_mcFlag(1),
6 m_type(0),
7 m_bgbins(15),
8 m_upperbg(5.0),
9 m_lowerbg(0.1),
10 m_cosbins(18),
11 m_uppercos(0.93),
12 m_lowercos(0.0)
13{}
14
15HadronPrep::HadronPrep( TString infile, int mcFlag, int type, int bgbins, double upperbg, double lowerbg, int cosbins, double uppercos, double lowercos ){
16 m_filename = infile;
17 m_mcFlag = mcFlag;
18 m_type = type;
19 m_bgbins = bgbins;
20 m_upperbg = upperbg;
21 m_lowerbg = lowerbg;
22 m_cosbins = cosbins;
23 m_uppercos = uppercos;
24 m_lowercos = lowercos;
25}
26
27void
28HadronPrep::FormatGraph( TGraphErrors* gr, int flag ){
29 if( flag == 0 ){
30 gr->SetTitle(";cos(#theta);#chi_{mean}");
31 gr->SetMarkerStyle(24);
32 gr->SetMarkerColor(2);
33 gr->SetMarkerSize(.7);
34 gr->SetMaximum(1);
35 gr->SetMinimum(-1);
36 gr->GetXaxis()->SetRangeUser(-1,1);
37 }
38 else if( flag == 1 ){
39 gr->SetMarkerStyle(22);
40 gr->SetMarkerColor(9);
41 gr->SetMarkerSize(.7);
42 }
43 else if( flag == 2 ){
44 gr->SetTitle(";cos(#theta);#sigma");
45 gr->SetMarkerStyle(24);
46 gr->SetMarkerColor(2);
47 gr->SetMarkerSize(.7);
48 gr->SetMaximum(2);
49 gr->SetMinimum(0);
50 gr->GetXaxis()->SetRangeUser(-1,1);
51 }
52}
53
54void
55HadronPrep::bgCosThetaFits( TString particle, TFile* outfile, bool correct, std::string paramfile ){
56
57 // supress TCanvas::Print messages
58 gErrorIgnoreLevel = kWarning;
59
60 m_gpar.setParameters(paramfile);
61
62 double mass = 0.0;
63 if( particle == "pion" ) mass = Widget::mpion;
64 else if( particle == "kaon" ) mass = Widget::mkaon;
65 else if( particle == "proton" ) mass = Widget::mproton;
66 else if( particle == "muon" ) mass = Widget::mmuon;
67 else if( particle == "electron" ) mass = Widget::melectron;
68 if( mass == 0.0 ) exit(1);
69
70 TFile* infile = new TFile(m_filename);
71 TTree* hadron = (TTree*)infile->Get("track");
72
73 std::cout << "HadronPrep: reading in file " << m_filename << std::endl;
74 std::cout << "\t beta-gamma range: " << m_lowerbg << " to " << m_upperbg << std::endl;
75 std::cout << "\t cos(theta) range: " << m_lowercos << " to " << m_uppercos << std::endl;
76
77 // --------------------------------------------------
78 // INITIALIZE CONTAINERS
79 // --------------------------------------------------
80
81 double dedxpub; // dE/dx without electron saturation correction
82 double dedx; // dE/dx with electron saturation correction
83 double p; // track momentum
84 double bg; // track beta-gamma
85 double costh; // cosine of track polar angle
86 double db; // the nearest distance to the IP in the xy plane
87 double dz; // the nearest distance to the IP in the z plane
88 double chiPi; // PID chi value
89 double eop; // energy over momentum in the calorimeter
90
91 // Belle II variables
92 int b2nhit; // number of hits on this track
93
94 // BES III variables
95 double b3nhit; // number of hits on this track
96
97 hadron->SetBranchAddress("dedx", &dedxpub);
98 hadron->SetBranchAddress("dedxsat", &dedx);
99 hadron->SetBranchAddress("pF", &p);
100 hadron->SetBranchAddress("costh", &costh);
101 hadron->SetBranchAddress("db", &db);
102 hadron->SetBranchAddress("dz", &dz);
103 hadron->SetBranchAddress("chiPi", &chiPi);
104 hadron->SetBranchAddress("eopst", &eop);
105
106 if( m_type == 0 )
107 hadron->SetBranchAddress("numGoodHits", &b3nhit);
108 else if( m_type == 1 )
109 hadron->SetBranchAddress("lNHitsUsed", &b2nhit);
110 // hadron->SetBranchAddress("numGoodLayerHits", &b2nhit);
111
112 double bgstep = (m_upperbg-m_lowerbg)/m_bgbins;
113 double cosstep = (m_uppercos-m_lowercos)/m_cosbins;
114
115 // Create the histograms to be fit (before correction)
116 TH1F* hdedx_costh[m_bgbins][m_cosbins];
117 TH1F* hdedx_costh_pub[m_bgbins][m_cosbins];
118 TH1F* hchi_costh[m_bgbins][m_cosbins];
119 TH1F* hchi_costh_pub[m_bgbins][m_cosbins];
120
121 // initialize the histograms
122 for( int i = 0; i < m_bgbins; ++i ){
123 for( int j = 0; j < m_cosbins; ++j ){
124 char histname[100],histname2[100],histname3[100],histname4[100],histname5[100];
125 sprintf(histname, "dedx_costh_p_%d_cos_%d",i,j);
126 sprintf(histname2, "chi_costh_p_%d_cos_%d",i,j);
127 sprintf(histname3, "chipub_costh_p_%d_cos_%d",i,j);
128 sprintf(histname4, "dedxpub_costh_p_%d_cos_%d",i,j);
129 sprintf(histname5, "dedxnew_costh_p_%d_cos_%d",i,j);
130
131
132 if( particle == "electron" ){
133 hdedx_costh[i][j] = new TH1F(histname,histname,200,0,2);
134 hchi_costh[i][j] = new TH1F(histname2,histname2,200,-2,2);
135 hchi_costh_pub[i][j] = new TH1F(histname3,histname3,200,-5,5);
136 hdedx_costh_pub[i][j] = new TH1F(histname4,histname4,200,0,2);
137 }
138 else if( particle == "proton" ){
139 if( i < 5 ){
140 hdedx_costh[i][j] = new TH1F(histname,histname,200,0,6);
141 hchi_costh[i][j] = new TH1F(histname2,histname2,200,-30,30);
142 hchi_costh_pub[i][j] = new TH1F(histname3,histname3,200,0,40);
143 hdedx_costh_pub[i][j] = new TH1F(histname4,histname4,200,0,6);
144 }
145 else if( i < 9 ){
146 hdedx_costh[i][j] = new TH1F(histname,histname,200,0,4);
147 hchi_costh[i][j] = new TH1F(histname2,histname2,200,-5,5);
148 hchi_costh_pub[i][j] = new TH1F(histname3,histname3,200,-20,20);
149 hdedx_costh_pub[i][j] = new TH1F(histname4,histname4,200,0,4);
150 }
151 else{
152 hdedx_costh[i][j] = new TH1F(histname,histname,200,0,2);
153 hchi_costh[i][j] = new TH1F(histname2,histname2,200,-5,5);
154 hchi_costh_pub[i][j] = new TH1F(histname3,histname3,200,-5,5);
155 hdedx_costh_pub[i][j] = new TH1F(histname4,histname4,200,0,2);
156 }
157 }
158 else if( i < 2 ){
159 hdedx_costh[i][j] = new TH1F(histname,histname,200,0,4);
160 hchi_costh[i][j] = new TH1F(histname2,histname2,200,-5,5);
161 hchi_costh_pub[i][j] = new TH1F(histname3,histname3,200,-5,5);
162 hdedx_costh_pub[i][j] = new TH1F(histname4,histname4,200,0,4);
163 }
164 else{
165 hdedx_costh[i][j] = new TH1F(histname,histname,200,0,2);
166 hchi_costh[i][j] = new TH1F(histname2,histname2,200,-5,5);
167 hchi_costh_pub[i][j] = new TH1F(histname3,histname3,200,-5,5);
168 hdedx_costh_pub[i][j] = new TH1F(histname4,histname4,200,0,2);
169 }
170 }
171 }
172
173 // Create some containers to calculate averages
174 double sumcos[m_bgbins][m_cosbins];
175 double sumsin[m_bgbins][m_cosbins];
176 double sumbg[m_bgbins][m_cosbins];
177 double sumressq[m_bgbins][m_cosbins];
178 int sumsize[m_bgbins][m_cosbins];
179 for( int i = 0; i < m_bgbins; ++i ){
180 for( int j = 0; j < m_cosbins; ++j ){
181 sumcos[i][j] = 0;
182 sumsin[i][j] = 0;
183 sumbg[i][j] = 0;
184 sumressq[i][j] = 0;
185 sumsize[i][j] = 0;
186 }
187 }
188
189 // Create some containers for means and errors
190 double means[m_bgbins][m_cosbins];
191 double errors[m_bgbins][m_cosbins];
192 double widths[m_bgbins][m_cosbins];
193 for( int i = 0; i < m_bgbins; ++i ){
194 for( int j = 0; j < m_cosbins; ++j ){
195 means[i][j] = 0;
196 errors[i][j] = 0;
197 widths[i][j] = 0;
198 }
199 }
200
201
202 // get the hadron saturation parameters
203 // if the parameters do not exist, use the values in the default constructor
204 m_hadsat.setParameters("sat-pars.fit.txt");
205
206 // --------------------------------------------------
207 // LOOP OVER EVENTS AND FILL CONTAINERS
208 // --------------------------------------------------
209
210 // Fill the histograms to be fitted
211 for( unsigned int index = 0; index < hadron->GetEntries(); ++index ){
212 hadron->GetEvent(index);
213
214 bg = fabs(p)/mass;
215 costh = fabs(costh);
216
217 int nhit;
218 if( m_type == 0 )
219 nhit = std::floor(b3nhit);
220 else if( m_type == 1 )
221 nhit = b2nhit;
222
223 // clean up bad events and restrict the momentum range
224 if( nhit < 0 || nhit > 100 || dedx <= 0 || costh != costh ||
225 bg <= m_lowerbg || bg >= m_upperbg || costh <= m_lowercos || costh >= m_uppercos )
226 continue;
227
228 // apply an E/p cut for pions
229 if( particle == "pion" && eop > 0.75 ) continue;
230
231 // use loose dE/dx restrictions to remove contaminations
232 if( particle == "proton" && (dedxpub < 0 || dedxpub > 100 || dedxpub < 0.85/fabs(p)) )
233 continue;
234 if( particle == "kaon" && dedxpub < 0.4/fabs(p) )
235 continue;
236
237 int bgBin = (int)((bg-m_lowerbg)/(m_upperbg-m_lowerbg) * m_bgbins);
238 int cosBin = (int)((costh-m_lowercos)/(m_uppercos-m_lowercos) * m_cosbins);
239
240 double dedx_pre = dedx;
241 double dedx_new = dedx_pre;
242 if( !m_mcFlag ) dedx_new = m_hadsat.D2I(costh,m_hadsat.I2D(costh,1.0)*dedx);
243 double dedx_cur = m_gpar.dedxPrediction(bg);
244 double dedx_res = m_gpar.sigmaPrediction(dedx_cur,nhit,sqrt(1-costh*costh));
245 if( dedx_res == 0 ){
246 std::cout << "RESOLUTION IS ZERO!!!" << std::endl;
247 continue;
248 }
249
250 double chi_new = (dedx_new-dedx_cur)/dedx_res;
251 if( particle == "electron" ) chi_new = (dedx_new-1)/dedx_res;
252 if( correct ) hdedx_costh[bgBin][cosBin]->Fill(dedx_new);
253 else hdedx_costh[bgBin][cosBin]->Fill(dedx_pre);
254
255 if( correct ) hchi_costh[bgBin][cosBin]->Fill(chi_new);
256 else hchi_costh[bgBin][cosBin]->Fill(chiPi);
257
258 hdedx_costh_pub[bgBin][cosBin]->Fill(dedxpub);
259 hchi_costh_pub[bgBin][cosBin]->Fill(chiPi);
260
261 // if( fabs(bg-(m_lowerbg+0.5*bgstep+bgBin*bgstep)) < 0.5*bgstep &&
262 // fabs(costh-(m_lowercos+0.5*cosstep+cosBin*cosstep)) < 0.5*cosstep ){
263 sumcos[bgBin][cosBin] += costh;
264 sumsin[bgBin][cosBin] += sqrt(1-costh*costh);
265 sumbg[bgBin][cosBin] += bg;
266 sumressq[bgBin][cosBin] += pow(dedx_res,2);
267 sumsize[bgBin][cosBin] += 1;
268 // }
269
270 }// end of event loop
271
272 // --------------------------------------------------
273 // FIT IN BINS OF BETA-GAMMA AND COS(THETA)
274 // --------------------------------------------------
275
276 // fit the histograms with Gaussian functions
277 // and extract the means and errors
278 TTree* satTree = new TTree(particle,"dE/dx means and errors");
279
280 double satbg; // beta-gamma value for this bin
281 double satcosth; // cos(theta) value for this bin
282 double satdedx; // mean dE/dx value for this bin
283 double saterror; // error on ^
284 double satbg_avg; // average beta-gamma value for this sample
285 double satcosth_avg; // average cos(theta) value for this sample
286 double satsinth_avg; // average sin(theta) value for this sample
287 double satdedxres_avg; // average dE/dx error squared for this sample
288 double satpubdedx; // mean "public" dE/dx value for this bin
289 double satpuberror; // error on ^
290 double satpubwidth; // width of ^ distribution
291 double satchi; // mean chi value for this bin
292 double satchierr; // error on ^
293 double satchiwidth; // width of ^ distribution
294 double satchipub; // mean "public" chi value for this bin
295 double satchipuberr; // error on ^
296 double satchipubwidth; // width of ^ distribution
297 double ratio; // ratio of the predicted mean to that of the average
298
299 satTree->Branch("bg",&satbg,"bg/D");
300 satTree->Branch("costh",&satcosth,"costh/D");
301 satTree->Branch("dedx",&satdedx,"dedx/D");
302 satTree->Branch("error",&saterror,"error/D");
303 satTree->Branch("bg_avg",&satbg_avg,"bg_avg/D");
304 satTree->Branch("costh_avg",&satcosth_avg,"costh_avg/D");
305 satTree->Branch("sinth_avg",&satsinth_avg,"sinth_avg/D");
306 satTree->Branch("dedxres_avg",&satdedxres_avg,"dedxres_avg/D");
307 satTree->Branch("dedx_pub",&satpubdedx,"dedx_pub/D");
308 satTree->Branch("dedxerr_pub",&satpuberror,"dedxerr_pub/D");
309 satTree->Branch("chi",&satchi,"chi/D");
310 satTree->Branch("chi_pub",&satchipub,"chi_pub/D");
311 satTree->Branch("sigma",&satchiwidth,"sigma/D");
312 satTree->Branch("sigma_pub",&satchipubwidth,"sigma_pub/D");
313 satTree->Branch("ratio",&ratio,"ratio/D");
314
315 // Fit the histograms
316 int fs1=0, fs2=0, fs3=0, fs4=0;
317 double dedxmax = 1.0, dedxmin = 1.0;
318 for( int i = 0; i < m_bgbins; ++i ){
319 for( int j = 0; j < m_cosbins; ++j ){
320
321 std::cout << "\t\t" << "i= " << i << " j= " << j << std::endl;
322 // fill some details for this bin
323 satbg = m_lowerbg+0.5*bgstep+i*bgstep;
324 satcosth = m_lowercos+0.5*cosstep+j*cosstep;
325
326 satbg_avg = sumbg[i][j]/sumsize[i][j];
327 satcosth_avg = sumcos[i][j]/sumsize[i][j];
328 satsinth_avg = sumsin[i][j]/sumsize[i][j];
329 satdedxres_avg = sumressq[i][j]/sumsize[i][j];
330
331 int fs;
332 // fit the dE/dx distribution in bins of beta-gamma and cosine
333 fs = hdedx_costh[i][j]->Fit("gaus","ql");
334 double mean = hdedx_costh[i][j]->GetFunction("gaus")->GetParameter(1);
335 double width = hdedx_costh[i][j]->GetFunction("gaus")->GetParameter(2);
336 // fs = hdedx_costh[i][j]->Fit("gaus","ql","",mean-2.5*width,mean+2.5*width);
337 fs = hdedx_costh[i][j]->Fit("gaus","l","",mean-2.5*width,mean+2.5*width);
338 if( fs != 0 ){
339 std::cout << "\t\t" <<"MEAN FIT STATUS " << fs << std::endl;
340 fs1++;
341 }
342
343 satdedx = hdedx_costh[i][j]->GetFunction("gaus")->GetParameter(1);
344 saterror = hdedx_costh[i][j]->GetFunction("gaus")->GetParError(1);
345 means[i][j] = satdedx;
346 errors[i][j] = saterror;
347 widths[i][j] = hdedx_costh[i][j]->GetFunction("gaus")->GetParameter(2);
348
349 if( satdedx > dedxmax ) dedxmax = satdedx;
350 if( satdedx < dedxmin ) dedxmin = satdedx;
351
352 // fit the dE/dx distribution without correction in bins of beta-gamma and cosine
353 fs = hdedx_costh_pub[i][j]->Fit("gaus","ql");
354 mean = hdedx_costh_pub[i][j]->GetFunction("gaus")->GetParameter(1);
355 width = hdedx_costh_pub[i][j]->GetFunction("gaus")->GetParameter(2);
356 // fs = hdedx_costh_pub[i][j]->Fit("gaus","ql","",mean-2.5*width,mean+2.5*width);
357 fs = hdedx_costh_pub[i][j]->Fit("gaus","ql","",mean-2.5*width,mean+2.5*width);
358 if( fs != 0 ){
359 std::cout << "\t\t" <<"MEAN PUB FIT STATUS " << fs << std::endl;
360 fs2++;
361 }
362
363 satpubdedx = hdedx_costh_pub[i][j]->GetFunction("gaus")->GetParameter(1);
364 satpuberror = hdedx_costh_pub[i][j]->GetFunction("gaus")->GetParError(1);
365 satpubwidth = hdedx_costh_pub[i][j]->GetFunction("gaus")->GetParameter(2);
366
367
368 // fit the chi distribution
369 fs = hchi_costh[i][j]->Fit("gaus","ql");
370 mean = hchi_costh[i][j]->GetFunction("gaus")->GetParameter(1);
371 width = hchi_costh[i][j]->GetFunction("gaus")->GetParameter(2);
372 fs = hchi_costh[i][j]->Fit("gaus","ql","",mean-2.5*width,mean+2.5*width);
373 if( fs != 0 ){
374 // std::cout << "\t\t" <<"CHI FIT STATUS " << fs << std::endl;
375 fs3++;
376 }
377
378 satchi = hchi_costh[i][j]->GetFunction("gaus")->GetParameter(1);
379 satchierr = hchi_costh[i][j]->GetFunction("gaus")->GetParError(1);
380 satchiwidth = hchi_costh[i][j]->GetFunction("gaus")->GetParameter(2);
381
382
383 // fit the chi distribution without correction
384 /*
385 fs = hchi_costh_pub[i][j]->Fit("gaus","ql");
386 mean = hchi_costh_pub[i][j]->GetFunction("gaus")->GetParameter(1);
387 width = hchi_costh_pub[i][j]->GetFunction("gaus")->GetParameter(2);
388 fs = hchi_costh_pub[i][j]->Fit("gaus","ql","",mean-2.5*width,mean+2.5*width);
389 if( fs != 0 ){
390 // std::cout << "\t\t" <<"CHI PUB FIT STATUS " << fs << std::endl;
391 fs4++;
392 }
393
394 satchipub = hchi_costh_pub[i][j]->GetFunction("gaus")->GetParameter(1);
395 satchipuberr = hchi_costh_pub[i][j]->GetFunction("gaus")->GetParError(1);
396 satchipubwidth = hchi_costh_pub[i][j]->GetFunction("gaus")->GetParameter(2);
397 */
398
399 // determine the ratio of the predicted mean at a given bg to that of the average
400 ratio = m_gpar.dedxPrediction(satbg_avg)/m_gpar.dedxPrediction(satbg);
401
402 // fill the tree for this bin
403 satTree->Fill();
404 }
405 }
406
407 if( fs1+fs2+fs3+fs4 != 0 ) std::cout << "FIT RESULTS: " << particle << std::endl;
408 if( fs1 != 0 ) std::cout << "\t\t" <<"MEAN FIT FAILS " << fs1 << std::endl;
409 if( fs2 != 0 ) std::cout << "\t\t" <<"MEAN PUB FIT FAILS " << fs2 << std::endl;
410 if( fs3 != 0 ) std::cout << "\t\t" <<"CHI FIT FAILS " << fs3 << std::endl;
411 if( fs4 != 0 ) std::cout << "\t\t" <<"CHI PUB FIT FAILS " << fs4 << std::endl;
412
413 std::string corname("uncorrected");
414 if( correct == true ) corname = "corrected";
415
416 // Print the histograms for quality control
417 TCanvas* ctmp = new TCanvas("tmp","tmp",900,900);
418 ctmp->Divide(3,3);
419 std::stringstream psname; psname << "plots/dedx_" << particle << "_" << corname << ".ps[";
420 ctmp->Print(psname.str().c_str());
421 psname.str(""); psname << "plots/dedx_" << particle << "_" << corname << ".ps";
422 for( int i = 0 ; i < m_bgbins; ++i ){
423 for( int j = 0; j < m_cosbins; ++j ){
424 ctmp->cd(j%9+1);
425 hdedx_costh[i][j]->Draw();
426 if((j+1)%9==0)
427 ctmp->Print(psname.str().c_str());
428 }
429 }
430 psname.str(""); psname << "plots/dedx_" << particle << "_" << corname << ".ps]";
431 ctmp->Print(psname.str().c_str());
432 delete ctmp;
433
434
435 // Print the histograms for quality control
436 TCanvas* ctmp2 = new TCanvas("tmp2","tmp2",900,900);
437 ctmp2->Divide(3,3);
438 psname.str(""); psname << "plots/dedx_" << particle << "_pub_" << corname << ".ps[";
439 ctmp2->Print(psname.str().c_str());
440 psname.str(""); psname << "plots/dedx_" << particle << "_pub_" << corname << ".ps";
441 for( int i = 0 ; i < m_bgbins; ++i ){
442 for( int j = 0; j < m_cosbins; ++j ){
443 ctmp2->cd(j%9+1);
444 hdedx_costh_pub[i][j]->Draw();
445 if((j+1)%9==0)
446 ctmp2->Print(psname.str().c_str());
447 }
448 }
449 psname.str(""); psname << "plots/dedx_" << particle << "_pub_" << corname << ".ps]";
450 ctmp2->Print(psname.str().c_str());
451 delete ctmp2;
452
453
454 // Print the histograms for quality control
455 TCanvas* ctmp3 = new TCanvas("tmp3","tmp3",900,900);
456 ctmp3->Divide(3,3);
457 psname.str(""); psname << "plots/chi_" << particle << "_" << corname << ".ps[";
458 ctmp3->Print(psname.str().c_str());
459 psname.str(""); psname << "plots/chi_" << particle << "_" << corname << ".ps";
460 for( int i = 0 ; i < m_bgbins; ++i ){
461 for( int j = 0; j < m_cosbins; ++j ){
462 ctmp3->cd(j%9+1);
463 hchi_costh[i][j]->Draw();
464 if((j+1)%9==0)
465 ctmp3->Print(psname.str().c_str());
466 }
467 }
468 psname.str(""); psname << "plots/chi_" << particle << "_" << corname << ".ps]";
469 ctmp3->Print(psname.str().c_str());
470 delete ctmp3;
471
472
473 // Print the histograms for quality control
474 TCanvas* ctmp4 = new TCanvas("tmp4","tmp4",900,900);
475 ctmp4->Divide(3,3);
476 psname.str(""); psname << "plots/chi_" << particle << "_pub_" << corname << ".ps[";
477 ctmp4->Print(psname.str().c_str());
478 psname.str(""); psname << "plots/chi_" << particle << "_pub_" << corname << ".ps";
479 for( int i = 0 ; i < m_bgbins; ++i ){
480 for( int j = 0; j < m_cosbins; ++j ){
481 ctmp4->cd(j%9+1);
482 hchi_costh_pub[i][j]->Draw();
483 if((j+1)%9==0)
484 ctmp4->Print(psname.str().c_str());
485 }
486 }
487 psname.str(""); psname << "plots/chi_" << particle << "_pub_" << corname << ".ps]";
488 ctmp4->Print(psname.str().c_str());
489 delete ctmp4;
490
491 double cbcenters[m_cosbins], cberrors[m_cosbins];
492 for( int j = 0; j < m_cosbins; ++j ){
493 cbcenters[j] = m_lowercos+0.5*cosstep+j*cosstep;
494 cberrors[j] = 0;
495 }
496
497 // Plot the dE/dx means vs. cos(theta) for validation
498 TGraphErrors gdedx_costh[m_bgbins];
499 for( int i = 0; i < m_bgbins; ++i ){
500 char graphname[100];
501 satbg = m_lowerbg+0.5*bgstep+i*bgstep;
502 sprintf(graphname, "dedx_costh_p_%d",i);
503 gdedx_costh[i] = TGraphErrors(m_cosbins,cbcenters,means[i],cberrors,errors[i]);
504 }
505 TH1F* base = new TH1F("base","dE/dx vs. cos(#theta);cos(#theta);dE/dx",100,0,1.0);
506 base->GetXaxis()->SetTitleOffset(1.5);
507 base->SetMaximum(dedxmax*1.1);
508 base->SetMinimum(dedxmin*0.9);
509 base->SetStats(0);
510
511 // Print the histograms for quality control
512 TCanvas* ctmp6 = new TCanvas("tmp6","tmp6",900,900);
513 base->DrawCopy();
514 for( int i = 0 ; i < m_bgbins; ++i ){
515 gdedx_costh[i].SetMarkerSize(0.8);
516 gdedx_costh[i].SetMarkerColor(50+i);
517 gdedx_costh[i].SetLineColor(50+i);
518 gdedx_costh[i].DrawClone("same");
519 }
520 psname.str(""); psname << "plots/costh_" << particle << "_" << corname << ".eps";
521 ctmp6->SaveAs(psname.str().c_str());
522 delete ctmp6;
523
524 std::cout << "HadronPrep: saving output to " << outfile->GetName() << std::endl;
525
526 // write out the data to file
527 outfile->cd();
528 satTree->Write();
529
530 delete base;
531 for( int i = 0; i < m_bgbins; ++i ){
532 for( int j = 0; j < m_cosbins; ++j ){
533 delete hdedx_costh[i][j];
534 delete hdedx_costh_pub[i][j];
535 delete hchi_costh[i][j];
536 delete hchi_costh_pub[i][j];
537 }
538 }
539 infile->Close();
540}
541
542
543void
544HadronPrep::bgFits( TString particle, TFile* outfile, bool correct, std::string paramfile ){
545
546 // supress TCanvas::Print messages
547 gErrorIgnoreLevel = kWarning;
548
549 m_gpar.setParameters(paramfile);
550
551 double mass = 0.0;
552 if( particle == "pion" ) mass = Widget::mpion;
553 else if( particle == "kaon" ) mass = Widget::mkaon;
554 else if( particle == "proton" ) mass = Widget::mproton;
555 else if( particle == "muon" ) mass = Widget::mmuon;
556 else if( particle == "electron" ) mass = Widget::melectron;
557 if( mass == 0.0 ) exit(1);
558
559 TFile* infile = new TFile(m_filename);
560 TTree* hadron = (TTree*)infile->Get("track");
561
562 std::cout << "HadronPrep: reading in file " << m_filename << std::endl;
563 std::cout << "\t beta-gamma range: " << m_lowerbg << " to " << m_upperbg << std::endl;
564
565 // --------------------------------------------------
566 // INITIALIZE CONTAINERS
567 // --------------------------------------------------
568
569 double dedxpub; // dE/dx without electron saturation correction
570 double dedx; // dE/dx with electron saturation correction
571 double p; // track momentum
572 double bg; // track beta-gamma
573 double costh; // cosine of track polar angle
574 double db; // the nearest distance to the IP in the xy plane
575 double dz; // the nearest distance to the IP in the z plane
576 double chiPi; // PID chi value
577 double eop; // energy over momentum in the calorimeter
578
579 // Belle II variables
580 int b2nhit; // number of hits on this track
581
582 // BES III variables
583 double b3nhit; // number of hits on this track
584
585 hadron->SetBranchAddress("dedx", &dedxpub);
586 hadron->SetBranchAddress("dedxsat", &dedx);
587 hadron->SetBranchAddress("pF", &p);
588 hadron->SetBranchAddress("costh", &costh);
589 hadron->SetBranchAddress("db", &db);
590 hadron->SetBranchAddress("dz", &dz);
591 hadron->SetBranchAddress("chiPi", &chiPi);
592 hadron->SetBranchAddress("eopst", &eop);
593
594 if( m_type == 0 )
595 hadron->SetBranchAddress("numGoodHits", &b3nhit);
596 else if( m_type == 1 )
597 hadron->SetBranchAddress("lNHitsUsed", &b2nhit);
598 // hadron->SetBranchAddress("numGoodLayerHits", &b2nhit);
599
600 double bgstep = (m_upperbg-m_lowerbg)/m_bgbins;
601
602 // Create the histograms to be fit (before correction)
603 TH1F* hdedx_bg[m_bgbins];
604 TH1F* hdedx_bg_pub[m_bgbins];
605 TH1F* hchi_bg[m_bgbins];
606 TH1F* hchi_bg_pub[m_bgbins];
607 TH1F* hsigma_bg[m_bgbins];
608
609 // initialize the histograms
610 for( int i = 0; i < m_bgbins; ++i ){
611 char histname[100], histname2[100], histname3[100];
612 char histname4[100], histname5[100], histname6[100];
613 sprintf(histname, "dedx_bg_%d",i);
614 sprintf(histname2, "chi_bg_%d",i);
615 sprintf(histname3, "chipub_bg_%d",i);
616 sprintf(histname4, "dedxpub_bg_%d",i);
617 sprintf(histname5, "dedxnew_bg_%d",i);
618 sprintf(histname6, "sigma_bg_%d",i);
619
620 if( particle == "electron" ){
621 hdedx_bg[i] = new TH1F(histname,histname,200,0,2);
622 hchi_bg[i] = new TH1F(histname2,histname2,200,-2,2);
623 hchi_bg_pub[i] = new TH1F(histname3,histname3,200,-5,5);
624 hdedx_bg_pub[i] = new TH1F(histname4,histname4,200,0,2);
625 }
626 else if( particle == "proton" ){
627 if( i < 5 ){
628 hdedx_bg[i] = new TH1F(histname,histname,200,0,6);
629 hchi_bg[i] = new TH1F(histname2,histname2,200,-30,30);
630 hchi_bg_pub[i] = new TH1F(histname3,histname3,200,0,40);
631 hdedx_bg_pub[i] = new TH1F(histname4,histname4,200,0,6);
632 }
633 else if( i < 9 ){
634 hdedx_bg[i] = new TH1F(histname,histname,200,0,4);
635 hchi_bg[i] = new TH1F(histname2,histname2,200,-5,5);
636 hchi_bg_pub[i] = new TH1F(histname3,histname3,200,-20,20);
637 hdedx_bg_pub[i] = new TH1F(histname4,histname4,200,0,4);
638 }
639 else{
640 hdedx_bg[i] = new TH1F(histname,histname,200,0,2);
641 hchi_bg[i] = new TH1F(histname2,histname2,200,-5,5);
642 hchi_bg_pub[i] = new TH1F(histname3,histname3,200,-5,5);
643 hdedx_bg_pub[i] = new TH1F(histname4,histname4,200,0,2);
644 }
645 }
646 else if( i < 2 ){
647 hdedx_bg[i] = new TH1F(histname,histname,200,0,4);
648 hchi_bg[i] = new TH1F(histname2,histname2,200,-5,5);
649 hchi_bg_pub[i] = new TH1F(histname3,histname3,200,-20,20);
650 hdedx_bg_pub[i] = new TH1F(histname4,histname4,200,0,4);
651 }
652 else{
653 hdedx_bg[i] = new TH1F(histname,histname,200,0,2);
654 hchi_bg[i] = new TH1F(histname2,histname2,200,-5,5);
655 hchi_bg_pub[i] = new TH1F(histname3,histname3,100,-5,5);
656 hdedx_bg_pub[i] = new TH1F(histname4,histname4,200,0,2);
657 }
658 hsigma_bg[i] = new TH1F(histname6,histname6,100,-3,3);
659 }
660
661 // create some containers to calculate averages
662 double sumcos[m_bgbins];
663 double sumsin[m_bgbins];
664 double sumbg[m_bgbins];
665 double sumressq[m_bgbins];
666 int sumsize[m_bgbins];
667 for( int i = 0; i < m_bgbins; ++i ){
668 sumcos[i] = 0;
669 sumsin[i] = 0;
670 sumbg[i] = 0;
671 sumressq[i] = 0;
672 sumsize[i] = 0;
673 }
674
675 // create some containers for means and errors
676 double means[m_bgbins];
677 double errors[m_bgbins];
678 double widths[m_bgbins];
679 for( int i = 0; i < m_bgbins; ++i ){
680 means[i] = 0;
681 errors[i] = 0;
682 widths[i] = 0;
683 }
684
685 // create some containers for checking the residual saturation vs. cos(theta)
686 const int cosbins = 20;
687 const double cosstep = 2.0 / cosbins;
688 double cosArray[cosbins], cosArrayErr[cosbins];
689 for( int i = 0; i < cosbins; ++i ){
690 cosArray[i] = -1 + (i*cosstep + cosstep/2.0);
691 cosArrayErr[i] = 0.0;
692 }
693
694 TH1F* hchi_costh_all[2][cosbins];
695 TH1F* hchi_costh_0[2][cosbins];
696 TH1F* hchi_costh_1[2][cosbins];
697 TH1F* hchi_costh_2[2][cosbins];
698
699 for( int i = 0; i < cosbins; ++i ){
700 char histname[100], histname0[100], histname1[100], histname2[100];
701 sprintf(histname, "chi_costh_pos_%d",i);
702 sprintf(histname0, "chi_costh0_pos_%d",i);
703 sprintf(histname1, "chi_costh1_pos_%d",i);
704 sprintf(histname2, "chi_costh2_pos_%d",i);
705
706 hchi_costh_all[0][i] = new TH1F(histname,histname,100,-5,5);
707 hchi_costh_0[0][i] = new TH1F(histname0,histname0,100,-5,5);
708 hchi_costh_1[0][i] = new TH1F(histname1,histname1,100,-5,5);
709 hchi_costh_2[0][i] = new TH1F(histname2,histname2,100,-5,5);
710
711 sprintf(histname, "chi_costh_neg_%d",i);
712 sprintf(histname0, "chi_costh0_neg_%d",i);
713 sprintf(histname1, "chi_costh1_neg_%d",i);
714 sprintf(histname2, "chi_costh2_neg_%d",i);
715
716 hchi_costh_all[1][i] = new TH1F(histname,histname,100,-5,5);
717 hchi_costh_0[1][i] = new TH1F(histname0,histname0,100,-5,5);
718 hchi_costh_1[1][i] = new TH1F(histname1,histname1,100,-5,5);
719 hchi_costh_2[1][i] = new TH1F(histname2,histname2,100,-5,5);
720 }
721
722 // get the hadron saturation parameters
723 // if the parameters do not exist, use the values in the default constructor
724 m_hadsat.setParameters("sat-pars.fit.txt");
725
726 // --------------------------------------------------
727 // LOOP OVER EVENTS AND FILL CONTAINERS
728 // --------------------------------------------------
729
730 // Fill the histograms to be fitted
731 for( unsigned int index = 0; index < hadron->GetEntries(); ++index ){
732 hadron->GetEvent(index);
733
734 int charge = (p < 0)? 1 : 0;
735 bg = fabs(p)/mass;
736
737 int nhit;
738 if( m_type == 0 )
739 nhit = std::floor(b3nhit);
740 else if( m_type == 1 )
741 nhit = b2nhit;
742
743 // clean up bad events and restrict the momentum range
744 if( nhit < 0 || nhit > 100 || dedx <= 0 || costh != costh ||
745 bg <= m_lowerbg || bg >= m_upperbg )
746 continue;
747
748 // apply an E/p cut for pions
749 if( particle == "pion" && eop > 0.75 ) continue;
750
751 // use loose dE/dx restrictions to remove contaminations
752 if( particle == "proton" && (dedxpub < 0 || dedxpub > 100 || dedxpub < 0.85/fabs(p)) )
753 continue;
754 if( particle == "kaon" && dedxpub < 0.4/fabs(p) )
755 continue;
756
757 int bgBin = (int)((bg-m_lowerbg)/(m_upperbg-m_lowerbg) * m_bgbins);
758
759 double dedx_pre = dedx;
760 double dedx_new = dedx_pre;
761 if( !m_mcFlag ) dedx_new = m_hadsat.D2I(costh,m_hadsat.I2D(costh,1.0)*dedx);
762 double dedx_cur = m_gpar.dedxPrediction(bg);
763 double dedx_res = m_gpar.sigmaPrediction(dedx_cur,nhit,sqrt(1-costh*costh));
764
765 double chi_new = (dedx_new-dedx_cur)/dedx_res;
766 if( particle == "electron" ) chi_new = (dedx_new-1)/dedx_res;
767 double res_cor = m_gpar.resPrediction(nhit,sqrt(1-costh*costh));
768 int i = (int) ( (fabs(bg)-m_lowerbg)/bgstep );
769 hsigma_bg[i]->Fill((dedx_new-dedx_cur)/res_cor);
770
771 if( correct ) hdedx_bg[bgBin]->Fill(dedx_new);
772 else hdedx_bg[bgBin]->Fill(dedx_pre);
773 hdedx_bg_pub[bgBin]->Fill(dedxpub);
774 if( correct ) hchi_bg[bgBin]->Fill(chi_new);
775 hchi_bg_pub[bgBin]->Fill(chiPi);
776
777 sumcos[bgBin] += costh;
778 sumsin[bgBin] += sqrt(1-costh*costh);
779 sumbg[bgBin] += bg;
780 sumressq[bgBin] += pow(dedx_res,2);
781 sumsize[bgBin] += 1;
782
783 // make histograms of dE/dx vs. cos(theta) for validation
784 int icos = (int)((costh+1)/cosstep);
785 hchi_costh_all[charge][icos]->Fill(chi_new);
786 if( bgBin < int(m_bgbins/3) )
787 hchi_costh_0[charge][icos]->Fill(chi_new);
788 else if( bgBin < 2*int(m_bgbins/3) )
789 hchi_costh_1[charge][icos]->Fill(chi_new);
790 else
791 hchi_costh_2[charge][icos]->Fill(chi_new);
792
793 }// end of event loop
794
795 // --------------------------------------------------
796 // FIT IN BINS OF BETA-GAMMA
797 // --------------------------------------------------
798
799 // fit the histograms with Gaussian functions
800 // and extract the means and errors
801 TTree* satTree = new TTree(particle,"dE/dx means and errors");
802
803 double satbg; // beta-gamma value for this bin
804 double satcosth; // cos(theta) value for this bin
805 double satdedx; // mean dE/dx value for this bin
806 double saterror; // error on ^
807 double satbg_avg; // average beta-gamma value for this sample
808 double satcosth_avg; // average cos(theta) value for this sample
809 double satsinth_avg; // average sin(theta) value for this sample
810 double satdedxres_avg; // average dE/dx error squared for this sample
811 double satpubdedx; // mean "public" dE/dx value for this bin
812 double satpuberror; // error on ^
813 double satpubwidth; // width of ^ distribution
814 double satchi; // mean chi value for this bin
815 double satchierr; // error on ^
816 double satchiwidth; // width of ^ distribution
817 double satchipub; // mean "public" chi value for this bin
818 double satchipuberr; // error on ^
819 double satchipubwidth; // width of ^ distribution
820 double ratio; // ratio of the predicted mean to that of the average
821
822 satTree->Branch("bg",&satbg,"bg/D");
823 satTree->Branch("costh",&satcosth,"costh/D");
824 satTree->Branch("dedx",&satdedx,"dedx/D");
825 satTree->Branch("error",&saterror,"error/D");
826 satTree->Branch("bg_avg",&satbg_avg,"bg_avg/D");
827 satTree->Branch("costh_avg",&satcosth_avg,"costh_avg/D");
828 satTree->Branch("sinth_avg",&satsinth_avg,"sinth_avg/D");
829 satTree->Branch("dedxres_avg",&satdedxres_avg,"dedxres_avg/D");
830 satTree->Branch("dedx_pub",&satpubdedx,"dedx_pub/D");
831 satTree->Branch("dedxerr_pub",&satpuberror,"dedxerr_pub/D");
832 satTree->Branch("chi",&satchi,"chi/D");
833 satTree->Branch("chi_pub",&satchipub,"chi_pub/D");
834 satTree->Branch("sigma",&satchiwidth,"sigma/D");
835 satTree->Branch("sigma_pub",&satchipubwidth,"sigma_pub/D");
836 satTree->Branch("ratio",&ratio,"ratio/D");
837
838 // Fit the histograms
839 for( int i = 0; i < m_bgbins; ++i ){
840
841 // fill some details for this bin
842 satbg = m_lowerbg+0.5*bgstep+i*bgstep;
843
844 satbg_avg = sumbg[i]/sumsize[i];
845 satcosth_avg = sumcos[i]/sumsize[i];
846 satsinth_avg = sumsin[i]/sumsize[i];
847 // satdedxres_avg = sumressq[i]/sumsize[i];
848
849 hsigma_bg[i]->Fit("gaus","ql");
850 satdedxres_avg = hsigma_bg[i]->GetFunction("gaus")->GetParameter(2);
851
852 int fs;
853
854 // fit the dE/dx distribution in bins of beta-gamma
855 fs = hdedx_bg[i]->Fit("gaus","ql");
856 if( fs != 0 ) std::cout << "\t\t" <<"MEAN FIT STATUS " << fs << std::endl;
857 double mean = hdedx_bg[i]->GetFunction("gaus")->GetParameter(1);
858 double width = hdedx_bg[i]->GetFunction("gaus")->GetParameter(2);
859 fs = hdedx_bg[i]->Fit("gaus","ql","",mean-2.5*width,mean+2.5*width);
860 if( fs != 0 ) std::cout << "\t\t" <<"MEAN FIT STATUS " << fs << std::endl;
861
862 satdedx = hdedx_bg[i]->GetFunction("gaus")->GetParameter(1);
863 saterror = hdedx_bg[i]->GetFunction("gaus")->GetParError(1);
864 means[i] = satdedx;
865 errors[i] = saterror;
866 widths[i] = hdedx_bg[i]->GetFunction("gaus")->GetParameter(2);
867
868
869 // fit the dE/dx distribution without correction in bins of beta-gamma
870 fs = hdedx_bg_pub[i]->Fit("gaus","ql");
871 if( fs != 0 ) std::cout << "\t\t" <<"MEAN PUB FIT STATUS " << fs << std::endl;
872 mean = hdedx_bg_pub[i]->GetFunction("gaus")->GetParameter(1);
873 width = hdedx_bg_pub[i]->GetFunction("gaus")->GetParameter(2);
874 fs = hdedx_bg_pub[i]->Fit("gaus","ql","",mean-2.5*width,mean+2.5*width);
875 if( fs != 0 ) std::cout << "\t\t" <<"MEAN PUB FIT STATUS " << fs << std::endl;
876
877 satpubdedx = hdedx_bg_pub[i]->GetFunction("gaus")->GetParameter(1);
878 satpuberror = hdedx_bg_pub[i]->GetFunction("gaus")->GetParError(1);
879 satpubwidth = hdedx_bg_pub[i]->GetFunction("gaus")->GetParameter(2);
880
881
882 // fit the chi distribution
883 fs = hchi_bg[i]->Fit("gaus","ql");
884 if( fs != 0 ) std::cout << "\t\t" <<"CHI FIT STATUS " << fs << std::endl;
885 mean = hchi_bg[i]->GetFunction("gaus")->GetParameter(1);
886 width = hchi_bg[i]->GetFunction("gaus")->GetParameter(2);
887 fs = hchi_bg[i]->Fit("gaus","ql","",mean-2.5*width,mean+2.5*width);
888 if( fs != 0 ) std::cout << "\t\t" <<"CHI FIT STATUS " << fs << std::endl;
889
890 satchi = hchi_bg[i]->GetFunction("gaus")->GetParameter(1);
891 satchierr = hchi_bg[i]->GetFunction("gaus")->GetParError(1);
892 satchiwidth = hchi_bg[i]->GetFunction("gaus")->GetParameter(2);
893
894
895 // fit the chi distribution without correction
896 /*
897 fs = hchi_bg_pub[i]->Fit("gaus","ql");
898 if( fs != 0 ) std::cout << "\t\t" <<"CHI PUB FIT STATUS " << fs << std::endl;
899 mean = hchi_bg_pub[i]->GetFunction("gaus")->GetParameter(1);
900 width = hchi_bg_pub[i]->GetFunction("gaus")->GetParameter(2);
901 fs = hchi_bg_pub[i]->Fit("gaus","ql","",mean-2.5*width,mean+2.5*width);
902 if( fs != 0 ) std::cout << "\t\t" <<"CHI PUB FIT STATUS " << fs << std::endl;
903
904 satchipub = hchi_bg_pub[i]->GetFunction("gaus")->GetParameter(1);
905 satchipuberr = hchi_bg_pub[i]->GetFunction("gaus")->GetParError(1);
906 satchipubwidth = hchi_bg_pub[i]->GetFunction("gaus")->GetParameter(2);
907 */
908
909 // determine the ratio of the predicted mean at a given bg to that of the average
910 ratio = m_gpar.dedxPrediction(satbg_avg)/m_gpar.dedxPrediction(satbg);
911
912 // fill the tree for this bin
913 satTree->Fill();
914 }
915
916 std::string corname("uncorrected");
917 if( correct == true ) corname = "corrected";
918
919 // Print the histograms for quality control
920 TCanvas* ctmp = new TCanvas("tmp","tmp",900,900);
921 ctmp->Divide(3,3);
922 std::stringstream psname; psname << "plots/dedx_" << particle << "_" << corname << ".ps[";
923 ctmp->Print(psname.str().c_str());
924 psname.str(""); psname << "plots/dedx_" << particle << "_" << corname << ".ps";
925 for( int i = 0 ; i < m_bgbins; ++i ){
926 ctmp->cd(i%9+1);
927 hdedx_bg[i]->Draw();
928 if((i+1)%9==0)
929 ctmp->Print(psname.str().c_str());
930 }
931 psname.str(""); psname << "plots/dedx_" << particle << "_" << corname << ".ps]";
932 ctmp->Print(psname.str().c_str());
933 delete ctmp;
934
935 std::cout << "HadronPrep: saving output to " << outfile->GetName() << std::endl;
936
937 // write out the data to file
938 outfile->cd();
939 satTree->Write();
940
941 std::cout << "making validation plots..." << std::endl;
942
943 // --------------------------------------------------
944 // FIT IN BINS OF COS(THETA) FOR VALIDATION
945 // --------------------------------------------------
946
947 double chicos[2][cosbins], sigmacos[2][cosbins];
948 double chicoserr[2][cosbins], sigmacoserr[2][cosbins];
949
950 double chicos0[2][cosbins], sigmacos0[2][cosbins];
951 double chicos0err[2][cosbins], sigmacos0err[2][cosbins];
952
953 double chicos1[2][cosbins], sigmacos1[2][cosbins];
954 double chicos1err[2][cosbins], sigmacos1err[2][cosbins];
955
956 double chicos2[2][cosbins], sigmacos2[2][cosbins];
957 double chicos2err[2][cosbins], sigmacos2err[2][cosbins];
958
959 for( int c = 0; c < 2; ++c ){
960 for( int i = 0; i < cosbins; ++i ){
961 if( hchi_costh_all[c][i]->Integral(1,100) == 0 ) continue;
962 hchi_costh_all[c][i]->Fit("gaus","ql");
963 chicos[c][i] = hchi_costh_all[c][i]->GetFunction("gaus")->GetParameter(1);
964 chicoserr[c][i] = hchi_costh_all[c][i]->GetFunction("gaus")->GetParError(1);
965 sigmacos[c][i] = hchi_costh_all[c][i]->GetFunction("gaus")->GetParameter(2);
966 sigmacoserr[c][i] = hchi_costh_all[c][i]->GetFunction("gaus")->GetParError(2);
967 }
968 for( int i = 0; i < cosbins; ++i ){
969 if( hchi_costh_0[c][i]->Integral(1,100) == 0 ) continue;
970 hchi_costh_0[c][i]->Fit("gaus","ql");
971 chicos0[c][i] = hchi_costh_0[c][i]->GetFunction("gaus")->GetParameter(1);
972 chicos0err[c][i] = hchi_costh_0[c][i]->GetFunction("gaus")->GetParError(1);
973 sigmacos0[c][i] = hchi_costh_0[c][i]->GetFunction("gaus")->GetParameter(2);
974 sigmacos0err[c][i] = hchi_costh_0[c][i]->GetFunction("gaus")->GetParError(2);
975 }
976 for( int i = 0; i < cosbins; ++i ){
977 if( hchi_costh_1[c][i]->Integral(1,100) == 0 ) continue;
978 hchi_costh_1[c][i]->Fit("gaus","ql");
979 chicos1[c][i] = hchi_costh_1[c][i]->GetFunction("gaus")->GetParameter(1);
980 chicos1err[c][i] = hchi_costh_1[c][i]->GetFunction("gaus")->GetParError(1);
981 sigmacos1[c][i] = hchi_costh_1[c][i]->GetFunction("gaus")->GetParameter(2);
982 sigmacos1err[c][i] = hchi_costh_1[c][i]->GetFunction("gaus")->GetParError(2);
983 }
984 for( int i = 0; i < cosbins; ++i ){
985 if( hchi_costh_2[c][i]->Integral(1,100) == 0 ) continue;
986 hchi_costh_2[c][i]->Fit("gaus","ql");
987 chicos2[c][i] = hchi_costh_2[c][i]->GetFunction("gaus")->GetParameter(1);
988 chicos2err[c][i] = hchi_costh_2[c][i]->GetFunction("gaus")->GetParError(1);
989 sigmacos2[c][i] = hchi_costh_2[c][i]->GetFunction("gaus")->GetParameter(2);
990 sigmacos2err[c][i] = hchi_costh_2[c][i]->GetFunction("gaus")->GetParError(2);
991 }
992 }
993
994 // Print the histograms for quality control
995 TCanvas* ctmp2 = new TCanvas("tmp2","tmp2",900,900);
996 ctmp2->Divide(3,3);
997 psname.str(""); psname << "plots/chivscos_fits_" << particle << ".ps[";
998 ctmp2->Print(psname.str().c_str());
999 psname.str(""); psname << "plots/chivscos_fits_" << particle << ".ps";
1000 for( int i = 0 ; i < cosbins; ++i ){
1001 ctmp2->cd(i%9+1);
1002 hchi_costh_all[0][i]->Draw();
1003 hchi_costh_all[1][i]->SetMarkerColor(kRed);
1004 hchi_costh_all[1][i]->Draw("same");
1005 if((i+1)%9==0)
1006 ctmp2->Print(psname.str().c_str());
1007 }
1008 psname.str(""); psname << "plots/chivscos_fits_" << particle << ".ps]";
1009 ctmp2->Print(psname.str().c_str());
1010 delete ctmp2;
1011
1012 TGraphErrors* grchicos = new TGraphErrors(cosbins,cosArray,chicos[0],cosArrayErr,chicoserr[0]);
1013 TGraphErrors* grchicos0 = new TGraphErrors(cosbins,cosArray,chicos0[0],cosArrayErr,chicos0err[0]);
1014 TGraphErrors* grchicos1 = new TGraphErrors(cosbins,cosArray,chicos1[0],cosArrayErr,chicos1err[0]);
1015 TGraphErrors* grchicos2 = new TGraphErrors(cosbins,cosArray,chicos2[0],cosArrayErr,chicos2err[0]);
1016
1017 TGraphErrors* grchicosn = new TGraphErrors(cosbins,cosArray,chicos[1],cosArrayErr,chicoserr[1]);
1018 TGraphErrors* grchicos0n = new TGraphErrors(cosbins,cosArray,chicos0[1],cosArrayErr,chicos0err[1]);
1019 TGraphErrors* grchicos1n = new TGraphErrors(cosbins,cosArray,chicos1[1],cosArrayErr,chicos1err[1]);
1020 TGraphErrors* grchicos2n = new TGraphErrors(cosbins,cosArray,chicos2[1],cosArrayErr,chicos2err[1]);
1021
1022 TGraphErrors* grsigmacos = new TGraphErrors(cosbins,cosArray,sigmacos[0],cosArrayErr,sigmacoserr[0]);
1023 TGraphErrors* grsigmacos0 = new TGraphErrors(cosbins,cosArray,sigmacos0[0],cosArrayErr,sigmacos0err[0]);
1024 TGraphErrors* grsigmacos1 = new TGraphErrors(cosbins,cosArray,sigmacos1[0],cosArrayErr,sigmacos1err[0]);
1025 TGraphErrors* grsigmacos2 = new TGraphErrors(cosbins,cosArray,sigmacos2[0],cosArrayErr,sigmacos2err[0]);
1026
1027 TGraphErrors* grsigmacosn = new TGraphErrors(cosbins,cosArray,sigmacos[1],cosArrayErr,sigmacoserr[1]);
1028 TGraphErrors* grsigmacos0n = new TGraphErrors(cosbins,cosArray,sigmacos0[1],cosArrayErr,sigmacos0err[1]);
1029 TGraphErrors* grsigmacos1n = new TGraphErrors(cosbins,cosArray,sigmacos1[1],cosArrayErr,sigmacos1err[1]);
1030 TGraphErrors* grsigmacos2n = new TGraphErrors(cosbins,cosArray,sigmacos2[1],cosArrayErr,sigmacos2err[1]);
1031
1032 TLine* line0 = new TLine(-1,0,1,0);
1033 line0->SetLineStyle(kDashed);
1034 line0->SetLineColor(kRed);
1035 TLine* line1 = new TLine(-1,1,1,1);
1036 line1->SetLineStyle(kDashed);
1037 line1->SetLineColor(kRed);
1038
1039 TString ptype("");
1040 if( particle == "pion" ) ptype += "#pi";
1041 else if( particle == "kaon" ) ptype += "K";
1042 else if( particle == "proton" ) ptype += "p";
1043 else if( particle == "muon" ) ptype += "#mu";
1044 else if( particle == "electron" ) ptype += "e";
1045
1046 TLegend* lchi = new TLegend(0.6,0.75,0.8,0.85);
1047 lchi->SetBorderSize(0);
1048 lchi->SetFillColor(0);
1049 lchi->AddEntry(grchicos,ptype+"^{+}","p");
1050 lchi->AddEntry(grchicosn,ptype+"^{-}","p");
1051 TCanvas* cchi = new TCanvas("cchi","cchi",700,600);
1052
1053 FormatGraph(grchicos,0);
1054 FormatGraph(grchicosn,1);
1055 grchicos->Draw("AP");
1056 grchicosn->Draw("P,same");
1057 line0->Draw("same");
1058 lchi->Draw("same");
1059 cchi->SaveAs("plots/chiall"+particle+".eps");
1060
1061 FormatGraph(grchicos0,0);
1062 FormatGraph(grchicos0n,1);
1063 grchicos0->Draw("AP");
1064 grchicos0n->Draw("P,same");
1065 line0->Draw("same");
1066 lchi->Draw("same");
1067 cchi->SaveAs("plots/chi0"+particle+".eps");
1068
1069 FormatGraph(grchicos1,0);
1070 FormatGraph(grchicos1n,1);
1071 grchicos1->Draw("AP");
1072 grchicos1n->Draw("P,same");
1073 line0->Draw("same");
1074 lchi->Draw("same");
1075 cchi->SaveAs("plots/chi1"+particle+".eps");
1076
1077 FormatGraph(grchicos2,0);
1078 FormatGraph(grchicos2n,1);
1079 grchicos2->Draw("AP");
1080 grchicos2n->Draw("P,same");
1081 line0->Draw("same");
1082 lchi->Draw("same");
1083 cchi->SaveAs("plots/chi2"+particle+".eps");
1084
1085 FormatGraph(grsigmacos,2);
1086 FormatGraph(grsigmacosn,1);
1087 grsigmacos->Draw("AP");
1088 grsigmacosn->Draw("P,same");
1089 line1->Draw("same");
1090 lchi->Draw("same");
1091 cchi->SaveAs("plots/sigmaall"+particle+".eps");
1092
1093 FormatGraph(grsigmacos0,2);
1094 FormatGraph(grsigmacos0n,1);
1095 grsigmacos0->Draw("AP");
1096 grsigmacos0n->Draw("P,same");
1097 line1->Draw("same");
1098 lchi->Draw("same");
1099 cchi->SaveAs("plots/sigma0"+particle+".eps");
1100
1101 FormatGraph(grsigmacos1,2);
1102 FormatGraph(grsigmacos1n,1);
1103 grsigmacos1->Draw("AP");
1104 grsigmacos1n->Draw("P,same");
1105 line1->Draw("same");
1106 lchi->Draw("same");
1107 cchi->SaveAs("plots/sigma1"+particle+".eps");
1108
1109 FormatGraph(grsigmacos2,2);
1110 FormatGraph(grsigmacos2n,1);
1111 grsigmacos2->Draw("AP");
1112 grsigmacos2n->Draw("P,same");
1113 line1->Draw("same");
1114 lchi->Draw("same");
1115 cchi->SaveAs("plots/sigma2"+particle+".eps");
1116
1117 for( int i = 0; i < 2; ++i ){
1118 for( int j = 0; j < cosbins; ++j ){
1119 delete hchi_costh_all[i][j];
1120 delete hchi_costh_0[i][j];
1121 delete hchi_costh_1[i][j];
1122 delete hchi_costh_2[i][j];
1123 }
1124 }
1125
1126 delete grchicos;
1127 delete grchicos0;
1128 delete grchicos1;
1129 delete grchicos2;
1130
1131 delete grchicosn;
1132 delete grchicos0n;
1133 delete grchicos1n;
1134 delete grchicos2n;
1135
1136 delete grsigmacos;
1137 delete grsigmacos0;
1138 delete grsigmacos1;
1139 delete grsigmacos2;
1140
1141 delete grsigmacosn;
1142 delete grsigmacos0n;
1143 delete grsigmacos1n;
1144 delete grsigmacos2n;
1145
1146 delete line0;
1147 delete line1;
1148
1149 delete lchi;
1150 delete cchi;
1151}
mass
double mass
Definition: CosmicGenerator.cxx:138
Fill
curve Fill()
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)
gr
TGraph * gr
Definition: DataBase/tau_mode.c:69
HadronPrep::FormatGraph
void FormatGraph(TGraphErrors *gr, int flag)
Definition: HadronPrep.cc:28
HadronPrep::bgCosThetaFits
void bgCosThetaFits(TString particle, TFile *outfile, bool correct, std::string paramfile)
Definition: HadronPrep.cc:55
HadronPrep::bgFits
void bgFits(TString particle, TFile *outfile, bool correct, std::string paramfile)
Definition: HadronPrep.cc:544
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::setParameters
void setParameters(std::string infile)
Definition: WidgetParameterization.cc:12
WidgetParameterization::resPrediction
double resPrediction(double nhit, double sin)
Definition: WidgetParameterization.cc:164
WidgetParameterization::sigmaPrediction
double sigmaPrediction(double dedx, double nhit, double sin)
Definition: WidgetParameterization.cc:130
charge
float charge
Definition: prepare_data.cxx:28
bg
float bg
Definition: prepare_data.cxx:28