demoFitContours.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 "Utilities/rootlogon.C"
18 #include "OscLib/IOscCalc.h"
19 
20 #include "CAFAna/nue/Ana2018/joint_fit_2018_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 void DrawHieTag(int hie, bool th13=true, bool vert=false);
44 
45 void FillGraphs(std::vector<TGraph*> g1,std::vector<TGraph*> g2,std::vector<TGraph*> g3,
46  const Int_t kFillColor1, const Int_t kFillColor2, const Int_t kFillColor3,
47  const Int_t kDarkColor, const TString surftype);
48 
49 void DrawReactorConstraint(double mean, double err);
50 
51 // kCMYK analog
53  {
54  Double_t white[9] = { 1, 1, 1, 1, 1, 1, 1, 1, 1 };
55  Double_t Red[9] = { 61./255., 99./255., 136./255., 181./255., 213./255., 225./255., 198./255., 99./255., 61./255.};
56  Double_t Green[9] = { 149./255., 140./255., 96./255., 83./255., 132./255., 178./255., 190./255., 140./255., 149./255};
57  Double_t Blue[9] = { 214./255., 203./255., 168./255., 135./255., 110./255., 100./255., 111./255., 203./255., 214./255.};
58  Double_t Length[9] = { 0.0000, 0.1250, 0.2500, 0.3750, 0.5000, 0.6250, 0.7500, 0.8750, 1.0000};
59 
60  TColor::CreateGradientColorTable(9, Length, Red, Green, Blue, 255, 1.0);
61  gStyle->SetNumberContours(255);
62  }
63 
64 /// Start the Demo script (it's a bit simplified nue / Ana2018 / FitandFC / joint_fit_2018_contours.C):
65 
66 void demoFitContours(bool createFile=false, /// Do you want to fit and write it into the file (true) or make plots (false)?
67  bool corrSysts = true, /// Stat only (false) or with systs (true)?
68  TString options="joint_both_onlyIH", /// Write exactly what fit is it? joint/nueOnly/numuOnly + both/FHCOnly/RHCOnly + onlyIH/onlyNH
69  bool dmsqSurf = false, /// What kind of plot do you want, is it dmsq vs th23 surface?
70  bool th13Surf = false, /// What kind of plot do you want, is it th13 vs delta CP surface?
71  /// If dmsqSurf = false && th13Surf = false then it will be th23 vs delta CP surface
72  bool zoomIn = true, /// Since it's pretty universal code axes are very extanded for cases of nueOnly || numuOnly || RHCOnly || FHCOnly.
73  /// Thus sometimes you'll want to zoom axis (especially in case of FHC&&RHC)
74  bool fillContour = false) /// Do you want to fill it with color or stay with border lines only?
75 {
76  /// Some general info about further fit
77 
78  bool nueOnly = options.Contains("nueOnly");
79  bool numuOnly = options.Contains("numuOnly");
80  bool joint = options.Contains("joint");
81  assert (nueOnly || numuOnly || joint);
82 
83  bool FHCOnly = options.Contains("FHCOnly");
84  bool RHCOnly = options.Contains("RHCOnly");
85  bool both = options.Contains("both");
86  assert (FHCOnly || RHCOnly || both);
87 
88  bool onlyNH = options.Contains("onlyNH");
89  bool onlyIH = options.Contains("onlyIH");
90 
91  auto suffix = options;
92  if(dmsqSurf) suffix+= "_dmsq";
93  if(th13Surf) suffix+= "_th13";
94  if(corrSysts) suffix+="_systs";
95  else suffix+="_stats";
96 
97  ///How do you want to organize and sort your outputs?
98 
99  /*TString outdir="";
100 
101  if(th13Surf) outdir += "th13_delta/";
102  else if (dmsqSurf) outdir += "th23_dmsq/";
103  else outdir += "delta_th23/";
104 
105  TString outfilename (outdir + (corrSysts?"syst/":"stat/") + "hist_contours_2018_" + suffix);*/
106 
107  /// For tute purposes save all plots in current dir
108  TString outdir="";
109  TString outfilename (outdir + "hist_contours_2018_" + suffix);
110 
111  /// Enter the fit and file creation part
112 
113  if(createFile){
114 
115 //////////////////////////////////////////////////
116 // Load Nue and Numu experiments
117 //////////////////////////////////////////////////
118 
119  /// Create vectors for all predictions we need - MC predictions, cosmic predictions, data or fake data
120  /// Then we'll need vectors of experiments for making a multiexperiment
121  /// And we'll need separate set of all items for making numu-only preseeds (since we follow 2018 ana procedure)
122 
123  std::vector <const IPrediction * > preds;
124  std::vector <const IPrediction * > preds_numu_only;
125  std::vector <std::pair <TH1D*, double > > cosmics;
126  std::vector <std::pair <TH1D*, double > > cosmics_numu_only;
127  std::vector <Spectrum * > data;
128  std::vector <Spectrum * > data_numu_only;
129  std::vector <const IExperiment * > expts;
130  std::vector <const IExperiment * > expts_numu_only;
131 
132  auto calc_fake = DefaultOscCalc();
133  SetFakeCalc(calc_fake, 0.58, 2.51e-3, 0.17); /// use 2018 BF for making fake data
134 
135  /// And fill all vectors we created, the structure is:
136  // nue FHC, nue RHC, 4 numu FHC quartiles, 4 numu RHC quartiles - for the joint_both option selected
137  // If you selected numuOnly|nueOnly|RHCOnly|FHCOnly then some of these predictions will be switched off
138 
139  if(!numuOnly) { /// fill nue part
140  if(FHCOnly || both ) {
141  preds.push_back(GetNuePrediction2018("combo", DefaultOscCalc(), corrSysts, "fhc", false)); // nue FHC
142  cosmics.push_back(GetNueCosmics2018("fhc"));
143  }
144  if(RHCOnly || both ) {
145  preds.push_back(GetNuePrediction2018("prop", DefaultOscCalc(), corrSysts, "rhc", false)); // nue FHC
146  cosmics.push_back(GetNueCosmics2018("rhc"));
147  }
148  }
149 
150  int nnumu = 4;
151  if(!nueOnly) { /// fill numu part and numu_only for preseeds
152  if(FHCOnly || both ) {
153 
154  auto numu_preds = GetNumuPredictions2018(nnumu, corrSysts, "fhc"); // numu FHC
155  preds.insert(preds.end(),numu_preds.begin(), numu_preds.end());
156  preds_numu_only.insert(preds_numu_only.end(),numu_preds.begin(), numu_preds.end());
157 
158  auto numu_cosmics = GetNumuCosmics2018(nnumu, "fhc");
159  cosmics.insert(cosmics.end(),numu_cosmics.begin(), numu_cosmics.end());
160  cosmics_numu_only.insert(cosmics_numu_only.end(),numu_cosmics.begin(), numu_cosmics.end());
161 
162  }
163  if(RHCOnly || both ) {
164 
165  auto numu_preds = GetNumuPredictions2018(nnumu, corrSysts, "rhc"); // numu RHC
166  preds.insert(preds.end(),numu_preds.begin(), numu_preds.end());
167  preds_numu_only.insert(preds_numu_only.end(),numu_preds.begin(), numu_preds.end());
168 
169  auto numu_cosmics = GetNumuCosmics2018(nnumu, "rhc");
170  cosmics.insert(cosmics.end(),numu_cosmics.begin(), numu_cosmics.end());
171  cosmics_numu_only.insert(cosmics_numu_only.end(),numu_cosmics.begin(), numu_cosmics.end());
172  }
173  }
174 
175  /// Now assign proper POT for each of these predictions following the same structure
176  /// nue FHC, nue RHC, 4 numu FHC and 4 numu RHC quartiles
177  /// end make fake data. Feed everything to the experiments
178 
179  cout<<"\n\n Predictions used for the fit:\n";
180 
181  for(int i = 0; i < int(preds.size()); ++i){
182  double POT;
183  if(FHCOnly) POT = kAna2018FHCPOT;
184  if(RHCOnly) POT = kAna2018RHCPOT;
185  if(both) {
186  if(nueOnly || both) {
187  if(i==0) POT = kAna2018FHCPOT;
188  if(i==1) POT = kAna2018RHCPOT;
189  }
190  if(numuOnly) {
191  if(i >= 0 && i < 4) POT = kAna2018FHCPOT;
192  if(i >= 4 && i < 8) POT = kAna2018RHCPOT;
193  }
194  if(both) {
195  if(i >= 2 && i < 6) POT = kAna2018FHCPOT;
196  if(i >= 6 && i < 10) POT = kAna2018RHCPOT;
197  }
198  }
199  /// done with POT, let's make fake data
200  auto thisdata = GetFakeData (preds[i], calc_fake, POT, cosmics[i].first);
201  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;
202  expts.push_back(new SingleSampleExperiment(preds[i],*thisdata, cosmics[i].first,cosmics[i].second)); // make an experiment with current predictions and data
203  /// for debug purposes print spectra to the file
204  auto debug=true;
205  if(debug){
206  new TCanvas();
207  DataMCComparison (*thisdata, preds[i], DefaultOscCalc());
208  cosmics[i].first->SetMarkerStyle(34);
209  cosmics[i].first->SetMarkerColor(kCosmicBackgroundColor);
210  cosmics[i].first->Draw("hist p same");
211  //gPad->Print(outdir + (corrSysts?"syst/":"stat/")+ "debug_predictions_" + suffix + "_" + std::to_string(i) + ".pdf");
212  gPad->Print(outdir + "debug_predictions_" + suffix + "_" + std::to_string(i) + ".pdf");
213  }//end debug
214  }
215 
216  /// Out main preporation part is done - we have 10 experiments for fit
217  /// one more step - do the same but with numu only predictions for pre-seeding
218 
219  cout<<"\n\nMake numu only experiment to get the seeds"<<endl;
220  //Follow the same scheme but with reduced number of predictions:
221  for(int i = 0; i < (int) preds_numu_only.size(); ++i){
222  double POT;
223  if(FHCOnly) POT = kAna2018FHCPOT;
224  if(RHCOnly) POT = kAna2018RHCPOT;
225  if(both) {
226  if(i >= 0 && i < 4) POT = kAna2018FHCPOT;
227  if(i >= 4 && i < 8) POT = kAna2018RHCPOT;
228  }
229  auto thisdata = GetFakeData (preds_numu_only[i],calc_fake, POT, cosmics_numu_only[i].first);
230  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;
231  expts_numu_only.push_back(new SingleSampleExperiment(preds_numu_only[i],*thisdata, cosmics_numu_only[i].first,cosmics_numu_only[i].second));
232  }
233 
234 //////////////////////////////////////////////////
235 /// Add constraints, make experiments
236 //////////////////////////////////////////////////
237 
238  /// Lets make a multiexperiment for simultaneous fit
239  /// For different contours/fits we will use different additional information
240 
241  std::cout << "\nAll experiments are done \nCreating multiexperiment\n" << std::endl;
242  if(!th13Surf){ /// We use reactor constraint as 11th "experiment" - th13 will be varied in PDG's error range with PDG central value
243  /// but we wouldn't use such constraint if we want to make our own measurement, right?
244  std::cout << "Adding WorldReactorConstraint2017\n";
245  expts.push_back(WorldReactorConstraint2017());
246  }
247  if(nueOnly) { /// If we don't use NOvA's result in the fit, then we will use PDG's constraint
248  std::cout << "Adding Dmsq32ConstraintPDG2017\n";
249  expts.push_back(&kDmsq32ConstraintPDG2017);
250  }
251  /// Let's make a multiexperiment
252  std::cout << "Creating Multiexperiment with a total of "
253  << expts.size() << " experiments\n\n" << std::endl;
254  auto exptThis = new MultiExperiment(expts);
255  /// and multiexperiment for numu pre-seeds
256  std::cout << "Creating Multiexperiment of numu only SimpleExp with a total of "
257  << expts_numu_only.size() << " experiments\n\n" << std::endl;
258  auto exptThis_numu_only = new MultiExperiment(expts_numu_only);
259 
260  std::cout<<"Clean verctors"<<std::endl;
261  preds.clear();
262  preds_numu_only.clear();
263 
264 ////////////////////////////////////////////////////////////
265 /// Systematics
266 ////////////////////////////////////////////////////////////
267 
268  /// Now, let's prepare systematics
269  std::vector< const ISyst * > systs = {};
270  std::vector< const ISyst * > systs_numu_only = {};
271  if (corrSysts){ // syst case
272  systs = {&kAnaCalibrationSyst}; //here will be your own syst, which you made previously
273  systs_numu_only = {&kAnaCalibrationSyst}; //here will be your own syst, which you made previously
274  }
275 
276  /// And that's how it looks like in the analysis - you need to set proper correlations for all simple experiments
277 
278  /*std::cout << "Systematics for the fit:\n";
279  auto systs = GetJointFitSystematicList(corrSysts, numuOnly, nueOnly, true, true, true);
280 
281  std::cout << "\n\nSystematic correlations...\n";
282  if(!nueOnly && ! numuOnly && corrSysts){
283  if(FHCOnly){
284  exptThis->SetSystCorrelations(0, GetCorrelations(true, true));
285  auto notfornumu = GetCorrelations(false, true);
286  for(int i =0; i < nnumu; ++i) exptThis->SetSystCorrelations(i+1, notfornumu);
287  }
288  if(RHCOnly){
289  exptThis->SetSystCorrelations(0, GetCorrelations(true, false));
290  auto notfornumu = GetCorrelations(false, true);
291  for(int i =0; i < nnumu; ++i) exptThis->SetSystCorrelations(i+1, notfornumu);
292  }
293  if(both){
294  exptThis->SetSystCorrelations(0, GetCorrelations(true, true));
295  exptThis->SetSystCorrelations(1, GetCorrelations(true, false));
296  auto notfornumu = GetCorrelations(false, true);
297  for(int i =0; i < 8; ++i) exptThis->SetSystCorrelations(i+2, notfornumu);
298  }
299  auto systs_numu_only = GetJointFitSystematicList(corrSysts, false, true, true, true, true);
300  }
301  */
302 
303 ////////////////////////////////////////////////////////////
304 /// Fit
305 ////////////////////////////////////////////////////////////
306  /// And finally let's start fitting
307 
308  std::cout << "Starting the fit" << std::endl;
309 
311 
312  SystShifts auxShifts = SystShifts::Nominal();
313  // In case some systs need different seeds
314  std::vector <SystShifts> seedShifts = {};
315 
316 //////////////////////////////////////////////////////////////////////
317 /////////////////////////// Seed controller ////////////////////////////
318 ////////////////////////////////////////////////////////////////////////
319  /// This section contains all seeds that will be used further
320  /// There is numu pre-seed part here as well. Idea is to choose th23 and dmsq23 seeds
321  /// from the numu only fast fit. It is more important for FC, but for similarity we put this
322  /// part in gaussian fit as well
323 
324  cout<<"------------------- Start prestage seeds --------------------------"<<endl;
325 
326  double minchi_numu = 1E20;
327  double pre_seed_th23;
328  double pre_seed_dmsq;
329  ResetOscCalcToDefault(calc);
330  auxShifts.ResetToNominal();
331 
332  double maxmix = 0.514; // from the numu results
333  double numu_pre_seedLONH, numu_pre_seedUONH, numu_pre_seedLOIH, numu_pre_seedUOIH, dmsq_numu_pre_seedNH, dmsq_numu_pre_seedIH;
334  /// it simply searches BF in NH and use BF th23 and dmsq values as seeds
335  for(int hie:{1}){
336  std::cout << "\n\nFinding test best fit " << (hie>0? "NH " : "IH ") << "\n\n";
337  Fitter fitnumu_only(exptThis_numu_only, {&kFitSinSqTheta23, &kFitDmSq32Scaled}, {});
338 
339  auto thisminchi = fitnumu_only.Fit(calc, auxShifts ,
340  {{&kFitDmSq32Scaled,{hie*2.5}},
341  {&kFitSinSqTheta23, {0.4} }}, {},// seedShifts
342  Fitter::kVerbose);
343 
344  if(thisminchi < minchi_numu) minchi_numu = thisminchi;
345  }
346  /// our BF values
347  pre_seed_th23 = kFitSinSqTheta23.GetValue(calc);
348  pre_seed_dmsq = kFitDmSq32Scaled.GetValue(calc);
349  /// make two seeds in LO and UO with the help of pre_seed_dmsq
350  numu_pre_seedLONH = ((pre_seed_th23>maxmix)?(2*maxmix-pre_seed_th23):pre_seed_th23);
351  numu_pre_seedUONH = ((pre_seed_th23>maxmix)?pre_seed_th23:(2*maxmix-pre_seed_th23));
352 
353  ResetOscCalcToDefault(calc);
354  auxShifts.ResetToNominal();
355 
356  cout<<"------------------- End prestage seeds --------------------------"<<endl;
357 
358  struct th23helper{
359  std::vector<double> seeds;
360  const IFitVar * var;
361  TString label;
362  };
363 
364  std::vector <th23helper> th23seeds;
365  //for NH:
366  /// make two seeds in both UO and LO (near maxmixing and th23 from pre-seeding stage)
367  th23seeds.push_back( { {0.499, numu_pre_seedLONH}, kFitSinSqTheta23LowerOctant, "LO"});
368  th23seeds.push_back( { {0.501, numu_pre_seedUONH}, kFitSinSqTheta23UpperOctant, "UO"});
369  /// Put three seeds - LO, UO, maxmixing
370  std::vector<double> th23_all_seeds_NH = {numu_pre_seedLONH, 0.5, numu_pre_seedUONH};
371  /// delta CP seeds
372  std::vector<double> delta_seeds = {0, 0.5, 1., 1.5};
373 
374 /////////////////////////////////////////////////////////////////////////
375 //////////////////////////// Different best fit searches //////////////////
376 ///////////////////////////////////////////////////////////////////////////
377  /// Now, let's find BF value, it searches BF in IH/NH and UO/LO and chooses the best chisq
378  double minchi23 = 1E20;
379  for(int hie:{-1,1}){
380  /// will use th23seeds which we defined previously for LO and UO
381  for (auto & thseed: th23seeds){
382  std::cout << "\n\nFinding best fit "
383  << thseed.label
384  << (hie>0? " NH " : " IH ")
385  << "\n\n";
386  for (auto const & s:thseed.seeds) std::cout << s << ", ";
387  std::cout << std::endl;
388  cout<<"pre seed from numu is "<<pre_seed_dmsq<<endl;
389  /// replace th23 variable by kFitSinSqTheta23LowerOctant/kFitSinSqTheta23UpperOctant
390  // We minimize deltaCP, th23, dmsq23, th13 (which is constrained)
391  std::vector <const IFitVar*> fitvars = {&kFitDeltaInPiUnits,
392  thseed.var,
395 
396  Fitter fit23(exptThis, fitvars, systs);
397  /// Minimize and use seeds for dmsq (the one from pre-seed stage * 1 or -1 for NH and IH)
398  /// th23 will use seeds from prestage for kFitSinSqTheta23LowerOctant/kFitSinSqTheta23UpperOctant
399  /// th13 seed is PDG's central value
400  /// deltaCP seeds at 0, pi/2, pi, 3pi/2
401  auto thisminchi = fit23.Fit(calc,auxShifts ,{
402  {&kFitDmSq32Scaled,{hie*pre_seed_dmsq}},
403  {thseed.var, thseed.seeds},
404  {&kFitDeltaInPiUnits, delta_seeds},
405  {&kFitSinSq2Theta13, {calc->GetTh13()}} },
406  {}, Fitter::kVerbose
407  );
408  /// compare with global BF
409  minchi23= min(minchi23, thisminchi);
410  ResetOscCalcToDefault(calc);
411  auxShifts.ResetToNominal();
412  }//end th23
413  }//end hie
414  std::cout << "\nFound overall minchi " << minchi23 << std::endl;
415 
416  //We'll save to this file, make pretty plots later
417  TFile * outfile = new TFile(outfilename+".root","recreate");
418 
419  outfile->cd();
420  TVectorD v(1);
421  v[0] = minchi23;
422  cout<<"minchi wrote into the file "<<minchi23<<endl;
423  v.Write("overall_minchi");
424 
425  /// How much point on the grid do you want? It will be steps*steps
426  int steps = 30;//10;// 5;
427  if(dmsqSurf){
428  std::cerr << "\n WARNING Using 20 bins for dmsq surface\n\n";
429  steps = 20;
430  }
431  //Now create all the surfaces
432  /// Generally need both IH and NH, but it's too long to run them all together
433  for(int hie: {-1, +1}){
434 
435  if((onlyNH && hie<0) || (onlyIH && hie>0)) continue; /// Will run only one hie which you chose at the beginning
436 
437  std::cout << "Starting surface " << (hie>0? "NH ": "IH") << "\n\n";
438  cout<<"the delta seeds are ";
439  for (auto const & s:delta_seeds) std::cout << s << ", ";
440  std::string hieStr = hie>0 ? "NH" : "IH";
441 
442  /// It will be deltaCP vs Th23 surface
443  if(!th13Surf && ! dmsqSurf){
444  ResetOscCalcToDefault(calc);
445  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
446  /// Set Y axis's low and up values
447  double low = (nueOnly||RHCOnly)?0:0.3;
448  double high = (nueOnly||RHCOnly)?1:0.7;
449  cout<<"low "<<low<<" high "<<high<<endl;
450  /// It will minimize dmsq23 and th13 only, no seeds here
451  FrequentistSurface surf23(exptThis, calc,
452  &kFitDeltaInPiUnits,steps, 0, 2,
453  &kFitSinSqTheta23, steps, low, high,
455  //{},seedShifts
456  );
457 
458  auto surf1=surf23.ToTH2(minchi23);
459  outfile->cd();
460  surf1->Write((TString)"delta_th23_"+hieStr);
461  surf23.SaveTo(outfile->mkdir((TString)"surf_delta_th23_"+hieStr));
462  }
463 
464  /// This part will make th23 vs dmsq surface
465  if(!th13Surf && dmsqSurf){
466  calc = DefaultOscCalc();
467  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
468  /// Put Y axis's low and up values
469  double low = (nueOnly||RHCOnly)?0.2:0.35;
470  double high = (nueOnly||RHCOnly)?1:0.7;
471  cout<<"low "<<low<<" high "<<high<<endl;
472  /// It will minimize deltaCP and th13, use seeds for deltaCP
473  /// also, there is some logic for x axis's values as well
474  /// for different cases of fit (RHC/FHC/both) contours can be larger or smaller. To much pain in these values ...
475  FrequentistSurface surf23m(exptThis, calc,
476  &kFitSinSqTheta23, steps,low, high,
477  &kFitDmSq32Scaled, steps,
478  (hie>0?(RHCOnly?2.0:2.2):(RHCOnly?-3.5:-2.8)),(hie>0?(RHCOnly?3.2:2.85):(RHCOnly?-2.0:-2.3)),
480  {{&kFitDeltaInPiUnits, delta_seeds}});
481  auto surf6=surf23m.ToTH2(minchi23);
482  outfile->cd();
483  surf6->Write((TString)"th23_dm32_"+(hie>0?"NH":"IH"));
484  surf23m.SaveTo(outfile->mkdir((TString)"surf_th23_dm32_"+(hie>0?"NH":"IH")));
485  }
486 
487  /// Also, we can do th13 vs deltaCP contour, why not
488  /// (we don't show them to public since there are some issues with FC for th13 measurements)
489  /// but gaussian fit is ok
490  if(th13Surf){
491  calc = DefaultOscCalc();
492  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
493  /// Minimize dmsq23, th23, put seeds for th23 as (smth in LO, 0.5, smth in UO)
494  FrequentistSurface surf13(exptThis, calc,
495  &kFitSinSq2Theta13, steps, 0, (RHCOnly?1:0.35),
496  &kFitDeltaInPiUnits,steps, 0, 2,
498  {{&kFitSinSqTheta23, th23_all_seeds_NH}});//, seedShifts);
499  auto surf4 = surf13.ToTH2(minchi23);
500  outfile->cd();
501  surf4->Write((TString)"th13_delta_"+(hie>0?"NH":"IH"));
502  surf13.SaveTo(outfile->mkdir((TString)"surf_th13_delta_"+(hie>0?"NH":"IH")));
503 
504  }
505  } //end hie
506 
507  //delete outfile;
508  std::cout << "\nSurfaces saved to " << outfilename << std::endl;
509  return;
510  }//end createfile
511 
512 /// Done with fit! Now play with style and make pretty plots.
513 
514 //////////////////////////////////////////////////////////////////////
515 // Make plots
516 //////////////////////////////////////////////////////////////////////
517 
518  else{ //plot only
519 
520  /// read file from the previous section
521  TFile * infile = new TFile (outfilename+".root","read");
522  /// Get BF
523  auto mins =* (TVectorD*)infile->Get("overall_minchi");
524  double minchi23 = mins[0];
525  /// Some logic to select contour which will you plot
526  std::vector <TString> surfNames;// = {"delta_th23_"};<
527  if(!th13Surf && ! dmsqSurf)surfNames = {"delta_th23_"};
528  if (dmsqSurf) surfNames.push_back("th23_dm32_");
529  if (th13Surf) surfNames.push_back("th13_delta_");
530  for(TString surfName:surfNames){
531  /// again, need two hierarchies, but all files were done with just one
532  for (int hie:{-1,+1}){
533 
534  if((onlyNH && hie<0) || (onlyIH && hie>0)) continue;
535 
536  //Setup canvas for each contour type, zoom if necessary
537  if(surfName.Contains("delta_th23")) {
538  if (zoomIn) DrawContourCanvas(surfName, 0., 2., (RHCOnly?0.15:0.275), (RHCOnly?0.8:0.725));
539  else DrawContourCanvas(surfName, 0., 2., 0., 1.);
540  }
541  if(surfName.Contains("th23_dm32")) {
542  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));
543  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));
544  }
545  if(surfName.Contains("th13_delta")) {
546  if (zoomIn) DrawContourCanvas(surfName, 0., 0.38, 0, 2);
547  else DrawContourCanvas(surfName, 0., 1.0, 0, 2);
548  }
549  /// Read your surface and set 1,2,3 sigma curves
550  auto surf = *FrequentistSurface::LoadFrom(infile->GetDirectory("surf_"+surfName+(hie>0?"NH":"IH")));
551  TH2 * hsurf = (TH2*)infile->Get(surfName+(hie>0?"NH":"IH"));
552  TH2 * surf_1Sigma = Gaussian68Percent2D(surf);
553  TH2 * surf_2Sigma = Gaussian2Sigma2D(surf);
554  TH2 * surf_3Sigma = Gaussian3Sigma2D(surf);
555  /// colors
556  Color_t k2SigmaConfidenceColor = ((hie>0)? k2SigmaConfidenceColorNH:k2SigmaConfidenceColorIH);
557  Int_t kFillColor1 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.75);
558  Int_t kFillColor2 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.38);
559  Int_t kFillColor3 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.18);
560  Int_t kFCFillColor1 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.75);
561  Int_t kFCFillColor2 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.38);
562  Int_t kFCFillColor3 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.18);
563  Int_t kDarkColor = kGray+3;
564 
565  if(fillContour){
566  auto g1 = surf.GetGraphs(surf_1Sigma,minchi23);
567  auto g2 = surf.GetGraphs(surf_2Sigma,minchi23);
568  auto g3 = surf.GetGraphs(surf_3Sigma,minchi23);
569  FillGraphs(g1, g2, g3,
570  kFillColor1, kFillColor2, kFillColor3,
571  kDarkColor, suffix + surfName + (hie>0?"NH":"IH"));
572  }
573  else{
574  cout<<"minchi from file is "<<minchi23<<endl;
575 
576  surf.DrawContour(surf_1Sigma, kSolid, kFillColor1,minchi23);
577  surf.DrawContour(surf_2Sigma, k2SigmaConfidenceStyle, kFillColor2,minchi23);
578  surf.DrawContour(surf_3Sigma, kSolid, kFillColor3,minchi23);
579  }
580 
581  gPad->RedrawAxis();
582 
583  /// Make pretty point at BF value
584  bool bestfit = 1;
585  double delta, th23, dmsq32, th13;
586  if(bestfit && both && corrSysts && hie>0){
587  if(!th13Surf){
588  delta = 0.16576;
589  th23 = 0.584509;
590  dmsq32 = 2.50515;
591  th13 = 0.0824663;
592  }
593  else{delta = 0.678516; th13 = 0.0957894;}
594  TMarker bf;
595  bf.SetMarkerColor(kDarkColor);
596  // bf.SetMarkerStyle(kFullStar);
597  bf.SetMarkerStyle(43);
598  bf.SetMarkerSize(2);
599  if(!dmsqSurf && ! th13Surf){
600  bf.SetX(delta);
601  bf.SetY(th23);
602  bf.DrawClone();
603  }
604  if(dmsqSurf){
605  bf.SetX(th23);
606  bf.SetY(dmsq32);
607  bf.DrawClone();
608  }
609  if(th13Surf){
610  bf.SetX(th13);
611  bf.SetY(delta);
612  bf.DrawClone();
613  }
614  }
615  /// Set legend
616  auto leg = ContourLegend(hie, fillContour, false,
617  kFillColor1,kFillColor2,kFillColor3,kDarkColor,
618  bestfit && both);
619  /// Put PDG's constraint line on the th13 vs deltaCP surface
620  if(th13Surf) {
621  DrawReactorConstraint(0.082,0.004);
622  leg->SetFillStyle(1001);
623  leg->SetNColumns(1);
624  leg->SetX1(0.7); leg->SetX2(0.8);
625  leg->SetY1(0.4); leg->SetY2(0.85);
626  leg->SetHeader("No FC");
627 
628  TGraph * dummy = new TGraph;
629  dummy->SetLineWidth(2);
630  dummy->SetLineColor(kGray+2);
631  dummy->SetFillColor(kGray);
632  leg->AddEntry(dummy,"#splitline{Reactor}{68%C.L.}","f");
633  }
634  leg->Draw();
635  Preliminary();
636  DrawHieTag(hie,false);
637 
638  /// and finally (!) save the plot to the file
639  for(TString ext: {".pdf",".root"}){
640  gPad->Print(outdir + "contour_" + suffix + "_" +
641  surfName +
642  (hie>0?"NH":"IH") +
643  (zoomIn?"_zoom":"") +
644  ext);
645  }//end ext
646 
647  /// in case if you want (of course you want) to see debug plots for each syst or osc parameter put true here
648  auto debug = true;
649  if(debug){
650  auto hists = surf.GetProfiledHists();
651  int nSysts = corrSysts ? 51 : 0;
652  for (int i = 0; i < nSysts; ++i){
653  hists.push_back((TH2*) infile->Get((TString) "surf_" + surfName +
654  (hie > 0 ? "NH" : "IH") +
655  "/profHists/hist" + std::to_string(i)));
656  }
657  hists.push_back(surf.ToTH2(minchi23));
658 
659  for(uint i = 0; i < hists.size(); ++i){
660  if(!hists[i]) continue;
661  if(TString(hists[i]->GetTitle()).Contains("delta")) {
663  hists[i]->GetZaxis()->SetRangeUser(0,2);
664  }
665  else gStyle->SetPalette(87);
666  if( i > 3) {
667  double minz = hists[i]->GetBinContent(hists[i]->GetMinimumBin());
668  double maxz = hists[i]->GetBinContent(hists[i]->GetMaximumBin());
669  double lim = std::max(fabs(minz), fabs(maxz));
670  hists[i]->GetZaxis()->SetRangeUser(-lim,+lim);
671  }
672  if(i==(hists.size()-1)) {
673  gStyle->SetPalette(56);
674  hists[i]->GetZaxis()->SetRangeUser(0.,11.83);
675  }
676 
677  hists[i]->Draw("colz");
678  hists[i]->GetYaxis()->SetTitleOffset(0.8);
679  gPad->SetRightMargin(0.14);
680  surf.DrawContour(surf_1Sigma, kSolid, kDarkColor,minchi23);
681  surf.DrawContour(surf_2Sigma, kSolid, kDarkColor,minchi23);
682  surf.DrawContour(surf_3Sigma, kSolid, kDarkColor,minchi23);
683  gPad->Print(outdir + "debug_contour_" + suffix + "_" +
684  surfName +
685  (hie > 0 ? "NH":"IH") +
686  (zoomIn ? "_zoom" : "") +
687  (std::to_string(i)) + ".pdf");
688 
689  }//end profiles
690  }//end debug
691  }//end hie
692  }//end surf
693  } // end plot only
694 }//end all
695 
696 
697 void DrawHieTag(int hie,bool th13, bool vert){
698  if(hie>0){
699  TLatex * ltxNH = new TLatex(0.16,0.89,"#bf{NH}");
700  ltxNH->SetTextColor(kPrimColorNH);
701  ltxNH->SetNDC();
702  ltxNH->SetTextAlign(22);
703  ltxNH->SetTextSize(kBlessedLabelFontSize);
704  if(th13) ltxNH->Draw();
705  else if(!vert) ltxNH->DrawLatex(.85,.20,"#bf{NH}");
706  else ltxNH->DrawLatex(.85,.85,"#bf{NH}");
707  }
708  else{
709  TLatex * ltxIH = new TLatex (0.16,0.46,"#bf{IH}");
710  ltxIH->SetTextColor(kPrimColorIH);
711  ltxIH->SetNDC();
712  ltxIH->SetTextAlign(22);
713  ltxIH->SetTextSize(kBlessedLabelFontSize);
714  if(th13) ltxIH->Draw();
715  else if(!vert) ltxIH->DrawLatex(.85,.20,"#bf{IH}");
716  else ltxIH->DrawLatex(.85,.85,"#bf{IH}");
717  }
718 }
719 
720 
721 void FillGraphs(std::vector<TGraph*> g1,std::vector<TGraph*> g2,std::vector<TGraph*> g3,
722  const Int_t kFillColor1, const Int_t kFillColor2, const Int_t kFillColor3,
723  const Int_t kDarkColor, const TString surftype)
724 {
725  bool isJoint = surftype.Contains("joint");
726  bool isBoth = surftype.Contains("both");
727  bool isRHC = surftype.Contains("RHCOnly");
728  bool isNH = surftype.Contains("NH");
729  bool isSysts = surftype.Contains("systs");
730  bool isFccorr = surftype.Contains("fccorr");
731 
732  if (surftype.Contains("delta_th23")){
733  cout<<"delta"<<endl;
734  if(isJoint && isNH ){
735  cout<<"joint nh"<<endl;
736  if(isRHC){
737  JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f");
738  JoinGraphs(g2[0], g2[1], kWhite)->Draw("f");
739  JoinGraphs(g2[1], g2[2], kWhite)->Draw("f");
740  JoinGraphs(g2[0], g2[1], kFillColor2)->Draw("f same");
741  JoinGraphs(g2[1], g2[2], kFillColor2)->Draw("f same");
742  }
743  else{
744  JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f");
745  JoinGraphs(g3[0], g3[1], kWhite)->Draw("l");
746  JoinGraphs(g2[0], g2[1], kWhite)->Draw("c f");
747  JoinGraphs(g2[0], g2[1], kFillColor2)->Draw("c f");
748  JoinGraphs(g2[0], g2[1], kWhite)->Draw("c");
749  }
750  if(isBoth){
751  for (auto & g:g1) { g->SetFillColor(kWhite); g->DrawClone("c f");}
752  for (auto & g:g1) { g->SetFillColor(kFillColor1); g->Draw("c f");}
753  for (auto & g:g1) { g->Draw("c");}
754  }
755  }
756 
757  if(isJoint && !isNH){
758  if(isRHC) JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f");
759  else{
760  JoinGraphs(g3[0], g3[2], kFillColor3)->Draw("f c");
761  g3[1]->SetFillColor(kFillColor3); g3[1]->Draw("f c");
762  g3[3]->SetFillColor(kFillColor3); g3[3]->Draw("f c");
763  for (auto &g:g3) {
764  g->SetLineWidth(3);
765  g->Draw("c");
766  }
767  }
768  for (auto &g:g2) {
769  g->SetFillColor(kWhite); g->DrawClone("f c");
770  g->SetFillColor(kFillColor2); g->Draw("f c");
771  g->SetLineColor(kDarkColor);
772  g->SetLineWidth(3);
773  g->Draw("c");
774  }
775  for (auto &g:g1) {
776  g->SetFillColor(kWhite); g->DrawClone("f c");
777  g->SetFillColor(kFillColor1); g->Draw("f c");
778  g->SetLineColor(kDarkColor);
779  g->SetLineWidth(3);
780  g->Draw("c");
781  }
782  for(auto & gs:{g3, g2, g1}){
783  for(auto & g:gs) {
784  g->SetLineColor(kDarkColor);
785  g->SetLineWidth(3);
786  g->Draw("c");
787  }
788  }
789  }
790  }
791 
792  if (surftype.Contains("th23_dm32")){
793  for (auto &g:g3) {
794  g->SetFillColor(kFillColor3); g->Draw("cf");
795  }
796  for (auto &g:g2) {
797  g->SetFillColor(kWhite); g->DrawClone("cf");
798  g->SetFillColor(kFillColor2); g->Draw("cf");
799  }
800  for (auto &g:g1) {
801  g->SetFillColor(kWhite); g->DrawClone("cf");
802  g->SetFillColor(kFillColor1); g->Draw("cf");
803  }
804  for (auto & gs:{g3,g2,g1}){
805  for (auto & g:gs) {
806  g->SetLineColor(kDarkColor);
807  g->SetLineWidth(3);
808  g->Draw("c");
809  }
810  }
811 
812  }
813 
814  if (surftype.Contains("th13_delta")) {
815  if(isBoth) {
816  JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f");
817  if(isNH) {
818  JoinGraphs(g2[0], g2[1], kWhite)->Draw("f");
819  JoinGraphs(g2[0], g2[1], kFillColor2)->Draw("f");
820  }
821  else{
822  for (auto & g:g2) { g->SetFillColor(kWhite); g->DrawClone("f");}
823  for (auto & g:g2) { g->SetFillColor(kFillColor2); g->Draw("f");}
824  }
825  for (auto & g:g1) { g->SetFillColor(kWhite); g->DrawClone("f");}
826  for (auto & g:g1) { g->SetFillColor(kFillColor1); g->Draw("f");}
827  }
828  else{
829  //for (auto &g:g3) {
830  // g->SetFillColor(kFillColor3); g->Draw("f");}
831  //JoinGraphs(g3[0], g3[1], kWhite)->Draw("f");
832  JoinGraphs(g3[0], g3[1], kFillColor3)->Draw("f c");
833  JoinGraphs(g2[0], g2[1], kWhite)->Draw("f c");
834  JoinGraphs(g2[0], g2[1], kFillColor2)->Draw("f c");
835  /*for (auto &g:g2) {
836  g->SetFillColor(kWhite); g->DrawClone("f");
837  g->SetFillColor(kFillColor2); g->Draw("f");
838  }*/
839  JoinGraphs(g1[0], g1[1], kWhite)->Draw("f c");
840  JoinGraphs(g1[0], g1[1], kFillColor1)->Draw("f c");
841  }
842 
843  for (auto & gs:{g3,g2,g1}){
844  for (auto & g:gs) {
845  g->SetLineColor(kDarkColor);
846  g->SetLineWidth(3);
847  g->Draw("l");
848  }
849  }
850  }
851 
852 }
853 
854 void DrawReactorConstraint(double mean, double err)
855 {
856  TBox* box = new TBox(mean-err, 0, mean+err, 2);
857  box->SetFillColorAlpha(kGray,0.7);
858 
859  TLine* linLo = new TLine(mean-err, 0, mean-err, 2);
860  linLo->SetLineColor(kGray+2);
861  linLo->SetLineWidth(2);
862 
863  TLine* linHi = new TLine(mean+err, 0, mean+err, 2);
864  linHi->SetLineColor(kGray+2);
865  linHi->SetLineWidth(2);
866 
867  box->Draw();
868  linLo->Draw();
869  linHi->Draw();
870  gPad->RedrawAxis();
871 }
872 
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
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
double GetValue(const osc::IOscCalcAdjustable *osc) const override
const Color_t kCosmicBackgroundColor
Definition: Style.h:41
double th23
Definition: runWimpSim.h:98
double delta
Definition: runWimpSim.h:98
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)
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
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.
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)
TString hists[nhists]
Definition: bdt_com.C:3
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
virtual void SetDmsq32(const T &dmsq32)=0
const double kAna2018RHCPOT
Definition: Exposures.h:208
const char * label
tuple white
Definition: rootlogon.py:71
const XML_Char const XML_Char * data
Definition: expat.h:268
void demoFitContours(bool createFile=false, bool corrSysts=true, TString options="joint_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):
Log-likelihood scan across two parameters.
const XML_Char * s
Definition: expat.h:262
string outfilename
knobs that need extra care
const DummyAnaSyst kAnaCalibrationSyst("Calibration","AbsCalib")
void SetPaletteDeltaCyclicNew()
void DrawHieTag(int hie, bool th13=true, bool vert=false)
void DrawReactorConstraint(double mean, double err)
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)
std::vector< std::pair< TH1D *, double > > GetNumuCosmics2018(const int nq=4, std::string beam="fhc", bool GetFromUPS=false, bool NERSC=false)
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
TH2 * ToTH2(double minchi=-1) const
Definition: ISurface.cxx:161
TLegend * ContourLegend(int hie, bool fillContour, bool fccorr, Int_t kFillColor1, Int_t kFillColor2, Int_t kFillColor3, Int_t kDarkColor, bool bestFit)
Combine multiple component experiments.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const Style_t k2SigmaConfidenceStyle
Definition: Style.h:71
std::pair< TH1D *, double > GetNueCosmics2018(std::string beam, bool GetFromUPS=false, bool NERSC=false)
double dmsq32
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
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
TH1 * DataMCComparison(const Spectrum &data, const Spectrum &mc, EBinType bintype)
Definition: Plots.cxx:35
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
Float_t e
Definition: plot.C:35
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
static std::unique_ptr< FrequentistSurface > LoadFrom(TDirectory *dir, const std::string &name)
double th13
Definition: runWimpSim.h:98
const std::string outdir
const Color_t k2SigmaConfidenceColorIH
Definition: Style.h:82
FILE * outfile
Definition: dump_event.C:13
Compare a single data spectrum to the MC + cosmics expectation.
surf
Definition: demo5.py:35
unsigned int uint