BaseParticlePropagator.h

Go to the documentation of this file.

/*Emacs: -*- C++ -*- */
#ifndef BASEPARTICLEPROPAGATOR_H
#define BASEPARTICLEPROPAGATOR_H

//FAMOS
#include "FastSimulation/Particle/interface/RawParticle.h"

class BaseParticlePropagator : public RawParticle {
  
public:

  BaseParticlePropagator();

  BaseParticlePropagator(const RawParticle& myPart, 
             double r, double z, double B);
  BaseParticlePropagator(const RawParticle& myPart, 
             double r, double z, double B, double t);

  void init();

  bool propagate();

  bool backPropagate();

  BaseParticlePropagator propagated() const;

  bool propagateToClosestApproach(double x0=0.,double y0=0,bool first=true);
  bool propagateToEcal(bool first=true);
  bool propagateToPreshowerLayer1(bool first=true);
  bool propagateToPreshowerLayer2(bool first=true);
  bool propagateToEcalEntrance(bool first=true);
  bool propagateToHcalEntrance(bool first=true);
  bool propagateToVFcalEntrance(bool first=true);
  bool propagateToHcalExit(bool first=true);
  bool propagateToHOLayer(bool first=true);
  bool propagateToNominalVertex(const XYZTLorentzVector& hit2=XYZTLorentzVector(0.,0.,0.,0.));
  bool propagateToBeamCylinder(const XYZTLorentzVector& v, double radius=0.); 
  void setPropagationConditions(double r, double z, bool firstLoop=true);

private:
  //  RawParticle   particle;
  double rCyl;
  double rCyl2;
  double zCyl;
  double bField;
  double properDecayTime;
  bool debug;

protected:
  int success;
  bool fiducial;

  inline double c_light() const { return 299.792458; }

private:
  bool firstLoop;
  bool decayed;
  double properTime;
  int propDir;

public:

  inline void setProperDecayTime(double t) { properDecayTime = t; }

  inline void increaseRCyl(double delta) {rCyl = rCyl + delta; rCyl2 = rCyl*rCyl; }

  double xyImpactParameter(double x0=0., double y0=0.) const;

  inline double zImpactParameter(double x0=0, double y0=0.) const {
    // Longitudinal impact parameter
    return Z() - Pz() * std::sqrt( ((X()-x0)*(X()-x0) + (Y()-y0)*(Y()-y0) ) / Perp2());
  }

  inline double helixRadius() const { 
    // The helix Radius
    //
    // The helix' Radius sign accounts for the orientation of the magnetic field 
    // (+ = along z axis) and the sign of the electric charge of the particle. 
    // It signs the rotation of the (charged) particle around the z axis: 
    // Positive means anti-clockwise, negative means clockwise rotation.
    //
    // The radius is returned in cm to match the units in RawParticle.
    return charge() == 0 ? 0.0 : - Pt() / ( c_light() * 1e-5 * bField * charge() );
  }

  inline double helixRadius(double pT) const { 
    // a faster version of helixRadius, once Perp() has been computed
    return charge() == 0 ? 0.0 : - pT / ( c_light() * 1e-5 * bField * charge() );
  }

  inline double helixStartPhi() const { 
    // The azimuth of the momentum at the vertex
    return Px() == 0.0 && Py() == 0.0 ? 0.0 : std::atan2(Py(),Px());
  }
  
  inline double helixCentreX() const { 
    // The x coordinate of the helix axis
    return X() - helixRadius() * std::sin ( helixStartPhi() );
  }

  inline double helixCentreX(double radius, double phi) const { 
    // Fast version of helixCentreX()
    return X() - radius * std::sin (phi);
  }

  inline double helixCentreY() const { 
    // The y coordinate of the helix axis
    return Y() + helixRadius() * std::cos ( helixStartPhi() );
}

  inline double helixCentreY(double radius, double phi) const { 
    // Fast version of helixCentreX()
    return Y() + radius * std::cos (phi);
  }

  inline double helixCentreDistToAxis() const { 
    // The distance between the cylinder and the helix axes
    double xC = helixCentreX();
    double yC = helixCentreY();
    return std::sqrt( xC*xC + yC*yC );
  }

  inline double helixCentreDistToAxis(double xC, double yC) const { 
    // Faster version of helixCentreDistToAxis
    return std::sqrt( xC*xC + yC*yC );
  }

  inline double helixCentrePhi() const { 
    // The azimuth if the vector joining the cylinder and the helix axes
    double xC = helixCentreX();
    double yC = helixCentreY();
    return xC == 0.0 && yC == 0.0 ? 0.0 : std::atan2(yC,xC);
  }
  
  inline double helixCentrePhi(double xC, double yC) const { 
    // Faster version of helixCentrePhi() 
    return xC == 0.0 && yC == 0.0 ? 0.0 : std::atan2(yC,xC);
  }

  inline bool inside() const {
    return (R2()<rCyl2-0.00001*rCyl && fabs(Z())<zCyl-0.00001);}

  inline bool inside(double rPos2) const {
    return (rPos2<rCyl2-0.00001*rCyl && fabs(Z())<zCyl-0.00001);}


  inline bool onSurface() const {
    return ( onBarrel() || onEndcap() ); 
  }

  inline bool onSurface(double rPos2) const {
    return ( onBarrel(rPos2) || onEndcap(rPos2) ); 
  }

  inline bool onBarrel() const {
    double rPos2 = R2();
    return ( fabs(rPos2-rCyl2) < 0.00001*rCyl && fabs(Z()) <= zCyl );
  }

  inline bool onBarrel(double rPos2) const {
    return ( fabs(rPos2-rCyl2) < 0.00001*rCyl && fabs(Z()) <= zCyl );
  }

  inline bool onEndcap() const {
    return ( fabs(fabs(Z())-zCyl) < 0.00001 && R2() <= rCyl2 ); 
  }
  
  inline bool onEndcap(double rPos2) const {
    return ( fabs(fabs(Z())-zCyl) < 0.00001 && rPos2 <= rCyl2 ); 
  }

  inline bool onFiducial() const { return fiducial; }

  inline bool hasDecayed() const { return decayed; }

  inline int  getSuccess() const { return success;  }

  inline void setMagneticField(double b) {  bField=b; }

  inline double getMagneticField() const {  return bField; }

  inline void setDebug() { debug = true; } 
  inline void resetDebug() { debug = false; }
    
};

#endif