#include <FastSimulation/MuonSimHitProducer/src/MuonSimHitProducer.cc>

Inheritance diagram for MuonSimHitProducer:

Public Member Functions
	MuonSimHitProducer (const edm::ParameterSet &)

	~MuonSimHitProducer () override

Public Member Functions inherited from edm::stream::EDProducer<>
	EDProducer ()=default

bool	hasAbilityToProduceInBeginLumis () const final

bool	hasAbilityToProduceInBeginRuns () const final

bool	hasAbilityToProduceInEndLumis () const final

bool	hasAbilityToProduceInEndRuns () const final

Private Member Functions
void	applyMaterialEffects (TrajectoryStateOnSurface &tsosWithdEdx, TrajectoryStateOnSurface &tsos, double radPath, RandomEngineAndDistribution const *, HepPDT::ParticleDataTable const &)
	Simulate material effects in iron (dE/dx, multiple scattering) More...

void	beginRun (edm::Run const &run, const edm::EventSetup &es) override

void	produce (edm::Event &, const edm::EventSetup &) override

void	readParameters (const edm::ParameterSet &, const edm::ParameterSet &, const edm::ParameterSet &)

Private Attributes
const CSCGeometry *	cscGeom

bool	doGL_

bool	doL1_

bool	doL3_

const DTGeometry *	dtGeom

double	fCSC

double	fDT

bool	fullPattern_

double	kCSC

double	kDT

const MagneticField *	magfield

const Propagator *	propagatorWithMaterial

Propagator *	propagatorWithoutMaterial

const RPCGeometry *	rpcGeom

edm::InputTag	simMuonLabel

edm::EDGetTokenT< std::vector< SimTrack > >	simMuonToken

edm::InputTag	simVertexLabel

edm::EDGetTokenT< std::vector< SimVertex > >	simVertexToken

Chi2MeasurementEstimator	theEstimator

MaterialEffects *	theMaterialEffects

MuonServiceProxy *	theService

Additional Inherited Members
Public Types inherited from edm::stream::EDProducer<>
typedef CacheContexts< T... >	CacheTypes

typedef CacheTypes::GlobalCache	GlobalCache

typedef AbilityChecker< T... >	HasAbility

typedef CacheTypes::LuminosityBlockCache	LuminosityBlockCache

typedef LuminosityBlockContextT< LuminosityBlockCache, RunCache, GlobalCache >	LuminosityBlockContext

typedef CacheTypes::LuminosityBlockSummaryCache	LuminosityBlockSummaryCache

typedef CacheTypes::RunCache	RunCache

typedef RunContextT< RunCache, GlobalCache >	RunContext

typedef CacheTypes::RunSummaryCache	RunSummaryCache

Detailed Description

Description: <one line="" class="" summary>="">

Implementation: <Notes on="" implementation>="">

Description: Fast simulation producer of Muon Sim Hits (to be used for realistic Muon reconstruction)

Implementation: <Notes on="" implementation>="">

Definition at line 58 of file MuonSimHitProducer.h.

Constructor & Destructor Documentation

◆ MuonSimHitProducer()

MuonSimHitProducer::MuonSimHitProducer ( const edm::ParameterSet & iConfig )

explicit

Definition at line 78 of file MuonSimHitProducer.cc.

     : theEstimator(iConfig.getParameter<double>("Chi2EstimatorCut")), propagatorWithoutMaterial(nullptr) {
   // Read relevant parameters
   readParameters(iConfig.getParameter<edm::ParameterSet>("MUONS"),
                  iConfig.getParameter<edm::ParameterSet>("TRACKS"),
                  iConfig.getParameter<edm::ParameterSet>("MaterialEffectsForMuons"));
  
   //
   //  register your products ... need to declare at least one possible product...
   //
   produces<edm::PSimHitContainer>("MuonCSCHits");
   produces<edm::PSimHitContainer>("MuonDTHits");
   produces<edm::PSimHitContainer>("MuonRPCHits");
  
   edm::ParameterSet serviceParameters = iConfig.getParameter<edm::ParameterSet>("ServiceParameters");
   theService = new MuonServiceProxy(serviceParameters, consumesCollector(), MuonServiceProxy::UseEventSetupIn::Run);
  
   // consumes
   simMuonToken = consumes<std::vector<SimTrack> >(simMuonLabel);
   simVertexToken = consumes<std::vector<SimVertex> >(simVertexLabel);
 }

References edm::ParameterSet::getParameter(), MuonServiceProxy_cff::MuonServiceProxy, readParameters(), MuonServiceProxy::Run, simMuonLabel, simMuonToken, simVertexLabel, simVertexToken, and theService.

◆ ~MuonSimHitProducer()

MuonSimHitProducer::~MuonSimHitProducer ( )

override

Definition at line 129 of file MuonSimHitProducer.cc.

                                         {
   // do anything here that needs to be done at destruction time
   // (e.g. close files, deallocate resources etc.)
  
   if (theMaterialEffects)
     delete theMaterialEffects;
   if (propagatorWithoutMaterial)
     delete propagatorWithoutMaterial;
 }

References propagatorWithoutMaterial, and theMaterialEffects.

Member Function Documentation

◆ applyMaterialEffects()

