ParticlePropagator.cc

Go to the documentation of this file.

//CMSSW headers
#include "DataFormats/GeometrySurface/interface/BoundDisk.h"
#include "DataFormats/GeometrySurface/interface/BoundCylinder.h"

//FAMOS headers
#include "FastSimulation/ParticlePropagator/interface/ParticlePropagator.h"
#include "FastSimulation/ParticlePropagator/interface/MagneticFieldMap.h"
#include "FastSimulation/TrackerSetup/interface/TrackerLayer.h"
#include "FastSimulation/Event/interface/FSimTrack.h"
#include "FastSimulation/Event/interface/FSimVertex.h"
#include "FastSimulation/Particle/interface/pdg_functions.h"
#include "FastSimulation/Particle/interface/makeParticle.h"
#include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"

ParticlePropagator::ParticlePropagator() : 
  BaseParticlePropagator(), random(nullptr) {;}

ParticlePropagator::ParticlePropagator(const RawParticle& myPart,
                       double RCyl, double ZCyl, 
                       const MagneticFieldMap* aFieldMap,
                       const RandomEngineAndDistribution* engine,
                                       const HepPDT::ParticleDataTable* table) :
  BaseParticlePropagator(myPart,RCyl,ZCyl,0.),
  theFieldMap(aFieldMap),
  random(engine),
  theTable(table)
{
  setMagneticField(fieldMap(particle().X(),particle().Y(),particle().Z()));
  initProperDecayTime();
}

ParticlePropagator::ParticlePropagator( const RawParticle& myPart,
                    const MagneticFieldMap* aFieldMap,
                    const RandomEngineAndDistribution* engine,
                                        const HepPDT::ParticleDataTable* table) :
  BaseParticlePropagator(myPart,0.,0.,0.),
  theFieldMap(aFieldMap),
  random(engine),
  theTable(table)
 
{
  setMagneticField(fieldMap(particle().X(),particle().Y(),particle().Z()));
  initProperDecayTime();
}

ParticlePropagator::ParticlePropagator(const XYZTLorentzVector& mom, 
                       const XYZTLorentzVector& vert, float q,
                       const MagneticFieldMap* aFieldMap,
                                       const HepPDT::ParticleDataTable* table) :
  BaseParticlePropagator(RawParticle(mom,vert,q),0.,0.,0.),
  theFieldMap(aFieldMap),
  random(nullptr),
  theTable(table)
{
  setMagneticField(fieldMap(particle().X(),particle().Y(),particle().Z()));
}

ParticlePropagator::ParticlePropagator(const XYZTLorentzVector& mom, 
                       const XYZVector& vert, float q,
                       const MagneticFieldMap* aFieldMap,
                                       const HepPDT::ParticleDataTable* table) :
  BaseParticlePropagator(
  RawParticle(mom,XYZTLorentzVector(vert.X(),vert.Y(),vert.Z(),0.0), q),0.,0.,0.),
  theFieldMap(aFieldMap),
  random(nullptr),
  theTable(table)
{
  setMagneticField(fieldMap(particle().X(),particle().Y(),particle().Z()));
}

ParticlePropagator::ParticlePropagator(const FSimTrack& simTrack,
                       const MagneticFieldMap* aFieldMap,
                       const RandomEngineAndDistribution* engine,
                                       const HepPDT::ParticleDataTable* table) :
  BaseParticlePropagator(makeParticle(table,
                                      simTrack.type(),
                                      simTrack.momentum(),
                                      simTrack.vertex().position()),
             0.,0.,0.),
  theFieldMap(aFieldMap),
  random(engine),
  theTable(table)
{
  setMagneticField(fieldMap(particle().X(),particle().Y(),particle().Z()));
  if ( simTrack.decayTime() < 0. ) { 
    if ( simTrack.nDaughters() ) 
      // This particle already decayed, don't decay it twice
      this->setProperDecayTime(1E99);
    else
      // This particle hasn't decayed yet. Decay time according to particle lifetime
      initProperDecayTime();
  } else {
    // Decay time pre-defined at generator level 
    this->setProperDecayTime(simTrack.decayTime());
  }
}

ParticlePropagator::ParticlePropagator(const ParticlePropagator& myPropPart) :
  BaseParticlePropagator(myPropPart),
  theFieldMap(myPropPart.theFieldMap)
{  
  //  setMagneticField(fieldMap(x(),y(),z()));
}

