Fork me on GitHub

source: svn/trunk/Utilities/Hector/include/H_BeamParticle.h@ 281

Last change on this file since 281 was 281, checked in by Xavier Rouby, 15 years ago

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