#include <string>
#include <sstream>
#include <TSystem.h>
#include <TFile.h>
#include <TTree.h>
#include <TArrayF.h>
#include "Algorithm/AlgConfigPool.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"

Macros
#define	MAX_COEF 4
	A simple program to illustrate how to use the GENIE event reweighting for use with the z-expansion axial form factor. More...

Functions
void	PrintSyntax ()

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

void	GetCommandLineArgs (int argc, char **argv)

int	GetNumberOfWeights (int *ntwk, int kmaxinc, int normtwk, bool donorm)

bool	IncrementCoefficients (int ntwk, int kmaxinc, int normtwk, bool donorm, float twkvals, GSystSet &syst)

GSyst_t	GetZExpSystematic (int ip)

int	main (int argc, char **argv)

Variables
string	gOptInpFilename

string	gOptOutFilename

Long64_t	gOptNEvt1

Long64_t	gOptNEvt2

int	gOptKmaxInc = 0

int	gOptNormTweaks = 0

bool	gOptDoNorm = false

bool	gOptSigmaDefined = false

int	gOptNTweaks [MAX_COEF] = {0 }

float	gOptSigMin [MAX_COEF] = {0.}

float	gOptSigMax [MAX_COEF] = {0.}

Macro Definition Documentation

#define MAX_COEF 4

A simple program to illustrate how to use the GENIE event reweighting for use with the z-expansion axial form factor.

gRwghtZExpAxFF

