#include <TrajectoryManager.h>

Public Types
typedef ROOT::Math::AxisAngle	Rotation
Public Member Functions
void	createPSimHits (const TrackerLayer &layer, const ParticlePropagator &P_before, std::map< double, PSimHit > &theHitMap, int trackID, int partID, const TrackerTopology *tTopo)
	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)
void	reconstruct (const TrackerTopology *tTopo)
	Does the real job.
const TrackerInteractionGeometry *	theGeometry ()
	Returns the pointer to geometry.
	TrajectoryManager ()
	Default Constructor.
	TrajectoryManager (FSimEvent aSimEvent, const edm::ParameterSet &matEff, const edm::ParameterSet &simHits, const edm::ParameterSet &decays, const RandomEngine engine)
	Constructor from a FSimEvent.
	~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, const TrackerTopology tTopo)
	and there
std::pair< double, PSimHit >	makeSinglePSimHit (const GeomDetUnit &det, const TrajectoryStateOnSurface &ts, int tkID, float el, float thick, int pID, const TrackerTopology *tTopo) const
	and there
TrajectoryStateOnSurface	makeTrajectoryState (const DetLayer layer, const ParticlePropagator &pp, const MagneticField field) const
	Teddy, you must put comments there.
void	moveAllDaughters (int fsimi, const Rotation &r, double rescale)
	Move, rescale and rotate all daughters after propagation, material effects and decay of the mother.
void	updateWithDaughters (ParticlePropagator &PP, int fsimi)
	Decay the particle and update the SimEvent with daughters.
Private Attributes
const MagneticFieldMap *	_theFieldMap
const TrackerInteractionGeometry *	_theGeometry
std::string	decayer
double	distCut
bool	firstLoop
PythiaDecays *	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
bool	use_hardcoded

Detailed Description

Definition at line 60 of file TrajectoryManager.h.

Member Typedef Documentation

typedef ROOT::Math::AxisAngle TrajectoryManager::Rotation

Definition at line 65 of file TrajectoryManager.h.

Constructor & Destructor Documentation

TrajectoryManager::TrajectoryManager ( ) [inline]

Default Constructor.

Definition at line 68 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 47 of file TrajectoryManager.cc.

References gather_cfg::cout, decayer, distCut, firstLoop, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), myDecayEngine, pTmin, random, AlCaHLTBitMon_QueryRunRegistry::string, theMaterialEffects, and use_hardcoded.

                                                                 : 
  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),
  use_hardcoded(1)

{  
  //std::cout << "TrajectoryManager.cc 1 use_hardcoded = " << use_hardcoded << std::endl;
  use_hardcoded = matEff.getParameter<bool>("use_hardcoded_geometry");

  // Initialize Bthe stable particle decay engine 
  if ( decays.getParameter<bool>("ActivateDecays") && ( decays.getParameter<std::string>("Decayer") == "pythia6" || decays.getParameter<std::string>("Decayer") == "pythia8" ) ) { 
    decayer = decays.getParameter<std::string>("Decayer");
    myDecayEngine = new PythiaDecays(decayer);
    distCut = decays.getParameter<double>("DistCut");
  } else if (! ( decays.getParameter<std::string>("Decayer") == "pythia6" || decays.getParameter<std::string>("Decayer") == "pythia8" ) )
    std::cout << "No valid decayer has been selected! No decay performed..." << std::endl;
  // Initialize the Material Effects updator, if needed
  if ( matEff.getParameter<bool>("PairProduction") || 
       matEff.getParameter<bool>("Bremsstrahlung") ||
       matEff.getParameter<bool>("MuonBremsstrahlung") ||
       matEff.getParameter<bool>("EnergyLoss") || 
       matEff.getParameter<bool>("MultipleScattering") || 
       matEff.getParameter<bool>("NuclearInteraction")
       )
       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("h302",1210,-121.,121.,1210,-121.,121.);
  myHistos->book("h300",1210,-121.,121.,1210,-121.,121.);
  myHistos->book("h301",1200,-300.,300.,1210,-121.,121.);  
  myHistos->book("h303",1200,-300.,300.,1210,-121.,121.);
  */
}

