starPlot.C

Go to the documentation of this file.

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Fit/MinuitFitter.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Fit/FrequentistSurface.h"
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Progress.h"
#include "CAFAna/Core/ReweightableSpectrum.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/TimingCuts.h"
#include "3FlavorAna/Decomp/NumuDecomp.h"
#include "CAFAna/Experiment/AtmConstraint.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Experiment/ReactorExperiment.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Experiment/SolarConstraints.h"
#include "CAFAna/Extrap/ModularExtrap.h"
#include "CAFAna/Extrap/ModularExtrapComponent.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "CAFAna/Prediction/PredictionInterp.h"
#include "CAFAna/Prediction/PredictionCombinePeriods.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "3FlavorAna/Systs/EnergySysts.h"
#include "3FlavorAna/Systs/NumuSysts.h"
#include "CAFAna/Vars/FitVars.h"
#include "3FlavorAna/Vars/FitVarsNumu.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "CAFAna/Vars/Vars.h"
#include "3FlavorAna/Prediction/PredictionSystNumu2017.h"
#include "3FlavorAna/Systs/JointAna2017Systs.h"
#include "OscLib/OscCalcPMNSOpt.h"

using namespace ana;

#include "TCanvas.h"
#include "TFile.h"
#include "TGraph.h"
#include "TGraphAsymmErrors.h"
#include "TH1.h"
#include "TH2.h"
#include "TLatex.h"
#include "TLegend.h"
#include "TStyle.h"

void Preliminary() // For some reason it doesn't include it automatically
{
  TLatex* prelim = new TLatex(.9, .95, "NOvA Preliminary");
  prelim->SetTextColor(kBlue);
  prelim->SetNDC();
  prelim->SetTextSize(2/30.);
  prelim->SetTextAlign(32);
  prelim->Draw();
}

void RestartCalculator( osc::IOscCalcAdjustable *calc  )
{
  ana::ResetOscCalcToDefault(calc);
  calc->SetDmsq32(2.44006e-3);
  calc->SetTh23(asin(sqrt(0.558965)));  
  calc->SetL(810);
  calc->SetRho(0); // No matter effects
  calc->SetDmsq21(7.59e-5);
  calc->SetTh12(.601);
  calc->SetTh13(.1567);
  calc->SetdCP(0);

}


// These two clamp th23 to the uo / lo
class FitSinSqTheta23UO: public IConstrainedFitVar
{
public:
  virtual double GetValue(const osc::IOscCalcAdjustable* osc) const
  {
    return util::sqr(sin(osc->GetTh23()));
  };
  virtual void SetValue(osc::IOscCalcAdjustable* osc, double val) const
  {
    osc->SetTh23(asin(sqrt(Clamp(val))));
  };
  virtual std::string ShortName() const {return "sinsqth23_UO";}
  virtual std::string LatexName() const {return "sin^{2}#theta_{23} Constrained to Upper Octant";}
 
  // And now set max/min values for upper octant boundaries
  virtual double LowLimit() const {return 0.514;}
  virtual double HighLimit() const {return 1.;}
};
const FitSinSqTheta23UO kFitSinSqTheta23UO = FitSinSqTheta23UO();

class FitSinSqTheta23LO: public IConstrainedFitVar
{
public:
  virtual double GetValue(const osc::IOscCalcAdjustable* osc) const
  {
    return util::sqr(sin(osc->GetTh23()));
  };
  virtual void SetValue(osc::IOscCalcAdjustable* osc, double val) const
  {
    osc->SetTh23(asin(sqrt(Clamp(val))));
  };
  virtual std::string ShortName() const {return "sinsqth23_LO";}
  virtual std::string LatexName() const {return "sin^{2}#theta_{23} Constraint to Lower Octant";}
 
  // And now set max/min values for lower octant boundaries
  virtual double LowLimit() const {return 0.;}
  //virtual double HighLimit() const {return 0.514;}
  virtual double HighLimit() const {return 0.49;}
};
const FitSinSqTheta23LO kFitSinSqTheta23LO = FitSinSqTheta23LO();