void MuonSimHitProducer::applyMaterialEffects	(	TrajectoryStateOnSurface &	tsosWithdEdx,
		TrajectoryStateOnSurface &	tsos,
		double	radPath,
		RandomEngineAndDistribution const *	random,
		HepPDT::ParticleDataTable const &	table
	)

private

Simulate material effects in iron (dE/dx, multiple scattering)

Definition at line 538 of file MuonSimHitProducer.cc.

                                                                                     {
   // The energy loss simulator
   EnergyLossSimulator* energyLoss = theMaterialEffects->energyLossSimulator();
  
   // The multiple scattering simulator
   MultipleScatteringSimulator* multipleScattering = theMaterialEffects->multipleScatteringSimulator();
  
   // The Muon Bremsstrahlung simulator
   MuonBremsstrahlungSimulator* bremsstrahlung = theMaterialEffects->muonBremsstrahlungSimulator();
  
   // Initialize the Particle position, momentum and energy
   const Surface& nextSurface = tsos.surface();
   GlobalPoint gPos = energyLoss ? tsos.globalPosition() : tsosWithdEdx.globalPosition();
   GlobalVector gMom = energyLoss ? tsos.globalMomentum() : tsosWithdEdx.globalMomentum();
   double mu = 0.1056583692;
   double en = std::sqrt(gMom.mag2() + mu * mu);
  
   // And now create the Particle
   XYZTLorentzVector position(gPos.x(), gPos.y(), gPos.z(), 0.);
   XYZTLorentzVector momentum(gMom.x(), gMom.y(), gMom.z(), en);
   float charge = (float)(tsos.charge());
   ParticlePropagator theMuon(rawparticle::makeMuon(charge < 1., momentum, position), nullptr, nullptr, &table);
  
   // Recompute the energy loss to get the fluctuations
   if (energyLoss) {
     // Difference between with and without dE/dx (average only)
     // (for corrections once fluctuations are applied)
     GlobalPoint gPosWithdEdx = tsosWithdEdx.globalPosition();
     GlobalVector gMomWithdEdx = tsosWithdEdx.globalMomentum();
     double enWithdEdx = std::sqrt(gMomWithdEdx.mag2() + mu * mu);
     XYZTLorentzVector deltaPos(
         gPosWithdEdx.x() - gPos.x(), gPosWithdEdx.y() - gPos.y(), gPosWithdEdx.z() - gPos.z(), 0.);
     XYZTLorentzVector deltaMom(
         gMomWithdEdx.x() - gMom.x(), gMomWithdEdx.y() - gMom.y(), gMomWithdEdx.z() - gMom.z(), enWithdEdx - en);
  
     // Energy loss in iron (+ fluctuations)
     energyLoss->updateState(theMuon, radPath, random);
  
     // Correcting factors to account for slight differences in material descriptions
     // (Material description is more accurate in the stepping helix propagator)
     radPath *= -deltaMom.E() / energyLoss->mostLikelyLoss();
     double fac = energyLoss->deltaMom().E() / energyLoss->mostLikelyLoss();
  
     // Particle momentum & position after energy loss + fluctuation
     XYZTLorentzVector theNewMomentum = theMuon.particle().momentum() + energyLoss->deltaMom() + fac * deltaMom;
     XYZTLorentzVector theNewPosition = theMuon.particle().vertex() + fac * deltaPos;
     fac = (theNewMomentum.E() * theNewMomentum.E() - mu * mu) / theNewMomentum.Vect().Mag2();
     fac = fac > 0. ? std::sqrt(fac) : 1E-9;
     theMuon.particle().setMomentum(
         theNewMomentum.Px() * fac, theNewMomentum.Py() * fac, theNewMomentum.Pz() * fac, theNewMomentum.E());
     theMuon.particle().setVertex(theNewPosition);
   }
  
   // Does the actual mutliple scattering
   if (multipleScattering) {
     // Pass the vector normal to the "next" surface
     GlobalVector normal = nextSurface.tangentPlane(tsos.globalPosition())->normalVector();
     multipleScattering->setNormalVector(normal);
     // Compute the amount of multiple scattering after a given path length
     multipleScattering->updateState(theMuon, radPath, random);
   }
  
   // Muon Bremsstrahlung
   if (bremsstrahlung) {
     // Compute the amount of Muon Bremsstrahlung after given path length
     bremsstrahlung->updateState(theMuon, radPath, random);
   }
  
   // Fill the propagated state
   GlobalPoint propagatedPosition(theMuon.particle().X(), theMuon.particle().Y(), theMuon.particle().Z());
   GlobalVector propagatedMomentum(theMuon.particle().Px(), theMuon.particle().Py(), theMuon.particle().Pz());
   GlobalTrajectoryParameters propagatedGtp(propagatedPosition, propagatedMomentum, (int)charge, magfield);
   tsosWithdEdx = TrajectoryStateOnSurface(propagatedGtp, nextSurface);
 }

