CMS 3D CMS Logo

List of all members | Public Types | Public Member Functions | Private Member Functions | Private Attributes
TrajectoryManager Class Reference

#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, RandomEngineAndDistribution const *)
 Propagate the particle through the calorimeters. More...
 
bool propagateToLayer (ParticlePropagator &PP, unsigned layer)
 
void reconstruct (const TrackerTopology *tTopo, RandomEngineAndDistribution const *)
 Does the real job. More...
 
const TrackerInteractionGeometrytheGeometry ()
 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)
 Constructor from a FSimEvent. More...
 
 ~TrajectoryManager ()
 Default Destructor. More...
 

Private Member Functions

const DetLayerdetLayer (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, PSimHitmakeSinglePSimHit (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, RandomEngineAndDistribution const *)
 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
 
PythiaDecaysmyDecayEngine
 
FSimEventmySimEvent
 
double pTmin
 
const GeometricSearchTrackertheGeomSearchTracker
 
const TrackerGeometrytheGeomTracker
 
std::vector< const DetLayer * > theLayerMap
 
MaterialEffectstheMaterialEffects
 
int theNegLayerOffset
 
std::map< unsigned, std::map< double, PSimHit > > thePSimHits
 
bool use_hardcoded
 

Detailed Description

Definition at line 59 of file TrajectoryManager.h.

Member Typedef Documentation

◆ Rotation

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

Definition at line 63 of file TrajectoryManager.h.

Constructor & Destructor Documentation

◆ TrajectoryManager() [1/2]

TrajectoryManager::TrajectoryManager ( )
inline

Default Constructor.

Definition at line 66 of file TrajectoryManager.h.

66 { ; }

◆ TrajectoryManager() [2/2]

TrajectoryManager::TrajectoryManager ( FSimEvent aSimEvent,
const edm::ParameterSet matEff,
const edm::ParameterSet simHits,
const edm::ParameterSet decays 
)

Constructor from a FSimEvent.

Definition at line 47 of file TrajectoryManager.cc.

References bJpsiMuMuTrigSettings_cff::decays, distCut, firstLoop, edm::ParameterSet::getParameter(), MaterialEffects_cfi::MaterialEffects, myDecayEngine, pTmin, FastTrackerRecHitCombiner_cfi::simHits, theMaterialEffects, and use_hardcoded.

51  : mySimEvent(aSimEvent),
52  _theGeometry(nullptr),
53  _theFieldMap(nullptr),
54  theMaterialEffects(nullptr),
55  myDecayEngine(nullptr),
56  theGeomTracker(nullptr),
57  theGeomSearchTracker(nullptr),
58  theLayerMap(56, static_cast<const DetLayer*>(nullptr)), // reserve space for layers here
60  // myHistos(0),
61  use_hardcoded(true)
62 
63 {
64  //std::cout << "TrajectoryManager.cc 1 use_hardcoded = " << use_hardcoded << std::endl;
65  use_hardcoded = matEff.getParameter<bool>("use_hardcoded_geometry");
66 
67  // Initialize Bthe stable particle decay engine
68  if (decays.getParameter<bool>("ActivateDecays")) {
70  distCut = decays.getParameter<double>("DistCut");
71  }
72  // Initialize the Material Effects updator, if needed
73  if (matEff.getParameter<bool>("PairProduction") || matEff.getParameter<bool>("Bremsstrahlung") ||
74  matEff.getParameter<bool>("MuonBremsstrahlung") || matEff.getParameter<bool>("EnergyLoss") ||
75  matEff.getParameter<bool>("MultipleScattering") || matEff.getParameter<bool>("NuclearInteraction"))
76  theMaterialEffects = new MaterialEffects(matEff);
77 
78  // Save SimHits according to Optiom
79  // Only the hits from first half loop is saved
80  firstLoop = simHits.getUntrackedParameter<bool>("firstLoop", true);
81  // Only if pT>pTmin are the hits saved
82  pTmin = simHits.getUntrackedParameter<double>("pTmin", 0.5);
83 
84  /*
85  // Get the Famos Histos pointer
86  myHistos = Histos::instance();
87 
88  // Initialize a few histograms
89 
90  myHistos->book("h302",1210,-121.,121.,1210,-121.,121.);
91  myHistos->book("h300",1210,-121.,121.,1210,-121.,121.);
92  myHistos->book("h301",1200,-300.,300.,1210,-121.,121.);
93  myHistos->book("h303",1200,-300.,300.,1210,-121.,121.);
94  */
95 }
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
PythiaDecays * myDecayEngine
const MagneticFieldMap * _theFieldMap
const TrackerGeometry * theGeomTracker
const GeometricSearchTracker * theGeomSearchTracker
const TrackerInteractionGeometry * _theGeometry
MaterialEffects * theMaterialEffects
std::vector< const DetLayer * > theLayerMap

◆ ~TrajectoryManager()

TrajectoryManager::~TrajectoryManager ( )

Default Destructor.

Definition at line 116 of file TrajectoryManager.cc.

References myDecayEngine, and theMaterialEffects.

116  {
117  if (myDecayEngine)
118  delete myDecayEngine;
119  if (theMaterialEffects)
120  delete theMaterialEffects;
121 
122  //Write the histograms
123  /*
124  myHistos->put("histos.root");
125  if ( myHistos ) delete myHistos;
126  */
127 }
PythiaDecays * myDecayEngine
MaterialEffects * theMaterialEffects

Member Function Documentation

◆ createPSimHits()

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

References anyDirection, GeometricSearchDet::compatibleDets(), detLayer(), MaterialEffects::energyLoss(), BaseParticlePropagator::getMagneticField(), mps_fire::i, nano_mu_digi_cff::layer, makePSimHits(), makeTrajectoryState(), BaseParticlePropagator::particle(), theMaterialEffects, MaterialEffects::thickness(), Calorimetry_cff::thickness, and RawParticle::Z().

Referenced by reconstruct().

525  {
526  // Propagate the particle coordinates to the closest tracker detector(s)
527  // in this layer and create the PSimHit(s)
528 
529  // const MagneticField& mf = MagneticFieldMap::instance()->magneticField();
530  // This solution is actually much faster !
531  LocalMagneticField mf(PP.getMagneticField());
532  AnalyticalPropagator alongProp(&mf, anyDirection);
534 
535  // std::cout << "PP.X() = " << PP.X() << std::endl;
536  // std::cout << "PP.Y() = " << PP.Y() << std::endl;
537  // std::cout << "PP.Z() = " << PP.Z() << std::endl;
538 
540  const DetLayer* tkLayer = detLayer(layer, PP.particle().Z());
541 
542  TrajectoryStateOnSurface trajState = makeTrajectoryState(tkLayer, PP, &mf);
544  float eloss = theMaterialEffects ? theMaterialEffects->energyLoss() : 0.;
545 
546  // Find, in the corresponding layers, the detectors compatible
547  // with the current track
548  std::vector<DetWithState> compat = tkLayer->compatibleDets(trajState, alongProp, est);
549 
550  // And create the corresponding PSimHits
551  std::map<double, PSimHit> theTrackHits;
552  for (std::vector<DetWithState>::const_iterator i = compat.begin(); i != compat.end(); i++) {
553  // Correct Eloss for last 3 rings of TEC (thick sensors, 0.05 cm)
554  // Disgusting fudge factor !
555  makePSimHits(i->first, i->second, theHitMap, trackID, eloss, thickness, partID, tTopo);
556  }
557 }
virtual std::vector< DetWithState > compatibleDets(const TrajectoryStateOnSurface &startingState, const Propagator &prop, const MeasurementEstimator &est) const
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
const DetLayer * detLayer(const TrackerLayer &layer, float zpos) const
Returns the DetLayer pointer corresponding to the FAMOS layer.
double energyLoss() const
Return the energy loss by ionization in the current layer.
std::pair< const GeomDet *, TrajectoryStateOnSurface > DetWithState
MaterialEffects * theMaterialEffects
double thickness() const
Return the thickness of the current layer.
TrajectoryStateOnSurface makeTrajectoryState(const DetLayer *layer, const ParticlePropagator &pp, const MagneticField *field) const
Teddy, you must put comments there.
std::pair< const GeomDet *, TrajectoryStateOnSurface > DetWithState

◆ detLayer()

const DetLayer * TrajectoryManager::detLayer ( const TrackerLayer layer,
float  zpos 
) const
private

Returns the DetLayer pointer corresponding to the FAMOS layer.

Definition at line 844 of file TrajectoryManager.cc.

References nano_mu_digi_cff::layer, theLayerMap, and theNegLayerOffset.

Referenced by createPSimHits().

844  {
845  if (zpos > 0 || !layer.forward())
846  return theLayerMap[layer.layerNumber()];
847  else
848  return theLayerMap[layer.layerNumber() + theNegLayerOffset];
849 }
std::vector< const DetLayer * > theLayerMap

◆ initializeLayerMap()

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

References _theGeometry, GeometricSearchTracker::barrelLayers(), TrackerInteractionGeometry::cylinderBegin(), TrackerInteractionGeometry::cylinderEnd(), newFWLiteAna::found, mps_fire::i, LogDebug, GeometricSearchTracker::negForwardLayers(), GeometricSearchTracker::posForwardLayers(), theGeomSearchTracker, theLayerMap, and theNegLayerOffset.

Referenced by initializeRecoGeometry().

757  {
758  // These are the BoundSurface&, the BoundDisk* and the BoundCylinder* for that layer
759  // const BoundSurface& theSurface = layer.surface();
760  // BoundDisk* theDisk = layer.disk(); // non zero for endcaps
761  // BoundCylinder* theCylinder = layer.cylinder(); // non zero for barrel
762  // int theLayer = layer.layerNumber(); // 1->3 PixB, 4->5 PixD,
763  // // 6->9 TIB, 10->12 TID,
764  // // 13->18 TOB, 19->27 TEC
765 
768 
769  const std::vector<const BarrelDetLayer*>& barrelLayers = theGeomSearchTracker->barrelLayers();
770  LogDebug("FastTracking") << "Barrel DetLayer dump: ";
771  for (auto bl = barrelLayers.begin(); bl != barrelLayers.end(); ++bl) {
772  LogDebug("FastTracking") << "radius " << (**bl).specificSurface().radius();
773  }
774 
775  const std::vector<const ForwardDetLayer*>& posForwardLayers = theGeomSearchTracker->posForwardLayers();
776  LogDebug("FastTracking") << "Positive Forward DetLayer dump: ";
777  for (auto fl = posForwardLayers.begin(); fl != posForwardLayers.end(); ++fl) {
778  LogDebug("FastTracking") << "Z pos " << (**fl).surface().position().z() << " radii "
779  << (**fl).specificSurface().innerRadius() << ", "
780  << (**fl).specificSurface().outerRadius();
781  }
782 
783  const float rTolerance = 1.5;
784  const float zTolerance = 3.;
785 
786  LogDebug("FastTracking") << "Dump of TrackerInteractionGeometry cylinders:";
787  for (std::list<TrackerLayer>::const_iterator i = _theGeometry->cylinderBegin(); i != _theGeometry->cylinderEnd();
788  ++i) {
789  const BoundCylinder* cyl = i->cylinder();
790  const BoundDisk* disk = i->disk();
791 
792  LogDebug("FastTracking") << "Famos Layer no " << i->layerNumber() << " is sensitive? " << i->sensitive() << " pos "
793  << i->surface().position();
794  if (!i->sensitive())
795  continue;
796 
797  if (cyl != nullptr) {
798  LogDebug("FastTracking") << " cylinder radius " << cyl->radius();
799  bool found = false;
800  for (auto bl = barrelLayers.begin(); bl != barrelLayers.end(); ++bl) {
801  if (fabs(cyl->radius() - (**bl).specificSurface().radius()) < rTolerance) {
802  theLayerMap[i->layerNumber()] = *bl;
803  found = true;
804  LogDebug("FastTracking") << "Corresponding DetLayer found with radius " << (**bl).specificSurface().radius();
805  break;
806  }
807  }
808  if (!found) {
809  edm::LogWarning("FastTracking") << " Trajectory manager FAILED to find a corresponding DetLayer!";
810  }
811  } else {
812  LogDebug("FastTracking") << " disk radii " << disk->innerRadius() << ", " << disk->outerRadius();
813  bool found = false;
814  for (auto fl = posForwardLayers.begin(); fl != posForwardLayers.end(); ++fl) {
815  if (fabs(disk->position().z() - (**fl).surface().position().z()) < zTolerance) {
816  theLayerMap[i->layerNumber()] = *fl;
817  found = true;
818  LogDebug("FastTracking") << "Corresponding DetLayer found with Z pos " << (**fl).surface().position().z()
819  << " and radii " << (**fl).specificSurface().innerRadius() << ", "
820  << (**fl).specificSurface().outerRadius();
821  break;
822  }
823  }
824  if (!found) {
825  edm::LogWarning("FastTracking") << "FAILED to find a corresponding DetLayer!";
826  }
827  }
828  }
829 
830  // Put the negative layers in the same map but with an offset
831  const std::vector<const ForwardDetLayer*>& negForwardLayers = theGeomSearchTracker->negForwardLayers();
832  for (auto nl = negForwardLayers.begin(); nl != negForwardLayers.end(); ++nl) {
833  for (int i = 0; i <= theNegLayerOffset; i++) {
834  if (theLayerMap[i] == nullptr)
835  continue;
836  if (fabs((**nl).surface().position().z() + theLayerMap[i]->surface().position().z()) < zTolerance) {
838  break;
839  }
840  }
841  }
842 }
std::vector< ForwardDetLayer const * > const & negForwardLayers() const
const GeometricSearchTracker * theGeomSearchTracker
const TrackerInteractionGeometry * _theGeometry
std::list< TrackerLayer >::const_iterator cylinderEnd() const
Returns the last pointer in the cylinder list.
std::vector< ForwardDetLayer const * > const & posForwardLayers() const
std::list< TrackerLayer >::const_iterator cylinderBegin() const
Returns the first pointer in the cylinder list.
std::vector< BarrelDetLayer const * > const & barrelLayers() const
std::vector< const DetLayer * > theLayerMap
Log< level::Warning, false > LogWarning
#define LogDebug(id)

◆ initializeRecoGeometry()

void TrajectoryManager::initializeRecoGeometry ( const GeometricSearchTracker geomSearchTracker,
const TrackerInteractionGeometry interactionGeometry,
const MagneticFieldMap aFieldMap 
)

Initialize the Reconstruction Geometry.

Definition at line 97 of file TrajectoryManager.cc.

References _theFieldMap, _theGeometry, initializeLayerMap(), and theGeomSearchTracker.

99  {
100  // Initialize the reco tracker geometry
101  theGeomSearchTracker = geomSearchTracker;
102 
103  // Initialize the simplified tracker geometry
104  _theGeometry = interactionGeometry;
105 
107 
108  // Initialize the magnetic field
109  _theFieldMap = aFieldMap;
110 }
void initializeLayerMap()
Initialize correspondence map between Famos interaction geometry and tracker reco geometry...
const MagneticFieldMap * _theFieldMap
const GeometricSearchTracker * theGeomSearchTracker
const TrackerInteractionGeometry * _theGeometry

◆ initializeTrackerGeometry()

void TrajectoryManager::initializeTrackerGeometry ( const TrackerGeometry geomTracker)

Initialize the full Tracker Geometry.

Definition at line 112 of file TrajectoryManager.cc.

References theGeomTracker.

112 { theGeomTracker = geomTracker; }
const TrackerGeometry * theGeomTracker

◆ loadSimHits()

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

Definition at line 851 of file TrajectoryManager.cc.

References HltBtagPostValidation_cff::c, mps_fire::end, ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::it, and thePSimHits.

851  {
852  std::map<unsigned, std::map<double, PSimHit> >::const_iterator itrack = thePSimHits.begin();
853  std::map<unsigned, std::map<double, PSimHit> >::const_iterator itrackEnd = thePSimHits.end();
854  for (; itrack != itrackEnd; ++itrack) {
855  std::map<double, PSimHit>::const_iterator it = (itrack->second).begin();
856  std::map<double, PSimHit>::const_iterator itEnd = (itrack->second).end();
857  for (; it != itEnd; ++it) {
858  /*
859  DetId theDetUnitId((it->second).detUnitId());
860  const GeomDet* theDet = theGeomTracker->idToDet(theDetUnitId);
861  std::cout << "Track/z/r after : "
862  << (it->second).trackId() << " "
863  << theDet->surface().toGlobal((it->second).localPosition()).z() << " "
864  << theDet->surface().toGlobal((it->second).localPosition()).perp() << std::endl;
865  */
866  // Keep only those hits that are on the physical volume of a module
867  // (The other hits have been assigned a negative <double> value.
868  if (it->first > 0.)
869  c.push_back(it->second);
870  }
871  }
872 }
std::map< unsigned, std::map< double, PSimHit > > thePSimHits

◆ makePSimHits()

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

References AlCaHLTBitMon_QueryRunRegistry::comp, GeomDet::components(), mps_fire::i, and makeSinglePSimHit().

Referenced by createPSimHits().

576  {
577  std::vector<const GeomDet*> comp = det->components();
578  if (!comp.empty()) {
579  for (std::vector<const GeomDet*>::const_iterator i = comp.begin(); i != comp.end(); i++) {
580  auto du = (*i);
581  if (du->isLeaf()) // not even needed (or it should iterate if really not leaf)
582  theHitMap.insert(theHitMap.end(), makeSinglePSimHit(*du, ts, tkID, el, thick, pID, tTopo));
583  }
584  } else {
585  auto du = (det);
586  theHitMap.insert(theHitMap.end(), makeSinglePSimHit(*du, ts, tkID, el, thick, pID, tTopo));
587  }
588 }
virtual std::vector< const GeomDet * > components() const
Returns direct components, if any.
Definition: GeomDet.h:73
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

◆ makeSinglePSimHit()

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

References anyDirection, PV3DBase< T, PVType, FrameType >::basicVector(), Surface::bounds(), FSimTrack::closestDaughterId(), gather_cfg::cout, mps_splice::entry, beamvalidation::exit(), GeomDet::geographicalId(), TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), FSimTrack::id(), listHistos::IP, Bounds::length(), TrajectoryStateOnSurface::localMomentum(), TrajectoryStateOnSurface::localPosition(), PV3DBase< T, PVType, FrameType >::mag(), FSimTrack::mother(), mySimEvent, FSimTrack::noMother(), castor_dqm_sourceclient_file_cfg::path, FSimVertex::position(), TrackerTopology::pxbLayer(), TrackerTopology::pxfDisk(), DetId::rawId(), nano_mu_digi_cff::rawId, GeomDet::surface(), TrackerTopology::tecRing(), TrackerTopology::tecWheel(), TrackerTopology::tibLayer(), TrackerTopology::tidRing(), TrackerTopology::tidWheel(), TrackerTopology::tobLayer(), TrackerTopology::tobStereo(), Surface::toGlobal(), GeomDet::toLocal(), FBaseSimEvent::track(), TrajectoryStateOnSurface::transverseCurvature(), Vector3DBase< T, FrameTag >::unit(), FSimTrack::vertex(), Bounds::width(), hit::x, hit::y, z, PV3DBase< T, PVType, FrameType >::z(), and hit::z.