void starPlot()
{

  gStyle->SetPadLeftMargin(0.12);
  gStyle->SetTitleOffset(1.15,"x");
  gStyle->SetTitleOffset(0.95,"y");


  osc::OscCalcPMNSOpt calc;
  ana::ResetOscCalcToDefault(&calc);

  // make syst list of the largest systs
  std::vector<const ISyst*> systs_XSec;
  std::vector<const ISyst*> systs_Beam;
  std::vector<const ISyst*> systs_FileSysts;
  std::vector<const ISyst*> systs_Others;
  AddJointAna2017XSecSysts(systs_XSec, kNumuAna2017, false);
  AddJointAna2017BeamSysts(systs_Beam, kNumuAna2017);
  AddJointAna2017FileSysts(systs_FileSysts, kNumuAna2017);
  AddJointAna2017OtherSysts(systs_Others, kNumuAna2017);

  std::map< TString, std::vector<const ISyst*> > mapSystLists;
  mapSystLists["XSec"] = systs_XSec;
  mapSystLists["Beam"] = systs_Beam;
  mapSystLists["FileSysts"] = systs_FileSysts;
  mapSystLists["Others"] = systs_Others;

  ENumu2017ExtrapType extrap = kNuMu;
  PredictionSystNumu2017 prediction_Q1(extrap, &calc, 1);
  PredictionSystNumu2017 prediction_Q2(extrap, &calc, 2);
  PredictionSystNumu2017 prediction_Q3(extrap, &calc, 3);
  PredictionSystNumu2017 prediction_Q4(extrap, &calc, 4);

  // get the contour
  TFile* fSurfSens = new TFile("resultCont__allSysts.root"); // result
  FrequentistSurface* surfaceSens = LoadFrom< FrequentistSurface >(fSurfSens, "surface").release();
  fSurfSens -> Close();
  double bestFitDm = surfaceSens -> fBestFitY;
  double bestFitSsq = surfaceSens -> fBestFitX;
  std::cout<< "Best fit from result contour: ssqth23: " << bestFitSsq 
           << ", dm^2_32: " << bestFitDm <<std::endl;
  delete fSurfSens;


  for(auto const& thisSystList:  mapSystLists){

    std::string NameZoom = "starPlotZoom" + (std::string)thisSystList.first;
    TCanvas* cZoom = new TCanvas("cZoom",(NameZoom).c_str());
    cZoom->SetBottomMargin(0.15);
    TH2* h2d = (TH2*) surfaceSens -> ToTH2();
    for(int binx = 1; binx<h2d->GetNbinsX(); binx++)
      for(int biny = 1; biny<h2d->GetNbinsY(); biny++)
        h2d -> SetBinContent(binx, biny, 0);
    h2d -> GetYaxis()->SetRangeUser(2.4,2.48);
    h2d -> GetXaxis()->SetRangeUser(0.455,0.48);
    h2d -> GetYaxis()->CenterTitle();
    h2d -> GetXaxis()->CenterTitle();
    h2d -> Draw();

    std::string Name = "starPlot" + (std::string)thisSystList.first;
    TCanvas* c = new TCanvas("c",(Name).c_str());
    c->SetBottomMargin(0.15);
    surfaceSens->DrawContour(Gaussian90Percent2D(*surfaceSens), kSolid, kBlack);
    surfaceSens->DrawBestFit(kBlack);


    TLegend *leg = new TLegend(0.15, 0.65, 0.60, 0.85);
    int nOfGraphs = 0;


    std::vector<const ISyst*> systList = thisSystList.second;

    for(const ISyst* s: systList){
      Progress prog("Making star plot for "+s->ShortName());
      TGraph* g = new TGraph;
      TGraph* gLO = new TGraph;
      for(double sigma = -1; sigma < +1.05; sigma += 0.2){

        SystShifts systListShift;
        systListShift.SetShift(s, sigma);
        RestartCalculator(&calc); // Restart &calculator after fitting

        // Systematically shifted fake data to fit to
        Spectrum fakeData_Q1 = prediction_Q1.PredictSyst(&calc, systListShift).FakeData(kAna2017POT);
        Spectrum fakeData_Q2 = prediction_Q2.PredictSyst(&calc, systListShift).FakeData(kAna2017POT);
        Spectrum fakeData_Q3 = prediction_Q3.PredictSyst(&calc, systListShift).FakeData(kAna2017POT);
        Spectrum fakeData_Q4 = prediction_Q4.PredictSyst(&calc, systListShift).FakeData(kAna2017POT);

        SingleSampleExperiment expt1(&prediction_Q1, fakeData_Q1);
        SingleSampleExperiment expt2(&prediction_Q2, fakeData_Q2);
        SingleSampleExperiment expt3(&prediction_Q3, fakeData_Q3);
        SingleSampleExperiment expt4(&prediction_Q4, fakeData_Q4);

        std::vector<const IExperiment*> vConstraints;
        vConstraints.push_back(&expt1);
        vConstraints.push_back(&expt2);
        vConstraints.push_back(&expt3);
        vConstraints.push_back(&expt4);
        vConstraints.push_back(&kSolarConstraintsPDG2017);
        vConstraints.push_back(WorldReactorConstraint2017());
        MultiExperiment exptSyst(vConstraints);

        MinuitFitter fitter(&exptSyst, {&kFitSinSqTheta23UO, &kFitDmSq32Scaled});
        fitter.Fit(&calc, IFitter::kQuiet);

        g->SetPoint(g->GetN(),
                    kFitSinSqTheta23UO.GetValue(&calc),
                    kFitDmSq32Scaled.GetValue(&calc));

        std::cout<< "\nsigma: " << sigma <<std::endl;
        std::cout << kFitSinSqTheta23.GetValue(&calc) << " " << kFitDmSq32Scaled.GetValue(&calc) << std::endl;


        MinuitFitter fitterLO(&exptSyst, {&kFitSinSqTheta23LO, &kFitDmSq32Scaled});
        fitterLO.Fit(&calc, IFitter::kQuiet);

        gLO->SetPoint(gLO->GetN(),
                      kFitSinSqTheta23LO.GetValue(&calc),
                      kFitDmSq32Scaled.GetValue(&calc));

        std::cout << kFitSinSqTheta23LO.GetValue(&calc) << " " << kFitDmSq32Scaled.GetValue(&calc) << std::endl;



        prog.SetProgress((sigma+1)/2.);

      }
      if(nOfGraphs<=12)
        {
          // This would still need some playing with the output to make it look pretty,
          // but at least it's not a bunch of graphs all looking the same
          g->SetLineColor(7-nOfGraphs/2); // Hopefully these are different enough
          g->SetLineStyle(nOfGraphs%2 + 1); // Alternating dashed and normal line

          gLO->SetLineColor(7-nOfGraphs/2); // Hopefully these are different enough
          gLO->SetLineStyle(nOfGraphs%2 + 1); // Alternating dashed and normal line
        }
      else
        {
          g->SetLineColor(nOfGraphs%10); 
          g->SetLineStyle((nOfGraphs-10)%10);

          gLO->SetLineColor(nOfGraphs%10); 
          gLO->SetLineStyle((nOfGraphs-10)%10);
        }

      nOfGraphs++;


      g->SetTitle((s->LatexName()).c_str());
      g->SetFillColor(10);
      g->SetMarkerStyle(0);
      g->SetLineWidth(3);
      gLO->SetTitle((s->LatexName()).c_str());
      gLO->SetFillColor(10);
      gLO->SetMarkerStyle(0);
      gLO->SetLineWidth(3);

      leg->AddEntry(g, (s->LatexName()).c_str(), "l");
      leg -> SetFillColor(0);
      leg -> SetFillStyle(0);

      c -> cd();
      g->Draw("l");
      gLO->Draw("l");
      cZoom -> cd();
      g->Draw("l");
      gLO->Draw("l");

    }

    c -> cd();
    Preliminary();
    leg->Draw();

    c->SaveAs(("plots/"+Name+".root").c_str());
    c->SaveAs(("plots/"+Name+".pdf").c_str());
    c->SaveAs(("plots/"+Name+".png").c_str());
    c->Close();


    cZoom -> cd();
    Preliminary();
    leg->Draw();

    cZoom->SaveAs(("plots/"+NameZoom+".root").c_str());
    cZoom->SaveAs(("plots/"+NameZoom+".pdf").c_str());
    cZoom->SaveAs(("plots/"+NameZoom+".png").c_str());
    cZoom->Close();

  }

} // End of function