References fastSimProducer_cff::bremsstrahlung, ALCARECOTkAlJpsiMuMu_cff::charge, TrajectoryStateOnSurface::charge(), fastSimProducer_cff::energyLoss, MaterialEffects::energyLossSimulator(), dqmMemoryStats::float, TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), PV3DBase< T, PVType, FrameType >::mag2(), magfield, rawparticle::makeMuon(), RawParticle::momentum(), amptDefaultParameters_cff::mu, fastSimProducer_cff::multipleScattering, MaterialEffects::multipleScatteringSimulator(), MaterialEffects::muonBremsstrahlungSimulator(), BaseParticlePropagator::particle(), position, RawParticle::Px(), RawParticle::Py(), RawParticle::Pz(), RawParticle::setMomentum(), RawParticle::setVertex(), mathSSE::sqrt(), TrajectoryStateOnSurface::surface(), TableParser::table, Surface::tangentPlane(), theMaterialEffects, RawParticle::vertex(), PV3DBase< T, PVType, FrameType >::x(), RawParticle::X(), PV3DBase< T, PVType, FrameType >::y(), RawParticle::Y(), PV3DBase< T, PVType, FrameType >::z(), and RawParticle::Z().

Referenced by produce().

◆ beginRun()

void MuonSimHitProducer::beginRun	(	edm::Run const &	run,
		const edm::EventSetup &	es
	)

overrideprivate

Definition at line 101 of file MuonSimHitProducer.cc.

                                                                             {
   //services
   edm::ESHandle<MagneticField> magField;
   edm::ESHandle<DTGeometry> dtGeometry;
   edm::ESHandle<CSCGeometry> cscGeometry;
   edm::ESHandle<RPCGeometry> rpcGeometry;
   edm::ESHandle<Propagator> propagator;
  
   es.get<IdealMagneticFieldRecord>().get(magField);
   es.get<MuonGeometryRecord>().get("MisAligned", dtGeometry);
   es.get<MuonGeometryRecord>().get("MisAligned", cscGeometry);
   es.get<MuonGeometryRecord>().get(rpcGeometry);
  
   magfield = &(*magField);
   dtGeom = &(*dtGeometry);
   cscGeom = &(*cscGeometry);
   rpcGeom = &(*rpcGeometry);
  
   bool duringEvent = false;
   theService->update(es, duringEvent);
  
   // A few propagators
   propagatorWithMaterial = &(*(theService->propagator("SteppingHelixPropagatorAny")));
   propagatorWithoutMaterial = propagatorWithMaterial->clone();
   SteppingHelixPropagator* SHpropagator = dynamic_cast<SteppingHelixPropagator*>(propagatorWithoutMaterial);  // Beuark!
   SHpropagator->setMaterialMode(true);  // switches OFF material effects;
 }

References Propagator::clone(), cscGeom, dtGeom, edm::EventSetup::get(), get, magfield, TrackCandidateProducer_cfi::propagator, MuonServiceProxy::propagator(), propagatorWithMaterial, propagatorWithoutMaterial, rpcGeom, SteppingHelixPropagator::setMaterialMode(), theService, and MuonServiceProxy::update().

◆ produce()

void MuonSimHitProducer::produce	(	edm::Event &	iEvent,
		const edm::EventSetup &	iSetup
	)

overrideprivate

