getTimePeak.C
Go to the documentation of this file.
1 #include "CAFAna/Cuts/Cuts.h"
7 #include "CAFAna/Core/Loaders.h"
9 #include "CAFAna/Core/Spectrum.h"
11 #include "CAFAna/Vars/Vars.h"
15 #include "OscLib/OscCalcPMNSOpt.h"
16 #include "CAFAna/Analysis/Calcs.h"
17 #include "CAFAna/Analysis/Plots.h"
18 
29 #include "CAFAna/Vars/FitVars.h"
30 
31 #include "TStyle.h"
32 #include "TCanvas.h"
33 #include "TGraph.h"
34 #include "TGraphAsymmErrors.h"
35 #include "TH1.h"
36 #include "TH2.h"
37 #include "TLegend.h"
38 #include "TLatex.h"
39 
40 using namespace ana;
41 
42 const int kNumQuantiles = 4;
43 
44 void Preliminary() // For some reason it doesn't include it automatically
45 {
46  TLatex* prelim = new TLatex(.9, .95, "NOvA Preliminary");
47  prelim->SetTextColor(kBlue);
48  prelim->SetNDC();
49  prelim->SetTextSize(2/30.);
50  prelim->SetTextAlign(32);
51  prelim->Draw();
52 }
53 
54 void Exposure()
55 {
56  TLatex* ExposureLabel = new TLatex(0.55,0.80,"6.05#times10^{20} POT-equiv." );
57 
58  ExposureLabel->SetNDC();
59  ExposureLabel->SetTextSize(0.05);
60  ExposureLabel->Draw();
61 
62 }
63 
64 const Var kSliceTimeShift(
65  [](const caf::SRProxy *sr)
66  {
67  const double t = sr->slc.meantime / 1000;
68  if(sr->hdr.ismc)
69  return t;
70  const int run = sr->hdr.run;
72  return t;
73  if (!util::IsInBeamWindow(run, t*1000) )
74  return t;
75  if (t < 250.)
76  return t;
77  else
78  return (t-64.);
79  });
80 
81 struct Selection
82 {
86 };
87 
89 {
90  std::cout << "Drawing plots..." << std::endl;
91 
92  std::string fd_dir = "/pnfs/nova/persistent/users/bzamoran/CheckEvents/";
93 
94  // std::string nd_data = "prod_sumdecaf_R17-03-01-prod3reco.d_nd_numi_fhc_period3_v1_goodruns_numu2017";
95  // std::string nd_nonswap = "prod_sumdecaf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_period3_v1_numu2017";
96  std::string fd_nonswap = fd_dir + "nonswap.root";
97  std::string fd_fluxswap = fd_dir + "fluxswap.root";
98  std::string fd_tau = fd_dir + "tau.root";
99 
102 
103  // loaders.SetLoaderPath( nd_data, caf::kNEARDET, Loaders::kData, ana::kBeam,
104  // Loaders::kNonSwap );
105  // loaders.SetLoaderPath( nd_nonswap, caf::kNEARDET, Loaders::kMC, ana::kBeam,
106  // Loaders::kNonSwap );
107  loaders.SetLoaderPath( fd_nonswap, caf::kFARDET, Loaders::kMC, ana::kBeam,
109  loaders.SetLoaderPath( fd_fluxswap, caf::kFARDET, Loaders::kMC, ana::kBeam,
113 
114  gStyle->SetStatFormat("6.3g"); // Significant figures. Fix if this is not satisfactory
115  gStyle->SetPaintTextFormat("6.3g");
116  gStyle->SetFitFormat("6.3g");
117 
120  calc.SetDmsq32(2.42993e-3);
121  calc.SetTh23(asin(sqrt(0.476033))); // preliminary Ana2017
122 
123  HistAxis kAxis( "#Deltat from t_{0} (#mus)", Binning::Simple(35,25,445), kSliceTimeShift );
124 
125  // Get FD files for HadEFrac vs. neutrino energy distribution
126  std::string quantpath = "/pnfs/nova/persistent/analysis/numu/Ana2017/all_FD_histo_for_quantile_cuts.root";
127  TFile* quantfile = TFile::Open(pnfs2xrootd(quantpath).c_str());
128  TH2* quanthist = (TH2*)quantfile->Get("dir_FDSpec2D/FDSpec2D");
129  std::vector<Cut> QuantCuts = QuantileCutsFromTH2(quanthist, kNumuCCOptimisedAxis, kHadEFracAxis, 4);
130 
131  const int kNumCuts = 5;
132 
133  Selection cuts[kNumCuts] = {
134  {"Quartile 1", "Q1", QuantCuts[0]},
135  {"Quartile 2", "Q2", QuantCuts[1]},
136  {"Quartile 3", "Q3", QuantCuts[2]},
137  {"Quartile 4", "Q4", QuantCuts[3]},
138  {"All quartiles", "AllQ", kNoCut}
139  };
140 
141  // make a vector with an entry for each quantiles
142  // std::vector<NumuExtrapGenerator*> NumuCCExtrapVec;
143  // std::vector<PredictionInterp*> predictionVec;
144  std::vector<PredictionNoExtrap*> predictionVec;
145 
146  for(int k = 0; k < kNumCuts; k++)
147  {
148  // NumuCCExtrapVec.push_back(new NumuExtrapGenerator(kAxis, cuts[k].cut && kNumuCutFD2017, cuts[k].cut && kNumuCutND2017, kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017));
149  // predictionVec.push_back(new PredictionInterp({}, &calc, *(NumuCCExtrapVec.back()), loaders));
150  predictionVec.push_back( new PredictionNoExtrap(loaders, kAxis, kNumuCutFD2017 && cuts[k].cut, kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017) );
151 
152  } // loop over quantiles
153 
154  loaders.Go();
155 
156  std::string OutFileName = "";
157  // save predictions for later
158  for(int k = 0; k < kNumCuts; k++)
159  {
160  auto prediction = predictionVec[k];
161  OutFileName = "Predictions/pred_" + cuts[k].fname + ".root";
162  std::cerr << "Saving cut " << (cuts[k].label).c_str() << " to file: " << OutFileName << std::endl;
163  TFile *f = TFile::Open(OutFileName.c_str(), "RECREATE");
164  prediction -> SaveTo(f, "prediction");
165  f->Close();
166  delete f;
167  }
168 
169  // now actually do the plots!
170 
171  std::cout << "All plots finished!" << std::endl;
172 
173 
174 } // End of function
Far Detector at Ash River.
Definition: SREnums.h:11
void Exposure()
Definition: getTimePeak.C:54
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Cut kNumuCutFD2017
Definition: NumuCuts2017.h:39
std::string fname
Definition: getTimePeak.C:84
const int kLastBadTimingRun
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2049
static bool IsInBeamWindow(const int run, const double time)
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2038
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
T sqrt(T number)
Definition: d0nt_math.hpp:156
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
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: Utilities.cxx:620
osc::OscCalcDumb calc
caf::Proxy< unsigned int > run
Definition: SRProxy.h:243
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
const Var kSliceTimeShift([](const caf::SRProxy *sr){const double t=sr->slc.meantime/1000;if(sr->hdr.ismc) return t;const int run=sr->hdr.run;if(run > util::kLastBadTimingRun) return t;if(!util::IsInBeamWindow(run, t *1000)) return t;if(t< 250.) return t;else return(t-64.);})
void SetSpillCut(const SpillCut &cut)
Definition: Loaders.cxx:121
if(dump)
std::string label
Definition: getTimePeak.C:83
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
const int kNumQuantiles
Definition: getTimePeak.C:42
void SetTh23(const T &th23) override
TLatex * prelim
Definition: Xsec_final.C:133
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
Definition: run.py:1
const SystShifts kNoShift
Definition: SystShifts.h:115
OStream cout
Definition: OStream.cxx:6
void SetDmsq32(const T &dmsq32) override
void Preliminary()
Definition: getTimePeak.C:44
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const Cut cut
Definition: exporter_fd.C:30
Template for Cut and SpillCut.
Definition: Cut.h:15
caf::Proxy< caf::SRSlice > slc
Definition: SRProxy.h:2054
std::vector< Loaders * > loaders
Definition: syst_header.h:386
caf::Proxy< float > meantime
Definition: SRProxy.h:1258
caf::Proxy< bool > ismc
Definition: SRProxy.h:237
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...
void getPredictions()
Definition: getTimePeak.C:88
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:46
Prediction that just uses FD MC, with no extrapolation.
const Var kXSecCVWgt2017
Definition: XsecTunes.h:37
Float_t e
Definition: plot.C:35
void SetLoaderPath(const std::string &path, caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Configure loader via wildcard path.
Definition: Loaders.cxx:25
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:38
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:117
static constexpr Double_t sr
Definition: Munits.h:164
T asin(T number)
Definition: d0nt_math.hpp:60
Template for Var and SpillVar.
Definition: Var.h:16