attenuation_calibration_comparisons.C
Go to the documentation of this file.
1 //Author: Prabhjot Singh (prabhjot@fnal.gov).
2 //Dated April 19 2016 03:00 am
3 //This program is a template macro to plot and compare
4 //a. uncorrected PE/cm profiles,
5 //b. threshold and shielding corrections,
6 //c. attenution fitting plots for any randomly chosen plane and cell,
7 //d. calibration constants, coeff_exp, attenlength, background and chisq,
8 //e. list and number of calibrated channels for the following cases
9 //1. ND MC,
10 //2. ND Data,
11 //3. FD MC,
12 //4. FD Data,
13 //5. Tricell vs. Trajectory (two cases in one time)
14 //6. Any two cases in one time e.g.
15 // a. different periods in same analysis,
16 // b. new calibration vs. old calibration,
17 // c. one analyis vs. another
18 //7. Many vs. one case, e.g. all 2nd analysis cases vs. 1st analysis e.g. period1, period2, epoch3b, epoch3c vs. 1st analysis.
19 //
20 //8. You have to set the booleans given at the start of the program to select one or more cases.
21 //
22
23 #include "HistogramAttr.h"
25
26 bool uncorrectedPEplots = false; // set true if you want plots of uncorrected PE/cm vs. W
27 bool thresholdshadowingplots = false; // set true if you want plots of threshold and shadowing correction factors.
28 bool correctedPEplots = true; // set true if you want plots of corrected PE/cm vs. W. That is, attenuation fitting plots and thier ratios.
29 bool calibrationconstantsplots = false; // set it to true if you wans to make plots of calibration constants
30
31 bool nd = true; // set true for the ND comparisons
32 bool fd = false; // set true for the FD comparisons
33 bool mc = false; // set true for the Monte Carlo based comparisons
34 bool data = true; // set true for the Data based comparisons
35 bool tricellvstrajectory = false; // set true for the tricell vs. trajectory comparisons
36 bool onetoone = true; // set true for one vs. one case e.g. point 6 above
37 bool manytoone = false; // set true for many vs. one case e..g. point 7 above
38 bool nddatamuoncatcher = true; // set true for the ND data muon catcher plots
39
41  gStyle->SetOptStat(" ");
42  //-----------------------------------------------------------------------------------------------------
43  //set of all .root files for uncorrected PE/cm vs. w plots
44  char* c1 = "/nova/ana/users/prabhjot/gridjob/Work_with_Adam_and_Brian/nd/NDPCHitData_S150224/mainbody_and_muoncatcher/mainboby/ConsolidateViewsMode_false/Fromgrid_aftercosmiccalib/histofiles/calibNDPCHitData_hist_mainbody_ConsolidateViewsMode_false_S150224.root"; // base file
45  char* c2 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/nddata/Period2_DDActivity/mainbody/calib_hist_2ndAna_fulldetector_nddata_mainbody_period2ddactivity.root"; // 2nd file
46
47  //.root file for uncorrected PE/cm for tricell vs trajectory study
48  char* c3 = "";
49
50  //set of .root files for threshold and shadowing corrections
51  char* cthreshold1 = ""; // base file
52  char* cthreshold2 = ""; // 2nd file
53
54  //set of .root files to plot attenuation fitting plots
55  char* cattenuation1 = "/nova/app/users/prabhjot/nova_test_S150224/AttenuationCalibration/nd/pchitdata/fromgrid/mainbody_and_muoncatcher/mainboby/ConsolidateViewsMode_false/Atten_Fit_Results.root"; // base file
56  char* cattenuation2 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/nddata/Period2_DDActivity/mainbody/FitAttenuation/Atten_Fit_Results.root"; // 2nd file
57
58  //set of .root files to plot attenuation fitting plots for the ND data muon catcher
59  char* cattenuation_muoncatcher1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/nddata/Period2_DDActivity/muoncatcher/FitAttenuation/Atten_Fit_Results.root"; // base file
60  char* cattenuation_muoncatcher2 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/nddata/Epoch3b_DDActivity/muoncatcher/FitAttenuation/Atten_Fit_Results.root"; // 2nd file
61
62  //set of attenuation constants .txt files for attenuation constants comparisons
63  char* cconstants1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/nddata/comparison_plots/calib_atten_consts_uniq_sort_1st_ana.txt"; // base file
64  char* cconstants2 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/nddata/comparison_plots/calib_atten_consts_uniq_sort_period2.txt"; // 2nd file
65  //---------------------------------------------------------------------------------------------------
66  //---------------------------------------------------------------------------------------------------
67
68  TString hist1; // base histogram from the base file. You want to compare other histograms wrt base histogram.
69  TString hist2; // histogram from the second file that you want to compare with the base histogram.
70  TString det; // name of the detector.
71  TString mcdata; // labeling for the Monte Carlo or the data .
72  TString analysis1; // labeling for the base file.
73  TString analysis2; // labelling for the second file.
74  TString var; // variable labelling on the x axis. It is W (cm).
75  TString plane; // plane number.
76  TString cell; // cell number.
77
78  analysis1 = "2nd_analysis_period2";
79  analysis2 = "2nd_analysis_epoch3b";
80
81  hist1 = "h_WPE_corr_xy";
82  hist2 = "h_WPE_corr_xy"; //_corr_xy_ for the tricell, _corr_traj_ for the trajectory
83
84  var = "W";
85  plane = "003";
86  cell = "023";
87
88  //ND MC two cases comparison
89  if(nd && onetoone){
90  det = "nd";
91  if(mc) mcdata = "mc";
92  if(data) mcdata = "data";
93
94  if(tricellvstrajectory && uncorrectedPEplots) uncorrected_PE_tricellvstrajectory(c3, det, mcdata, analysis1, analysis2, hist1, hist2);
95  else if(uncorrectedPEplots) uncorrected_PE(c1, c2, det, mcdata, analysis1, analysis2, hist1, hist2);
96
98
99  if(data && nddatamuoncatcher && correctedPEplots)corrected_PE(cattenuation_muoncatcher1, cattenuation_muoncatcher2, det, mcdata, analysis1, analysis2, plane, cell);
100  else if(correctedPEplots) corrected_PE(cattenuation1, cattenuation2, det, mcdata, analysis1, analysis2, plane, cell);
101
102  if(calibrationconstantsplots) calibrationconstants(cconstants1, cconstants2, det, mcdata, analysis1, analysis2);
103  }//end of if(nd && mc && onetoone)
104
105  //FD two case comparison
106  if(fd && onetoone){
107  det = "fd";
108  if(mc) mcdata = "mc";
109  if(data) mcdata = "data";
110
111  if(tricellvstrajectory && uncorrectedPEplots) uncorrected_PE_tricellvstrajectory(c3, det, mcdata, analysis1, analysis2, hist1, hist2);
112  else if(uncorrectedPEplots) uncorrected_PE(c1, c2, det, mcdata, analysis1, analysis2, hist1, hist2);
113
115
116  if(correctedPEplots) corrected_PE(cattenuation1, cattenuation2, det, mcdata, analysis1, analysis2, plane, cell);
117
118  if(calibrationconstantsplots) calibrationconstants(cconstants1, cconstants2, det, mcdata, analysis1, analysis2);
119  }//end of if(fd && mc && onetoone)
120 }//attenuation_calibration_comparisons()
121 /////////////////////////////////////////////////////////////////////////////////
122
123 //----------------------------------------------------------------------------------------------------------------------
124 uncorrected_PE(char* c1, char* c2, TString det, TString mcdata, TString analysis1, TString analysis2, TString hist1, TString hist2){
125  TFile* f1 = new TFile(c1, "read");
126  TFile* f2 = new TFile(c2, "read");
127  TString viewstr;
128  TString viewStr;
129  for(int view = 0; view < 2; ++view){
130  if(nddatamuoncatcher) viewstr = view ? "muy" : "mux";
131  else if(!nddatamuoncatcher) viewstr = view ? "y" : "x";
132  if(nddatamuoncatcher) viewStr = view ? "muY" : "muX";
133  else if(!nddatamuoncatcher) viewStr = view ? "Y" : "X";
134
135  f1->cd();
136  if(f1->GetListOfKeys()->Contains("calib"))
137  TH2F* h2D_1 = (TH2F*)f1->Get("calib/"+hist1+"_in"+viewStr);
138  if(f1->GetListOfKeys()->Contains("make"))
139  TH2F* h2D_1 = (TH2F*)f1->Get("make/"+hist1+"_in"+viewStr);
140  TProfile* ProfileX_1 = h2D_1->ProfileX();
141  if(det=="nd" && mcdata=="data") ProfileX_1->Rebin();
142
143  f2->cd();
144  if(f2->GetListOfKeys()->Contains("calib"))
145  TH2F* h2D_2 = (TH2F*)f2->Get("calib/"+hist2+"_in"+viewStr);
146  else if(f2->GetListOfKeys()->Contains("make"))
147  TH2F* h2D_2 = (TH2F*)f2->Get("make/"+hist2+"_in"+viewStr);
148  TProfile* ProfileX_2 = h2D_2->ProfileX();
149  if(det=="nd" && mcdata=="data") ProfileX_2->Rebin();
150
151  new TCanvas;
153  HistogramAttr1D(ProfileX_1, "W (cm)", "Uncorrected PE/cm", ProfileX_1->GetXaxis()->GetBinCenter(1), ProfileX_1->GetXaxis()->GetBinCenter(ProfileX_1->GetXaxis()->GetNbins()), 15, 80, kBlue);
154  ProfileX_1->Draw();
155  ProfileX_1->SetLineColor(kBlue);
156  ProfileX_2->SetLineColor(kRed);
157  ProfileX_2->Draw("same");
158
159  TLegend *legend = new TLegend(.12,.62,.64,.90);
160  legend->SetBorderSize(0);
161  legend->SetFillStyle(0);
165  legend->Draw();
166
168
169  //ratios
170  new TCanvas;
172  const TString hrationame = Form("hratio%s_%s_%s", det.Data(), mcdata.Data(), viewstr.Data());
173  TH1D* hratio = new TH1D(hrationame, "",ProfileX_1->GetXaxis()->GetNbins(), ProfileX_1->GetBinCenter(0), ProfileX_1->GetBinCenter(ProfileX_1->GetXaxis()->GetNbins()));
174  Ratiohist(ProfileX_1, ProfileX_2, hratio, kRed, 24, "W (cm)", "Ratio", 0.85, 1.10);
175
176  TLegend *legend = new TLegend(.12,.62,.64,.90);
177  legend->SetBorderSize(0);
178  legend->SetFillStyle(0);
179  legend->SetTextSize(0.04);
181  legend->AddEntry(hratio,analysis2+" to "+analysis1+" "+viewStr+" view");
182  legend->Draw();
184  }//end of loop over views
185 } // end of uncorrected_PE()
186 //-------------------------------------------------------------------------------------------------------------------------------
187 void uncorrected_PE_tricellvstrajectory(char* c1, TString det, TString mcdata, TString analysis1, TString analysis2, TString hist1, TString hist2){
188  uncorrected_PE(c1, c1, det, mcdata, analysis1, analysis2, hist1, hist2);
189 }//end of uncorrected_PE_tricellvstrajectory()
190 //-------------------------------------------------------------------------------------------------------------------------------
191 void thresholdshadowing(char* c1, char* c2, TString det, TString mcdata, TString analysis1, TString analysis2, TString var){
192  TFile* f1 = new TFile(c1, "read");
193  TFile* f2 = new TFile(c2, "read");
194  for(int view = 0; view < 2; ++view){
195  const TString viewstr = view ? "y" : "x";
196  const TString viewStr = view ? "Y" : "X";
197
198  f1->cd();
199  TH1D* correction_1 = (TH1D*)f1->Get("a"+viewstr+"view");
200  correction_1->SetLineColor(kBlack);
201
202  f2->cd();
203  TH1D* correction_2 = (TH1D*)f2->Get("a"+viewstr+"view");
204  correction_2->SetLineColor(kRed);
205
206  new TCanvas;
208  if(correction_1->FindObject("stats")){
209  TPaveStats* pave1 = (TPaveStats*)correction_1->FindObject("stats");
210  pave1->Delete();
211  }
212  if(correction_2->FindObject("stats")){
213  TPaveStats* pave2 = (TPaveStats*)correction_2->FindObject("stats");
214  pave2->Delete();
215  }
216  correction_1->Draw();
217  correction_2->Draw("same");
218  if(var=="W") TString labelx = "W (cm)";
219  HistogramAttr1D(correction_1, labelx.Data(), "correction factor", correction_1->GetXaxis()->GetBinCenter(1), correction_1->GetXaxis()->GetBinCenter(correction_1->GetXaxis()->GetNbins()), 0.5, 2.0, kBlack);
220
221  TLegend *legend = new TLegend(.12,.62,.64,.90);
222  legend->SetBorderSize(0);
223  legend->SetFillStyle(0);
227  legend->Draw();
229
230  //ratios
231  new TCanvas;
233  const TString hrationame = Form("hratio%s_%s_%s_%s", det.Data(), mcdata.Data(), viewstr.Data(), var.Data());
234  TH1D* hratio = new TH1D(hrationame, "",correction_1->GetXaxis()->GetNbins(), correction_1->GetBinCenter(0), correction_1->GetBinCenter(correction_1->GetXaxis()->GetNbins()));
235  Ratiohist(correction_1, correction_2, hratio, kRed, 24, labelx.Data(), "Correction factor ratio", 1.8, 1.2);
236
237  TLegend *legend = new TLegend(.12,.62,.64,.90);
238  legend->SetBorderSize(0);
239  legend->SetFillStyle(0);
240  legend->SetTextSize(0.04);
242  legend->AddEntry(hratio,analysis2+" to "+analysis1+" "+viewStr+" view");
243  legend->Draw();
245  }//end of loop over views
247 //-------------------------------------------------------------------------------------------------------------------------------
248 void corrected_PE(char* c1, char* c2, TString det, TString mcdata, TString analysis1, TString analysis2, TString plane, TString cell){
249
250  if(nddatamuoncatcher){
251  if(plane.Atoi() > 3 || plane.Atoi() < 2){
252  std::cout << "It is a ND muon catcher data file and has only 002, 003 planes. Plane " << plane << " is out of range. " << std::endl;
253  abort();
254  }
255  if(cell.Atoi() > 95 || cell.Atoi() < 0){
256  std::cout << "It is a Data file. Cell " << cell << " is out of range. " << std::endl;
257  abort();
258  }
259  }
260  if(det.Data()=="nd"){
261  if(mcdata.Data()=="mc"){
262  if(plane.Atoi() > 3 || plane.Atoi() < 0){
263  std::cout << "It is a Monte Carlo file. Plane " << plane << " out of range. " << std::endl;
264  abort();
265  }
266  if(cell.Atoi() > 95 || cell.Atoi() < 0){
267  std::cout << "It is a Monte Carlo file. Cell " << cell << " out of range. " << std::endl;
268  abort();
269  }
270  }//end of MC condition
271  if(mcdata.Data()=="data"){
272  if(plane.Atoi() > 191 || plane.Atoi() < 0){
273  std::cout << "It is a ND data file and has only 191 main body planes. Plane " << plane << " out of range. " << std::endl;
274  abort();
275  }
276  if(cell.Atoi() > 95 || cell.Atoi() < 0){
277  std::cout << "It is a Data file. Cell " << cell << " out of range. " << std::endl;
278  abort();
279  }
280  }//end of data condition
281  }//end of the ND condition
282
283  if(det.Data()=="fd"){
284  if(mcdata.Data()=="mc"){
285  if(plane.Atoi() > 1 || plane.Atoi() < 0){
286  std::cout << "It is a Monte Carlo file. Plane " << plane << " out of range. " << std::endl;
287  abort();
288  }
289  if(cell.Atoi() > 383 || cell.Atoi() < 0){
290  std::cout << "It is a Monte Carlo file. Cell " << cell << " out of range. " << std::endl;
291  abort();
292  }
293  }//end of MC condition
294  else if(mcdata.Data()=="data"){
295  if(plane.Atoi() > 895 || plane.Atoi() < 0){
296  std::cout << "It is a Data file. Plane " << plane << " out of range. " << std::endl;
297  abort();
298  }
299  if(cell.Atoi() > 383 || cell.Atoi() < 0){
300  std::cout << "It is a Data file. Cell " << cell << " out of range. " << std::endl;
301  abort();
302  }
303  }//end of data condition
304  }//end of FD condition
305
306  TFile* f1 = new TFile(c1, "read");
307  TFile* f2 = new TFile(c2, "read");
308
309  TF1* fit1;
310  TF1* fit2;
311
312  f1->cd();
313  TH1F* h1 = (TH1F*)f1->Get("plane_"+plane+"/Atten_Fit_"+plane+"_"+cell);
314  Color_t color1 = kBlack;
315  if(mcdata=="mc"){
316  if(plane.Atoi()==0) fit1 = (TF1*)h1->FindObject("fitX");
317  else if(plane.Atoi()==1) fit1 = (TF1*)h1->FindObject("fitY");
318  else if(plane.Atoi()==2) fit1 = (TF1*)h1->FindObject("fitMuX");
319  else if(plane.Atoi()==3) fit1 = (TF1*)h1->FindObject("fitMuY");
320  if(fit1) fit1->SetLineColor(color1);
321  }//end of MC condition
322  else if(mcdata=="data"){
323  if(nddatamuoncatcher){
324  if(plane.Atoi()==2) fit1 = (TF1*)h1->FindObject("fitMuX");
325  else if(plane.Atoi()==3) fit1 = (TF1*)h1->FindObject("fitMuY");
326  if(fit1) fit1->SetLineColor(color1);
327  }//end of muon catcher condition
328  else if(plane.Atoi()%2==0) fit1 = (TF1*)h1->FindObject("fitY");
329  else if(plane.Atoi()%2!=0) fit1 = (TF1*)h1->FindObject("fitX");
330  if(fit1) fit1->SetLineColor(color1);
331  }//end of data condition
332
333  f2->cd();
334  TH1F* h2 = (TH1F*)f2->Get("plane_"+plane+"/Atten_Fit_"+plane+"_"+cell);
335  Color_t color2 = kRed;
336  if(mcdata=="mc"){
337  if(plane.Atoi()==0) fit2 = (TF1*)h2->FindObject("fitX");
338  else if(plane.Atoi()==1) fit2 = (TF1*)h2->FindObject("fitY");
339  else if(plane.Atoi()==2) fit2 = (TF1*)h2->FindObject("fitMuX");
340  else if(plane.Atoi()==3) fit2 = (TF1*)h2->FindObject("fitMuY");
341  if(fit2) fit2->SetLineColor(color2);
342  }//end of MC condition
343  else if(mcdata=="data"){
344  if(nddatamuoncatcher){
345  if(plane.Atoi()==2) fit2 = (TF1*)h2->FindObject("fitMuX");
346  else if(plane.Atoi()==3) fit2 = (TF1*)h2->FindObject("fitMuY");
347  if(fit2) fit2->SetLineColor(color2);
348  }//end of muon catcher condition
349  else if(plane.Atoi()%2==0) fit2 = (TF1*)h2->FindObject("fitY");
350  else if(plane.Atoi()%2!=0) fit2 = (TF1*)h2->FindObject("fitX");
351  if(fit2) fit2->SetLineColor(color2);
352  }//end of data condition
353
354  new TCanvas;
356  h1->Draw();
357  HistogramAttr1D(h1, "W (cm)", "Mean PE/cm", h1->GetXaxis()->GetBinCenter(1), h1->GetXaxis()->GetBinCenter(h1->GetXaxis()->GetNbins()), 0., 150., color1);
358  h1->SetLineColor(color1);
359  h1->SetMarkerColor(color1);
360  h1->SetMarkerStyle(24);
361  h2->Draw("same");
362  h2->SetLineColor(color2);
363  h2->SetMarkerColor(color2);
364  h2->SetMarkerStyle(24);
365
366  TLegend *legend = new TLegend(.12,.62,.64,.90);
367  legend->SetBorderSize(0);
368  legend->SetFillStyle(0);
369  legend->SetTextSize(0.04);
373  legend->Draw();
375
376  //ratio of data points
377  new TCanvas;
379  const TString hrationame = Form("hratio%s_%s_%s_%s", det.Data(), mcdata.Data(), plane.Data(), cell.Data());
380  TH1D* hratio = new TH1D(hrationame, "",h1->GetXaxis()->GetNbins(), h1->GetBinCenter(0), h1->GetBinCenter(h1->GetXaxis()->GetNbins()));
381  Ratiohist(h1, h2, hratio, kRed, 24, "W (cm)", "Ratio of Mean PE/cm data points", 0.3, 1.0);
382
383  TLegend *legend = new TLegend(.12,.62,.64,.90);
384  legend->SetBorderSize(0);
385  legend->SetFillStyle(0);
386  legend->SetTextSize(0.04);
387  legend->AddEntry((TObject*)0, det+" "+mcdata+" : Plane "+plane+" Cell "+cell," ");
389  legend->Draw();
391
392  //ratio of fitting curve
393  new TCanvas;
395  const TString hratiofitname = Form("hratiofit%s_%s_%s_%s", det.Data(), mcdata.Data(), plane.Data(), cell.Data());
396  if(det=="nd"){
397  float x_limit = 250.;
398  int x_bins = 2*x_limit;
399  }
400  else if(det=="fd"){
401  float x_limit = 800.;
402  int x_bins = 2*x_limit;
403  }
404  TH1D* hfit = new TH1D("hfit", "", x_bins, -x_limit, +x_limit);
405  HistogramAttr1D(hfit, "W (cm)", "Ratio of Mean PE/cm fitting curves", -x_limit, +x_limit, 0.55, 0.68, kWhite);
406  hfit->Draw();
407  TH1D* hratiofit = new TH1D(hratiofitname, "", x_bins, -x_limit, +x_limit);
408  Ratiofit(fit1, fit2, hratiofit, plane.Data(), cell.Data(), color2, 24, "W (cm)", "Ratio of Mean PE/cm fitting curves", 0.55, 0.68);
409  TLegend *legend = new TLegend(.12,.62,.64,.90);
410  legend->SetBorderSize(0);
411  legend->SetFillStyle(0);
412  legend->SetTextSize(0.04);
413  legend->AddEntry((TObject*)0, det+" "+mcdata+" : Plane "+plane+" Cell "+cell," ");
415  legend->Draw();
417 }//end of corrected_PE()
418 //---------------------------------------------------------------------------------------------------
419 void calibrationconstants(char* c1, char* c2, TString det, TString mcdata, TString analysis1, TString analysis2){
420  ifstream in_1;
421  ifstream in_2;
422
423  ofstream out;
424
425  FILE* fout = fopen("uncalibrated_channels_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+".csv","w");
426  fprintf(fout, "%s%50s%s\n", analysis1.Data(), " ", analysis2.Data());
427  fprintf(fout, "channel,coeff_exp,atten_length,background,chisq, %9s, channel,coeff_exp,atten_length,background,chisq\n", " ");
428
429  in_1.open(c1);
430  in_2.open(c2);
431
432
433  //x=coeff_exp, y=atten_length, z=background, a=edgx_low, b=edgx_high, c=coeff_low
434  Float_t x_1, y_1, z_1, a_1, b_1, c_1, d_1, e_1, g_1;
435  Float_t x_2, y_2, z_2, a_2, b_2, c_2, d_2, e_2, g_2;
436  Float_t chisq = 0.2;
437
438  Int_t channels;
439  Int_t totalchannels1 = 0;
440  Int_t totalchannels2 = 0;
441  Int_t calibratedchannels1 = 0;
442  Int_t calibratedchannels2 = 0;
443  Int_t f_1;
444  Int_t f_2;
445  Int_t bins1;
446
447  Double_t limit1D = 3.;
448  Double_t limit2D = 1000.;
449
450  if(det=="nd" && mcdata=="mc" ) { channels = 3063; }
451  if(det=="nd" && mcdata=="data" ) { channels = 1003063; }
452  if(det=="fd" && mcdata=="mc" ) { channels = 1383; }
453  if(det=="fd" && mcdata=="data" ) { channels = 895383; }
454
455
456  TH1F* hcoeffexp = new TH1F("hcoeffexp", "coeffexp difference distribution", 120, -limit1D, +limit1D );
457  TH1F* hattenlength = new TH1F("hattenlength", "attenlength difference distribution", 600, -limit1D, +limit1D );
458  TH1F* hbackground = new TH1F("hbackground", "background difference distribution", 120, -limit1D, +limit1D );
459  TH1F* hchisq = new TH1F("hchisq", "chisq difference distribution", 600, -limit1D, +limit1D );
460
461  TH2F* h2coeffexp = new TH2F("h2coeffexp", " " , 2*(Int_t)limit2D, -limit2D, limit2D, 2*(Int_t)limit2D, -limit2D, limit2D );
462  TH2F* h2attenlength = new TH2F("h2attenlength"," " , 2*(Int_t)limit2D, -limit2D, limit2D, 2*(Int_t)limit2D, -limit2D, limit2D );
463  TH2F* h2background = new TH2F("h2background", " " , 2*(Int_t)limit2D, -limit2D, limit2D, 2*(Int_t)limit2D, -limit2D, limit2D );
464  TH2F* h2chisq = new TH2F("h2chisq"," ", 200,0.,2., 200,0.,2. );
465
466  while(1){
467  in_1 >> x_1 >> y_1 >> z_1 >> a_1 >> b_1 >> c_1 >> d_1 >> e_1 >> f_1 >> g_1;
468  in_2 >> x_2 >> y_2 >> z_2 >> a_2 >> b_2 >> c_2 >> d_2 >> e_2 >> f_2 >> g_2;
469  if(!in_1.good()) break;
470  if(!in_2.good()) break;
471
472  if(f_1<=channels && e_1 < chisq && e_2 < chisq){
473  hcoeffexp->Fill((x_1 - x_2)/x_1);
474  hattenlength->Fill((y_1 - y_2)/y_1);
475  hbackground->Fill((z_1 - z_2)/z_1);
476
477  h2coeffexp->Fill(x_2, x_1);
478  h2attenlength->Fill(y_2, y_1);
479  h2background->Fill(z_2, z_1);
480  }//end of loop over channels and chisq conditions.
481  if(f_1<=channels){
482  totalchannels1++;
483  totalchannels2++;
484  if(e_1<chisq) calibratedchannels1++;
485  if(e_2<chisq) calibratedchannels2++;
486
487  hchisq->Fill((e_1 - e_2)/e_1);
488  h2chisq->Fill(e_2, e_1);
489
490  if(e_1>=chisq || e_2>=chisq){
491  fprintf(fout, "%i,%g,%g,%g,%g,%20s,%i,%g,%g,%g,%g\n",f_1,x_1,y_1,z_1,e_1, " ", f_2,x_2,y_2,z_2,e_2);
492  out << f_1 << setw(22) << x_1 << setw(22) << y_1 << setw(22) << z_1 << setw(22) << e_1 << f_2 << setw(22) << x_2 << setw(22) << y_2 << setw(22) << z_2 << setw(22) << e_2 << endl;
493  }
494  }//end of only channels condition for chisq plots. For chisq plots we have to consider all channels.
495  }//end of while condition.
496
497  cout << " " << endl;
498  cout << "Total channels in "+analysis1+" = " << totalchannels1 << " and calibrated channels = " << calibratedchannels1 << " = " << ((float)calibratedchannels1/(float)totalchannels1)*100. << " %." << endl;
499  cout << "Total channels in "+analysis2+" = " << totalchannels2 << " and calibrated channels = " << calibratedchannels2 << " = " << ((float)calibratedchannels2/(float)totalchannels2)*100. << " %." << endl;
500  cout << " " << endl;
501
502  //coeffexp plots
503  new TCanvas;
505  HistogramAttr1D(hcoeffexp, "#Deltacoeffexp/coeffexp", "# of channels", 0.5*hcoeffexp->GetXaxis()->GetBinCenter(1), 0.5*hcoeffexp->GetXaxis()->GetBinCenter(hcoeffexp->GetXaxis()->GetNbins()), 0., 7e5, kBlack);
506  hcoeffexp->Draw();
507  TLegend *legend = new TLegend(0.00, 0.60, 0.47, 0.93);
508  legend->SetBorderSize(0);
509  legend->SetFillStyle(0);
510  legend->SetTextSize(0.05);
513  legend->Draw();
515
516  new TCanvas;
519  HistogramAttr2D(h2coeffexp, "coeffexp_{"+analysis2+"}", "coeffexp_{"+analysis1+"}", " ", 0, 100, 0, 100);
520  h2coeffexp->Draw("colz");
522
523  //attenlength plots
524  new TCanvas;
526  HistogramAttr1D(hattenlength, "#Deltaattenlength/attenlength", "# of channels", 0.5*hattenlength->GetXaxis()->GetBinCenter(1), 0.5*hattenlength->GetXaxis()->GetBinCenter(hattenlength->GetXaxis()->GetNbins()), 0., 166500, kBlack);
527  hattenlength->Draw();
528  TLegend *legend = new TLegend(0.00, 0.60, 0.47, 0.93);
529  legend->SetBorderSize(0);
530  legend->SetFillStyle(0);
531  legend->SetTextSize(0.05);
534  legend->Draw();
536
537  new TCanvas;
540  HistogramAttr2D(h2attenlength, "attenlength_{"+analysis2+"} (cm)", "attenlength_{"+analysis1+"} (cm)", " ", 0, 1000, 0, 1000);
541  h2attenlength->Draw("colz");
543
544  //backgrounds plots
545  new TCanvas;
547  HistogramAttr1D(hbackground, "#Deltabackground/background", "# of channels", 0.5*hbackground->GetXaxis()->GetBinCenter(1), 0.5*hbackground->GetXaxis()->GetBinCenter(hbackground->GetXaxis()->GetNbins()), 0., 3e4, kBlack);
548  hbackground->Draw();
549  TLegend *legend = new TLegend(0.00, 0.60, 0.47, 0.93);
550  legend->SetBorderSize(0);
551  legend->SetFillStyle(0);
552  legend->SetTextSize(0.05);
555  legend->Draw();
557
558  new TCanvas;
561  HistogramAttr2D(h2background, "background_{"+analysis2+"}", "background_{"+analysis1+"}", " ", -50, 50, -50, 50);
562  h2background->Draw("colz");
564
565  //chisq plots
566  new TCanvas;
568  HistogramAttr1D(hchisq, "#Deltachisq/chisq", "# of channels", 0.5*hchisq->GetXaxis()->GetBinCenter(1), 0.5*hchisq->GetXaxis()->GetBinCenter(hchisq->GetXaxis()->GetNbins()), 0., 4200., kBlack);
569  hchisq->Draw();
570  TLegend *legend = new TLegend(0.00, 0.60, 0.47, 0.93);
571  legend->SetBorderSize(0);
572  legend->SetFillStyle(0);
573  legend->SetTextSize(0.05);
576  legend->Draw();
578
579  new TCanvas;
582  HistogramAttr2D(h2chisq, "#chi^{2}_{"+analysis2+"}", "#chi^{2}_{"+analysis1+"}", " ", 0, 0.5, 0, 0.5);
583  h2chisq->Draw("colz");
585
586  in_1.close();
587  in_2.close();
588  out.close();
589  fclose(fout);
590 }//end of calibrationconstants()
591 //---------------------------------------------------------------------------------------------------
enum BeamMode kRed
void Ratiohist(TH1 *h1, TH1 *h2, TH1 *hratio, Color_t color, Style_t mstyle, TString labelx, TString labely, Double_t lowy, Double_t highy)
TH1D * hist2
Definition: plotHisto.C:12
var_value< double > var
Definition: StanTypedefs.h:14
Float_t f2
static constexpr Double_t c_2
Definition: Munits.h:245
void corrected_PE(char *c1, char *c2, TString det, TString mcdata, TString analysis1, TString analysis2, TString plane, TString cell)
const XML_Char const XML_Char * data
Definition: expat.h:268
c2
Definition: demo5.py:33
fclose(fg1)
std::map< ToFCounter, std::vector< unsigned int > > channels
void HistogramAttr2D(TH2 *h, char *titlx_x, char *titlx_y, char *titlx_z, Double_t binlowx, Double_t binhighx, Double_t binlowy, Double_t binhighy)
uncorrected_PE(char *c1, char *c2, TString det, TString mcdata, TString analysis1, TString analysis2, TString hist1, TString hist2)
TH1D * hist1
Definition: plotHisto.C:9
Float_t f1
OStream cout
Definition: OStream.cxx:6
TH1F * h2
Definition: plot.C:45
void Ratiofit(TF1 *f1, TF1 *f2, TH1 *hratio, TString plane, TString cell, Color_t color, Style_t mstyle, TString labelx, TString labely, Double_t lowy, Double_t highy)
TH1F * h1
void thresholdshadowing(char *c1, char *c2, TString det, TString mcdata, TString analysis1, TString analysis2, TString var)
void uncorrected_PE_tricellvstrajectory(char *c1, TString det, TString mcdata, TString analysis1, TString analysis2, TString hist1, TString hist2)
c1
Definition: demo5.py:24
HistogramAttr1D(TH1 *h, char *title_x, char *title_y, Double_t binlowx, Double_t binhighx, Double_t binlowy, Double_t binhighy, Color_t lcolor)
void calibrationconstants(char *c1, char *c2, TString det, TString mcdata, TString analysis1, TString analysis2)
enum BeamMode kBlue