saveContours_complete.C

Go to the documentation of this file.

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

#include "../WorkshopIncludes.h"

void saveContours_complete()
{
  // Period 1
  TFile* fPred_p1 = new TFile("../Predictions/Prediction_period1_systs_kCCE.root");
  PredictionInterp* prediction_ugly_p1 = LoadFrom<PredictionInterp>(fPred_p1->GetDirectory("pred_ugly")).release();
  PredictionInterp* prediction_bad_p1  = LoadFrom<PredictionInterp>(fPred_p1->GetDirectory("pred_bad")).release();
  PredictionInterp* prediction_good_p1 = LoadFrom<PredictionInterp>(fPred_p1->GetDirectory("pred_good")).release();

  PredictionInterp* prediction_std_p1 = LoadFrom<PredictionInterp>(fPred_p1->GetDirectory("pred_std")).release();

  fPred_p1->Close();
  // Period 2
  TFile* fPred_p2 = new TFile("../Predictions/Prediction_period2_systs_kCCE.root");
  PredictionInterp* prediction_ugly_p2 = LoadFrom<PredictionInterp>(fPred_p2->GetDirectory("pred_ugly")).release();
  PredictionInterp* prediction_bad_p2  = LoadFrom<PredictionInterp>(fPred_p2->GetDirectory("pred_bad")).release();
  PredictionInterp* prediction_good_p2 = LoadFrom<PredictionInterp>(fPred_p2->GetDirectory("pred_good")).release();

  PredictionInterp* prediction_std_p2 = LoadFrom<PredictionInterp>(fPred_p2->GetDirectory("pred_std")).release();

  fPred_p2->Close();
  // Epoch 3b
  TFile* fPred_3b = new TFile("../Predictions/Prediction_epoch3b_systs_kCCE.root");
  PredictionInterp* prediction_ugly_3b = LoadFrom<PredictionInterp>(fPred_3b->GetDirectory("pred_ugly")).release();
  PredictionInterp* prediction_bad_3b  = LoadFrom<PredictionInterp>(fPred_3b->GetDirectory("pred_bad")).release();
  PredictionInterp* prediction_good_3b = LoadFrom<PredictionInterp>(fPred_3b->GetDirectory("pred_good")).release();

  PredictionInterp* prediction_std_3b = LoadFrom<PredictionInterp>(fPred_3b->GetDirectory("pred_std")).release();

  fPred_3b->Close();
  // Epoch 3c
  TFile* fPred_3c = new TFile("../Predictions/Prediction_period1_systs_kCCE.root");
  PredictionInterp* prediction_ugly_3c = LoadFrom<PredictionInterp>(fPred_3c->GetDirectory("pred_ugly")).release();
  PredictionInterp* prediction_bad_3c  = LoadFrom<PredictionInterp>(fPred_3c->GetDirectory("pred_bad")).release();
  PredictionInterp* prediction_good_3c = LoadFrom<PredictionInterp>(fPred_3c->GetDirectory("pred_good")).release();

  PredictionInterp* prediction_std_3c = LoadFrom<PredictionInterp>(fPred_3c->GetDirectory("pred_std")).release();

  fPred_3c->Close();

  // Now let's combine the different periods

  PredictionCombinePeriods* prediction_ugly = new PredictionCombinePeriods({ 
       {prediction_ugly_p1, kSecondAnaPeriod1POT},
       {prediction_ugly_p2, kSecondAnaPeriod2POT},
       {prediction_ugly_3b, kSecondAnaEpoch3bPOT},
       {prediction_ugly_3c, kSecondAnaEpoch3cPOT + kSecondAnaEpoch3dPOT}  });

  PredictionCombinePeriods* prediction_bad = new PredictionCombinePeriods({ 
       {prediction_bad_p1, kSecondAnaPeriod1POT},
       {prediction_bad_p2, kSecondAnaPeriod2POT},
       {prediction_bad_3b, kSecondAnaEpoch3bPOT},
       {prediction_bad_3c, kSecondAnaEpoch3cPOT + kSecondAnaEpoch3dPOT}  });

  PredictionCombinePeriods* prediction_good = new PredictionCombinePeriods({ 
       {prediction_good_p1, kSecondAnaPeriod1POT},
       {prediction_good_p2, kSecondAnaPeriod2POT},
       {prediction_good_3b, kSecondAnaEpoch3bPOT},
       {prediction_good_3c, kSecondAnaEpoch3cPOT + kSecondAnaEpoch3dPOT}  });

  PredictionCombinePeriods* prediction_std = new PredictionCombinePeriods({ 
       {prediction_std_p1, kSecondAnaPeriod1POT},
       {prediction_std_p2, kSecondAnaPeriod2POT},
       {prediction_std_3b, kSecondAnaEpoch3bPOT},
       {prediction_std_3c, kSecondAnaEpoch3cPOT + kSecondAnaEpoch3dPOT}  });

  // And we're ready to go

  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();

  std::vector<const IFitVar*> oscpar = getNumuMarginalisedOscParam();
  std::vector<const IFitVar*> margPars = getNumuMarginalisedOscParam();
  margPars.push_back(&kFitDmSq32Scaled);

  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, oscpar, systs);

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

  // 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/Complete.root");
  c->SaveAs("../Contours/Complete.pdf");

  c->Clear();

  // Get 1D marginalised plots
  TH1* profileSinsq = Profile(&ExptWconstr, &calc, 
    &kFitSinSqTheta23, 40, 0.3, 0.7, -1,
    margPars, systs );

  TH1* profileSinsqMulti = Profile(&multiExpt, &calc, 
    &kFitSinSqTheta23, 40, 0.3, 0.7, -1,
    margPars, systs );

  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_Complete.root");
  c->SaveAs("../Contours/1D_Complete.pdf");
  c->Close();
  c->Close();

} // End of function


#endif