demoStarPlots.C

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////////////
//
// Another method of evaluating the impact of a systematic
//
// Taken largely from macros by Luke, and Diana
// I've trimmed it down a bit for simplicity
// if you want further examples like how to include cosmics
// they can be found here:
// CAFAna/numu/Analysis2018/BoxOpening/make_starplots.C
//
// The idea here is to apply a systematic shift then run a fit that is ignorant
// of the systematics. If we didn't account for it in our fit but it was there
// in data where would it pull the fit.
//
////////////////////////////////////////////////////////////////////////////////

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Fit/Fit.h"
#include "CAFAna/Analysis/Plots.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 "3FlavorAna/Prediction/PredictionSystJoint2018.h"
//#include "CAFAna/Prediction/PredictionSystJointDemo.h"
#include "3FlavorAna/Systs/EnergySysts.h"
#include "3FlavorAna/Systs/NumuSysts.h"
#include "CAFAna/Systs/Systs.h"
#include "CAFAna/Vars/FitVars.h"
#include "OscLib/IOscCalc.h"

#include "TCanvas.h"
#include "TGraph.h"
#include "TLegend.h"
#include "TStyle.h"



#include <fstream>
#include <iostream>

using namespace ana;

// style choices
std::vector<Style_t> styleline = {kSolid, kSolid, kSolid, kSolid, kSolid, kSolid, kSolid, kSolid, kSolid};
std::vector<Color_t> color = {kMagenta, kGreen+1, kOrange+1, kBlue-7, kCyan+1, kRed-4, kBlue+1, kGray+1, kOrange+3};

// Helper function to organize the syst groups
std::map< std::string, std::vector<const ISyst*> > GetSystMap(){
  std::map< std::string, std::vector<const ISyst*> > mapSystVecs;

  std::vector<const ISyst*> muESysts;
  muESysts.push_back(&kMuEScaleSyst2017);
  muESysts.push_back(&kRelMuEScaleSyst2017);
  std::vector<const ISyst*> neutronSysts;
  neutronSysts.push_back(&kNeutronVisEScalePrimariesSyst2018);

  mapSystVecs["MuESysts"] = muESysts;
  mapSystVecs["NeutronSysts"] = neutronSysts;

  return mapSystVecs;
}

// And another helper to reset the calculator
void ResetCalc (osc::IOscCalcAdjustable * calc, double ssth23 = -5, double dmsq32 = -5)
{
  ana::ResetOscCalcToDefault(calc);
  calc->SetTh23(asin(sqrt(ssth23)));
  calc->SetDmsq32(dmsq32);
}

