---

TrajectoryManager Class Reference

#include <FastSimulation/TrajectoryManager/interface/TrajectoryManager.h>

List of all members.

Public Member Functions
void	createPSimHits (const TrackerLayer &layer, const ParticlePropagator &P_before, std::map< double, PSimHit > &theHitMap, int trackID, int partID)
	Create a vector of PSimHits.
void	initializeRecoGeometry (const GeometricSearchTracker *geomSearchTracker, const TrackerInteractionGeometry *interactionGeometry, const MagneticFieldMap *aFieldMap)
	Initialize the Reconstruction Geometry.
void	initializeTrackerGeometry (const TrackerGeometry *geomTracker)
	Initialize the full Tracker Geometry.
void	loadSimHits (edm::PSimHitContainer &c) const
void	propagateToCalorimeters (ParticlePropagator &PP, int fsimi)
	Propagate the particle through the calorimeters.
bool	propagateToLayer (ParticlePropagator &PP, unsigned layer)
	Propagate a particle to a given tracker layer (for electron pixel matching mostly).
void	reconstruct ()
	Does the real job.
const TrackerInteractionGeometry *	theGeometry ()
	Returns the pointer to geometry.
	TrajectoryManager (FSimEvent *aSimEvent, const edm::ParameterSet &matEff, const edm::ParameterSet &simHits, const edm::ParameterSet &decays, const RandomEngine *engine)
	Constructor from a FSimEvent.
	TrajectoryManager ()
	Default Constructor.
	~TrajectoryManager ()
	Default Destructor.
Private Member Functions
const DetLayer *	detLayer (const TrackerLayer &layer, float zpos) const
	Returns the DetLayer pointer corresponding to the FAMOS layer.
void	initializeLayerMap ()
	Initialize correspondence map between Famos interaction geometry and tracker reco geometry.
void	makePSimHits (const GeomDet *det, const TrajectoryStateOnSurface &ts, std::map< double, PSimHit > &theHitMap, int tkID, float el, float thick, int pID)
	and there
std::pair< double, PSimHit >	makeSinglePSimHit (const GeomDetUnit &det, const TrajectoryStateOnSurface &ts, int tkID, float el, float thick, int pID) const
	and there
TrajectoryStateOnSurface	makeTrajectoryState (const DetLayer *layer, const ParticlePropagator &pp, const MagneticField *field) const
	Teddy, you must put comments there.
void	updateWithDaughters (ParticlePropagator &PP, int fsimi)
	Decay the particle and update the SimEvent with daughters.
Private Attributes
const MagneticFieldMap *	_theFieldMap
const TrackerInteractionGeometry *	_theGeometry
double	distCut
bool	firstLoop
Pythia6Decays *	myDecayEngine
FSimEvent *	mySimEvent
double	pTmin
const RandomEngine *	random
const GeometricSearchTracker *	theGeomSearchTracker
const TrackerGeometry *	theGeomTracker
std::vector< const DetLayer * >	theLayerMap
MaterialEffects *	theMaterialEffects
int	theNegLayerOffset
std::map< unsigned, std::map < double, PSimHit > >	thePSimHits

---

Detailed Description

Definition at line 57 of file TrajectoryManager.h.

---

Constructor & Destructor Documentation

TrajectoryManager::TrajectoryManager

(

)

[inline]

Default Constructor.

Definition at line 63 of file TrajectoryManager.h.

{;}

TrajectoryManager::TrajectoryManager	(	FSimEvent *	aSimEvent,
		const edm::ParameterSet &	matEff,
		const edm::ParameterSet &	simHits,
		const edm::ParameterSet &	decays,
		const RandomEngine *	engine
	)

Constructor from a FSimEvent.

Definition at line 50 of file TrajectoryManager.cc.

References distCut, firstLoop, RandomEngine::flatShoot(), edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), int, myDecayEngine, pTmin, random, and theMaterialEffects.

                                                                 : 
  mySimEvent(aSimEvent), 
  _theGeometry(0),
  _theFieldMap(0),
  theMaterialEffects(0), 
  myDecayEngine(0), 
  theGeomTracker(0),
  theGeomSearchTracker(0),
  theLayerMap(56, static_cast<const DetLayer*>(0)), // reserve space for layers here
  theNegLayerOffset(27),
  //  myHistos(0),
  random(engine)

{
  
  // Initialize Bthe stable particle decay engine 
  if ( decays.getParameter<bool>("ActivateDecays") ) { 
    int seed = (int) ( 900000000. * random->flatShoot() );
    double comE = decays.getParameter<double>("comEnergy");
    myDecayEngine = new Pythia6Decays(seed,comE);
    distCut = decays.getParameter<double>("DistCut");
  }

  // Initialize the Material Effects updator, if needed
  if ( matEff.getParameter<bool>("PairProduction") || 
       matEff.getParameter<bool>("Bremsstrahlung") ||
       matEff.getParameter<bool>("EnergyLoss") || 
       matEff.getParameter<bool>("MultipleScattering") )
       theMaterialEffects = new MaterialEffects(matEff,random);

  // Save SimHits according to Optiom
  // Only the hits from first half loop is saved
  firstLoop = simHits.getUntrackedParameter<bool>("firstLoop",true);
  // Only if pT>pTmin are the hits saved
  pTmin = simHits.getUntrackedParameter<double>("pTmin",0.5);

  // Get the Famos Histos pointer
  //  myHistos = Histos::instance();

  // Initialize a few histograms
  /* 
  myHistos->book("h300",1210,-121.,121.,1210,-121.,121.);
  myHistos->book("h301",1200,-300.,300.,1210,-121.,121.);
  */

  
}

TrajectoryManager::~TrajectoryManager

(

)

Default Destructor.

Definition at line 133 of file TrajectoryManager.cc.