Definition at line 145 of file MuonSimHitProducer.cc.

                                                                               {
   // using namespace edm;
   // using namespace std;
   edm::ESHandle<HepPDT::ParticleDataTable> pdg;
   iSetup.getData(pdg);
  
   RandomEngineAndDistribution random(iEvent.streamID());
  
   MuonPatternRecoDumper dumper;
  
   edm::Handle<std::vector<SimTrack> > simMuons;
   edm::Handle<std::vector<SimVertex> > simVertices;
   std::vector<PSimHit> theCSCHits;
   std::vector<PSimHit> theDTHits;
   std::vector<PSimHit> theRPCHits;
  
   DirectMuonNavigation navigation(theService->detLayerGeometry());
   iEvent.getByToken(simMuonToken, simMuons);
   iEvent.getByToken(simVertexToken, simVertices);
  
   for (unsigned int itrk = 0; itrk < simMuons->size(); itrk++) {
     const SimTrack& mySimTrack = (*simMuons)[itrk];
     math::XYZTLorentzVector mySimP4(
         mySimTrack.momentum().x(), mySimTrack.momentum().y(), mySimTrack.momentum().z(), mySimTrack.momentum().t());
  
     // Decaying hadrons are now in the list, and so are their muon daughter
     // Ignore the hadrons here.
     int pid = mySimTrack.type();
     if (abs(pid) != 13 && abs(pid) != 1000024)
       continue;
     double t0 = 0;
     GlobalPoint initialPosition;
     int ivert = mySimTrack.vertIndex();
     if (ivert >= 0) {
       t0 = (*simVertices)[ivert].position().t();
       GlobalPoint xyzzy((*simVertices)[ivert].position().x(),
                         (*simVertices)[ivert].position().y(),
                         (*simVertices)[ivert].position().z());
       initialPosition = xyzzy;
     }
     //
     //  Presumably t0 has dimensions of cm if not mm?
     //  Convert to ns for internal calculations.
     //  I wonder where we should get c from?
     //
     double tof = t0 / 29.98;
  
 #ifdef FAMOS_DEBUG
     std::cout << " ===> MuonSimHitProducer::reconstruct() found SIMTRACK - pid = " << pid;
     std::cout << " : pT = " << mySimP4.Pt() << ", eta = " << mySimP4.Eta() << ", phi = " << mySimP4.Phi() << std::endl;
 #endif
  
     //
     //  Produce muons sim hits starting from undecayed simulated muons
     //
  
     GlobalPoint startingPosition(mySimTrack.trackerSurfacePosition().x(),
                                  mySimTrack.trackerSurfacePosition().y(),
                                  mySimTrack.trackerSurfacePosition().z());
     GlobalVector startingMomentum(mySimTrack.trackerSurfaceMomentum().x(),
                                   mySimTrack.trackerSurfaceMomentum().y(),
                                   mySimTrack.trackerSurfaceMomentum().z());
     //
     //  Crap... there's no time-of-flight to the trackerSurfacePosition()...
     //  So, this will be wrong when the curvature can't be neglected, but that
     //  will be rather seldom...  May as well ignore the mass too.
     //
     GlobalVector dtracker = startingPosition - initialPosition;
     tof += dtracker.mag() / 29.98;
  
 #ifdef FAMOS_DEBUG
     std::cout << " the Muon START position " << startingPosition << std::endl;
     std::cout << " the Muon START momentum " << startingMomentum << std::endl;
 #endif
  
     //
     //  Some magic to define a TrajectoryStateOnSurface
     //
     PlaneBuilder pb;
     GlobalVector zAxis = startingMomentum.unit();
     GlobalVector yAxis(zAxis.y(), -zAxis.x(), 0);
     GlobalVector xAxis = yAxis.cross(zAxis);
     Surface::RotationType rot = Surface::RotationType(xAxis, yAxis, zAxis);
     PlaneBuilder::ReturnType startingPlane = pb.plane(startingPosition, rot);
     GlobalTrajectoryParameters gtp(startingPosition, startingMomentum, (int)mySimTrack.charge(), magfield);
     TrajectoryStateOnSurface startingState(gtp, *startingPlane);
  
     std::vector<const DetLayer*> navLayers;
     if (fabs(startingState.globalMomentum().eta()) > 4.5) {
       navLayers = navigation.compatibleEndcapLayers(*(startingState.freeState()), alongMomentum);
     } else {
       navLayers = navigation.compatibleLayers(*(startingState.freeState()), alongMomentum);
     }
     /*
     edm::ESHandle<Propagator> propagator =
       theService->propagator("SteppingHelixPropagatorAny");
     */
  
     if (navLayers.empty())
       continue;
  
 #ifdef FAMOS_DEBUG
     std::cout << "Found " << navLayers.size() << " compatible DetLayers..." << std::endl;
 #endif
  
     TrajectoryStateOnSurface propagatedState = startingState;
     for (unsigned int ilayer = 0; ilayer < navLayers.size(); ilayer++) {
 #ifdef FAMOS_DEBUG
       std::cout << "Propagating to layer " << ilayer << " " << dumper.dumpLayer(navLayers[ilayer]) << std::endl;
 #endif
  
       std::vector<DetWithState> comps =
           navLayers[ilayer]->compatibleDets(propagatedState, *propagatorWithMaterial, theEstimator);
       if (comps.empty())
         continue;
  
 #ifdef FAMOS_DEBUG
       std::cout << "Propagating " << propagatedState << std::endl;
 #endif
  
       // Starting momentum
       double pi = propagatedState.globalMomentum().mag();
  
       // Propagate with material effects (dE/dx average only)
       SteppingHelixStateInfo shsStart(*(propagatedState.freeTrajectoryState()));
       SteppingHelixStateInfo shsDest;
       ((const SteppingHelixPropagator*)propagatorWithMaterial)
           ->propagate(shsStart, navLayers[ilayer]->surface(), shsDest);
       std::pair<TrajectoryStateOnSurface, double> next(shsDest.getStateOnSurface(navLayers[ilayer]->surface()),
                                                        shsDest.path());
       // No need to continue if there is no valid propagation available.
       // This happens rarely (~0.1% of ttbar events)
       if (!next.first.isValid())
         continue;
       // This is the estimate of the number of radiation lengths traversed,
       // together with the total path length
       double radPath = shsDest.radPath();
       double pathLength = next.second;
  
       // Now propagate without dE/dx (average)
       // [To add the dE/dx fluctuations to the actual dE/dx]
       std::pair<TrajectoryStateOnSurface, double> nextNoMaterial =
           propagatorWithoutMaterial->propagateWithPath(propagatedState, navLayers[ilayer]->surface());
  
       // Update the propagated state
       propagatedState = next.first;
       double pf = propagatedState.globalMomentum().mag();
  
       // Insert dE/dx fluctuations and multiple scattering
       // Skip this step if nextNoMaterial.first is not valid
       // This happens rarely (~0.02% of ttbar events)
       if (theMaterialEffects && nextNoMaterial.first.isValid())
         applyMaterialEffects(propagatedState, nextNoMaterial.first, radPath, &random, *pdg);
       // Check that the 'shaken' propagatedState is still valid, otherwise continue
       if (!propagatedState.isValid())
         continue;
       // (No evidence that this ever happens)
       //
       //  Consider this... 1 GeV muon has a velocity that is only 0.5% slower than c...
       //  We probably can safely ignore the mass for anything that makes it out to the
       //  muon chambers.
       //
       double pavg = 0.5 * (pi + pf);
       double m = mySimP4.M();
       double rbeta = sqrt(1 + m * m / (pavg * pavg)) / 29.98;
       double dtof = pathLength * rbeta;
  
 #ifdef FAMOS_DEBUG
       std::cout << "Propagated to next surface... path length = " << pathLength << " cm, dTOF = " << dtof << " ns"
                 << std::endl;
 #endif
  
       tof += dtof;
  
       for (unsigned int icomp = 0; icomp < comps.size(); icomp++) {
         const GeomDet* gd = comps[icomp].first;
         if (gd->subDetector() == GeomDetEnumerators::DT) {
           DTChamberId id(gd->geographicalId());
           const DTChamber* chamber = dtGeom->chamber(id);
           std::vector<const DTSuperLayer*> superlayer = chamber->superLayers();
           for (unsigned int isl = 0; isl < superlayer.size(); isl++) {
             std::vector<const DTLayer*> layer = superlayer[isl]->layers();
             for (unsigned int ilayer = 0; ilayer < layer.size(); ilayer++) {
               DTLayerId lid = layer[ilayer]->id();
 #ifdef FAMOS_DEBUG
               std::cout << "    Extrapolated to DT (" << lid.wheel() << "," << lid.station() << "," << lid.sector()
                         << "," << lid.superlayer() << "," << lid.layer() << ")" << std::endl;
 #endif
  
               const GeomDetUnit* det = dtGeom->idToDetUnit(lid);
  
               HelixArbitraryPlaneCrossing crossing(propagatedState.globalPosition().basicVector(),
                                                    propagatedState.globalMomentum().basicVector(),
                                                    propagatedState.transverseCurvature(),
                                                    anyDirection);
               std::pair<bool, double> path = crossing.pathLength(det->surface());
               if (!path.first)
                 continue;
               LocalPoint lpos = det->toLocal(GlobalPoint(crossing.position(path.second)));
               if (!det->surface().bounds().inside(lpos))
                 continue;
               const DTTopology& dtTopo = layer[ilayer]->specificTopology();
               int wire = dtTopo.channel(lpos);
               if (wire - dtTopo.firstChannel() < 0 || wire - dtTopo.lastChannel() > 0)
                 continue;
               // no drift cell here (on the chamber edge or just outside)
               // this hit would otherwise be discarded downstream in the digitizer
  
               DTWireId wid(lid, wire);
               double thickness = det->surface().bounds().thickness();
               LocalVector lmom = det->toLocal(GlobalVector(crossing.direction(path.second)));
               lmom = lmom.unit() * propagatedState.localMomentum().mag();
  
               // Factor that takes into account the (rec)hits lost because of delta's, etc.:
               // (Not fully satisfactory patch, but it seems to work...)
               double pmu = lmom.mag();
               double theDTHitIneff = pmu > 0 ? exp(kDT * log(pmu) + fDT) : 0.;
               if (random.flatShoot() < theDTHitIneff)
                 continue;
  
               double eloss = 0;
               double pz = fabs(lmom.z());
               LocalPoint entry = lpos - 0.5 * thickness * lmom / pz;
               LocalPoint exit = lpos + 0.5 * thickness * lmom / pz;
               double dtof = path.second * rbeta;
               int trkid = mySimTrack.trackId();
               unsigned int id = wid.rawId();
               short unsigned int processType = 2;
               PSimHit hit(
                   entry, exit, lmom.mag(), tof + dtof, eloss, pid, id, trkid, lmom.theta(), lmom.phi(), processType);
               theDTHits.push_back(hit);
             }
           }
         } else if (gd->subDetector() == GeomDetEnumerators::CSC) {
           CSCDetId id(gd->geographicalId());
           const CSCChamber* chamber = cscGeom->chamber(id);
           std::vector<const CSCLayer*> layers = chamber->layers();
           for (unsigned int ilayer = 0; ilayer < layers.size(); ilayer++) {
             CSCDetId lid = layers[ilayer]->id();
  
 #ifdef FAMOS_DEBUG
             std::cout << "    Extrapolated to CSC (" << lid.endcap() << "," << lid.ring() << "," << lid.station() << ","
                       << lid.layer() << ")" << std::endl;
 #endif
  
             const GeomDetUnit* det = cscGeom->idToDetUnit(lid);
             HelixArbitraryPlaneCrossing crossing(propagatedState.globalPosition().basicVector(),
                                                  propagatedState.globalMomentum().basicVector(),
                                                  propagatedState.transverseCurvature(),
                                                  anyDirection);
             std::pair<bool, double> path = crossing.pathLength(det->surface());
             if (!path.first)
               continue;
             LocalPoint lpos = det->toLocal(GlobalPoint(crossing.position(path.second)));
             // For CSCs the Bounds are for chamber frames not gas regions
             //      if ( ! det->surface().bounds().inside(lpos) ) continue;
             // New function knows where the 'active' volume is:
             const CSCLayerGeometry* laygeom = layers[ilayer]->geometry();
             if (!laygeom->inside(lpos))
               continue;
             //double thickness = laygeom->thickness(); gives number which is about 20 times too big
             double thickness = det->surface().bounds().thickness();  // this one works much better...
             LocalVector lmom = det->toLocal(GlobalVector(crossing.direction(path.second)));
             lmom = lmom.unit() * propagatedState.localMomentum().mag();
  
             // Factor that takes into account the (rec)hits lost because of delta's, etc.:
             // (Not fully satisfactory patch, but it seems to work...)
             double pmu = lmom.mag();
             double theCSCHitIneff = pmu > 0 ? exp(kCSC * log(pmu) + fCSC) : 0.;
             // Take into account the different geometry in ME11:
             if (id.station() == 1 && id.ring() == 1)
               theCSCHitIneff = theCSCHitIneff * 0.442;
             if (random.flatShoot() < theCSCHitIneff)
               continue;
  
             double eloss = 0;
             double pz = fabs(lmom.z());
             LocalPoint entry = lpos - 0.5 * thickness * lmom / pz;
             LocalPoint exit = lpos + 0.5 * thickness * lmom / pz;
             double dtof = path.second * rbeta;
             int trkid = mySimTrack.trackId();
             unsigned int id = lid.rawId();
             short unsigned int processType = 2;
             PSimHit hit(
                 entry, exit, lmom.mag(), tof + dtof, eloss, pid, id, trkid, lmom.theta(), lmom.phi(), processType);
             theCSCHits.push_back(hit);
           }
         } else if (gd->subDetector() == GeomDetEnumerators::RPCBarrel ||
                    gd->subDetector() == GeomDetEnumerators::RPCEndcap) {
           RPCDetId id(gd->geographicalId());
           const RPCChamber* chamber = rpcGeom->chamber(id);
           std::vector<const RPCRoll*> roll = chamber->rolls();
           for (unsigned int iroll = 0; iroll < roll.size(); iroll++) {
             RPCDetId rid = roll[iroll]->id();
  
 #ifdef FAMOS_DEBUG
             std::cout << "    Extrapolated to RPC (" << rid.ring() << "," << rid.station() << "," << rid.sector() << ","
                       << rid.subsector() << "," << rid.layer() << "," << rid.roll() << ")" << std::endl;
 #endif
  
             const GeomDetUnit* det = rpcGeom->idToDetUnit(rid);
             HelixArbitraryPlaneCrossing crossing(propagatedState.globalPosition().basicVector(),
                                                  propagatedState.globalMomentum().basicVector(),
                                                  propagatedState.transverseCurvature(),
                                                  anyDirection);
             std::pair<bool, double> path = crossing.pathLength(det->surface());
             if (!path.first)
               continue;
             LocalPoint lpos = det->toLocal(GlobalPoint(crossing.position(path.second)));
             if (!det->surface().bounds().inside(lpos))
               continue;
             double thickness = det->surface().bounds().thickness();
             LocalVector lmom = det->toLocal(GlobalVector(crossing.direction(path.second)));
             lmom = lmom.unit() * propagatedState.localMomentum().mag();
             double eloss = 0;
             double pz = fabs(lmom.z());
             LocalPoint entry = lpos - 0.5 * thickness * lmom / pz;
             LocalPoint exit = lpos + 0.5 * thickness * lmom / pz;
             double dtof = path.second * rbeta;
             int trkid = mySimTrack.trackId();
             unsigned int id = rid.rawId();
             short unsigned int processType = 2;
             PSimHit hit(
                 entry, exit, lmom.mag(), tof + dtof, eloss, pid, id, trkid, lmom.theta(), lmom.phi(), processType);
             theRPCHits.push_back(hit);
           }
         } else {
           std::cout << "Extrapolated to unknown subdetector '" << gd->subDetector() << "'..." << std::endl;
         }
       }
     }
   }
  
   std::unique_ptr<edm::PSimHitContainer> pcsc(new edm::PSimHitContainer);
   int n = 0;
   for (std::vector<PSimHit>::const_iterator i = theCSCHits.begin(); i != theCSCHits.end(); i++) {
     pcsc->push_back(*i);
     n += 1;
   }
   iEvent.put(std::move(pcsc), "MuonCSCHits");
  
   std::unique_ptr<edm::PSimHitContainer> pdt(new edm::PSimHitContainer);
   n = 0;
   for (std::vector<PSimHit>::const_iterator i = theDTHits.begin(); i != theDTHits.end(); i++) {
     pdt->push_back(*i);
     n += 1;
   }
   iEvent.put(std::move(pdt), "MuonDTHits");
  
   std::unique_ptr<edm::PSimHitContainer> prpc(new edm::PSimHitContainer);
   n = 0;
   for (std::vector<PSimHit>::const_iterator i = theRPCHits.begin(); i != theRPCHits.end(); i++) {
     prpc->push_back(*i);
     n += 1;
   }
   iEvent.put(std::move(prpc), "MuonRPCHits");
 }

