gRwghtZExpDirect.cxx File Reference

#include <string>
#include <sstream>
#include <TSystem.h>
#include <TFile.h>
#include <TTree.h>
#include <TArrayF.h>
#include "Algorithm/AlgConfigPool.h"
#include "Algorithm/AlgFactory.h"
#include "Base/XSecAlgorithmI.h"
#include "Conventions/Controls.h"
#include "EVGCore/EventRecord.h"
#include "Ntuple/NtpMCFormat.h"
#include "Ntuple/NtpMCTreeHeader.h"
#include "Ntuple/NtpMCEventRecord.h"
#include "Messenger/Messenger.h"
#include "ReWeight/GReWeightI.h"
#include "ReWeight/GSystSet.h"
#include "ReWeight/GReWeight.h"
#include "ReWeight/GReWeightNuXSecCCQE.h"
#include "ReWeight/GReWeightNuXSecCCQEvec.h"
#include "ReWeight/GReWeightNuXSecCCRES.h"
#include "ReWeight/GReWeightNuXSecNCRES.h"
#include "ReWeight/GReWeightNuXSecDIS.h"
#include "ReWeight/GReWeightNuXSecCOH.h"
#include "ReWeight/GReWeightNonResonanceBkg.h"
#include "ReWeight/GReWeightFGM.h"
#include "ReWeight/GReWeightDISNuclMod.h"
#include "ReWeight/GReWeightResonanceDecay.h"
#include "ReWeight/GReWeightFZone.h"
#include "ReWeight/GReWeightINuke.h"
#include "ReWeight/GReWeightAGKY.h"
#include "ReWeight/GSystUncertainty.h"
#include "Utils/CmdLnArgParser.h"
#include "Utils/StringUtils.h"

Include dependency graph for gRwghtZExpDirect.cxx:

Go to the source code of this file.

Macros
#define	MAX_COEF   4
	A simple program to reweight the GENIE z-expansion axial form factor from one set of z-expansion parameters directly to another. More...

Functions
void	PrintSyntax ()
void	GetEventRange (Long64_t nev_in_file, Long64_t &nfirst, Long64_t &nlast)
void	GetCommandLineArgs (int argc, char **argv)
GSyst_t	GetZExpSystematic (int ip)
int	main (int argc, char **argv)

Variables
string	gOptInpFilename
string	gOptOutFilename
Long64_t	gOptNEvt1
Long64_t	gOptNEvt2
bool	gOptDoNorm = false
double	gOptNormValue = 1.
double	gOptParameterValue [MAX_COEF] = {0.}

Macro Definition Documentation

#define MAX_COEF   4

A simple program to reweight the GENIE z-expansion axial form factor from one set of z-expansion parameters directly to another.

gRwghtZExpDirect

grwghtzexpaxff -f filename -v val1,val2,val3,val4 [-n nev] [-o fileOutName] [-m norm]

     where 
     [] is an optional argument
     -f specifies a GENIE event file (GHEP format)
     -o specifies a GENIE output filename
     -n specifies the number of events to process (default: all)
     -v z-expansion values to reweight to
     -m reweight normalization of form factor (default: 1)

Author

Aaron Meyer <asmeyer2012 uchicago.edu> University of Chicago, Fermi National Accelerator Laboratory

based on gtestRewght by

Costas Andreopoulos <c.andreopoulos cern.ch> University of Liverpool

Created:

Mar 14, 2016

License:

Copyright (c) 2003-2024, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org

Definition at line 74 of file gRwghtZExpDirect.cxx.

Referenced by main().

Function Documentation

void GetCommandLineArgs	(	int	argc,
		char **	argv
	)

void GetEventRange	(	Long64_t	nev_in_file,
		Long64_t &	nfirst,
		Long64_t &	nlast
	)

GSyst_t GetZExpSystematic

(

int

ip

)

int main	(	int	argc,
		char **	argv
	)

Definition at line 95 of file gRwghtZExpDirect.cxx.

References genie::AlgFactory::AdoptAlgorithm(), genie::Algorithm::AdoptSubstructure(), genie::NtpMCEventRecord::Clear(), genie::NtpMCEventRecord::event, genie::gAbortingInErr, GetCommandLineArgs(), genie::Algorithm::GetConfig(), genie::Registry::GetDouble(), genie::Registry::GetDoubleDef(), GetEventRange(), GetZExpSystematic(), gOptDoNorm, gOptInpFilename, gOptNormValue, gOptOutFilename, gOptParameterValue, genie::AlgFactory::Instance(), genie::AlgConfigPool::Instance(), LOG, MAX_COEF, pFATAL, pINFO, pNOTICE, and PrintSyntax().