References myDecayEngine, and theMaterialEffects.

                                      {

  if ( myDecayEngine ) delete myDecayEngine;
  if ( theMaterialEffects ) delete theMaterialEffects;

  //Write the histograms
  //myHistos->put("histos.root");
  //  if ( myHistos ) delete myHistos;

}

---

Member Function Documentation

void TrajectoryManager::createPSimHits	(	const TrackerLayer &	layer,
		const ParticlePropagator &	P_before,
		std::map< double, PSimHit > &	theHitMap,
		int	trackID,
		int	partID
	)

Create a vector of PSimHits.

Definition at line 478 of file TrajectoryManager.cc.

References anyDirection, GeometricSearchDet::compatibleDets(), detLayer(), MaterialEffects::energyLoss(), BaseParticlePropagator::getMagneticField(), i, makePSimHits(), makeTrajectoryState(), theMaterialEffects, MaterialEffects::thickness(), and RawParticle::Z().

Referenced by reconstruct().

                                                           {

  // Propagate the particle coordinates to the closest tracker detector(s) 
  // in this layer and create the PSimHit(s)

  //  const MagneticField& mf = MagneticFieldMap::instance()->magneticField();
  // This solution is actually much faster !
  LocalMagneticField mf(PP.getMagneticField());
  AnalyticalPropagator alongProp(&mf, anyDirection);
  InsideBoundsMeasurementEstimator est;

  typedef GeometricSearchDet::DetWithState   DetWithState;
  const DetLayer* tkLayer = detLayer(layer,PP.Z());

  TrajectoryStateOnSurface trajState = makeTrajectoryState( tkLayer, PP, &mf);
  float thickness = theMaterialEffects ? theMaterialEffects->thickness() : 0.;
  float eloss = theMaterialEffects ? theMaterialEffects->energyLoss() : 0.;

  // Find, in the corresponding layers, the detectors compatible 
  // with the current track 
  std::vector<DetWithState> compat 
    = tkLayer->compatibleDets( trajState, alongProp, est);

  // And create the corresponding PSimHits
  std::map<double,PSimHit> theTrackHits;
  for (std::vector<DetWithState>::const_iterator i=compat.begin(); i!=compat.end(); i++) {
    // Correct Eloss for last 3 rings of TEC (thick sensors, 0.05 cm)
    // Disgusting fudge factor ! 
      makePSimHits( i->first, i->second, theHitMap, trackID, eloss, thickness, partID);
  }
}

const DetLayer * TrajectoryManager::detLayer	(	const TrackerLayer &	layer,
		float	zpos
	)			const [private]

Returns the DetLayer pointer corresponding to the FAMOS layer.

Definition at line 842 of file TrajectoryManager.cc.

References TrackerLayer::forward(), TrackerLayer::layerNumber(), theLayerMap, and theNegLayerOffset.

Referenced by createPSimHits().

{
  if (zpos > 0 || !layer.forward() ) return theLayerMap[layer.layerNumber()];
  else return theLayerMap[layer.layerNumber()+theNegLayerOffset];
}

void TrajectoryManager::initializeLayerMap

(

)

[private]

Initialize correspondence map between Famos interaction geometry and tracker reco geometry.

ATTENTION: HARD CODED LOGIC! If Famos layer numbering changes this logic needs to be adapted to the new numbering!

Definition at line 736 of file TrajectoryManager.cc.

References _theGeometry, GeometricSearchTracker::barrelLayers(), TrackerInteractionGeometry::cylinderBegin(), TrackerInteractionGeometry::cylinderEnd(), muonGeometry::disk, i, BoundDisk::innerRadius(), LogDebug, GeometricSearchTracker::negForwardLayers(), BoundDisk::outerRadius(), GeometricSearchTracker::posForwardLayers(), GloballyPositioned< T >::position(), Cylinder::radius(), theGeomSearchTracker, theLayerMap, and theNegLayerOffset.

Referenced by initializeRecoGeometry().