// This should  probably go in a header
// 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 demoStarPlots()
{
  /////////////////////////////////////////////////////////////
  //
  // all the setup stuff
  //
  /////////////////////////////////////////////////////////////
  double pot = kAna2018FHCPOT;
  double livetime = kAna2018FHCLivetime;

  // some fit points
  const double ssth23 = 0.585;
  const double dmsq32 = 2.501e-3;

  // for later use when plotting
  float minx      = 0.35, maxx      = 0.65;
  float minx_zoom = 0.575, maxx_zoom = 0.595;
  float miny      = 2.30, maxy      = 2.70;
  float miny_zoom = 2.44, maxy_zoom = 2.52;

  auto calc = DefaultOscCalc();
  ResetCalc(calc, ssth23, dmsq32);

  std::vector< const ISyst* > systs = {};
  auto systmap = GetSystMap();
  for(auto thissyst : systmap){
    for(unsigned int sysVecIt = 0; sysVecIt < thissyst.second.size(); ++sysVecIt){
      systs.push_back(thissyst.second[sysVecIt]);
    }
  }

  /////////////////////////////////////////////////////////////
  //
  // Loading in predictions and surface if you want it
  //
  /////////////////////////////////////////////////////////////
  std::vector< PredictionSystJoint2018* > predictions;
  for(int quant = 1; quant < 5; ++quant){
    predictions.push_back(new PredictionSystJoint2018(kNuMu, calc, "fhc", quant, systs, "/nova/ana/nu_mu_ana/Ana2018/Predictions/pred_numuconcat_fhc__numu2018.root"));
  }

/*
  // get contour and best fit points
  TFile* fSurfSens = new TFile("/nova/ana/nu_mu_ana/Ana2018/Results/surfprof_fakedata_cosmics_fhc__extrap_stats_2018calc_nh__numu2018.root");
  Surface* surfaceSens = LoadFrom< Surface >(fSurfSens -> GetDirectory("surface")).release();
  fSurfSens -> Close();

  double bestFitDm = surfaceSens-> GetMinY();
  double bestFitSsq = surfaceSens-> GetMinX();
  delete fSurfSens;
*/

  /////////////////////////////////////////////////////////////
  //
  // actually do something
  //
  /////////////////////////////////////////////////////////////
  for(auto const& thisSystList:  systmap){
    std::string NameZoom = "sysfit_starplot_" + (std::string)thisSystList.first + "__zoom";
    TCanvas* canvas_zoom = new TCanvas("canvas_zoom", (NameZoom).c_str());
    TString out_zoom = NameZoom.c_str();

    std::string Name = "sysfit_starplot_"+(std::string)thisSystList.first;
    TCanvas* canvas = new TCanvas("canvas", (Name).c_str());
    TString out = Name.c_str();

    // lots of beautification stuff that you can skip over
    TH2F *axes_zoom = new TH2F("setaxes_zoom","", 30, minx_zoom, maxx_zoom, 50, miny_zoom, maxy_zoom);
    axes_zoom->GetXaxis()->SetTitle("sin^{2}#theta_{23}");
    axes_zoom->GetYaxis()->SetTitle("#Deltam^{2}_{32} (10^{-3} eV^{2})");
    axes_zoom->GetXaxis()->SetTitleOffset(0.95);
    axes_zoom->GetYaxis()->SetTitleOffset(0.85);
    axes_zoom->GetXaxis()->SetTitleSize(0.055);
    axes_zoom->GetYaxis()->SetTitleSize(0.055);
    axes_zoom->GetXaxis()->SetLabelSize(0.04);
    axes_zoom->GetYaxis()->SetLabelSize(0.04);
    axes_zoom->GetXaxis()->CenterTitle();
    axes_zoom->GetYaxis()->CenterTitle();
    axes_zoom->SetTitle("");
    axes_zoom->GetYaxis()->SetDecimals();

    TH2F *axes = new TH2F("setaxes","" ,38 ,minx, maxx, 52, miny, maxy);
    axes->GetXaxis()->SetTitle("sin^{2}#theta_{23}");
    axes->GetYaxis()->SetTitle("#Deltam^{2}_{32} (10^{-3} eV^{2})");
    axes->GetXaxis()->SetTitleOffset(0.95);
    axes->GetYaxis()->SetTitleOffset(0.85);
    axes->GetXaxis()->SetTitleSize(0.055);
    axes->GetYaxis()->SetTitleSize(0.055);
    axes->GetXaxis()->SetLabelSize(0.04);
    axes->GetYaxis()->SetLabelSize(0.04);
    axes->GetXaxis()->CenterTitle();
    axes->GetYaxis()->CenterTitle();
    axes->SetTitle("");
    axes->GetYaxis()->SetDecimals();
    axes->Draw();
    axes->Draw();


    canvas_zoom->cd();
    gPad->SetFillStyle(0);
    gPad->SetTopMargin(0.08);
    gPad->SetRightMargin(0.08);
    gPad->SetBottomMargin(0.12);
    gPad->SetLeftMargin(0.12);
    gStyle->SetTitleOffset(0.95,"x");
    gStyle->SetTitleOffset(0.85,"y");
    axes_zoom->Draw();
//    surfaceSens->DrawContour(Gaussian90Percent2D(*surfaceSens), kSolid, kBlack);
//    surfaceSens->DrawBestFit(kBlack, kFullStar);

    canvas->cd();
    gPad->SetFillStyle(0);
    gPad->SetTopMargin(0.08);
    gPad->SetRightMargin(0.08);
    gPad->SetBottomMargin(0.12);
    gPad->SetLeftMargin(0.12);
    gStyle->SetTitleOffset(0.95,"x");
    gStyle->SetTitleOffset(0.85,"y");
    axes->Draw();
//    surfaceSens->DrawContour(Gaussian90Percent2D(*surfaceSens), kSolid, kBlack);
//    surfaceSens->DrawBestFit(kBlack, kFullStar);
//    DrawLegendBFNull(surfaceSens -> GetMinX(), surfaceSens -> GetMinY(), kBlack, kFullStar);

    TLegend *leg = new TLegend(0.125, 0.45, 0.60, 0.85);
    TLegendEntry *entry = leg->AddEntry("NULL","NOvA NH #nu_{#mu}+#bar{#nu}_{#mu} Systematic","h");
    int graphId = 0;

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


    /////////////////////////////////////////////////////////////
    //
    // now really do something
    //
    /////////////////////////////////////////////////////////////
    for(const ISyst* s: systList){
      // Two graphs one for upper octant and one for lower
      TGraph* gUO = new TGraph;
      TGraph* gLO = new TGraph;
      for(double sigma = -1; sigma < +1.05; sigma += 0.2){
        SystShifts systListShift;
        systListShift.SetShift(s, sigma);
        ResetCalc(calc, ssth23, dmsq32);

        // fake data and experiments
        std::vector<Spectrum> s_fakedata;
        std::vector <const IExperiment*> experiments;
        for(int quant = 0; quant < 4; quant++){
          s_fakedata.push_back( predictions[quant]->PredictSyst(calc, systListShift).FakeData(pot) );
          experiments.push_back(new SingleSampleExperiment(predictions[quant], s_fakedata.back()));
        }
        experiments.push_back(&kSolarConstraintsPDG2017);
        experiments.push_back(WorldReactorConstraint2017());
        MultiExperiment exptSyst(experiments);

        std::vector <const IFitVar*> fitvarsUO = {&kFitSinSqTheta23UO, &kFitDmSq32Scaled};
        std::vector <const IFitVar*> fitvarsLO = {&kFitSinSqTheta23LO, &kFitDmSq32Scaled};

        // Fit the UO and put in graph
        ResetCalc(calc, ssth23, dmsq32);
        Fitter fitter(&exptSyst, fitvarsUO/*, {}, Fitter::kCareful*/);
        fitter.Fit(calc, Fitter::kQuiet);
        gUO->SetPoint(gUO->GetN(), kFitSinSqTheta23UO.GetValue(calc), kFitDmSq32Scaled.GetValue(calc));

        // Fit the LO and put in graph
        ResetCalc(calc, ssth23, dmsq32);
        Fitter fitterLO(&exptSyst, fitvarsLO/*, {}, Fitter::kCareful*/);
        fitterLO.Fit(calc, Fitter::kQuiet);
        gLO->SetPoint(gLO->GetN(), kFitSinSqTheta23LO.GetValue(calc), kFitDmSq32Scaled.GetValue(calc));


      }// end of loop through sigmas

      // More beautification
      gUO->SetLineColor(color[graphId]);
      gUO->SetLineStyle(styleline[graphId]);
      gLO->SetLineColor(color[graphId]);
      gLO->SetLineStyle(styleline[graphId]);
      graphId++;

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

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

      canvas-> cd();
      gUO  ->Draw("l");
      gLO->Draw("l");
//      surfaceSens->DrawBestFit(kBlack,kFullStar);

      canvas_zoom-> cd();
      gUO  ->Draw("l");
      gLO->Draw("l");
//      surfaceSens->DrawBestFit(kBlack,kFullStar);

    }// end of loop over individual systs

    // Draw and save
    canvas-> cd();
    leg->Draw();
    canvas->SaveAs((Name+".png").c_str());
    canvas->Close();

    canvas_zoom-> cd();
    leg->Draw();
    canvas_zoom->SaveAs((NameZoom+".png").c_str());
    canvas_zoom->Close();

  }// end of loop over systmap

}