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
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:40
const double kAna2018RHCPOT
Definition: Exposures.h:208
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic >::Index size_type
Definition: typedefs.hpp:11
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.
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
const SystShifts kNoShift
Definition: SystShifts.cxx:21
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:49
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