purity.C
Go to the documentation of this file.
1 // Author: Prabhjot Singh (prabhjot@fnal.gov; prabhjot.singh@qmul.ac.uk)
2 // Date: 04 March 2021
3 // Purpose: To make plots of purity for the numubar CC low hadronic activity analysis
4 //
5 // Command: cafe -b -q -l 1 --numubarccinc numubarcc_low_hadronic/Analysis/Preselection/PurityEfficiency/purity.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/Core/Ratio.h"
15 #include "CAFAna/Core/SystShifts.h"
16 #include "CAFAna/Vars/Vars.h"
17 #include "CAFAna/Vars/XsecTunes.h"
19 #include "CAFAna/Cuts/TruthCuts.h"
21 
22 // NDAna includes
25 
26 // Standard Record
28 
29 // ROOT includes
30 #include "TH1.h"
31 #include "TFile.h"
32 
33 using namespace ana;
34 
35 void purity()
36 {
37 
38  // Prod 5.1. Flatsumdecaf RHC ND MC dataset. This dataset has 400 files
39  const std::string dataset_name = "prod_flatsumdecaf_development_nd_genie_N1810j0211a_nonswap_rhc_nova_v08_full_ndphysics_contain_v1";
40 
41  // Start spectrum loader
42  ana::SpectrumLoader loader(dataset_name);
43 
44  // Binning: let's use some simple binning
45  const ana::Binning bins = Binning::Simple(100, -1.0, 1.0);
46 
47  // Muon ID variable to plot
49 
50  // Signal and background definitions including the true fiducial volumes
51  const ana::Cut kTrueSignal = xsec::numubarcc::kIsNumubarCC && xsec::numubarcc::kTrueVtxCut; // True numubar CC signal
52  const ana::Cut kTrueWrongSign = xsec::numubarcc::kIsNumuCC && xsec::numubarcc::kTrueVtxCut; // True numu CC for the wrong sign contamination
53  const ana::Cut kTrueNC = xsec::numubarcc::kIsNC && xsec::numubarcc::kTrueVtxCut; // True NC
54  const ana::Cut kTrueBkgd = xsec::numubarcc::kIsNueorbarCC && xsec::numubarcc::kTrueVtxCut; // True NuE or the true NuE bar
55 
56  // Selections
57  // Qualty, Containment and Fiducial Selections
58  const ana::Cut kSelections = ana::xsec::numubarcc::kQualityCut &&
61 
62  // XSec and PPFX weights
64  const ana::Var xsec_wgt = ana::VarFromNuTruthVar(xsec_wgt_ST, 1);
66  const ana::Var flux_wgt = ana::VarFromNuTruthVar(flux_wgt_ST, 1);
67  const ana::NuTruthVar std_wgt_ST = flux_wgt_ST*xsec_wgt_ST;
68  const ana::Var std_wgt = ana::VarFromNuTruthVar(std_wgt_ST, 1);
69 
70  // Make a spectrum
71  Spectrum sTrueSignal ("MuonID", bins, loader, kMuonID, kTrueSignal && kSelections, ana::kNoShift, std_wgt);
72  Spectrum sTrueAllEvents ("MuonID", bins, loader, kMuonID, (kTrueSignal || kTrueWrongSign || kTrueNC || kTrueBkgd) && kSelections, ana::kNoShift, std_wgt);
73 
74  // Go ahead with spectrum loader
75  loader.Go();
76 
77  // POT for scaling histograms
78  const double kPOTRHC = ana::kAna2020RHCPOT; // = 12.5003e20 RHC POT
79 
80  //Histograms
81  TH1D* h_TrueSignal = (TH1D*)sTrueSignal .ToTH1(kPOTRHC, 1 );
82  TH1D* h_AllEvents = (TH1D*)sTrueAllEvents .ToTH1(kPOTRHC, 2 );
83  TH1D* h_SignalPurity = (TH1D*)h_TrueSignal ->Clone("h_SignalPurity");
84  h_SignalPurity ->Divide(h_AllEvents);
85  h_SignalPurity ->SetLineColor(8);
86 
87  // Set Names of the histograms
88  h_TrueSignal ->SetName("TrueNuMubarCCSignal" );
89  h_AllEvents ->SetName("TrueAllEvents" );
90  h_SignalPurity ->SetName("NuMuBarCCSignalPurity" );
91 
92  // Make a root file and save histogram in it
93  TFile* f_output = new TFile("purity.root", "recreate");
94 
95  // Save histograms
96  h_TrueSignal ->Write();
97  h_AllEvents ->Write();
98  h_SignalPurity ->Write();
99 
100  // Close output root file
101  f_output ->Close();
102 
103 }// 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'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
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);})
virtual void Go() override
Load all the registered spectra.
void purity()
Definition: purity.C:35
loader
Definition: demo0.py:10
const double kAna2020RHCPOT
Definition: Exposures.h:235
const SystShifts kNoShift
Definition: SystShifts.cxx:22
const ana::NuTruthVar xsec_wgt_ST
Definition: muonid_opt.C:24
const Binning bins
Definition: NumuCC_CPiBin.h:8
Var VarFromNuTruthVar(const NuTruthVar &stv, double _default)
Definition: Var.cxx:8
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 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
const ana::NuTruthVar std_wgt_ST
Definition: muonid_opt.C:28
enum BeamMode string