FD_plots.C
Go to the documentation of this file.
1 // Fills spectra for systematic error band
2 
3 #include "CAFAna/Cuts/Cuts.h"
7 #include "CAFAna/Core/Binning.h"
8 #include "CAFAna/Core/Loaders.h"
9 #include "CAFAna/Core/Spectrum.h"
10 #include "CAFAna/Core/Var.h"
13 #include "CAFAna/Systs/Systs.h"
19 #include "CAFAna/Cuts/SpillCuts.h"
20 #include "CAFAna/Analysis/Calcs.h"
22 #include "CAFAna/Analysis/Plots.h"
23 #include "CAFAna/Analysis/Style.h"
24 #include "CAFAna/Vars/HistAxes.h"
25 
26 #include "OscLib/OscCalcPMNSOpt.h"
27 
28 #include "TCanvas.h"
29 #include "TFile.h"
30 #include "TGraph.h"
31 #include "TGraphAsymmErrors.h"
32 #include "TH1.h"
33 #include "TH2.h"
34 #include "TMath.h"
35 #include "TGaxis.h"
36 #include "TMultiGraph.h"
37 #include "TLegend.h"
38 #include "TLatex.h"
39 #include "TStyle.h"
40 #include "THStack.h"
41 #include "TPaveText.h"
42 #include "TList.h"
43 #include "TGaxis.h"
44 #include "TAttLine.h"
45 #include "TAttMarker.h"
46 
47 #include <cmath>
48 #include <iostream>
49 #include <vector>
50 #include <list>
51 #include <sstream>
52 #include <string>
53 #include <sstream>
54 #include <fstream>
55 #include <iomanip>
56 
57 using namespace ana;
58 
60 {
61  TGraph* gr1 = new TGraph();
62  gr1->SetPoint(0, 1.5, -1e90);
63  gr1->SetPoint(1, 1.5, 1e90);
64 
65  TGraph* gr2 = new TGraph();
66  gr2->SetPoint(0, 1.75, -1e90);
67  gr2->SetPoint(1, 1.75, 1e90);
68 
69  gr1->SetLineColor(kMagenta);
70  gr2->SetLineColor(kMagenta);
71 
72  gr1->Draw("L");
73  gr2->Draw("L");
74 
75 }
76 
78 {
79  TGraph* gr1 = new TGraph();
80  gr1->SetPoint(0, -1e90, 1.5);
81  gr1->SetPoint(1, 1e90, 1.5);
82 
83  TGraph* gr2 = new TGraph();
84  gr2->SetPoint(0, -1e90, 1.75);
85  gr2->SetPoint(1, 1e90, 1.75);
86 
87  gr1->SetLineColor(kMagenta);
88  gr2->SetLineColor(kMagenta);
89 
90  gr1->Draw("L");
91  gr2->Draw("L");
92 
93 }
94 
95 void FD_plots(int JustOne = -1)
96 {
97 
98  std::string fd_dir = "/pnfs/nova/persistent/users/bzamoran/CheckEvents/";
99 
100  std::string fd_nonswap = fd_dir + "nonswap.root";
101  std::string fd_fluxswap = fd_dir + "fluxswap.root";
102  std::string fd_tau = fd_dir + "tau.root";
103 
105 
106  loaders.SetLoaderPath( fd_nonswap, caf::kFARDET, Loaders::kMC, ana::kBeam,
108  loaders.SetLoaderPath( fd_fluxswap, caf::kFARDET, Loaders::kMC, ana::kBeam,
112 
114 
117  calc.SetDmsq32(2.43071e-3);
118  calc.SetTh23(asin(sqrt(0.471974)));
119  calc.SetdCP(M_PI*0.915869); // preliminary Ana2017
120 
121  struct Plot
122  {
125  Binning bins;
126  Var var;
127  };
128 
129  // And some binning definitions
130  const Binning kXYBins = Binning::Simple(36,-9,9);
131  const Binning kZBins = Binning::Simple(30,0,60);
132  const Binning kEnergyBinning = Binning::Simple(50,0,5);
133 
134  // Vars
135  const Var kVisHadE(
136  [](const caf::SRProxy* sr)
137  {
138  return sr->energy.numu.hadtrkE + sr->energy.numu.hadcalE;
139  }
140  );
141 
142  const Var kSAMuE(
143  [](const caf::SRProxy* sr)
144  {
145  float muE;
146  if(sr->trk.kalman.ntracks < 1) muE = -5.f;
147  else muE = SAMuE(sr->trk.kalman.tracks[0].len);
148  return muE;
149  }
150  );
151 
152  const Var kRun = SIMPLEVAR(hdr.run);
153  const Var kSAHadE(
154  [kVisHadE, kRun](const caf::SRProxy* sr)
155  {
156  float hadE = SAHadEFD(kVisHadE(sr), kRun(sr));
157  return hadE;
158  }
159  );
160 
161  const Var kSACCE = kSAMuE + kSAHadE;
162 
163  const int kNumPlots = 31;
164 
165  Plot plots[kNumPlots] = {
166  {"SA neutrino energy (GeV)", "SACCE", Binning::Simple(50, 0 ,5), kSACCE},
167  {"Cosmic rejection BDT output", "bdt", Binning::Simple(60, 0.5, 0.8), kNumuContPID},
168  {"CVNm", "cvnm", Binning::Simple(50,0.5, 1), kCVNm},
169  {"Track start X position (m)", "trkStartX", kXYBins, kTrkStartX},
170  {"Track start Y position (m)", "trkStartY", kXYBins, kTrkStartY},
171  {"Track start Z position (m)", "trkStartZ", kZBins, kTrkStartZ},
172  {"cos#theta_{NuMI}", "cosNumi", Binning::Simple(40, 0.3, 1), kCosNumi},
173  {"Muon ID", "remid", Binning::Simple(41,0.50-1./160,1.+1./160), kRemID},
174  {"Reconstructed muon energy (GeV)","muonE", Binning::Simple(50,0,5), kMuE},
175  {"Hadronic energy fraction", "hadEfrac",Binning::Simple(50,0,1), kHadEFrac},
176  {"Calorimetric energy (GeV))", "calE", kEnergyBinning, SIMPLEVAR(slc.calE)},
177  {"Hadronic energy (GeV))", "hadE", kEnergyBinning, kHadE},
178  {"Calorimetric hadronic energy (GeV))", "hadcalE", Binning::Simple(30,0,3), kNumuHadCalE},
179  // {"Maximum Y (m)", "maxy", kXYBins, kSlcMaxY},
180  // {"Maximum X (m)", "maxx", kXYBins, kSlcMaxX},
181  // {"Maximum Z (m)", "maxz", kZBins, kSlcMaxZ},
182  {"Number of tracks in slice", "nkal", Binning::Simple(10, 0.5, 10.5), kNKalman},
183  {"Total number of hits", "nhit", Binning::Simple(40, 0, 400), kNHit},
184  {"Cosmic rejection BDT output", "bdt", Binning::Simple(60, 0.5, 0.8), kNumuContPID},
185  {"Track end X position (m)", "trkEndX", kXYBins, kTrkEndX},
186  {"Track end Y position (m)", "trkEndY", kXYBins, kTrkEndY},
187  {"Track end Z position (m)", "trkEndZ", kZBins, kTrkEndZ},
188  {"Number of planes to front", "nPlanesFront", Binning::Simple(30,0,900), kNplanesToFront},
189  {"Number of planes to back", "nPlanesBack", Binning::Simple(30,0,900), kNplanesToBack},
190  {"Number of cells to edge", "nCellsFromEdge", Binning::Simple(36,0,180), kNcellsFromEdge},
191  {"Basic track forward distance to edge (m)", "cosFwdDist", kZBins, kCosFwdDist},
192  {"Basic track backward distance to edge (m)", "cosBackDist", kZBins, kCosBackDist},
193  {"Kalman track forward distance to edge (m)", "kalFwdDist", kZBins, kKalmanFwdDist},
194  {"Kalman track backward distance to edge (m)", "kalBackDist", kZBins, kKalmanBackDist},
195  {"Number of hits in main track", "trkNhits", Binning::Simple(40,0,400), kTrkNhits},
196  {"Number of hits not in main track", "hadNhit", Binning::Simple(30, 0, 150), kHadNHit},
197  {"Length of main track (m)", "trkLength", Binning::Simple(50,0,25), kTrkLength},
198  {"Scattering log-likelihood", "scatLL", Binning::Simple(40,-0.2,0.3), kReMIdScatLLH},
199  {"dE/dx log-likelihood", "dedxLL", Binning::Simple(40,-0.1,0.5), kReMIdDEDxLLH}
200  };
201 
202  PredictionNoExtrap* predictionVec[kNumPlots][4];
203 
204  std::string quantpath = "/pnfs/nova/persistent/analysis/numu/Ana2017/all_FD_histo_for_quantile_cuts.root";
205  TFile* quantfile = TFile::Open(pnfs2xrootd(quantpath).c_str());
206  TH2* quanthist = (TH2*)quantfile->Get("dir_FDSpec2D/FDSpec2D");
207  std::vector<Cut> QuantCuts = QuantileCutsFromTH2(quanthist, kNumuCCOptimisedAxis, kHadEFracAxis, 4);
208 
209  for(int i = 0; i < kNumPlots; i++)
210  {
211  if(JustOne >= 0 && i != JustOne) continue;
212 
213  HistAxis kAxis("Reconstructed neutrino energy (GeV)", Binning::Simple(50,0,5), kCCE, plots[i].label, plots[i].bins, plots[i].var);
214  for(int j = 0; j < 4; j++){
215  predictionVec[i][j] = new PredictionNoExtrap(loaders, kAxis, kNumuCutFD2017 && QuantCuts[j], kNoShift, kPPFXFluxCVWgt * kXSecCVWgt2017);
216  } // loop over quantiles
217  }
218 
219  ////////////////////
220  // GOOOOOOOO!!!!!
221  loaders.Go();
222  ////////////////////
223 
224  for(int i = 0; i < kNumPlots; i++)
225  {
226 
227  if(JustOne >= 0 && i != JustOne) continue;
228 
229  for(int j = 0; j < 4; j++)
230  {
231  TGraph* gr = MakeGraph("/pnfs/nova/persistent/users/bzamoran/CheckEvents/fd_data.root", kNumuCutFD2017 && QuantCuts[j] && kInBeamSpill, kCCE, plots[i].var, &kStandardSpillCuts);
232 
233  TGraph* gr2 = MakeGraph("/pnfs/nova/persistent/users/bzamoran/CheckEvents/fd_data.root", kNumuCutFD2017 && QuantCuts[j] && kInTimingSideband, kCCE, plots[i].var, &kStandardSpillCuts);
234 
235  TH2* h2 = predictionVec[i][j]->Predict(&calc).ToTH2(kAna2017POT);
236  double max1 = h2->GetBinContent(h2->GetMaximumBin());
237  h2->Scale(1./max1);
238  h2->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
239  h2->GetYaxis()->SetTitle( (plots[i].label).c_str() );
240  h2->GetYaxis()->CenterTitle();
241  h2->GetXaxis()->CenterTitle();
242 
243  TCanvas* c1 = new TCanvas("c1","c1");
244 
245  h2->Draw("COLZ");
246 
247  gr2->SetMarkerStyle(34);
248  gr2->SetMarkerColor(kRed);
249  gr2->Draw("P");
250 
251  TGraph* grIn = (TGraph*) gr->Clone();
252  gr->SetMarkerStyle(24);
253  gr->SetMarkerColor(kBlack);
254  grIn->SetMarkerStyle(20);
255  grIn->SetMarkerColor(kYellow);
256  gr->Draw("P");
257  grIn->Draw("P");
258 
260  if(i == 0) DrawDipBoundariesY();
261 
262  std::string OutFileName = "plots/" + plots[i].fname + (std::string)Form("_quant%d_", j+1 );
263 
264  c1->SaveAs( (OutFileName + ".root").c_str() );
265  c1->SaveAs( (OutFileName + ".pdf").c_str() );
266  }
267  }
268 
269 }
const Var kHadE
Definition: NumuVars.h:23
const Var kHadNHit([](const caf::SRProxy *sr){unsigned int nought=0;if(sr->trk.kalman.ntracks< 1) return nought;return sr->slc.nhit-sr->trk.kalman.tracks[0].nhit;})
Definition: NumuVars.h:61
const Var kNKalman
Definition: NumuVars.cxx:540
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
Definition: Spectrum.cxx:165
Far Detector at Ash River.
Definition: SREnums.h:11
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
const Var kNcellsFromEdge
Definition: NumuVars.cxx:560
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
const Var kReMIdScatLLH
Definition: NumuVars.cxx:555
const Var kVisHadE([](const caf::SRProxy *sr){int ibesttrk=kBestTrack(sr);if(sr->trk.kalman.ntracks< 1) return-1000.f;if(ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return-1000.f;float ExtraHadE=sr->trk.kalman.tracks[ibesttrk].overlapE;float CalhadE=sr->slc.calE-sr->trk.kalman.tracks[ibesttrk].calE;float vishadE=CalhadE+ExtraHadE;return vishadE;})
Definition: NumuCCIncVars.h:57
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Cut kNumuCutFD2017
Definition: NumuCuts2017.h:39
caf::Proxy< size_t > ntracks
Definition: SRProxy.h:1733
const Var kNumuContPID
Definition: NumuVars.cxx:553
caf::Proxy< caf::SRNumuEnergy > numu
Definition: SRProxy.h:213
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2085
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
const Var kTrkStartY([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].start.Y()/100;})
Definition: NumuVars.h:52
caf::Proxy< float > hadtrkE
Definition: SRProxy.h:170
void FD_plots(int JustOne=-1)
Definition: FD_plots.C:95
const Var kCosBackDist([](const caf::SRProxy *sr){return sr->sel.contain.cosbakdist/100.;})
Definition: NumuVars.h:71
const double kAna2017POT
Definition: Exposures.h:174
const Var kCosFwdDist([](const caf::SRProxy *sr){return sr->sel.contain.cosfwddist/100.;})
Definition: NumuVars.h:72
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: UtilsExt.cxx:237
caf::Proxy< caf::SREnergyBranch > energy
Definition: SRProxy.h:2095
osc::OscCalcDumb calc
const Var kTrkLength([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-5.f;return sr->trk.kalman.tracks[0].len/100;})
Definition: NumuVars.h:65
#define M_PI
Definition: SbMath.h:34
TGraph * gr1
Definition: compare.C:42
const Var kNumuHadCalE
Definition: NumuVars.cxx:538
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
const char * label
TGraph * gr2
Definition: compare.C:43
const Cut kInTimingSideband([](const caf::SRProxy *sr){if(sr->spill.run > util::kLastBadTimingRun) return(kInTimingSideband_before(sr)|| kInTimingSideband_after(sr));else return(kInTimingSideband_before(sr)|| kInTimingSideband_afterA(sr)|| kInTimingSideband_afterB(sr));}, [](const caf::SRSpillProxy *spill){if(spill->run > util::kLastBadTimingRun) return(kInTimingSideband_before.Livetime(spill)+ kInTimingSideband_after.Livetime(spill));else return(kInTimingSideband_before.Livetime(spill)+ kInTimingSideband_afterA.Livetime(spill)+ kInTimingSideband_afterB.Livetime(spill));}, [](const caf::SRSpillProxy *spill){return 0;})
Definition: TimingCuts.h:12
void SetSpillCut(const SpillCut &cut)
Definition: Loaders.cxx:121
#define SIMPLEVAR(CAFNAME)
For Vars where literally all you need is a single CAF variable.
Definition: Var.h:88
====================================================================== ///
Definition: CutFlow_Data.C:28
const Var kKalmanFwdDist([](const caf::SRProxy *sr){return sr->sel.contain.kalfwddist/100.;})
Definition: NumuVars.h:74
const Var kRemID
PID
Definition: Vars.cxx:81
caf::Proxy< caf::SRTrackBranch > trk
Definition: SRProxy.h:2104
Spectrum Predict(osc::IOscCalc *calc) const override
const Var kHadEFrac
Definition: NumuVars.h:24
const int kNumPlots
Definition: GetSpectra.h:1
const Var kNHit
Definition: Vars.cxx:71
const std::vector< Plot > plots
const Var kTrkStartZ([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].start.Z()/100;})
Definition: NumuVars.h:53
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
const Var kCCE
Definition: NumuVars.h:21
const Binning kEnergyBinning
const Binning kXYBins
Definition: VarsAndCuts.h:95
caf::StandardRecord * sr
const double j
Definition: BetheBloch.cxx:29
caf::Proxy< float > hadcalE
Definition: SRProxy.h:168
void SetTh23(const T &th23) override
const Var kKalmanBackDist([](const caf::SRProxy *sr){return sr->sel.contain.kalbakdist/100.;})
Definition: NumuVars.h:75
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 Var kTrkStartX([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].start.X()/100;})
Definition: NumuVars.h:51
void DrawDipBoundaries()
Definition: FD_plots.C:59
const SystShifts kNoShift
Definition: SystShifts.cxx:21
const Var kTrkEndZ([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].stop.Z()/100;})
Definition: NumuVars.h:57
TH1F * h2
Definition: plot.C:45
void SetDmsq32(const T &dmsq32) override
float SAMuE(double trkLen)
Definition: NumuEFxs.cxx:647
const Binning bins
Definition: NumuCC_CPiBin.h:8
std::string fname
Definition: CutFlow_Data.C:30
caf::Proxy< caf::SRKalman > kalman
Definition: SRProxy.h:1752
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void DrawDipBoundariesY()
Definition: FD_plots.C:77
float SAHadEFD(double hadVisE, int run)
Definition: NumuEFxs.cxx:718
const Var kTrkEndY([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].stop.Y()/100;})
Definition: NumuVars.h:56
const Var kTrkEndX([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].stop.X()/100;})
Definition: NumuVars.h:55
const Var kNplanesToBack
Definition: NumuVars.cxx:559
void SetdCP(const T &dCP) override
std::vector< Loaders * > loaders
Definition: syst_header.h:386
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
c1
Definition: demo5.py:24
const Var kReMIdDEDxLLH
Definition: NumuVars.cxx:556
const Var kNplanesToFront
Definition: NumuVars.cxx:558
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
caf::Proxy< std::vector< caf::SRKalmanTrack > > tracks
Definition: SRProxy.h:1735
Prediction that just uses FD MC, with no extrapolation.
const Var kRun
Definition: Vars.cxx:20
const Var kXSecCVWgt2017
Definition: XsecTunes.h:37
Float_t e
Definition: plot.C:35
const Var kMuE
Definition: NumuVars.h:22
const Binning kZBins
Definition: VarsAndCuts.h:96
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
const Var kCVNm
PID
Definition: Vars.cxx:39
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:114
TGraph * MakeGraph(const std::string &wildcard, const Cut &cut, const Var &VarX, const Var &VarY, const SpillCut *spillCut)
Get a graph with two variables when a histogram is not the best choice.
Definition: GraphDrawer.cxx:13
const Var kCosNumi([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks > 0 &&sr->trk.kalman.idxremid!=999){if(sr->hdr.det==1){return sr->trk.kalman.tracks[0].dir.Dot(beamDirND);}if(sr->hdr.det==2){return sr->trk.kalman.tracks[0].dir.Dot(beamDirFD);}}return-5.f;})
Definition: NumuVars.h:43
const Var kTrkNhits([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return 65535;return int(sr->trk.kalman.tracks[0].nhit);})
Definition: NumuVars.h:59
T asin(T number)
Definition: d0nt_math.hpp:60