SystsPPFXLoad17.C

Go to the documentation of this file.

#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/SystShifts.h"
#include "CAFAna/Core/Utilities.h"
#include "CAFAna/Cuts/Cuts.h"
#include "NuXAna/Cuts/NusCuts.h"
#include "NuXAna/Cuts/NusCuts17.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Decomp/CheatDecomp.h"
#include "CAFAna/Decomp/ProportionalDecomp.h"
#include "CAFAna/Extrap/ExtrapSterile.h"
#include "CAFAna/Prediction/PredictionGenerator.h"
#include "NuXAna/Prediction/PredictionSterile.h"
#include "CAFAna/Systs/Systs.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "NuXAna/Vars/HistAxes.h"
#include "NuXAna/macros/NuSLoadMacroDefs.h"
#include "NuXAna/macros/NuSSystFunctions.h"
#include "CAFAna/Vars/PPFXWeights.h"

using namespace ana;


   std::string StringFromInt1(int i)
  {
    char pot[3];
    sprintf(pot, "%d", i);
    std::string potStr = std::string(pot);

    return potStr;
  }


void SystsPPFXLoad17()
{
  TH1::AddDirectory(0);


   const Var kWgtP3 = kXSecCVWgt2017*kPPFXFluxCVWgt;   
  

    const std::string fnamenear_conc="defname: prod_caf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1";
    const std::string fnameneardata_conc="defname: prod_caf_R17-03-01-prod3reco.d_nd_numi_fhc_full_v1_goodruns";
  const std::string fnamefar_con="defname: prod_caf_R17-03-01-prod3reco.l_fd_genie_nonswap_fhc_nova_v08_full_v1";
  const std::string fnameswap_con="defname: prod_caf_R17-03-01-prod3reco.l_fd_genie_fluxswap_fhc_nova_v08_full_v1";
  const std::string  fnametau_conc="defname: prod_caf_R17-03-01-prod3reco.l_fd_genie_tau_fhc_nova_v08_full_v1";


  SpectrumLoader floadnear(fnamenear_conc);
  SpectrumLoader floadfar(fnamefar_con);
  SpectrumLoader floadswap(fnameswap_con);
  SpectrumLoader floadtau(fnametau_conc);
  SpectrumLoader floadneardata(fnameneardata_conc);

  floadnear    .SetSpillCut(kStandardSpillCuts);
  floadfar     .SetSpillCut(kStandardSpillCuts);
  floadswap    .SetSpillCut(kStandardSpillCuts);
  floadtau     .SetSpillCut(kStandardSpillCuts);
  floadneardata.SetSpillCut(kStandardSpillCuts);

  // A map of systematic objects, indexed by the same label as above
/*  std::map<std::string, std::vector<SpectrumLoader*> > systematics;

  for(const auto& shiftlabel : shiftlabels) {
    systematics[shiftlabel.first] = std::vector<SpectrumLoader*>();
    for(const auto& samdef : shiftlabel.second){
      systematics[shiftlabel.first].push_back(new SpectrumLoader(samdef));
      systematics[shiftlabel.first].back()->SetSpillCut(kStandardSpillCuts);
    }
  }
*/
 //  std::map<std::string, std::string> systematics;
 // systematics["PPFX"] = "PPFX Weight";

  // Define a map of samples and cuts/axes to set them up
  // Each sample, essentially, is a way to set up the analysis
  std::string labelRecoE = "Energy Deposited in Scintillator (GeV)";
  std::map<std::string, std::pair<HistAxis*, std::pair<Cut*, Cut*> > > cut_samples;
  std::pair<Cut*, Cut*> Ana01(new Cut(kNus17ND), new Cut(kNus17FDAlt));
  cut_samples["Ana01"] = std::pair<HistAxis*, std::pair<Cut*, Cut*> >(
    new HistAxis(labelRecoE, kNCDisappearanceEnergyBinning, kNusNDEnergy17),
    Ana01
  );

  // Set up maps to the decompositions/predictions (each flavor component)
  // Nominal maps are indexed purely by sample label
  // Shifted maps are indexed by sample label, systematic label, then sigma of the shift
   std::map<std::string,                                     IDecomp*>     decomps_nominal;
   std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > > decomps_shifted;
   std::map<std::string,                                     PredictionNoExtrap*>     predsNE_nominal;
   std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > > predsNE_shifted;
   std::map<std::string,                                     PredictionSterile*>      predsSt_nominal;
   std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > >  predsSt_shifted;

  

  std::vector<Var> kPPFXFluxUnivWgts;
  for(unsigned int UnivIdx = 0; UnivIdx < 100; UnivIdx++){
        const Var tempPPFXWgt(
                  [UnivIdx](const caf::SRProxy* sr){
                    if (!sr->hdr.ismc)return 1.f;
                    if (sr->mc.nnu != 1)return 1.f;
                    if(sr->mc.nu[0].rwgt.ppfx.vuniv[UnivIdx] <= 0)return 1.f;
                    return (float)sr->mc.nu[0].rwgt.ppfx.vuniv[UnivIdx];
                  });
        kPPFXFluxUnivWgts.push_back(tempPPFXWgt);
  }



  //std::cout << "Hey I am here "<<std::cout;

  // Set up the actual decompositions/predictions
  // sample.first = the sample label
  // shiftlabel.first = the systematic shift label
  // sample.second.first = the sample HistAxis
  // sample.second.second = the sample Cut
  for(const auto& sample : cut_samples) {
    
     decomps_nominal[sample.first] = new CheatDecomp(
      floadnear,
      *sample.second.first, *sample.second.second.first,
      kNoShift, kWgtP3
    );

     predsNE_nominal[sample.first] = new PredictionNoExtrap(
      floadfar, floadswap, floadtau,
      labelRecoE, kNCDisappearanceEnergyBinning,kNusNDEnergy17,*sample.second.second.second,
      kNoShift, kWgtP3
    );



    ProportionalDecomp* decompnc_n   = new ProportionalDecomp(
      floadnear, floadneardata,
      *sample.second.first, *sample.second.second.first,
      kNoShift, kNoShift, kWgtP3
    );
    ProportionalDecomp* decompnumu_n = new ProportionalDecomp(
      floadnear, floadneardata,
      kNCBinsNumuCCAxis, kNumuND,
      kNoShift, kNoShift, kWgtP3
    );
    ModularExtrapSterile* extrap_n = new ModularExtrapSterile(ModularExtrapSterile::NCDisappearance(
      floadnear, floadswap, floadfar, floadtau, *decompnc_n, *decompnumu_n,
      *sample.second.first, kNCBinsNumuCCAxis,
      *sample.second.second.second, *sample.second.second.first, kNumuND,
      kNoShift, kWgtP3
    ));
    predsSt_nominal[sample.first] = new PredictionSterile(extrap_n);


     for(unsigned int Idx = 0; Idx < 100; Idx++){
          const Var kPPFXFluxWgt = kPPFXFluxUnivWgts[Idx];
         const Var kMCUnivWeight = kPPFXFluxWgt*kXSecCVWgt2017;

       decomps_shifted[sample.first][StringFromInt1(Idx)][1] = new CheatDecomp(
          floadnear,
          *sample.second.first, *sample.second.second.first,
          kNoShift, kMCUnivWeight
        );
       

        predsNE_shifted[sample.first][StringFromInt1(Idx)][1] = new PredictionNoExtrap(
        floadfar, floadswap, floadtau,
        labelRecoE, kNCDisappearanceEnergyBinning ,
        kNusNDEnergy17, *sample.second.second.second,
        kNoShift, kMCUnivWeight
      );

      ProportionalDecomp* decompnc_s   = new ProportionalDecomp(
        floadnear, floadneardata,
        *sample.second.first, *sample.second.second.first,
        kNoShift, kNoShift, kMCUnivWeight
      );
      ProportionalDecomp* decompnumu_s = new ProportionalDecomp(
        floadnear, floadneardata,
        kNCBinsNumuCCAxis, kNumuND,
        kNoShift, kNoShift, kMCUnivWeight
      );
      ModularExtrapSterile* extrap_s = new ModularExtrapSterile(ModularExtrapSterile::NCDisappearance(
        floadnear, floadswap, floadfar, floadtau, *decompnc_s, *decompnumu_s,
        *sample.second.first, kNCBinsNumuCCAxis,
        *sample.second.second.second, *sample.second.second.first, kNumuND,
        kNoShift, kMCUnivWeight
      ));
      predsSt_shifted[sample.first][StringFromInt1(Idx)][1] = new PredictionSterile(
        extrap_s
      );
    }
  }

 /* for(const auto& systlabel : systematics) {
    for(const auto& loader : systlabel.second) {
      loader->Go();
    }
  }*/
  floadnear.Go();
  floadfar .Go();
  floadswap.Go();
  floadtau .Go();
  floadneardata.Go();
 /* floadnear_ns.Go();
  floadfar_ns .Go();
  floadswap_ns.Go();
  floadtau_ns.Go();
*/
  std::string folder = "./";
  std::string filenm = "ppfx_debug";

  std::string fullLocation = folder + filenm + ".root";
  TFile* rootF = new TFile(fullLocation.c_str(), "RECREATE");

  SaveMaps(rootF,
    decomps_nominal, decomps_shifted,
    predsNE_nominal, predsNE_shifted,
    predsSt_nominal, predsSt_shifted
  );
}