Fork me on GitHub

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

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

update to Hector_1_5_2
File size: 8.7 KB
RevLine 
[3]1#ifndef _H_BeamParticle_
2#define _H_BeamParticle_
3
[281]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 * * * * * * * * * * * * * * * * * * * * * * * * * * * */
[3]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"
[281]36#include "TRandom3.h"
[3]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&);
[216]55 ~H_BeamParticle() {if(stop_position) delete stop_position; if(!stop_element) delete stop_element; positions.clear();}
[3]56 //@}
57 /// Smears the (x,y) coordinates of the particle [\f$ \mu m \f$]
[281]58 void smearPos(const double dx=SX,const double dy=SY, TRandom* r=gRandom);
[3]59 /// Smears the (x,y) angular coordinates of the particle [\f$ \mu rad \f$]
[281]60 void smearAng(const double tx=STX, const double ty=STY, TRandom* r=gRandom);
[3]61 /// Smears the Energy of the particle [GeV]
[281]62 void smearE(const double erre=SBE, TRandom* r=gRandom);
[3]63 /// Smears the longitudinal position of the particle [\f$ \mu m \f$]
[281]64 void smearS(const double errs=SS, TRandom* r=gRandom);
[3]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]
[281]74 const double getM() const {return mp;};
[3]75 /// Returns the particle charge [e]
[281]76 const double getQ() const {return qp;};
[3]77 /// Returns the current x coordinate [\f$ \mu \f$m]
[281]78 const double getX() const {return fx;};
[3]79 /// Returns the current y coordinate [\f$ \mu \f$m]
[281]80 const double getY() const {return fy;};
[3]81 /// Returns the current s coordinate [m]
[281]82 const double getS() const {return fs;};
[3]83 /// Returns the current \f$ \theta_x \f$ angular coordinate [\f$ \mu \f$rad]
[281]84 const double getTX() const {return thx;};
[3]85 /// Returns the current \f$ \theta_y \f$ angular coordinate [\f$ \mu \f$rad]
[281]86 const double getTY() const {return thy;};
[3]87 /// Returns the current particle energy [GeV]
[281]88 const double getE() const {return energy;};
[3]89 /// Returns all the positions
90 vector<TVectorD> getPositions() const {return positions;};
[281]91 const bool isPhysical() const {return isphysical;};
[3]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 //@{
[281]100 void emitGamma(const double gee, const double gq2, const double phimin=0, const double phimax=2*TMath::Pi());
[3]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
[216]112 void propagate(const H_AbstractBeamLine *, const string&);
[3]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
[281]120 void printProperties() const {cout << *this; return;};
[3]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
[281]126 const bool stopped(const H_AbstractBeamLine *);
[3]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;
[400]134 void getPath(const int x_or_y, const string& filename) const;
[3]135 /// Computes the position of the particle at the end of each element of the beam.
[281]136 void computePath(const H_AbstractBeamLine * beam, const bool NonLinear=true);
[3]137 /// Computes the position of the particle at the end of each element of the beam.
[281]138 void computePath(const H_AbstractBeamLine & beam, const bool NonLinear=true);
[3]139 /// Clears H_BeamParticle::positions but keeps the initial vector.
140 void resetPath();
141
142 private:
[281]143 void init();
[3]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 );
[237]174
[281]175 friend std::ostream& operator<< (std::ostream& os, const H_BeamParticle& p);
[3]176};
177#endif
Note: See TracBrowser for help on using the repository browser.