Tutorial2019FitContours.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 "Utilities/rootlogon.C"
18 #include "OscLib/IOscCalc.h"
19 
20 #include "CAFAna/nue/Ana2019/joint_fit_2019_loader_tools.h"
21 
22 #include "TCanvas.h"
23 #include "TMarker.h"
24 #include "TBox.h"
25 #include "TLatex.h"
26 #include "TColor.h"
27 #include "TGraph.h"
28 #include "TVectorD.h"
29 #include "TF1.h"
30 #include "TLegend.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 std::vector <Spectrum * > GetMockData(TString ana_options);
44 
45 void DrawHieTag(int hie, bool th13=true, bool vert=false);
46 
47 void FillGraphs(std::vector<TGraph*> g1,std::vector<TGraph*> g2,std::vector<TGraph*> g3,
48  const Int_t kFillColor1, const Int_t kFillColor2, const Int_t kFillColor3,
49  const Int_t kDarkColor, const TString surftype);
50 
51 void DrawReactorConstraint(double mean, double err);
52 
53 // kCMYK analog
55  {
56  Double_t white[9] = { 1, 1, 1, 1, 1, 1, 1, 1, 1 };
57  Double_t Red[9] = { 61./255., 99./255., 136./255., 181./255., 213./255., 225./255., 198./255., 99./255., 61./255.};
58  Double_t Green[9] = { 149./255., 140./255., 96./255., 83./255., 132./255., 178./255., 190./255., 140./255., 149./255};
59  Double_t Blue[9] = { 214./255., 203./255., 168./255., 135./255., 110./255., 100./255., 111./255., 203./255., 214./255.};
60  Double_t Length[9] = { 0.0000, 0.1250, 0.2500, 0.3750, 0.5000, 0.6250, 0.7500, 0.8750, 1.0000};
61 
62  TColor::CreateGradientColorTable(9, Length, Red, Green, Blue, 255, 1.0);
63  gStyle->SetNumberContours(255);
64  }
65 
66 /// Start the Demo script (it's a bit simplified nue / Ana2018 / FitandFC / joint_fit_2018_contours.C):
67 
68 void Tutorial2019FitContours(bool createFile=false, /// Do you want to fit and write it into the file (true) or make plots (false)?
69  bool corrSysts = true, /// Stat only (false) or with systs (true)?
70  TString options="joint_realData_both_onlyIH", /// Write exactly what fit is it? joint/nueOnly/numuOnly + both/FHCOnly/RHCOnly + onlyIH/onlyNH
71  bool dmsqSurf = false, /// What kind of plot do you want, is it dmsq vs th23 surface?
72  bool th13Surf = false, /// What kind of plot do you want, is it th13 vs delta CP surface?
73  /// If dmsqSurf = false && th13Surf = false then it will be th23 vs delta CP surface
74  bool zoomIn = true, /// Since it's pretty universal code axes are very extanded for cases of nueOnly || numuOnly || RHCOnly || FHCOnly.
75  /// Thus sometimes you'll want to zoom axis (especially in case of FHC&&RHC)
76  bool fillContour = false) /// Do you want to fill it with color or stay with border lines only?
77 {
78  /// Some general info about further fit
79 
80  bool nueOnly = options.Contains("nueOnly");
81  bool numuOnly = options.Contains("numuOnly");
82  bool joint = options.Contains("joint");
83  assert (nueOnly || numuOnly || joint);
84 
85  bool FHCOnly = options.Contains("FHCOnly");
86  bool RHCOnly = options.Contains("RHCOnly");
87  bool both = options.Contains("both");
88  assert (FHCOnly || RHCOnly || both);
89 
90  bool fake2019 = options.Contains("fake2019");
91  bool realData = options.Contains("realData");
92 
93  bool onlyNH = options.Contains("onlyNH");
94  bool onlyIH = options.Contains("onlyIH");
95 
96  auto suffix = options;
97  if(dmsqSurf) suffix+= "_dmsq";
98  if(th13Surf) suffix+= "_th13";
99  if(corrSysts) suffix+="_systs";
100  else suffix+="_stats";
101 
102  ///How do you want to organize and sort your outputs?
103 
104  /*TString outdir="";
105 
106  if(th13Surf) outdir += "th13_delta/";
107  else if (dmsqSurf) outdir += "th23_dmsq/";
108  else outdir += "delta_th23/";
109 
110  TString outfilename (outdir + (corrSysts?"syst/":"stat/") + "hist_contours_2018_" + suffix);*/
111 
112  /// For tute purposes save all plots in current dir
113  TString outdir="";
114  TString outfilename (outdir + "hist_contours_2019_" + suffix);
115 
116  /// Enter the fit and file creation part
117 
118  if(createFile){
119 
120 //////////////////////////////////////////////////
121 // Load Nue and Numu experiments
122 //////////////////////////////////////////////////
123 
124  /// Create vectors for all predictions we need - MC predictions, cosmic predictions, data or fake data
125  /// Then we'll need vectors of experiments for making a multiexperiment
126  /// And we'll need separate set of all items for making numu-only preseeds (since we follow 2018 ana procedure)
127 
128  struct predictions {
129  const string name;
130  const IPrediction * pred;
131  std::pair <TH1D*, double> cos;
132  double pot;
133  };
134 
135  std::vector <predictions> preds;
136  std::vector <Spectrum*> data;
137  std::vector <const IExperiment * > expts;
138 
139  auto calc_fake = DefaultOscCalc();
140  if(fake2019) SetFakeCalc(calc_fake, 0.565, 2.48e-3, 1.99); /// use 2019 BF for making fake data
141 
142  if(!numuOnly) {
143  if(FHCOnly || both ) {
144  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
145  }
146  if(RHCOnly || both ) {
147  preds.push_back({"nue_fhc", GetNuePrediction2019("prop", DefaultOscCalc(), corrSysts, "rhc", false), GetNueCosmics2019("rhc"), GetPOT(false)});
148  }
149  }
150 
151  int nnumu = 4;
152  if(!nueOnly) {
153  if(FHCOnly || both ) {
154  auto numu_preds = GetNumuPredictions2019(nnumu, corrSysts, "fhc");
155  auto numu_cosmics = GetNumuCosmics2019(nnumu, "fhc");
156  for (int i = 0; i< nnumu; i++ ){
157  preds.push_back({"numu_fhc_"+std::to_string(i+1), numu_preds[i], numu_cosmics[i], GetPOT()});
158  }
159  }
160  if(RHCOnly || both ) {
161  auto numu_preds = GetNumuPredictions2019(nnumu, corrSysts, "rhc");
162  auto numu_cosmics = GetNumuCosmics2019(nnumu, "rhc");
163  for (int i = 0; i< nnumu; i++ ){
164  preds.push_back({"numu_rhc_"+std::to_string(i+1), numu_preds[i], numu_cosmics[i], GetPOT(false)});
165  }
166  }
167  }
168 
169  /// Now assign proper POT for each of these predictions following the same structure
170  /// nue FHC, nue RHC, 4 numu FHC and 4 numu RHC quartiles
171  /// end make fake data. Feed everything to the experiments
172 
173  cout<<"\n\n Predictions used for the fit:\n";
174 
175  if(realData){
176  if(!numuOnly){
177  if(FHCOnly || both ) {
178  auto temp = GetMockData("nue_fhc");
179  data.insert(data.end(),temp.begin(), temp.end());
180  }
181  if(RHCOnly || both ) {
182  auto temp = GetMockData("nue_rhc");
183  data.insert(data.end(),temp.begin(), temp.end());
184  }
185  }
186  if(!nueOnly){
187  if(FHCOnly || both ){
188  auto temp = GetMockData("numu_fhc");
189  data.insert(data.end(), temp.begin(), temp.end());}
190  if(RHCOnly || both ){
191  auto temp = GetMockData("numu_rhc");
192  data.insert(data.end(), temp.begin(), temp.end());}
193  }
194  }
195 
196  for(int i = 0; i < int(preds.size()); ++i){
197  double POT = preds[i].pot;
198  auto thisdata = GetFakeData (preds[i].pred, calc_fake, POT, preds[i].cos.first);
199  if(realData) thisdata = data[i];
200  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;
201  expts.push_back(new SingleSampleExperiment(preds[i].pred, *thisdata, preds[i].cos.first, preds[i].cos.second));
202 
203  }
204 
205 
206 //////////////////////////////////////////////////
207 /// Add constraints, make experiments
208 //////////////////////////////////////////////////
209 
210  /// Lets make a multiexperiment for simultaneous fit
211  /// For different contours/fits we will use different additional information
212 
213  std::cout << "\nAll experiments are done \nCreating multiexperiment\n" << std::endl;
214  if(!th13Surf){ /// We use reactor constraint as 11th "experiment" - th13 will be varied in PDG's error range with PDG central value
215  /// but we wouldn't use such constraint if we want to make our own measurement, right?
216  std::cout << "Adding WorldReactorConstraint2017\n";
217  expts.push_back(WorldReactorConstraint2017());
218  }
219  if(nueOnly) { /// If we don't use NOvA's result in the fit, then we will use PDG's constraint
220  std::cout << "Adding Dmsq32ConstraintPDG2017\n";
221  expts.push_back(&kDmsq32ConstraintPDG2017);
222  }
223  /// Let's make a multiexperiment
224  std::cout << "Creating Multiexperiment with a total of "
225  << expts.size() << " experiments\n\n" << std::endl;
226  auto exptThis = new MultiExperiment(expts);
227 
228  std::cout<<"Clean verctors"<<std::endl;
229  preds.clear();
230 
231 ////////////////////////////////////////////////////////////
232 /// Systematics
233 ////////////////////////////////////////////////////////////
234 
235  /// Now, let's prepare systematics
236  std::vector< const ISyst * > systs = {};
237  if (corrSysts) systs = {&kAnaCalibrationSyst};
238 
239 
240 ////////////////////////////////////////////////////////////
241 /// Fit
242 ////////////////////////////////////////////////////////////
243  /// And finally let's start fitting
244 
245  std::cout << "Starting the fit" << std::endl;
246 
248 
249  SystShifts auxShifts = SystShifts::Nominal();
250  // In case some systs need different seeds
251  std::vector <SystShifts> seedShifts = {};
252 
253 //////////////////////////////////////////////////////////////////////
254 /////////////////////////// Seed controller ////////////////////////////
255 ////////////////////////////////////////////////////////////////////////
256  /// This section contains all seeds that will be used further
257 
258  struct th23helper{
259  std::vector<double> seeds;
260  const IFitVar * var;
261  TString label;
262  };
263 
264  std::vector <th23helper> th23seeds;
265  //for NH:
266  /// make two seeds in both UO and LO (near maxmixing and th23 from pre-seeding stage)
267  th23seeds.push_back( { {0.499, 0.45}, kFitSinSqTheta23LowerOctant, "LO"});
268  th23seeds.push_back( { {0.501, 0.55}, kFitSinSqTheta23UpperOctant, "UO"});
269 
270  /// Put three seeds - LO, UO, maxmixing
271  std::vector<double> th23_all_seeds_NH = {0.45, 0.5, 0.55};
272  /// delta CP seeds
273  std::vector<double> delta_seeds = {0, 0.5, 1., 1.5};
274  double seed_dmsq = 2.5e-3;
275 
276 /////////////////////////////////////////////////////////////////////////
277 //////////////////////////// Different best fit searches //////////////////
278 ///////////////////////////////////////////////////////////////////////////
279  /// Now, let's find BF value, it searches BF in IH/NH and UO/LO and chooses the best chisq
280  double minchi23 = 1E20;
281  for(int hie:{-1,1}){
282  /// will use th23seeds which we defined previously for LO and UO
283  for (auto & thseed: th23seeds){
284  std::cout << "\n\nFinding best fit "
285  << thseed.label
286  << (hie>0? " NH " : " IH ")
287  << "\n\n";
288  for (auto const & s:thseed.seeds) std::cout << s << ", ";
289  std::cout << std::endl;
290  std::vector <const IFitVar*> fitvars = {&kFitDeltaInPiUnits,
291  thseed.var,
294 
295  Fitter fit23(exptThis, fitvars, systs);
296  auto thisminchi = fit23.Fit(calc,auxShifts ,{
297  {&kFitDmSq32Scaled,{hie*seed_dmsq}},
298  {thseed.var, thseed.seeds},
299  {&kFitDeltaInPiUnits, delta_seeds},
300  {&kFitSinSq2Theta13, {calc->GetTh13()}} },
301  {}, Fitter::kVerbose
302  );
303  /// compare with global BF
304  minchi23= min(minchi23, thisminchi);
305  ResetOscCalcToDefault(calc);
306  auxShifts.ResetToNominal();
307  }//end th23
308  }//end hie
309  std::cout << "\nFound overall minchi " << minchi23 << std::endl;
310 
311  //We'll save to this file, make pretty plots later
312  TFile * outfile = new TFile(outfilename+".root","recreate");
313 
314  outfile->cd();
315  TVectorD v(1);
316  v[0] = minchi23;
317  cout<<"minchi wrote into the file "<<minchi23<<endl;
318  v.Write("overall_minchi");
319 
320  /// How much point on the grid do you want? It will be steps*steps
321  int steps = 30;//10;// 5;
322  if(dmsqSurf){
323  std::cerr << "\n WARNING Using 20 bins for dmsq surface\n\n";
324  steps = 20;
325  }
326  //Now create all the surfaces
327  /// Generally need both IH and NH, but it's too long to run them all together
328  for(int hie: {-1, +1}){
329 
330  if((onlyNH && hie<0) || (onlyIH && hie>0)) continue; /// Will run only one hie which you chose at the beginning
331 
332  std::cout << "Starting surface " << (hie>0? "NH ": "IH") << "\n\n";
333  cout<<"the delta seeds are ";
334  for (auto const & s:delta_seeds) std::cout << s << ", ";
335  std::string hieStr = hie>0 ? "NH" : "IH";
336 
337  /// It will be deltaCP vs Th23 surface
338  if(!th13Surf && ! dmsqSurf){
339  ResetOscCalcToDefault(calc);
340  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
341  /// Set Y axis's low and up values
342  double low = (nueOnly||RHCOnly)?0:0.3;
343  double high = (nueOnly||RHCOnly)?1:0.7;
344  cout<<"low "<<low<<" high "<<high<<endl;
345  /// It will minimize dmsq23 and th13 only, no seeds here
346  Surface surf23(exptThis, calc,
347  &kFitDeltaInPiUnits,steps, 0, 2,
348  &kFitSinSqTheta23, steps, low, high,
350  //{},seedShifts
351  );
352 
353  auto surf1=surf23.ToTH2(minchi23);
354  outfile->cd();
355  surf1->Write((TString)"delta_th23_"+hieStr);
356  surf23.SaveTo(outfile->mkdir((TString)"surf_delta_th23_"+hieStr));
357  }
358 
359  /// This part will make th23 vs dmsq surface
360  if(!th13Surf && dmsqSurf){
361  calc = DefaultOscCalc();
362  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
363  /// Put Y axis's low and up values
364  double low = (nueOnly||RHCOnly)?0.2:0.35;
365  double high = (nueOnly||RHCOnly)?1:0.7;
366  cout<<"low "<<low<<" high "<<high<<endl;
367  /// It will minimize deltaCP and th13, use seeds for deltaCP
368  /// also, there is some logic for x axis's values as well
369  /// for different cases of fit (RHC/FHC/both) contours can be larger or smaller. To much pain in these values ...
370  Surface surf23m(exptThis, calc,
371  &kFitSinSqTheta23, steps,low, high,
372  &kFitDmSq32Scaled, steps,
373  (hie>0?(RHCOnly?2.0:2.2):(RHCOnly?-3.5:-2.8)),(hie>0?(RHCOnly?3.2:2.85):(RHCOnly?-2.0:-2.3)),
375  {{&kFitDeltaInPiUnits, delta_seeds}});
376  auto surf6=surf23m.ToTH2(minchi23);
377  outfile->cd();
378  surf6->Write((TString)"th23_dm32_"+(hie>0?"NH":"IH"));
379  surf23m.SaveTo(outfile->mkdir((TString)"surf_th23_dm32_"+(hie>0?"NH":"IH")));
380  }
381 
382  /// Also, we can do th13 vs deltaCP contour, why not
383  /// (we don't show them to public since there are some issues with FC for th13 measurements)
384  /// but gaussian fit is ok
385  if(th13Surf){
386  calc = DefaultOscCalc();
387  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
388  /// Minimize dmsq23, th23, put seeds for th23 as (smth in LO, 0.5, smth in UO)
389  Surface surf13(exptThis, calc,
390  &kFitSinSq2Theta13, steps, 0, (RHCOnly?1:0.35),
391  &kFitDeltaInPiUnits,steps, 0, 2,
393  {{&kFitSinSqTheta23, th23_all_seeds_NH}});//, seedShifts);
394  auto surf4 = surf13.ToTH2(minchi23);
395  outfile->cd();
396  surf4->Write((TString)"th13_delta_"+(hie>0?"NH":"IH"));
397  surf13.SaveTo(outfile->mkdir((TString)"surf_th13_delta_"+(hie>0?"NH":"IH")));
398 
399  }
400  } //end hie
401 
402  //delete outfile;
403  std::cout << "\nSurfaces saved to " << outfilename << std::endl;
404  return;
405  }//end createfile
406 
407 /// Done with fit! Now play with style and make pretty plots.
408 
409 //////////////////////////////////////////////////////////////////////
410 // Make plots
411 //////////////////////////////////////////////////////////////////////
412 
413  else{ //plot only
414 
415  /// read file from the previous section
416  TFile * infile = new TFile (outfilename+".root","read");
417  /// Get BF
418  auto mins =* (TVectorD*)infile->Get("overall_minchi");
419  double minchi23 = mins[0];
420  /// Some logic to select contour which you will plot
421  std::vector <TString> surfNames;
422  if(!th13Surf && ! dmsqSurf)surfNames = {"delta_th23_"};
423  if (dmsqSurf) surfNames.push_back("th23_dm32_");
424  if (th13Surf) surfNames.push_back("th13_delta_");
425  for(TString surfName:surfNames){
426  /// again, need two hierarchies, but all files were done with just one type
427  for (int hie:{-1,+1}){
428 
429  if((onlyNH && hie<0) || (onlyIH && hie>0)) continue;
430 
431  //Setup canvas for each contour type, zoom if necessary
432  if(surfName.Contains("delta_th23")) {
433  if (zoomIn) DrawContourCanvas(surfName, 0., 2., (RHCOnly?0.15:0.275), (RHCOnly?0.8:0.725));
434  else DrawContourCanvas(surfName, 0., 2., 0., 1.);
435  }
436  if(surfName.Contains("th23_dm32")) {
437  if (hie>0) DrawContourCanvas(surfName, ((nueOnly||RHCOnly)?0.2:0.3), ((nueOnly||RHCOnly)?0.8:0.7), (RHCOnly?1.8:2.05), (RHCOnly?3.3:2.85));
438  else DrawContourCanvas(surfName, ((nueOnly||RHCOnly)?0.2:0.3), ((nueOnly||RHCOnly)?0.8:0.7), (RHCOnly?-3.5:-2.95), (RHCOnly?-2.0:-2.25));
439  }
440  if(surfName.Contains("th13_delta")) {
441  if (zoomIn) DrawContourCanvas(surfName, 0., 0.38, 0, 2);
442  else DrawContourCanvas(surfName, 0., 1.0, 0, 2);
443  }
444  /// Read your surface and set 1,2,3 sigma curves
445  auto surf = *Surface::LoadFrom(infile->GetDirectory("surf_"+surfName+(hie>0?"NH":"IH")));
446  TH2 * hsurf = (TH2*)infile->Get(surfName+(hie>0?"NH":"IH"));
447  TH2 * surf_1Sigma = Gaussian68Percent2D(surf);
448  TH2 * surf_2Sigma = Gaussian2Sigma2D(surf);
449  TH2 * surf_3Sigma = Gaussian3Sigma2D(surf);
450  /// colors
451  Color_t k2SigmaConfidenceColor = ((hie>0)? k2SigmaConfidenceColorNH:k2SigmaConfidenceColorIH);
452  Int_t kFillColor1 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.75);
453  Int_t kFillColor2 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.38);
454  Int_t kFillColor3 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.18);
455  Int_t kDarkColor = kGray+3;
456 
457  if(fillContour){
458  auto g1 = surf.GetGraphs(surf_1Sigma,minchi23);
459  auto g2 = surf.GetGraphs(surf_2Sigma,minchi23);
460  auto g3 = surf.GetGraphs(surf_3Sigma,minchi23);
461  FillGraphs(g1, g2, g3,
462  kFillColor1, kFillColor2, kFillColor3,
463  kDarkColor, suffix + surfName + (hie>0?"NH":"IH"));
464  }
465  else{
466  cout<<"minchi from file is "<<minchi23<<endl;
467 
468  surf.DrawContour(surf_1Sigma, kSolid, kFillColor1,minchi23);
469  surf.DrawContour(surf_2Sigma, k2SigmaConfidenceStyle, kFillColor2,minchi23);
470  surf.DrawContour(surf_3Sigma, kSolid, kFillColor3,minchi23);
471  }
472 
473  gPad->RedrawAxis();
474 
475  /// Set legend
476  auto leg = ContourLegend(hie, fillContour, false,
477  kFillColor1,kFillColor2,kFillColor3,kDarkColor,
478  false);
479  /// Put PDG constraint line on the th13 vs deltaCP surface
480  if(th13Surf) {
481  DrawReactorConstraint(0.082,0.004);
482  leg->SetFillStyle(1001);
483  leg->SetNColumns(1);
484  leg->SetX1(0.7); leg->SetX2(0.8);
485  leg->SetY1(0.4); leg->SetY2(0.85);
486  leg->SetHeader("No FC");
487 
488  TGraph * dummy = new TGraph;
489  dummy->SetLineWidth(2);
490  dummy->SetLineColor(kGray+2);
491  dummy->SetFillColor(kGray);
492  leg->AddEntry(dummy,"#splitline{Reactor}{68%C.L.}","f");
493  }
494  leg->Draw();
495  Preliminary();
496  DrawHieTag(hie,false);
497 
498  /// and finally (!) save the plot to the file
499  for(TString ext: {".pdf",".root"}){
500  gPad->Print(outdir + "contour_" + suffix + "_" +
501  surfName +
502  (hie>0?"NH":"IH") +
503  (zoomIn?"_zoom":"") +
504  ext);
505  }//end ext
506 
507  /// in case if you want (of course you want) to see debug plots for each syst or osc parameter put true here
508  auto debug = true;
509  if(debug){
510  auto hists = surf.GetProfiledHists();
511  int nSysts = corrSysts ? 51 : 0;
512  for (int i = 0; i < nSysts; ++i){
513  hists.push_back((TH2*) infile->Get((TString) "surf_" + surfName +
514  (hie > 0 ? "NH" : "IH") +
515  "/profHists/hist" + std::to_string(i)));
516  }
517  hists.push_back(surf.ToTH2(minchi23));
518 
519  for(uint i = 0; i < hists.size(); ++i){
520  if(!hists[i]) continue;
521  if(TString(hists[i]->GetTitle()).Contains("delta")) {
523  hists[i]->GetZaxis()->SetRangeUser(0,2);
524  }
525  else gStyle->SetPalette(87);
526  if( i > 3) {
527  double minz = hists[i]->GetBinContent(hists[i]->GetMinimumBin());
528  double maxz = hists[i]->GetBinContent(hists[i]->GetMaximumBin());
529  double lim = std::max(fabs(minz), fabs(maxz));
530  hists[i]->GetZaxis()->SetRangeUser(-lim,+lim);
531  }
532  if(i==(hists.size()-1)) {
533  gStyle->SetPalette(56);
534  hists[i]->GetZaxis()->SetRangeUser(0.,11.83);
535  }
536 
537  hists[i]->Draw("colz");
538  hists[i]->GetYaxis()->SetTitleOffset(0.8);
539  gPad->SetRightMargin(0.14);
540  surf.DrawContour(surf_1Sigma, kSolid, kDarkColor,minchi23);
541  surf.DrawContour(surf_2Sigma, kSolid, kDarkColor,minchi23);
542  surf.DrawContour(surf_3Sigma, kSolid, kDarkColor,minchi23);
543  gPad->Print(outdir + "debug_contour_" + suffix + "_" +
544  surfName +
545  (hie > 0 ? "NH":"IH") +
546  (zoomIn ? "_zoom" : "") +
547  (std::to_string(i)) + ".pdf");
548 
549  }//end profiles
550  }//end debug
551  }//end hie
552  }//end surf
553  } // end plot only
554 }//end all
555 
556 
557 void DrawHieTag(int hie,bool th13, bool vert){
558  if(hie>0){
559  TLatex * ltxNH = new TLatex(0.16,0.89,"#bf{NH}");
560  ltxNH->SetTextColor(kPrimColorNH);
561  ltxNH->SetNDC();
562  ltxNH->SetTextAlign(22);
563  ltxNH->SetTextSize(kBlessedLabelFontSize);
564  if(th13) ltxNH->Draw();
565  else if(!vert) ltxNH->DrawLatex(.85,.20,"#bf{NH}");
566  else ltxNH->DrawLatex(.85,.85,"#bf{NH}");
567  }
568  else{
569  TLatex * ltxIH = new TLatex (0.16,0.46,"#bf{IH}");
570  ltxIH->SetTextColor(kPrimColorIH);
571  ltxIH->SetNDC();
572  ltxIH->SetTextAlign(22);
573  ltxIH->SetTextSize(kBlessedLabelFontSize);
574  if(th13) ltxIH->Draw();
575  else if(!vert) ltxIH->DrawLatex(.85,.20,"#bf{IH}");
576  else ltxIH->DrawLatex(.85,.85,"#bf{IH}");
577  }
578 }
579 
580 /// This one needs some work by hands with putting JoinGraphs properly for each contour
581 /// I didn't make it for the tutorial cases
582 void FillGraphs(std::vector<TGraph*> g1,std::vector<TGraph*> g2,std::vector<TGraph*> g3,
583  const Int_t kFillColor1, const Int_t kFillColor2, const Int_t kFillColor3,
584  const Int_t kDarkColor, const TString surftype)
585 {
586  bool isJoint = surftype.Contains("joint");
587  bool isBoth = surftype.Contains("both");
588  bool isRHC = surftype.Contains("RHCOnly");
589  bool isNH = surftype.Contains("NH");
590  bool isSysts = surftype.Contains("systs");
591  bool isFccorr = surftype.Contains("fccorr");
592 
593  if (surftype.Contains("delta_th23")){
594  cout<<"delta"<<endl;
595  if(isJoint && isNH ){
596  cout<<"joint nh"<<endl;
597  if(isRHC){
598  JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f");
599  JoinGraphs(g2[0], g2[1], kWhite)->Draw("f");
600  JoinGraphs(g2[1], g2[2], kWhite)->Draw("f");
601  JoinGraphs(g2[0], g2[1], kFillColor2)->Draw("f same");
602  JoinGraphs(g2[1], g2[2], kFillColor2)->Draw("f same");
603  }
604  else{
605  JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f");
606  JoinGraphs(g3[0], g3[1], kWhite)->Draw("l");
607  JoinGraphs(g2[0], g2[1], kWhite)->Draw("c f");
608  JoinGraphs(g2[0], g2[1], kFillColor2)->Draw("c f");
609  JoinGraphs(g2[0], g2[1], kWhite)->Draw("c");
610  }
611  if(isBoth){
612  for (auto & g:g1) { g->SetFillColor(kWhite); g->DrawClone("c f");}
613  for (auto & g:g1) { g->SetFillColor(kFillColor1); g->Draw("c f");}
614  for (auto & g:g1) { g->Draw("c");}
615  }
616  }
617 
618  if(isJoint && !isNH){
619  if(isRHC) JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f");
620  else{
621  JoinGraphs(g3[0], g3[2], kFillColor3)->Draw("f c");
622  g3[1]->SetFillColor(kFillColor3); g3[1]->Draw("f c");
623  g3[3]->SetFillColor(kFillColor3); g3[3]->Draw("f c");
624  for (auto &g:g3) {
625  g->SetLineWidth(3);
626  g->Draw("c");
627  }
628  }
629  for (auto &g:g2) {
630  g->SetFillColor(kWhite); g->DrawClone("f c");
631  g->SetFillColor(kFillColor2); g->Draw("f c");
632  g->SetLineColor(kDarkColor);
633  g->SetLineWidth(3);
634  g->Draw("c");
635  }
636  for (auto &g:g1) {
637  g->SetFillColor(kWhite); g->DrawClone("f c");
638  g->SetFillColor(kFillColor1); g->Draw("f c");
639  g->SetLineColor(kDarkColor);
640  g->SetLineWidth(3);
641  g->Draw("c");
642  }
643  for(auto & gs:{g3, g2, g1}){
644  for(auto & g:gs) {
645  g->SetLineColor(kDarkColor);
646  g->SetLineWidth(3);
647  g->Draw("c");
648  }
649  }
650  }
651  }
652 
653  if (surftype.Contains("th23_dm32")){
654  for (auto &g:g3) {
655  g->SetFillColor(kFillColor3); g->Draw("cf");
656  }
657  for (auto &g:g2) {
658  g->SetFillColor(kWhite); g->DrawClone("cf");
659  g->SetFillColor(kFillColor2); g->Draw("cf");
660  }
661  for (auto &g:g1) {
662  g->SetFillColor(kWhite); g->DrawClone("cf");
663  g->SetFillColor(kFillColor1); g->Draw("cf");
664  }
665  for (auto & gs:{g3,g2,g1}){
666  for (auto & g:gs) {
667  g->SetLineColor(kDarkColor);
668  g->SetLineWidth(3);
669  g->Draw("c");
670  }
671  }
672 
673  }
674 
675  if (surftype.Contains("th13_delta")) {
676  if(isBoth) {
677  JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f");
678  if(isNH) {
679  JoinGraphs(g2[0], g2[1], kWhite)->Draw("f");
680  JoinGraphs(g2[0], g2[1], kFillColor2)->Draw("f");
681  }
682  else{
683  for (auto & g:g2) { g->SetFillColor(kWhite); g->DrawClone("f");}
684  for (auto & g:g2) { g->SetFillColor(kFillColor2); g->Draw("f");}
685  }
686  for (auto & g:g1) { g->SetFillColor(kWhite); g->DrawClone("f");}
687  for (auto & g:g1) { g->SetFillColor(kFillColor1); g->Draw("f");}
688  }
689  else{
690  //for (auto &g:g3) {
691  // g->SetFillColor(kFillColor3); g->Draw("f");}
692  //JoinGraphs(g3[0], g3[1], kWhite)->Draw("f");
693  JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f c");
694  JoinGraphs(g2[0], g2[1], kWhite)->Draw("f c");
695  JoinGraphs(g2[0], g2[1], kFillColor2)->Draw("f c");
696  /*for (auto &g:g2) {
697  g->SetFillColor(kWhite); g->DrawClone("f");
698  g->SetFillColor(kFillColor2); g->Draw("f");
699  }*/
700  JoinGraphs(g1[0], g1[1], kWhite)->Draw("f c");
701  JoinGraphs(g1[0], g1[1], kFillColor1)->Draw("f c");
702  }
703 
704  for (auto & gs:{g3,g2,g1}){
705  for (auto & g:gs) {
706  g->SetLineColor(kDarkColor);
707  g->SetLineWidth(3);
708  g->Draw("l");
709  }
710  }
711  }
712 
713 }
714 
715 void DrawReactorConstraint(double mean, double err)
716 {
717  TBox* box = new TBox(mean-err, 0, mean+err, 2);
718  box->SetFillColorAlpha(kGray,0.7);
719 
720  TLine* linLo = new TLine(mean-err, 0, mean-err, 2);
721  linLo->SetLineColor(kGray+2);
722  linLo->SetLineWidth(2);
723 
724  TLine* linHi = new TLine(mean+err, 0, mean+err, 2);
725  linHi->SetLineColor(kGray+2);
726  linHi->SetLineWidth(2);
727 
728  box->Draw();
729  linLo->Draw();
730  linHi->Draw();
731  gPad->RedrawAxis();
732 }
733 
734 
735 std::vector <Spectrum * > GetMockData(TString ana_options){
736 
737  cout<<"\n\n ============= You're loading Mock data from the file ============="<<endl;
738  std::vector <Spectrum *> spec;
739  std::string Dir = "./";
740  bool RHCmode = ana_options.Contains("rhc");
741  bool FHCmode = ana_options.Contains("fhc");
742  bool NueData = ana_options.Contains("nue");
743  bool NumuData = ana_options.Contains("numu");
744  if (NueData){
745  if (FHCmode) spec.push_back(LoadFromFile<Spectrum>(Dir+"mock_data_for_yn_tutorial.root", "data_exp_0").release());
746  if (RHCmode) spec.push_back(LoadFromFile<Spectrum>(Dir+"mock_data_for_yn_tutorial.root", "data_exp_1").release());
747  }
748  if (NumuData){
749  if (FHCmode){
750  for (int i = 2; i < 6; ++i) {
751  auto temp = LoadFromFile<Spectrum>(Dir+"mock_data_for_yn_tutorial.root", "data_exp_"+std::to_string(i)).release();
752  spec.push_back(temp);
753  }
754  }
755  if (RHCmode){
756  for (int i = 6; i < 10; ++i) {
757  auto temp = LoadFromFile<Spectrum>(Dir+"mock_data_for_yn_tutorial.root", "data_exp_"+std::to_string(i)).release();
758  spec.push_back(temp);
759  }
760  }
761 
762  }
763  return spec;
764 
765 }
Spectrum * GetFakeData(const IPrediction *pred, osc::IOscCalc *calc, const double pot, TH1D *cosmics=0)
const Color_t kPrimColorIH
Definition: Style.h:64
T max(const caf::Proxy< T > &a, T b)
const Color_t kPrimColorNH
Definition: Style.h:61
void FillGraphs(std::vector< TGraph * > g1, std::vector< TGraph * > g2, std::vector< TGraph * > g3, const Int_t kFillColor1, const Int_t kFillColor2, const Int_t kFillColor3, const Int_t kDarkColor, const TString surftype)
void Tutorial2019FitContours(bool createFile=false, bool corrSysts=true, TString options="joint_realData_both_onlyIH", bool dmsqSurf=false, bool th13Surf=false, bool zoomIn=true, bool fillContour=false)
Start the Demo script (it&#39;s a bit simplified nue / Ana2018 / FitandFC / joint_fit_2018_contours.C):
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
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)
const Dmsq32Constraint kDmsq32ConstraintPDG2017(2.45e-3, 0.05e-3, 2.52e-3, 0.05e-3)
const FitSinSqTheta23UpperOctant kFitSinSqTheta23UpperOctant
Definition: FitVars.cxx:16
TH2 * Gaussian68Percent2D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 2D in gaussian approximation.
OStream cerr
Definition: OStream.cxx:7
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.
TString hists[nhists]
Definition: bdt_com.C:3
static SystShifts Nominal()
Definition: SystShifts.h:34
std::vector< Spectrum * > GetMockData(TString ana_options)
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
std::unique_ptr< T > LoadFrom(TDirectory *dir, const std::string &label)
Definition: LoadFromFile.h:17
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
tuple white
Definition: rootlogon.py:71
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
TH2 * Gaussian3Sigma2D(const FrequentistSurface &s)
Up-value surface for 3 sigma confidence in 2D in gaussian approximation.
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
virtual T GetDmsq32() const
Definition: IOscCalc.h:86
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
void SetFakeCalc(osc::IOscCalcAdjustable *calc, double ssth23=-5, double dmsq32=-5, double dCP_Pi=-5)
void Preliminary()
void DrawContourCanvas(TString surfName, double minx=0, double maxx=2, double miny=0, double maxy=1)
TGraph * JoinGraphs(TGraph *a, TGraph *b, int fillcol)
Join graphs and set a fill color. Useful for contours.
Definition: Plots.cxx:1529
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
virtual T GetTh13() const
Definition: IOscCalc.h:88
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.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void DrawReactorConstraint(double mean, double err)
const Style_t k2SigmaConfidenceStyle
Definition: Style.h:71
std::vector< std::pair< TH1D *, double > > GetNumuCosmics2019(const int nq=4, std::string beam="fhc", bool GetFromUPS=false, bool NERSC=false)
const FitSinSqTheta23LowerOctant kFitSinSqTheta23LowerOctant
Definition: FitVars.cxx:17
void ResetToNominal()
Definition: SystShifts.cxx:141
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
assert(nhit_max >=nhit_nbins)
Interface definition for fittable variables.
Definition: IFitVar.h:16
double GetPOT(bool isFHC=true)
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const Float_t kBlessedLabelFontSize
Definition: Style.h:90
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
std::pair< Spectrum *, double > cos
void DrawHieTag(int hie, bool th13=true, bool vert=false)
Float_t e
Definition: plot.C:35
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
double th13
Definition: runWimpSim.h:98
void SetPaletteDeltaCyclicNew()
const std::string outdir
const Color_t k2SigmaConfidenceColorIH
Definition: Style.h:82
FILE * outfile
Definition: dump_event.C:13
std::pair< TH1D *, double > GetNueCosmics2019(std::string beam, bool GetFromUPS=false, bool NERSC=false)
Compare a single data spectrum to the MC + cosmics expectation.
surf
Definition: demo5.py:35
unsigned int uint