getCrossSectionAnalysis_Spectra.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/XSec/Flux.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", Binning::Simple(40,0,1), kcvnProd3}},
  //{"lid", {"LID", Binning::Simple(40,-0.1,1.1), klid}},
  {"prongCVN", {"Electron Prong Selector", pidbins, kprongCVN}},
  //{"twoviewCVN", {"Electron Prong Selector 2 View", 
  //              Binning::Simple(40,0,1), ktwoviewCVN}},
  //{"flatfluxCVN", {"Electron Prong Selector Flat Flux",
  //               Binning::Simple(40,0,1), kflatfluxCVN}},
  //{"genielikeCVN", {"Electron Prong Selector GENIE Like", 
  //                Binning::Simple(40,0,1), kgenielikeCVN}},
};


const Var kMySignalWgt([](const caf::SRProxy* sr)->float
                         {      
                           if(sr->mc.nnu == 0) return 0.;       
                           assert(sr->mc.nnu == 1);
                           if( (sr->mc.nu[0].iscc) 
                               && (abs(sr->mc.nu[0].pdg) == 12)) return 1.06;
                           else return 1.;                
                         });

const Var kQ2Wgt([](const caf::SRProxy* sr)->float
                    { 
                      if(sr->mc.nnu == 0) return 1;
                      
                      double q2 = sr->mc.nu[0].q2;
                      if(q2 < 0) return 1.0;
                      double wt = 0.8 + (q2*0.08);
                      if (wt <= 0) wt = 1.0;
                      
                      return wt;
                    });

const NuTruthVar kMySignalWgtNT([](const caf::SRNeutrinoProxy* sr)->float
                         {      
                           if( (sr->iscc) 
                               && (abs(sr->pdg) == 12)) return 1.06;
                           else return 1.;                
                         });

const NuTruthVar kQ2WgtNT([](const caf::SRNeutrinoProxy* sr)->float
                    { 
                      double q2 = sr->q2;
                      if(q2 < 0) return 1.0;
                      double wt = 0.8 + (q2*0.08);
                      if (wt <= 0) wt = 1.0;
                      
                      return wt;
                    });