TrajectoryManager::~TrajectoryManager ( )

Default Destructor.

Definition at line 136 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,
		const TrackerTopology *	tTopo
	)

Create a vector of PSimHits.

Definition at line 554 of file TrajectoryManager.cc.

References AnalyticalPropagator_cfi::AnalyticalPropagator, 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;

//   std::cout << "PP.X() = " << PP.X() << std::endl;
//   std::cout << "PP.Y() = " << PP.Y() << std::endl;
//   std::cout << "PP.Z() = " << PP.Z() << std::endl;
  
  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,tTopo);
  }

}

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

Returns the DetLayer pointer corresponding to the FAMOS layer.

Definition at line 928 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 820 of file TrajectoryManager.cc.

References _theGeometry, GeometricSearchTracker::barrelLayers(), TrackerInteractionGeometry::cylinderBegin(), TrackerInteractionGeometry::cylinderEnd(), newFWLiteAna::found, i, LogDebug, GeometricSearchTracker::negForwardLayers(), GeometricSearchTracker::posForwardLayers(), 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("FastTracking") << "Barrel DetLayer dump: ";
  for (std::vector< BarrelDetLayer*>::const_iterator bl=barrelLayers.begin();
       bl != barrelLayers.end(); ++bl) {
    LogDebug("FastTracking")<< "radius " << (**bl).specificSurface().radius(); 
  }

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

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

  LogDebug("FastTracking")<< "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("FastTracking") << "Famos Layer no " << i->layerNumber()
                            << " is sensitive? " << i->sensitive()
                            << " pos " << i->surface().position();
    if (!i->sensitive()) continue;

    if (cyl != 0) {
      LogDebug("FastTracking") << " 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("FastTracking")<< "Corresponding DetLayer found with radius "
                                 << (**bl).specificSurface().radius();
          break;
        }
      }
      if (!found) {
        edm::LogWarning("FastTracking") << " Trajectory manager FAILED to find a corresponding DetLayer!";
      }
    }
    else {
      LogDebug("FastTracking") << " 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("FastTracking") << "Corresponding DetLayer found with Z pos "
                                  << (**fl).surface().position().z()
                                  << " and radii " 
                                  << (**fl).specificSurface().innerRadius() 
                                  << ", " 
                                  << (**fl).specificSurface().outerRadius(); 
          break;
        }
      }
      if (!found) {
        edm::LogWarning("FastTracking") << "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 106 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 125 of file TrajectoryManager.cc.

References theGeomTracker.

Referenced by FamosManager::setupGeometryAndField().

                                                                               { 
  
  theGeomTracker = geomTracker;

}

void TrajectoryManager::loadSimHits ( edm::PSimHitContainer & c ) const

Definition at line 935 of file TrajectoryManager.cc.

References begin, end, and thePSimHits.

{

  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,
		const TrackerTopology *	tTopo
	)		`[private]`

and there

Definition at line 607 of file TrajectoryManager.cc.

References AlCaHLTBitMon_QueryRunRegistry::comp, 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,tTopo));
    }
  }
  else {
    const GeomDetUnit* du = dynamic_cast<const GeomDetUnit*>(det);
    if (du != 0)
      theHitMap.insert(theHitMap.end(),makeSinglePSimHit( *du, ts, tkID, el, thick, pID,tTopo));
  }

}

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

and there

