gRWIOExample1.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3 
4 \program gRWIOExample1
5 
6 \brief Simple example to demostrate use of GENIE ReWeight I/O infrastructure.
7  It is a 2-step procedure.
8  First, it processes ALL events in the input file, but performs re-weighting
9  only where it is applicable.
10  An "empty" RW IO record will still be generated and written out for those
11  events where re-weighting does not apply.
12  Second, it opens the GENIE file again, this time in the "READ" mode, and
13  extracts and prints the RW information.
14  WARNING-1: this example is NOT entirely "fool proof" !!!
15  WARNING-2: this example can be quite verbose, due to all the printouts at step-2 !!!
16 
17  Syntax :
18  gRWIOExample1 -f input_ghep_file
19 
20  -f : Input GENIE event file
21 
22 \author(s) Julia Yarba (FNAL)
23 
24 \created May 11, 2016
25 
26 \cpright Copyright (c) 2003-2016, GENIE Neutrino MC Generator Collaboration
27  All rights reserved.
28  For the licensing terms see $GENIE/USER_LICENSE.
29 */
30 //____________________________________________________________________________
31 
32 #include <string>
33 #include <sstream>
34 #include <cassert>
35 
36 #include <TSystem.h>
37 #include <TFile.h>
38 #include <TTree.h>
39 
40 #include "EVGCore/EventRecord.h"
42 
44 
45 #include "ReWeight/GReWeightI.h"
46 #include "ReWeight/GSystSet.h"
47 #include "ReWeight/GSyst.h"
48 #include "ReWeight/GReWeight.h"
49 #include "ReWeight/GSystUncertainty.h"
50 
51 // Modules to calc weights/uncertainties
52 // by varying specific variables
53 //
54 // This one right below is for MaCCQE
55 //
56 #include "ReWeight/GReWeightNuXSecCCQE.h"
57 
58 // I/O for re-weighting info
59 //
60 #include "ReWeight/GReWeightIOBranchDesc.h"
61 #include "ReWeight/GReWeightIORecord.h"
62 
63 #include "Utils/CmdLnArgParser.h"
64 
65 int main( int argc, char ** argv )
66 {
67 
68  // GENIE event file on input
69  //
70  std::string isample = "";
71 
72  genie::CmdLnArgParser* lp = new genie::CmdLnArgParser( argc, argv );
73  if ( lp->OptionExists('f') )
74  {
75  isample = lp->ArgAsString('f');
76  }
77 
78  if ( isample.empty() )
79  {
80  // Nothing to re-weight, bail out... may also a warning will be useful...
81  std::cout << " Missing GENIE event file on input " << std::endl;
82  return 1;
83  }
84 
85  // Define what parameter we want to vary/tweak
86  //
87  // In principle, it should be run-time configurable
88  // But for simplicity we will define it here in this particular case
89  //
90  genie::rew::GSyst_t param_to_tweak = genie::rew::kXSecTwkDial_MaCCQE ;
91 
92  // Create a GReWeight object and add to it a set of weight calculators
93  //
94  genie::rew::GReWeight rw;
95 
96  // Add weight calculator for MaCCQE
97  // NOTE: will add other weight calculators later
98  //
99  rw.AdoptWghtCalc( "xsec_ccqe", new genie::rew::GReWeightNuXSecCCQE() );
100 
101  // Get GSystSet and include the (single) input systematic parameter
102  // NOTE: in this case we use kXSecTwkDial_MaCCQE (for "MaCCQE")
103  //
104  genie::rew::GSystSet& syst = rw.Systematics();
105  // syst.Init( genie::rew::kXSecTwkDial_MaCCQE );
106  syst.Init( param_to_tweak );
107 
108  // By default GReWeightNuXSecCCQE is in `NormAndMaShape' mode
109  // where Ma affects the shape of dsigma/dQ2 and a different param affects the normalization
110  // If the input is MaCCQE, switch the weight calculator to `Ma' mode
111  //
112  genie::rew::GReWeightNuXSecCCQE* rwccqe = dynamic_cast<genie::rew::GReWeightNuXSecCCQE*>( rw.WghtCalc("xsec_ccqe") );
113  rwccqe->SetMode( genie::rew::GReWeightNuXSecCCQE::kModeMa );
114 
115  // Now open input GENIE sample (fetched at the beginning of the job)
116  //
117  TFile* file = new TFile( isample.c_str(), "UPDATE" );
118 
119  // if invalid input file, bail out
120  //
121  if ( !file )
122  {
123  std::cout << " Can NOT open input GENIE file " << isample << std::endl;
124  return 1;
125  }
126 
127  //
128  // Fetch or create a tree for RW records
129  //
130  TTree* rwtree = 0;
131  rwtree = dynamic_cast<TTree*>( file->Get("reweighting") );
132  if ( !rwtree )
133  {
134  rwtree = new TTree( "reweighting", "GENIE weights tree" );
135  TTree::SetBranchStyle(1);
136  rwtree->SetAutoSave( 200000000 ); // autosave when 0.2 Gbyte written
137  // - it's the same for "gtree" but I need to double check
138  // how to *get* autosave from "gtree"
139  }
140 
141  // Now create a branch to correspond to a specific parameter/variable to vary
142  //
143  // FIXME !!!
144  // In principle, one should also check if a branch, and the corresponding metadata already exist, etc.
145  //
146  genie::rew::GReWeightIORecord* rwrec = 0;
147  std::string param_name = genie::rew::GSyst::AsString( param_to_tweak );
148  TBranch* rwbr = rwtree->Branch( param_name.c_str(),
149  "genie::rew::GReWeightIORecord",
150  &rwrec, 32000, 1 ); // FIXME !!! also check more "sophisticated" options
151  assert(rwbr);
152  rwbr->SetAutoDelete(kFALSE);
153 
154  // Add meta-data (UserInfo) to the RW tree
155  //
156  // MaCCQE=0.99 has been extracted in the course of run under gdb from GReWeightNuXSecCCQE class
157  // (member data fMaDef).
158  // In principle, it depends on the physics model - in the case of GReWeightNuXSecCCQE the model
159  // is based on the "LwlynSmithQELCCPXSec" algorithm (see GReWeightNuXSecCCQE::Init() method)
160  // A more uniform machinery to access such information would be useful, but details of it need
161  // to be discussed additionally.
162  //
163  // Sigma's (+/-) can be extracted from GSystUncertainty
164  //
165  genie::rew::GSystUncertainty* syser = genie::rew::GSystUncertainty::Instance();
166  double sigpls = syser->OneSigmaErr( param_to_tweak, 1 );
167  double sigmin = syser->OneSigmaErr( param_to_tweak, -1 );
168  //
169  rwtree->GetUserInfo()->Add( new genie::rew::GReWeightIOBranchDesc( param_name, 0.99, sigpls, sigmin ) );
170 
171  // Fetch the Evt tree
172  //
173  TTree* evtree = dynamic_cast<TTree*>( file->Get("gtree") );
174 
175  // Connect Evt record (branch)
176  //
177  genie::NtpMCEventRecord* mcrec = 0;
178  evtree->SetBranchAddress( "gmcrec", &mcrec );
179 
180  // "Tie" together these trees, Evt & RW !
181  //
182  evtree->AddFriend( rwtree );
183 
184  // now loop over events and see what needs to be re-weighted
185  //
186  int nevt_total = evtree->GetEntries();
187 
188  double twk = 0.;
189  double wt = 1.;
190 
191  for ( int iev=0; iev<nevt_total; ++iev )
192  {
193 
194  evtree->GetEntry(iev);
195  genie::EventRecord& evt = *(mcrec->event);
196 
197  //
198  // Select events to be re-weighted
199  //
200  // Specifically, check if it's QEL && WeakCC process
201  // because that's what we want to reweight (MaCCQE).
202  // Also skip charm events (although if those are quite rare in this case)
203  //
204  genie::Interaction* interaction = evt.Summary();
205  const genie::ProcessInfo& prinfo = interaction->ProcInfo();
206  const genie::XclsTag& xclsv = interaction->ExclTag();
207  bool accept = ( prinfo.IsQuasiElastic() && prinfo.IsWeakCC() && !xclsv.IsCharmEvent() );
208  if ( !accept )
209  {
210  rwrec = new genie::rew::GReWeightIORecord();
211  rwrec->SetOriginalEvtNumber(iev);
212  rwtree->Fill();
213  delete rwrec;
214  rwrec=0;
215  continue ;
216  }
217 
218  rwrec = new genie::rew::GReWeightIORecord();
219 
220  rwrec->SetOriginalEvtNumber( iev );
221 
222  twk = -0.5; // in the units of MaCCQE SIGMA !!! (that's how the weigh calculator "understands" it)
223  syst.Set( param_to_tweak, twk );
224  rw.Reconfigure();
225  wt = rw.CalcWeight(evt);
226  rwrec->Insert( twk, wt );
227 
228  twk = 0.5;
229  syst.Set( param_to_tweak, twk );
230  rw.Reconfigure();
231  wt = rw.CalcWeight(evt);
232  rwrec->Insert( twk, wt );
233 
234  rwtree->Fill();
235 
236  // Clear mc evt record before the next one
237  //
238  mcrec->Clear();
239  if ( rwrec ) delete rwrec;
240  rwrec = 0;
241 
242  }
243 
244  file->cd();
245  rwtree->Write("",TObject::kOverwrite);
246 
247  delete rwtree;
248  rwtree=0;
249 
250  file->Close();
251 
252  delete lp; // destroy input args line parser
253 
254 
255  // RE_TEST NOW !!!
256  //
257  // Open up Genie event file (now also with the RW tree in it),
258  // this time in the READ mode
259  //
260  TFile* tfile = new TFile( isample.c_str(), "READ" );
261 
262  //
263  // Fetch the RW tree
264  //
265  TTree* rwtree_test = dynamic_cast<TTree*>( tfile->Get("reweighting") );
266 
267  TList* hdr = rwtree_test->GetUserInfo();
268  assert(hdr);
269 
270  int nentries = rwtree_test->GetEntries();
271  int nrw = 0;
272  std::cout << " num of entries of re-test: " << nentries << std::endl;
273  genie::rew::GReWeightIORecord* rwrec_test = 0;
274  rwtree_test->SetBranchAddress( param_name.c_str(), &rwrec_test );
275  for ( int i=0; i<nentries; ++i )
276  {
277  rwtree_test->GetEntry(i);
278  // genie::EventRecord& evt = *(mcrec->event);
279  int nres = rwrec_test->GetNumOfRWResults();
280  if ( nres <= 0 ) continue;
281  for ( int ir=0; ir<nres; ++ir )
282  {
283  double twk_test = rwrec_test->GetTweak( ir );
284  double wt_test = rwrec_test->GetWeight( ir );
285  std::cout << " twk = " << twk_test << " wt = " << wt_test << std::endl;
286  }
287  nrw++;
288  }
289 
290  std::cout << " num of re-weighted results of re-test: " << nrw << std::endl;
291 
292  // now print meta-data
293  //
294  int nhdr = hdr->GetEntries();
295  for ( int i=0; i<nhdr; ++i )
296  {
297  genie::rew::GReWeightIOBranchDesc* brdesc = dynamic_cast<genie::rew::GReWeightIOBranchDesc*>( hdr->At(i) );
298  std::cout << " branch name: " << brdesc->GetParameterName() << std::endl;
299  std::cout << " parameter: " << brdesc->GetParameterMean() << " "
300  << brdesc->GetParameterSigmaPlus() << " " << brdesc->GetParameterSigmaMinus() << std::endl;
301  }
302 
303  return 0;
304 
305 }
306 
int iev
Definition: runWimpSim.h:118
bool IsWeakCC(void) const
string ArgAsString(char opt)
virtual Interaction * Summary(void) const
Definition: GHepRecord.cxx:135
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.
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:67
Contains minimal information for tagging exclusive processes.
Definition: XclsTag.h:37
bool IsCharmEvent(void) const
Definition: XclsTag.h:48
Summary information for an interaction.
Definition: Interaction.h:56
Long64_t nentries
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:44
int evt
int main(int argc, char **argv)
Simple example to demostrate use of GENIE ReWeight I/O infrastructure. It is a 2-step procedure...
OStream cout
Definition: OStream.cxx:6
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
Definition: EventRecord.h:38
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
TFile * file
Definition: cellShifts.C:17
assert(nhit_max >=nhit_nbins)
const char * AsString(Resonance_t res)
resonance id -> string
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
Command line argument parser.
void Clear(Option_t *opt="")
bool OptionExists(char opt)
was option set?
EventRecord * event
event
enum BeamMode string