Referenced by makePSimHits().

596  {
597  const float onSurfaceTolarance = 0.01; // 10 microns
598 
599  LocalPoint lpos;
600  LocalVector lmom;
601  if (fabs(det.toLocal(ts.globalPosition()).z()) < onSurfaceTolarance) {
602  lpos = ts.localPosition();
603  lmom = ts.localMomentum();
604  } else {
607  std::pair<bool, double> path = crossing.pathLength(det.surface());
608  if (!path.first) {
609  // edm::LogWarning("FastTracking") << "TrajectoryManager ERROR: crossing with det failed, skipping PSimHit";
610  return std::pair<double, PSimHit>(0., PSimHit());
611  }
612  lpos = det.toLocal(GlobalPoint(crossing.position(path.second)));
613  lmom = det.toLocal(GlobalVector(crossing.direction(path.second)));
614  lmom = lmom.unit() * ts.localMomentum().mag();
615  }
616 
617  // The module (half) thickness
618  const BoundPlane& theDetPlane = det.surface();
619  float halfThick = 0.5 * theDetPlane.bounds().thickness();
620  // The Energy loss rescaled to the module thickness
621  float eloss = el;
622  if (thick > 0.) {
623  // Total thickness is in radiation lengths, 1 radlen = 9.36 cm
624  // Sensitive module thickness is about 30 microns larger than
625  // the module thickness itself
626  eloss *= (2. * halfThick - 0.003) / (9.36 * thick);
627  }
628  // The entry and exit points, and the time of flight
629  float pZ = lmom.z();
630  LocalPoint entry = lpos + (-halfThick / pZ) * lmom;
631  LocalPoint exit = lpos + halfThick / pZ * lmom;
632  float tof = ts.globalPosition().mag() / 30.; // in nanoseconds, FIXME: very approximate
633 
634  // If a hadron suffered a nuclear interaction, just assign the hits of the closest
635  // daughter to the mother's track. The same applies to a charged particle decay into
636  // another charged particle.
637  int localTkID = tkID;
638  if (!mySimEvent->track(tkID).noMother() && mySimEvent->track(tkID).mother().closestDaughterId() == tkID)
639  localTkID = mySimEvent->track(tkID).mother().id();
640 
641  // FIXME: fix the track ID and the particle ID
642  PSimHit hit(
643  entry, exit, lmom.mag(), tof, eloss, pID, det.geographicalId().rawId(), localTkID, lmom.theta(), lmom.phi());
644 
645  // Check that the PSimHit is physically on the module!
646  unsigned subdet = DetId(hit.detUnitId()).subdetId();
647  double boundX = theDetPlane.bounds().width() / 2.;
648  double boundY = theDetPlane.bounds().length() / 2.;
649 
650  // Special treatment for TID and TEC trapeziodal modules
651  if (subdet == 4 || subdet == 6)
652  boundX *= 1. - hit.localPosition().y() / theDetPlane.position().perp();
653 
654 #ifdef FAMOS_DEBUG
655  unsigned detid = DetId(hit.detUnitId()).rawId();
656  unsigned stereo = 0;
657  unsigned theLayer = 0;
658  unsigned theRing = 0;
659  switch (subdet) {
660  case 1: {
661  theLayer = tTopo->pxbLayer(detid);
662  std::cout << "\tPixel Barrel Layer " << theLayer << std::endl;
663  stereo = 1;
664  break;
665  }
666  case 2: {
667  theLayer = tTopo->pxfDisk(detid);
668  std::cout << "\tPixel Forward Disk " << theLayer << std::endl;
669  stereo = 1;
670  break;
671  }
672  case 3: {
673  theLayer = tTopo->tibLayer(detid);
674  std::cout << "\tTIB Layer " << theLayer << std::endl;
675  stereo = module.stereo();
676  break;
677  }
678  case 4: {
679  theLayer = tTopo->tidWheel(detid);
680  theRing = tTopo->tidRing(detid);
681  unsigned int theSide = module.side();
682  if (theSide == 1)
683  std::cout << "\tTID Petal Back " << std::endl;
684  else
685  std::cout << "\tTID Petal Front" << std::endl;
686  std::cout << "\tTID Layer " << theLayer << std::endl;
687  std::cout << "\tTID Ring " << theRing << std::endl;
688  stereo = module.stereo();
689  break;
690  }
691  case 5: {
692  theLayer = tTopo->tobLayer(detid);
693  stereo = tTopo->tobStereo(detid);
694  std::cout << "\tTOB Layer " << theLayer << std::endl;
695  break;
696  }
697  case 6: {
698  theLayer = tTopo->tecWheel(detid);
699  theRing = tTopo->tecRing(detid);
700  unsigned int theSide = module.petal()[0];
701  if (theSide == 1)
702  std::cout << "\tTEC Petal Back " << std::endl;
703  else
704  std::cout << "\tTEC Petal Front" << std::endl;
705  std::cout << "\tTEC Layer " << theLayer << std::endl;
706  std::cout << "\tTEC Ring " << theRing << std::endl;
707  stereo = module.stereo();
708  break;
709  }
710  default: {
711  stereo = 0;
712  break;
713  }
714  }
715 
716  std::cout << "Thickness = " << 2. * halfThick - 0.003 << "; " << thick * 9.36 << std::endl
717  << "Length = " << det.surface().bounds().length() << std::endl
718  << "Width = " << det.surface().bounds().width() << std::endl;
719 
720  std::cout << "Hit position = " << hit.localPosition().x() << " " << hit.localPosition().y() << " "
721  << hit.localPosition().z() << std::endl;
722 #endif
723 
724  // Check if the hit is on the physical volume of the module
725  // (It happens that it is not, in the case of double sided modules,
726  // because the envelope of the gluedDet is larger than each of
727  // the mono and the stereo modules)
728 
729  double dist = 0.;
730  GlobalPoint IP(mySimEvent->track(localTkID).vertex().position().x(),
731  mySimEvent->track(localTkID).vertex().position().y(),
732  mySimEvent->track(localTkID).vertex().position().z());
733 
734  dist = (fabs(hit.localPosition().x()) > boundX || fabs(hit.localPosition().y()) > boundY)
735  ?
736  // Will be used later as a flag to reject the PSimHit!
737  -(det.surface().toGlobal(hit.localPosition()) - IP).mag2()
738  :
739  // These hits are kept!
740  (det.surface().toGlobal(hit.localPosition()) - IP).mag2();
741 
742  // Fill Histos (~poor man event display)
743  /*
744  GlobalPoint gpos( det.toGlobal(hit.localPosition()));
745  // std::cout << "gpos.x() = " << gpos.x() << std::endl;
746  // std::cout << "gpos.y() = " << gpos.y() << std::endl;
747 
748  myHistos->fill("h300",gpos.x(),gpos.y());
749  if ( sin(gpos.phi()) > 0. )
750  myHistos->fill("h301",gpos.z(),gpos.perp());
751  else
752  myHistos->fill("h301",gpos.z(),-gpos.perp());
753  */
754  return std::pair<double, PSimHit>(dist, hit);
755 }
bool noMother() const
no mother particle
FSimTrack & track(int id) const
Return track with given Id.
unsigned int tobLayer(const DetId &id) const
unsigned int pxbLayer(const DetId &id) const
virtual float length() const =0
LocalPoint toLocal(const GlobalPoint &gp) const
Conversion to the R.F. of the GeomDet.
Definition: GeomDet.h:58
T z() const
Definition: PV3DBase.h:61
unsigned int tidWheel(const DetId &id) const
unsigned int tecWheel(const DetId &id) const
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
const FSimVertex vertex() const
Origin vertex.
unsigned int tecRing(const DetId &id) const
ring id
uint32_t tobStereo(const DetId &id) const
GlobalPoint globalPosition() const
T mag() const
Definition: PV3DBase.h:64
unsigned int pxfDisk(const DetId &id) const
DetId geographicalId() const
The label of this GeomDet.
Definition: GeomDet.h:64
const BasicVectorType & basicVector() const
Definition: PV3DBase.h:53
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:79
Definition: DetId.h:17
int closestDaughterId() const
Get the index of the closest charged daughter.
Definition: FSimTrack.h:206
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:37
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:57
const math::XYZTLorentzVector & position() const
Temporary (until CMSSW moves to Mathcore) - No ! Actually very useful.
Definition: FSimVertex.h:48
GlobalVector globalMomentum() const
unsigned int tidRing(const DetId &id) const
unsigned int tibLayer(const DetId &id) const
Vector3DBase unit() const
Definition: Vector3DBase.h:54
virtual float width() const =0
int id() const
the index in FBaseSimEvent and other vectors
Definition: FSimTrack.h:96
const FSimTrack & mother() const
mother
Global3DVector GlobalVector
Definition: GlobalVector.h:10
def exit(msg="")
const Bounds & bounds() const
Definition: Surface.h:87

