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