GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
gtestRewght.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3 
4 \program gtestRewght
5 
6 \brief A simple program to illustrate how to use the GENIE event reweighting.
7 
8 \syntax gtestRewght -f filename [-n nev]
9 
10  where
11  [] is an optional argument
12  -f specifies a GENIE event file (GHEP format)
13  -n specifies the number of events to process (default: all)
14 
15 \author Costas Andreopoulos <c.andreopoulos \at cern.ch>
16  University of Liverpool
17 
18 \created May 19, 2010
19 
20 \cpright Copyright (c) 2003-2024, The GENIE Collaboration
21  For the full text of the license visit http://copyright.genie-mc.org
22 
23 */
24 //____________________________________________________________________________
25 
26 
27 #include <string>
28 
29 #include <TSystem.h>
30 #include <TFile.h>
31 #include <TTree.h>
32 
38 #include "Tools/ReWeight/GReWeightI.h"
39 #include "Tools/ReWeight/GSystSet.h"
40 #include "Tools/ReWeight/GReWeight.h"
41 #include "Tools/ReWeight/GReWeightNuXSecCCQE.h"
42 #include "Tools/ReWeight/GReWeightNuXSecCCQEvec.h"
43 #include "Tools/ReWeight/GReWeightNuXSecCCRES.h"
44 #include "Tools/ReWeight/GReWeightNuXSecNCRES.h"
45 #include "Tools/ReWeight/GReWeightNuXSecDIS.h"
46 #include "Tools/ReWeight/GReWeightNuXSecCOH.h"
47 #include "Tools/ReWeight/GReWeightNonResonanceBkg.h"
48 #include "Tools/ReWeight/GReWeightFGM.h"
49 #include "Tools/ReWeight/GReWeightDISNuclMod.h"
50 #include "Tools/ReWeight/GReWeightResonanceDecay.h"
51 #include "Tools/ReWeight/GReWeightFZone.h"
52 #include "Tools/ReWeight/GReWeightINuke.h"
53 #include "Tools/ReWeight/GReWeightAGKY.h"
55 
56 using std::string;
57 
58 using namespace genie;
59 using namespace genie::rew;
60 
61 void GetCommandLineArgs (int argc, char ** argv);
62 
65 
66 //___________________________________________________________________
67 int main(int argc, char ** argv)
68 {
69  GetCommandLineArgs (argc, argv);
70 
71  // open the ROOT file and get the TTree & its header
72  TTree * tree = 0;
73  NtpMCTreeHeader * thdr = 0;
74 
75  TFile file(gOptInpFilename.c_str(),"READ");
76 
77  tree = dynamic_cast <TTree *> ( file.Get("gtree") );
78  thdr = dynamic_cast <NtpMCTreeHeader *> ( file.Get("header") );
79 
80  if(!tree) return 1;
81 
82  LOG("test", pNOTICE) << "Input tree header: " << *thdr;
83 
84  int nev = (gOptNEvt > 0) ?
85  TMath::Min(gOptNEvt, (int)tree->GetEntries()) :
86  (int) tree->GetEntries();
87 
88  LOG("test", pNOTICE) << "Will process " << nev << " events";
89 
90  //
91  // Create a GReWeight object and add to it a set of
92  // weight calculators
93  //
94 
95  GReWeight rw;
96 
97  rw.AdoptWghtCalc( "xsec_ccqe", new GReWeightNuXSecCCQE );
98  rw.AdoptWghtCalc( "xsec_ccqe_vec", new GReWeightNuXSecCCQEvec );
99  rw.AdoptWghtCalc( "xsec_ccres", new GReWeightNuXSecCCRES );
100  rw.AdoptWghtCalc( "xsec_ncres", new GReWeightNuXSecNCRES );
101  rw.AdoptWghtCalc( "xsec_nonresbkg", new GReWeightNonResonanceBkg );
102  rw.AdoptWghtCalc( "xsec_dis", new GReWeightNuXSecDIS );
103  rw.AdoptWghtCalc( "xsec_coh", new GReWeightNuXSecCOH );
104  rw.AdoptWghtCalc( "nuclear_qe", new GReWeightFGM );
105  rw.AdoptWghtCalc( "nuclear_dis", new GReWeightDISNuclMod );
106  rw.AdoptWghtCalc( "hadro_res_decay", new GReWeightResonanceDecay );
107  rw.AdoptWghtCalc( "hadro_fzone", new GReWeightFZone );
108  rw.AdoptWghtCalc( "hadro_intranuke", new GReWeightINuke );
109  rw.AdoptWghtCalc( "hadro_agky", new GReWeightAGKY );
110 
111  //
112  // Create a list of systematic params (more to be found at GSyst.h)
113  // set non-default values and re-configure.
114  // Weight calculators included above must be able to handle the tweaked params.
115  // Each tweaking dial t modifies a physics parameter p as:
116  // p_{tweaked} = p_{default} ( 1 + t * dp/p )
117  // So setting a tweaking dial to +/-1 modifies a physics quantity
118  // by +/- 1sigma.
119  // Default fractional errors are defined in GSystUncertainty
120  // and can be overriden.
121  //
122 
123  GSystSet & syst = rw.Systematics();
124 
125  syst.Set(kXSecTwkDial_NormCCQE, +1.0);
126  syst.Set(kXSecTwkDial_MaCCQEshape, +1.0);
127  syst.Set(kXSecTwkDial_NormCCRES, -1.0);
128  syst.Set(kXSecTwkDial_VecFFCCQEshape, -1.0);
129  syst.Set(kXSecTwkDial_MaCCRESshape, -1.0);
130  syst.Set(kXSecTwkDial_MvCCRESshape, +0.5);
131  syst.Set(kXSecTwkDial_NormNCRES, +1.0);
132  syst.Set(kXSecTwkDial_MaNCRESshape, -0.7);
133  syst.Set(kXSecTwkDial_MvNCRESshape, +0.3);
134  syst.Set(kXSecTwkDial_RvpCC1pi, +0.5);
135  syst.Set(kXSecTwkDial_RvnCC1pi, +0.5);
136  syst.Set(kXSecTwkDial_MaCOHpi, -0.5);
137  syst.Set(kINukeTwkDial_MFP_pi, +1.0);
138  syst.Set(kINukeTwkDial_MFP_N, -1.0);
139  syst.Set(kINukeTwkDial_FrPiProd_pi, -0.7);
140  syst.Set(kHadrAGKYTwkDial_xF1pi, -1.0);
141  syst.Set(kHadrAGKYTwkDial_pT1pi, +1.0);
142  syst.Set(kHadrNuclTwkDial_FormZone, +1.0);
143  syst.Set(kRDcyTwkDial_Theta_Delta2Npi, +1.0);
144 
145  rw.Reconfigure();
146 
147  //
148  // Concrete weight calculators can be retrieved and fine-tuned.
149  // For example:
150 
151  GReWeightNuXSecCCQE * rwccqe =
152  dynamic_cast<GReWeightNuXSecCCQE *> (rw.WghtCalc("xsec_ccqe"));
153  rwccqe -> RewNue (false);
154  rwccqe -> RewNuebar (false);
155  rwccqe -> RewNumubar(false);
156 
157  //
158  // Event loop
159  //
160 
161  NtpMCEventRecord * mcrec = 0;
162  tree->SetBranchAddress("gmcrec", &mcrec);
163 
164  for(int i = 0; i < nev; i++) {
165  tree->GetEntry(i);
166 
167  EventRecord & event = *(mcrec->event);
168  LOG("test", pNOTICE) << event;
169 
170  double wght = rw.CalcWeight(event);
171  LOG("test", pNOTICE) << "Overall weight = " << wght;
172 
173  mcrec->Clear();
174  }
175 
176  file.Close();
177 
178  LOG("test", pNOTICE) << "Done!";
179  return 0;
180 }
181 //___________________________________________________________________
182 void GetCommandLineArgs(int argc, char ** argv)
183 {
184  LOG("test", pINFO) << "*** Parsing command line arguments";
185 
186  CmdLnArgParser parser(argc,argv);
187 
188  // get GENIE event sample
189  if( parser.OptionExists('f') ) {
190  LOG("test", pINFO) << "Reading event sample filename";
191  gOptInpFilename = parser.ArgAsString('f');
192  } else {
193  LOG("test", pFATAL)
194  << "Unspecified input filename - Exiting";
195  exit(1);
196  }
197 
198  // number of events:
199  if( parser.OptionExists('n') ) {
200  LOG("test", pINFO) << "Reading number of events to analyze";
201  gOptNEvt = parser.ArgAsInt('n');
202  } else {
203  LOG("test", pINFO)
204  << "Unspecified number of events to analyze - Use all";
205  gOptNEvt = -1;
206  }
207 }
208 //_________________________________________________________________________________
int gOptNEvt
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
int main(int argc, char **argv)
Definition: gAtmoEvGen.cxx:327
#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
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
Command line argument parser.
#define pNOTICE
Definition: Messenger.h:61
void GetCommandLineArgs(int argc, char **argv)
Definition: gAtmoEvGen.cxx:563
void Clear(Option_t *opt="")
EventRecord * event
event