grwghtzexpaxff -f filename -t NTwk1,NTwk2,... [-n nev] [-o fileOutName] [-s SigmaLo1,SigmaHi1,SigmaLo2,SigmaHi2,...] [-m NTwkN]

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) -t specify number of tweaks on each z-expansion coefficient values are comma separated (# < 2 are ignored) -s specify +- one-sigma bounds on all coefficients up to max values are comma separated, given as percentages requires 2x number of fields from -t option default value is 10% on all coefficients -m number of tweaks on normalization puts reweighting into norm+shape mode

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

based on gtestRewght by

Costas Andreopoulos <constantinos.andreopoulos cern.ch> STFC, Rutherford Appleton Laboratory

Dec 26, 2014

Definition at line 79 of file gRwghtZExpAxFF.cxx.

Function Documentation

void GetCommandLineArgs	(	int	argc,
		char **	argv
	)

Definition at line 324 of file gRwghtZExpAxFF.cxx.

 {
   LOG("rwghtzexpaxff", pINFO) << "*** Parsing command line arguments";
 
   CmdLnArgParser parser(argc,argv);
 
   // get GENIE event sample
   if( parser.OptionExists('f') ) {  
     LOG("rwghtzexpaxff", pINFO) << "Reading event sample filename";
     gOptInpFilename = parser.ArgAsString('f');
   } else {
     LOG("rwghtzexpaxff", pFATAL) 
       << "Unspecified input filename - Exiting";
     PrintSyntax();
     exit(1);
   }
 
   // output weight file
   if(parser.OptionExists('o')) {
     LOG("rwghtzexpaxff", pINFO) << "Reading requested output filename";
     gOptOutFilename = parser.ArgAsString('o');
   } else {
     LOG("rwghtzexpaxff", pINFO) << "Setting default output filename";
     gOptOutFilename = "test_rw_zexp_axff.root";
   }
 
   if ( parser.OptionExists('n') ) {
     //
     LOG("grwghtzexpaxff", pINFO) << "Reading number of events to analyze";
     string nev =  parser.ArgAsString('n');
     if (nev.find(",") != string::npos) {
       vector<long> vecn = parser.ArgAsLongTokens('n',",");
       if(vecn.size()!=2) {
          LOG("grwghtzexpaxff", pFATAL) << "Invalid syntax";
          gAbortingInErr = true;
          PrintSyntax();
          exit(1);
       }
       // User specified a comma-separated set of values n1,n2.
       // Use [n1,n2] as the event range to process.
       gOptNEvt1 = vecn[0];
       gOptNEvt2 = vecn[1];
     } else {
       // User specified a single number n.
       // Use [0,n] as the event range to process.
       gOptNEvt1 = -1;
       gOptNEvt2 = parser.ArgAsLong('n');
     }
   } else {
     LOG("grwghtzexpaxff", pINFO)
       << "Unspecified number of events to analyze - Use all";
     gOptNEvt1 = -1;
     gOptNEvt2 = -1;
   }
   LOG("grwghtzexpaxff", pDEBUG)
     << "Input event range: " << gOptNEvt1 << ", " << gOptNEvt2;
 
   // number of tweaks:
   if( parser.OptionExists('t') ) {
     LOG("rwghtzexpaxff", pINFO) << "Reading number of tweaks";
     string coef = parser.ArgAsString('t');
     
     // split into sections of min,max,inc(rement)
     vector<string> coefrange = utils::str::Split(coef, ",");
     gOptKmaxInc = coefrange.size();
     if (gOptKmaxInc > MAX_COEF)
     {
       LOG("rwghtzexpaxff", pFATAL) 
         << "Too many coefficients: increase MAX_COEF in source code and recompile";
       exit(1);
     }
     //AlgConfigPool * confp = AlgConfigPool::Instance();
     //const Registry * gc = confp->GlobalParameterList();
     //int kDef = fConfig->GetIntDef("QEL-Kmax",gc->GetInt("QEL-Kmax"));
     //else if (gOptKmaxInc > kDef)
     //{
     //  LOG("rwghtzexpaxff", pFATAL) 
     //    << "Too many coefficients:
     //        requested number of coefficients is more than defined in UserPhysicsOptions.xml";
     //  exit(1);
     //}
 
     LOG("rwghtzexpaxff", pINFO) << "Largest coefficient to tweak : " << gOptKmaxInc;
     for (int ik = 0;ik<gOptKmaxInc;ik++)
     {
       gOptNTweaks[ik] = atof(coefrange[ik].c_str());
     }
     for (int ik = 0;ik<gOptKmaxInc;ik++)
     {
       LOG("rwghtzexpaxff",pINFO)<<"Number of tweaks on coefficient "<<ik+1<<" : "<< gOptNTweaks[ik];
     }
   } else {
     LOG("rwghtzexpaxff", pFATAL) 
       << "Unspecified tweaks for parameters - Exiting";
     PrintSyntax();
     exit(1);
   }
 
   // lower/upper sigma:
   if( parser.OptionExists('s') ) {
     LOG("rwghtzexpaxff", pINFO) << "Reading specified parameter uncertainties";
     string coef = parser.ArgAsString('s');
     
     // split into sections of min,max
     vector<string> sigrange = utils::str::Split(coef, ",");
     // gOptKmaxInc defined by number of tweaks (-t)
     assert(sigrange.size() == (unsigned int) 2*gOptKmaxInc);
     gOptSigmaDefined = true;
     for (int ik = 0;ik<gOptKmaxInc;ik++)
     {
       gOptSigMin[ik] = atof(sigrange[ik*2  ].c_str());
       gOptSigMax[ik] = atof(sigrange[ik*2+1].c_str());
     }
     for (int ik = 0;ik<gOptKmaxInc;ik++)
     {
       LOG("rwghtzexpaxff",pINFO)<<ik+1<<": "<< gOptSigMin[ik] <<","<< gOptSigMax[ik];
     }
   }
 
   // number of norm tweaks:
   if( parser.OptionExists('m') ) {
     LOG("rwghtzexpaxff", pINFO) << "Reading number of tweaks on normalization";
     string coef = parser.ArgAsString('m');
     gOptDoNorm = true;
     gOptNormTweaks   = atof(coef.c_str());
     LOG("rwghtzexpaxff",pINFO)<<"Number of tweaks on norm : "<< gOptNormTweaks;
   }
 }

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

Definition at line 453 of file gRwghtZExpAxFF.cxx.

 {
   nfirst = 0;
   nlast  = 0;
 
   if(gOptNEvt1>=0 && gOptNEvt2>=0) {
     // Input was `-n N1,N2'.
     // Process events [N1,N2].
     // Note: Incuding N1 and N2.
     nfirst = gOptNEvt1;
     nlast  = TMath::Min(nev_in_file-1, gOptNEvt2);
   }
   else
   if(gOptNEvt1<0 && gOptNEvt2>=0) {
     // Input was `-n N'.
     // Process first N events [0,N). 
     // Note: Event N is not included.
     nfirst = 0;
     nlast  = TMath::Min(nev_in_file-1, gOptNEvt2-1);
   }
   else
   if(gOptNEvt1<0 && gOptNEvt2<0) {
     // No input. Process all events.
     nfirst = 0;
     nlast  = nev_in_file-1;
   }
 
   assert(nfirst <= nlast && nfirst >= 0 && nlast <= nev_in_file-1);
 }

int GetNumberOfWeights	(	int *	ntwk,
		int	kmaxinc,
		int	normtwk,
		bool	donorm
	)

Definition at line 547 of file gRwghtZExpAxFF.cxx.

 {
   int  num_pts = 1;
   for (int i=0;i<kmaxinc;i++)
   {
     if (ntwk[i] > 1) num_pts *= ntwk[i];
   }
   if (donorm)
   {
     if (normtwk > 1) num_pts *= normtwk;
   }
   return num_pts;
 }

GSyst_t GetZExpSystematic ( int ip )

Definition at line 561 of file gRwghtZExpAxFF.cxx.

 {
     switch(ip){
       case 1: return kXSecTwkDial_ZExpA1CCQE; break;
       case 2: return kXSecTwkDial_ZExpA2CCQE; break;
       case 3: return kXSecTwkDial_ZExpA3CCQE; break;
       case 4: return kXSecTwkDial_ZExpA4CCQE; break;
       default:
         LOG("rwghtzexpaxff", pFATAL) 
           << "Cannot find systematic corresponding to parameter " << ip;
         exit(0);
         break;
     }
 }

bool IncrementCoefficients	(	int *	ntwk,
		int	kmaxinc,
		int	normtwk,
		bool	donorm,
		float *	twkvals,
		GSystSet &	syst
	)

Definition at line 483 of file gRwghtZExpAxFF.cxx.

 {
   if (kmaxinc < 2 && ! donorm)
   {
     LOG("grwghtzexpaxff",pERROR) << "No coefficients to increment";
     return false;
   } else {
 
   int ip = -1;
   bool stopflag = false;
   GSyst_t gsyst = kXSecTwkDial_ZNormCCQE;
   do
   {
     if (ip > 0 || (ip == 0 && donorm))
     { // fully incremented a coefficient; reset it and continue to the next
       if (ip == 0) { twkvals[0 ] = (normtwk    > 1 ? -1. : 0.); }
       else         { twkvals[ip] = (ntwk[ip-1] > 1 ? -1. : 0.); }
 
       // set the value manually
       syst.Set(gsyst, twkvals[ip]);
       LOG("rwghtzexpaxff", pNOTICE) << "Setting z-expansion tweak for param " 
         <<ip<<" : " << twkvals[ip];
     }
     stopflag = true;
 
     ip++;                                 // increment index
     if (ip == kmaxinc) { return false; }  // done with incrementing
     if (ip == 0 && ! donorm)
     {
       stopflag = false;
       continue;  // skip when not doing norm
     }
     // set to next systematic
     if (ip == 0) { gsyst = kXSecTwkDial_ZNormCCQE; }
     else         { gsyst = GetZExpSystematic(ip); }
 
     // increment systematic
     if (ip == 0)
     {
       if (normtwk > 1) { twkvals[0] += 2./float(normtwk-1); }
       else             { stopflag = false; continue; }
     }
     else
     {
       if (ntwk[ip-1] > 1) { twkvals[ip] += 2./float(ntwk[ip-1]-1); }
       else                { stopflag = false; continue; }
     }
 
     // set the systematic to the new tweak value
     // actual updating of this will be handled by GReWeight::Reconfigure
     syst.Set(gsyst, twkvals[ip]);
     LOG("rwghtzexpaxff", pNOTICE) << "Setting z-expansion tweak for param " 
       <<ip<<" : " << twkvals[ip];
     if (twkvals[ip] > 1. + controls::kASmallNum) { stopflag=false; } // went over
 
   } while (! stopflag); // loop
 
   return true;
   } // if kmaxinc >= 1
 
   return false;
 }

int main	(	int	argc,
		char **	argv
	)

Definition at line 108 of file gRwghtZExpAxFF.cxx.

 {
   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();
 
   // 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) (gOptKmaxInc + 1); // +1 for norm
   const int n_points =
    (const int) GetNumberOfWeights(gOptNTweaks,gOptKmaxInc,gOptNormTweaks,gOptDoNorm);
   // if segfaulting here, may need to increase MAX_COEF
   // copied from gtestRwght... seems inefficient
   // -- couldn't we just load straight to the tree? doing so would prevent segfaults
   float weights  [n_events][n_points];
   float twkvals  [n_points][n_params];
 
   // Initialize
   for (int ipt = 0; ipt < n_points; ipt++)
   {
     for (int iev = 0; iev < nev; iev++) { weights[iev][ipt] = 1.; }
     twkvals[ipt][0] = (gOptDoNorm && (gOptNormTweaks > 1)) ? -1 : 0;
     for (int ipr = 1; ipr < n_params; ipr++)
     {
       twkvals[ipt][ipr] = (gOptNTweaks[ipr-1] > 1 ? -1 : 0);
     }
   }
 
   // set first values for weighting
   if (gOptDoNorm)
   {
     syst.Set(kXSecTwkDial_ZNormCCQE, twkvals[0][0]);
     LOG("rwghtzexpaxff", pNOTICE) << "Setting z-expansion tweak for norm : "
       << twkvals[0][0];
   }
   GSyst_t gsyst;
   for (int ipr = 1; ipr < n_params; ipr++)
   {
     gsyst = GetZExpSystematic(ipr);
     syst.Set(gsyst, twkvals[0][ipr]);
     LOG("rwghtzexpaxff", pNOTICE) << "Setting z-expansion tweak for param " 
       <<ipr<<" : " << twkvals[0][ipr];
     if (gOptSigmaDefined)
     {
       unc->SetUncertainty(gsyst,gOptSigMin[ipr-1],gOptSigMax[ipr-1]);
       LOG("rwghtzexpaxff", pNOTICE) << "Setting z-expansion sigma for param " 
         <<ipr<<" : " << gOptSigMin[ipr-1] <<","<< gOptSigMax[ipr-1];
     }
   }
 
   // point loop (number of parameter combinations)
   for (int ipt = 0; ipt < n_points; ipt++) {
     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][ipt] = wght;
   
       mcrec->Clear();
     } // events
 
     // set the next set of coefficients to match the previous, then increment
     if (ipt < n_points-1) {
       for (int ipr=0;ipr<n_params;ipr++)
       {
         twkvals[ipt+1][ipr] = twkvals[ipt][ipr];
       }
       IncrementCoefficients(gOptNTweaks,n_params,gOptNormTweaks,gOptDoNorm,
         twkvals[ipt+1],syst);
     }
   }   // points
 
   // 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;
   TArrayF *  branch_weight_array   = new TArrayF(n_points);
   TArrayF ** branch_twkdials_array = new TArrayF* [n_params];
 
   wght_tree->Branch("eventnum", &branch_eventnum);
   wght_tree->Branch("weights",  &branch_weight_array);
 
   // create and add branches for each z-expansion coefficient
   ostringstream twk_dial_brnch_name;
   for (int ipr = 0; ipr < n_params; ipr++) {
     if (!gOptDoNorm && ipr == 0) { continue; } // skip norm if not requested
     twk_dial_brnch_name.str("");
     if (ipr == 0) { twk_dial_brnch_name << "twk_dial_param_norm";      }
     else          { twk_dial_brnch_name << "twk_dial_param_" << ipr; }
     LOG("rwghtzexpaxff", pWARN) << "Branch name = " << twk_dial_brnch_name.str();
     branch_twkdials_array[ipr] = new TArrayF(n_points);
     wght_tree->Branch(twk_dial_brnch_name.str().c_str(), branch_twkdials_array[ipr]);
   }
 
   ostringstream str_wght;
   for(int iev = nfirst; iev <= nlast; iev++) {
     branch_eventnum = iev;
 
     for(int ipt = 0; ipt < n_points; ipt++){
 
        // printout
        str_wght.str("");
        str_wght << ", tweaked parameter values : ";
        for (int ipr = 0; ipr < n_params; ipr++) {
           if (ipr > 0) str_wght << ", ";
           str_wght << ipr << " -> " << twkvals[ipt][ipr];
        }
        LOG("grwghtzexpaxff", pNOTICE)
           << "Filling tree with wght = " << weights[iev - nfirst][ipt] << str_wght.str();
 
        // fill tree
        branch_weight_array   -> AddAt (weights [iev - nfirst][ipt], ipt);
        for (int ipr = (gOptDoNorm ? 0 : 1); ipr < n_params; ipr++) // skip norm if not requested
          { branch_twkdials_array[ipr] -> AddAt (twkvals[ipt][ipr], ipt); }
 
     } // twk_dial loop
     wght_tree->Fill();
   }
 
   wght_file->cd();
   wght_tree->Write();
   wght_tree = 0;
   wght_file->Close();
 
   // free memory
   delete wght_tree;
   delete wght_file;
   for (int ipr = 0; ipr < n_params; ipr++) {
     if (!gOptDoNorm && ipr == 0) { continue; }
     delete branch_twkdials_array[ipr];
   }
   delete branch_twkdials_array;
   delete branch_weight_array;
 
   LOG("rwghtzexpaxff", pNOTICE)  << "Done!";
   return 0;
 }

