demoFitSlices.C
Go to the documentation of this file.
3 #include "CAFAna/Fit/Fit.h"
11 #include "CAFAna/FC/FCSurface.h"
16 #include "CAFAna/Vars/FitVars.h"
17 #include "OscLib/IOscCalc.h"
18 #include "CAFAna/nue/Ana2018/joint_fit_2018_tools.h"
19 #include "Utilities/rootlogon.C"
20 
21 #include "TCanvas.h"
22 #include "TBox.h"
23 #include "TColor.h"
24 #include "TGraph.h"
25 #include "TVectorD.h"
26 #include "TF1.h"
27 #include "TLegend.h"
28 #include "TText.h"
29 #include "TLatex.h"
30 #include "TPad.h"
31 #include "TLine.h"
32 #include "TMarker.h"
33 #include "TStyle.h"
34 #include "TSystem.h"
35 #include "TGaxis.h"
36 
37 #include <algorithm>
38 #include <vector>
39 #include <string>
40 
41 using namespace ana;
42 
43 void POTtag(bool both, bool FHCOnly, bool RHCOnly)
44 {
45  TLatex* ltx;
46  auto ltx1 = new TLatex(.16, 0.93, TString::Format("NOvA FD"));
47  if (FHCOnly){ ltx = new TLatex(.42, 0.93, TString::Format("%g#times10^{20} POT equiv. of #nu data only", kAna2017FullDetEquivPOT/1E20).Data());}
48  if (RHCOnly){ ltx = new TLatex(.52, 0.93, TString::Format("%1.2g#times10^{20} POT of #bar{#nu} data only", kAna2018RHCPOT/1E20).Data());}
49  if(both) {ltx = new TLatex(.35, 0.93, TString::Format("%g#times10^{20} POT equiv #nu + %1.2g#times10^{20} POT #bar{#nu}", kAna2017FullDetEquivPOT/1E20, kAna2018RHCPOT/1E20).Data());}
50  ltx1->SetNDC();
51  ltx1->SetTextAlign(11);
52  ltx1->SetTextSize(0.8*kBlessedLabelFontSize);
53  ltx1->Draw();
54  ltx->SetNDC();
55  ltx->SetTextAlign(11);
56  ltx->SetTextSize(0.8*kBlessedLabelFontSize);
57  ltx->Draw("same");
58 }
59 
60 
61 void MakeMaps ( std::vector <const IFitVar*> profvars,
62  std::map<const IFitVar*, TGraph*> &profMap,
63  std::vector <const ISyst* > systs,
64  std::map<const ISyst*, TGraph*> &systMap)
65 {
66  for (const IFitVar * var:profvars)
67  profMap.insert(std::pair<const IFitVar*, TGraph*> (var, new TGraph()));
68  for (const ISyst* syst : systs)
69  systMap.insert(std::pair<const ISyst*, TGraph*> (syst, new TGraph()));
70 }
71 void SaveMaps ( TDirectory * dir,
72  TString prefix,
73  std::vector <const IFitVar*> profvars,
74  std::vector <TString > profnames,
75  std::map<const IFitVar*, TGraph*> &profMap,
76  std::vector <const ISyst* > systs,
77  std::map<const ISyst*, TGraph*> &systMap)
78 {
79  TDirectory *tmp = gDirectory;
80  dir->cd();
81  for (int i = 0; i < (int) profvars.size(); ++i){
82  profMap[profvars[i]]->Write((prefix + "_" + profnames[i]));
83  }
84  for (const ISyst* s : systs)
85  systMap[s]->Write((prefix + "_" + s->ShortName()));
86  tmp->cd();
87 }
88 
89 
90 void demoFitSlices(bool createFile = false, /// Do you want to make fit and write it into the file or just plot existing slice
91  bool corrSysts = true, /// syst or stat only?
92  TString options="joint_both", /// what kind of analysis? numuOnly||nueOnly||joint ana for RHCOnly||RHCOnly||both modes
93  bool th23slice = false, /// is it th23 slice?
94  bool octantSlice = true, /// dou you want x2 curves split by octant?
95  bool dmsqSlice = false, /// is it dmsq32 slice?
96  bool th13Slice = false) /// is it th13 slice?
97  /// in case of th23slice = false && dmsqSlice = false && th13Slice = false it will be deltaCP slice
98 
99 {
100  ///First part follows demoFitContour logic, so I won't comment the loading and BF searches much
101 
102  bool nueOnly = options.Contains("nueOnly");
103  bool numuOnly = options.Contains("numuOnly");
104  bool joint = options.Contains("joint");
105  assert (nueOnly || numuOnly || joint);
106 
107  bool FHCOnly = options.Contains("FHCOnly");
108  bool RHCOnly = options.Contains("RHCOnly");
109  bool both = options.Contains("both");
110  assert (FHCOnly || RHCOnly || both);
111 
112  auto suffix = options;
113  if(corrSysts) suffix+="_systs";
114  else suffix+="_stats";
115 
116 /* TString tag;
117  if(FHCOnly) tag="FHCOnly/";
118  if(RHCOnly) tag="RHCOnly/";
119  if(both) tag = "RHC_and_FHC/";*/
120 
121  if(th23slice) suffix+="_th23";
122  if(dmsqSlice) suffix+="_dmsq";
123  if(th13Slice) suffix+="_th13";
124  if(!octantSlice)suffix+="_noOct";
125 
126  TString outdir = "";
127  TString outfilename (outdir +"hist_slices_2017_" + suffix);
128  TString outFCfilename (outdir + "slices_FCInput_2018_" + suffix); /// need for syst. debug profiles
129  if(joint && corrSysts && !th23slice && ! dmsqSlice && !th13Slice)
130  outFCfilename = outdir + "slices_FCInput_2018_"+ suffix;
131 
132 
133  if(createFile){
134 
135  //////////////////////////////////////////////////
136  // Load Nue and Numu experiments
137  //////////////////////////////////////////////////
138  //need numu only for prestage seeds
139  std::vector <const IPrediction * > preds;
140  std::vector <const IPrediction * > preds_numu_only;
141  std::vector <std::pair <TH1D*, double > > cosmics;
142  std::vector <std::pair <TH1D*, double > > cosmics_numu_only;
143  std::vector <Spectrum * > data;
144  std::vector <Spectrum * > data_numu_only;
145  std::vector <const IExperiment * > expts;
146  std::vector <const IExperiment * > expts_numu_only;
147 
148  auto calc_fake = DefaultOscCalc();
149  SetFakeCalc(calc_fake, 0.58, 2.51e-3, 0.17);
150 
151  if(!numuOnly) {
152  if(FHCOnly || both ) {
153  preds.push_back(GetNuePrediction2018("combo", DefaultOscCalc(), corrSysts, "fhc", false));
154  cosmics.push_back(GetNueCosmics2018("fhc"));
155  }
156  if(RHCOnly || both ) {
157  preds.push_back(GetNuePrediction2018("prop", DefaultOscCalc(), corrSysts, "rhc", false));
158  cosmics.push_back(GetNueCosmics2018("rhc"));
159  }
160  }
161 
162  int nnumu = 4;
163  if(!nueOnly) {
164  if(FHCOnly || both ) {
165 
166  auto numu_preds = GetNumuPredictions2018(nnumu, corrSysts, "fhc");
167  preds.insert(preds.end(),numu_preds.begin(), numu_preds.end());
168  preds_numu_only.insert(preds_numu_only.end(),numu_preds.begin(), numu_preds.end());
169 
170  auto numu_cosmics = GetNumuCosmics2018(nnumu, "fhc");
171  cosmics.insert(cosmics.end(),numu_cosmics.begin(), numu_cosmics.end());
172  cosmics_numu_only.insert(cosmics_numu_only.end(),numu_cosmics.begin(), numu_cosmics.end());
173  }
174  if(RHCOnly || both ) {
175 
176  auto numu_preds = GetNumuPredictions2018(nnumu, corrSysts, "rhc");
177  preds.insert(preds.end(),numu_preds.begin(), numu_preds.end());
178  preds_numu_only.insert(preds_numu_only.end(),numu_preds.begin(), numu_preds.end());
179 
180  auto numu_cosmics = GetNumuCosmics2018(nnumu, "rhc");
181  cosmics.insert(cosmics.end(),numu_cosmics.begin(), numu_cosmics.end());
182  cosmics_numu_only.insert(cosmics_numu_only.end(),numu_cosmics.begin(), numu_cosmics.end());
183  }
184  }
185 
186  cout<<"\n\n Predictions used for the fit:\n";
187 
188  for(int i = 0; i < (int) preds.size(); ++i){
189  double POT;
190  if(FHCOnly) POT = kAna2018FHCPOT;
191  if(RHCOnly) POT = kAna2018RHCPOT;
192  if(both) {
193  if(nueOnly || both) {
194  if(i==0) POT = kAna2018FHCPOT;
195  if(i==1) POT = kAna2018RHCPOT;
196  }
197  if(numuOnly) {
198  if(i >= 0 && i < 4) POT = kAna2018FHCPOT;
199  if(i >= 4 && i < 8) POT = kAna2018RHCPOT;
200  }
201  if(both) {
202  if(i >= 2 && i < 6) POT = kAna2018FHCPOT;
203  if(i >= 6 && i < 10) POT = kAna2018RHCPOT;
204  }
205  }
206  auto thisdata = GetFakeData (preds[i],calc_fake, POT, cosmics[i].first);
207  cout<<i<<" POT "<<POT<<" tot MC "<<preds[i]->Predict(calc_fake).Integral(POT)<<" cos "<<cosmics[i].first->Integral()<<" cos er "<<cosmics[i].second<<" analyze data "<<thisdata->Integral(POT)<<endl;
208  expts.push_back(new SingleSampleExperiment(preds[i],*thisdata, cosmics[i].first,cosmics[i].second));
209  }
210 
211  cout<<"Make numu only experiment to get the seeds"<<endl;
212  //make numu only experiment for seeds:
213  for(int i = 0; i < (int) preds_numu_only.size(); ++i){
214  double POT;
215  if(FHCOnly) POT = kAna2018FHCPOT;
216  if(RHCOnly) POT = kAna2018RHCPOT;
217  if(both) {
218  if(i >= 0 && i < 4) POT = kAna2018FHCPOT;
219  if(i >= 4 && i < 8) POT = kAna2018RHCPOT;
220  }
221  auto thisdata = GetFakeData (preds_numu_only[i],calc_fake, POT, cosmics_numu_only[i].first);
222  cout<<i<<" POT "<<POT<<" tot MC "<<preds_numu_only[i]->Predict(calc_fake).Integral(POT)<<" cos "<<cosmics_numu_only[i].first->Integral()<<" cos er "<<cosmics_numu_only[i].second<<" analyze data "<<thisdata->Integral(POT)<<endl;
223  expts_numu_only.push_back(new SingleSampleExperiment(preds_numu_only[i],*thisdata, cosmics_numu_only[i].first,cosmics_numu_only[i].second));
224  }
225 
226  ////////////////////////////////////////////////////////////
227  // Add constraints, make experiments
228  ////////////////////////////////////////////////////////////
229  std::cout << "\nCreating multiexperiment\n" << std::endl;
230 
231  if (!th13Slice) {expts.push_back(WorldReactorConstraint2017()); std::cout << "Adding WorldReactorConstraint2017\n";}
232  else {std::cout << "No reactor constraint, that's th13 slice\n";}
233  if(nueOnly) {
234  std::cout << "Adding Dmsq32ConstraintPDG2017\n";
235  expts.push_back(&kDmsq32ConstraintPDG2017);
236  }
237  std::cout << "Creating Multiexperiment with a total of "
238  << expts.size() << " experiments\n\n" << std::endl;
239  auto exptThis = new MultiExperiment(expts);
240 
241  std::cout << "Creating Multiexperiment of numu only SimpleExp with a total of "
242  << expts_numu_only.size() << " experiments\n\n" << std::endl;
243  auto exptThis_numu_only = new MultiExperiment(expts_numu_only);
244 
245  ////////////////////////////////////////////////////////////
246  // Systematics
247  ////////////////////////////////////////////////////////////
248  std::cout << "Systematics for the fit:\n";
249 
250  std::vector< const ISyst * > systs = {};
251  std::vector< const ISyst * > systs_numu_only = {};
252  if (corrSysts){ // syst case
253  systs = {&kAnaCalibrationSyst}; //here will be your own syst, which you made previously
254  systs_numu_only = {&kAnaCalibrationSyst}; //here will be your own syst, which you made previously
255  }
256 
257  ////////////////////////////////////////////////////////////
258  // Fit
259  ////////////////////////////////////////////////////////////
260  std::cout << "Starting the fit" << std::endl;
261 
263 
264  SystShifts auxShifts = SystShifts::Nominal();
265 
266  // In case some systs need different seeds
267  // If you run script with 2018 syst list you'll need seeds for some systematics
268  // Comment for the workshop, but leave here for the education purposes
269 
270  std::vector <SystShifts> seedShifts = {};
271  /*if(corrSysts && !th23slice && !dmsqSlice && !th13Slice && octantSlice && (both || FHCOnly)){ //Cherenkov seed for deltaCP and neutron as well :(
272  cout<<"Add Cherenkov seeds \n";
273  for (double systshift:{-0.5, +0.5}){ //for the FHC Only or FHC+RHC ana
274  SystShifts tempShifts (&kAnaCherenkovSyst, systshift);
275  seedShifts.push_back(tempShifts);
276  }
277  if(both){
278  cout<<"Add neutron systematic seeds \n"; //For the joint RHC+FHC ana only
279  for (double systshift:{-0.5, +0.5}){
280  SystShifts tempShifts (&kNeutronVisEScalePrimariesSyst2018, systshift);
281  seedShifts.push_back(tempShifts);
282  }
283  }
284  }
285  */
286 
287 
288 //////////////////////////////////////////////////////////////////////
289 ///////////////////////// Seed controller ////////////////////////////
290 //////////////////////////////////////////////////////////////////////
291 
292  cout<<"------------------- Start prestage seeds --------------------------"<<endl;
293 
294  double minchi_numu = 1E20;
295  double pre_seed_th23;
296  double pre_seed_dmsq;
297  ResetOscCalcToDefault(calc);
298  auxShifts.ResetToNominal();
299 
300  double maxmix = 0.514; // from the numu results
301  double numu_pre_seedLONH, numu_pre_seedUONH, numu_pre_seedLOIH, numu_pre_seedUOIH, dmsq_numu_pre_seedNH, dmsq_numu_pre_seedIH;
302 
303  for(int hie:{1}){
304  std::cout << "\n\nFinding test best fit " << (hie>0? "NH " : "IH ") << "\n\n";
305  Fitter fitnumu_only(exptThis_numu_only, {&kFitSinSqTheta23, &kFitDmSq32}, {});
306 
307  auto thisminchi = fitnumu_only.Fit(calc, auxShifts ,
308  {{&kFitDmSq32,{hie*2.5e-3}},
309  {&kFitSinSqTheta23, {0.4} }}, seedShifts,
310  Fitter::kVerbose);
311 
312  if(thisminchi < minchi_numu) minchi_numu = thisminchi;
313 
314  pre_seed_th23 = kFitSinSqTheta23.GetValue(calc);
315  pre_seed_dmsq = kFitDmSq32.GetValue(calc);
316  }
317 
318  numu_pre_seedLONH = ((pre_seed_th23>maxmix)?(2*maxmix-pre_seed_th23):pre_seed_th23);
319  numu_pre_seedUONH = ((pre_seed_th23>maxmix)?pre_seed_th23:(2*maxmix-pre_seed_th23));
320 
321  ResetOscCalcToDefault(calc);
322  auxShifts.ResetToNominal();
323 
324  cout<<"------------------- End prestage seeds --------------------------"<<endl;
325 
326 
327  struct th23helper{
328  std::vector<double> seeds;
329  const IFitVar * var;
330  TString label;
331  };
332 
333 
334  std::vector <th23helper> th23seeds;
335 
336  if(octantSlice) {
337  th23seeds.push_back( { {0.499, numu_pre_seedLONH}, kFitSinSqTheta23LowerOctant, "LO"});
338  th23seeds.push_back( { {0.501, numu_pre_seedUONH}, kFitSinSqTheta23UpperOctant, "UO"});
339  }
340  else th23seeds.push_back({ {numu_pre_seedLONH, 0.5, numu_pre_seedUONH}, &kFitSinSqTheta23, ""});
341 
342  std::vector<double> delta_seeds = {0, 0.5, 1., 1.5};
343  std::vector<double> th23_old_seeds = {0.45, 0.5, 0.55};
344 
345 /////////////////////////////////////////////////////////////////////////
346 ////////////////////////// Different best fit searches //////////////////
347 /////////////////////////////////////////////////////////////////////////
348 
349  // Find the best fit points
350  double minchi23 = 1E20;
351  for(int hie:{-1, 1}){
352  for (auto & thseed:th23seeds){
353 
354  std::cout << "\n\nFinding best fit " << (hie > 0 ? "NH " : "IH ")
355  << thseed.label << ", "
356  << "ssth23 seeds ";
357  for (auto const & s:thseed.seeds) std::cout << s << ", ";
358  std::cout << std::endl;
359 
360  std::vector <const IFitVar*> fitvars = {&kFitDeltaInPiUnits,
361  thseed.var,
362  &kFitDmSq32,
364 
365  Fitter fit23(exptThis, fitvars, systs);
366  cout<<" pre dmsq seed is "<<pre_seed_dmsq<<endl;
367  auto thisminchi = fit23.Fit(calc, auxShifts,
368  {
369  { &kFitDmSq32, {hie*pre_seed_dmsq} },
370  { thseed.var, thseed.seeds },
371  { &kFitDeltaInPiUnits, delta_seeds }
372  }, seedShifts,
373  Fitter::kVerbose);
374  minchi23= min(minchi23, thisminchi);
375 
376  if(corrSysts){
377  // Plot the systematic pulls and label with BFV - "dot" plot for systematics
378  auto shifts = PlotSystShifts(auxShifts);
379  TString str = "Best fit " ;
380  for (auto &v:fitvars){
381  str += TString::Format(" %s=%.3f ",v->LatexName().c_str(),v->GetValue(calc));
382  }
383  str+= TString::Format(" LL=%.6f", thisminchi);
384  shifts->SetTitle(str);
385  gPad->Update();
386  TLine *l=new TLine(gPad->GetUxmin(),0,gPad->GetUxmax(),0);
387  l->Draw("same");
388  gPad->Print(outdir + "debug_slice_shifts_bestfit_" + suffix +
389  (hie>0? "NH": "IH") + thseed.label + ".pdf");
390  }
391 
392  ResetOscCalcToDefault(calc);
393  auxShifts.ResetToNominal();
394  }
395  }
396  std::cout << "\nFound overall minchi " << minchi23 << "\n\n";
397 /////////////////////////////////////////////////////////////////////////////
398 /////////////////////// Now, slices! ////////////////////////////////////////
399 ////////////////////////////////////////////////////////////////////////////
400 
401  ///We'll save to this file, make pretty plots later
402  TFile * outfile = new TFile(outfilename+".root","recreate");
403  TFile * outFCfile = new TFile(outFCfilename+".root","recreate");
404  /// again, save BF into the file
405  outfile->cd();
406  TVectorD v(1);
407  v[0] = minchi23;
408  v.Write("overall_minchi");
409  // how many points do you want?
410  int steps = 60;// 80;
411 
412  ///Now create all the slices
413  for(int hie: {-1, +1}){
414 
415  ResetOscCalcToDefault(calc);
416  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
417 
418  std::cout << "Starting profile " << (hie>0? "NH ": "IH") << "\n\n";
419  ///Structure defining the slice
420  struct ProfDef{
421  const IFitVar * fitvar; // which var to fix in each point?
422  double minX;
423  double maxX;
424  std::vector <const IFitVar * > profvars; /// which vars to profile?
425  std::vector <TString> profvarnames; ///names of these vars
426  std::map <const IFitVar *, std::vector <double> >profseeds; /// seeds for osc. params
427  TString shortname;
428  };
429 
430  /// will loop once if you selected w/o octant split and twice if you want LO and UO curves
431  for(auto const & thseed:th23seeds){
432  ProfDef profdef;
433  const IFitVar* kCurTheta23 = thseed.var;
434 
435  /// deltaCP slice
436  if(!th23slice && !dmsqSlice && !th13Slice){
437  profdef.fitvar = &kFitDeltaInPiUnits;
438  profdef.minX = 0;
439  profdef.maxX = 2;
440  profdef.profvars = {&kFitDmSq32, &kFitSinSq2Theta13, kCurTheta23}; ///takes th23 vars from th23seeds (defined in seed controller) - it can be th23 varconstrained in LO or UO, or ordinary th23
441  profdef.profvarnames = {"DmSq32","SinSq2Theta13","SinSqTheta23"};
442  profdef.profseeds = {{kCurTheta23, thseed.seeds}, { &kFitDmSq32, {hie*pre_seed_dmsq}} }; /// th23 seed from the seed controller
443  profdef.shortname="delta";
444  }
445 
446  ///th23 slice - no split into octants
447  if(th23slice){
448  profdef.fitvar = &kFitSinSqTheta23;
449  profdef.minX = (nueOnly ? 0:0.3);
450  profdef.maxX = (nueOnly ? 1:0.7);
451  profdef.profvars = {&kFitDmSq32, &kFitSinSq2Theta13, &kFitDeltaInPiUnits};
452  profdef.profvarnames = {"DmSq32","SinSq2Theta13","DeltaCPpi"};
453  profdef.profseeds = { {&kFitDeltaInPiUnits,delta_seeds}};
454  profdef.shortname="th23";
455  }
456  ///th13 slice
457  if(th13Slice){
458  profdef.fitvar = &kFitSinSq2Theta13;
459  profdef.minX = 0.01;
460  profdef.maxX = 0.2;
461  profdef.profvars = {&kFitDmSq32, kCurTheta23, &kFitDeltaInPiUnits};
462  profdef.profvarnames = {"DmSq32","SinSqTheta23","DeltaCPpi"};
463  profdef.profseeds = { {kCurTheta23, thseed.seeds}, {&kFitDeltaInPiUnits,delta_seeds}, { &kFitDmSq32, {hie*pre_seed_dmsq}} };
464  profdef.shortname="th13";
465  }
466  ///dmsq32 slice
467  if(dmsqSlice){
468  profdef.fitvar = &kFitDmSq32;
469  profdef.minX = hie > 0 ? 2e-3 : -3e-3;
470  profdef.maxX = hie > 0 ? 3e-3 : -2e-3;
471  profdef.profvars = {kCurTheta23, &kFitSinSq2Theta13, &kFitDeltaInPiUnits};
472  profdef.profvarnames = {"SinSqTheta23","SinSq2Theta13","DeltaCPpi"};
473  profdef.profseeds = { {kCurTheta23, thseed.seeds},
474  {&kFitDeltaInPiUnits,delta_seeds}};
475  profdef.shortname="dmsq";
476  }
477  // these three lines are arguments for SqrtProfile, they are used for FC
478  std::map<const IFitVar*, TGraph*> profVarsMap;
479  std::map<const ISyst*, TGraph*> systMap;
480  MakeMaps (profdef.profvars, profVarsMap, systs, systMap);
481 
482 
483  std::cout << " Profile "<< thseed.label << "\n";
484  ///make a profile
485  auto slice = SqrtProfile(exptThis, calc,
486  profdef.fitvar, steps, profdef.minX, profdef.maxX,
487  minchi23,
488  profdef.profvars,
489  systs,
490  profdef.profseeds,
491  seedShifts,
492  profVarsMap, systMap);
493  //Save!
494  TString hieroctstr = (hie > 0 ? "NH" : "IH") + thseed.label;
495  //save to the file
496  outfile->cd();
497  slice->Write((TString ) "slice_" + profdef.shortname + "_" + hieroctstr);
498  SaveMaps (outFCfile, hieroctstr,
499  profdef.profvars, profdef.profvarnames, profVarsMap, systs, systMap);
500 
501  }//end seeds
502  }//end hie
503 
504  delete outfile;
505  std::cout << "\nProfiles saved to " << outfilename << std::endl;
506  return;
507  }//end createfile
508 
509 
510 ///Done with fitting!
511 /// Now stile work
512 
513  //////////////////////////////////////////////////////////////////////
514  // Make plots
515  //////////////////////////////////////////////////////////////////////
516  else{
517  std::cout << "Making plots from " << outfilename+".root" << std::endl;
518  TFile * infile = new TFile (outfilename+".root","read");
519 
520  //read minichi
521  auto mins =* (TVectorD*)infile->Get("overall_minchi");
522  double minchi23 = mins[0];
523 
524  std::vector <TString> sliceNames = {};
525  if (!th23slice && ! dmsqSlice && !th13Slice) sliceNames = {"slice_delta_"};
526  if (th23slice) sliceNames = {"slice_th23_"};
527  if (dmsqSlice) sliceNames = {"slice_dmsq_"};
528  if (th13Slice) sliceNames = {"slice_th13_"};
529 
530  // it shouldn't play any role, but copy old style seeding here for sake of naming, everywhere is just .lable
531 
532  struct th23helper{
533  std::vector<double> seeds;
534  const IFitVar * var;
535  TString label;
536  };
537  std::vector <th23helper> th23seeds;
538  if(octantSlice) {
539  th23seeds.push_back( { {0.45}, kFitSinSqTheta23LowerOctant, "LO"});
540  th23seeds.push_back( { {0.55}, kFitSinSqTheta23UpperOctant, "UO"});
541  }
542  else th23seeds.push_back({ {0.45, 0.5, 0.55}, &kFitSinSqTheta23, ""});
543 
544  //end of old-style-seeding
545 
546  /// loop over all slices in this ploe (IH, NH, LO, UO)
547  for(TString sliceName:sliceNames){
548  std::vector < std::pair <TGraph*, TString> > graphs;
549 
550  if(sliceName.Contains("delta")) {
551  for (TString hie:{"NH","IH"}){
552  for(auto const & thseed:th23seeds){
553  auto h = (TH1D*) infile ->Get( sliceName + hie + thseed.label);
554  if(!h) {
555  std::cerr << "Problem " << sliceName + hie + thseed.label << "\n";
556  exit(1);
557  }
558  graphs.push_back({new TGraph(h), hie + " " + thseed.label});
559  // graphs.push_back({HistoToCyclicGraph(h),hie + " " + thseed.label});
560  }
561  }
562  DrawSliceCanvas(sliceName, 5, 0, 2);
563  }
564  else {
565  for (TString hie:{"NH","IH"}) {
566  for(auto const & thseed:th23seeds){
567  auto h = (TH1D*) infile ->Get( sliceName + hie + thseed.label);
568  if(!h) {std::cerr << "Problem " << sliceName + hie + thseed.label << "\n"; }
569  graphs.push_back({new TGraph(h),hie + " " + thseed.label});
570  }
571  }
572 
573  if(sliceName.Contains("th23")) {
574  DrawSliceCanvas(sliceName, 3,
575  (nueOnly ? 0 : 0.32), (nueOnly ? 1 : 0.72));
576  }
577 
578  if(sliceName.Contains("dmsq")) {
579  DrawSliceCanvas(sliceName, 3, 2.15, 2.75);
580  }
581  if(sliceName.Contains("th13")) {
582  DrawSliceCanvas(sliceName, 3, 0.05, 0.2);
583  }
584  }
585 
586  //loop over all slices and set colors
587  for (auto &gr:graphs) {
588  if(gr.second.Contains("IH")) gr.first->SetLineColor(kPrimColorIH);
589  if(gr.second.Contains("NH")) gr.first->SetLineColor(k3SigmaConfidenceColorNH);
590  if(gr.second.Contains("LO")) gr.first->SetLineStyle(7);
591  gr.first->SetLineWidth(3);
592  if(dmsqSlice){ //HACK
593  for (int i =0; i < gr.first->GetN(); ++i){
594  gr.first->SetPoint(i,1000*fabs(gr.first->GetX()[i]),gr.first->GetY()[i]);
595  }
596  }
597  gr.first->Draw("l same");
598  }
599 
600  POTtag(both, FHCOnly, RHCOnly);
601  PreliminarySide();
602  gPad->RedrawAxis();
603 
604  auto leg = SliceLegend(graphs, !dmsqSlice && !th23slice && !th13Slice, 0.7, dmsqSlice);
605  leg->Draw();
606 
607  // outdir = "test_delta_slice/";
608  //gPad->SetGrid();
609  for(TString ext: {".pdf",".root"}) {
610  gPad->Print(outdir + "slice_" + suffix +
611  ext);
612  }//end ext
613 
614  }//end slicename
615 
616 
617  bool debug = true;
618 
619  if(debug){
620 
621  TFile * inFCfile = new TFile(outFCfilename+".root","read");
622  std::cout << "\n\nOpen profiled " << outFCfilename+".root\n" << std::endl;
623 
624  std::vector <TString> strprof = {"DmSq32","SinSq2Theta13","SinSqTheta23"};
625  if(th23slice) strprof = {"DmSq32","SinSq2Theta13","DeltaCPpi"};
626  if(th13Slice) strprof = {"DmSq32","DeltaCPpi", "SinSqTheta23"};
627  if(dmsqSlice) strprof = {"SinSqTheta23","SinSq2Theta13","DeltaCPpi"};
628 
629  std::vector< const ISyst * > systs = {};
630  if (corrSysts){ // syst case
631  systs = {&kAnaCalibrationSyst}; // replace by your new syst
632  }
633 
634  for (const ISyst* s : systs){
635  strprof.push_back(s->ShortName());
636  }
637 
638  for(auto const & str:strprof){
639  new TCanvas;
640  double miny=-1.5, maxy=1.5;
641  if(str=="DeltaCPpi") {miny=0; maxy=2;}
642  if(str=="DmSq32") {miny=2.2; maxy=2.9;}
643  if(str=="SinSq2Theta13") {miny=0.080; maxy=0.086;}
644  if(str=="SinSqTheta23") {miny=0.3; maxy=0.8;}
645 
646  auto axes = new TH2F("",";#delta_{CP}/#pi;" + str, 100, 0, 2, 100, miny, maxy);
647  if(dmsqSlice) axes = new TH2F("",";|#Deltam^{2}_{32}| (10^{-3} eV^{2});" + str, 100, 2.2, 2.8, 100, miny, maxy);
648  if(th23slice) axes = new TH2F("",";sin^{2}#theta_{23};" + str, 100, 0.3, 0.7, 100, miny, maxy);
649  if(th13Slice) axes = new TH2F("",";sin^{2}#theta_{13};" + str, 100, 0.01, 0.2, 100, miny, maxy);
650  axes->Draw("hist");
651 
652  for (TString hie:{"NH","IH"}){
653  for(auto const & thseed:th23seeds){
654  auto gr = (TGraph*) inFCfile->Get(hie + thseed.label + "_" + str);
655  if(!gr) {
656  std::cerr << "Problem " << hie + thseed.label + "_" + str << "\n";
657  exit(1);
658  }
659  if(str=="DmSq32"){
660  for (int i = 0; i < (int) gr->GetN(); ++i){
661  gr->SetPoint(i, gr->GetX()[i], 1000* fabs(gr->GetY()[i]));
662  }
663  }
664  if(dmsqSlice){
665  for (int i = 0; i < (int) gr->GetN(); ++i){
666  gr->SetPoint(i, 1000*fabs(gr->GetX()[i]), gr->GetY()[i]);
667  }
668  }
669  gr->SetLineColor(hie=="NH"?k3SigmaConfidenceColorNH:kPrimColorIH);
670  if(thseed.label.Contains("LO")) gr->SetLineStyle(7);
671  gr->SetLineWidth(3);
672  gr->SetMarkerColor(gr->GetLineColor());
673  gr->Draw("l same");
674  }//end th23
675  }//end hie
676  gPad->SetGrid();
677  gPad->Print(outdir + "debug_slice_prof_" + suffix + "_"+ str + ".pdf");
678  }//end profile vars
679  }//end debug
680 
681  }//end plots
682 }//end all
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
TH1 * PlotSystShifts(const SystShifts &shifts, bool sortName)
Definition: Plots.cxx:1495
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
double GetValue(const osc::IOscCalcAdjustable *osc) const override
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
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
void MakeMaps(std::vector< const IFitVar * > profvars, std::map< const IFitVar *, TGraph * > &profMap, std::vector< const ISyst * > systs, std::map< const ISyst *, TGraph * > &systMap)
Definition: demoFitSlices.C:61
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
std::vector< const IPrediction * > GetNumuPredictions2018(const int nq=4, bool useSysts=true, std::string beam="fhc", bool GetFromUPS=false, ENu2018ExtrapType numuExtrap=kNuMu, bool minimizeMemory=false, 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
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
virtual void SetDmsq32(const T &dmsq32)=0
void PreliminarySide()
Definition: rootlogon.C:114
const double kAna2018RHCPOT
Definition: Exposures.h:208
const char * label
const XML_Char const XML_Char * data
Definition: expat.h:268
void demoFitSlices(bool createFile=false, bool corrSysts=true, TString options="joint_both", bool th23slice=false, bool octantSlice=true, bool dmsqSlice=false, bool th13Slice=false)
Definition: demoFitSlices.C:90
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
void POTtag(bool both, bool FHCOnly, bool RHCOnly)
Definition: demoFitSlices.C:43
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)
std::vector< std::pair< TH1D *, double > > GetNumuCosmics2018(const int nq=4, std::string beam="fhc", bool GetFromUPS=false, bool NERSC=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::pair< TH1D *, double > GetNueCosmics2018(std::string beam, bool GetFromUPS=false, bool NERSC=false)
const double kAna2017FullDetEquivPOT
Definition: Exposures.h:195
TDirectory * dir
Definition: macro.C:5
const FitSinSqTheta23LowerOctant kFitSinSqTheta23LowerOctant
Definition: FitVars.cxx:17
void ResetToNominal()
Definition: SystShifts.cxx:143
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
exit(0)
assert(nhit_max >=nhit_nbins)
Interface definition for fittable variables.
Definition: IFitVar.h:16
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
const Float_t kBlessedLabelFontSize
Definition: Style.h:90
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
const IPrediction * GetNuePrediction2018(std::string decomp, osc::IOscCalc *calc, bool corrSysts, std::string beam, bool isFakeData, bool GetFromUPS=false, bool minimizeMemory=false, bool NERSC=false)
const double kAna2018FHCPOT
Definition: Exposures.h:207
Float_t e
Definition: plot.C:35
#define for
Definition: msvc_pragmas.h:3
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
const std::string outdir
FILE * outfile
Definition: dump_event.C:13
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