getData.C
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////////
2 // NOvA Muon Neutrino Appearance Analysis 2017 //
3 // //
4 // Author: Diana Patricia Mendez d.p.mendez@sussex.ac.uk //
5 // //
6 //////////////////////////////////////////////////////////////////////////
7 // //
8 // getData.C //
9 // Gets FD data spectra for each quantile //
10 // //
11 //////////////////////////////////////////////////////////////////////////
12 
14 #include "CAFAna/Core/Spectrum.h"
15 
16 #include "CAFAna/Cuts/Cuts.h"
17 #include "CAFAna/Cuts/SpillCuts.h"
18 #include "CAFAna/Cuts/TimingCuts.h"
21 
22 #include "CAFAna/Vars/Vars.h"
23 #include "CAFAna/Vars/HistAxes.h"
27 
28 #include "TFile.h"
29 #include "varsandcuts.h"
30 
31 using namespace ana;
32 
33 void getData(int sampleCut = 9, std::string period = "full", std::string anaCut = "3A", bool energy3a = true, bool optBin = true){
34 
35 
36 
37  // sample
38  std::string SampleName = "";
39  std::cout << "\nSample:" << std::endl;
40  SpillCut runSpillCut = (const_cast<SpillCut&>(kStandardSpillCuts));
41 
42  if(sampleCut == 9){
43  std::cout << "9e20 POT" << std::endl;
44  SampleName = "9e20";
45  }
46  if(sampleCut == 6){
47  std::cout << "6e20 POT" << std::endl;
48  SampleName = "6e20";
49  runSpillCut = (const_cast<SpillCut&>(kIs6e20));
50  }
51  if(sampleCut == 3){
52  std::cout << "3e20 POT" << std::endl;
53  SampleName = "3e20";
54  runSpillCut = (const_cast<SpillCut&>(kIs3e20));
55  }
56 
57  const SpillCut kSampleSpillCuts = (const SpillCut)runSpillCut;
58 
59 
60  std::string CutName = "";
61  std::cout << "Cut:" << std::endl;
62  Cut NDCutAna = (const_cast<Cut&>(kNumuND));
63  Cut FDCutAna = (const_cast<Cut&>(kNumuFD));
64 
65  if(anaCut == "SA"){
66  std::cout << "SA\n" << std::endl;
67  CutName = "SACut";
68  NDCutAna = (const_cast<Cut&>(kNumuND_Prod3));
69  FDCutAna = (const_cast<Cut&>(kNumuFD_Prod3));
70  }
71  if(anaCut == "3A"){
72  std::cout<< "3A\n" << std::endl;
73  CutName = "3ACut";
74  NDCutAna = (const_cast<Cut&>(kNumuCutND2017));
75  FDCutAna = (const_cast<Cut&>(kNumuCutFD2017));
76  }
77  const Cut kNDCut_Here = (const Cut)NDCutAna;
78  const Cut kFDCut_Here = (const Cut)FDCutAna;
79 
80 
81 
82 
83  Var EnergyAna = (const_cast<Var&>(kCCE));
84  std::string EnergyName = "";
85  std::cout << "Energy Estimator:" << std::endl;
86  if(!energy3a){
87  std::cout << "SA" << std::endl;
88  EnergyName = "SAEnergy";
89  EnergyAna = (const_cast<Var&>(kSAEnergy));
90  }
91  if(energy3a){
92  std::cout << "3A" << std::endl;
93  EnergyName = "3AEnergy";
94  EnergyAna = (const_cast<Var&>(k3AEnergy));
95  }
96  const Var kNumuEnergy_Here = (const Var)EnergyAna;
97 
98 
99 
100 
101  // binning
102  Binning BinAna = (const_cast<Binning&>(kNumuEnergyBinning));
103  std::string BinName = "SimpleBinning";
104  std::string BinNameDir = "SimpleBin";
105  std::cout << "Binning:" << std::endl;
106  if(!optBin){
107  std::cout << "Simple binning" << std::endl;
108  }
109  if(optBin){
110  std::cout << "Optimised binning" << std::endl;
111  BinName = "OptBinning";
112  BinNameDir = "OptBin";
113  BinAna = (const_cast<Binning&>(kNumuCCEOptimisedBinning));
114  }
115  const Binning kBinning_Here = (const Binning)BinAna;
116  const HistAxis kNumuCCAxis_Here( "Reconstructed Neutrino Energy (GeV)", kBinning_Here, kNumuEnergy_Here );
117 
118 
119 
120 
121  // quantiles
122  //Read and get quantile cuts
123  std::cout << "Get quantile cuts from file" << std::endl;
124  std::string FDSpecFile = "/nova/ana/nu_mu_ana/Ana2017/Quantiles/all_FD_histo_for_quantile_cuts.root";
125  TFile inFile(FDSpecFile.c_str());
126  TH2 *FDSpec2D = (TH2*)inFile.FindObjectAny("FDSpec2D");
127  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, 4);
128 
129  // input fd numi data
130  std::string fd_numi = "prod_sumrestricteddecaf_R17-03-01-prod3reco.k_fd_numi_fhc_" + period + "_v1_goodruns_numu2017";
131  SpectrumLoader loader(fd_numi);
132  loader.SetSpillCut(kStandardSpillCuts && kSampleSpillCuts);
133 
135  // make a vector with an entry for each quantiles
136  std::vector<Spectrum*> sDataVec;
137  for(auto thisCut : HadEFracQuantCuts){
138  sDataVec.push_back(new Spectrum(loader, kNumuCCOptimisedAxis, kNumuCutFD2017 && kInBeamSpill && thisCut));
139  } // loop over quantiles
140  loader.Go();
141 
142 
143  // Now we should save the predictions
144  std::cerr << "Saving Spectra" << std::endl;
145  std::string OutFileName = "fd_data_3ACut_3AEnergy_OptBinning_full_" + SampleName + "_Quant0.root";
146  TFile *fAQ = TFile::Open(OutFileName.c_str(), "RECREATE");
147  sData->SaveTo(fAQ, "fd_data");
148  fAQ->Close();
149  delete fAQ;
150 
151  for(size_t quant=1; quant <= sDataVec.size(); quant++){
152  auto fd_data = sDataVec[quant-1];
153  std::string OutFileName = "fd_data_3ACut_3AEnergy_OptBinning_full_" + SampleName + "_" + (std::string)Form("Quant%d", quant) + ".root";
154  std::cerr << "Saving Quantile " << quant << " to Spectrum file: " << OutFileName << std::endl;
155  TFile *f = TFile::Open(OutFileName.c_str(), "RECREATE");
156  fd_data -> SaveTo(f, "fd_data");
157  f->Close();
158  delete f;
159  }
160 
161 
162 }
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 getData(int sampleCut=9, std::string period="full", std::string anaCut="3A", bool energy3a=true, bool optBin=true)
Definition: getData.C:33
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
const Cut kNumuCutFD2017
Definition: NumuCuts2017.h:39
const SpillCut kIs6e20([](const caf::SRSpillProxy *s){int RunNumber=s->run;int lastFDSARun=22900;int lastNDSARun=11383;if(s->det==caf::kFARDET) return(RunNumber<=lastFDSARun);if(s->det==caf::kNEARDET) return(RunNumber<=lastNDSARun);})
const Var k3AEnergy
Definition: varsandcuts.h:16
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
OStream cerr
Definition: OStream.cxx:7
const Cut kNumuFD
Definition: NumuCuts.h:53
void SetSpillCut(const SpillCut &cut)
const Var kSAEnergy([](const caf::SRProxy *sr){return sr->energy.numu.E;})
_Var< caf::SRProxy > Var
Representation of a variable to be retrieved from a caf::StandardRecord object.
Definition: Var.h:74
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
ifstream inFile
Definition: AnaPlotMaker.h:34
const SpillCut kIs3e20([](const caf::SRSpillProxy *s){int RunNumber=s->run;int lastFDSARun=22900;int lastNDSARun=11383;if(s->det==caf::kFARDET) return(RunNumber >lastFDSARun);if(s->det==caf::kNEARDET) return(RunNumber >lastNDSARun);})
const Binning kNumuEnergyBinning
Definition: Binnings.cxx:13
const Cut kNumuND
Definition: NumuCuts.h:55
const Var kCCE
Definition: NumuVars.h:21
_Cut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
Definition: Cut.h:96
virtual void Go() override
Load all the registered spectra.
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:506
const Cut kNumuCutND2017
Definition: NumuCuts2017.h:41
loader
Definition: demo0.py:10
std::vector< float > Spectrum
Definition: Constants.h:610
_Cut< caf::SRSpillProxy > SpillCut
Equivalent of Cut acting on caf::SRSpill. For use in spill-by-spill data quality cuts.
Definition: Cut.h:100
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
OStream cout
Definition: OStream.cxx:6
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
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 Cut kNumuND_Prod3
Definition: varsandcuts.h:35
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 Cut kInBeamSpill([](const caf::SRProxy *sr){if(sr->spill.run > util::kLastBadTimingRun) return kInBeamSpill_main(sr);else return kInBeamSpill_main(sr)||kInBeamSpill_shifted(sr);}, [](const caf::SRSpillProxy *spill){if(spill->run > util::kLastBadTimingRun) return kInBeamSpill_main.Livetime(spill);else return kInBeamSpill_main.Livetime(spill)+kInBeamSpill_shifted.Livetime(spill);}, [](const caf::SRSpillProxy *spill) -> double{return spill->spillpot;})
Does the event fall inside the window we call the beam spill?
Definition: TimingCuts.h:8
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
Template for Cut and SpillCut.
Definition: Cut.h:15
const Cut kNumuFD_Prod3
Definition: varsandcuts.h:36
enum BeamMode string