helperFunctions.cxx
Go to the documentation of this file.
2 
3 
4 //-------------------------------------
5 
6 // TODO I think maybe this whole thing needs to be restructured as a class
7 
8 
9 uint GetTrajPointIdx( TLorentzVector daughter_pos,
10  std::vector<TLorentzVector> traj_points)
11 // Returns index of traj_points vector whose corresponding element
12 // is smallest distance from daughter_pos.This is a way of using daughter
13 // process info to set primary processes.
14 //
15 // We set a threshold. If the separation between the points is less than the
16 // threshold, the association is made.
17 
18 // This is hacky, but the reason it is necessary is that (for some reason)
19 // particle->Trajectory()->TrajectoryProcesses() is not filled in our
20 // simulation
21 {
22 
23  double thresh = 0.001;
24 
25  TLorentzVector elem;
26  double elem_X;
27  double elem_Y;
28  double elem_Z;
29 
30  double daughter_X = daughter_pos(0);
31  double daughter_Y = daughter_pos(1);
32  double daughter_Z = daughter_pos(2);
33 
34  double deltaX;
35  double deltaY;
36  double deltaZ;
37 
38  double length;
39 
40  uint numPoints = traj_points.size();
41 
42  uint traj_idx = 12345; // Some ridiculous sentinel value
43  // Can't use -5 as it is unsigned int
44 
45  int min_len = 100.0; // Some large initial value
46 
47  for(uint point_idx=0; point_idx < numPoints; point_idx++){
48 
49  elem = traj_points.at(point_idx);
50  elem_X = elem(0);
51  elem_Y = elem(1);
52  elem_Z = elem(2);
53 
54  deltaX = elem_X - daughter_X;
55  deltaY = elem_Y - daughter_Y;
56  deltaZ = elem_Z - daughter_Z;
57 
58  length = len(deltaX, deltaY, deltaZ);
59 
60  if(length < min_len){
61  min_len = length;
62  }
63  else{
64  continue;
65  }
66 
67  if(length < thresh){
68  /*
69  if(traj_idx != 12345){ // Has traj_idx already been set
70 
71  std::cout << "***********************************" << std::endl;
72  std::cout << "WARNING - two points below threshold" << std::endl;
73  std::cout << "Ambigous association in GetTrajPointIdx() " << std::endl;
74  std::cout << "***********************************" << std::endl;
75 
76  }
77  */
78  traj_idx = point_idx;
79 
80  }
81  }
82 
83  return traj_idx;
84 }
85 
86 double len(double x1, double x2, double x3)
87 {
88 
89  // Returns length of vector
90  double length = sqrt(
91  std::pow( x1, 2 ) +
92  std::pow( x2, 2 ) +
93  std::pow( x3, 2 )
94  );
95 
96  return length;
97 
98 }
99 
101 // does this particle exit the nucleus and get handed to geant4
102 {
103  std::string process = part->Process();
104 
105  if(strcmp( process.c_str(), "Primary") == 0){
106  return true;
107  }
108 
109  return false;
110 }
111 
112 
113 G4ReweightTraj makeTraj(const sim::Particle* part,
114  art::Event& evt,
116 // Constructs trajectory object geant4reweight needs
117 {
118 
119 
120  int event = evt.id().event();
121 
122  // ~~~ Assign values to variables ~~~
123 
124  int prim_pdg = part->PdgCode();
125 
126  int prim_id = part->TrackId();
127 
128  int prim_nDaughters = part->NumberDaughters();
129 
130  // Don't know why these two are needed, but G4ReweightTraj needs them
131  // The values I'm passing I assume do nothing
132  int parent_id = 0;
133  std::pair<int,int> range = std::make_pair(0,0);
134 
135  double mass;
136  mass = part->Mass();
137  mass *= 1.e3; // geant4reweight wants mass in MeV
138 
139 
140  //Make a G4ReweightTraj -- This is the reweightable object
141  G4ReweightTraj theTraj(prim_id, prim_pdg, parent_id,
142  event, range);
143 
144  uint nTrajPoints = part->NumberTrajectoryPoints();
145  std::vector<TLorentzVector> trajPoints;
146 
147  // Process string at each point in trajectory. Initialise values to "none"
148  std::string processes[nTrajPoints];
149  std::fill_n(processes, nTrajPoints, "none");
150 
151  std::vector<double> prim_X;
152  std::vector<double> prim_Y;
153  std::vector<double> prim_Z;
154 
155  std::vector<double> prim_Px;
156  std::vector<double> prim_Py;
157  std::vector<double> prim_Pz;
158 
159  //art::ServiceHandle<geo::Geometry> geom; //TODO probs bad to declare
160  // service handle in a
161  // function which will
162  // be called multiple times
163 
164  for( uint step_idx = 0; step_idx < nTrajPoints;
165  step_idx++){
166 
167  trajPoints.push_back( part->Position(step_idx) );
168 
169  double X = part->Position(step_idx).X();
170  double Y = part->Position(step_idx).Y();
171  double Z = part->Position(step_idx).Z();
172 
173  /*
174  const TGeoMaterial* testmaterial1 = geom->Material( X,Y,Z );
175 
176  std::cout << "At point (" << X << ", " << Y << ", " << Z << ") " << testmaterial1->GetName()
177  << " " << testmaterial1->GetZ() << " " << testmaterial1->GetA()
178  << " " << testmaterial1->GetDensity() << std::endl;
179  */
180 
181  prim_X.push_back(X);
182  prim_Y.push_back(Y);
183  prim_Z.push_back(Z);
184 
185  prim_Px.push_back(part->Px(step_idx));
186  prim_Py.push_back(part->Py(step_idx));
187  prim_Pz.push_back(part->Pz(step_idx));
188 
189  } // End loop over trajectory points
190 
191  const TLorentzVector prim_start = part->Position(0);
192  const TLorentzVector prim_end = part->Position(nTrajPoints -1);
193 
194  for(int i=0; i < prim_nDaughters; i++){
195 
196  const int daughter_idx = part->Daughter(i);
197 
198  const sim::Particle* daughter =
199  bt->TrackIDToParticle(daughter_idx);
200 
201  const int daughter_pdg = daughter->PdgCode();
202  std::string daughter_proc = daughter->Process();
203 
204  theTraj.AddChild(new G4ReweightTraj(daughter_idx, daughter_pdg,
205  prim_id, event, range) );
206 
207 
208  const TLorentzVector daughter_start = daughter->Position(0);
209 
210  //****************
211  uint trajPointIdx = GetTrajPointIdx(daughter_start, trajPoints);
212  // TODO: There is redundancy here as e.g. at the end of the trajectory
213  // where there may be many daughters at the same point we are
214  // setting the process at that point once for each daughter.
215  if(trajPointIdx != 12345){
216  processes[trajPointIdx] = daughter->Process();
217  }
218  //****************
219 
220 
221  }// end loop over pion daughters
222 
223  std::string endProc = processes[nTrajPoints -1];
224 
225  if( strcmp( endProc.c_str(), "hadElastic") == 0 ){
226  // What is anyone supposed to do if they see this warning ??
227  //std::cout << "WARNING: last step and elastic" << std::endl;
228  }
229 
230 
231  // Loop where we add steps to traj
232  // This could possibly go in an existing loop for the sake of
233  // efficiency ??
234  for(uint idx=1; idx < nTrajPoints; idx++){
235 
236  double deltaX = ( prim_X.at(idx) - prim_X.at(idx -1) );
237  double deltaY = ( prim_Y.at(idx) - prim_Y.at(idx -1) );
238  double deltaZ = ( prim_Z.at(idx) - prim_Z.at(idx -1) );
239 
240  double length = len(deltaX, deltaY, deltaZ);
241 
242  double preStepP[3] = {
243  prim_Px.at(idx-1)*1.e3, // geant4reweight wants Momentum in MeV
244  prim_Py.at(idx-1)*1.e3,
245  prim_Pz.at(idx-1)*1.e3
246  };
247 
248  double postStepP[3] = {
249  prim_Px.at(idx)*1.e3,
250  prim_Py.at(idx)*1.e3,
251  prim_Pz.at(idx)*1.e3
252  };
253 
254 
255  if(idx == 1){
256 
257  double momentum = sqrt( std::pow(preStepP[0], 2) +
258  std::pow(preStepP[1], 2) +
259  std::pow(preStepP[2], 2)
260  );
261 
262 
263  // Set initial energy
264  double energy = sqrt( std::pow(momentum, 2) +
265  std::pow(mass, 2)
266  );
267 
268  theTraj.SetEnergy(energy);
269  }
270 
271  std::string proc = processes[idx];
272 
273 
274 
275  // Events where a single geant4 step is huge
276  // seem to spit out huge weights, so lets skip those
277  //
278  // TODO this could be an area of further investigation
279  // re: non-unitarity
280  if(length > 250.0) continue;
281 
282 
283  G4ReweightStep * theStep = new G4ReweightStep( prim_id, prim_pdg,
284  parent_id, event,
285  preStepP, postStepP,
286  length, proc );
287 
288  theStep->SetDeltaX(deltaX);
289  theStep->SetDeltaY(deltaY);
290  theStep->SetDeltaZ(deltaZ);
291 
292  theTraj.AddStep(theStep);
293 
294  } // end loop over trajPoints
295 
296  return theTraj;
297 }
298 
299 bool RockFilter_A(std::vector<const sim::Particle*> parts,
301 // Return true if any of the particles enter the detector,
302 // otherwise return false
303 {
304  for(uint i=0; i<parts.size(); i++){
305  const sim::Particle *part = parts.at(i);
306  uint nTrajPoints = part->NumberTrajectoryPoints();
307  for(uint j=0; j<nTrajPoints; j++){
308  if(geom->isInsideDetectorBigBox( part->Position(j).X(),
309  part->Position(j).Y(),
310  part->Position(j).Z() )
311  ){
312  return true;
313  }
314 
315  }// end for(trajPoints)
316 
317  }// end for(parts)
318 
319  return false;
320 }
321 
322 
323 bool RockFilter_B(std::vector<const sim::Particle*> parts,
325 // Return true if any of the particles have endpoint in the detector,
326 // otherwise return false
327 {
328  for(uint i=0; i<parts.size(); i++){
329  const sim::Particle *part = parts.at(i);
330  uint nTrajPoints = part->NumberTrajectoryPoints();
331 
332  if(geom->isInsideDetectorBigBox( part->Position(nTrajPoints -1).X(),
333  part->Position(nTrajPoints -1).Y(),
334  part->Position(nTrajPoints -1).Z() )
335  ){
336  return true;
337  }
338 
339 
340 
341  }// end for(parts)
342 
343  return false;
344 }
345 
348 // Return true if interaction vertex is inside detector
349 {
350  const TVector3 vtx = mctruth->GetNeutrino().Nu().Position().Vect();
351 
352  if(geom->isInsideDetectorBigBox(vtx.X(), vtx.Y(), vtx.Z())){
353  return true;
354  }
355  return false;
356 
357 }
358 
359 bool RockFilter_D(std::vector<const sim::Particle*> parts,
361 // Return true if all pions stay within detector,
362 // otherwise return false
363 {
364  for(uint i=0; i<parts.size(); i++){
365  const sim::Particle *part = parts.at(i);
366 
367  if(abs(part->PdgCode()) != 211) continue;
368 
369  uint nTrajPoints = part->NumberTrajectoryPoints();
370  for(uint j=0; j<nTrajPoints; j++){
371  if(!geom->isInsideDetectorBigBox( part->Position(j).X(),
372  part->Position(j).Y(),
373  part->Position(j).Z() )
374  ){
375  return false;
376  }
377 
378  }// end for(trajPoints)
379 
380  }// end for(parts)
381 
382  return true;
383 }
384 
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:217
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:218
int PdgCode() const
Definition: MCParticle.h:211
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:76
double Py(const int i=0) const
Definition: MCParticle.h:230
G4ReweightTraj makeTraj(const sim::Particle *part, art::Event &evt, art::ServiceHandle< cheat::BackTracker > bt)
bool isPrimary(const sim::Particle *part)
Float_t x1[n_points_granero]
Definition: compare.C:5
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
T sqrt(T number)
Definition: d0nt_math.hpp:156
double Mass() const
Definition: MCParticle.h:238
constexpr T pow(T x)
Definition: pow.h:75
double Px(const int i=0) const
Definition: MCParticle.h:229
void abs(TH1 *hist)
bool RockFilter_C(art::Ptr< simb::MCTruth > mctruth, art::ServiceHandle< geo::Geometry > geom)
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
std::string Process() const
Definition: MCParticle.h:214
int NumberDaughters() const
Definition: MCParticle.h:216
int TrackId() const
Definition: MCParticle.h:209
int Daughter(const int i) const
Definition: MCParticle.cxx:112
Float_t Y
Definition: plot.C:38
Float_t Z
Definition: plot.C:38
bool RockFilter_D(std::vector< const sim::Particle * > parts, art::ServiceHandle< geo::Geometry > geom)
TString part[npart]
Definition: Style.C:32
length
Definition: demo0.py:21
bool RockFilter_B(std::vector< const sim::Particle * > parts, art::ServiceHandle< geo::Geometry > geom)
double len(double x1, double x2, double x3)
int evt
double energy
Definition: plottest35.C:25
bool isInsideDetectorBigBox(const double x_cm, const double y_cm, const double z_cm) const
Is the particle inside the detector Big Box?
const double j
Definition: BetheBloch.cxx:29
bool RockFilter_A(std::vector< const sim::Particle * > parts, art::ServiceHandle< geo::Geometry > geom)
const sim::Particle * TrackIDToParticle(int const &id) const
Returns a pointer to the sim::Particle object corresponding to the given TrackID. ...
EventNumber_t event() const
Definition: EventID.h:116
void geom(int which=0)
Definition: geom.C:163
double Pz(const int i=0) const
Definition: MCParticle.h:231
uint GetTrajPointIdx(TLorentzVector daughter_pos, std::vector< TLorentzVector > traj_points)
Float_t X
Definition: plot.C:38
EventID id() const
Definition: Event.h:56
enum BeamMode string
unsigned int uint