getFitTemplates.C

Go to the documentation of this file.

#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/ISyst.h"
#include "CAFAna/Core/SystShifts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/TruthCuts.h"
#include "CAFAna/Systs/XSecSystLists.h"
#include "3FlavorAna/Cuts/NueCutsSecondAna.h"
#include "3FlavorAna/Vars/NueVars.h"
#include "3FlavorAna/Vars/NueVarsExtra.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "CAFAna/XSec/GenieMultiverseSyst.h"
#include "CAFAna/XSec/FluxMultiverseSyst.h"

#include "NDAna/nuebarcc_inc/NuebarCCIncVars.h"
#include "NDAna/nuebarcc_inc/NuebarCCIncCuts.h"
#include "NDAna/nuebarcc_inc/NuebarCCIncExtra.h"

using namespace ana;

const Var kcvn ([](const caf::SRProxy* sr){
    return (sr->sel.cvn.nueid);
  });

const Var kcvn2017([] (const caf::SRProxy* sr){
    return (sr->sel.cvn2017.nueid);
  });

const Var kcvnProd3([] (const caf::SRProxy* sr){
    return (sr->sel.cvnProd3Train.nueid);
  });

const Var klid([](const caf::SRProxy* sr){
    return sr->sel.lid.ann;
  });

//sr->vtx.elastic[0].fuzzyk.png[iProng].cvnpart.photonid
const Var kprongCVN([](const caf::SRProxy* sr){
    if (sr->vtx.nelastic < 1) return -1000.f;
    if (sr->vtx.elastic[0].fuzzyk.npng < 1) return -1000.f;
    
    int nProngs = sr->vtx.elastic[0].fuzzyk.npng;
    float maxProngCVNe = -5.;
    for(int iProng = 0; iProng < nProngs; iProng++){
      if(sr->vtx.elastic[0].fuzzyk.png[iProng].cvnpart.electronid >
         maxProngCVNe) {
        maxProngCVNe = 
          sr->vtx.elastic[0].fuzzyk.png[iProng].cvnpart.electronid;
      }
    }
    return maxProngCVNe;
  });

const int kNumPIDVars = 1;
const mHistAxisDef pid_vars[kNumPIDVars] = {
  //{"cvne", {"CVNe", Binning::Simple(40,0,1), kcvn}},
  //{"cvne2017", {"CVNe 2017", pidbins, kcvn2017}},
  //{"cvneProd3", {"CVNe Prod 3", pidbins, kcvnProd3}},
  //{"lid", {"LID", pidbins, klid}},
  {"prongCVN", {"Electron Prong Selector", pidbins, kprongCVN}},
  //{"twoviewCVN", {"Electron Prong Selector 2 View", 
  //              pidbins, ktwoviewCVN}},
  //{"flatfluxCVN", {"Electron Prong Selector Flat Flux",
  //               pidbins, kflatfluxCVN}},
  //{"genielikeCVN", {"Electron Prong Selector GENIE Like", 
  //                pidbins, kgenielikeCVN}},
};



