source: svn/trunk/external/Hector/H_BeamParticle.h@ 1379

#ifndef _H_BeamParticle_
#define _H_BeamParticle_

/*
---- Hector the simulator ----
   A fast simulator of particles through generic beamlines.
   J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be

        http://www.fynu.ucl.ac.be/hector.html

   Centre de Physique des Particules et de Phénoménologie (CP3)
   Université Catholique de Louvain (UCL)
*/

/// \file H_BeamParticle.h
/// \brief Class aiming at simulating a particle in the beamline

// from IP to RP, with emission of a photon of defined energy and Q.
// Units : angles [rad], distances [m], energies [GeV], masses [GeV], c=[1].
// !!! no comment statement at the end of a #define line !!!

// c++ #includes
#include <vector>

// ROOT #includes
#include "TMatrixD.h"
#include "TVectorD.h"
//#include "TPolyLine3D.h"
#include "TRandom.h"

// local #includes
#include "H_Parameters.h"
#include "H_AbstractBeamLine.h"
#include "H_OpticalElement.h"

using namespace std;

// local defines
#define LENGTH_VEC 5
#define INDEX_X 0
#define INDEX_TX 1
#define INDEX_Y 2
#define INDEX_TY 3
#define INDEX_S 4
// (x,theta_x,y,theta_y,s)

/// Defines a particle from the beam and its transport through the beamline
class H_BeamParticle {

        public:
                void init();
                /// Constructors and Destructor
                //@{
                H_BeamParticle();
                H_BeamParticle(const H_BeamParticle&);
                H_BeamParticle(const double, const double);
                H_BeamParticle& operator=(const H_BeamParticle&);
                ~H_BeamParticle() {delete stop_position; if(!stop_element) delete stop_element; positions.clear(); return; }
                //@}
                /// Smears the (x,y) coordinates of the particle [\f$ \mu m \f$]
                void smearPos(const double dx=SX,const double dy=SY, TRandom* r=gRandom);
                /// Smears the (x,y) angular coordinates of the particle [\f$ \mu rad \f$]
                void smearAng(const double tx=STX, const double ty=STY, TRandom* r=gRandom);
                /// Smears the Energy of the particle [GeV]
                void smearE(const double erre=SBE, TRandom* r=gRandom);
                /// Smears the longitudinal position of the particle [\f$ \mu m \f$]
                void smearS(const double errs=SS, TRandom* r=gRandom);
                /// Sets the energy [GeV].
                void setE(const double);
                /// Sets the particle 4-momentum \f$ P^\mu \f$
                void set4Momentum(const double, const double, const double, const double);      
                /// Clears H_BeamParticle::positions and sets the initial one.
                void setPosition(const double , const double , const double , const double , const double );
                /// Returns the particle mass [GeV]
                double getM() const {return mp;};
                /// Returns the particle charge [e]
                double getQ() const {return qp;};
                /// Returns the current x coordinate [\f$ \mu \f$m]
                double getX() const {return fx;};
                /// Returns the current y coordinate [\f$ \mu \f$m]
                double getY() const {return fy;};
                /// Returns the current s coordinate [m]
                inline double getS() const {return fs;};
                /// Returns the current \f$ \theta_x \f$ angular coordinate [\f$ \mu \f$rad]
                inline double getTX() const {return thx;};
                /// Returns the current \f$ \theta_y \f$ angular coordinate [\f$ \mu \f$rad]
                inline double getTY() const {return thy;};
                /// Returns the current particle energy [GeV]
                inline double getE() const {return energy;};
                /// Returns all the positions 
                vector<TVectorD> getPositions() const {return positions;};
                bool isPhysical() const {return isphysical;};
                /// \brief Simulates the emission of a photon in a random direction
                ///
                /// For \f$ p_{1} \rightarrow p_{2} \gamma \f$, kinematics imposes that
                /// \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>
                /// 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>
                /// As \f$ E^{2}_{\gamma} - (p_{1}-p_{2})^{2} \f$ could be numerically instable, we use here another form of this formula : <BR>
                /// \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$
                //@{
                void emitGamma(const double, const double, const double, const double);
                void emitGamma(const double, const double);
                //@}
                /// uses Pythia to generate some inelastic pp->pX collision as background
                void doInelastic();
                /// \brief Propagates the particule across the beamline until the s coordinate is reached
                ///
                /// Caution : "computePath" should be used before any "propagate" call <BR>
                /// Caution : "stopped" is not included in "propagate" : please run it afterward if needed 
                void propagate(const double ) ;
                /// Propagates the particle accross the beamline until a given element
                void propagate(const H_AbstractBeamLine *, const H_OpticalElement *);
                /// Propagates the particle accross the beamline until a given element
                void propagate(const H_AbstractBeamLine *, const string);
                /// Propagates the particle until the end of the beamline
                void propagate(const H_AbstractBeamLine *);
                /// Returns the phase vector of the particle
                const TMatrixD * getV() const;
                /// Returns the current phase vector of the particle (in H_BeamParticle::positions)
                const TVectorD * getPosition(const int ) const;
                /// Prints the properties of the particle
                void printProperties() const;
                /// Prints the phase vector of the particle
                void printV() const;
                /// Returns the element where the particle has been stopped
                const H_OpticalElement * getStoppingElement() const;
                /// Checks if the particle has been stopped in any element of the beamline
                bool stopped(const H_AbstractBeamLine *);
                /// Returns the StopPosition vector
                inline const TVectorD * getStopPosition() const { return stop_position; };
                // returns (-1,-1,-1,-1,-1) if not stopped (and then hasstopped is false)
                /// Shows all the vectors \f$ (x, \theta_x, y, \theta_y ,s) \f$ in H_BeamParticle::positions
                void showPositions() const; 
                /// Returns the particle path in the beamline
                ////TGraph * getPath(const int , const int ) const;
                /// Draws the particle path in the beamline in 3D
                ////TPolyLine3D * getPath3D(const H_AbstractBeamLine *, const bool, const int, const int) const;
                /// Computes the position of the particle at the end of each element of the beam, without non linear effects
                void computePath(const H_AbstractBeamLine *);
                /// Computes the position of the particle at the end of each element of the beam.
                void computePath(const H_AbstractBeamLine *, const bool);
                /// Computes the position of the particle at the end of each element of the beam.
                void computePath(const H_AbstractBeamLine &, const bool);
                /// Clears H_BeamParticle::positions but keeps the initial vector.
                void resetPath();

