joint_fit_2017_make_fc_surf.C

Go to the documentation of this file.

#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/Utilities.h"

#include "CAFAna/Fit/MinuitFitter.h"
#include "CAFAna/Analysis/Exposures.h"

#include "CAFAna/Prediction/IPrediction.h"
#include "3FlavorAna/Prediction/PredictionSystNue2017.h"

#include "CAFAna/Vars/FitVars.h"

#include "CAFAna/FC/FCPoint.h"
#include "CAFAna/FC/FCCollection.h"
#include "CAFAna/FC/FCSurface.h"

#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Experiment/ReactorExperiment.h"

#include "3FlavorAna/Ana2017/nue/joint_fit_2017_tools.h"

#include "CAFAna/Analysis/Calcs.h"
#include "OscLib/IOscCalc.h"

#include "TRandom3.h"
#include "TH1.h"

using namespace ana;

// Needed for running on the OSG in a tag prior to r20638
bool SetHierarchy(osc::IOscCalcAdjustable* calc, bool hier)
{
  double dmsq32 = hier ? 2.44e-3 : -2.44e-3;
  calc->SetDmsq32(dmsq32);
}


bool RunFitter(MinuitFitter fit, osc::IOscCalcAdjustable* calc,
               double& bestChi, double& bestX, double& bestY,
               SystShifts seed=kNoShift)
{
  double currentChiSq = fit.Fit(calc, seed, IFitter::kQuiet)->EvalMetricVal();
  if (currentChiSq < bestChi){
    bestChi = currentChiSq;
    bestX = kFitDeltaInPiUnits.GetValue(calc);
    bestY = kFitSinSqTheta23.  GetValue(calc);
    return true;
  }
  return false;
}


