getPlotsIncludingAllCC.C
Go to the documentation of this file.
1 //Run using: root -b -q -l getPlots.C
2 
3 //This script should be run first in one's journey to find new fit constants.
4 //This script makes the base histograms that you will need to fit.
5 //Fitting is fast and takes tinkering at times, so it is separate from this stage, which can take a long time
6 
7 //This needs to be run to determine the muon reco E.
8 //THEN run it again to get the hadronic colz plots.
9 //DON'T do these simultaneously because it won't be consistent!
10 //You need to re-create files after you committing the final FD muon fit values!
11 
12 #include <iostream.h>
13 #include <fstream>
14 #include <iomanip.h>
15 #include <sstream>
16 #include <string>
17 #include <cstring>
18 
19 Double_t MuonEFromTrackLength(double trkLen)
20 {
21 
22  //This function returns muon energy, in GeV, given the muon track length in cm
23 
24  double muonE = 0.0;
25 
26  if (trkLen <= 0.0) return muonE;
27 
28  double stitch1 = 317.721; // cm
29  double stitch2 = 539.019; // cm
30  double stitch3 = 1012.83; // cm
31  double offset = 0.150468; // GeV
32  double slope1 = 0.00190560; // GeV/cm
33  double slope2 = 0.00198726; // GeV/cm
34  double slope3 = 0.00203830; // GeV/cm
35  double slope4 = 0.00212812; // GeV/cm
36 
37  if (trkLen < stitch1){
38  muonE = slope1*trkLen + offset;
39  }
40  else if (trkLen < stitch2){
41  muonE = slope2*trkLen + ((slope1-slope2)*stitch1 + offset);
42  }
43  else if (trkLen <stitch3){
44  muonE = slope3*trkLen + ((slope2-slope3)*stitch2 +(slope1-slope2)*stitch1 + offset);
45  }
46  else{
47  muonE = slope4*trkLen + ((slope3-slope4)*stitch3 +(slope2-slope3)*stitch2 +(slope1-slope2)*stitch1 + offset);
48  }
49 
50  return muonE;
51 
52 }//End of MuonEFromTrackLength
53 
54 Double_t ECCFromCalE(double E)
55 {
56  // This function returns the energy, in GeV, of the hadronic shower, given the calorimetric energy
57  double HadE = 0.0;
58 
59  if(E <= 0.0) return HadE;
60 
61  double stitch1 = 3.75525e-02; // GeV
62  double stitch2 = 2.05000e-01; // GeV
63  double stitch3 = 1.22411e+00; // GeV
64  double offset = 2.71346e-01; // GeV
65  double slope1 = 1.26967e+00; // unitless
66  double slope2 = 9.31211e-01; // unitless
67  double slope3 = 1.56865e+00; // unitless
68  double slope4 = 2.30327e+00; // unitless
69 
70  if (E < stitch1){
71  HadE = slope1*E + offset;
72  }
73  else if (E < stitch2){
74  HadE = slope2*E + ((slope1-slope2)*stitch1 + offset);
75  }
76  else if (E < stitch3){
77  HadE = slope3*E + ((slope2-slope3)*stitch2 +(slope1-slope2)*stitch1 + offset);
78  }
79  else{
80  HadE = slope4*E + ((slope3-slope4)*stitch3 +(slope2-slope3)*stitch2 +(slope1-slope2)*stitch1 + offset);
81  }
82 
83  return HadE;
84 } // End of ECCFromCalE
85 
86 
87 void getPlots()
88 {
89  std::cout << "Hiya!" << std::endl;
90 
91  // Energy estimators
92  float CalCCEnergy ;
93  float QEEnergy ;
94  float NonQEEnergy ;
95  float AngleQEEnergy ;
96  float AngleQEError ;
97  float HadCalEnergy ;
98  float HadTrkEnergy ;
99  float ndTrkLenAct ;
100  float ndTrkLenCat ;
101  float ndTrkCalAct ;
102  float ndTrkCalTran ;
103  float ndTrkCalCat ;
104  float ndHadCalAct ;
105  float ndHadCalTran ;
106  float ndHadCalCat ;
107 
108  bool fND; ///< is detector ND?
109  bool fFD; ///< is detector FD?
110 
111  bool isCC;
112  bool TrueQEInt; ///< is the interaction QE?
113  bool trackMatchesMuon;
114  bool contained;
115 
116  int Run;
117  int SubRun;
118  int Evt;
119  int SubEvt;
120  int goodMuon;
121  int bestTrack;
122  int PDGCode;
123  int nKalman;
124  int firstplane;
125  int lastplane;
126  int nplanestofront;
127  int nplanestoback;
128  unsigned int mincellsfromedge;
129 
130  double NuPurity;
131  double NuEfficiency;
132  double trackLength;
133  double bestPID;
134  double trackMuonEffHits;
135  double trackMuonPurHits;
136  double start_x_reco;
137  double start_y_reco;
138  double start_z_reco;
139  double end_x_reco;
140  double end_y_reco;
141  double end_z_reco;
142  double TrueNuE;
143  double TrueMuonE;
144 
145  TChain* tree = new TChain("numuana/NumuETree","");
146  tree->Add("../*.root");
147 
148  std::cout << "Reading " << tree->GetEntries() << " entries..." << std::endl;
149 
150  tree->SetBranchAddress("nKalman",&nKalman);
151  tree->SetBranchAddress("bestTrack",&bestTrack);
152  tree->SetBranchAddress("CalCCEnergy",&CalCCEnergy);
153  tree->SetBranchAddress("NuEfficiency",&NuEfficiency);
154  tree->SetBranchAddress("NuPurity",&NuPurity);
155  tree->SetBranchAddress("isCC",&isCC);
156  tree->SetBranchAddress("PDGCode",&PDGCode);
157  tree->SetBranchAddress("trackMatchesMuon",&trackMatchesMuon);
158  tree->SetBranchAddress("contained",&contained);
159  tree->SetBranchAddress("TrueQEInt",&TrueQEInt);
160  tree->SetBranchAddress("goodMuon",&goodMuon);
161  tree->SetBranchAddress("bestTrack",&bestTrack);
162  tree->SetBranchAddress("trackMuonEffHits",&trackMuonEffHits);
163  tree->SetBranchAddress("trackMuonPurHits",&trackMuonPurHits);
164  tree->SetBranchAddress("TrueNuE",&TrueNuE);
165  tree->SetBranchAddress("TrueMuonE",&TrueMuonE);
166  tree->SetBranchAddress("HadCalEnergy",&HadCalEnergy);
167  tree->SetBranchAddress("HadTrkEnergy",&HadTrkEnergy);
168  tree->SetBranchAddress("trackLength",&trackLength);
169  tree->SetBranchAddress("bestPID",&bestPID);
170  tree->SetBranchAddress("Run",&Run);
171  tree->SetBranchAddress("SubRun",&SubRun);
172  tree->SetBranchAddress("Evt",&Evt);
173  tree->SetBranchAddress("SubEvt",&SubEvt);
174  tree->SetBranchAddress("QEEnergy",&QEEnergy);
175  tree->SetBranchAddress("NonQEEnergy",&NonQEEnergy);
176 
177  //Making array for variable binning
178  double hadronQEBins[78];
179  double hadronQEAxis = 0.0;
180  for(int i = 0; i < 78; ++i){
181  hadronQEBins[i] = hadronQEAxis;
182  if (hadronQEAxis < 0.5){hadronQEAxis = hadronQEAxis + 0.01;}
183  else if (hadronQEAxis < 1.0){hadronQEAxis = hadronQEAxis + 0.025;}
184  else if (hadronQEAxis < 1.5){hadronQEAxis = hadronQEAxis + 0.1;}
185  else {hadronQEAxis = hadronQEAxis + 0.25;}
186  }
187 
188  double hadronNonQEBins[131];
189  double hadronNonQEAxis = 0.0;
190  for(int i = 0; i < 131; ++i){
191  hadronNonQEBins[i] = hadronNonQEAxis;
192  if (hadronNonQEAxis < 1.0){hadronNonQEAxis = hadronNonQEAxis + 0.01;}
193  else if (hadronNonQEAxis < 1.5){hadronNonQEAxis = hadronNonQEAxis + 0.025;}
194  else {hadronNonQEAxis = hadronNonQEAxis + 0.05;}
195  }
196 
197  double hadronCCBins[131];
198  double hadronCCAxis = 0.0;
199  for(int i = 0; i < 131; ++i){
200  hadronCCBins[i] = hadronCCAxis;
201  if (hadronCCAxis < 1.0){hadronCCAxis = hadronCCAxis + 0.01;}
202  else if (hadronCCAxis < 1.5){hadronCCAxis = hadronCCAxis + 0.025;}
203  else {hadronCCAxis = hadronCCAxis + 0.05;}
204  }
205 
206  TH2D* f1 = new TH2D("1",";Muon track length (rec) [cm];E_{#mu} (true) [GeV]",150, 0.0, 2000.0, 150, 0.0, 5.0);
207  TH2D* f2 = new TH2D("2",";Visible Hadronic E [GeV];E_{#nu} (true) - E_{#mu} (rec) [GeV]",77, hadronQEBins, 150, 0.0, 5.0);
208  TH2D* f3 = new TH2D("3",";Visible Hadronic E [GeV];E_{#nu} (true) - E_{#mu} (rec) [GeV]",130, hadronNonQEBins, 150, 0.0, 5.0);
209  f1->GetXaxis()->SetNoExponent(kTRUE);
210 
211  TH1D* f4 = new TH1D("4",";(Reco - True E_{#mu})/True E_{#mu};Slices",150, -1.0, 1.0);
212  TH2D* f5 = new TH2D("5",";Muon track length (rec) [cm];(Reco - True E_{#mu})/True E_{#mu}",150, 0.0, 2000.0, 150, -1.0, 1.0);
213  f5->GetXaxis()->SetNoExponent(kTRUE);
214 
215  TH1D* f6 = new TH1D("6",";(Reco - True Hadron Energy)/True Hadron Energy;Slices",150, -1.0, 1.0);
216  TH2D* f7 = new TH2D("7",";Reco Hadron Energy [GeV];(Reco - True Hadron E)/True Hadron E",150, 0.0, 5.0, 150, -1.0, 1.0);
217  f7->GetXaxis()->SetNoExponent(kTRUE);
218 
219  TH1D* f8 = new TH1D("8",";(Reco - True Hadron Energy)/True Hadron Energy;Slices",150, -1.0, 1.0);
220  TH2D* f9 = new TH2D("9",";Reco Hadron Energy [GeV];(Reco - True Hadron E)/True Hadron E",150, 0.0, 5.0, 150, -1.0, 1.0);
221  f9->GetXaxis()->SetNoExponent(kTRUE);
222 
223  TH1D* f10 = new TH1D("10",";(Reco - True E_{#nu})/True E_{#nu};Slices",150, -1.0, 1.0);
224  TH2D* f11 = new TH2D("11",";E_{#nu} (rec) [GeV];(Reco - True E_{#nu}/True E_{#nu}",150, 0.0, 5.0, 150, -1.0, 1.0);
225  f11->GetXaxis()->SetNoExponent(kTRUE);
226 
227  TH1D* f12 = new TH1D("12",";(Reco - True E_{#nu})/True E_{#nu};Slices",150, -1.0, 1.0);
228  TH2D* f13 = new TH2D("13",";E_{#nu} (rec) [GeV];(Reco - True E_{#nu})/True E_{#nu}",150, 0.0, 5.0, 150, -1.0, 1.0);
229  f13->GetXaxis()->SetNoExponent(kTRUE);
230 
231  TH2D* f14 = new TH2D("14",";E_{#mu} (rec) [GeV];E_{#mu} (true)",150, 0.0, 5.0, 150, 0.0, 5.0);
232  TH2D* f15 = new TH2D("15",";Reco Hadronic Energy [GeV];True Hadronic Energy [GeV]",150, 0.0, 5.0, 150, 0.0, 5.0);
233  TH2D* f16 = new TH2D("16",";Reco Hadronic Energy [GeV];True Hadronic Energy [GeV]",150, 0.0, 5.0, 150, 0.0, 5.0);
234  TH2D* f17 = new TH2D("17",";E_{#nu} (rec) [GeV];E_{#nu} (true) [GeV]",150, 0.0, 5.0, 150, 0.0, 5.0);
235  TH2D* f18 = new TH2D("18",";E_{#nu} (rec) [GeV];E_{#nu} (true) [GeV]",150, 0.0, 5.0, 150, 0.0, 5.0);
236 
237  // All CC
238  TH2D* f19 = new TH2D("19",";Visible Hadronic E [GeV];E_{#nu} (true) - E_{#mu} (rec) [GeV]",130, hadronCCBins, 150, 0.0, 5.0);
239  TH1D* f20 = new TH1D("20",";(Reco - True Hadron Energy)/True Hadron Energy;Slices",150, -1.0, 1.0);
240  TH2D* f21 = new TH2D("21",";Reco Hadron Energy [GeV];(Reco - True Hadron E)/True Hadron E",150, 0.0, 5.0, 150, -1.0, 1.0);
241  TH1D* f22 = new TH1D("22",";(Reco - True E_{#nu})/True E_{#nu};Slices",150, -1.0, 1.0);
242  TH2D* f23 = new TH2D("23",";E_{#nu} (rec) [GeV];(Reco - True E_{#nu})/True E_{#nu}",150, 0.0, 5.0, 150, -1.0, 1.0);
243  TH2D* f24 = new TH2D("24",";Reco Hadronic Energy [GeV];True Hadronic Energy [GeV]",150, 0.0, 5.0, 150, 0.0, 5.0);
244  TH2D* f25 = new TH2D("25",";E_{#nu} (rec) [GeV];E_{#nu} (true) [GeV]",150, 0.0, 5.0, 150, 0.0, 5.0);
245  f21->GetXaxis()->SetNoExponent(kTRUE);
246  f23->GetXaxis()->SetNoExponent(kTRUE);
247  //gStyle->SetOptStat(0);
248 
249  for (int i=0; i<tree->GetEntries(); ++i){
250  //for (int i=0; i<1000; ++i){ //Use this instead if you just want to do a quick test
251 
252  if (i%50000==0)
253  std::cout << "Entry: " << i << std::endl;
254 
255  tree->GetEntry(i);
256 
257  bool tracksExist = nKalman>0;
258  bool decentTrackExists = bestPID!=999;
259  bool energyExists = CalCCEnergy!=-5;
260  bool lessThan5GeV = TrueNuE < 5;
261 
262  if (tracksExist && decentTrackExists && energyExists && trackMatchesMuon){
263 
264  bool neutrinoEffPur = NuPurity == 1;
265  bool numuCC = isCC && (abs(PDGCode) == 14);
266  bool isQE = TrueQEInt;
267  bool muonEff = trackMuonEffHits >= 0.8;
268  bool muonPur = trackMuonPurHits >= 0.8;
269 
270  //Use the first method if have re-run files with muon E; otherwise use second to do on fly
271  //double muonRecoE = rec->energy.mutrkE[best3DRemid].E;
272  double muonRecoE = MuonEFromTrackLength(trackLength);
273 
274  double recoHadQE = QEEnergy - muonRecoE;
275  double recoHadNonQE = NonQEEnergy - muonRecoE;
276  // double recoHadCC = CalCCEnergy - muonRecoE; // Adjust to the real estimator
277  double recoHadCC = ECCFromCalE(HadCalEnergy);
278  double CCEnergy = recoHadCC + muonRecoE;
279 
280  if (neutrinoEffPur && numuCC && goodMuon && bestTrack && muonEff && muonPur && contained){
281  //For muon, going to use QE and nonQE events together!
282  f1->Fill(trackLength,TrueMuonE);
283  if (lessThan5GeV){
284  f4->Fill((muonRecoE-TrueMuonE)/TrueMuonE);
285  }
286  f5->Fill(trackLength,(muonRecoE-TrueMuonE)/TrueMuonE);
287  f14->Fill(muonRecoE,TrueMuonE);
288 
289  // First fill the All-CC estimator histograms
290  f19->Fill(HadCalEnergy,(TrueNuE-muonRecoE));
291  if (lessThan5GeV){
292  f20->Fill((recoHadCC - (TrueNuE-muonRecoE))/(TrueNuE-muonRecoE));
293  }
294  f21->Fill(recoHadCC,(recoHadCC - (TrueNuE-muonRecoE))/(TrueNuE-muonRecoE));
295  if (lessThan5GeV){
296  f22->Fill((CCEnergy-TrueNuE)/TrueNuE);
297  }
298  f23->Fill(CCEnergy,(CCEnergy-TrueNuE)/TrueNuE);
299  f24->Fill(recoHadCC,TrueNuE-muonRecoE);
300  f25->Fill(CCEnergy,TrueNuE);
301 
302  if (isQE){
303  f2->Fill(HadCalEnergy,(TrueNuE-muonRecoE));
304  if (lessThan5GeV){
305  f6->Fill((recoHadQE - (TrueNuE-muonRecoE))/(TrueNuE-muonRecoE));
306  }
307  f7->Fill(recoHadQE,(recoHadQE - (TrueNuE-muonRecoE))/(TrueNuE-muonRecoE));
308  if (lessThan5GeV){
309  f10->Fill((QEEnergy-TrueNuE)/TrueNuE);
310  }
311  f11->Fill(QEEnergy,(QEEnergy-TrueNuE)/TrueNuE);
312  f15->Fill(recoHadQE,TrueNuE-muonRecoE);
313  f17->Fill(QEEnergy,TrueNuE);
314  }
315  else{
316  f3->Fill(HadCalEnergy,(TrueNuE-muonRecoE));
317  if (lessThan5GeV){
318  f8->Fill((recoHadNonQE - (TrueNuE-muonRecoE))/(TrueNuE-muonRecoE));
319  }
320  f9->Fill(recoHadNonQE,(recoHadNonQE - (TrueNuE-muonRecoE))/(TrueNuE-muonRecoE));
321  if (lessThan5GeV){
322  f12->Fill((NonQEEnergy-TrueNuE)/TrueNuE);
323  }
324  f13->Fill(NonQEEnergy,(NonQEEnergy-TrueNuE)/TrueNuE);
325  f16->Fill(recoHadNonQE,TrueNuE-muonRecoE);
326  f18->Fill(NonQEEnergy,TrueNuE);
327  }
328 
329  }//End of loop over things that pass all basic super pure cuts
330  }//End of loop over things with tracks and energy
331  }//End of loop over entries
332 
333 
334  TCanvas* c1 = new TCanvas("c1","c1");
335  f1->Draw("colz");
336  c1->Print("muonTrkLenVTrueEnergy.pdf");
337 
338  TCanvas* c2 = new TCanvas("c2","c2");
339  f2->Draw("colz");
340  c2->Print("QEHadronEVTrueEnergy.pdf");
341 
342  TCanvas* c3 = new TCanvas("c3","c3");
343  f3->Draw("colz");
344  c3->Print("NonQEHadronEVTrueEnergy.pdf");
345 
346  TCanvas* c4 = new TCanvas("c4","c4");
347  f4->Draw();
348  c4->Print("muon1dRes.pdf");
349 
350  TCanvas* c5 = new TCanvas("c5","c5");
351  f5->Draw("colz");
352  c5->Print("muon2dRes.pdf");
353 
354  TCanvas* c6 = new TCanvas("c6","c6");
355  f6->Draw();
356  c6->Print("hadqe1dRes.pdf");
357 
358  TCanvas* c7 = new TCanvas("c7","c7");
359  f7->Draw("colz");
360  c7->Print("hadqe2dRes.pdf");
361 
362  TCanvas* c8 = new TCanvas("c8","c8");
363  f8->Draw();
364  c8->Print("hadnonqe1dRes.pdf");
365 
366  TCanvas* c9 = new TCanvas("c9","c9");
367  f9->Draw("colz");
368  c9->Print("hadnonqe2dRes.pdf");
369 
370  TCanvas* c10 = new TCanvas("c10","c10");
371  f10->Draw();
372  c10->Print("qe1dRes.pdf");
373 
374  TCanvas* c11 = new TCanvas("c11","c11");
375  f11->Draw("colz");
376  c11->Print("qe2dRes.pdf");
377 
378  TCanvas* c12 = new TCanvas("c12","c12");
379  f12->Draw();
380  c12->Print("nonqe1dRes.pdf");
381 
382  TCanvas* c13 = new TCanvas("c13","c13");
383  f13->Draw("colz");
384  c13->Print("nonqe2dRes.pdf");
385 
386  TCanvas* c14 = new TCanvas("c14","c14");
387  f14->Draw("colz");
388  c14->Print("muonReco.pdf");
389 
390  TCanvas* c15 = new TCanvas("c15","c15");
391  f15->Draw("colz");
392  c15->Print("hadqeReco.pdf");
393 
394  TCanvas* c16 = new TCanvas("c16","c16");
395  f16->Draw("colz");
396  c16->Print("hadnonqeReco.pdf");
397 
398  TCanvas* c17 = new TCanvas("c17","c17");
399  f17->Draw("colz");
400  c17->Print("neutrinoqeReco.pdf");
401 
402  TCanvas* c18 = new TCanvas("c18","c18");
403  f18->Draw("colz");
404  c18->Print("neutrinononqeReco.pdf");
405 
406  // All CC
407  TCanvas* c19 = new TCanvas("c19","c19");
408  f19->Draw("colz");
409  c19->Print("CCHadronEVTrueEnergy.pdf");
410 
411  TCanvas* c20 = new TCanvas("c20","c20");
412  f20->Draw();
413  c20->Print("hadcc1dRes.pdf");
414 
415  TCanvas* c21 = new TCanvas("c21","c21");
416  f21->Draw("colz");
417  c21->Print("hadcc2dRes.pdf");
418 
419  TCanvas* c22 = new TCanvas("c22","c22");
420  f22->Draw();
421  c22->Print("cc1dRes.pdf");
422 
423  TCanvas* c23 = new TCanvas("c23","c23");
424  f23->Draw("colz");
425  c23->Print("cc2dRes.pdf");
426 
427  TCanvas* c24 = new TCanvas("c24","c24");
428  f24->Draw("colz");
429  c24->Print("hadccReco.pdf");
430 
431  TCanvas* c25 = new TCanvas("c25","c25");
432  f25->Draw("colz");
433  c25->Print("neutrinoccReco.pdf");
434 
435 
436  TFile f("2DPlotsForFitting.root","RECREATE");
437  f1->Write();
438  f2->Write();
439  f3->Write();
440  f4->Write();
441  f5->Write();
442  f6->Write();
443  f7->Write();
444  f8->Write();
445  f9->Write();
446  f10->Write();
447  f11->Write();
448  f12->Write();
449  f13->Write();
450  f14->Write();
451  f15->Write();
452  f16->Write();
453  f17->Write();
454  f18->Write();
455  f19->Write();
456  f20->Write();
457  f21->Write();
458  f22->Write();
459  f23->Write();
460  f24->Write();
461  f25->Write();
462 
463 }
Float_t f4
void abs(TH1 *hist)
Float_t f2
c2
Definition: demo5.py:33
Float_t f3
Float_t E
Definition: plot.C:20
Float_t f1
OStream cout
Definition: OStream.cxx:6
void getPlots()
c1
Definition: demo5.py:24
Double_t ECCFromCalE(double E)
Double_t MuonEFromTrackLength(double trkLen)