getHists_FNEX.C
Go to the documentation of this file.
13 #include "CAFAna/Systs/NumuSysts.h"
15 #include "CAFAna/Vars/FitVars.h"
16 
17 #include "Utilities/rootlogon.C"
19 
20 #include "TVectorD.h"
21 
22 #include <algorithm>
23 
24 using namespace ana;
25 
26 void SaveHistogramsND (TFile*, TFile*);
27 void SaveHistogramsFD (TFile*, TFile*, bool, double&);
28 void SaveOscPrediction (TDirectory*, IPrediction * , osc::IOscCalculator *, double);
29 
30 void getHists_FNEX(const bool statsOnly = true,
31  const bool nueOnly = true)
32 {
33 
34  TString suffix = (nueOnly? "nueOnly":"jointFit");
35  suffix += (statsOnly? "_stats":"_systs");
36 
37  TString outfilename ("histograms_forFNEX_CVN_"+suffix+".root");
38  TFile * outFile = new TFile(outfilename,"recreate");
39 
40  //////////////////////////////////////////////////
41  // Load Nue and Numu experiments
42  //////////////////////////////////////////////////
43 
44  // Created with joint_fit_make_experiment
45  TFile * exptFile = new TFile ("expts_sa_numu_nue_cvn_prop.root","read");
46  if(exptFile->IsZombie()) return;
47 
48  auto exptNue = SingleSampleExperiment::LoadFrom(exptFile->GetDirectory("exptNue"))
49  .release();
50  auto exptNumu = SingleSampleExperiment::LoadFrom(exptFile->GetDirectory("exptNumu"))
51  .release();
52 
53  // Funny story: the nue experiment doesn't save things in the format I expected
54  // So load prediction files as well
55 
56  TString predDir = "/nova/ana/nu_e_ana/SecondAna/Predictions/";
57  TFile * predFileNue = new TFile(predDir+
58  "CombinedPeriods/NueSA_prediction_NominalF.root", "READ");
59  TFile * predFileNumu = new TFile(predDir+
60  "pred2_numu_sa_123b_3c.root", "READ");
61 
62  ////////////////////////////////////////////////////////////
63  // Save some spectra
64  ////////////////////////////////////////////////////////////
65 
66  // Note: these functions expect certain file structure
67  double fdpot;
68  SaveHistogramsND (predFileNue, outFile);
69  SaveHistogramsFD (exptFile, outFile, nueOnly, fdpot);
70 
71  //////////////////////////////////////////////////
72  // Add constraints, make experiments
73  //////////////////////////////////////////////////
74 
75  const ReactorExperiment react(0.085, 0.005);
76  MultiExperiment exptJoint({exptNue,exptNumu,&react});
77  MultiExperiment exptNueOnly({exptNue,&react,&kDmsq32ConstraintPDG2015});
78 
79  MultiExperiment * exptThis = &exptJoint;
80  if(nueOnly) exptThis = &exptNueOnly;
81 
82  ////////////////////////////////////////////////////////////
83  // Systematics
84  ////////////////////////////////////////////////////////////
85 
86  std::vector<const ISyst*> systs;
87 
88  // Some definitions for the joint analysis
89  GenieRwgtTableSyst kNumuReweightMaNCEL(rwgt::fReweightMaNCEL);
90  GenieRwgtTableSyst kNumuReweightMaCCQEshape(rwgt::fReweightMaCCQEshape);
91  GenieRwgtTableSyst kNumuReweightNormCCQE(rwgt::fReweightNormCCQE);
92  GenieRwgtTableSyst kNumuReweightMaNCRES(rwgt::fReweightMaNCRES);
93  GenieRwgtTableSyst kNumuReweightMvNCRES(rwgt::fReweightMvNCRES);
94  GenieRwgtTableSyst kNumuReweightCCQEPauliSupViaKF(rwgt::fReweightCCQEPauliSupViaKF);
95  GenieRwgtTableSyst kNumuReweightMaCCRES(rwgt::fReweightMaCCRES);
96  GenieRwgtTableSyst kNumuReweightMvCCRES(rwgt::fReweightMvCCRES);
97 
98  systs={
109  &kSAMuEnergyScaleSyst, &kFDSAMuEnergyScaleSyst,
110  &kNumuSummedSmallGENIESyst, &kNumuReweightMaNCEL
111  };
112 
113  // Correlations: nue systematics to their numu conterparts
114  // Some systematics only apply to numu
115  exptThis->SetSystCorrelations(0, {
116  {&kBeamAll, &kNueNA49Syst},
118  {&kNumuRelNormSyst, NULL},
119  {&kNumuNCScaleSyst, NULL},
120  {&kSAMuEnergyScaleSyst, NULL},
121  {&kFDSAMuEnergyScaleSyst, NULL},
122  {&kNumuSummedSmallGENIESyst, NULL},
123  {&kNumuReweightMaNCEL , NULL}
124  });
125 
126  // Correlations: numu systematics to their nue counterparts
127  // Some systematics only apply to nue
128  if(!nueOnly )
129  exptThis->SetSystCorrelations(1, {
130  {&kNueMaCCQEShapeSyst, &kNumuReweightMaCCQEshape},
131  {&kNueCCQENormSyst, &kNumuReweightNormCCQE},
132  {&kNuePauliSyst, &kNumuReweightCCQEPauliSupViaKF},
133  {&kNueMaRESSyst, &kNumuReweightMaCCRES},
134  {&kNueMaRESSyst, &kNumuReweightMaNCRES},
135  {&kNueMvRESSyst, &kNumuReweightMvCCRES},
136  {&kNueMvRESSyst, &kNumuReweightMvNCRES},
138  {&kNueRelCalibSyst, &kSARelHadESyst},
139  {&kNueMECSyst, &kMECScaleSyst},
140  {&kNueRPASyst, &kRPACCQESyst},
144  {&kNueMaCOHSyst, NULL},
145  {&kNueAllPiSyst, NULL},
146  {&kNueShiftXSyst, NULL},
147  {&kNueNormSystSig, NULL},
148  {&kNueNormSystBkg, NULL }
149  });
150 
151  if(statsOnly) systs={};
152 
153  ////////////////////////////////////////////////////////////
154  // Fit
155  ////////////////////////////////////////////////////////////
156 
157  // Seeds for oscillation parameters and systshifts
159  SystShifts sysTemp = SystShifts::Nominal();
160 
161  // Find the best fit
162  MinuitFitter fit23 (exptThis,
165  systs);
166 
167  // Seed both hierarchies, both octants, diff values of dCP
168  auto minchi23 = fit23.Fit (calc, sysTemp, {
169  {&kFitDmSq32,{-fabs(calc->GetDmsq32()),fabs(calc->GetDmsq32())}},
170  {&kFitSinSqTheta23,{0.3,0.7}},
171  {&kFitDeltaInPiUnits,{0.,0.5,1.,1.5}}
172  });
173 
174  // Store in the file
175  auto fitDir = outFile->mkdir("Fit");
176  fitDir->cd();
177  TVectorD bests(4);
178  bests[0]= calc->GetdCP();
179  bests[1]= calc->GetTh13();
180  bests[2]= calc->GetTh23();
181  bests[3]= calc->GetDmsq32();
182  bests.Write("best_dcp_th13_th23_dmsq32");
183 
184  TVectorD bestchi(1);
185  bestchi[0] = minchi23;
186  bestchi.Write("best_chisq");
187 
188  // Save FD spectra at these values
189 
191  (predFileNue->GetDirectory("prediction_Nominal/nue_pred_CVN_EnergyPID2D")).release();
193  (predFileNue->GetDirectory("prediction_Nominal/nue_pred_noextrap_CVN_EnergyPID2D"))
194  .release();
195 
196  SaveOscPrediction(outFile->mkdir("FD_Nue_Extrap_Oscil_BF"),xpNue,calc,fdpot);
197  SaveOscPrediction(outFile->mkdir("FD_Nue_NoExtrap_Oscil_BF"),nxpNue,calc,fdpot);
198 
199  if(!nueOnly){
200  auto predNumu1 = PredictionInterp::LoadFrom(predFileNumu->GetDirectory
201  ("numu_pred_p123b")).release();
202  auto predNumu2 = PredictionInterp::LoadFrom(predFileNumu->GetDirectory
203  ("numu_pred_p3c")).release();
204  auto xpNumu = new PredictionCombinePeriods (
207  SaveOscPrediction(outFile->mkdir("FD_Numu_Extrap_Oscil_BF"),xpNumu,calc,fdpot);
208  }
209 
210  // Now create the dCP th23 surfaces
211  for(int hie: {-1, +1}){
212  calc = DefaultOscCalc();
213  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
214 
215  FrequentistSurface surf23(exptThis, calc,
216  &kFitDeltaInPiUnits,30, 0, 2,
217  &kFitSinSqTheta23, 30, (nueOnly?0:0.3), (nueOnly?1:0.7),
218  {&kFitDmSq32,&kFitSinSq2Theta13},systs);
219 
220  auto surf1=surf23.ToTH2(minchi23);
221  fitDir->cd();
222  surf1->Write((TString)"delta_th23_"+(hie>0?"NH":"IH"));
223  surf23.SaveTo(fitDir->mkdir((TString)"surf_delta_th23_"+(hie>0?"NH":"IH")));
224  } //end hie
225 
226  // We're done! good job, people
227  std::cout << "Cleaning up... this might take a while\n";
228  delete exptFile, predFileNue, predFileNumu, outFile;
229  std::cout << "Spectra and surfaces saved to " << outfilename << std::endl;
230 }
231 
232 ////////////////////////////////////////////////////////////
233 
234 void SaveHistogramsND (TFile* inFile, TFile* outFile)
235 {
236  TString nueDir = "prediction_Nominal/nue_pred_CVN_EnergyPID2D/pred0/extrap/";
237  TString nueDecomp = nueDir +"EEextrap/Decomp/";
238  TString numuDecomp = nueDir +"MEextrap/Decomp/";
239 
240  for(auto isNue:{true,false}){
241 
242  TString dirDecomp = (isNue? nueDecomp:numuDecomp);
243  TString outDir = (isNue? "ND_Nue":"ND_Numu");
244  auto od = outFile->mkdir(outDir);
245 
246  auto data = Spectrum::LoadFrom(inFile->GetDirectory(dirDecomp+"data_comp"));
247  auto ndpot = data->POT();
248 
249  od->cd();
250  TVectorD pot(1);
251  pot[0]={ndpot};
252  pot.Write("nd_pot");
253  data->ToTH1(ndpot)->Write("nd_data");
254 
255  for (TString comp:{"nue_comp","antinue_comp",
256  "numu_comp","antinumu_comp",
257  "nc_comp"} ) {
258  auto hsp = Spectrum::LoadFrom(inFile->GetDirectory(dirDecomp+comp))->ToTH1(ndpot);
259  od->cd();
260  hsp->Write(comp);
261 
262  } //end mc componets
263  }//end nue or numu decomp
264 }
265 
266 ////////////////////////////////////////////////////////////
267 
268 void SaveHistogramsFD (TFile* exptFile, TFile* outFile, bool nueOnly, double &fdpot)
269 {
270  auto dataNue = Spectrum::LoadFrom(exptFile->GetDirectory("exptNue/data"));
271  auto cosmNue = (TH1D*) exptFile->Get("exptNue/cosmic");
272 
273  fdpot = dataNue->POT();
274 
275  auto od = outFile->mkdir("FD_Data");
276  od->cd();
277 
278  TVectorD pot(1);
279  pot[0]={fdpot};
280  pot.Write("fd_pot");
281 
282  cosmNue->Write("fd_cosmics_nue");
283  dataNue->ToTH1(fdpot)->Write("fd_data_nue");
284 
285  if(!nueOnly){
286  auto dataNumu = Spectrum::LoadFrom(exptFile->GetDirectory("exptNumu/data"));
287  auto cosmNumu = (TH1D*) exptFile->Get("exptNumu/cosmic");
288  od->cd();
289  cosmNumu->Write("fd_cosmics_numu");
290  dataNumu->ToTH1(fdpot)->Write("fd_data_numu");
291  }
292 
293 }
294 
295 void SaveOscPrediction (TDirectory * dir, IPrediction * pred,
296  osc::IOscCalculator *calc, double fdpot)
297 {
298  dir->cd();
299 
301  .ToTH1(fdpot)->Write("nue_app");
303  .ToTH1(fdpot)->Write("nue_surv");
305  .ToTH1(fdpot)->Write("numu_surv");
307  .ToTH1(fdpot)->Write("numu_app");
308 
310  .ToTH1(fdpot)->Write("anue_app");
312  .ToTH1(fdpot)->Write("anue_surv");
314  .ToTH1(fdpot)->Write("anumu_surv");
316  .ToTH1(fdpot)->Write("anumu_app");
317 
319  .ToTH1(fdpot)->Write("tau_cc_all");
320 
322  .ToTH1(fdpot)->Write("nc");
323 }
const DummyNueSyst kNueAllTransportSyst("nueSAAllTransport","All transport",{-1, 0,+1})
const Dmsq32Constraint kDmsq32ConstraintPDG2015(2.44e-3, 0.06e-3, 2.49e-3, 0.06e-3)
const NumuNCScaleSyst kNumuNCScaleSyst
Definition: NumuSysts.cxx:36
Cuts and Vars for the 2020 FD DiF Study.
Definition: CutFlow_header.h:5
osc::IOscCalculatorAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
Antineutrinos-only.
Definition: IPrediction.h:50
std::vector< SystGroupDef > systs
Definition: syst_header.h:384
const DummyNueSyst kNueCCQENormSyst("nueSANormCCQE","CC QE norm",{-2,-1, 0,+1,+2}, kGenieRWNominal)
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:564
const NumuNormSyst kNumuNormSyst
Definition: NumuSysts.cxx:38
virtual Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
const DummyNueSyst kNueAbsCalibSyst("nueSAShiftXYAbs","Abs. calib",{-1, 0,+1}, kNoTauNominal)
(&#39;beam &#39;)
Definition: IPrediction.h:15
const double kSecondAnaEpoch3dPOT
Definition: Exposures.h:77
const DummyNueSyst kNueAllPiSyst("nueSAallpi","All pi",{-2,-1, 0,+1,+2}, kGenieRWNominal)
const DummyNueSyst kNueBirksCSyst("nueSABirksC","Birks C",{0,+1}, kIsAltSample|kNoRockNominal)
const DummyNueSyst kNueNA49Syst("nueSANA49Nue","NA49",{-1, 0,+1})
const double kSecondAnaEpoch3bPOT
Definition: Exposures.h:75
const DummyNueSyst kNueNormSystBoth("nueSAnormBothSA","Sig+Bkg Norm.",{0,+1})
static std::unique_ptr< PredictionCombinePeriods > LoadFrom(TDirectory *dir)
const DummyNueSyst kNueMvRESSyst("nueSAMvRES","M_{V} RES",{-2,-1, 0,+1,+2}, kGenieRWNominal)
void SaveHistogramsND(TFile *, TFile *)
const NumuRelNormSyst kNumuRelNormSyst
Definition: NumuSysts.cxx:39
const double kSecondAnaPeriod1POT
Definition: Exposures.h:73
static SystShifts Nominal()
Definition: SystShifts.h:34
std::string outDir
const DummyNueSyst kNueNormSystSig("nueSAnormSigSA","Sig. Norm.",{0,+1})
const DummyNueSyst kNueNormSystBkg("nueSAnormBkgSA","Bkg. Norm.",{0,+1})
const XML_Char const XML_Char * data
Definition: expat.h:268
const SABirksSyst kSABirksSyst
const DummyNueSyst kNueMaCCQEShapeSyst("nueSAMaCCQEshape","M_{A} CCQE shape",{-2,-1, 0,+1,+2}, kGenieRWNominal)
const DummyNueSyst kNueMaCOHSyst("nueSAMaCOH","M_{A} COH",{-2,-1, 0,+1,+2}, kGenieRWNominal)
Log-likelihood scan across two parameters.
ifstream inFile
Definition: AnaPlotMaker.h:34
const DummyNueSyst kNueRPASyst("nueSARPACCQE","RPA CCQE",{0,+1}, kGenieRWNominal)
string outfilename
knobs that need extra care
Charged-current interactions.
Definition: IPrediction.h:39
const DummyNueSyst kNueMaRESSyst("nueSAMaRES","M_{A} RES",{-2,-1, 0,+1,+2}, kGenieRWNominal)
const DummyNueSyst kNueRelCalibSyst("nueSAShiftXYRel","Rel. calib",{-1, 0,+1}, kNoTauNominal)
TFile * outFile
Definition: PlotXSec.C:135
const DummyNueSyst kNueMECSyst("nueSAMECScale","MEC scale",{-2,-1, 0,+1,+2}, kGenieRWNominal)
#define pot
static const BeamSyst & kBeamAll
Definition: BeamSysts.h:127
Very simple model allowing inclusion of reactor constraints.
osc::OscCalculatorDumb calc
Neutrinos-only.
Definition: IPrediction.h:49
static std::unique_ptr< PredictionInterp > LoadFrom(TDirectory *dir)
void SaveHistogramsFD(TFile *, TFile *, bool, double &)
OStream cout
Definition: OStream.cxx:6
(&#39; survival&#39;)
Definition: IPrediction.h:19
const DummyNueSyst kNueShiftYSyst("nueSAShiftY","Calib y shift",{0,+1}, kIsAltSample|kNoTauNominal)
const SACalibXYSyst kSACalibXYSyst
TH2 * ToTH2(double minchi=-1) const
Definition: ISurface.cxx:170
Combine multiple component experiments.
const SACalibYFuncSyst kSACalibYFuncSyst
const DummyNueSyst kNuePauliSyst("nueSACCQEPauliSupViaKF","Pauli Sup",{-2,-1, 0,+1,+2}, kGenieRWNominal)
virtual T GetDmsq32() const
TDirectory * dir
Definition: macro.C:5
void SaveOscPrediction(TDirectory *, IPrediction *, osc::IOscCalculator *, double)
const FitSinSqTheta23 kFitSinSqTheta23
Definition: FitVars.cxx:15
const NumuSummedSmallGENIESyst kNumuSummedSmallGENIESyst(FindCAFAnaDir()+"/data/numu/summed_genie_systs_nd.root", FindCAFAnaDir()+"/data/numu/summed_genie_systs_fd.root")
Definition: NumuSysts.h:62
Neutral-current interactions.
Definition: IPrediction.h:40
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
void getHists_FNEX()
Definition: getHists_FNEX.C:40
const double kSecondAnaEpoch3cPOT
Definition: Exposures.h:76
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
Sum MC predictions from different periods scaled according to data POT targets.
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir)
Definition: Spectrum.cxx:1066
All neutrinos, any flavor.
Definition: IPrediction.h:26
(&#39; appearance&#39;)
Definition: IPrediction.h:16
const DummyNueSyst kNueShiftXSyst("nueSAShiftX","Calib x shift",{0,+1}, kIsAltSample|kNoTauNominal)
virtual double Fit(osc::IOscCalculatorAdjustable *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:56
const FitDeltaInPiUnits kFitDeltaInPiUnits
Definition: FitVars.cxx:14
static std::unique_ptr< SingleSampleExperiment > LoadFrom(TDirectory *dir)
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:15