Tutorial2019FitSlices.C
Go to the documentation of this file.
3 #include "CAFAna/Analysis/Fit.h"
5 #include "CAFAna/Analysis/NuePlotStyle.h"
6 #include "CAFAna/Analysis/Surface.h"
12 #include "CAFAna/FC/FCSurface.h"
13 #include "CAFAna/Prediction/PredictionSystJoint2018.h"
15 #include "CAFAna/Systs/NumuSysts.h"
16 #include "CAFAna/Vars/FitVars.h"
17 #include "OscLib/IOscCalc.h"
18 #include "Utilities/rootlogon.C"
19 
20 #include "CAFAna/nue/Ana2019/joint_fit_2019_loader_tools.h"
21 
22 #include "TCanvas.h"
23 #include "TBox.h"
24 #include "TColor.h"
25 #include "TGraph.h"
26 #include "TVectorD.h"
27 #include "TF1.h"
28 #include "TLegend.h"
29 #include "TText.h"
30 #include "TLatex.h"
31 #include "TPad.h"
32 #include "TLine.h"
33 #include "TMarker.h"
34 #include "TStyle.h"
35 #include "TSystem.h"
36 #include "TGaxis.h"
37 
38 #include <algorithm>
39 #include <vector>
40 #include <string>
41 
42 using namespace ana;
43 
44 std::vector <Spectrum * > GetMockData(TString ana_options);
45 
46 void POTtag(bool both, bool FHCOnly, bool RHCOnly);
47 
48 void MakeMaps ( std::vector <const IFitVar*> profvars,
49  std::map<const IFitVar*, TGraph*> &profMap,
50  std::vector <const ISyst* > systs,
51  std::map<const ISyst*, TGraph*> &systMap);
52 
53 void SaveMaps ( TDirectory * dir,
54  TString prefix,
55  std::vector <const IFitVar*> profvars,
56  std::vector <TString > profnames,
57  std::map<const IFitVar*, TGraph*> &profMap,
58  std::vector <const ISyst* > systs,
59  std::map<const ISyst*, TGraph*> &systMap);
60 
61 
62 void Tutorial2019FitSlices(bool createFile = false, /// Do you want to make fit and write it into the file or just plot existing slice
63  bool corrSysts = true, /// syst or stat only?
64  TString options="joint_realData_both", /// what kind of analysis? numuOnly||nueOnly||joint ana for RHCOnly||RHCOnly||both modes
65  bool th23slice = false, /// is it th23 slice?
66  bool octantSlice = true, /// dou you want x2 curves split by octant?
67  bool dmsqSlice = false, /// is it dmsq32 slice?
68  bool th13Slice = false) /// is it th13 slice?
69  /// in case of th23slice = false && dmsqSlice = false && th13Slice = false it will be deltaCP slice
70 
71 {
72 
73  bool nueOnly = options.Contains("nueOnly");
74  bool numuOnly = options.Contains("numuOnly");
75  bool joint = options.Contains("joint");
76  assert (nueOnly || numuOnly || joint);
77 
78  bool FHCOnly = options.Contains("FHCOnly");
79  bool RHCOnly = options.Contains("RHCOnly");
80  bool both = options.Contains("both");
81  assert (FHCOnly || RHCOnly || both);
82 
83  bool fake2019 = options.Contains("fake2019");
84  bool realData = options.Contains("realData");
85 
86  auto suffix = options;
87  if(corrSysts) suffix+="_systs";
88  else suffix+="_stats";
89 
90 /* TString tag;
91  if(FHCOnly) tag="FHCOnly/";
92  if(RHCOnly) tag="RHCOnly/";
93  if(both) tag = "RHC_and_FHC/";*/
94 
95  if(th23slice) suffix+="_th23";
96  if(dmsqSlice) suffix+="_dmsq";
97  if(th13Slice) suffix+="_th13";
98  if(!octantSlice)suffix+="_noOct";
99 
100  TString outdir = "";
101  TString outfilename (outdir +"hist_slices_2019_" + suffix);
102  TString outFCfilename (outdir + "slices_FCInput_2019_" + suffix); /// need for syst. debug profiles
103 
104  if(createFile){
105 
106  //////////////////////////////////////////////////
107  // Load Nue and Numu experiments
108  //////////////////////////////////////////////////
109  //need numu only for prestage seeds
110 
111  struct predictions {
112  const string name;
113  const IPrediction * pred;
114  std::pair <TH1D*, double> cos;
115  double pot;
116  };
117 
118  std::vector <predictions> preds;
119  std::vector <Spectrum*> data;
120  std::vector <const IExperiment * > expts;
121 
122  auto calc_fake = DefaultOscCalc();
123  if(fake2019) SetFakeCalc(calc_fake, 0.565, 2.48e-3, 1.99); /// use 2019 BF for making fake data
124 
125  if(!numuOnly) {
126  if(FHCOnly || both ) {
127  preds.push_back({"nue_fhc", GetNuePrediction2019("combo", DefaultOscCalc(), corrSysts, "fhc", false), GetNueCosmics2019("fhc"), GetPOT()}); /// prediction downloading functions are defined in the ..._tools.h header, GetPOT () returns the POT exposure value for rhc or fhc mode
128  }
129  if(RHCOnly || both ) {
130  preds.push_back({"nue_fhc", GetNuePrediction2019("prop", DefaultOscCalc(), corrSysts, "rhc", false), GetNueCosmics2019("rhc"), GetPOT(false)});
131  }
132  }
133 
134  int nnumu = 4;
135  if(!nueOnly) {
136  if(FHCOnly || both ) {
137  auto numu_preds = GetNumuPredictions2019(nnumu, corrSysts, "fhc");
138  auto numu_cosmics = GetNumuCosmics2019(nnumu, "fhc");
139  for (int i = 0; i< nnumu; i++ ){
140  preds.push_back({"numu_fhc_"+std::to_string(i+1), numu_preds[i], numu_cosmics[i], GetPOT()});
141  }
142  }
143  if(RHCOnly || both ) {
144  auto numu_preds = GetNumuPredictions2019(nnumu, corrSysts, "rhc");
145  auto numu_cosmics = GetNumuCosmics2019(nnumu, "rhc");
146  for (int i = 0; i< nnumu; i++ ){
147  preds.push_back({"numu_rhc_"+std::to_string(i+1), numu_preds[i], numu_cosmics[i], GetPOT(false)});
148  }
149  }
150  }
151 
152  cout<<"\n\n Predictions used for the fit:\n";
153 
154  if(realData){
155  if(!numuOnly){
156  if(FHCOnly || both ) {
157  auto temp = GetMockData("nue_fhc");
158  data.insert(data.end(),temp.begin(), temp.end());
159  }
160  if(RHCOnly || both ) {
161  auto temp = GetMockData("nue_rhc");
162  data.insert(data.end(),temp.begin(), temp.end());
163  }
164  }
165  if(!nueOnly){
166  if(FHCOnly || both ){
167  auto temp = GetMockData("numu_fhc");
168  data.insert(data.end(), temp.begin(), temp.end());}
169  if(RHCOnly || both ){
170  auto temp = GetMockData("numu_rhc");
171  data.insert(data.end(), temp.begin(), temp.end());}
172  }
173  }
174 
175 
176  for(int i = 0; i < int(preds.size()); ++i){
177  double POT = preds[i].pot;
178  auto thisdata = GetFakeData (preds[i].pred, calc_fake, POT, preds[i].cos.first);
179  if(realData) thisdata = data[i];
180  cout<<i<<" "<<preds[i].name<<" POT "<<POT<<" tot MC "<<preds[i].pred->Predict(calc_fake).Integral(POT)<<" cos "<<preds[i].cos.first->Integral()<<" cos er "<<preds[i].cos.second<<" analyze data "<<thisdata->Integral(POT)<<endl;
181  expts.push_back(new SingleSampleExperiment(preds[i].pred, *thisdata, preds[i].cos.first, preds[i].cos.second));
182 
183  }
184 
185  ////////////////////////////////////////////////////////////
186  // Add constraints, make experiments
187  ////////////////////////////////////////////////////////////
188  std::cout << "\nCreating multiexperiment\n" << std::endl;
189 
190  if (!th13Slice) {expts.push_back(WorldReactorConstraint2017()); std::cout << "Adding WorldReactorConstraint2017\n";}
191  else {std::cout << "No reactor constraint, that's th13 slice\n";}
192  if(nueOnly) {
193  std::cout << "Adding Dmsq32ConstraintPDG2017\n";
194  expts.push_back(&kDmsq32ConstraintPDG2017);
195  }
196  std::cout << "Creating Multiexperiment with a total of "
197  << expts.size() << " experiments\n\n" << std::endl;
198  auto exptThis = new MultiExperiment(expts);
199 
200 
201  ////////////////////////////////////////////////////////////
202  // Systematics
203  ////////////////////////////////////////////////////////////
204  std::cout << "Systematics for the fit:\n";
205 
206  /// Now, let's prepare systematics
207  std::vector< const ISyst * > systs = {};
208  if (corrSysts) systs = {&kAnaCalibrationSyst};
209 
210  ////////////////////////////////////////////////////////////
211  // Fit
212  ////////////////////////////////////////////////////////////
213  std::cout << "Starting the fit" << std::endl;
214 
216 
217  SystShifts auxShifts = SystShifts::Nominal();
218 
219  // In case some systs need different seeds
220  // If you run ana script for 2019 you will see seeds for some systematics for paricular plots
221  // Comment for the workshop, but leave here for the education purposes
222 
223  std::vector <SystShifts> seedShifts = {};
224  /*if(corrSysts && !th23slice && !dmsqSlice && !th13Slice && octantSlice && (both || FHCOnly)){ //Cherenkov seed for deltaCP and neutron as well :(
225  cout<<"Add Cherenkov seeds \n";
226  for (double systshift:{-0.5, +0.5}){ //for the FHC Only or FHC+RHC ana
227  SystShifts tempShifts (&kAnaCherenkovSyst, systshift);
228  seedShifts.push_back(tempShifts);
229  }
230  if(both){
231  cout<<"Add neutron systematic seeds \n"; //For the joint RHC+FHC ana only
232  for (double systshift:{-0.5, +0.5}){
233  SystShifts tempShifts (&kNeutronVisEScalePrimariesSyst2018, systshift);
234  seedShifts.push_back(tempShifts);
235  }
236  }
237  }
238  */
239 
240 
241 //////////////////////////////////////////////////////////////////////
242 ///////////////////////// Seed controller ////////////////////////////
243 //////////////////////////////////////////////////////////////////////
244 
245  struct th23helper{
246  std::vector<double> seeds;
247  const IFitVar * var;
248  TString label;
249  };
250 
251  std::vector <th23helper> th23seeds;
252 
253  if(octantSlice) {
254  th23seeds.push_back( { {0.499, 0.45}, kFitSinSqTheta23LowerOctant, "LO"});
255  th23seeds.push_back( { {0.501, 0.55}, kFitSinSqTheta23UpperOctant, "UO"});
256  }
257  else th23seeds.push_back({ {0.45, 0.5, 0.55}, &kFitSinSqTheta23, ""});
258 
259  std::vector<double> delta_seeds = {0, 0.5, 1., 1.5};
260  std::vector<double> th23_old_seeds = {0.45, 0.5, 0.55};
261  double seed_dmsq = 2.5e-3;
262 
263 /////////////////////////////////////////////////////////////////////////
264 ////////////////////////// Different best fit searches //////////////////
265 /////////////////////////////////////////////////////////////////////////
266 
267  // Find the best fit points
268  double minchi23 = 1E20;
269  for(int hie:{-1, 1}){
270  for (auto & thseed:th23seeds){
271 
272  std::cout << "\n\nFinding best fit " << (hie > 0 ? "NH " : "IH ")
273  << thseed.label << ", "
274  << "ssth23 seeds ";
275  for (auto const & s:thseed.seeds) std::cout << s << ", ";
276  std::cout << std::endl;
277 
278  std::vector <const IFitVar*> fitvars = {&kFitDeltaInPiUnits,
279  thseed.var,
280  &kFitDmSq32,
282 
283  Fitter fit23(exptThis, fitvars, systs);
284  auto thisminchi = fit23.Fit(calc, auxShifts,
285  {
286  { &kFitDmSq32, {hie*seed_dmsq} },
287  { thseed.var, thseed.seeds },
288  { &kFitDeltaInPiUnits, delta_seeds }
289  }, seedShifts,
290  Fitter::kVerbose);
291  minchi23= min(minchi23, thisminchi);
292 
293  ResetOscCalcToDefault(calc);
294  auxShifts.ResetToNominal();
295  }
296  }
297  std::cout << "\nFound overall minchi " << minchi23 << "\n\n";
298 
299 /////////////////////////////////////////////////////////////////////////////
300 /////////////////////// Now, slices! ////////////////////////////////////////
301 ////////////////////////////////////////////////////////////////////////////
302 
303  ///We'll save to this file, make pretty plots later
304  TFile * outfile = new TFile(outfilename+".root","recreate");
305  TFile * outFCfile = new TFile(outFCfilename+".root","recreate");
306  /// again, save BF into the file
307  outfile->cd();
308  TVectorD v(1);
309  v[0] = minchi23;
310  v.Write("overall_minchi");
311  // how many points do you want?
312  int steps = 60;// 80;
313 
314  ///Now create all the slices
315  for(int hie: {-1, +1}){
316 
317  ResetOscCalcToDefault(calc);
318  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
319 
320  std::cout << "Starting profile " << (hie>0? "NH ": "IH") << "\n\n";
321  ///Structure defining the slice
322  struct ProfDef{
323  const IFitVar * fitvar; // which var to fix in each point? Or what kind of slice it is?
324  double minX;
325  double maxX;
326  std::vector <const IFitVar * > profvars; /// which vars to profile?
327  std::vector <TString> profvarnames; ///names of these vars
328  std::map <const IFitVar *, std::vector <double> >profseeds; /// seeds for osc. params
329  TString shortname;
330  };
331 
332  /// will loop once if you selected w/o octant split and twice if you want LO and UO curves
333  for(auto const & thseed:th23seeds){
334  ProfDef profdef;
335  const IFitVar* kCurTheta23 = thseed.var;
336 
337  /// deltaCP slice
338  if(!th23slice && !dmsqSlice && !th13Slice){
339  profdef.fitvar = &kFitDeltaInPiUnits;
340  profdef.minX = 0;
341  profdef.maxX = 2;
342  profdef.profvars = {&kFitDmSq32, &kFitSinSq2Theta13, kCurTheta23}; ///takes th23 vars from th23seeds (defined in seed controller) - it can be th23 var constrained in LO or UO, or ordinary th23
343  profdef.profvarnames = {"DmSq32","SinSq2Theta13","SinSqTheta23"};
344  profdef.profseeds = {{kCurTheta23, thseed.seeds}, { &kFitDmSq32, {hie*seed_dmsq}} }; /// th23 seed from the seed controller
345  profdef.shortname="delta";
346  }
347 
348  ///th23 slice - no split into octants
349  if(th23slice){
350  profdef.fitvar = &kFitSinSqTheta23;
351  profdef.minX = (nueOnly ? 0:0.3);
352  profdef.maxX = (nueOnly ? 1:0.7);
353  profdef.profvars = {&kFitDmSq32, &kFitSinSq2Theta13, &kFitDeltaInPiUnits};
354  profdef.profvarnames = {"DmSq32","SinSq2Theta13","DeltaCPpi"};
355  profdef.profseeds = { {&kFitDeltaInPiUnits,delta_seeds}};
356  profdef.shortname="th23";
357  }
358  ///th13 slice
359  if(th13Slice){
360  profdef.fitvar = &kFitSinSq2Theta13;
361  profdef.minX = 0.01;
362  profdef.maxX = 0.2;
363  profdef.profvars = {&kFitDmSq32, kCurTheta23, &kFitDeltaInPiUnits};
364  profdef.profvarnames = {"DmSq32","SinSqTheta23","DeltaCPpi"};
365  profdef.profseeds = { {kCurTheta23, thseed.seeds}, {&kFitDeltaInPiUnits,delta_seeds}, { &kFitDmSq32, {hie*seed_dmsq}} };
366  profdef.shortname="th13";
367  }
368  ///dmsq32 slice
369  if(dmsqSlice){
370  profdef.fitvar = &kFitDmSq32;
371  profdef.minX = hie > 0 ? 2e-3 : -3e-3;
372  profdef.maxX = hie > 0 ? 3e-3 : -2e-3;
373  profdef.profvars = {kCurTheta23, &kFitSinSq2Theta13, &kFitDeltaInPiUnits};
374  profdef.profvarnames = {"SinSqTheta23","SinSq2Theta13","DeltaCPpi"};
375  profdef.profseeds = { {kCurTheta23, thseed.seeds},
376  {&kFitDeltaInPiUnits,delta_seeds}};
377  profdef.shortname="dmsq";
378  }
379  // these three lines are arguments for SqrtProfile, they are used for FC
380  std::map<const IFitVar*, TGraph*> profVarsMap;
381  std::map<const ISyst*, TGraph*> systMap;
382  MakeMaps (profdef.profvars, profVarsMap, systs, systMap);
383 
384 
385  std::cout << " Profile "<< thseed.label << "\n";
386  ///make a profile
387  auto slice = SqrtProfile(exptThis, calc,
388  profdef.fitvar, steps, profdef.minX, profdef.maxX,
389  minchi23,
390  profdef.profvars,
391  systs,
392  profdef.profseeds,
393  seedShifts,
394  profVarsMap, systMap);
395  //Save!
396  TString hieroctstr = (hie > 0 ? "NH" : "IH") + thseed.label;
397  //save to the file
398  outfile->cd();
399  slice->Write((TString ) "slice_" + profdef.shortname + "_" + hieroctstr);
400  SaveMaps (outFCfile, hieroctstr,
401  profdef.profvars, profdef.profvarnames, profVarsMap, systs, systMap);
402 
403  }//end seeds
404  }//end hie
405 
406  delete outfile;
407  std::cout << "\nProfiles saved to " << outfilename << std::endl;
408  return;
409  }//end createfile
410 
411 
412 ///Done with fitting!
413 /// Now style work
414 
415  //////////////////////////////////////////////////////////////////////
416  // Make plots
417  //////////////////////////////////////////////////////////////////////
418  else{
419  std::cout << "Making plots from " << outfilename+".root" << std::endl;
420  TFile * infile = new TFile (outfilename+".root","read");
421 
422  //read minichi
423  auto mins =* (TVectorD*)infile->Get("overall_minchi");
424  double minchi23 = mins[0];
425 
426  std::vector <TString> sliceNames = {};
427  if (!th23slice && ! dmsqSlice && !th13Slice) sliceNames = {"slice_delta_"};
428  if (th23slice) sliceNames = {"slice_th23_"};
429  if (dmsqSlice) sliceNames = {"slice_dmsq_"};
430  if (th13Slice) sliceNames = {"slice_th13_"};
431 
432  // it shouldn't play any role, but copy old style seeding here for sake of naming, everywhere is just .lable
433 
434  struct th23helper{
435  std::vector<double> seeds;
436  const IFitVar * var;
437  TString label;
438  };
439  std::vector <th23helper> th23seeds;
440  if(octantSlice) {
441  th23seeds.push_back( { {0.45}, kFitSinSqTheta23LowerOctant, "LO"});
442  th23seeds.push_back( { {0.55}, kFitSinSqTheta23UpperOctant, "UO"});
443  }
444  else th23seeds.push_back({ {0.45, 0.5, 0.55}, &kFitSinSqTheta23, ""});
445 
446  //end of old-style-seeding
447 
448  /// loop over all slices in this file (IH, NH, LO, UO)
449  for(TString sliceName:sliceNames){
450  std::vector < std::pair <TGraph*, TString> > graphs;
451 
452  if(sliceName.Contains("delta")) {
453  for (TString hie:{"NH","IH"}){
454  for(auto const & thseed:th23seeds){
455  auto h = (TGraph*) infile ->Get( sliceName + hie + thseed.label);
456  if(!h) {
457  std::cerr << "Problem " << sliceName + hie + thseed.label << "\n";
458  exit(1);
459  }
460  graphs.push_back({h, hie + " " + thseed.label});
461  }
462  }
463  DrawSliceCanvas(sliceName, 7, 0, 2);
464  }
465  else {
466  for (TString hie:{"NH","IH"}) {
467  for(auto const & thseed:th23seeds){
468  auto h = (TGraph*) infile ->Get( sliceName + hie + thseed.label);
469  if(!h) {std::cerr << "Problem " << sliceName + hie + thseed.label << "\n"; }
470  graphs.push_back({h,hie + " " + thseed.label});
471  }
472  }
473 
474  if(sliceName.Contains("th23")) {
475  DrawSliceCanvas(sliceName, 5,
476  (nueOnly ? 0 : 0.32), (nueOnly ? 1 : 0.72));
477  }
478 
479  if(sliceName.Contains("dmsq")) {
480  DrawSliceCanvas(sliceName, 5, 2.15, 2.75);
481  }
482  if(sliceName.Contains("th13")) {
483  DrawSliceCanvas(sliceName, 5, 0.05, 0.2);
484  }
485  }
486 
487  //loop over all slices and set colors
488  for (auto &gr:graphs) {
489  if(gr.second.Contains("IH")) gr.first->SetLineColor(kPrimColorIH);
490  if(gr.second.Contains("NH")) gr.first->SetLineColor(k3SigmaConfidenceColorNH);
491  if(gr.second.Contains("LO")) gr.first->SetLineStyle(7);
492  gr.first->SetLineWidth(3);
493  if(dmsqSlice){
494  for (int i =0; i < gr.first->GetN(); ++i){
495  gr.first->SetPoint(i,1000*fabs(gr.first->GetX()[i]),gr.first->GetY()[i]);
496  }
497  }
498  gr.first->Draw("l same");
499  }
500 
501  //POTtag(both, FHCOnly, RHCOnly);
502  PreliminarySide();
503  gPad->RedrawAxis();
504 
505  auto leg = SliceLegend(graphs, !dmsqSlice && !th23slice && !th13Slice, 0.7, dmsqSlice);
506  leg->Draw();
507 
508  // outdir = "test_delta_slice/";
509  //gPad->SetGrid();
510  for(TString ext: {".pdf",".root"}) {
511  gPad->Print(outdir + "slice_" + suffix +
512  ext);
513  }//end ext
514 
515  }//end slicename
516 
517 
518  bool debug = true;
519 
520  if(debug){
521 
522  TFile * inFCfile = new TFile(outFCfilename+".root","read");
523  std::cout << "\n\nOpen profiled " << outFCfilename+".root\n" << std::endl;
524 
525  std::vector <TString> strprof = {"DmSq32","SinSq2Theta13","SinSqTheta23"};
526  if(th23slice) strprof = {"DmSq32","SinSq2Theta13","DeltaCPpi"};
527  if(th13Slice) strprof = {"DmSq32","DeltaCPpi", "SinSqTheta23"};
528  if(dmsqSlice) strprof = {"SinSqTheta23","SinSq2Theta13","DeltaCPpi"};
529 
530  std::vector< const ISyst * > systs = {};
531  if (corrSysts){ // syst case
532  systs = {&kAnaCalibrationSyst}; // replace the only one syst
533  }
534 
535  for (const ISyst* s : systs){
536  strprof.push_back(s->ShortName());
537  }
538 
539  for(auto const & str:strprof){
540  new TCanvas;
541  double miny=-1.5, maxy=1.5;
542  if(str=="DeltaCPpi") {miny=0; maxy=2;}
543  if(str=="DmSq32") {miny=2.2; maxy=2.9;}
544  if(str=="SinSq2Theta13") {miny=0.080; maxy=0.086;}
545  if(str=="SinSqTheta23") {miny=0.3; maxy=0.8;}
546 
547  auto axes = new TH2F("",";#delta_{CP}/#pi;" + str, 100, 0, 2, 100, miny, maxy);
548  if(dmsqSlice) axes = new TH2F("",";|#Deltam^{2}_{32}| (10^{-3} eV^{2});" + str, 100, 2.2, 2.8, 100, miny, maxy);
549  if(th23slice) axes = new TH2F("",";sin^{2}#theta_{23};" + str, 100, 0.3, 0.7, 100, miny, maxy);
550  if(th13Slice) axes = new TH2F("",";sin^{2}#theta_{13};" + str, 100, 0.01, 0.2, 100, miny, maxy);
551  axes->Draw("hist");
552 
553  for (TString hie:{"NH","IH"}){
554  for(auto const & thseed:th23seeds){
555  auto gr = (TGraph*) inFCfile->Get(hie + thseed.label + "_" + str);
556  if(!gr) {
557  std::cerr << "Problem " << hie + thseed.label + "_" + str << "\n";
558  exit(1);
559  }
560  if(str=="DmSq32"){
561  for (int i = 0; i < (int) gr->GetN(); ++i){
562  gr->SetPoint(i, gr->GetX()[i], 1000* fabs(gr->GetY()[i]));
563  }
564  }
565  if(dmsqSlice){
566  for (int i = 0; i < (int) gr->GetN(); ++i){
567  gr->SetPoint(i, 1000*fabs(gr->GetX()[i]), gr->GetY()[i]);
568  }
569  }
570  gr->SetLineColor(hie=="NH"?k3SigmaConfidenceColorNH:kPrimColorIH);
571  if(thseed.label.Contains("LO")) gr->SetLineStyle(7);
572  gr->SetLineWidth(3);
573  gr->SetMarkerColor(gr->GetLineColor());
574  gr->Draw("l same");
575  }//end th23
576  }//end hie
577  gPad->SetGrid();
578  gPad->Print(outdir + "debug_slice_prof_" + suffix + "_"+ str + ".pdf");
579  }//end profile vars
580  }//end debug
581 
582  }//end plots
583 }//end all
584 
585 
586 void MakeMaps ( std::vector <const IFitVar*> profvars,
587  std::map<const IFitVar*, TGraph*> &profMap,
588  std::vector <const ISyst* > systs,
589  std::map<const ISyst*, TGraph*> &systMap)
590 {
591  for (const IFitVar * var:profvars)
592  profMap.insert(std::pair<const IFitVar*, TGraph*> (var, new TGraph()));
593  for (const ISyst* syst : systs)
594  systMap.insert(std::pair<const ISyst*, TGraph*> (syst, new TGraph()));
595 }
596 void SaveMaps ( TDirectory * dir,
597  TString prefix,
598  std::vector <const IFitVar*> profvars,
599  std::vector <TString > profnames,
600  std::map<const IFitVar*, TGraph*> &profMap,
601  std::vector <const ISyst* > systs,
602  std::map<const ISyst*, TGraph*> &systMap)
603 {
604  TDirectory *tmp = gDirectory;
605  dir->cd();
606  for (int i = 0; i < (int) profvars.size(); ++i){
607  profMap[profvars[i]]->Write((prefix + "_" + profnames[i]));
608  }
609  for (const ISyst* s : systs)
610  systMap[s]->Write((prefix + "_" + s->ShortName()));
611  tmp->cd();
612 }
613 
614 
615 std::vector <Spectrum * > GetMockData(TString ana_options){
616 
617  cout<<"\n\n ============= You're loading Mock data from the file ============="<<endl;
618  std::vector <Spectrum *> spec;
619  std::string Dir = "./";
620  bool RHCmode = ana_options.Contains("rhc");
621  bool FHCmode = ana_options.Contains("fhc");
622  bool NueData = ana_options.Contains("nue");
623  bool NumuData = ana_options.Contains("numu");
624  if (NueData){
625  if (FHCmode) spec.push_back(LoadFromFile<Spectrum>(Dir+"mock_data_for_yn_tutorial.root", "data_exp_0").release());
626  if (RHCmode) spec.push_back(LoadFromFile<Spectrum>(Dir+"mock_data_for_yn_tutorial.root", "data_exp_1").release());
627  }
628  if (NumuData){
629  if (FHCmode){
630  for (int i = 2; i < 6; ++i) {
631  auto temp = LoadFromFile<Spectrum>(Dir+"mock_data_for_yn_tutorial.root", "data_exp_"+std::to_string(i)).release();
632  spec.push_back(temp);
633  }
634  }
635  if (RHCmode){
636  for (int i = 6; i < 10; ++i) {
637  auto temp = LoadFromFile<Spectrum>(Dir+"mock_data_for_yn_tutorial.root", "data_exp_"+std::to_string(i)).release();
638  spec.push_back(temp);
639  }
640  }
641 
642  }
643  return spec;
644 
645 }
646 
647 
Spectrum * GetFakeData(const IPrediction *pred, osc::IOscCalc *calc, const double pot, TH1D *cosmics=0)
const Color_t kPrimColorIH
Definition: Style.h:64
void DrawSliceCanvas(TString slicename, const double ymax, const double xmin=0, const double xmax=2.)
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
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
std::vector< const IPrediction * > GetNumuPredictions2019(const int nq=4, bool useSysts=true, std::string beam="fhc", bool GetFromUPS=false, ENu2018ExtrapType numuExtrap=kNuMu, bool minimizeMemory=false, bool NERSC=false)
double maxy
const Dmsq32Constraint kDmsq32ConstraintPDG2017(2.45e-3, 0.05e-3, 2.52e-3, 0.05e-3)
const FitSinSqTheta23UpperOctant kFitSinSqTheta23UpperOctant
Definition: FitVars.cxx:16
void SaveMaps(TDirectory *out, std::map< std::string, IDecomp * > decomps_nominal, std::map< std::string, std::map< std::string, std::map< int, IDecomp * > > > decomps_shifted, std::map< std::string, PredictionNoExtrap * > predsNE_nominal, std::map< std::string, std::map< std::string, std::map< int, PredictionNoExtrap * > > > predsNE_shifted, std::map< std::string, PredictionSterile * > predsSt_nominal, std::map< std::string, std::map< std::string, std::map< int, PredictionSterile * > > > predsSt_shifted)
Save all of the objects in the input maps to the out directory/file.
Definition: PPFXHelper.h:1077
OStream cerr
Definition: OStream.cxx:7
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
Float_t tmp
Definition: plot.C:36
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
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
const IPrediction * GetNuePrediction2019(std::string decomp, osc::IOscCalc *calc, bool corrSysts, std::string beam, bool isFakeData, bool GetFromUPS=false, bool minimizeMemory=false, bool NERSC=false)
virtual void SetDmsq32(const T &dmsq32)=0
const char * label
const XML_Char const XML_Char * data
Definition: expat.h:268
const XML_Char * s
Definition: expat.h:262
string outfilename
knobs that need extra care
const DummyAnaSyst kAnaCalibrationSyst("Calibration","AbsCalib")
string infile
const Color_t k3SigmaConfidenceColorNH
Definition: Style.h:78
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
virtual T GetDmsq32() const
Definition: IOscCalc.h:91
TLegend * SliceLegend(std::vector< std::pair< TGraph *, TString > > &graphs, bool isDelta)
void SetFakeCalc(osc::IOscCalcAdjustable *calc, double ssth23=-5, double dmsq32=-5, double dCP_Pi=-5)
void Tutorial2019FitSlices(bool createFile=false, bool corrSysts=true, TString options="joint_realData_both", bool th23slice=false, bool octantSlice=true, bool dmsqSlice=false, bool th13Slice=false)
const XML_Char * prefix
Definition: expat.h:380
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
Combine multiple component experiments.
std::vector< std::pair< TH1D *, double > > GetNumuCosmics2019(const int nq=4, std::string beam="fhc", bool GetFromUPS=false, bool NERSC=false)
TDirectory * dir
Definition: macro.C:5
const FitSinSqTheta23LowerOctant kFitSinSqTheta23LowerOctant
Definition: FitVars.cxx:17
void ResetToNominal()
Definition: SystShifts.cxx:144
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
exit(0)
assert(nhit_max >=nhit_nbins)
Interface definition for fittable variables.
Definition: IFitVar.h:16
double GetPOT(bool isFHC=true)
void MakeMaps(std::vector< const IFitVar * > profvars, std::map< const IFitVar *, TGraph * > &profMap, std::vector< const ISyst * > systs, std::map< const ISyst *, TGraph * > &systMap)
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
std::pair< Spectrum *, double > cos
void PreliminarySide()
Definition: rootlogon.C:114
Float_t e
Definition: plot.C:35
#define for
Definition: msvc_pragmas.h:3
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
std::vector< Spectrum * > GetMockData(TString ana_options)
const std::string outdir
FILE * outfile
Definition: dump_event.C:13
void POTtag(bool both, bool FHCOnly, bool RHCOnly)
std::pair< TH1D *, double > GetNueCosmics2019(std::string beam, bool GetFromUPS=false, bool NERSC=false)
Compare a single data spectrum to the MC + cosmics expectation.
TGraph * SqrtProfile(const IExperiment *expt, osc::IOscCalcAdjustable *calc, const IFitVar *v, int nbinsx, double minx, double maxx, double minchi, std::vector< const IFitVar * > profVars, std::vector< const ISyst * > profSysts, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts, std::map< const IFitVar *, TGraph * > &profVarsMap, std::map< const ISyst *, TGraph * > &systsMap, MinuitFitter::FitOpts opts)
Forward to Profile but sqrt the result for a crude significance.
Definition: Fit.cxx:143
enum BeamMode string