void getFitTemplates()
{
  std::vector<std::string> dataset_labels = {"nominal","lightdown",
                                             "lightup", "ckv", "calibneg",
                                             "calibpos","calibshape"};

  std::vector<std::string> datasets = 
    {
      NominalMC_entire_RHC_prod4,
      DWNLightFiles_full_RHC_prod4,
      UPLightFiles_full_RHC_prod4, 
      CherenkovFiles_full_RHC_prod4,
      DWNCalibFiles_full_RHC_prod4, 
      UPCalibFiles_full_RHC_prod4,
      ShapeCalibFiles_full_RHC_prod4
    };

  //Loaders, Chns 
  std::vector<std::vector<ana::Spectrum*>> spec_double_templates;
  std::vector<std::vector<ana::Spectrum*>> spec_double_xsec_templates;
  std::vector<std::vector<ana::Spectrum*>> spec_double_ppfx_templates;

  std::vector<ana::SpectrumLoaderBase*> loaders;
  for(uint idataset = 0; idataset < datasets.size(); ++idataset){
    SpectrumLoaderBase* loader;
    loader = new SpectrumLoader(datasets[idataset]);
    loader->SetSpillCut(kStandardSpillCuts);
    loaders.push_back(loader);
  }//end loop over loaders

  //Get All the Alternative flux and genie universe weights and shifts
  std::vector<ana::Var> ppfx_univ = ana::GetkPPFXFluxUnivWgt();
  GenieMultiverseParameters verse(100,
                                  ana::getAllXsecSysts_2018());
  std::vector<ana::SystShifts> genie_shifts = verse.GetSystShifts();
  
  for(uint iload = 0; iload < loaders.size(); ++iload){
    std::vector<ana::Spectrum*> vspec;
    
    for(int ichn = 0; ichn < kcNumChns; ichn++){
      vspec.push_back(new Spectrum(*loaders[iload],
                                   analysis_vars[1].axis, // costheta
                                   analysis_vars[2].axis, // electronE
                                   pid_vars[0].axis, // prongcvn
                                   kcPresel && chns[ichn].cut,
                                   kNoShift,
                                   kPPFXFluxCVWgt * kXSecCVWgt2018));
    }//loop over interaction channels

    //loop over all the flux and genie universes
    //there are 100 of them so we can do it together
    if(iload == 0){
      std::vector<ana::Spectrum*> vspec_xsec;
      std::vector<ana::Spectrum*> vspec_flux;

      for(int iuniv = 0; iuniv < 100; iuniv++){
        vspec_xsec.push_back(new Spectrum(*loaders[iload],
                                          analysis_vars[1].axis, //costheta
                                          analysis_vars[2].axis, //electronE
                                          pid_vars[0].axis, // prongcvn
                                          kcPresel,
                                          genie_shifts[iuniv],
                                          kPPFXFluxCVWgt * kXSecCVWgt2018));    
        vspec_flux.push_back(new Spectrum(*loaders[iload],
                                          analysis_vars[1].axis, //costheta
                                          analysis_vars[2].axis, //electronE
                                          pid_vars[0].axis, //prongcnv
                                          kcPresel,
                                          kNoShift,
                                          ppfx_univ[iuniv]* kXSecCVWgt2018));   
      }//loop over universes

      spec_double_xsec_templates.push_back(vspec_xsec);
      spec_double_ppfx_templates.push_back(vspec_flux);
    }//if statement
    
    spec_double_templates.push_back(vspec);
  }//loader


  for(uint iload = 0; iload < loaders.size(); ++iload){
    loaders[iload]->Go();
  }

  TFile* outf = 
    new TFile("fFitTemplates_CovarianceMatrix_newPIDBins_prongcvn.root", 
              "recreate");

  for(uint iload = 0; iload < loaders.size(); ++iload){    
    for(int ichn = 0; ichn < kcNumChns; ichn++){
      char name[50];
      sprintf(name, "%s_%s_%s_%s", dataset_labels[iload].c_str(),
              pid_vars[0].name.c_str(), "double", chns[ichn].name.c_str());

      spec_double_templates[iload][ichn]->SaveTo(outf, name);
    }//loop over interaction channels

    //loop over all the flux and genie universes
    //there are 100 of them so we can do it together
    if(iload == 0){
      for(int iuniv = 0; iuniv < 100; iuniv++){
        char name[50];
        sprintf(name, "%s_%i_%s_%s", "genie", iuniv,
                pid_vars[0].name.c_str(), "double");
        spec_double_xsec_templates[iload][iuniv]->SaveTo(outf, name);   

        sprintf(name, "%s_%i_%s_%s", "ppfx", iuniv,
                pid_vars[0].name.c_str(), "double");
        spec_double_ppfx_templates[iload][iuniv]->SaveTo(outf, name);   
      }//loop over universes
    }//if statement
  }//loader
  outf->Close();
}