ComparisonPlots_MC.C
Go to the documentation of this file.
1 #include "CAFAna/Analysis/Prod3Loaders.h"
2 
3 #include "CAFAna/Cuts/Cuts.h"
4 #include "3FlavorAna/Cuts/NumuCuts.h"
5 #include "3FlavorAna/Cuts/NumuCuts2017.h"
6 #include "3FlavorAna/Cuts/QuantileCuts.h"
7 #include "CAFAna/Cuts/TimingCuts.h"
8 #include "CAFAna/Core/Spectrum.h"
9 #include "CAFAna/Core/SpectrumLoader.h"
10 #include "CAFAna/Core/Loaders.h"
11 #include "CAFAna/nus/NuSPlotFunctions.h"
12 
13 #include "CAFAna/Vars/GenieWeights.h"
14 #include "CAFAna/Vars/HistAxes.h"
15 #include "3FlavorAna/Vars/HistAxes.h"
16 #include "CAFAna/Vars/PPFXWeights.h"
17 
18 #include "TCanvas.h"
19 #include "TH2.h"
20 #include "TProfile.h"
21 
22 #include "OscLib/OscCalcPMNSOpt.h"
23 #include "CAFAna/Prediction/PredictionNoExtrap.h"
24 
25 #include "3FlavorAna/Ana2017/numu/Supporting_plots/DataMCComp_Comparison/ComparisonHints.h"
26 #include "3FlavorAna/Vars/NumuVars.h"
27 
28 using namespace ana;
29 
30 void ComparisonPlots_MC() {
31  std::cout << "Please provide;"
32  << "A Detector;\n\t1 For ND.\n\t2 for FD. \n"
33  << "An input samweb definition. \n"
34  << "An output file name. \n"
35  << std::endl;
36 }
37 
38 void ComparisonPlots_MC( int WhichDetector, std::string OutputName,
39  std::string InputDef = "full",
40  unsigned int NumCuts=1, bool Debug=false )
41 {
42  // The variables I am going to use for WhichDetector are;
43  // WhichDetector == 1 -- Near Detector
44  // WhichDetector == 2 -- Far Detector
45 
46  // --- Set my loader
47  ECAFType CAFType = ECAFType::kFullCAF;
48  if ( (InputDef.find("prod_sumdecaf_") != std::string::npos ) &&
49  (InputDef.find("numu2017") != std::string::npos ) ) {
50  CAFType = ECAFType::kNumuConcat;
51  }
52  Loaders *loader = new Prod3NomLoaders( CAFType, "full");
53  loader->SetSpillCut(kStandardSpillCuts);
54 
55  /* // To use Bruno's files.
56  // --- Define my loader, and pass it the file / dataset which it will load...
57  std::string fd_dir = "/pnfs/nova/persistent/users/bzamoran/CheckEvents/";
58  std::string fd_nonswap = fd_dir + "nonswap.root";
59  std::string fd_fluxswap = fd_dir + "fluxswap.root";
60  std::string fd_tau = fd_dir + "tau.root";
61  Loaders loader;
62  loader.SetSpillCut(kStandardSpillCuts);
63  loader.SetLoaderPath( fd_nonswap, caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap );
64  loader.SetLoaderPath( fd_fluxswap, caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kFluxSwap);
65  loader.SetLoaderPath( fd_tau, caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kTauSwap );
66  */
67  // --- Define my broad cuts - this is taken from Michael - Devs Repo /mbaird42/tracker_comparisons/FourSquare_plots.C
68  // Add the detector specific ones.
69  Cut kDetCut = kNoCut;
70  if ( WhichDetector == 1 ) {
71  kDetCut = kNumuContainND2017;
72  } else if (WhichDetector == 2) {
73  kDetCut = kNumuContainFD2017;
74  }
75 
76  // ----- Define my cuts....
77  // ----- Define my cuts....
78  std::vector<Cut> kDefinedCut;
79  std::string CutNames[NumCuts];
80  if (NumCuts == 4) { // --- NumCuts == 4
81  // --- Cut
82  kDefinedCut.push_back( kInBeamSpill && kNumuQuality ); // Good quality NuMu events
83  kDefinedCut.push_back( kInBeamSpill && kNumuQuality && kDetCut ); // That the event is contained.
84  kDefinedCut.push_back( kInBeamSpill && kNumuQuality && kDetCut ); // Reject cosmics
85  kDefinedCut.push_back( kInBeamSpill && kNumuQuality && kDetCut && kNumuPID2017 ); // Reject Neutral current events
86  // --- CutNames
87  CutNames[0] = "_QualCuts";
88  CutNames[1] = "_DetCuts" ;
89  CutNames[2] = "_CosmCuts";
90  CutNames[3] = "_NeuCCuts";
91  // --- CosmicRej
92  if ( WhichDetector == 2 ) {
93  kDefinedCut[2] = kDefinedCut[2] && kNumuCosmicRej2017;
94  kDefinedCut[3] = kDefinedCut[3] && kNumuCosmicRej2017;
95  }
96  }
97  else if (NumCuts == 1) { // --- NumCuts == 1
98  // --- Cut
99  kDefinedCut.push_back( kInBeamSpill && kNumuQuality && kDetCut && kNumuPID2017 ); // Reject Neutral current events
100  // --- CutNames
101  CutNames[0] = "_NeuCCuts";
102  // --- CosmicRej
103  if ( WhichDetector == 2 ) {
104  kDefinedCut[0] = kDefinedCut[0] && kNumuCosmicRej2017;
105  }
106  }
107 
108  // --- Load Michaels QuantileCuts.
109  std::string WhichPeriod = "full";
110  if (InputDef.find("period1") != std::string::npos ) WhichPeriod = "period1";
111  else if (InputDef.find("period2") != std::string::npos ) WhichPeriod = "period2";
112  else if (InputDef.find("period3") != std::string::npos ) WhichPeriod = "period3";
113  else if (InputDef.find("period5") != std::string::npos ) WhichPeriod = "period5";
114  TString fdspecfile = "/nova/app/users/karlwarb/Workspace/NuMuSystematics/FDQuantileHists/final_"+WhichPeriod+"_FD_histo_for_quantile_cuts.root";
115  TFile inFile(fdspecfile);
116  gDirectory->cd("dir_FDSpec2D");
117  TH2 *FDSpec2D = (TH2*)inFile.FindObjectAny("FDSpec2D");
118  // Add get the cuts....
119  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, 4);
120 
121  // ----- Define my quantiles.
122  const unsigned int NumQuants = 5;
123  Cut kQuantCut[NumQuants] = { kNoCut, HadEFracQuantCuts[0], HadEFracQuantCuts[1], HadEFracQuantCuts[2], HadEFracQuantCuts[3] };
124  const std::string QuantNames[NumQuants] = {
125  "_AllQuants",
126  "_Quant1",
127  "_Quant2",
128  "_Quant3",
129  "_Quant4"
130  };
131 
132  // ---- Write some output.
133  std::cout << "Before I declare anything, the following is set;\n"
134  << "\tWhichDetector is " << WhichDetector << "\n"
135  << "\tInput definit is " << InputDef << "\n"
136  << "\tOutput file name " << OutputName << "\n"
137  << "\tNumber of Cuts " << NumCuts << "\n"
138  << "\tIs debug mode is " << Debug << "\n"
139  << std::endl;
140 
141 
142  // --- Declare all of my spectra, before setting them in a loop...
143  // *** Some key ones first.
144  PredictionNoExtrap *sReconstEnergy[NumCuts][NumQuants]; // --- Reconstructed energy --- Vars.h ===> kCCE ---> sr->energy.numu.E
145  PredictionNoExtrap *sSliceTimeFull[NumCuts][NumQuants]; // --- Mean time of slice Full --- NumuVars.h ===> kSliceTime ---> sr->slc.meantime/1000
146  PredictionNoExtrap *sSliceTimeZoom[NumCuts][NumQuants]; // --- Mean time of slice Zoom --- NumuVars.h ===> kSliceTime ---> sr->slc.meantime/1000
147  PredictionNoExtrap *sMuonEnergy [NumCuts][NumQuants]; // --- Muon energy in the slice --- Vars.h ===> kMuE --> energy.numu.recomuonE
148  PredictionNoExtrap *sMuonEnPerHit [NumCuts][NumQuants]; // --- Muon energy per hit --- NumuVars.h ===> kTrkEPerNHit --> sr->trk.kalman.tracks[0].calE / sr->trk.kalman.tracks[0].nhit
149  PredictionNoExtrap *sHadronEnergy [NumCuts][NumQuants]; // --- Hadronic energy --- NumuVars.h ===> kHadE --> sr->energy.numu.E - sr->energy.numu.recomuonE
150  PredictionNoExtrap *sHadroEnPerHit[NumCuts][NumQuants]; // --- Hadronic energy per hit --- NumuVars.h ===> kHadEPerNHit --> function.
151  PredictionNoExtrap *sHadFracEnergy[NumCuts][NumQuants]; // --- Hadronic energy fraction --- NumuVars.h ===> kHadEFrac --> kHadE / kCCE
152  PredictionNoExtrap *sHitsPerSlice [NumCuts][NumQuants]; // --- Number of hits per slice --- Vars.h ===> kNHit --> slc.nhit
153  PredictionNoExtrap *sRemIDScore [NumCuts][NumQuants]; // --- The RemID score --- Vars.h ===> kRemID --> sr->sel.remid.pid
154  PredictionNoExtrap *sCVNCosmicScor[NumCuts][NumQuants]; // --- CVN Cosmic score --- Vars.h ===> sr->sel.cvn.output[14]
155  PredictionNoExtrap *sCVNNuMuIDScor[NumCuts][NumQuants]; // --- CVN NuMuID score --- Vars.h ===> sr->sel.cvn.numuid
156  PredictionNoExtrap *sNuMuContPID [NumCuts][NumQuants]; // --- Third analysis BDT --- NuMuVars.h ===> sr->sel.cosrej.numucontpid
157  PredictionNoExtrap *sSANuMuContPID[NumCuts][NumQuants]; // --- Second analysis BDT --- Above.
158  // *** Some important Truth based ones.
159  PredictionNoExtrap *sTrueNuEnergy [NumCuts][NumQuants]; // --- True neutrino energy --- Vars.h ===> kTrueE ---> sr->mc.nu[0].E
160  PredictionNoExtrap *sReTrOverTrNuE[NumCuts][NumQuants]; // --- Reco - True / True Nu E --- Above.
161 
162  // *** Basic positional plots.
163  PredictionNoExtrap *sTrkStartX[NumCuts][NumQuants]; // --- Start Pos of Track --- NumuVars.h ===> kTrkStartX --> trk.kalman.tracks[0].start.X()/100
164  PredictionNoExtrap *sTrkStartY[NumCuts][NumQuants]; // --- Start Pos of Track --- NumuVars.h ===> kTrkStartY --> trk.kalman.tracks[0].start.Y()/100
165  PredictionNoExtrap *sTrkStartZ[NumCuts][NumQuants]; // --- Start Pos of Track --- NumuVars.h ===> kTrkStartZ --> trk.kalman.tracks[0].start.Z()/100
166  PredictionNoExtrap *sTrkEndX [NumCuts][NumQuants]; // --- End Pos of Track --- NumuVars.h ===> kTrkEndX --> trk.kalman.tracks[0].stop.X()/100
167  PredictionNoExtrap *sTrkEndY [NumCuts][NumQuants]; // --- End Pos of Track --- NumuVars.h ===> kTrkEndY --> trk.kalman.tracks[0].stop.Y()/100
168  PredictionNoExtrap *sTrkEndZ [NumCuts][NumQuants]; // --- End Pos of Track --- NumuVars.h ===> kTrkEndZ --> trk.kalman.tracks[0].stop.Z()/100
169  PredictionNoExtrap *sTrkLenXY [NumCuts][NumQuants]; // --- Diff in Track Len X Y --- Above. ===> kTrkLenXY --> (tracks[0].start.X - tracks[0].stop.X) / 100
170  //*
171  // *** Track based ones...
172  PredictionNoExtrap *sNumKalTracks [NumCuts][NumQuants]; // --- Number of Kal trk in slc --- NuMuVars.h ===> sr->trk.kalman.ntracks
173  PredictionNoExtrap *sNumKalTrHits [NumCuts][NumQuants]; // --- Number of hits in Kal tr --- NuMuVars.h ===> sr->trk.kalman.tracks[0].nhit
174  PredictionNoExtrap *sRatKalHitSlc [NumCuts][NumQuants]; // --- Ratio of Kal hits in slc --- Above.
175  PredictionNoExtrap *sKalTrLength [NumCuts][NumQuants]; // --- Length of Kalman track --- NuMuVars.h ===> sr->trk.kalman.tracks[0].len / 100
176  PredictionNoExtrap *sKalTrBeamAng [NumCuts][NumQuants]; // --- Angle Kal track & beam --- NuMuVars.h ===> kCosNumi --> sr->trk.kalman.tracks[0].dir.Dot(beamDirND)
177  PredictionNoExtrap *sKalMostFwdCel[NumCuts][NumQuants]; // --- Most forward Kal tr hit --- Above.
178  PredictionNoExtrap *sKalMostBakCel[NumCuts][NumQuants]; // --- Most backward Kal tr hit --- Above.
179  PredictionNoExtrap *sKalTrVer_MaxY[NumCuts][NumQuants]; // --- Highest Y Kalman tr vert --- Above.
180  PredictionNoExtrap *sKalTrVer_MaxZ[NumCuts][NumQuants]; // --- Highest Z Kalman tr vert --- Above.
181  PredictionNoExtrap *sKalTrDir_Y [NumCuts][NumQuants]; // --- Y direction of Kal track --- NuMuVars.h ===> sr->trk.kalman.tracks[0].dir.Y()
182  PredictionNoExtrap *sScattKalTrLen[NumCuts][NumQuants]; // --- Scattering over tr len --- Above.
183 
184  // *** Hit based ones...
185  PredictionNoExtrap *sFirstHitCell [NumCuts][NumQuants]; // --- The first cell hit --- Above.
186  PredictionNoExtrap *sLastHitCell [NumCuts][NumQuants]; // --- The last cell hit --- Above.
187  PredictionNoExtrap *sMaxActivity_Y[NumCuts][NumQuants]; // --- Maximum Y pos of activty --- NuMuVars.h ===> sr->slc.boxmax.Y() / 100.
188  PredictionNoExtrap *sMinActivity_Y[NumCuts][NumQuants]; // --- Minimum Y pos of activty --- NuMuVars.h ===> sr->slc.boxmin.Y() / 100.
189  PredictionNoExtrap *sMaxActivity_Z[NumCuts][NumQuants]; // --- Maximum Z pos of activty --- NuMuVars.h ===> sr->slc.boxmax.Z() / 100.
190  PredictionNoExtrap *sMinActivity_Z[NumCuts][NumQuants]; // --- Minimum Z pos of activty --- NuMuVars.h ===> sr->slc.boxmin.Z() / 100.
191  PredictionNoExtrap *sMinCellToEdge[NumCuts][NumQuants]; // --- Minimum cells to edge --- Above.
192 
193  // *** Visible energy quanta.
194  // Slice
195  PredictionNoExtrap *sVisEInSlice [NumCuts][NumQuants]; // --- The Vis En in the slice --- Above. ---> kVisESlc ---> sr->slc.calE
196  PredictionNoExtrap *sNHitInSlice [NumCuts][NumQuants]; // --- The Num.Hits in the slic --- Above. ---> kNHitSlc ---> sr->slc.nhit
197  PredictionNoExtrap *sVisEPerHitSlc[NumCuts][NumQuants]; // --- Visibile E per hit in sl --- Above. ---> kEnPHSlc ---> kVisESlc / (1.78 * kNHitSlc)
198  // Track
199  PredictionNoExtrap *sVisEInTrack [NumCuts][NumQuants]; // --- The Vis En of the track --- Above. ---> kVisETrk ---> sr->trk.kalman.tracks[0].calE
200  PredictionNoExtrap *sNHitInTrack [NumCuts][NumQuants]; // --- The NumHits in the track --- Above. ---> kNHitTrk ---> sr->trk.kalman.tracks[0].nhit
201  PredictionNoExtrap *sVisEPerHitTrk[NumCuts][NumQuants]; // --- Visibile E PerHit in trk --- Above. ---> kEnPHTrk ---> kVisESlc / (1.78 * kNHitTrk)
202  // Hadronic
203  PredictionNoExtrap *sVisEInHadron [NumCuts][NumQuants]; // --- The Hadronic Vis Energy --- Above. ---> kVisEHad ---> (kVisESlc - kVisETrk)
204  PredictionNoExtrap *sNHitInHadron [NumCuts][NumQuants]; // --- The Num of hadronic Hits --- Above. ---> kNHitHad ---> (kNHitHad - kNHitTrk)
205  PredictionNoExtrap *sVisEPerHitHad[NumCuts][NumQuants]; // --- Visibile E PerHit hadron --- Above. ---> kEnPHHad ---> kVisEHad / (1.78 * kNHitHad)
206 
207  // *** RemID score / inputs
208  PredictionNoExtrap *sRemIDScatLLH [NumCuts][NumQuants]; // --- RemID scattering angle --- NuMuVars.h ===> kReMIdScatLLH --> sel.remid.scatllh
209  PredictionNoExtrap *sRemIDdEdxLLH [NumCuts][NumQuants]; // --- RemID dE/dx --- NuMuVars.h ===> kReMIdDEDxLLH --> sel.remid.dedxllh
210  PredictionNoExtrap *sRemIDMeasFrac[NumCuts][NumQuants]; // --- RemID fraction of planes --- NuMuVars.h ===> kReMIdMeasFrac -> sel.remid.measfrac
211 
212  // *** Some calibration quantities
213  PredictionNoExtrap *sCor_RecEn[NumCuts][NumQuants][4];
214  PredictionNoExtrap *sCor_MuoEn[NumCuts][NumQuants][4];
215  PredictionNoExtrap *sCor_HadEn[NumCuts][NumQuants][4];
216  PredictionNoExtrap *sCor_HFrEn[NumCuts][NumQuants][4];
217  // Some special truth based ones.
218  PredictionNoExtrap *sCor_RecTr[NumCuts][NumQuants][4];
219  PredictionNoExtrap *sCor_MuoTr[NumCuts][NumQuants][4];
220  PredictionNoExtrap *sCor_HadTr[NumCuts][NumQuants][4];
221  PredictionNoExtrap *sCor_HFrTr[NumCuts][NumQuants][4];
222  PredictionNoExtrap *sCor_ReTr [NumCuts][NumQuants][4];
223  //*/
224  const Cut kCorner[4] = { kW_TL, kW_TR, kW_BL, kW_BR };
225  const std::string CutName[4] = { "_TopL", "_TopR", "_BotL", "_BotR" };
226 
227  // --- And now set what is going to go into each of my spectra.
228  for (unsigned int cc=0; cc<NumCuts; ++cc) {
229  // --- Loop through my quantiles.
230  for (unsigned int qq=0; qq<NumQuants; ++qq) {
231  // What is the cut this time?
232  const Cut kThisCut = kDefinedCut[cc] && kQuantCut[qq];
233 
234  // *** Some key ones first...
235  sReconstEnergy[cc][qq] = new PredictionNoExtrap( *loader, "The reconstructed energy", EnergyBins, kCCE , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
236  sSliceTimeFull[cc][qq] = new PredictionNoExtrap( *loader, "The time of slice - Full", TimeFBins , kSliceTime , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
237  sSliceTimeZoom[cc][qq] = new PredictionNoExtrap( *loader, "The time of slice - Zoom", TimeZBins , kSliceTime , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
238  sMuonEnergy [cc][qq] = new PredictionNoExtrap( *loader, "The muon energy " , EnergyBins, kMuE , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
239  sMuonEnPerHit [cc][qq] = new PredictionNoExtrap( *loader, "The muon energy per hit" , EnPHitBins, kTrkEPerNHit , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
240  sHadronEnergy [cc][qq] = new PredictionNoExtrap( *loader, "The hadronic energy" , EnergyBins, kHadE , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
241  sHadroEnPerHit[cc][qq] = new PredictionNoExtrap( *loader, "The hadronic En per hit ", EnPHitBins, kHadEPerNHit , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
242  sHadFracEnergy[cc][qq] = new PredictionNoExtrap( *loader, "The hadronic energy frac", RatioBins , kHadEFrac , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
243  sHitsPerSlice [cc][qq] = new PredictionNoExtrap( *loader, "Number of hits per slice", NumHitBins, kNHit , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
244  sRemIDScore [cc][qq] = new PredictionNoExtrap( *loader, "The RemID score" , RatioBins , kRemID , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
245  sCVNCosmicScor[cc][qq] = new PredictionNoExtrap( *loader, "CVN Cosmic score" , RatioBins , kCVNcos2017 , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
246  sCVNNuMuIDScor[cc][qq] = new PredictionNoExtrap( *loader, "CVN NuMu score" , RatioBins , kCVNm , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
247  sNuMuContPID [cc][qq] = new PredictionNoExtrap( *loader, "The NuMu containment PID", RatioBins , kNumuContPID , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
248  sSANuMuContPID[cc][qq] = new PredictionNoExtrap( *loader, "SA NuMu containment PID" , RatioBins , kSANumuContPID, kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
249  // *** Some important truth based ones.
250  sTrueNuEnergy [cc][qq] = new PredictionNoExtrap( *loader, "The true neutrino energy", EnergyBins, kTrueE , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
251  sReTrOverTrNuE[cc][qq] = new PredictionNoExtrap( *loader, "The R-T ov T Nu energy" , RangeBins , kRTOvTNuE, kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
252 
253  // *** Positional plots
254  sTrkStartX[cc][qq] = new PredictionNoExtrap( *loader, "Starting X Pos of Track", YRangeBins, kTrkStartX, kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
255  sTrkStartY[cc][qq] = new PredictionNoExtrap( *loader, "Starting Y Pos of Track", YRangeBins, kTrkStartY, kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
256  sTrkStartZ[cc][qq] = new PredictionNoExtrap( *loader, "Starting Z Pos of Track", ZRangeBins, kTrkStartZ, kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
257  sTrkEndX[cc][qq] = new PredictionNoExtrap( *loader, "Ending X Pos of Track" , YRangeBins, kTrkEndX , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
258  sTrkEndY[cc][qq] = new PredictionNoExtrap( *loader, "Ending Y Pos of Track" , YRangeBins, kTrkEndY , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
259  sTrkEndZ[cc][qq] = new PredictionNoExtrap( *loader, "Ending Z Pos of Track" , ZRangeBins, kTrkEndZ , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
260  sTrkLenXY[cc][qq] = new PredictionNoExtrap( *loader, "Difference in XY Tr Len", ZRangeBins, kTrkLenXY , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
261  //*
262  // *** Track based ones...
263  sNumKalTracks [cc][qq] = new PredictionNoExtrap( *loader, "Number of Kalman Tracks" , NumTrkBins, kNKalman , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
264  sNumKalTrHits [cc][qq] = new PredictionNoExtrap( *loader, "Number of Kal Track hits", NumHitBins, kTrkNhits , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
265  sRatKalHitSlc [cc][qq] = new PredictionNoExtrap( *loader, "Ratio of KalHit in Slice", RatioBins , kRatOfKalHi, kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
266  sKalTrLength [cc][qq] = new PredictionNoExtrap( *loader, "Length of Kalman Track" , LengthBins, kTrkLength , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
267  sKalTrBeamAng [cc][qq] = new PredictionNoExtrap( *loader, "Angle of KalTr to Beam" , RatioBins , kCosNumi , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
268  sKalMostFwdCel[cc][qq] = new PredictionNoExtrap( *loader, "Most forward KalTr Cell" , CellsBins , kKalFwdCell, kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
269  sKalMostBakCel[cc][qq] = new PredictionNoExtrap( *loader, "Most backward KalTr Cell", CellsBins , kKalBakCell, kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
270  sKalTrVer_MaxY[cc][qq] = new PredictionNoExtrap( *loader, "Most pos Y KalTr Vertex" , YRangeBins, kMaxKalYPos, kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
271  sKalTrVer_MaxZ[cc][qq] = new PredictionNoExtrap( *loader, "Most pos Z KalTr Vertex" , ZRangeBins, kMaxKalZPos, kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
272  sKalTrDir_Y [cc][qq] = new PredictionNoExtrap( *loader, "Y Direction of KalTr" , CompBins , kDirY , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
273  sScattKalTrLen[cc][qq] = new PredictionNoExtrap( *loader, "Scatt over KalTr length" , TinyBins , kScattTrLen, kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
274 
275  // *** Hit based ones...
276  sFirstHitCell [cc][qq] = new PredictionNoExtrap( *loader, "The first cell hit in sl", CellsBins , kFirstCell , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
277  sLastHitCell [cc][qq] = new PredictionNoExtrap( *loader, "The last cell hit in sl" , CellsBins , kLastCell , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
278  sMaxActivity_Y[cc][qq] = new PredictionNoExtrap( *loader, "Largest Y loc with a hit", YRangeBins, kSlcMaxY , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
279  sMinActivity_Y[cc][qq] = new PredictionNoExtrap( *loader, "Smalest Y loc with a hit", YRangeBins, kSlcMinY , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
280  sMaxActivity_Z[cc][qq] = new PredictionNoExtrap( *loader, "Largest Z loc with a hit", ZRangeBins, kSlcMaxZ , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
281  sMinActivity_Z[cc][qq] = new PredictionNoExtrap( *loader, "Largest Z loc with a hit", ZRangeBins, kSlcMinZ , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
282  sMinCellToEdge[cc][qq] = new PredictionNoExtrap( *loader, "Min cells hit to edge" , CellsBins , kMinCellEdg, kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
283 
284  // *** Visible energy quanta.
285  // Slice
286  sVisEInSlice [cc][qq] = new PredictionNoExtrap( *loader, "Visible Energy in slice" , EnergyBins, kVisESlc, kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
287  sNHitInSlice [cc][qq] = new PredictionNoExtrap( *loader, "Number of hits in slice" , NumHitBins, kNHitSlc, kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
288  sVisEPerHitSlc[cc][qq] = new PredictionNoExtrap( *loader, "Visible E per hit in slc", EnPHitBins, kEnPHSlc, kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
289  // Track
290  sVisEInTrack [cc][qq] = new PredictionNoExtrap( *loader, "Visible Energy in track" , EnergyBins, kVisETrk, kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
291  sNHitInTrack [cc][qq] = new PredictionNoExtrap( *loader, "Number of hits in track" , NumHitBins, kNHitTrk, kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
292  sVisEPerHitTrk[cc][qq] = new PredictionNoExtrap( *loader, "Visible E per hit in trk", EnPHitBins, kEnPHTrk, kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
293  // Hadronic
294  sVisEInHadron [cc][qq] = new PredictionNoExtrap( *loader, "Visible Hadronic Energy" , EnergyBins, kVisEHad, kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
295  sNHitInHadron [cc][qq] = new PredictionNoExtrap( *loader, "Number of hadronic hits" , NumHitBins, kNHitHad, kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
296  sVisEPerHitHad[cc][qq] = new PredictionNoExtrap( *loader, "Visible E per hadron hit", EnPHitBins, kEnPHHad, kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
297 
298  // *** RemID score / inputs
299  sRemIDScatLLH [cc][qq] = new PredictionNoExtrap( *loader, "RemID Scattering angle" , CompZBins , kReMIdScatLLH , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
300  sRemIDdEdxLLH [cc][qq] = new PredictionNoExtrap( *loader, "RemID dE/dx" , CompBins , kReMIdDEDxLLH , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
301  sRemIDMeasFrac[cc][qq] = new PredictionNoExtrap( *loader, "RemID fraction of planes", RatioBins , kReMIdMeasFrac , kThisCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017 );
302 
303  // *** Calibration quantities
304  for (int cor=0; cor<4; ++cor) {
305  const Cut kCorCut = kThisCut && kCorner[cor];
306  sCor_RecEn[cc][qq][cor] = new PredictionNoExtrap( *loader, "Reconstructed Energy with corner cut", EnergyBins, kCCE , kCorCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017);
307  sCor_MuoEn[cc][qq][cor] = new PredictionNoExtrap( *loader, "Muon energy with corner cut" , EnergyBins, kMuE , kCorCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017);
308  sCor_HadEn[cc][qq][cor] = new PredictionNoExtrap( *loader, "Hadronic energy with corner cut" , EnergyBins, kHadE , kCorCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017);
309  sCor_HFrEn[cc][qq][cor] = new PredictionNoExtrap( *loader, "Hadronic energy frac with corner cut", RatioBins , kHadEFrac, kCorCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017);
310  // Additional truth based ones.
311  sCor_RecTr[cc][qq][cor] = new PredictionNoExtrap( *loader, "True Energy with corner cut" , EnergyBins, kTrueE , kCorCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017);
312  sCor_MuoTr[cc][qq][cor] = new PredictionNoExtrap( *loader, "True Mu En with corner cut" , EnergyBins, kTrueMuonE, kCorCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017);
313  sCor_HadTr[cc][qq][cor] = new PredictionNoExtrap( *loader, "True Had En with corner cut" , EnergyBins, kTrueHadE , kCorCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017);
314  sCor_HFrTr[cc][qq][cor] = new PredictionNoExtrap( *loader, "True HadEFrac with corner cut", RatioBins , kTrueHFrE , kCorCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017);
315  sCor_ReTr [cc][qq][cor] = new PredictionNoExtrap( *loader, "Reco-True/True w. corner cut" , RangeBins , kRTOvTNuE , kCorCut, kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2017);
316  }
317  //*/
318  } // Quantile loop
319  } // Cut loop
320 
321  // --- Do it!
322  loader->Go();
323 
324  // --- Open my outfile so that I save histograms!
325  TFile *outFile = new TFile(OutputName.c_str(),"RECREATE");
326  std::cerr << "Made my output file...." << std::endl;
327 
328  osc::OscCalcPMNSOpt calc;
329  ana::ResetOscCalcToDefault(&calc);
330  calc.SetDmsq32(2.43071e-3);
331  calc.SetTh23(asin(sqrt(0.471974)));
332  calc.SetdCP(M_PI*0.915869); // preliminary Ana2017
333 
334  // --- Finally, write them to disk!
335  for (unsigned int cc=0; cc<NumCuts; ++cc) {
336  // --- Loop through my quantiles.
337  for (unsigned int qq=0; qq<NumQuants; ++qq) {
338  std::string MyAppend = CutNames[cc] + QuantNames[qq];
339  // *** Some key ones first...
340  sReconstEnergy[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sReconstEnergy")+TString(MyAppend) ) ;
341  sSliceTimeFull[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sSliceTimeFull")+TString(MyAppend) ) ;
342  sSliceTimeZoom[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sSliceTimeZoom")+TString(MyAppend) ) ;
343  sMuonEnergy [cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sMuonEnergy") +TString(MyAppend) ) ;
344  sMuonEnPerHit [cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sMuonEnPerHit") +TString(MyAppend) ) ;
345  sHadronEnergy [cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sHadronEnergy") +TString(MyAppend) ) ;
346  sHadroEnPerHit[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sHadroEnPerHit")+TString(MyAppend) ) ;
347  sHadFracEnergy[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sHadFracEnergy")+TString(MyAppend) ) ;
348  sHitsPerSlice [cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sHitsPerSlice") +TString(MyAppend) ) ;
349  sRemIDScore [cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sRemIDScore") +TString(MyAppend) ) ;
350  sCVNCosmicScor[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sCVNCosmicScor")+TString(MyAppend) ) ;
351  sCVNNuMuIDScor[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sCVNNuMuIDScor")+TString(MyAppend) ) ;
352  sNuMuContPID [cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sNuMuContPID") +TString(MyAppend) ) ;
353  sSANuMuContPID[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sSANuMuContPID")+TString(MyAppend) ) ;
354  // *** Some important truth based ones.
355  sTrueNuEnergy [cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sTrueNuEnergy") +TString(MyAppend) ) ;
356  sReTrOverTrNuE[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sReTrOverTrNuE")+TString(MyAppend) ) ;
357 
358  // *** Positional plots
359  sTrkStartX[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sTrkStartX")+TString(MyAppend) ) ;
360  sTrkStartY[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sTrkStartY")+TString(MyAppend) ) ;
361  sTrkStartZ[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sTrkStartZ")+TString(MyAppend) ) ;
362  sTrkEndX[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sTrkEndX") +TString(MyAppend) ) ;
363  sTrkEndY[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sTrkEndY") +TString(MyAppend) ) ;
364  sTrkEndZ[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sTrkEndZ") +TString(MyAppend) ) ;
365  sTrkLenXY[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sTrkLenXY") +TString(MyAppend) ) ;
366  //*
367  // *** Track based ones...
368  sNumKalTracks [cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sNumKalTracks") +TString(MyAppend) ) ;
369  sNumKalTrHits [cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sNumKalTrHits") +TString(MyAppend) ) ;
370  sRatKalHitSlc [cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sRatKalHitSlc") +TString(MyAppend) ) ;
371  sKalTrLength [cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sKalTrLength") +TString(MyAppend) ) ;
372  sKalTrBeamAng [cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sKalTrBeamAng") +TString(MyAppend) ) ;
373  sKalMostFwdCel[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sKalMostFwdCel")+TString(MyAppend) ) ;
374  sKalMostBakCel[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sKalMostBakCel")+TString(MyAppend) ) ;
375  sKalTrVer_MaxY[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sKalTrVer_MaxY")+TString(MyAppend) ) ;
376  sKalTrVer_MaxZ[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sKalTrVer_MaxZ")+TString(MyAppend) ) ;
377  sKalTrDir_Y [cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sKalTrDir_Y") +TString(MyAppend) ) ;
378  sScattKalTrLen[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sScattKalTrLen")+TString(MyAppend) ) ;
379 
380  // *** Hit based ones...
381  sFirstHitCell [cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sFirstHitCell") +TString(MyAppend) ) ;
382  sLastHitCell [cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sLastHitCell") +TString(MyAppend) ) ;
383  sMaxActivity_Y[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sMaxActivity_Y")+TString(MyAppend) ) ;
384  sMinActivity_Y[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sMinActivity_Y")+TString(MyAppend) ) ;
385  sMinCellToEdge[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sMinCellToEdge")+TString(MyAppend) ) ;
386  sMaxActivity_Z[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sMaxActivity_Z")+TString(MyAppend) ) ;
387  sMinActivity_Z[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sMinActivity_Z")+TString(MyAppend) ) ;
388 
389  // *** Visible energy quanta.
390  // Slice
391  sVisEInSlice [cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sVisEInSlice") +TString(MyAppend) ) ;
392  sNHitInSlice [cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sNHitInSlice") +TString(MyAppend) ) ;
393  sVisEPerHitSlc[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sVisEPerHitSlc")+TString(MyAppend) ) ;
394  // Track
395  sVisEInTrack [cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sVisEInTrack") +TString(MyAppend) ) ;
396  sNHitInTrack [cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sNHitInTrack") +TString(MyAppend) ) ;
397  sVisEPerHitTrk[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sVisEPerHitTrk")+TString(MyAppend) ) ;
398  // Hadronic
399  sVisEInHadron [cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sVisEInHadron") +TString(MyAppend) ) ;
400  sNHitInHadron [cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sNHitInHadron") +TString(MyAppend) ) ;
401  sVisEPerHitHad[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sVisEPerHitHad")+TString(MyAppend) ) ;
402 
403  // *** RemID score / inputs
404  sRemIDScatLLH [cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sRemIDScatLLH") +TString(MyAppend) ) ;
405  sRemIDdEdxLLH [cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sRemIDdEdxLLH") +TString(MyAppend) ) ;
406  sRemIDMeasFrac[cc][qq] -> Predict(&calc).SaveTo( outFile, TString("sRemIDMeasFrac")+TString(MyAppend) ) ;
407 
408  // *** Calibration quantities
409  for (int cor=0; cor<4; ++cor) {
410  sCor_RecEn[cc][qq][cor] -> Predict(&calc).SaveTo( outFile, TString("sCor_RecEn")+TString(MyAppend)+TString(CutName[cor]) ) ;
411  sCor_MuoEn[cc][qq][cor] -> Predict(&calc).SaveTo( outFile, TString("sCor_MuoEn")+TString(MyAppend)+TString(CutName[cor]) ) ;
412  sCor_HadEn[cc][qq][cor] -> Predict(&calc).SaveTo( outFile, TString("sCor_HadEn")+TString(MyAppend)+TString(CutName[cor]) ) ;
413  sCor_HFrEn[cc][qq][cor] -> Predict(&calc).SaveTo( outFile, TString("sCor_HFrEn")+TString(MyAppend)+TString(CutName[cor]) ) ;
414  // Additional truth based ones.
415  sCor_RecTr[cc][qq][cor] -> Predict(&calc).SaveTo( outFile, TString("sCor_RecTr")+TString(MyAppend)+TString(CutName[cor]) ) ;
416  sCor_MuoTr[cc][qq][cor] -> Predict(&calc).SaveTo( outFile, TString("sCor_MuoTr")+TString(MyAppend)+TString(CutName[cor]) ) ;
417  sCor_HadTr[cc][qq][cor] -> Predict(&calc).SaveTo( outFile, TString("sCor_HadTr")+TString(MyAppend)+TString(CutName[cor]) ) ;
418  sCor_HFrTr[cc][qq][cor] -> Predict(&calc).SaveTo( outFile, TString("sCor_HFrTr")+TString(MyAppend)+TString(CutName[cor]) ) ;
419  sCor_ReTr [cc][qq][cor] -> Predict(&calc).SaveTo( outFile, TString("sCor_ReTr" )+TString(MyAppend)+TString(CutName[cor]) ) ;
420  }
421  //*/
422  std::cerr << "\tWritten quantile set " << qq << " for cut set " << cc << " to disk!!!" << std::endl;
423  }
424  std::cerr << "\n\tWritten cut set " << cc << " to disk!!!" << std::endl;
425  }
426 
427 }
ana::kHadE
const Var kHadE
Definition: NumuVars.h:23
ana::Prod3NomLoaders
For nominal spectra and reweighting systs (xsec/flux)
Definition: Prod3Loaders.h:40
ana::ZRangeBins
const Binning ZRangeBins
Definition: ComparisonHints.h:36
ana::kNKalman
const Var kNKalman
Definition: NumuVars.cxx:540
ana
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
ana::kEnPHHad
const Var kEnPHHad([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return(float)(sr->slc.calE-sr->trk.kalman.tracks[0].calE)/(float)(1.78 *(sr->slc.nhit-sr->trk.kalman.tracks[0].nhit));})
ana::kMaxKalYPos
const Var kMaxKalYPos([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return std::max(sr->trk.kalman.tracks[0].start.Y()/100, sr->trk.kalman.tracks[0].stop.Y()/100);})
ana::ResetOscCalcToDefault
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
ana::kReMIdScatLLH
const Var kReMIdScatLLH
Definition: NumuVars.cxx:555
ana::kVisESlc
const Var kVisESlc([](const caf::SRProxy *sr){return sr->slc.calE;})
ana::kNHitSlc
const Var kNHitSlc([](const caf::SRProxy *sr){return sr->slc.nhit;})
ana::kPPFXFluxCVWgt
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
ana::kSlcMaxZ
const Var kSlcMaxZ([](const caf::SRProxy *sr){return sr->slc.boxmax.Z()/100.;})
Definition: NumuVars.h:93
ana::kSlcMaxY
const Var kSlcMaxY([](const caf::SRProxy *sr){return sr->slc.boxmax.Y()/100.;})
Definition: NumuVars.h:92
ana::kNumuContPID
const Var kNumuContPID
Definition: NumuVars.cxx:553
ana::kNumuContainND2017
const Cut kNumuContainND2017([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;for(unsigned int i=0;i< sr->vtx.elastic.fuzzyk.nshwlid;++i){TVector3 start=sr->vtx.elastic.fuzzyk.png[i].shwlid.start;TVector3 stop=sr->vtx.elastic.fuzzyk.png[i].shwlid.stop;if(std::min(start.X(), stop.X())< -180.0) return false;if(std::max(start.X(), stop.X()) > 180.0) return false;if(std::min(start.Y(), stop.Y())< -180.0) return false;if(std::max(start.Y(), stop.Y()) > 180.0) return false;if(std::min(start.Z(), stop.Z())< 20.0) return false;if(std::max(start.Z(), stop.Z()) > 1525.0) return false;}if(sr->trk.kalman.ntracks< 1) return false;for(unsigned int i=0;i< sr->trk.kalman.ntracks;++i){if(i==sr->trk.kalman.idxremid) continue;else if(sr->trk.kalman.tracks[i].start.Z() > 1275||sr->trk.kalman.tracks[i].stop.Z() > 1275) return false;}return(sr->trk.kalman.ntracks > sr->trk.kalman.idxremid &&sr->slc.firstplane > 1 &&sr->slc.lastplane< 212 &&sr->trk.kalman.tracks[0].start.Z()< 1100 &&(sr->trk.kalman.tracks[0].stop.Z()< 1275 ||sr->sel.contain.kalyposattrans< 55) &&sr->sel.contain.kalfwdcellnd > 5 &&sr->sel.contain.kalbakcellnd > 10);})
Definition: NumuCuts2017.h:11
ana::kEnPHSlc
const Var kEnPHSlc([](const caf::SRProxy *sr){if(sr->slc.nhit< 1) return-10.0f;return(float)(sr->slc.calE)/(float)(1.78 *sr->slc.nhit);})
ana::Loaders
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
ana::YRangeBins
const Binning YRangeBins
Definition: ComparisonHints.h:35
ana::kHadEFracAxis
const HistAxis kHadEFracAxis("E_{had.} / E_{#nu}", Binning::Simple(200, 0, 1), kHadEFrac)
HistAxis that implements Hadronic Energy fraction binning used by L.Vinton to derive Hadronic Energy ...
Definition: HistAxes.h:30
ana::EnPHitBins
const Binning EnPHitBins
Definition: ComparisonHints.h:29
ana::kNHitHad
const Var kNHitHad([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return(float)(sr->slc.nhit-sr->trk.kalman.tracks[0].nhit);})
ana::kMinCellEdg
const Var kMinCellEdg([](const caf::SRProxy *sr){return std::min((sr->sel.contain.kalfwdcell+sr->sel.contain.kalbakcell), (sr->sel.contain.cosfwdcell+sr->sel.contain.cosbakcell));})
ana::kTrkStartY
const Var kTrkStartY([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].start.Y()/100;})
Definition: NumuVars.h:52
ana::kFirstCell
const Var kFirstCell([](const caf::SRProxy *sr){return sr->slc.firstcell;})
ana::kSliceTime
const Var kSliceTime([](const caf::SRProxy *sr){return sr->slc.meantime/1000;})
Definition: NumuVars.h:34
om::cerr
OStream cerr
Definition: OStream.cxx:7
ana::kTrueHFrE
const Var kTrueHFrE([](const caf::SRProxy *sr){float HadE=kTrueHadE(sr);float TruE=kTrueE(sr);return(HadE/TruE);})
ana::kLastCell
const Var kLastCell([](const caf::SRProxy *sr){return sr->slc.lastcell;})
ana::TimeZBins
const Binning TimeZBins
Definition: ComparisonHints.h:41
kTrkEPerNHit
const Var kTrkEPerNHit([](const caf::SRProxy *sr){return(sr->trk.kalman.tracks[0].calE/sr->trk.kalman.tracks[0].nhit);})
ana::kReMIdMeasFrac
const Var kReMIdMeasFrac
Definition: NumuVars.cxx:557
ana::kKalFwdCell
const Var kKalFwdCell([](const caf::SRProxy *sr){return sr->sel.contain.kalfwdcell;})
calc
osc::OscCalcDumb calc
Definition: make_nueFDprediction_kinematics.C:70
ana::kTrkLength
const Var kTrkLength([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-5.f;return sr->trk.kalman.tracks[0].len/100;})
Definition: NumuVars.h:65
ana::kScattTrLen
const Var kScattTrLen([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return(float)((float) sr->sel.cosrej.scatt/(float) sr->trk.kalman.tracks[0].len);})
ana::kHadEPerNHit
const Var kHadEPerNHit([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return 0.0f;int nHit=sr->slc.nhit-sr->trk.kalman.tracks[0].nhit;if(nHit<=0) return 0.0f;float hadE=sr->energy.numu.hadcalE;return hadE/nHit;})
Definition: NumuVars.h:63
M_PI
#define M_PI
Definition: SbMath.h:34
ana::kDirY
const Var kDirY([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-5.f;return sr->trk.kalman.tracks[0].dir.Y();})
Definition: NumuVars.h:38
ana::Loaders::Go
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
ana::kSlcMinY
const Var kSlcMinY([](const caf::SRProxy *sr){return sr->slc.boxmin.Y()/100.;})
Definition: NumuVars.h:88
inFile
ifstream inFile
Definition: AnaPlotMaker.h:34
ana::Loaders::SetSpillCut
void SetSpillCut(const SpillCut &cut)
Definition: Loaders.cxx:121
ana::kTrueE
const Var kTrueE([](const caf::SRProxy *sr){assert(sr->mc.nnu==1);return sr->mc.nu[0].E;})
Definition: Vars.cxx:85
ana::kRatOfKalHi
const Var kRatOfKalHi([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return(float)((float) sr->trk.kalman.tracks[0].nhit/(float) sr->slc.nhit);})
ana::EnergyBins
const Binning EnergyBins
Definition: ComparisonHints.h:27
outFile
TFile * outFile
Definition: PlotXSec.C:135
ana::kNumuContainFD2017
const Cut kNumuContainFD2017
Definition: NumuCuts2017.h:21
ana::kRemID
const Var kRemID
PID
Definition: Vars.cxx:81
ana::kW_TR
const Cut kW_TR([](const caf::SRProxy *sr){bool IsTrue(sr->trk.kalman.tracks[0].start.X() > XCent && sr->trk.kalman.tracks[0].start.Y() > YCent);return IsTrue;})
ana::kSANumuContPID
const Var kSANumuContPID([](const caf::SRProxy *sr){return sr->sel.cosrej.numuSAcontpid;})
ana::kHadEFrac
const Var kHadEFrac
Definition: NumuVars.h:24
ana::kNHit
const Var kNHit
Definition: Vars.cxx:71
ana::kW_BR
const Cut kW_BR([](const caf::SRProxy *sr){bool IsTrue(sr->trk.kalman.tracks[0].start.X() > XCent && sr->trk.kalman.tracks[0].start.Y()< YCent);return IsTrue;})
ana::kTrkStartZ
const Var kTrkStartZ([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].start.Z()/100;})
Definition: NumuVars.h:53
osc::_OscCalcPMNSOpt
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
ana::kCCE
const Var kCCE
Definition: NumuVars.h:21
ana::kNHitTrk
const Var kNHitTrk([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return(float) sr->trk.kalman.tracks[0].nhit;})
osc::_OscCalcPMNSOpt::SetTh23
void SetTh23(const T &th23) override
Definition: OscCalcPMNSOpt.h:42
ana::CompBins
const Binning CompBins
Definition: ComparisonHints.h:24
QuantNames
const std::string QuantNames[NumQuant]
Definition: DoThePlots.C:207
demo0.loader
loader
Definition: demo0.py:10
ana::CompZBins
const Binning CompZBins
Definition: ComparisonHints.h:25
ana::kNumuPID2017
const Cut kNumuPID2017([](const caf::SRProxy *sr){return(sr->sel.remid.pid > 0.5 &&sr->sel.cvn.numuid > 0.5);})
Definition: NumuCuts2017.h:27
ana::PredictionNoExtrap::SaveTo
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
Definition: PredictionNoExtrap.cxx:79
ana::kNumuCCOptimisedAxis
const HistAxis kNumuCCOptimisedAxis("Reconstructed Neutrino Energy (GeV)", kNumuCCEOptimisedBinning, kCCE)
HistAxis that implements optimised numuCCE from L. Vinton. See docdb 16332. This was close to &#39;custC&#39;...
Definition: HistAxes.h:24
ana::kW_TL
const Cut kW_TL([](const caf::SRProxy *sr){bool IsTrue=(sr->trk.kalman.tracks[0].start.X()< XCent && sr->trk.kalman.tracks[0].start.Y() > YCent);return IsTrue;})
ana::kNumuCosmicRej2017
const Cut kNumuCosmicRej2017([](const caf::SRProxy *sr){return(sr->sel.cosrej.anglekal > 0.5 && sr->sel.cosrej.numucontpid2020 > 0.5 && sr->slc.nhit< 400);})
Definition: NumuCuts2017.h:24
ana::CellsBins
const Binning CellsBins
Definition: ComparisonHints.h:38
ana::kTrkStartX
const Var kTrkStartX([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].start.X()/100;})
Definition: NumuVars.h:51
ComparisonPlots_MC
void ComparisonPlots_MC()
Definition: ComparisonPlots_MC.C:30
ana::kNoShift
const SystShifts kNoShift
Definition: SystShifts.cxx:21
ana::kTrkEndZ
const Var kTrkEndZ([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].stop.Z()/100;})
Definition: NumuVars.h:57
om::cout
OStream cout
Definition: OStream.cxx:6
ana::kTrkLenXY
const Var kTrkLenXY([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;float XLen=(sr->trk.kalman.tracks[0].start.X()-sr->trk.kalman.tracks[0].stop.X())/100;float YLen=(sr->trk.kalman.tracks[0].start.Y()-sr->trk.kalman.tracks[0].stop.Y())/100;return XLen-YLen;})
osc::_OscCalcPMNSOpt::SetDmsq32
void SetDmsq32(const T &dmsq32) override
Definition: OscCalcPMNSOpt.h:39
ana::kKalBakCell
const Var kKalBakCell([](const caf::SRProxy *sr){return sr->sel.contain.kalbakcell;})
ana::kW_BL
const Cut kW_BL([](const caf::SRProxy *sr){bool IsTrue(sr->trk.kalman.tracks[0].start.X()< XCent && sr->trk.kalman.tracks[0].start.Y()< YCent);return IsTrue;})
ana::kRTOvTNuE
const Var kRTOvTNuE
Definition: ComparisonHints.h:86
ana::kTrkEndY
const Var kTrkEndY([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].stop.Y()/100;})
Definition: NumuVars.h:56
ana::NumHitBins
const Binning NumHitBins
Definition: ComparisonHints.h:22
ana::kTrkEndX
const Var kTrkEndX([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].stop.X()/100;})
Definition: NumuVars.h:55
ana::kVisETrk
const Var kVisETrk([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return(float) sr->trk.kalman.tracks[0].calE;})
ana::LengthBins
const Binning LengthBins
Definition: ComparisonHints.h:37
cc
void cc()
Definition: test_ana.C:28
ana::kTrueMuonE
const Var kTrueMuonE([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return 0.f;if(sr->mc.nu[0].prim.empty()) return 0.f;if(std::abs(sr->mc.nu[0].prim[0].pdg)!=13) return 0.f;return float(sr->mc.nu[0].prim[0].p.E);})
Definition: NumuVars.h:107
ana::kTrueHadE
const Var kTrueHadE
Definition: ComparisonHints.h:162
osc::_OscCalcPMNSOpt::SetdCP
void SetdCP(const T &dCP) override
Definition: OscCalcPMNSOpt.h:43
ana::QuantileCutsFromTH2
std::vector< Cut > QuantileCutsFromTH2(TH2 *quantileHist, const HistAxis &independentAxis, const HistAxis &quantileAxis, const unsigned int &numQuantiles, const bool verbose)
: Do the same as the QuantileCuts function but taking in the TH2 instead of making it...
Definition: QuantileCuts.cxx:198
ana::kInBeamSpill
const Cut kInBeamSpill([](const caf::SRProxy *sr){if(sr->spill.run > util::kLastBadTimingRun) return kInBeamSpill_main(sr);else return kInBeamSpill_main(sr)||kInBeamSpill_shifted(sr);}, [](const caf::SRSpillProxy *spill){if(spill->run > util::kLastBadTimingRun) return kInBeamSpill_main.Livetime(spill);else return kInBeamSpill_main.Livetime(spill)+kInBeamSpill_shifted.Livetime(spill);}, [](const caf::SRSpillProxy *spill) -> double{return spill->spillpot;})
Does the event fall inside the window we call the beam spill?
Definition: TimingCuts.h:8
ana::kVisEHad
const Var kVisEHad([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return(float)(sr->slc.calE-sr->trk.kalman.tracks[0].calE);})
ana::RangeBins
const Binning RangeBins
Definition: ComparisonHints.h:30
ana::NumTrkBins
const Binning NumTrkBins
Definition: ComparisonHints.h:34
ana::kReMIdDEDxLLH
const Var kReMIdDEDxLLH
Definition: NumuVars.cxx:556
ana::kStandardSpillCuts
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
ana::TimeFBins
const Binning TimeFBins
Definition: ComparisonHints.h:40
ana::kMaxKalZPos
const Var kMaxKalZPos([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return std::max(sr->trk.kalman.tracks[0].start.Z()/100, sr->trk.kalman.tracks[0].stop.Z()/100);})
ana::PredictionNoExtrap
Prediction that just uses FD MC, with no extrapolation.
Definition: PredictionNoExtrap.h:17
ana::_Cut< caf::SRProxy >
ana::kNumuQuality
const Cut kNumuQuality
Definition: NumuCuts.h:18
ana::kXSecCVWgt2017
const Var kXSecCVWgt2017
Definition: XsecTunes.h:37
e
Float_t e
Definition: plot.C:35
ana::kMuE
const Var kMuE
Definition: NumuVars.h:22
ana::kCVNm
const Var kCVNm
PID
Definition: Vars.cxx:39
ana::kNoCut
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
ComparisonHints.h
ana::kCosNumi
const Var kCosNumi([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks > 0 &&sr->trk.kalman.idxremid!=999){if(sr->hdr.det==1){return sr->trk.kalman.tracks[0].dir.Dot(beamDirND);}if(sr->hdr.det==2){return sr->trk.kalman.tracks[0].dir.Dot(beamDirFD);}}return-5.f;})
Definition: NumuVars.h:43
ana::kSlcMinZ
const Var kSlcMinZ([](const caf::SRProxy *sr){return sr->slc.boxmin.Z()/100.;})
Definition: NumuVars.h:89
ana::kTrkNhits
const Var kTrkNhits([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return 65535;return int(sr->trk.kalman.tracks[0].nhit);})
Definition: NumuVars.h:59
std::asin
T asin(T number)
Definition: d0nt_math.hpp:60
ana::TinyBins
const Binning TinyBins
Definition: ComparisonHints.h:32
ana::RatioBins
const Binning RatioBins
Definition: ComparisonHints.h:23
ana::kEnPHTrk
const Var kEnPHTrk([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return(float)(sr->trk.kalman.tracks[0].calE)/(float)(1.78 *sr->trk.kalman.tracks[0].nhit);})
CutNames
std::vector< std::string > CutNames
Definition: MakeCutFlow.C:49
ana::ECAFType
ECAFType
Definition: Loaders.h:19
string
enum BeamMode string
Generated on Mon Nov 30 2020 05:13:12 for NOvA by   doxygen 1.8.11