◆ makeTrajectoryState()

TrajectoryStateOnSurface TrajectoryManager::makeTrajectoryState ( const DetLayer layer,
const ParticlePropagator pp,
const MagneticField field 
) const
private

Teddy, you must put comments there.

Definition at line 559 of file TrajectoryManager.cc.

References nano_mu_digi_cff::layer, and createTree::pp.

Referenced by createPSimHits().

561  {
562  GlobalPoint pos(pp.particle().X(), pp.particle().Y(), pp.particle().Z());
563  GlobalVector mom(pp.particle().Px(), pp.particle().Py(), pp.particle().Pz());
564  auto plane = layer->surface().tangentPlane(pos);
565  return TrajectoryStateOnSurface(GlobalTrajectoryParameters(pos, mom, TrackCharge(pp.particle().charge()), field),
566  *plane);
567 }
int TrackCharge
Definition: TrackCharge.h:4

◆ moveAllDaughters()

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

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

Referenced by updateWithDaughters().

502  {
503  //
504  for (unsigned idaugh = 0; idaugh < (unsigned)(mySimEvent->track(fsimi).nDaughters()); ++idaugh) {
505  // Initial momentum of the daughter
506  XYZTLorentzVector daughMomentum(mySimEvent->track(fsimi).daughter(idaugh).momentum());
507  // Rotate and rescale
508  XYZVector newMomentum(r * daughMomentum.Vect());
509  newMomentum *= rescale;
510  double newEnergy = std::sqrt(newMomentum.mag2() + daughMomentum.mag2());
511  // Set the new momentum
512  mySimEvent->track(fsimi).setMomentum(
513  XYZTLorentzVector(newMomentum.X(), newMomentum.Y(), newMomentum.Z(), newEnergy));
514  // Watch out : recursive call to get all grand-daughters
515  int fsimDaug = mySimEvent->track(fsimi).daughter(idaugh).id();
516  moveAllDaughters(fsimDaug, r, rescale);
517  }
518 }
FSimTrack & track(int id) const
Return track with given Id.
void moveAllDaughters(int fsimi, const Rotation &r, double rescale)
Move, rescale and rotate all daughters after propagation, material effects and decay of the mother...
const FSimTrack & daughter(int i) const
Ith daughter.
const XYZTLorentzVector & momentum() const
Temporary (until move of SimTrack to Mathcore) - No! Actually very useful.
Definition: FSimTrack.h:209
T sqrt(T t)
Definition: SSEVec.h:19
int nDaughters() const
Number of daughters.
math::XYZVector XYZVector
Definition: RawParticle.h:26
void setMomentum(const math::XYZTLorentzVector &newMomentum)
Reset the momentum (to be used with care)
Definition: FSimTrack.h:212
int id() const
the index in FBaseSimEvent and other vectors
Definition: FSimTrack.h:96
math::XYZTLorentzVector XYZTLorentzVector
Definition: RawParticle.h:25