{

// These are the BoundSurface&, the BoundDisk* and the BoundCylinder* for that layer
//   const BoundSurface& theSurface = layer.surface();
//   BoundDisk* theDisk = layer.disk();  // non zero for endcaps
//   BoundCylinder* theCylinder = layer.cylinder(); // non zero for barrel
//   int theLayer = layer.layerNumber(); // 1->3 PixB, 4->5 PixD, 
//                                       // 6->9 TIB, 10->12 TID, 
//                                       // 13->18 TOB, 19->27 TEC


  std::vector< BarrelDetLayer*>   barrelLayers = 
    theGeomSearchTracker->barrelLayers();
  LogDebug("FastTracker") << "Barrel DetLayer dump: ";
  for (std::vector< BarrelDetLayer*>::const_iterator bl=barrelLayers.begin();
       bl != barrelLayers.end(); ++bl) {
    LogDebug("FastTracker")<< "radius " << (**bl).specificSurface().radius(); 
  }

  std::vector< ForwardDetLayer*>  posForwardLayers = 
    theGeomSearchTracker->posForwardLayers();
  LogDebug("FastTracker") << "Positive Forward DetLayer dump: ";
  for (std::vector< ForwardDetLayer*>::const_iterator fl=posForwardLayers.begin();
       fl != posForwardLayers.end(); ++fl) {
    LogDebug("FastTracker") << "Z pos "
                            << (**fl).surface().position().z()
                            << " radii " 
                            << (**fl).specificSurface().innerRadius() 
                            << ", " 
                            << (**fl).specificSurface().outerRadius(); 
  }

  const float rTolerance = 1.5;
  const float zTolerance = 3.;

  LogDebug("FastTracker")<< "Dump of TrackerInteractionGeometry cylinders:";
  for( std::list<TrackerLayer>::const_iterator i=_theGeometry->cylinderBegin();
       i!=_theGeometry->cylinderEnd(); ++i) {
    const BoundCylinder* cyl = i->cylinder();
    const BoundDisk* disk = i->disk();

    LogDebug("FastTracker") << "Famos Layer no " << i->layerNumber()
                            << " is sensitive? " << i->sensitive()
                            << " pos " << i->surface().position();
    if (!i->sensitive()) continue;

    if (cyl != 0) {
      LogDebug("FastTracker") << " cylinder radius " << cyl->radius();
      bool found = false;
      for (std::vector< BarrelDetLayer*>::const_iterator 
             bl=barrelLayers.begin(); bl != barrelLayers.end(); ++bl) {
        if (fabs( cyl->radius() - (**bl).specificSurface().radius()) < rTolerance) {
          theLayerMap[i->layerNumber()] = *bl;
          found = true;
          LogDebug("FastTracker")<< "Corresponding DetLayer found with radius "
                                 << (**bl).specificSurface().radius();
          break;
        }
      }
      if (!found) {
        edm::LogError("FastTracker") << "FAILED to find a corresponding DetLayer!";
      }
    }
    else {
      LogDebug("FastTracker") << " disk radii " << disk->innerRadius() 
                 << ", " << disk->outerRadius();
      bool found = false;
      for (std::vector< ForwardDetLayer*>::const_iterator fl=posForwardLayers.begin();
           fl != posForwardLayers.end(); ++fl) {
        if (fabs( disk->position().z() - (**fl).surface().position().z()) < zTolerance) {
          theLayerMap[i->layerNumber()] = *fl;
          found = true;
          LogDebug("FastTracker") << "Corresponding DetLayer found with Z pos "
                                  << (**fl).surface().position().z()
                                  << " and radii " 
                                  << (**fl).specificSurface().innerRadius() 
                                  << ", " 
                                  << (**fl).specificSurface().outerRadius(); 
          break;
        }
      }
      if (!found) {
        edm::LogError("FastTracker") << "FAILED to find a corresponding DetLayer!";
      }
    }
  }

  // Put the negative layers in the same map but with an offset
  std::vector< ForwardDetLayer*>  negForwardLayers = theGeomSearchTracker->negForwardLayers();
  for (std::vector< ForwardDetLayer*>::const_iterator nl=negForwardLayers.begin();
       nl != negForwardLayers.end(); ++nl) {
    for (int i=0; i<=theNegLayerOffset; i++) {
      if (theLayerMap[i] == 0) continue;
      if ( fabs( (**nl).surface().position().z() +theLayerMap[i]-> surface().position().z()) < zTolerance) {
        theLayerMap[i+theNegLayerOffset] = *nl;
        break;
      }
    }
  }  

}

void TrajectoryManager::initializeRecoGeometry	(	const GeometricSearchTracker *	geomSearchTracker,
		const TrackerInteractionGeometry *	interactionGeometry,
		const MagneticFieldMap *	aFieldMap
	)

Initialize the Reconstruction Geometry.

Definition at line 103 of file TrajectoryManager.cc.

References _theFieldMap, _theGeometry, initializeLayerMap(), and theGeomSearchTracker.

Referenced by FamosManager::setupGeometryAndField().

{
  
  // Initialize the reco tracker geometry
  theGeomSearchTracker = geomSearchTracker;
  
  // Initialize the simplified tracker geometry
  _theGeometry = interactionGeometry;

  initializeLayerMap();

  // Initialize the magnetic field
  _theFieldMap = aFieldMap;

}

void TrajectoryManager::initializeTrackerGeometry

(

const TrackerGeometry *

geomTracker

)

Initialize the full Tracker Geometry.

Definition at line 122 of file TrajectoryManager.cc.

References theGeomTracker.

Referenced by FamosManager::setupGeometryAndField().

                                                                               { 
  
  theGeomTracker = geomTracker;

}

void TrajectoryManager::loadSimHits

(

edm::PSimHitContainer &

c

)

const

Definition at line 849 of file TrajectoryManager.cc.

References begin, end, it, and thePSimHits.

Referenced by FamosProducer::produce().

{

  std::map<unsigned,std::map<double,PSimHit> >::const_iterator itrack = thePSimHits.begin();
  std::map<unsigned,std::map<double,PSimHit> >::const_iterator itrackEnd = thePSimHits.end();
  for ( ; itrack != itrackEnd; ++itrack ) {
    std::map<double,PSimHit>::const_iterator it = (itrack->second).begin();
    std::map<double,PSimHit>::const_iterator itEnd = (itrack->second).end();
    for( ; it!= itEnd; ++it) { 
      /*
      DetId theDetUnitId((it->second).detUnitId());
      const GeomDet* theDet = theGeomTracker->idToDet(theDetUnitId);
      std::cout << "Track/z/r after : "
                << (it->second).trackId() << " " 
                << theDet->surface().toGlobal((it->second).localPosition()).z() << " " 
                << theDet->surface().toGlobal((it->second).localPosition()).perp() << std::endl;
      */
      // Keep only those hits that are on the physical volume of a module
      // (The other hits have been assigned a negative <double> value. 
      if ( it->first > 0. ) c.push_back(it->second); 
    }
  }

}

void TrajectoryManager::makePSimHits	(	const GeomDet *	det,
		const TrajectoryStateOnSurface &	ts,
		std::map< double, PSimHit > &	theHitMap,
		int	tkID,
		float	el,
		float	thick,
		int	pID
	)			[private]

and there

Definition at line 526 of file TrajectoryManager.cc.

References GeomDet::components(), i, and makeSinglePSimHit().

Referenced by createPSimHits().

