#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)
	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 ()
	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)
	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	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

Detailed Description

Definition at line 59 of file TrajectoryManager.h.

Member Typedef Documentation

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

Definition at line 64 of file TrajectoryManager.h.

Constructor & Destructor Documentation

TrajectoryManager::TrajectoryManager ( ) [inline]

Default Constructor.

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

References gather_cfg::cout, decayer, distCut, firstLoop, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), myDecayEngine, pTmin, random, AlCaHLTBitMon_QueryRunRegistry::string, 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") && ( 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("h300",1210,-121.,121.,1210,-121.,121.);
  myHistos->book("h301",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
	)

Create a vector of PSimHits.

Definition at line 540 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;

  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 904 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 798 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 911 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
	)		`[private]`

and there

Definition at line 588 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));
    }
  }
  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 613 of file TrajectoryManager.cc.

References anyDirection, PV3DBase< T, PVType, FrameType >::basicVector(), BoundSurface::bounds(), FSimTrack::closestDaughterId(), gather_cfg::cout, cond::rpcobgas::detid, PXFDetId::disk(), 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(), PV3DBase< T, PVType, FrameType >::mag2(), python::rootplot::argparse::module, FSimTrack::mother(), mySimEvent, getHLTPrescaleColumns::path, TECDetId::petal(), FSimVertex::position(), GloballyPositioned< T >::position(), DetId::rawId(), TIDDetId::ring(), TECDetId::ring(), TIDDetId::side(), SiStripDetId::stereo(), GeomDet::surface(), Bounds::thickness(), Surface::toGlobal(), GeomDet::toLocal(), FBaseSimEvent::track(), TrajectoryStateOnSurface::transverseCurvature(), Vector3DBase< T, FrameTag >::unit(), FSimTrack::vertex(), TIDDetId::wheel(), TECDetId::wheel(), Bounds::width(), hit::x, hit::y, PV3DBase< T, PVType, FrameType >::z(), hit::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::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: 
    {
      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 576 of file TrajectoryManager.cc.

References RawParticle::charge(), pos, 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::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 521 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 366 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 422 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 ( )

Does the real job.

Definition at line 148 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(), and updateWithDaughters().

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 ( ( 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

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

          }
        }

        // 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, 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 451 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);
    }

  }

}