◆ propagateToCalorimeters()

void TrajectoryManager::propagateToCalorimeters ( ParticlePropagator PP,
int  fsimi,
RandomEngineAndDistribution const *  random 
)

Propagate the particle through the calorimeters.

Definition at line 361 of file TrajectoryManager.cc.

References BaseParticlePropagator::getSuccess(), BaseParticlePropagator::hasDecayed(), RawParticle::momentum(), mySimEvent, FSimTrack::notYetToEndVertex(), BaseParticlePropagator::particle(), 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().

363  {
364  FSimTrack& myTrack = mySimEvent->track(fsimi);
365 
366  // Set the position and momentum at the end of the tracker volume
367  myTrack.setTkPosition(PP.particle().vertex().Vect());
368  myTrack.setTkMomentum(PP.particle().momentum());
369 
370  // Propagate to Preshower Layer 1
371  PP.propagateToPreshowerLayer1(false);
372  if (PP.hasDecayed()) {
373  updateWithDaughters(PP, fsimi, random);
374  return;
375  }
376  if (myTrack.notYetToEndVertex(PP.particle().vertex()) && PP.getSuccess() > 0)
377  myTrack.setLayer1(PP.particle(), PP.getSuccess());
378 
379  // Propagate to Preshower Layer 2
380  PP.propagateToPreshowerLayer2(false);
381  if (PP.hasDecayed()) {
382  updateWithDaughters(PP, fsimi, random);
383  return;
384  }
385  if (myTrack.notYetToEndVertex(PP.particle().vertex()) && PP.getSuccess() > 0)
386  myTrack.setLayer2(PP.particle(), PP.getSuccess());
387 
388  // Propagate to Ecal Endcap
389  PP.propagateToEcalEntrance(false);
390  if (PP.hasDecayed()) {
391  updateWithDaughters(PP, fsimi, random);
392  return;
393  }
394  if (myTrack.notYetToEndVertex(PP.particle().vertex()))
395  myTrack.setEcal(PP.particle(), PP.getSuccess());
396 
397  // Propagate to HCAL entrance
398  PP.propagateToHcalEntrance(false);
399  if (PP.hasDecayed()) {
400  updateWithDaughters(PP, fsimi, random);
401  return;
402  }
403  if (myTrack.notYetToEndVertex(PP.particle().vertex()))
404  myTrack.setHcal(PP.particle(), PP.getSuccess());
405 
406  // Propagate to VFCAL entrance
407  PP.propagateToVFcalEntrance(false);
408  if (PP.hasDecayed()) {
409  updateWithDaughters(PP, fsimi, random);
410  return;
411  }
412  if (myTrack.notYetToEndVertex(PP.particle().vertex()))
413  myTrack.setVFcal(PP.particle(), PP.getSuccess());
414 }
FSimTrack & track(int id) const
Return track with given Id.
bool propagateToPreshowerLayer1(bool first=true)
void setEcal(const RawParticle &pp, int success)
Set the ecal variables.
Definition: FSimTrack.cc:118
void setLayer2(const RawParticle &pp, int success)
Set the preshower layer2 variables.
Definition: FSimTrack.cc:112
const XYZTLorentzVector & momentum() const
the momentum fourvector
Definition: RawParticle.h:321
bool propagateToVFcalEntrance(bool first=true)
void setLayer1(const RawParticle &pp, int success)
Set the preshower layer1 variables.
Definition: FSimTrack.cc:106
bool hasDecayed() const
Has the particle decayed while propagated ?
RawParticle const & particle() const
The particle being propagated.
bool notYetToEndVertex(const XYZTLorentzVector &pos) const
Compare the end vertex position with another position.
Definition: FSimTrack.cc:85
bool propagateToEcalEntrance(bool first=true)
int getSuccess() const
Has propagation been performed and was barrel or endcap reached ?
void setTkPosition(const math::XYZVectorD &pos)
Definition: SimTrack.h:44
bool propagateToHcalEntrance(bool first=true)
void setTkMomentum(const math::XYZTLorentzVectorD &mom)
Definition: SimTrack.h:46
void setVFcal(const RawParticle &pp, int success)
Set the hcal variables.
Definition: FSimTrack.cc:130
void setHcal(const RawParticle &pp, int success)
Set the hcal variables.
Definition: FSimTrack.cc:124
void updateWithDaughters(ParticlePropagator &PP, int fsimi, RandomEngineAndDistribution const *)
Decay the particle and update the SimEvent with daughters.
bool propagateToPreshowerLayer2(bool first=true)
const XYZTLorentzVector & vertex() const
the vertex fourvector
Definition: RawParticle.h:320

