ppfx_smooth_weights_save.C

Go to the documentation of this file.

#include "StandardRecord/Proxy/SRProxy.h"

#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/SpectrumLoaderBase.h"
#include "CAFAna/Core/HistAxis.h"
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Cut.h"
#include "CAFAna/Core/SystShifts.h"
#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "CAFAna/Cuts/TruthCuts.h"
#include "CAFAna/Cuts/NueCutsSecondAna.h"
#include "CAFAna/Cuts/TimingCuts.h"
#include "CAFAna/Systs/Systs.h"

#include "TStopwatch.h"
#include "CAFAna/Core/Utilities.h"
#include "Utilities/func/MathUtil.h"

#include "TFile.h"
#include "TRandom.h"

#include <fstream>
#include <iostream>
#include <cmath>


using namespace ana;
  //----------------------------------------------------------------------

void ppfx_smooth_weights_save()
{
  const Cut kIsNumu(
            [](const caf::SRProxy* sr)
            {
              if(sr->mc.nnu == 0) return false;
              assert(sr->mc.nnu == 1);
              return (abs(sr->mc.nu[0].pdg) == 14);
            });
  
  const Cut kIsNue(
            [](const caf::SRProxy* sr)
            {
              if(sr->mc.nnu == 0) return false;
              assert(sr->mc.nnu == 1);
              return (abs(sr->mc.nu[0].pdg) == 12);
            });

  const Cut kIsRock(
            [](const caf::SRProxy* sr)
            {
              if (sr->mc.nnu == 0)
              return false;
              if(fabs(sr->mc.nu[0].vtx.X())>180||fabs(sr->mc.nu[0].vtx.Y())>180||sr->mc.nu[0].vtx.Z()<0||sr->mc.nu[0].vtx.Z()>1250)return true;
              return false;
            });
  

  const Var kTrueE(
            [](const caf::SRProxy* sr){
              if(sr->mc.nnu!=1)return (float)-5.0;
              return (float)sr->mc.nu[0].E;
            });


  std::vector<Var> kPPFXFluxUnivWgt;
  
  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];
                  });

        kPPFXFluxUnivWgt.push_back(tempPPFXWgt);
  }

  const Binning bins = Binning::Simple(200, 0, 10);
  const HistAxis TrueEAxis("True Energy (GeV)", bins, kTrueE);

  struct CompDef
  {
    std::string name;
    Cut cut;
  };
  
  std::vector <CompDef> comps = {
    {"numu", kIsNumu && !kIsAntiNu},
    {"nue", kIsNue && !kIsAntiNu},
    {"anumu", kIsNumu && kIsAntiNu},
    {"anue", kIsNue && kIsAntiNu},
  };

  std::map<std::string, Spectrum*> spects;

  const std::string nearNonSFHC = "prod_caf_R17-11-14-prod4reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1";
  const std::string farNonSFHC = "prod_caf_R17-11-14-prod4reco.d_fd_genie_nonswap_fhc_nova_v08_full_v1";
  const std::string nearNonSRHC = "prod_caf_R17-11-14-prod4reco.e_nd_genie_nonswap_rhc_nova_v08_full_v1";
  const std::string farNonSRHC = "prod_caf_R17-11-14-prod4reco.e_fd_genie_nonswap_rhc_nova_v08_full_v1";

  SpectrumLoader loaderNDMCFHC(nearNonSFHC);
  SpectrumLoader loaderFDMCFHC(farNonSFHC);
  SpectrumLoader loaderNDMCRHC(nearNonSRHC);
  SpectrumLoader loaderFDMCRHC(farNonSRHC);

  loaderNDMCFHC.SetSpillCut(kStandardDQCuts);
  loaderFDMCFHC.SetSpillCut(kStandardDQCuts);
  loaderNDMCRHC.SetSpillCut(kStandardDQCuts);
  loaderFDMCRHC.SetSpillCut(kStandardDQCuts);

  int N_PPFX = 100;
  for(int compIdx = 0; compIdx < (int) comps.size(); ++compIdx){
    Spectrum* temp_unweighted_ndFHC = new Spectrum(loaderNDMCFHC, TrueEAxis, comps[compIdx].cut && !kIsRock, kNoShift, kUnweighted);
    Spectrum* temp_unweighted_fdFHC = new Spectrum(loaderFDMCFHC, TrueEAxis, comps[compIdx].cut, kNoShift, kUnweighted);
    Spectrum* temp_unweighted_ndRHC = new Spectrum(loaderNDMCRHC, TrueEAxis, comps[compIdx].cut && !kIsRock, kNoShift, kUnweighted);
    Spectrum* temp_unweighted_fdRHC = new Spectrum(loaderFDMCRHC, TrueEAxis, comps[compIdx].cut, kNoShift, kUnweighted);

    spects.insert(std::make_pair("unweighted_FHC_ND_"+comps[compIdx].name, temp_unweighted_ndFHC));
    spects.insert(std::make_pair("unweighted_FHC_FD_"+comps[compIdx].name, temp_unweighted_fdFHC));
    spects.insert(std::make_pair("unweighted_RHC_ND_"+comps[compIdx].name, temp_unweighted_ndRHC));
    spects.insert(std::make_pair("unweighted_RHC_FD_"+comps[compIdx].name, temp_unweighted_fdRHC));

    Spectrum* temp_ppfx_cv_ndFHC = new Spectrum(loaderNDMCFHC, TrueEAxis, comps[compIdx].cut && !kIsRock, kNoShift, kPPFXFluxCVWgt);
    Spectrum* temp_ppfx_cv_fdFHC = new Spectrum(loaderFDMCFHC, TrueEAxis, comps[compIdx].cut, kNoShift, kPPFXFluxCVWgt);
    Spectrum* temp_ppfx_cv_ndRHC = new Spectrum(loaderNDMCRHC, TrueEAxis, comps[compIdx].cut && !kIsRock, kNoShift, kPPFXFluxCVWgt);
    Spectrum* temp_ppfx_cv_fdRHC = new Spectrum(loaderFDMCRHC, TrueEAxis, comps[compIdx].cut, kNoShift, kPPFXFluxCVWgt);

    spects.insert(std::make_pair("ppfx_cv_FHC_ND_"+comps[compIdx].name, temp_ppfx_cv_ndFHC));
    spects.insert(std::make_pair("ppfx_cv_FHC_FD_"+comps[compIdx].name, temp_ppfx_cv_fdFHC));
    spects.insert(std::make_pair("ppfx_cv_RHC_ND_"+comps[compIdx].name, temp_ppfx_cv_ndRHC));
    spects.insert(std::make_pair("ppfx_cv_RHC_FD_"+comps[compIdx].name, temp_ppfx_cv_fdRHC));

    for(int UnivIdx = 0; UnivIdx < N_PPFX; UnivIdx++){
      const Var kPPFXFluxWgt = kPPFXFluxUnivWgt[UnivIdx];
      Spectrum* temp_ppfx_univ_ndFHC = new Spectrum(loaderNDMCFHC, TrueEAxis, comps[compIdx].cut && !kIsRock, kNoShift, kPPFXFluxWgt);
      Spectrum* temp_ppfx_univ_fdFHC = new Spectrum(loaderFDMCFHC, TrueEAxis, comps[compIdx].cut, kNoShift, kPPFXFluxWgt);
      Spectrum* temp_ppfx_univ_ndRHC = new Spectrum(loaderNDMCRHC, TrueEAxis, comps[compIdx].cut && !kIsRock, kNoShift, kPPFXFluxWgt);
      Spectrum* temp_ppfx_univ_fdRHC = new Spectrum(loaderFDMCRHC, TrueEAxis, comps[compIdx].cut, kNoShift, kPPFXFluxWgt);

      int dig0 = UnivIdx % 10;
      int dig1 = (UnivIdx/10) % 10;
      spects.insert(std::make_pair("ppfx_univ"+std::to_string(dig1)+std::to_string(dig0)+"_FHC_ND_"+comps[compIdx].name, temp_ppfx_univ_ndFHC));
      spects.insert(std::make_pair("ppfx_univ"+std::to_string(dig1)+std::to_string(dig0)+"_FHC_FD_"+comps[compIdx].name, temp_ppfx_univ_fdFHC));
      spects.insert(std::make_pair("ppfx_univ"+std::to_string(dig1)+std::to_string(dig0)+"_RHC_ND_"+comps[compIdx].name, temp_ppfx_univ_ndRHC));
      spects.insert(std::make_pair("ppfx_univ"+std::to_string(dig1)+std::to_string(dig0)+"_RHC_FD_"+comps[compIdx].name, temp_ppfx_univ_fdRHC));

    }
  }
  
  TStopwatch gtimer;
  gtimer.Start();
  loaderNDMCFHC.Go();
  loaderFDMCFHC.Go();
  loaderNDMCRHC.Go();
  loaderFDMCRHC.Go();
  gtimer.Stop();
  std::cout<<"loading time "<<gtimer.CpuTime()<<std::endl;


  std::string filename = "ppfx_weights.root";
  TFile* fout = new TFile(filename.c_str(), "RECREATE");
  std::for_each(spects.begin(), spects.end(),
                  [&](std::pair<std::string, Spectrum*> spect){
                   spect.second->SaveTo(fout, spect.first.c_str());
               });
  delete fout;


}