ParticlePropagator::ParticlePropagator(const BaseParticlePropagator& myPropPart,
                       const MagneticFieldMap* aFieldMap,
                                       const HepPDT::ParticleDataTable* table) :
  BaseParticlePropagator(myPropPart),
  theFieldMap(aFieldMap),
  theTable(table)
{  
  setMagneticField(fieldMap(particle().X(),particle().Y(),particle().Z()));
}


void 
ParticlePropagator::initProperDecayTime() {

  // And this is the proper time at which the particle will decay
  double properDecayTime = 
    (particle().pid()==0||particle().pid()==22||abs(particle().pid())==11||abs(particle().pid())==2112||abs(particle().pid())==2212||
     !random) ?
    1E99 : -pdg::cTau(particle().pid(), theTable) * std::log(random->flatShoot());

  this->setProperDecayTime(properDecayTime);

}


bool
ParticlePropagator::propagateToClosestApproach(double x0, double y0, bool first) {
  setMagneticField(fieldMap(0.,0.,0.));
  return BaseParticlePropagator::propagateToClosestApproach(x0,y0,first);
}

bool
ParticlePropagator::propagateToNominalVertex(const XYZTLorentzVector& v) {
  setMagneticField(fieldMap(0.,0.,0.));
  return BaseParticlePropagator::propagateToNominalVertex(v);
}

ParticlePropagator
ParticlePropagator::propagated() const {
  return ParticlePropagator(BaseParticlePropagator::propagated(),theFieldMap,theTable);
}

double
ParticlePropagator::fieldMap(double xx,double yy, double zz) {
  // Arguments now passed in cm.
  //  return MagneticFieldMap::instance()->inTesla(GlobalPoint(xx/10.,yy/10.,zz/10.)).z();
  // Return a dummy value for neutral particles!
  return particle().charge() == 0.0 || theFieldMap == nullptr ? 
    4. : theFieldMap->inTeslaZ(GlobalPoint(xx,yy,zz));
}

double
ParticlePropagator::fieldMap(const TrackerLayer& layer, double coord, int success) {
  // Arguments now passed in cm.
  //  return MagneticFieldMap::instance()->inTesla(GlobalPoint(xx/10.,yy/10.,zz/10.)).z();
  // Return a dummy value for neutral particles!
  return particle().charge() == 0.0 || theFieldMap == nullptr ? 
    4. : theFieldMap->inTeslaZ(layer,coord,success);
}

bool
ParticlePropagator::propagateToBoundSurface(const TrackerLayer& layer) {


  fiducial = true;
  BoundDisk const* disk = layer.disk();
  //  bool disk = layer.forward();
  //  double innerradius=-999;
  double innerradius = disk ? layer.diskInnerRadius() : -999. ;

  //  if( disk ) {
    //    const Surface& surface = layer.surface();
    //    const BoundDisk & myDisk = dynamic_cast<const BoundDisk&>(surface);
    //    innerradius=myDisk.innerRadius();   
  //    innerradius=myDisk->innerRadius();    
  //  }

  bool done = propagate();

  // Set the magnetic field at the new location (if succesfully propagated)
  if ( done && !hasDecayed() ) { 
    if ( success == 2 ) 
      setMagneticField(fieldMap(layer,particle().r(),success));
    else if ( success == 1 )
      setMagneticField(fieldMap(layer,particle().z(),success));
  } 
       
  // There is some real material here
  fiducial = !(!disk &&  success!=1) &&
    !( disk && (success!=2  || particle().r()<innerradius));

  return done;
}

void 
ParticlePropagator::setPropagationConditions(const TrackerLayer& layer, 
                         bool firstLoop) { 
  // Set the magentic field
  // setMagneticField(fieldMap(x(),y(),z()));

  // Set R and Z according to the Tracker Layer characteristics.
  //  const Surface& surface = layer.surface();

  if( layer.forward() ) {

    //    const BoundDisk & myDisk = dynamic_cast<const BoundDisk&>(surface);
    // ParticlePropagator works in mm, whereas the detector geometry is in cm
    BaseParticlePropagator::setPropagationConditions(
                                  layer.diskOuterRadius(),
                  fabs(layer.disk()->position().z()),
                  firstLoop);       

  // ... or if it is a cylinder barrel 
  } else {

    //    const BoundCylinder & myCylinder = dynamic_cast<const BoundCylinder &>(surface);
    // ParticlePropagator works now in cm
    BaseParticlePropagator::setPropagationConditions(
                     layer.cylinder()->bounds().width()/2.,
                     layer.cylinder()->bounds().length()/2.,
                     firstLoop);
  }

}