Analyse_Data2DataComp_HigherEnergyCuts.C
Go to the documentation of this file.
1 /******************************************************
2  * A SCRIPT TO MAKE SPECTRA OF KEY VARIABLES
3  * IN THE NUMU ANALYSIS FOR DATA SETS TAKEN
4  * AT DIFFERENT POINTS IN TIME. OUT FILE
5  * FORMATTED IN SUCH A WAY THAT MAKING RATIOS
6  * IS STRAIGHTFORWARD.
7  * USAGE:
8  * Analyse_Data2DataComp_HigherEnergyCuts.C <bool isFHC [true,falese]>
9  * <std::string name of cut to be applied. [OPTIONS: "BEN", "Rock", "NC", "highE", "NCUncont"]>
10  * <samdef for first period's dataset>
11  * <samdef for second period's dataset>
12  * <Out file that will contain spectra ["/Full/Path/To/OutFile.root"]>
13  ******************************************************/
14 
15 #include "CAFAna/Cuts/Cuts.h"
21 #include "CAFAna/Cuts/SpillCuts.h"
22 #include "CAFAna/Core/Var.h"
24 #include "CAFAna/Vars/HistAxes.h"
28 
29 #include "TFile.h"
30 #include "TSystem.h"
31 
32 
33 using namespace ana;
34 
35 const Cut AssignCut(std::string const & s_CutName)
36 {
37  if(s_CutName=="BEN")
38  {
39  std::cout << "\n\n\nAPPLYING BEN CUT\n\n\n" << std::endl;
40  return kBENKaNumuCut;
41  }
42  else if(s_CutName=="Rock")
43  {
44  std::cout << "\n\n\nAPPLYING Rock CUT\n\n\n" << std::endl;
46  }
47  else if(s_CutName=="NC")
48  {
49  std::cout << "\n\n\nAPPLYING NC CUT\n\n\n" << std::endl;
51  }
52  else if(s_CutName=="NCUncont")
53  {
54  std::cout << "\n\n\nAPPLYING NCUncont CUT\n\n\n" << std::endl;
56  }
57  else if(s_CutName=="highE")
58  {
59  std::cout << "\n\n\nAPPLYING highE CUT\n\n\n" << std::endl;
61  }
62  else
63  {
64  std::cout << "\n\n\nTHE CUT NAME WAS NOT RECOGNISED\n\n\n" << std::endl;
65  std::exit(1);
66  }
67 }
68 
69 
70 void Analyse_Data2DataComp_HigherEnergyCuts(bool isFHC = true, std::string s_CutName = "BEN", std::string data1 = "", std::string data2 = "", std::string s_OutFile = "") {
71 
72  const Cut kChosenCut = AssignCut(s_CutName);
73 
74  //CUTS TO DIVIDE DETECTOR INTO QUADRANTS.
75  const Cut kUpperWest([](const caf::SRProxy* sr){
76  return ( sr->slc.meanpos.x > 0.0 ) && ( sr->slc.meanpos.y > 0.0 );
77  });
78 
79  const Cut kLowerEast([](const caf::SRProxy* sr){
80  return ( sr->slc.meanpos.x < 0.0 ) && ( sr->slc.meanpos.y < 0.0 );
81  });
82 
83  //DEFINE SOME VARIABLES TO PLOT.
84  const Var kTrkCalE ([](const caf::SRProxy* sr) { return (sr->trk.kalman.ntracks < 1 ? 0. : (float)sr->trk.kalman.tracks[0].calE); });
85  const Var kTrkCalEPerNHit([](const caf::SRProxy* sr){
86  float TrkCalE = (sr->trk.kalman.ntracks < 1 ? 0. : (float)sr->trk.kalman.tracks[0].calE);
87  float MuNHit = (sr->trk.kalman.ntracks < 1 ? 0. : (float)sr->trk.kalman.tracks[0].nhit);
88  return (MuNHit <= 0 ? 0. : TrkCalE/MuNHit);
89  });
90 
91  //DEFINE ANY SPECIALISED BINNING.
92  const Binning kBinningHadNHit = Binning::Custom({0, 2, 4, 6, 8, 10, 12, 14, 18, 22, 26, 30, 34, 38, 48, 60, 72, 85, 100, 120, 140, 160, 180, 200, 220, 250});
93 
94  std::vector<SpectrumLoader*> loaders;
95  loaders.push_back(new SpectrumLoader(data1));
96  loaders.push_back(new SpectrumLoader(data2));
97 
98  //GET HADE QUANTILES.
99  std::string sFHC_RHC = "rhc";
100  if(isFHC) sFHC_RHC = "fhc";
101 
102  std::string fdspecfile = (std::string)std::getenv("NUMUDATA_DIR")+"ana2018/Quantiles/quantiles__"+sFHC_RHC+"_full__numu2018.root";
103 
104  std::cout << "\n\nACCESSING QUANTILE BOUNDARY INFORMATION FROM FILE:\n" << fdspecfile << "\n\n\n";
105 
106  TFile* f_In = TFile::Open( fdspecfile.c_str() );
107  TH2 *FDSpec2D = (TH2*)f_In->FindObjectAny( "FDSpec2D" );
108  const int NHadEFracQuantiles = 4;
109  std::vector<Cut> QuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, NHadEFracQuantiles);
110  QuantCuts.push_back(kNoCut);
111 
112  //DEFINE ARRAY OF LOADERS AND SPECTRA.
113  unsigned int nLoader = 2;
114  unsigned int nQuant = 5;
115  std::vector<Spectrum*> mySpecs[nLoader][nQuant];
116 
117  for(unsigned int Q = 0; Q < nQuant; ++Q) {
118  for(unsigned int L = 0; L < nLoader; ++L){
119 
120  // DEFINE ALL SPECTRA HERE:
121  mySpecs[L][Q].push_back(new Spectrum("Nu E (all slices)", Binning::Simple(50, 0, 15), *loaders[L], kCCE, kChosenCut && QuantCuts[Q]));
122  mySpecs[L][Q].push_back(new Spectrum("Nu E (upper-West)", Binning::Simple(50, 0, 15), *loaders[L], kCCE, kChosenCut && QuantCuts[Q] && kUpperWest));
123  mySpecs[L][Q].push_back(new Spectrum("Nu E (lower-East)", Binning::Simple(50, 0, 15), *loaders[L], kCCE, kChosenCut && QuantCuts[Q] && kLowerEast));
124 
125  mySpecs[L][Q].push_back(new Spectrum("NHit in Slice (all slices)", Binning::Simple(50, 0, 600), *loaders[L], kNHit, kChosenCut && QuantCuts[Q]));
126  mySpecs[L][Q].push_back(new Spectrum("NHit in Slice (upper-West)", Binning::Simple(50, 0, 600), *loaders[L], kNHit, kChosenCut && QuantCuts[Q] && kUpperWest));
127  mySpecs[L][Q].push_back(new Spectrum("NHit in Slice (lower-East)", Binning::Simple(50, 0, 600), *loaders[L], kNHit, kChosenCut && QuantCuts[Q] && kLowerEast));
128 
129  mySpecs[L][Q].push_back(new Spectrum("Mu Trk Length, (m) (all slices)", Binning::Simple(100, 0, 30), *loaders[L], kTrkLength, kChosenCut && QuantCuts[Q]));
130  mySpecs[L][Q].push_back(new Spectrum("Mu Trk Length, (m) (upper-West)", Binning::Simple(100, 0, 30), *loaders[L], kTrkLength, kChosenCut && QuantCuts[Q] && kUpperWest));
131  mySpecs[L][Q].push_back(new Spectrum("Mu Trk Length, (m) (lower-East)", Binning::Simple(100, 0, 30), *loaders[L], kTrkLength, kChosenCut && QuantCuts[Q] && kLowerEast));
132 
133  mySpecs[L][Q].push_back(new Spectrum("Mu RecoE, (GeV) (all slices)", Binning::Simple(100, 0, 10.), *loaders[L], kCCE - kHadE, kChosenCut && QuantCuts[Q]));
134  mySpecs[L][Q].push_back(new Spectrum("Mu RecoE, (GeV) (upper-West)", Binning::Simple(100, 0, 10.), *loaders[L], kCCE - kHadE, kChosenCut && QuantCuts[Q] && kUpperWest));
135  mySpecs[L][Q].push_back(new Spectrum("Mu RecoE, (GeV) (lower-East)", Binning::Simple(100, 0, 10.), *loaders[L], kCCE - kHadE, kChosenCut && QuantCuts[Q] && kLowerEast));
136 
137  mySpecs[L][Q].push_back(new Spectrum("Mu CalE, (GeV) (all slices)", Binning::Simple(100, 0, 10), *loaders[L], kTrkCalE, kChosenCut && QuantCuts[Q]));
138  mySpecs[L][Q].push_back(new Spectrum("Mu CalE, (GeV) (upper-West)", Binning::Simple(100, 0, 10), *loaders[L], kTrkCalE, kChosenCut && QuantCuts[Q] && kUpperWest));
139  mySpecs[L][Q].push_back(new Spectrum("Mu CalE, (GeV) (lower-East)", Binning::Simple(100, 0, 10), *loaders[L], kTrkCalE, kChosenCut && QuantCuts[Q] && kLowerEast));
140 
141  mySpecs[L][Q].push_back(new Spectrum("Mu NHit (all slices)", Binning::Simple(50, 0, 500), *loaders[L], kTrkNhits, kChosenCut && QuantCuts[Q]));
142  mySpecs[L][Q].push_back(new Spectrum("Mu NHit (upper-West)", Binning::Simple(50, 0, 500), *loaders[L], kTrkNhits, kChosenCut && QuantCuts[Q] && kUpperWest));
143  mySpecs[L][Q].push_back(new Spectrum("Mu NHit (lower-East)", Binning::Simple(50, 0, 500), *loaders[L], kTrkNhits, kChosenCut && QuantCuts[Q] && kLowerEast));
144 
145  mySpecs[L][Q].push_back(new Spectrum("Muon CalE / Hit, (GeV) (all slices)", Binning::Simple(50, 0., 0.03), *loaders[L], kTrkCalEPerNHit, kChosenCut && QuantCuts[Q]));
146  mySpecs[L][Q].push_back(new Spectrum("Muon CalE / Hit, (GeV) (upper-West)", Binning::Simple(50, 0., 0.03), *loaders[L], kTrkCalEPerNHit, kChosenCut && QuantCuts[Q] && kUpperWest));
147  mySpecs[L][Q].push_back(new Spectrum("Muon CalE / Hit, (GeV) (lower-East)", Binning::Simple(50, 0., 0.03), *loaders[L], kTrkCalEPerNHit, kChosenCut && QuantCuts[Q] && kLowerEast));
148 
149  mySpecs[L][Q].push_back(new Spectrum("Had RecoE, (GeV) (all slices)", Binning::Simple(100, 0, 8), *loaders[L], kHadE, kChosenCut && QuantCuts[Q]));
150  mySpecs[L][Q].push_back(new Spectrum("Had RecoE, (GeV) (upper-West)", Binning::Simple(100, 0, 8), *loaders[L], kHadE, kChosenCut && QuantCuts[Q] && kUpperWest));
151  mySpecs[L][Q].push_back(new Spectrum("Had RecoE, (GeV) (lower-East)", Binning::Simple(100, 0, 8), *loaders[L], kHadE, kChosenCut && QuantCuts[Q] && kLowerEast));
152 
153  mySpecs[L][Q].push_back(new Spectrum("Had CalE, (GeV) (all slices)", Binning::Simple(100, 0, 8), *loaders[L], kNumuHadCalE, kChosenCut && QuantCuts[Q]));
154  mySpecs[L][Q].push_back(new Spectrum("Had CalE, (GeV) (upper-West)", Binning::Simple(100, 0, 8), *loaders[L], kNumuHadCalE, kChosenCut && QuantCuts[Q] && kUpperWest));
155  mySpecs[L][Q].push_back(new Spectrum("Had CalE, (GeV) (lower-East)", Binning::Simple(100, 0, 8), *loaders[L], kNumuHadCalE, kChosenCut && QuantCuts[Q] && kLowerEast));
156 
157  mySpecs[L][Q].push_back(new Spectrum("Had NHit (all slices)", kBinningHadNHit, *loaders[L], kHadNHit, kChosenCut && QuantCuts[Q]));
158  mySpecs[L][Q].push_back(new Spectrum("Had NHit (upper-West)", kBinningHadNHit, *loaders[L], kHadNHit, kChosenCut && QuantCuts[Q] && kUpperWest));
159  mySpecs[L][Q].push_back(new Spectrum("Had NHit (lower-East)", kBinningHadNHit, *loaders[L], kHadNHit, kChosenCut && QuantCuts[Q] && kLowerEast));
160 
161  mySpecs[L][Q].push_back(new Spectrum("Had CalE / Hit, (GeV) (all slices)", Binning::Simple(50, 0., 0.03), *loaders[L], kHadEPerNHit, kChosenCut && QuantCuts[Q]));
162  mySpecs[L][Q].push_back(new Spectrum("Had CalE / Hit, (GeV) (upper-West)", Binning::Simple(50, 0., 0.03), *loaders[L], kHadEPerNHit, kChosenCut && QuantCuts[Q] && kUpperWest));
163  mySpecs[L][Q].push_back(new Spectrum("Had CalE / Hit, (GeV) (lower-East)", Binning::Simple(50, 0., 0.03), *loaders[L], kHadEPerNHit, kChosenCut && QuantCuts[Q] && kLowerEast));
164 
165  mySpecs[L][Q].push_back(new Spectrum("CVN Numu 2017 Score (all slices)", Binning::Simple(25, 0, 1.01), *loaders[L], kCVNm2017, kChosenCut && QuantCuts[Q]));
166  mySpecs[L][Q].push_back(new Spectrum("CVN Numu 2017 Score (upper-West)", Binning::Simple(25, 0, 1.01), *loaders[L], kCVNm2017, kChosenCut && QuantCuts[Q] && kUpperWest));
167  mySpecs[L][Q].push_back(new Spectrum("CVN Numu 2017 Score (lower-East)", Binning::Simple(25, 0, 1.01), *loaders[L], kCVNm2017, kChosenCut && QuantCuts[Q] && kLowerEast));
168 
169  mySpecs[L][Q].push_back(new Spectrum("CVN Numu 2018 Score (all slices)", Binning::Simple(25, 0, 1.01), *loaders[L], kCVNm, kChosenCut && QuantCuts[Q]));
170  mySpecs[L][Q].push_back(new Spectrum("CVN Numu 2018 Score (upper-West)", Binning::Simple(25, 0, 1.01), *loaders[L], kCVNm, kChosenCut && QuantCuts[Q] && kUpperWest));
171  mySpecs[L][Q].push_back(new Spectrum("CVN Numu 2018 Score (lower-East)", Binning::Simple(25, 0, 1.01), *loaders[L], kCVNm, kChosenCut && QuantCuts[Q] && kLowerEast));
172 
173  mySpecs[L][Q].push_back(new Spectrum("Remid (all slices)", kRemidBinning, *loaders[L], kRemID, kChosenCut && QuantCuts[Q]));
174  mySpecs[L][Q].push_back(new Spectrum("Remid (upper-West)", kRemidBinning, *loaders[L], kRemID, kChosenCut && QuantCuts[Q] && kUpperWest));
175  mySpecs[L][Q].push_back(new Spectrum("Remid (lower-East)", kRemidBinning, *loaders[L], kRemID, kChosenCut && QuantCuts[Q] && kLowerEast));
176  }
177  }
178 
179  loaders[0]->SetSpillCut(kStandardSpillCuts);
180  loaders[0]->Go();
181  loaders[1]->SetSpillCut(kStandardSpillCuts);
182  loaders[1]->Go();
183 
184  //LOOP OVER ARRAY AND WRITE OUT SPECTRA TO FILE.
185  char dirName[256];
186  TFile *outFile = new TFile(s_OutFile.c_str(), "RECREATE");
187 
188  for(unsigned int L = 0; L < nLoader; ++L){
189  for(unsigned int Q = 0; Q < nQuant; ++Q) {
190  for(unsigned int S = 0; S < mySpecs[L][Q].size(); ++S) {
191 
192  sprintf(dirName, "dir_Var%d_Q%d_DataSet%d", S, Q, L);
193  mySpecs[L][Q][S]->SaveTo(outFile, dirName);
194  }
195  }
196  }
197 
198  std::cout << "\n\n\nALL DONE!\n\n\n";
199 
200 }
const Var kHadE
Definition: NumuVars.h:23
const Var kHadNHit([](const caf::SRProxy *sr){unsigned int nought=0;if(sr->trk.kalman.ntracks< 1) return nought;return sr->slc.nhit-sr->trk.kalman.tracks[0].nhit;})
Definition: NumuVars.h:61
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 Cut kBENKaNumuCut
Definition: BeamNueCuts.h:18
#define S(x, n)
const Binning kRemidBinning
Binning for plotting remid attractively.
Definition: Binning.cxx:80
caf::Proxy< size_t > ntracks
Definition: SRProxy.h:1756
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
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2108
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
const Cut AssignCut(std::string const &s_CutName)
const Var kTrkCalEPerNHit([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return 0.f;return sr->trk.kalman.tracks[0].calE/sr->trk.kalman.tracks[0].nhit;})
Definition: NumuVars.h:68
const Var kTrkLength([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-5.f;return sr->trk.kalman.tracks[0].len/100;})
Definition: NumuVars.h:65
const Var kHadEPerNHit([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return 0.0f;int nHit=sr->slc.nhit-sr->trk.kalman.tracks[0].nhit;if(nHit<=0) return 0.0f;float hadE=sr->energy.numu.hadcalE;return hadE/nHit;})
Definition: NumuVars.h:63
const Var kNumuHadCalE
Definition: NumuVars.cxx:538
TFile * outFile
Definition: PlotXSec.C:135
const Var kRemID
PID
Definition: Vars.cxx:81
caf::Proxy< caf::SRTrackBranch > trk
Definition: SRProxy.h:2127
static constexpr double L
caf::Proxy< float > x
Definition: SRProxy.h:105
std::string getenv(std::string const &name)
const Var kNHit
Definition: Vars.cxx:71
const Var kCCE
Definition: NumuVars.h:21
const Var kTrkCalE([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return 0.f;return float(sr->trk.kalman.tracks[0].calE);})
Definition: NumuVars.h:67
caf::StandardRecord * sr
static Binning Custom(const std::vector< double > &edges)
Definition: Binning.cxx:161
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
static bool isFHC
OStream cout
Definition: OStream.cxx:6
caf::Proxy< caf::SRKalman > kalman
Definition: SRProxy.h:1775
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void Analyse_Data2DataComp_HigherEnergyCuts(bool isFHC=true, std::string s_CutName="BEN", std::string data1="", std::string data2="", std::string s_OutFile="")
std::string dirName
Definition: PlotSpectra.h:47
const Cut kNumuBasicQuality([](const caf::SRProxy *sr){return(sr->energy.numu.trkccE > 0 && sr->sel.remid.pid > 0 && sr->slc.nhit > 20 && sr->slc.ncontplanes > 4 && sr->trk.cosmic.ntracks > 0);})
Definition: NumuCuts.h:14
caf::Proxy< caf::SRSlice > slc
Definition: SRProxy.h:2124
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
caf::Proxy< float > y
Definition: SRProxy.h:106
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 int nQuant
Definition: varsandcuts.h:4
caf::Proxy< caf::SRVector3D > meanpos
Definition: SRProxy.h:1289
std::vector< Loaders * > loaders
Definition: syst_header.h:386
exit(0)
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...
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
caf::Proxy< std::vector< caf::SRKalmanTrack > > tracks
Definition: SRProxy.h:1758
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:114
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
const Var kTrkNhits([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return 65535;return int(sr->trk.kalman.tracks[0].nhit);})
Definition: NumuVars.h:59