Analyse_Data2DataComp.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.C <bool isFHC [true,falese]>
9  * <samdef for first period's dataset>
10  * <samdef for second period's dataset>
11  * <Out file that will contain spectra ["/Full/Path/To/OutFile.root"]>
12  ******************************************************/
13 
14 #include "CAFAna/Cuts/Cuts.h"
19 #include "CAFAna/Cuts/SpillCuts.h"
20 #include "CAFAna/Core/Var.h"
22 #include "CAFAna/Vars/HistAxes.h"
26 
27 #include "TFile.h"
28 #include "TSystem.h"
29 
30 
31 using namespace ana;
32 
33 void Analyse_Data2DataComp(bool isFHC = true, std::string data1 = "", std::string data2 = "", std::string s_OutFile = "") {
34 
35  //CUTS TO DIVIDE DETECTOR INTO QUADRANTS.
36  const Cut kUpperWest([](const caf::SRProxy* sr){
37  return ( sr->slc.meanpos.x > 0.0 ) && ( sr->slc.meanpos.y > 0.0 );
38  });
39 
40  const Cut kLowerEast([](const caf::SRProxy* sr){
41  return ( sr->slc.meanpos.x < 0.0 ) && ( sr->slc.meanpos.y < 0.0 );
42  });
43 
44  //DEFINE SOME VARIABLES TO PLOT.
45  const Var kTrkCalE ([](const caf::SRProxy* sr) { return (sr->trk.kalman.ntracks < 1 ? 0. : (float)sr->trk.kalman.tracks[0].calE); });
46  const Var kTrkCalEPerNHit([](const caf::SRProxy* sr){
47  float TrkCalE = (sr->trk.kalman.ntracks < 1 ? 0. : (float)sr->trk.kalman.tracks[0].calE);
48  float MuNHit = (sr->trk.kalman.ntracks < 1 ? 0. : (float)sr->trk.kalman.tracks[0].nhit);
49  return (MuNHit <= 0 ? 0. : TrkCalE/MuNHit);
50  });
51 
52  //DEFINE ANY SPECIALISED BINNING.
53  const Binning kBinningHadNHit = Binning::Custom({0, 2, 4, 6, 8, 10, 12, 14, 18, 22, 26, 30, 34, 38, 45, 55, 70, 85, 100});
54 
55  std::vector<SpectrumLoader*> loaders;
56  loaders.push_back(new SpectrumLoader(data1));
57  loaders.push_back(new SpectrumLoader(data2));
58 
59  //GET HADE QUANTILES.
60  std::string sFHC_RHC = "rhc";
61  if(isFHC) sFHC_RHC = "fhc";
62 
63  std::string fdspecfile = (std::string)std::getenv("NUMUDATA_DIR")+"ana2018/Quantiles/quantiles__"+sFHC_RHC+"_full__numu2018.root";
64 
65  std::cout << "\n\nACCESSING QUANTILE BOUNDARY INFORMATION FROM FILE:\n" << fdspecfile << "\n\n\n";
66 
67  TFile* f_In = TFile::Open( fdspecfile.c_str() );
68  TH2 *FDSpec2D = (TH2*)f_In->FindObjectAny( "FDSpec2D" );
69  const int NHadEFracQuantiles = 4;
70  std::vector<Cut> QuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, NHadEFracQuantiles);
71  QuantCuts.push_back(kNoCut);
72 
73  //DEFINE ARRAY OF LOADERS AND SPECTRA.
74  unsigned int nLoader = 2;
75  unsigned int nQuant = 5;
76  std::vector<Spectrum*> mySpecs[nLoader][nQuant];
77 
78  for(unsigned int Q = 0; Q < nQuant; ++Q) {
79  for(unsigned int L = 0; L < nLoader; ++L){
80 
81  // DEFINE ALL SPECTRA HERE:
82  mySpecs[L][Q].push_back(new Spectrum("Nu E (all slices)", kNumuCCEOptimisedBinning, *loaders[L], kCCE, kNumuCutND2018 && QuantCuts[Q]));
83  mySpecs[L][Q].push_back(new Spectrum("Nu E (upper-West)", kNumuCCEOptimisedBinning, *loaders[L], kCCE, kNumuCutND2018 && QuantCuts[Q] && kUpperWest));
84  mySpecs[L][Q].push_back(new Spectrum("Nu E (lower-East)", kNumuCCEOptimisedBinning, *loaders[L], kCCE, kNumuCutND2018 && QuantCuts[Q] && kLowerEast));
85 
86  mySpecs[L][Q].push_back(new Spectrum("NHit in Slice (all slices)", Binning::Simple(25, 0, 300), *loaders[L], kNHit, kNumuCutND2018 && QuantCuts[Q]));
87  mySpecs[L][Q].push_back(new Spectrum("NHit in Slice (upper-West)", Binning::Simple(25, 0, 300), *loaders[L], kNHit, kNumuCutND2018 && QuantCuts[Q] && kUpperWest));
88  mySpecs[L][Q].push_back(new Spectrum("NHit in Slice (lower-East)", Binning::Simple(25, 0, 300), *loaders[L], kNHit, kNumuCutND2018 && QuantCuts[Q] && kLowerEast));
89 
90  mySpecs[L][Q].push_back(new Spectrum("Mu Trk Length, (m) (all slices)", Binning::Simple(50, 0, 16), *loaders[L], kTrkLength, kNumuCutND2018 && QuantCuts[Q]));
91  mySpecs[L][Q].push_back(new Spectrum("Mu Trk Length, (m) (upper-West)", Binning::Simple(50, 0, 16), *loaders[L], kTrkLength, kNumuCutND2018 && QuantCuts[Q] && kUpperWest));
92  mySpecs[L][Q].push_back(new Spectrum("Mu Trk Length, (m) (lower-East)", Binning::Simple(50, 0, 16), *loaders[L], kTrkLength, kNumuCutND2018 && QuantCuts[Q] && kLowerEast));
93 
94  mySpecs[L][Q].push_back(new Spectrum("Mu RecoE, (GeV) (all slices)", Binning::Simple(50, 0, 4.5), *loaders[L], kCCE - kHadE, kNumuCutND2018 && QuantCuts[Q]));
95  mySpecs[L][Q].push_back(new Spectrum("Mu RecoE, (GeV) (upper-West)", Binning::Simple(50, 0, 4.5), *loaders[L], kCCE - kHadE, kNumuCutND2018 && QuantCuts[Q] && kUpperWest));
96  mySpecs[L][Q].push_back(new Spectrum("Mu RecoE, (GeV) (lower-East)", Binning::Simple(50, 0, 4.5), *loaders[L], kCCE - kHadE, kNumuCutND2018 && QuantCuts[Q] && kLowerEast));
97 
98  mySpecs[L][Q].push_back(new Spectrum("Mu CalE, (GeV) (all slices)", Binning::Simple(25, 0, 4), *loaders[L], kTrkCalE, kNumuCutND2018 && QuantCuts[Q]));
99  mySpecs[L][Q].push_back(new Spectrum("Mu CalE, (GeV) (upper-West)", Binning::Simple(25, 0, 4), *loaders[L], kTrkCalE, kNumuCutND2018 && QuantCuts[Q] && kUpperWest));
100  mySpecs[L][Q].push_back(new Spectrum("Mu CalE, (GeV) (lower-East)", Binning::Simple(25, 0, 4), *loaders[L], kTrkCalE, kNumuCutND2018 && QuantCuts[Q] && kLowerEast));
101 
102  mySpecs[L][Q].push_back(new Spectrum("Mu NHit (all slices)", Binning::Simple(25, 0, 220), *loaders[L], kTrkNhits, kNumuCutND2018 && QuantCuts[Q]));
103  mySpecs[L][Q].push_back(new Spectrum("Mu NHit (upper-West)", Binning::Simple(25, 0, 220), *loaders[L], kTrkNhits, kNumuCutND2018 && QuantCuts[Q] && kUpperWest));
104  mySpecs[L][Q].push_back(new Spectrum("Mu NHit (lower-East)", Binning::Simple(25, 0, 220), *loaders[L], kTrkNhits, kNumuCutND2018 && QuantCuts[Q] && kLowerEast));
105 
106  mySpecs[L][Q].push_back(new Spectrum("Muon CalE / Hit, (GeV) (all slices)", Binning::Simple(25, 0., 0.03), *loaders[L], kTrkCalEPerNHit, kNumuCutND2018 && QuantCuts[Q]));
107  mySpecs[L][Q].push_back(new Spectrum("Muon CalE / Hit, (GeV) (upper-West)", Binning::Simple(25, 0., 0.03), *loaders[L], kTrkCalEPerNHit, kNumuCutND2018 && QuantCuts[Q] && kUpperWest));
108  mySpecs[L][Q].push_back(new Spectrum("Muon CalE / Hit, (GeV) (lower-East)", Binning::Simple(25, 0., 0.03), *loaders[L], kTrkCalEPerNHit, kNumuCutND2018 && QuantCuts[Q] && kLowerEast));
109 
110  mySpecs[L][Q].push_back(new Spectrum("Had RecoE, (GeV) (all slices)", Binning::Simple(25, 0, 3), *loaders[L], kHadE, kNumuCutND2018 && QuantCuts[Q]));
111  mySpecs[L][Q].push_back(new Spectrum("Had RecoE, (GeV) (upper-West)", Binning::Simple(50, 0, 3), *loaders[L], kHadE, kNumuCutND2018 && QuantCuts[Q] && kUpperWest));
112  mySpecs[L][Q].push_back(new Spectrum("Had RecoE, (GeV) (lower-East)", Binning::Simple(50, 0, 3), *loaders[L], kHadE, kNumuCutND2018 && QuantCuts[Q] && kLowerEast));
113 
114  mySpecs[L][Q].push_back(new Spectrum("Had CalE, (GeV) (all slices)", Binning::Simple(50, 0, 2), *loaders[L], kNumuHadCalE, kNumuCutND2018 && QuantCuts[Q]));
115  mySpecs[L][Q].push_back(new Spectrum("Had CalE, (GeV) (upper-West)", Binning::Simple(50, 0, 2), *loaders[L], kNumuHadCalE, kNumuCutND2018 && QuantCuts[Q] && kUpperWest));
116  mySpecs[L][Q].push_back(new Spectrum("Had CalE, (GeV) (lower-East)", Binning::Simple(50, 0, 2), *loaders[L], kNumuHadCalE, kNumuCutND2018 && QuantCuts[Q] && kLowerEast));
117 
118  mySpecs[L][Q].push_back(new Spectrum("Had NHit (all slices)", kBinningHadNHit, *loaders[L], kHadNHit, kNumuCutND2018 && QuantCuts[Q]));
119  mySpecs[L][Q].push_back(new Spectrum("Had NHit (upper-West)", kBinningHadNHit, *loaders[L], kHadNHit, kNumuCutND2018 && QuantCuts[Q] && kUpperWest));
120  mySpecs[L][Q].push_back(new Spectrum("Had NHit (lower-East)", kBinningHadNHit, *loaders[L], kHadNHit, kNumuCutND2018 && QuantCuts[Q] && kLowerEast));
121 
122  mySpecs[L][Q].push_back(new Spectrum("Had CalE / Hit, (GeV) (all slices)", Binning::Simple(50, 0., 0.03), *loaders[L], kHadEPerNHit, kNumuCutND2018 && QuantCuts[Q]));
123  mySpecs[L][Q].push_back(new Spectrum("Had CalE / Hit, (GeV) (upper-West)", Binning::Simple(50, 0., 0.03), *loaders[L], kHadEPerNHit, kNumuCutND2018 && QuantCuts[Q] && kUpperWest));
124  mySpecs[L][Q].push_back(new Spectrum("Had CalE / Hit, (GeV) (lower-East)", Binning::Simple(50, 0., 0.03), *loaders[L], kHadEPerNHit, kNumuCutND2018 && QuantCuts[Q] && kLowerEast));
125 
126  mySpecs[L][Q].push_back(new Spectrum("CVN Numu 2017 Score (all slices)", Binning::Simple(25, 0, 1.01), *loaders[L], kCVNm2017, kNumuCutND2018 && QuantCuts[Q]));
127  mySpecs[L][Q].push_back(new Spectrum("CVN Numu 2017 Score (upper-West)", Binning::Simple(25, 0, 1.01), *loaders[L], kCVNm2017, kNumuCutND2018 && QuantCuts[Q] && kUpperWest));
128  mySpecs[L][Q].push_back(new Spectrum("CVN Numu 2017 Score (lower-East)", Binning::Simple(25, 0, 1.01), *loaders[L], kCVNm2017, kNumuCutND2018 && QuantCuts[Q] && kLowerEast));
129 
130  mySpecs[L][Q].push_back(new Spectrum("CVN Numu 2018 Score (all slices)", Binning::Simple(25, 0, 1.01), *loaders[L], kCVNm, kNumuCutND2018 && QuantCuts[Q]));
131  mySpecs[L][Q].push_back(new Spectrum("CVN Numu 2018 Score (upper-West)", Binning::Simple(25, 0, 1.01), *loaders[L], kCVNm, kNumuCutND2018 && QuantCuts[Q] && kUpperWest));
132  mySpecs[L][Q].push_back(new Spectrum("CVN Numu 2018 Score (lower-East)", Binning::Simple(25, 0, 1.01), *loaders[L], kCVNm, kNumuCutND2018 && QuantCuts[Q] && kLowerEast));
133 
134  mySpecs[L][Q].push_back(new Spectrum("Remid (all slices)", kRemidBinning, *loaders[L], kRemID, kNumuCutND2018 && QuantCuts[Q]));
135  mySpecs[L][Q].push_back(new Spectrum("Remid (upper-West)", kRemidBinning, *loaders[L], kRemID, kNumuCutND2018 && QuantCuts[Q] && kUpperWest));
136  mySpecs[L][Q].push_back(new Spectrum("Remid (lower-East)", kRemidBinning, *loaders[L], kRemID, kNumuCutND2018 && QuantCuts[Q] && kLowerEast));
137  }
138  }
139 
140  loaders[0]->SetSpillCut(kStandardSpillCuts);
141  loaders[0]->Go();
142  loaders[1]->SetSpillCut(kStandardSpillCuts);
143  loaders[1]->Go();
144 
145  //LOOP OVER ARRAY AND WRITE OUT SPECTRA TO FILE.
146  char dirName[256];
147  TFile *outFile = new TFile(s_OutFile.c_str(), "RECREATE");
148 
149  for(unsigned int L = 0; L < nLoader; ++L){
150  for(unsigned int Q = 0; Q < nQuant; ++Q) {
151  for(unsigned int S = 0; S < mySpecs[L][Q].size(); ++S) {
152 
153  sprintf(dirName, "dir_Var%d_Q%d_DataSet%d", S, Q, L);
154  mySpecs[L][Q][S]->SaveTo(outFile, dirName);
155  }
156  }
157  }
158 
159  std::cout << "\n\n\nALL DONE!\n\n\n";
160 
161 }
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
#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 kNumuCutND2018
Definition: NumuCuts2018.h:41
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
void Analyse_Data2DataComp(bool isFHC=true, std::string data1="", std::string data2="", std::string s_OutFile="")
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:573
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
std::string dirName
Definition: PlotSpectra.h:47
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 int nQuant
Definition: varsandcuts.h:4
caf::Proxy< caf::SRVector3D > meanpos
Definition: SRProxy.h:1289
std::vector< Loaders * > loaders
Definition: syst_header.h:386
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
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