References funct::abs(), alongMomentum, anyDirection, applyMaterialEffects(), PV3DBase< T, PVType, FrameType >::basicVector(), Surface::bounds(), relativeConstraints::chamber, RPCGeometry::chamber(), CSCGeometry::chamber(), DTGeometry::chamber(), DTTopology::channel(), CoreSimTrack::charge(), gather_cfg::cout, GeomDetEnumerators::CSC, cscGeom, MuonServiceProxy::detLayerGeometry(), GeomDetEnumerators::DT, dtGeom, CSCDetId::endcap(), mps_splice::entry, beamvalidation::exit(), JetChargeProducer_cfi::exp, fCSC, fDT, DTTopology::firstChannel(), TrajectoryStateOnSurface::freeTrajectoryState(), GeomDet::geographicalId(), edm::EventSetup::getData(), SteppingHelixStateInfo::getStateOnSurface(), TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), mps_fire::i, triggerObjects_cff::id, RPCGeometry::idToDetUnit(), DTGeometry::idToDetUnit(), CSCGeometry::idToDetUnit(), iEvent, Bounds::inside(), CSCLayerGeometry::inside(), TrajectoryStateOnSurface::isValid(), kCSC, kDT, DTTopology::lastChannel(), DTLayerId::layer(), CSCDetId::layer(), hgcalTopologyTester_cfi::layers, TrajectoryStateOnSurface::localMomentum(), dqm-mbProfile::log, visualization-live-secondInstance_cfg::m, PV3DBase< T, PVType, FrameType >::mag(), magfield, CoreSimTrack::momentum(), eostools::move(), dqmiodumpmetadata::n, GetRecoTauVFromDQM_MC_cff::next, castor_dqm_sourceclient_file_cfg::path, SteppingHelixStateInfo::path(), dileptonTrigSettings_cff::pdg, packedPFCandidateRefMixer_cfi::pf, PV3DBase< T, PVType, FrameType >::phi(), pi, PlaneBuilder::plane(), position, Propagator::propagateWithPath(), propagatorWithMaterial, propagatorWithoutMaterial, SteppingHelixStateInfo::radPath(), DetId::rawId(), bril_dqm_clientPB-live_cfg::rid, CSCDetId::ring(), relativeConstraints::ring, makeMuonMisalignmentScenario::rot, GeomDetEnumerators::RPCBarrel, GeomDetEnumerators::RPCEndcap, rpcGeom, DTChamberId::sector(), simMuonToken, simVertexToken, HGCalValidator_cfi::simVertices, mathSSE::sqrt(), DTChamberId::station(), relativeConstraints::station, CSCDetId::station(), GeomDet::subDetector(), DTSuperLayerId::superlayer(), GeomDet::surface(), FrontierCondition_GT_autoExpress_cfi::t0, theEstimator, theMaterialEffects, theService, PV3DBase< T, PVType, FrameType >::theta(), Bounds::thickness(), Calorimetry_cff::thickness, GeomDet::toLocal(), SimTrack::trackerSurfaceMomentum(), SimTrack::trackerSurfacePosition(), CoreSimTrack::trackId(), TrajectoryStateOnSurface::transverseCurvature(), CoreSimTrack::type(), Vector3DBase< T, FrameTag >::unit(), SimTrack::vertIndex(), DTChamberId::wheel(), x, HLT_2018_cff::xAxis, y, HLT_2018_cff::yAxis, z, PV3DBase< T, PVType, FrameType >::z(), and HLT_2018_cff::zAxis.

