BaseParticlePropagator/interface/BaseParticlePropagator.h

0001 /*Emacs: -*- C++ -*- */
0002 #ifndef CommonTools_BaseParticlePropagator_BaseParticlePropagator_h
0003 #define CommonTools_BaseParticlePropagator_BaseParticlePropagator_h
0004
0005 /**
0006 * This class is aimed at propagating charged and neutral
0007 * particles (yet under the form of a RawParticle) from
0008 * a given vertex to a cylinder, defined by a radius rCyl
0009 * and a length 2*zCyl, centered in (0,0,0) and whose axis
0010 * is parallel to the B field (B is oriented along z, by
0011 * definition of the z axis).
0012 *
0013 * The main method
0014 *
0015 *     bool propagate()
0016 *
0017 * returns true if an intersection between the propagated
0018 * RawParticle and the cylinder is found. The location of
0019 * the intersection is given by the value of success:
0020 *
0021 *   * success = 2 : an intersection is found forward in
0022 *                   the cylinder endcaps. The RawParticle
0023 *                   vertex and momentum are overwritten by
0024 *                   the values at this intersection.
0025 *
0026 *   * success = 1 : an intersection is found forward in
0027 *                   the cylinder barrel. The RawParticle
0028 *                   vertex and momentum are overwritten by
0029 *                   the values at this intersection.
0030 *
0031 *   * success =-2 : an intersection is found backward in
0032 *                   the cylinder endcaps. The RawParticle
0033 *                   vertex and momentum are NOT overwritten
0034 *                   by the values at this intersection.
0035 *
0036 *   * success =-1 : an intersection is found backward in
0037 *                   the cylinder barrel. The RawParticle
0038 *                   vertex and momentum are NOT overwritten
0039 *                   by the values at this intersection.
0040 *
0041 *   * success = 0 : No intersection is found after half a
0042 *                   loop (only possible if firstLoop=true).
0043 *                   The vertex and momentum are overwritten
0044 *                   by the values after this half loop.
0045 *
0046 *
0047 * The method
0048 *
0049 *     propagated()
0050 *
0051 * returns a new RawParticle with the propagated coordinates,
0052 * if overwriting is not considered an advantage by the user.
0053 *
0054 * Particles can be propagated backward with the method
0055 *
0056 *     backPropagate()
0057 *
0058 * Member functions
0059 *
0060 *  o propagateToPreshowerLayer1(),
0061 *  o propagateToPreshowerLayer2(),
0062 *  o propagateToEcalEntrance(),
0063 *  o propagateToHcalEntrance(),
0064 *  o propagateToHcalExit(),
0065 *  o propagateToClosestApproach(),
0066 *
0067 * only need a momentum, a vertex and an electric charge
0068 * to operate. Radii, half-lengths and default B field (4T)
0069 * are set therein by default.
0070 *
0071 * As of today, no average loss of energy (dE/dx, Brems) is
0072 * considered in the propagation. No uncertainty (multiple
0073 * scattering, dE/dx, Brems) is yet implemented.
0074 *
0075 * \author Patrick Janot
0076 * $Date : 19-Aug-2002, with subsequent modification for FAMOS
0077 * \version 15-Dec-2003 */
0078
0079 //FAMOS
0080 #include "CommonTools/BaseParticlePropagator/interface/RawParticle.h"
0081
0082 class BaseParticlePropagator {
0083 public:
0084   /// Default c'tor
0085   BaseParticlePropagator();
0086
0087   /** Constructors taking as arguments a RawParticle, as well as the radius,
0088       half-height and magnetic field defining the cylinder for which
0089       propagation is to be performed, and optionally, the proper decay time */
0090   BaseParticlePropagator(const RawParticle& myPart, double r, double z, double B);
0091   BaseParticlePropagator(const RawParticle& myPart, double r, double z, double B, double t);
0092
0093   /// Initialize internal switches and quantities
0094   void init();
0095
0096   /** Update the current instance, after the propagation of the particle
0097       to the surface of the cylinder */
0098   bool propagate();
0099
0100   /** Update the current instance, after the back-propagation of the particle
0101       to the surface of the cylinder */
0102   bool backPropagate();
0103
0104   /** Return a new instance, corresponding to the particle propagated
0105       to the surface of the cylinder */
0106   BaseParticlePropagator propagated() const;
0107
0108   /** Update the particle after propagation to the closest approach from
0109       Z axis, to the preshower layer 1 & 2, to the ECAL entrance, to the
0110       HCAL entrance, the HCAL 2nd and 3rd layer (not coded yet), the VFCAL
0111       entrance, or any BoundSurface(disk or cylinder)*/
0112
0113   bool propagateToClosestApproach(double x0 = 0., double y0 = 0, bool first = true);
0114   bool propagateToEcal(bool first = true);
0115   bool propagateToPreshowerLayer1(bool first = true);
0116   bool propagateToPreshowerLayer2(bool first = true);
0117   bool propagateToEcalEntrance(bool first = true);
0118   bool propagateToHcalEntrance(bool first = true);
0119   bool propagateToVFcalEntrance(bool first = true);
0120   bool propagateToHcalExit(bool first = true);
0121   bool propagateToHOLayer(bool first = true);
0122   bool propagateToNominalVertex(const XYZTLorentzVector& hit2 = XYZTLorentzVector(0., 0., 0., 0.));
0123   bool propagateToBeamCylinder(const XYZTLorentzVector& v, double radius = 0.);
0124   /// Set the propagation characteristics (rCyl, zCyl and first loop only)
0125   void setPropagationConditions(double r, double z, bool firstLoop = true);
0126
0127 private:
0128   RawParticle particle_;
0129   /// Simulated particle that is to be resp has been propagated
0130   //  RawParticle   particle;
0131   /// Radius of the cylinder (centred at 0,0,0) to which propagation is done
0132   double rCyl;
0133   double rCyl2;
0134   /// Half-height of the cylinder (centred at 0,0,0) to which propagation is done
0135   double zCyl;
0136   /// Magnetic field in the cylinder, oriented along the Z axis
0137   double bField;
0138   /// The proper decay time of the particle
0139   double properDecayTime;
0140   /// The debug level
0141   bool debug;
0142
0143 protected:
0144   /// 0:propagation still be done, 1:reached 'barrel', 2:reached 'endcaps'
0145   int success;
0146   /// The particle traverses some real material
0147   bool fiducial;
0148
0149   /// The speed of light in mm/ns (!) without clhep (yeaaahhh!)
0150   inline double c_light() const { return 299.792458; }
0151
0152 private:
0153   /// Do only the first half-loop
0154   bool firstLoop;
0155   /// The particle decayed while propagated !
0156   bool decayed;
0157   /// The proper time of the particle
0158   double properTime;
0159   /// The propagation direction
0160   int propDir;
0161
0162 public:
0163   /// The particle being propagated
0164   inline RawParticle const& particle() const { return particle_; }
0165   inline RawParticle& particle() { return particle_; }
0166   void setParticle(RawParticle const& iParticle) { particle_ = iParticle; }
0167
0168   /// Set the proper decay time
0169   inline void setProperDecayTime(double t) { properDecayTime = t; }
0170
0171   /// Just an internal trick
0172   inline void increaseRCyl(double delta) {
0173     rCyl = rCyl + delta;
0174     rCyl2 = rCyl * rCyl;
0175   }
0176
0177   /// Transverse impact parameter
0178   double xyImpactParameter(double x0 = 0., double y0 = 0.) const;
0179
0180   /// Longitudinal impact parameter
0181   inline double zImpactParameter(double x0 = 0, double y0 = 0.) const {
0182     // Longitudinal impact parameter
0183     return particle_.Z() - particle_.Pz() * std::sqrt(((particle_.X() - x0) * (particle_.X() - x0) +
0184                                                        (particle_.Y() - y0) * (particle_.Y() - y0)) /
0185                                                       particle_.Perp2());
0186   }
0187
0188   /// The helix Radius
0189   inline double helixRadius() const {
0190     // The helix Radius
0191     //
0192     // The helix' Radius sign accounts for the orientation of the magnetic field
0193     // (+ = along z axis) and the sign of the electric charge of the particle.
0194     // It signs the rotation of the (charged) particle around the z axis:
0195     // Positive means anti-clockwise, negative means clockwise rotation.
0196     //
0197     // The radius is returned in cm to match the units in RawParticle.
0198     return particle_.charge() == 0 ? 0.0 : -particle_.Pt() / (c_light() * 1e-5 * bField * particle_.charge());
0199   }
0200
0201   inline double helixRadius(double pT) const {
0202     // a faster version of helixRadius, once Perp() has been computed
0203     return particle_.charge() == 0 ? 0.0 : -pT / (c_light() * 1e-5 * bField * particle_.charge());
0204   }
0205
0206   /// The azimuth of the momentum at the vertex
0207   inline double helixStartPhi() const {
0208     // The azimuth of the momentum at the vertex
0209     return particle_.Px() == 0.0 && particle_.Py() == 0.0 ? 0.0 : std::atan2(particle_.Py(), particle_.Px());
0210   }
0211
0212   /// The x coordinate of the helix axis
0213   inline double helixCentreX() const {
0214     // The x coordinate of the helix axis
0215     return particle_.X() - helixRadius() * std::sin(helixStartPhi());
0216   }
0217
0218   inline double helixCentreX(double radius, double phi) const {
0219     // Fast version of helixCentreX()
0220     return particle_.X() - radius * std::sin(phi);
0221   }
0222
0223   /// The y coordinate of the helix axis
0224   inline double helixCentreY() const {
0225     // The y coordinate of the helix axis
0226     return particle_.Y() + helixRadius() * std::cos(helixStartPhi());
0227   }
0228
0229   inline double helixCentreY(double radius, double phi) const {
0230     // Fast version of helixCentreX()
0231     return particle_.Y() + radius * std::cos(phi);
0232   }
0233
0234   /// The distance between the cylinder and the helix axes
0235   inline double helixCentreDistToAxis() const {
0236     // The distance between the cylinder and the helix axes
0237     double xC = helixCentreX();
0238     double yC = helixCentreY();
0239     return std::sqrt(xC * xC + yC * yC);
0240   }
0241
0242   inline double helixCentreDistToAxis(double xC, double yC) const {
0243     // Faster version of helixCentreDistToAxis
0244     return std::sqrt(xC * xC + yC * yC);
0245   }
0246
0247   /// The azimuth if the vector joining the cylinder and the helix axes
0248   inline double helixCentrePhi() const {
0249     // The azimuth if the vector joining the cylinder and the helix axes
0250     double xC = helixCentreX();
0251     double yC = helixCentreY();
0252     return xC == 0.0 && yC == 0.0 ? 0.0 : std::atan2(yC, xC);
0253   }
0254
0255   inline double helixCentrePhi(double xC, double yC) const {
0256     // Faster version of helixCentrePhi()
0257     return xC == 0.0 && yC == 0.0 ? 0.0 : std::atan2(yC, xC);
0258   }
0259
0260   /// Is the vertex inside the cylinder ? (stricly inside : true)
0261   inline bool inside() const {
0262     return (particle_.R2() < rCyl2 - 0.00001 * rCyl && fabs(particle_.Z()) < zCyl - 0.00001);
0263   }
0264
0265   inline bool inside(double rPos2) const {
0266     return (rPos2 < rCyl2 - 0.00001 * rCyl && fabs(particle_.Z()) < zCyl - 0.00001);
0267   }
0268
0269   /// Is the vertex already on the cylinder surface ?
0270   inline bool onSurface() const { return (onBarrel() || onEndcap()); }
0271
0272   inline bool onSurface(double rPos2) const { return (onBarrel(rPos2) || onEndcap(rPos2)); }
0273
0274   /// Is the vertex already on the cylinder barrel ?
0275   inline bool onBarrel() const {
0276     double rPos2 = particle_.R2();
0277     return (fabs(rPos2 - rCyl2) < 0.00001 * rCyl && fabs(particle_.Z()) <= zCyl);
0278   }
0279
0280   inline bool onBarrel(double rPos2) const {
0281     return (fabs(rPos2 - rCyl2) < 0.00001 * rCyl && fabs(particle_.Z()) <= zCyl);
0282   }
0283
0284   /// Is the vertex already on the cylinder endcap ?
0285   inline bool onEndcap() const { return (fabs(fabs(particle_.Z()) - zCyl) < 0.00001 && particle_.R2() <= rCyl2); }
0286
0287   inline bool onEndcap(double rPos2) const { return (fabs(fabs(particle_.Z()) - zCyl) < 0.00001 && rPos2 <= rCyl2); }
0288
0289   /// Is the vertex on some material ?
0290   inline bool onFiducial() const { return fiducial; }
0291
0292   /// Has the particle decayed while propagated ?
0293   inline bool hasDecayed() const { return decayed; }
0294
0295   /// Has propagation been performed and was barrel or endcap reached ?
0296   inline int getSuccess() const { return success; }
0297
0298   /// Set the magnetic field
0299   inline void setMagneticField(double b) { bField = b; }
0300
0301   /// Get the magnetic field
0302   inline double getMagneticField() const { return bField; }
0303
0304   /// Set the debug leve;
0305   inline void setDebug() { debug = true; }
0306   inline void resetDebug() { debug = false; }
0307 };
0308
0309 #endif