runCheatDecomp.C

Go to the documentation of this file.

#include "3FlavorAna/Cuts/NumuCuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/HistAxis.h"
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Decomp/CheatDecomp.h"
#include "CAFAna/Core/ISyst.h"
#include "CAFAna/Systs/Systs.h"
#include "CAFAna/Core/Utilities.h"

#include "nugen/NuReweight/ReweightLabels.h"



#include <iostream>
#include <sstream>
#include <map>
#include <utility>

#include "TBrowser.h"
#include "TFile.h"
#include "3FlavorAna/Vars/NumuVars.h"

using namespace ana;


void runCheatDecomp(TString fileListTextFileName)
{

   
    std::map<std::string, rwgt::ReweightLabel_t> parameters;       
     // NCEL tweaking parameters:
     parameters["ReweightMaNCEL"]        =   rwgt::fReweightMaNCEL;              ///< tweak Ma NCEL, affects dsigma(NCEL)/dQ2 both in shape and normalization
     parameters["ReweightEtaNCEL"]        =   rwgt::fReweightEtaNCEL;            ///< tweak NCEL strange axial form factor eta, affects dsigma(NCEL)/dQ2 both in shape and normalization
     // CCQE tweaking parameters:
     parameters["ReweightNormCCQE"]        =   rwgt::fReweightNormCCQE;          ///< tweak CCQE normalization (energy independent)
     parameters["ReweightNormCCQEenu"]        =   rwgt::fReweightNormCCQEenu;    ///< tweak CCQE normalization (maintains dependence on neutrino energy)
     parameters["ReweightMaCCQEshape"]        =   rwgt::fReweightMaCCQEshape;    ///< tweak Ma CCQE, affects dsigma(CCQE)/dQ2 in shape only (normalized to constant integral)
     parameters["ReweightMaCCQE"]        =   rwgt::fReweightMaCCQE;              ///< tweak Ma CCQE, affects dsigma(CCQE)/dQ2 both in shape and normalization
     parameters["ReweightVecCCQEshape"]        =   rwgt::fReweightVecCCQEshape;///< tweak elastic nucleon form factors (BBA/default -> dipole) - shape only effect of dsigma(CCQE)/dQ2
     // Resonance neutrino-production tweaking parameters:
     parameters["ReweightNormCCRES"]        =   rwgt::fReweightNormCCRES;         ///< tweak CCRES normalization
     parameters["ReweightMaCCRESshape"]        =   rwgt::fReweightMaCCRESshape;   ///< tweak Ma CCRES, affects d2sigma(CCRES)/dWdQ2 in shape only (normalized to constant integral)
     parameters["ReweightMvCCRESshape"]        =   rwgt::fReweightMvCCRESshape;   ///< tweak Mv CCRES, affects d2sigma(CCRES)/dWdQ2 in shape only (normalized to constant integral)
     parameters["ReweightMaCCRES"]        =   rwgt::fReweightMaCCRES;             ///< tweak Ma CCRES, affects d2sigma(CCRES)/dWdQ2 both in shape and normalization
     parameters["ReweightMvCCRES"]        =   rwgt::fReweightMvCCRES;             ///< tweak Mv CCRES, affects d2sigma(CCRES)/dWdQ2 both in shape and normalization
     
     parameters["ReweightNormNCRES"]        =   rwgt::fReweightNormNCRES;         ///< tweak NCRES normalization
     parameters["ReweightMaNCRESshape"]        =   rwgt::fReweightMaNCRESshape;   ///< tweak Ma NCRES, affects d2sigma(NCRES)/dWdQ2 in shape only (normalized to constant integral)
     parameters["ReweightMvNCRESshape"]        =   rwgt::fReweightMvNCRESshape;   ///< tweak Mv NCRES, affects d2sigma(NCRES)/dWdQ2 in shape only (normalized to constant integral)
     parameters["ReweightMaNCRES"]        =   rwgt::fReweightMaNCRES;             ///< tweak Ma NCRES, affects d2sigma(NCRES)/dWdQ2 both in shape and normalization
     parameters["ReweightMvNCRES"]        =   rwgt::fReweightMvNCRES;             ///< tweak Mv NCRES, affects d2sigma(NCRES)/dWdQ2 both in shape and normalization
     
     // Coherent pion production tweaking parameters:
     parameters["ReweightMaCOHpi"]        =   rwgt::fReweightMaCOHpi;             ///< tweak Ma for COH pion production
     parameters["ReweightR0COHpi"]        =   rwgt::fReweightR0COHpi;             ///< tweak R0 for COH pion production
     // Non-resonance background tweaking parameters:
     parameters["ReweightRvpCC1pi"]        =   rwgt::fReweightRvpCC1pi;           ///< tweak the 1pi non-RES bkg in the RES region, for v+p CC
     parameters["ReweightRvpCC2pi"]        =   rwgt::fReweightRvpCC2pi;           ///< tweak the 2pi non-RES bkg in the RES region, for v+p CC
     parameters["ReweightRvpNC1pi"]        =   rwgt::fReweightRvpNC1pi;           ///< tweak the 1pi non-RES bkg in the RES region, for v+p NC
     parameters["ReweightRvpNC2pi"]        =   rwgt::fReweightRvpNC2pi;           ///< tweak the 2pi non-RES bkg in the RES region, for v+p NC
     parameters["ReweightRvnCC1pi"]        =   rwgt::fReweightRvnCC1pi;           ///< tweak the 1pi non-RES bkg in the RES region, for v+n CC
     parameters["ReweightRvnCC2pi"]        =   rwgt::fReweightRvnCC2pi;           ///< tweak the 2pi non-RES bkg in the RES region, for v+n CC
     parameters["ReweightRvnNC1pi"]        =   rwgt::fReweightRvnNC1pi;           ///< tweak the 1pi non-RES bkg in the RES region, for v+n NC
     parameters["ReweightRvnNC2pi"]        =   rwgt::fReweightRvnNC2pi;           ///< tweak the 2pi non-RES bkg in the RES region, for v+n NC
     parameters["ReweightRvbarpCC1pi"]        =   rwgt::fReweightRvbarpCC1pi;     ///< tweak the 1pi non-RES bkg in the RES region, for vbar+p CC
     parameters["ReweightRvbarpCC2pi"]        =   rwgt::fReweightRvbarpCC2pi;     ///< tweak the 2pi non-RES bkg in the RES region, for vbar+p CC
     parameters["ReweightRvbarpNC1pi"]        =   rwgt::fReweightRvbarpNC1pi;     ///< tweak the 1pi non-RES bkg in the RES region, for vbar+p NC
     parameters["ReweightRvbarpNC2pi"]        =   rwgt::fReweightRvbarpNC2pi;     ///< tweak the 2pi non-RES bkg in the RES region, for vbar+p NC
     parameters["ReweightRvbarnCC1pi"]        =   rwgt::fReweightRvbarnCC1pi;     ///< tweak the 1pi non-RES bkg in the RES region, for vbar+n CC
     parameters["ReweightRvbarnCC2pi"]        =   rwgt::fReweightRvbarnCC2pi;     ///< tweak the 2pi non-RES bkg in the RES region, for vbar+n CC
     parameters["ReweightRvbarnNC1pi"]        =   rwgt::fReweightRvbarnNC1pi;     ///< tweak the 1pi non-RES bkg in the RES region, for vbar+n NC
     parameters["ReweightRvbarnNC2pi"]        =   rwgt::fReweightRvbarnNC2pi;     ///< tweak the 2pi non-RES bkg in the RES region, for vbar+n NC
     // DIS tweaking parameters - applied for DIS events with (Q2>Q2o, W>Wo), 
     // typically Q2okReweight =1GeV^2, WokReweight =1.7-2.0GeV
     parameters["ReweightAhtBY"]        =   rwgt::fReweightAhtBY;                 ///< tweak the Bodek-Yang model parameter A_{ht} - incl. both shape and normalization effect
     parameters["ReweightBhtBY"]        =   rwgt::fReweightBhtBY;                 ///< tweak the Bodek-Yang model parameter B_{ht} - incl. both shape and normalization effect 
     parameters["ReweightCV1uBY"]        =   rwgt::fReweightCV1uBY;               ///< tweak the Bodek-Yang model parameter CV1u - incl. both shape and normalization effect 
     parameters["ReweightCV2uBY"]        =   rwgt::fReweightCV2uBY;               ///< tweak the Bodek-Yang model parameter CV2u - incl. both shape and normalization effect 
     parameters["ReweightAhtBYshape"]        =   rwgt::fReweightAhtBYshape;       ///< tweak the Bodek-Yang model parameter A_{ht} - shape only effect to d2sigma(DIS)/dxdy
     parameters["ReweightBhtBYshape"]        =   rwgt::fReweightBhtBYshape;       ///< tweak the Bodek-Yang model parameter B_{ht} - shape only effect to d2sigma(DIS)/dxdy
     parameters["ReweightCV1uBYshape"]        =   rwgt::fReweightCV1uBYshape;     ///< tweak the Bodek-Yang model parameter CV1u - shape only effect to d2sigma(DIS)/dxdy
     parameters["ReweightCV2uBYshape"]        =   rwgt::fReweightCV2uBYshape;     ///< tweak the Bodek-Yang model parameter CV2u - shape only effect to d2sigma(DIS)/dxdy
     parameters["ReweightNormDISCC"]        =   rwgt::fReweightNormDISCC;         ///< tweak the inclusive DIS CC normalization (not currently working in genie)
     parameters["ReweightRnubarnuCC"]        =   rwgt::fReweightRnubarnuCC;       ///< tweak the ratio of \sigma(\bar\nu CC) / \sigma(\nu CC) (not currently working in genie)
     parameters["ReweightDISNuclMod"]        =   rwgt::fReweightDISNuclMod;       ///< tweak DIS nuclear modification (shadowing, anti-shadowing, EMC).  Does not appear to be working in GENIE at the moment
     //
     parameters["ReweightNC"]        =   rwgt::fReweightNC;                ///< 
     
     
     //
     // Hadronization (free nucleon target)
     // 
     
     parameters["ReweightAGKY_xF1pi"]        =   rwgt::fReweightAGKY_xF1pi;         ///< tweak xF distribution for low multiplicity (N + pi) DIS f/s produced by AGKY
     parameters["ReweightAGKY_pT1pi"]        =   rwgt::fReweightAGKY_pT1pi;         ///< tweak pT distribution for low multiplicity (N + pi) DIS f/s produced by AGKY
     
     
     //
     // Medium-effects to hadronization
     // 
     
     parameters["ReweightFormZone"]        =   rwgt::fReweightFormZone;         ///< tweak formation zone
     
     
     //
     // Intranuclear rescattering systematics.
     // There are 2 sets of parameters:
     // - parameters that control the total rescattering probability, P(total)
     // - parameters that control the fraction of each process (`fate'), given a total rescat. prob., P(fate|total)
     // These parameters are considered separately for pions and nucleons.
     //
     
     parameters["ReweightMFP_pi"]        =   rwgt::fReweightMFP_pi;           ///< tweak mean free path for pions
     parameters["ReweightMFP_N"]        =   rwgt::fReweightMFP_N;             ///< tweak mean free path for nucleons
     parameters["ReweightFrCEx_pi"]        =   rwgt::fReweightFrCEx_pi;       ///< tweak charge exchange probability for pions, for given total rescattering probability
     parameters["ReweightFrElas_pi"]        =   rwgt::fReweightFrElas_pi;     ///< tweak elastic   probability for pions, for given total rescattering probability
     parameters["ReweightFrInel_pi"]        =   rwgt::fReweightFrInel_pi;     ///< tweak inelastic probability for pions, for given total rescattering probability
     parameters["ReweightFrAbs_pi"]        =   rwgt::fReweightFrAbs_pi;       ///< tweak absorption probability for pions, for given total rescattering probability
     parameters["ReweightFrPiProd_pi"]        =   rwgt::fReweightFrPiProd_pi; ///< tweak pion production probability for pions, for given total rescattering probability
     parameters["ReweightFrCEx_N"]        =   rwgt::fReweightFrCEx_N;         ///< tweak charge exchange probability for nucleons, for given total rescattering probability
     parameters["ReweightFrElas_N"]        =   rwgt::fReweightFrElas_N;       ///< tweak elastic    probability for nucleons, for given total rescattering probability
     parameters["ReweightFrInel_N"]        =   rwgt::fReweightFrInel_N;       ///< tweak inelastic  probability for nucleons, for given total rescattering probability
     parameters["ReweightFrAbs_N"]        =   rwgt::fReweightFrAbs_N;         ///< tweak absorption probability for nucleons, for given total rescattering probability
     parameters["ReweightFrPiProd_N"]        =   rwgt::fReweightFrPiProd_N;   ///< tweak pion production probability for nucleons, for given total rescattering probability
     
     //
     // Nuclear model
     // 
     
     parameters["ReweightCCQEPauliSupViaKF"]        =   rwgt::fReweightCCQEPauliSupViaKF;   ///<
     parameters["ReweightCCQEMomDistroFGtoSF"]        =   rwgt::fReweightCCQEMomDistroFGtoSF; ///<
     
     //
     // Resonance decays
     // 
     
     parameters["ReweightBR1gamma"]        =   rwgt::fReweightBR1gamma;               ///< tweak Resonance -> X + gamma branching ratio, eg Delta+(1232) -> p gamma
     parameters["ReweightBR1eta"]        =   rwgt::fReweightBR1eta;                   ///< tweak Resonance -> X + eta   branching ratio, eg N+(1440) -> p eta
     parameters["ReweightTheta_Delta2Npi"]        =   rwgt::fReweightTheta_Delta2Npi; ///< distort pi angular distribution in Delta -> N + pi



  int nParameters = 0;
  for(auto parameter:parameters) nParameters++;
  std::cout << "N Parameters : " << nParameters << std::endl;

  std::string fileListTextFile =  (std::string)fileListTextFileName;
  assert(fileListTextFile.find(".txt") != std::string::npos && "Failed to find .txt in argument.  We're going to to replace it with .root for the output file, so it needs to be there.");
  std::vector<std::string> fileList = LoadFileList(fileListTextFile);
  SpectrumLoader loader (fileList);
  

   
  std::map<std::string, std::map<int, SystShifts*> > shifts; 
  std::vector<int> sigmas; 
  sigmas.push_back(-2);
  sigmas.push_back(-1);
  sigmas.push_back(1);
  sigmas.push_back(2);

  for (const auto & parameter:parameters)
  { 
    for(const auto & sigma:sigmas)
    {

      shifts[parameter.first][sigma] = new SystShifts(GetGenieKnobSyst(parameter.second), sigma);

    }
  }

  const Binning bins = Binning::Simple(50, 0, 5);

  std::map<std::string, std::pair<Cut*, HistAxis*> > samples;
//  samples["qe"] = std::pair<Cut*, HistAxis* >(new Cut(kNumuNCRej && kNumuContainND && kNumuQuality && kQE), new HistAxis("QE Energy [GeV]", bins, kQEE));
//  samples["nonqe"] = std::pair<Cut*, HistAxis* >(new Cut(kNumuNCRej && kNumuContainND && kNumuQuality && !kQE), new HistAxis("Non-QE Energy [GeV]", bins, kCCE));
  samples["allcc"] = std::pair<Cut*, HistAxis* >(new Cut(kNumuNCRej && kNumuContainND && kNumuQuality), new HistAxis("Reconstructed Neutrino Energy [GeV]", bins, kCCE));

  std::map<std::string, CheatDecomp*> decomps;

  std::map<std::string, std::map<std::string, std::map<int, CheatDecomp*> > > shiftedDecomps;


  
  for(const auto & sample:samples){
    std::cout << "Sample: " << sample.first << std::endl;
    decomps[sample.first] = new CheatDecomp(loader, *sample.second.second, *sample.second.first);  
    for(const auto & parameter:parameters)
    {
      for(const auto & shift:shifts[parameter.first])
      {
        shiftedDecomps[sample.first][parameter.first][shift.first] = new CheatDecomp(loader, *sample.second.second, *sample.second.first, *shift.second);  

      }
    }
  }
  
  loader.Go();


  std::map<std::string, TDirectory*> dirs;

  TFile outputFile(fileListTextFile.replace(fileListTextFile.find(".txt"), 4, ".root").c_str(), "RECREATE");
  for(const auto & decomp:decomps){

    
    dirs[decomp.first] = outputFile.mkdir(decomp.first.c_str());
    decomp.second->SaveTo(dirs[decomp.first]);

  }

  for(const auto & parameter:shiftedDecomps)
  {
    for(const auto & shift:parameter.second)
    {
      for(const auto & shiftedDecomp:shift.second)
      {
        std::stringstream dirName("");
        dirName << parameter.first << "_" << shift.first << "_" << shiftedDecomp.first ;  
        dirs[dirName.str()] = outputFile.mkdir(dirName.str().c_str());
        shiftedDecomp.second->SaveTo(dirs[dirName.str()]);        
      }
    }

  }
}