CAF_makeCAFSensitivities_for_FNEX.C
Go to the documentation of this file.
1 
2 #include "CAFAna/Core/Loaders.h"
5 #include "CAFAna/Core/Spectrum.h"
7 #include "CAFAna/Cuts/Cuts.h"
10 
11 #include "CAFAna/Analysis/Calcs.h"
15 #include "CAFAna/Analysis/Style.h"
26 #include "CAFAna/FC/FCSurface.h"
29 #include "CAFAna/Vars/FitVars.h"
31 
32 #include "CAFAna/Analysis/Plots.h"
33 
34 using namespace ana;
35 
36 #include "Utilities/rootlogon.C"
37 #include "OscLib/IOscCalc.h"
38 
39 #include "TCanvas.h"
40 #include "TVectorD.h"
41 #include "TBox.h"
42 #include "TH2.h"
43 #include "TColor.h"
44 #include "TGraph.h"
45 #include "TLegend.h"
46 #include "TLine.h"
47 #include "TStyle.h"
48 #include "TFile.h"
49 
50 
52 {
53 
54  // Set the analysis
55  std::string decomp = "prop";
56  bool corrSysts = false;
57  bool debug = true;
58  bool th13Surf = false;
59  bool dmsqSurf = false;
60  bool numuOnly = false;
61  bool nueOnly = false;
62  bool joint = true;
63  bool onlyNH = true;
64  bool onlyIH = false;
65  bool fakedata = false;
66  bool realdata = true;
67 
68  // Set the name of the output files
69  TString savefile;
70  if(fakedata && joint && onlyNH) savefile = "sensitivity_numu_nue_joint_NH";
71  if(fakedata && joint && onlyIH) savefile = "sensitivity_numu_nue_joint_IH";
72  if(fakedata && nueOnly && onlyNH) savefile = "sensitivity_nueOnly_NH";
73  if(fakedata && nueOnly && onlyIH) savefile = "sensitivity_nueOnly_IH";
74  if(fakedata && numuOnly && onlyNH) savefile = "sensitivity_numuOnly_NH";
75  if(fakedata && numuOnly && onlyIH) savefile = "sensitivity_numuOnly_IH";
76  if(realdata && joint && onlyNH) savefile = "contours_numu_nue_joint_NH";
77  if(realdata && joint && onlyIH) savefile = "contours_numu_nue_joint_IH";
78  if(realdata && nueOnly && onlyNH) savefile = "contours_nueOnly_NH";
79  if(realdata && nueOnly && onlyIH) savefile = "contours_nueOnly_IH";
80  if(realdata && numuOnly && onlyNH) savefile = "contours_numuOnly_NH";
81  if(realdata && numuOnly && onlyIH) savefile = "contours_numuOnly_IH";
82  //else savefile = "";
83 
84  // The output directory
85  TString outdir = "/nova/app/users/prabhjot/ReadCAFMakeEventList/";
86  TString outfilename (outdir + "hist_contours_2017_" + decomp + "_" + savefile + "_2017");
87 
88  //We'll save to this file, make pretty plots later
89  TFile * outfile = new TFile(outfilename+".root", "recreate");
90 
91  //////////////////////////////////////////////////
92  // Load Nue experiments
93  //////////////////////////////////////////////////
94 
95  std::vector <const IPrediction * > preds;
96  std::vector <std::pair <TH1D*, double > > cosmics;
97  std::vector <Spectrum * > data;
98  std::vector <const IExperiment * > expts;
99 
100  auto calc_fake = DefaultOscCalc();
101 
102  if(fakedata && onlyNH) SetFakeCalc(calc_fake, 0.477, 2.43e-3, 1.432);
103  else if(fakedata && onlyIH) SetFakeCalc(calc_fake, 0.564422, -2.498e-3, 1.45913);
104 
105  if(!numuOnly){
106  preds .push_back(GetNuePrediction2017(decomp, calc_fake, corrSysts));
107  cosmics .push_back(GetNueCosmics2017());
108  if(realdata){
109  auto nue_data = GetNueData2017();
110  data.insert(data.end(), nue_data);
111  }//end of condition that it is a real data
112  }//end of condition that we have nue or joint fit
113 
114  int nnumu = 4;
115 
116  if(!nueOnly){
117  auto numu_preds = GetNumuPredictions2017(nnumu, corrSysts);
118  preds.insert(preds.end(), numu_preds.begin(), numu_preds.end());
119  auto numu_cosmics = GetNumuCosmics2017(nnumu);
120  cosmics.insert(cosmics.end(), numu_cosmics.begin(), numu_cosmics.end());
121  if(realdata){
122  auto numu_data = GetNumuData2017(nnumu);
123  data.insert(data.end(), numu_data.begin(), numu_data.end());
124  }//end of condition that it is a real data
125  }//end of condition that it is not nue only case
126 
127 
128  //get the data
129  for(int i = 0; i < int(preds.size()); ++i){
130  Spectrum* thisdata;
131  if(fakedata) thisdata = GetFakeData (preds[i], calc_fake, kAna2017POT, cosmics[i].first);
132  if(realdata) thisdata = data[i];
133  expts.push_back(new SingleSampleExperiment(preds[i],
134  *thisdata,
135  cosmics[i].first,
136  cosmics[i].second));
137 
138  if(debug){
139  new TCanvas();
140  DataMCComparison (*thisdata, preds[i], calc_fake);
141  cosmics[i].first->SetMarkerStyle(34);
142  cosmics[i].first->SetMarkerColor(kCosmicBackgroundColor);
143  cosmics[i].first->Draw("hist p same");
144  gPad->Print(outdir + "debug_predictions.pdf");
145  }//end of debugging
146  }//end of loop over predictions
147 
148  //////////////////////////////////////////////////
149  // Add constraints, make experiments
150  //////////////////////////////////////////////////
151 
152  std::cout << "\nCreating multiexperiment\n" << std::endl;
153  if(!th13Surf){
154  std::cout << "Adding WorldReactorConstraint2017\n";
155  expts.push_back(WorldReactorConstraint2017());
156  }//end of condition to add th13 world constrainet
157 
158  if(nueOnly){
159  std::cout << "Adding Dmsq32ConstraintPDG2017\n";
160  expts.push_back(&kDmsq32ConstraintPDG2017);
161  }//end of condition that it is a nue only analysis
162 
163  std::cout << "Creating Multiexperiment with a total of " << expts.size() << " experiments\n\n" << std::endl;
164  auto exptThis = new MultiExperiment(expts);
165 
166  ////////////////////////////////////////////////////////////
167  // Systematics - STATS only case
168  ////////////////////////////////////////////////////////////
169  std::vector<const ISyst*> systs = {}; //No systematics
170 
171  SystShifts auxShifts = SystShifts::Nominal(); //Empty systematics - http://nusoft.fnal.gov/nova/novasoft/doxygen/html/classana_1_1SystShifts.html#a23fd2bdae34aaadf646f422d66513403
172 
173  std::vector <SystShifts> seedShifts = {}; //Empty SystShifts seeds
174 
175  ////////////////////////////////////////////////////////////
176  // Fit
177  ////////////////////////////////////////////////////////////
178 
179  std::cout << "Starting the fit" << std::endl;
180 
182 
183  //Find the best fit
185 
186  double minchi23 = 1E20;
187 
188  for(double thseed : {0.3, 0.7}){
189 
190  for(int hie:{-1, 1}){
191 
192  std::cout << "\n\nFinding best fit " << (thseed < 0.5 ? "LO " : "UO ") << (hie > 0 ? "NH " : "IH ") << "\n\n";
193 
194  //find minimum chisq
195  auto thisminchi = fit23.Fit(calc, auxShifts ,{
196  {&kFitDmSq32, {hie*fabs(calc->GetDmsq32())} },
197  {&kFitSinSqTheta23, {thseed}},
198  {&kFitDeltaInPiUnits,{0.33, 1., 1.66}}
199  },
200  seedShifts);
201 
202  minchi23 = min(minchi23, thisminchi);
203 
204  ResetOscCalcToDefault(calc);
205  auxShifts.ResetToNominal();
206 
207  }//end of loop over NH and IH hierarchies
208  }//end of loop over sinsqth23 values as seed
209 
210  std::cout << "\nFound overall minchi " << minchi23 << std::endl;
211 
212  outfile->cd();
213 
214  TVectorD v(1);
215  v[0] = minchi23;
216  v.Write("overall_minchi");
217 
218  TVectorD bests(4);
219  bests[0]= calc_fake->GetdCP();
220  bests[1]= calc_fake->GetTh13();
221  bests[2]= calc_fake->GetTh23();
222  bests[3]= calc_fake->GetDmsq32();
223  bests.Write("best_dcp_th13_th23_dmsq32");
224 
225  int steps = 30;
226 
227  std::cout << "Using 30 bins for dmsq surface\n\n";
228 
229  //Now create all the surfaces
230  for(int hie: {-1, +1}){
231 
232  if((onlyNH && hie < 0) || (onlyIH && hie > 0)) continue;
233 
234  std::cout << "Starting surface " << (hie > 0? "NH ": "IH") << "\n\n";
235 
236  //Draw surfaces only if you are not making surfaces for th13 and Dmsq32
237  if(!th13Surf && ! dmsqSurf){
238 
239  ResetOscCalcToDefault(calc);
240  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
241 
242  FrequentistSurface surf23(exptThis, calc,
243  &kFitDeltaInPiUnits, steps, 0, 2,
244  &kFitSinSqTheta23, steps, (nueOnly ? 0 : 0.3), (nueOnly ? 1 : 0.7),
245  {&kFitDmSq32, &kFitSinSq2Theta13}, systs,
246  {}, seedShifts);
247 
248  auto surf1 = surf23.ToTH2(minchi23);
249 
250  std::string hieStr = hie > 0 ? "NH" : "IH";
251 
252  outfile->cd();
253  surf1->Write((TString)"delta_th23_"+hieStr);
254  surf23.SaveTo(outfile, (TString)"surf_delta_th23_"+hieStr);
255 
256  std::vector<TH2*> profhists = surf23.GetProfiledHists();
257 
258  //Osc parameters are special case
259  profhists[0]->Write((hieStr+"_DmSq32").c_str());
260  profhists[1]->Write((hieStr+"_SinSq2Theta13").c_str());
261  }//end of condition that we don't want surfaces for th13 and Dmsq32
262 
263  //Draw surfaces if you are making surfaces for dmsq but not for th13
264  if(!th13Surf && dmsqSurf){
265 
266  calc = DefaultOscCalc();
267  calc->SetDmsq32(hie * fabs(calc->GetDmsq32()));
268 
269  FrequentistSurface surf23m(exptThis, calc,
270  &kFitSinSqTheta23, steps, 0.3, 0.7,
271  &kFitDmSq32, steps, (hie > 0 ? 2e-3 : -3.0e-3 ), (hie > 0 ? 3.0e-3 : -2e-3),
273  {{&kFitDeltaInPiUnits, {0.5,1.5}}}, seedShifts);
274 
275  auto surf6=surf23m.ToTH2(minchi23);
276 
277  outfile->cd();
278 
279  surf6->Write((TString)"th23_dm32_"+(hie > 0 ? "NH" : "IH"));
280 
281  surf23m.SaveTo(outfile, (TString)"surf_th23_dm32_"+(hie > 0 ? "NH" : "IH"));
282  }//end end of condition that we are making surfaces for dmsq but not for th13
283 
284  //Draw surfaces for th13 (or draw contours for th13)
285  if(th13Surf){
286 
287  calc = DefaultOscCalc();
288  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
289 
290  FrequentistSurface surf13(exptThis, calc,
291  &kFitSinSq2Theta13, steps, 0, 0.35,
292  &kFitDeltaInPiUnits, steps, 0, 2,
294  {{&kFitSinSqTheta23, {0.3, 0.7}}}, seedShifts);
295 
296  auto surf4 = surf13.ToTH2(-1);
297 
298  outfile->cd();
299 
300  surf4->Write((TString)"th13_delta_"+(hie > 0 ? "NH": "IH"));
301  surf13.SaveTo(outfile, (TString)"surf_th13_delta_"+(hie > 0 ? "NH" : "IH"));
302  }//end of condition to make surfaces for th13
303 
304  }//end of loop over NH and IH hierarchies
305 
306  std::cout << "\nSurfaces saved to " << outfilename << std::endl;
307 
308  //Close file
309  outfile->Close();
310  std::cout << "Name of the outfile is " << outfilename + ".root" << std::endl;
311 }//end of main function CAF_makeCAFSensitivities_for_FNEX
312 
Spectrum * GetFakeData(const IPrediction *pred, osc::IOscCalc *calc, const double pot, TH1D *cosmics=0)
std::vector< std::pair< TH1D *, double > > GetNumuCosmics2017(const int nq=4, bool GetFromUPS=false)
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)
const IPrediction * GetNuePrediction2017(std::string decomp, osc::IOscCalc *calc, bool corrSysts, bool GetFromUPS=false)
const double kAna2017POT
Definition: Exposures.h:174
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 XML_Char const XML_Char * data
Definition: expat.h:268
Log-likelihood scan across two parameters.
string outfilename
knobs that need extra care
std::pair< TH1D *, double > GetNueCosmics2017(bool GetFromUPS=false)
std::vector< const IPrediction * > GetNumuPredictions2017(const int nq=4, bool useSysts=true, bool GetFromUPS=false)
std::vector< Spectrum * > GetNumuData2017(const int nq=4)
Spectrum * GetNueData2017()
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)
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
void ResetToNominal()
Definition: SystShifts.cxx:141
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
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
void CAF_makeCAFSensitivities_for_FNEX()
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