{
  GetCommandLineArgs (argc, argv);

  // open the ROOT file and get the TTree & its header
  TTree *           tree = 0;
  NtpMCTreeHeader * thdr = 0;
  TFile file(gOptInpFilename.c_str(),"READ");
  tree = dynamic_cast <TTree *>           ( file.Get("gtree")  );
  thdr = dynamic_cast <NtpMCTreeHeader *> ( file.Get("header") );
  LOG("rwghtzexpaxff", pNOTICE) << "Input tree header: " << *thdr;
  if(!tree){
    LOG("grwghtzexpaxff", pFATAL)
      << "Can't find a GHEP tree in input file: "<< file.GetName();
    gAbortingInErr = true;
    PrintSyntax();
    exit(1);
  }
  NtpMCEventRecord * mcrec = 0;
  tree->SetBranchAddress("gmcrec", &mcrec);

  Long64_t nev_in_file = tree->GetEntries();
  Long64_t nfirst = 0;
  Long64_t nlast  = 0;
  GetEventRange(nev_in_file, nfirst, nlast);
  int nev = int(nlast - nfirst + 1);

  LOG("rwghtzexpaxff", pNOTICE) << "Will process " << nev << " events";

  //
  // Create a GReWeight object and add to it a set of 
  // weight calculators
  //
  // If seg-faulting here, need to change
  // AxialFormFactorModel in UserPhysicsOptions.xml and LwlynSmithFFCC.xml
  //

  GReWeight rw;
  rw.AdoptWghtCalc( "xsec_ccqe",       new GReWeightNuXSecCCQE      );

  //
  // Create a list of systematic params (more to be found at GSyst.h)
  // set non-default values and re-configure.
  // Weight calculators included above must be able to handle the tweaked params.
  // Each tweaking dial t modifies a physics parameter p as:
  // p_{tweaked} = p_{default} ( 1 + t * dp/p )
  // So setting a tweaking dial to +/-1 modifies a physics quantity
  // by +/- 1sigma.
  // Default fractional errors are defined in GSystUncertainty
  // and can be overriden.
  //

  GSystSet & syst = rw.Systematics();

  // Create a concrete weight calculator to fine-tune
  GReWeightNuXSecCCQE * rwccqe = 
    dynamic_cast<GReWeightNuXSecCCQE *> (rw.WghtCalc("xsec_ccqe"));
  rwccqe->SetMode(GReWeightNuXSecCCQE::kModeZExp);
  // In case uncertainties need to be altered
  GSystUncertainty * unc = GSystUncertainty::Instance();

  // Set up algorithm for loading central values
  AlgFactory * algf = AlgFactory::Instance();
  AlgId id("genie::LwlynSmithQELCCPXSec","ZExp");

  Algorithm * alg = algf->AdoptAlgorithm(id);
  XSecAlgorithmI* fXSecModel = dynamic_cast<XSecAlgorithmI*>(alg);
  fXSecModel->AdoptSubstructure();

  Registry * fXSecModelConfig = new Registry(fXSecModel->GetConfig());
  AlgConfigPool * confp = AlgConfigPool::Instance();
  const Registry * gc = confp->GlobalParameterList();

  // further optional fine-tuning
  //rwccqe -> RewNue    (false); 
  //rwccqe -> RewNuebar (false); 
  //rwccqe -> RewNumubar(false); 

  // Declare the weights and twkvals arrays 
  const int n_events = (const int) nev;
  const int n_params = (const int) MAX_COEF;
  // if segfaulting here, may need to increase MAX_COEF
  float weights [n_events];

  // Initialize
  for (int iev = 0; iev < nev; iev++) { weights[iev] = 1.; }

  // set first values for weighting
  if (gOptDoNorm)
  {
    //LOG("rwghtzexpaxff", pNOTICE) << "Setting z-expansion tweak for norm : "
    //  << (gOptNormValue-1.)
    //  /unc->OneSigmaErr(kXSecTwkDial_ZNormCCQE,TMath::Sign(1,(int)(gOptNormValue-1.)));
    syst.Set(kXSecTwkDial_ZNormCCQE, (gOptNormValue-1.)
     /unc->OneSigmaErr(kXSecTwkDial_ZNormCCQE,TMath::Sign(1,(int)(gOptNormValue-1.))));
  }
  GSyst_t gsyst;
  double cval = 0.; // central value for parameter
  ostringstream zval_name; // string to use to extract central value
  for (int ipr = 0; ipr < n_params; ipr++)
  {
    gsyst = GetZExpSystematic(ipr+1); // get tweak dial
    zval_name.str("");
    zval_name << "QEL-Z_A" << ipr+1;

    cval = fXSecModelConfig->GetDoubleDef(zval_name.str(),gc->GetDouble(zval_name.str()));
    unc->SetUncertainty(gsyst,1.,1.); // easier to deal with
    
    //LOG("rwghtzexpaxff", pNOTICE) << "Setting z-expansion tweak for param " 
    //  <<ipr<<" : " << (gOptParameterValue[ipr]-cval)/cval;
    syst.Set(gsyst, (gOptParameterValue[ipr]-cval)/cval);
  }

  rw.Reconfigure();
  // Event loop
  for(int iev = nfirst; iev <= nlast; iev++) {
    tree->GetEntry(iev);
  
    EventRecord & event = *(mcrec->event);
    LOG("rwghtzexpaxff", pNOTICE) << "Event number   = " << iev;
    LOG("rwghtzexpaxff", pNOTICE) << event;
  
    double wght = rw.CalcWeight(event);

    LOG("rwghtzexpaxff", pNOTICE) << "Overall weight = " << wght;

    // add to arrays
    weights[iev - nfirst] = wght;
  
    mcrec->Clear();
  } // events

  // Close event file
  file.Close();

  //
  // Save weights 
  //

  // Make an output tree for saving the weights.
  TFile * wght_file = new TFile(gOptOutFilename.c_str(), "RECREATE");
  TTree * wght_tree = new TTree("ZExpCCQE","GENIE weights tree");
  // objects to pass elements into tree
  int branch_eventnum = 0;
  float branch_weight = 1.;
  float branch_norm_val = gOptNormValue;
  float branch_zparam_val[MAX_COEF] = {0.};

  wght_tree->Branch("eventnum", &branch_eventnum);
  wght_tree->Branch("weights",  &branch_weight);
  wght_tree->Branch("norm", &branch_norm_val);

  // create and add branches for each z-expansion coefficient
  ostringstream zparam_brnch_name;
  for (int ipr = 0; ipr < n_params; ipr++) {
    zparam_brnch_name.str("");
    zparam_brnch_name << "param_" << ipr+1;
    LOG("rwghtzexpaxff", pINFO) << "Branch name = " << zparam_brnch_name.str();
    branch_zparam_val[ipr] = gOptParameterValue[ipr];
    LOG("rwghtzexpaxff", pINFO) << "Setting parameter value = " << gOptParameterValue[ipr];
    wght_tree->Branch(zparam_brnch_name.str().c_str(), &branch_zparam_val[ipr]);
  }

  ostringstream str_wght;
  for(int iev = nfirst; iev <= nlast; iev++) {
    branch_eventnum = iev;

    // printout
    LOG("grwghtzexpaxff", pNOTICE)
       << "Filling tree with wght = " << weights[iev - nfirst];
    branch_weight = weights[iev - nfirst];
    wght_tree->Fill();
  }

  wght_file->cd();
  wght_tree->Write();
  wght_tree = 0;
  wght_file->Close();

  // free memory
  delete wght_tree;
  delete wght_file;

  LOG("rwghtzexpaxff", pNOTICE)  << "Done!";
  return 0;
}

void PrintSyntax

(

void

)

Variable Documentation

bool gOptDoNorm = false

Definition at line 90 of file gRwghtZExpDirect.cxx.

string gOptInpFilename

Definition at line 86 of file gRwghtZExpDirect.cxx.

Long64_t gOptNEvt1

Definition at line 88 of file gRwghtZExpDirect.cxx.

Long64_t gOptNEvt2

Definition at line 89 of file gRwghtZExpDirect.cxx.

double gOptNormValue = 1.

Definition at line 91 of file gRwghtZExpDirect.cxx.

Referenced by main().

string gOptOutFilename

Definition at line 87 of file gRwghtZExpDirect.cxx.

double gOptParameterValue[MAX_COEF] = {0.}

Definition at line 92 of file gRwghtZExpDirect.cxx.

Referenced by main().