{

  std::vector< const GeomDet*> comp = det->components();
  if (!comp.empty()) {
    for (std::vector< const GeomDet*>::const_iterator i = comp.begin();
         i != comp.end(); i++) {
      const GeomDetUnit* du = dynamic_cast<const GeomDetUnit*>(*i);
      if (du != 0)
        theHitMap.insert(theHitMap.end(),makeSinglePSimHit( *du, ts, tkID, el, thick, pID));
    }
  }
  else {
    const GeomDetUnit* du = dynamic_cast<const GeomDetUnit*>(det);
    if (du != 0)
      theHitMap.insert(theHitMap.end(),makeSinglePSimHit( *du, ts, tkID, el, thick, pID));
  }


}

std::pair< double, PSimHit > TrajectoryManager::makeSinglePSimHit	(	const GeomDetUnit &	det,
		const TrajectoryStateOnSurface &	ts,
		int	tkID,
		float	el,
		float	thick,
		int	pID
	)			const [private]

and there

Definition at line 551 of file TrajectoryManager.cc.

References anyDirection, PV3DBase< T, PVType, FrameType >::basicVector(), BoundSurface::bounds(), FSimTrack::closestDaughterId(), GenMuonPlsPt100GeV_cfg::cout, PXFDetId::disk(), dist(), lat::endl(), cmsRelvalreport::exit, GeomDet::geographicalId(), TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), FSimTrack::id(), TIBDetId::layer(), TOBDetId::layer(), PXBDetId::layer(), Bounds::length(), TrajectoryStateOnSurface::localMomentum(), TrajectoryStateOnSurface::localPosition(), PV3DBase< T, PVType, FrameType >::mag(), module(), FSimTrack::mother(), mySimEvent, path(), TECDetId::petal(), FSimVertex::position(), GloballyPositioned< T >::position(), DetId::rawId(), TIDDetId::ring(), TECDetId::ring(), TIDDetId::side(), SiStripDetId::stereo(), GeomDet::surface(), theLayer, Bounds::thickness(), Surface::toGlobal(), GeomDet::toLocal(), FBaseSimEvent::track(), TrajectoryStateOnSurface::transverseCurvature(), Vector3DBase< T, FrameTag >::unit(), FSimTrack::vertex(), TIDDetId::wheel(), TECDetId::wheel(), Bounds::width(), PV3DBase< T, PVType, FrameType >::z(), and z.

Referenced by makePSimHits().

{

  const float onSurfaceTolarance = 0.01; // 10 microns

  LocalPoint lpos;
  LocalVector lmom;
  if ( fabs( det.toLocal(ts.globalPosition()).z()) < onSurfaceTolarance) {
    lpos = ts.localPosition();
    lmom = ts.localMomentum();
  }
  else {
    HelixArbitraryPlaneCrossing crossing( ts.globalPosition().basicVector(),
                                          ts.globalMomentum().basicVector(),
                                          ts.transverseCurvature(),
                                          anyDirection);
    std::pair<bool,double> path = crossing.pathLength(det.surface());
    if (!path.first) {
      // edm::LogError("FastTracker") << "TrajectoryManager ERROR: crossing with det failed, skipping PSimHit";
      return  std::pair<double,PSimHit>(0.,PSimHit());
    }
    lpos = det.toLocal( GlobalPoint( crossing.position(path.second)));
    lmom = det.toLocal( GlobalVector( crossing.direction(path.second)));
    lmom = lmom.unit() * ts.localMomentum().mag();
  }

  // The module (half) thickness 
  const BoundPlane& theDetPlane = det.surface();
  float halfThick = 0.5*theDetPlane.bounds().thickness();
  // The Energy loss rescaled to the module thickness
  float eloss = el;
  if ( thick > 0. ) {
    // Total thickness is in radiation lengths, 1 radlen = 9.36 cm
    // Sensitive module thickness is about 30 microns larger than 
    // the module thickness itself
    eloss *= (2.* halfThick - 0.003) / (9.36 * thick);
  }
  // The entry and exit points, and the time of flight
  float pZ = lmom.z();
  LocalPoint entry = lpos + (-halfThick/pZ) * lmom;
  LocalPoint exit = lpos + halfThick/pZ * lmom;
  float tof = ts.globalPosition().mag() / 30. ; // in nanoseconds, FIXME: very approximate

  // If a hadron suffered a nuclear interaction, just assign the hits of the closest 
  // daughter to the mother's track. The same applies to a charged particle decay into
  // another charged particle.
  int localTkID = tkID;
  if ( mySimEvent->track(tkID).mother().closestDaughterId() == tkID )
    localTkID = mySimEvent->track(tkID).mother().id();

  // FIXME: fix the track ID and the particle ID
  PSimHit hit( entry, exit, lmom.mag(), tof, eloss, pID,
                  det.geographicalId().rawId(), localTkID,
                  lmom.theta(),
                  lmom.phi());

  // Check that the PSimHit is physically on the module!
  unsigned subdet = DetId(hit.detUnitId()).subdetId(); 
  double boundX = theDetPlane.bounds().width()/2.;
  double boundY = theDetPlane.bounds().length()/2.;

  // Special treatment for TID and TEC trapeziodal modules
  if ( subdet == 4 || subdet == 6 ) 
    boundX *=  1. - hit.localPosition().y()/theDetPlane.position().perp();

#ifdef FAMOS_DEBUG
  unsigned detid  = DetId(hit.detUnitId()).rawId();
  unsigned stereo = 0;
  unsigned theLayer = 0;
  unsigned theRing = 0;
  switch (subdet) { 
  case 1: 
    {
      PXBDetId module(detid);
      theLayer = module.layer();
      std::cout << "\tPixel Barrel Layer " << theLayer << std::endl;
      stereo = 1;
      break;
    }
  case 2: 
    {
      PXFDetId module(detid);
      theLayer = module.disk();
      std::cout << "\tPixel Forward Disk " << theLayer << std::endl;
      stereo = 1;
      break;
    }
  case 3:
    {
      TIBDetId module(detid);
      theLayer  = module.layer();
      std::cout << "\tTIB Layer " << theLayer << std::endl;
      stereo = module.stereo();
      break;
    }
  case 4:
    {
      TIDDetId module(detid);
      theLayer = module.wheel();
      theRing  = module.ring();
      unsigned int theSide = module.side();
      if ( theSide == 1 ) 
        std::cout << "\tTID Petal Back " << std::endl; 
      else
        std::cout << "\tTID Petal Front" << std::endl; 
      std::cout << "\tTID Layer " << theLayer << std::endl;
      std::cout << "\tTID Ring " << theRing << std::endl;
      stereo = module.stereo();
      break;
    }
  case 5:
    {
      TOBDetId module(detid);
      theLayer  = module.layer();
      stereo = module.stereo();
      std::cout << "\tTOB Layer " << theLayer << std::endl;
      break;
    }
  case 6:
    {
      TECDetId module(detid);
      theLayer = module.wheel();
      theRing  = module.ring();
      unsigned int theSide = module.petal()[0];
      if ( theSide == 1 ) 
        std::cout << "\tTEC Petal Back " << std::endl; 
      else
        std::cout << "\tTEC Petal Front" << std::endl; 
      std::cout << "\tTEC Layer " << theLayer << std::endl;
      std::cout << "\tTEC Ring " << theRing << std::endl;
      stereo = module.stereo();
      break;
    }
  default:
    {
      stereo = 0;
      break;
    }
  }
  
  std::cout << "Thickness = " << 2.*halfThick-0.003 << "; " << thick * 9.36 << std::endl
            << "Length    = " << det.surface().bounds().length() << std::endl
            << "Width     = " << det.surface().bounds().width() << std::endl;
    
  std::cout << "Hit position = " 
            << hit.localPosition().x() << " " 
            << hit.localPosition().y() << " " 
            << hit.localPosition().z() << std::endl;
#endif

  // Check if the hit is on the physical volume of the module
  // (It happens that it is not, in the case of double sided modules,
  //  because the envelope of the gluedDet is larger than each of 
  //  the mono and the stereo modules)

  double dist = 0.;
  GlobalPoint IP (mySimEvent->track(localTkID).vertex().position().x(),
                  mySimEvent->track(localTkID).vertex().position().y(),
                  mySimEvent->track(localTkID).vertex().position().z());

  dist = ( fabs(hit.localPosition().x()) > boundX  || 
           fabs(hit.localPosition().y()) > boundY ) ?  
    // Will be used later as a flag to reject the PSimHit!
    -( det.surface().toGlobal(hit.localPosition()) - IP ).mag2() 
    : 
    // These hits are kept!
     ( det.surface().toGlobal(hit.localPosition()) - IP ).mag2();

  // Fill Histos (~poor man event display)
  /* 
     GlobalPoint gpos( det.toGlobal(hit.localPosition()));
     myHistos->fill("h300",gpos.x(),gpos.y());
     if ( sin(gpos.phi()) > 0. ) 
     myHistos->fill("h301",gpos.z(),gpos.perp());
     else
     myHistos->fill("h301",gpos.z(),-gpos.perp());
  */
  
  return std::pair<double,PSimHit>(dist,hit);

}

