Toy_analyses.C
Go to the documentation of this file.
6 #include "CAFAna/Cuts/Cuts.h"
11 #include "CAFAna/Core/EventList.h"
12 #include "CAFAna/Core/Loaders.h"
16 #include "CAFAna/Core/Spectrum.h"
32 #include "CAFAna/Vars/FitVars.h"
34 #include "CAFAna/Vars/Vars.h"
37 #include "OscLib/OscCalcPMNSOpt.h"
38 
39 using namespace ana;
40 
41 #include "TCanvas.h"
42 #include "TFile.h"
43 #include "TGraph.h"
44 #include "TGraphAsymmErrors.h"
45 #include "TH1.h"
46 #include "TH2.h"
47 #include "TF1.h"
48 #include "TLatex.h"
49 #include "TLegend.h"
50 #include "TStyle.h"
51 #include "TRandom3.h"
52 #include "TTree.h"
53 
54 TRandom3 r(0); // make this global so I don't have to pass it across functions
55 
56 struct FitResult
57 {
58 double bestX;
59 double bestY;
60 double LL;
61 double sigma;
62 };
63 
64 FitResult getBestFit(IPrediction* prediction, TH1F* h1, TH1D* hCosmic);
65 FitResult getBestFit(IPrediction* quant1, IPrediction* quant2, IPrediction* quant3, IPrediction* quant4,
66  TH1F* h1, TH1F* h2, TH1F* h3, TH1F* h4,
67  TH1D* hCosmic1, TH1D* hCosmic2, TH1D* hCosmic3, TH1D* hCosmic4);
68 
69 std::vector<std::vector<double>> QuantileVectorFromTH2(TH2* quantileHist,
70  const HistAxis& independentAxis,
71  const HistAxis& quantileAxis,
72  const unsigned int& numQuantiles);
73 
74 unsigned int quantile(std::vector<std::vector<double>> v, float E, float hadEfrac);
75 void getPoissonRandom(TH1D* s, TH1F* h);
76 void getSArebin(TH1D* s, TH1D* h);
77 void getSArebin(TH1F* s, TH1F* h);
78 
79 void Toy_analyses(unsigned int nExperiments, unsigned int volume = 1)
80 {
81  std::string inputDir = "/pnfs/nova/persistent/users/bzamoran/ToyPvalue/";
82 
83  TFile* fInputTree = TFile::Open(pnfs2xrootd(inputDir + "inputTree.root").c_str());
84  TTree* inputTree = (TTree*) fInputTree->Get("inputTree");
85  float TrueE, CCE, CCESA, HadEfrac, PPFX, CXW;
86  int isSA, isA17;
87  inputTree->SetBranchAddress("TrueE", &TrueE);
88  inputTree->SetBranchAddress("isSA", &isSA);
89  inputTree->SetBranchAddress("isA17", &isA17);
90  inputTree->SetBranchAddress("CCE", &CCE);
91  inputTree->SetBranchAddress("CCESA", &CCESA);
92  inputTree->SetBranchAddress("HadEfrac", &HadEfrac);
93  inputTree->SetBranchAddress("PPFX", &PPFX);
94  inputTree->SetBranchAddress("CXW", &CXW);
95 
96  std::cout << "Tree" << std::endl;
97 
98  TFile* fPred = TFile::Open(pnfs2xrootd(inputDir + "Predictions/pred_3A_prod3_AllQuants.root").c_str());
99  PredictionInterp* prediction_OR = LoadFrom<PredictionInterp>(fPred, "prediction_OR" ).release();
100  PredictionInterp* prediction_SA = LoadFrom<PredictionInterp>(fPred, "prediction_SA" ).release();
101  fPred->Close();
102 
103  std::cout << "Predictions 1" << std::endl;
104 
105  TFile* fPred_Q1 = TFile::Open(pnfs2xrootd(inputDir + "Predictions/pred_3A_prod3_Quant1.root").c_str());
106  TFile* fPred_Q2 = TFile::Open(pnfs2xrootd(inputDir + "Predictions/pred_3A_prod3_Quant2.root").c_str());
107  TFile* fPred_Q3 = TFile::Open(pnfs2xrootd(inputDir + "Predictions/pred_3A_prod3_Quant3.root").c_str());
108  TFile* fPred_Q4 = TFile::Open(pnfs2xrootd(inputDir + "Predictions/pred_3A_prod3_Quant4.root").c_str());
109 
110  PredictionInterp* prediction_Q1 = LoadFrom<PredictionInterp>(fPred_Q1, "prediction").release();
111  PredictionInterp* prediction_Q2 = LoadFrom<PredictionInterp>(fPred_Q2, "prediction").release();
112  PredictionInterp* prediction_Q3 = LoadFrom<PredictionInterp>(fPred_Q3, "prediction").release();
113  PredictionInterp* prediction_Q4 = LoadFrom<PredictionInterp>(fPred_Q4, "prediction").release();
114 
115  std::cout << "Predictions quant" << std::endl;
116 
117  TFile* fCosmics = TFile::Open(pnfs2xrootd(inputDir + "cos_bkg_hists_v2.root").c_str());
118  TH1D* hSpecCosmics_Q1 = (TH1D*) fCosmics->Get("q1all");
119  TH1D* hSpecCosmics_Q2 = (TH1D*) fCosmics->Get("q2all");
120  TH1D* hSpecCosmics_Q3 = (TH1D*) fCosmics->Get("q3all");
121  TH1D* hSpecCosmics_Q4 = (TH1D*) fCosmics->Get("q4all");
122 
123  hSpecCosmics_Q1->Scale(kSecondAnaLivetime / kAna2017Livetime ); // scale to 6e20
124  hSpecCosmics_Q2->Scale(kSecondAnaLivetime / kAna2017Livetime );
125  hSpecCosmics_Q3->Scale(kSecondAnaLivetime / kAna2017Livetime );
126  hSpecCosmics_Q4->Scale(kSecondAnaLivetime / kAna2017Livetime );
127 
128  TF1* fSA = new TF1("fSA","TMath::Exp(1.976*(x-2.26))/(1+TMath::Exp(1.976*(x-2.26)))",0,5);
129  TF1* fA17 = new TF1("fA17","TMath::Exp(-1.907*(x-1.967))/(1+TMath::Exp(-1.907*(x-1.967)))",0,5);
130  TF1* fBoth = new TF1("fBoth","1.-fA17-fSA",0,5);
131 
132  TF1* fAdjustSA = new TF1("fAdjustSA","TMath::Min(fSA/(1.-fSA),5.0)", 0.5, 5); // don't let it blow
133  TF1* fAdjustA17 = new TF1("fAdjustA17","fA17/(1.-fSA)", 0.5, 5);
134  TF1* fAdjustBoth = new TF1("fAdjustBoth","fBoth/(1.-fSA)", 0.5, 5);
135  TF1* fAdjustSAandBoth = new TF1("fAdjustSAandBoth","TMath::Min((1.-fA17)/(1.-fSA),5.0)", 0.5, 5);
136 
137  TH1D* hSpecCosmics_SA = (TH1D*) hSpecCosmics_Q1->Clone("hSpecCosmics_SA");
138  hSpecCosmics_SA->Add(hSpecCosmics_Q2);
139  hSpecCosmics_SA->Add(hSpecCosmics_Q3);
140  hSpecCosmics_SA->Add(hSpecCosmics_Q4);
141  hSpecCosmics_SA->Multiply(fAdjustSA);
142  TH1D* hSpecCosmics_SA_rebin = new TH1D("hSpecCosmics_SA_rebin","",20,0,5);
143  getSArebin(hSpecCosmics_SA, hSpecCosmics_SA_rebin);
144 
145  TH1D* hSpecCosmics_SA_and_both = (TH1D*) hSpecCosmics_Q1->Clone("hSpecCosmics_SA_and_both");
146  hSpecCosmics_SA_and_both->Add(hSpecCosmics_Q2);
147  hSpecCosmics_SA_and_both->Add(hSpecCosmics_Q3);
148  hSpecCosmics_SA_and_both->Add(hSpecCosmics_Q4);
149  hSpecCosmics_SA_and_both->Multiply(fAdjustSAandBoth);
150  TH1D* hSpecCosmics_SA_and_both_rebin = new TH1D("hSpecCosmics_SA_and_both_rebin","",20,0,5);
151  getSArebin(hSpecCosmics_SA_and_both, hSpecCosmics_SA_and_both_rebin);
152 
153  TH1D* hSpecCosmics_Both_Q1 = (TH1D*) hSpecCosmics_Q1->Clone("hSpecCosmics_Both_Q1");
154  TH1D* hSpecCosmics_Both_Q2 = (TH1D*) hSpecCosmics_Q2->Clone("hSpecCosmics_Both_Q2");
155  TH1D* hSpecCosmics_Both_Q3 = (TH1D*) hSpecCosmics_Q3->Clone("hSpecCosmics_Both_Q3");
156  TH1D* hSpecCosmics_Both_Q4 = (TH1D*) hSpecCosmics_Q4->Clone("hSpecCosmics_Both_Q4");
157  hSpecCosmics_Both_Q1->Multiply(fAdjustBoth);
158  hSpecCosmics_Both_Q2->Multiply(fAdjustBoth);
159  hSpecCosmics_Both_Q3->Multiply(fAdjustBoth);
160  hSpecCosmics_Both_Q4->Multiply(fAdjustBoth);
161 
162  TH1D* hSpecCosmics_A17_Q1 = (TH1D*) hSpecCosmics_Q1->Clone("hSpecCosmics_A17_Q1");
163  TH1D* hSpecCosmics_A17_Q2 = (TH1D*) hSpecCosmics_Q2->Clone("hSpecCosmics_A17_Q2");
164  TH1D* hSpecCosmics_A17_Q3 = (TH1D*) hSpecCosmics_Q3->Clone("hSpecCosmics_A17_Q3");
165  TH1D* hSpecCosmics_A17_Q4 = (TH1D*) hSpecCosmics_Q4->Clone("hSpecCosmics_A17_Q4");
166  hSpecCosmics_A17_Q1->Multiply(fAdjustA17);
167  hSpecCosmics_A17_Q2->Multiply(fAdjustA17);
168  hSpecCosmics_A17_Q3->Multiply(fAdjustA17);
169  hSpecCosmics_A17_Q4->Multiply(fAdjustA17);
170 
171 
172  std::cout << "All's well with loading" << std::endl;
173 
174  // Get FD files for HadEFrac vs. neutrino energy distribution
175  std::string quantpath = "/pnfs/nova/persistent/analysis/numu/Ana2017/all_FD_histo_for_quantile_cuts.root";
176  TFile* quantfile = TFile::Open(pnfs2xrootd(quantpath).c_str());
177  TH2* quanthist = (TH2*)quantfile->Get("dir_FDSpec2D/FDSpec2D");
178  std::vector<std::vector<double>> QuantBins = QuantileVectorFromTH2(quanthist, kNumuCCOptimisedAxis, kHadEFracAxis, 4);
179 
180  // Now the fittting
181 
184  // calc.SetDmsq32(2.67e-3);
185  // calc.SetTh23(asin(sqrt(0.404))); // PRL-published Second Ana
186  calc.SetDmsq32(2.43071e-3);
187  calc.SetTh23(asin(sqrt(0.471974)));
188 
189  TFile* fOutput = new TFile(Form("Toy_Output_vol%i.root", volume),"RECREATE");
190  TTree* fTree = new TTree("T","Summary tree from p-value study");
191 
192  TH1F* hSA = new TH1F("hSA","",20,0,5);
193  TH1F* hSA_cosmic = new TH1F("hSA_cosmic","",20,0,5);
194 
195  float A17Ebins[20] = { 0, 0.75, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8,
196  1.9, 2, 2.25, 2.5, 2.75, 3, 3.5, 4, 5 };
197 
198  TH1F* hA17 = new TH1F("hA17","",19, A17Ebins);
199  TH1F* hA17_Q1 = new TH1F("hA17_Q1","",19, A17Ebins);
200  TH1F* hA17_Q2 = new TH1F("hA17_Q2","",19, A17Ebins);
201  TH1F* hA17_Q3 = new TH1F("hA17_Q3","",19, A17Ebins);
202  TH1F* hA17_Q4 = new TH1F("hA17_Q4","",19, A17Ebins);
203 
204  TH1F* hA17_cosmic = new TH1F("hA17_cosmic","",19, A17Ebins);
205  TH1F* hA17_Q1_cosmic = new TH1F("hA17_Q1_cosmic","",19, A17Ebins);
206  TH1F* hA17_Q2_cosmic = new TH1F("hA17_Q2_cosmic","",19, A17Ebins);
207  TH1F* hA17_Q3_cosmic = new TH1F("hA17_Q3_cosmic","",19, A17Ebins);
208  TH1F* hA17_Q4_cosmic = new TH1F("hA17_Q4_cosmic","",19, A17Ebins);
209 
210  TH1F* hBoth = new TH1F("hBoth","",19, A17Ebins);
211 
212  TH1F* hBoth_cosmic_Q1 = new TH1F("hBoth_cosmic_Q1","",19, A17Ebins);
213  TH1F* hBoth_cosmic_Q2 = new TH1F("hBoth_cosmic_Q2","",19, A17Ebins);
214  TH1F* hBoth_cosmic_Q3 = new TH1F("hBoth_cosmic_Q3","",19, A17Ebins);
215  TH1F* hBoth_cosmic_Q4 = new TH1F("hBoth_cosmic_Q4","",19, A17Ebins);
216  TH1F* hBoth_cosmic = new TH1F("hBoth_cosmic","",19, A17Ebins);
217  TH1F* hBoth_cosmic_rebin = new TH1F("hA17_Both_cosmic_rebin","",20,0,5);
218 
219  double sigmaMaxSA, sigmaMaxA17, LLSA, LLA17;
220  double dmsq32SA, dmsq32A17, ssth23SA, ssth23A17;
221  unsigned int nPoints;
222 
223  fTree->Branch("sigmaMaxSA", &sigmaMaxSA);
224  fTree->Branch("sigmaMaxA17", &sigmaMaxA17);
225  fTree->Branch("LLSA", &LLSA);
226  fTree->Branch("LLA17", &LLA17);
227  fTree->Branch("dmsq32SA", &dmsq32SA);
228  fTree->Branch("dmsq32A17", &dmsq32A17);
229  fTree->Branch("ssth23SA", &ssth23SA);
230  fTree->Branch("ssth23A17", &ssth23A17);
231  fTree->Branch("Spectrum_SA", &hSA);
232  fTree->Branch("Cosmics_SA", &hSA_cosmic);
233  fTree->Branch("Spectrum_A17", &hA17);
234  fTree->Branch("Spectrum_A17_Q1", &hA17_Q1);
235  fTree->Branch("Spectrum_A17_Q2", &hA17_Q2);
236  fTree->Branch("Spectrum_A17_Q3", &hA17_Q3);
237  fTree->Branch("Spectrum_A17_Q4", &hA17_Q4);
238  fTree->Branch("Cosmics_A17", &hA17_cosmic);
239  fTree->Branch("Cosmics_A17_Q1", &hA17_Q1_cosmic);
240  fTree->Branch("Cosmics_A17_Q2", &hA17_Q2_cosmic);
241  fTree->Branch("Cosmics_A17_Q3", &hA17_Q3_cosmic);
242  fTree->Branch("Cosmics_A17_Q4", &hA17_Q4_cosmic);
243  fTree->Branch("Spectrum_Both", &hBoth);
244  fTree->Branch("Cosmics_Both", &hBoth_cosmic);
245 
246  unsigned int nEntries = inputTree->GetEntries();
247 
248  double MCpoints = prediction_OR->Predict(&calc).ToTH1(kSecondAnaPOT)->Integral();
249 
250  for(unsigned int i = 0; i < nExperiments; i++)
251  {
252  hSA->Reset();
253  hSA_cosmic->Reset();
254  hA17->Reset();
255  hA17_Q1->Reset();
256  hA17_Q2->Reset();
257  hA17_Q3->Reset();
258  hA17_Q4->Reset();
259  hA17_cosmic->Reset();
260  hA17_Q1_cosmic->Reset();
261  hA17_Q2_cosmic->Reset();
262  hA17_Q3_cosmic->Reset();
263  hA17_Q4_cosmic->Reset();
264  hBoth->Reset();
265  hBoth_cosmic_Q1->Reset();
266  hBoth_cosmic_Q2->Reset();
267  hBoth_cosmic_Q3->Reset();
268  hBoth_cosmic_Q4->Reset();
269  hBoth_cosmic->Reset();
270  hBoth_cosmic_rebin->Reset();
271 
272  nPoints = r.Poisson(MCpoints);
273 
274  unsigned int k = 0;
275  while(k <= nPoints)
276  {
277  unsigned int iEntry = r.Integer(nEntries);
278  inputTree->GetEntry(iEntry);
279 
280  double p = r.Uniform();
281  double prob = calc.P(14, 14, TrueE) * TMath::Min(CXW * PPFX * 0.5, 1.); // adjusting for weights
282 
283  if(p < prob && TrueE > 0.0 && TrueE < 5.0)
284  {
285  if(isSA)
286  {
287  hSA->Fill(CCE);
288  }
289  if(isA17)
290  {
291  int qi = quantile(QuantBins, CCE, HadEfrac);
292  if(qi == 0) hA17_Q1->Fill(CCE);
293  else if(qi == 1) hA17_Q2->Fill(CCE);
294  else if(qi == 2) hA17_Q3->Fill(CCE);
295  else hA17_Q4->Fill(CCE);
296  }
297  if(isA17 && isSA) hBoth->Fill(CCE);
298  if(isSA || isA17) k++;
299  }
300  } // loop over points
301 
302  // Protection against errors reading the tree
303  if(hSA->FindFirstBinAbove(1) == hSA->FindLastBinAbove(1)) continue;
304  if(hBoth->FindFirstBinAbove(1) == hBoth->FindLastBinAbove(1)) continue;
305  if(hA17_Q1->FindFirstBinAbove(1) == hA17_Q1->FindLastBinAbove(1)) continue;
306  if(hA17_Q2->FindFirstBinAbove(1) == hA17_Q2->FindLastBinAbove(1)) continue;
307  if(hA17_Q3->FindFirstBinAbove(1) == hA17_Q3->FindLastBinAbove(1)) continue;
308  if(hA17_Q4->FindFirstBinAbove(1) == hA17_Q4->FindLastBinAbove(1)) continue;
309 
310  // Now handle the cosmics
311  getPoissonRandom(hSpecCosmics_Both_Q1, hBoth_cosmic_Q1);
312  getPoissonRandom(hSpecCosmics_Both_Q2, hBoth_cosmic_Q2);
313  getPoissonRandom(hSpecCosmics_Both_Q3, hBoth_cosmic_Q3);
314  getPoissonRandom(hSpecCosmics_Both_Q4, hBoth_cosmic_Q4);
315 
316  hBoth_cosmic->Add(hBoth_cosmic_Q1);
317  hBoth_cosmic->Add(hBoth_cosmic_Q2);
318  hBoth_cosmic->Add(hBoth_cosmic_Q3);
319  hBoth_cosmic->Add(hBoth_cosmic_Q4);
320 
321  hBoth->Add(hBoth_cosmic);
322 
323  getSArebin(hBoth_cosmic, hBoth_cosmic_rebin);
324 
325  getPoissonRandom(hSpecCosmics_SA_rebin, hSA_cosmic);
326  hSA_cosmic->Add(hBoth_cosmic_rebin);
327  hSA->Add(hSA_cosmic);
328 
329  getPoissonRandom(hSpecCosmics_A17_Q1, hA17_Q1_cosmic);
330  getPoissonRandom(hSpecCosmics_A17_Q2, hA17_Q2_cosmic);
331  getPoissonRandom(hSpecCosmics_A17_Q3, hA17_Q3_cosmic);
332  getPoissonRandom(hSpecCosmics_A17_Q4, hA17_Q4_cosmic);
333 
334  hA17_Q1_cosmic->Add(hBoth_cosmic_Q1);
335  hA17_Q2_cosmic->Add(hBoth_cosmic_Q2);
336  hA17_Q3_cosmic->Add(hBoth_cosmic_Q3);
337  hA17_Q4_cosmic->Add(hBoth_cosmic_Q4);
338 
339  hA17_Q1->Add(hA17_Q1_cosmic);
340  hA17_Q2->Add(hA17_Q2_cosmic);
341  hA17_Q3->Add(hA17_Q3_cosmic);
342  hA17_Q4->Add(hA17_Q4_cosmic);
343 
344  hA17->Add(hA17_Q1);
345  hA17->Add(hA17_Q2);
346  hA17->Add(hA17_Q3);
347  hA17->Add(hA17_Q4);
348 
349  hA17_cosmic->Add(hA17_Q1_cosmic);
350  hA17_cosmic->Add(hA17_Q2_cosmic);
351  hA17_cosmic->Add(hA17_Q3_cosmic);
352  hA17_cosmic->Add(hA17_Q4_cosmic);
353 
354  // And finally make the fits
355  FitResult fitSA = getBestFit(prediction_SA, hSA, hSpecCosmics_SA_and_both_rebin);
356 
357  ssth23SA = fitSA.bestX;
358  dmsq32SA = fitSA.bestY;
359  LLSA = fitSA.LL;
360  sigmaMaxSA = fitSA.sigma;
361 
362  std::cout << "SA. rejection for expt " << i+1 << " " << sigmaMaxSA << std::endl;
363 
364  FitResult fitA17 = getBestFit(prediction_Q1, prediction_Q2, prediction_Q3, prediction_Q4,
365  hA17_Q1, hA17_Q2, hA17_Q3, hA17_Q4,
366  hSpecCosmics_Q1, hSpecCosmics_Q2, hSpecCosmics_Q3, hSpecCosmics_Q4);
367 
368  ssth23A17 = fitA17.bestX;
369  dmsq32A17 = fitA17.bestY;
370  LLA17 = fitA17.LL;
371  sigmaMaxA17 = fitA17.sigma;
372 
373  std::cout << "--- A17 rejection for expt " << i+1 << " " << sigmaMaxA17 << std::endl;
374 
375  fTree->Fill();
376 
377  } // loop over entries
378 
379  fTree->Write();
380  fOutput->Close();
381 
382 } // End of function
383 
385  TH1F* h1, TH1F* h2, TH1F* h3, TH1F* h4,
386  TH1D* hCosmic1, TH1D* hCosmic2, TH1D* hCosmic3, TH1D* hCosmic4)
387 {
388  osc::OscCalcPMNSOpt calc2;
390 
391  Spectrum fakeData1(h1, kSecondAnaPOT, kSecondAnaLivetime);
392  Spectrum fakeData2(h2, kSecondAnaPOT, kSecondAnaLivetime);
393  Spectrum fakeData3(h3, kSecondAnaPOT, kSecondAnaLivetime);
394  Spectrum fakeData4(h4, kSecondAnaPOT, kSecondAnaLivetime);
395 
396  SingleSampleExperiment fakeExpt1 (quant1, fakeData1, hCosmic1);
397  SingleSampleExperiment fakeExpt2 (quant2, fakeData2, hCosmic2);
398  SingleSampleExperiment fakeExpt3 (quant3, fakeData3, hCosmic3);
399  SingleSampleExperiment fakeExpt4 (quant4, fakeData4, hCosmic4);
400 
401  std::vector<const IExperiment*> vExpt;
402  vExpt.push_back(&fakeExpt1);
403  vExpt.push_back(&fakeExpt2);
404  vExpt.push_back(&fakeExpt3);
405  vExpt.push_back(&fakeExpt4);
406  vExpt.push_back(&kSolarConstraintsPDG2017);
407  vExpt.push_back(WorldReactorConstraint2017());
408  MultiExperiment fakeExpt(vExpt);
409 
410  std::vector<const IFitVar*> margPars_ssth23 = getNumuMarginalisedOscParam();
411  margPars_ssth23.push_back(&kFitDmSq32Scaled);
412 
413  MinuitFitter fitMaxMix(&fakeExpt, margPars_ssth23);
414  double MMLL = fitMaxMix.Fit(&calc2, IFitter::kQuiet)->EvalMetricVal();
415 
416  margPars_ssth23.push_back(&kFitSinSqTheta23);
417  MinuitFitter fit(&fakeExpt, margPars_ssth23);
418  double bestLL = fit.Fit(&calc2, IFitter::kQuiet)->EvalMetricVal();
419 
420  double bestX = pow(TMath::Sin(calc2.GetTh23()),2);
421  double bestY = 1e3 * calc2.GetDmsq32();
422 
423  double MMrej = 0;
424  if(MMLL > bestLL) MMrej = sqrt(MMLL-bestLL); // this should always happen, but just in case
425 
426  FitResult result = {bestX, bestY, bestLL, MMrej};
427 
428  return result;
429 
430 }
431 
432 FitResult getBestFit(IPrediction* prediction, TH1F* h1, TH1D* hCosmic)
433 {
434  osc::OscCalcPMNSOpt calc2;
436 
438 
439  // Experiment
440  SingleSampleExperiment fakeExpt1(prediction, fakeData, hCosmic);
441 
442  std::vector<const IExperiment*> vExpt;
443  vExpt.push_back(&fakeExpt1);
444  vExpt.push_back(&kSolarConstraintsPDG2017);
445  vExpt.push_back(WorldReactorConstraint2017());
446  MultiExperiment fakeExpt(vExpt);
447 
448  std::vector<const IFitVar*> margPars_ssth23 = getNumuMarginalisedOscParam();
449  margPars_ssth23.push_back(&kFitDmSq32Scaled);
450 
451  MinuitFitter fitMaxMix(&fakeExpt, margPars_ssth23);
452  double MMLL = fitMaxMix.Fit(&calc2, IFitter::kQuiet)->EvalMetricVal();
453 
454  margPars_ssth23.push_back(&kFitSinSqTheta23);
455  MinuitFitter fit(&fakeExpt, margPars_ssth23);
456  double bestLL = fit.Fit(&calc2, IFitter::kQuiet)->EvalMetricVal();
457 
458  double bestX = pow(TMath::Sin(calc2.GetTh23()),2);
459  double bestY = 1e3 * calc2.GetDmsq32();
460 
461  double MMrej = 0;
462  if(MMLL > bestLL) MMrej = sqrt(MMLL-bestLL); // this should always happen, but just in case
463 
464  FitResult result = {bestX, bestY, bestLL, MMrej};
465 
466  return result;
467 }
468 
469 std::vector<std::vector<double>> QuantileVectorFromTH2(TH2* quantileHist,
470  const HistAxis& independentAxis,
471  const HistAxis& quantileAxis,
472  const unsigned int& numQuantiles){
473 
474  std::vector<std::vector<double>> quantileBins = GetQuantileBins(quantileHist,
475  independentAxis, quantileAxis, numQuantiles);
476 
477  return quantileBins;
478 }
479 
480 unsigned int quantile(std::vector<std::vector<double>> v, float E, float hadEfrac)
481 {
482  float A17Ebins[20] = { 0, 0.75, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8,
483  1.9, 2, 2.25, 2.5, 2.75, 3, 3.5, 4, 5 };
484  unsigned int quant = 0;
485 
486  unsigned int k = 0;
487 
488  while (A17Ebins[k] < E && k < 20)
489  // for(unsigned int k = 0; k < 20; k++)
490  {
491  unsigned int l = 0;
492  while (v.at(k)[l] < hadEfrac && l < v.at(k).size())
493  // for(unsigned int l = 0; l < v.at(k).size(); l++)
494  {
495  quant = l;
496  l++;
497  }
498  k++;
499  }
500 
501  return quant;
502 }
503 
504 void getPoissonRandom(TH1D* s, TH1F* h)
505 {
506  for(int i = 1; i <= s->GetNbinsX(); i++)
507  {
508  float fX = s->GetBinCenter(i);
509  float fY = r.Poisson(s->GetBinContent(i));
510 
511  h->Fill(fX, fY);
512  }
513 }
514 
515 void getSArebin(TH1D* s, TH1D* h)
516 {
517  for(int i = 1; i <= s->GetNbinsX(); i++)
518  {
519  double fX = s->GetBinCenter(i);
520  double fY = s->GetBinContent(i);
521 
522  h->Fill(fX, fY);
523  }
524 }
525 
526 void getSArebin(TH1F* s, TH1F* h)
527 {
528  for(int i = 1; i <= s->GetNbinsX(); i++)
529  {
530  float fX = s->GetBinCenter(i);
531  float fY = s->GetBinContent(i);
532 
533  h->Fill(fX, fY);
534  }
535 }
double LL
Definition: Toy_analyses.C:60
Implements systematic errors by interpolation between shifted templates.
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
size_t nPoints
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
void getPoissonRandom(TH1D *s, TH1F *h)
Definition: Toy_analyses.C:504
TH1F * h3
Definition: berger.C:36
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
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
virtual T GetTh23() const
Definition: IOscCalc.h:94
std::vector< std::vector< double > > GetQuantileBins(TH2 *quantileHist, const HistAxis &independentAxis, const HistAxis &quantileAxis, const unsigned int &numQuantiles, const bool verbose)
Returns a 2D vector First index is the independentAxis bin number Second index is the high bin edge f...
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: UtilsExt.cxx:237
osc::OscCalcDumb calc
double sigma
Definition: Toy_analyses.C:61
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const XML_Char * s
Definition: expat.h:262
FitResult getBestFit(IPrediction *prediction, TH1F *h1, TH1D *hCosmic)
Definition: Toy_analyses.C:432
const double kAna2017Livetime
Definition: Exposures.h:200
Float_t E
Definition: plot.C:20
virtual std::unique_ptr< IFitSummary > Fit(osc::IOscCalcAdjustable *seed, SystShifts &bestSysts=junkShifts, const SeedList &seedPts=SeedList(), const std::vector< SystShifts > &systSeedPts={}, Verbosity verb=kVerbose) const
Master fitting method. Depends on FitHelper and FitHelperSeeded.
Definition: IFitter.cxx:69
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
Definition: structs.h:1
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
T P(int flavBefore, int flavAfter, double E) override
E in GeV; flavors as PDG codes (so, neg==>antinu)
double bestX
Definition: Toy_analyses.C:58
void Toy_analyses(unsigned int nExperiments, unsigned int volume=1)
Definition: Toy_analyses.C:79
void SetTh23(const T &th23) override
virtual T GetDmsq32() const
Definition: IOscCalc.h:91
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
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
unsigned int quantile(std::vector< std::vector< double >> v, float E, float hadEfrac)
Definition: Toy_analyses.C:480
OStream cout
Definition: OStream.cxx:6
TH1F * h2
Definition: plot.C:45
void SetDmsq32(const T &dmsq32) override
Spectrum Predict(osc::IOscCalc *calc) const override
TH1F * h1
Combine multiple component experiments.
double bestY
Definition: Toy_analyses.C:59
const SolarConstraints kSolarConstraintsPDG2017(7.53e-5, 0.18e-5, 0.851, 0.020)
void getSArebin(TH1D *s, TH1D *h)
Definition: Toy_analyses.C:515
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
TRandom3 r(0)
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
std::vector< std::vector< double > > QuantileVectorFromTH2(TH2 *quantileHist, const HistAxis &independentAxis, const HistAxis &quantileAxis, const unsigned int &numQuantiles)
Definition: Toy_analyses.C:469
Float_t e
Definition: plot.C:35
const double kSecondAnaPOT
Definition: Exposures.h:87
T asin(T number)
Definition: d0nt_math.hpp:60
const double kSecondAnaLivetime
Definition: Exposures.h:113
Compare a single data spectrum to the MC + cosmics expectation.
std::vector< const IFitVar * > getNumuMarginalisedOscParam()
Definition: FitVarsNumu.cxx:9
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17
enum BeamMode string