prepare_mcnp_root_for_g4nova.C
Go to the documentation of this file.
1 // Simplifies an MCNP-derived root file (convert_mcnp_text_to_root.C) into a
2 // TTree of primaries and their associated immediate daughters, in the format
3 // that g4nova/NeutronSubstitutionPhysics accepts.
4 
5 #include "TChain.h"
6 #include "TFile.h"
7 #include "TTree.h"
8 #include "TLorentzVector.h"
9 #include "TVector3.h"
10 
11 #include <iostream>
12 
13 // NB make sure to protect wildcard from the shell
14 void prepare_mcnp_root_for_g4nova(TString wildcard, bool byParentID = false)
15 {
16  TChain ch("tr");
17  ch.Add(wildcard);
18 
19  // Setup to read fields in
20  float x, y, z;
21  ch.SetBranchAddress("x", &x);
22  ch.SetBranchAddress("y", &y);
23  ch.SetBranchAddress("z", &z);
24 
25  float px, py, pz;
26  ch.SetBranchAddress("px", &px);
27  ch.SetBranchAddress("py", &py);
28  ch.SetBranchAddress("pz", &pz);
29 
30  float t;
31  ch.SetBranchAddress("t", &t);
32 
33  int pdg, evt, trk, parent;
34  ch.SetBranchAddress("pdg", &pdg);
35  ch.SetBranchAddress("evt", &evt);
36  ch.SetBranchAddress("trk", &trk);
37  ch.SetBranchAddress("parent", &parent);
38 
39 
40  TFile fout("fates.root", "RECREATE");
41  TTree trout("tr", "tr");
42 
43  // Setup to write fields out
44  float p, dt, stepz, stepx;
45  trout.Branch("p", &p);
46  trout.Branch("dt", &dt);
47  trout.Branch("stepz", &stepz);
48  trout.Branch("stepx", &stepx);
49 
50  // Child fields
51  int nchildren;
52  trout.Branch("nchildren", &nchildren);
53  int childpdg[10000];
54  trout.Branch("childpdg", &childpdg, "childpdg[nchildren]/I");
55  float childx[10000], childy[10000], childz[10000], childt[10000];
56  trout.Branch("childx", &childx, "childx[nchildren]/F");
57  trout.Branch("childy", &childy, "childy[nchildren]/F");
58  trout.Branch("childz", &childz, "childz[nchildren]/F");
59  trout.Branch("childt", &childt, "childt[nchildren]/F");
60  float childpx[10000], childpy[10000], childpz[10000];
61  trout.Branch("childpx", &childpx, "childpx[nchildren]/F");
62  trout.Branch("childpy", &childpy, "childpy[nchildren]/F");
63  trout.Branch("childpz", &childpz, "childpz[nchildren]/F");
64 
65 
66  long row = 0;
67  const long maxRow = ch.GetEntries();
68 
69  int nfalsepos = 0;
70 
71  while(true){
72  if(row%1000 == 0)
73  std::cout << row/1000000 << "m / " << maxRow/1000000 << "m" << std::endl;
74 
75  // Read the initial information for the parent neutron
76  ch.GetEntry(row++);
77  if(row == maxRow) goto end;
78  const TVector3 parentStart(x, y, z); // mm
79  const TVector3 parentMom(px, py, pz); // MeV
80  const float parentT = t; // ns
81  const int curEvt = evt;
82 
83  // We'll loop through lines looking for the first that isn't trackID #1,
84  // and call the details from the last line before that the position/time
85  // the neutron interacted.
86  TVector3 parentStop = parentStart;
87  float parentStopT = parentT;
88 
89  std::vector<TLorentzVector> parentTraj;
90 
91  while(true){
92  ch.GetEntry(row++);
93  if(row == maxRow) goto end;
94  if(trk != 1) break; // new track
95  if(evt != curEvt) break; // new event (no daughters)
96  parentTraj.emplace_back(x, y, z, t);
97  // Still the primary, so update the stop details
98  parentStop = TVector3(x, y, z);
99  parentStopT = t;
100  } // end search for end of primary
101 
102  --row; // back up so we'll re-read the first row belonging to a child next
103 
104  const TVector3 step = parentStop-parentStart;
105 
106  // z is defined as the initial momentum direction
107  const TVector3 zUnit = parentMom.Unit();
108  // x is defined as the part of the motion orthogonal to that
109  TVector3 xUnit = (step - step.Dot(zUnit)*zUnit).Unit();
110  if(fabs(xUnit.Mag()-1) > .01){
111  // Probably the motion transverse is zero, we need a different way to
112  // construct a basis.
113  xUnit = TVector3(1, 0, 0).Cross(zUnit).Unit();
114  }
115  // and y is orthogonal to both
116  const TVector3 yUnit = zUnit.Cross(xUnit);
117 
118  // Setup parent part of the output
119  p = parentMom.Mag();
120  dt = parentStopT - parentT;
121  stepz = step.Dot(zUnit);
122  stepx = step.Dot(xUnit);
123 
124  // Loop through following rows accumulating daughter tracks of this initial
125  // interaction.
126  int curTrk = 0;
127  nchildren = 0;
128  while(true){
129  ch.GetEntry(row++);
130  if(row == maxRow) goto end;
131  if(evt != curEvt) break; // run out of daughters, move on to next event
132 
133  // New track, check if it's directly produced by the primary
134  if(trk != curTrk){
135  const TLorentzVector childStart(x, y, z, t);
136 
137  bool match = false;
138  for(const TLorentzVector& v: parentTraj){
139  if(childStart == v) match = true;
140  }
141 
142  // If this is a g4 file where we have an independent record of the
143  // parentage, do some cross-checks.
144  if(parent != -1){
145  if(parent == 1 && !match){
146  std::cout << "Warning: particle known to be a child of the neutron had nonzero seperations to all trajectory points!" << std::endl;
147  }
148  if(parent != 1 && match){
149  ++nfalsepos;
150  // std::cout << "Warning: particle known not to be a child of the neutron had nonzero seperation to a trajectory point! (" << nfalsepos << " / " << trout.GetEntries() << ", parent = " << parent << ")" << std::endl;
151  }
152  }
153 
154  // Override the matching decision by just looking at the parent ID if
155  // requested (this will only work for g4).
156  if(byParentID) match = (parent == 1);
157 
158  if(match){
159  childpdg[nchildren] = pdg;
160 
161  const TVector3 childMom(px, py, pz);
162  childpx[nchildren] = childMom.Dot(xUnit);
163  childpy[nchildren] = childMom.Dot(yUnit);
164  childpz[nchildren] = childMom.Dot(zUnit);
165 
166  childx[nchildren] = childStart.Vect().Dot(xUnit);
167  childy[nchildren] = childStart.Vect().Dot(yUnit);
168  childz[nchildren] = childStart.Vect().Dot(zUnit);
169  childt[nchildren] = childStart.T() - parentT;
170 
171  ++nchildren;
172  }
173 
174  // Move on to looking for the next track
175  curTrk = trk;
176  } // end processing new track
177  } // event loop over tracks
178 
179  trout.Fill();
180 
181  --row; // back up one so we will read the parent particle again
182  } // end loop over events
183 
184  end:
185  trout.Write();
186 }
void prepare_mcnp_root_for_g4nova(TString wildcard, bool byParentID=false)
const char * p
Definition: xmltok.h:285
Track finder for cosmic rays.
Definition: Cand.cxx:23
std::vector< std::string > wildcard(const std::string &wildcardString)
Definition: convert.C:9
int evt
TVector3 Unit() const
z
Definition: test.py:28
OStream cout
Definition: OStream.cxx:6