        private:
                /// Particle mass [GeV]
                double mp;
                /// Particle charge [e]
                double qp;
                /// Longitudinal co-moving coordinate [m]
                double fs;
                /// Transverse (horizontal) coordinate [m]
                double fx;
                /// Transverse (vertical) coordinate [m]
                double fy;
                /// Direction of the 3-momentum in the horizontal plane [rad]
                double thx;
                /// Direction of the 3-momentum in the vertical plane [rad]
                double thy;
                /// Kinetic energy of the particle [GeV]
                double energy;
                /// 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.
                bool hasstopped;
                /// True if the particle has lost some (E,Q), i.e. if H_BeamParticle::emitGamma was used. Default = false.
                bool hasemitted;
                /// False if the particle has emitted a photon with impossible (E,Q). Default = true.
                bool isphysical;
                /// Vector (x,tx,y,ty,s) where the particle has stopped.
                TVectorD * stop_position; 
                /// Element of the beamline (H_OpticalElement) where the particle has stopped. 
                H_OpticalElement * stop_element;
                /// List of (x,tx,y,ty,s) vectors, after each optical element of the beam
                vector<TVectorD> positions; // vector (x,tx,y,ty,s) after each optical element of the beam ([m],[rad],[m],[rad],[m])
                /// Adds a new vector (x,tx,y,ty,s) at the end of H_BeamParticle::positions
                void addPosition(const double , const double , const double , const double , const double );
};
#endif