CutFlow_NearDet.C
Go to the documentation of this file.
1 #include "CAFAna/Analysis/Exposures.h"
2 #include "CAFAna/Analysis/Prod4Loaders.h"
3 
4 #include "CAFAna/Cuts/Cuts.h"
5 #include "3FlavorAna/Cuts/NumuCuts.h"
6 #include "3FlavorAna/Cuts/NumuCuts2017.h"
7 #include "3FlavorAna/Cuts/NumuCuts2018.h"
8 #include "3FlavorAna/Cuts/QuantileCuts.h"
9 #include "CAFAna/Cuts/SpillCuts.h"
10 #include "CAFAna/Cuts/TimingCuts.h"
11 
12 #include "CAFAna/Core/SpectrumLoader.h"
13 #include "CAFAna/Core/Spectrum.h"
14 #include "CAFAna/Core/ReweightableSpectrum.h"
15 
16 #include "CAFAna/Vars/GenieWeights.h"
17 #include "CAFAna/Vars/HistAxes.h"
18 #include "3Flavor/Vars/HistAxes.h"
19 #include "CAFAna/Vars/Vars.h"
20 #include "3FlavorAna/Vars/Binnings.h"
21 #include "3FlavorAna/Vars/NumuVars.h"
22 #include "CAFAna/Vars/PPFXWeights.h"
23 #include "CAFAna/Systs/XSecSystLists.h"
24 
25 using namespace ana;
26 
27 #include "TStyle.h"
28 #include "TCanvas.h"
29 #include "TH1.h"
30 #include "TLegend.h"
31 #include "TLatex.h"
32 
33 #include <fstream>
34 
35 void CutFlow_NearDet( bool IsMC, bool IsFHC, bool IsTest=false )
36 {
37  // ---- Figure out what files I'm loading / what I'm calling things....
38  std::string sFHC_RHC = (IsFHC == true ? "fhc" : "rhc" );
39  std::string sIsMC = (IsMC == true ? "Mont" : "Data" );
40  std::string OutName = "NearDet_CutFlow_"+sIsMC+".root";
41  std::string MyLoader = "";
42  if (IsTest) {
43  if ( IsMC && IsFHC) { MyLoader = "karlwarb_ND_FHC_MC"; }
44  else if ( IsMC && !IsFHC) { MyLoader = "karlwarb_ND_RHC_MC"; }
45  else if (!IsMC && IsFHC) { MyLoader = "karlwarb_ND_FHC_Data"; }
46  else if (!IsMC && !IsFHC) { MyLoader = "karlwarb_ND_RHC_Data"; }
47  } else {
48  if ( IsMC && IsFHC) { MyLoader = "prod_caf_R17-11-14-prod4reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1"; }
49  else if ( IsMC && !IsFHC) { MyLoader = "prod_caf_R17-11-14-prod4reco.e_nd_genie_nonswap_rhc_nova_v08_full_v1"; }
50  else if (!IsMC && IsFHC) { MyLoader = "prod_caf_R17-09-05-prod4recopreview.f_nd_numi_fhc_full_v1_addShortSimpleCVN_goodruns"; }
51  else if (!IsMC && !IsFHC) { MyLoader = "prod_caf_R17-09-05-prod4recopreview.f_nd_numi_rhc_full_v1_addShortSimpleCVN_goodruns"; }
52  }
53  // --- Determine what file I'm loading for the quantile cuts.
54  std::string QuantFile = "/pnfs/nova/persistent/analysis/numu/Ana2018/provisional/quantiles/quantiles__"+sFHC_RHC+"_full__numu2018.root";
55 
56  std::cout << "\n\n My macro is configured as follows; IsMC " << IsMC << ", IsFHC " << IsFHC << "\nTherefore;"
57  << "\n\t sFHC_RHC = " << sFHC_RHC
58  << "\n\t sIsMC = " << sIsMC
59  << "\n\t OutName = " << OutName
60  << "\n\t MyLoader = " << MyLoader
61  << "\n\t QuantFile = " << QuantFile
62  << std::endl;
63 
64  // --- Set my loader.
65  SpectrumLoader loader( MyLoader );
66  loader.SetSpillCut(kStandardSpillCuts);
67 
68  // --- Define my cuts and their names.
69  std::vector<Cut> TieredCuts;
70  std::vector<std::string> CutNames;
71  CutNames.emplace_back("NoCut"); TieredCuts.emplace_back(kNoCut);
72  CutNames.emplace_back("Quality"); TieredCuts.emplace_back(kNumuQuality);
73  CutNames.emplace_back("Containment"); TieredCuts.emplace_back(kNumuQuality && kNumuContainND2017);
74  CutNames.emplace_back("Particle_ID"); TieredCuts.emplace_back(kNumuQuality && kNumuContainND2017 && kNumuPID2018 );
75 
76  size_t NCuts = TieredCuts.size();
77 
78  // Get stuff for the hadEFrac Cut:
79  TFile* inFile = TFile::Open( pnfs2xrootd(QuantFile).c_str() );
80  TH2 *FDSpec2D = (TH2*)inFile->FindObjectAny( "FDSpec2D" );
81  // How many quantiles?
82  const int NHadEFracQuantiles = 4;
83  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, NHadEFracQuantiles);
84 
85  // --- Define my objects.
86  Spectrum* FullMC [NCuts][NHadEFracQuantiles+1];
87  Spectrum* NuMuCC [NCuts][NHadEFracQuantiles+1];
88  Spectrum* AntiNuMuCC[NCuts][NHadEFracQuantiles+1];
89  Spectrum* AllOtherCC[NCuts][NHadEFracQuantiles+1];
90  Spectrum* AllNC [NCuts][NHadEFracQuantiles+1];
91 
92  // --- Set my objects.
93  for ( size_t cut = 0; cut < NCuts; ++cut ) {
94  for ( int quant=0; quant<NHadEFracQuantiles; ++quant ) {
95 
96  // What Cut and Weight am I using?
97  Cut ThisCut = TieredCuts[cut] && HadEFracQuantCuts[quant] && kInBeamSpill;
98  Var ThisWgt = kPPFXFluxCVWgt*kXSecCVWgt2018;
99  if (!IsMC) ThisWgt = kUnweighted;
100 
101  // --- Now declare my spectrum.
102  FullMC [cut][quant] = new Spectrum( "Reconstructed energy [GeV]", kNumuCCEOptimisedBinning, loader, kCCE, ThisCut , kNoShift, ThisWgt );
103  NuMuCC [cut][quant] = new Spectrum( "Reconstructed energy [GeV]", kNumuCCEOptimisedBinning, loader, kCCE, ThisCut && kIsNumuCC && kIsNu, kNoShift, ThisWgt );
104  AntiNuMuCC[cut][quant] = new Spectrum( "Reconstructed energy [GeV]", kNumuCCEOptimisedBinning, loader, kCCE, ThisCut && kIsNumuCC && !kIsNu, kNoShift, ThisWgt );
105  AllOtherCC[cut][quant] = new Spectrum( "Reconstructed energy [GeV]", kNumuCCEOptimisedBinning, loader, kCCE, ThisCut && !kIsNumuCC && !kIsNC, kNoShift, ThisWgt );
106  AllNC [cut][quant] = new Spectrum( "Reconstructed energy [GeV]", kNumuCCEOptimisedBinning, loader, kCCE, ThisCut && kIsNC, kNoShift, ThisWgt );
107 
108  // --- Finally, push back an all quantile spectrum
109  if (quant == NHadEFracQuantiles-1) {
110  // Change the cut...
111  ThisCut = TieredCuts[cut] && kInBeamSpill;
112  // Declare the spects.
113  FullMC [cut][quant+1] = new Spectrum( "Reconstructed energy [GeV]", kNumuCCEOptimisedBinning, loader, kCCE, ThisCut , kNoShift, ThisWgt );
114  NuMuCC [cut][quant+1] = new Spectrum( "Reconstructed energy [GeV]", kNumuCCEOptimisedBinning, loader, kCCE, ThisCut && kIsNumuCC && kIsNu, kNoShift, ThisWgt );
115  AntiNuMuCC[cut][quant+1] = new Spectrum( "Reconstructed energy [GeV]", kNumuCCEOptimisedBinning, loader, kCCE, ThisCut && kIsNumuCC && !kIsNu, kNoShift, ThisWgt );
116  AllOtherCC[cut][quant+1] = new Spectrum( "Reconstructed energy [GeV]", kNumuCCEOptimisedBinning, loader, kCCE, ThisCut && !kIsNumuCC && !kIsNC, kNoShift, ThisWgt );
117  AllNC [cut][quant+1] = new Spectrum( "Reconstructed energy [GeV]", kNumuCCEOptimisedBinning, loader, kCCE, ThisCut && kIsNC, kNoShift, ThisWgt );
118  }
119  }
120  }
121 
122  // --- Set the loader off!
123  loader.Go();
124 
125  // --- Define an outfile.
126  TFile *OutFile = new TFile(OutName.c_str(), "RECREATE");
127  OutFile -> cd();
128  for(size_t cut = 0; cut < NCuts; ++cut){
129  for ( int quant=0; quant<NHadEFracQuantiles; ++quant ) {
130  // --- Now figure out a name for my spectrum.
131  std::string ThisDir = CutNames[cut]+"_Quant"+std::to_string(quant+1);
132  FullMC [cut][quant] -> SaveTo( OutFile, TString("FullMC_" )+TString(ThisDir) ) ;
133  NuMuCC [cut][quant] -> SaveTo( OutFile, TString("NuMuCC_" )+TString(ThisDir) ) ;
134  AntiNuMuCC[cut][quant] -> SaveTo( OutFile, TString("AntiNuMuCC_")+TString(ThisDir) ) ;
135  AllOtherCC[cut][quant] -> SaveTo( OutFile, TString("AllOtherCC_")+TString(ThisDir) ) ;
136  AllNC [cut][quant] -> SaveTo( OutFile, TString("AllNC_" )+TString(ThisDir) ) ;
137 
138  // --- Finally, the last spectra is an all quantile spectrum
139  if (quant == NHadEFracQuantiles-1) {
140  ThisDir = CutNames[cut]+"_QuantAll";
141  FullMC [cut][quant] -> SaveTo( OutFile, TString("FullMC_" )+TString(ThisDir) ) ;
142  NuMuCC [cut][quant] -> SaveTo( OutFile, TString("NuMuCC_" )+TString(ThisDir) ) ;
143  AntiNuMuCC[cut][quant] -> SaveTo( OutFile, TString("AntiNuMuCC_")+TString(ThisDir) ) ;
144  AllOtherCC[cut][quant] -> SaveTo( OutFile, TString("AllOtherCC_")+TString(ThisDir) ) ;
145  AllNC [cut][quant] -> SaveTo( OutFile, TString("AllNC_" )+TString(ThisDir) ) ;
146  }
147  }
148  }
149 } // End of function
ana::kIsNumuCC
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
NCuts
size_t NCuts
Definition: MakeCutFlow.C:50
ana
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
ana::SaveTo
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
Definition: LoadFromFile.cxx:105
ana::kPPFXFluxCVWgt
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
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
IsFHC
bool IsFHC
Definition: ND_DataMC_Comp_Systs.C:4
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
CutFlow_NearDet
void CutFlow_NearDet(bool IsMC, bool IsFHC, bool IsTest=false)
Definition: CutFlow_NearDet.C:35
ana::SpectrumLoaderBase::SetSpillCut
void SetSpillCut(const SpillCut &cut)
Definition: SpectrumLoaderBase.cxx:64
ana::pnfs2xrootd
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: UtilsExt.cxx:237
ana::kIsNu
const Cut kIsNu([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return sr->mc.nu[0].pdg > 0;})
Definition: TruthCuts.h:54
ana::Spectrum
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
inFile
ifstream inFile
Definition: AnaPlotMaker.h:34
ana::kIsNC
const Cut kIsNC([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return!sr->mc.nu[0].iscc;})
Is this a Neutral Current event?
Definition: TruthCuts.h:8
ana::kCCE
const Var kCCE
Definition: NumuVars.h:21
ana::SpectrumLoader::Go
virtual void Go() override
Load all the registered spectra.
Definition: SpectrumLoader.cxx:66
demo0.loader
loader
Definition: demo0.py:10
cmf::Spectrum
std::vector< float > Spectrum
Definition: Constants.h:728
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::kNoShift
const SystShifts kNoShift
Definition: SystShifts.cxx:22
om::cout
OStream cout
Definition: OStream.cxx:6
cut
const Cut cut
Definition: exporter_fd.C:30
ana::SpectrumLoader
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
Definition: SpectrumLoader.h:21
ana::kNumuPID2018
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
ana::kXSecCVWgt2018
const Var kXSecCVWgt2018
Definition: XsecTunes.h:48
ana::kNumuCCEOptimisedBinning
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
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::kStandardSpillCuts
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
OutFile
TFile * OutFile
Definition: NuMu2020_BasicPIDPlots_Plot.C:70
art::to_string
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
ana::_Cut< caf::SRProxy >
ana::kNumuQuality
const Cut kNumuQuality
Definition: NumuCuts.h:18
ana::_Var< caf::SRProxy >
ana::kUnweighted
const Var kUnweighted
The simplest possible Var, always 1. Used as a default weight.
Definition: Var.h:96
cd
c cd(1)
ana::kNoCut
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
CutNames
std::vector< std::string > CutNames
Definition: MakeCutFlow.C:49
string
enum BeamMode string
Generated on Sun Feb 28 2021 14:11:55 for NOvA by   doxygen 1.8.11