BaseParticlePropagator.h

Go to the documentation of this file.

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

//FAMOS
#include "CommonTools/BaseParticlePropagator/interface/RawParticle.h"

class BaseParticlePropagator {
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_;
  //  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 RawParticle const& particle() const { return particle_; }
  inline RawParticle& particle() { return particle_; }
  void setParticle(RawParticle const& iParticle) { particle_ = iParticle; }

  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 particle_.Z() - particle_.Pz() * std::sqrt(((particle_.X() - x0) * (particle_.X() - x0) +
                                                       (particle_.Y() - y0) * (particle_.Y() - y0)) /
                                                      particle_.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 particle_.charge() == 0 ? 0.0 : -particle_.Pt() / (c_light() * 1e-5 * bField * particle_.charge());
  }

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

  inline double helixStartPhi() const {
    // The azimuth of the momentum at the vertex
    return particle_.Px() == 0.0 && particle_.Py() == 0.0 ? 0.0 : std::atan2(particle_.Py(), particle_.Px());
  }

  inline double helixCentreX() const {
    // The x coordinate of the helix axis
    return particle_.X() - helixRadius() * std::sin(helixStartPhi());
  }

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

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

  inline double helixCentreY(double radius, double phi) const {
    // Fast version of helixCentreX()
    return particle_.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 (particle_.R2() < rCyl2 - 0.00001 * rCyl && fabs(particle_.Z()) < zCyl - 0.00001);
  }

  inline bool inside(double rPos2) const {
    return (rPos2 < rCyl2 - 0.00001 * rCyl && fabs(particle_.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 = particle_.R2();
    return (fabs(rPos2 - rCyl2) < 0.00001 * rCyl && fabs(particle_.Z()) <= zCyl);
  }

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

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

  inline bool onEndcap(double rPos2) const { return (fabs(fabs(particle_.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