◆ propagateToLayer()

bool TrajectoryManager::propagateToLayer ( ParticlePropagator PP,
unsigned  layer 
)

Propagate a particle to a given tracker layer (for electron pixel matching mostly)

Definition at line 416 of file TrajectoryManager.cc.

References _theGeometry, TrackerInteractionGeometry::cylinderBegin(), TrackerInteractionGeometry::cylinderEnd(), fileCollector::done, BaseParticlePropagator::getSuccess(), nano_mu_digi_cff::layer, BaseParticlePropagator::onFiducial(), ParticlePropagator::propagateToBoundSurface(), and ParticlePropagator::setPropagationConditions().

416  {
417  std::list<TrackerLayer>::const_iterator cyliter;
418  bool done = false;
419 
420  // Get the geometry elements
421  cyliter = _theGeometry->cylinderBegin();
422 
423  // Find the layer to propagate to.
424  for (; cyliter != _theGeometry->cylinderEnd(); ++cyliter) {
425  if (layer != cyliter->layerNumber())
426  continue;
427 
428  PP.setPropagationConditions(*cyliter);
429 
430  done = PP.propagateToBoundSurface(*cyliter) && PP.getSuccess() > 0 && PP.onFiducial();
431 
432  break;
433  }
434 
435  return done;
436 }
bool propagateToBoundSurface(const TrackerLayer &)
bool onFiducial() const
Is the vertex on some material ?
const TrackerInteractionGeometry * _theGeometry
void setPropagationConditions(const TrackerLayer &, bool firstLoop=true)
std::list< TrackerLayer >::const_iterator cylinderEnd() const
Returns the last pointer in the cylinder list.
std::list< TrackerLayer >::const_iterator cylinderBegin() const
Returns the first pointer in the cylinder list.
int getSuccess() const
Has propagation been performed and was barrel or endcap reached ?

◆ reconstruct()

