#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. More...

void	initializeRecoGeometry (const GeometricSearchTracker geomSearchTracker, const TrackerInteractionGeometry interactionGeometry, const MagneticFieldMap *aFieldMap)
	Initialize the Reconstruction Geometry. More...

void	initializeTrackerGeometry (const TrackerGeometry *geomTracker)
	Initialize the full Tracker Geometry. More...

void	loadSimHits (edm::PSimHitContainer &c) const

void	propagateToCalorimeters (ParticlePropagator &PP, int fsimi)
	Propagate the particle through the calorimeters. More...

bool	propagateToLayer (ParticlePropagator &PP, unsigned layer)

void	reconstruct (const TrackerTopology *tTopo)
	Does the real job. More...

const TrackerInteractionGeometry *	theGeometry ()
	Returns the pointer to geometry. More...

	TrajectoryManager ()
	Default Constructor. More...

	TrajectoryManager (FSimEvent aSimEvent, const edm::ParameterSet &matEff, const edm::ParameterSet &simHits, const edm::ParameterSet &decays, const RandomEngine engine)
	Constructor from a FSimEvent. More...

	~TrajectoryManager ()
	Default Destructor. More...

Private Member Functions
const DetLayer *	detLayer (const TrackerLayer &layer, float zpos) const
	Returns the DetLayer pointer corresponding to the FAMOS layer. More...

void	initializeLayerMap ()
	Initialize correspondence map between Famos interaction geometry and tracker reco geometry. More...

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 More...

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 More...

TrajectoryStateOnSurface	makeTrajectoryState (const DetLayer layer, const ParticlePropagator &pp, const MagneticField field) const
	Teddy, you must put comments there. More...

void	moveAllDaughters (int fsimi, const Rotation &r, double rescale)
	Move, rescale and rotate all daughters after propagation, material effects and decay of the mother. More...

void	updateWithDaughters (ParticlePropagator &PP, int fsimi)
	Decay the particle and update the SimEvent with daughters. More...

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.

68 {;}

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(), 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 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(), FSimTrack::momentum(), RawParticle::momentum(), moveAllDaughters(), myDecayEngine, mySimEvent, FSimTrack::nDaughters(), PythiaDecays::particleDaughtersPy6(), PythiaDecays::particleDaughtersPy8(), alignCSCRings::r, FSimTrack::setClosestDaughterId(), FSimVertex::setPosition(), mathSSE::sqrt(), FBaseSimEvent::track(), FBaseSimEvent::vertex(), and RawParticle::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);
     }
 
   }
 
 }