NuMu2019_BasicPIDPlots_ND.C
Go to the documentation of this file.
4 
5 #include "CAFAna/Core/Loaders.h"
7 #include "CAFAna/Core/Var.h"
8 
14 
15 #include "TH2.h"
16 #include "TFile.h"
18 
19 #include <iostream>
20 
21 using namespace ana;
22 
23 // ---------------------------------------------- //
24 // This has to be ran out of R19-04-17-2019ana.a //
25 // ---------------------------------------------- //
26 
27 void NuMu2019_BasicPIDPlots_ND( bool isFHC, bool isData, std::string period = "full" ) {
28 
29  std::string polarity = "fhc"; if (!isFHC ) polarity = "rhc" ;
30  std::string data = "mc" ; if ( isData) data = "data";
31  std::string detector = "nd" ;
32 
33  std::cout << "\n================================= \n"
34  << "\n\t isFHC --> " << isFHC << " --> " << polarity
35  << "\n\t isData --> " << isData << " --> " << data
36  << "\n\t Period --> " << period
37  << "\n\t Far detector isFD --> " << detector
38  << "\n================================= \n"
39  << std::endl;
40 
41  // Make my vector of cuts.
42  std::vector<std::pair<Cut, std::string> > TrueCuts, SelCuts, MyCuts;
43  TrueCuts.clear(); SelCuts.clear(); MyCuts.clear();
44  // Create my true cuts...
45  TrueCuts. push_back( std::make_pair( kNoCut , "_NoTrue" ) );
46  if (!isData) {
47  TrueCuts.push_back( std::make_pair( kHasNeutrino , "_TrueNu" ) );
48  TrueCuts.push_back( std::make_pair( kIsNumuCC && kIsNu , "_TrueNuMu" ) );
49  TrueCuts.push_back( std::make_pair( kIsNumuCC && kIsAntiNu, "_TrueAntiNuMu" ) );
50  }
51  // Create my selection cuts...
52  SelCuts.push_back( std::make_pair( kNumuQuality&&kNumuContainND2017, "_QualCont" ) );
53  SelCuts.push_back( std::make_pair( kNumuCutND2018 , "_FullCut" ) );
54  // Put them together...
55  for (size_t sel=0; sel<SelCuts.size(); ++sel) {
56  for (size_t tr=0; tr<TrueCuts.size(); ++tr) {
57  Cut kThisCut = SelCuts[sel].first && TrueCuts[tr].first;
58  std::string kThisName = SelCuts[sel].second + TrueCuts[tr].second;
59  MyCuts.push_back( std::make_pair( kThisCut, kThisName ) );
60  }
61  }
62  const size_t NCuts = MyCuts.size();
63 
64  // What XSec and PPFX weights am I using?
65  Var kMyWeight = kXSecCVWgt2018 * kPPFXFluxCVWgt;
66  if (isData) kMyWeight = kUnweighted;
67 
68  // Figure out my definition, and then declare my spectrum loader.
69  std::string MyDef = "";
70  if ( isFHC && !isData) MyDef ="prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_v1";
71  else if ( isFHC && isData) MyDef ="prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_numi_fhc_full_v1_addShortSimpleCVN_goodruns";
72  else if (!isFHC && !isData) MyDef ="prod_caf_R17-11-14-prod4reco.neutron-respin.b_nd_genie_nonswap_rhc_nova_v08_full_v1";
73  else if (!isFHC && isData) MyDef ="prod_caf_R17-11-14-prod4reco.neutron-respin.b_nd_numi_rhc_full_v1_goodruns";
74  SpectrumLoader loader( MyDef );
76 
77  // Declare my spectra
78  Spectrum *sNuRecoE [NCuts];
79  Spectrum *sReMIdScore [NCuts];
80  Spectrum *sProngScoremuon[NCuts];
81  Spectrum *sCosRejScoreP4 [NCuts];
82  Spectrum *sCVNScoremuon [NCuts];
83 
84  // Declare some binnings.
86  const Binning kSimpBins = Binning::Simple(102, -0.01, 1.01 );
87 
88  for (size_t cc=0; cc < NCuts; ++cc) {
89  Cut kThisCut = MyCuts[cc].first;
90 
91  sNuRecoE [cc] = new Spectrum( "Reco Nu E" , kOptEnBins , loader, kCCE , kThisCut, kNoShift, kMyWeight );
92  sReMIdScore [cc] = new Spectrum( "ReMId Score" , kRemidBinning, loader, kRemID , kThisCut, kNoShift, kMyWeight );
93  sProngScoremuon[cc] = new Spectrum( "Prong CVN Muon Score" , kSimpBins , loader, kCVNBestMuonScore, kThisCut, kNoShift, kMyWeight );
94  sCosRejScoreP4 [cc] = new Spectrum( "Prod4 CosRej Score" , kSimpBins , loader, kNumuContPID , kThisCut, kNoShift, kMyWeight );
95  sCVNScoremuon [cc] = new Spectrum( "CVN Muon Score" , kSimpBins , loader, kCVNm , kThisCut, kNoShift, kMyWeight );
96  }
97 
98  // Set my loader going.
99  loader.Go();
100 
101  // Make my Output file
102  std::string OutName = "BasicPIDPlots_2019_"+detector+"_"+data+"_"+polarity+"_"+period+".root";
103  TFile *OutFile = TFile::Open(OutName.c_str(), "RECREATE");
104  for (size_t cc=0; cc < NCuts; ++cc) {
105  std::string ThisCutNa = MyCuts[cc].second;
106 
107  sNuRecoE [cc] -> SaveTo( OutFile, TString("RecoNuE" )+TString(ThisCutNa) ) ;
108  sReMIdScore [cc] -> SaveTo( OutFile, TString("ReMIdScore" )+TString(ThisCutNa) ) ;
109  sProngScoremuon[cc] -> SaveTo( OutFile, TString("ProngCVNMuonScore" )+TString(ThisCutNa) ) ;
110  sCosRejScoreP4 [cc] -> SaveTo( OutFile, TString("Prod4CosRejScore" )+TString(ThisCutNa) ) ;
111  sCVNScoremuon [cc] -> SaveTo( OutFile, TString("CVNMuonScore" )+TString(ThisCutNa) ) ;
112  }
113  OutFile -> Close();
114 }
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
size_t NCuts
Definition: MakeCutFlow.C:50
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
const Binning kRemidBinning
Binning for plotting remid attractively.
Definition: Binning.cxx:80
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
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 Cut kNumuCutND2018
Definition: NumuCuts2018.h:41
const Cut kIsAntiNu([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return sr->mc.nu[0].pdg< 0;})
Is this truly an antineutrino?
Definition: TruthCuts.h:53
void SetSpillCut(const SpillCut &cut)
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
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
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const XML_Char const XML_Char * data
Definition: expat.h:268
const Var kRemID
PID
Definition: Vars.cxx:81
base_types push_back(int_type())
const Var kCCE
Definition: NumuVars.h:21
const Binning kSimpBins
virtual void Go() override
Load all the registered spectra.
loader
Definition: demo0.py:10
std::vector< float > Spectrum
Definition: Constants.h:610
c1 Close()
static bool isFHC
const SystShifts kNoShift
Definition: SystShifts.cxx:22
OStream cout
Definition: OStream.cxx:6
const Cut kHasNeutrino([](const caf::SRProxy *sr){return(sr->mc.nnu!=0);})
Check if MC slice has neutrino information (useful for in-and-out tests)
Definition: TruthCuts.h:61
const Binning kOptEnBins
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
void cc()
Definition: test_ana.C:28
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
void NuMu2019_BasicPIDPlots_ND(bool isFHC, bool isData, std::string period="full")
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 SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
TFile * OutFile
const Cut kNumuQuality
Definition: NumuCuts.h:18
const Var kCVNBestMuonScore([](const caf::SRProxy *sr){float muonScore=-5.0;if(kCVNMuonIdx(sr)< 0.0) return muonScore;if(!sr->vtx.elastic.IsValid) return muonScore;muonScore=sr->vtx.elastic.fuzzyk.png[(unsigned int) kCVNMuonIdx(sr)].cvnpart.muonid;if(sr->vtx.elastic.fuzzyk.png[(unsigned int) kCVNMuonIdx(sr)].len > 500.0) muonScore=1.0;return muonScore;})
: Muon score for best muon prong by CVN score & length
Definition: CVNProngVars.h:27
const Var kCVNm
PID
Definition: Vars.cxx:39
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
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
enum BeamMode string