void TrajectoryManager::reconstruct ( const TrackerTopology tTopo,
RandomEngineAndDistribution const *  random 
)

Does the real job.

Definition at line 129 of file TrajectoryManager.cc.

References _theFieldMap, _theGeometry, KineParticleFilter::acceptParticle(), FBaseSimEvent::addSimVertex(), createPSimHits(), pdg::cTau(), TrackerInteractionGeometry::cylinderBegin(), TrackerInteractionGeometry::cylinderEnd(), FSimVertexType::END_VERTEX, FBaseSimEvent::filter(), firstLoop, createfilelist::int, MaterialEffects::interact(), heppy_loop::loop, mySimEvent, FSimTrack::nDaughters(), FSimTrack::notYetToEndVertex(), FSimEvent::nTracks(), propagateToCalorimeters(), pTmin, MaterialEffects::save(), FSimTrack::setPropagate(), Validation_hcalonly_cfi::sign, summarizeEdmComparisonLogfiles::success, theGeomTracker, theMaterialEffects, thePSimHits, FBaseSimEvent::theTable(), FBaseSimEvent::track(), CoreSimTrack::type(), updateWithDaughters(), and use_hardcoded.

129  {
130  // Clear the hits of the previous event
131  // thePSimHits->clear();
132  thePSimHits.clear();
133 
134  // The new event
135  XYZTLorentzVector myBeamPipe = XYZTLorentzVector(0., 2.5, 9999999., 0.);
136 
137  std::list<TrackerLayer>::const_iterator cyliter;
138 
139  // bool debug = mySimEvent->id().event() == 8;
140 
141  // Loop over the particles (watch out: increasing upper limit!)
142  for (int fsimi = 0; fsimi < (int)mySimEvent->nTracks(); ++fsimi) {
143  // If the particle has decayed inside the beampipe, or decays
144  // immediately, there is nothing to do
145  //if ( debug ) std::cout << mySimEvent->track(fsimi) << std::endl;
146  //if ( debug ) std::cout << "Not yet at end vertex ? " << mySimEvent->track(fsimi).notYetToEndVertex(myBeamPipe) << std::endl;
147  if (!mySimEvent->track(fsimi).notYetToEndVertex(myBeamPipe))
148  continue;
149  mySimEvent->track(fsimi).setPropagate();
150 
151  // Get the geometry elements
152  cyliter = _theGeometry->cylinderBegin();
153  // Prepare the propagation
155  //The real work starts here
156  int success = 1;
157  int sign = +1;
158  int loop = 0;
159  int cyl = 0;
160 
161  // Find the initial cylinder to propagate to.
162  for (; cyliter != _theGeometry->cylinderEnd(); ++cyliter) {
163  PP.setPropagationConditions(*cyliter);
164  if (PP.inside() && !PP.onSurface())
165  break;
166  ++cyl;
167  }
168 
169  // The particle has a pseudo-rapidity (position or momentum direction)
170  // in excess of 3.0. Just simply go to the last tracker layer
171  // without bothering with all the details of the propagation and
172  // material effects.
173  // 08/02/06 - pv: increase protection from 0.99 (eta=2.9932) to 0.9998 (eta=4.9517)
174  // to simulate material effects at large eta
175  // if above 0.99: propagate to the last tracker cylinder where the material is concentrated!
176  double ppcos2T = PP.particle().cos2Theta();
177  double ppcos2V = PP.particle().cos2ThetaV();
178 
179  if (use_hardcoded) {
180  if ((ppcos2T > 0.99 && ppcos2T < 0.9998) && (cyl == 0 || (ppcos2V > 0.99 && ppcos2V < 0.9998))) {
181  if (cyliter != _theGeometry->cylinderEnd()) {
182  cyliter = _theGeometry->cylinderEnd();
183  --cyliter;
184  }
185  // if above 0.9998: don't propagate at all (only to the calorimeters directly)
186  } else if (ppcos2T > 0.9998 && (cyl == 0 || ppcos2V > 0.9998)) {
187  cyliter = _theGeometry->cylinderEnd();
188  }
189  } else {
190  if (ppcos2T > 0.9998 && (cyl == 0 || ppcos2V > 0.9998)) {
191  cyliter = _theGeometry->cylinderEnd();
192  }
193  }
194 
195  // Loop over the cylinders
196  while (cyliter != _theGeometry->cylinderEnd() && loop < 100 && // No more than 100 loops
197  mySimEvent->track(fsimi).notYetToEndVertex(PP.particle().vertex())) { // The particle decayed
198 
199  // Skip layers with no material (kept just for historical reasons)
200  if (cyliter->surface().mediumProperties().radLen() < 1E-10) {
201  ++cyliter;
202  ++cyl;
203  continue;
204  }
205 
206  // Pathological cases:
207  // To prevent from interacting twice in a row with the same layer
208  // bool escapeBarrel = (PP.getSuccess() == -1 && success == 1);
209  bool escapeBarrel = PP.getSuccess() == -1;
210  bool escapeEndcap = (PP.getSuccess() == -2 && success == 1);
211  // To break the loop
212  bool fullPropagation = (PP.getSuccess() <= 0 && success == 0) || escapeEndcap;
213 
214  if (escapeBarrel) {
215  ++cyliter;
216  ++cyl;
217  while (cyliter != _theGeometry->cylinderEnd() && cyliter->forward()) {
218  sign = 1;
219  ++cyliter;
220  ++cyl;
221  }
222 
223  if (cyliter == _theGeometry->cylinderEnd()) {
224  --cyliter;
225  --cyl;
226  fullPropagation = true;
227  }
228  }
229 
230  // Define the propagation conditions
231  PP.setPropagationConditions(*cyliter, !fullPropagation);
232  if (escapeEndcap)
233  PP.increaseRCyl(0.0005);
234 
235  // Remember last propagation outcome
236  success = PP.getSuccess();
237 
238  // Propagation was not successful :
239  // Change the sign of the cylinder increment and count the loops
240  if (!PP.propagateToBoundSurface(*cyliter) || PP.getSuccess() <= 0) {
241  sign = -sign;
242  ++loop;
243  }
244 
245  // The particle may have decayed on its way... in which the daughters
246  // have to be added to the event record
247  if (PP.hasDecayed() ||
248  (!mySimEvent->track(fsimi).nDaughters() && pdg::cTau(PP.particle().pid(), mySimEvent->theTable()) < 1E-3)) {
249  updateWithDaughters(PP, fsimi, random);
250  break;
251  }
252 
253  // Exit by the endcaps or innermost cylinder :
254  // Positive cylinder increment
255  if (PP.getSuccess() == 2 || cyliter == _theGeometry->cylinderBegin())
256  sign = +1;
257 
258  // Successful propagation to a cylinder, with some Material :
259  if (PP.getSuccess() > 0 && PP.onFiducial()) {
260  bool saveHit = ((loop == 0 && sign > 0) || !firstLoop) && // Save only first half loop
261  PP.particle().charge() != 0. && // Consider only charged particles
262  cyliter->sensitive() && // Consider only sensitive layers
263  PP.particle().Perp2() > pTmin * pTmin; // Consider only pT > pTmin
264 
265  // Material effects are simulated there
266  if (theMaterialEffects)
267  theMaterialEffects->interact(*mySimEvent, *cyliter, PP, fsimi, random);
268 
269  // There is a PP.setXYZT=(0,0,0,0) if bremss fails
270  saveHit &= PP.particle().E() > 1E-6;
271 
272  if (saveHit) {
273  // Consider only active layers
274  if (cyliter->sensitive()) {
275  // Add information to the FSimTrack (not yet available)
276  // myTrack.addSimHit(PP.particle(),layer);
277 
278  // Return one or two (for overlap regions) PSimHits in the full
279  // tracker geometry
280  if (theGeomTracker)
281  createPSimHits(*cyliter, PP, thePSimHits[fsimi], fsimi, mySimEvent->track(fsimi).type(), tTopo);
282 
283  /*
284  myHistos->fill("h302",PP.particle().X() ,PP.particle().Y());
285  if ( sin(PP.particle().vertex().Phi()) > 0. )
286  myHistos->fill("h303",PP.particle().Z(),PP.particle().R());
287  else
288  myHistos->fill("h303",PP.Z(),-PP.particle().R());
289  */
290  }
291  }
292 
293  // Fill Histos (~poor man event display)
294  /*
295  myHistos->fill("h300",PP.particle().x(),PP.particle().y());
296  if ( sin(PP.particle().vertex().phi()) > 0. )
297  myHistos->fill("h301",PP.particle().z(),sqrt(PP.particle().vertex().Perp2()));
298  else
299  myHistos->fill("h301",PP.particle().z(),-sqrt(PP.particle().vertex().Perp2()));
300  */
301 
302  //The particle may have lost its energy in the material
303  if (mySimEvent->track(fsimi).notYetToEndVertex(PP.particle().vertex()) &&
304  !mySimEvent->filter().acceptParticle(PP.particle()))
305  mySimEvent->addSimVertex(PP.particle().vertex(), fsimi, FSimVertexType::END_VERTEX);
306  }
307 
308  // Stop here if the particle has reached an end
309  if (mySimEvent->track(fsimi).notYetToEndVertex(PP.particle().vertex())) {
310  // Otherwise increment the cylinder iterator
311  // do {
312  if (sign == 1) {
313  ++cyliter;
314  ++cyl;
315  } else {
316  --cyliter;
317  --cyl;
318  }
319 
320  // Check if the last surface has been reached
321  if (cyliter == _theGeometry->cylinderEnd()) {
322  // Try to propagate to the ECAL in half a loop
323  // Note: Layer1 = ECAL Barrel entrance, or Preshower
324  // entrance, or ECAL Endcap entrance (in the corner)
325  PP.propagateToEcal();
326  // PP.propagateToPreshowerLayer1();
327 
328  // If it is not possible, try go back to the last cylinder
329  if (PP.getSuccess() == 0) {
330  --cyliter;
331  --cyl;
332  sign = -sign;
333  PP.setPropagationConditions(*cyliter);
334  PP.propagateToBoundSurface(*cyliter);
335 
336  // If there is definitely no way, leave it here.
337  if (PP.getSuccess() < 0) {
338  ++cyliter;
339  ++cyl;
340  }
341  }
342 
343  // Check if the particle has decayed on the way to ECAL
344  if (PP.hasDecayed())
345  updateWithDaughters(PP, fsimi, random);
346  }
347  }
348  }
349 
350  // Propagate all particles without a end vertex to the Preshower,
351  // theECAL and the HCAL.
352  if (mySimEvent->track(fsimi).notYetToEndVertex(PP.particle().vertex()))
353  propagateToCalorimeters(PP, fsimi, random);
354  }
355 
356  // Save the information from Nuclear Interaction Generation
357  if (theMaterialEffects)
359 }
int addSimVertex(const XYZTLorentzVector &decayVertex, int im=-1, FSimVertexType::VertexType type=FSimVertexType::ANY)
Add a new vertex to the Event and to the various lists.
FSimTrack & track(int id) const
Return track with given Id.
void propagateToCalorimeters(ParticlePropagator &PP, int fsimi, RandomEngineAndDistribution const *)
Propagate the particle through the calorimeters.
double cTau(int pdgID, const HepPDT::ParticleDataTable *pdt)
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.
const MagneticFieldMap * _theFieldMap
const TrackerGeometry * theGeomTracker
void setPropagate()
The particle has been propgated through the tracker.
Definition: FSimTrack.cc:103
bool acceptParticle(const RawParticle &p) const
unsigned int nTracks() const
Number of tracks.
Definition: FSimEvent.cc:24
const HepPDT::ParticleDataTable * theTable() const
Get the pointer to the particle data table.
Definition: FBaseSimEvent.h:54
int type() const
particle type (HEP PDT convension)
Definition: CoreSimTrack.h:22
const TrackerInteractionGeometry * _theGeometry
MaterialEffects * theMaterialEffects
void save()
Save nuclear interaction information.
std::list< TrackerLayer >::const_iterator cylinderEnd() const
Returns the last pointer in the cylinder list.
bool notYetToEndVertex(const XYZTLorentzVector &pos) const
Compare the end vertex position with another position.
Definition: FSimTrack.cc:85
std::map< unsigned, std::map< double, PSimHit > > thePSimHits
int nDaughters() const
Number of daughters.
void interact(FSimEvent &simEvent, const TrackerLayer &layer, ParticlePropagator &PP, unsigned i, RandomEngineAndDistribution const *)
std::list< TrackerLayer >::const_iterator cylinderBegin() const
Returns the first pointer in the cylinder list.
const KineParticleFilter & filter() const
void updateWithDaughters(ParticlePropagator &PP, int fsimi, RandomEngineAndDistribution const *)
Decay the particle and update the SimEvent with daughters.
math::XYZTLorentzVector XYZTLorentzVector
Definition: RawParticle.h:25

