SystematicComp.C
Go to the documentation of this file.
1 #include "CAFAna/Cuts/Cuts.h"
7 #include "CAFAna/Core/Spectrum.h"
9 #include "CAFAna/Core/Loaders.h"
10 //#include "CAFAna/nus/NuSPlotFunctions.h"
11 
15 
17 
18 #include "TCanvas.h"
19 #include "TH2.h"
20 #include "TProfile.h"
21 
22 #include "TFile.h"
23 
24 #include <vector>
25 #include <sstream>
26 #include <string>
27 
28 using namespace ana;
29 
31  std::cout << "Please provide;"
32  << "An input samweb definition. \n"
33  << "Whether you want to run all cuts or not. I will assume full cuts.\n"
34  << std::endl;
35 }
36 
37 void SystematicComp( std::string InputDef, bool IsData = false, bool IsCosmics = false, bool RunProd4Cuts=true, bool RunAllCuts=false, bool Debug=false )
38 {
39  // --- Set my loader
40  SpectrumLoader loader( InputDef );
42 
43  // Add the detector specific ones.
44  Cut kDetCut = kNoCut;
45  if ( InputDef.find("nd") != std::string::npos) {
46  kDetCut = kNumuContainND2017;
47  } else {
48  kDetCut = kNumuContainFD2017;
49  }
50 
51  // What PID and Cos Rej Cut am I using?
52  Cut kMyPIDCut = kNumuPID2018; Cut kMyCosRejCut = kNumuCosmicRej2018;
53  if (!RunProd4Cuts) {
54  Cut kMyPIDCut = kNumuPID2017; Cut kMyCosRejCut = kNumuCosmicRej2017;
55  }
56 
57  // What weights do I want to use?
58  Var MyWeight = kUnweighted;
59  if ( IsData && IsCosmics) { std::cout << "Cosmics --> Sideband" << std::endl; MyWeight = kTimingSidebandWeight; }
60  if (!IsData && RunProd4Cuts) { std::cout << "Ana18MC --> XSec2018" << std::endl; MyWeight = kPPFXFluxCVWgt*kXSecCVWgt2018; }
61  if (!IsData && !RunProd4Cuts) { std::cout << "Ana17MC --> XSec2017" << std::endl; MyWeight = kPPFXFluxCVWgt*kXSecCVWgt2017; }
62 
63  // --- The timing cut.
64  Cut TimingCut = kInBeamSpill;
65  if (IsCosmics) { TimingCut = kInTimingSideband; std::cout << "Cosmics --> SidebandTime" << std::endl; }
66 
67  // ----- Define my cuts....
68  std::vector<Cut> kDefinedCut;
69  unsigned int NumCuts = 1;
70  if (RunAllCuts) NumCuts=4;
71  std::string CutNames[NumCuts];
72  if (NumCuts == 4) { // --- NumCuts == 4
73  // --- Cut
74  kDefinedCut.push_back( TimingCut && kNumuQuality ); // Good quality NuMu events
75  kDefinedCut.push_back( TimingCut && kNumuQuality && kDetCut ); // That the event is contained.
76  kDefinedCut.push_back( TimingCut && kNumuQuality && kDetCut ); // Reject cosmics
77  kDefinedCut.push_back( TimingCut && kNumuQuality && kDetCut && kMyPIDCut ); // Reject Neutral current events
78  // --- CutNames
79  CutNames[0] = "_QualCuts";
80  CutNames[1] = "_DetCuts" ;
81  CutNames[2] = "_CosmCuts";
82  CutNames[3] = "_PIDCuts";
83  // --- CosmicRej
84  if ( InputDef.find("fd") != std::string::npos ) {
85  kDefinedCut[2] = kDefinedCut[2] && kMyCosRejCut;
86  kDefinedCut[3] = kDefinedCut[3] && kMyCosRejCut;
87  }
88  }
89  else if (NumCuts == 1) { // --- NumCuts == 1
90  // --- Cut
91  kDefinedCut.push_back( TimingCut && kNumuQuality && kDetCut && kMyPIDCut ); // Go straight to full cuts.
92  // --- CutNames
93  CutNames[0] = "_AllCuts";
94  // --- CosmicRej
95  if ( InputDef.find("fd") != std::string::npos ) {
96  kDefinedCut[0] = kDefinedCut[0] && kMyCosRejCut;
97  }
98 
99  if (IsData) kDefinedCut[0] = kDefinedCut[0] && kNumu2017Veto;
100  }
101 
102  // --- Load the Quantile Cuts.
103  // Figure out the file name should be...
104  std::string sFHC_RHC = "fhc";
105  if (InputDef.find("rhc") != std::string::npos ) sFHC_RHC = "rhc";
106  std::string fdspecfile = "/pnfs/nova/persistent/analysis/numu/Ana2018/official/Quantiles/quantiles__"+sFHC_RHC+"_full__numu2018.root";
107  std::cout << "\nLoading quantile cuts from " << fdspecfile << "\n" << std::endl;
108  // Load the file...
109  TFile* inFile = TFile::Open( pnfs2xrootd(fdspecfile).c_str() );
110  TH2 *FDSpec2D = (TH2*)inFile->FindObjectAny( "FDSpec2D" );
111  // Add get the cuts....
112  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, 4);
113 
114  // ----- Define my quantiles.
115  const unsigned int NumQuants = 5;
116  Cut kQuantCut[NumQuants] = { kNoCut, HadEFracQuantCuts[0], HadEFracQuantCuts[1], HadEFracQuantCuts[2], HadEFracQuantCuts[3] };
117  const std::string QuantNames[NumQuants] = {
118  "_AllQuants",
119  "_Quant1",
120  "_Quant2",
121  "_Quant3",
122  "_Quant4"
123  };
124 
125  // --- Decide which kind of sample I am looking at....
126  // If Nominal
127  std::string SystName = "Nominal";
128  // If Calib shift
129  if (InputDef.find("calib-shift") != std::string::npos) {
130  if (InputDef.find("neg") != std::string::npos) SystName = "Calib_Neg";
131  else if (InputDef.find("pos") != std::string::npos) SystName = "Calib_Pos";
132  else if (InputDef.find("func") != std::string::npos) SystName = "Calib_Func";
133  } else if (InputDef.find("light") != std::string::npos) {
134  // If looking at Light shifts
135  if (InputDef.find("lightdown") != std::string::npos) SystName = "Light_LightDown";
136  else if (InputDef.find("lightup") != std::string::npos) SystName = "Light_LightUp";
137  else if (InputDef.find("noshift") != std::string::npos) SystName = "Light_NoShift";
138  } else if (InputDef.find("ckv-proton") != std::string::npos) SystName = "Light_CkvProton";
139 
140  // --- I need to add the period to the output name.
141  std::string TPer = "full_";
142  if (InputDef.find("period1") != std::string::npos ) TPer = "Period1_";
143  else if (InputDef.find("period2") != std::string::npos ) TPer = "Period2_";
144  else if (InputDef.find("period3") != std::string::npos ) TPer = "Period3_";
145  else if (InputDef.find("period4") != std::string::npos ) TPer = "Period4_";
146  else if (InputDef.find("period5") != std::string::npos ) TPer = "Period5_";
147  else if (InputDef.find("period6") != std::string::npos ) TPer = "Period6_";
148 
149  // --- If data need to flag that, and if MC I need to state the swap.
150  std::string TSwap = "";
151  if ( IsData && !IsCosmics) TSwap = "Data_";
152  if ( IsData && IsCosmics) TSwap = "Cosmics_";
153  if (!IsData) {
154  std::string TSwap = "nonswap_";
155  if (InputDef.find("fluxswap") != std::string::npos ) TSwap = "fluxswap_";
156  if (InputDef.find("tau" ) != std::string::npos ) TSwap = "tau_";
157  }
158 
159  // --- Make the name of my output file.
160  std::string OutputName = "SystematicComps_";
161  if (InputDef.find("sumdecaf") != std::string::npos ) OutputName += "concat_";
162  else if (InputDef.find("restricted") != std::string::npos ) OutputName += "concat_";
163  else OutputName += "caf_";
164  OutputName += ( InputDef.find("fhc") != std::string::npos ? "fhc_" :"rhc_" );
165  OutputName += TPer + TSwap;
166  OutputName += ( InputDef.find("fd") != std::string::npos ? "FD_" :"ND_" );
167  OutputName += ( RunProd4Cuts == true ? "Ana18Cuts_":"Ana17Cuts_");
168  OutputName += SystName;
169  OutputName += ".root";
170 
171  // ---- Write some output.
172  std::cout << "Before I declare anything, the following is set;\n"
173  << "\tInput definit is " << InputDef << "\n"
174  << "\tOutput file name " << OutputName << "\n"
175  << "\tNumber of Cuts " << NumCuts << "\n"
176  << "\tIs debug mode is " << Debug << "\n"
177  << std::endl;
178 
179 
180  // --- Declare all of my spectra, before setting them in a loop...
181  // *** Some key ones first.
182  Spectrum *sReconstEnergy[NumCuts][NumQuants]; // --- Reconstructed energy --- Vars.h ===> kCCE ---> sr->energy.numu.E
183  Spectrum *sSliceTimeFull[NumCuts][NumQuants]; // --- Mean time of slice Full --- NumuVars.h ===> kSliceTime ---> sr->slc.meantime/1000
184  Spectrum *sSliceTimeZoom[NumCuts][NumQuants]; // --- Mean time of slice Zoom --- NumuVars.h ===> kSliceTime ---> sr->slc.meantime/1000
185  Spectrum *sMuonEnergy [NumCuts][NumQuants]; // --- Muon energy in the slice --- Vars.h ===> kMuE --> energy.numu.recomuonE
186  Spectrum *sMuonEnPerHit [NumCuts][NumQuants]; // --- Muon energy per hit --- NumuVars.h ===> kTrkEPerNHit --> sr->trk.kalman.tracks[0].calE / sr->trk.kalman.tracks[0].nhit
187  Spectrum *sHadronEnergy [NumCuts][NumQuants]; // --- Hadronic energy --- NumuVars.h ===> kHadE --> sr->energy.numu.E - sr->energy.numu.recomuonE
188  Spectrum *sHadroEnPerHit[NumCuts][NumQuants]; // --- Hadronic energy per hit --- NumuVars.h ===> kHadEPerNHit --> function.
189  Spectrum *sHadFracEnergy[NumCuts][NumQuants]; // --- Hadronic energy fraction --- NumuVars.h ===> kHadEFrac --> kHadE / kCCE
190  Spectrum *sHitsPerSlice [NumCuts][NumQuants]; // --- Number of hits per slice --- Vars.h ===> kNHit --> slc.nhit
191  Spectrum *sRemIDScore [NumCuts][NumQuants]; // --- The RemID score --- Vars.h ===> kRemID --> sr->sel.remid.pid
192  Spectrum *sCVNCosmicScor[NumCuts][NumQuants]; // --- CVN Cosmic score --- Vars.h ===> sr->sel.cvn.output[390]
193  Spectrum *sCVNNuMuIDScor[NumCuts][NumQuants]; // --- CVN NuMuID score --- Vars.h ===> sr->sel.cvnProd3Train.numuid
194  Spectrum *sCVNNuMuID2017[NumCuts][NumQuants]; // --- CVN NuMuID score --- Vars.h ===> sr->sel.cvn2017.numuid
195  Spectrum *sNuMuContPID [NumCuts][NumQuants]; // --- Third analysis BDT --- NuMuVars.h ===> sr->sel.cosrej.numucontpid
196  Spectrum *sSANuMuContPID[NumCuts][NumQuants]; // --- Second analysis BDT --- Above.
197  // *** Some important Truth based ones.
198  Spectrum *sTrueNuEnergy [NumCuts][NumQuants]; // --- True neutrino energy --- Vars.h ===> kTrueE ---> sr->mc.nu[0].E
199  Spectrum *sReTrOverTrNuE[NumCuts][NumQuants]; // --- Reco - True / True Nu E --- Above.
200 
201  // *** Basic positional plots.
202  Spectrum *sTrkStartX[NumCuts][NumQuants]; // --- Start Pos of Track --- NumuVars.h ===> kTrkStartX --> trk.kalman.tracks[0].start.X()/100
203  Spectrum *sTrkStartY[NumCuts][NumQuants]; // --- Start Pos of Track --- NumuVars.h ===> kTrkStartY --> trk.kalman.tracks[0].start.Y()/100
204  Spectrum *sTrkStartZ[NumCuts][NumQuants]; // --- Start Pos of Track --- NumuVars.h ===> kTrkStartZ --> trk.kalman.tracks[0].start.Z()/100
205  Spectrum *sTrkEndX [NumCuts][NumQuants]; // --- End Pos of Track --- NumuVars.h ===> kTrkEndX --> trk.kalman.tracks[0].stop.X()/100
206  Spectrum *sTrkEndY [NumCuts][NumQuants]; // --- End Pos of Track --- NumuVars.h ===> kTrkEndY --> trk.kalman.tracks[0].stop.Y()/100
207  Spectrum *sTrkEndZ [NumCuts][NumQuants]; // --- End Pos of Track --- NumuVars.h ===> kTrkEndZ --> trk.kalman.tracks[0].stop.Z()/100
208  Spectrum *sTrkLenXY [NumCuts][NumQuants]; // --- Diff in Track Len X Y --- Above. ===> kTrkLenXY --> (tracks[0].start.X - tracks[0].stop.X) / 100
209  //*
210  // *** Track based ones...
211  Spectrum *sNumKalTracks [NumCuts][NumQuants]; // --- Number of Kal trk in slc --- NuMuVars.h ===> sr->trk.kalman.ntracks
212  Spectrum *sNumKalTrHits [NumCuts][NumQuants]; // --- Number of hits in Kal tr --- NuMuVars.h ===> sr->trk.kalman.tracks[0].nhit
213  Spectrum *sRatKalHitSlc [NumCuts][NumQuants]; // --- Ratio of Kal hits in slc --- Above.
214  Spectrum *sKalTrLength [NumCuts][NumQuants]; // --- Length of Kalman track --- NuMuVars.h ===> sr->trk.kalman.tracks[0].len / 100
215  Spectrum *sKalTrBeamAng [NumCuts][NumQuants]; // --- Angle Kal track & beam --- NuMuVars.h ===> kCosNumi --> sr->trk.kalman.tracks[0].dir.Dot(beamDirND)
216  Spectrum *sKalMostFwdCel[NumCuts][NumQuants]; // --- Most forward Kal tr hit --- Above.
217  Spectrum *sKalMostBakCel[NumCuts][NumQuants]; // --- Most backward Kal tr hit --- Above.
218  Spectrum *sKalTrVer_MaxY[NumCuts][NumQuants]; // --- Highest Y Kalman tr vert --- Above.
219  Spectrum *sKalTrVer_MaxZ[NumCuts][NumQuants]; // --- Highest Z Kalman tr vert --- Above.
220  Spectrum *sKalTrDir_Y [NumCuts][NumQuants]; // --- Y direction of Kal track --- NuMuVars.h ===> sr->trk.kalman.tracks[0].dir.Y()
221  Spectrum *sScattKalTrLen[NumCuts][NumQuants]; // --- Scattering over tr len --- Above.
222 
223  // *** Hit based ones...
224  Spectrum *sFirstHitCell [NumCuts][NumQuants]; // --- The first cell hit --- Above.
225  Spectrum *sLastHitCell [NumCuts][NumQuants]; // --- The last cell hit --- Above.
226  Spectrum *sMaxActivity_Y[NumCuts][NumQuants]; // --- Maximum Y pos of activty --- NuMuVars.h ===> sr->slc.boxmax.Y() / 100.
227  Spectrum *sMinActivity_Y[NumCuts][NumQuants]; // --- Minimum Y pos of activty --- NuMuVars.h ===> sr->slc.boxmin.Y() / 100.
228  Spectrum *sMaxActivity_Z[NumCuts][NumQuants]; // --- Maximum Z pos of activty --- NuMuVars.h ===> sr->slc.boxmax.Z() / 100.
229  Spectrum *sMinActivity_Z[NumCuts][NumQuants]; // --- Minimum Z pos of activty --- NuMuVars.h ===> sr->slc.boxmin.Z() / 100.
230  Spectrum *sMinCellToEdge[NumCuts][NumQuants]; // --- Minimum cells to edge --- Above.
231 
232  // *** Visible energy quanta.
233  // Slice
234  Spectrum *sVisEInSlice [NumCuts][NumQuants]; // --- The Vis En in the slice --- Above. ---> kVisESlc ---> sr->slc.calE
235  Spectrum *sNHitInSlice [NumCuts][NumQuants]; // --- The Num.Hits in the slic --- Above. ---> kNHitSlc ---> sr->slc.nhit
236  Spectrum *sVisEPerHitSlc[NumCuts][NumQuants]; // --- Visibile E per hit in sl --- Above. ---> kEnPHSlc ---> kVisESlc / (1.78 * kNHitSlc)
237  // Track
238  Spectrum *sVisEInTrack [NumCuts][NumQuants]; // --- The Vis En of the track --- Above. ---> kVisETrk ---> sr->trk.kalman.tracks[0].calE
239  Spectrum *sNHitInTrack [NumCuts][NumQuants]; // --- The NumHits in the track --- Above. ---> kNHitTrk ---> sr->trk.kalman.tracks[0].nhit
240  Spectrum *sVisEPerHitTrk[NumCuts][NumQuants]; // --- Visibile E PerHit in trk --- Above. ---> kEnPHTrk ---> kVisESlc / (1.78 * kNHitTrk)
241  // Hadronic
242  Spectrum *sVisEInHadron [NumCuts][NumQuants]; // --- The Hadronic Vis Energy --- Above. ---> kVisEHad ---> (kVisESlc - kVisETrk)
243  Spectrum *sNHitInHadron [NumCuts][NumQuants]; // --- The Num of hadronic Hits --- Above. ---> kNHitHad ---> (kNHitHad - kNHitTrk)
244  Spectrum *sVisEPerHitHad[NumCuts][NumQuants]; // --- Visibile E PerHit hadron --- Above. ---> kEnPHHad ---> kVisEHad / (1.78 * kNHitHad)
245 
246  // *** RemID score / inputs
247  Spectrum *sRemIDScatLLH [NumCuts][NumQuants]; // --- RemID scattering angle --- NuMuVars.h ===> kReMIdScatLLH --> sel.remid.scatllh
248  Spectrum *sRemIDdEdxLLH [NumCuts][NumQuants]; // --- RemID dE/dx --- NuMuVars.h ===> kReMIdDEDxLLH --> sel.remid.dedxllh
249  Spectrum *sRemIDMeasFrac[NumCuts][NumQuants]; // --- RemID fraction of planes --- NuMuVars.h ===> kReMIdMeasFrac -> sel.remid.measfrac
250 
251  // *** Some calibration quantities
252  Spectrum *sCor_RecEn[NumCuts][NumQuants][4];
253  Spectrum *sCor_MuoEn[NumCuts][NumQuants][4];
254  Spectrum *sCor_HadEn[NumCuts][NumQuants][4];
255  Spectrum *sCor_HFrEn[NumCuts][NumQuants][4];
256  // Some special truth based ones.
257  Spectrum *sCor_RecTr[NumCuts][NumQuants][4];
258  Spectrum *sCor_MuoTr[NumCuts][NumQuants][4];
259  Spectrum *sCor_HadTr[NumCuts][NumQuants][4];
260  Spectrum *sCor_HFrTr[NumCuts][NumQuants][4];
261  Spectrum *sCor_ReTr [NumCuts][NumQuants][4];
262  //*/
263  const Cut kCorner[4] = { kW_TL, kW_TR, kW_BL, kW_BR };
264  const std::string CutName[4] = { "_TopL", "_TopR", "_BotL", "_BotR" };
265 
266  // --- And now set what is going to go into each of my spectra.
267  for (unsigned int cc=0; cc<NumCuts; ++cc) {
268  // --- Loop through my quantiles.
269  for (unsigned int qq=0; qq<NumQuants; ++qq) {
270  // What is the cut this time?
271  const Cut kThisCut = kDefinedCut[cc] && kQuantCut[qq];
272 
273  // *** Some key ones first...
274  sReconstEnergy[cc][qq] = new Spectrum( "The reconstructed energy", OptEnBins, loader, kCCE , kThisCut, kNoShift, MyWeight );
275  sSliceTimeFull[cc][qq] = new Spectrum( "The time of slice - Full", TimeFBins , loader, kSliceTime , kThisCut, kNoShift, MyWeight );
276  sSliceTimeZoom[cc][qq] = new Spectrum( "The time of slice - Zoom", TimeZBins , loader, kSliceTime , kThisCut, kNoShift, MyWeight );
277  sMuonEnergy [cc][qq] = new Spectrum( "The muon energy " , LowEnBins , loader, kMuE , kThisCut, kNoShift, MyWeight );
278  sMuonEnPerHit [cc][qq] = new Spectrum( "The muon energy per hit" , EnPHitBins, loader, kTrkEPerNHit , kThisCut, kNoShift, MyWeight );
279  sHadronEnergy [cc][qq] = new Spectrum( "The hadronic energy" , LowEnBins , loader, kHadE , kThisCut, kNoShift, MyWeight );
280  sHadroEnPerHit[cc][qq] = new Spectrum( "The hadronic En per hit ", EnPHitBins, loader, kHadEPerNHit , kThisCut, kNoShift, MyWeight );
281  sHadFracEnergy[cc][qq] = new Spectrum( "The hadronic energy frac", RatioBins , loader, kHadEFrac , kThisCut, kNoShift, MyWeight );
282  sHitsPerSlice [cc][qq] = new Spectrum( "Number of hits per slice", NumHitBins, loader, kNHit , kThisCut, kNoShift, MyWeight );
283  sRemIDScore [cc][qq] = new Spectrum( "The RemID score" , RemIDBins , loader, kRemID , kThisCut, kNoShift, MyWeight );
284  sCVNCosmicScor[cc][qq] = new Spectrum( "CVN Cosmic score" , RatioBins , loader, kCVNcos , kThisCut, kNoShift, MyWeight );
285  sCVNNuMuIDScor[cc][qq] = new Spectrum( "CVN NuMu score" , RatioBins , loader, kCVNm , kThisCut, kNoShift, MyWeight );
286  sCVNNuMuID2017[cc][qq] = new Spectrum( "2017 CVN NuMu score" , RatioBins , loader, kCVNm2017 , kThisCut, kNoShift, MyWeight );
287  sNuMuContPID [cc][qq] = new Spectrum( "The NuMu containment PID", RatioBins , loader, kNumuContPID , kThisCut, kNoShift, MyWeight );
288  sSANuMuContPID[cc][qq] = new Spectrum( "SA NuMu containment PID" , RatioBins , loader, kSANumuContPID, kThisCut, kNoShift, MyWeight );
289  // *** Some important truth based ones.
290  sTrueNuEnergy [cc][qq] = new Spectrum( "The true neutrino energy", EnergyBins, loader, kTrueE , kThisCut, kNoShift, MyWeight );
291  sReTrOverTrNuE[cc][qq] = new Spectrum( "The R-T ov T Nu energy" , RangeBins , loader, kRTOvTNuE, kThisCut, kNoShift, MyWeight );
292 
293  // *** Positional plots
294  sTrkStartX[cc][qq] = new Spectrum( "Starting X Pos of Track", YRangeBins, loader, kTrkStartX, kThisCut, kNoShift, MyWeight );
295  sTrkStartY[cc][qq] = new Spectrum( "Starting Y Pos of Track", YRangeBins, loader, kTrkStartY, kThisCut, kNoShift, MyWeight );
296  sTrkStartZ[cc][qq] = new Spectrum( "Starting Z Pos of Track", ZRangeBins, loader, kTrkStartZ, kThisCut, kNoShift, MyWeight );
297  sTrkEndX[cc][qq] = new Spectrum( "Ending X Pos of Track" , YRangeBins, loader, kTrkEndX , kThisCut, kNoShift, MyWeight );
298  sTrkEndY[cc][qq] = new Spectrum( "Ending Y Pos of Track" , YRangeBins, loader, kTrkEndY , kThisCut, kNoShift, MyWeight );
299  sTrkEndZ[cc][qq] = new Spectrum( "Ending Z Pos of Track" , ZRangeBins, loader, kTrkEndZ , kThisCut, kNoShift, MyWeight );
300  sTrkLenXY[cc][qq] = new Spectrum( "Difference in XY Tr Len", ZRangeBins, loader, kTrkLenXY , kThisCut, kNoShift, MyWeight );
301  //*
302  // *** Track based ones...
303  sNumKalTracks [cc][qq] = new Spectrum( "Number of Kalman Tracks" , NumTrkBins, loader, kNKalman , kThisCut, kNoShift, MyWeight );
304  sNumKalTrHits [cc][qq] = new Spectrum( "Number of Kal Track hits", NumHitBins, loader, kTrkNhits , kThisCut, kNoShift, MyWeight );
305  sRatKalHitSlc [cc][qq] = new Spectrum( "Ratio of KalHit in Slice", RatioBins , loader, kRatOfKalHi, kThisCut, kNoShift, MyWeight );
306  sKalTrLength [cc][qq] = new Spectrum( "Length of Kalman Track" , LengthBins, loader, kTrkLength , kThisCut, kNoShift, MyWeight );
307  sKalTrBeamAng [cc][qq] = new Spectrum( "Angle of KalTr to Beam" , RatioBins , loader, kCosNumi , kThisCut, kNoShift, MyWeight );
308  sKalMostFwdCel[cc][qq] = new Spectrum( "Most forward KalTr Cell" , CellsBins , loader, kKalFwdCell, kThisCut, kNoShift, MyWeight );
309  sKalMostBakCel[cc][qq] = new Spectrum( "Most backward KalTr Cell", CellsBins , loader, kKalBakCell, kThisCut, kNoShift, MyWeight );
310  sKalTrVer_MaxY[cc][qq] = new Spectrum( "Most pos Y KalTr Vertex" , YRangeBins, loader, kMaxKalYPos, kThisCut, kNoShift, MyWeight );
311  sKalTrVer_MaxZ[cc][qq] = new Spectrum( "Most pos Z KalTr Vertex" , ZRangeBins, loader, kMaxKalZPos, kThisCut, kNoShift, MyWeight );
312  sKalTrDir_Y [cc][qq] = new Spectrum( "Y Direction of KalTr" , CompBins , loader, kDirY , kThisCut, kNoShift, MyWeight );
313  sScattKalTrLen[cc][qq] = new Spectrum( "Scatt over KalTr length" , TinyBins , loader, kScattTrLen, kThisCut, kNoShift, MyWeight );
314 
315  // *** Hit based ones...
316  sFirstHitCell [cc][qq] = new Spectrum( "The first cell hit in sl", CellsBins , loader, kFirstCell , kThisCut, kNoShift, MyWeight );
317  sLastHitCell [cc][qq] = new Spectrum( "The last cell hit in sl" , CellsBins , loader, kLastCell , kThisCut, kNoShift, MyWeight );
318  sMaxActivity_Y[cc][qq] = new Spectrum( "Largest Y loc with a hit", YRangeBins, loader, kSlcMaxY , kThisCut, kNoShift, MyWeight );
319  sMinActivity_Y[cc][qq] = new Spectrum( "Smalest Y loc with a hit", YRangeBins, loader, kSlcMinY , kThisCut, kNoShift, MyWeight );
320  sMaxActivity_Z[cc][qq] = new Spectrum( "Largest Z loc with a hit", ZRangeBins, loader, kSlcMaxZ , kThisCut, kNoShift, MyWeight );
321  sMinActivity_Z[cc][qq] = new Spectrum( "Largest Z loc with a hit", ZRangeBins, loader, kSlcMinZ , kThisCut, kNoShift, MyWeight );
322  sMinCellToEdge[cc][qq] = new Spectrum( "Min cells hit to edge" , CellsBins , loader, kMinCellEdg, kThisCut, kNoShift, MyWeight );
323 
324  // *** Visible energy quanta.
325  // Slice
326  sVisEInSlice [cc][qq] = new Spectrum( "Visible Energy in slice" , EnergyBins, loader, kVisESlc, kThisCut, kNoShift, MyWeight );
327  sNHitInSlice [cc][qq] = new Spectrum( "Number of hits in slice" , NumHitBins, loader, kNHitSlc, kThisCut, kNoShift, MyWeight );
328  sVisEPerHitSlc[cc][qq] = new Spectrum( "Visible E per hit in slc", EnPHitBins, loader, kEnPHSlc, kThisCut, kNoShift, MyWeight );
329  // Track
330  sVisEInTrack [cc][qq] = new Spectrum( "Visible Energy in track" , EnergyBins, loader, kVisETrk, kThisCut, kNoShift, MyWeight );
331  sNHitInTrack [cc][qq] = new Spectrum( "Number of hits in track" , NumHitBins, loader, kNHitTrk, kThisCut, kNoShift, MyWeight );
332  sVisEPerHitTrk[cc][qq] = new Spectrum( "Visible E per hit in trk", EnPHitBins, loader, kEnPHTrk, kThisCut, kNoShift, MyWeight );
333  // Hadronic
334  sVisEInHadron [cc][qq] = new Spectrum( "Visible Hadronic Energy" , EnergyBins, loader, kVisEHad, kThisCut, kNoShift, MyWeight );
335  sNHitInHadron [cc][qq] = new Spectrum( "Number of hadronic hits" , NumHitBins, loader, kNHitHad, kThisCut, kNoShift, MyWeight );
336  sVisEPerHitHad[cc][qq] = new Spectrum( "Visible E per hadron hit", EnPHitBins, loader, kEnPHHad, kThisCut, kNoShift, MyWeight );
337 
338  // *** RemID score / inputs
339  sRemIDScatLLH [cc][qq] = new Spectrum( "RemID Scattering angle" , CompZBins , loader, kReMIdScatLLH , kThisCut, kNoShift, MyWeight );
340  sRemIDdEdxLLH [cc][qq] = new Spectrum( "RemID dE/dx" , CompBins , loader, kReMIdDEDxLLH , kThisCut, kNoShift, MyWeight );
341  sRemIDMeasFrac[cc][qq] = new Spectrum( "RemID fraction of planes", RatioBins , loader, kReMIdMeasFrac , kThisCut, kNoShift, MyWeight );
342 
343  // *** Calibration quantities
344  for (int cor=0; cor<4; ++cor) {
345  const Cut kCorCut = kThisCut && kCorner[cor];
346  sCor_RecEn[cc][qq][cor] = new Spectrum( "Reconstructed Energy with corner cut", EnergyBins, loader, kCCE , kCorCut, kNoShift, MyWeight);
347  sCor_MuoEn[cc][qq][cor] = new Spectrum( "Muon energy with corner cut" , EnergyBins, loader, kMuE , kCorCut, kNoShift, MyWeight);
348  sCor_HadEn[cc][qq][cor] = new Spectrum( "Hadronic energy with corner cut" , EnergyBins, loader, kHadE , kCorCut, kNoShift, MyWeight);
349  sCor_HFrEn[cc][qq][cor] = new Spectrum( "Hadronic energy frac with corner cut", RatioBins , loader, kHadEFrac, kCorCut, kNoShift, MyWeight);
350  // Additional truth based ones.
351  sCor_RecTr[cc][qq][cor] = new Spectrum( "True Energy with corner cut" , EnergyBins, loader, kTrueE , kCorCut, kNoShift, MyWeight);
352  sCor_MuoTr[cc][qq][cor] = new Spectrum( "True Mu En with corner cut" , EnergyBins, loader, kTrueMuonE, kCorCut, kNoShift, MyWeight);
353  sCor_HadTr[cc][qq][cor] = new Spectrum( "True Had En with corner cut" , EnergyBins, loader, kTrueHadE , kCorCut, kNoShift, MyWeight);
354  sCor_HFrTr[cc][qq][cor] = new Spectrum( "True HadEFrac with corner cut", RatioBins , loader, kTrueHFrE , kCorCut, kNoShift, MyWeight);
355  sCor_ReTr [cc][qq][cor] = new Spectrum( "Reco-True/True w. corner cut" , RangeBins , loader, kRTOvTNuE , kCorCut, kNoShift, MyWeight);
356  }
357  //*/
358  } // Quantile loop
359  } // Cut loop
360 
361  // --- Do it!
362  loader.Go();
363 
364  // --- Open my outfile so that I save histograms!
365  TFile *outFile = new TFile(OutputName.c_str(),"RECREATE");
366  std::cerr << "Made my output file...." << std::endl;
367 
368  // --- Finally, write them to disk!
369  for (unsigned int cc=0; cc<NumCuts; ++cc) {
370  // --- Loop through my quantiles.
371  for (unsigned int qq=0; qq<NumQuants; ++qq) {
372  std::string MyAppend = CutNames[cc] + QuantNames[qq];
373  // *** Some key ones first...
374  sReconstEnergy[cc][qq] -> SaveTo( outFile, TString("sReconstEnergy")+TString(MyAppend) ) ;
375  sSliceTimeFull[cc][qq] -> SaveTo( outFile, TString("sSliceTimeFull")+TString(MyAppend) ) ;
376  sSliceTimeZoom[cc][qq] -> SaveTo( outFile, TString("sSliceTimeZoom")+TString(MyAppend) ) ;
377  sMuonEnergy [cc][qq] -> SaveTo( outFile, TString("sMuonEnergy") +TString(MyAppend) ) ;
378  sMuonEnPerHit [cc][qq] -> SaveTo( outFile, TString("sMuonEnPerHit") +TString(MyAppend) ) ;
379  sHadronEnergy [cc][qq] -> SaveTo( outFile, TString("sHadronEnergy") +TString(MyAppend) ) ;
380  sHadroEnPerHit[cc][qq] -> SaveTo( outFile, TString("sHadroEnPerHit")+TString(MyAppend) ) ;
381  sHadFracEnergy[cc][qq] -> SaveTo( outFile, TString("sHadFracEnergy")+TString(MyAppend) ) ;
382  sHitsPerSlice [cc][qq] -> SaveTo( outFile, TString("sHitsPerSlice") +TString(MyAppend) ) ;
383  sRemIDScore [cc][qq] -> SaveTo( outFile, TString("sRemIDScore") +TString(MyAppend) ) ;
384  sCVNCosmicScor[cc][qq] -> SaveTo( outFile, TString("sCVNCosmicScor")+TString(MyAppend) ) ;
385  sCVNNuMuIDScor[cc][qq] -> SaveTo( outFile, TString("sCVNNuMuIDScor")+TString(MyAppend) ) ;
386  sCVNNuMuID2017[cc][qq] -> SaveTo( outFile, TString("sCVNNuMuID2017")+TString(MyAppend) ) ;
387  sNuMuContPID [cc][qq] -> SaveTo( outFile, TString("sNuMuContPID") +TString(MyAppend) ) ;
388  sSANuMuContPID[cc][qq] -> SaveTo( outFile, TString("sSANuMuContPID")+TString(MyAppend) ) ;
389  // *** Some important truth based ones.
390  sTrueNuEnergy [cc][qq] -> SaveTo( outFile, TString("sTrueNuEnergy") +TString(MyAppend) ) ;
391  sReTrOverTrNuE[cc][qq] -> SaveTo( outFile, TString("sReTrOverTrNuE")+TString(MyAppend) ) ;
392 
393  // *** Positional plots
394  sTrkStartX[cc][qq] -> SaveTo( outFile, TString("sTrkStartX")+TString(MyAppend) ) ;
395  sTrkStartY[cc][qq] -> SaveTo( outFile, TString("sTrkStartY")+TString(MyAppend) ) ;
396  sTrkStartZ[cc][qq] -> SaveTo( outFile, TString("sTrkStartZ")+TString(MyAppend) ) ;
397  sTrkEndX[cc][qq] -> SaveTo( outFile, TString("sTrkEndX") +TString(MyAppend) ) ;
398  sTrkEndY[cc][qq] -> SaveTo( outFile, TString("sTrkEndY") +TString(MyAppend) ) ;
399  sTrkEndZ[cc][qq] -> SaveTo( outFile, TString("sTrkEndZ") +TString(MyAppend) ) ;
400  sTrkLenXY[cc][qq] -> SaveTo( outFile, TString("sTrkLenXY") +TString(MyAppend) ) ;
401  //*
402  // *** Track based ones...
403  sNumKalTracks [cc][qq] -> SaveTo( outFile, TString("sNumKalTracks") +TString(MyAppend) ) ;
404  sNumKalTrHits [cc][qq] -> SaveTo( outFile, TString("sNumKalTrHits") +TString(MyAppend) ) ;
405  sRatKalHitSlc [cc][qq] -> SaveTo( outFile, TString("sRatKalHitSlc") +TString(MyAppend) ) ;
406  sKalTrLength [cc][qq] -> SaveTo( outFile, TString("sKalTrLength") +TString(MyAppend) ) ;
407  sKalTrBeamAng [cc][qq] -> SaveTo( outFile, TString("sKalTrBeamAng") +TString(MyAppend) ) ;
408  sKalMostFwdCel[cc][qq] -> SaveTo( outFile, TString("sKalMostFwdCel")+TString(MyAppend) ) ;
409  sKalMostBakCel[cc][qq] -> SaveTo( outFile, TString("sKalMostBakCel")+TString(MyAppend) ) ;
410  sKalTrVer_MaxY[cc][qq] -> SaveTo( outFile, TString("sKalTrVer_MaxY")+TString(MyAppend) ) ;
411  sKalTrVer_MaxZ[cc][qq] -> SaveTo( outFile, TString("sKalTrVer_MaxZ")+TString(MyAppend) ) ;
412  sKalTrDir_Y [cc][qq] -> SaveTo( outFile, TString("sKalTrDir_Y") +TString(MyAppend) ) ;
413  sScattKalTrLen[cc][qq] -> SaveTo( outFile, TString("sScattKalTrLen")+TString(MyAppend) ) ;
414 
415  // *** Hit based ones...
416  sFirstHitCell [cc][qq] -> SaveTo( outFile, TString("sFirstHitCell") +TString(MyAppend) ) ;
417  sLastHitCell [cc][qq] -> SaveTo( outFile, TString("sLastHitCell") +TString(MyAppend) ) ;
418  sMaxActivity_Y[cc][qq] -> SaveTo( outFile, TString("sMaxActivity_Y")+TString(MyAppend) ) ;
419  sMinActivity_Y[cc][qq] -> SaveTo( outFile, TString("sMinActivity_Y")+TString(MyAppend) ) ;
420  sMinCellToEdge[cc][qq] -> SaveTo( outFile, TString("sMinCellToEdge")+TString(MyAppend) ) ;
421  sMaxActivity_Z[cc][qq] -> SaveTo( outFile, TString("sMaxActivity_Z")+TString(MyAppend) ) ;
422  sMinActivity_Z[cc][qq] -> SaveTo( outFile, TString("sMinActivity_Z")+TString(MyAppend) ) ;
423 
424  // *** Visible energy quanta.
425  // Slice
426  sVisEInSlice [cc][qq] -> SaveTo( outFile, TString("sVisEInSlice") +TString(MyAppend) ) ;
427  sNHitInSlice [cc][qq] -> SaveTo( outFile, TString("sNHitInSlice") +TString(MyAppend) ) ;
428  sVisEPerHitSlc[cc][qq] -> SaveTo( outFile, TString("sVisEPerHitSlc")+TString(MyAppend) ) ;
429  // Track
430  sVisEInTrack [cc][qq] -> SaveTo( outFile, TString("sVisEInTrack") +TString(MyAppend) ) ;
431  sNHitInTrack [cc][qq] -> SaveTo( outFile, TString("sNHitInTrack") +TString(MyAppend) ) ;
432  sVisEPerHitTrk[cc][qq] -> SaveTo( outFile, TString("sVisEPerHitTrk")+TString(MyAppend) ) ;
433  // Hadronic
434  sVisEInHadron [cc][qq] -> SaveTo( outFile, TString("sVisEInHadron") +TString(MyAppend) ) ;
435  sNHitInHadron [cc][qq] -> SaveTo( outFile, TString("sNHitInHadron") +TString(MyAppend) ) ;
436  sVisEPerHitHad[cc][qq] -> SaveTo( outFile, TString("sVisEPerHitHad")+TString(MyAppend) ) ;
437 
438  // *** RemID score / inputs
439  sRemIDScatLLH [cc][qq] -> SaveTo( outFile, TString("sRemIDScatLLH") +TString(MyAppend) ) ;
440  sRemIDdEdxLLH [cc][qq] -> SaveTo( outFile, TString("sRemIDdEdxLLH") +TString(MyAppend) ) ;
441  sRemIDMeasFrac[cc][qq] -> SaveTo( outFile, TString("sRemIDMeasFrac")+TString(MyAppend) ) ;
442 
443  // *** Calibration quantities
444  for (int cor=0; cor<4; ++cor) {
445  sCor_RecEn[cc][qq][cor] -> SaveTo( outFile, TString("sCor_RecEn")+TString(MyAppend)+TString(CutName[cor]) ) ;
446  sCor_MuoEn[cc][qq][cor] -> SaveTo( outFile, TString("sCor_MuoEn")+TString(MyAppend)+TString(CutName[cor]) ) ;
447  sCor_HadEn[cc][qq][cor] -> SaveTo( outFile, TString("sCor_HadEn")+TString(MyAppend)+TString(CutName[cor]) ) ;
448  sCor_HFrEn[cc][qq][cor] -> SaveTo( outFile, TString("sCor_HFrEn")+TString(MyAppend)+TString(CutName[cor]) ) ;
449  // Additional truth based ones.
450  sCor_RecTr[cc][qq][cor] -> SaveTo( outFile, TString("sCor_RecTr")+TString(MyAppend)+TString(CutName[cor]) ) ;
451  sCor_MuoTr[cc][qq][cor] -> SaveTo( outFile, TString("sCor_MuoTr")+TString(MyAppend)+TString(CutName[cor]) ) ;
452  sCor_HadTr[cc][qq][cor] -> SaveTo( outFile, TString("sCor_HadTr")+TString(MyAppend)+TString(CutName[cor]) ) ;
453  sCor_HFrTr[cc][qq][cor] -> SaveTo( outFile, TString("sCor_HFrTr")+TString(MyAppend)+TString(CutName[cor]) ) ;
454  sCor_ReTr [cc][qq][cor] -> SaveTo( outFile, TString("sCor_ReTr" )+TString(MyAppend)+TString(CutName[cor]) ) ;
455  }
456  //*/
457  std::cerr << "\tWritten quantile set " << qq << " for cut set " << cc << " to disk!!!" << std::endl;
458  }
459  std::cerr << "\n\tWritten cut set " << cc << " to disk!!!" << std::endl;
460  }
461 
462 }
const Var kHadE
Definition: NumuVars.h:23
const Binning ZRangeBins
const Var kNKalman
Definition: NumuVars.cxx:540
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
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));})
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);})
const Var kReMIdScatLLH
Definition: NumuVars.cxx:555
const Var kVisESlc([](const caf::SRProxy *sr){return sr->slc.calE;})
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
const Var kNHitSlc([](const caf::SRProxy *sr){return sr->slc.nhit;})
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Var kSlcMaxZ([](const caf::SRProxy *sr){return sr->slc.boxmax.Z()/100.;})
Definition: NumuVars.h:93
const Var kSlcMaxY([](const caf::SRProxy *sr){return sr->slc.boxmax.Y()/100.;})
Definition: NumuVars.h:92
const Var kNumuContPID
Definition: NumuVars.cxx:553
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
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);})
const Binning YRangeBins
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
const Binning EnPHitBins
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);})
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));})
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
const Var kFirstCell([](const caf::SRProxy *sr){return sr->slc.firstcell;})
const Var kSliceTime([](const caf::SRProxy *sr){return sr->slc.meantime/1000;})
Definition: NumuVars.h:34
OStream cerr
Definition: OStream.cxx:7
const Binning LowEnBins
const Var kTrueHFrE([](const caf::SRProxy *sr){float HadE=kTrueHadE(sr);float TruE=kTrueE(sr);return(HadE/TruE);})
const Var kLastCell([](const caf::SRProxy *sr){return sr->slc.lastcell;})
const Binning TimeZBins
void SetSpillCut(const SpillCut &cut)
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: UtilsExt.cxx:237
const Var kTrkEPerNHit([](const caf::SRProxy *sr){return(sr->trk.kalman.tracks[0].calE/sr->trk.kalman.tracks[0].nhit);})
const Var kReMIdMeasFrac
Definition: NumuVars.cxx:557
const Var kKalFwdCell([](const caf::SRProxy *sr){return sr->sel.contain.kalfwdcell;})
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
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);})
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
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
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const Var kSlcMinY([](const caf::SRProxy *sr){return sr->slc.boxmin.Y()/100.;})
Definition: NumuVars.h:88
const Cut kInTimingSideband([](const caf::SRProxy *sr){if(sr->spill.run > util::kLastBadTimingRun) return(kInTimingSideband_before(sr)|| kInTimingSideband_after(sr));else return(kInTimingSideband_before(sr)|| kInTimingSideband_afterA(sr)|| kInTimingSideband_afterB(sr));}, [](const caf::SRSpillProxy *spill){if(spill->run > util::kLastBadTimingRun) return(kInTimingSideband_before.Livetime(spill)+ kInTimingSideband_after.Livetime(spill));else return(kInTimingSideband_before.Livetime(spill)+ kInTimingSideband_afterA.Livetime(spill)+ kInTimingSideband_afterB.Livetime(spill));}, [](const caf::SRSpillProxy *spill){return 0;})
Definition: TimingCuts.h:12
const Cut kNumuCosmicRej2018([](const caf::SRProxy *sr){return(sr->sel.cosrej.anglekal > 0.5 && sr->sel.cosrej.numucontpid2019 > 0.53 && sr->slc.nhit< 400 && sr->sel.nuecosrej.pngptp< 0.9 );})
Definition: NumuCuts2018.h:19
const Binning RemIDBins
ifstream inFile
Definition: AnaPlotMaker.h:34
const Cut kNumu2017Veto([](const caf::SRProxy *sr){return sr->sel.veto.keep;})
cosmic veto cut for FD
Definition: NumuCuts2017.h:30
const Var kTrueE([](const caf::SRProxy *sr){assert(sr->mc.nnu==1);return sr->mc.nu[0].E;})
Definition: Vars.cxx:85
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);})
const Binning EnergyBins
TFile * outFile
Definition: PlotXSec.C:135
const Cut kNumuContainFD2017
Definition: NumuCuts2017.h:21
const Var kRemID
PID
Definition: Vars.cxx:81
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;})
const Var kSANumuContPID([](const caf::SRProxy *sr){return sr->sel.cosrej.numuSAcontpid;})
const Var kHadEFrac
Definition: NumuVars.h:24
const Var kNHit
Definition: Vars.cxx:71
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;})
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
const Var kCCE
Definition: NumuVars.h:21
const Var kNHitTrk([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return(float) sr->trk.kalman.tracks[0].nhit;})
virtual void Go() override
Load all the registered spectra.
const Binning CompBins
const std::string QuantNames[NumQuant]
Definition: DoThePlots.C:207
loader
Definition: demo0.py:10
const Binning CompZBins
const Cut kNumuPID2017([](const caf::SRProxy *sr){return(sr->sel.remid.pid > 0.5 &&sr->sel.cvn.numuid > 0.5);})
Definition: NumuCuts2017.h:27
std::vector< float > Spectrum
Definition: Constants.h:610
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
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;})
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
const Binning CellsBins
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
const SystShifts kNoShift
Definition: SystShifts.cxx:21
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
OStream cout
Definition: OStream.cxx:6
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;})
const Var kKalBakCell([](const caf::SRProxy *sr){return sr->sel.contain.kalbakcell;})
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;})
const Var kRTOvTNuE
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
const Binning NumHitBins
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
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Var kVisETrk([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return(float) sr->trk.kalman.tracks[0].calE;})
const Binning LengthBins
void cc()
Definition: test_ana.C:28
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
const Cut kNumuPID2018([](const caf::SRProxy *sr){std::cout<< "ERROR::kNumuPID2018, cutting on both cvnProd3Train and cvn2017."<< " Neither branch exists anymore. Returning False."<< std::endl;abort();return false;})
Definition: NumuCuts2018.h:22
const Var kTrueHadE
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
const Binning OptEnBins
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...
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
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);})
const Binning RangeBins
const Binning NumTrkBins
const Var kReMIdDEDxLLH
Definition: NumuVars.cxx:556
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
const Binning TimeFBins
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);})
const Cut kNumuQuality
Definition: NumuCuts.h:18
const Var kXSecCVWgt2017
Definition: XsecTunes.h:37
const Var kMuE
Definition: NumuVars.h:22
void SystematicComp()
const Var kCVNm
PID
Definition: Vars.cxx:39
const Var kUnweighted
The simplest possible Var, always 1. Used as a default weight.
Definition: Var.h:96
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
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
const Var kSlcMinZ([](const caf::SRProxy *sr){return sr->slc.boxmin.Z()/100.;})
Definition: NumuVars.h:89
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
const Binning TinyBins
const Binning RatioBins
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);})
std::vector< std::string > CutNames
Definition: MakeCutFlow.C:49
enum BeamMode string