37 bool corrSysts =
false;
39 bool th13Surf =
false;
40 bool dmsqSurf =
false;
42 bool numuOnly =
false;
45 bool fakedata = true ;
46 bool realdata =
false;
53 if(!onlyNH) onlyIH =
true;
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";
72 TString
outdir =
"/nova/app/users/prabhjot/ReadCAFMakeEventList/development/results/sensitivities/my_macro/";
73 TString
outfilename (outdir +
"hist_contours_2018_" + decomp +
"_" + savefile +
"_2018");
76 TFile *
outfile =
new TFile(outfilename+
".root",
"recreate");
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;
90 if(fakedata && onlyNH)
SetFakeCalc(calc_fake, 0.594, 2.49134
e-3, 0.805 );
91 else if(fakedata && onlyIH)
SetFakeCalc(calc_fake, 0.588, -2.51776
e-3, 1.507 );
95 if(FHCOnly || FHCRHC){
99 if(RHCOnly || FHCRHC){
105 if(FHCOnly || FHCRHC){
107 data.insert(data.end(), nue_data);
109 if(RHCOnly || FHCRHC){
111 data.insert(data.end(), nue_data);
120 if(FHCOnly || FHCRHC){
122 preds.insert(preds.end(), numu_preds.begin(), numu_preds.end());
124 cosmics.insert(cosmics.end(), numu_cosmics.begin(), numu_cosmics.end());
126 if(RHCOnly || FHCRHC){
128 preds.insert(preds.end(), numu_preds.begin(), numu_preds.end());
130 cosmics.insert(cosmics.end(), numu_cosmics.begin(), numu_cosmics.end());
133 if(FHCOnly || FHCRHC){
135 data.insert(data.end(), numu_data.begin(), numu_data.end());
137 if(RHCOnly || FHCRHC){
139 data.insert(data.end(), numu_data.begin(), numu_data.end());
146 for(
int i = 0;
i <
int(preds.size()); ++
i){
155 if(fakedata) thisdata =
GetFakeData(preds[
i], calc_fake, kPOT, cosmics[i].first);
156 if(realdata) thisdata = data[
i];
165 cosmics[
i].first->SetMarkerStyle(34);
167 cosmics[
i].first->Draw(
"hist p same");
168 gPad->Print(outdir +
"debug_predictions.pdf");
178 std::cout <<
"Adding WorldReactorConstraint2017\n";
183 std::cout <<
"Adding Dmsq32ConstraintPDG2017\n";
187 std::cout <<
"Creating Multiexperiment with a total of " << expts.size() <<
" experiments\n\n" <<
std::endl;
193 std::vector<const ISyst*>
systs = {};
197 std::vector <SystShifts> seedShifts = {};
213 double minchi23 = 1E20;
215 for(
double thseed : {0.3, 0.7}){
217 for(
int hie:{-1, 1}){
219 std::cout <<
"\n\nFinding best fit " << (thseed < 0.5 ?
"LO " :
"UO ") << (hie > 0 ?
"NH " :
"IH ") <<
"\n\n";
222 auto thisminchi = fit23.
Fit(calc, auxShifts ,{
227 seedShifts)->EvalMetricVal();
229 minchi23 =
min(minchi23, thisminchi);
243 v.Write(
"overall_minchi");
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");
254 std::cout <<
"Using 30 bins for dmsq surface\n\n";
257 for(
int hie: {-1, +1}){
259 if((onlyNH && hie < 0) || (onlyIH && hie > 0))
continue;
261 std::cout <<
"Starting surface " << (hie > 0?
"NH ":
"IH") <<
"\n\n";
264 if(!th13Surf && ! dmsqSurf){
275 auto surf1 = surf23.
ToTH2(minchi23);
280 surf1->Write((TString)
"delta_th23_"+hieStr);
281 surf23.SaveTo(outfile, (TString)
"surf_delta_th23_"+hieStr);
283 std::vector<TH2*> profhists = surf23.GetProfiledHists();
286 profhists[0]->Write((hieStr+
"_DmSq32").c_str());
287 profhists[1]->Write((hieStr+
"_SinSq2Theta13").c_str());
291 if(!th13Surf && dmsqSurf){
298 &kFitDmSq32, steps, (hie > 0 ? 2
e-3 : -3.0
e-3 ), (hie > 0 ? 3.0
e-3 : -2
e-3),
302 auto surf6=surf23m.
ToTH2(minchi23);
306 surf6->Write((TString)
"th23_dm32_"+(hie > 0 ?
"NH" :
"IH"));
308 surf23m.SaveTo(outfile, (TString)
"surf_th23_dm32_"+(hie > 0 ?
"NH" :
"IH"));
318 &kFitSinSq2Theta13, steps, 0, 0.35,
323 auto surf4 = surf13.
ToTH2(-1);
327 surf4->Write((TString)
"th13_delta_"+(hie > 0 ?
"NH":
"IH"));
328 surf13.SaveTo(outfile, (TString)
"surf_th13_delta_"+(hie > 0 ?
"NH" :
"IH"));
Cuts and Vars for the 2020 FD DiF Study.
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
fvar< T > fabs(const fvar< T > &x)
const Color_t kCosmicBackgroundColor
const FitDmSq32 kFitDmSq32
Simple record of shifts applied to systematic parameters.
const Dmsq32Constraint kDmsq32ConstraintPDG2017(2.45e-3, 0.05e-3, 2.52e-3, 0.05e-3)
static SystShifts Nominal()
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
virtual void SetDmsq32(const T &dmsq32)=0
Representation of a spectrum in any variable, with associated POT.
const double kAna2018RHCPOT
const XML_Char const XML_Char * data
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.
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
virtual T GetDmsq32() const
static float min(const float a, const float b, const float c)
TH2 * ToTH2(double minchi=-1) const
Combine multiple component experiments.
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
void MakeCAFSensitivities_for_FNEX()
const FitSinSq2Theta13 kFitSinSq2Theta13
const double kAna2018FHCPOT
TH1 * DataMCComparison(const Spectrum &data, const Spectrum &mc, EBinType bintype)
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
Compare a single data spectrum to the MC + cosmics expectation.
Perform MINUIT fits in one or two dimensions.