make_nue_pidtuning.C
Go to the documentation of this file.
1 #include "OscLib/OscCalcDumb.h"
3 
4 #include "CAFAna/Core/Spectrum.h"
6 #include "CAFAna/Systs/Systs.h"
9 
10 
11 
12 #include "CAFAna/Cuts/Cuts.h"
13 #include "CAFAna/Cuts/TimingCuts.h"
15 #include "CAFAna/Cuts/SpillCuts.h"
16 
17 
18 #include "CAFAna/Vars/Vars.h"
19 #include "CAFAna/Vars/FitVars.h"
23 #include "CAFAna/Vars/HistAxes.h"
24 
26 
27 #include "Utilities/func/MathUtil.h"
28 
29 #include "TCanvas.h"
30 #include "TH2.h"
31 #include "TH1.h"
32 #include "TFile.h"
33 #include "TLegend.h"
34 
35 #include <iostream>
36 #include <cmath>
37 
38 using namespace ana;
39 
40 
41 
42 // Quick numbers, POT, livetime, and osc parameters
43 // livetime/pot from Lyudmila's checks
44 
45 double pot = 8.0e20;
46 double lt = 199;
47 
48 // Change these here to look at cut stability if you want
49 double dCP = 1.2*TMath::Pi();
50 double dmsq32 = 0.00244;
51 double ssth23 = 0.56;
52 double ss2th13= 0.082;
53 
54 
55 
56 // Need a way to translate CVN-BDT 1D plot into a TH2
57 TH2 *TH1ToTH2(TH1 *h)
58 {
59  TH2 *ret = new TH2F(UniqueName().c_str(),"",50,0.5,1,50,0.35,0.85);
60  for (int i = 1; i <= 2500; i++){
61  int xbin = (i-1)%50 + 1;
62  int ybin = (i-1)/50 + 1;
63  ret->SetBinContent(xbin,ybin,h->GetBinContent(i));
64  }
65  return ret;
66 }
67 
69 {
70  std::string fname = "out_fdpred.root";
72  LoadFromFile<PredictionNoExtrap>(fname,"pred"+name).release();
73  Spectrum *cos = LoadFromFile<Spectrum>(fname,"cos"+name).release();
74 
75  // Setup oscillations
77 
78  kFitDeltaInPiUnits.SetValue(calc,dCP/TMath::Pi());
82 
83  // Make signal, background spectra
86  Spectrum bb = pred->Predict(calc);
87  bb -= sig;
88 
89  TH2 *hbb = TH1ToTH2(bb .ToTH1(pot));
90  TH2 *hsig = TH1ToTH2(sig.ToTH1(pot));
91  TH2 *hcos = TH1ToTH2(cos->ToTH1(lt,kLivetime));
92 
93  TH2 *htot = TH1ToTH2(cos->ToTH1(lt,kLivetime));
94  htot->Add(hbb);
95  htot->Add(hsig);
96 
97  std::vector<double> sigs(50,0);
98  std::vector<double> tots(50,0);
99 
100  // scan over CVN bins
101  for (int i = 1; i <= 50; i++){
102  double fom = 0;
103  double bestY = 0;
104  // Calculate a best BDT cut
105  for (int j = 1; j <= 50; j++){
106  double sig = hsig->Integral(i,i,j,50);
107  double tot = htot->Integral(i,i,j,50);
108  if (tot == 0) continue;
109  double cfom = sig/sqrt(tot);
110  if (cfom > fom){
111  fom = cfom;
112  bestY = 0.35 + 0.01*(j-1);
113  }
114  }
115  if (i >=46)
116  std::cout << 0.49 + 0.01*i << " " << bestY << " " << fom << std::endl;
117  //foms.push_back(fom);
118  for (int j = 0; j < i; j++){
119  // Accumalate sig / tot events above this CVN cut, and BDT cuts
120  sigs[j] += hsig->Integral(i,i,(bestY-0.35)/0.01+1,50);
121  tots[j] += htot->Integral(i,i,(bestY-0.35)/0.01+1,50);
122  }
123  }
124 
125  TH1 *h = new TH1F(UniqueName().c_str(),"",50,0.5,1); // FOM vs CVN plot
126  for (int i = 50; i >= 1; i--){
127  h->SetBinContent(i,sigs[i-1]/sqrt(tots[i-1]));
128  }
129 
130  h->GetXaxis()->SetTitle("CVN-SS");
131  h->GetYaxis()->SetTitle("Total FOM");
132  h->GetXaxis()->CenterTitle();
133  h->GetYaxis()->CenterTitle();
134 
135  TCanvas *can = new TCanvas();
136  h->Draw();
137 
138 }
139 
141 {
142  std::string fname = "out_fdpred.root";
144  LoadFromFile<PredictionNoExtrap>(fname,"pred"+name).release();
145  Spectrum *cos = LoadFromFile<Spectrum>(fname,"cos"+name).release();
146 
148 
149  kFitDeltaInPiUnits.SetValue(calc,dCP/TMath::Pi());
150  kFitDmSq32.SetValue(calc,dmsq32);
153 
154 
157  Spectrum bb = pred->Predict(calc);
158  bb -= sig;
159 
160  TH1 *hbb = bb.ToTH1(pot);
161  TH1 *hcos = cos->ToTH1(lt,kLivetime);
162  TH1 *hsig = sig.ToTH1(pot);
163 
164  TH1 *htot = cos->ToTH1(lt,kLivetime);
165  htot->Add(hbb);
166  htot->Add(hsig);
167 
168  double fom = -1;
169  // Bin boundaries for CVN
170  double best1 = -1;
171  double best2 = -1;
172  double best3 = -1;
173 
174  std::cout << "Finding optimal FOM for 1 CVN bin analysis" << std::endl;
175  for (int i = 1; i <= 50; i++){
176  double sig1 = hsig->Integral(i,50);
177  double tot1 = htot->Integral(i,50);
178  double fom1 = sig1/sqrt(tot1);
179  //double fom1 = sig1/sqrt(tot1);
180  if (fom1 > fom){
181  fom = fom1;
182  best1 = 0.5 + 0.01*(i-1);
183  }
184  }
185  std::cout << fom << ": " << best1 << "-" << "1" << std::endl;
186 
187 
188  std::cout << "Finding optimal FOM for 2 CVN bin analysis" << std::endl;
189  fom = -1; // reset fom
190  for (int i = 1; i <= 49; i++){ // Can't go all the way to 1
191  for (int j = i+1; j <= 50; j++){ // Must start above 1 bin max
192  double sig1 = hsig->Integral(i,j-1);
193  double tot1 = htot->Integral(i,j-1);
194  double fom1 = sig1/sqrt(tot1);
195  double sig2 = hsig->Integral(j,50);
196  double tot2 = htot->Integral(j,50);
197  double fom2 = sig2/sqrt(tot2);
198  if (sqrt(fom1*fom1+fom2*fom2) > fom){
199  fom = sqrt(fom1*fom1+fom2*fom2);
200  best1 = 0.5 + 0.01*(i-1);
201  best2 = 0.5 + 0.01*(j-1);
202  }
203  }
204  }
205  std::cout << fom << ": " << best1 << "-" << best2 << "-1" << std::endl;
206 
207 
208  std::cout << "Finding optimal FOM for 2 CVN bin analysis" << std::endl;
209  fom = -1; // Reset fom
210  for (int i = 1; i <= 48; i++){
211  for (int j = i+1; j <= 49; j++){
212  for (int k = j+1; k <= 50; k++){
213  double sig1 = hsig->Integral(i,j-1);
214  double tot1 = htot->Integral(i,j-1);
215  double fom1 = sig1/sqrt(tot1);
216 
217  double sig2 = hsig->Integral(j,k-1);
218  double tot2 = htot->Integral(j,k-1);
219  double fom2 = sig2/sqrt(tot2);
220 
221  double sig3 = hsig->Integral(k,50);
222  double tot3 = htot->Integral(k,50);
223  double fom3 = sig2/sqrt(tot2);
224 
225  double tfom = sqrt(fom1*fom1+fom2*fom2+fom3*fom3);
226  if (tfom>fom){
227  fom = tfom;
228  best1 = 0.5 + 0.01*(i-1);
229  best2 = 0.5 + 0.01*(j-1);
230  best3 = 0.5 + 0.01*(k-1);
231  }
232  }
233  }
234  }
235  std::cout << fom << ": " << best1 << "-" << best2 << "-" << best3 << "-1" << std::endl;
236 
237 }
238 
239 
240 
241 
242 
243 
244 // TODO: move these cuts to a header, once PID cuts are finalized
245 
246 const Cut kNue2018CVNVsCosPID([](const caf::SRProxy *sr){
247  double cvn = sr->sel.cvnProd3Train.nueid;
248  double bdt = sr->sel.nuecosrej.cospidcontain;
249  if (cvn > 0.99 || (cvn>0.95 && bdt > 0.57)) return true;
250  return false;
251  });
252 
255 
256 const Var kNue2018SelectionBin([](const caf::SRProxy *sr){
257  bool isPeri = kNue2018FDPeripheral(sr);
258  double cvn = sr->sel.cvnProd3Train.nueid;
259  if (isPeri)
260  return float(3.5);
261  if (cvn > 0.75 && cvn < 0.87)
262  return float(0.5);
263  if (cvn > 0.87 && cvn < 0.95)
264  return float(1.5);
265  if (cvn > 0.95)
266  return float(2.5);
267  return float(-0.5);
268  });
269 
270 
271 const Var kNue2018AnaBin([](const caf::SRProxy *sr){
272  int selBin = kNue2018SelectionBin(sr);
273 
274  float nuE = kNueEnergy2017(sr);
275  int nuEBin = nuE/0.5;
276  assert(nuEBin <= 8 && "An event with nuE > 4.5 should never happen");
277 
278  int anaBin = 9*selBin + nuEBin;
279 
280  if (anaBin > 36)
281  anaBin = -1;
282 
283  return anaBin;
284  });
285 
286 
287 const Var kCVNSSe([](const caf::SRProxy* sr){
288  return sr->sel.cvnProd3Train.nueid;
289  });
290 
291 const Cut kNue2018FDCore = kCVNSSe > 0.75 && kNue2017CorePresel;
292 const Cut kNue2018FD = kNue2018FDCore || kNue2018FDPeripheral;
293 
294 
295 
296 // Macro for optimizing PID cuts / binning in the FD. Just works on RHC now
297 // TODO: incorperate FHC. Update presel cuts if we decide to change them
298 void make_nue_pidtuning(bool makePlots=false)
299 {
300 
301  if (!makePlots){
302 
303  SpectrumLoader lSwap("dpershey_prod4-cosrej_rhc_mc_fluxswap");
304  SpectrumLoader lNonS("dpershey_prod4-cosrej_rhc_mc_nonswap");
305  SpectrumLoader lTauS("dpershey_prod4-cosrej_rhc_mc_tauswap");
306 
307  SpectrumLoader lCosmic("dpershey_prod4-cosrej_cosmicdata_rhctune");
308 
309  lSwap.SetSpillCut(kStandardSpillCuts);
310  lNonS.SetSpillCut(kStandardSpillCuts);
311  lTauS.SetSpillCut(kStandardSpillCuts);
312  lCosmic.SetSpillCut(kStandardSpillCuts);
313 
314  Cut selCore = kNue2017CorePresel;
315  Cut selPeri = kNue2017PeripheralPresel;
316  Var wei = kXSecCVWgt2017*kPPFXFluxCVWgt;
317 
318 
319 
320  // Peripheral cuts are 2D PID cuts, set up 2D vars now
321  Var kBDTLight = SIMPLEVAR(sel.nuecosrej.cospidlight);
322  Var kBDTCont = SIMPLEVAR(sel.nuecosrej.cospidcontain);
323 
324  Binning bdtbins = Binning::Simple(50,0.35,0.85);
325  Binning cvnbins = Binning::Simple(50,0.5,1.);
326 
327  Var kLightCVN = Var2D(kBDTLight,bdtbins,
328  kCVNSSe,cvnbins);
329  Var kContCVN = Var2D(kBDTCont,bdtbins,
330  kCVNSSe,cvnbins);
331 
332 
333  // Prediction as a function of CVN, for main PID tuning
334  // Be sure you're using the right core preselection, it may change
335  PredictionNoExtrap predCVN(lNonS,lSwap,lTauS,"",
336  Binning::Simple(50,0.5,1),
337  kCVNSSe, selCore, kNoShift, wei);
338  Spectrum cosCVN(lCosmic,HistAxis("",Binning::Simple(50,0.5,1),kCVNSSe),
339  kInCosmicTimingWindow&&selCore,kNoShift,wei);
340 
341  // 2D predictions for determining peripheral PID cuts
342  PredictionNoExtrap predLightPeri(lNonS,lSwap,lTauS,
343  "",Binning::Simple(2500,0,2500),kLightCVN,
344  selPeri,kNoShift,wei); // BDTLight
345  Spectrum cosLightPeri(lCosmic,
346  HistAxis("",Binning::Simple(2500,0,2500),kLightCVN),
347  selPeri&&kInCosmicTimingWindow,kNoShift,wei);
348  PredictionNoExtrap predContPeri (lNonS,lSwap,lTauS,
349  "",Binning::Simple(2500,0,2500),kContCVN,
350  selPeri,kNoShift,wei); // BDTContain
351  Spectrum cosContPeri (lCosmic,
352  HistAxis("",Binning::Simple(2500,0,2500),kContCVN),
353  selPeri&&kInCosmicTimingWindow,kNoShift,wei);
354 
355  // Prediction with 2017 bin boundaries, just for fun
356  PredictionNoExtrap pred2017(lNonS,lSwap,lTauS,"",
357  Binning::Simple(36,0,36),
358  kNue2018AnaBin,
359  kNue2018FD,kNoShift,wei);
360  Spectrum cos2017(lCosmic,
361  HistAxis("",Binning::Simple(36,0,36),kNue2018AnaBin),
362  kInCosmicTimingWindow&&kNue2018FD,kNoShift,wei);
363 
364 
365 
366 
367 
368 
369 
370 
371 
372  lSwap.Go();
373  lNonS.Go();
374  lTauS.Go();
375  lCosmic.Go();
376 
377  TFile* outFile = new TFile("out_fdcosrej.root","RECREATE");
378  outFile->cd();
379 
380  pred2017.SaveTo(outFile, "pred2017");
381  cos2017.SaveTo(outFile, "cos2017");
382 
383  predCVN.SaveTo(outFile, "predCVN");
384  cosCVN.SaveTo(outFile, "cosCVN");
385 
386  predLightPeri.SaveTo(outFile, "predLightPeri");
387  cosLightPeri .SaveTo(outFile, "cosLightPeri");
388 
389  predContPeri.SaveTo(outFile, "predContPeri");
390  cosContPeri .SaveTo(outFile, "cosContPeri");
391 
392 
393  outFile->Close();
394  } // End if !makePlots
395  else {
396 
397  PeripheralCuts("LightPeri");
398  CVNCuts("CVN");
399 
400  }
401 }
Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
const XML_Char * name
Definition: expat.h:151
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
double lt
double ssth23
Antineutrinos-only.
Definition: IPrediction.h:50
(&#39; appearance&#39;)
Definition: IPrediction.h:18
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:148
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Definition: FitVars.cxx:177
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
T sqrt(T number)
Definition: d0nt_math.hpp:156
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Forward to wrapped Var&#39;s SetValue()
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:249
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Definition: FitVars.cxx:48
Defines an enumeration for prong classification.
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:40
const Var kNue2018AnaBin([](const caf::SRProxy *sr){int selBin=kNue2018SelectionBin(sr);float nuE=kNueEnergy2018(sr);int nuEBin=nuE/0.5;assert(nuEBin<=8 &&"An event with nuE > 4.5 should never happen");int anaBin=9 *selBin+nuEBin;return anaBin;})
Use this Analysis Binning for Ana2018, official Binning.
Definition: NueCuts2018.h:171
caf::Proxy< caf::SRNueCosRej > nuecosrej
Definition: SRProxy.h:1265
TH2 * TH1ToTH2(TH1 *h)
Charged-current interactions.
Definition: IPrediction.h:39
const Var kNueEnergy2017([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;return NueRecoE_2017FDFit(kCVNemE(sr), kCVNhadE(sr));})
Definition: NueEnergy2017.h:11
Spectrum Predict(osc::IOscCalc *calc) const override
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
void CVNCuts(std::string name)
double dCP
double pot
caf::StandardRecord * sr
const double j
Definition: BetheBloch.cxx:29
const Var kNue2018SelectionBin([](const caf::SRProxy *sr){bool isPeri=kNue2018FDPeripheral(sr);if(isPeri) return float(2.5);std::cout<< "ERROR::kNue2018SelectionBin. Looking for cvnProd3Train. Branch no longer exists."<< std::endl;abort();return float(-5.0);})
Definition: NueCuts2018.h:168
const Cut kNue2017PeripheralPresel
Definition: NueCuts2017.h:110
const Cut kNue2018FDPeripheral(kNue2018FDPeripheralFunc)
OStream cout
Definition: OStream.cxx:6
const Cut kNue2018CVNVsCosPID(kNue2018CVNVsCosPIDFunc)
double dmsq32
double ss2th13
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
T cos(T number)
Definition: d0nt_math.hpp:78
assert(nhit_max >=nhit_nbins)
caf::Proxy< float > cospidcontain
Definition: SRProxy.h:1035
def CVN(num_classes)
Definition: models.py:37
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
Prediction that just uses FD MC, with no extrapolation.
caf::Proxy< caf::SRIDBranch > sel
Definition: SRProxy.h:2141
Template for Cut and SpillCut.
Definition: Cut.h:15
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
void PeripheralCuts(std::string name)
enum BeamMode string