TrajectoryStateOnSurface TrajectoryManager::makeTrajectoryState	(	const DetLayer *	layer,
		const ParticlePropagator &	pp,
		const MagneticField *	field
	)			const [private]

Teddy, you must put comments there.

Definition at line 514 of file TrajectoryManager.cc.

References RawParticle::charge(), GeometricSearchDet::surface(), Surface::tangentPlane(), RawParticle::X(), RawParticle::Y(), and RawParticle::Z().

Referenced by createPSimHits().

{
  GlobalPoint  pos( pp.X(), pp.Y(), pp.Z());
  GlobalVector mom( pp.Px(), pp.Py(), pp.Pz());
  ReferenceCountingPointer<TangentPlane> plane = layer->surface().tangentPlane(pos);
  return TrajectoryStateOnSurface
    (GlobalTrajectoryParameters( pos, mom, TrackCharge( pp.charge()), field), *plane);
}

void TrajectoryManager::propagateToCalorimeters	(	ParticlePropagator &	PP,
		int	fsimi
	)

Propagate the particle through the calorimeters.

Definition at line 356 of file TrajectoryManager.cc.

References BaseParticlePropagator::getSuccess(), BaseParticlePropagator::hasDecayed(), RawParticle::momentum(), mySimEvent, FSimTrack::notYetToEndVertex(), BaseParticlePropagator::propagateToEcalEntrance(), BaseParticlePropagator::propagateToHcalEntrance(), BaseParticlePropagator::propagateToPreshowerLayer1(), BaseParticlePropagator::propagateToPreshowerLayer2(), BaseParticlePropagator::propagateToVFcalEntrance(), FSimTrack::setEcal(), FSimTrack::setHcal(), FSimTrack::setLayer1(), FSimTrack::setLayer2(), SimTrack::setTkMomentum(), SimTrack::setTkPosition(), FSimTrack::setVFcal(), FBaseSimEvent::track(), updateWithDaughters(), and RawParticle::vertex().

