[5b822e5] | 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 |
|
---|
| 43 | using namespace std;
|
---|
| 44 |
|
---|
| 45 | /// Defines a particle from the beam and its transport through the beamline
|
---|
| 46 | class 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 | void getPath(const int x_or_y, const string& filename) const;
|
---|
| 135 | /// Computes the position of the particle at the end of each element of the beam.
|
---|
| 136 | void computePath(const H_AbstractBeamLine * beam, const bool NonLinear=true);
|
---|
| 137 | /// Computes the position of the particle at the end of each element of the beam.
|
---|
| 138 | void computePath(const H_AbstractBeamLine & beam, const bool NonLinear=true);
|
---|
| 139 | /// Clears H_BeamParticle::positions but keeps the initial vector.
|
---|
| 140 | void resetPath();
|
---|
| 141 |
|
---|
| 142 | private:
|
---|
| 143 | void init();
|
---|
| 144 | /// Particle mass [GeV]
|
---|
| 145 | double mp;
|
---|
| 146 | /// Particle charge [e]
|
---|
| 147 | double qp;
|
---|
| 148 | /// Longitudinal co-moving coordinate [m]
|
---|
| 149 | double fs;
|
---|
| 150 | /// Transverse (horizontal) coordinate [m]
|
---|
| 151 | double fx;
|
---|
| 152 | /// Transverse (vertical) coordinate [m]
|
---|
| 153 | double fy;
|
---|
| 154 | /// Direction of the 3-momentum in the horizontal plane [rad]
|
---|
| 155 | double thx;
|
---|
| 156 | /// Direction of the 3-momentum in the vertical plane [rad]
|
---|
| 157 | double thy;
|
---|
| 158 | /// Kinetic energy of the particle [GeV]
|
---|
| 159 | double energy;
|
---|
| 160 | /// 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.
|
---|
| 161 | bool hasstopped;
|
---|
| 162 | /// True if the particle has lost some (E,Q), i.e. if H_BeamParticle::emitGamma was used. Default = false.
|
---|
| 163 | bool hasemitted;
|
---|
| 164 | /// False if the particle has emitted a photon with impossible (E,Q). Default = true.
|
---|
| 165 | bool isphysical;
|
---|
| 166 | /// Vector (x,tx,y,ty,s) where the particle has stopped.
|
---|
| 167 | TVectorD * stop_position;
|
---|
| 168 | /// Element of the beamline (H_OpticalElement) where the particle has stopped.
|
---|
| 169 | H_OpticalElement * stop_element;
|
---|
| 170 | /// List of (x,tx,y,ty,s) vectors, after each optical element of the beam
|
---|
| 171 | vector<TVectorD> positions; // vector (x,tx,y,ty,s) after each optical element of the beam ([m],[rad],[m],[rad],[m])
|
---|
| 172 | /// Adds a new vector (x,tx,y,ty,s) at the end of H_BeamParticle::positions
|
---|
| 173 | void addPosition(const double , const double , const double , const double , const double );
|
---|
| 174 |
|
---|
| 175 | friend std::ostream& operator<< (std::ostream& os, const H_BeamParticle& p);
|
---|
| 176 | };
|
---|
| 177 | #endif
|
---|