◆ readParameters()

void MuonSimHitProducer::readParameters	(	const edm::ParameterSet &	fastMuons,
		const edm::ParameterSet &	fastTracks,
		const edm::ParameterSet &	matEff
	)

private

Definition at line 503 of file MuonSimHitProducer.cc.

                                                                        {
   // Muons
   std::string _simModuleLabel = fastMuons.getParameter<std::string>("simModuleLabel");
   std::string _simModuleProcess = fastMuons.getParameter<std::string>("simModuleProcess");
   simMuonLabel = edm::InputTag(_simModuleLabel, _simModuleProcess);
   simVertexLabel = edm::InputTag(_simModuleLabel);
  
   std::vector<double> simHitIneffDT = fastMuons.getParameter<std::vector<double> >("simHitDTIneffParameters");
   std::vector<double> simHitIneffCSC = fastMuons.getParameter<std::vector<double> >("simHitCSCIneffParameters");
   kDT = simHitIneffDT[0];
   fDT = simHitIneffDT[1];
   kCSC = simHitIneffCSC[0];
   fCSC = simHitIneffCSC[1];
  
   // Tracks
   fullPattern_ = fastTracks.getUntrackedParameter<bool>("FullPatternRecognition");
  
   // The following should be on LogInfo
   //  std::cout << " MUON SIM HITS: FastSimulation parameters " << std::endl;
   //  std::cout << " ============================================== " << std::endl;
   //  if ( fullPattern_ )
   //    std::cout << " The FULL pattern recognition option is turned ON" << std::endl;
   //  else
   //    std::cout << " The FAST tracking option is turned ON" << std::endl;
  
   // Material Effects
   theMaterialEffects = nullptr;
   if (matEff.getParameter<bool>("PairProduction") || matEff.getParameter<bool>("Bremsstrahlung") ||
       matEff.getParameter<bool>("MuonBremsstrahlung") || matEff.getParameter<bool>("EnergyLoss") ||
       matEff.getParameter<bool>("MultipleScattering"))
     theMaterialEffects = new MaterialEffects(matEff);
 }