◆ theGeometry()

const TrackerInteractionGeometry * TrajectoryManager::theGeometry ( )

Returns the pointer to geometry.

Definition at line 114 of file TrajectoryManager.cc.

References _theGeometry.

114 { return _theGeometry; }
const TrackerInteractionGeometry * _theGeometry

◆ updateWithDaughters()

void TrajectoryManager::updateWithDaughters ( ParticlePropagator PP,
int  fsimi,
RandomEngineAndDistribution const *  random 
)
private

Decay the particle and update the SimEvent with daughters.

Definition at line 438 of file TrajectoryManager.cc.

References FBaseSimEvent::addSimTrack(), FBaseSimEvent::addSimVertex(), angle(), RawParticle::charge(), BPhysicsValidation_cfi::daughters, FSimVertexType::DECAY_VERTEX, distCut, muonTagProbeFilters_cff::distMin, FSimTrack::endVertex(), FSimVertex::id(), FSimTrack::momentum(), RawParticle::momentum(), moveAllDaughters(), myDecayEngine, mySimEvent, FSimTrack::nDaughters(), BaseParticlePropagator::particle(), PythiaDecays::particleDaughters(), alignCSCRings::r, FSimTrack::setClosestDaughterId(), FSimVertex::setPosition(), mathSSE::sqrt(), RandomEngineAndDistribution::theEngine(), FBaseSimEvent::track(), RawParticle::Vect(), FBaseSimEvent::vertex(), and RawParticle::vertex().

Referenced by propagateToCalorimeters(), and reconstruct().

