makeNUSYS.cxx
Go to the documentation of this file.
1 #include "NUSYS.C"
2 #include "TRandom3.h"
3 #include "TFile.h"
4 #include "TTree.h"
5 #include "TVector3.h"
6 #include "TLorentzVector.h"
8 #include "EVGCore/EventRecord.h"
9 #include "nusystematics/artless/response_helper.hh"
10 #include <stdio.h>
11 
12 TRandom3 * rando;
13 
14 struct params {
15  double OA_xcoord;
16  bool fhc, grid;
18 };
19 
20 // main loop function
21 void loop( CAF &caf, params &par, TTree * tree, std::string ghepdir, std::string fhicl_filename )
22 {
23 
24  // Get GHEP file for genie::EventRecord from other file
25  int ievt, ifileNo;
26  tree->SetBranchAddress( "ievt", &ievt );
27  tree->SetBranchAddress( "ifileNo", &ifileNo );
28  int current_file = -1;
29  TFile * ghep_file = NULL;
30  TTree * gtree = NULL;
31  std::string mode = ( par.fhc ? "neutrino" : "antineutrino" );
32 
33  // DUNE reweight getter
34  nusyst::response_helper rh( fhicl_filename );
35  // Get list of variations, and make tree branch for each one
36  std::vector<unsigned int> parIds = rh.GetParameters();
37  for( unsigned int i = 0; i < parIds.size(); ++i ) {
38  systtools::SystParamHeader head = rh.GetHeader(parIds[i]);
39  printf( "Adding reweight branch %u for %s with %lu shifts\n", parIds[i], head.prettyName.c_str(), head.paramVariations.size() );
40  bool is_wgt = head.isWeightSystematicVariation;
41  std::string wgt_var = ( is_wgt ? "wgt" : "var" );
42  caf.addRWbranch( parIds[i], head.prettyName, wgt_var, head.paramVariations );
43  caf.iswgt[parIds[i]] = is_wgt;
44  }
45 
46 
47  // Main event loop
48  int N = tree->GetEntries();
49  if( par.n > 0 && par.n < N ) N = par.n + par.first;
50  for( int ii = par.first; ii < N; ++ii ) {
51 
52  tree->GetEntry(ii);
53  if( ii % 100 == 0 ) printf( "Event %d of %d...\n", ii, N );
54 
55 
56  // set defaults
57  for( int j = 0; j < 100; ++j ) {
58  caf.nwgt[j] = 7;
59  caf.cvwgt[j] = ( caf.iswgt[j] ? 1. : 0. );
60  for( unsigned int k = 0; k < 100; ++k ) {
61  caf.wgt[j][k] = ( caf.iswgt[j] ? 1. : 0. );
62  }
63  }
64 
65  // make sure ghep file matches the current one, otherwise update to the current ghep file
66  if( ifileNo != current_file ) {
67  // close the previous file
68  if( ghep_file ) ghep_file->Close();
69 
70  if( par.grid ) ghep_file = new TFile( Form("genie.%d.root", ifileNo) );
71  else ghep_file = new TFile( Form("%s/GAr.%s.%d.ghep.root", ghepdir.c_str(), mode.c_str(), ifileNo) );
72 
73  gtree = (TTree*) ghep_file->Get( "gtree" );
74 
75 
76  gtree->SetBranchAddress( "gmcrec", &caf.mcrec );
77  current_file = ifileNo;
78  }
79 
80 
81 
82  // get GENIE event record
83 
84  gtree->GetEntry( ievt );
85  genie::EventRecord * event = caf.mcrec->event;
86  genie::Interaction * in = event->Summary();
87 
88  // Get truth stuff out of GENIE ghep record
89 
90  TLorentzVector lepP4;
91  TLorentzVector nuP4nuc = *(in->InitState().GetProbeP4(genie::kRfHitNucRest));
92  TLorentzVector nuP4 = *(in->InitState().GetProbeP4(genie::kRfLab));
93 
94 
95 
96 
97  // Add DUNErw weights to tree
98 
99  systtools::event_unit_response_w_cv_t resp = rh.GetEventVariationAndCVResponse(*event);
100  for( systtools::event_unit_response_w_cv_t::iterator it = resp.begin(); it != resp.end(); ++it ) {
101  caf.nwgt[(*it).pid] = (*it).responses.size();
102  caf.cvwgt[(*it).pid] = (*it).CV_response;
103  for( unsigned int i = 0; i < (*it).responses.size(); ++i ) {
104  caf.wgt[(*it).pid][i] = (*it).responses[i];
105  }
106  }
107 
108 
109  caf.fill();
110  }
111 
112  // set POT
113  caf.meta_run = par.run;
114  caf.meta_subrun = par.subrun;
115 
116 }
117 
118 int main( int argc, char const *argv[] )
119 {
120 
121  if( (argc == 2) && ((std::string("--help") == argv[1]) || (std::string("-h") == argv[1])) ) {
122  std::cout << "Help yourself by looking at the source code to see what the options are." << std::endl;
123  return 0;
124  }
125 
126  // get command line options
127  std::string ghepdir;
129  std::string edepfile;
130  std::string fhicl_filename;
131 
132  // Make parameter object and set defaults
133  params par;
134  par.OA_xcoord = 0.; // on-axis by default
135  par.fhc = true;
136  par.grid = false;
137  par.seed = 7; // a very random number
138  par.run = 1; // CAFAna doesn't like run number 0
139  par.subrun = 0;
140  par.n = -1;
141  par.nfiles = 1;
142  par.first = 0;
143 
144  int i = 0;
145  while( i < argc ) {
146  if( argv[i] == std::string("--edepfile") ) {
147  edepfile = argv[i+1];
148  i += 2;
149  } else if( argv[i] == std::string("--ghepdir") ) {
150  ghepdir = argv[i+1];
151  i += 2;
152  } else if( argv[i] == std::string("--outfile") ) {
153  outfile = argv[i+1];
154  i += 2;
155  } else if( argv[i] == std::string("--fhicl") ) {
156  fhicl_filename = argv[i+1];
157  i += 2;
158  } else if( argv[i] == std::string("--seed") ) {
159  par.seed = atoi(argv[i+1]);
160  par.run = par.seed;
161  i += 2;
162  } else if( argv[i] == std::string("--nevents") ) {
163  par.n = atoi(argv[i+1]);
164  i += 2;
165  } else if( argv[i] == std::string("--nfiles") ) {
166  par.nfiles = atoi(argv[i+1]);
167  i += 2;
168  } else if( argv[i] == std::string("--first") ) {
169  par.first = atoi(argv[i+1]);
170  i += 2;
171  } else if( argv[i] == std::string("--oa") ) {
172  par.OA_xcoord = atof(argv[i+1]);
173  i += 2;
174  } else if( argv[i] == std::string("--grid") ) {
175  par.grid = true;
176  i += 1;
177  } else if( argv[i] == std::string("--rhc") ) {
178  par.fhc = false;
179  i += 1;
180  } else i += 1; // look for next thing
181  }
182 
183  printf( "Making CAF from edep-sim tree dump: %s\n", edepfile.c_str() );
184  printf( "Searching for GENIE ghep files here: %s\n", ghepdir.c_str() );
185  if( par.fhc ) printf( "Running neutrino mode (FHC)\n" );
186  else printf( "Running antineutrino mode (RHC)\n" );
187  if( par.grid ) printf( "Running grid mode\n" );
188  else printf( "Running test mode\n" );
189  printf( "Output CAF file: %s\n", outfile.c_str() );
190 
191  rando = new TRandom3( par.seed );
192  CAF caf( outfile );
193 
194 
195  TFile * tf = new TFile( edepfile.c_str() );
196  TTree * tree = (TTree*) tf->Get( "tree" );
197 
198  loop( caf, par, tree, ghepdir, fhicl_filename );
199 
200  caf.version = 4;
201  printf( "Run %d POT %g\n", caf.meta_run, caf.pot );
202  caf.fillPOT();
203  caf.write();
204 
205  caf.cafFile->Close();
206 
207  printf( "-30-\n" );
208 
209 
210 }
intermediate_table::iterator iterator
void addRWbranch(int parId, std::string name, std::string wgt_var, std::vector< double > &vars)
int first
Definition: makeNUSYS.cxx:17
int meta_subrun
Definition: NUSYS.h:42
int seed
Definition: makeNUSYS.cxx:17
std::string string
Definition: nybbler.cc:12
bool fhc
Definition: makeNUSYS.cxx:16
bool iswgt[100]
Definition: NUSYS.h:35
void fill()
genie::NtpMCEventRecord * mcrec
Definition: NUSYS.h:38
double cvwgt[100]
Definition: NUSYS.h:33
double OA_xcoord
Definition: makeNUSYS.cxx:15
void fillPOT()
Definition: tf_graph.h:23
int subrun
Definition: makeNUSYS.cxx:17
int nfiles
Definition: makeNUSYS.cxx:17
void write()
int meta_run
Definition: NUSYS.h:42
Summary information for an interaction.
Definition: Interaction.h:56
TFile * cafFile
The output TFile pointer */.
Definition: CAF.h:53
int run
Definition: makeNUSYS.cxx:17
bool grid
Definition: makeNUSYS.cxx:16
int version
Definition: NUSYS.h:43
Definition: CAF.h:12
int main(int argc, char const *argv[])
Definition: makeNUSYS.cxx:118
int n
Definition: makeNUSYS.cxx:17
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
Definition: EventRecord.h:37
int nwgt[100]
Definition: NUSYS.h:32
const InitialState & InitState(void) const
Definition: Interaction.h:69
Common Analysis Files.
Definition: SRGAr.h:13
double wgt[100][100]
Definition: NUSYS.h:34
void loop(CAF &caf, params &par, TTree *tree, std::string ghepdir, std::string fhicl_filename)
Definition: makeNUSYS.cxx:21
EventRecord * event
event
QTextStream & endl(QTextStream &s)
TLorentzVector * GetProbeP4(RefFrame_t rf=kRfHitNucRest) const
TRandom3 * rando
Definition: makeNUSYS.cxx:12
Event finding and building.
double pot
Definition: NUSYS.h:41