getPlotsFromNtupleND.C
Go to the documentation of this file.
1 //Run using: root -b -q -l getPlotsFromNtuple.C
2 
3 //This script should be run first in one's journey to find new fit constants. After you have a ton of stats, please
4 //hadd the histogram files for one input file. Then 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 //to loop over all the ntuple entries.
7 
8 //This needs to be run to determine the muon reco E. THEN run it again to get the hadronic colz plots. Do NOT
9 //do these simultaneously because it won't be consistant!!!!
10 //This means you need to re-create histogram files after you have committed 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 
25  double muonE = 0.0;
26 
27  if (trkLen <= 0.0) return muonE;
28 
29  //New tuning values rounded for errors
30  double stitch1 = 334; // cm
31  double stitch2 = 539; // cm
32  double stitch3 = 1064; // cm
33  double offset = 0.1503; // GeV
34  double slope1 = 0.001910; // GeV/cm
35  double slope2 = 0.001987; // GeV/cm
36  double slope3 = 0.002039; // GeV/cm
37  double slope4 = 0.002159; // GeV/cm
38 
39  if (trkLen < stitch1){
40  muonE = slope1*trkLen + offset;
41  }
42  else if (trkLen < stitch2){
43  muonE = slope2*trkLen + ((slope1-slope2)*stitch1 + offset);
44  }
45  else if (trkLen <stitch3){
46  muonE = slope3*trkLen + ((slope2-slope3)*stitch2 +(slope1-slope2)*stitch1 + offset);
47  }
48  else{
49  muonE = slope4*trkLen + ((slope3-slope4)*stitch3 +(slope2-slope3)*stitch2 +(slope1-slope2)*stitch1 + offset);
50  }
51 
52  return muonE;
53 
54 }//End of MuonEFromTrackLength
55 
56 double ActiveTrackLengthFromMuonE(double muonE)
57 {
58  //This function returns effective active-region track length, in cm, given the muon energy in GeV.
59  //It reverses a spline fit from May 2015.
60 
61  double trkLen = 0.0;
62 
63  //New tuning values rounded for errors
64  double FDMstitch1 = 334; // cm
65  double FDMstitch2 = 539; // cm
66  double FDMstitch3 = 1064; // cm
67  double FDMoffset = 0.1503; // GeV
68  double FDMslope1 = 0.001910; // GeV/cm
69  double FDMslope2 = 0.001987; // GeV/cm
70  double FDMslope3 = 0.002039; // GeV/cm
71  double FDMslope4 = 0.002159; // GeV/cm
72 
73  if (muonE <= 0.0) return trkLen;
74 
75  if (muonE <= FDMoffset ) return trkLen;
76 
77  else if (muonE < (FDMslope1*FDMstitch1 + FDMoffset)){
78  trkLen = (muonE - FDMoffset)/FDMslope1;
79  }
80 
81  else if (muonE < (FDMslope2*FDMstitch2 + ((FDMslope1-FDMslope2)*FDMstitch1 + FDMoffset))){
82  trkLen = (muonE - ((FDMslope1-FDMslope2)*FDMstitch1 + FDMoffset))/FDMslope2;
83  }
84 
85  else if (muonE < (FDMslope3*FDMstitch3 + ((FDMslope2-FDMslope3)*FDMstitch2 + (FDMslope1-FDMslope2)*FDMstitch1 + FDMoffset))){
86  trkLen = (muonE - ((FDMslope2-FDMslope3)*FDMstitch2 +(FDMslope1-FDMslope2)*FDMstitch1 + FDMoffset))/FDMslope3;
87  }
88 
89  else {
90  trkLen = (muonE - ((FDMslope3-FDMslope4)*FDMstitch3 + (FDMslope2-FDMslope3)*FDMstitch2 + (FDMslope1-FDMslope2)*FDMstitch1 + FDMoffset))/FDMslope4;
91  }
92 
93  return trkLen;
94 
95 }//End of ActiveTrackLengthFromMuonE
96 
97 double NDMuonEFromTrackLength(double actTrkLen, double catTrkLen, double trkCalAct, double trkCalTran, double trkCalCat)
98 {
99  //This function returns muon energy, in GeV, given the muon track length in cm.
100  //It distinguishes between the four types of muon tracks and calls the
101  //correct energy estimator for each. This fitting work was done in May 2015.
102 
103  double muonE = 0.0;
104 
105  if ((catTrkLen <= 0.0) && (actTrkLen <= 0.0)) return muonE;
106 
107  if (((trkCalAct > 0.0)||(trkCalTran > 0.0))&&(trkCalCat == 0.0)){
108  muonE = NDMuonEInActiveOnly(actTrkLen);
109  }//Entirely contained in active region and transition plane
110 
111  if (((trkCalAct > 0.0)||(trkCalTran > 0.0))&&(trkCalCat > 0.0)){
112  muonE = NDMuonEInActiveAndCatcher(actTrkLen,catTrkLen);
113  }//Length in active and catcher region
114 
115  return muonE;
116 
117 }//End of NDMuonEFromTrackLength
118 
119 double NDMuonEInActiveOnly(double actTrkLen)
120 {
121  //This function returns muon energy, in GeV, given the muon track length in cm.
122  //It is for ND muons that are entirely contained in the active region of the
123  //detector. It uses the same fit used for FD muons, since these are also
124  //entirely contained in the active material.
125 
126  double muonE = 0.0;
127 
128  if (actTrkLen <= 0.0) return muonE;
129 
130  muonE = MuonEFromTrackLength(actTrkLen);
131 
132  return muonE;
133 
134 }//End of NDMuonEInActiveOnly
135 
136 double NDMuonEInActiveAndCatcher(double actTrkLen, double catTrkLen)
137 {
138  //This function returns muon energy, in GeV, given the muon track length in cm.
139  //It is for ND muons that have energy in both the active and catcher region.
140  //It uses a linear fit to determine the muon energy when it enters the muon catcher.
141  //This fit is from May 2015..
142  //This is then converted into an effective active track length, which is added
143  //to the existing active track length for a total active track length.
144  //This length is then fed into the FD muon fit.
145 
146  double muonE = 0.0;
147 
148  if ((catTrkLen <= 0.0) && (actTrkLen <= 0.0)) return muonE;
149 
150  double offset = 0.152; // GeV
151  double slope = 0.00536; // GeV/cm
152 
153  double muonECat = slope*catTrkLen + offset;
154 
155  double catcherActLen = ActiveTrackLengthFromMuonE(muonECat);
156 
157  double totalActLen = catcherActLen + actTrkLen;
158 
159  muonE = MuonEFromTrackLength(totalActLen);
160 
161  return muonE;
162 
163 }//End of NDMuonEInActiveAndCatcher
164 
165 double NDMuonEInActiveAndCatcherJustCatcherE(double actTrkLen, double catTrkLen)
166 {
167  //This function returns muon energy, in GeV, given the muon track length in cm.
168  //It is for ND muons that have energy in both the active and catcher region.
169  //It uses a linear fit to determine the muon energy when it enters the muon catcher.
170  //This fit is from May 2015..
171  //This is then converted into an effective active track length, which is added
172  //to the existing active track length for a total active track length.
173  //This length is then fed into the FD muon fit.
174 
175  double muonE = 0.0;
176 
177  if ((catTrkLen <= 0.0) && (actTrkLen <= 0.0)) return muonE;
178 
179  double offset = 0.152; // GeV
180  double slope = 0.00536; // GeV/cm
181 
182  double muonECat = slope*catTrkLen + offset;
183 
184  return muonECat;
185 
186 }//End of NDMuonEInActiveAndCatcher
187 
188 double NDMuonEInActiveAndCatcherEffLen(double actTrkLen, double catTrkLen)
189 {
190  //This function returns muon energy, in GeV, given the muon track length in cm.
191  //It is for ND muons that have energy in both the active and catcher region.
192  //It uses a linear fit to determine the muon energy when it enters the muon catcher.
193  //This fit is from May 2015..
194  //This is then converted into an effective active track length, which is added
195  //to the existing active track length for a total active track length.
196  //This length is then fed into the FD muon fit.
197 
198  double muonE = 0.0;
199 
200  if ((catTrkLen <= 0.0) && (actTrkLen <= 0.0)) return muonE;
201 
202  double offset = 0.152; // GeV
203  double slope = 0.00536; // GeV/cm
204 
205  double muonECat = slope*catTrkLen + offset;
206 
207  double catcherActLen = ActiveTrackLengthFromMuonE(muonECat);
208 
209  double totalActLen = catcherActLen + actTrkLen;
210 
211  return totalActLen;
212 
213 }//End of NDMuonEInActiveAndCatcher
214 
215 
217 {
218 
219  std::cout<<"Hello!"<<std::endl;
220 
221  TFile* file = new TFile("/local/nova/users/lein/NumuE/checkND.root","READ");
222  TNtuple* tree = (TNtuple*)file->Get("numuana/NumuETree");
223 
224  float CalCCEnergy;
225  float QEEnergy;
226  float NonQEEnergy;
227  float CCEnergy;
228  float AngleQEEnergy;
229  float AngleQEError;
230  float HadCalEnergy;
231  float HadTrkEnergy;
232  float ndTrkLenAct;
233  float ndTrkLenCat;
234  float ndTrkCalAct;
235  float ndTrkCalTran;
236  float ndTrkCalCat;
237  float ndHadCalAct;
238  float ndHadCalTran;
239  float ndHadCalCat;
240  float ndHadTrkAct;
241  float ndHadTrkTran;
242  float ndHadTrkCat;
243  float ndTrkTranX;
244  float ndTrkTranY;
245  bool ND;
246  bool FD;
247  bool isCC;
248  bool TrueQEInt;
249  bool trackMatchesMuon;
250  bool contained;
251  int Run;
252  int SubRun;
253  int Evt;
254  int SubEvt;
255  int goodMuon;
256  int bestTrack;
257  int PDGCode;
258  int nKalman;
259  int nCosmic;
260  int qepidntrk;
261  int firstplane;
262  int lastplane;
263  int nplanestofront;
264  int nplanestoback;
265  unsigned int mincellsfromedge;
266  int nslicehit;
267  int nslicecontplanes;
268  double NuPurity;
269  double NuEfficiency;
270  double trackLength;
271  double catcherE;
272  double bestPID;
273  double qepidPID;
274  double cosrejanglekal;
275  double cosrejnumucontpid;
276  double trackMuonEffHits;
277  double trackMuonPurHits;
278  double start_x_reco;
279  double start_y_reco;
280  double start_z_reco;
281  double end_x_reco;
282  double end_y_reco;
283  double end_z_reco;
284  double TrueMuonE;
285  double RecoMuonE;
286  double TrueNuE;
287  double TrueTrackLength;
288 
289 
290 
291  tree->SetBranchAddress("CalCCEnergy", &CalCCEnergy);
292  tree->SetBranchAddress("QEEnergy", &QEEnergy);
293  tree->SetBranchAddress("NonQEEnergy", &NonQEEnergy);
294  tree->SetBranchAddress("CCEnergy", &CCEnergy);
295  tree->SetBranchAddress("AngleQEEnergy", &AngleQEEnergy);
296  tree->SetBranchAddress("AngleQEError", &AngleQEError);
297  tree->SetBranchAddress("HadCalEnergy", &HadCalEnergy);
298  tree->SetBranchAddress("HadTrkEnergy", &HadTrkEnergy);
299  tree->SetBranchAddress("ndTrkLenAct", &ndTrkLenAct);
300  tree->SetBranchAddress("ndTrkLenCat", &ndTrkLenCat);
301  tree->SetBranchAddress("ndTrkCalAct", &ndTrkCalAct);
302  tree->SetBranchAddress("ndTrkCalTran", &ndTrkCalTran);
303  tree->SetBranchAddress("ndTrkCalCat", &ndTrkCalCat);
304  tree->SetBranchAddress("ndHadCalAct", &ndHadCalAct);
305  tree->SetBranchAddress("ndHadCalTran", &ndHadCalTran);
306  tree->SetBranchAddress("ndHadCalCat", &ndHadCalCat);
307  tree->SetBranchAddress("ndHadTrkAct", &ndHadTrkAct);
308  tree->SetBranchAddress("ndHadTrkTran", &ndHadTrkTran);
309  tree->SetBranchAddress("ndHadTrkCat", &ndHadTrkCat);
310  tree->SetBranchAddress("ndTrkTranX", &ndTrkTranX);
311  tree->SetBranchAddress("ndTrkTranY", &ndTrkTranY);
312  tree->SetBranchAddress("ND", &ND);
313  tree->SetBranchAddress("FD", &FD);
314  tree->SetBranchAddress("isCC", &isCC);
315  tree->SetBranchAddress("TrueQEInt", &TrueQEInt);
316  tree->SetBranchAddress("trackMatchesMuon", &trackMatchesMuon);
317  tree->SetBranchAddress("contained", &contained);
318  tree->SetBranchAddress("Run", &Run);
319  tree->SetBranchAddress("SubRun", &SubRun);
320  tree->SetBranchAddress("Evt", &Evt);
321  tree->SetBranchAddress("SubEvt", &SubEvt);
322  tree->SetBranchAddress("goodMuon", &goodMuon);
323  tree->SetBranchAddress("bestTrack", &bestTrack);
324  tree->SetBranchAddress("PDGCode", &PDGCode);
325  tree->SetBranchAddress("nKalman", &nKalman);
326  tree->SetBranchAddress("nCosmic", &nCosmic);
327  tree->SetBranchAddress("qepidntrk", &qepidntrk);
328  tree->SetBranchAddress("firstplane", &firstplane);
329  tree->SetBranchAddress("lastplane", &lastplane);
330  tree->SetBranchAddress("nplanestofront", &nplanestofront);
331  tree->SetBranchAddress("nplanestoback", &nplanestoback);
332  tree->SetBranchAddress("mincellsfromedge", &mincellsfromedge);
333  tree->SetBranchAddress("nslicehit", &nslicehit);
334  tree->SetBranchAddress("nslicecontplanes", &nslicecontplanes);
335  tree->SetBranchAddress("NuPurity", &NuPurity);
336  tree->SetBranchAddress("NuEfficiency", &NuEfficiency);
337  tree->SetBranchAddress("trackLength", &trackLength);
338  tree->SetBranchAddress("catcherE", &catcherE);
339  tree->SetBranchAddress("bestPID", &bestPID);
340  tree->SetBranchAddress("qepidPID", &qepidPID);
341  tree->SetBranchAddress("cosrejanglekal", &cosrejanglekal);
342  tree->SetBranchAddress("cosrejnumucontpid", &cosrejnumucontpid);
343  tree->SetBranchAddress("trackMuonEffHits", &trackMuonEffHits);
344  tree->SetBranchAddress("trackMuonPurHits", &trackMuonPurHits);
345  tree->SetBranchAddress("start_x_reco", &start_x_reco);
346  tree->SetBranchAddress("start_y_reco", &start_y_reco);
347  tree->SetBranchAddress("start_z_reco", &start_z_reco);
348  tree->SetBranchAddress("end_x_reco", &end_x_reco);
349  tree->SetBranchAddress("end_y_reco", &end_y_reco);
350  tree->SetBranchAddress("end_z_reco", &end_z_reco);
351  tree->SetBranchAddress("TrueMuonE", &TrueMuonE);
352  tree->SetBranchAddress("RecoMuonE", &RecoMuonE);
353  tree->SetBranchAddress("TrueNuE", &TrueNuE);
354  tree->SetBranchAddress("TrueTrackLength", &TrueTrackLength);
355 
356  //Making array for variable binning
357  double hadronQEBins[58];
358  double hadronQEAxis = 0.0;
359  for(int i = 0; i < 58; ++i){
360  hadronQEBins[i] = hadronQEAxis;
361  if (hadronQEAxis < 0.5){hadronQEAxis = hadronQEAxis + 0.01;}
362  else if (hadronQEAxis < 0.75){hadronQEAxis = hadronQEAxis + 0.05;}
363  else if (hadronQEAxis < 1.25){hadronQEAxis = hadronQEAxis + 0.5;}
364  else {hadronQEAxis = hadronQEAxis + 0.75;}
365  }
366 
367  double hadronNonQEBins[116];
368  double hadronNonQEAxis = 0.0;
369  for(int i = 0; i < 116; ++i){
370  hadronNonQEBins[i] = hadronNonQEAxis;
371  if (hadronNonQEAxis < 1.0){hadronNonQEAxis = hadronNonQEAxis + 0.01;}
372  else if (hadronNonQEAxis < 1.5){hadronNonQEAxis = hadronNonQEAxis + 0.05;}
373  else {hadronNonQEAxis = hadronNonQEAxis + 0.1;}
374  }
375 
376  //All Active Case
377  TH2D* f1 = new TH2D("1",";Reco Muon Track Length (cm);True Muon Energy (GeV)",150, 0.0, 1500.0, 150, 0.0, 5.0);
378  //All Catcher Case
379  TH2D* f2 = new TH2D("2",";Reco Muon Track Length (cm);True Muon Energy (GeV)",150, 0.0, 400.0, 150, 0.0, 3.0);
380  //Catcher Energy for Active + Catcher Case
381  TH2D* f3 = new TH2D("3",";Reco Muon Track Length (cm);True Muon Energy (GeV)",150, 0.0, 400.0, 150, 0.0, 3.0);
382  //Active + Catcher Case
383  TH2D* f4 = new TH2D("4",";Reco Muon Track Length (cm);True Muon Energy (GeV)",150, 0.0, 1500.0, 150, 0.0, 5.0);
384  //Hadronic QE
385  TH2D* f6 = new TH2D("6",";Visible Hadronic E (GeV);True Neutrino E - Reco Muon E (GeV)",57, hadronQEBins, 150, 0.0, 5.0);
386  //Hadronic NonQE
387  TH2D* f7 = new TH2D("7",";Visible Hadronic E (GeV);True Neutrino E - Reco Muon E (GeV)",115, hadronNonQEBins, 150, 0.0, 5.0);
388  //Hadronic CC
389  TH2D* f8 = new TH2D("8",";Visible Hadronic E (GeV);True Neutrino E - Reco Muon E (GeV)",115, hadronNonQEBins, 150, 0.0, 5.0);
390 
391  f1->GetXaxis()->SetNoExponent(kTRUE);
392  f4->GetXaxis()->SetNoExponent(kTRUE);
393 
394  //Validation plots for all-acitve muons
395  TH1D* f9 = new TH1D("9",";(Reco - True Muon Energy)/True Muon Energy;Slices",150, -1.0, 1.0);
396  TH2D* f10 = new TH2D("10",";Reco Muon Track Length (cm);(Reco - True Muon E)/True Muon E",150, 0.0, 1500.0, 150, -1.0, 1.0);
397  f10->GetXaxis()->SetNoExponent(kTRUE);
398  TH2D* f11 = new TH2D("11",";Reco Muon Energy (GeV);True Muon Energy (GeV)",150, 0.0, 5.0, 150, 0.0, 5.0);
399  //Validation plots for act+cat muons
400  TH1D* f12 = new TH1D("12",";(Reco - True Muon Energy)/True Muon Energy;Slices",150, -1.0, 1.0);
401  TH2D* f13 = new TH2D("13",";Reco Muon Track Length (cm);(Reco - True Muon E)/True Muon E",150, 0.0, 400.0, 150, -1.0, 1.0);
402  f13->GetXaxis()->SetNoExponent(kTRUE);
403  TH2D* f14 = new TH2D("14",";Reco Muon Energy (GeV);True Muon Energy (GeV)",150, 0.0, 3.0, 150, 0.0, 3.0);
404  TH1D* f15 = new TH1D("15",";(Reco - True Muon Energy)/True Muon Energy;Slices",150, -1.0, 1.0);
405  TH2D* f16 = new TH2D("16",";Reco Muon Track Length (cm);(Reco - True Muon E)/True Muon E",150, 0.0, 1500.0, 150, -1.0, 1.0);
406  f16->GetXaxis()->SetNoExponent(kTRUE);
407  TH2D* f17 = new TH2D("17",";Reco Muon Energy (GeV);True Muon Energy (GeV)",150, 0.0, 5.0, 150, 0.0, 5.0);
408  //Validation plots for all muons
409  TH1D* f18 = new TH1D("18",";(Reco - True Muon Energy)/True Muon Energy;Slices",150, -1.0, 1.0);
410  TH2D* f19 = new TH2D("19",";Reco Muon Track Length (cm);(Reco - True Muon E)/True Muon E",150, 0.0, 1500.0, 150, -1.0, 1.0);
411  f19->GetXaxis()->SetNoExponent(kTRUE);
412  TH2D* f20 = new TH2D("20",";Reco Muon Energy (GeV);True Muon Energy (GeV)",150, 0.0, 5.0, 150, 0.0, 5.0);
413  //Validation plots for hadronic QE
414  TH1D* f21 = new TH1D("21",";(Reco - Desired Hadronic E)/Desired Hadronic E;Slices",150, -1.0, 1.0);
415  TH2D* f22 = new TH2D("22",";Reco Hadronic Energy (GeV);(Reco - Desired Had E)/Desired Had E",150, 0.0, 5.0, 150, -1.0, 1.0);
416  TH2D* f23 = new TH2D("23",";Reco Hadronic Energy (GeV);Desired Hadronic Energy (GeV)",150, 0.0, 5.0, 150, 0.0, 5.0);
417  //Validation plots for hadronic NonQE
418  TH1D* f24 = new TH1D("24",";(Reco - Desired Hadronic E)/Desired Hadronic E;Slices",150, -1.0, 1.0);
419  TH2D* f25 = new TH2D("25",";Reco Hadronic Energy (GeV);(Reco - Desired Had E)/Desired Had E",150, 0.0, 5.0, 150, -1.0, 1.0);
420  TH2D* f26 = new TH2D("26",";Reco Hadronic Energy (GeV);Desired Hadronic Energy (GeV)",150, 0.0, 5.0, 150, 0.0, 5.0);
421  //Validation plots for hadronic CC
422  TH1D* f27 = new TH1D("27",";(Reco - Desired Hadronic E)/Desired Hadronic E;Slices",150, -1.0, 1.0);
423  TH2D* f28 = new TH2D("28",";Reco Hadronic Energy (GeV);(Reco - Desired Had E)/Desired Had E",150, 0.0, 5.0, 150, -1.0, 1.0);
424  TH2D* f29 = new TH2D("29",";Reco Hadronic Energy (GeV);Desired Hadronic Energy (GeV)",150, 0.0, 5.0, 150, 0.0, 5.0);
425  //Validation plots for QE neutrinos
426  TH1D* f30 = new TH1D("30",";(Reco - True Neutrino Energy)/True Neutrino Energy;Slices",150, -1.0, 1.0);
427  TH2D* f31 = new TH2D("31",";Reco Neutrino Energy (GeV);(Reco - True Neutrino E)/True Neutrino E",150, 0.0, 5.0, 150, -1.0, 1.0);
428  TH2D* f32 = new TH2D("32",";Reco Neutrino Energy (GeV);True Neutrino Energy (GeV)",150, 0.0, 5.0, 150, 0.0, 5.0);
429  //Validation plots for NonQE neutrinos
430  TH1D* f33 = new TH1D("33",";(Reco - True Neutrino Energy)/True Neutrino Energy;Slices",150, -1.0, 1.0);
431  TH2D* f34 = new TH2D("34",";Reco Neutrino Energy (GeV);(Reco - True Neutrino E)/True Neutrino E",150, 0.0, 5.0, 150, -1.0, 1.0);
432  TH2D* f35 = new TH2D("35",";Reco Neutrino Energy (GeV);True Neutrino Energy (GeV)",150, 0.0, 5.0, 150, 0.0, 5.0);
433  //Validation plots for CC neutrinos
434  TH1D* f36 = new TH1D("36",";(Reco - True Neutrino Energy)/True Neutrino Energy;Slices",150, -1.0, 1.0);
435  TH2D* f37 = new TH2D("37",";Reco Neutrino Energy (GeV);(Reco - True Neutrino E)/True Neutrino E",150, 0.0, 5.0, 150, -1.0, 1.0);
436  TH2D* f38 = new TH2D("38",";Reco Neutrino Energy (GeV);True Neutrino Energy (GeV)",150, 0.0, 5.0, 150, 0.0, 5.0);
437 
438  /*
439  TH2D* f1 = new TH2D("1",";Reco Muon Track Length (cm);True Muon Energy (GeV)",150, 0.0, 2000.0, 150, 0.0, 5.0);
440  TH2D* f2 = new TH2D("2",";Visible Hadronic E (GeV);True Neutrino E - Reco Muon E (GeV)",77, hadronQEBins, 150, 0.0, 5.0);
441  TH2D* f3 = new TH2D("3",";Visible Hadronic E (GeV);True Neutrino E - Reco Muon E (GeV)",130, hadronNonQEBins, 150, 0.0, 5.0);
442  TH2D* f19 = new TH2D("19",";Visible Hadronic E (GeV);True Neutrino E - Reco Muon E (GeV)",130, hadronNonQEBins, 150, 0.0, 5.0);
443  f1->GetXaxis()->SetNoExponent(kTRUE);
444 
445  TH1D* f4 = new TH1D("4",";(Reco - True Muon Energy)/True Muon Energy;Slices",150, -1.0, 1.0);
446  TH2D* f5 = new TH2D("5",";Reco Muon Track Length (cm);(Reco - True Muon E)/True Muon E",150, 0.0, 2000.0, 150, -1.0, 1.0);
447  f5->GetXaxis()->SetNoExponent(kTRUE);
448 
449  TH1D* f6 = new TH1D("6",";(Reco - Desired Hadronic E)/Desired Hadronic E;Slices",150, -1.0, 1.0);
450  TH2D* f7 = new TH2D("7",";Reco Hadronic Energy (GeV);(Reco - Desired Had E)/Desired Had E",150, 0.0, 5.0, 150, -1.0, 1.0);
451  f7->GetXaxis()->SetNoExponent(kTRUE);
452 
453  TH1D* f8 = new TH1D("8",";(Reco - Desired Hadronic E)/Desired Hadronic E;Slices",150, -1.0, 1.0);
454  TH2D* f9 = new TH2D("9",";Reco Hadronic Energy (GeV);(Reco - Desired Had E)/Desired Had E",150, 0.0, 5.0, 150, -1.0, 1.0);
455  f9->GetXaxis()->SetNoExponent(kTRUE);
456 
457  TH1D* f20 = new TH1D("20",";(Reco - Desired Hadronic E)/Desired Hadronic E;Slices",150, -1.0, 1.0);
458  TH2D* f21 = new TH2D("21",";Reco Hadronic Energy (GeV);(Reco - Desired Had E)/Desired Had E",150, 0.0, 5.0, 150, -1.0, 1.0);
459  f21->GetXaxis()->SetNoExponent(kTRUE);
460 
461  TH1D* f10 = new TH1D("10",";(Reco - True Neutrino Energy)/True Neutrino Energy;Slices",150, -1.0, 1.0);
462  TH2D* f11 = new TH2D("11",";Reco Neutrino Energy (GeV);(Reco - True Neutrino E)/True Neutrino E",150, 0.0, 5.0, 150, -1.0, 1.0);
463  f11->GetXaxis()->SetNoExponent(kTRUE);
464 
465  TH1D* f12 = new TH1D("12",";(Reco - True Neutrino Energy)/True Neutrino Energy;Slices",150, -1.0, 1.0);
466  TH2D* f13 = new TH2D("13",";Reco Neutrino Energy (GeV);(Reco - True Neutrino E)/True Neutrino E",150, 0.0, 5.0, 150, -1.0, 1.0);
467  f13->GetXaxis()->SetNoExponent(kTRUE);
468 
469  TH1D* f22 = new TH1D("22",";(Reco - True Neutrino Energy)/True Neutrino Energy;Slices",150, -1.0, 1.0);
470  TH2D* f23 = new TH2D("23",";Reco Neutrino Energy (GeV);(Reco - True Neutrino E)/True Neutrino E",150, 0.0, 5.0, 150, -1.0, 1.0);
471  f23->GetXaxis()->SetNoExponent(kTRUE);
472 
473  TH2D* f14 = new TH2D("14",";Reco Muon Energy (GeV);True Muon Energy (GeV)",150, 0.0, 5.0, 150, 0.0, 5.0);
474  TH2D* f15 = new TH2D("15",";Reco Hadronic Energy (GeV);Desired Hadronic Energy (GeV)",150, 0.0, 5.0, 150, 0.0, 5.0);
475  TH2D* f16 = new TH2D("16",";Reco Hadronic Energy (GeV);Desired Hadronic Energy (GeV)",150, 0.0, 5.0, 150, 0.0, 5.0);
476  TH2D* f24 = new TH2D("24",";Reco Hadronic Energy (GeV);Desired Hadronic Energy (GeV)",150, 0.0, 5.0, 150, 0.0, 5.0);
477  TH2D* f17 = new TH2D("17",";Reco Neutrino Energy (GeV);True Neutrino Energy (GeV)",150, 0.0, 5.0, 150, 0.0, 5.0);
478  TH2D* f18 = new TH2D("18",";Reco Neutrino Energy (GeV);True Neutrino Energy (GeV)",150, 0.0, 5.0, 150, 0.0, 5.0);
479  TH2D* f25 = new TH2D("25",";Reco Neutrino Energy (GeV);True Neutrino Energy (GeV)",150, 0.0, 5.0, 150, 0.0, 5.0);
480  */
481  gStyle->SetOptStat(0);
482 
483  for (int i=0; i<tree->GetEntries(); ++i){
484  //for (int i=0; i<1000; ++i){ //Use this instead if you just want to do a quick test
485  tree->GetEntry(i);
486 
487  bool tracksExist = nKalman>0;
488  bool energyExists = CalCCEnergy!=-5;
489 
490  if (tracksExist&&energyExists){
491 
492  bool neutrinoEffPur = ((NuPurity>=0.8)&&(NuEfficiency>=0.8));
493  bool numuCC = (isCC)&&(abs(PDGCode) == 14);
494  bool isQE = TrueQEInt;
495  bool bestTrackBool = (bestTrack != -1);
496  bool muonEff = trackMuonEffHits >= 0.8;
497  bool muonPur = trackMuonPurHits >= 0.8;
498  bool lessThan5GeV = TrueNuE <= 5.0;
499 
500  bool kNumuQuality = ((NonQEEnergy > 0)&&(bestPID > 0)&&(nslicehit > 20)&&(nslicecontplanes > 4)&&(nCosmic > 0));
501  bool kNumuNCRej = (bestPID > 0.75);
502 
503  double hadActive = ndHadCalAct + ndHadTrkAct;
504  double hadTran = ndHadCalTran + ndHadTrkTran;
505  double hadCatcher = ndHadCalCat + ndHadTrkCat;
506  double hadAll = hadActive + hadTran + hadCatcher;
507 
508  //Use the first method if have re-run files with muon E; otherwise use second to do on fly
509  double muonRecoE = RecoMuonE;
510  //double muonRecoE = NDMuonEFromTrackLength(ndTrkLenAct,ndTrkLenCat,ndTrkCalAct,ndTrkCalTran,ndTrkCalCat);
511  double muonRecoEInCatcher = NDMuonEInActiveAndCatcherJustCatcherE(ndTrkLenAct,ndTrkLenCat);
512  double muonActCatEffLen = NDMuonEInActiveAndCatcherEffLen(ndTrkLenAct,ndTrkLenCat);
513 
514  double recoHadQE = QEEnergy - muonRecoE;
515  double recoHadNonQE = NonQEEnergy - muonRecoE;
516  double recoHadCC = CCEnergy - muonRecoE;
517 
518  if (numuCC&&goodMuon&&bestTrackBool&&contained&&kNumuQuality&&kNumuNCRej){
519 
520  //For muon, going to use QE and nonQE events together!
521  if (((ndTrkCalAct > 0.0)||(ndTrkCalTran > 0.0))&&(ndTrkCalCat == 0.0)){ //All Active
522  f1->Fill(ndTrkLenAct,TrueMuonE);
523  if (lessThan5GeV){
524  f9->Fill((muonRecoE-TrueMuonE)/TrueMuonE);
525  }
526  f10->Fill(trackLength,(muonRecoE-TrueMuonE)/TrueMuonE);
527  f11->Fill(muonRecoE,TrueMuonE);
528  f19->Fill(trackLength,(muonRecoE-TrueMuonE)/TrueMuonE);
529  }
530 
531  if ((ndTrkCalAct == 0.0)&&(ndTrkCalTran == 0.0)&&(ndTrkCalCat > 0.0)){ //All Catcher
532  f2->Fill(ndTrkLenCat,TrueMuonE);
533  }
534 
535  if (((ndTrkCalAct > 0.0)||(ndTrkCalTran > 0.0))&&(ndTrkCalCat > 0.0)){ //Active + Catcher
536  f3->Fill(ndTrkLenCat,catcherE);
537  //Plots to evaluate just the catcher portion of track
538  if (lessThan5GeV){
539  f12->Fill((muonRecoEInCatcher-catcherE)/catcherE);
540  }
541  f13->Fill(ndTrkLenCat,(muonRecoEInCatcher-catcherE)/catcherE);
542  f14->Fill(muonRecoEInCatcher,catcherE);
543  if (lessThan5GeV){
544  f15->Fill((muonRecoE-TrueMuonE)/TrueMuonE);
545  }
546  f16->Fill(muonActCatEffLen,(muonRecoE-TrueMuonE)/TrueMuonE);
547  f17->Fill(muonRecoE,TrueMuonE);
548  f19->Fill(muonActCatEffLen,(muonRecoE-TrueMuonE)/TrueMuonE);
549  }
550 
551  if (lessThan5GeV){
552  f18->Fill((muonRecoE-TrueMuonE)/TrueMuonE);
553  }
554  f20->Fill(muonRecoE,TrueMuonE);
555 
556  f8->Fill(hadAll,(TrueNuE-muonRecoE));
557  if (lessThan5GeV){
558  f27->Fill((recoHadCC-(TrueNuE-muonRecoE))/(TrueNuE-muonRecoE));
559  }
560  f28->Fill(recoHadCC,(recoHadCC-(TrueNuE-muonRecoE))/(TrueNuE-muonRecoE));
561  f29->Fill(recoHadCC,TrueNuE-muonRecoE);
562  if (lessThan5GeV){
563  f36->Fill((CCEnergy-TrueNuE)/TrueNuE);
564  }
565  f37->Fill(CCEnergy,(CCEnergy-TrueNuE)/TrueNuE);
566  f38->Fill(CCEnergy,TrueNuE);
567 
568  if (isQE){
569  f6->Fill(hadAll,(TrueNuE-muonRecoE));
570  if (lessThan5GeV){
571  f21->Fill((recoHadQE-(TrueNuE-muonRecoE))/(TrueNuE-muonRecoE));
572  }
573  f22->Fill(recoHadQE,(recoHadQE-(TrueNuE-muonRecoE))/(TrueNuE-muonRecoE));
574  f23->Fill(recoHadQE,TrueNuE-muonRecoE);
575  if (lessThan5GeV){
576  f30->Fill((QEEnergy-TrueNuE)/TrueNuE);
577  }
578  f31->Fill(QEEnergy,(QEEnergy-TrueNuE)/TrueNuE);
579  f32->Fill(QEEnergy,TrueNuE);
580  }
581  else{
582  f7->Fill(hadAll,(TrueNuE-muonRecoE));
583  if (lessThan5GeV){
584  f24->Fill((recoHadNonQE-(TrueNuE-muonRecoE))/(TrueNuE-muonRecoE));
585  }
586  f25->Fill(recoHadNonQE,(recoHadNonQE-(TrueNuE-muonRecoE))/(TrueNuE-muonRecoE));
587  f26->Fill(recoHadNonQE,TrueNuE-muonRecoE);
588  if (lessThan5GeV){
589  f33->Fill((NonQEEnergy-TrueNuE)/TrueNuE);
590  }
591  f34->Fill(NonQEEnergy,(NonQEEnergy-TrueNuE)/TrueNuE);
592  f35->Fill(NonQEEnergy,TrueNuE);
593  }
594 
595 
596  /*
597 
598 
599  f19->Fill(hadE,(TrueNuE-muonRecoE));
600  if (lessThan5GeV){
601  f20->Fill((recoHadCC - (TrueNuE-muonRecoE))/(TrueNuE-muonRecoE));
602  }
603  f21->Fill(recoHadCC,(recoHadCC - (TrueNuE-muonRecoE))/(TrueNuE-muonRecoE));
604  if (lessThan5GeV){
605  f22->Fill((CCEnergy-TrueNuE)/TrueNuE);
606  }
607  f23->Fill(CCEnergy,(CCEnergy-TrueNuE)/TrueNuE);
608  f24->Fill(recoHadCC,TrueNuE-muonRecoE);
609  f25->Fill(CCEnergy,TrueNuE);
610 
611  if (isQE){
612  f2->Fill(hadE,(TrueNuE-muonRecoE));
613  if (lessThan5GeV){
614  f6->Fill((recoHadQE - (TrueNuE-muonRecoE))/(TrueNuE-muonRecoE));
615  }
616  f7->Fill(recoHadQE,(recoHadQE - (TrueNuE-muonRecoE))/(TrueNuE-muonRecoE));
617  if (lessThan5GeV){
618  f10->Fill((QEEnergy-TrueNuE)/TrueNuE);
619  }
620  f11->Fill(QEEnergy,(QEEnergy-TrueNuE)/TrueNuE);
621  f15->Fill(recoHadQE,TrueNuE-muonRecoE);
622  f17->Fill(QEEnergy,TrueNuE);
623  }
624  else{
625  f3->Fill(hadE,(TrueNuE-muonRecoE));
626  if (lessThan5GeV){
627  f8->Fill((recoHadNonQE - (TrueNuE-muonRecoE))/(TrueNuE-muonRecoE));
628  }
629  f9->Fill(recoHadNonQE,(recoHadNonQE - (TrueNuE-muonRecoE))/(TrueNuE-muonRecoE));
630  if (lessThan5GeV){
631  f12->Fill((NonQEEnergy-TrueNuE)/TrueNuE);
632  }
633  f13->Fill(NonQEEnergy,(NonQEEnergy-TrueNuE)/TrueNuE);
634  f16->Fill(recoHadNonQE,TrueNuE-muonRecoE);
635  f18->Fill(NonQEEnergy,TrueNuE);
636  }
637  */
638  }//End of loop over things that pass all basic numu reco cuts
639  }//End of loop over things with tracks and energy
640  }//End of loop over entries
641 
642 
643  TCanvas* c1 = new TCanvas("c1","c1");
644  f1->Draw("colz");
645  c1->Print("muonTrkLenVTrueEnergy_ActiveOnly.pdf");
646 
647  TCanvas* c2 = new TCanvas("c2","c2");
648  f2->Draw("colz");
649  c2->Print("muonTrkLenVTrueEnergy_CatcherOnly.pdf");
650 
651  TCanvas* c3 = new TCanvas("c3","c3");
652  f3->Draw("colz");
653  c3->Print("muonTrkLenVTrueEnergy_CatcherForActiveCatcher.pdf");
654 
655  TCanvas* c6 = new TCanvas("c6","c6");
656  f6->Draw("colz");
657  c6->Print("QEHadronEVTrueEnergy.pdf");
658 
659  TCanvas* c7 = new TCanvas("c7","c7");
660  f7->Draw("colz");
661  c7->Print("NonQEHadronEVTrueEnergy.pdf");
662 
663  TCanvas* c8 = new TCanvas("c8","c8");
664  f8->Draw("colz");
665  c8->Print("CCHadronEVTrueEnergy.pdf");
666 
667  TCanvas* c9 = new TCanvas("c9","c9");
668  f9->Draw();
669  c9->Print("NDmuon1dRes_allAct.pdf");
670 
671  TCanvas* c10 = new TCanvas("c10","c10");
672  f10->Draw("colz");
673  c10->Print("NDmuon2dRes_allAct.pdf");
674 
675  TCanvas* c11 = new TCanvas("c11","c11");
676  f11->Draw("colz");
677  c11->Print("NDmuonReco_allAct.pdf");
678 
679  TCanvas* c12 = new TCanvas("c12","c12");
680  f12->Draw();
681  c12->Print("NDmuon1dRes_actCat_CatOnly.pdf");
682 
683  TCanvas* c13 = new TCanvas("c13","c13");
684  f13->Draw("colz");
685  c13->Print("NDmuon2dRes_actCat_CatOnly.pdf");
686 
687  TCanvas* c14 = new TCanvas("c14","c14");
688  f14->Draw("colz");
689  c14->Print("NDmuonReco_actCat_CatOnly.pdf");
690 
691  TCanvas* c15 = new TCanvas("c15","c15");
692  f15->Draw();
693  c15->Print("NDmuon1dRes_actCat_tog.pdf");
694 
695  TCanvas* c16 = new TCanvas("c16","c16");
696  f16->Draw("colz");
697  c16->Print("NDmuon2dRes_actCat_tog.pdf");
698 
699  TCanvas* c17 = new TCanvas("c17","c17");
700  f17->Draw("colz");
701  c17->Print("NDmuonReco_actCat_tog.pdf");
702 
703  TCanvas* c18 = new TCanvas("c18","c18");
704  f18->Draw();
705  c18->Print("NDmuon1dRes.pdf");
706 
707  TCanvas* c19 = new TCanvas("c19","c19");
708  f19->Draw("colz");
709  c19->Print("NDmuon2dRes.pdf");
710 
711  TCanvas* c20 = new TCanvas("c20","c20");
712  f20->Draw("colz");
713  c20->Print("NDmuonReco.pdf");
714 
715  TCanvas* c21 = new TCanvas("c21","c21");
716  f21->Draw();
717  c21->Print("NDhadqe1dRes.pdf");
718 
719  TCanvas* c22 = new TCanvas("c22","c22");
720  f22->Draw("colz");
721  c22->Print("NDhadqe2dRes.pdf");
722 
723  TCanvas* c23 = new TCanvas("c23","c23");
724  f23->Draw("colz");
725  c23->Print("NDhadqeReco.pdf");
726 
727  TCanvas* c24 = new TCanvas("c24","c24");
728  f24->Draw();
729  c24->Print("NDhadnonqe1dRes.pdf");
730 
731  TCanvas* c25 = new TCanvas("c25","c25");
732  f25->Draw("colz");
733  c25->Print("NDhadnonqe2dRes.pdf");
734 
735  TCanvas* c26 = new TCanvas("c26","c26");
736  f26->Draw("colz");
737  c26->Print("NDhadnonqeReco.pdf");
738 
739  TCanvas* c27 = new TCanvas("c27","c27");
740  f27->Draw();
741  c27->Print("NDhadcc1dRes.pdf");
742 
743  TCanvas* c28 = new TCanvas("c28","c28");
744  f28->Draw("colz");
745  c28->Print("NDhadcc2dRes.pdf");
746 
747  TCanvas* c29 = new TCanvas("c29","c29");
748  f29->Draw("colz");
749  c29->Print("NDhadccReco.pdf");
750 
751  TCanvas* c30 = new TCanvas("c30","c30");
752  f30->Draw();
753  c30->Print("NDneutrinoqe1dRes.pdf");
754 
755  TCanvas* c31 = new TCanvas("c31","c31");
756  f31->Draw("colz");
757  c31->Print("NDneutrinoqe2dRes.pdf");
758 
759  TCanvas* c32 = new TCanvas("c32","c32");
760  f32->Draw("colz");
761  c32->Print("NDneutrinoqeReco.pdf");
762 
763  TCanvas* c33 = new TCanvas("c33","c33");
764  f33->Draw();
765  c33->Print("NDneutrinononqe1dRes.pdf");
766 
767  TCanvas* c34 = new TCanvas("c34","c34");
768  f34->Draw("colz");
769  c34->Print("NDneutrinononqe2dRes.pdf");
770 
771  TCanvas* c35 = new TCanvas("c35","c35");
772  f35->Draw("colz");
773  c35->Print("NDneutrinononqeReco.pdf");
774 
775  TCanvas* c36 = new TCanvas("c36","c36");
776  f36->Draw();
777  c36->Print("NDneutrinocc1dRes.pdf");
778 
779  TCanvas* c37 = new TCanvas("c37","c37");
780  f37->Draw("colz");
781  c37->Print("NDneutrinocc2dRes.pdf");
782 
783  TCanvas* c38 = new TCanvas("c38","c38");
784  f38->Draw("colz");
785  c38->Print("NDneutrinoccReco.pdf");
786 
787  TFile f("2DPlotsForFitting.root","RECREATE");
788  f1->Write();
789  f2->Write();
790  f3->Write();
791  f4->Write();
792  f6->Write();
793  f7->Write();
794  f8->Write();
795  f9->Write();
796  f10->Write();
797  f11->Write();
798  f12->Write();
799  f13->Write();
800  f14->Write();
801  f15->Write();
802  f16->Write();
803  f17->Write();
804  f18->Write();
805  f19->Write();
806  f20->Write();
807  f21->Write();
808  f22->Write();
809  f23->Write();
810  f24->Write();
811  f25->Write();
812  f26->Write();
813  f27->Write();
814  f28->Write();
815  f29->Write();
816  f30->Write();
817  f31->Write();
818  f32->Write();
819  f33->Write();
820  f34->Write();
821  f35->Write();
822  f36->Write();
823  f37->Write();
824  f38->Write();
825 }
Float_t f4
Double_t MuonEFromTrackLength(double trkLen)
void abs(TH1 *hist)
void getPlotsFromNtupleND()
Float_t f2
c2
Definition: demo5.py:33
Float_t f3
double NDMuonEInActiveOnly(double actTrkLen)
double NDMuonEInActiveAndCatcherJustCatcherE(double actTrkLen, double catTrkLen)
Float_t f1
double NDMuonEFromTrackLength(double actTrkLen, double catTrkLen, double trkCalAct, double trkCalTran, double trkCalCat)
const Cut kNumuNCRej([](const caf::SRProxy *sr){return(sr->sel.remid.pid >0.75);})
Definition: NumuCuts.h:24
OStream cout
Definition: OStream.cxx:6
TFile * file
Definition: cellShifts.C:17
c1
Definition: demo5.py:24
double NDMuonEInActiveAndCatcher(double actTrkLen, double catTrkLen)
const Cut kNumuQuality
Definition: NumuCuts.h:18
double ActiveTrackLengthFromMuonE(double muonE)
double NDMuonEInActiveAndCatcherEffLen(double actTrkLen, double catTrkLen)