void joint_fit_2017_make_fc_surf(int NPts, int bin, bool nh, int N,
                                 std::string plot)
{
  assert(bin >= 0 && "Invalid bin number");

  std::string hierStr = nh ? "NH" : "IH";

  FCCollection fccol;
  osc::IOscCalcAdjustable* calc = DefaultOscCalc();

  // Relevant exposures, to be used later
  double pot = kAna2017POT;
  double livetime = kAna2017Livetime;

  // Grab the systs
  std::vector<const ISyst*> systs = GetJointFitSystematicList(true, true,
                                                              true, true);
  std::vector<const ISyst*> snue  = GetJointFitSystematicList(true, true,
                                                              false,true);
  std::vector<const ISyst*> snumu = GetJointFitSystematicList(true, false,
                                                              true, true);

  // We're setting up the requisite inputs to the nue experiment
  // Nue Prediction
  const IPrediction *p = GetNuePrediction2017("combo",calc,true,true);
  std::pair<TH1D*,double> hcosmicnue = GetNueCosmics2017(true);
  // We need to store cosmic hists as spectra to add into the our fake data
  Spectrum cosmicnue(hcosmicnue.first,pot,livetime);

  // Construct numu inputs
  std::vector<const IPrediction*> ps_numu =
    GetNumuPredictions2017(4,true,true);
  std::vector< std::pair<TH1D*,double> > hcosmics_nm =
    GetNumuCosmics2017(4,true);
  std::vector<Spectrum> cosmics_nm;
  for (std::pair<TH1D*,double> h : hcosmics_nm)
    cosmics_nm.push_back(Spectrum(h.first,pot,livetime));



  // Find the bin center for this job, probably should be arguments
  // dCP / pi
  double minX = 0;
  double maxX = 2;
  if (plot == "ssth23dmsq32"){
    minX = 0.3;
    maxX = 0.7;
  }
  int    NX   = 20;
  double widthX = (maxX-minX)/NX;

  // ss theta23
  double minY = 0.3;
  double maxY = 0.7;
  if (plot == "ssth23dmsq32"){
    minY = nh ? 0.002 : -0.003;
    maxY = nh ? 0.003 : -0.002;
  }
  double NY   = 20;
  double widthY = (maxY-minY)/NY;

  assert(bin < NX*NY && "Invalid bin number");

  int binY = bin/NX;
  int binX = bin - NX*binY;

  double centerX = minX + (binX+0.5)*widthX;
  double centerY = minY + (binY+0.5)*widthY;


  // Grab the FCHelper2017 file
  std::string helperupsname = std::string(getenv("FCHELPERANA2017_LIB_PATH"));
  std::string helpername =
    helperupsname+"/contours_"+plot+"_FCInput_2017.root";
  TFile *fchelp = TFile::Open(helpername.c_str());
  fchelp->cd();

  // Pull out hists and get values at current bin
  double seeddmsq32 = 0;
  if (plot == "deltassth23"){
    TH2* hdmsq = (TH2F*)fchelp->Get((hierStr+"_DmSq32").c_str());
    seeddmsq32 = hdmsq->GetBinContent(binX+1,binY+1); // root bins off-by-one
    delete hdmsq;
  }
  double seeddelta = 0;
  if (plot == "ssth23dmsq32"){
    TH2* hdelta = (TH2F*)fchelp->Get((hierStr+"_DeltaInPiUnits").c_str());
    seeddelta = hdelta->GetBinContent(binX+1,binY+1);
    delete hdelta;
  }
  TH2* hss2th13 = (TH2F*)fchelp->Get((hierStr+"_SinSq2Theta13").c_str());
  double seedss2th13 = hss2th13->GetBinContent(binX+1,binY+1);
  delete hss2th13;
  std::map<std::string,double> seedsyst;
  for (const ISyst* syst : systs){
    TH2* h = (TH2F*)fchelp->Get((hierStr+"_"+syst->ShortName()).c_str());
    if (!h){
      std::cout << "Don't see a prof hist for " << syst->ShortName() << ".  Continuing, but check your ups version." << std::endl;
      continue;
    }
    double seedval = h->GetBinContent(binX+1,binY+1);
    seedsyst.insert(std::pair<std::string,double>(syst->ShortName(),seedval));
  }

  // Need random numbers to sample the bin
  TRandom3 rnd(0);

  for (int i = 0; i < NPts; i++){
    ResetOscCalcToDefault(calc);

    if (100*i % NPts == 0){
      std::cout << i << " / " << NPts << "  " << std::endl;
    }
    // Throw a true X,Y value, randomly accross current bin
    double trueX = -1;
    // Edge bins centered on a physical edge  :(
    while (trueX < minX || trueX > maxX)
      trueX = rnd.Uniform(centerX-widthX/2,centerX+widthX/2);
    double trueY = -1;
    while (trueY < minY || trueY > maxY)
      trueY = rnd.Uniform(centerY-widthY/2,centerY+widthY/2);

    kFitDeltaInPiUnits.SetValue(calc, plot=="deltassth23"?trueX:seeddelta);
    kFitSinSqTheta23.  SetValue(calc, plot=="deltassth23"?trueY:trueX);
    kFitDmSq32.        SetValue(calc, plot=="deltassth23"?seeddmsq32:trueY);
    kFitSinSq2Theta13.SetValue(calc,seedss2th13);

    SystShifts seednue;
    for (const ISyst *syst : snue){
      if (seedsyst.find(syst->ShortName()) == seedsyst.end()) continue;
      seednue.SetShift(syst,seedsyst[syst->ShortName()]);
    }
    SystShifts seednumu;
    for (const ISyst *syst : snumu){
      if (seedsyst.find(syst->ShortName()) == seedsyst.end()) continue;
      seednumu.SetShift(syst,seedsyst[syst->ShortName()]);
    }


    // Make a mock data spectrum for nue
    // Use mock data + cosmic as predictions
    Spectrum fakenue = p->PredictSyst(calc,seednue);
    // Add in cosmic prediction here, make sure it gets Poisson fluctuated
    fakenue += cosmicnue;
    Spectrum mocknue = fakenue.MockData(pot);
    // Save this, for posterity's sake
    int mock_integral = mocknue.Integral(pot);

    // Make a mock data spectrum for numu
    std::vector<Spectrum> mocknumu;
    assert(hcosmics_nm.size()==4 && "Issue with N quantiles in cosmics");
    assert(ps_numu.size()==4 && "Issue with N quantiles in predictions");
    for (int i = 0; i < 4; i++){
      Spectrum fakenm = ps_numu[i]->PredictSyst(calc,seednumu);
      fakenm += cosmics_nm[i];
      mocknumu.push_back(fakenm.MockData(pot));
    }

    // Set up experiment, nova nue + nova numu + 2017 reactor
    // Get true chi-square before fit
    std::vector<const IExperiment*> expts;
    expts.push_back(new SingleSampleExperiment(p, mocknue, hcosmicnue.first,
                                               hcosmicnue.second)); // Nue
    for (int i = 0; i < 4; i++)
      expts.push_back(new SingleSampleExperiment(ps_numu[i], mocknumu[i],
                                                 hcosmics_nm[i].first,
                                                 hcosmics_nm[i].second));
    expts.push_back( WorldReactorConstraint2017() );
    MultiExperiment expt(expts);

    // Correlate your systs
    expt.SetSystCorrelations(0, GetCorrelations(true));
    auto notfornumu = GetCorrelations(false);
    for(int i = 0; i < 4; ++i) expt.SetSystCorrelations(i+1, notfornumu);


    std::vector<const IFitVar*> trueVars;
    if (plot=="deltassth23"){
      trueVars.push_back(&kFitDmSq32);
      trueVars.push_back(&kFitSinSq2Theta13);
    }
    if (plot=="ssth23dmsq32"){
      trueVars.push_back(&kFitDeltaInPiUnits);
      trueVars.push_back(&kFitSinSq2Theta13);
    }

    // Finally, we're get to do some FC
    // Set up fitter for #chi^{2}_{true}
    MinuitFitter fittrue(&expt, trueVars, systs);
    double trueChi = 1e6;
    SystShifts mecq0seed;
    for (int i: {0,1}){
      for (double mecq0: {-1,+1}){
        mecq0seed.SetShift(&kMECq0ShapeSyst2017,mecq0);
        if (plot != "deltassth23") kFitDeltaInPiUnits.SetValue(calc,i);
        trueChi = std::min(trueChi,fittrue.Fit(calc,mecq0seed,IFitter::kQuiet)->EvalMetricVal());
      }
    }

    // Set up for best fit
    double dcp = calc->GetdCP();
    double bestX = kFitDeltaInPiUnits.GetValue(calc);
    double bestY = kFitSinSqTheta23.  GetValue(calc);

    // Set up fitter
    MinuitFitter fit(&expt, {&kFitDeltaInPiUnits, &kFitSinSqTheta23, &kFitDmSq32, &kFitSinSq2Theta13},systs);

    double bestChi = 1e6;
    // Save best fit xy for seeding with #chi^{2}_{best}
    int bestHie = 0; int bestDel = 0; int bestOct = 0;
    // Seed fit at dif delta / oct values to deal with local minima
    for (bool hier : {false, true}){
      for (bool oct : {false, true}){
        for (int i = 0; i <=1; i++){
          for (double mecq0 : {-1,+1}){
            mecq0seed.SetShift(&kMECq0ShapeSyst2017,mecq0);
            SetHierarchy(calc, hier);
            // Test at best fit from fitter, and the value reflected about 0.5
            kFitSinSqTheta23.SetValue(calc,oct?bestY:2*0.5-bestY);
            calc->SetdCP(dcp + i*M_PI);
            bool isBest = RunFitter(fit, calc, bestChi, bestX, bestY,
                                    mecq0seed);
            // Set fit location if this is global best chisquare
            if (isBest){
              bestHie=hier; bestDel=i; bestOct=oct;
            }
          }
        }
      }
    }
    // If trueChi is lower than bestChi, we know bestChi is lower than we found
    if (trueChi < bestChi){
      std::cout << "Warning:" << std::endl;
      std::cout << "Pseudo-experiment bestChi, " << bestChi << " is less than the trueChi, " << trueChi << std::endl;
      bestChi = trueChi;
    }


    // These are all useful numbers, encode them and save as FCPoint seed
    double fitloc = mock_integral*1e3+bestOct*1e2+bestHie*1e1+bestDel;


    // Make pseudoexperiment dst
    FCPoint pt(trueX, trueY, trueChi, bestX, bestY, bestChi, fitloc);
    // Every once in a while, print out values we're getting
    if (100*i % NPts == 0){
      std::cout << i << " / " << NPts << std::endl;
      std::cout << "Thrown parameters:  " << trueX << "  " << trueY << std::endl;
      std::cout << "Best fit parameters:  " << bestX << "  " << bestY << std::endl;
      std::cout << "ChiSqaures:  chitrue = " << trueChi << ", chibest =  " << bestChi << ", deltachisquare = " << trueChi-bestChi << std::endl;
    }
    fccol.AddPoint(pt);
  }

  // Save to file
  std::string binname = std::to_string(bin);
  std::string hier = nh ? "_nh_" : "_ih_";
  std::string number = std::to_string(N);
  fccol.SaveToFile("FCCol_"+plot+"_"+binname+hier+number+".root");
}