sensitivity2020.C
Go to the documentation of this file.
6 #include "CAFAna/Core/Spectrum.h"
7 #include "CAFAna/Core/Registry.h"
12 #include "CAFAna/FC/FCSurface.h"
13 #include "CAFAna/Fit/IFitter.h"
17 #include "CAFAna/Systs/Systs.h"
23 #include "CAFAna/Vars/FitVars.h"
24 
26 #include "OscLib/IOscCalc.h"
27 #include "Utilities/rootlogon.C"
28 
29 #include "TFile.h"
30 #include "TStyle.h"
31 #include "TGraph.h"
32 #include "TCanvas.h"
33 #include "TSystem.h"
34 #include "TH2.h"
35 
36 using namespace ana;
37 
38 // MakeFile: Do you need to make the intermediate files, or just load from them?
39 // MakeSurf: true false
40 // whatsurf: "th23dcp", "dmsqth23"
41 // MakeSens: true false
42 // whatplot: "mh", "dcp", "maxmix", "oct"
43 // isSyst : use syst, not only nominals
44 // NuMuOnly: true false --> Run over only numu, or include nue too?
45 // horn : "both", "fhc", "rhc"
46 // isNH : true false
47 void sensitivity2020(bool MakeFile = false,
48  bool MakeSurf = false,
49  std::string whatsurf = "th23dcp",
50  bool MakeSens = false,
51  bool Varyth23 = false,
52  std::string whatplot = "mh",
53  bool isSyst = false,
54  bool isPtExtrap = false,
55  bool NuMuOnly = false,
56  std::string horn = "both",
57  bool isNH = true,
58  std::string option = "2019BestFit",
59  bool OnGrid = false,
60  bool isFutureSens = false ) {
61 
62  // --- Set some constants that I'll want to use.
63  const int nnumu = 4;
64  const std::vector<double> delta_seeds = {0, 0.5, 1., 1.5};
65  const std::vector<double> th23_seeds = {0.45, 0.55};
66  ////////////////////////////////////////////////////////////////
67  /////////////////////// Load nue and numu stuff ///////////////
68  ///////////////////////////////////////////////////////////////
69 
70  std::string SystTag = isSyst ?"with_syst":"just_stat";
71  std::string sNuMu = NuMuOnly ?"numuonly" :"numu_nue";
72  std::string sExtrap = isPtExtrap?"PtExtrap" :"nonPtExtrap";
73  std::string sHier = isNH ?"NH" :"IH";
74  std::string TagName = sNuMu+"_"+horn+"_"+sHier+"_"+SystTag+"_"+sExtrap+"_"+option;
75  std::string outfilename = "run_sensitivity_"+TagName;
76 
77  std::cout << "\n Greetings! You requested to " << (MakeFile?" make files." : " only plot." )
78  << "\nIt will use " << horn << " files and " << SystTag << " for " << sNuMu << ", plotting: "
79  << (MakeSurf?"contours ":"no contours ")
80  << (isNH ?" NH " :" IH " ) << " ; " << whatsurf
81  << (MakeSens?" slices " :" no slices" )
82  <<" in case of slices it will be "<< whatplot
83  << "\nAll using option " << option
84  << "\n\t TagName is thus " << TagName
85  << std::endl;
86 
88  // Some options for different best fits.
89  double BASE_dcp = 0.82*M_PI;
90  double BASE_th23 = 0.568;
91  double BASE_dmsq32 = 2.41e-3;
92  std::cout << "\n Default calc is set to 2020 BF." << std::endl;
93  // --- 2019 best fit
94  if ( option.find("2019BestFit") != std::string::npos ) {
95  BASE_dcp = 2*M_PI;
96  BASE_th23 = 0.565;
97  BASE_dmsq32 = 2.48e-3;
98  std::cout << "\t Updated calc to 2019 BF." << std::endl;
99  }
100  // --- 2019 Best Fit - NH, LO
101  if ( option.find("NHLO") != std::string::npos ) {
102  BASE_dcp = 1.9*M_PI;
103  BASE_th23 = 0.480;
104  BASE_dmsq32 = 2.49e-3;
105  std::cout << "\t Updated calc to NH LO 2019 BF." << std::endl;
106  }
107  // --- 2019 Best Fit - IH UO
108  if ( option.find("IHUO") != std::string::npos ) {
109  BASE_dcp = 1.5*M_PI;
110  BASE_th23 = 0.560;
111  BASE_dmsq32 = -2.54e-3;
112  std::cout << "\t Updated calc to IH UO 2019 BF." << std::endl;
113  }
114  calc ->SetdCP ( BASE_dcp );
115  calc ->SetTh23 ( asin(sqrt(BASE_th23)) );
116  calc ->SetDmsq32( BASE_dmsq32 );
117  // --- 2018 best fit
118  //calc->SetdCP(1.21*M_PI);
119  //calc->SetTh23(asin(sqrt(0.558)));
120  //calc->SetDmsq32(2.44e-3);
121 
122  std::vector<double> pot;
123  std::vector<double> livetime;
124 
125  double temp_potFD_rhc = kAna2020RHCPOT, temp_cosFD_rhc = kAna2020RHCLivetime;
126  double temp_potFD_fhc = kAna2020FHCPOT, temp_cosFD_fhc = kAna2020FHCLivetime;
127  if (isFutureSens) {
128  temp_potFD_rhc = kFutureRHCPOT; temp_cosFD_rhc = kFutureRHCLivetime;
129  temp_potFD_fhc = kFutureFHCPOT; temp_cosFD_fhc = kFutureFHCLivetime;
130  }
131  const double potFD_rhc = temp_potFD_rhc;
132  const double potFD_fhc = temp_potFD_fhc;
133  const double cosFD_rhc = temp_cosFD_rhc;
134  const double cosFD_fhc = temp_cosFD_fhc;
135 
136  std::vector<const IPrediction*> predsvector;
137  std::vector<std::pair <Spectrum*, double>> cosmicsvector;
138 
139  std::vector<const IPrediction*> predictionsNue;
140  std::vector<const IPrediction*> predictionsNumu;
141  std::vector<std::pair <Spectrum*, double>> cosmicsNue;
142  std::vector<std::pair <Spectrum*, double>> cosmicsNumu;
143 
144  string decomp_fhc = "combo"; //should be combo
145  string decomp_rhc = "prop";
146  if (!isPtExtrap) {
147  decomp_fhc += "_noPt";
148  decomp_rhc += "_noPt";
149  outfilename += "_noPtExtrap";
150  }
151 
152  bool FakeNDData = false;
153 
154  // --- Load the Nu_e things.
155  if (!NuMuOnly) {
156  std::cout << "\n\t Not NuMuOnly, so pushing back the nue stuff...\n" << std::endl;
157  if (horn=="fhc" || horn=="both") {
158  predictionsNue .push_back( GetNuePrediction2020(decomp_fhc, DefaultOscCalc(), isSyst, "fhc", FakeNDData) );
159  cosmicsNue .push_back( GetNueCosmics2020("fhc") );
160  }
161  if(horn=="rhc" || horn=="both") {
162  predictionsNue.push_back( GetNuePrediction2020(decomp_rhc, DefaultOscCalc(), isSyst, "rhc", FakeNDData) );
163  cosmicsNue .push_back( GetNueCosmics2020("rhc") );
164  }
165 
166  // --- Push the Nue stuff into the Cosmics and Preds vectors
167  cosmicsvector.insert(cosmicsvector.end(), cosmicsNue .begin(), cosmicsNue .end());
168  predsvector .insert(predsvector .end(), predictionsNue.begin(), predictionsNue.end());
169  }
170 
171  // --- Load the Nu_mu things
172  if (horn=="fhc" || horn=="both") {
173  auto numu_preds_fhc = GetNumuPredictions2020(nnumu, decomp_fhc, DefaultOscCalc(), isSyst, "fhc", FakeNDData);
174  auto numu_cosmics_fhc = GetNumuCosmics2020 (nnumu, "fhc");
175  predictionsNumu.insert(predictionsNumu.end(), numu_preds_fhc .begin(), numu_preds_fhc .end() );
176  cosmicsNumu .insert(cosmicsNumu .end(), numu_cosmics_fhc.begin(), numu_cosmics_fhc.end() );
177  }
178  if (horn=="rhc" || horn=="both") {
179  auto numu_preds_rhc = GetNumuPredictions2020(nnumu, decomp_rhc, DefaultOscCalc(), isSyst, "rhc", FakeNDData);
180  auto numu_cosmics_rhc = GetNumuCosmics2020 (nnumu, "rhc");
181  predictionsNumu.insert( predictionsNumu.end(), numu_preds_rhc .begin(), numu_preds_rhc .end() );
182  cosmicsNumu .insert( cosmicsNumu .end(), numu_cosmics_rhc.begin(), numu_cosmics_rhc.end() );
183  }
184  // --- Push the Numu stuff into the Cosmics and Preds vectors
185  cosmicsvector.insert(cosmicsvector.end(), cosmicsNumu .begin(), cosmicsNumu .end());
186  predsvector .insert(predsvector .end(), predictionsNumu.begin(), predictionsNumu.end());
187 
188  std::cout<<"\n\n Good job, you've loaded everything! \n\n"<< std::endl;
189 
190  //////////////////////////////////////////////////////////////
191  ///////////////////// Check Predictions at first /////////////
192  /////////////////////////////////////////////////////////////
193  if (horn=="fhc" || horn=="both") {
194  pot.push_back(potFD_fhc);
195  livetime.push_back(cosFD_fhc);
196  }
197  if (horn=="rhc" || horn=="both") {
198  pot.push_back(potFD_rhc);
199  livetime.push_back(cosFD_rhc);
200  }
201 
202  std::vector<TH1*> hNue; // nue signal
203  std::vector<TH1*> hMC; //total MC
204  std::vector<TH1*> hTotbkg; //total sum of all bkg
205  std::vector<TH1*> hBkg; // all bkg except cosmics
206  std::vector<TH1*> hNC; //1
207  std::vector<TH1*> hCC; //2
208  std::vector<TH1*> hBeam; //3
209  std::vector<TH1*> hTau; //4
210  std::vector<TH1*> hCos; //5
211  std::vector<TH1*> hRock; //6 rock events
212  std::vector<TH1*> hNue_bg; //7 wrong sign nue
213 
214  std::vector<TH1*> hNumu;
215  std::vector<TH1*> hCosNumu;
216  std::vector<TH1*> hMCNumu;
217 
218  std::cout<<">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> Here are the predictions for 2020 nue ana with 2019 Best fit values: "<< std::endl;
219  for(int i = 0; i < (int)predictionsNue.size(); i++){
221  if (i==0){
222  sFHC = "FHC";
223  hNue .push_back(predictionsNue[i]->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kNu ).ToTH1(pot[i])); //signal FHC
224  hNue_bg.push_back(predictionsNue[i]->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kAntiNu).ToTH1(pot[i])); //wrong sign bkg for FHC
225  }
226  if (i==1){
227  sFHC = "RHC";
228  hNue .push_back(predictionsNue[i]->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kAntiNu).ToTH1(pot[i])); //signal RHC
229  hNue_bg.push_back(predictionsNue[i]->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kNu ).ToTH1(pot[i])); //wrong sign bkg for RHC
230  }
231 
232  hTau .push_back(predictionsNue[i]->PredictComponent(calc, Flavors::kAllNuTau, Current::kCC, Sign::kBoth).ToTH1(pot[i])); // tau CC bkg (2)
233  hBeam.push_back(predictionsNue[i]->PredictComponent(calc, Flavors::kNuEToNuE, Current::kCC, Sign::kBoth).ToTH1(pot[i])); // beam bkg (3)
234  hNC .push_back(predictionsNue[i]->PredictComponent(calc, Flavors::kAll , Current::kNC, Sign::kBoth).ToTH1(pot[i])); // NC bkg (4)
235  hCC .push_back(predictionsNue[i]->PredictComponent(calc, Flavors::kAllNuMu , Current::kCC, Sign::kBoth).ToTH1(pot[i])); // numu CC bkg (5)
236  hMC .push_back(predictionsNue[i]->Predict(calc).ToTH1(pot[i]));
237  TH1D* temp = (TH1D*)hNue_bg[i]->Clone(UniqueName().c_str());
238  temp->Add(hTau [i]);
239  temp->Add(hBeam[i]);
240  temp->Add(hNC [i]);
241  temp->Add(hCC [i]);
242  hBkg.push_back( (TH1D*)temp->Clone(UniqueName().c_str()));
243  TH1D *tempCos = cosmicsNue[i].first->ToTH1(livetime[i], kBlack, kSolid, kLivetime);
244  hCos.push_back( tempCos ); //Clone(UniqueName().c_str())); //cosmic bkg (6)
245  std::cerr << "Done the cosmics..." << std::endl;
246  temp ->Add(hCos[i]);
247  hMC[i]->Add(hCos[i]);
248  hTotbkg.push_back(temp);
249 
250  std::cout << sFHC << " mode, with " << pot[i] <<" POT, and " << livetime[i] << " s livetime."
251  <<"\n\t Total MC " << hMC [i]->Integral()
252  <<"\n\t Nu_e Signal " << hNue [i]->Integral()
253  <<"\n\t Total Bkg " << hTotbkg[i]->Integral()
254  <<"\n\t Total Beam bkg " << hBkg [i]->Integral()
255  <<"\n\t Wrong Sign bkg " << hNue_bg[i]->Integral()
256  <<"\n\t Beam Nu_e " << hBeam [i]->Integral()
257  <<"\n\t NC " << hNC [i]->Integral()
258  <<"\n\t Nu_mu CC " << hCC [i]->Integral()
259  <<"\n\t Nu_tau CC " << hTau [i]->Integral()
260  <<"\n\t Cosmic " << hCos [i]->Integral()
261  << std::endl;
262  }
263 
264  std::cout<<">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> Here are the predictions for 2020 numu ana with 2019 Best fit values: "<< std::endl;
265  double AllQ_NuMu_Full, AllQ_Beam_Full, AllQ_Cosm_Full;
266  double AllQ_NuMu_Dip , AllQ_Beam_Dip , AllQ_Cosm_Dip;
267  AllQ_NuMu_Full = AllQ_Beam_Full = AllQ_Cosm_Full = AllQ_NuMu_Dip = AllQ_Beam_Dip = AllQ_Cosm_Dip = 0;
268  for(int i = 0; i < (int)predictionsNumu.size(); i++){
269  std::string sFHC = "FHC";
270  double pot = potFD_fhc;
271  double live = cosFD_fhc;
272  int quant = i+1;
273  if (i>=4) {
274  sFHC = "RHC";
275  pot = potFD_rhc;
276  live = cosFD_rhc;
277  quant= quant-4;
278  }
279 
280  hNumu .push_back(predictionsNumu[i]->PredictComponent(calc, Flavors::kNuMuToNuMu, Current::kCC, Sign::kBoth).ToTH1(pot));
281  hMCNumu .push_back(predictionsNumu[i]->Predict(calc).ToTH1(pot));
282  TH1D *tempCos = cosmicsNumu[i].first->ToTH1(live, kBlack, kSolid, kLivetime);
283  hCosNumu.push_back( tempCos ); //Clone(UniqueName().c_str()));
284 
285  double NuMu_Full = hNumu [i]->Integral();
286  double Beam_Full = hMCNumu [i]->Integral() - NuMu_Full;
287  double Cosm_Full = hCosNumu[i]->Integral();
288  double FOM_Full = NuMu_Full / pow( NuMu_Full+Beam_Full+Cosm_Full, 0.5 );
289  int bmin = hNumu[i]->GetXaxis()->FindBin(1.);
290  int bmax = hNumu[i]->GetXaxis()->FindBin(2.);
291  double NuMu_Dip = hNumu [i]->Integral( bmin, bmax );
292  double Beam_Dip = hMCNumu [i]->Integral( bmin, bmax ) - NuMu_Dip;
293  double Cosm_Dip = hCosNumu[i]->Integral( bmin, bmax );
294  double FOM_Dip = NuMu_Dip / pow( NuMu_Dip+Beam_Dip+Cosm_Dip, 0.5 );
295  std::cout << sFHC << " mode, Quartile " << quant << " with " << pot <<" POT, and " << live << " s livetime."
296  << "\n\t Nu_mu Signal Full: " << NuMu_Full << "\t Dip: " << NuMu_Dip
297  << "\n\t Beam Bkg Full: " << Beam_Full << "\t Dip: " << Beam_Dip
298  << "\n\t Cosmic Full: " << Cosm_Full << "\t Dip: " << Cosm_Dip
299  << "\n\t FOM_Full Full: " << FOM_Full << "\t Dip: " << FOM_Dip
300  << std::endl;
301  AllQ_NuMu_Full += NuMu_Full; AllQ_NuMu_Dip += NuMu_Dip;
302  AllQ_Beam_Full += Beam_Full; AllQ_Beam_Dip += Beam_Dip;
303  AllQ_Cosm_Full += Cosm_Full; AllQ_Cosm_Dip += Cosm_Dip;
304  if (quant%4==0) {
305  double AllQ_FOM_Full = AllQ_NuMu_Full / pow( AllQ_NuMu_Full+AllQ_Beam_Full+AllQ_Cosm_Full, 0.5 );
306  double AllQ_FOM_Dip = AllQ_NuMu_Dip / pow( AllQ_NuMu_Dip +AllQ_Beam_Dip +AllQ_Cosm_Dip , 0.5 );
307  std::cout << sFHC << " mode, All Quartile with " << pot <<" POT, and " << live << " s livetime."
308  << "\n\t Nu_mu Signal Full: " << AllQ_NuMu_Full << "\t Dip: " << AllQ_NuMu_Dip
309  << "\n\t Beam Bkg Full: " << AllQ_Beam_Full << "\t Dip: " << AllQ_Beam_Dip
310  << "\n\t Cosmic Full: " << AllQ_Cosm_Full << "\t Dip: " << AllQ_Cosm_Dip
311  << "\n\t FOM_Full Full: " << AllQ_FOM_Full << "\t Dip: " << AllQ_FOM_Dip
312  << "\n" << std::endl;
313  AllQ_NuMu_Full = AllQ_Beam_Full = AllQ_Cosm_Full = AllQ_NuMu_Dip = AllQ_Beam_Dip = AllQ_Cosm_Dip = 0;
314  }
315  }
316 
317 
318  std::cout<<"\n\n Good job, you've got your integrals, now to make sensitivites! \n\n"<< std::endl;
319  /////////////////////////////////////////////////////////////
320  /////////////////// Sensitivities ///////////////////////////
321  ////////////////////////////////////////////////////////////
322 
323  // TODO Extrap, Extrap + Syst; also nue, numu, both; RHC, FHC, both
324  // species from 2019 TN: th23 vs delta; th23 vs dmsq23; slices: hie, dcp, octant, max mixing
325 
326  //make fakedata at 2019 BF(make choosable)
327  //ResetOscCalcToDefault(calc);
328  //calc->SetdCP(1.21*M_PI);
329  //calc->SetTh23(asin(sqrt(0.558)));
330  //calc->SetDmsq32(2.44e-3);
331 
332  std::vector <const IExperiment * > expts;
333  std::cout<<"Create fake data with systematic pred. and multiexp"<< std::endl;
334  for(int i = 0; i < (int) predsvector.size(); ++i){
335  double pot = potFD_fhc;
336  double live = cosFD_fhc;
337  // Cases for when I need RHC pot/livetime.
338  if(horn=="both" && ( (!NuMuOnly && (i==1 || i>5) ) || // Nue+NuMu then vector is (nueFHC, nueRHC, numuFHC*4, numuRHC*4)
339  ( NuMuOnly && i>3 ) // NuMu only then vector is (numuFHC*4, numuRHC*4)
340  ) ) {
341  pot = potFD_rhc;
342  live = cosFD_rhc;
343  }
344  std::cout << "\t My POT is " << pot << ", My Livetime is " << live << std::endl;
345  auto thisdata = GetFakeData ( predsvector[i], calc, pot, cosmicsvector[i].first );
346  std::cout << "\t Fake data integral is: "<< thisdata->Integral(pot) << std::endl;
347  expts.push_back(new SingleSampleExperiment( predsvector[i], *thisdata, *cosmicsvector[i].first, cosmicsvector[i].second, true) );
348  }
349 
350  expts.push_back(WorldReactorConstraint2019());
351  auto exptsAll = new MultiExperiment(expts);
352  std::cout<<"Experiments are ready, number of items in vector is "<<expts.size()<<" (the above plus WordReactorConstraint)" << std::endl;
353 
354  // --- Set systematic correlations if requested --- blindly copying out for now...
355  std::cout<<"\n Now I want to set your systematics"<< std::endl;
356  std::vector< const ISyst * > systs = {};
357  SystShifts auxShifts = SystShifts::Nominal();
358  std::vector <SystShifts> seedShifts = {};
359 
360  if(isSyst){
361  std::cout << "\t Getting syst list and the correlations." << std::endl;
362  if (NuMuOnly) {
363  std::cout << "\t\t NuMuOnly, so get a reduced list of systs..." << std::endl;
364  systs = GetJointFitSystematicList(isSyst, false, true, true, true, true, isPtExtrap);
365  auto notfornumuFHC = GetCorrelations(false, true , isPtExtrap);
366  auto notfornumuRHC = GetCorrelations(false, false, isPtExtrap);
367  for (size_t i=0; i<expts.size(); ++i) {
368  if (i<4) {
369  exptsAll->SetSystCorrelations(i, notfornumuFHC);
370  } else {
371  if (i==8) continue;
372  exptsAll->SetSystCorrelations(i, notfornumuRHC);
373  }
374  }
375  } else {
376  std::cout << "\t\t NuMu and Nue, so get all of the systs..." << std::endl;
377  systs = GetJointFitSystematicList(isSyst, false, false, true, true, true, isPtExtrap);
378  auto notfornueFHC = GetCorrelations(true , true , isPtExtrap);
379  auto notfornueRHC = GetCorrelations(true , false, isPtExtrap);
380  auto notfornumuFHC = GetCorrelations(false, true , isPtExtrap);
381  auto notfornumuRHC = GetCorrelations(false, false, isPtExtrap);
382  if(horn=="fhc"){
383  exptsAll->SetSystCorrelations(0, notfornueFHC);
384  for(int i=1; i<5;i++) {
385  exptsAll->SetSystCorrelations(i, notfornumuFHC);
386  }
387  }
388  if(horn=="both"){
389  exptsAll->SetSystCorrelations(0, notfornueFHC);
390  exptsAll->SetSystCorrelations(1, notfornueRHC);
391  for(int i=0; i < 4; ++i) exptsAll->SetSystCorrelations(i+2, notfornumuFHC);
392  for(int i=0; i < 4; ++i) exptsAll->SetSystCorrelations(i+6, notfornumuRHC);
393  }
394 
395  auxShifts = SystShifts::Nominal();
396  //Registry::Print();
397  }
398  }
399  else {
400  std::cout<<"\t You're not setting any systematics"<< std::endl;
401  systs={};
402  //auxShifts = junkShifts;
403  auxShifts = SystShifts();
404  }
405 
406  // --- find BF here:
407  std::cout << "\n Now to find the Best Fit." << std::endl;
408  std::vector <const IFitVar *> fitvars = {&kFitDeltaInPiUnits, &kFitSinSqTheta23, &kFitDmSq32, &kFitSinSq2Theta13};
409  MinuitFitter fit(exptsAll, fitvars, systs);
410 
411  double minchi23 = 1E20;
412  ResetOscCalcToDefault(calc);
413  for(auto & th23seed:th23_seeds){
414  for(int hie:{-1,1}){
415  std::cout << "Finding best fit for " << (th23seed<0.5?"LO ":"UO ") << (hie>0? "NH " : "IH ") << std::endl;
416  auto thisminchi = fit.Fit( calc, auxShifts ,
417  SeedList({
418  { &kFitDmSq32 , {hie*fabs(calc->GetDmsq32())} },
419  { &kFitSinSqTheta23 , {th23seed} },
420  { &kFitDeltaInPiUnits, delta_seeds }
421  }),
422  seedShifts, IFitter::kVerbose
423  );
424  minchi23= min(minchi23, thisminchi);
425  ResetOscCalcToDefault(calc);
426  if(isSyst){
427  auxShifts = SystShifts::Nominal();
428  std::cout<<"reset syst as well"<< std::endl;
429  }
430  }//end hie
431  }//end th23
432  std::cout << "Found overall min chi " << minchi23 << "\n" << std::endl;
433 
434  std::string FileDir = "";
435  //std::string FileDir = "/nova/ana/3flavor/Ana2020/FutureSensitivities/";
436  std::string PlotDir = "";
437  if (!OnGrid) {
438  FileDir = "files/";
439  PlotDir = "plots/";
440  gSystem->MakeDirectory( FileDir.c_str() );
441  gSystem->MakeDirectory( PlotDir.c_str() );
442  }
443  // -----------------------------------
444  // contours from 2019 TN --> Make Surf
445  if(!MakeSurf){
446  std::cout << "\n NOT making surfaces.\n" << std::endl;
447  } else {
448 
449  int hie = 1;
450  if(!isNH) { hie = -1; }
451 
452  std::string FileName = FileDir+outfilename+"_Surfonly_"+whatsurf+".root";
453  std::string PlotName;
454  if (whatsurf == "th23dcp" ) PlotName = TagName+"-surf_delta_th23";
455  if (whatsurf == "dmsqth23") PlotName = TagName+"-surf_th23_dmsq23";
456 
457  std::cout << "\n You now want to make the surfaces called " << PlotName << ". They will be saved/loaded from " << FileName << std::endl;
458 
459  if(MakeFile){
460  std::cout << "You are going to make the file containing the surfaces." << std::endl;
461  TFile * outfile = new TFile( FileName.c_str(), "recreate" );
462  outfile->cd();
463  TVectorD v(1);
464  v[0] = minchi23;
465  v.Write("overall_minchi");
466 
467  std::cout<<"\t Clean the vectors"<< std::endl;
468  predsvector.clear();
469 
470  ResetOscCalcToDefault(calc);
471  calc->SetDmsq32(hie*fabs(BASE_dmsq32));
472  if(isSyst) auxShifts = SystShifts::Nominal();
473 
474  // --- th23dcp
475  if (whatsurf=="th23dcp") {
476  int Ysteps = 30; int Xsteps = 30;
477  std::cout<<"\t Doing th23dcp --> minimize with grid "<< Ysteps<<" * "<<Xsteps<<" points"<< std::endl;
478 
479  FrequentistSurface surf_th23_dcp (exptsAll, calc,
480  &kFitDeltaInPiUnits, Xsteps, 0 , 2,
481  &kFitSinSqTheta23 , Ysteps, 0.2, 0.7,
482  {&kFitDmSq32Scaled, &kFitSinSq2Theta13},
483  systs,
484  SeedList({{&kFitSinSq2Theta13, {calc->GetTh13()}}})
485  //,seedShifts
486  );
487  outfile->cd();
488  surf_th23_dcp.SaveTo(outfile, PlotName.c_str() ) ;
489  }
490  // --- dmsqth23
491  else if (whatsurf=="dmsqth23") {
492  int Ysteps = 26; int Xsteps = 26;
493  std::cout<<"\t Doing dmsqth23 --> minimize with grid "<< Ysteps<<" * "<<Xsteps<<" points"<< std::endl;
494 
495  FrequentistSurface surf_th23_dmsq23 (exptsAll, calc,
496  &kFitSinSqTheta23, Xsteps, 0.3, 0.7,
497  &kFitDmSq32Scaled, Ysteps, (hie>0?2.1:-2.8), (hie>0?2.8:-2.1),
499  systs,
500  SeedList({{&kFitDeltaInPiUnits, delta_seeds}})
501  //,seedShifts
502  );
503  outfile->cd();
504  auto surf_test=surf_th23_dmsq23.ToTH2(minchi23);
505  surf_test ->Write( TString(PlotName)+TString("_2DTest") );
506  surf_th23_dmsq23.SaveTo( outfile, PlotName.c_str() ) ;
507  }
508  outfile->Close();
509  std::cout<<"Done with surfaces!"<< std::endl;
510  }
511 
512  // ----> If not Make File, or once finished the above loop.
513  std::cout << "You are going to load premade surfaces." << std::endl;
514  TFile * infile = new TFile( FileName.c_str(),"read");
515 
516  auto mins =* (TVectorD*)infile->Get("overall_minchi");
517  double minchi23 = mins[0];
518  std::cout<<"Min chi: "<<minchi23<< std::endl;
519 
520  Int_t kFillColor1 = TColor::GetColorTransparent(k1SigmaConfidenceColorNH, 0.70);
521  Int_t kFillColor2 = TColor::GetColorTransparent(k2SigmaConfidenceColorNH, 0.35);
522  Int_t kFillColor3 = TColor::GetColorTransparent(k3SigmaConfidenceColorNH, 0.20);
523  if (!isNH) {
524  kFillColor1 = TColor::GetColorTransparent(k1SigmaConfidenceColorIH, 0.70);
525  kFillColor2 = TColor::GetColorTransparent(k2SigmaConfidenceColorIH, 0.35);
526  kFillColor3 = TColor::GetColorTransparent(k3SigmaConfidenceColorIH, 0.20);
527  }
528  auto leg = ContourLegend(hie, false, false, kFillColor1, kFillColor2, kFillColor3, kBlack, false);
529 
530  // --- th23dcp
531  if (whatsurf=="th23dcp"){
532  auto surf_th23_dcp = *FrequentistSurface::LoadFrom(infile, PlotName.c_str() ) ;
533 
534  TH2 * surf_th23_dcp_1Sigma = Gaussian68Percent2D(surf_th23_dcp);
535  TH2 * surf_th23_dcp_2Sigma = Gaussian2Sigma2D (surf_th23_dcp);
536  TH2 * surf_th23_dcp_3Sigma = Gaussian3Sigma2D (surf_th23_dcp);
537  DrawContourCanvas("delta_th23", 0., 2., 0.225, 0.725);
538  surf_th23_dcp.DrawContour(surf_th23_dcp_1Sigma, kSolid , kFillColor1, minchi23);
539  surf_th23_dcp.DrawContour(surf_th23_dcp_2Sigma, k2SigmaConfidenceStyle, kFillColor2, minchi23);
540  surf_th23_dcp.DrawContour(surf_th23_dcp_3Sigma, kSolid , kFillColor3, minchi23);
541 
542  leg->Draw();
543  Simulation();
544 
545  gPad->Print((TString)PlotDir+PlotName+".pdf");
546  }
547  // --- dmsqth23
548  else if (whatsurf=="dmsqth23") {
549  auto surf_th23_dmsq23 = *FrequentistSurface::LoadFrom(infile, PlotName.c_str() );
550 
551  TH2 * surf_th23_dmsq23_1Sigma = Gaussian68Percent2D(surf_th23_dmsq23);
552  TH2 * surf_th23_dmsq23_2Sigma = Gaussian2Sigma2D (surf_th23_dmsq23);
553  TH2 * surf_th23_dmsq23_3Sigma = Gaussian3Sigma2D (surf_th23_dmsq23);
554 
555  if (hie>0) DrawContourCanvas("th23_dm32", 0.3, 0.7, 2.05, 2.75);
556  else DrawContourCanvas("th23_dm32", 0.3, 0.7, -2.95, -2.25);
557 
558  surf_th23_dmsq23.DrawContour(surf_th23_dmsq23_1Sigma, kSolid , kFillColor1, minchi23);
559  surf_th23_dmsq23.DrawContour(surf_th23_dmsq23_2Sigma, k2SigmaConfidenceStyle, kFillColor2, minchi23);
560  surf_th23_dmsq23.DrawContour(surf_th23_dmsq23_3Sigma, kSolid , kFillColor3, minchi23);
561 
562  leg->Draw();
563  Simulation();
564 
565  gPad->Print((TString)PlotDir+PlotName+".pdf");
566  }
567  } //MakeSurf end
568  // contours from 2019 TN --> Make Surf
569  // -----------------------------------
570 
571  // -----------------------------------
572  // Now for the Make Sens loop
573  if(!MakeSens){
574  std::cout << "\n NOT making any sensitivities.\n" << std::endl;
575  } else {
576  std::vector<std::pair<TGraph*, TString>> slice;
577  std::vector<std::pair<TGraph*, TString>> slice2019;
578  slice.clear();
579  slice2019.clear();
580 
581  struct th23helper{
582  std::vector<double> seeds;
583  const IFitVar * var;
584  TString label;
585  };
586  std::vector <th23helper> th23seeds;
587  th23seeds.push_back( { {0.4,0.6}, &kFitSinSqTheta23 , "2019" } );
588  th23seeds.push_back( { {0.45} , &kFitSinSqTheta23LowerOctant, "LO SA"} );
589  th23seeds.push_back( { {0.55} , &kFitSinSqTheta23UpperOctant, "UO SA"} );
590 
591  int id, revid;
592  std::string sVary = Varyth23 ? "Varyth23":"VarydCP";
593  std::string FileName = FileDir+outfilename+"_"+sVary+"_Slicesonly_"+whatplot+".root";
594 
595  if(MakeFile){
596  std::cout << "\nYou are going to make the file containing the sensitivites." << std::endl;
597  TFile * outfile = new TFile( FileName.c_str(), "recreate" );
598  outfile->cd();
599 
600  TVectorD v(1);
601  v[0] = minchi23;
602  v.Write("overall_minchi");
603 
604  th23helper th23seed = th23seeds[0];
605 
606  double th23 = BASE_th23;
607  double delta = BASE_dcp;
608 
609  for(int hie : {-1,1}) {
610  std::cout << "\nNow looking at hie = " << hie << std::endl;
611  int steps = 25; double chi2r[steps], chi2w[steps], chi[steps]; double Fits[steps];
612  double step = 2.0/(steps-1);
613  for (int i =0; i<steps; i++){
614  if ( Varyth23 ) { th23 = (i*step*0.1 ) + 0.4; }
615  else { delta = (i*step*M_PI); }
616 
617  ResetOscCalcToDefault(calc);
618  calc ->SetdCP ( delta );
619  calc ->SetTh23 ( asin(sqrt(th23)) );
620  calc ->SetDmsq32( hie*fabs(BASE_dmsq32) );
621 
622  std::vector <const IExperiment * > expts_temp;
623  for(int i = 0; i < (int) predsvector.size(); ++i){
624 
625  double pot = potFD_fhc;
626  if(horn=="both" && ( (!NuMuOnly && (i==1 || i>5) ) || // Nue+NuMu then vector is (nueFHC, nueRHC, numuFHC*4, numuRHC*4)
627  ( NuMuOnly && i>3 ) // NuMu only then vector is (numuFHC*4, numuRHC*4)
628  ) ) {
629  pot = potFD_rhc;
630  }
631  auto thisdata = GetFakeData (predsvector[i],calc, pot, cosmicsvector[i].first);
632  expts_temp.push_back(new SingleSampleExperiment(predsvector[i], *thisdata, *cosmicsvector[i].first, cosmicsvector[i].second, true));
633  }
634  expts_temp.push_back(WorldReactorConstraint2019());
635  auto exptsAll_temp = new MultiExperiment(expts_temp);
636  auxShifts = SystShifts::Nominal();
637 
638  std::cout << "\tStep " << i << " of " << steps << ", th23 = " << th23 << ", delta = " << delta << " ("<<delta/M_PI<<"*M_PI)" << std::endl;
639 
640  if(isSyst){
641  std::cout << "\t Getting syst list and the correlations." << std::endl;
642  if (NuMuOnly) {
643  std::cout << "\t\t NuMuOnly, so get a reduced list of systs..." << std::endl;
644  systs = GetJointFitSystematicList(isSyst, false, true, true, true, true, isPtExtrap);
645  auto notfornumuFHC = GetCorrelations(false, true , isPtExtrap);
646  auto notfornumuRHC = GetCorrelations(false, false, isPtExtrap);
647  for (size_t i=0; i<expts.size(); ++i) {
648  if (i<4) {
649  exptsAll_temp->SetSystCorrelations(i, notfornumuFHC);
650  } else {
651  if (i==8) continue;
652  exptsAll_temp->SetSystCorrelations(i, notfornumuRHC);
653  }
654  }
655  } else {
656  systs = GetJointFitSystematicList(isSyst, false, false, true, true, true, isPtExtrap);
657  auto notfornueFHC = GetCorrelations(true , true , isPtExtrap);
658  auto notfornueRHC = GetCorrelations(true , false, isPtExtrap);
659  auto notfornumuFHC = GetCorrelations(false, true , isPtExtrap);
660  auto notfornumuRHC = GetCorrelations(false, false, isPtExtrap);
661  if(horn=="fhc"){
662  exptsAll_temp->SetSystCorrelations(0, notfornueFHC);
663  for(int i=1; i<5;i++) {
664  exptsAll_temp->SetSystCorrelations(i, notfornumuFHC);
665  }
666  }
667  if(horn=="both"){
668  exptsAll_temp->SetSystCorrelations(0, notfornueFHC);
669  exptsAll_temp->SetSystCorrelations(1, notfornueRHC);
670  for(int i=0; i < 4; ++i) exptsAll_temp->SetSystCorrelations(i+2, notfornumuFHC);
671  for(int i=0; i < 4; ++i) exptsAll_temp->SetSystCorrelations(i+6, notfornumuRHC);
672  }
673  //Registry::Print();
674  }
675  }
676 
677  // ----- Now to do the MH sensitivities
678  if( whatplot=="mh" ){
679  std::cout << "\t\t Now to make the MH sensitivities...First the Normal Hierarchy." << std::endl;
680 
681  std::vector <const IFitVar *> fitvars = {&kFitDeltaInPiUnits, &kFitSinSqTheta23, &kFitDmSq32, &kFitSinSq2Theta13};
682  MinuitFitter fit(exptsAll_temp, fitvars, systs);
683 
684  ResetOscCalcToDefault(calc);
685  calc->SetDmsq32(hie*fabs(BASE_dmsq32));
686  chi2r[i] = fit.Fit(calc, auxShifts,
687  SeedList({
688  { &kFitDeltaInPiUnits, delta_seeds },
689  { &kFitSinSqTheta23 , th23_seeds }
690  }), {}, IFitter::kQuiet
691  );
692  std::cout << "\t\t And now the Inverse Hierarchy." << std::endl;
693  ResetOscCalcToDefault(calc);
694  calc->SetDmsq32(-hie*fabs(BASE_dmsq32));
695  chi2w[i] = fit.Fit(calc, auxShifts,
696  SeedList({
697  { &kFitDeltaInPiUnits, delta_seeds },
698  { &kFitSinSqTheta23 , th23_seeds }
699  }), {}, IFitter::kQuiet);
700  }
701 
702  // --- Now to do the Delta CP sensitivities
703  if( whatplot=="dcp" ){
704  std::cout << "\t\t Now to make the dCP sensitivities...First do DeltaCP = delta = " << delta << std::endl;
705 
706  std::vector <const IFitVar *> fitvars = {&kFitSinSqTheta23, &kFitDmSq32, &kFitSinSq2Theta13};
707  MinuitFitter fit(exptsAll_temp, fitvars, systs);
708 
709  ResetOscCalcToDefault(calc);
710  calc->SetDmsq32(hie*fabs(BASE_dmsq32));
711  calc->SetdCP(delta);
712  chi2r[i] = fit.Fit(calc, auxShifts,
713  SeedList({
714  { &kFitSinSqTheta23, th23_seeds},
715  { &kFitDmSq32 , {fabs(calc->GetDmsq32()), -fabs(calc->GetDmsq32())} }
716  }), {}, IFitter::kQuiet);
717  std::cout << "\t\t And now do DeltaCP = 0" << std::endl;
718  ResetOscCalcToDefault(calc);
719  calc->SetDmsq32(hie*fabs(BASE_dmsq32));
720  calc->SetdCP(0);
721  double chi2w0 = fit.Fit(calc, auxShifts,
722  SeedList({
723  { &kFitSinSqTheta23, th23_seeds},
724  { &kFitDmSq32 , {fabs(calc->GetDmsq32()), -fabs(calc->GetDmsq32())} }
725  }), {}, IFitter::kQuiet);
726  std::cout << "\t\t And now do DeltaCP = PI" << std::endl;
727  ResetOscCalcToDefault(calc);
728  calc->SetdCP(M_PI);
729  calc->SetDmsq32(hie*fabs(BASE_dmsq32));
730  double chi2wpi = fit.Fit(calc, auxShifts,
731  SeedList({
732  { &kFitSinSqTheta23, th23_seeds},
733  { &kFitDmSq32 , {fabs(calc->GetDmsq32()), -fabs(calc->GetDmsq32())} }
734  }), {}, IFitter::kQuiet);
735  // And find the minimum value.
736  chi2w[i] = std::min(chi2w0, chi2wpi);
737  std::cout << "The smallest chi2 was " << chi2w[i] << std::endl;
738  }
739 
740  // --- Now to do the Maximal Mixing sensitivities
741  if( whatplot=="maxmix" ){
742  std::cout << "\t\t Now to make the Maximal Mixing sensitivities...First do theta23 = " << th23 << std::endl;
743 
744  std::vector <const IFitVar *> fitvars = {&kFitDeltaInPiUnits, &kFitDmSq32, &kFitSinSq2Theta13};
745  std::vector <const IFitVar *> fitvars23 = {&kFitDeltaInPiUnits, &kFitDmSq32, &kFitSinSq2Theta13, &kFitSinSqTheta23 };
746  MinuitFitter fit (exptsAll_temp, fitvars , systs);
747  MinuitFitter fit23(exptsAll_temp, fitvars23, systs);
748 
749  ResetOscCalcToDefault(calc);
750  calc->SetDmsq32(hie*fabs(BASE_dmsq32));
751  calc->SetTh23 (asin(sqrt(th23)));
752  chi2r[i] = fit23.Fit(calc, auxShifts,
753  SeedList({
754  { &kFitDeltaInPiUnits, delta_seeds },
755  { &kFitDmSq32 , {fabs(calc->GetDmsq32()), -fabs(calc->GetDmsq32())} },
756  { &kFitSinSqTheta23 , th23_seeds }
757  }), {}, IFitter::kQuiet);
758  std::cout << "\t\t And now do theta23 = 0.5" << std::endl;
759  ResetOscCalcToDefault(calc);
760  calc->SetDmsq32(hie*fabs(BASE_dmsq32));
761  calc->SetTh23(asin(sqrt(0.5)));
762  chi2w[i] = fit.Fit(calc, auxShifts,
763  SeedList({
764  { &kFitDeltaInPiUnits, delta_seeds},
765  { &kFitDmSq32 ,{fabs(calc->GetDmsq32()), -fabs(calc->GetDmsq32())} }
766  }), {}, IFitter::kQuiet);
767  }
768 
769  // --- Now to do the Octant sensitivities
770  if ( whatplot=="oct" ) {
771  std::cout << "\t\t Now to make the octant sensitivities...First use the correct octant for th23 " << th23 << std::endl;
772  // Primitive flip between octants: minimize LO with kFitSinSqTheta23LowerOctant,
773  // then use wrong octant and minimize there with kFitSinSqTheta23UpperOctant
774  id=2;
775  revid=1;
776 
777  std::vector <const IFitVar *> fitvars1 = {&kFitDeltaInPiUnits, &kFitDmSq32, &kFitSinSq2Theta13, th23seeds[id ].var};
778  std::vector <const IFitVar *> fitvars2 = {&kFitDeltaInPiUnits, &kFitDmSq32, &kFitSinSq2Theta13, th23seeds[revid].var};
779 
780  MinuitFitter fit1(exptsAll_temp, fitvars1, systs); // Correct oct
781  MinuitFitter fit2(exptsAll_temp, fitvars2, systs); // Wrong oct
782 
783  ResetOscCalcToDefault(calc);
784  calc->SetDmsq32(hie*fabs(BASE_dmsq32));
785  calc->SetTh23(asin(sqrt(th23)));
786  chi2r[i] = fit1.Fit(calc, auxShifts,
787  SeedList({
788  { &kFitDeltaInPiUnits, delta_seeds },
789  { &kFitDmSq32 ,{fabs(calc->GetDmsq32()), -fabs(calc->GetDmsq32())} },
790  { th23seeds[id].var , th23seeds[id].seeds }
791  }), {}, IFitter::kQuiet);
792  std::cout << "\t\t And now, use the wrong octant " << 1-th23 << std::endl;
793  ResetOscCalcToDefault(calc);
794  calc->SetDmsq32(hie*fabs(BASE_dmsq32));
795  calc->SetTh23(asin(sqrt(1-th23)));
796  chi2w[i] = fit2.Fit(calc, auxShifts,
797  SeedList({
798  { &kFitDeltaInPiUnits, delta_seeds },
799  { &kFitDmSq32 ,{fabs(calc->GetDmsq32()), -fabs(calc->GetDmsq32())} },
800  { th23seeds[revid].var , th23seeds[revid].seeds }
801  }), {}, IFitter::kQuiet);
802  }
803 
804  if (Varyth23) Fits[i] = th23;
805  else Fits[i] = delta/M_PI;
806  chi[i] = sqrt(fabs(chi2w[i]-chi2r[i]));
807 
808  //std::cout<<delt<<" chi2 "<<chi[i]<<" right "<<chi2r[i]<<" wrong assumption :"<<chi2w[i]<< std::endl;
809  }
810  TString name = th23seed.label+" "+(hie>0?"NH":"IH");
811  slice.push_back({new TGraph(steps, Fits, chi), name}); //0 is ih, 1 is nh
812 
813  } // end hie
814  outfile->cd();
815 
816  TVectorD ss(1);
817  ss[0] = slice.size();
818  ss.Write("numsl");
819 
820  for(int i=0; i<(int) slice.size(); i++){
821  slice[i].first->Write(slice[i].second);
822  }
823 
824  outfile->Close();
825  std::cout<<"Done with slices"<< std::endl;
826  }//end MakeFile
827 
828  // ----> If not Make File, or when that loop finishes.
829  std::cout << "You are going to read the file containing the sensitivities -- " << FileName << std::endl;
830  TFile * infile = new TFile (FileName.c_str(),"read");
831  auto mins =* (TVectorD*)infile->Get("overall_minchi");
832  double minchi23 = mins[0];
833  std::cout<<"Min chi: "<<minchi23<< std::endl;
834 
835  auto ss =* (TVectorD*)infile->Get("numsl");
836  double numsl = ss[0];
837  std::cout<<"Num slices: "<<numsl<< std::endl;
838 
839  // Clear the slice vectors.
840  slice.clear();
841  slice2019.clear();
842 
843  for (double th23 : {0.565}){ //2019 BF
844 
845  th23helper th23seed;
846  if(th23 == 0.565) {th23seed = th23seeds[0];}
847  if(th23 == 0.404) {th23seed = th23seeds[1];}
848  if(th23 == 0.623) {th23seed = th23seeds[2];}
849 
850  for(int hie : {-1,1}){
851  TString name = th23seed.label+" "+(hie>0?"NH":"IH");
852  TGraph* get_slice;
853  infile->GetObject(name, get_slice);
854  if(name.Contains("2019")) slice2019.push_back({get_slice,name});
855  else slice.push_back({get_slice,name});
856  }
857  }
858 
859  new TCanvas;
860  double max = 5.1, xmin = 0, xmax = 2;
861  if (Varyth23) {
862  xmin = 0.4;
863  xmax = 0.6;
864  if(whatplot=="oct" ) max = 4.1;
865  if(whatplot=="maxmix" && isFutureSens) max = 6.1;
866  } else {
867  if(whatplot=="dcp" ) max = 3.1;
868  if(whatplot=="maxmix") max = 1.5;
869  if(whatplot=="oct" && NuMuOnly ) max = 0.5;
870  if(whatplot=="mh" && isFutureSens) max = 6.1;
871  }
872  DrawSliceCanvas(whatplot, max, xmin, xmax, false, false, isFutureSens );
873  for(unsigned int i =0; i < slice2019.size(); i++){
874 
875  if(slice2019[i].second.Contains("IH")) slice2019[i].first->SetLineColor(kPrimColorIH);
876  if(slice2019[i].second.Contains("NH")) slice2019[i].first->SetLineColor(k3SigmaConfidenceColorNH);
877 
878  slice2019[i].first->SetLineWidth(3);
879  slice2019[i].first->Draw("same");
880  }
881 
882  auto leg2 = SliceLegend(slice2019, (whatplot!="dcp"&&whatplot!="mh"), 0.53, whatplot=="mh", false, isSyst, isFutureSens);
883  leg2->Draw();
884  gPad->SetGrid();
885  if (isFutureSens) {
886  latex(0.138, 0.94, "NOvA 2020 Best Fit" , 0, 12, 1.3/30.);
887  latex(0.9 , 0.94, "#nu 31.5#times10^{20} + #bar{#nu} 31.5#times10^{20} POT", 0, 32, 1.3/30.);
888  SimulationSide();
889  } else {
890  Simulation();
891  }
892  gPad->Print((TString)PlotDir+TagName+"-Sens-"+sVary+"-"+whatplot+".pdf");
893 
894  }//end MakeSens
895 }
Spectrum * GetFakeData(const IPrediction *pred, osc::IOscCalc *calc, const double pot, TH1D *cosmics=0)
::xsd::cxx::tree::id< char, ncname > id
Definition: Database.h:165
const Color_t kPrimColorIH
Definition: Style.h:64
void Simulation()
Definition: tools.h:16
void DrawSliceCanvas(TString slicename, const double ymax, const double xmin=0, const double xmax=2.)
const XML_Char * name
Definition: expat.h:151
void latex(double x, double y, std::string txt, double ang=0, int align=22, double size=0.042)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::map< std::string, double > xmax
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
Antineutrinos-only.
Definition: IPrediction.h:50
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
(&#39; appearance&#39;)
Definition: IPrediction.h:18
double th23
Definition: runWimpSim.h:98
double delta
Definition: runWimpSim.h:98
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
(&#39;beam &#39;)
Definition: IPrediction.h:15
T sqrt(T number)
Definition: d0nt_math.hpp:156
const FitSinSqTheta23UpperOctant kFitSinSqTheta23UpperOctant
Definition: FitVars.cxx:16
constexpr T pow(T x)
Definition: pow.h:75
TH2 * Gaussian68Percent2D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 2D in gaussian approximation.
Float_t ss
Definition: plot.C:24
OStream cerr
Definition: OStream.cxx:7
const double kFutureRHCPOT
Definition: Exposures.h:245
const Color_t k1SigmaConfidenceColorIH
Definition: Style.h:81
const Color_t k2SigmaConfidenceColorNH
Definition: Style.h:77
TH2 * Gaussian2Sigma2D(const FrequentistSurface &s)
Up-value surface for 2 sigma confidence in 2D in gaussian approximation.
const double kFutureFHCPOT
Definition: Exposures.h:244
std::vector< const IPrediction * > GetNumuPredictions2020(const int nq=4, std::string decomp="noPt", osc::IOscCalc *calc=DefaultOscCalc(), bool useSysts=true, std::string beam="fhc", bool isFakeData=false, bool GetFromUPS=false, bool minimizeMemory=true, bool NERSC=false)
static SystShifts Nominal()
Definition: SystShifts.h:34
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
const Color_t k1SigmaConfidenceColorNH
Definition: Style.h:75
virtual void SetDmsq32(const T &dmsq32)=0
#define M_PI
Definition: SbMath.h:34
std::pair< Spectrum *, double > GetNueCosmics2020(std::string beam, bool GetFromUPS=false, bool NERSC=false)
const char * label
Log-likelihood scan across two parameters.
string outfilename
knobs that need extra care
std::string sFHC
Definition: MakeCutFlow.C:35
Charged-current interactions.
Definition: IPrediction.h:39
void sensitivity2020(bool MakeFile=false, bool MakeSurf=false, std::string whatsurf="th23dcp", bool MakeSens=false, bool Varyth23=false, std::string whatplot="mh", bool isSyst=false, bool isPtExtrap=false, bool NuMuOnly=false, std::string horn="both", bool isNH=true, std::string option="2019BestFit", bool OnGrid=false, bool isFutureSens=false)
string infile
const Color_t k3SigmaConfidenceColorNH
Definition: Style.h:78
const double kFutureFHCLivetime
Definition: Exposures.h:247
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
#define pot
TH2 * Gaussian3Sigma2D(const FrequentistSurface &s)
Up-value surface for 3 sigma confidence in 2D in gaussian approximation.
const IPrediction * GetNuePrediction2020(std::string decomp, osc::IOscCalc *calc, bool corrSysts, std::string beam, bool isFakeData, bool GetFromUPS=false, bool minimizeMemory=true, bool NERSC=false)
const double kAna2020FHCLivetime
Definition: Exposures.h:237
std::string PlotDir
virtual T GetDmsq32() const
Definition: IOscCalc.h:91
TLegend * SliceLegend(std::vector< std::pair< TGraph *, TString > > &graphs, bool isDelta)
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
const double kAna2020FHCPOT
Definition: Exposures.h:233
const double kFutureRHCLivetime
Definition: Exposures.h:248
const double kAna2020RHCPOT
Definition: Exposures.h:235
Neutrinos-only.
Definition: IPrediction.h:49
void DrawContourCanvas(TString surfName, double minx=0, double maxx=2, double miny=0, double maxy=1)
OStream cout
Definition: OStream.cxx:6
(&#39; survival&#39;)
Definition: IPrediction.h:19
std::vector< std::pair< Spectrum *, double > > GetNumuCosmics2020(const int nq=4, std::string beam="fhc", bool GetFromUPS=false, bool NERSC=false)
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
virtual T GetTh13() const
Definition: IOscCalc.h:93
const Color_t k3SigmaConfidenceColorIH
Definition: Style.h:83
TH2 * ToTH2(double minchi=-1) const
Definition: ISurface.cxx:161
TLegend * ContourLegend(int hie, bool fillContour, bool fccorr, Int_t kFillColor1, Int_t kFillColor2, Int_t kFillColor3, Int_t kDarkColor, bool bestFit)
Combine multiple component experiments.
const Style_t k2SigmaConfidenceStyle
Definition: Style.h:71
void SimulationSide()
Definition: rootlogon.C:137
double livetime
Definition: saveFDMCHists.C:21
void SaveTo(TDirectory *dir, const std::string &name) const
const FitSinSqTheta23LowerOctant kFitSinSqTheta23LowerOctant
Definition: FitVars.cxx:17
const double kAna2020RHCLivetime
Definition: Exposures.h:238
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
virtual void SetTh23(const T &th23)=0
Neutral-current interactions.
Definition: IPrediction.h:40
Interface definition for fittable variables.
Definition: IFitVar.h:16
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
std::vector< std::pair< const ISyst *, const ISyst * > > GetCorrelations(bool isNue, bool isFHC, bool ptExtrap)
const ReactorExperiment * WorldReactorConstraint2019()
Reactor constraint from PDG 2019.
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
T min(const caf::Proxy< T > &a, T b)
All neutrinos, any flavor.
Definition: IPrediction.h:26
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
static std::unique_ptr< FrequentistSurface > LoadFrom(TDirectory *dir, const std::string &name)
std::vector< const ISyst * > GetJointFitSystematicList(bool corrSysts, bool nueExclusive=false, bool numuExclusive=false, bool isFHC=true, bool isRHC=true, bool intersection=true, bool ptExtrap=true)
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
const Color_t k2SigmaConfidenceColorIH
Definition: Style.h:82
FILE * outfile
Definition: dump_event.C:13
virtual void SetdCP(const T &dCP)=0
T asin(T number)
Definition: d0nt_math.hpp:60
Compare a single data spectrum to the MC + cosmics expectation.
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17
enum BeamMode string