Referenced by reconstruct().

                                                                            {

  FSimTrack& myTrack = mySimEvent->track(fsimi);

  // Set the position and momentum at the end of the tracker volume
  myTrack.setTkPosition(PP.vertex().Vect());
  myTrack.setTkMomentum(PP.momentum());

  // Propagate to Preshower Layer 1 
  PP.propagateToPreshowerLayer1(false);
  if ( PP.hasDecayed() ) {
    updateWithDaughters(PP,fsimi);
    return;
  }
  if ( myTrack.notYetToEndVertex(PP.vertex()) && PP.getSuccess() > 0 )
    myTrack.setLayer1(PP,PP.getSuccess());
  
  // Propagate to Preshower Layer 2 
  PP.propagateToPreshowerLayer2(false);
  if ( PP.hasDecayed() ) { 
    updateWithDaughters(PP,fsimi);
    return;
  }
  if ( myTrack.notYetToEndVertex(PP.vertex()) && PP.getSuccess() > 0 )
    myTrack.setLayer2(PP,PP.getSuccess());

  // Propagate to Ecal Endcap
  PP.propagateToEcalEntrance(false);
  if ( PP.hasDecayed() ) { 
    updateWithDaughters(PP,fsimi);
    return;
  }
  if ( myTrack.notYetToEndVertex(PP.vertex()) )
    myTrack.setEcal(PP,PP.getSuccess());

  // Propagate to HCAL entrance
  PP.propagateToHcalEntrance(false);
  if ( PP.hasDecayed() ) { 
    updateWithDaughters(PP,fsimi);
    return;
  }
  if ( myTrack.notYetToEndVertex(PP.vertex()) )
    myTrack.setHcal(PP,PP.getSuccess());

  // Propagate to VFCAL entrance
  PP.propagateToVFcalEntrance(false);
  if ( PP.hasDecayed() ) { 
    updateWithDaughters(PP,fsimi);
    return;
  }
  if ( myTrack.notYetToEndVertex(PP.vertex()) )
    myTrack.setVFcal(PP,PP.getSuccess());
    
}

bool TrajectoryManager::propagateToLayer	(	ParticlePropagator &	PP,
		unsigned	layer
	)

Propagate a particle to a given tracker layer (for electron pixel matching mostly).

Definition at line 412 of file TrajectoryManager.cc.

References _theGeometry, TrackerInteractionGeometry::cylinderBegin(), TrackerInteractionGeometry::cylinderEnd(), BaseParticlePropagator::getSuccess(), BaseParticlePropagator::onFiducial(), ParticlePropagator::propagateToBoundSurface(), and ParticlePropagator::setPropagationConditions().

                                                                          {

  std::list<TrackerLayer>::const_iterator cyliter;
  bool done = false;

  // Get the geometry elements 
  cyliter = _theGeometry->cylinderBegin();

  // Find the layer to propagate to.      
  for ( ; cyliter != _theGeometry->cylinderEnd() ; ++cyliter ) {

    if ( layer != cyliter->layerNumber() ) continue;
      
    PP.setPropagationConditions(*cyliter);

    done =  
      PP.propagateToBoundSurface(*cyliter) &&
      PP.getSuccess() > 0 && 
      PP.onFiducial();

    break;
    
  }

  return done;

}

void TrajectoryManager::reconstruct

(

)

Does the real job.

Definition at line 145 of file TrajectoryManager.cc.

References _theFieldMap, _theGeometry, BaseRawParticleFilter::accept(), FBaseSimEvent::addSimVertex(), createPSimHits(), TrackerInteractionGeometry::cylinderBegin(), TrackerInteractionGeometry::cylinderEnd(), FBaseSimEvent::filter(), firstLoop, Pythia6Decays::getRandom(), int, MaterialEffects::interact(), python::cmstools::loop(), myDecayEngine, mySimEvent, FSimTrack::notYetToEndVertex(), FSimEvent::nTracks(), propagateToCalorimeters(), pTmin, random, MaterialEffects::save(), Pythia6Decays::saveRandom(), FSimTrack::setPropagate(), theGeomTracker, theMaterialEffects, thePSimHits, FBaseSimEvent::track(), CoreSimTrack::type(), and updateWithDaughters().

Referenced by FamosManager::reconstruct().

