bdtcvn.C
Go to the documentation of this file.
2 #include "CAFAna/Core/Spectrum.h"
5 
6 #include "CAFAna/Cuts/Cuts.h"
10 
11 #include "CAFAna/Vars/Vars.h"
12 #include "CAFAna/Vars/HistAxes.h"
15 
17 
18 #include "CAFAna/Vars/FitVars.h"
19 
21 
22 #include "CAFAna/Analysis/Calcs.h"
23 #include "OscLib/IOscCalc.h"
24 
25 #include "TH2.h"
26 #include "TH1.h"
27 #include "TCanvas.h"
28 #include "TLegend.h"
29 #include "TLatex.h"
30 #include "TColor.h"
31 #include "TStyle.h"
32 #include "TLine.h"
33 #include "TROOT.h"
34 
35 using namespace ana;
36 
37 // Put a "NOvA Preliminary" tag in the corner
39 {
40  TLatex* prelim = new TLatex(.9, .95, "NOvA Preliminary");
41  prelim->SetTextColor(kBlue);
42  prelim->SetNDC();
43  prelim->SetTextSize(2/30.);
44  prelim->SetTextAlign(32);
45  prelim->Draw();
46 }
47 
48 
49 // Custom function to translate 1D axis into 2D
50 TH2* TH1ToTH2(TH1* h)
51 {
52  TH2* ret = new TH2F(UniqueName().c_str(),";CVN;BDT",25,0.75,1,45,0.26,0.71);
53 
54  for (int i = 1; i <= h->GetNbinsX(); i++){
55  int xbin = (i-1)/45 + 1; // CVN bin
56  int ybin = (i-1)%45 + 1; // BDT bin
57  ret->SetBinContent(xbin,ybin,h->GetBinContent(i));
58  }
59  return ret;
60 }
61 
62 // Makes a plot of signal purity in bdt-cvn space for peripheral pre-selected
63 // events
65 {
66  double pot = isRHC ? kAna2018RHCPOT : kAna2018FHCPOT;
67  double lt = isRHC ? kAna2018RHCLivetime : kAna2018FHCLivetime;
68  std::string appen = isRHC ? "RHC" : "FHC";
69  std::string fname = "Ana2018_cvnbdt.root";
70 
71  // Set up an osc calculator, now the 2017 best fit
73  kFitDeltaInPiUnits.SetValue(calc, 0.805);
74  kFitSinSqTheta23. SetValue(calc, 0.594);
75  kFitSinSq2Theta13. SetValue(calc, 0.082);
76  kFitDmSq32. SetValue(calc, 2.49e-3);
77 
78  // Pull our peripheral prediciton in bdt-cvn
80  LoadFromFile<PredictionNoExtrap>(fname,"pred_"+appen).release();
81  TH1 *htot1D = pred->Predict(calc).ToTH1(kAna2017POT);
82  TH2 *htot = TH1ToTH2(htot1D); // translate 1D axis to a 2D grid
83  TH1 *hsig1D = pred->PredictComponent(calc,Flavors::kNuMuToNuE,
85  if (isRHC)
86  hsig1D = pred->PredictComponent(calc,Flavors::kNuMuToNuE,
88  TH2 *hsig = TH1ToTH2(hsig1D); // Hist used in plot for overlaying boxes
89  TH2 *hfrac= TH1ToTH2(hsig1D); // Hist used in plot for purity color
90 
91  // Also need to pull the cosmic spectrum to work out signal purity
92  Spectrum *cos = LoadFromFile<Spectrum>(fname,"cos_"+appen).release();
93  TH2 *hcos = TH1ToTH2(cos->ToTH1(lt,kLivetime));
94 
95  htot->Add(hcos);
96  htot->RebinY(3);
97  hfrac->RebinY(3);
98  hsig->RebinY(3);
99  hfrac->Divide(htot); // Purity with all beam + cosmics in denominator
100  // We have the two histograms used in our plot, now make the canvas
101 
102  hfrac->GetYaxis()->SetTitle("Cosmic Rejection BDT");
103  hfrac->GetYaxis()->CenterTitle();
104  hfrac->GetXaxis()->SetTitle("CVN");
105  hfrac->GetXaxis()->CenterTitle();
106  hfrac->GetXaxis()->SetRangeUser(0.90,1.0);
107 
108  // Set up a mono-chromatic (white-purple) color for our event purity
109  Double_t Red[3] = {1.0,0.4};
110  Double_t Green[3] = {1.0,0.4};
111  Double_t Blue[3] = {1.0,1.0};
112  Double_t Length[3] = {0.0,1.0};
113  Int_t colors[255];
114  Int_t FI = TColor::CreateGradientColorTable(2,Length,Red,Green,Blue,255);
115  for (int i = 0; i < 255; i++) colors[i] = FI+i;
116  gStyle->SetPalette(255, colors);
117 
118 
119  TCanvas *can = new TCanvas();
120 
121  hsig->SetLineColor(kBlack);
122 
123  // Adjust to make boxes the size you want
124  double boxScale = isRHC ? 10 : 1;
125  hsig->Scale(boxScale);
126  hfrac->GetZaxis()->SetRangeUser(0,1);
127 
128  hfrac->Draw("colz");
129  hsig->Draw("same,box");
130 
131  // And now draw the cut regions
132  if (isRHC){
133  TLine *l1 = new TLine(0.98,0.71,0.98,0.57);
134  TLine *l2 = new TLine(0.98,0.57,1.00,0.57);
135  l1->SetLineColor(kRed);
136  l2->SetLineColor(kRed);
137  l1->SetLineWidth(5);
138  l2->SetLineWidth(5);
139  l1->Draw();
140  l2->Draw();
141  }
142  else{
143  TLine *l1 = new TLine(0.96,0.71,0.96,0.53);
144  TLine *l2 = new TLine(0.96,0.53,0.99,0.53);
145  TLine *l3 = new TLine(0.99,0.53,0.99,0.26);
146  l1->SetLineColor(kRed);
147  l2->SetLineColor(kRed);
148  l3->SetLineColor(kRed);
149  l1->SetLineWidth(5);
150  l2->SetLineWidth(5);
151  l3->SetLineWidth(5);
152  l1->Draw();
153  l2->Draw();
154  l3->Draw();
155  }
156 
157  // Comment off if including in a paper
159 
160  can->SaveAs(("Fig_PeripheralSel_"+appen+".png").c_str());
161  can->SaveAs(("Fig_PeripheralSel_"+appen+".pdf").c_str());
162 }
163 
165 {
166 
167  // Basic MC weights, along with our appropriate CVN / BDT axes and vars
169  const HistAxis kBDTAxis("CosRej BDT",Binning::Simple(45,0.26,0.71),
171  const HistAxis kCVNAxis("CVN",Binning::Simple(50,0.75,1),kCVNSSe);
172  Var kCVNBDT = Var2D(kCVNe,Binning::Simple(25,0.75,1),
173  kCosPIDContain,Binning::Simple(45,0.26,0.71));
174 
175  // Set up cosmic histogram from timing sideband
176  // Must set up both FHC and RHC
177  SpectrumLoader lCosFHC("prod4_sumdecaf_fd_numi_fhc_full_goodruns_nue2018");
179  Spectrum cos_FHC(lCosFHC, kCVNAxis, kBDTAxis, // Must apply s.b. timing cut
181  SpectrumLoader lCosRHC("prod4_sumdecaf_fd_numi_rhc_full_goodruns_nue2018");
183  Spectrum cos_RHC(lCosRHC, kCVNAxis, kBDTAxis,
185 
186  // Set up MC prediction from concats
187  SpectrumLoader lSwapRHC("prod_sumdecaf_R17-11-14-prod4reco.e_fd_genie_fluxswap_rhc_nova_v08_full_v1_nue2018");
188  SpectrumLoader lNonSRHC("prod_sumdecaf_R17-11-14-prod4reco.e_fd_genie_nonswap_rhc_nova_v08_full_v1_nue2018");
189  SpectrumLoader lTauSRHC("prod_sumdecaf_R17-11-14-prod4reco.e_fd_genie_tau_rhc_nova_v08_full_v1_nue2018");
190 
194 
195  PredictionNoExtrap pred_RHC(lNonSRHC,lSwapRHC,lTauSRHC,"CVN/BDT",
196  Binning::Simple(1125,0,1125),
197  kCVNBDT,kNue2018PeripheralPresel,
198  kNoShift,wei);
199 
200  SpectrumLoader lSwapFHC("prod_sumdecaf_R17-11-14-prod4reco.d_fd_genie_fluxswap_fhc_nova_v08_full_v1_nue2018");
201  SpectrumLoader lNonSFHC("prod_sumdecaf_R17-11-14-prod4reco.d_fd_genie_nonswap_fhc_nova_v08_full_v1_nue2018");
202  SpectrumLoader lTauSFHC("prod_sumdecaf_R17-11-14-prod4reco.d_fd_genie_tau_fhc_nova_v08_full_v1_nue2018");
203 
207 
208  PredictionNoExtrap pred_FHC(lNonSFHC,lSwapFHC,lTauSFHC,"CVN/BDT",
209  Binning::Simple(1125,0,1125),
210  kCVNBDT,kNue2018PeripheralPresel,
211  kNoShift,wei);
212 
213  // Press go
214  lCosRHC.Go();
215  lCosFHC.Go();
216 
217  lSwapRHC.Go();
218  lNonSRHC.Go();
219  lTauSRHC.Go();
220  lSwapFHC.Go();
221  lNonSFHC.Go();
222  lTauSFHC.Go();
223 
224  TFile *out = new TFile("Ana2018_cvnbdt.root","recreate");
225  cos_FHC.SaveTo(out, "cos_FHC");
226  cos_RHC.SaveTo(out, "cos_RHC");
227 
228  pred_FHC.SaveTo(out, "pred_FHC");
229  pred_RHC.SaveTo(out, "pred_RHC");
230  // Done with filling the histograms
231 }
232 
233 void bdtcvn(bool fillSpectra)
234 {
235  if (fillSpectra){
236  FillSpectra();
237  }
238  else{
239  MakePeriCutPlot(true);
240  MakePeriCutPlot(false);
241  }
242 }
243 
Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
double lt
Antineutrinos-only.
Definition: IPrediction.h:50
const Var kCVNe
PID
Definition: Vars.cxx:35
(&#39; appearance&#39;)
Definition: IPrediction.h:18
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:176
void MakePeriCutPlot(bool isRHC)
Definition: bdtcvn.C:64
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
const Cut kNue2018PeripheralPresel(kNue2018PeripheralPreselFunc)
const Var kCosPIDContain
Definition: NueVars.cxx:110
const double kAna2017POT
Definition: Exposures.h:174
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Forward to wrapped Var&#39;s SetValue()
void SetSpillCut(const SpillCut &cut)
TH1F * hsig
Definition: Xdiff_gwt.C:59
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:39
const double kAna2018RHCPOT
Definition: Exposures.h:208
void FillSpectra()
Definition: bdtcvn.C:164
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
Charged-current interactions.
Definition: IPrediction.h:39
int colors[6]
Definition: tools.h:1
Spectrum Predict(osc::IOscCalc *calc) const override
void PrintPreliminary()
Definition: bdtcvn.C:38
#define pot
virtual void Go() override
Load all the registered spectra.
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:534
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
TLatex * prelim
Definition: Xsec_final.C:133
const Var kCVNSSe([](const caf::SRProxy *sr){throw std::runtime_error("kCVNSSe is no longer available. Fix your macro so you don't use it.");return-5.;})
2018 nue PID
Definition: Vars.h:52
Neutrinos-only.
Definition: IPrediction.h:49
const SystShifts kNoShift
Definition: SystShifts.cxx:22
TH2 * TH1ToTH2(TH1 *h)
Definition: bdtcvn.C:50
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Var kXSecCVWgt2018
Definition: XsecTunes.h:48
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
T cos(T number)
Definition: d0nt_math.hpp:78
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
Prediction that just uses FD MC, with no extrapolation.
const double kAna2018FHCPOT
Definition: Exposures.h:207
Float_t e
Definition: plot.C:35
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
enum BeamMode kBlue
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
void bdtcvn(bool fillSpectra)
Definition: bdtcvn.C:233
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:28
const double kAna2018FHCLivetime
Definition: Exposures.h:213
const double kAna2018RHCLivetime
Definition: Exposures.h:214
enum BeamMode string