void getCrossSectionAnalysis_Spectra()
{
  std::vector<std::string> dataset_labels = {"nominal","fakedata",
                                             "fakedata_q2", "fakedata_sig"};

  std::vector<std::string> datasets = 
    {
      NominalMC_entire_RHC_prod4, NominalMC_entire_RHC_prod4 //, NominalFakeData_2_full_RHC
    };

  std::vector<Var> weights = { (kPPFXFluxCVWgt * kXSecCVWgt2018),
                               (kPPFXFluxCVWgt * kXSecCVWgt2018),
                               (kPPFXFluxCVWgt * kXSecCVWgt2018 * kQ2Wgt),
                               (kPPFXFluxCVWgt * kXSecCVWgt2018 * kMySignalWgt)
  };

  std::vector<NuTruthVar> weights_nt = { (kPPFXFluxCVWgtST * 
                                          kXSecCVWgt2018_smallerDISScale_NT),
                                         (kPPFXFluxCVWgtST * 
                                          kXSecCVWgt2018_smallerDISScale_NT),
                                         (kPPFXFluxCVWgtST * 
                                          kXSecCVWgt2018_smallerDISScale_NT *
                                          kQ2WgtNT),
                                         (kPPFXFluxCVWgtST * 
                                          kXSecCVWgt2018_smallerDISScale_NT * 
                                          kMySignalWgtNT)
  };


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

  std::vector<std::vector<ana::Spectrum*>> spec_analysis_3d;
  std::vector<ana::Spectrum*> spec_template;
  std::vector<ana::Spectrum*> spec_unfolding_1d;
  std::vector<ana::Spectrum*> spec_unfolding_2d;
  std::vector<ana::Spectrum*> spec_efficiency_num_1d;
  std::vector<ana::Spectrum*> spec_efficiency_num_2d;
  std::vector<ana::Spectrum*> spec_efficiency_denom_1d;
  std::vector<ana::Spectrum*> spec_efficiency_denom_2d;

  std::vector<ana::Var> ppfx_univ = ana::GetkPPFXFluxUnivWgt();
  std::vector<ana::NuTruthVar> ppfx_univ_st = ana::GetkPPFXFluxUnivWgtST();

  for(uint iload = 0; iload < loaders.size(); ++iload){
    
    std::vector<ana::Spectrum*> vspec_analysis_3d;
    for(int ichn = 0; ichn < kcNumChns; ichn++){
      vspec_analysis_3d.push_back(new Spectrum(*loaders[iload],
                                               analysis_vars[1].axis, // costheta
                                               analysis_vars[2].axis, // electron E
                                               analysis_vars[0].axis, // recoE
                                               kcPresel && chns[ichn].cut,
                                               kNoShift,
                                               weights[iload]));
    }//loop over interaction channels
    spec_analysis_3d.push_back(vspec_analysis_3d);


    spec_template.push_back(new Spectrum(*loaders[iload],
                                         analysis_vars[1].axis, // costheta
                                         analysis_vars[2].axis, // electron E
                                         pid_vars[0].axis, // prongCVN
                                         kcPresel,
                                         kNoShift,
                                         weights[iload]));

    spec_unfolding_1d.push_back(new Spectrum(*loaders[iload],
                                             analysis_vars[0].axis, // recoE
                                             analysis_vars[5].axis, // true E
                                             kcPresel &&
                                             CutFromNuTruthCut(kIsTrueSigST),
                                             kNoShift,
                                             weights[iload]));
    spec_unfolding_2d.push_back(new Spectrum(*loaders[iload],
                                             kRecoElecEVsCosStandardAxis,
                                             kTrueElecEVsCosStandardAxis,
                                             kcPresel &&
                                             CutFromNuTruthCut(kIsTrueSigST),
                                             kNoShift,
                                             weights[iload]));
    spec_efficiency_num_1d.push_back(new Spectrum(*loaders[iload],
                      HistAxisFromNuTruthHistAxis(analysis_vars_truth[0].axis),
                                              kcPresel && 
                                              CutFromNuTruthCut(kIsTrueSigST),
                                              kNoShift,
                                              weights[iload]));
    spec_efficiency_num_2d.push_back(new Spectrum(*loaders[iload],
                                             kTrueElecEVsCosStandardAxis,
                                             kcPresel &&
                                             CutFromNuTruthCut(kIsTrueSigST),
                                             kNoShift,
                                             weights[iload]));
    spec_efficiency_denom_1d.push_back(new Spectrum(*loaders[iload],
                                            analysis_vars_truth[0].axis,
                                            kIsTrueSigST,
                                            kNoShift,
                                            weights_nt[iload]));
    spec_efficiency_denom_2d.push_back(new Spectrum(*loaders[iload],
                                                kTrueElecEVsCosStandardAxisST,
                                                kIsTrueSigST,
                                                kNoShift,
                                                weights_nt[iload]));
  }

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


  TFile* outf =
    new TFile("fCrossSection_Nominal_FakeData_prongcvn_new.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", "mc", dataset_labels[iload].c_str(),
              "analysis", chns[ichn].name.c_str());
      spec_analysis_3d[iload][ichn]->SaveTo(outf, name);
    }//loop over interaction channels

    char name[50];
    sprintf(name, "%s_%s_%s_%s", "mc", dataset_labels[iload].c_str(),
            "template", "tot");
    spec_template[iload]->SaveTo(outf, name);
    sprintf(name, "%s_%s_%s_%s", "mc", dataset_labels[iload].c_str(),
            "unfolding", "1d");
    spec_unfolding_1d[iload]->SaveTo(outf, name);
    sprintf(name, "%s_%s_%s_%s", "mc", dataset_labels[iload].c_str(),
            "unfolding", "2d");
    spec_unfolding_2d[iload]->SaveTo(outf, name);
    sprintf(name, "%s_%s_%s_%s", "mc", dataset_labels[iload].c_str(),
            "efficiency", "num_1d");
    spec_efficiency_num_1d[iload]->SaveTo(outf, name);
    sprintf(name, "%s_%s_%s_%s", "mc", dataset_labels[iload].c_str(),
            "efficiency", "num_2d");
    spec_efficiency_num_2d[iload]->SaveTo(outf, name);
    sprintf(name, "%s_%s_%s_%s", "mc", dataset_labels[iload].c_str(),
            "efficiency", "denom_1d");
    spec_efficiency_denom_1d[iload]->SaveTo(outf, name);
    sprintf(name, "%s_%s_%s_%s", "mc", dataset_labels[iload].c_str(),
            "efficiency", "denom_2d");
    spec_efficiency_denom_2d[iload]->SaveTo(outf, name);
  }//loader

  outf->Close();



                                         
}