{

  // Get pythia random state for decay (after saving current pythia state)
  if ( myDecayEngine ) myDecayEngine->getRandom();
  
  // Clear the hits of the previous event
  //  thePSimHits->clear();
  thePSimHits.clear();

  // The new event
  XYZTLorentzVector myBeamPipe = XYZTLorentzVector(0.,25., 9999999.,0.);

  std::list<TrackerLayer>::const_iterator cyliter;

  //  bool debug = mySimEvent->id().event() == 3;

  // Loop over the particles (watch out: increasing upper limit!)
  for( int fsimi=0; fsimi < (int) mySimEvent->nTracks(); ++fsimi) {

    // If the particle has decayed inside the beampipe, or decays 
    // immediately, there is nothing to do
    if( !mySimEvent->track(fsimi).notYetToEndVertex(myBeamPipe) ) continue;
    mySimEvent->track(fsimi).setPropagate();

    // Get the geometry elements 
    cyliter = _theGeometry->cylinderBegin();
    // Prepare the propagation  
    ParticlePropagator PP(mySimEvent->track(fsimi),_theFieldMap,random);
    //The real work starts here
    int success = 1;
    int sign = +1;
    int loop = 0;
    int cyl = 0;

    // Find the initial cylinder to propagate to.      
    for ( ; cyliter != _theGeometry->cylinderEnd() ; ++cyliter ) {
      
      PP.setPropagationConditions(*cyliter);
      if ( PP.inside() && !PP.onSurface() ) break;
      ++cyl;

    }

    // The particle has a pseudo-rapidity (position or momentum direction) 
    // in excess of 3.0. Just simply go to the last tracker layer
    // without bothering with all the details of the propagation and 
    // material effects.
    // 08/02/06 - pv: increase protection from 0.99 (eta=2.9932) to 0.9998 (eta=4.9517)
    //                to simulate material effects at large eta 
    // if above 0.99: propagate to the last tracker cylinder where the material is concentrated!
    double ppcos2T =  PP.cos2Theta();
    double ppcos2V =  PP.cos2ThetaV();
    if ( ( ppcos2T > 0.99 && ppcos2T < 0.9998 ) && ( cyl == 0 || ( ppcos2V > 0.99 && ppcos2V < 0.9998 ) ) ){ 
      if ( cyliter != _theGeometry->cylinderEnd() ) { 
        cyliter = _theGeometry->cylinderEnd(); 
        --cyliter;
      }
    // if above 0.9998: don't propagate at all (only to the calorimeters directly)
    } else if ( ppcos2T > 0.9998 && ( cyl == 0 || ppcos2V > 0.9998 ) ) { 
      cyliter = _theGeometry->cylinderEnd();
    }
        
    // Loop over the cylinders
    while ( cyliter != _theGeometry->cylinderEnd() &&
            loop<100 &&                            // No more than 100 loops
            mySimEvent->track(fsimi).notYetToEndVertex(PP.vertex())) { // The particle decayed
      
      // Pathological cases:
      // To prevent from interacting twice in a row with the same layer
      //      bool escapeBarrel    = (PP.getSuccess() == -1 && success == 1);
      bool escapeBarrel    = PP.getSuccess() == -1;
      bool escapeEndcap    = (PP.getSuccess() == -2 && success == 1);
      // To break the loop
      bool fullPropagation = 
        (PP.getSuccess() <= 0 && success==0) || escapeEndcap;

      if ( escapeBarrel ) {
        ++cyliter; ++cyl;
        while (cyliter != _theGeometry->cylinderEnd() && cyliter->forward() ) {
          sign=1; ++cyliter; ++cyl;
        }

        if ( cyliter == _theGeometry->cylinderEnd()  ) {
          --cyliter; --cyl; fullPropagation=true; 
        }

      }

      // Define the propagation conditions
      PP.setPropagationConditions(*cyliter,!fullPropagation);
      if ( escapeEndcap ) PP.increaseRCyl(0.0005);

      // Remember last propagation outcome
      success = PP.getSuccess();

      // Propagation was not successful :
      // Change the sign of the cylinder increment and count the loops
      if ( !PP.propagateToBoundSurface(*cyliter) || 
           PP.getSuccess()<=0) {
        sign = -sign;
        ++loop;
      }

      // The particle may have decayed on its way... in which the daughters
      // have to be added to the event record
      if ( PP.hasDecayed() || PP.PDGcTau()<1E-3 ) updateWithDaughters(PP,fsimi);
      if ( PP.hasDecayed() || PP.PDGcTau()<1E-3 ) break;

      // Exit by the endcaps or innermost cylinder :
      // Positive cylinder increment
      if ( PP.getSuccess()==2 || cyliter==_theGeometry->cylinderBegin() ) 
        sign = +1; 
          
      // Successful propagation to a cylinder, with some Material :
      if( PP.getSuccess() > 0 && PP.onFiducial() ) {

        bool saveHit = 
          ( (loop==0 && sign>0) || !firstLoop ) &&   // Save only first half loop
          PP.charge()!=0. &&                         // Consider only charged particles
          cyliter->sensitive() &&                    // Consider only sensitive layers
          PP.Perp2()>pTmin*pTmin;                    // Consider only pT > pTmin

        // Material effects are simulated there
        if ( theMaterialEffects ) 
          theMaterialEffects->interact(*mySimEvent,*cyliter,PP,fsimi); 
        
        if ( saveHit ) { 

          // Consider only active layers
          if ( cyliter->sensitive() ) {
            // Add information to the FSimTrack (not yet available)
            //      myTrack.addSimHit(PP,layer);

            // Return one or two (for overlap regions) PSimHits in the full 
            // tracker geometry
            if ( theGeomTracker ) 
              createPSimHits(*cyliter, PP, thePSimHits[fsimi], fsimi,mySimEvent->track(fsimi).type());

          }
        }

        // Fill Histos (~poor man event display)
        /* 
        myHistos->fill("h300",PP.x(),PP.y());
        if ( sin(PP.vertex().phi()) > 0. ) 
          myHistos->fill("h301",PP.z(),PP.vertex().perp());
        else
          myHistos->fill("h301",PP.z(),-PP.vertex().perp());
        */

        //The particle may have lost its energy in the material
        if ( mySimEvent->track(fsimi).notYetToEndVertex(PP.vertex()) && 
             !mySimEvent->filter().accept(PP)  ) 
          mySimEvent->addSimVertex(PP.vertex(),fsimi);
          
      }

      // Stop here if the particle has reached an end
      if ( mySimEvent->track(fsimi).notYetToEndVertex(PP.vertex()) ) {

        // Otherwise increment the cylinder iterator
        //      do { 
        if (sign==1) {++cyliter;++cyl;}
        else         {--cyliter;--cyl;}

        // Check if the last surface has been reached 
        if( cyliter==_theGeometry->cylinderEnd()) {

          // Try to propagate to the ECAL in half a loop
          // Note: Layer1 = ECAL Barrel entrance, or Preshower
          // entrance, or ECAL Endcap entrance (in the corner)
          PP.propagateToEcal();
          // PP.propagateToPreshowerLayer1();

          // If it is not possible, try go back to the last cylinder
          if(PP.getSuccess()==0) {
            --cyliter; --cyl; sign = -sign;
            PP.setPropagationConditions(*cyliter);
            PP.propagateToBoundSurface(*cyliter);

            // If there is definitely no way, leave it here.
            if(PP.getSuccess()<0) {++cyliter; ++cyl;}

          }

          // Check if the particle has decayed on the way to ECAL
          if ( PP.hasDecayed() )
            updateWithDaughters(PP,fsimi);

        }
      }

    }

    // Propagate all particles without a end vertex to the Preshower, 
    // theECAL and the HCAL.
    if ( mySimEvent->track(fsimi).notYetToEndVertex(PP.vertex()) )
      propagateToCalorimeters(PP,fsimi);

  }

  // Save the information from Nuclear Interaction Generation
  if ( theMaterialEffects ) theMaterialEffects->save();

  // Save pythia random state for decay (and put pythia state for event generation back on the stack)
  if ( myDecayEngine ) myDecayEngine->saveRandom();
  
}

