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

96 {
97  GetCommandLineArgs (argc, argv);
98 
99  // open the ROOT file and get the TTree & its header
100  TTree * tree = 0;
101  NtpMCTreeHeader * thdr = 0;
102  TFile file(gOptInpFilename.c_str(),"READ");
103  tree = dynamic_cast <TTree *> ( file.Get("gtree") );
104  thdr = dynamic_cast <NtpMCTreeHeader *> ( file.Get("header") );
105  LOG("rwghtzexpaxff", pNOTICE) << "Input tree header: " << *thdr;
106  if(!tree){
107  LOG("grwghtzexpaxff", pFATAL)
108  << "Can't find a GHEP tree in input file: "<< file.GetName();
109  gAbortingInErr = true;
110  PrintSyntax();
111  exit(1);
112  }
113  NtpMCEventRecord * mcrec = 0;
114  tree->SetBranchAddress("gmcrec", &mcrec);
115 
116  Long64_t nev_in_file = tree->GetEntries();
117  Long64_t nfirst = 0;
118  Long64_t nlast = 0;
119  GetEventRange(nev_in_file, nfirst, nlast);
120  int nev = int(nlast - nfirst + 1);
121 
122  LOG("rwghtzexpaxff", pNOTICE) << "Will process " << nev << " events";
123 
124  //
125  // Create a GReWeight object and add to it a set of
126  // weight calculators
127  //
128  // If seg-faulting here, need to change
129  // AxialFormFactorModel in UserPhysicsOptions.xml and LwlynSmithFFCC.xml
130  //
131 
132  GReWeight rw;
133  rw.AdoptWghtCalc( "xsec_ccqe", new GReWeightNuXSecCCQE );
134 
135  //
136  // Create a list of systematic params (more to be found at GSyst.h)
137  // set non-default values and re-configure.
138  // Weight calculators included above must be able to handle the tweaked params.
139  // Each tweaking dial t modifies a physics parameter p as:
140  // p_{tweaked} = p_{default} ( 1 + t * dp/p )
141  // So setting a tweaking dial to +/-1 modifies a physics quantity
142  // by +/- 1sigma.
143  // Default fractional errors are defined in GSystUncertainty
144  // and can be overriden.
145  //
146 
147  GSystSet & syst = rw.Systematics();
148 
149  // Create a concrete weight calculator to fine-tune
150  GReWeightNuXSecCCQE * rwccqe =
151  dynamic_cast<GReWeightNuXSecCCQE *> (rw.WghtCalc("xsec_ccqe"));
152  rwccqe->SetMode(GReWeightNuXSecCCQE::kModeZExp);
153  // In case uncertainties need to be altered
154  GSystUncertainty * unc = GSystUncertainty::Instance();
155 
156  // Set up algorithm for loading central values
157  AlgFactory * algf = AlgFactory::Instance();
158  AlgId id("genie::LwlynSmithQELCCPXSec","ZExp");
159 
160  Algorithm * alg = algf->AdoptAlgorithm(id);
161  XSecAlgorithmI* fXSecModel = dynamic_cast<XSecAlgorithmI*>(alg);
162  fXSecModel->AdoptSubstructure();
163 
164  Registry * fXSecModelConfig = new Registry(fXSecModel->GetConfig());
165  AlgConfigPool * confp = AlgConfigPool::Instance();
166  const Registry * gc = confp->GlobalParameterList();
167 
168  // further optional fine-tuning
169  //rwccqe -> RewNue (false);
170  //rwccqe -> RewNuebar (false);
171  //rwccqe -> RewNumubar(false);
172 
173  // Declare the weights and twkvals arrays
174  const int n_events = (const int) nev;
175  const int n_params = (const int) MAX_COEF;
176  // if segfaulting here, may need to increase MAX_COEF
177  float weights [n_events];
178 
179  // Initialize
180  for (int iev = 0; iev < nev; iev++) { weights[iev] = 1.; }
181 
182  // set first values for weighting
183  if (gOptDoNorm)
184  {
185  //LOG("rwghtzexpaxff", pNOTICE) << "Setting z-expansion tweak for norm : "
186  // << (gOptNormValue-1.)
187  // /unc->OneSigmaErr(kXSecTwkDial_ZNormCCQE,TMath::Sign(1,(int)(gOptNormValue-1.)));
188  syst.Set(kXSecTwkDial_ZNormCCQE, (gOptNormValue-1.)
189  /unc->OneSigmaErr(kXSecTwkDial_ZNormCCQE,TMath::Sign(1,(int)(gOptNormValue-1.))));
190  }
191  GSyst_t gsyst;
192  double cval = 0.; // central value for parameter
193  ostringstream zval_name; // string to use to extract central value
194  for (int ipr = 0; ipr < n_params; ipr++)
195  {
196  gsyst = GetZExpSystematic(ipr+1); // get tweak dial
197  zval_name.str("");
198  zval_name << "QEL-Z_A" << ipr+1;
199 
200  cval = fXSecModelConfig->GetDoubleDef(zval_name.str(),gc->GetDouble(zval_name.str()));
201  unc->SetUncertainty(gsyst,1.,1.); // easier to deal with
202 
203  //LOG("rwghtzexpaxff", pNOTICE) << "Setting z-expansion tweak for param "
204  // <<ipr<<" : " << (gOptParameterValue[ipr]-cval)/cval;
205  syst.Set(gsyst, (gOptParameterValue[ipr]-cval)/cval);
206  }
207 
208  rw.Reconfigure();
209  // Event loop
210  for(int iev = nfirst; iev <= nlast; iev++) {
211  tree->GetEntry(iev);
212 
213  EventRecord & event = *(mcrec->event);
214  LOG("rwghtzexpaxff", pNOTICE) << "Event number = " << iev;
215  LOG("rwghtzexpaxff", pNOTICE) << event;
216 
217  double wght = rw.CalcWeight(event);
218 
219  LOG("rwghtzexpaxff", pNOTICE) << "Overall weight = " << wght;
220 
221  // add to arrays
222  weights[iev - nfirst] = wght;
223 
224  mcrec->Clear();
225  } // events
226 
227  // Close event file
228  file.Close();
229 
230  //
231  // Save weights
232  //
233 
234  // Make an output tree for saving the weights.
235  TFile * wght_file = new TFile(gOptOutFilename.c_str(), "RECREATE");
236  TTree * wght_tree = new TTree("ZExpCCQE","GENIE weights tree");
237  // objects to pass elements into tree
238  int branch_eventnum = 0;
239  float branch_weight = 1.;
240  float branch_norm_val = gOptNormValue;
241  float branch_zparam_val[MAX_COEF] = {0.};
242 
243  wght_tree->Branch("eventnum", &branch_eventnum);
244  wght_tree->Branch("weights", &branch_weight);
245  wght_tree->Branch("norm", &branch_norm_val);
246 
247  // create and add branches for each z-expansion coefficient
248  ostringstream zparam_brnch_name;
249  for (int ipr = 0; ipr < n_params; ipr++) {
250  zparam_brnch_name.str("");
251  zparam_brnch_name << "param_" << ipr+1;
252  LOG("rwghtzexpaxff", pINFO) << "Branch name = " << zparam_brnch_name.str();
253  branch_zparam_val[ipr] = gOptParameterValue[ipr];
254  LOG("rwghtzexpaxff", pINFO) << "Setting parameter value = " << gOptParameterValue[ipr];
255  wght_tree->Branch(zparam_brnch_name.str().c_str(), &branch_zparam_val[ipr]);
256  }
257 
258  ostringstream str_wght;
259  for(int iev = nfirst; iev <= nlast; iev++) {
260  branch_eventnum = iev;
261 
262  // printout
263  LOG("grwghtzexpaxff", pNOTICE)
264  << "Filling tree with wght = " << weights[iev - nfirst];
265  branch_weight = weights[iev - nfirst];
266  wght_tree->Fill();
267  }
268 
269  wght_file->cd();
270  wght_tree->Write();
271  wght_tree = 0;
272  wght_file->Close();
273 
274  // free memory
275  delete wght_tree;
276  delete wght_file;
277 
278  LOG("rwghtzexpaxff", pNOTICE) << "Done!";
279  return 0;
280 }
Cross Section Calculation Interface.
RgDbl GetDoubleDef(RgKey key, RgDbl def_opt, bool set_def=true)
Definition: Registry.cxx:535
GSyst_t GetZExpSystematic(int ip)
#define MAX_COEF
A simple program to reweight the GENIE z-expansion axial form factor from one set of z-expansion para...
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.
A singleton class holding all configuration registries built while parsing all loaded XML configurati...
Definition: AlgConfigPool.h:40
#define pFATAL
Definition: Messenger.h:56
Algorithm abstract base class.
Definition: Algorithm.h:54
RgDbl GetDouble(RgKey key) const
Definition: Registry.cxx:474
virtual const Registry & GetConfig(void) const
Definition: Algorithm.cxx:246
double gOptParameterValue[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.
#define pINFO
Definition: Messenger.h:62
void AdoptSubstructure(void)
Definition: Algorithm.cxx:400
void GetEventRange(Long64_t nev, Long64_t &n1, Long64_t &n2)
Definition: gEvDump.cxx:217
Algorithm * AdoptAlgorithm(const AlgId &algid) const
Definition: AlgFactory.cxx:116
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
Algorithm ID (algorithm name + configuration set name)
Definition: AlgId.h:34
double gOptNormValue
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
#define pNOTICE
Definition: Messenger.h:61
void GetCommandLineArgs(int argc, char **argv)
Definition: gAtmoEvGen.cxx:563
The GENIE Algorithm Factory.
Definition: AlgFactory.h:39
void Clear(Option_t *opt="")
bool gAbortingInErr
Definition: Messenger.cxx:34
void PrintSyntax(void)
EventRecord * event
event
string gOptOutFilename
Definition: gEvScan.cxx:93
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().