convert_flugg.C
Go to the documentation of this file.
1 //
2 // Convert a flugg ntuple to the new format
3 //
4 // this script can be run using:
5 // $ root $(DK2NU)/snippets/load_dk2nu.C convert_flugg.C("file.root")
6 //
7 // rhatcher@fnal.gov 2012-11-06
8 //====================================================================
9 
10 #include <iostream>
11 #include <iomanip>
12 #include <string>
13 using namespace std;
14 
15 #include "TMath.h"
16 
17 #include "convert/common_convert.C"
18 #include "convert/flugg/flugg.C"
19 
20 void copy_flugg_to_dk2nu(const flugg& fluggObj, string volumeFilePath);
21 void fluggCrossChecks(const flugg& fluggObj, string inputloc);
22 
23 void convert_flugg(string ifname="../fluxfiles/generic_flugg.root",
24  int jobnum=42,
25  string volumeFilePath="../fluxfiles/g4numi001_Volumes_index.inp", // Fluka Volumes File
26  string inputloc="MINOS", // "MINOS"or "NOvA" for xcheck
27  Long64_t maxentries=-1,
28  Long64_t moddump=-1) // modulo for dump
29 {
30  // set globals
31  myjob = jobnum; // allow override because flugg files forgot to set this
32  pots = 0;
33  // allowance in location energy/weight cross-check
34  frac_diff_tolerance = 2.5e-4;
35 
36  int highest_potnum = 0;
37 
38  // open input ntuple
39  TFile* fin = TFile::Open(ifname.c_str());
40  if ( ! fin ) {
41  cout << "couldn't open input file: " << ifname << endl;
42  return;
43  }
44  const char* objName = "h10";
45  TTree* tin = 0;
46  fin->GetObject(objName,tin);
47 
48  if ( ! tin ) {
49  cout << "couldn't find " << objName << endl;
50  return;
51  }
52  cout << endl << "Input file: " << ifname << endl;
53 
54  // use makeclass interface for flugg
55  flugg fluggObj(tin);
56 
57  // construct output filename
58  string ofname = construct_outfilename(ifname);
59  cout << "Output file: " << ofname << endl << endl;
60 
61  ConvertInit();
62  ConvertBookNtuple(ofname);
63 
64  ///
65  /// loop over entries
66  ///
67  Long64_t nentries = fluggObj.fChain->GetEntriesFast();
68  if ( maxentries > 0 ) {
69  cout << "limit # of entries to " << nentries
70  << " / " << maxentries << endl;
71  nentries = TMath::Min(nentries,maxentries);
72  }
73  cout << endl;
74 
75  Long64_t nbytes = 0, nb = 0;
76  for (Long64_t jentry = 0; jentry < nentries; ++jentry) {
77  Long64_t ientry = fluggObj.LoadTree(jentry);
78  if (ientry < 0) break;
79  nb = fluggObj.fChain->GetEntry(jentry);
80  nbytes += nb;
81 
82  // always clear the dk2nu object
83  dk2nu->clear(); // !!! important !!! always do this
84 
85  // fill the dk2nu object from the flugg entry
86  copy_flugg_to_dk2nu(fluggObj, volumeFilePath);
87 
88  // fill location specific p3, energy and weights
89  // locations to fill are in the metadata
90  // assumes that prior copying filled the first entry w/ random decay
92 
93  // keep track of potnum
94  if ( dk2nu->potnum > highest_potnum )
95  highest_potnum = dk2nu->potnum;
96 
97  // push entry out to tree
98  dk2nuTree->Fill();
99 
100  // just for fun print every n entries
101  if ( moddump > 0 && jentry%moddump == 0 ) cout << endl << *dk2nu << endl;
102 
103  fluggCrossChecks(fluggObj,inputloc);
104 
105  }
106  cout << endl;
107 
108  pots = estimate_pots(highest_potnum);
109  cout << "estimated pots " << pots << " (" << highest_potnum << ")" << endl;
110  cout << "maximum observed fractional difference in cross checks: "
111  << obs_frac_diff_max << endl
112  << "print ## message when exceeded " << frac_diff_tolerance << endl;
113 
114  dkmeta->job = myjob;
115  dkmeta->pots = pots;
116  dkmeta->beamsim = "convert_flugg.C";
117  dkmeta->physics = "bogus";
118 
119  // print last entries
120  cout << endl << "========== last entries" << endl << endl;
121  cout << *dk2nu << endl << endl;
122  cout << *dkmeta << endl;
123 
124  // finish off ntuples (including meta data)
125  ConvertFinish();
126 
127  // close input file
128  fin->Close();
129 }
130 
131 void copy_flugg_to_dk2nu(const flugg& fluggObj, string volumeFilePath)
132 {
133  // fill the global dk2nu object
134 
135  dk2nu->job = myjob;
136  dk2nu->potnum = fluggObj.evtno;
137 
138  // calcLocationWeights needs random decay entries filled first
139  // but don't copy any other (i.e. near/far) values as they'll
140  // be recalculated for the whole set of locations
141  double pzrndm = fluggObj.Npz;
142  double pxrndm = fluggObj.Ndxdz * pzrndm;
143  double pyrndm = fluggObj.Ndydz * pzrndm;
144  bsim::NuRay nuray(pxrndm,pyrndm,pzrndm,fluggObj.Nenergy,1.0);
145  dk2nu->nuray.push_back(nuray);
146 
147  // calcLocationWeights needs these filled if it isn't going assert()
148  // really need to fill the other bits at this point as well:
149  // ntype, ptype, vx, vy, vz, pdpx, pdpy, pdpz, necm,
150  // ppenergy, ppdxdz, ppdydz, pppz,
151  // muparpx, muparpy, muparpz, mupare
152 
153  dk2nu->decay.norig = fluggObj.Norig;
154  dk2nu->decay.ndecay = fluggObj.Ndecay;
155  dk2nu->decay.ntype = Convert5xToPdg(fluggObj.Ntype);
156  dk2nu->decay.vx = fluggObj.Vx;
157  dk2nu->decay.vy = fluggObj.Vy;
158  dk2nu->decay.vz = fluggObj.Vz;
159  dk2nu->decay.pdpx = fluggObj.pdPx;
160  dk2nu->decay.pdpy = fluggObj.pdPy;
161  dk2nu->decay.pdpz = fluggObj.pdPz;
162  dk2nu->decay.ppdxdz = fluggObj.ppdxdz;
163  dk2nu->decay.ppdydz = fluggObj.ppdydz;
164  dk2nu->decay.pppz = fluggObj.pppz;
165  dk2nu->decay.ppenergy = fluggObj.ppenergy;
166  dk2nu->decay.ppmedium = fluggObj.ppmedium;
167  dk2nu->decay.ptype = ConvertGeantToPdg(fluggObj.ptype,"ptype");
168  dk2nu->decay.muparpx = fluggObj.muparpx;
169  dk2nu->decay.muparpy = fluggObj.muparpy;
170  dk2nu->decay.muparpz = fluggObj.muparpz;
171  dk2nu->decay.mupare = fluggObj.mupare;
172 
173  dk2nu->decay.necm = fluggObj.Necm;
174  dk2nu->decay.nimpwt = fluggObj.Nimpwt;
175 
176  dk2nu->ppvx = fluggObj.ppvx;
177  dk2nu->ppvy = fluggObj.ppvy;
178  dk2nu->ppvz = fluggObj.ppvz;
179 
180  //not-in-new//dk2nu->xpoint = fluggObj.xpoint;
181  //not-in-new//dk2nu->ypoint = fluggObj.ypoint;
182  //not-in-new//dk2nu->zpoint = fluggObj.zpoint;
183 
184  dk2nu->tgtexit.tvx = fluggObj.tvx;
185  dk2nu->tgtexit.tvy = fluggObj.tvy;
186  dk2nu->tgtexit.tvz = fluggObj.tvz;
187  dk2nu->tgtexit.tpx = fluggObj.tpx;
188  dk2nu->tgtexit.tpy = fluggObj.tpy;
189  dk2nu->tgtexit.tpz = fluggObj.tpz;
190  dk2nu->tgtexit.tptype = ConvertGeantToPdg(fluggObj.tptype,"tptype");
191  dk2nu->tgtexit.tgen = fluggObj.tgen;
192 
193  //no-equiv// dk2nu->tgptype = ConvertGeantToPdg(fluggObj.tgptype,"tgptype");
194  //no-equiv// dk2nu->tgppx = fluggObj.tgppx;
195  //no-equiv// dk2nu->tgppy = fluggObj.tgppy;
196  //no-equiv// dk2nu->tgppz = fluggObj.tgppz;
197  //no-equiv// dk2nu->tprivx = fluggObj.tprivx;
198  //no-equiv// dk2nu->tprivy = fluggObj.tprivy;
199  //no-equiv// dk2nu->tprivz = fluggObj.tprivz;
200 
201  //not-in-new//dk2nu->beamx = fluggObj.protonX;
202  //not-in-new//dk2nu->beamy = fluggObj.protonY;
203  //not-in-new//dk2nu->beamz = fluggObj.protonZ;
204  //not-in-new//dk2nu->beampx = fluggObj.protonPx;
205  //not-in-new//dk2nu->beampy = fluggObj.protonPy;
206  //not-in-new//dk2nu->beampz = fluggObj.protonPz;
207 
208  //no-equiv// dk2nu->tgptype = fluggObj.tgptype;
209 
210 
211 
212  // now copy ancestor history - added by Marco on 4-29-2015
213 
214  int nmx = fluggObj.ancestor; // number of intermediate levels
215  for (int ian=0; ian < nmx; ++ian) {
216 
218 
219  ancestor.pdg = ConvertGeantToPdg(fluggObj.track_pdg[ian],"ancestor_track_pdg");
220  ancestor.SetStartXYZT(fluggObj.startx[ian], // cm
221  fluggObj.starty[ian], // cm
222  fluggObj.startz[ian], // cm
223  fluggObj.startt[ian]); // cm
224  ancestor.SetStartP(fluggObj.startpx[ian], // GeV
225  fluggObj.startpy[ian], // GeV
226  fluggObj.startpz[ian]); // GeV
227  ancestor.SetStopP(fluggObj.stoppx[ian],
228  fluggObj.stoppy[ian],
229  fluggObj.stoppz[ian]);
230  //ancestor.SetPProdP(0,
231  // 0,
232  // 0);
233 
234  // ancestor.polx = ancestor.poly = ancestor.polz = ?
235  // ancestor.nucleau = ?
236 
237  TString procString, volumeString;
238  procString = ConvertFlukaInteractionCodeToString(fluggObj.proc[ian]);
239  ancestor.proc = procString.Data();
240  volumeString = ConvertVolumeCodeToString(fluggObj.ivol[ian],volumeFilePath);
241  ancestor.ivol = volumeString.Data();
242 
243  ancestor.parIndex = fluggObj.parIndex[ian];
244 
245  dk2nu->ancestor.push_back(ancestor);
246  }
247 
248 
249 }
250 
251 void fluggCrossChecks(const flugg& fluggObj, string inputloc)
252 {
253 
254  static bool first = true;
255  static size_t indxNear = -1, indxFar = -1;
256  if ( first ) {
257  first = false;
258  /// find indices for new weights so we can compare later
259  indxNear = find_loc_index(inputloc+" NearDet");
260  indxFar = find_loc_index(inputloc+" FarDet");
261  cout << "indx " << inputloc << " Near " << indxNear
262  << " Far " << indxFar << endl;
263  }
264  // cross check location energy/weights
265  //
266  histCompare(fluggObj.Nenergyn,dk2nu->nuray[indxNear].E, "near energy");
267  histCompare(fluggObj.Nwtnear, dk2nu->nuray[indxNear].wgt,"near wgt");
268  histCompare(fluggObj.Nenergyf,dk2nu->nuray[indxFar].E, "far energy");
269  histCompare(fluggObj.Nwtfar, dk2nu->nuray[indxFar].wgt, "far wgt");
270 
271 }
272 
Double_t Nwtfar
Definition: flugg.h:34
Double_t stoppz[10000]
Definition: flugg.h:113
TString fin
Definition: Style.C:24
Int_t ppmedium
Definition: flugg.h:48
Int_t ancestor
Definition: flugg.h:101
Double_t ppdydz
% direction of nu parent at its production point
Definition: dk2nu.h:139
Double_t startpx[10000]
Definition: flugg.h:108
Double_t muparpx
%
Definition: dk2nu.h:160
Int_t job
identifying job #
Definition: dk2nu.h:323
TString ConvertVolumeCodeToString(int volumeCode, string volumeFilePath)
Int_t parIndex
particle index, from nu (0), parent (1) ... to proton (n)
Definition: dk2nu.h:222
Double_t pppz
% z momentum of nu parent at its production point
Definition: dk2nu.h:140
void SetStopP(Double_t px, Double_t py, Double_t pz)
TTree * dk2nuTree
std::string ivol
name of the volume where the particle starts
Definition: dk2nu.h:225
void copy_flugg_to_dk2nu(const flugg &fluggObj, string volumeFilePath)
Double_t pdpx
% px momentum of nu parent at (vx,vy,vz)
Definition: dk2nu.h:133
void clear(const std::string &opt="")
reset everything
Double_t ppenergy
Definition: flugg.h:47
Int_t track_pdg[10000]
Definition: flugg.h:102
std::vector< bsim::NuRay > nuray
rays through detector fixed points
Definition: dk2nu.h:328
Int_t evtno
Definition: flugg.h:22
Double_t startz[10000]
Definition: flugg.h:106
Double_t Npz
Definition: flugg.h:25
Int_t potnum
proton # processed by simulation
Definition: dk2nu.h:324
std::string proc
name of the process that creates this particle
Definition: dk2nu.h:224
int pots
Int_t tgen
Definition: dk2nu.h:276
Double_t nimpwt
% cumulative importance weight prod to decay
Definition: dk2nu.h:166
Double_t Nenergyf
Definition: flugg.h:33
Double_t tpy
Definition: flugg.h:66
Double_t ppdxdz
Definition: flugg.h:44
Double_t ppvy
Definition: flugg.h:51
Int_t proc[10000]
Definition: flugg.h:115
Double_t necm
% nu energy in center-of-mass frame
Definition: dk2nu.h:165
double estimate_pots(int highest_potnum)
Double_t muparpx
Definition: flugg.h:53
bsim::Decay decay
basic decay information
Definition: dk2nu.h:327
int myjob
Double_t ppvz
production vertex z of nu parent
Definition: dk2nu.h:336
int ConvertGeantToPdg(int geant_code, std::string tag="?")
Double_t startpy[10000]
Definition: flugg.h:109
Double_t tvy
y position of nu ancestor as it exits target
Definition: dk2nu.h:270
string construct_outfilename(string infilename)
Int_t job
identifying job # (keep files distinct)
Definition: dkmeta.h:86
void SetStartXYZT(Double_t x, Double_t y, Double_t z, Double_t t)
set triplets
Double_t pdPz
Definition: flugg.h:43
Double_t pdPx
Definition: flugg.h:41
bsim::Dk2Nu * dk2nu
Double_t ppdxdz
% direction of nu parent at its production point
Definition: dk2nu.h:138
Double_t ppvy
production vertex y of nu parent
Definition: dk2nu.h:335
Double_t vz
% neutrino production vertex z
Definition: dk2nu.h:132
Double_t Ndxdz
Definition: flugg.h:23
bsim::TgtExit tgtexit
info about leaving the target
Definition: dk2nu.h:338
Double_t pots
protons-on-target
Definition: dkmeta.h:87
Int_t tptype
species of ancestor exiting the target
Definition: dk2nu.h:275
Int_t pdg
ancestor species
Definition: dk2nu.h:188
Double_t tvy
Definition: flugg.h:63
void ConvertFinish()
double frac_diff_tolerance
Double_t pdpz
% pz momentum of nu parent at (vx,vy,vz)
Definition: dk2nu.h:135
Double_t muparpy
%
Definition: dk2nu.h:161
Double_t ppvx
production vertex x of nu parent
Definition: dk2nu.h:334
Double_t muparpz
%
Definition: dk2nu.h:162
Double_t stoppx[10000]
Definition: flugg.h:111
Int_t ptype
Definition: flugg.h:49
Double_t tpz
z momentum of nu ancestor as it exits target
Definition: dk2nu.h:274
Long64_t nentries
void SetStartP(Double_t px, Double_t py, Double_t pz)
std::string physics
e.g. "fluka08", "g4.9.3p01"
Definition: dkmeta.h:96
virtual Long64_t LoadTree(Long64_t entry)
Double_t ppenergy
% energy of nu parent at its production point
Definition: dk2nu.h:141
Double_t tvx
Definition: flugg.h:62
Double_t Nenergy
Definition: flugg.h:26
Double_t startt[10000]
Definition: flugg.h:107
std::string beamsim
e.g. "flugg" or "g4numi/<tag>"
Definition: dkmeta.h:95
Double_t tpy
y momentum of nu ancestor as it exits target
Definition: dk2nu.h:273
Int_t Ndecay
Definition: flugg.h:36
Int_t ntype
% neutrino flavor (PDG? code)
Definition: dk2nu.h:128
Double_t Nwtnear
Definition: flugg.h:30
void ConvertBookNtuple(std::string ofname="test_dk2nu.root")
Double_t pppz
Definition: flugg.h:46
Double_t Vy
Definition: flugg.h:39
Definition: flugg.h:15
Double_t ppvz
Definition: flugg.h:52
Double_t muparpz
Definition: flugg.h:55
Double_t startpz[10000]
Definition: flugg.h:110
Double_t starty[10000]
Definition: flugg.h:105
double obs_frac_diff_max
void convert_flugg(string ifname="../fluxfiles/generic_flugg.root", int jobnum=42, string volumeFilePath="../fluxfiles/g4numi001_Volumes_index.inp", string inputloc="MINOS", Long64_t maxentries=-1, Long64_t moddump=-1)
Definition: convert_flugg.C:23
Double_t Necm
Definition: flugg.h:57
Int_t tptype
Definition: flugg.h:68
OStream cout
Definition: OStream.cxx:6
Double_t mupare
Definition: flugg.h:56
Double_t Nenergyn
Definition: flugg.h:29
Int_t Norig
Definition: flugg.h:35
Double_t Vx
Definition: flugg.h:38
Double_t tpx
Definition: flugg.h:65
Double_t Ndydz
Definition: flugg.h:24
Double_t tvz
z position of nu ancestor as it exits target
Definition: dk2nu.h:271
Double_t stoppy[10000]
Definition: flugg.h:112
bool histCompare(double a, double b, string s)
bsim::DkMeta * dkmeta
Double_t startx[10000]
Definition: flugg.h:104
Long64_t nbytes
Double_t ppdydz
Definition: flugg.h:45
Double_t ppvx
Definition: flugg.h:50
int Convert5xToPdg(int old_ntype)
Int_t parIndex[10000]
Definition: flugg.h:114
size_t find_loc_index(string match)
void fluggCrossChecks(const flugg &fluggObj, string inputloc)
Double_t tpx
x momentum of nu ancestor as it exits target
Definition: dk2nu.h:272
Double_t muparpy
Definition: flugg.h:54
Double_t pdPy
Definition: flugg.h:42
Double_t vx
% neutrino production vertex x
Definition: dk2nu.h:130
Double_t tvz
Definition: flugg.h:64
void calcLocationWeights(const bsim::DkMeta *dkmeta, bsim::Dk2Nu *dk2nu)
user interface
Int_t ndecay
decay process (see dkproc_t)
Definition: dk2nu.h:127
Double_t mupare
% energy of nu grandparent
Definition: dk2nu.h:163
Double_t vy
% neutrino production vertex y
Definition: dk2nu.h:131
Int_t tgen
Definition: flugg.h:69
Int_t Ntype
Definition: flugg.h:37
static const double nb
Definition: Units.h:89
Int_t ivol[10000]
Definition: flugg.h:116
Double_t tpz
Definition: flugg.h:67
Int_t norig
not used?
Definition: dk2nu.h:126
Int_t ppmedium
material nu parent was produced in
Definition: dk2nu.h:143
TString ConvertFlukaInteractionCodeToString(int interactionCode)
Double_t Vz
Definition: flugg.h:40
Double_t tvx
x position of nu ancestor as it exits target
Definition: dk2nu.h:269
std::vector< bsim::Ancestor > ancestor
chain from proton to neutrino
Definition: dk2nu.h:329
Double_t pdpy
% py momentum of nu parent at (vx,vy,vz)
Definition: dk2nu.h:134
TTree * fChain
Definition: flugg.h:17
Definition: nuray.py:1
Int_t ptype
% nu parent species (PDG? code)
Definition: dk2nu.h:144
void ConvertInit(std::string locfilename="$(DK2NU)/etc/locations.txt")
Double_t Nimpwt
Definition: flugg.h:58