gRwghtZExpAxFF.cxx File Reference

#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"

Include dependency graph for gRwghtZExpAxFF.cxx:

Go to the source code of this file.

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 <c.andreopoulos cern.ch> University of Liverpool

Created:

Dec 26, 2014

License:

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

Definition at line 79 of file gRwghtZExpAxFF.cxx.

Function Documentation

void GetCommandLineArgs	(	int	argc,
		char **	argv
	)

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

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

Definition at line 547 of file gRwghtZExpAxFF.cxx.

Referenced by main().

{
  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.

References LOG, and pFATAL.

Referenced by IncrementCoefficients(), and main().

{
    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.

References GetZExpSystematic(), genie::controls::kASmallNum, LOG, pERROR, and pNOTICE.

{
  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.

References genie::NtpMCEventRecord::Clear(), genie::NtpMCEventRecord::event, genie::gAbortingInErr, GetCommandLineArgs(), GetEventRange(), GetNumberOfWeights(), GetZExpSystematic(), gOptDoNorm, gOptInpFilename, gOptKmaxInc, gOptNormTweaks, gOptNTweaks, gOptOutFilename, gOptSigmaDefined, gOptSigMax, gOptSigMin, IncrementCoefficients(), LOG, pFATAL, pNOTICE, PrintSyntax(), and pWARN.

{
  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

)

Variable Documentation

bool gOptDoNorm = false

Definition at line 101 of file gRwghtZExpAxFF.cxx.

Referenced by main().

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.

Referenced by main().

int gOptNTweaks[MAX_COEF] = {0 }

Definition at line 103 of file gRwghtZExpAxFF.cxx.

Referenced by main().

string gOptOutFilename

Definition at line 95 of file gRwghtZExpAxFF.cxx.

bool gOptSigmaDefined = false

Definition at line 102 of file gRwghtZExpAxFF.cxx.

Referenced by main().

float gOptSigMax[MAX_COEF] = {0.}

Definition at line 105 of file gRwghtZExpAxFF.cxx.

Referenced by main().

float gOptSigMin[MAX_COEF] = {0.}

Definition at line 104 of file gRwghtZExpAxFF.cxx.

Referenced by main().