Jpp  18.2.1
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JLorentzBoost.hh
Go to the documentation of this file.
1 #ifndef __JAANET__JLORENTZBOOST__
2 #define __JAANET__JLORENTZBOOST__
3 
4 #include <vector>
5 
10 
11 #include "JLang/JException.hh"
12 
13 #include "JGeometry3D/JVector3D.hh"
14 #include "JGeometry3D/JVertex3D.hh"
15 #include "JGeometry3D/JTrack3E.hh"
16 
17 #include "JAAnet/JAAnetToolkit.hh"
18 
19 #include "JPhysics/JConstants.hh"
20 
21 #include "Math/Vector4D.h"
22 #include "Math/GenVector/Boost.h"
23 #include "Math/GenVector/LorentzVector.h"
24 
25 
26 /**
27  * \file
28  *
29  * Auxiliary methods for calculating Lorentz boosts.
30  * \author bjjung
31  */
32 
33 namespace JAANET {}
34 namespace JPP { using namespace JAANET; }
35 
36 namespace JAANET {
37 
38  using ROOT::Math::Boost;
39  using ROOT::Math::LorentzVector;
40 
41 
42  /**
43  * Auxiliary class for performing Lorentz boosts.
44  */
45  struct JLorentzBoost :
46  public Boost
47  {
48  /**
49  * Default constructor.
50  */
52  {}
53 
54 
55  /**
56  * Copy constructor.
57  *
58  * \param boost Lorentz boost
59  */
60  JLorentzBoost(const Boost& boost) :
61  Boost(boost)
62  {}
63 
64 
65  /**
66  * Retrieve gamma factor.
67  *
68  * \return gamma factor
69  */
70  double getGamma() const
71  {
72  const double beta2 = this->BetaVector().Mag2();
73 
74  return 1.0 / sqrt(1.0 - beta2);
75  }
76 
77 
78  /**
79  * Lorentz boost operator.
80  *
81  * \param x4 4-vector
82  * \return boosted 4-vector
83  */
84  template<class CoordSystem>
85  LorentzVector<CoordSystem> operator()(const LorentzVector<CoordSystem>& x4) const
86  {
87  return static_cast<const Boost&>(*this)(x4);
88  }
89 
90 
91  /**
92  * Lorentz boost operator.
93  *
94  * \param vertex event vertex
95  */
96  void operator()(JVertex3D& vertex) const
97  {
98  using namespace JPP;
99  using namespace ROOT::Math;
100 
101  const XYZTVector x0(vertex.getX(), vertex.getY(), vertex.getZ(), vertex.getT() * getSpeedOfLight());
102  const XYZTVector x1 = (*this)(x0);
103 
104  const JPosition3D pos1(x1.X(), x1.Y(), x1.Z());
105 
106  vertex.setPosition(pos1);
107  vertex.setT(x1.T() / getSpeedOfLight());
108  }
109 
110 
111  /**
112  * Lorentz boost operator.
113  *
114  * \param track track
115  * \param mass track particle mass
116  */
117  void operator()(JTrack3E& track,
118  const double mass) const
119  {
120  using namespace JPP;
121  using namespace ROOT::Math;
122 
123  const double Ekin = getKineticEnergy(track.getE(), mass);
124 
125  const XYZTVector x0(track.getX(), track.getY(), track.getZ(), track.getT() * getSpeedOfLight());
126  const XYZTVector x1 = (*this)(x0);
127 
128  const PxPyPzEVector p0(Ekin * track.getDX(), Ekin * track.getDY(), Ekin * track.getDZ(), track.getE());
129  const PxPyPzEVector p1 = (*this)(p0);
130 
131  const JVector3D pos1(x1.X(), x1.Y(), x1.Z());
132  const JVersor3D dir1(p1.X(), p1.Y(), p1.Z());
133 
134  track.setPosition (pos1);
135  track.setDirection(dir1);
136  track.setT(x1.T());
137  track.setE(p1.E());
138  }
139 
140 
141  /**
142  * Lorentz boost operator.
143  * N.B.: Only initial and final state tracks are boosted.
144  *
145  * \param track track
146  */
147  void operator()(Trk& track) const
148  {
149  using namespace std;
150  using namespace JPP;
151  using namespace ROOT::Math;
152 
153  // Boost momentum 4-vector
154 
155  if (is_finalstate(track) || is_initialstate(track)) {
156 
157  const Vec Ekin = (track.status != TRK_ST_ININUCLEI ?
158  getKineticEnergy(track) * track.dir : Vec(0.0, 0.0, 0.0));
159 
160  const PxPyPzEVector p0(Ekin.x, Ekin.y, Ekin.z, track.E);
161  const PxPyPzEVector p1 = (*this)(p0);
162 
163  track.E = p1.E();
164  track.dir = Vec(p1.Px(), p1.Py(), p1.Pz());
165  if (track.dir.len() > 0) { track.dir.normalize(); }
166 
167  // Boost track vertex and compute new track length,
168  // assuming a relativistic particle traveling at the speed of light
169 
170  JVertex3D x0(getPosition(track.pos), track.t); // Vertex
171  JVertex3D x1(getPosition(track.pos + track.len * track.dir),
172  track.t + track.len / getSpeedOfLight()); // Vertex + travel distance
173 
174  (*this)(x0); // Boost vertices
175  (*this)(x1);
176 
177  track.t = x0.getT() / getSpeedOfLight();
178  track.pos = Vec(x0.getX(), x0.getY(), x0.getZ());
179 
180  track.len = x0.getDistance(x1);
181  }
182  }
183 
184 
185  /**
186  * Lorentz boost operator.
187  * N.B.: The time-over-threshold is not boosted.
188  *
189  * \param hit hit
190  */
191  void operator()(Hit& hit) const
192  {
193  using namespace std;
194  using namespace JPP;
195 
196  const double dt = hit.t - hit.tdc;
197 
198  JVertex3D vertex(getPosition(hit.pos), hit.t);
199  (*this)(vertex);
200 
201  hit.pos = Vec(vertex.getX(), vertex.getY(), vertex.getZ());
202  hit.t = vertex.getT();
203 
204  hit.tdc = hit.t - getGamma() * dt;
205  }
206 
207 
208  /**
209  * Lorentz boost operator.
210  *
211  * \param __begin begin of data
212  * \param __end end of data
213  */
214  template<class T>
215  void operator()(T __begin, T __end) const
216  {
217  using namespace std;
218 
219  for (T i = __begin; i != __end; ++i) {
220  (*this)(*i);
221  }
222  }
223 
224 
225  /**
226  * Lorentz boost operator.
227  * N.B.: Intermediate state tracks as well as\n
228  * the event timestamps and MC-time variable are not boosted.
229  *
230  * \param event event
231  */
232  void operator()(Evt& event) const
233  {
234  (*this)(event.mc_trks.begin(), event.mc_trks.end());
235  (*this)(event.mc_hits.begin(), event.mc_hits.end());
236  (*this)(event.trks.begin(), event.trks.end());
237  (*this)(event.hits.begin(), event.hits.end());
238  }
239  };
240 
241 
242  /**
243  * Get Lorentz boost to the Center of Mass (COM) frame for a given neutrino interaction.
244  *
245  * \param event event
246  */
248  {
249  using namespace JPP;
250 
251  if (has_neutrino(event)) {
252 
253  const Trk& nu = get_neutrino(event);
254 
255  const Vec beta = -nu.E / getE0(event) * nu.dir;
256 
257  return Boost(beta.x, beta.y, beta.z);
258 
259  } else {
260  THROW(JValueOutOfRange, "getBoostToCOM(): Given event does not correspond to a neutrino interaction.");
261  }
262  }
263 
264 
265  /**
266  * Boost event to the Center of Mass (COM) frame.
267  *
268  * \param event event
269  */
270  void boostToCOM(Evt& event)
271  {
272  const JLorentzBoost boost = getBoostToCOM(event);
273  boost(event);
274  }
275 }
276 
277 #endif
Exceptions.
Vec pos
hit position
Definition: Hit.hh:25
TPaveText * p1
void operator()(Hit &hit) const
Lorentz boost operator.
void boostToCOM(Evt &event)
Boost event to the Center of Mass (COM) frame.
double t
track time [ns] (when the particle is at pos )
Definition: Trk.hh:19
double z
Definition: Vec.hh:14
double getE0(const Evt &evt)
Get initial state energy of a neutrino interaction.
bool has_neutrino(const Evt &evt)
Test whether given event has an incoming neutrino.
Vec dir
track direction
Definition: Trk.hh:18
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:712
JLorentzBoost getBoostToCOM(const Evt &event)
Get Lorentz boost to the Center of Mass (COM) frame for a given neutrino interaction.
unsigned int tdc
hit tdc (=time in ns)
Definition: Hit.hh:16
double getDistance(const JVector3D &pos) const
Get distance to point.
Definition: JVector3D.hh:270
double E
Energy [GeV] (either MC truth or reconstructed)
Definition: Trk.hh:20
double y
Definition: Vec.hh:14
3D track with energy.
Definition: JTrack3E.hh:30
void operator()(JTrack3E &track, const double mass) const
Lorentz boost operator.
double len() const
Get length.
Definition: Vec.hh:145
JLorentzBoost()
Default constructor.
bool is_finalstate(const Trk &track)
Test whether given track corresponds to a final state particle.
LorentzVector< CoordSystem > operator()(const LorentzVector< CoordSystem > &x4) const
Lorentz boost operator.
double x
Definition: Vec.hh:14
The Vec class is a straightforward 3-d vector, which also works in pyroot.
Definition: Vec.hh:12
void operator()(Evt &event) const
Lorentz boost operator.
double getKineticEnergy(const double E, const double m)
Get kinetic energy of particle with given mass.
double getE() const
Get energy.
Definition: JTrack3E.hh:93
double getDY() const
Get y direction.
Definition: JVersor3D.hh:106
double getDX() const
Get x direction.
Definition: JVersor3D.hh:95
do set_variable OUTPUT_DIRECTORY $WORKDIR T
void setE(const double E)
Set energy.
Definition: JTrack3E.hh:104
JPosition3D getPosition(const Vec &pos)
Get position.
Definition of hit and track types and auxiliary methods for handling Monte Carlo data.
double len
length, if applicable [m]
Definition: Trk.hh:22
Physics constants.
JLorentzBoost(const Boost &boost)
Copy constructor.
static const int TRK_ST_ININUCLEI
Initial state nuclei (gseagen)
Definition: trkmembers.hh:19
double getY() const
Get y position.
Definition: JVector3D.hh:104
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
Definition: JTrack3D.hh:126
Auxiliary class for performing Lorentz boosts.
bool is_initialstate(const Trk &track)
Test whether given track corresponds to an initial state particle.
void setT(const double time)
Set time.
Vec & normalize()
Normalise this vector.
Definition: Vec.hh:159
int status
MC status code, see km3net-dataformat/definitions/trkmembers.csv for values.
Definition: Trk.hh:28
Vec pos
postion [m] of the track at time t
Definition: Trk.hh:17
void setDirection(const JDirection3D &dir)
Set direction.
Definition: Hit.hh:8
const double getSpeedOfLight()
Get speed of light.
double getGamma() const
Retrieve gamma factor.
void operator()(Trk &track) const
Lorentz boost operator.
double t
hit time (from tdc+calibration or MC truth)
Definition: Hit.hh:23
void operator()(T __begin, T __end) const
Lorentz boost operator.
double getX() const
Get x position.
Definition: JVector3D.hh:94
Data structure for position in three dimensions.
Definition: JPosition3D.hh:36
Exception for accessing a value in a collection that is outside of its range.
Definition: JException.hh:178
void operator()(JVertex3D &vertex) const
Lorentz boost operator.
const Trk & get_neutrino(const Evt &evt)
Get incoming neutrino.
double getT(const JVector3D &pos) const
Get arrival time of Cherenkov light at given position.
Definition: JVertex3D.hh:147
The Trk class represents a Monte Carlo (MC) particle as well as a reconstructed track/shower.
Definition: Trk.hh:14
double getZ() const
Get z position.
Definition: JVector3D.hh:115
double getDZ() const
Get z direction.
Definition: JVersor3D.hh:117
The Evt class respresent a Monte Carlo (MC) event as well as an offline event.
Definition: Evt.hh:20
void setPosition(const JVector3D &pos)
Set position.
Definition: JPosition3D.hh:152