References fCSC, fDT, fullPattern_, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), HLT_2018_cff::InputTag, kCSC, kDT, TrajectoryFactories_cff::MaterialEffects, simMuonLabel, simVertexLabel, AlCaHLTBitMon_QueryRunRegistry::string, and theMaterialEffects.

Referenced by MuonSimHitProducer().

Member Data Documentation

◆ cscGeom

const CSCGeometry* MuonSimHitProducer::cscGeom

private

Definition at line 69 of file MuonSimHitProducer.h.

Referenced by beginRun(), and produce().

◆ doGL_

bool MuonSimHitProducer::doGL_

private

Definition at line 95 of file MuonSimHitProducer.h.

◆ doL1_

bool MuonSimHitProducer::doL1_

private

Definition at line 95 of file MuonSimHitProducer.h.

◆ doL3_

bool MuonSimHitProducer::doL3_

private

Definition at line 95 of file MuonSimHitProducer.h.

◆ dtGeom

const DTGeometry* MuonSimHitProducer::dtGeom

private

Definition at line 68 of file MuonSimHitProducer.h.

Referenced by beginRun(), and produce().

◆ fCSC

double MuonSimHitProducer::fCSC

private

Definition at line 84 of file MuonSimHitProducer.h.

Referenced by produce(), and readParameters().

