energyResolution.C
Go to the documentation of this file.
1 // to get histograms of cosmics from cosmic trigger data and numi stream, and GENIE prediction, by period and by quantile
2 // argument is period: 1, 2, 3, or 5
3 // author: Kirk Bays (bays@caltech.edu) Sep 7 2017
4 // Histograms made: cosrej BDT inputs, neutrino E, and optimized neutrino E
5 // suitable to run on the grid (set bool grid to true to keep track of pot/livetime and not scale per job)
6 // another important option is 'quantiletype': 0 means use quantiles determined per period, 1 means use quantiles determined from all periods combined
7 
10 #include "CAFAna/Core/Var.h"
11 #include "CAFAna/Analysis/Plots.h"
12 
13 #include <cmath>
14 #include <vector>
15 #include <sstream>
16 #include <string>
17 #include <iomanip>
18 
19 
20 using namespace ana;
21 
22 void energyResolution(bool fhc = true, bool isFirstRun = false)
23 {
24 
25  gStyle->SetStatFormat("6.3g");
26  gStyle->SetPaintTextFormat("6.3g");
27  gStyle->SetFitFormat("6.3g");
28  gStyle->SetTitleOffset(1.15,"x");
29  gStyle->SetTitleOffset(0.95,"y");
30  gStyle->SetPadLeftMargin(0.12);
31  gStyle->SetPadBottomMargin(0.15);
32  gStyle->SetHistLineWidth(3);
33 
34  double livetime, pot;
35  std::string numifiles, cosfiles, geniefiles, outfile, quantpath, horn, prodgen;
36  if(fhc){
37  livetime = kAna2018FHCLivetime;
38  pot = kAna2018FHCPOT;
39  horn = "fhc";
40  prodgen = ".d";
41  }
42  if(!fhc){
43  livetime = kAna2018RHCLivetime;
44  pot = kAna2018RHCPOT;
45  horn = "rhc";
46  prodgen = ".e";
47  }
48 
49  geniefiles = "prod_sumdecaf_R17-11-14-prod4reco" + prodgen + "_fd_genie_nonswap_" + horn + "_nova_v08_full_v1_numu2018";
50  outfile = horn + "_energyResolution_numu2018.root";
51 
52  const Var kRMTOT([](const caf::SRProxy *sr)
53  {
54  if(sr->mc.nnu == 0) return -10.0;
55  return (kCCE(sr)-kTrueE(sr))/kTrueE(sr);
56  });
57 
58  std::cout << "\n\nGet quantile cuts from file" << std::endl;
59  quantpath = "/nova/ana/nu_mu_ana/Ana2018/Quantiles/quantiles__" + horn + "_full__numu2018.root";
60  TFile* quantfile = TFile::Open(pnfs2xrootd(quantpath).c_str());
61  TH2* quanthist = (TH2*)quantfile->FindObjectAny("FDSpec2D");
62  std::vector<Cut> QuantCuts = QuantileCutsFromTH2(quanthist, kNumuCCOptimisedAxis, kHadEFracAxis, 4);
63  quantfile->Close();
64  delete quantfile;
65 
66 
67  if(isFirstRun){
68  SpectrumLoader loaderGenie(geniefiles);
69  loaderGenie.SetSpillCut(kStandardSpillCuts);
70 
71  HistAxis axisRMTOT("(E_{reco}-E_{true})/E_{true}", Binning::Simple(200,-1,1), kRMTOT);
72 
73  Spectrum* specTotal = new Spectrum(loaderGenie, axisRMTOT, kNumuCutFD2018,
75  std::map< TString, Spectrum* > specMap;
76  // for(auto thisCut : QuantCuts){
77  for(std::vector<Cut>::size_type i = 0; i != QuantCuts.size(); i++) {
78  std::cout<< "i: " << i <<std::endl;
79  specMap[Form("rmtot_q%d",(i+1))] = new Spectrum(loaderGenie, axisRMTOT, QuantCuts[i] && kNumuCutFD2018,
80  kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2018);
81  }
82 
83 
84  loaderGenie.Go();
85 
86 
87  TFile fout(outfile.c_str(), "recreate");
88  TH1D* histTotal = specTotal -> ToTH1(pot);
89  histTotal -> Write("rmtot_comb");
90  for(std::vector<Cut>::size_type i = 0; i != QuantCuts.size(); i++) {
91  std::cout<< "i: " << i <<std::endl;
92  TH1D* hist = specMap[Form("rmtot_q%d",(i+1))] -> ToTH1(pot);
93  hist -> Write(Form("rmtot_q%d",(i+1)));
94  }
95  fout.Close();
96 
97  }
98 
99 
100  //////////////////////////////
101  // Draw hists
102  DontAddDirectory guard;
103  TFile fout(outfile.c_str());
104 
105  TH1D* histTotal = (TH1D*)fout.Get("rmtot_comb");
106 
107  TCanvas* c = new TCanvas();
108  TLegend *leg = new TLegend(0.15,0.45,0.48,0.85);
109  histTotal -> GetYaxis() -> SetRangeUser(0, 1.3 * histTotal->GetMaximum());
110  histTotal -> Draw("hist");
111  histTotal -> GetXaxis() -> CenterTitle();
112  double rms = histTotal -> GetRMS();
113  double mean = histTotal -> GetMean();
114  leg -> AddEntry(histTotal, "Combined", "l");
115  leg -> AddEntry(histTotal, Form("Mean:%.03lf, RMS:%.03lf",mean,rms), "l");
116 
117  for(std::vector<Cut>::size_type i = 0; i != QuantCuts.size(); i++) {
118  TH1D* hist = (TH1D*) fout.Get(Form("rmtot_q%d",(i+1)));
119  hist -> SetLineColor(i+2);
120  if((i+2) == 5) hist -> SetLineColor(6);
121  hist -> Draw("histsame");
122  double rms = hist -> GetRMS();
123  double mean = hist -> GetMean();
124  leg -> AddEntry(hist, Form("Q%d",(i+1)), "l");
125  leg -> AddEntry(hist, Form("Mean:%.03lf, RMS:%.03lf",mean,rms), "l");
126  }
127  leg -> SetFillStyle(0);
128  leg -> SetFillColor(0);
129  leg -> Draw();
130  c -> SaveAs((horn+"_energyResSpec.pdf").c_str());
131 
132  fout.Close();
133 
134 
135 
136 }
tree Draw("slc.nhit")
bin1_2sigma SetFillColor(3)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
correl_xv GetYaxis() -> SetDecimals()
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
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
std::vector< double > Spectrum
Definition: Constants.h:746
caf::Proxy< short int > nnu
Definition: SRProxy.h:617
void SetSpillCut(const SpillCut &cut)
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: UtilsExt.cxx:237
ntuple SetFillStyle(1001)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:39
const double kAna2018RHCPOT
Definition: Exposures.h:208
const Var kTrueE([](const caf::SRProxy *sr){assert(sr->mc.nnu==1);return sr->mc.nu[0].E;})
Definition: Vars.cxx:85
correl_xv GetXaxis() -> SetDecimals()
leg AddEntry(GRdata,"data","p")
TH1D * histTotal
Definition: PlotSpectra.h:53
hmean SetLineColor(4)
const Var kCCE
Definition: NumuVars.h:21
#define pot
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
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
const SystShifts kNoShift
Definition: SystShifts.cxx:22
OStream cout
Definition: OStream.cxx:6
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
double livetime
Definition: saveFDMCHists.C:21
const Var kXSecCVWgt2018
Definition: XsecTunes.h:48
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 kNumuCutFD2018
Definition: NumuCuts2018.h:39
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
const double kAna2018FHCPOT
Definition: Exposures.h:207
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
cosmicTree SaveAs("cosmicTree.root")
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
const double kAna2018FHCLivetime
Definition: Exposures.h:213
FILE * outfile
Definition: dump_event.C:13
void energyResolution(bool fhc=true, bool isFirstRun=false)
const double kAna2018RHCLivetime
Definition: Exposures.h:214
gm Write()
enum BeamMode string