Definition at line 632 of file TrajectoryManager.cc.

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::LogWarning("FastTracking") << "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: 
    {
      
      theLayer = tTopo->pxbLayer(detid);
      std::cout << "\tPixel Barrel Layer " << theLayer << std::endl;
      stereo = 1;
      break;
    }
  case 2: 
    {
      
      theLayer = tTopo->pxfDisk(detid);
      std::cout << "\tPixel Forward Disk " << theLayer << std::endl;
      stereo = 1;
      break;
    }
  case 3:
    {
      
      theLayer  = tTopo->tibLayer(detid);
      std::cout << "\tTIB Layer " << theLayer << std::endl;
      stereo = module.stereo();
      break;
    }
  case 4:
    {
      
      theLayer = tTopo->tidWheel(detid);
      theRing  = tTopo->tidRing(detid);
      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:
    {
      
      theLayer  = tTopo->tobLayer(detid);
      stereo = tTopo->tobStereo(detid);
      std::cout << "\tTOB Layer " << theLayer << std::endl;
      break;
    }
  case 6:
    {
      
      theLayer = tTopo->tecWheel(detid);
      theRing  = tTopo->tecRing(detid);
      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()));
//      std::cout << "gpos.x() = " << gpos.x() << std::endl;
//      std::cout << "gpos.y() = " << gpos.y() << std::endl;

     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 595 of file TrajectoryManager.cc.

References RawParticle::charge(), pos, GeometricSearchDet::surface(), 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::moveAllDaughters	(	int	fsimi,
		const Rotation &	r,
		double	rescale
	)		`[private]`

Move, rescale and rotate all daughters after propagation, material effects and decay of the mother.

Definition at line 535 of file TrajectoryManager.cc.

References FSimTrack::daughter(), FSimTrack::id(), FSimTrack::momentum(), mySimEvent, FSimTrack::nDaughters(), FSimTrack::setMomentum(), mathSSE::sqrt(), and FBaseSimEvent::track().

Referenced by updateWithDaughters().

                                                                                { 

  //
  for ( unsigned idaugh=0; idaugh < (unsigned)(mySimEvent->track(fsimi).nDaughters()); ++idaugh) { 
    // Initial momentum of the daughter
    XYZTLorentzVector daughMomentum (mySimEvent->track(fsimi).daughter(idaugh).momentum()); 
    // Rotate and rescale
    XYZVector newMomentum (r * daughMomentum.Vect()); 
    newMomentum *= rescale;
    double newEnergy = std::sqrt(newMomentum.mag2() + daughMomentum.mag2());
    // Set the new momentum
    mySimEvent->track(fsimi).setMomentum(XYZTLorentzVector(newMomentum.X(),newMomentum.Y(),newMomentum.Z(),newEnergy));
    // Watch out : recursive call to get all grand-daughters
    int fsimDaug = mySimEvent->track(fsimi).daughter(idaugh).id();
    moveAllDaughters(fsimDaug,r,rescale);
  }
}

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

Propagate the particle through the calorimeters.

Definition at line 381 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 437 of file TrajectoryManager.cc.

References _theGeometry, TrackerInteractionGeometry::cylinderBegin(), TrackerInteractionGeometry::cylinderEnd(), run_regression::done, 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 ( const TrackerTopology * tTopo )

Does the real job.

Definition at line 149 of file TrajectoryManager.cc.

References _theFieldMap, _theGeometry, BaseRawParticleFilter::accept(), FBaseSimEvent::addSimVertex(), createPSimHits(), TrackerInteractionGeometry::cylinderBegin(), TrackerInteractionGeometry::cylinderEnd(), FSimVertexType::END_VERTEX, FBaseSimEvent::filter(), firstLoop, MaterialEffects::interact(), python::cmstools::loop(), mySimEvent, FSimTrack::nDaughters(), FSimTrack::notYetToEndVertex(), FSimEvent::nTracks(), propagateToCalorimeters(), pTmin, random, MaterialEffects::save(), FSimTrack::setPropagate(), summarizeEdmComparisonLogfiles::success, theGeomTracker, theMaterialEffects, thePSimHits, FBaseSimEvent::track(), CoreSimTrack::type(), updateWithDaughters(), and use_hardcoded.

Referenced by FamosManager::reconstruct().

