saveContours_systs.C

Go to the documentation of this file.

#ifdef __CINT__
void saveContours_systs(){
  std::cout << "Sorry, you must run in compiled mode" << std::endl;
}
#else

#include "../WorkshopIncludes.h"

void saveContours_systs()
{
  TFile* fPred_addExpt = new TFile("../Predictions/Prediction_epoch3b_someSysts_kCCE.root");
  PredictionInterp* prediction_ugly = LoadFrom<PredictionInterp>(fPred_addExpt->GetDirectory("pred_ugly")).release();
  PredictionInterp* prediction_bad  = LoadFrom<PredictionInterp>(fPred_addExpt->GetDirectory("pred_bad")).release();
  PredictionInterp* prediction_good = LoadFrom<PredictionInterp>(fPred_addExpt->GetDirectory("pred_good")).release();

  PredictionInterp* prediction_std = LoadFrom<PredictionInterp>(fPred_addExpt->GetDirectory("pred_std")).release();

  fPred_addExpt->Close();

  osc::OscCalculatorPMNSOpt calc;
  ana::ResetOscCalcToDefault(&calc); // This seeds to default (in particular max-mixing)

  calc.SetTh23(TMath::ASin(sqrt(0.403))); // Second analysis best fit                             
  calc.SetDmsq32(2.67e-3); // Same

  //std::vector<const ISyst*> systs = getAllNumuSysts(); -> contains all the numu systematics
  std::vector<const ISyst*> someSysts = {&kNumuNCScaleSyst, &kSAMuEnergyScaleSyst, 
                                     &kFDSAMuEnergyScaleSyst, &kSARelHadESyst }; // short subset to run faster

  // std::vector<const IFitVar*> oscpar = getNumuMarginalisedOscParam(); --> This gives all the osc param
  std::vector<const IFitVar*> someOscpar = {&kFitSinSq2Theta13}; // to seeep up, th13

  Spectrum fakePred = prediction_std->Predict(&calc);
  Spectrum fakeData = fakePred.FakeData(kSecondAnaPOT);

  // Split experiments data
  Spectrum fakePred_ugly = prediction_ugly->Predict(&calc);
  Spectrum fakeData_ugly = fakePred_ugly.FakeData(kSecondAnaPOT);

  Spectrum fakePred_bad = prediction_bad->Predict(&calc);
  Spectrum fakeData_bad = fakePred_bad.FakeData(kSecondAnaPOT);

  Spectrum fakePred_good = prediction_good->Predict(&calc);
  Spectrum fakeData_good = fakePred_good.FakeData(kSecondAnaPOT);

  // Experiments with constraints

  SingleSampleExperiment fakeExpt(prediction_std, fakeData);
  std::vector<const IExperiment*> vConstraints;
  vConstraints.push_back(&fakeExpt);
  vConstraints.push_back(new SolarConstraints());
  vConstraints.push_back(WorldReactorConstraint2016());
  MultiExperiment ExptWconstr(vConstraints);


  // Multiple exprts
  SingleSampleExperiment fakeExpt_ugly(prediction_ugly, fakeData_ugly);
  SingleSampleExperiment fakeExpt_bad(prediction_bad, fakeData_bad);
  SingleSampleExperiment fakeExpt_good(prediction_good, fakeData_good);
  /// Combine experiments
  std::vector<const IExperiment*> vExpts;
  vExpts.push_back(&fakeExpt_ugly);
  vExpts.push_back(&fakeExpt_bad);
  vExpts.push_back(&fakeExpt_good);
  vExpts.push_back(new SolarConstraints());
  vExpts.push_back(WorldReactorConstraint2016());

  MultiExperiment multiExpt(vExpts);

  // 2D Fits

  Surface* surface = new Surface( &ExptWconstr, &calc,
  &kFitSinSqTheta23, 40, 0.3, 0.7,
  &kFitDmSq32Scaled, 35, 2.3, 3.0, someOscpar, someSysts);

  Surface* surfaceExtra = new Surface( &multiExpt, &calc,
  &kFitSinSqTheta23, 40, 0.3, 0.7,
  &kFitDmSq32Scaled, 35, 2.3, 3.0, someOscpar, someSysts);

  // Plot
  TCanvas* c = new TCanvas();
  c->SetLeftMargin(0.12);
  c->SetBottomMargin(0.15);

  surface->DrawBestFit(kBlue);
  surface->DrawContour( Gaussian90Percent2D(*surface), kSolid, kBlue);

  surfaceExtra->DrawBestFit(kRed);
  surfaceExtra->DrawContour( Gaussian90Percent2D(*surfaceExtra), kDashed, kRed);

  c->SaveAs("../Contours/ExtrapOscparAndSysts.root");
  c->SaveAs("../Contours/ExtrapOscparAndSysts.pdf");

  c->Clear();

  // Get 1D marginalised plots
  TH1* profileSinsq = Profile(&ExptWconstr, &calc, 
    &kFitSinSqTheta23, 40, 0.3, 0.7, -1,
    {&kFitDmSq32Scaled, &kFitSinSq2Theta13}, someSysts );

  TH1* profileSinsqMulti = Profile(&multiExpt, &calc, 
    &kFitSinSqTheta23, 40, 0.3, 0.7, -1,
    {&kFitDmSq32Scaled, &kFitSinSq2Theta13}, someSysts );

  profileSinsq->SetLineColor(kBlue);
  profileSinsq->Draw("hist");
  profileSinsqMulti->SetLineColor(kRed);
  profileSinsqMulti->SetLineStyle(7);
  profileSinsqMulti->Draw("hist same");

  std::cout << "Max-mix rej 1    : " << sqrt(profileSinsq->Interpolate(0.5)) << std::endl;
  std::cout << "Max-mix rej Multi: " << sqrt(profileSinsqMulti->Interpolate(0.5)) << std::endl;

  c->SaveAs("../Contours/1D_SomeSysts.root");
  c->SaveAs("../Contours/1D_SomeSysts.pdf");
  c->Close();
  c->Close();

} // End of function


#endif