muonID.C
Go to the documentation of this file.
1 // Author: Prabhjot Singh (prabhjot@fnal.gov; prabhjot.singh@qmul.ac.uk)
2 // Date: 26 Feb 2021
3 // Purpose: To make plots of muonID for the numubar CC low hadronic activity analysis
4 //
5 // Command: cafe -b -q -l 1 --numubarccinc numubarcc_low_hadronic/Analysis/Preselection/Analysis_Variables/muonID.C
6 // To run it over all dataset remove -l 1 option
7 //
8 // Includes
9 // CAFAna specific
10 #include "CAFAna/Cuts/Cuts.h"
11 #include "CAFAna/Core/Binning.h"
12 #include "CAFAna/Core/Spectrum.h"
14 #include "CAFAna/Vars/Vars.h"
15 #include "CAFAna/Vars/XsecTunes.h"
17 #include "CAFAna/Cuts/TruthCuts.h"
19 #include "CAFAna/Core/SystShifts.h"
20 
21 // NDAna includes
26 
27 // Standard Record
29 
30 // ROOT includes
31 #include "TCanvas.h"
32 #include "TH1.h"
33 #include "TH2.h"
34 #include "TFile.h"
35 
36 // C++ includes
37 #include <iostream>
38 
39 using namespace ana;
40 
41 void muonID()
42 {
43 
44  // Prod 5.1. Flatsumdecaf RHC ND MC dataset. This dataset has 400 files
45  const std::string dataset_name = "prod_flatsumdecaf_development_nd_genie_N1810j0211a_nonswap_rhc_nova_v08_full_ndphysics_contain_v1";
46 
47  // Start spectrum loader
48  ana::SpectrumLoader loader(dataset_name);
49 
50  // Binning: let's use some simple binning
51  const ana::Binning bins = Binning::Simple(100, -1.0, 1.0);
52  const ana::Binning bins_NewMuonID = Binning::Simple(50, 0.0, 1.0);
53  const ana::Binning bins_prong_length = Binning::Simple(400, 0.0, 16);
54 
55  // variables to plot
57 
59 
60  // Signal and background definitions including the true fiducial volumes
61  const ana::Cut kTrueNumubarCC = xsec::numubarcc::kIsNumubarCC && xsec::numubarcc::kTrueVtxCut; // True numubar CC signal
62  const ana::Cut kTrueNumuCC = xsec::numubarcc::kIsNumuCC && xsec::numubarcc::kTrueVtxCut; // True numu CC for the wrong sign contamination
63  const ana::Cut kTrueNC = xsec::numubarcc::kIsNC && xsec::numubarcc::kTrueVtxCut; // True NC
64  const ana::Cut kTrueNueorNuebarCC = xsec::numubarcc::kIsNueorbarCC && xsec::numubarcc::kTrueVtxCut; // True NuE or the true NuE bar
65 
66  // Selections
67  // Qualty, Containment and Loose Fiducial Selections
68  const ana::Cut kSelections = ana::xsec::numubarcc::kQualityCut &&
71 
72  // XSec and PPFX weights
74  const ana::Var xsec_wgt = ana::VarFromNuTruthVar(xsec_wgt_ST, 1);
76  const ana::Var flux_wgt = ana::VarFromNuTruthVar(flux_wgt_ST, 1);
77  const ana::NuTruthVar std_wgt_ST = flux_wgt_ST*xsec_wgt_ST;
78  const ana::Var std_wgt = ana::VarFromNuTruthVar(std_wgt_ST, 1);
79 
80  // Make a spectrum
81  // 1D spectra
82  Spectrum sTrueNumubarCC ("MuonID", bins, loader, kMuonID, kTrueNumubarCC && kSelections, ana::kNoShift, std_wgt);
83  Spectrum sTrueNumuCC ("MuonID", bins, loader, kMuonID, kTrueNumuCC && kSelections, ana::kNoShift, std_wgt);
84  Spectrum sTrueNC ("MuonID", bins, loader, kMuonID, kTrueNC && kSelections, ana::kNoShift, std_wgt);
85  Spectrum sTrueNueorNuebarCC ("MuonID", bins, loader, kMuonID, kTrueNueorNuebarCC && kSelections, ana::kNoShift, std_wgt);
86 
87  // 2D spectra
88  Spectrum skMuonTrueNuMuBarCC ("Prong Length (m)", "MuonID", loader, bins_prong_length, kTrkLenActandCat, bins, kMuonID, kTrueNumubarCC && kSelections, ana::kNoShift, std_wgt);
89  Spectrum skMuonTrueNuMuCC ("Prong Length (m)", "MuonID", loader, bins_prong_length, kTrkLenActandCat, bins, kMuonID, kTrueNumuCC && kSelections, ana::kNoShift, std_wgt);
90  Spectrum skMuonTrueNC ("Prong Length (m)", "MuonID", loader, bins_prong_length, kTrkLenActandCat, bins, kMuonID, kTrueNC && kSelections, ana::kNoShift, std_wgt);
91  Spectrum skMuonTrueNueorNuebarCC ("Prong Length (m)", "MuonID", loader, bins_prong_length, kTrkLenActandCat, bins, kMuonID, kTrueNueorNuebarCC && kSelections, ana::kNoShift, std_wgt);
92 
93  // Spectra for new MuonID: CVN5LabelMuonID
94  // 1D Spectra
95  Spectrum sTrueNumubarCC_NewMuonID ("CVN5LabelMuonID", bins_NewMuonID, loader, ana::CVN5LabelMuonID, kTrueNumubarCC && kSelections, ana::kNoShift, std_wgt);
96  Spectrum sTrueNumuCC_NewMuonID ("CVN5LabelMuonID", bins_NewMuonID, loader, ana::CVN5LabelMuonID, kTrueNumuCC && kSelections, ana::kNoShift, std_wgt);
97  Spectrum sTrueNC_NewMuonID ("CVN5LabelMuonID", bins_NewMuonID, loader, ana::CVN5LabelMuonID, kTrueNC && kSelections, ana::kNoShift, std_wgt);
98  Spectrum sTrueNueorNuebarCC_NewMuonID ("CVN5LabelMuonID", bins_NewMuonID, loader, ana::CVN5LabelMuonID, kTrueNueorNuebarCC && kSelections, ana::kNoShift, std_wgt);
99 
100  // 2D spectra
101  Spectrum sCVNMuonTrueNuMuBarCC ("Prong Length (m)", "CVN5LabelMuonID", loader, bins_prong_length, kTrkLenActandCat, bins_NewMuonID, ana::CVN5LabelMuonID, kTrueNumubarCC && kSelections, ana::kNoShift, std_wgt);
102  Spectrum sCVNMuonTrueNuMuCC ("Prong Length (m)", "CVN5LabelMuonID", loader, bins_prong_length, kTrkLenActandCat, bins_NewMuonID, ana::CVN5LabelMuonID, kTrueNumuCC && kSelections, ana::kNoShift, std_wgt);
103  Spectrum sCVNMuonTrueNC ("Prong Length (m)", "CVN5LabelMuonID", loader, bins_prong_length, kTrkLenActandCat, bins_NewMuonID, ana::CVN5LabelMuonID, kTrueNC && kSelections, ana::kNoShift, std_wgt);
104  Spectrum sCVNMuonTrueNueorNuebarCC ("Prong Length (m)", "CVN5LabelMuonID", loader, bins_prong_length, kTrkLenActandCat, bins_NewMuonID, ana::CVN5LabelMuonID, kTrueNueorNuebarCC && kSelections, ana::kNoShift, std_wgt);
105 
106  // Go ahead with spectrum loader
107  loader.Go();
108 
109  // POT for scaling histograms
110  const double kPOTRHC = ana::kAna2020RHCPOT; // = 12.5003e20 RHC POT
111 
112  // Histograms scaled to the RHC 2020 POT
113  TH1D* h_TrueNumubarCC = (TH1D*)sTrueNumubarCC .ToTH1(kPOTRHC, 4 );
114  TH1D* h_TrueNumuCC = (TH1D*)sTrueNumuCC .ToTH1(kPOTRHC, 38 );
115  TH1D* h_TrueNC = (TH1D*)sTrueNC .ToTH1(kPOTRHC, 2 );
116  TH1D* h_TrueNueorNuebarCC = (TH1D*)sTrueNueorNuebarCC .ToTH1(kPOTRHC, 3 );
117 
118  TH1D* h_TrueNumubarCC_NewMuonID = (TH1D*)sTrueNumubarCC_NewMuonID .ToTH1(kPOTRHC, 4 );
119  TH1D* h_TrueNumuCC_NewMuonID = (TH1D*)sTrueNumuCC_NewMuonID .ToTH1(kPOTRHC, 38 );
120  TH1D* h_TrueNC_NewMuonID = (TH1D*)sTrueNC_NewMuonID .ToTH1(kPOTRHC, 2 );
121  TH1D* h_TrueNueorNuebarCC_NewMuonID = (TH1D*)sTrueNueorNuebarCC_NewMuonID .ToTH1(kPOTRHC, 3 );
122 
123  TH2D* h_kMuonTrueNuMuBarCC = (TH2D*)skMuonTrueNuMuBarCC .ToTH2(kPOTRHC);
124  TH2D* h_kMuonTrueNuMuCC = (TH2D*)skMuonTrueNuMuCC .ToTH2(kPOTRHC);
125  TH2D* h_kMuonTrueNC = (TH2D*)skMuonTrueNC .ToTH2(kPOTRHC);
126  TH2D* h_kMuonTrueNueorNuebarCC = (TH2D*)skMuonTrueNueorNuebarCC .ToTH2(kPOTRHC);
127 
128  TH2D* h_CVNMuonTrueNuMuBarCC = (TH2D*)sCVNMuonTrueNuMuBarCC .ToTH2(kPOTRHC);
129  TH2D* h_CVNMuonTrueNuMuCC = (TH2D*)sCVNMuonTrueNuMuCC .ToTH2(kPOTRHC);
130  TH2D* h_CVNMuonTrueNC = (TH2D*)sCVNMuonTrueNC .ToTH2(kPOTRHC);
131  TH2D* h_CVNMuonTrueNueorNuebarCC = (TH2D*)sCVNMuonTrueNueorNuebarCC .ToTH2(kPOTRHC);
132 
133  // Set Names of the histograms
134  h_TrueNumubarCC ->SetName("TrueNuMuBarCC" );
135  h_TrueNumuCC ->SetName("TrueNuMurCC" );
136  h_TrueNC ->SetName("TrueNC" );
137  h_TrueNueorNuebarCC ->SetName("TrueNueorNuebarCC" );
138 
139  h_TrueNumubarCC_NewMuonID ->SetName("TrueNuMuBarCC_NewMuonID" );
140  h_TrueNumuCC_NewMuonID ->SetName("TrueNuMurCC_NewMuonID" );
141  h_TrueNC_NewMuonID ->SetName("TrueNC_NewMuonID" );
142  h_TrueNueorNuebarCC_NewMuonID ->SetName("TrueNueorNuebarCC_NewMuonID" );
143 
144  h_kMuonTrueNuMuBarCC ->SetName("ProngLength_MuonID_TrueNuMuBarCC");
145  h_kMuonTrueNuMuCC ->SetName("ProngLength_MuonID_TrueNuMuCC");
146  h_kMuonTrueNC ->SetName("ProngLength_MuonID_TrueNC");
147  h_kMuonTrueNueorNuebarCC ->SetName("ProngLength_MuonID_TrueNueorNuebarCC");
148 
149  h_CVNMuonTrueNuMuBarCC ->SetName("ProngLength_CVNMuonID_TrueNuMuBarCC");
150  h_CVNMuonTrueNuMuCC ->SetName("ProngLength_CVNMuonID_TrueNuMuCC");
151  h_CVNMuonTrueNC ->SetName("ProngLength_CVNMuonID_TrueNC");
152  h_CVNMuonTrueNueorNuebarCC ->SetName("ProngLength_CVNMuonID_TrueNueorNuebarCC");
153 
154  // Can we just get the variable and put it in the tree. Let's try this.
155  //
156 
157  // Make a root file and save histogram in it
158  TFile* f_output = new TFile("muonID.root", "recreate");
159 
160  // Save histograms
161  h_TrueNumubarCC ->Write();
162  h_TrueNumuCC ->Write();
163  h_TrueNC ->Write();
164  h_TrueNueorNuebarCC ->Write();
165 
166  h_TrueNumubarCC_NewMuonID ->Write();
167  h_TrueNumuCC_NewMuonID ->Write();
168  h_TrueNC_NewMuonID ->Write();
169  h_TrueNueorNuebarCC_NewMuonID ->Write();
170 
171  h_kMuonTrueNuMuBarCC ->Write();
172  h_kMuonTrueNuMuCC ->Write();
173  h_kMuonTrueNC ->Write();
174  h_kMuonTrueNueorNuebarCC ->Write();
175 
176  h_CVNMuonTrueNuMuBarCC ->Write();
177  h_CVNMuonTrueNuMuCC ->Write();
178  h_CVNMuonTrueNC ->Write();
179  h_CVNMuonTrueNueorNuebarCC ->Write();
180 
181  // Close output root file
182  f_output ->Close();
183 
184  std::cout << " " << std::endl;
185  std::cout << " Output file is " << f_output ->GetName() << std::endl;
186  std::cout << " " << std::endl;
187 
188 }// end of main function
const ana::NuTruthVar flux_wgt_ST
Definition: muonid_opt.C:26
const ana::Var std_wgt
Definition: muonid_opt.C:29
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
const NuTruthVar kXSecCVWgt2020_NT
Final 2020 xsec tune. See documentation in NOvARwgt.
Definition: XsecTunes.h:104
void muonID()
Definition: muonID.C:41
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:39
const Var kMuonID
Definition: NDXSecMuonPID.h:32
const Cut kQualityCut([](const caf::SRProxy *sr){return(sr->trk.kalman.ntracks > 0 &&sr->slc.nhit > 20 &&sr->slc.ncontplanes > 4);})
Preselection.
const Cut kContainmentCut([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;int ibesttrk=kBestMuonIDIndex(sr);if(ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return false;for(const caf::SRFuzzyKProngProxy &prong:sr->vtx.elastic.fuzzyk.png) if(!VtxInBounds(&prong.shwlid.start, containLow, containHigh)||!VtxInBounds(&prong.shwlid.stop, containLow, containHigh)) return false;if(sr->trk.kalman.ntracks< 1) return false;const unsigned short muon_catcher_edge=1275;for(unsigned int i=0;i< sr->trk.kalman.ntracks;++i){if(int(i)==ibesttrk) continue;if(sr->trk.kalman.tracks[i].start.Z() > muon_catcher_edge|| sr->trk.kalman.tracks[i].stop.Z() > muon_catcher_edge) return false;}const caf::SRKalmanTrackProxy &besttrack=sr->trk.kalman.tracks[ibesttrk];return((besttrack.stop.Z()< muon_catcher_edge||besttrack.trkyposattrans< 55) &&besttrack.trkfwdcellnd > 5 &&besttrack.trkbakcellnd > 10);})
const Var CVN5LabelMuonID([](const caf::SRProxy *sr){if(!(sr->vtx.elastic.IsValid)) return-999.9f;if(sr->vtx.elastic.fuzzyk.npng==0) return-999.9f;float maxscore=-999.9f;for(unsigned int i=0;i< sr->vtx.elastic.fuzzyk.npng;i++) if(maxscore< sr->vtx.elastic.fuzzyk.png[i].spprongcvnpart5label.muonid) maxscore=sr->vtx.elastic.fuzzyk.png[i].spprongcvnpart5label.muonid;return(float(maxscore));})
virtual void Go() override
Load all the registered spectra.
const Var kTrkLenCat([](const caf::SRProxy *sr){int ibesttrk=kBestMuonIDIndex(sr);if(sr->trk.kalman.ntracks< 1||ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return-1000.f;const caf::SRKalmanTrackProxy &bestmuon=sr->trk.kalman.tracks[ibesttrk]; if(bestmuon.leninact > 0 && bestmuon.lenincat<=0) return 0.f; if(bestmuon.leninact > 0 && bestmuon.lenincat > 0) return float(bestmuon.lenincat/100.); if(bestmuon.leninact<=0 && bestmuon.lenincat > 0) return float((bestmuon.leninact/100.) +(bestmuon.lenincat/100.));return-1000.f;})
loader
Definition: demo0.py:10
const double kAna2020RHCPOT
Definition: Exposures.h:235
const SystShifts kNoShift
Definition: SystShifts.cxx:22
OStream cout
Definition: OStream.cxx:6
const ana::NuTruthVar xsec_wgt_ST
Definition: muonid_opt.C:24
const Var kTrkLenActandCat
Definition: NumuCCIncVars.h:77
const Binning bins
Definition: NumuCC_CPiBin.h:8
Var VarFromNuTruthVar(const NuTruthVar &stv, double _default)
Definition: Var.cxx:8
const Var kTrkLenAct([](const caf::SRProxy *sr){int ibesttrk=kBestMuonIDIndex(sr);if(sr->trk.kalman.ntracks< 1||ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return-1000.f;const caf::SRKalmanTrackProxy &bestmuon=sr->trk.kalman.tracks[ibesttrk]; if(bestmuon.leninact > 0 && bestmuon.lenincat<=0) return float((bestmuon.leninact/100.) +(bestmuon.lenincat/100.)); if(bestmuon.leninact > 0 && bestmuon.lenincat > 0) return float(bestmuon.leninact/100.); if(bestmuon.leninact<=0 && bestmuon.lenincat > 0) return 0.f;return-1000.f;})
float MuonID(const caf::SRProxy *sr)
Definition: MuonID.cxx:74
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const ana::Var xsec_wgt
Definition: muonid_opt.C:25
const Cut kTrueNumuCC
Definition: OverlayCuts.cxx:65
const ana::Var flux_wgt
Definition: muonid_opt.C:27
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
const NuTruthVar kPPFXFluxCVWgtST([](const caf::SRNeutrinoProxy *nu){ if(nu->rwgt.ppfx.cv!=nu->rwgt.ppfx.cv){return 1.f;}if(nu->rwgt.ppfx.cv >90){return 1.f;}return float(nu->rwgt.ppfx.cv);})
weight events with the flux PPFX Central value correction.
Definition: PPFXWeights.h:12
TH2 * ToTH2(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, ana::EBinType bintype)
For use with Var2D.
Definition: UtilsExt.cxx:115
const ana::NuTruthVar std_wgt_ST
Definition: muonid_opt.C:28
enum BeamMode string