neutKEsyst.C
Go to the documentation of this file.
1 #pragma once
2 
3 #include "CAFAna/Cuts/Cuts.h"
10 
11 #include "CAFAna/Core/Binning.h"
12 #include "CAFAna/Core/Spectrum.h"
14 #include "CAFAna/Core/SystShifts.h"
15 
19 #include "CAFAna/Vars/Vars.h"
20 #include "CAFAna/Core/Var.h"
21 #include "CAFAna/Core/MultiVar.h"
22 
26 
30 
31 #include "Utilities/rootlogon.C"
32 
33 #include "TCanvas.h"
34 #include "TH2.h"
35 #include "TProfile.h"
36 #include "TSystem.h"
37 
38 using namespace ana;
39 
40 void neutKEsyst( bool IsFHC = false, bool IsFarDet = true)
41 {
42 
43  std::string sFHC_RHC = ( IsFHC == true ? "fhc" : "rhc" );
44  std::string sFD_ND = ( IsFarDet == true ? "FD" : "ND" );
45 
46  // --- Define my loader, and the file / dataset which it will load...
47  std::string fname = "";
48  if ( IsFHC ) {
49  if (IsFarDet) {
50  fname = "prod_sumdecaf_R17-11-14-prod4reco.neutron-respin.b_fd_genie_nonswap_fhc_nova_v08_full_v1_numu2018"; // numu concats
51 // fname = "prod_caf_R17-11-14-prod4reco.d_fd_genie_nonswap_fhc_nova_v08_full_v1"; // full cafs
52 // fname = "howard_neutronstudy_respin_birks_caf_fd_genie_nonswap_fhc_nova_v08_full_v1";
53  }
54  else {
55  fname = "prod_sumdecaf_R17-11-14-prod4reco.neutron-respin.b_nd_genie_nonswap_fhc_nova_v08_full_v1_numu2018"; // numu concats
56 // fname = "prod_caf_R17-11-14-prod4reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1"; // full cafs
57 // fname = "howard_neutronstudy_respin_birks_caf_nd_genie_nonswap_fhc_nova_v08_full_v1";
58  }
59  }
60  else {
61  if (IsFarDet) {
62  fname = "prod_sumdecaf_R17-11-14-prod4reco.neutron-respin.b_fd_genie_nonswap_rhc_nova_v08_full_v1_numu2018"; // numu concats
63 // fname = "prod_caf_R17-11-14-prod4reco.e_fd_genie_nonswap_rhc_nova_v08_full_v1"; // full cafs
64 // fname = "howard_neutronstudy_respin_birks_caf_fd_genie_nonswap_rhc_nova_v08_full_v1";
65  }
66  else {
67  fname = "prod_sumdecaf_R17-11-14-prod4reco.neutron-respin.b_nd_genie_nonswap_rhc_nova_v08_full_v1_numu2018"; // numu concats
68 // fname = "prod_caf_R17-11-14-prod4reco.e_nd_genie_nonswap_rhc_nova_v08_full_v1"; // full cafs
69 // fname = "howard_neutronstudy_respin_birks_caf_nd_genie_nonswap_rhc_nova_v08_full_v1";
70  }
71  }
72 
73  SpectrumLoader loader(fname);
75 
76  // weights, shifts
79 
80  // --- Define my output directory.
81  std::string sOutFile = "Files_"+sFD_ND+"_"+sFHC_RHC+".root";
82 
83  // --- Define the number of GENIE interaction types I'm looking at....
84  const std::vector<Cut> GENIECuts = { kNoCut };
85  const std::vector<std::string> GENIEStr = { "All" };
86  //const std::vector<Cut> GENIECuts = { kNoCut, kIsDytmanMEC, kIsDIS, kIsRes, kIsCoh, kIsQE };
87  //const std::vector<std::string> GENIEStr = { "All", "MEC", "DIS", "Resonance", "Coherent", "QuasiElactic" };
88  const unsigned int NumGENIE = GENIECuts.size();
89 
90  // --- Define my broad cuts
91  std::vector<std::string> CutNames; std::vector<Cut> kDetCuts;
92  if (IsFarDet) {
93  CutNames.push_back( "FD_2018full" ); kDetCuts.push_back( kNumuQuality && kNumuContainFD2017 && kNumuPID2018 && kNumuCosmicRej2018 );
94  } else {
95  CutNames.push_back( "ND_2018full" ); kDetCuts.push_back( kNumuQuality && kNumuContainND2017 && kNumuPID2018 );
96  }
97  unsigned int nDetCuts = kDetCuts.size();
98 
99  // --- Figure out where to load my quantile boundaries from.
100  std::string fdspecfile = "/pnfs/nova/persistent/analysis/numu/Ana2018/provisional/quantiles/quantiles__" + sFHC_RHC + "_full__numu2018.root";
101 
102  TFile* inFile = TFile::Open( pnfs2xrootd(fdspecfile.c_str()).c_str() );
103  TH2 *FDSpec2D = (TH2*)inFile->FindObjectAny( "FDSpec2D" ); //TString("FD_full_") + TString(sFHC_RHC) );
104  // How many quantiles?
105  const int NHadEFracQuantiles = 4; // defines how many divisions of hadEFrac are used
106  // Make my cuts.
107  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, NHadEFracQuantiles);
108  HadEFracQuantCuts.push_back( kNoCut );
109  std::vector<std::string> QuantNames = { "Quant1", "Quant2", "Quant3", "Quant4", "AllQuant" };
110  const unsigned int nQuants = QuantNames.size();
111 
112 
113  // Extra Vars
114  // ----- Neutron multiplicity
115  const Var kNNeutrons([](const caf::SRProxy* sr){
116  if (sr->mc.nnu == 0) return -5.; //if (sr->mc.nu.size() == 0) return false;
117  return double(sr->mc.nu[0].nneutron);
118  });
119 
120  // ----- Visible energy from neutron daughters
121  const Var kVisEsum_Ndau([](const caf::SRProxy* sr){
122  double visE = -5.0;
123  if(sr->mc.nnu == 0) return visE;
124  visE = 0.0;
125  // loop over primary list
126  for (unsigned int PrimL=0; PrimL < sr->mc.nu[0].prim.size(); ++PrimL) {
127  // is it a neutron?
128  if (sr->mc.nu[0].prim[PrimL].pdg != 2112) continue;
129  visE+=sr->mc.nu[0].prim[PrimL].daughterVisEinslc;
130  }
131  return visE;
132  });
133 
134  const MultiVar kVisE_Ndau([](const caf::SRProxy* sr){
135  std::vector<double> visE;
136  if(sr->mc.nnu == 0) return visE;
137  for (unsigned int PrimL=0; PrimL < sr->mc.nu[0].prim.size(); ++PrimL){
138  if (sr->mc.nu[0].prim[PrimL].pdg != 2112) continue;
139  visE.push_back(sr->mc.nu[0].prim[PrimL].daughterVisEinslc);
140  }
141  return visE;
142  });
143 
144  // ----- Number of 2D prongs
145  const Var kNPng2D([](const caf::SRProxy* sr){
146  double npng = -5.0;
147  if(!sr->vtx.elastic.IsValid) return -5.0;
148  npng = sr->vtx.elastic.fuzzyk.npng2d;
149  return npng;
150  });
151 
152  // ----- Number of 3D prongs
153  const Var kNPng3D([](const caf::SRProxy* sr){
154  double npng = -5.0;
155  if(!sr->vtx.elastic.IsValid) return -5.0;
156  npng = sr->vtx.elastic.fuzzyk.npng;
157  return npng;
158  });
159 
160 
161  // ----- Cut -----
162  const Cut kRemoveNoTruth([](const caf::SRProxy* sr) {
163  if (sr->mc.nnu == 0) return false; //if (sr->mc.nu.size() == 0) return false;
164  return true;
165  });
166 
167 
168  // --- 1D Spectra
169  // Energy plots
170  Spectrum *NuE_nom [NumGENIE][nDetCuts][nQuants];
171  Spectrum *HadE_nom [NumGENIE][nDetCuts][nQuants];
172 
173  Spectrum *NuE_scale_p [NumGENIE][nDetCuts][nQuants];
174  Spectrum *HadE_scale_p [NumGENIE][nDetCuts][nQuants];
175 
176  Spectrum *NuE_scale_m [NumGENIE][nDetCuts][nQuants];
177  Spectrum *HadE_scale_m [NumGENIE][nDetCuts][nQuants];
178 
179  Spectrum *NNeu [NumGENIE][nDetCuts][nQuants];
180  Spectrum *VisE_sum [NumGENIE][nDetCuts][nQuants];
181  Spectrum *VisE [NumGENIE][nDetCuts][nQuants];
182 
183  // --- Define my Spectra which are to be filled from the loader
184  for (unsigned int gen=0; gen<NumGENIE; ++gen) {
185  for (unsigned int det=0; det<nDetCuts; ++det) {
186  for (unsigned int quant=0; quant<nQuants; ++quant) {
187  // Define my new cut.
188  Cut GENIECut = kRemoveNoTruth && GENIECuts[gen] && kDetCuts[det] && HadEFracQuantCuts[quant];
189  std::string MyStr = " for GENIE Int. "+GENIEStr[gen]+", Cut "+CutNames[det]+", "+QuantNames[quant];
190  // --- Now fill our 1D spectra!
191  // Energy spectra
192  NuE_nom [gen][det][quant] = new Spectrum("Reco #nu E (GeV), No Shift" + MyStr, kNumuCCEOptimisedBinning, loader, kCCE, GENIECut, kNoShift, weight);
193  HadE_nom[gen][det][quant] = new Spectrum("Reco Had E (Gev), No Shift" + MyStr, Binning::Simple(50, 0, 3), loader, kHadE, GENIECut, kNoShift, weight);
194 
195 
196  NuE_scale_p [gen][det][quant] = new Spectrum("Reco #nu E (GeV), Neutron Scale Up" + MyStr, kNumuCCEOptimisedBinning, loader, kCCE, GENIECut, SystShifts(scale, +1), weight);
197  HadE_scale_p[gen][det][quant] = new Spectrum("Reco Had E (Gev), Neutron Scale Up" + MyStr, Binning::Simple(50, 0, 3), loader, kHadE, GENIECut, SystShifts(scale, +1), weight);
198 
199  NuE_scale_m [gen][det][quant] = new Spectrum("Reco #nu E (GeV), Neutron Scale Down" + MyStr, kNumuCCEOptimisedBinning, loader, kCCE, GENIECut, SystShifts(scale, -1), weight);
200  HadE_scale_m[gen][det][quant] = new Spectrum("Reco Had E (Gev), Neutron Scale Down" + MyStr, Binning::Simple(50, 0, 3), loader, kHadE, GENIECut, SystShifts(scale, -1), weight);
201 
202  NNeu [gen][det][quant] = new Spectrum("Neutron Multiplicity" +MyStr, Binning::Simple(5,0,5), loader, kNNeutrons, GENIECut, kNoShift, weight);
203  VisE_sum [gen][det][quant] = new Spectrum("Summed visible energy from neutron daughters (GeV)" + MyStr, Binning::Simple(50, 0, 3), loader, kVisEsum_Ndau, GENIECut, kNoShift, weight);
204  VisE [gen][det][quant] = new Spectrum("Visible energy from neutron daughters (GeV)" + MyStr, Binning::Simple(50, 0, 3), loader, kVisE_Ndau, GENIECut, kNoShift, weight);
205 
206  } // Loop through my Spectrum arrays.
207  }
208  }
209 
210  // --- Do it!
211  loader.Go();
212 
213  // --- Open my outfile so that I save my histograms!
214  TFile *OutFile = new TFile(sOutFile.c_str(),"RECREATE");
215 
216  // --- Save plots
217  for (unsigned int gen=0; gen<NumGENIE; ++gen) {
218  for (unsigned int det=0; det<nDetCuts; ++det) {
219  for (unsigned int quant=0; quant<nQuants; ++quant) {
220  std::string MyStr = GENIEStr[gen]+"_"+CutNames[det]+"_"+QuantNames[quant];
221  // --- Save my 1D histograms!
222  // Energy spectra.
223  NuE_nom [gen][det][quant] -> SaveTo( OutFile, TString("NuE_nom_") + TString(MyStr) ) ;
224  HadE_nom [gen][det][quant] -> SaveTo( OutFile, TString("HadE_nom_") + TString(MyStr) ) ;
225  NuE_scale_p [gen][det][quant] -> SaveTo( OutFile, TString("NuE_scale_p_") + TString(MyStr) ) ;
226  HadE_scale_p[gen][det][quant] -> SaveTo( OutFile, TString("HadE_scale_p_") + TString(MyStr) ) ;
227  NuE_scale_m [gen][det][quant] -> SaveTo( OutFile, TString("NuE_scale_m_") + TString(MyStr) ) ;
228  HadE_scale_m[gen][det][quant] -> SaveTo( OutFile, TString("HadE_scale_m_") + TString(MyStr) ) ;
229 
230  NNeu[gen][det][quant] -> SaveTo( OutFile, TString("NNeu_") + TString(MyStr) ) ;
231  VisE_sum[gen][det][quant] -> SaveTo( OutFile, TString("VisE_sum_") + TString(MyStr) ) ;
232  VisE[gen][det][quant] -> SaveTo( OutFile, TString("VisE_") + TString(MyStr) ) ;
233 
234  }
235  }
236  }
237 
238 }
const Var kHadE
Definition: NumuVars.h:23
caf::Proxy< size_t > npng
Definition: SRProxy.h:1997
caf::Proxy< caf::SRFuzzyK > fuzzyk
Definition: SRProxy.h:2018
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var kNPng2D([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-5;return(int) sr->vtx.elastic.fuzzyk.npng2d;})
Definition: NumuVarsExtra.h:57
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Var weight
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
bool IsFHC
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2085
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:574
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
caf::Proxy< short int > nnu
Definition: SRProxy.h:573
void SetSpillCut(const SpillCut &cut)
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: UtilsExt.cxx:237
const NeutronVisEScaleSyst2018 kNeutronVisEScalePrimariesSyst2018(false)
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const Var kNPng3D([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-5;return(int) sr->vtx.elastic.fuzzyk.npng;})
Definition: NumuVarsExtra.h:56
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2077
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
ifstream inFile
Definition: AnaPlotMaker.h:34
Double_t scale
Definition: plot.C:25
const Cut kNumuContainFD2017
Definition: NumuCuts2017.h:21
const SpillCut kStandardDQCuts([](const caf::SRSpillProxy *spill){if(spill->dcmedgematchfrac==0 &&spill->fracdcm3hits==0 &&spill->nmissingdcmslg==0) return bool(spill->isgoodspill); if(spill->det==caf::kNEARDET && (spill->fracdcm3hits > 0.45|| spill->nmissingdcms > 0)) return false; if(spill->eventincomplete) return false; if(spill->det==caf::kFARDET && spill->nmissingdcmslg > 0) return false; if(spill->det==caf::kFARDET && !spill->ismc && spill->dcmedgematchfrac<=0.2) return false;return true;})
Cut out events with a noisy detector or with parts missing.
Definition: SpillCuts.h:16
const Var kCCE
Definition: NumuVars.h:21
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
const std::string QuantNames[NumQuant]
Definition: DoThePlots.C:207
loader
Definition: demo0.py:10
const Var kNNeutrons([](const caf::SRProxy *sr){float NPar=0.;if(DEBUGGING){std::cout<< "\nLooking at event "<< sr->hdr.evt<< ", there are "<< sr->mc.nu.size()<< ", the 0th neutrino is a "<< sr->mc.nu[0].pdg<< ", with Energy "<< sr->mc.nu[0].E<< ". There were a total of "<< sr->mc.nu[0].prim.size()<< " other primaries associated with that neutrino, "<< sr->mc.nu[0].nneutron<< " were neutrons."<< std::endl;if(kIsDytmanMEC(sr)) std::cout<< " The neutrino was a MEC"<< std::endl;else if(kIsDIS(sr)) std::cout<< " The neutrino was a DIS"<< std::endl;else if(kIsRes(sr)) std::cout<< " The neutrino was a Resonance"<< std::endl;else if(kIsCoh(sr)) std::cout<< " The neutrino was a Coherent"<< std::endl;float TotEn=0.0;for(unsigned int PrimL=0;PrimL< sr->mc.nu[0].prim.size();++PrimL){std::cout<< "\t\tLooking at Prim "<< PrimL<< " of "<< sr->mc.nu[0].prim.size()<< ", it was a "<< sr->mc.nu[0].prim[PrimL].pdg<< ", VisE "<< sr->mc.nu[0].prim[PrimL].visE<< ", VisEInslc "<< sr->mc.nu[0].prim[PrimL].visEinslc<< ", E0 "<< sr->mc.nu[0].prim[PrimL].p.E<< ", enteringE "<< sr->mc.nu[0].prim[PrimL].enteringE<< ", totEscE "<< sr->mc.nu[0].prim[PrimL].totEscE<< std::endl;TotEn+=sr->mc.nu[0].prim[PrimL].totEscE;}std::cout<< "\tThe total uncontained energy was "<< TotEn<< std::endl;NPar=sr->mc.nu[0].nneutron;}return NPar;})
std::vector< float > Spectrum
Definition: Constants.h:570
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 SystShifts kNoShift
Definition: SystShifts.cxx:21
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2097
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
caf::Proxy< bool > IsValid
Definition: SRProxy.h:2017
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
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 kXSecCVWgt2018
Definition: XsecTunes.h:49
void neutKEsyst(bool IsFHC=false, bool IsFarDet=true)
Definition: neutKEsyst.C:40
const Binning kNumuCCEOptimisedBinning
Optimised binning for numuCCE from L. Vinton. See docdb 16332. This was close to &#39;custC&#39; in that talk...
Definition: Binnings.cxx:28
const std::vector< std::string > GENIEStr
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...
TFile * OutFile
caf::Proxy< size_t > npng2d
Definition: SRProxy.h:1998
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2105
const Cut kNumuQuality
Definition: NumuCuts.h:18
#define for
Definition: msvc_pragmas.h:3
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:114
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
std::vector< std::string > CutNames
Definition: MakeCutFlow.C:49