Belle II Software  light-2309-munchkin
MCInitialParticles.h
1 /**************************************************************************
2  * basf2 (Belle II Analysis Software Framework) *
3  * Author: The Belle II Collaboration *
4  * *
5  * See git log for contributors and copyright holders. *
6  * This file is licensed under LGPL-3.0, see LICENSE.md. *
7  **************************************************************************/
8 
9 #pragma once
10 
11 #include <framework/gearbox/Const.h>
12 
13 #include <TMath.h>
14 #include <TVector3.h>
15 #include <Math/AxisAngle.h>
16 #include <Math/Boost.h>
17 #include <Math/LorentzRotation.h>
18 #include <Math/Vector3D.h>
19 #include <Math/Vector4D.h>
20 
21 namespace Belle2 {
32  class MCInitialParticles: public TObject {
33 
34  public:
35 
39  c_generateCMS = 1 << 0,
47  c_smearVertex = 1 << 3,
50  };
51 
53  MCInitialParticles(): TObject() {}
54 
58 
61  {
62  m_her = b.m_her; m_ler = b.m_ler; m_vertex = b.m_vertex;
63  m_validFlag = b.m_validFlag; m_generationFlags = b.m_generationFlags;
64  resetBoost();
65  return *this;
66  }
67 
69  bool operator==(const MCInitialParticles& b) const
70  {
71  // FIXME: sin(x) returns slightly different values on
72  // different platforms in some cases so we cannot just do an equality
73  // comparison. We need to do this more elegantly, this is just for
74  // testing if it solves all problems
75 #if defined(MCP_DBL_CMP) || defined(MCP_VEC3_CMP) || defined(MCP_VEC4_CMP)
76 #error Macro already defined, cannot continue
77 #endif
78 #define MCP_DBL_CMP(a,b,x) ((a.X()==b.X())||(std::abs(a.X()-b.X())<1e-10))
79 #define MCP_VEC3_CMP(a,b) (MCP_DBL_CMP(a,b,X) && MCP_DBL_CMP(a,b,Y) && MCP_DBL_CMP(a,b,Z))
80 #define MCP_VEC4_CMP(a,b) (MCP_VEC3_CMP(a,b) && MCP_DBL_CMP(a,b,E))
81  return MCP_VEC4_CMP(m_her, b.m_her) && MCP_VEC4_CMP(m_ler, b.m_ler) && MCP_VEC3_CMP(m_vertex, b.m_vertex)
82  && m_validFlag == b.m_validFlag &&
83  m_generationFlags == b.m_generationFlags;
84 #undef MCP_DBL_CMP
85 #undef MCP_VEC3_CMP
86 #undef MCP_VEC4_CMP
87  }
88 
91 
98  void set(const ROOT::Math::PxPyPzEVector& her, const ROOT::Math::PxPyPzEVector& ler, const ROOT::Math::XYZVector& vertex)
99  {
100  m_her = her;
101  m_ler = ler;
102  m_vertex = vertex;
103  m_validFlag = true;
104  resetBoost();
105  }
106 
117  void setByLorentzTransformation(double Ecms, double bX, double bY, double bZ, double angleXZ, double angleYZ,
118  const ROOT::Math::XYZVector& vertex)
119  {
120  if (m_labToCMS) delete m_labToCMS;
121  if (m_CMSToLab) delete m_CMSToLab;
122 
123  m_invariantMass = Ecms;
124  m_CMSToLab = new ROOT::Math::LorentzRotation();
125  m_labToCMS = new ROOT::Math::LorentzRotation();
126  *m_CMSToLab = cmsToLab(bX, bY, bZ, angleXZ, angleYZ);
127  *m_labToCMS = m_CMSToLab->Inverse();
128 
129  const double me = Const::electron.getMass();
130  double p = sqrt(Ecms * Ecms / 4 - me * me);
131  m_her = (*m_CMSToLab) * ROOT::Math::PxPyPzEVector(0.0, 0.0, p, Ecms / 2);
132  m_ler = (*m_CMSToLab) * ROOT::Math::PxPyPzEVector(0.0, 0.0, -p, Ecms / 2);
133 
134  m_vertex = vertex;
135  m_validFlag = true;
136  }
137 
139  void setHER(const ROOT::Math::PxPyPzEVector& her)
140  {
141  m_her = her;
142  resetBoost();
143  }
144 
146  void setLER(const ROOT::Math::PxPyPzEVector& ler)
147  {
148  m_ler = ler;
149  resetBoost();
150  }
151 
153  void setVertex(const ROOT::Math::XYZVector& vertex)
154  {
155  m_vertex = vertex;
156  }
157 
159  void setTime(double time) {m_time = time;}
160 
162  virtual void setGenerationFlags(int flags) { m_generationFlags = flags; }
163 
165  const ROOT::Math::PxPyPzEVector& getHER() const { return m_her; }
166 
168  const ROOT::Math::PxPyPzEVector& getLER() const { return m_ler; }
169 
171  const ROOT::Math::XYZVector& getVertex() const { return m_vertex; }
172 
174  double getTime() const {return m_time;}
175 
177  double getEnergy() const { return (m_her + m_ler).E(); }
178 
180  double getMass() const { calculateBoost(); return m_invariantMass; }
181 
183  const ROOT::Math::LorentzRotation& getLabToCMS() const
184  {
185  calculateBoost(); return *m_labToCMS;
186  }
187 
189  const ROOT::Math::LorentzRotation& getCMSToLab() const
190  {
191  calculateBoost(); return *m_CMSToLab;
192  }
193 
195  bool getValidFlag() const { return m_validFlag; }
196 
198  int getGenerationFlags() const { return m_generationFlags; }
199 
201  bool hasGenerationFlags(int flags) const { return (m_generationFlags & flags) == flags; }
202 
205  std::string getGenerationFlagString(const std::string& separator = " ") const;
206 
207 
215  static ROOT::Math::LorentzRotation cmsToLab(double bX, double bY, double bZ, double angleXZ, double angleYZ);
216 
217  private:
218 
220  void calculateBoost() const;
222  void resetBoost();
224  ROOT::Math::PxPyPzEVector m_her;
226  ROOT::Math::PxPyPzEVector m_ler;
228  ROOT::Math::XYZVector m_vertex;
230  double m_time = 0;
232  mutable ROOT::Math::LorentzRotation* m_labToCMS{nullptr};
234  mutable ROOT::Math::LorentzRotation* m_CMSToLab{nullptr};
236  mutable double m_invariantMass{0.0};
238  bool m_validFlag = false;
239 
240  protected:
241 
244 
245  private:
246 
249  };
250 
252  {
253  if (m_labToCMS)
254  return;
255 
256  ROOT::Math::PxPyPzEVector beam = m_her + m_ler;
257  // Save the invariant mass because it's used very often in analysis
258  m_invariantMass = beam.M();
259 
260  // If we generate events in CMS we already are in CMS and there is no
261  // transformation so let's use the identity
263  m_labToCMS = new ROOT::Math::LorentzRotation();
264  m_CMSToLab = new ROOT::Math::LorentzRotation();
265  return;
266  }
267 
268  // Transformation from Lab system to CMS system
269  m_labToCMS = new ROOT::Math::LorentzRotation(ROOT::Math::Boost(beam.BoostToCM()));
270  // boost HER e- from Lab system to CMS system
271  const ROOT::Math::PxPyPzEVector electronCMS = (*m_labToCMS) * m_her;
272  // now rotate CMS such that incoming e- is parallel to z-axis
273  const ROOT::Math::XYZVector zaxis(0., 0., 1.);
274  ROOT::Math::XYZVector rotaxis = zaxis.Cross(electronCMS.Vect()) / electronCMS.P();
275  double rotangle = TMath::ASin(rotaxis.R());
276  const ROOT::Math::LorentzRotation rotation(ROOT::Math::AxisAngle(rotaxis, -rotangle));
277  *m_labToCMS = rotation * (*m_labToCMS);
278 
279  //cache derived quantities
280  m_CMSToLab = new ROOT::Math::LorentzRotation(m_labToCMS->Inverse());
281  }
282 
283 
284 
285 
286 
288  {
289  delete m_labToCMS;
290  delete m_CMSToLab;
291  m_labToCMS = nullptr;
292  m_CMSToLab = nullptr;
293  m_invariantMass = 0.0;
294  }
295 
297 } //Belle2 namespace
double getMass() const
Particle mass.
Definition: UnitConst.cc:356
static const ChargedStable electron
electron particle
Definition: Const.h:650
This class contains the initial state for the given event.
const ROOT::Math::PxPyPzEVector & getLER() const
Get 4vector of the low energy beam.
int m_generationFlags
Flags to be used when generating events.
ROOT::Math::LorentzRotation * m_CMSToLab
transient
bool operator==(const MCInitialParticles &b) const
Equality operator.
int getGenerationFlags() const
Get the generation flags to be used for event generation (ORed combination of EGenerationFlags)
bool getValidFlag() const
Get the flag to check if a valid MCInitialParticles object was already generated and filled in an eve...
void setVertex(const ROOT::Math::XYZVector &vertex)
Set the vertex position.
const ROOT::Math::PxPyPzEVector & getHER() const
Get 4vector of the high energy beam.
double m_time
collision time
ROOT::Math::PxPyPzEVector m_ler
LER 4vector.
void setTime(double time)
Set collison time.
double getEnergy() const
Get the the actual collision energy (in lab system)
ROOT::Math::PxPyPzEVector m_her
HER 4vector.
EGenerationFlags
Possible Flags for initial event generation.
@ c_smearBeam
smear the full beam momentum (energy and direction)
@ c_generateCMS
generate initial event in CMS instead of lab
@ c_smearBeamEnergy
smear energy of HER and LER (but not direction)
@ c_smearBeamDirection
smear direction of HER and LER (but not energy)
ROOT::Math::LorentzRotation * m_labToCMS
Boost from Lab into CMS.
const ROOT::Math::XYZVector & getVertex() const
Get the position of the collision.
MCInitialParticles & operator=(const MCInitialParticles &b)
Assignment operator.
ROOT::Math::XYZVector m_vertex
collision position
std::string getGenerationFlagString(const std::string &separator=" ") const
Return string representation of all active flags for printing.
virtual void setGenerationFlags(int flags)
Set the generation flags to be used for event generation (ORed combination of EGenerationFlags)
static ROOT::Math::LorentzRotation cmsToLab(double bX, double bY, double bZ, double angleXZ, double angleYZ)
Return the LorentzRotation from CMS to LAB based on the following parameters.
MCInitialParticles()
Default constructor.
const ROOT::Math::LorentzRotation & getLabToCMS() const
Return the LorentzRotation to convert from lab to CMS frame.
void setByLorentzTransformation(double Ecms, double bX, double bY, double bZ, double angleXZ, double angleYZ, const ROOT::Math::XYZVector &vertex)
Initialize the event values from CMS energy and parameters of the Lorentz transformation between LAB ...
bool hasGenerationFlags(int flags) const
Check if a certain set of EGenerationFlags is set.
void setHER(const ROOT::Math::PxPyPzEVector &her)
Set the High Energy Beam 4-momentum.
MCInitialParticles(const MCInitialParticles &b)
Copy constructor.
void setLER(const ROOT::Math::PxPyPzEVector &ler)
Set the Low Energy Beam 4-momentum.
virtual ~MCInitialParticles()
Free memory of the LorentzRotation if it was created.
const ROOT::Math::LorentzRotation & getCMSToLab() const
Return the LorentzRotation to convert from CMS to lab frame.
double getMass() const
Get the invariant mass of the collision (= energy in CMS)
double getTime() const
Get collison time.
void set(const ROOT::Math::PxPyPzEVector &her, const ROOT::Math::PxPyPzEVector &ler, const ROOT::Math::XYZVector &vertex)
Set the initial event values, i.e.
ClassDef(MCInitialParticles, 5)
ROOT Dictionary.
void resetBoost()
Reset cached transformations after changing parameters.
void calculateBoost() const
Calculate the boost if necessary.
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24