warp_gsimple.C
Go to the documentation of this file.
1 //========================================================================
2 //
3 // This is a script for creating a modified GSimpleNtpFlux file from
4 // an original file, but modifying entries. Assumes it will be passed
5 // the name of a single file.
6 //
7 // User must supply two functions "warp_entry()" and "warp_meta()"
8 //
9 //========================================================================
10 #include "GENIE/Framework/Numerical/RandomGen.h"
11 #include "GENIE/Tools/Flux/GSimpleNtpFlux.h"
12 #include "GENIE/Framework/Utils/UnitUtils.h"
13 
14 #include "TSystem.h"
15 #include "TStopwatch.h"
16 #include "TLorentzVector.h"
17 #include "TNtuple.h"
18 #include "TFile.h"
19 #include "TChain.h"
20 #include "TF1.h"
21 #include "TH1F.h"
22 #include "TRandom.h"
23 
24 #include <iostream>
25 #include <iomanip>
26 #include <string>
27 #include <sstream>
28 #include <set>
29 #include <algorithm>
30 
31 using namespace std;
32 using namespace genie;
33 using namespace genie::flux;
34 
35 TF1* fcn;
36 // function prototype
37 double pick_energy(double, double);
38 
39 bool gDoWarp = true;
40 
41 //========================================================================
42 // this is what the USER must supply
43 
45  GSimpleNtpAux* aux,
46  GSimpleNtpNuMI* numi)
47 {
48  if ( ! gDoWarp ) return;
49  // this is a dumb weighting scheme
50  //if ( entry->E > 1.5 ) {
51  //// don't set absolute weight in case we started with weighted entries
52  //// scale whatever is there from GSimpleNtpFlux driver
53  //entry->wgt *= 0.5;
54  //}
55  //// change some flavors
56  //if ( gRandom->Rndm() < 0.25 ) {
57  //entry->pdg = ( entry->pdg > 0 ) ? 16 : -16;
58  //}
59 
60  // record the 4 momentum before warp
61  TLorentzVector fmb4 = TLorentzVector(entry->px, entry->py, entry->pz,
62  entry->E);
63  // sample an energy value from the 1/E distribution
64  // min energy 0.01 GeV, where GENIE spline starts
65  // max energy 100 GeV, by eye
66  entry->E = pick_energy(0.01, 50);
67  // scale momentum so that the new momentum is in the same direction as the
68  // old one but obeys the energy-momentum relation
69  double k = sqrt(1+(entry->E*entry->E-fmb4.E()*fmb4.E())/fmb4.Vect().Mag2());
70  entry->px *= k;
71  entry->py *= k;
72  entry->pz *= k;
73 
74  // use all branches so there are no complaints from compiler
75  static double trash = 0;
76  if ( numi ) trash = numi->entryno;
77  if ( aux ) trash = aux->auxint.size();
78  if ( trash == -1 ) cout << "trash was -1" << endl;
79 }
80 
81 void warp_meta(GSimpleNtpMeta* meta, string fnamein)
82 {
83  meta->infiles.push_back("THIS IS MDC WARPED FLUX FROM:");
84  meta->infiles.push_back(fnamein);
85  if ( ! gDoWarp ) {
86  meta->infiles.push_back("NO ACTUAL WARPING APPLIED");
87  } else {
88  string msg ="USING WARP FUNCTION CODENAMED WHISKEY-TANGO-FOXTROT";
89  meta->infiles.push_back(msg);
90  }
91 }
92 
93 double pick_energy(double emin, double emax)
94 {
95  // pick a distribution of shape 1/E
96  // between the values [emin, emax]
97  // must have low energy cut to avoid infrared divergence
98  // double r = gRandom->Rndm(); // flat variate between (0,1]
99  // // double c = log(emin);
100  // // double A = 1. / (log(emax) - c);
101  // // return exp(r/A + c);
102  // double c = emin;
103  // double A = 1./(emax - c);
104  // return r/A + c;
105  return fcn->GetRandom(emin, emax);
106 
107 }
108 
109 double fit_fcn(double* x, double* par){
110  double energy = x[0];
111  auto profile = [](double y, double* p){
112  return p[0]+p[1]*TMath::Power(y, -p[2])+p[3]*TMath::Power(y, -p[4])+p[5]*y;
113  };
114  double thresholdE = 0.1;
115  if(energy > thresholdE)
116  return profile(energy, par);
117  else
118  return profile(100., par);
119 }
120 
121 
122 void update_flavors(GSimpleNtpMeta* meta, std::set<int>& gFlavors);
123 
124 //========================================================================
125 // main routine
126 
127 // assumes a single input file (meta-data handling)
128 
129 
130 void warp_gsimple(string fnameout="gsimple_output.root",
131  string fnamein="gsimple_input.root",
132  bool dowarp=true)
133 {
134 
135  gRandom->SetSeed(0);
136  gDoWarp = dowarp;
137 
138  fcn = new TF1("fcn", fit_fcn, 0., 50, 6);
139  // fcn->SetParameters(-8.10907e+01,-2.83834e+02,-3.21121e+00,1.24471e+04,9.42167e-01);
140  fcn->SetParameters(-2.43864e+02, 4.42211e+02, 3.00555e+00, 1.23288e+04, 8.95128e-01, 2.38448e+00);
141  fcn->SetNpx(100000);
142 
143  cout << "-----------" << endl;
144  string fnameinlong(gSystem->ExpandPathName(fnamein.c_str()));
145  cout << fnameinlong << endl;
146  TFile* fin = TFile::Open(fnameinlong.c_str(),"READONLY");
147  TTree* etreein = (TTree*)fin->Get("flux");
148  TTree* mtreein = (TTree*)fin->Get("meta");
153 
154  long int nentries = etreein->GetEntries();
155 
156  int sba_status[4] = { -999, -999, -999, -999 };
157  sba_status[0] = etreein->SetBranchAddress("entry",&entry_in);
158  sba_status[1] = etreein->SetBranchAddress("numi",&numi_in);
159  sba_status[2] = etreein->SetBranchAddress("aux",&aux_in);
160  sba_status[3] = mtreein->SetBranchAddress("meta",&meta_in);
161  cout << "SetBranchAddress results "
162  << sba_status[0] << ","
163  << sba_status[1] << ","
164  << sba_status[2] << ","
165  << sba_status[3]
166  << endl;
167  bool donumi = ( sba_status[1] == 0 );
168  bool doaux = ( sba_status[2] == 0 );
169 
170  int nindices = mtreein->BuildIndex("metakey"); // tie metadata to entry
171  cout << "saw " << nindices << " metakey indices" << endl;
172 
173  cout << "Creating: " << fnameout << endl;
174  cout << "Input file: " << fnamein << endl;
175  cout << "Branches: entry "
176  << ( (doaux) ? "aux ":"" )
177  << ( (donumi) ? "numi ":"" )
178  << endl;
179  if ( ! gDoWarp ) cout << "++++++++ NO ACTUAL WARP APPLIED" << endl;
180 
181  TFile* fout = TFile::Open(fnameout.c_str(),"RECREATE");
182  TTree* fluxntp = new TTree("flux","a simple flux n-tuple");
183  TTree* metantp = new TTree("meta","metadata for flux n-tuple");
188 
189  fluxntp->Branch("entry",&fentry);
190  if ( doaux ) fluxntp->Branch("aux",&faux);
191  if ( donumi ) fluxntp->Branch("numi",&fnumi);
192  metantp->Branch("meta",&fmeta);
193 
194  cout << "=========================== Start " << endl;
195 
196  TStopwatch sw;
197  sw.Start();
198 
199  const double large = 1.0e300;
200  double minwgt = large, maxwgt = -large, maxenergy = -large;
201  std::set<int> gFlavors;
202 
203  for ( long int ientry = 0; ientry < nentries; ++ientry ) {
204 
205  // reset what's been read in anticipation of a new entry
206  entry_in->Reset();
207  aux_in->Reset();
208  numi_in->Reset();
209 
210  // read the next entry, get metadata if it's different
211  //int nbytes =
212  etreein->GetEntry(ientry);
213  UInt_t metakey = entry_in->metakey;
214  if ( fmeta->metakey != metakey ) {
215  // UInt_t oldkey = meta_in->metakey;
216  int nbmeta = mtreein->GetEntryWithIndex(metakey);
217  cout << "on entry " << ientry << " fetched metadata w/ key "
218  << metakey << " read " << nbmeta << " bytes" << endl;
219  }
220 
221  // reset what we're writing
222  fentry->Reset();
223  fnumi->Reset();
224  faux->Reset();
225 
226  // copy read in objects to output objects
227  *fentry = *entry_in;
228  *fnumi = *numi_in;
229  *faux = *aux_in;
230 
231  // mess with the values to your hearts content
232  warp_entry(fentry,faux,fnumi);
233  // keep track of any weights for the meta data
234  if ( fentry->wgt < minwgt ) minwgt = fentry->wgt;
235  if ( fentry->wgt > maxwgt ) maxwgt = fentry->wgt;
236  // user might have changed an energy larger than any in original file
237  if ( fentry->E > maxenergy ) maxenergy = fentry->E;
238  // user might have changed a flavor keep track of all we see
239  gFlavors.insert(fentry->pdg);
240 
241  // process currently held metadata after transition to metadata key
242  if ( fmeta->metakey != meta_in->metakey ) {
243  // new meta data found
244  cout << "new meta found " << fmeta->metakey << "/" << meta_in->metakey << endl;
245  if ( fmeta->metakey != 0 ) {
246  // have metadata that needs adjustment and writing
247  // before processing new metadata
248  warp_meta(fmeta,fnamein);
249  fmeta->minWgt = minwgt;
250  fmeta->maxWgt = maxwgt;
251  fmeta->maxEnergy = maxenergy;
252  update_flavors(fmeta,gFlavors);
253  gFlavors.clear();
254  cout << "metantp->Fill() " << *fmeta << endl;
255  metantp->Fill();
256  // next meta-data restarts weight range (and energy max)
257  minwgt = large;
258  maxwgt = -large;
259  maxenergy = -large;
260  }
261  // get a copy of metadata for accumulating adjustments
262  *fmeta = *meta_in;
263  //cout << " copy meta_in " << *meta_in << endl << *fmeta << endl;
264  }
265 
266  fluxntp->Fill();
267  }
268 
269  // write last set of meta-data after all is done
270  if ( fmeta->metakey != 0 ) {
271  cout << "final meta flush " << fmeta->metakey << endl;
272  // have metadata that needs adjustment and writing
273  warp_meta(fmeta,fnamein);
274  fmeta->minWgt = minwgt;
275  fmeta->maxWgt = maxwgt;
276  fmeta->maxEnergy = maxenergy;
277  update_flavors(fmeta,gFlavors);
278  gFlavors.clear();
279  cout << "metantp->Fill() " << *fmeta << endl;
280  metantp->Fill();
281  }
282 
283  cout << "=========================== Complete " << endl;
284 
285  sw.Stop();
286  cout << "Generated " << nentries
287  << endl
288  << "Time to generate: " << endl;
289  sw.Print();
290 
291  fout->cd();
292 
293  // write ntuples out
294  fluxntp->Write();
295  metantp->Write();
296  fout->Close();
297 
298  cout << endl << endl;
299 }
300 
301 void update_flavors(GSimpleNtpMeta* meta, std::set<int>& gFlavors)
302 {
303  // add any new flavors here, if necessary
304  std::set<int>::const_iterator flitr = gFlavors.begin();
305  std::vector<Int_t>& flvlist = meta->pdglist;
306  for ( ; flitr != gFlavors.end(); ++flitr ) {
307  int seen_pdg = *flitr;
308  std::vector<Int_t>::iterator knwitr =
309  std::find(flvlist.begin(),flvlist.end(),seen_pdg);
310  bool known_pdg = ( knwitr != flvlist.end() );
311  if ( ! known_pdg ) meta->pdglist.push_back(seen_pdg);
312  }
313 }
Double_t E
energy in lab frame
TString fin
Definition: Style.C:24
void warp_gsimple(string fnameout="gsimple_output.root", string fnamein="gsimple_input.root", bool dowarp=true)
Definition: warp_gsimple.C:130
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
Double_t px
x momentum in lab frame
TF1 * fcn
Definition: warp_gsimple.C:35
Double_t maxEnergy
maximum energy
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
Int_t par
Definition: SimpleIterate.C:24
std::vector< Int_t > auxint
additional ints associated w/ entry
const double emin
GENIE flux drivers.
std::vector< Int_t > pdglist
list of neutrino flavors
void update_flavors(GSimpleNtpMeta *meta, std::set< int > &gFlavors)
Definition: warp_gsimple.C:301
void warp_entry(GSimpleNtpEntry *entry, GSimpleNtpAux *aux, GSimpleNtpNuMI *numi)
Definition: warp_gsimple.C:44
const double emax
Long64_t nentries
double energy
Definition: plottest35.C:25
void warp_meta(GSimpleNtpMeta *meta, string fnamein)
Definition: warp_gsimple.C:81
Double_t minWgt
minimum weight
Class for parsing *.fcl files for their metadata information.
double pick_energy(double, double)
Definition: warp_gsimple.C:93
std::vector< std::string > infiles
list of input files
OStream cout
Definition: OStream.cxx:6
Double_t maxWgt
maximum weight
Float_t sw
Definition: plot.C:20
bool gDoWarp
Definition: warp_gsimple.C:39
UInt_t metakey
index key to tie to individual entries
double fit_fcn(double *x, double *par)
Definition: warp_gsimple.C:109
Double_t pz
z momentum in lab frame
Float_t e
Definition: plot.C:35
UInt_t metakey
key to meta data
Double_t py
y momentum in lab frame
def entry(str)
Definition: HTMLTools.py:26