440  {
441  // The particle was already decayed in the GenEvent, but still the particle was
442  // allowed to propagate (for magnetic field bending, for material effects, etc...)
443  // Just modify the momentum of the daughters in that case
444  unsigned nDaugh = mySimEvent->track(fsimi).nDaughters();
445  if (nDaugh) {
446  // Move the vertex
447  unsigned vertexId = mySimEvent->track(fsimi).endVertex().id();
448  mySimEvent->vertex(vertexId).setPosition(PP.particle().vertex());
449 
450  // Before-propagation and after-propagation momentum and vertex position
451  XYZTLorentzVector momentumBefore = mySimEvent->track(fsimi).momentum();
452  const XYZTLorentzVector& momentumAfter = PP.particle().momentum();
453  double magBefore = std::sqrt(momentumBefore.Vect().mag2());
454  double magAfter = std::sqrt(momentumAfter.Vect().mag2());
455  // Rotation to be applied
456  XYZVector axis = momentumBefore.Vect().Cross(momentumAfter.Vect());
457  double angle = std::acos(momentumBefore.Vect().Dot(momentumAfter.Vect()) / (magAfter * magBefore));
458  Rotation r(axis, angle);
459  // Rescaling to be applied
460  double rescale = magAfter / magBefore;
461 
462  // Move, rescale and rotate daugthers, grand-daughters, etc.
463  moveAllDaughters(fsimi, r, rescale);
464 
465  // The particle is not decayed in the GenEvent, decay it with PYTHIA
466  } else {
467  // Decays are not activated : do nothing
468  if (!myDecayEngine)
469  return;
470 
471  // Invoke PYDECY (Pythia6) or Pythia8 to decay the particle and get the daughters
472  const DaughterParticleList& daughters = myDecayEngine->particleDaughters(PP, &random->theEngine());
473 
474  // Update the FSimEvent with an end vertex and with the daughters
475  if (!daughters.empty()) {
476  double distMin = 1E99;
477  int theClosestChargedDaughterId = -1;
478  DaughterParticleIterator daughter = daughters.begin();
479 
480  int ivertex = mySimEvent->addSimVertex(daughter->vertex(), fsimi, FSimVertexType::DECAY_VERTEX);
481 
482  if (ivertex != -1) {
483  for (; daughter != daughters.end(); ++daughter) {
484  int theDaughterId = mySimEvent->addSimTrack(&(*daughter), ivertex);
485  // Find the closest charged daughter (if charged mother)
486  if (PP.particle().charge() * daughter->charge() > 1E-10) {
487  double dist = (daughter->Vect().Unit().Cross(PP.particle().Vect().Unit())).R();
488  if (dist < distCut && dist < distMin) {
489  distMin = dist;
490  theClosestChargedDaughterId = theDaughterId;
491  }
492  }
493  }
494  }
495  // Attach mother and closest daughter sp as to cheat tracking ;-)
496  if (theClosestChargedDaughterId >= 0)
497  mySimEvent->track(fsimi).setClosestDaughterId(theClosestChargedDaughterId);
498  }
499  }
500 }
int addSimVertex(const XYZTLorentzVector &decayVertex, int im=-1, FSimVertexType::VertexType type=FSimVertexType::ANY)
Add a new vertex to the Event and to the various lists.
FSimTrack & track(int id) const
Return track with given Id.
PythiaDecays * myDecayEngine
void setPosition(const math::XYZTLorentzVector &newPosition)
Reset the position (to be used with care)
Definition: FSimVertex.h:51
int addSimTrack(const RawParticle *p, int iv, int ig=-1, const HepMC::GenVertex *ev=nullptr)
Add a new track to the Event and to the various lists.
void moveAllDaughters(int fsimi, const Rotation &r, double rescale)
Move, rescale and rotate all daughters after propagation, material effects and decay of the mother...
const DaughterParticleList & particleDaughters(ParticlePropagator &particle, CLHEP::HepRandomEngine *)
Definition: PythiaDecays.cc:34
DaughterParticleList::const_iterator DaughterParticleIterator
Definition: PythiaDecays.h:26
void setClosestDaughterId(int id)
Set the index of the closest charged daughter.
Definition: FSimTrack.h:203
const XYZTLorentzVector & momentum() const
the momentum fourvector
Definition: RawParticle.h:321
ROOT::Math::AxisAngle Rotation
double charge() const
get the MEASURED charge
Definition: RawParticle.h:294
const FSimVertex & endVertex() const
end vertex
const XYZTLorentzVector & momentum() const
Temporary (until move of SimTrack to Mathcore) - No! Actually very useful.
Definition: FSimTrack.h:209
int id() const
the index in FBaseSimEvent
Definition: FSimVertex.h:43
T sqrt(T t)
Definition: SSEVec.h:19
RawParticle const & particle() const
The particle being propagated.
int nDaughters() const
Number of daughters.
XYZVector Vect() const
the momentum threevector
Definition: RawParticle.h:323
FSimVertex & vertex(int id) const
Return vertex with given Id.
math::XYZVector XYZVector
Definition: RawParticle.h:26
math::XYZTLorentzVector XYZTLorentzVector
Definition: RawParticle.h:25
std::vector< RawParticle > DaughterParticleList
Definition: PythiaDecays.h:25
const XYZTLorentzVector & vertex() const
the vertex fourvector
Definition: RawParticle.h:320
T angle(T x1, T y1, T z1, T x2, T y2, T z2)
Definition: angle.h:11

Member Data Documentation

◆ _theFieldMap

const MagneticFieldMap* TrajectoryManager::_theFieldMap
private

Definition at line 149 of file TrajectoryManager.h.

Referenced by initializeRecoGeometry(), and reconstruct().

◆ _theGeometry

const TrackerInteractionGeometry* TrajectoryManager::_theGeometry
private

◆ decayer

std::string TrajectoryManager::decayer
private

Definition at line 154 of file TrajectoryManager.h.

◆ distCut

double TrajectoryManager::distCut
private

Definition at line 155 of file TrajectoryManager.h.

Referenced by TrajectoryManager(), and updateWithDaughters().

◆ firstLoop

bool TrajectoryManager::firstLoop
private

Definition at line 158 of file TrajectoryManager.h.

Referenced by reconstruct(), and TrajectoryManager().

◆ myDecayEngine

PythiaDecays* TrajectoryManager::myDecayEngine
private

Definition at line 153 of file TrajectoryManager.h.

Referenced by TrajectoryManager(), updateWithDaughters(), and ~TrajectoryManager().

◆ mySimEvent

FSimEvent* TrajectoryManager::mySimEvent
private

◆ pTmin

double TrajectoryManager::pTmin
private

Definition at line 157 of file TrajectoryManager.h.

Referenced by reconstruct(), and TrajectoryManager().

◆ theGeomSearchTracker

const GeometricSearchTracker* TrajectoryManager::theGeomSearchTracker
private

Definition at line 162 of file TrajectoryManager.h.

Referenced by initializeLayerMap(), and initializeRecoGeometry().

◆ theGeomTracker

const TrackerGeometry* TrajectoryManager::theGeomTracker
private

Definition at line 161 of file TrajectoryManager.h.

Referenced by initializeTrackerGeometry(), and reconstruct().

◆ theLayerMap

std::vector<const DetLayer*> TrajectoryManager::theLayerMap
private

Definition at line 163 of file TrajectoryManager.h.

Referenced by detLayer(), and initializeLayerMap().

◆ theMaterialEffects

MaterialEffects* TrajectoryManager::theMaterialEffects
private

◆ theNegLayerOffset

int TrajectoryManager::theNegLayerOffset
private

Definition at line 164 of file TrajectoryManager.h.

Referenced by detLayer(), and initializeLayerMap().

◆ thePSimHits

std::map<unsigned, std::map<double, PSimHit> > TrajectoryManager::thePSimHits
private

Definition at line 159 of file TrajectoryManager.h.

Referenced by loadSimHits(), and reconstruct().

◆ use_hardcoded

bool TrajectoryManager::use_hardcoded
private

Definition at line 168 of file TrajectoryManager.h.

Referenced by reconstruct(), and TrajectoryManager().