demo_numuOnly.C

Go to the documentation of this file.

#ifdef __CINT__
void demo_numuOnly(const TString options = "demo1")
{
  std::cout << "Sorry, you must run in compiled mode" << std::endl;
}
#else


#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Fit.h"
#include "CAFAna/Analysis/Surface.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/FC/FCSurface.h"
#include "CAFAna/Prediction/IPrediction.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Experiment/ReactorExperiment.h"
#include "CAFAna/Experiment/SolarConstraints.h"
#include "CAFAna/Vars/FitVars.h"
#include "OscLib/IOscCalc.h"
#include "TFile.h"
#include "TCanvas.h"
#include "TH2.h"

using namespace ana;


void PlotNumuPredData (IPrediction*, osc::IOscCalcAdjustable*, double,
                       Spectrum&,
                       TString, bool);

void demo_numuOnly(const TString options = "demo1")
{

  TString predFileName = "pred_tutorial_nue_numu.root";
  TFile * predFile = new TFile(predFileName,"read");
  auto predictionNumu = ana::LoadFrom<IPrediction>(predFile->GetDirectory("prediction_numu")); 
  osc::IOscCalcAdjustable* calc = DefaultOscCalc();
  const double potFD = 6E20;
   
  if(options.Contains("demo1")) {
    
    calc->SetTh23(asin(sqrt(0.4)));
    auto fakeDataNumu = predictionNumu->Predict(calc).FakeData(potFD);
    
    ResetOscCalcToDefault(calc);

    PlotNumuPredData(predictionNumu.get(), calc, potFD, fakeDataNumu, "demo1", false);
  }

  if(options.Contains("demo2")){
    
    for(int potfactor:{1,10,100}){
      double pot = potFD * potfactor;
      //Create different types of fake data
      ResetOscCalcToDefault(calc);
      auto fakeData_max = predictionNumu->Predict(calc).FakeData(pot);       
      auto mockData_max = predictionNumu->Predict(calc).MockData(pot);       
      calc->SetTh23(asin(sqrt(0.4)));
      auto fakeData_nonMax = predictionNumu->Predict(calc).FakeData(pot);       
      
      //reduce copy paste
      std::vector <TString> labels = {"fakeMax","mockMax","fakeNonMax"};
      std::vector <Spectrum> data = {fakeData_max,mockData_max,fakeData_nonMax};
      
      for(int i=0;i<3;++i){
        ResetOscCalcToDefault(calc);
        PlotNumuPredData(predictionNumu.get(), calc, pot, data[i], "demo2_"+labels[i]+"_potX"+std::to_string(potfactor), true);
        auto exptNumu = new SingleSampleExperiment(predictionNumu.get(), data[i]);
        
        std::cout << "Fit: prediction maximal mixing, data " << labels[i] <<  std::endl;
        ResetOscCalcToDefault(calc);
        //We are not fitting any parameters
        //Whatever is in the calculator, is used for the prediction
        Fitter fitter (exptNumu,{},{});    
        fitter.Fit(calc);
      }
    }
  }

  ////////////////////////////////////////////////////////////  
  if(options.Contains("demo3")){
    
    ResetOscCalcToDefault(calc);
    calc->SetTh23(asin(sqrt(0.4)));
    auto fakeData = predictionNumu->Predict(calc).FakeData(potFD);       
    auto exptNumu = new SingleSampleExperiment(predictionNumu.get(), fakeData);

    ResetOscCalcToDefault(calc);
    PlotNumuPredData(predictionNumu.get(), calc, potFD, fakeData, "demo3_before", true);

    Fitter fitter (exptNumu,{&kFitSinSqTheta23},{});
    fitter.Fit(calc);
    //Fit result is printed on screen and saved in calc :)
    PlotNumuPredData(predictionNumu.get(), calc, potFD, fakeData, "demo3_after_th23", true);

   //Crazy alternative: fit everything but th23
    ResetOscCalcToDefault(calc);
    Fitter fitter2 (exptNumu,{&kFitSinSq2Theta13, &kFitDeltaInPiUnits, &kFitDmSq32, &kFitSinSq2Theta12, &kFitDmSq21},{});
    fitter2.Fit(calc);
    PlotNumuPredData(predictionNumu.get(), calc, potFD, fakeData, "demo3_after_other", true);
  }
  
  ////////////////////////////////////////////////////////////
  if(options.Contains("demo4")){
    
    ResetOscCalcToDefault(calc);
    calc->SetTh23(asin(sqrt(0.4)));
    auto fakeData = predictionNumu->Predict(calc).FakeData(potFD);       
    auto exptNumu = new SingleSampleExperiment(predictionNumu.get(), fakeData);

    auto reactor = WorldReactorConstraint2015();
    auto solar = new SolarConstraints();

    auto exptNumuConstraint = new MultiExperiment({exptNumu,reactor,solar});

    ResetOscCalcToDefault(calc);
    PlotNumuPredData(predictionNumu.get(), calc, potFD, fakeData, "demo4_before", true);

    Fitter fitter1 (exptNumuConstraint,{&kFitSinSq2Theta13, &kFitDeltaInPiUnits, &kFitDmSq32, &kFitSinSq2Theta12, &kFitDmSq21},{});
    fitter1.Fit(calc);
    PlotNumuPredData(predictionNumu.get(), calc, potFD, fakeData, "demo4_constraints_noth23", true);

    ResetOscCalcToDefault(calc);
    Fitter fitter2 (exptNumuConstraint,{&kFitSinSqTheta23, &kFitSinSq2Theta13, &kFitDeltaInPiUnits, &kFitDmSq32, &kFitSinSq2Theta12, &kFitDmSq21},{});
    fitter2.Fit(calc);
    PlotNumuPredData(predictionNumu.get(), calc, potFD, fakeData, "demo4_constraints_all", true);

    ResetOscCalcToDefault(calc);
    Fitter fitter3 (exptNumu,{ &kFitSinSqTheta23, &kFitSinSq2Theta13, &kFitDeltaInPiUnits, &kFitDmSq32, &kFitSinSq2Theta12, &kFitDmSq21},{});
    fitter3.Fit(calc);
    PlotNumuPredData(predictionNumu.get(), calc, potFD, fakeData, "demo4_noconstraints_all", true);
  }

  if(options.Contains("demo5")){

    ResetOscCalcToDefault(calc);
    calc->SetTh23(asin(sqrt(0.4)));
    auto fakeData = predictionNumu->Predict(calc).FakeData(potFD);       
    auto exptNumu = new SingleSampleExperiment(predictionNumu.get(), fakeData);

    ResetOscCalcToDefault(calc);

    auto prof_th23 = Profile (exptNumu, calc,  &kFitSinSqTheta23, 30, 0, 1,
                              -1, {});

    prof_th23->SetMaximum(450);
    prof_th23->Draw("hist");
    gPad->Print("plot_demo5_profile_th23.pdf");
    
    ResetOscCalcToDefault(calc);
    auto prof_th23_dmsq =  Profile (exptNumu, calc,  &kFitSinSqTheta23, 30, 0, 1,
                                    -1, {&kFitDmSq32});
    
    prof_th23_dmsq->SetMaximum(450);
    prof_th23_dmsq->Draw("hist");
    gPad->Print("plot_demo5_profile_th23_dmsq.pdf");

    ResetOscCalcToDefault(calc);

    auto prof_dCP = Profile (exptNumu, calc,  &kFitDeltaInPiUnits, 30, 0, 2,
                             -1, {});

    prof_dCP->SetMaximum(20);
    prof_dCP->Draw("hist");
    gPad->Print("plot_demo5_profile_dCP.pdf");
    
    ResetOscCalcToDefault(calc);
    auto prof_dCP_th23_dmsq =  Profile (exptNumu, calc,  &kFitDeltaInPiUnits, 30, 0, 2,
                                        -1, {&kFitDmSq32,&kFitSinSqTheta23});
    
    prof_dCP_th23_dmsq->SetMaximum(20);
    prof_dCP_th23_dmsq->Draw("hist");
    gPad->Print("plot_demo5_profile_dCP_th23_dmsq.pdf");
  }

  if(options.Contains("demo6")){

    ResetOscCalcToDefault(calc);
    calc->SetTh23(asin(sqrt(0.4)));
    auto fakeData = predictionNumu->Predict(calc).FakeData(potFD);       
    auto exptNumu = new SingleSampleExperiment(predictionNumu.get(), fakeData);
    auto reactor = WorldReactorConstraint2015();

    auto exptNumuConstraint = new MultiExperiment({exptNumu,reactor});
    ResetOscCalcToDefault(calc);

    Surface surf (exptNumu, calc, 
                  &kFitSinSqTheta23, 30, 0, 1,
                  &kFitDmSq32Scaled, 30, 2,3.5,
                  {&kFitDeltaInPiUnits, &kFitSinSq2Theta13});
    new TCanvas;
    surf.Draw();
    gPad->Print("plot_demo6_surface.pdf");
    auto margs = surf.GetProfiledHists();
    margs[0]->Draw("colz");
    gPad->Print("plot_demo6_surface_prof_dCP.pdf");

    margs[1]->Draw("colz");
    gPad->Print("plot_demo6_surface_prof_th13.pdf");


    ResetOscCalcToDefault(calc);
    
    Surface surf2 (exptNumuConstraint, calc, 
                   &kFitSinSqTheta23, 30, 0, 1,
                   &kFitDmSq32Scaled, 30, 2,3.5,
                   {&kFitDeltaInPiUnits, &kFitSinSq2Theta13},{},
                   {{&kFitDeltaInPiUnits,{0.,0.5,1.,1.5}}});
    new TCanvas;
    surf2.Draw();
    gPad->Print("plot_demo6_surface2.pdf");
    auto margs2 = surf2.GetProfiledHists();
    margs2[0]->Draw("colz");
    gPad->Print("plot_demo6_surface2_prof_dCP.pdf");

    margs2[1]->Draw("colz");
    gPad->Print("plot_demo6_surface2_prof_th13.pdf");
  }
  if(options.Contains("demo7")){

    ResetOscCalcToDefault(calc);
    calc->SetTh23(asin(sqrt(0.4)));
    auto fakeData = predictionNumu->Predict(calc).FakeData(potFD);       
    auto exptNumu = new SingleSampleExperiment(predictionNumu.get(), fakeData);
    ResetOscCalcToDefault(calc);
    Surface surf (exptNumu, calc, 
                  &kFitSinSqTheta23, 30, 0, 1,
                  &kFitDmSq32Scaled, 30, 2,3.5);
    new TCanvas;
    TH2 * surf_1Sigma = Gaussian68Percent2D(surf);
    TH2 * surf_2Sigma = Gaussian2Sigma2D(surf);
    TH2 * surf_3Sigma = Gaussian3Sigma2D(surf);
    surf.DrawContour(surf_1Sigma, kSolid, kRed);
    surf.DrawContour(surf_2Sigma, kSolid, kViolet);
    surf.DrawContour(surf_3Sigma, kSolid, kBlue);
    gPad->Print("plot_demo7_contours.pdf");


  }


  
}


  
void PlotNumuPredData (IPrediction * predictionNumu, osc::IOscCalcAdjustable* calc, 
                       double potFD,
                       Spectrum & fakeDataNumu, 
                       TString suffix, bool sameCanvas=false){
  new TCanvas;

  auto hpred = predictionNumu->Predict(calc).ToTH1(potFD,kRed);
  auto hfake = fakeDataNumu.ToTH1(potFD);
  
  hpred->SetMaximum(1.2*max(hpred->GetMaximum(),hfake->GetMaximum()));
  hfake->SetMaximum(1.2*max(hpred->GetMaximum(),hfake->GetMaximum()));
  
  hfake->SetMarkerStyle(kFullCircle);
  hfake->SetTitle("Fake Data");

  if(sameCanvas){
    hpred->Draw("hist");
    hfake->Draw("same");
    gPad->Print("plot_"+suffix+"_overlay.pdf");
  }
  else{
    hpred->SetTitle("Prediction");
    hpred->Draw("hist");
    gPad->Print("plot_"+suffix+"_pred.pdf");
    hfake->Draw("p");
    gPad->Print("plot_"+suffix+"_data.pdf");
  }
}




#endif