{

  // Clear the hits of the previous event
  //  thePSimHits->clear();
  thePSimHits.clear();

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

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

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

  // 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 ( debug ) std::cout << mySimEvent->track(fsimi) << std::endl;
    //if ( debug ) std::cout << "Not yet at end vertex ? " << mySimEvent->track(fsimi).notYetToEndVertex(myBeamPipe) << std::endl;
    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(use_hardcoded){
      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();
      } 
    }
    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

      // Skip layers with no material (kept just for historical reasons)
      if ( cyliter->surface().mediumProperties().radLen() < 1E-10 ) { 
        ++cyliter; ++cyl;
        continue;
      }
      
      // 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() || (!mySimEvent->track(fsimi).nDaughters() && PP.PDGcTau()<1E-3 ) ) { 
        updateWithDaughters(PP,fsimi);
        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); 

        // There is a PP.setXYZT=(0,0,0,0) if bremss fails
        saveHit &= PP.E()>1E-6;

        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(), tTopo);

            /*
            myHistos->fill("h302",PP.X() ,PP.Y());
            if ( sin(PP.vertex().Phi()) > 0. ) 
              myHistos->fill("h303",PP.Z(),PP.R());
            else
              myHistos->fill("h303",PP.Z(),-PP.R());
            */

          }
        }

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

        //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, FSimVertexType::END_VERTEX);
          
      }

      // 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();

}

const TrackerInteractionGeometry * TrajectoryManager::theGeometry ( )

Returns the pointer to geometry.

Definition at line 132 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 466 of file TrajectoryManager.cc.

References FBaseSimEvent::addSimTrack(), FBaseSimEvent::addSimVertex(), angle(), RawParticle::charge(), FSimVertexType::DECAY_VERTEX, decayer, distCut, FSimTrack::endVertex(), FSimVertex::id(), RawParticle::momentum(), FSimTrack::momentum(), moveAllDaughters(), myDecayEngine, mySimEvent, FSimTrack::nDaughters(), PythiaDecays::particleDaughtersPy6(), PythiaDecays::particleDaughtersPy8(), alignCSCRings::r, FSimTrack::setClosestDaughterId(), FSimVertex::setPosition(), mathSSE::sqrt(), FBaseSimEvent::track(), RawParticle::vertex(), and FBaseSimEvent::vertex().

Referenced by propagateToCalorimeters(), and reconstruct().

                                                                        {

  // The particle was already decayed in the GenEvent, but still the particle was 
  // allowed to propagate (for magnetic field bending, for material effects, etc...)
  // Just modify the momentum of the daughters in that case 
  unsigned nDaugh = mySimEvent->track(fsimi).nDaughters();
  if ( nDaugh ) {

    // Move the vertex
    unsigned vertexId = mySimEvent->track(fsimi).endVertex().id();
    mySimEvent->vertex(vertexId).setPosition(PP.vertex());

    // Before-propagation and after-propagation momentum and vertex position
    XYZTLorentzVector momentumBefore = mySimEvent->track(fsimi).momentum();
    XYZTLorentzVector momentumAfter = PP.momentum();
    double magBefore = std::sqrt(momentumBefore.Vect().mag2());
    double magAfter = std::sqrt(momentumAfter.Vect().mag2());
    // Rotation to be applied
    XYZVector axis = momentumBefore.Vect().Cross(momentumAfter.Vect());
    double angle = std::acos(momentumBefore.Vect().Dot(momentumAfter.Vect())/(magAfter*magBefore));
    Rotation r(axis,angle);
    // Rescaling to be applied
    double rescale = magAfter/magBefore;

    // Move, rescale and rotate daugthers, grand-daughters, etc. 
    moveAllDaughters(fsimi,r,rescale);

  // The particle is not decayed in the GenEvent, decay it with PYTHIA 
  } else { 

    // Decays are not activated : do nothing
    if ( !myDecayEngine ) return;
    
    // Invoke PYDECY (Pythia6) or Pythia8 to decay the particle and get the daughters
    const DaughterParticleList& daughters = (decayer == "pythia6") ? myDecayEngine->particleDaughtersPy6(PP) : myDecayEngine->particleDaughtersPy8(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, 
                                             FSimVertexType::DECAY_VERTEX);
      
      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);
    }

  }

}