Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:11:25 

0001 //Framework Headers
0002 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0003 
0004 //CMSSW Headers
0005 #include "DataFormats/GeometrySurface/interface/BoundDisk.h"
0006 #include "DataFormats/GeometrySurface/interface/BoundCylinder.h"
0007 #include "DataFormats/GeometrySurface/interface/Surface.h"
0008 #include "DataFormats/GeometrySurface/interface/TangentPlane.h"
0009 
0010 // Tracker reco geometry headers
0011 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0012 #include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
0013 #include "TrackingTools/DetLayers/interface/ForwardDetLayer.h"
0014 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0015 #include "FastSimulation/TrajectoryManager/interface/InsideBoundsMeasurementEstimator.h"
0016 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0017 #include "TrackingTools/GeomPropagators/interface/HelixArbitraryPlaneCrossing.h"
0018 #include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
0019 //#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0020 
0021 //FAMOS Headers
0022 #include "FastSimulation/TrajectoryManager/interface/TrajectoryManager.h"
0023 #include "FastSimulation/TrajectoryManager/interface/LocalMagneticField.h"
0024 #include "FastSimulation/ParticlePropagator/interface/ParticlePropagator.h"
0025 #include "FastSimulation/TrackerSetup/interface/TrackerInteractionGeometry.h"
0026 #include "FastSimulation/ParticleDecay/interface/PythiaDecays.h"
0027 #include "FastSimulation/Event/interface/FSimEvent.h"
0028 #include "FastSimulation/Event/interface/FSimVertex.h"
0029 #include "FastSimulation/Event/interface/KineParticleFilter.h"
0030 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0031 #include "FastSimulation/Particle/interface/pdg_functions.h"
0032 
0033 //#include "FastSimulation/Utilities/interface/Histos.h"
0034 //#include "FastSimulation/Utilities/interface/FamosLooses.h"
0035 // Numbering scheme
0036 
0037 //#define FAMOS_DEBUG
0038 #ifdef FAMOS_DEBUG
0039 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0040 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0041 #endif
0042 
0043 #include <list>
0044 
0045 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0046 
0047 TrajectoryManager::TrajectoryManager(FSimEvent* aSimEvent,
0048                                      const edm::ParameterSet& matEff,
0049                                      const edm::ParameterSet& simHits,
0050                                      const edm::ParameterSet& decays)
0051     : mySimEvent(aSimEvent),
0052       _theGeometry(nullptr),
0053       _theFieldMap(nullptr),
0054       theMaterialEffects(nullptr),
0055       myDecayEngine(nullptr),
0056       theGeomTracker(nullptr),
0057       theGeomSearchTracker(nullptr),
0058       theLayerMap(56, static_cast<const DetLayer*>(nullptr)),  // reserve space for layers here
0059       theNegLayerOffset(27),
0060       //  myHistos(0),
0061       use_hardcoded(true)
0062 
0063 {
0064   //std::cout << "TrajectoryManager.cc 1 use_hardcoded = " << use_hardcoded << std::endl;
0065   use_hardcoded = matEff.getParameter<bool>("use_hardcoded_geometry");
0066 
0067   // Initialize Bthe stable particle decay engine
0068   if (decays.getParameter<bool>("ActivateDecays")) {
0069     myDecayEngine = new PythiaDecays();
0070     distCut = decays.getParameter<double>("DistCut");
0071   }
0072   // Initialize the Material Effects updator, if needed
0073   if (matEff.getParameter<bool>("PairProduction") || matEff.getParameter<bool>("Bremsstrahlung") ||
0074       matEff.getParameter<bool>("MuonBremsstrahlung") || matEff.getParameter<bool>("EnergyLoss") ||
0075       matEff.getParameter<bool>("MultipleScattering") || matEff.getParameter<bool>("NuclearInteraction"))
0076     theMaterialEffects = new MaterialEffects(matEff);
0077 
0078   // Save SimHits according to Optiom
0079   // Only the hits from first half loop is saved
0080   firstLoop = simHits.getUntrackedParameter<bool>("firstLoop", true);
0081   // Only if pT>pTmin are the hits saved
0082   pTmin = simHits.getUntrackedParameter<double>("pTmin", 0.5);
0083 
0084   /*
0085   // Get the Famos Histos pointer
0086   myHistos = Histos::instance();
0087 
0088   // Initialize a few histograms
0089    
0090   myHistos->book("h302",1210,-121.,121.,1210,-121.,121.);
0091   myHistos->book("h300",1210,-121.,121.,1210,-121.,121.);
0092   myHistos->book("h301",1200,-300.,300.,1210,-121.,121.);  
0093   myHistos->book("h303",1200,-300.,300.,1210,-121.,121.);
0094   */
0095 }
0096 
0097 void TrajectoryManager::initializeRecoGeometry(const GeometricSearchTracker* geomSearchTracker,
0098                                                const TrackerInteractionGeometry* interactionGeometry,
0099                                                const MagneticFieldMap* aFieldMap) {
0100   // Initialize the reco tracker geometry
0101   theGeomSearchTracker = geomSearchTracker;
0102 
0103   // Initialize the simplified tracker geometry
0104   _theGeometry = interactionGeometry;
0105 
0106   initializeLayerMap();
0107 
0108   // Initialize the magnetic field
0109   _theFieldMap = aFieldMap;
0110 }
0111 
0112 void TrajectoryManager::initializeTrackerGeometry(const TrackerGeometry* geomTracker) { theGeomTracker = geomTracker; }
0113 
0114 const TrackerInteractionGeometry* TrajectoryManager::theGeometry() { return _theGeometry; }
0115 
0116 TrajectoryManager::~TrajectoryManager() {
0117   if (myDecayEngine)
0118     delete myDecayEngine;
0119   if (theMaterialEffects)
0120     delete theMaterialEffects;
0121 
0122   //Write the histograms
0123   /*
0124     myHistos->put("histos.root");
0125     if ( myHistos ) delete myHistos;
0126   */
0127 }
0128 
0129 void TrajectoryManager::reconstruct(const TrackerTopology* tTopo, RandomEngineAndDistribution const* random) {
0130   // Clear the hits of the previous event
0131   //  thePSimHits->clear();
0132   thePSimHits.clear();
0133 
0134   // The new event
0135   XYZTLorentzVector myBeamPipe = XYZTLorentzVector(0., 2.5, 9999999., 0.);
0136 
0137   std::list<TrackerLayer>::const_iterator cyliter;
0138 
0139   // bool debug = mySimEvent->id().event() == 8;
0140 
0141   // Loop over the particles (watch out: increasing upper limit!)
0142   for (int fsimi = 0; fsimi < (int)mySimEvent->nTracks(); ++fsimi) {
0143     // If the particle has decayed inside the beampipe, or decays
0144     // immediately, there is nothing to do
0145     //if ( debug ) std::cout << mySimEvent->track(fsimi) << std::endl;
0146     //if ( debug ) std::cout << "Not yet at end vertex ? " << mySimEvent->track(fsimi).notYetToEndVertex(myBeamPipe) << std::endl;
0147     if (!mySimEvent->track(fsimi).notYetToEndVertex(myBeamPipe))
0148       continue;
0149     mySimEvent->track(fsimi).setPropagate();
0150 
0151     // Get the geometry elements
0152     cyliter = _theGeometry->cylinderBegin();
0153     // Prepare the propagation
0154     ParticlePropagator PP(mySimEvent->track(fsimi), _theFieldMap, random, mySimEvent->theTable());
0155     //The real work starts here
0156     int success = 1;
0157     int sign = +1;
0158     int loop = 0;
0159     int cyl = 0;
0160 
0161     // Find the initial cylinder to propagate to.
0162     for (; cyliter != _theGeometry->cylinderEnd(); ++cyliter) {
0163       PP.setPropagationConditions(*cyliter);
0164       if (PP.inside() && !PP.onSurface())
0165         break;
0166       ++cyl;
0167     }
0168 
0169     // The particle has a pseudo-rapidity (position or momentum direction)
0170     // in excess of 3.0. Just simply go to the last tracker layer
0171     // without bothering with all the details of the propagation and
0172     // material effects.
0173     // 08/02/06 - pv: increase protection from 0.99 (eta=2.9932) to 0.9998 (eta=4.9517)
0174     //                to simulate material effects at large eta
0175     // if above 0.99: propagate to the last tracker cylinder where the material is concentrated!
0176     double ppcos2T = PP.particle().cos2Theta();
0177     double ppcos2V = PP.particle().cos2ThetaV();
0178 
0179     if (use_hardcoded) {
0180       if ((ppcos2T > 0.99 && ppcos2T < 0.9998) && (cyl == 0 || (ppcos2V > 0.99 && ppcos2V < 0.9998))) {
0181         if (cyliter != _theGeometry->cylinderEnd()) {
0182           cyliter = _theGeometry->cylinderEnd();
0183           --cyliter;
0184         }
0185         // if above 0.9998: don't propagate at all (only to the calorimeters directly)
0186       } else if (ppcos2T > 0.9998 && (cyl == 0 || ppcos2V > 0.9998)) {
0187         cyliter = _theGeometry->cylinderEnd();
0188       }
0189     } else {
0190       if (ppcos2T > 0.9998 && (cyl == 0 || ppcos2V > 0.9998)) {
0191         cyliter = _theGeometry->cylinderEnd();
0192       }
0193     }
0194 
0195     // Loop over the cylinders
0196     while (cyliter != _theGeometry->cylinderEnd() && loop < 100 &&                // No more than 100 loops
0197            mySimEvent->track(fsimi).notYetToEndVertex(PP.particle().vertex())) {  // The particle decayed
0198 
0199       // Skip layers with no material (kept just for historical reasons)
0200       if (cyliter->surface().mediumProperties().radLen() < 1E-10) {
0201         ++cyliter;
0202         ++cyl;
0203         continue;
0204       }
0205 
0206       // Pathological cases:
0207       // To prevent from interacting twice in a row with the same layer
0208       //      bool escapeBarrel    = (PP.getSuccess() == -1 && success == 1);
0209       bool escapeBarrel = PP.getSuccess() == -1;
0210       bool escapeEndcap = (PP.getSuccess() == -2 && success == 1);
0211       // To break the loop
0212       bool fullPropagation = (PP.getSuccess() <= 0 && success == 0) || escapeEndcap;
0213 
0214       if (escapeBarrel) {
0215         ++cyliter;
0216         ++cyl;
0217         while (cyliter != _theGeometry->cylinderEnd() && cyliter->forward()) {
0218           sign = 1;
0219           ++cyliter;
0220           ++cyl;
0221         }
0222 
0223         if (cyliter == _theGeometry->cylinderEnd()) {
0224           --cyliter;
0225           --cyl;
0226           fullPropagation = true;
0227         }
0228       }
0229 
0230       // Define the propagation conditions
0231       PP.setPropagationConditions(*cyliter, !fullPropagation);
0232       if (escapeEndcap)
0233         PP.increaseRCyl(0.0005);
0234 
0235       // Remember last propagation outcome
0236       success = PP.getSuccess();
0237 
0238       // Propagation was not successful :
0239       // Change the sign of the cylinder increment and count the loops
0240       if (!PP.propagateToBoundSurface(*cyliter) || PP.getSuccess() <= 0) {
0241         sign = -sign;
0242         ++loop;
0243       }
0244 
0245       // The particle may have decayed on its way... in which the daughters
0246       // have to be added to the event record
0247       if (PP.hasDecayed() ||
0248           (!mySimEvent->track(fsimi).nDaughters() && pdg::cTau(PP.particle().pid(), mySimEvent->theTable()) < 1E-3)) {
0249         updateWithDaughters(PP, fsimi, random);
0250         break;
0251       }
0252 
0253       // Exit by the endcaps or innermost cylinder :
0254       // Positive cylinder increment
0255       if (PP.getSuccess() == 2 || cyliter == _theGeometry->cylinderBegin())
0256         sign = +1;
0257 
0258       // Successful propagation to a cylinder, with some Material :
0259       if (PP.getSuccess() > 0 && PP.onFiducial()) {
0260         bool saveHit = ((loop == 0 && sign > 0) || !firstLoop) &&  // Save only first half loop
0261                        PP.particle().charge() != 0. &&             // Consider only charged particles
0262                        cyliter->sensitive() &&                     // Consider only sensitive layers
0263                        PP.particle().Perp2() > pTmin * pTmin;      // Consider only pT > pTmin
0264 
0265         // Material effects are simulated there
0266         if (theMaterialEffects)
0267           theMaterialEffects->interact(*mySimEvent, *cyliter, PP, fsimi, random);
0268 
0269         // There is a PP.setXYZT=(0,0,0,0) if bremss fails
0270         saveHit &= PP.particle().E() > 1E-6;
0271 
0272         if (saveHit) {
0273           // Consider only active layers
0274           if (cyliter->sensitive()) {
0275             // Add information to the FSimTrack (not yet available)
0276             //      myTrack.addSimHit(PP.particle(),layer);
0277 
0278             // Return one or two (for overlap regions) PSimHits in the full
0279             // tracker geometry
0280             if (theGeomTracker)
0281               createPSimHits(*cyliter, PP, thePSimHits[fsimi], fsimi, mySimEvent->track(fsimi).type(), tTopo);
0282 
0283             /*
0284           myHistos->fill("h302",PP.particle().X() ,PP.particle().Y());
0285           if ( sin(PP.particle().vertex().Phi()) > 0. ) 
0286           myHistos->fill("h303",PP.particle().Z(),PP.particle().R());
0287           else
0288           myHistos->fill("h303",PP.Z(),-PP.particle().R());
0289         */
0290           }
0291         }
0292 
0293         // Fill Histos (~poor man event display)
0294         /*   
0295          myHistos->fill("h300",PP.particle().x(),PP.particle().y());
0296          if ( sin(PP.particle().vertex().phi()) > 0. ) 
0297          myHistos->fill("h301",PP.particle().z(),sqrt(PP.particle().vertex().Perp2()));
0298          else
0299          myHistos->fill("h301",PP.particle().z(),-sqrt(PP.particle().vertex().Perp2()));
0300     */
0301 
0302         //The particle may have lost its energy in the material
0303         if (mySimEvent->track(fsimi).notYetToEndVertex(PP.particle().vertex()) &&
0304             !mySimEvent->filter().acceptParticle(PP.particle()))
0305           mySimEvent->addSimVertex(PP.particle().vertex(), fsimi, FSimVertexType::END_VERTEX);
0306       }
0307 
0308       // Stop here if the particle has reached an end
0309       if (mySimEvent->track(fsimi).notYetToEndVertex(PP.particle().vertex())) {
0310         // Otherwise increment the cylinder iterator
0311         //  do {
0312         if (sign == 1) {
0313           ++cyliter;
0314           ++cyl;
0315         } else {
0316           --cyliter;
0317           --cyl;
0318         }
0319 
0320         // Check if the last surface has been reached
0321         if (cyliter == _theGeometry->cylinderEnd()) {
0322           // Try to propagate to the ECAL in half a loop
0323           // Note: Layer1 = ECAL Barrel entrance, or Preshower
0324           // entrance, or ECAL Endcap entrance (in the corner)
0325           PP.propagateToEcal();
0326           // PP.propagateToPreshowerLayer1();
0327 
0328           // If it is not possible, try go back to the last cylinder
0329           if (PP.getSuccess() == 0) {
0330             --cyliter;
0331             --cyl;
0332             sign = -sign;
0333             PP.setPropagationConditions(*cyliter);
0334             PP.propagateToBoundSurface(*cyliter);
0335 
0336             // If there is definitely no way, leave it here.
0337             if (PP.getSuccess() < 0) {
0338               ++cyliter;
0339               ++cyl;
0340             }
0341           }
0342 
0343           // Check if the particle has decayed on the way to ECAL
0344           if (PP.hasDecayed())
0345             updateWithDaughters(PP, fsimi, random);
0346         }
0347       }
0348     }
0349 
0350     // Propagate all particles without a end vertex to the Preshower,
0351     // theECAL and the HCAL.
0352     if (mySimEvent->track(fsimi).notYetToEndVertex(PP.particle().vertex()))
0353       propagateToCalorimeters(PP, fsimi, random);
0354   }
0355 
0356   // Save the information from Nuclear Interaction Generation
0357   if (theMaterialEffects)
0358     theMaterialEffects->save();
0359 }
0360 
0361 void TrajectoryManager::propagateToCalorimeters(ParticlePropagator& PP,
0362                                                 int fsimi,
0363                                                 RandomEngineAndDistribution const* random) {
0364   FSimTrack& myTrack = mySimEvent->track(fsimi);
0365 
0366   // Set the position and momentum at the end of the tracker volume
0367   myTrack.setTkPosition(PP.particle().vertex().Vect());
0368   myTrack.setTkMomentum(PP.particle().momentum());
0369 
0370   // Propagate to Preshower Layer 1
0371   PP.propagateToPreshowerLayer1(false);
0372   if (PP.hasDecayed()) {
0373     updateWithDaughters(PP, fsimi, random);
0374     return;
0375   }
0376   if (myTrack.notYetToEndVertex(PP.particle().vertex()) && PP.getSuccess() > 0)
0377     myTrack.setLayer1(PP.particle(), PP.getSuccess());
0378 
0379   // Propagate to Preshower Layer 2
0380   PP.propagateToPreshowerLayer2(false);
0381   if (PP.hasDecayed()) {
0382     updateWithDaughters(PP, fsimi, random);
0383     return;
0384   }
0385   if (myTrack.notYetToEndVertex(PP.particle().vertex()) && PP.getSuccess() > 0)
0386     myTrack.setLayer2(PP.particle(), PP.getSuccess());
0387 
0388   // Propagate to Ecal Endcap
0389   PP.propagateToEcalEntrance(false);
0390   if (PP.hasDecayed()) {
0391     updateWithDaughters(PP, fsimi, random);
0392     return;
0393   }
0394   if (myTrack.notYetToEndVertex(PP.particle().vertex()))
0395     myTrack.setEcal(PP.particle(), PP.getSuccess());
0396 
0397   // Propagate to HCAL entrance
0398   PP.propagateToHcalEntrance(false);
0399   if (PP.hasDecayed()) {
0400     updateWithDaughters(PP, fsimi, random);
0401     return;
0402   }
0403   if (myTrack.notYetToEndVertex(PP.particle().vertex()))
0404     myTrack.setHcal(PP.particle(), PP.getSuccess());
0405 
0406   // Propagate to VFCAL entrance
0407   PP.propagateToVFcalEntrance(false);
0408   if (PP.hasDecayed()) {
0409     updateWithDaughters(PP, fsimi, random);
0410     return;
0411   }
0412   if (myTrack.notYetToEndVertex(PP.particle().vertex()))
0413     myTrack.setVFcal(PP.particle(), PP.getSuccess());
0414 }
0415 
0416 bool TrajectoryManager::propagateToLayer(ParticlePropagator& PP, unsigned layer) {
0417   std::list<TrackerLayer>::const_iterator cyliter;
0418   bool done = false;
0419 
0420   // Get the geometry elements
0421   cyliter = _theGeometry->cylinderBegin();
0422 
0423   // Find the layer to propagate to.
0424   for (; cyliter != _theGeometry->cylinderEnd(); ++cyliter) {
0425     if (layer != cyliter->layerNumber())
0426       continue;
0427 
0428     PP.setPropagationConditions(*cyliter);
0429 
0430     done = PP.propagateToBoundSurface(*cyliter) && PP.getSuccess() > 0 && PP.onFiducial();
0431 
0432     break;
0433   }
0434 
0435   return done;
0436 }
0437 
0438 void TrajectoryManager::updateWithDaughters(ParticlePropagator& PP,
0439                                             int fsimi,
0440                                             RandomEngineAndDistribution const* random) {
0441   // The particle was already decayed in the GenEvent, but still the particle was
0442   // allowed to propagate (for magnetic field bending, for material effects, etc...)
0443   // Just modify the momentum of the daughters in that case
0444   unsigned nDaugh = mySimEvent->track(fsimi).nDaughters();
0445   if (nDaugh) {
0446     // Move the vertex
0447     unsigned vertexId = mySimEvent->track(fsimi).endVertex().id();
0448     mySimEvent->vertex(vertexId).setPosition(PP.particle().vertex());
0449 
0450     // Before-propagation and after-propagation momentum and vertex position
0451     XYZTLorentzVector momentumBefore = mySimEvent->track(fsimi).momentum();
0452     const XYZTLorentzVector& momentumAfter = PP.particle().momentum();
0453     double magBefore = std::sqrt(momentumBefore.Vect().mag2());
0454     double magAfter = std::sqrt(momentumAfter.Vect().mag2());
0455     // Rotation to be applied
0456     XYZVector axis = momentumBefore.Vect().Cross(momentumAfter.Vect());
0457     double angle = std::acos(momentumBefore.Vect().Dot(momentumAfter.Vect()) / (magAfter * magBefore));
0458     Rotation r(axis, angle);
0459     // Rescaling to be applied
0460     double rescale = magAfter / magBefore;
0461 
0462     // Move, rescale and rotate daugthers, grand-daughters, etc.
0463     moveAllDaughters(fsimi, r, rescale);
0464 
0465     // The particle is not decayed in the GenEvent, decay it with PYTHIA
0466   } else {
0467     // Decays are not activated : do nothing
0468     if (!myDecayEngine)
0469       return;
0470 
0471     // Invoke PYDECY (Pythia6) or Pythia8 to decay the particle and get the daughters
0472     const DaughterParticleList& daughters = myDecayEngine->particleDaughters(PP, &random->theEngine());
0473 
0474     // Update the FSimEvent with an end vertex and with the daughters
0475     if (!daughters.empty()) {
0476       double distMin = 1E99;
0477       int theClosestChargedDaughterId = -1;
0478       DaughterParticleIterator daughter = daughters.begin();
0479 
0480       int ivertex = mySimEvent->addSimVertex(daughter->vertex(), fsimi, FSimVertexType::DECAY_VERTEX);
0481 
0482       if (ivertex != -1) {
0483         for (; daughter != daughters.end(); ++daughter) {
0484           int theDaughterId = mySimEvent->addSimTrack(&(*daughter), ivertex);
0485           // Find the closest charged daughter (if charged mother)
0486           if (PP.particle().charge() * daughter->charge() > 1E-10) {
0487             double dist = (daughter->Vect().Unit().Cross(PP.particle().Vect().Unit())).R();
0488             if (dist < distCut && dist < distMin) {
0489               distMin = dist;
0490               theClosestChargedDaughterId = theDaughterId;
0491             }
0492           }
0493         }
0494       }
0495       // Attach mother and closest daughter sp as to cheat tracking ;-)
0496       if (theClosestChargedDaughterId >= 0)
0497         mySimEvent->track(fsimi).setClosestDaughterId(theClosestChargedDaughterId);
0498     }
0499   }
0500 }
0501 
0502 void TrajectoryManager::moveAllDaughters(int fsimi, const Rotation& r, double rescale) {
0503   //
0504   for (unsigned idaugh = 0; idaugh < (unsigned)(mySimEvent->track(fsimi).nDaughters()); ++idaugh) {
0505     // Initial momentum of the daughter
0506     XYZTLorentzVector daughMomentum(mySimEvent->track(fsimi).daughter(idaugh).momentum());
0507     // Rotate and rescale
0508     XYZVector newMomentum(r * daughMomentum.Vect());
0509     newMomentum *= rescale;
0510     double newEnergy = std::sqrt(newMomentum.mag2() + daughMomentum.mag2());
0511     // Set the new momentum
0512     mySimEvent->track(fsimi).setMomentum(
0513         XYZTLorentzVector(newMomentum.X(), newMomentum.Y(), newMomentum.Z(), newEnergy));
0514     // Watch out : recursive call to get all grand-daughters
0515     int fsimDaug = mySimEvent->track(fsimi).daughter(idaugh).id();
0516     moveAllDaughters(fsimDaug, r, rescale);
0517   }
0518 }
0519 
0520 void TrajectoryManager::createPSimHits(const TrackerLayer& layer,
0521                                        const ParticlePropagator& PP,
0522                                        std::map<double, PSimHit>& theHitMap,
0523                                        int trackID,
0524                                        int partID,
0525                                        const TrackerTopology* tTopo) {
0526   // Propagate the particle coordinates to the closest tracker detector(s)
0527   // in this layer and create the PSimHit(s)
0528 
0529   //  const MagneticField& mf = MagneticFieldMap::instance()->magneticField();
0530   // This solution is actually much faster !
0531   LocalMagneticField mf(PP.getMagneticField());
0532   AnalyticalPropagator alongProp(&mf, anyDirection);
0533   InsideBoundsMeasurementEstimator est;
0534 
0535   //   std::cout << "PP.X() = " << PP.X() << std::endl;
0536   //   std::cout << "PP.Y() = " << PP.Y() << std::endl;
0537   //   std::cout << "PP.Z() = " << PP.Z() << std::endl;
0538 
0539   typedef GeometricSearchDet::DetWithState DetWithState;
0540   const DetLayer* tkLayer = detLayer(layer, PP.particle().Z());
0541 
0542   TrajectoryStateOnSurface trajState = makeTrajectoryState(tkLayer, PP, &mf);
0543   float thickness = theMaterialEffects ? theMaterialEffects->thickness() : 0.;
0544   float eloss = theMaterialEffects ? theMaterialEffects->energyLoss() : 0.;
0545 
0546   // Find, in the corresponding layers, the detectors compatible
0547   // with the current track
0548   std::vector<DetWithState> compat = tkLayer->compatibleDets(trajState, alongProp, est);
0549 
0550   // And create the corresponding PSimHits
0551   std::map<double, PSimHit> theTrackHits;
0552   for (std::vector<DetWithState>::const_iterator i = compat.begin(); i != compat.end(); i++) {
0553     // Correct Eloss for last 3 rings of TEC (thick sensors, 0.05 cm)
0554     // Disgusting fudge factor !
0555     makePSimHits(i->first, i->second, theHitMap, trackID, eloss, thickness, partID, tTopo);
0556   }
0557 }
0558 
0559 TrajectoryStateOnSurface TrajectoryManager::makeTrajectoryState(const DetLayer* layer,
0560                                                                 const ParticlePropagator& pp,
0561                                                                 const MagneticField* field) const {
0562   GlobalPoint pos(pp.particle().X(), pp.particle().Y(), pp.particle().Z());
0563   GlobalVector mom(pp.particle().Px(), pp.particle().Py(), pp.particle().Pz());
0564   auto plane = layer->surface().tangentPlane(pos);
0565   return TrajectoryStateOnSurface(GlobalTrajectoryParameters(pos, mom, TrackCharge(pp.particle().charge()), field),
0566                                   *plane);
0567 }
0568 
0569 void TrajectoryManager::makePSimHits(const GeomDet* det,
0570                                      const TrajectoryStateOnSurface& ts,
0571                                      std::map<double, PSimHit>& theHitMap,
0572                                      int tkID,
0573                                      float el,
0574                                      float thick,
0575                                      int pID,
0576                                      const TrackerTopology* tTopo) {
0577   std::vector<const GeomDet*> comp = det->components();
0578   if (!comp.empty()) {
0579     for (std::vector<const GeomDet*>::const_iterator i = comp.begin(); i != comp.end(); i++) {
0580       auto du = (*i);
0581       if (du->isLeaf())  // not even needed (or it should iterate if really not leaf)
0582         theHitMap.insert(theHitMap.end(), makeSinglePSimHit(*du, ts, tkID, el, thick, pID, tTopo));
0583     }
0584   } else {
0585     auto du = (det);
0586     theHitMap.insert(theHitMap.end(), makeSinglePSimHit(*du, ts, tkID, el, thick, pID, tTopo));
0587   }
0588 }
0589 
0590 std::pair<double, PSimHit> TrajectoryManager::makeSinglePSimHit(const GeomDetUnit& det,
0591                                                                 const TrajectoryStateOnSurface& ts,
0592                                                                 int tkID,
0593                                                                 float el,
0594                                                                 float thick,
0595                                                                 int pID,
0596                                                                 const TrackerTopology* tTopo) const {
0597   const float onSurfaceTolarance = 0.01;  // 10 microns
0598 
0599   LocalPoint lpos;
0600   LocalVector lmom;
0601   if (fabs(det.toLocal(ts.globalPosition()).z()) < onSurfaceTolarance) {
0602     lpos = ts.localPosition();
0603     lmom = ts.localMomentum();
0604   } else {
0605     HelixArbitraryPlaneCrossing crossing(
0606         ts.globalPosition().basicVector(), ts.globalMomentum().basicVector(), ts.transverseCurvature(), anyDirection);
0607     std::pair<bool, double> path = crossing.pathLength(det.surface());
0608     if (!path.first) {
0609       // edm::LogWarning("FastTracking") << "TrajectoryManager ERROR: crossing with det failed, skipping PSimHit";
0610       return std::pair<double, PSimHit>(0., PSimHit());
0611     }
0612     lpos = det.toLocal(GlobalPoint(crossing.position(path.second)));
0613     lmom = det.toLocal(GlobalVector(crossing.direction(path.second)));
0614     lmom = lmom.unit() * ts.localMomentum().mag();
0615   }
0616 
0617   // The module (half) thickness
0618   const BoundPlane& theDetPlane = det.surface();
0619   float halfThick = 0.5 * theDetPlane.bounds().thickness();
0620   // The Energy loss rescaled to the module thickness
0621   float eloss = el;
0622   if (thick > 0.) {
0623     // Total thickness is in radiation lengths, 1 radlen = 9.36 cm
0624     // Sensitive module thickness is about 30 microns larger than
0625     // the module thickness itself
0626     eloss *= (2. * halfThick - 0.003) / (9.36 * thick);
0627   }
0628   // The entry and exit points, and the time of flight
0629   float pZ = lmom.z();
0630   LocalPoint entry = lpos + (-halfThick / pZ) * lmom;
0631   LocalPoint exit = lpos + halfThick / pZ * lmom;
0632   float tof = ts.globalPosition().mag() / 30.;  // in nanoseconds, FIXME: very approximate
0633 
0634   // If a hadron suffered a nuclear interaction, just assign the hits of the closest
0635   // daughter to the mother's track. The same applies to a charged particle decay into
0636   // another charged particle.
0637   int localTkID = tkID;
0638   if (!mySimEvent->track(tkID).noMother() && mySimEvent->track(tkID).mother().closestDaughterId() == tkID)
0639     localTkID = mySimEvent->track(tkID).mother().id();
0640 
0641   // FIXME: fix the track ID and the particle ID
0642   PSimHit hit(
0643       entry, exit, lmom.mag(), tof, eloss, pID, det.geographicalId().rawId(), localTkID, lmom.theta(), lmom.phi());
0644 
0645   // Check that the PSimHit is physically on the module!
0646   unsigned subdet = DetId(hit.detUnitId()).subdetId();
0647   double boundX = theDetPlane.bounds().width() / 2.;
0648   double boundY = theDetPlane.bounds().length() / 2.;
0649 
0650   // Special treatment for TID and TEC trapeziodal modules
0651   if (subdet == 4 || subdet == 6)
0652     boundX *= 1. - hit.localPosition().y() / theDetPlane.position().perp();
0653 
0654 #ifdef FAMOS_DEBUG
0655   unsigned detid = DetId(hit.detUnitId()).rawId();
0656   unsigned stereo = 0;
0657   unsigned theLayer = 0;
0658   unsigned theRing = 0;
0659   switch (subdet) {
0660     case 1: {
0661       theLayer = tTopo->pxbLayer(detid);
0662       std::cout << "\tPixel Barrel Layer " << theLayer << std::endl;
0663       stereo = 1;
0664       break;
0665     }
0666     case 2: {
0667       theLayer = tTopo->pxfDisk(detid);
0668       std::cout << "\tPixel Forward Disk " << theLayer << std::endl;
0669       stereo = 1;
0670       break;
0671     }
0672     case 3: {
0673       theLayer = tTopo->tibLayer(detid);
0674       std::cout << "\tTIB Layer " << theLayer << std::endl;
0675       stereo = module.stereo();
0676       break;
0677     }
0678     case 4: {
0679       theLayer = tTopo->tidWheel(detid);
0680       theRing = tTopo->tidRing(detid);
0681       unsigned int theSide = module.side();
0682       if (theSide == 1)
0683         std::cout << "\tTID Petal Back " << std::endl;
0684       else
0685         std::cout << "\tTID Petal Front" << std::endl;
0686       std::cout << "\tTID Layer " << theLayer << std::endl;
0687       std::cout << "\tTID Ring " << theRing << std::endl;
0688       stereo = module.stereo();
0689       break;
0690     }
0691     case 5: {
0692       theLayer = tTopo->tobLayer(detid);
0693       stereo = tTopo->tobStereo(detid);
0694       std::cout << "\tTOB Layer " << theLayer << std::endl;
0695       break;
0696     }
0697     case 6: {
0698       theLayer = tTopo->tecWheel(detid);
0699       theRing = tTopo->tecRing(detid);
0700       unsigned int theSide = module.petal()[0];
0701       if (theSide == 1)
0702         std::cout << "\tTEC Petal Back " << std::endl;
0703       else
0704         std::cout << "\tTEC Petal Front" << std::endl;
0705       std::cout << "\tTEC Layer " << theLayer << std::endl;
0706       std::cout << "\tTEC Ring " << theRing << std::endl;
0707       stereo = module.stereo();
0708       break;
0709     }
0710     default: {
0711       stereo = 0;
0712       break;
0713     }
0714   }
0715 
0716   std::cout << "Thickness = " << 2. * halfThick - 0.003 << "; " << thick * 9.36 << std::endl
0717             << "Length    = " << det.surface().bounds().length() << std::endl
0718             << "Width     = " << det.surface().bounds().width() << std::endl;
0719 
0720   std::cout << "Hit position = " << hit.localPosition().x() << " " << hit.localPosition().y() << " "
0721             << hit.localPosition().z() << std::endl;
0722 #endif
0723 
0724   // Check if the hit is on the physical volume of the module
0725   // (It happens that it is not, in the case of double sided modules,
0726   //  because the envelope of the gluedDet is larger than each of
0727   //  the mono and the stereo modules)
0728 
0729   double dist = 0.;
0730   GlobalPoint IP(mySimEvent->track(localTkID).vertex().position().x(),
0731                  mySimEvent->track(localTkID).vertex().position().y(),
0732                  mySimEvent->track(localTkID).vertex().position().z());
0733 
0734   dist = (fabs(hit.localPosition().x()) > boundX || fabs(hit.localPosition().y()) > boundY)
0735              ?
0736              // Will be used later as a flag to reject the PSimHit!
0737              -(det.surface().toGlobal(hit.localPosition()) - IP).mag2()
0738              :
0739              // These hits are kept!
0740              (det.surface().toGlobal(hit.localPosition()) - IP).mag2();
0741 
0742   // Fill Histos (~poor man event display)
0743   /*  
0744       GlobalPoint gpos( det.toGlobal(hit.localPosition()));
0745       //      std::cout << "gpos.x() = " << gpos.x() << std::endl;
0746       //      std::cout << "gpos.y() = " << gpos.y() << std::endl;
0747 
0748       myHistos->fill("h300",gpos.x(),gpos.y());
0749       if ( sin(gpos.phi()) > 0. ) 
0750       myHistos->fill("h301",gpos.z(),gpos.perp());
0751       else
0752       myHistos->fill("h301",gpos.z(),-gpos.perp());
0753   */
0754   return std::pair<double, PSimHit>(dist, hit);
0755 }
0756 
0757 void TrajectoryManager::initializeLayerMap() {
0758   // These are the BoundSurface&, the BoundDisk* and the BoundCylinder* for that layer
0759   //   const BoundSurface& theSurface = layer.surface();
0760   //   BoundDisk* theDisk = layer.disk();  // non zero for endcaps
0761   //   BoundCylinder* theCylinder = layer.cylinder(); // non zero for barrel
0762   //   int theLayer = layer.layerNumber(); // 1->3 PixB, 4->5 PixD,
0763   //                                       // 6->9 TIB, 10->12 TID,
0764   //                                       // 13->18 TOB, 19->27 TEC
0765 
0766   /// ATTENTION: HARD CODED LOGIC! If Famos layer numbering changes this logic needs to
0767   /// be adapted to the new numbering!
0768 
0769   const std::vector<const BarrelDetLayer*>& barrelLayers = theGeomSearchTracker->barrelLayers();
0770   LogDebug("FastTracking") << "Barrel DetLayer dump: ";
0771   for (auto bl = barrelLayers.begin(); bl != barrelLayers.end(); ++bl) {
0772     LogDebug("FastTracking") << "radius " << (**bl).specificSurface().radius();
0773   }
0774 
0775   const std::vector<const ForwardDetLayer*>& posForwardLayers = theGeomSearchTracker->posForwardLayers();
0776   LogDebug("FastTracking") << "Positive Forward DetLayer dump: ";
0777   for (auto fl = posForwardLayers.begin(); fl != posForwardLayers.end(); ++fl) {
0778     LogDebug("FastTracking") << "Z pos " << (**fl).surface().position().z() << " radii "
0779                              << (**fl).specificSurface().innerRadius() << ", "
0780                              << (**fl).specificSurface().outerRadius();
0781   }
0782 
0783   const float rTolerance = 1.5;
0784   const float zTolerance = 3.;
0785 
0786   LogDebug("FastTracking") << "Dump of TrackerInteractionGeometry cylinders:";
0787   for (std::list<TrackerLayer>::const_iterator i = _theGeometry->cylinderBegin(); i != _theGeometry->cylinderEnd();
0788        ++i) {
0789     const BoundCylinder* cyl = i->cylinder();
0790     const BoundDisk* disk = i->disk();
0791 
0792     LogDebug("FastTracking") << "Famos Layer no " << i->layerNumber() << " is sensitive? " << i->sensitive() << " pos "
0793                              << i->surface().position();
0794     if (!i->sensitive())
0795       continue;
0796 
0797     if (cyl != nullptr) {
0798       LogDebug("FastTracking") << " cylinder radius " << cyl->radius();
0799       bool found = false;
0800       for (auto bl = barrelLayers.begin(); bl != barrelLayers.end(); ++bl) {
0801         if (fabs(cyl->radius() - (**bl).specificSurface().radius()) < rTolerance) {
0802           theLayerMap[i->layerNumber()] = *bl;
0803           found = true;
0804           LogDebug("FastTracking") << "Corresponding DetLayer found with radius " << (**bl).specificSurface().radius();
0805           break;
0806         }
0807       }
0808       if (!found) {
0809         edm::LogWarning("FastTracking") << " Trajectory manager FAILED to find a corresponding DetLayer!";
0810       }
0811     } else {
0812       LogDebug("FastTracking") << " disk radii " << disk->innerRadius() << ", " << disk->outerRadius();
0813       bool found = false;
0814       for (auto fl = posForwardLayers.begin(); fl != posForwardLayers.end(); ++fl) {
0815         if (fabs(disk->position().z() - (**fl).surface().position().z()) < zTolerance) {
0816           theLayerMap[i->layerNumber()] = *fl;
0817           found = true;
0818           LogDebug("FastTracking") << "Corresponding DetLayer found with Z pos " << (**fl).surface().position().z()
0819                                    << " and radii " << (**fl).specificSurface().innerRadius() << ", "
0820                                    << (**fl).specificSurface().outerRadius();
0821           break;
0822         }
0823       }
0824       if (!found) {
0825         edm::LogWarning("FastTracking") << "FAILED to find a corresponding DetLayer!";
0826       }
0827     }
0828   }
0829 
0830   // Put the negative layers in the same map but with an offset
0831   const std::vector<const ForwardDetLayer*>& negForwardLayers = theGeomSearchTracker->negForwardLayers();
0832   for (auto nl = negForwardLayers.begin(); nl != negForwardLayers.end(); ++nl) {
0833     for (int i = 0; i <= theNegLayerOffset; i++) {
0834       if (theLayerMap[i] == nullptr)
0835         continue;
0836       if (fabs((**nl).surface().position().z() + theLayerMap[i]->surface().position().z()) < zTolerance) {
0837         theLayerMap[i + theNegLayerOffset] = *nl;
0838         break;
0839       }
0840     }
0841   }
0842 }
0843 
0844 const DetLayer* TrajectoryManager::detLayer(const TrackerLayer& layer, float zpos) const {
0845   if (zpos > 0 || !layer.forward())
0846     return theLayerMap[layer.layerNumber()];
0847   else
0848     return theLayerMap[layer.layerNumber() + theNegLayerOffset];
0849 }
0850 
0851 void TrajectoryManager::loadSimHits(edm::PSimHitContainer& c) const {
0852   std::map<unsigned, std::map<double, PSimHit> >::const_iterator itrack = thePSimHits.begin();
0853   std::map<unsigned, std::map<double, PSimHit> >::const_iterator itrackEnd = thePSimHits.end();
0854   for (; itrack != itrackEnd; ++itrack) {
0855     std::map<double, PSimHit>::const_iterator it = (itrack->second).begin();
0856     std::map<double, PSimHit>::const_iterator itEnd = (itrack->second).end();
0857     for (; it != itEnd; ++it) {
0858       /*
0859     DetId theDetUnitId((it->second).detUnitId());
0860     const GeomDet* theDet = theGeomTracker->idToDet(theDetUnitId);
0861     std::cout << "Track/z/r after : "
0862     << (it->second).trackId() << " " 
0863     << theDet->surface().toGlobal((it->second).localPosition()).z() << " " 
0864     << theDet->surface().toGlobal((it->second).localPosition()).perp() << std::endl;
0865       */
0866       // Keep only those hits that are on the physical volume of a module
0867       // (The other hits have been assigned a negative <double> value.
0868       if (it->first > 0.)
0869         c.push_back(it->second);
0870     }
0871   }
0872 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!