gRwghtZExpDirect.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3 
4 \program gRwghtZExpDirect
5 
6 \brief A simple program to reweight the GENIE z-expansion axial form factor
7  from one set of z-expansion parameters directly to another
8 
9 \syntax grwghtzexpaxff -f filename -v val1,val2,val3,val4 [-n nev] [-o fileOutName] [-m norm]
10 
11  where
12  [] is an optional argument
13  -f specifies a GENIE event file (GHEP format)
14  -o specifies a GENIE output filename
15  -n specifies the number of events to process (default: all)
16  -v z-expansion values to reweight to
17  -m reweight normalization of form factor (default: 1)
18 
19 \author Aaron Meyer <asmeyer2012 \at uchicago.edu>
20  University of Chicago, Fermi National Accelerator Laboratory
21 
22  based on gtestRewght by
23 
24  Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
25  STFC, Rutherford Appleton Laboratory
26 
27 \created Mar 14, 2016
28 
29 \cpright Copyright (c) 2003-2013, GENIE Neutrino MC Generator Collaboration
30  For the full text of the license visit http://copyright.genie-mc.org
31  or see $GENIE/LICENSE
32 */
33 //____________________________________________________________________________
34 
35 
36 #include <string>
37 #include <sstream>
38 
39 #include <TSystem.h>
40 #include <TFile.h>
41 #include <TTree.h>
42 #include <TArrayF.h>
43 
45 #include "Algorithm/AlgFactory.h"
46 #include "Base/XSecAlgorithmI.h"
47 #include "Conventions/Controls.h"
48 #include "EVGCore/EventRecord.h"
49 #include "Ntuple/NtpMCFormat.h"
50 #include "Ntuple/NtpMCTreeHeader.h"
52 #include "Messenger/Messenger.h"
53 #include "ReWeight/GReWeightI.h"
54 #include "ReWeight/GSystSet.h"
55 #include "ReWeight/GReWeight.h"
56 #include "ReWeight/GReWeightNuXSecCCQE.h"
57 #include "ReWeight/GReWeightNuXSecCCQEvec.h"
58 #include "ReWeight/GReWeightNuXSecCCRES.h"
59 #include "ReWeight/GReWeightNuXSecNCRES.h"
60 #include "ReWeight/GReWeightNuXSecDIS.h"
61 #include "ReWeight/GReWeightNuXSecCOH.h"
62 #include "ReWeight/GReWeightNonResonanceBkg.h"
63 #include "ReWeight/GReWeightFGM.h"
64 #include "ReWeight/GReWeightDISNuclMod.h"
65 #include "ReWeight/GReWeightResonanceDecay.h"
66 #include "ReWeight/GReWeightFZone.h"
67 #include "ReWeight/GReWeightINuke.h"
68 #include "ReWeight/GReWeightAGKY.h"
69 #include "ReWeight/GSystUncertainty.h"
70 #include "Utils/CmdLnArgParser.h"
71 #include "Utils/StringUtils.h"
72 
73 // number of parameter tweaks which are coded into reweighting
74 #define MAX_COEF 4
75 
76 using namespace genie;
77 using namespace genie::rew;
78 using std::string;
79 using std::ostringstream;
80 
81 void PrintSyntax();
82 void GetEventRange (Long64_t nev_in_file, Long64_t & nfirst, Long64_t & nlast);
83 void GetCommandLineArgs (int argc, char ** argv);
84 GSyst_t GetZExpSystematic(int ip);
85 
88 Long64_t gOptNEvt1;
89 Long64_t gOptNEvt2;
90 bool gOptDoNorm = false;
91 double gOptNormValue = 1.;
93 
94 //___________________________________________________________________
95 int main(int argc, char ** argv)
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());
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 }
281 //___________________________________________________________________
282 void GetCommandLineArgs(int argc, char ** argv)
283 {
284  LOG("rwghtzexpaxff", pINFO) << "*** Parsing command line arguments";
285 
286  CmdLnArgParser parser(argc,argv);
287 
288  // get GENIE event sample
289  if( parser.OptionExists('f') ) {
290  LOG("rwghtzexpaxff", pINFO) << "Reading event sample filename";
291  gOptInpFilename = parser.ArgAsString('f');
292  } else {
293  LOG("rwghtzexpaxff", pFATAL)
294  << "Unspecified input filename - Exiting";
295  PrintSyntax();
296  exit(1);
297  }
298 
299  // output weight file
300  if(parser.OptionExists('o')) {
301  LOG("rwghtzexpaxff", pINFO) << "Reading requested output filename";
302  gOptOutFilename = parser.ArgAsString('o');
303  } else {
304  LOG("rwghtzexpaxff", pINFO) << "Setting default output filename";
305  gOptOutFilename = "test_rw_zexp_axff.root";
306  }
307 
308  if ( parser.OptionExists('n') ) {
309  //
310  LOG("grwghtzexpaxff", pINFO) << "Reading number of events to analyze";
311  string nev = parser.ArgAsString('n');
312  if (nev.find(",") != string::npos) {
313  vector<long> vecn = parser.ArgAsLongTokens('n',",");
314  if(vecn.size()!=2) {
315  LOG("grwghtzexpaxff", pFATAL) << "Invalid syntax";
316  gAbortingInErr = true;
317  PrintSyntax();
318  exit(1);
319  }
320  // User specified a comma-separated set of values n1,n2.
321  // Use [n1,n2] as the event range to process.
322  gOptNEvt1 = vecn[0];
323  gOptNEvt2 = vecn[1];
324  } else {
325  // User specified a single number n.
326  // Use [0,n] as the event range to process.
327  gOptNEvt1 = -1;
328  gOptNEvt2 = parser.ArgAsLong('n');
329  }
330  } else {
331  LOG("grwghtzexpaxff", pINFO)
332  << "Unspecified number of events to analyze - Use all";
333  gOptNEvt1 = -1;
334  gOptNEvt2 = -1;
335  }
336  LOG("grwghtzexpaxff", pDEBUG)
337  << "Input event range: " << gOptNEvt1 << ", " << gOptNEvt2;
338 
339  // values to reweight to:
340  if( parser.OptionExists('v') ) {
341  LOG("rwghtzexpaxff", pINFO) << "Reading requested parameter values";
342  string coef = parser.ArgAsString('v');
343 
344  vector<string> zvals = utils::str::Split(coef, ",");
345  // MAX_COEF should be set to the number of z-expansion tweaks which exist
346  assert(zvals.size() == (unsigned int) MAX_COEF);
347  for (int ik = 0;ik<MAX_COEF;ik++)
348  {
349  gOptParameterValue[ik] = atof(zvals[ik].c_str());
350  }
351  for (int ik = 0;ik<MAX_COEF;ik++)
352  {
353  LOG("rwghtzexpaxff",pINFO)<<"Parameter value "<<ik+1<<": "<<
354  gOptParameterValue[ik];
355  }
356  }
357 
358  // norm to reweight to:
359  if( parser.OptionExists('m') ) {
360  LOG("rwghtzexpaxff", pINFO) << "Reading requested normalization";
361  string coef = parser.ArgAsString('m');
362  gOptDoNorm = true;
363  gOptNormValue = atof(coef.c_str());
364  LOG("rwghtzexpaxff",pINFO)<<"Requested normalization: "<< gOptNormValue;
365  }
366 
367 }
368 //_________________________________________________________________________________
369 void GetEventRange(Long64_t nev_in_file, Long64_t & nfirst, Long64_t & nlast)
370 {
371  nfirst = 0;
372  nlast = 0;
373 
374  if(gOptNEvt1>=0 && gOptNEvt2>=0) {
375  // Input was `-n N1,N2'.
376  // Process events [N1,N2].
377  // Note: Incuding N1 and N2.
378  nfirst = gOptNEvt1;
379  nlast = TMath::Min(nev_in_file-1, gOptNEvt2);
380  }
381  else
382  if(gOptNEvt1<0 && gOptNEvt2>=0) {
383  // Input was `-n N'.
384  // Process first N events [0,N).
385  // Note: Event N is not included.
386  nfirst = 0;
387  nlast = TMath::Min(nev_in_file-1, gOptNEvt2-1);
388  }
389  else
390  if(gOptNEvt1<0 && gOptNEvt2<0) {
391  // No input. Process all events.
392  nfirst = 0;
393  nlast = nev_in_file-1;
394  }
395 
396  assert(nfirst <= nlast && nfirst >= 0 && nlast <= nev_in_file-1);
397 }
398 //_________________________________________________________________________________
399 GSyst_t GetZExpSystematic(int ip)
400 {
401  switch(ip){
402  case 1: return kXSecTwkDial_ZExpA1CCQE; break;
403  case 2: return kXSecTwkDial_ZExpA2CCQE; break;
404  case 3: return kXSecTwkDial_ZExpA3CCQE; break;
405  case 4: return kXSecTwkDial_ZExpA4CCQE; break;
406  default:
407  LOG("rwghtzexpaxff", pFATAL)
408  << "Cannot find systematic corresponding to parameter " << ip;
409  exit(0);
410  break;
411  }
412 }
413 //_________________________________________________________________________________
414 void PrintSyntax(void)
415 {
416  LOG("grwghtzexpaxff", pFATAL)
417  << "\n\n"
418  << "grwghtzexpaxff \n"
419  << " -f input_event_file \n"
420  << " -v val1,val2,val3,val4 \n"
421  << " [-n nev] \n"
422  << " [-o output_weights_file]";
423 }
424 //_________________________________________________________________________________
int iev
Definition: runWimpSim.h:118
::xsd::cxx::tree::id< char, ncname > id
Definition: Database.h:165
Cross Section Calculation Interface.
RgDbl GetDoubleDef(RgKey key, RgDbl def_opt, bool set_def=true)
Definition: Registry.cxx:550
long ArgAsLong(char opt)
void PrintSyntax()
string ArgAsString(char opt)
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define MAX_COEF
A simple program to reweight the GENIE z-expansion axial form factor from one set of z-expansion para...
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:41
#define pFATAL
Definition: Messenger.h:57
vector< long > ArgAsLongTokens(char opt, string delimeter)
int main(int argc, char **argv)
Algorithm abstract base class.
Definition: Algorithm.h:54
TString ip
Definition: loadincs.C:5
RgDbl GetDouble(RgKey key) const
Definition: Registry.cxx:489
void GetEventRange(Long64_t nev_in_file, Long64_t &nfirst, Long64_t &nlast)
virtual const Registry & GetConfig(void) const
Definition: Algorithm.cxx:254
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:97
void GetCommandLineArgs(int argc, char **argv)
GSyst_t GetZExpSystematic(int ip)
MINOS-style Ntuple Class to hold an output MC Tree Header.
Var weights
Long64_t gOptNEvt2
#define pINFO
Definition: Messenger.h:63
void AdoptSubstructure(void)
Definition: Algorithm.cxx:408
Algorithm * AdoptAlgorithm(const AlgId &algid) const
Definition: AlgFactory.cxx:127
bool gOptDoNorm
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
Definition: EventRecord.h:38
Algorithm ID (algorithm name + configuration set name)
Definition: AlgId.h:35
string gOptOutFilename
vector< string > Split(string input, string delim)
Definition: StringUtils.cxx:42
static AlgFactory * Instance()
Definition: AlgFactory.cxx:75
double gOptNormValue
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
string gOptInpFilename
exit(0)
TFile * file
Definition: cellShifts.C:17
assert(nhit_max >=nhit_nbins)
Command line argument parser.
#define pNOTICE
Definition: Messenger.h:62
The GENIE Algorithm Factory.
Definition: AlgFactory.h:40
Long64_t gOptNEvt1
void Clear(Option_t *opt="")
bool gAbortingInErr
Definition: Messenger.cxx:56
int n_events
Definition: test_Tier.py:22
bool OptionExists(char opt)
was option set?
EventRecord * event
event
static AlgConfigPool * Instance()
#define pDEBUG
Definition: Messenger.h:64
enum BeamMode string