void PrintSyntax ( void )

Definition at line 576 of file gRwghtZExpAxFF.cxx.

 {
   LOG("grwghtzexpaxff", pFATAL)
      << "\n\n"
      << "grwghtzexpaxff               \n"
      << "     -f input_event_file     \n"
      << "     -t ntwk1[,ntwk2[,...]]  \n"
      << "    [-n nev]                 \n"
      << "    [-s sigLo1,sigHi1[,...]] \n"
      << "    [-o output_weights_file] \n"
      << "    [-m ntwkNorm]" ;
 }

Variable Documentation

bool gOptDoNorm = false

Definition at line 101 of file gRwghtZExpAxFF.cxx.

string gOptInpFilename

Definition at line 94 of file gRwghtZExpAxFF.cxx.

int gOptKmaxInc = 0

Definition at line 99 of file gRwghtZExpAxFF.cxx.

Long64_t gOptNEvt1

Definition at line 97 of file gRwghtZExpAxFF.cxx.

Long64_t gOptNEvt2

Definition at line 98 of file gRwghtZExpAxFF.cxx.

int gOptNormTweaks = 0

Definition at line 100 of file gRwghtZExpAxFF.cxx.

int gOptNTweaks[MAX_COEF] = {0 }

Definition at line 103 of file gRwghtZExpAxFF.cxx.

string gOptOutFilename

Definition at line 95 of file gRwghtZExpAxFF.cxx.

bool gOptSigmaDefined = false

Definition at line 102 of file gRwghtZExpAxFF.cxx.

float gOptSigMax[MAX_COEF] = {0.}

Definition at line 105 of file gRwghtZExpAxFF.cxx.

float gOptSigMin[MAX_COEF] = {0.}

Definition at line 104 of file gRwghtZExpAxFF.cxx.

Macros

Functions

Variables

Macro Definition Documentation

Function Documentation

Variable Documentation