[3] | 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 "TLorentzVector.h"
|
---|
| 29 |
|
---|
| 30 | // local #includes
|
---|
| 31 | #include "H_Parameters.h"
|
---|
| 32 | #include "H_AbstractBeamLine.h"
|
---|
| 33 | #include "H_OpticalElement.h"
|
---|
| 34 |
|
---|
| 35 | using namespace std;
|
---|
| 36 |
|
---|
| 37 | /// Defines a particle from the beam and its transport through the beamline
|
---|
| 38 | class H_BeamParticle {
|
---|
| 39 |
|
---|
| 40 | public:
|
---|
| 41 | void init();
|
---|
| 42 | /// Constructors and Destructor
|
---|
| 43 | //@{
|
---|
| 44 | H_BeamParticle();
|
---|
| 45 | H_BeamParticle(const H_BeamParticle&);
|
---|
| 46 | H_BeamParticle(const double, const double);
|
---|
| 47 | H_BeamParticle& operator=(const H_BeamParticle&);
|
---|
[216] | 48 | ~H_BeamParticle() {if(stop_position) delete stop_position; if(!stop_element) delete stop_element; positions.clear();}
|
---|
[3] | 49 | //@}
|
---|
| 50 | /// Smears the (x,y) coordinates of the particle [\f$ \mu m \f$]
|
---|
| 51 | void smearPos(const double, const double);
|
---|
| 52 | /// Smears the (x,y) coordinates of the particle (using default values in parameters.h)
|
---|
| 53 | void smearPos();
|
---|
| 54 | /// Smears the (x,y) angular coordinates of the particle [\f$ \mu rad \f$]
|
---|
| 55 | void smearAng(const double, const double);
|
---|
| 56 | /// Smears the (x,y) angular coordinates of the particle (using default values in parameters.h)
|
---|
| 57 | void smearAng();
|
---|
| 58 | /// Smears the Energy of the particle [GeV]
|
---|
| 59 | void smearE(const double);
|
---|
| 60 | /// Smears the Energy of the particle (using default values in parameters.h)
|
---|
| 61 | void smearE();
|
---|
| 62 | /// Smears the longitudinal position of the particle [\f$ \mu m \f$]
|
---|
| 63 | void smearS(const double);
|
---|
| 64 | /// Smears the longitudinal position of the particle (using default values in parameters.h)
|
---|
| 65 | void smearS();
|
---|
| 66 | /// Sets the energy [GeV].
|
---|
| 67 | void setE(const double);
|
---|
| 68 | /// Sets the particle 4-momentum \f$ P^\mu \f$
|
---|
| 69 | void set4Momentum(const double, const double, const double, const double);
|
---|
| 70 | /// Sets the particle 4-momentum \f$ P^\mu \f$
|
---|
| 71 | void set4Momentum(const TLorentzVector&);
|
---|
| 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 | double getTX() const {return thx;};
|
---|
| 86 | /// Returns the current \f$ \theta_y \f$ angular coordinate [\f$ \mu \f$rad]
|
---|
| 87 | 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
|
---|
[216] | 114 | void propagate(const H_AbstractBeamLine *, const string&);
|
---|
[3] | 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 | /// Computes the position of the particle at the end of each element of the beam, without non linear effects
|
---|
| 137 | void computePath(const H_AbstractBeamLine *);
|
---|
| 138 | /// Computes the position of the particle at the end of each element of the beam.
|
---|
| 139 | void computePath(const H_AbstractBeamLine *, const bool);
|
---|
| 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 | /// Clears H_BeamParticle::positions but keeps the initial vector.
|
---|
| 143 | void resetPath();
|
---|
| 144 |
|
---|
| 145 | private:
|
---|
| 146 | /// Particle mass [GeV]
|
---|
| 147 | double mp;
|
---|
| 148 | /// Particle charge [e]
|
---|
| 149 | double qp;
|
---|
| 150 | /// Longitudinal co-moving coordinate [m]
|
---|
| 151 | double fs;
|
---|
| 152 | /// Transverse (horizontal) coordinate [m]
|
---|
| 153 | double fx;
|
---|
| 154 | /// Transverse (vertical) coordinate [m]
|
---|
| 155 | double fy;
|
---|
| 156 | /// Direction of the 3-momentum in the horizontal plane [rad]
|
---|
| 157 | double thx;
|
---|
| 158 | /// Direction of the 3-momentum in the vertical plane [rad]
|
---|
| 159 | double thy;
|
---|
| 160 | /// Kinetic energy of the particle [GeV]
|
---|
| 161 | double energy;
|
---|
| 162 | /// 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.
|
---|
| 163 | bool hasstopped;
|
---|
| 164 | /// True if the particle has lost some (E,Q), i.e. if H_BeamParticle::emitGamma was used. Default = false.
|
---|
| 165 | bool hasemitted;
|
---|
| 166 | /// False if the particle has emitted a photon with impossible (E,Q). Default = true.
|
---|
| 167 | bool isphysical;
|
---|
| 168 | /// Vector (x,tx,y,ty,s) where the particle has stopped.
|
---|
| 169 | TVectorD * stop_position;
|
---|
| 170 | /// Element of the beamline (H_OpticalElement) where the particle has stopped.
|
---|
| 171 | H_OpticalElement * stop_element;
|
---|
| 172 | /// List of (x,tx,y,ty,s) vectors, after each optical element of the beam
|
---|
| 173 | vector<TVectorD> positions; // vector (x,tx,y,ty,s) after each optical element of the beam ([m],[rad],[m],[rad],[m])
|
---|
| 174 | /// Adds a new vector (x,tx,y,ty,s) at the end of H_BeamParticle::positions
|
---|
| 175 | void addPosition(const double , const double , const double , const double , const double );
|
---|
| 176 | };
|
---|
| 177 | #endif
|
---|