◆ fDT

double MuonSimHitProducer::fDT

private

Definition at line 82 of file MuonSimHitProducer.h.

Referenced by produce(), and readParameters().

◆ fullPattern_

bool MuonSimHitProducer::fullPattern_

private

Definition at line 94 of file MuonSimHitProducer.h.

Referenced by readParameters().

◆ kCSC

double MuonSimHitProducer::kCSC

private

Definition at line 83 of file MuonSimHitProducer.h.

Referenced by produce(), and readParameters().

◆ kDT

double MuonSimHitProducer::kDT

private

Definition at line 81 of file MuonSimHitProducer.h.

Referenced by produce(), and readParameters().

◆ magfield

const MagneticField* MuonSimHitProducer::magfield

private

Definition at line 67 of file MuonSimHitProducer.h.

Referenced by applyMaterialEffects(), beginRun(), and produce().

◆ propagatorWithMaterial

const Propagator* MuonSimHitProducer::propagatorWithMaterial

private

Definition at line 71 of file MuonSimHitProducer.h.

Referenced by beginRun(), and produce().

◆ propagatorWithoutMaterial

Propagator* MuonSimHitProducer::propagatorWithoutMaterial

private

Definition at line 72 of file MuonSimHitProducer.h.

Referenced by beginRun(), produce(), and ~MuonSimHitProducer().

◆ rpcGeom

const RPCGeometry* MuonSimHitProducer::rpcGeom

private

Definition at line 70 of file MuonSimHitProducer.h.

Referenced by beginRun(), and produce().

◆ simMuonLabel

edm::InputTag MuonSimHitProducer::simMuonLabel

private

Definition at line 98 of file MuonSimHitProducer.h.

Referenced by MuonSimHitProducer(), and readParameters().

◆ simMuonToken

edm::EDGetTokenT<std::vector<SimTrack> > MuonSimHitProducer::simMuonToken

private

Definition at line 102 of file MuonSimHitProducer.h.

Referenced by MuonSimHitProducer(), and produce().

◆ simVertexLabel

edm::InputTag MuonSimHitProducer::simVertexLabel

private

Definition at line 99 of file MuonSimHitProducer.h.

Referenced by MuonSimHitProducer(), and readParameters().

◆ simVertexToken

edm::EDGetTokenT<std::vector<SimVertex> > MuonSimHitProducer::simVertexToken

private

Definition at line 103 of file MuonSimHitProducer.h.

Referenced by MuonSimHitProducer(), and produce().

◆ theEstimator

Chi2MeasurementEstimator MuonSimHitProducer::theEstimator

private

Definition at line 65 of file MuonSimHitProducer.h.

Referenced by produce().

◆ theMaterialEffects

MaterialEffects* MuonSimHitProducer::theMaterialEffects

private

Definition at line 74 of file MuonSimHitProducer.h.

Referenced by applyMaterialEffects(), produce(), readParameters(), and ~MuonSimHitProducer().

◆ theService

MuonServiceProxy* MuonSimHitProducer::theService

private

Definition at line 64 of file MuonSimHitProducer.h.

Referenced by beginRun(), MuonSimHitProducer(), and produce().

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ MuonSimHitProducer()

◆ ~MuonSimHitProducer()

Member Function Documentation

◆ applyMaterialEffects()

◆ beginRun()

◆ produce()

◆ readParameters()

Member Data Documentation

◆ cscGeom

◆ doGL_

◆ doL1_

◆ doL3_

◆ dtGeom

◆ fCSC

◆ fDT

◆ fullPattern_

◆ kCSC

◆ kDT

◆ magfield

◆ propagatorWithMaterial

◆ propagatorWithoutMaterial

◆ rpcGeom

◆ simMuonLabel

◆ simMuonToken

◆ simVertexLabel

◆ simVertexToken

◆ theEstimator

◆ theMaterialEffects

◆ theService