gnumi2simple.C
Go to the documentation of this file.
2 #include "Tools/Flux/GNuMIFlux.h"
5 
6 #include "TSystem.h"
7 #include "TStopwatch.h"
8 #include "TLorentzVector.h"
9 #include "TNtuple.h"
10 #include "TFile.h"
11 
12 #include <iostream>
13 #include <iomanip>
14 #include <string>
15 #include <set>
16 
17 using namespace std;
18 using namespace genie;
19 using namespace genie::flux;
20 
21 // worker routine
22 void fill_simple(string fname="fluxntgenie.root",
23  string cfg="MINOS-Near",
24  int flxset=0,
25  long int nentries=0,
26  double pots=1.0e30,
27  double enumin=0.0, bool doaux=false);
28 
29 
30 // main routine
31 void gnumi2simple(string fname, string det, string flux,
32  long int nentries=5000, double pots=1.0e30,
33  double enumin=0.0, bool doaux=false)
34 {
35 
36  cout << "output file: " << fname << endl;
37  cout << "request " << nentries << " entries, or " << pots << " POTs" << endl;
38 
39  string cfg = "MINOS-NearDet";
40  if ( det == "minos" ) cfg = "MINOS-NearDet";
41  if ( det == "novaipnd" ) cfg = "NOvA-IPNDShed";
42  if ( det == "uboone" ) cfg = "microboone-numi";
43 
44  cout << "det \"" << det << "\" ==> cfg " << cfg << endl;
45 
46  int flxset = 0;
47  if ( flux == "old" ) flxset = 0;
48  if ( flux == "old1" ) flxset = 10;
49  if ( flux == "old2" ) flxset = 20;
50  if ( flux == "lowcut" ) flxset = 1;
51  if ( flux == "lowcut1" ) flxset = 11;
52  if ( flux == "lowcut2" ) flxset = 21;
53  if ( flux == "flugg" ) flxset = 2;
54  if ( flux == "flugg1" ) flxset = 12;
55  if ( flux == "flugg2" ) flxset = 22;
56  if ( flux == "flugglowth" ) flxset = 3;
57  if ( flux == "flugglowth1" ) flxset = 13;
58  if ( flux == "flugglowth2" ) flxset = 23;
59 
60  cout << "flux \"" << flux << "\" ==> flxset " << flxset << endl;
61 
62  fill_simple(fname,cfg,flxset,nentries,pots,enumin,doaux);
63 }
64 
65 void fill_simple(string fname, string cfg, int flxset,
66  long int nentries, double pots, double enumin, bool doaux)
67 {
68  string fluxpatt = "$FLUXPATH";
69 
70  int flxver = flxset % 10;
71  switch ( flxver ) {
72  case 0:
73  fluxpatt += "/v19/fluka05_le010z185i/root/fluka05_le010z185i_";
74  break;
75  case 1:
76  fluxpatt += "/v19/fluka05_le010z185i_lowcut/root/fluka05_le010z185i_";
77  break;
78  case 2:
79  fluxpatt += "/flugg/flugg_le010z185i_run1/root/flugg_le010z185i_run1_";
80  break;
81  case 3:
82  fluxpatt += "/flugg_lowth/flugg_le010z185i_run1_lowth/root/flugg_le010z185i_run1_lowth_";
83  break;
84  default:
85  cout << " not a legal flux set " << endl;
86  return;
87  }
88 
89  if ( flxset/10 == 1 ) {
90  if ( flxver == 2 || flxver == 3 ) fluxpatt += "001";
91  else fluxpatt += "1";
92  } else if ( flxset/10 == 2 ) {
93  if ( flxver == 2 || flxver == 3 ) fluxpatt += "00[12]";
94  else fluxpatt += "[12]";
95  } else fluxpatt += "*";
96  fluxpatt += ".root";
97 
98  string fluxfname(gSystem->ExpandPathName(fluxpatt.c_str()));
99  flux::GNuMIFlux* gnumi = new GNuMIFlux();
100  gnumi->LoadBeamSimData(fluxfname,cfg);
101  gnumi->SetEntryReuse(1); // don't reuse entries when reformatting
102 
103  //gnumi->PrintConfig();
104  //cout << " change to cm:" << endl;
105  //double cm_unit = genie::utils::units::UnitFromString("cm");
106  //gnumi->SetLengthUnits(cm_unit);
107  //gnumi->PrintConfig();
108  //cout << " change to m:" << endl;
109  ///////RWH: the following 2 lines should work ...
110  //double m_unit = genie::utils::units::UnitFromString("m");
111  //gnumi->SetLengthUnits(m_unit);
112  gnumi->PrintConfig();
113  //double scaleusr = 100.;
114 
115  gnumi->SetEntryReuse(1); // reuse entries
116  gnumi->SetUpstreamZ(-3e38); // leave ray on flux window
117 
118  GFluxI* fdriver = dynamic_cast<GFluxI*>(gnumi);
119 
120  if (nentries == 0) nentries = 100000*20;
121 
122  // so as not to include scan time generate 1 nu
123  fdriver->GenerateNext();
124 
125  TFile* file = TFile::Open(fname.c_str(),"RECREATE");
126  TTree* fluxntp = new TTree("flux","a simple flux n-tuple");
127  TTree* metantp = new TTree("meta","metadata for flux n-tuple");
132  fluxntp->Branch("entry",&fentry);
133  fluxntp->Branch("numi",&fnumi);
134  if (doaux) fluxntp->Branch("aux",&faux);
135  metantp->Branch("meta",&fmeta);
136 
137  TLorentzVector p4u, x4u;
138  //TLorentzVector p4b, x4b;
139  long int ngen = 0;
140 
141  set<int> pdglist;
142  double maxe = 0;
143  double minwgt = +1.0e10;
144  double maxwgt = -1.0e10;
145 
146  UInt_t metakey = TString::Hash(fname.c_str(),strlen(fname.c_str()));
147  cout << "metakey " << metakey << endl;
148  // ensure that we don't get smashed by UInt_t vs Int_t
149  metakey &= 0x7FFFFFFF;
150  cout << "metakey " << metakey << " after 0x7FFFFFFF" << endl;
151  cout << "=========================== Start " << endl;
152 
153  TStopwatch sw;
154  sw.Start();
155  int nok = 0, nwrite = 0;
156  while ( nok < nentries && gnumi->UsedPOTs() < pots ) {
157  //for (int i=0; i<nentries; ++i) {
158  fdriver->GenerateNext();
159  ngen++;
160  fentry->Reset();
161  fnumi->Reset();
162  faux->Reset();
163 
164  fentry->metakey = metakey;
165  fentry->pdg = fdriver->PdgCode();
166  fentry->wgt = gnumi->Weight();
167  x4u = fdriver->Position();
168  fentry->vtxx = x4u.X();
169  fentry->vtxy = x4u.Y();
170  fentry->vtxz = x4u.Z();
171  fentry->dist = gnumi->GetDecayDist();
172  p4u = fdriver->Momentum();
173  fentry->px = p4u.Px();
174  fentry->py = p4u.Py();
175  fentry->pz = p4u.Pz();
176  fentry->E = p4u.E();
177 
178  fnumi->run = gnumi->PassThroughInfo().run;
179  fnumi->evtno = gnumi->PassThroughInfo().evtno;
180  fnumi->entryno = gnumi->GetEntryNumber();
181 
182  fnumi->tpx = gnumi->PassThroughInfo().tpx;
183  fnumi->tpy = gnumi->PassThroughInfo().tpy;
184  fnumi->tpz = gnumi->PassThroughInfo().tpz;
185  fnumi->vx = gnumi->PassThroughInfo().vx;
186  fnumi->vy = gnumi->PassThroughInfo().vy;
187  fnumi->vz = gnumi->PassThroughInfo().vz;
188 
189  fnumi->ndecay = gnumi->PassThroughInfo().ndecay;
190  fnumi->ppmedium = gnumi->PassThroughInfo().ppmedium;
191  fnumi->tptype = gnumi->PassThroughInfo().tptype;
192 
193 
194  if ( doaux ) {
195  faux->auxint.push_back(gnumi->PassThroughInfo().tptype);
196  faux->auxdbl.push_back(gnumi->PassThroughInfo().tpx);
197  faux->auxdbl.push_back(gnumi->PassThroughInfo().tpy);
198  faux->auxdbl.push_back(gnumi->PassThroughInfo().tpz);
199  }
200 
201  if ( fentry->E > enumin ) ++nok;
202  fluxntp->Fill();
203  ++nwrite;
204 
205  // accumulate meta data
206  pdglist.insert(fentry->pdg);
207  minwgt = TMath::Min(minwgt,fentry->wgt);
208  maxwgt = TMath::Max(maxwgt,fentry->wgt);
209  maxe = TMath::Max(maxe,fentry->E);
210 
211  }
212  cout << "=========================== Complete " << endl;
213 
214  //fmeta->nflavors = pdglist.size();
215  //if ( fmeta->nflavors > 0 ) {
216  // fmeta->flavor = new Int_t[fmeta->nflavors];
217  // set<int>::const_iterator itr = pdglist.begin();
218  // int indx = 0;
219  // for ( ; itr != pdglist.end(); ++itr, ++indx ) fmeta->flavor[indx] = *itr;
220  //}
221  fmeta->pdglist.clear();
222  set<int>::const_iterator setitr = pdglist.begin();
223  for ( ; setitr != pdglist.end(); ++setitr) fmeta->pdglist.push_back(*setitr);
224 
225  fmeta->maxEnergy = maxe;
226  fmeta->minWgt = minwgt;
227  fmeta->maxWgt = maxwgt;
228  fmeta->protons = gnumi->UsedPOTs();
229  TVector3 p0, p1, p2;
230  gnumi->GetFluxWindow(p0,p1,p2);
231  TVector3 d1 = p1 - p0;
232  TVector3 d2 = p2 - p0;
233  fmeta->windowBase[0] = p0.X();
234  fmeta->windowBase[1] = p0.Y();
235  fmeta->windowBase[2] = p0.Z();
236  fmeta->windowDir1[0] = d1.X();
237  fmeta->windowDir1[1] = d1.Y();
238  fmeta->windowDir1[2] = d1.Z();
239  fmeta->windowDir2[0] = d2.X();
240  fmeta->windowDir2[1] = d2.Y();
241  fmeta->windowDir2[2] = d2.Z();
242  if ( doaux ) {
243  fmeta->auxintname.push_back("tptype");
244  fmeta->auxdblname.push_back("tpx");
245  fmeta->auxdblname.push_back("tpy");
246  fmeta->auxdblname.push_back("tpz");
247  }
248  string fname = "flux_pattern=" + fluxfname;
249  fmeta->infiles.push_back(fluxfname); // should get this expanded from gnumi
250  vector<string> flist = gnumi->GetFileList();
251  for (size_t i = 0; i < flist.size(); ++i) {
252  fname = "flux_infile=" + flist[i];
253  fmeta->infiles.push_back(fname);
254  }
255 
256  fmeta->seed = RandomGen::Instance()->GetSeed();
257  fmeta->metakey = metakey;
258 
259  metantp->Fill();
260 
261  sw.Stop();
262  cout << "Generated " << nwrite << " (" << nok << " w/ E >" << enumin << " ) "
263  << " ( request " << nentries << " ) "
264  << endl
265  << gnumi->UsedPOTs() << " POTs " << " ( request " << pots << " )"
266  << endl
267  << " pulled NFluxNeutrinos " << gnumi->NFluxNeutrinos()
268  << endl
269  << "Time to generate: " << endl;
270  sw.Print();
271 
272  cout << "===================================================" << endl;
273 
274  cout << "Last GSimpleNtpEntry: " << endl
275  << *fentry
276  << endl
277  << "Last GSimpleNtpMeta: " << endl
278  << *fmeta
279  << endl;
280 
281  file->cd();
282 
283  // write meta as simple object at top of file
284  //fmeta->Write();
285 
286  // write ntuples out
287  fluxntp->Write();
288  metantp->Write();
289  file->Close();
290 
291  cout << endl << endl;
292  gnumi->PrintConfig();
293 }
294 
295 /////////////////////// old interface /////////////////////////////
296 void make_simple_grid(string fname, string det, string flux,
297  long int nentries=5000, double pots=1.0e30,
298  double enumin=0.0, bool doaux=false)
299 {
300  gnumi2simple(fname,det,flux,nentries,pots,enumin,doaux);
301 }
Double_t E
energy in lab frame
void GetFluxWindow(TVector3 &p1, TVector3 &p2, TVector3 &p3) const
3 points define a plane in beam coordinate
Definition: GNuMIFlux.cxx:1060
void SetEntryReuse(long int nuse=1)
of times to use entry before moving to next
Definition: GNuMIFlux.cxx:924
void fill_simple(string fname="fluxntgenie.root", string cfg="MINOS-Near", int flxset=0, long int nentries=0, double pots=1.0e30, double enumin=0.0, bool doaux=false)
Definition: gnumi2simple.C:65
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
Double_t px
x momentum in lab frame
Long64_t GetEntryNumber()
index in chain
Definition: GNuMIFlux.h:253
int pots
Double_t vtxy
y position in lab frame
Double_t vx
vertex position of hadron/muon decay
Double_t maxEnergy
maximum energy
virtual const TLorentzVector & Position(void)=0
returns the flux neutrino 4-position (note: expect SI rather than physical units) ...
std::vector< std::string > auxintname
tagname of aux ints associated w/ entry
Int_t tptype
parent particle type at target exit
void gnumi2simple(string fname, string det, string flux, long int nentries=5000, double pots=1.0e30, double enumin=0.0, bool doaux=false)
Definition: gnumi2simple.C:31
std::vector< Double_t > auxdbl
additional doubles associated w/ entry
Double_t windowDir1[3]
dx,dy,dz of window direction 1
std::vector< Int_t > auxint
additional ints associated w/ entry
Loaders::FluxType flux
Double_t protons
represented number of protons-on-target
Int_t seed
random seed used in generation
virtual void SetUpstreamZ(double z0)
Double_t windowDir2[3]
dx,dy,dz of window direction 2
GENIE flux drivers.
A GENIE flux driver encapsulating the NuMI neutrino flux. It reads-in the official GNUMI neutrino flu...
Definition: GNuMIFlux.h:217
double Weight(void)
returns the flux neutrino weight (if any)
Definition: GNuMIFlux.h:234
virtual void LoadBeamSimData(const std::vector< std::string > &filenames, const std::string &det_loc)
Definition: GNuMIFlux.cxx:640
std::vector< Int_t > pdglist
list of neutrino flavors
Double_t windowBase[3]
x,y,z position of window base point
std::vector< std::string > GetFileList()
list of files currently part of chain
Definition: GNuMIFlux.cxx:2592
const GNuMIFluxPassThroughInfo & PassThroughInfo(void)
GNuMIFluxPassThroughInfo.
Definition: GNuMIFlux.h:252
Long64_t nentries
Double_t vtxz
z position in lab frame
Double_t minWgt
minimum weight
Int_t ppmedium
tracking medium where parent was produced
std::vector< std::string > infiles
list of input files
Double_t tpx
parent particle px at target exit
OStream cout
Definition: OStream.cxx:6
virtual bool GenerateNext(void)=0
generate the next flux neutrino (return false in err)
double GetDecayDist() const
dist (user units) from dk to current pos
Definition: GNuMIFlux.cxx:551
const double enumin
Definition: runWimpSim.h:36
Double_t maxWgt
maximum weight
Float_t sw
Definition: plot.C:20
std::size_t Hash(const std::vector< std::size_t > &vals)
Generate a unique hash from a sequence of integers.
Definition: Utilities.cxx:384
UInt_t metakey
index key to tie to individual entries
Double_t vtxx
x position in lab frame
Double_t pz
z momentum in lab frame
virtual int PdgCode(void)=0
returns the flux neutrino pdg code
virtual const TLorentzVector & Momentum(void)=0
returns the flux neutrino 4-momentum
TFile * file
Definition: cellShifts.C:17
void make_simple_grid(string fname, string det, string flux, long int nentries=5000, double pots=1.0e30, double enumin=0.0, bool doaux=false)
Definition: gnumi2simple.C:296
void PrintConfig()
print the current configuration
Definition: GNuMIFlux.cxx:2492
Double_t dist
distance from hadron decay
virtual long int NFluxNeutrinos() const
of rays generated
Definition: GNuMIFlux.cxx:288
UInt_t metakey
key to meta data
Double_t py
y momentum in lab frame
std::vector< std::string > auxdblname
tagname of aux doubles associated w/ entry
GENIE Interface for user-defined flux classes.
Definition: GFluxI.h:37
double UsedPOTs(void) const
of protons-on-target used
Definition: GNuMIFlux.cxx:620