getHists_FNEX.C
Go to the documentation of this file.
1 
17 #include "CAFAna/FC/FCSurface.h"
20 #include "CAFAna/Vars/FitVars.h"
22 
23 using namespace ana;
24 
25 #include "Utilities/rootlogon.C"
26 #include "OscLib/IOscCalc.h"
27 
28 #include "TCanvas.h"
29 #include "TVectorD.h"
30 #include "TBox.h"
31 #include "TH2.h"
32 #include "TColor.h"
33 #include "TGraph.h"
34 #include "TLegend.h"
35 #include "TLine.h"
36 #include "TStyle.h"
37 
38 
40 {
41 
42  // Select the analyis and the hierarchy
43  // nue NH or nue IH
44  // joint NH or joint IH
45  bool nueOnly = false;
46  bool joint = true;
47  bool onlyNH = false;
48  bool onlyIH = true;
49 
50  // Output root file
51  TString outfilename;
52  outfilename = "histograms_forFNEX_CVN_prop_nueOnly_stats_no_rock_muons_realdata_IH.root";
53 
54  // The Input CAF prediction file - Proportional decomposition
55  TFile * infile = new TFile("/nova/ana/nu_e_ana/Ana2017/Predictions/provisional_singles/latest/other/pred_xp_prop_EnergyCVN2D_full_nueconcat_Nominal.root", "read");
56 
57  TFile * outfile = new TFile(outfilename,"recreate");
58 
59  // Write some information about the analysis
60  TString info;
61  info += "2017 analysis";
62  info += "Nue prediction uses CVN and proportional decomposition";
63  info += "ND MC scaled to data POT";
64  info += "Cosmic spectrum is scaled to data livetime \n";
65  info += "FD predicted components are oscillated and scaled to FD POT \n";
66  info += "Reactor constraint 0.085 \\pm 0.05\n";
67  info += "PDG dmsq32 constraint 2.44 \\pm 0 .06 e-3 and -2.49 \\pm 0.06 e-3 \n";
68  info += "Fit delta, ssth23, dmsq32, ss2th13, no systematics \n";
69  info += "NH IH surfaces wrt best fit point overall \n";
70  info += "POT and Best fit values stored in TVectorD's \n";
71  TNamed readme ("README",info.Data());
72  readme.Write();
73 
74  //ND Data/MC: Nue
75  auto data = Spectrum::LoadFrom(infile, "pred_xp_prop_Nominal/noShift/nue_pred_EnergyCVN2D/extrap/EEextrap/Decomp/data_comp");
76  double ndpot = data->POT();
77 
78  outfile->cd();
79 
80  data->ToTH1(ndpot)->Write("nd_data");
81 
82  outfile->mkdir("nd_mc");
83 
84  outfile->cd("nd_mc");
85 
86  auto nd_mc_nue_comp = Spectrum::LoadFrom(infile, "pred_xp_prop_Nominal/noShift/nue_pred_EnergyCVN2D/extrap/EEextrap/Decomp/nue_comp" );
87  auto nd_mc_antinue_comp = Spectrum::LoadFrom(infile, "pred_xp_prop_Nominal/noShift/nue_pred_EnergyCVN2D/extrap/EEextrap/Decomp/antinue_comp" );
88  auto nd_mc_numu_comp = Spectrum::LoadFrom(infile, "pred_xp_prop_Nominal/noShift/nue_pred_EnergyCVN2D/extrap/EEextrap/Decomp/numu_comp" );
89  auto nd_mc_antinumu_comp = Spectrum::LoadFrom(infile, "pred_xp_prop_Nominal/noShift/nue_pred_EnergyCVN2D/extrap/EEextrap/Decomp/antinumu_comp" );
90  auto nd_mc_nc_comp = Spectrum::LoadFrom(infile, "pred_xp_prop_Nominal/noShift/nue_pred_EnergyCVN2D/extrap/EEextrap/Decomp/nc_comp" );
91  auto nd_mc_total_comp = Spectrum::LoadFrom(infile, "pred_xp_prop_Nominal/noShift/nue_pred_EnergyCVN2D/extrap/EEextrap/Decomp/total_comp" );
92 
93  TH1D* h_nd_mc_nue = nd_mc_nue_comp ->ToTH1(ndpot);
94  TH1D* h_nd_mc_anue = nd_mc_antinue_comp ->ToTH1(ndpot);
95  TH1D* h_nd_mc_numu = nd_mc_numu_comp ->ToTH1(ndpot);
96  TH1D* h_nd_mc_anumu = nd_mc_antinumu_comp ->ToTH1(ndpot);
97  TH1D* h_nd_mc_nc = nd_mc_nc_comp ->ToTH1(ndpot);
98  TH1D* h_nd_mc_total = (TH1D*)h_nd_mc_nue ->Clone("h_nd_mc_total");
99  h_nd_mc_total ->Add(h_nd_mc_anue);
100  h_nd_mc_total ->Add(h_nd_mc_numu);
101  h_nd_mc_total ->Add(h_nd_mc_anumu);
102  h_nd_mc_total ->Add(h_nd_mc_nc);
103 
104  h_nd_mc_nue ->Write("nue" );
105  h_nd_mc_anue ->Write("anue" );
106  h_nd_mc_numu ->Write("numu" );
107  h_nd_mc_anumu ->Write("anumu" );
108  h_nd_mc_nc ->Write("nc" );
109  h_nd_mc_total ->Write("total" );
110 
111 
112  //FD Data
113  Spectrum* intime = LoadFromFile<Spectrum>("/nova/ana/nu_e_ana/Ana2017/data/data_spectra_fd.root", "All/spect_merged_4bin").release();
114  Spectrum* outtime = LoadFromFile<Spectrum>("/nova/ana/nu_e_ana/Ana2017/Predictions/cosmic/cosmic_spectra_prediction.root", "All/cosm_merged_4bin").release();
115 
116  double outN = outtime->Integral(outtime->Livetime(), 0, kLivetime);
117  double innN = outtime->Integral(intime->Livetime() , 0, kLivetime);
118  double cosmicscale = innN/outN;
119 
120  double fdpot = intime->POT();
121 
122  outfile->cd();
123 
124  intime ->ToTH1(fdpot) ->Write("fd_data");
125  outtime ->ToTH1(intime->Livetime(), kLivetime)->Write("fd_cosmic");
126 
127  //FD MC extrapolated predictions
128  std::string decomp = "prop";
129  std::cout
130  << "Calling GetNuePrediction2017 function in the getHists_FNEX.C CAF macro for - "
131  << decomp
132  << " decompositon"
133  << std::endl;
134  auto pred = GetNuePrediction2017(decomp, DefaultOscCalc(), false);
135  std::cout << decomp ;
136 
137  // Get the experiments
138  const SingleSampleExperiment exptNue(pred, *intime, *outtime, 1./sqrt(outN));
139  const ReactorExperiment th13Constr(0.085, 0.005);
140  MultiExperiment exptPDG ({&exptNue, &th13Constr, &kDmsq32ConstraintPDG2015});
141 
142  IExperiment* exptThis = &exptPDG;
143 
144  // Start the fitting to get the best fit values
146  double minchi23 = 1E20;
148  auto bestcalc = calc->Copy();
149  for(double ssth = 0.4; ssth < 0.6; ssth += 0.05){
150  for(int hie : {+1, -1}){
151 // int hie = -1;
152  ResetOscCalcToDefault(calc);
153  kFitSinSqTheta23.SetValue(calc, ssth);
154  calc->SetDmsq32(hie * fabs(calc->GetDmsq32()));
155  double thisminchi = fit23.Fit(calc)->EvalMetricVal();
156  if(thisminchi < minchi23){
157  minchi23 = thisminchi;
158  bestcalc = calc->Copy();
159  }
160  const double dcp = calc->GetdCP();
161  for(int n = 1; n <= 3; ++n){
162  calc->SetdCP(dcp+n*M_PI/2);
163  thisminchi = fit23.Fit(calc)->EvalMetricVal();
164  if(thisminchi<minchi23){
165  minchi23 = thisminchi;
166  bestcalc = calc->Copy();
167  }
168  }// end of loop over delta cp
169  } // end of loop over hierarchies
170  }// end of loop over sinsqth23 values
171 
172  std::cout << "Found best fit point with LL= " << minchi23 << std::endl;
173 
174  for(int hie: {+1, -1}){
175 // int hie = -1;
176  calc = DefaultOscCalc();
177  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
178 
179  double ymin = 0, ymax = 1;
180  FrequentistSurface surf23(exptThis, calc,
181  &kFitDeltaInPiUnits, 30, 0, 2,
182  &kFitSinSqTheta23, 30, ymin, ymax,
184  auto surf = surf23.ToTH2(minchi23);
185  outfile->cd();
186  surf->Write((TString)"dchisq_"+(hie > 0 ? "NH" : "IH"));
187  } // end of loop over hierarchies
188 
189  //Set desired oscillation parameters
190 // bestcalc->SetdCP(4.52);
191 // bestcalc->SetTh13(0.148);
192 // bestcalc->SetTh23(0.761);
193 // bestcalc->SetDmsq32(2.430e-3);
194 
195  //FD extrapolated MC- all components
196  outfile->mkdir("fd_extrap_oscil_BF");
197  outfile->cd("fd_extrap_oscil_BF");
198  auto fd_mc_nue_app = pred->PredictComponent(bestcalc, Flavors::kNuMuToNuE, Current::kCC, Sign::kNu) ;
199  auto fd_mc_nue_surv = pred->PredictComponent(bestcalc, Flavors::kNuEToNuE, Current::kCC, Sign::kNu) ;
200  auto fd_mc_numu_surv = pred->PredictComponent(bestcalc, Flavors::kNuMuToNuMu, Current::kCC, Sign::kNu) ;
201  auto fd_mc_numu_app = pred->PredictComponent(bestcalc, Flavors::kNuEToNuMu, Current::kCC, Sign::kNu) ;
202  auto fd_mc_anue_app = pred->PredictComponent(bestcalc, Flavors::kNuMuToNuE, Current::kCC, Sign::kAntiNu);
203  auto fd_mc_anue_surv = pred->PredictComponent(bestcalc, Flavors::kNuEToNuE, Current::kCC, Sign::kAntiNu);
204  auto fd_mc_anumu_surv = pred->PredictComponent(bestcalc, Flavors::kNuMuToNuMu, Current::kCC, Sign::kAntiNu);
205  auto fd_mc_anumu_app = pred->PredictComponent(bestcalc, Flavors::kNuEToNuMu, Current::kCC, Sign::kAntiNu);
206  auto fd_mc_tau_cc_all = pred->PredictComponent(bestcalc, Flavors::kAllNuTau, Current::kCC, Sign::kBoth) ;
207  auto fd_mc_nc = pred->PredictComponent(bestcalc, Flavors::kAll, Current::kNC, Sign::kBoth) ;
208 
209  TH1D* h_fd_mc_nue_app = fd_mc_nue_app . ToTH1(fdpot);
210  TH1D* h_fd_mc_nue_surv = fd_mc_nue_surv . ToTH1(fdpot);
211  TH1D* h_fd_mc_numu_surv = fd_mc_numu_surv . ToTH1(fdpot);
212  TH1D* h_fd_mc_numu_app = fd_mc_numu_app . ToTH1(fdpot);
213  TH1D* h_fd_mc_anue_app = fd_mc_anue_app . ToTH1(fdpot);
214  TH1D* h_fd_mc_anue_surv = fd_mc_anue_surv . ToTH1(fdpot);
215  TH1D* h_fd_mc_anumu_surv = fd_mc_anumu_surv . ToTH1(fdpot);
216  TH1D* h_fd_mc_anumu_app = fd_mc_anumu_app . ToTH1(fdpot);
217  TH1D* h_fd_mc_tau_cc_all = fd_mc_tau_cc_all . ToTH1(fdpot);
218  TH1D* h_fd_mc_nc = fd_mc_nc . ToTH1(fdpot);
219 
220  TH1D* h_fd_mc_total = (TH1D*)h_fd_mc_nue_app ->Clone("h_fd_mc_total");
221  h_fd_mc_total ->Add(h_fd_mc_nue_surv );
222  h_fd_mc_total ->Add(h_fd_mc_numu_surv );
223  h_fd_mc_total ->Add(h_fd_mc_numu_app );
224  h_fd_mc_total ->Add(h_fd_mc_anue_app );
225  h_fd_mc_total ->Add(h_fd_mc_anue_surv );
226  h_fd_mc_total ->Add(h_fd_mc_anumu_surv);
227  h_fd_mc_total ->Add(h_fd_mc_anumu_app );
228  h_fd_mc_total ->Add(h_fd_mc_tau_cc_all);
229  h_fd_mc_total ->Add(h_fd_mc_nc );
230 
231  h_fd_mc_nue_app ->Write("nue_app" );
232  h_fd_mc_nue_surv ->Write("nue_surv" );
233  h_fd_mc_numu_surv ->Write("numu_surv" );
234  h_fd_mc_numu_app ->Write("numu_app" );
235  h_fd_mc_anue_app ->Write("anue_app" );
236  h_fd_mc_anue_surv ->Write("anue_surv" );
237  h_fd_mc_anumu_surv ->Write("anumu_surv" );
238  h_fd_mc_anumu_app ->Write("anumu_app" );
239  h_fd_mc_tau_cc_all ->Write("tau_cc_all" );
240  h_fd_mc_nc ->Write("nc" );
241  h_fd_mc_total ->Write("total" );
242 
243  //FD un-extrapolated MC- all components
244  decomp = "noextrap";
245  auto pnxp = GetNuePrediction2017(decomp, DefaultOscCalc(), false);
246  std::cout << decomp ;
247 
248  outfile->mkdir("fd_noextrap_oscil_BF");
249  outfile->cd("fd_noextrap_oscil_BF");
250  auto noextrap_fd_mc_nue_app = pnxp->PredictComponent(bestcalc, Flavors::kNuMuToNuE, Current::kCC, Sign::kNu) ;
251  auto noextrap_fd_mc_nue_surv = pnxp->PredictComponent(bestcalc, Flavors::kNuEToNuE, Current::kCC, Sign::kNu) ;
252  auto noextrap_fd_mc_numu_surv = pnxp->PredictComponent(bestcalc, Flavors::kNuMuToNuMu, Current::kCC, Sign::kNu) ;
253  auto noextrap_fd_mc_numu_app = pnxp->PredictComponent(bestcalc, Flavors::kNuEToNuMu, Current::kCC, Sign::kNu) ;
254  auto noextrap_fd_mc_anue_app = pnxp->PredictComponent(bestcalc, Flavors::kNuMuToNuE, Current::kCC, Sign::kAntiNu);
255  auto noextrap_fd_mc_anue_surv = pnxp->PredictComponent(bestcalc, Flavors::kNuEToNuE, Current::kCC, Sign::kAntiNu);
256  auto noextrap_fd_mc_anumu_surv = pnxp->PredictComponent(bestcalc, Flavors::kNuMuToNuMu, Current::kCC, Sign::kAntiNu);
257  auto noextrap_fd_mc_anumu_app = pnxp->PredictComponent(bestcalc, Flavors::kNuEToNuMu, Current::kCC, Sign::kAntiNu);
258  auto noextrap_fd_mc_tau_cc_all = pnxp->PredictComponent(bestcalc, Flavors::kAllNuTau, Current::kCC, Sign::kBoth) ;
259  auto noextrap_fd_mc_nc = pnxp->PredictComponent(bestcalc, Flavors::kAll, Current::kNC, Sign::kBoth) ;
260 
261  TH1D* noextrap_h_fd_mc_nue_app = noextrap_fd_mc_nue_app . ToTH1(fdpot);
262  TH1D* noextrap_h_fd_mc_nue_surv = noextrap_fd_mc_nue_surv . ToTH1(fdpot);
263  TH1D* noextrap_h_fd_mc_numu_surv = noextrap_fd_mc_numu_surv . ToTH1(fdpot);
264  TH1D* noextrap_h_fd_mc_numu_app = noextrap_fd_mc_numu_app . ToTH1(fdpot);
265  TH1D* noextrap_h_fd_mc_anue_app = noextrap_fd_mc_anue_app . ToTH1(fdpot);
266  TH1D* noextrap_h_fd_mc_anue_surv = noextrap_fd_mc_anue_surv . ToTH1(fdpot);
267  TH1D* noextrap_h_fd_mc_anumu_surv = noextrap_fd_mc_anumu_surv . ToTH1(fdpot);
268  TH1D* noextrap_h_fd_mc_anumu_app = noextrap_fd_mc_anumu_app . ToTH1(fdpot);
269  TH1D* noextrap_h_fd_mc_tau_cc_all = noextrap_fd_mc_tau_cc_all . ToTH1(fdpot);
270  TH1D* noextrap_h_fd_mc_nc = noextrap_fd_mc_nc . ToTH1(fdpot);
271 
272  TH1D* noextrap_h_fd_mc_total = (TH1D*)noextrap_h_fd_mc_nue_app ->Clone("noextrap_h_fd_mc_total");
273  noextrap_h_fd_mc_total ->Add(noextrap_h_fd_mc_nue_surv );
274  noextrap_h_fd_mc_total ->Add(noextrap_h_fd_mc_numu_surv );
275  noextrap_h_fd_mc_total ->Add(noextrap_h_fd_mc_numu_app );
276  noextrap_h_fd_mc_total ->Add(noextrap_h_fd_mc_anue_app );
277  noextrap_h_fd_mc_total ->Add(noextrap_h_fd_mc_anue_surv );
278  noextrap_h_fd_mc_total ->Add(noextrap_h_fd_mc_anumu_surv);
279  noextrap_h_fd_mc_total ->Add(noextrap_h_fd_mc_anumu_app );
280  noextrap_h_fd_mc_total ->Add(noextrap_h_fd_mc_tau_cc_all);
281  noextrap_h_fd_mc_total ->Add(noextrap_h_fd_mc_nc );
282 
283  noextrap_h_fd_mc_nue_app ->Write("nue_app" );
284  noextrap_h_fd_mc_nue_surv ->Write("nue_surv" );
285  noextrap_h_fd_mc_numu_surv ->Write("numu_surv" );
286  noextrap_h_fd_mc_numu_app ->Write("numu_app" );
287  noextrap_h_fd_mc_anue_app ->Write("anue_app" );
288  noextrap_h_fd_mc_anue_surv ->Write("anue_surv" );
289  noextrap_h_fd_mc_anumu_surv ->Write("anumu_surv" );
290  noextrap_h_fd_mc_anumu_app ->Write("anumu_app" );
291  noextrap_h_fd_mc_tau_cc_all ->Write("tau_cc_all" );
292  noextrap_h_fd_mc_nc ->Write("nc" );
293  noextrap_h_fd_mc_total ->Write("total" );
294 
295  // Save the ND and FD POT values
296  TVectorD pots(2);
297  pots[0] = ndpot;
298  pots[1] = fdpot;
299  outfile->cd();
300  pots.Write("pot_nd_fd");
301 
302  // Save the oscillation values
303  TVectorD bests(4);
304  bests[0] = bestcalc->GetdCP();
305  bests[1] = bestcalc->GetTh13();
306  bests[2] = bestcalc->GetTh23();
307  bests[3] = bestcalc->GetDmsq32();
308  bests.Write("best_dcp_th13_th23_dmsq32");
309 
310  // Save the chisq values
311  TVectorD bestchi(1);
312  bestchi[0] = minchi23;
313  bestchi.Write("best_chisq");
314 
315 }//end of main function getHists_FNEX
316 
const Dmsq32Constraint kDmsq32ConstraintPDG2015(2.44e-3, 0.06e-3, 2.49e-3, 0.06e-3)
const XML_Char XML_Encoding * info
Definition: expat.h:530
TCut intime("tave>=217.0 && tave<=227.8")
virtual _IOscCalcAdjustable< T > * Copy() const =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
Antineutrinos-only.
Definition: IPrediction.h:50
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
(&#39; appearance&#39;)
Definition: IPrediction.h:18
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:148
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
int pots
(&#39;beam &#39;)
Definition: IPrediction.h:15
T sqrt(T number)
Definition: d0nt_math.hpp:156
const IPrediction * GetNuePrediction2017(std::string decomp, osc::IOscCalc *calc, bool corrSysts, bool GetFromUPS=false)
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:249
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
#define M_PI
Definition: SbMath.h:34
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const XML_Char const XML_Char * data
Definition: expat.h:268
Double_t ymax
Definition: plot.C:25
Log-likelihood scan across two parameters.
string outfilename
knobs that need extra care
Charged-current interactions.
Definition: IPrediction.h:39
string infile
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
virtual T GetDmsq32() const
Definition: IOscCalc.h:91
Very simple model allowing inclusion of reactor constraints.
void getHists_FNEX()
Definition: getHists_FNEX.C:39
double POT() const
Definition: Spectrum.h:227
Neutrinos-only.
Definition: IPrediction.h:49
virtual T GetdCP() const
Definition: IOscCalc.h:95
OStream cout
Definition: OStream.cxx:6
(&#39; survival&#39;)
Definition: IPrediction.h:19
Combine multiple component experiments.
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
Base class defining interface for experiments.
Definition: IExperiment.h:14
Neutral-current interactions.
Definition: IPrediction.h:40
TCut outtime("(tave>0.0&&tave<200.0) || (tave>250.0&&tave<450.0)")
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Double_t ymin
Definition: plot.C:24
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
All neutrinos, any flavor.
Definition: IPrediction.h:26
(&#39; appearance&#39;)
Definition: IPrediction.h:16
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:230
FILE * outfile
Definition: dump_event.C:13
virtual void SetdCP(const T &dCP)=0
Compare a single data spectrum to the MC + cosmics expectation.
gm Write()
surf
Definition: demo5.py:35
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17
enum BeamMode string