const TrackerInteractionGeometry * TrajectoryManager::theGeometry

(

)

Returns the pointer to geometry.

Definition at line 129 of file TrajectoryManager.cc.

References _theGeometry.

                               {
  return _theGeometry;
}

void TrajectoryManager::updateWithDaughters	(	ParticlePropagator &	PP,
		int	fsimi
	)			[private]

Decay the particle and update the SimEvent with daughters.

Definition at line 441 of file TrajectoryManager.cc.

References FBaseSimEvent::addSimTrack(), FBaseSimEvent::addSimVertex(), RawParticle::charge(), dist(), distCut, myDecayEngine, mySimEvent, Pythia6Decays::particleDaughters(), FSimTrack::setClosestDaughterId(), and FBaseSimEvent::track().

Referenced by propagateToCalorimeters(), and reconstruct().

                                                                        {

  // Decays are not activated : do nothing
  if ( !myDecayEngine ) return;

  // Invoke PYDECY to decay the particle and get the daughters
  const DaughterParticleList& daughters = myDecayEngine->particleDaughters(PP);
  
  // Update the FSimEvent with an end vertex and with the daughters
  if ( daughters.size() ) { 
    double distMin = 1E99;
    int theClosestChargedDaughterId = -1;
    DaughterParticleIterator daughter = daughters.begin();
    
    int ivertex = mySimEvent->addSimVertex(daughter->vertex(),fsimi);

    if ( ivertex != -1 ) {
      for ( ; daughter != daughters.end(); ++daughter) {
        int theDaughterId = mySimEvent->addSimTrack(&(*daughter), ivertex);
        // Find the closest charged daughter (if charged mother)
        if ( PP.charge() * daughter->charge() > 1E-10 ) {
          double dist = (daughter->Vect().Unit().Cross(PP.Vect().Unit())).R();
          if ( dist < distCut && dist < distMin ) { 
            distMin = dist;
            theClosestChargedDaughterId = theDaughterId;
          }
        }
      }
    }
    // Attach mother and closest daughter sp as to cheat tracking ;-)
    if ( theClosestChargedDaughterId >=0 ) 
      mySimEvent->track(fsimi).setClosestDaughterId(theClosestChargedDaughterId);
  }
}

---

Member Data Documentation

const MagneticFieldMap* TrajectoryManager::_theFieldMap [private]

Definition at line 137 of file TrajectoryManager.h.

Referenced by initializeRecoGeometry(), and reconstruct().

const TrackerInteractionGeometry* TrajectoryManager::_theGeometry [private]

Definition at line 136 of file TrajectoryManager.h.

Referenced by initializeLayerMap(), initializeRecoGeometry(), propagateToLayer(), reconstruct(), and theGeometry().

double TrajectoryManager::distCut [private]

Definition at line 142 of file TrajectoryManager.h.

Referenced by TrajectoryManager(), and updateWithDaughters().

bool TrajectoryManager::firstLoop [private]

Definition at line 145 of file TrajectoryManager.h.

Referenced by reconstruct(), and TrajectoryManager().

Pythia6Decays* TrajectoryManager::myDecayEngine [private]

Definition at line 141 of file TrajectoryManager.h.

Referenced by reconstruct(), TrajectoryManager(), updateWithDaughters(), and ~TrajectoryManager().

FSimEvent* TrajectoryManager::mySimEvent [private]

Definition at line 134 of file TrajectoryManager.h.

Referenced by makeSinglePSimHit(), propagateToCalorimeters(), reconstruct(), and updateWithDaughters().

double TrajectoryManager::pTmin [private]

Definition at line 144 of file TrajectoryManager.h.

Referenced by reconstruct(), and TrajectoryManager().

const RandomEngine* TrajectoryManager::random [private]

Definition at line 155 of file TrajectoryManager.h.

Referenced by reconstruct(), and TrajectoryManager().

const GeometricSearchTracker* TrajectoryManager::theGeomSearchTracker [private]

Definition at line 149 of file TrajectoryManager.h.

Referenced by initializeLayerMap(), and initializeRecoGeometry().

const TrackerGeometry* TrajectoryManager::theGeomTracker [private]

Definition at line 148 of file TrajectoryManager.h.

Referenced by initializeTrackerGeometry(), and reconstruct().

std::vector<const DetLayer*> TrajectoryManager::theLayerMap [private]

Definition at line 150 of file TrajectoryManager.h.

Referenced by detLayer(), and initializeLayerMap().

MaterialEffects* TrajectoryManager::theMaterialEffects [private]

Definition at line 139 of file TrajectoryManager.h.

Referenced by createPSimHits(), reconstruct(), TrajectoryManager(), and ~TrajectoryManager().

int TrajectoryManager::theNegLayerOffset [private]

Definition at line 151 of file TrajectoryManager.h.

Referenced by detLayer(), and initializeLayerMap().

std::map<unsigned,std::map<double,PSimHit> > TrajectoryManager::thePSimHits [private]

Definition at line 146 of file TrajectoryManager.h.

Referenced by loadSimHits(), and reconstruct().

---

The documentation for this class was generated from the following files:

• FastSimulation/TrajectoryManager/interface/TrajectoryManager.h
• FastSimulation/TrajectoryManager/src/TrajectoryManager.cc

---