get_fd_dataspectrum.C
Go to the documentation of this file.
2 #include "CAFAna/Core/Spectrum.h"
3 
4 #include "CAFAna/Cuts/Cuts.h"
9 
10 #include "CAFAna/Vars/Vars.h"
13 
14 #include "TFile.h"
15 
16 
17 using namespace ana;
18 
19 
20 void get_fd_dataspectrum(std::string polarity = "rhc",
21  std::string period = "full")
22 {
23 
24  std::cout << "get fd data spectrum" << std::endl;
25  std::cout << "polarity: " << polarity << ", period: " << period << std::endl;
26 
27  // read and get quantile cuts
28  std::cout << "\ngeting quantile cuts from file" << std::endl;
29  std::string outdir = "";
30  std::string indir_quant = "/nova/ana/nu_mu_ana/Ana2018/Quantiles";
31  std::string infile_quant = indir_quant + "/quantiles__" + polarity + "_full__numu2018.root";
32  TFile* infile = TFile::Open( pnfs2xrootd(infile_quant).c_str() );
33  TH2* FDSpec2D = (TH2*)infile->FindObjectAny("FDSpec2D");
34  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, 4);
35  infile->Close();
36  delete infile;
37 
38 
39  // input fd data
40  std::cout << "\ngeting fd numi files" << std::endl;
41  std::string fd_numi = "prod4_sumdecaf_fd_numi_" + polarity + "_" + period + "_goodruns_numu2018";
42  //std::string fd_numi = "prod4_sumrestricteddecaf_fd_numi_" + polarity + "_" + period + "_goodruns_numu2018";
43  SpectrumLoader loader(fd_numi);
45 
46 
47  std::vector<Spectrum*> s_data;
48  s_data.push_back(new Spectrum(loader, kNumuCCOptimisedAxis, kNumuCutFD2018 && kInBeamSpill));
49  for(auto thisCut : HadEFracQuantCuts){
50  s_data.push_back(new Spectrum(loader, kNumuCCOptimisedAxis, kNumuCutFD2018 && kInBeamSpill && thisCut));
51  } // loop over quantiles
52  loader.Go();
53 
54 
55  // save the fd data spectrums
56  std::string outname = outdir + "fd_dataspectrum_" + polarity + "_" + period + "__numu2018.root";
57  std::cerr << "all quantiles to file " << outname << std::endl;
58  TFile *fout = TFile::Open(outname.c_str(), "RECREATE");
59 
60  const int totquant = 4;
61  for(int quant=0; quant<totquant+1; quant++){
62  std::cout << "saving quantile " << quant << " spectrum" << std::endl;
63  auto fd_data = s_data[quant];
64  fd_data-> SaveTo(fout, Form("fd_data_q%d", quant));
65  }
66 
67  fout->Close();
68  delete fout;
69 
70 
71 } // end get_fd_dataspectrum
72 
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
int totquant
Definition: saveFDMCHists.C:40
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
void SetSpillCut(const SpillCut &cut)
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: UtilsExt.cxx:237
string infile
void get_fd_dataspectrum(std::string polarity="rhc", std::string period="full")
virtual void Go() override
Load all the registered spectra.
loader
Definition: demo0.py:10
std::vector< float > Spectrum
Definition: Constants.h:610
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.
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 Cut kNumuCutFD2018
Definition: NumuCuts2018.h:39
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
const std::string outdir
enum BeamMode string