MakeCAFSensitivities_for_FNEX.C
Go to the documentation of this file.
1 
2 #include "CAFAna/Core/Spectrum.h"
16 #include "CAFAna/Vars/FitVars.h"
18 
19 #include "CAFAna/Analysis/Plots.h"
20 
21 using namespace ana;
22 
23 #include "Utilities/rootlogon.C"
24 #include "OscLib/IOscCalc.h"
25 
26 #include "TCanvas.h"
27 #include "TVectorD.h"
28 #include "TH2.h"
29 #include "TStyle.h"
30 #include "TFile.h"
31 
33 {
34 
35  // Set the conditions we want sensitivites for
36  std::string decomp = "prop"; // Extrapolation tye
37  bool corrSysts = false; // Systematics
38  bool debug = true; // Data/MC debug plots
39  bool th13Surf = false; // Th13 surface
40  bool dmsqSurf = false; // Dmsq32 surface
41  bool nueOnly = true ; // NuE only analysis
42  bool numuOnly = false; // NuMu only analysis
43  bool joint = false; // NuMu-NuE joint analyses
44  bool onlyNH = true; // Normal hierarchy
45  bool fakedata = true ; // Fake FD Data
46  bool realdata = false; // Real FD Data
47  bool FHCOnly = false; // Neutrino mode only
48  bool RHCOnly = false; // Anti-neutrino mode only
49  bool FHCRHC = true ; // Both neutrino and anti-neutrino modes
50 
51  bool onlyIH = false; // Inverted hierarchy
52 
53  if(!onlyNH) onlyIH = true;
54 
55  // Set the name of the output files
56  TString savefile;
57  if(fakedata && joint && onlyNH) savefile = "sensitivity_numu_nue_joint_NH";
58  if(fakedata && joint && onlyIH) savefile = "sensitivity_numu_nue_joint_IH";
59  if(fakedata && nueOnly && onlyNH) savefile = "sensitivity_nueOnly_NH";
60  if(fakedata && nueOnly && onlyIH) savefile = "sensitivity_nueOnly_IH";
61  if(fakedata && numuOnly && onlyNH) savefile = "sensitivity_numuOnly_NH";
62  if(fakedata && numuOnly && onlyIH) savefile = "sensitivity_numuOnly_IH";
63  if(realdata && joint && onlyNH) savefile = "contours_numu_nue_joint_NH";
64  if(realdata && joint && onlyIH) savefile = "contours_numu_nue_joint_IH";
65  if(realdata && nueOnly && onlyNH) savefile = "contours_nueOnly_NH";
66  if(realdata && nueOnly && onlyIH) savefile = "contours_nueOnly_IH";
67  if(realdata && numuOnly && onlyNH) savefile = "contours_numuOnly_NH";
68  if(realdata && numuOnly && onlyIH) savefile = "contours_numuOnly_IH";
69  //else savefile = "NO_ANALYSIS_MENTIONED";
70 
71  // The output directory
72  TString outdir = "/nova/app/users/prabhjot/ReadCAFMakeEventList/development/results/sensitivities/my_macro/";
73  TString outfilename (outdir + "hist_contours_2018_" + decomp + "_" + savefile + "_2018");
74 
75  //We'll save to this file, make pretty plots later
76  TFile * outfile = new TFile(outfilename+".root", "recreate");
77 
78  //////////////////////////////////////////////////
79  // Load Nue experiments
80  //////////////////////////////////////////////////
81 
82  std::vector <const IPrediction * > preds;
83  std::vector <std::pair <TH1D*, double > > cosmics;
84  std::vector <Spectrum * > data;
85  std::vector <const IExperiment * > expts;
86 
87  auto calc_fake = DefaultOscCalc();
88 
89  // Set the values of oscillation paraeters for fake FD Data
90  if(fakedata && onlyNH) SetFakeCalc(calc_fake, 0.594, 2.49134e-3, 0.805 );
91  else if(fakedata && onlyIH) SetFakeCalc(calc_fake, 0.588, -2.51776e-3, 1.507 );
92 
93  // Predictions for the NuE only or joint analyses
94  if(!numuOnly){
95  if(FHCOnly || FHCRHC){
96  preds .push_back(GetNuePrediction2018(decomp, calc_fake, corrSysts, "fhc", false));
97  cosmics .push_back(GetNueCosmics2018("fhc"));
98  }
99  if(RHCOnly || FHCRHC){
100  preds .push_back(GetNuePrediction2018(decomp, calc_fake, corrSysts, "rhc", false));
101  cosmics .push_back(GetNueCosmics2018("rhc"));
102  }
103 
104  if(realdata){
105  if(FHCOnly || FHCRHC){
106  auto nue_data = GetNueData2018("fhc");
107  data.insert(data.end(), nue_data);
108  }
109  if(RHCOnly || FHCRHC){
110  auto nue_data = GetNueData2018("rhc");
111  data.insert(data.end(), nue_data);
112  }
113  }//end of condition that it is a real data
114  }//end of condition that we have nue or joint fit
115 
116  int nnumu = 4;
117 
118  // Predictions for the NuMu only or joint analyses
119  if(!nueOnly){
120  if(FHCOnly || FHCRHC){
121  auto numu_preds = GetNumuPredictions2018(nnumu, corrSysts, "fhc");
122  preds.insert(preds.end(), numu_preds.begin(), numu_preds.end());
123  auto numu_cosmics = GetNumuCosmics2018(nnumu, "fhc");
124  cosmics.insert(cosmics.end(), numu_cosmics.begin(), numu_cosmics.end());
125  }
126  if(RHCOnly || FHCRHC){
127  auto numu_preds = GetNumuPredictions2018(nnumu, corrSysts, "rhc");
128  preds.insert(preds.end(), numu_preds.begin(), numu_preds.end());
129  auto numu_cosmics = GetNumuCosmics2018(nnumu, "rhc");
130  cosmics.insert(cosmics.end(), numu_cosmics.begin(), numu_cosmics.end());
131  }
132  if(realdata){
133  if(FHCOnly || FHCRHC){
134  auto numu_data = GetNumuData2018(nnumu, "fhc");
135  data.insert(data.end(), numu_data.begin(), numu_data.end());
136  }
137  if(RHCOnly || FHCRHC){
138  auto numu_data = GetNumuData2018(nnumu, "rhc");
139  data.insert(data.end(), numu_data.begin(), numu_data.end());
140  }
141  }//end of condition that it is a real data
142  }//end of condition that it is not nue only case
143 
144 
145  //get the data
146  for(int i = 0; i < int(preds.size()); ++i){
147 
148  // Get the POTs
149  double kPOT = 1.0;
150  if(FHCOnly) kPOT = kAna2018FHCPOT;
151  else if(RHCOnly) kPOT = kAna2018RHCPOT;
152  else if(FHCRHC) kPOT = kAna2018FHCPOT + kAna2018RHCPOT;
153 
154  Spectrum* thisdata;
155  if(fakedata) thisdata = GetFakeData(preds[i], calc_fake, kPOT, cosmics[i].first);
156  if(realdata) thisdata = data[i];
157  expts.push_back(new SingleSampleExperiment(preds[i],
158  *thisdata,
159  cosmics[i].first,
160  cosmics[i].second));
161 
162  if(debug){
163  new TCanvas();
164  DataMCComparison (*thisdata, preds[i], calc_fake);
165  cosmics[i].first->SetMarkerStyle(34);
166  cosmics[i].first->SetMarkerColor(kCosmicBackgroundColor);
167  cosmics[i].first->Draw("hist p same");
168  gPad->Print(outdir + "debug_predictions.pdf");
169  }//end of debugging
170  }//end of loop over predictions
171 
172  //////////////////////////////////////////////////
173  // Add constraints, make experiments
174  //////////////////////////////////////////////////
175 
176  std::cout << "\nCreating multiexperiment\n" << std::endl;
177  if(!th13Surf){
178  std::cout << "Adding WorldReactorConstraint2017\n";
179  expts.push_back(WorldReactorConstraint2017());
180  }//end of condition to add th13 world constrainet
181 
182  if(nueOnly){
183  std::cout << "Adding Dmsq32ConstraintPDG2017\n";
184  expts.push_back(&kDmsq32ConstraintPDG2017);
185  }//end of condition that it is a nue only analysis
186 
187  std::cout << "Creating Multiexperiment with a total of " << expts.size() << " experiments\n\n" << std::endl;
188  auto exptThis = new MultiExperiment(expts);
189 
190  ////////////////////////////////////////////////////////////
191  // Systematics - STATS only case
192  ////////////////////////////////////////////////////////////
193  std::vector<const ISyst*> systs = {}; //No systematics
194 
195  SystShifts auxShifts = SystShifts::Nominal(); //Empty systematics - http://nusoft.fnal.gov/nova/novasoft/doxygen/html/classana_1_1SystShifts.html#a23fd2bdae34aaadf646f422d66513403
196 
197  std::vector <SystShifts> seedShifts = {}; //Empty SystShifts seeds
198 
199  ////////////////////////////////////////////////////////////
200  // Fit
201  ////////////////////////////////////////////////////////////
202 
203  std::cout << "Starting the fit" << std::endl;
204 
206 
207  //Find the best fit
208  std::vector <const IFitVar *> fitvars = {&kFitDeltaInPiUnits, &kFitSinSqTheta23, &kFitDmSq32, &kFitSinSq2Theta13 };
209  if(th13Surf) fitvars = {&kFitDeltaInPiUnits, &kFitSinSqTheta23, &kFitDmSq32 };
210 
211  MinuitFitter fit23(exptThis, fitvars, systs);
212 
213  double minchi23 = 1E20;
214 
215  for(double thseed : {0.3, 0.7}){
216 
217  for(int hie:{-1, 1}){
218 
219  std::cout << "\n\nFinding best fit " << (thseed < 0.5 ? "LO " : "UO ") << (hie > 0 ? "NH " : "IH ") << "\n\n";
220 
221  //find minimum chisq
222  auto thisminchi = fit23.Fit(calc, auxShifts ,{
223  {&kFitDmSq32, {hie*fabs(calc->GetDmsq32())} },
224  {&kFitSinSqTheta23, {thseed}},
225  {&kFitDeltaInPiUnits,{0, 0.5, 1, 1.5}}
226  },
227  seedShifts)->EvalMetricVal();
228 
229  minchi23 = min(minchi23, thisminchi);
230 
231  ResetOscCalcToDefault(calc);
232  auxShifts.ResetToNominal();
233 
234  }//end of loop over NH and IH hierarchies
235  }//end of loop over sinsqth23 values as seed
236 
237  std::cout << "\nFound overall minchi " << minchi23 << std::endl;
238 
239  outfile->cd();
240 
241  TVectorD v(1);
242  v[0] = minchi23;
243  v.Write("overall_minchi");
244 
245  TVectorD bests(4);
246  bests[0]= calc_fake->GetdCP();
247  bests[1]= calc_fake->GetTh13();
248  bests[2]= calc_fake->GetTh23();
249  bests[3]= calc_fake->GetDmsq32();
250  bests.Write("best_dcp_th13_th23_dmsq32");
251 
252  int steps = 30;
253 
254  std::cout << "Using 30 bins for dmsq surface\n\n";
255 
256  //Now create all the surfaces
257  for(int hie: {-1, +1}){
258 
259  if((onlyNH && hie < 0) || (onlyIH && hie > 0)) continue;
260 
261  std::cout << "Starting surface " << (hie > 0? "NH ": "IH") << "\n\n";
262 
263  //Draw surfaces only if you are not making surfaces for th13 and Dmsq32
264  if(!th13Surf && ! dmsqSurf){
265 
266  ResetOscCalcToDefault(calc);
267  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
268 
269  FrequentistSurface surf23(exptThis, calc,
270  &kFitDeltaInPiUnits, steps, 0, 2,
271  &kFitSinSqTheta23, steps, (nueOnly ? 0 : 0.3), (nueOnly ? 1 : 0.7),
272  {&kFitDmSq32, &kFitSinSq2Theta13}, systs,
273  {}, seedShifts);
274 
275  auto surf1 = surf23.ToTH2(minchi23);
276 
277  std::string hieStr = hie > 0 ? "NH" : "IH";
278 
279  outfile->cd();
280  surf1->Write((TString)"delta_th23_"+hieStr);
281  surf23.SaveTo(outfile, (TString)"surf_delta_th23_"+hieStr);
282 
283  std::vector<TH2*> profhists = surf23.GetProfiledHists();
284 
285  //Osc parameters are special case
286  profhists[0]->Write((hieStr+"_DmSq32").c_str());
287  profhists[1]->Write((hieStr+"_SinSq2Theta13").c_str());
288  }//end of condition that we don't want surfaces for th13 and Dmsq32
289 
290  //Draw surfaces if you are making surfaces for dmsq but not for th13
291  if(!th13Surf && dmsqSurf){
292 
293  calc = DefaultOscCalc();
294  calc->SetDmsq32(hie * fabs(calc->GetDmsq32()));
295 
296  FrequentistSurface surf23m(exptThis, calc,
297  &kFitSinSqTheta23, steps, 0.3, 0.7,
298  &kFitDmSq32, steps, (hie > 0 ? 2e-3 : -3.0e-3 ), (hie > 0 ? 3.0e-3 : -2e-3),
300  {{&kFitDeltaInPiUnits, {0.5,1.5}}}, seedShifts);
301 
302  auto surf6=surf23m.ToTH2(minchi23);
303 
304  outfile->cd();
305 
306  surf6->Write((TString)"th23_dm32_"+(hie > 0 ? "NH" : "IH"));
307 
308  surf23m.SaveTo(outfile, (TString)"surf_th23_dm32_"+(hie > 0 ? "NH" : "IH"));
309  }//end end of condition that we are making surfaces for dmsq but not for th13
310 
311  //Draw surfaces for th13 (or draw contours for th13)
312  if(th13Surf){
313 
314  calc = DefaultOscCalc();
315  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
316 
317  FrequentistSurface surf13(exptThis, calc,
318  &kFitSinSq2Theta13, steps, 0, 0.35,
319  &kFitDeltaInPiUnits, steps, 0, 2,
320  {&kFitSinSqTheta23, &kFitDmSq32}, systs,
321  {{&kFitSinSqTheta23, {0.3, 0.7}}}, seedShifts);
322 
323  auto surf4 = surf13.ToTH2(-1);
324 
325  outfile->cd();
326 
327  surf4->Write((TString)"th13_delta_"+(hie > 0 ? "NH": "IH"));
328  surf13.SaveTo(outfile, (TString)"surf_th13_delta_"+(hie > 0 ? "NH" : "IH"));
329  }//end of condition to make surfaces for th13
330 
331  }//end of loop over NH and IH hierarchies
332 
333  std::cout << "\nSurfaces saved to " << outfilename << std::endl;
334 
335  //Close file
336  outfile->Close();
337  std::cout << "Name of the outfile is " << outfilename + ".root" << std::endl;
338 }//end of main function MakeCAFSensitivities_for_FNEX
339 
Spectrum * GetFakeData(const IPrediction *pred, osc::IOscCalc *calc, const double pot, TH1D *cosmics=0)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
const Color_t kCosmicBackgroundColor
Definition: Style.h:41
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
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)
std::vector< Spectrum * > GetNumuData2018(const int nq=4, std::string beam="fhc")
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
virtual void SetDmsq32(const T &dmsq32)=0
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
const double kAna2018RHCPOT
Definition: Exposures.h:208
const XML_Char const XML_Char * data
Definition: expat.h:268
Log-likelihood scan across two parameters.
string outfilename
knobs that need extra care
virtual std::unique_ptr< IFitSummary > Fit(osc::IOscCalcAdjustable *seed, SystShifts &bestSysts=junkShifts, const SeedList &seedPts=SeedList(), const std::vector< SystShifts > &systSeedPts={}, Verbosity verb=kVerbose) const
Master fitting method. Depends on FitHelper and FitHelperSeeded.
Definition: IFitter.cxx:64
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
virtual T GetDmsq32() const
Definition: IOscCalc.h:86
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)
OStream cout
Definition: OStream.cxx:6
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
TH2 * ToTH2(double minchi=-1) const
Definition: ISurface.cxx:161
Combine multiple component experiments.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::pair< TH1D *, double > GetNueCosmics2018(std::string beam, bool GetFromUPS=false, bool NERSC=false)
void ResetToNominal()
Definition: SystShifts.cxx:141
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
void MakeCAFSensitivities_for_FNEX()
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
Float_t e
Definition: plot.C:35
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
Spectrum * GetNueData2018(std::string beam)
const std::string outdir
FILE * outfile
Definition: dump_event.C:13
Compare a single data spectrum to the MC + cosmics expectation.
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17