CMSSW_4_4_3_patch1/src/FastSimulation/TrajectoryManager/src/TrajectoryManager.cc

Go to the documentation of this file.

//Framework Headers
#include "FWCore/ParameterSet/interface/ParameterSet.h"

//CMSSW Headers
#include "DataFormats/GeometrySurface/interface/BoundDisk.h"
#include "DataFormats/GeometrySurface/interface/BoundCylinder.h"
#include "DataFormats/GeometrySurface/interface/Surface.h"
#include "DataFormats/GeometrySurface/interface/TangentPlane.h"

// Tracker reco geometry headers 
#include "TrackingTools/DetLayers/interface/DetLayer.h"
#include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
#include "TrackingTools/DetLayers/interface/ForwardDetLayer.h"
#include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
#include "FastSimulation/TrajectoryManager/interface/InsideBoundsMeasurementEstimator.h"
#include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
#include "TrackingTools/GeomPropagators/interface/HelixArbitraryPlaneCrossing.h"
#include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
//#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"

//FAMOS Headers
#include "FastSimulation/TrajectoryManager/interface/TrajectoryManager.h"
#include "FastSimulation/TrajectoryManager/interface/LocalMagneticField.h"
#include "FastSimulation/ParticlePropagator/interface/ParticlePropagator.h"
#include "FastSimulation/TrackerSetup/interface/TrackerInteractionGeometry.h"
#include "FastSimulation/ParticleDecay/interface/Pythia6Decays.h"
#include "FastSimulation/Event/interface/FSimEvent.h"
#include "FastSimulation/Event/interface/FSimVertex.h"
#include "FastSimulation/Event/interface/KineParticleFilter.h"

#include "FastSimulation/Utilities/interface/RandomEngine.h"

//#include "FastSimulation/Utilities/interface/Histos.h"
//#include "FastSimulation/Utilities/interface/FamosLooses.h"
// Numbering scheme

//#define FAMOS_DEBUG
#ifdef FAMOS_DEBUG
#include "DataFormats/SiStripDetId/interface/TIBDetId.h"
#include "DataFormats/SiStripDetId/interface/TIDDetId.h"
#include "DataFormats/SiStripDetId/interface/TOBDetId.h"
#include "DataFormats/SiStripDetId/interface/TECDetId.h"
#include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
#include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
#endif

#include <list>

#include "FWCore/MessageLogger/interface/MessageLogger.h"

TrajectoryManager::TrajectoryManager(FSimEvent* aSimEvent, 
                                     const edm::ParameterSet& matEff,
                                     const edm::ParameterSet& simHits,
                                     const edm::ParameterSet& decays,
                                     const RandomEngine* engine) : 
  mySimEvent(aSimEvent), 
  _theGeometry(0),
  _theFieldMap(0),
  theMaterialEffects(0), 
  myDecayEngine(0), 
  theGeomTracker(0),
  theGeomSearchTracker(0),
  theLayerMap(56, static_cast<const DetLayer*>(0)), // reserve space for layers here
  theNegLayerOffset(27),
  //  myHistos(0),
  random(engine)

{
  
  // Initialize Bthe stable particle decay engine 
  if ( decays.getParameter<bool>("ActivateDecays") ) { 
//    int seed = (int) ( 900000000. * random->flatShoot() );
//    double comE = decays.getParameter<double>("comEnergy");
//    myDecayEngine = new Pythia6Decays(seed,comE);
    myDecayEngine = new Pythia6Decays();
    distCut = decays.getParameter<double>("DistCut");
  }
   // new improvement: Muon brem effects 27-Fev-2011-S.Fonseca UERJ/Brazil
  // 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.);
  */

  
}

void 
TrajectoryManager::initializeRecoGeometry(const GeometricSearchTracker* geomSearchTracker,
                                          const TrackerInteractionGeometry* interactionGeometry,
                                          const MagneticFieldMap* aFieldMap)
{
  
  // 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) { 
  
  theGeomTracker = geomTracker;

}

const TrackerInteractionGeometry*
TrajectoryManager::theGeometry() {
  return _theGeometry;
}

TrajectoryManager::~TrajectoryManager() {

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

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

}

void
TrajectoryManager::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();

}

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

  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) {

  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::updateWithDaughters(ParticlePropagator& PP, int fsimi) {


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

  }

}


void
TrajectoryManager::moveAllDaughters(int fsimi, const Rotation& r, double rescale) { 

  //
  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::createPSimHits(const TrackerLayer& layer,
                                  const ParticlePropagator& PP,
                                  std::map<double,PSimHit>& theHitMap,
                                  int trackID, int partID) {

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

TrajectoryStateOnSurface 
TrajectoryManager::makeTrajectoryState( const DetLayer* layer, 
                                        const ParticlePropagator& pp,
                                        const MagneticField* field) const
{
  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::makePSimHits( const GeomDet* det, 
                                 const TrajectoryStateOnSurface& ts,
                                 std::map<double,PSimHit>& theHitMap,
                                 int tkID, float el, float thick, int pID ) 
{

  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
{

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

}

void 
TrajectoryManager::initializeLayerMap()
{

// 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;
      }
    }
  }  

}

const DetLayer*  
TrajectoryManager::detLayer( const TrackerLayer& layer, float zpos) const
{
  if (zpos > 0 || !layer.forward() ) return theLayerMap[layer.layerNumber()];
  else return theLayerMap[layer.layerNumber()+theNegLayerOffset];
}

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

  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); 
    }
  }

}