Fork me on GitHub

source: git/external/Hector/H_BeamParticle.h@ f98efbf

ImprovedOutputFile Timing dual_readout llp
Last change on this file since f98efbf was 3c40083, checked in by pavel <pavel@…>, 11 years ago

switch to a more stable Hector version
  • Property mode set to 100644
File size: 8.1 KB
Line 
1#ifndef _H_BeamParticle_
2#define _H_BeamParticle_
3
4/*
5---- Hector the simulator ----
6 A fast simulator of particles through generic beamlines.
7 J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be
8
9 http://www.fynu.ucl.ac.be/hector.html
10
11 Centre de Physique des Particules et de Phénoménologie (CP3)
12 Université Catholique de Louvain (UCL)
13*/
14
15/// \file H_BeamParticle.h
16/// \brief Class aiming at simulating a particle in the beamline
17
18// from IP to RP, with emission of a photon of defined energy and Q.
19// Units : angles [rad], distances [m], energies [GeV], masses [GeV], c=[1].
20// !!! no comment statement at the end of a #define line !!!
21
22// c++ #includes
23#include <vector>
24
25// ROOT #includes
26#include "TMatrixD.h"
27#include "TVectorD.h"
28//#include "TPolyLine3D.h"
29#include "TRandom.h"
30
31// local #includes
32#include "H_Parameters.h"
33#include "H_AbstractBeamLine.h"
34#include "H_OpticalElement.h"
35
36using namespace std;
37
38// local defines
39#define LENGTH_VEC 5
40#define INDEX_X 0
41#define INDEX_TX 1
42#define INDEX_Y 2
43#define INDEX_TY 3
44#define INDEX_S 4
45// (x,theta_x,y,theta_y,s)
46
47/// Defines a particle from the beam and its transport through the beamline
48class H_BeamParticle {
49
50 public:
51 void init();
52 /// Constructors and Destructor
53 //@{
54 H_BeamParticle();
55 H_BeamParticle(const H_BeamParticle&);
56 H_BeamParticle(const double, const double);
57 H_BeamParticle& operator=(const H_BeamParticle&);
58 ~H_BeamParticle() {delete stop_position; if(!stop_element) delete stop_element; positions.clear(); return; }
59 //@}
60 /// Smears the (x,y) coordinates of the particle [\f$ \mu m \f$]
61 void smearPos(const double dx=SX,const double dy=SY, TRandom* r=gRandom);
62 /// Smears the (x,y) angular coordinates of the particle [\f$ \mu rad \f$]
63 void smearAng(const double tx=STX, const double ty=STY, TRandom* r=gRandom);
64 /// Smears the Energy of the particle [GeV]
65 void smearE(const double erre=SBE, TRandom* r=gRandom);
66 /// Smears the longitudinal position of the particle [\f$ \mu m \f$]
67 void smearS(const double errs=SS, TRandom* r=gRandom);
68 /// Sets the energy [GeV].
69 void setE(const double);
70 /// Sets the particle 4-momentum \f$ P^\mu \f$
71 void set4Momentum(const double, const double, const double, const double);
72 /// Clears H_BeamParticle::positions and sets the initial one.
73 void setPosition(const double , const double , const double , const double , const double );
74 /// Returns the particle mass [GeV]
75 double getM() const {return mp;};
76 /// Returns the particle charge [e]
77 double getQ() const {return qp;};
78 /// Returns the current x coordinate [\f$ \mu \f$m]
79 double getX() const {return fx;};
80 /// Returns the current y coordinate [\f$ \mu \f$m]
81 double getY() const {return fy;};
82 /// Returns the current s coordinate [m]
83 inline double getS() const {return fs;};
84 /// Returns the current \f$ \theta_x \f$ angular coordinate [\f$ \mu \f$rad]
85 inline double getTX() const {return thx;};
86 /// Returns the current \f$ \theta_y \f$ angular coordinate [\f$ \mu \f$rad]
87 inline double getTY() const {return thy;};
88 /// Returns the current particle energy [GeV]
89 inline double getE() const {return energy;};
90 /// Returns all the positions
91 vector<TVectorD> getPositions() const {return positions;};
92 bool isPhysical() const {return isphysical;};
93 /// \brief Simulates the emission of a photon in a random direction
94 ///
95 /// For \f$ p_{1} \rightarrow p_{2} \gamma \f$, kinematics imposes that
96 /// \f$ Q^{2} = E^{2}_{\gamma} -p^{2}_{1} -p^{2}_{2} + 2p_{1}p_{2} cos(\theta) \f$ where \f$ \theta \f$ is the particle scattering angle and \f$ p_{i} = \|\vec{p_{i}}\| \f$. <BR>
97 /// So, \f$ Q^{2}_{min} = E^{2}_{\gamma} - (p_{1}+p_{2})^{2} \f$ and \f$ Q^{2}_{max} = E^{2}_{\gamma} - (p_{1}-p_{2})^{2} \f$.<BR>
98 /// As \f$ E^{2}_{\gamma} - (p_{1}-p_{2})^{2} \f$ could be numerically instable, we use here another form of this formula : <BR>
99 /// \f$ Q^{2}_{max} = -2 * \big( \frac{M_{p} E_{\gamma}}{p_{1}+p_{2}} \big) \big[ 1 + \frac{E^{2}_{1} + E^{2}_{2} - M^{2}_{p} }{ E_{1} E_{2} + p_{1} p_{2}} \big] \f$
100 //@{
101 void emitGamma(const double, const double, const double, const double);
102 void emitGamma(const double, const double);
103 //@}
104 /// uses Pythia to generate some inelastic pp->pX collision as background
105 void doInelastic();
106 /// \brief Propagates the particule across the beamline until the s coordinate is reached
107 ///
108 /// Caution : "computePath" should be used before any "propagate" call <BR>
109 /// Caution : "stopped" is not included in "propagate" : please run it afterward if needed
110 void propagate(const double ) ;
111 /// Propagates the particle accross the beamline until a given element
112 void propagate(const H_AbstractBeamLine *, const H_OpticalElement *);
113 /// Propagates the particle accross the beamline until a given element
114 void propagate(const H_AbstractBeamLine *, const string);
115 /// Propagates the particle until the end of the beamline
116 void propagate(const H_AbstractBeamLine *);
117 /// Returns the phase vector of the particle
118 const TMatrixD * getV() const;
119 /// Returns the current phase vector of the particle (in H_BeamParticle::positions)
120 const TVectorD * getPosition(const int ) const;
121 /// Prints the properties of the particle
122 void printProperties() const;
123 /// Prints the phase vector of the particle
124 void printV() const;
125 /// Returns the element where the particle has been stopped
126 const H_OpticalElement * getStoppingElement() const;
127 /// Checks if the particle has been stopped in any element of the beamline
128 bool stopped(const H_AbstractBeamLine *);
129 /// Returns the StopPosition vector
130 inline const TVectorD * getStopPosition() const { return stop_position; };
131 // returns (-1,-1,-1,-1,-1) if not stopped (and then hasstopped is false)
132 /// Shows all the vectors \f$ (x, \theta_x, y, \theta_y ,s) \f$ in H_BeamParticle::positions
133 void showPositions() const;
134 /// Returns the particle path in the beamline
135 ////TGraph * getPath(const int , const int ) const;
136 /// Draws the particle path in the beamline in 3D
137 ////TPolyLine3D * getPath3D(const H_AbstractBeamLine *, const bool, const int, const int) const;
138 /// Computes the position of the particle at the end of each element of the beam, without non linear effects
139 void computePath(const H_AbstractBeamLine *);
140 /// Computes the position of the particle at the end of each element of the beam.
141 void computePath(const H_AbstractBeamLine *, const bool);
142 /// Computes the position of the particle at the end of each element of the beam.
143 void computePath(const H_AbstractBeamLine &, const bool);
144 /// Clears H_BeamParticle::positions but keeps the initial vector.
145 void resetPath();
146
147 private:
148 /// Particle mass [GeV]
149 double mp;
150 /// Particle charge [e]
151 double qp;
152 /// Longitudinal co-moving coordinate [m]
153 double fs;
154 /// Transverse (horizontal) coordinate [m]
155 double fx;
156 /// Transverse (vertical) coordinate [m]
157 double fy;
158 /// Direction of the 3-momentum in the horizontal plane [rad]
159 double thx;
160 /// Direction of the 3-momentum in the vertical plane [rad]
161 double thy;
162 /// Kinetic energy of the particle [GeV]
163 double energy;
164 /// True if the particle has stopped (i.e. : if the particle transverse position has been out of any optics element aperture). <BR> See H_BeamParticle::stopped method. <BR>Default = false.
165 bool hasstopped;
166 /// True if the particle has lost some (E,Q), i.e. if H_BeamParticle::emitGamma was used. Default = false.
167 bool hasemitted;
168 /// False if the particle has emitted a photon with impossible (E,Q). Default = true.
169 bool isphysical;
170 /// Vector (x,tx,y,ty,s) where the particle has stopped.
171 TVectorD * stop_position;
172 /// Element of the beamline (H_OpticalElement) where the particle has stopped.
173 H_OpticalElement * stop_element;
174 /// List of (x,tx,y,ty,s) vectors, after each optical element of the beam
175 vector<TVectorD> positions; // vector (x,tx,y,ty,s) after each optical element of the beam ([m],[rad],[m],[rad],[m])
176 /// Adds a new vector (x,tx,y,ty,s) at the end of H_BeamParticle::positions
177 void addPosition(const double , const double , const double , const double , const double );
178};
179#endif
Note: See TracBrowser for help on using the repository browser.