GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Macros | Functions | Variables
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().

548 {
549  int num_pts = 1;
550  for (int i=0;i<kmaxinc;i++)
551  {
552  if (ntwk[i] > 1) num_pts *= ntwk[i];
553  }
554  if (donorm)
555  {
556  if (normtwk > 1) num_pts *= normtwk;
557  }
558  return num_pts;
559 }
GSyst_t GetZExpSystematic ( int  ip)

Definition at line 561 of file gRwghtZExpAxFF.cxx.

References LOG, and pFATAL.

Referenced by IncrementCoefficients(), and main().

562 {
563  switch(ip){
564  case 1: return kXSecTwkDial_ZExpA1CCQE; break;
565  case 2: return kXSecTwkDial_ZExpA2CCQE; break;
566  case 3: return kXSecTwkDial_ZExpA3CCQE; break;
567  case 4: return kXSecTwkDial_ZExpA4CCQE; break;
568  default:
569  LOG("rwghtzexpaxff", pFATAL)
570  << "Cannot find systematic corresponding to parameter " << ip;
571  exit(0);
572  break;
573  }
574 }
#define pFATAL
Definition: Messenger.h:56
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
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.

485 {
486  if (kmaxinc < 2 && ! donorm)
487  {
488  LOG("grwghtzexpaxff",pERROR) << "No coefficients to increment";
489  return false;
490  } else {
491 
492  int ip = -1;
493  bool stopflag = false;
494  GSyst_t gsyst = kXSecTwkDial_ZNormCCQE;
495  do
496  {
497  if (ip > 0 || (ip == 0 && donorm))
498  { // fully incremented a coefficient; reset it and continue to the next
499  if (ip == 0) { twkvals[0 ] = (normtwk > 1 ? -1. : 0.); }
500  else { twkvals[ip] = (ntwk[ip-1] > 1 ? -1. : 0.); }
501 
502  // set the value manually
503  syst.Set(gsyst, twkvals[ip]);
504  LOG("rwghtzexpaxff", pNOTICE) << "Setting z-expansion tweak for param "
505  <<ip<<" : " << twkvals[ip];
506  }
507  stopflag = true;
508 
509  ip++; // increment index
510  if (ip == kmaxinc) { return false; } // done with incrementing
511  if (ip == 0 && ! donorm)
512  {
513  stopflag = false;
514  continue; // skip when not doing norm
515  }
516  // set to next systematic
517  if (ip == 0) { gsyst = kXSecTwkDial_ZNormCCQE; }
518  else { gsyst = GetZExpSystematic(ip); }
519 
520  // increment systematic
521  if (ip == 0)
522  {
523  if (normtwk > 1) { twkvals[0] += 2./float(normtwk-1); }
524  else { stopflag = false; continue; }
525  }
526  else
527  {
528  if (ntwk[ip-1] > 1) { twkvals[ip] += 2./float(ntwk[ip-1]-1); }
529  else { stopflag = false; continue; }
530  }
531 
532  // set the systematic to the new tweak value
533  // actual updating of this will be handled by GReWeight::Reconfigure
534  syst.Set(gsyst, twkvals[ip]);
535  LOG("rwghtzexpaxff", pNOTICE) << "Setting z-expansion tweak for param "
536  <<ip<<" : " << twkvals[ip];
537  if (twkvals[ip] > 1. + controls::kASmallNum) { stopflag=false; } // went over
538 
539  } while (! stopflag); // loop
540 
541  return true;
542  } // if kmaxinc >= 1
543 
544  return false;
545 }
GSyst_t GetZExpSystematic(int ip)
#define pERROR
Definition: Messenger.h:59
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
static const double kASmallNum
Definition: Controls.h:40
#define pNOTICE
Definition: Messenger.h:61
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.

