nue_fnex_vs_caf_energycomparison.C
Go to the documentation of this file.
1 //Author: Prabhjot Singh (prabhjot@fnal.gov).
2 //Dated: Sep 29, 00:00
3 //Purpose: Compare Hadronic, Slice and Neutrino energies between CAF and FNEX to figure out the energy differences for the selected events for
4 //1. Nue Selected FD data
5 
6 #include <Header.h>
7 #include <HistogramAttr.h>
8 
9 void DrawHistograms(TH1D*, TH1D*, TString, TCanvas*, TString);
10 
11 //----------------------------------------------------------//
13  gStyle->SetOptStat(" ");
14 
15  int i = 0;
16 
17  //first file is caf file as base file
18  //second file is fnex file as secondary file
19  char* c_caf = "/nova/app/users/prabhjot/art/mytestrel_caf/run_jobs/cvn_numi_sorted.txt"; //sorted acording to the run number
20  char* c_fnex = "/nova/ana/users/prabhjot/FNEX/nue_analysis/FD_data_energies_info_sorted.txt"; //sorted acording to the run number
21 
22  std::ifstream in_caf;
23  std::ifstream in_fnex;
24 
25  in_caf.open(c_caf);
26  in_fnex.open(c_fnex);
27 
28  //checks to see if the CAF or/and FNEX files exist or not.
29  if(!in_caf){
30  std::cout << "CAF File does not exists. ABORTING..." << std::endl;
31  abort();
32  }
33  if(!in_fnex){
34  std::cout << "FNEX File does not exists. ABORTING..." << std::endl;
35  abort();
36  }
37 
38  //Variable to grab run, subrun, event, hadE, sliceE and neutrinoE
39  int kRun_caf;
40  int kSubrun_caf;
41  int kEvent_caf;
42  int kSlice_caf;
43 
44  double kSliceMeanTime_caf;
45  double kLID_caf;
46  double kLEM_caf;
47  double kCVN_caf;
48  double kNueEnergy_caf;
49  double kShowerE_caf;
50  double kHadE_caf;
51  double kSliceE_caf;
52 
53  int kDet_fnex;
54  int kFileType_fnex;
55  int kRun_fnex;
56  int kSubrun_fnex;
57  int kEvent_fnex;
58 
59  double kLID_fnex;
60  double kLEM_fnex;
61  double kCVN_fnex;
62  double kNueEnergy_fnex;
63  double kShowerE_fnex;
64  double kHadE_fnex;
65  double kSliceE_fnex;
66 
67  std::string sUnknown_fnex;
68  std::string sLID_fnex;
69  std::string sLEM_fnex;
70  std::string sCVN_fnex;
71  std::string sNueEnergy_fnex;
72  std::string sShowerE_fnex;
73  std::string sHadE_fnex;
74 
75  TH1D* h_neutrinoE_caf = new TH1D("h_neutrinoE_caf", " ", 33, 0, 33);
76  TH1D* h_neutrinoE_fnex = new TH1D("h_neutrinoE_fnex", " ", 33, 0, 33);
77 
78  TH1D* h_hadE_caf = new TH1D("h_hadE_caf", " ", 33, 0, 33);
79  TH1D* h_hadE_fnex = new TH1D("h_hadE_fnex", " ", 33, 0, 33);
80 
81  TH1D* h_showerE_caf = new TH1D("h_showerE_caf", " ", 33, 0, 33);
82  TH1D* h_showerE_fnex = new TH1D("h_showerE_fnex", " ", 33, 0, 33);
83 
84  TH1D* h_sliceE_caf = new TH1D("h_sliceE_caf", " ", 33, 0, 33);
85  TH1D* h_sliceE_fnex = new TH1D("h_sliceE_fnex", " ", 33, 0, 33);
86 
87  //read the variable
88  while(1){
89  in_caf
90  >> kRun_caf
91  >> kSubrun_caf
92  >> kEvent_caf
93  >> kSlice_caf
94  >> kSliceMeanTime_caf
95  >> kLID_caf
96  >> kLEM_caf
97  >> kCVN_caf
98  >> kNueEnergy_caf
99  >> kShowerE_caf
100  >> kHadE_caf
101  >> kSliceE_caf;
102 
103  in_fnex
104  >> kDet_fnex
105  >> kFileType_fnex
106  >> kRun_fnex
107  >> kSubrun_fnex
108  >> kEvent_fnex
109  >> sUnknown_fnex
110  >> kSliceE_fnex
111  >> sLID_fnex
112  >> kLID_fnex
113  >> sLEM_fnex
114  >> kLEM_fnex
115  >> sCVN_fnex
116  >> kCVN_fnex
117  >> sHadE_fnex
118  >> kHadE_fnex
119  >> sShowerE_fnex
120  >> kShowerE_fnex
121  >> sNueEnergy_fnex
122  >> kNueEnergy_fnex;
123 
124  if(!in_caf.good()) break;
125  if(!in_fnex.good()) break;
126 
127  if(kRun_caf != kRun_fnex){
128  std::cout << "NOT COMPARING same events. SORT by Run number" << std::endl;
129  std::cout << "CAF: Run = " << kRun_caf << " Subrun = " << kSubrun_caf << " and CVN = " << kCVN_caf << std::endl;
130  std::cout << "FNEX: Run = " << kRun_fnex << " Subrun = " << kSubrun_fnex << " and CVN = " << kCVN_fnex << std::endl;
131  std::cout << "........... ABORTING ............... " << std::endl;
132  abort();
133  }
134 
135  h_neutrinoE_caf ->SetBinContent(i+1, kNueEnergy_caf);
136  h_neutrinoE_caf ->GetXaxis()->SetBinLabel(i+1, Form("%d",kRun_caf) );
137  h_neutrinoE_fnex ->SetBinContent(i+1, kNueEnergy_fnex);
138 
139  h_hadE_caf ->SetBinContent(i+1, kHadE_caf);
140  h_hadE_caf ->GetXaxis()->SetBinLabel(i+1, Form("%d",kRun_caf) );
141  h_hadE_fnex ->SetBinContent(i+1, kHadE_fnex);
142 
143  h_showerE_caf ->SetBinContent(i+1, kShowerE_caf);
144  h_showerE_caf ->GetXaxis()->SetBinLabel(i+1, Form("%d",kRun_caf) );
145  h_showerE_fnex ->SetBinContent(i+1, kShowerE_fnex);
146 
147  h_sliceE_caf ->SetBinContent(i+1, kSliceE_caf);
148  h_sliceE_caf ->GetXaxis()->SetBinLabel(i+1, Form("%d",kRun_caf) );
149  h_sliceE_fnex ->SetBinContent(i+1, kSliceE_fnex);
150 
151  /* std::cout
152  << "i = "
153  << i
154  << " "
155  << "Run = "
156  << kRun_caf
157  << " "
158  << "neutrinoE: caf = "
159  << kNueEnergy_caf
160  << " "
161  << "neutrinoE: fnex = "
162  << kNueEnergy_fnex
163  << " and fraction (caf - fnex)/caf = "
164  << (kNueEnergy_caf - kNueEnergy_fnex)/ kNueEnergy_caf
165  << std::endl;
166  */
167  i++;
168  }//end of while loop i.e. reading of variables
169 
170 
171  //Draw histograms
172  TCanvas* c_neutrinoE = new TCanvas("c_neutrinoE", " ", 800, 800);
173  DrawHistograms(h_neutrinoE_caf, h_neutrinoE_fnex, "Neutrino energy (GeV)", c_neutrinoE, "FD_data_neutrinoE");
174 
175  TCanvas* c_hadE = new TCanvas("c_hadE", " ", 800, 800);
176  DrawHistograms(h_hadE_caf, h_hadE_fnex, "Hadronic energy (GeV)", c_hadE, "FD_data_hadE");
177 
178  TCanvas* c_showerE = new TCanvas("c_showerE", " ", 800, 800);
179  DrawHistograms(h_showerE_caf, h_showerE_fnex, "Shower energy (GeV)", c_showerE, "FD_data_showerE");
180 
181  TCanvas* c_sliceE = new TCanvas("c_sliceE", " ", 800, 800);
182  DrawHistograms(h_sliceE_caf, h_sliceE_fnex, "Slice energy (GeV)", c_sliceE, "FD_data_sliceE" );
183 }//end of main function nue_fnex_vs_caf_energycomparison
184 //----------------------------------------------------------//
185 //
186 //Draw histograms and their ratios
187 void DrawHistograms(TH1D* h1, TH1D* h2, TString Ylabel, TCanvas* c1, TString plotname){
188  //h1 is base histogram and h2 is secondary
189 
190  double ratio = 0.;
191  double ratio_error = 0.;
192 
193  TCanvas* c1;
194 
195  TPad* pad1;
196  TPad* pad2;
197 
198  TLegend* legend;
199  TLegend* legend_ratio;
200 
201  TLine* one;
202  TLine* one_main;
203 
204  TH1F* hist_ratio;
205 
206  pad1 = new TPad("pad1","",0.05,0.20,0.95,1);
207  pad2 = new TPad("pad2","",0.05,0.05,0.95,0.28);
208  pad1->Draw();
209  pad2->Draw();
210  pad1->cd();
211  gPad->SetLeftMargin(0.15);
212 
213  h1->Draw("P");
214  HistogramAttr1D(h1, "Run number", Ylabel);
215  h1->SetMarkerStyle(20);
216  h1->SetMarkerSize(1.50);
217  h1->SetMarkerColor(kRed);
218 
219  h2->Draw("P same");
220  h2->SetMarkerStyle(20);
221  h2->SetMarkerColor(kBlue);
222  h2->SetMarkerSize(1.00);
223  if(Ylabel=="Neutrino energy (GeV)"){
224  one_main = new TLine(0, 2, 33, 2);
225  one_main->SetLineWidth(2);
226  one_main->SetLineStyle(7);
227  one_main->Draw();
228  }
229 
230  //position of legends at right place for different cases
231  if(Ylabel=="Neutrino energy (GeV)") legend = new TLegend(0.31, 0.14, 0.69, 0.42);
232  else if(Ylabel=="Hadronic energy (GeV)") legend = new TLegend(0.63, 0.61, 0.95, 0.88);
233  else if(Ylabel=="Shower energy (GeV)") legend = new TLegend(0.66, 0.12, 0.94, 0.36);
234  else if(Ylabel=="Slice energy (GeV)") legend = new TLegend(0.35, 0.13, 0.64, 0.39);
235  else legend = new TLegend(0.63, 0.61, 0.95, 0.88);
236  legend->SetBorderSize(0);
237  legend->SetFillStyle(0);
238  legend->AddEntry((TObject*)0, "FD data" ," ");
239  legend->AddEntry(h1, "CAF", "p");
240  legend->AddEntry(h2, "FNEX", "p");
241  legend->Draw();
242 
243  hist_ratio = new TH1F("hist_ratio_"+plotname, " ", h1->GetNbinsX(), h1->GetXaxis()->GetXmin(), h1->GetXaxis()->GetXmax());
244  for(int i = 0; i < h1->GetNbinsX() + 1; i++){
245 
246  if(h1->GetBinContent(i)==0) ratio = 0.;
247 
248  else ratio = h2->GetBinContent(i)/h1->GetBinContent(i);
249 
250  if(h1->GetBinContent(i)==0) ratio_error = 0.;
251  //else ratio_error = error_division(h2->GetBinContent(i), h1->GetBinContent(i));
252  hist_ratio->SetBinContent(i, ratio);
253  hist_ratio->SetBinError(i, ratio_error);
254 
255  hist_ratio->GetXaxis()->SetBinLabel(i+1, h1->GetXaxis()->GetBinLabel(i+1));
256  }//end of bin by bin ratio
257 
258  pad2->cd();
259  gPad->SetBottomMargin(0.50);
260  gPad->SetLeftMargin(0.15);
261  hist_ratio->Draw("P");
262  hist_ratio->LabelsOption("v", "X");
263  hist_ratio->SetMarkerColor(kBlack);
264  hist_ratio->SetMarkerStyle(20);
265  hist_ratio->GetYaxis()->SetTitle("Ratio");
266  hist_ratio->GetYaxis()->CenterTitle();
267  hist_ratio->GetYaxis()->SetTitleSize(0.20);
268  hist_ratio->GetYaxis()->SetTitleOffset(.20);
269  hist_ratio->GetYaxis()->SetLabelSize(0.10);
270  hist_ratio->GetYaxis()->SetDecimals();
271  hist_ratio->GetYaxis()->SetRangeUser(0.95, 1.05);
272  hist_ratio->GetXaxis()->SetTitle("Run number");
273  hist_ratio->GetXaxis()->CenterTitle();
274  hist_ratio->GetXaxis()->SetTitleSize(0.15);
275  hist_ratio->GetXaxis()->SetTitleOffset(1.4);
276  hist_ratio->GetXaxis()->SetLabelSize(0.17);
277  hist_ratio->GetXaxis()->SetLabelOffset(0.015);
278 
279  legend->AddEntry(hist_ratio, "FNEX/CAF", "p");
280  legend_ratio = new TLegend(0.65, 0.45, 0.93, 0.97);
281  legend_ratio->SetBorderSize(0);
282  legend_ratio->SetFillStyle(0);
283  legend_ratio->Draw();
284 
285  one = new TLine(0, 1, 33, 1);
286  one->SetLineWidth(2);
287  one->SetLineStyle(7);
288  one->Draw();
289 
290  c1->cd();
291  c1->SaveAs("Images/fnexvscaf/"+plotname+".pdf");
292 
293 }//end of DrawHistograms()
294 //-----------------------------------------------------------------------------//
enum BeamMode kRed
void nue_fnex_vs_caf_energycomparison()
TH1 * ratio(TH1 *h1, TH1 *h2)
void DrawHistograms(TH1D *, TH1D *, TString, TCanvas *, TString)
OStream cout
Definition: OStream.cxx:6
TH1F * h2
Definition: plot.C:45
TH1F * h1
c1
Definition: demo5.py:24
auto one()
Definition: PMNS.cxx:49
TPad * pad2
Definition: analysis.C:13
HistogramAttr1D(TH1 *h, char *title_x, char *title_y, Double_t binlowx, Double_t binhighx, Double_t binlowy, Double_t binhighy, Color_t lcolor)
enum BeamMode kBlue
TPad * pad1
Definition: analysis.C:13
enum BeamMode string