109 {
110  GetCommandLineArgs (argc, argv);
111 
112  // open the ROOT file and get the TTree & its header
113  TTree * tree = 0;
114  NtpMCTreeHeader * thdr = 0;
115  TFile file(gOptInpFilename.c_str(),"READ");
116  tree = dynamic_cast <TTree *> ( file.Get("gtree") );
117  thdr = dynamic_cast <NtpMCTreeHeader *> ( file.Get("header") );
118  LOG("rwghtzexpaxff", pNOTICE) << "Input tree header: " << *thdr;
119  if(!tree){
120  LOG("grwghtzexpaxff", pFATAL)
121  << "Can't find a GHEP tree in input file: "<< file.GetName();
122  gAbortingInErr = true;
123  PrintSyntax();
124  exit(1);
125  }
126  NtpMCEventRecord * mcrec = 0;
127  tree->SetBranchAddress("gmcrec", &mcrec);
128 
129  Long64_t nev_in_file = tree->GetEntries();
130  Long64_t nfirst = 0;
131  Long64_t nlast = 0;
132  GetEventRange(nev_in_file, nfirst, nlast);
133  int nev = int(nlast - nfirst + 1);
134 
135  LOG("rwghtzexpaxff", pNOTICE) << "Will process " << nev << " events";
136 
137  //
138  // Create a GReWeight object and add to it a set of
139  // weight calculators
140  //
141  // If seg-faulting here, need to change
142  // AxialFormFactorModel in UserPhysicsOptions.xml and LwlynSmithFFCC.xml
143  //
144 
145  GReWeight rw;
146  rw.AdoptWghtCalc( "xsec_ccqe", new GReWeightNuXSecCCQE );
147 
148  //
149  // Create a list of systematic params (more to be found at GSyst.h)
150  // set non-default values and re-configure.
151  // Weight calculators included above must be able to handle the tweaked params.
152  // Each tweaking dial t modifies a physics parameter p as:
153  // p_{tweaked} = p_{default} ( 1 + t * dp/p )
154  // So setting a tweaking dial to +/-1 modifies a physics quantity
155  // by +/- 1sigma.
156  // Default fractional errors are defined in GSystUncertainty
157  // and can be overriden.
158  //
159 
160  GSystSet & syst = rw.Systematics();
161 
162  // Create a concrete weight calculator to fine-tune
163  GReWeightNuXSecCCQE * rwccqe =
164  dynamic_cast<GReWeightNuXSecCCQE *> (rw.WghtCalc("xsec_ccqe"));
165  rwccqe->SetMode(GReWeightNuXSecCCQE::kModeZExp);
166  // In case uncertainties need to be altered
167  GSystUncertainty * unc = GSystUncertainty::Instance();
168 
169  // further optional fine-tuning
170  //rwccqe -> RewNue (false);
171  //rwccqe -> RewNuebar (false);
172  //rwccqe -> RewNumubar(false);
173 
174  // Declare the weights and twkvals arrays
175  const int n_events = (const int) nev;
176  const int n_params = (const int) (gOptKmaxInc + 1); // +1 for norm
177  const int n_points =
179  // if segfaulting here, may need to increase MAX_COEF
180  // copied from gtestRwght... seems inefficient
181  // -- couldn't we just load straight to the tree? doing so would prevent segfaults
182  float weights [n_events][n_points];
183  float twkvals [n_points][n_params];
184 
185  // Initialize
186  for (int ipt = 0; ipt < n_points; ipt++)
187  {
188  for (int iev = 0; iev < nev; iev++) { weights[iev][ipt] = 1.; }
189  twkvals[ipt][0] = (gOptDoNorm && (gOptNormTweaks > 1)) ? -1 : 0;
190  for (int ipr = 1; ipr < n_params; ipr++)
191  {
192  twkvals[ipt][ipr] = (gOptNTweaks[ipr-1] > 1 ? -1 : 0);
193  }
194  }
195 
196  // set first values for weighting
197  if (gOptDoNorm)
198  {
199  syst.Set(kXSecTwkDial_ZNormCCQE, twkvals[0][0]);
200  LOG("rwghtzexpaxff", pNOTICE) << "Setting z-expansion tweak for norm : "
201  << twkvals[0][0];
202  }
203  GSyst_t gsyst;
204  for (int ipr = 1; ipr < n_params; ipr++)
205  {
206  gsyst = GetZExpSystematic(ipr);
207  syst.Set(gsyst, twkvals[0][ipr]);
208  LOG("rwghtzexpaxff", pNOTICE) << "Setting z-expansion tweak for param "
209  <<ipr<<" : " << twkvals[0][ipr];
210  if (gOptSigmaDefined)
211  {
212  unc->SetUncertainty(gsyst,gOptSigMin[ipr-1],gOptSigMax[ipr-1]);
213  LOG("rwghtzexpaxff", pNOTICE) << "Setting z-expansion sigma for param "
214  <<ipr<<" : " << gOptSigMin[ipr-1] <<","<< gOptSigMax[ipr-1];
215  }
216  }
217 
218  // point loop (number of parameter combinations)
219  for (int ipt = 0; ipt < n_points; ipt++) {
220  rw.Reconfigure();
221  // Event loop
222  for(int iev = nfirst; iev <= nlast; iev++) {
223  tree->GetEntry(iev);
224 
225  EventRecord & event = *(mcrec->event);
226  LOG("rwghtzexpaxff", pNOTICE) << "Event number = " << iev;
227  LOG("rwghtzexpaxff", pNOTICE) << event;
228 
229  double wght = rw.CalcWeight(event);
230 
231  LOG("rwghtzexpaxff", pNOTICE) << "Overall weight = " << wght;
232 
233  // add to arrays
234  weights[iev - nfirst][ipt] = wght;
235 
236  mcrec->Clear();
237  } // events
238 
239  // set the next set of coefficients to match the previous, then increment
240  if (ipt < n_points-1) {
241  for (int ipr=0;ipr<n_params;ipr++)
242  {
243  twkvals[ipt+1][ipr] = twkvals[ipt][ipr];
244  }
246  twkvals[ipt+1],syst);
247  }
248  } // points
249 
250  // Close event file
251  file.Close();
252 
253  //
254  // Save weights
255  //
256 
257  // Make an output tree for saving the weights.
258  TFile * wght_file = new TFile(gOptOutFilename.c_str(), "RECREATE");
259  TTree * wght_tree = new TTree("ZExpCCQE","GENIE weights tree");
260  // objects to pass elements into tree
261  int branch_eventnum = 0;
262  TArrayF * branch_weight_array = new TArrayF(n_points);
263  TArrayF ** branch_twkdials_array = new TArrayF* [n_params];
264 
265  wght_tree->Branch("eventnum", &branch_eventnum);
266  wght_tree->Branch("weights", &branch_weight_array);
267 
268  // create and add branches for each z-expansion coefficient
269  ostringstream twk_dial_brnch_name;
270  for (int ipr = 0; ipr < n_params; ipr++) {
271  if (!gOptDoNorm && ipr == 0) { continue; } // skip norm if not requested
272  twk_dial_brnch_name.str("");
273  if (ipr == 0) { twk_dial_brnch_name << "twk_dial_param_norm"; }
274  else { twk_dial_brnch_name << "twk_dial_param_" << ipr; }
275  LOG("rwghtzexpaxff", pWARN) << "Branch name = " << twk_dial_brnch_name.str();
276  branch_twkdials_array[ipr] = new TArrayF(n_points);
277  wght_tree->Branch(twk_dial_brnch_name.str().c_str(), branch_twkdials_array[ipr]);
278  }
279 
280  ostringstream str_wght;
281  for(int iev = nfirst; iev <= nlast; iev++) {
282  branch_eventnum = iev;
283 
284  for(int ipt = 0; ipt < n_points; ipt++){
285 
286  // printout
287  str_wght.str("");
288  str_wght << ", tweaked parameter values : ";
289  for (int ipr = 0; ipr < n_params; ipr++) {
290  if (ipr > 0) str_wght << ", ";
291  str_wght << ipr << " -> " << twkvals[ipt][ipr];
292  }
293  LOG("grwghtzexpaxff", pNOTICE)
294  << "Filling tree with wght = " << weights[iev - nfirst][ipt] << str_wght.str();
295 
296  // fill tree
297  branch_weight_array -> AddAt (weights [iev - nfirst][ipt], ipt);
298  for (int ipr = (gOptDoNorm ? 0 : 1); ipr < n_params; ipr++) // skip norm if not requested
299  { branch_twkdials_array[ipr] -> AddAt (twkvals[ipt][ipr], ipt); }
300 
301  } // twk_dial loop
302  wght_tree->Fill();
303  }
304 
305  wght_file->cd();
306  wght_tree->Write();
307  wght_tree = 0;
308  wght_file->Close();
309 
310  // free memory
311  delete wght_tree;
312  delete wght_file;
313  for (int ipr = 0; ipr < n_params; ipr++) {
314  if (!gOptDoNorm && ipr == 0) { continue; }
315  delete branch_twkdials_array[ipr];
316  }
317  delete branch_twkdials_array;
318  delete branch_weight_array;
319 
320  LOG("rwghtzexpaxff", pNOTICE) << "Done!";
321  return 0;
322 }
int gOptNormTweaks
int gOptKmaxInc
GSyst_t GetZExpSystematic(int ip)
float gOptSigMax[MAX_COEF]
bool gOptDoNorm
MINOS-style ntuple record. Each such ntuple record holds a generated EventRecord object. Ntuples of this type are intended for feeding GENIE events into other applications (for example the GEANT4 based MC generation framework of an experiment) if no direct interface exists.
#define pFATAL
Definition: Messenger.h:56
float gOptSigMin[MAX_COEF]
int gOptNTweaks[MAX_COEF]
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
MINOS-style Ntuple Class to hold an output MC Tree Header.
int GetNumberOfWeights(int *ntwk, int kmaxinc, int normtwk, bool donorm)
bool gOptSigmaDefined
void GetEventRange(Long64_t nev, Long64_t &n1, Long64_t &n2)
Definition: gEvDump.cxx:217
#define pWARN
Definition: Messenger.h:60
string gOptInpFilename
Definition: gEvDump.cxx:76
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
Definition: EventRecord.h:37
#define pNOTICE
Definition: Messenger.h:61
void GetCommandLineArgs(int argc, char **argv)
Definition: gAtmoEvGen.cxx:563
void Clear(Option_t *opt="")
bool gAbortingInErr
Definition: Messenger.cxx:34
bool IncrementCoefficients(double *coefmin, double *coefmax, double *coefinc, int kmaxinc, tdata params, AlgFactory *algf, Registry *r, const Registry *gc)
void PrintSyntax(void)
EventRecord * event
event
string gOptOutFilename
Definition: gEvScan.cxx:93
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().