Inheritance diagram for MuonSimHitProducer:

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

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

	EDProducer (const EDProducer &)=delete

bool	hasAbilityToProduceInBeginLumis () const final

bool	hasAbilityToProduceInBeginProcessBlocks () const final

bool	hasAbilityToProduceInBeginRuns () const final

bool	hasAbilityToProduceInEndLumis () const final

bool	hasAbilityToProduceInEndProcessBlocks () const final

bool	hasAbilityToProduceInEndRuns () const final

const EDProducer &	operator= (const EDProducer &)=delete

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

const edm::ESGetToken< CSCGeometry, MuonGeometryRecord >	cscGeometryESToken_

bool	doGL_

bool	doL1_

bool	doL3_

const DTGeometry *	dtGeom

const edm::ESGetToken< DTGeometry, MuonGeometryRecord >	dtGeometryESToken_

double	fCSC

double	fDT

bool	fullPattern_

double	kCSC

double	kDT

const MagneticField *	magfield

const edm::ESGetToken< MagneticField, IdealMagneticFieldRecord >	magneticFieldESToken_

const edm::ESGetToken< HepPDT::ParticleDataTable, edm::DefaultRecord >	particleDataTableESToken_

const Propagator *	propagatorWithMaterial

std::unique_ptr< Propagator >	propagatorWithoutMaterial

const RPCGeometry *	rpcGeom

const edm::ESGetToken< RPCGeometry, MuonGeometryRecord >	rpcGeometryESToken_

edm::InputTag	simMuonLabel

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

edm::InputTag	simVertexLabel

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

Chi2MeasurementEstimator	theEstimator

std::unique_ptr< MaterialEffects >	theMaterialEffects

MuonServiceProxy *	theService

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

using	GlobalCache = typename CacheTypes::GlobalCache

using	HasAbility = AbilityChecker< T... >

using	InputProcessBlockCache = typename CacheTypes::InputProcessBlockCache

using	LuminosityBlockCache = typename CacheTypes::LuminosityBlockCache

using	LuminosityBlockContext = LuminosityBlockContextT< LuminosityBlockCache, RunCache, GlobalCache >

using	LuminosityBlockSummaryCache = typename CacheTypes::LuminosityBlockSummaryCache

using	RunCache = typename CacheTypes::RunCache

using	RunContext = RunContextT< RunCache, GlobalCache >

using	RunSummaryCache = typename CacheTypes::RunSummaryCache

Detailed Description

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

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

Definition at line 51 of file MuonSimHitProducer.cc.

Constructor & Destructor Documentation

◆ MuonSimHitProducer()

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

explicit

Definition at line 110 of file MuonSimHitProducer.cc.

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

     : theEstimator(iConfig.getParameter<double>("Chi2EstimatorCut")),
       propagatorWithoutMaterial(nullptr),
       magneticFieldESToken_(esConsumes<edm::Transition::BeginRun>()),
       dtGeometryESToken_(esConsumes<edm::Transition::BeginRun>(edm::ESInputTag("", "MisAligned"))),
       cscGeometryESToken_(esConsumes<edm::Transition::BeginRun>(edm::ESInputTag("", "MisAligned"))),
       rpcGeometryESToken_(esConsumes<edm::Transition::BeginRun>()),
       particleDataTableESToken_(esConsumes()) {
   // 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);
 }

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 547 of file MuonSimHitProducer.cc.

References fastSimProducer_cff::bremsstrahlung, ALCARECOTkAlJpsiMuMu_cff::charge, TrajectoryStateOnSurface::charge(), fastSimProducer_cff::energyLoss, dqmMemoryStats::float, TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), PV3DBase< T, PVType, FrameType >::mag2(), magfield, rawparticle::makeMuon(), RawParticle::momentum(), amptDefaultParameters_cff::mu, fastSimProducer_cff::multipleScattering, normal, 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().

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

◆ beginRun()

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

overrideprivate

Definition at line 139 of file MuonSimHitProducer.cc.

References Propagator::clone(), cscGeom, cscGeometryESToken_, dtGeom, dtGeometryESToken_, edm::EventSetup::getData(), magfield, magneticFieldESToken_, MuonServiceProxy::propagator(), propagatorWithMaterial, propagatorWithoutMaterial, rpcGeom, rpcGeometryESToken_, SteppingHelixPropagator::setMaterialMode(), theService, and MuonServiceProxy::update().

                                                                             {
   //services
   magfield = &es.getData(magneticFieldESToken_);
   dtGeom = &es.getData(dtGeometryESToken_);
   cscGeom = &es.getData(cscGeometryESToken_);
   rpcGeom = &es.getData(rpcGeometryESToken_);
 
   bool duringEvent = false;
   theService->update(es, duringEvent);
 
   // A few propagators
   propagatorWithMaterial = &(*(theService->propagator("SteppingHelixPropagatorAny")));
   propagatorWithoutMaterial.reset(propagatorWithMaterial->clone());
   SteppingHelixPropagator* SHpropagator =
       dynamic_cast<SteppingHelixPropagator*>(propagatorWithoutMaterial.get());  // Beuark!
   SHpropagator->setMaterialMode(true);                                          // switches OFF material effects;
 }

◆ produce()

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

overrideprivate

Definition at line 163 of file MuonSimHitProducer.cc.

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::getHandle(), SteppingHelixStateInfo::getStateOnSurface(), TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), mps_fire::i, l1ctLayer2EG_cff::id, RPCGeometry::idToDetUnit(), DTGeometry::idToDetUnit(), CSCGeometry::idToDetUnit(), iEvent, Bounds::inside(), CSCLayerGeometry::inside(), TrajectoryStateOnSurface::isValid(), kCSC, kDT, DTTopology::lastChannel(), DTLayerId::layer(), CSCDetId::layer(), pixelTopology::layer, hgcalTBTopologyTester_cfi::layers, TrajectoryStateOnSurface::localMomentum(), dqm-mbProfile::log, visualization-live-secondInstance_cfg::m, PV3DBase< T, PVType, FrameType >::mag(), magfield, CoreSimTrack::momentum(), eostools::move(), GetRecoTauVFromDQM_MC_cff::next, particleDataTableESToken_, castor_dqm_sourceclient_file_cfg::path, SteppingHelixStateInfo::path(), packedPFCandidateRefMixer_cfi::pf, PV3DBase< T, PVType, FrameType >::phi(), pi, PlaneBuilder::plane(), position, 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_2023v12_cff::xAxis, y, HLT_2023v12_cff::yAxis, z, PV3DBase< T, PVType, FrameType >::z(), and HLT_2023v12_cff::zAxis.

                                                                               {
   auto pdg = iSetup.getHandle(particleDataTableESToken_);
 
   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);
   for (std::vector<PSimHit>::const_iterator i = theCSCHits.begin(); i != theCSCHits.end(); i++) {
     pcsc->push_back(*i);
   }
   iEvent.put(std::move(pcsc), "MuonCSCHits");
 
   std::unique_ptr<edm::PSimHitContainer> pdt(new edm::PSimHitContainer);
   for (std::vector<PSimHit>::const_iterator i = theDTHits.begin(); i != theDTHits.end(); i++) {
     pdt->push_back(*i);
   }
   iEvent.put(std::move(pdt), "MuonDTHits");
 
   std::unique_ptr<edm::PSimHitContainer> prpc(new edm::PSimHitContainer);
   for (std::vector<PSimHit>::const_iterator i = theRPCHits.begin(); i != theRPCHits.end(); i++) {
     prpc->push_back(*i);
   }
   iEvent.put(std::move(prpc), "MuonRPCHits");
 }

◆ readParameters()

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

private

Definition at line 512 of file MuonSimHitProducer.cc.

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

Referenced by MuonSimHitProducer().

                                                                        {
   // 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 = std::make_unique<MaterialEffects>(matEff);
 }

Member Data Documentation

◆ cscGeom

const CSCGeometry* MuonSimHitProducer::cscGeom

private

Definition at line 61 of file MuonSimHitProducer.cc.

Referenced by beginRun(), and produce().

◆ cscGeometryESToken_

const edm::ESGetToken<CSCGeometry, MuonGeometryRecord> MuonSimHitProducer::cscGeometryESToken_

private

Definition at line 99 of file MuonSimHitProducer.cc.

Referenced by beginRun().

◆ doGL_

bool MuonSimHitProducer::doGL_

private

Definition at line 87 of file MuonSimHitProducer.cc.

◆ doL1_

bool MuonSimHitProducer::doL1_

private

Definition at line 87 of file MuonSimHitProducer.cc.

◆ doL3_

bool MuonSimHitProducer::doL3_

private

Definition at line 87 of file MuonSimHitProducer.cc.

◆ dtGeom

const DTGeometry* MuonSimHitProducer::dtGeom

private

Definition at line 60 of file MuonSimHitProducer.cc.

Referenced by beginRun(), and produce().

◆ dtGeometryESToken_

const edm::ESGetToken<DTGeometry, MuonGeometryRecord> MuonSimHitProducer::dtGeometryESToken_

private

Definition at line 98 of file MuonSimHitProducer.cc.

Referenced by beginRun().

◆ fCSC

double MuonSimHitProducer::fCSC

private

Definition at line 76 of file MuonSimHitProducer.cc.

Referenced by produce(), and readParameters().

◆ fDT

double MuonSimHitProducer::fDT

private

Definition at line 74 of file MuonSimHitProducer.cc.

Referenced by produce(), and readParameters().

◆ fullPattern_

bool MuonSimHitProducer::fullPattern_

private

Definition at line 86 of file MuonSimHitProducer.cc.

Referenced by readParameters().

◆ kCSC

double MuonSimHitProducer::kCSC

private

Definition at line 75 of file MuonSimHitProducer.cc.

Referenced by produce(), and readParameters().

◆ kDT

double MuonSimHitProducer::kDT

private

Definition at line 73 of file MuonSimHitProducer.cc.

Referenced by produce(), and readParameters().

◆ magfield

const MagneticField* MuonSimHitProducer::magfield

private

Definition at line 59 of file MuonSimHitProducer.cc.

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

◆ magneticFieldESToken_

const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> MuonSimHitProducer::magneticFieldESToken_

private

Definition at line 97 of file MuonSimHitProducer.cc.

Referenced by beginRun().

◆ particleDataTableESToken_

const edm::ESGetToken<HepPDT::ParticleDataTable, edm::DefaultRecord> MuonSimHitProducer::particleDataTableESToken_

private

Definition at line 101 of file MuonSimHitProducer.cc.

Referenced by produce().

◆ propagatorWithMaterial

const Propagator* MuonSimHitProducer::propagatorWithMaterial

private

Definition at line 63 of file MuonSimHitProducer.cc.

Referenced by beginRun(), and produce().

◆ propagatorWithoutMaterial

std::unique_ptr<Propagator> MuonSimHitProducer::propagatorWithoutMaterial

private

Definition at line 64 of file MuonSimHitProducer.cc.

Referenced by beginRun(), and produce().

◆ rpcGeom

const RPCGeometry* MuonSimHitProducer::rpcGeom

private

Definition at line 62 of file MuonSimHitProducer.cc.

Referenced by beginRun(), and produce().

◆ rpcGeometryESToken_

const edm::ESGetToken<RPCGeometry, MuonGeometryRecord> MuonSimHitProducer::rpcGeometryESToken_

private

Definition at line 100 of file MuonSimHitProducer.cc.

Referenced by beginRun().

◆ simMuonLabel

edm::InputTag MuonSimHitProducer::simMuonLabel

private

Definition at line 90 of file MuonSimHitProducer.cc.

Referenced by MuonSimHitProducer(), and readParameters().

◆ simMuonToken

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

private

Definition at line 94 of file MuonSimHitProducer.cc.

Referenced by MuonSimHitProducer(), and produce().

◆ simVertexLabel

edm::InputTag MuonSimHitProducer::simVertexLabel

private

Definition at line 91 of file MuonSimHitProducer.cc.

Referenced by MuonSimHitProducer(), and readParameters().

◆ simVertexToken

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

private

Definition at line 95 of file MuonSimHitProducer.cc.

Referenced by MuonSimHitProducer(), and produce().

◆ theEstimator

Chi2MeasurementEstimator MuonSimHitProducer::theEstimator

private

Definition at line 57 of file MuonSimHitProducer.cc.

Referenced by produce().

◆ theMaterialEffects

std::unique_ptr<MaterialEffects> MuonSimHitProducer::theMaterialEffects

private

Definition at line 66 of file MuonSimHitProducer.cc.

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

◆ theService

MuonServiceProxy* MuonSimHitProducer::theService

private

Definition at line 56 of file MuonSimHitProducer.cc.

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

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ MuonSimHitProducer()

Member Function Documentation

◆ applyMaterialEffects()

◆ beginRun()

◆ produce()

◆ readParameters()

Member Data Documentation

◆ cscGeom

◆ cscGeometryESToken_

◆ doGL_

◆ doL1_

◆ doL3_

◆ dtGeom

◆ dtGeometryESToken_

◆ fCSC

◆ fDT

◆ fullPattern_

◆ kCSC

◆ kDT

◆ magfield

◆ magneticFieldESToken_

◆ particleDataTableESToken_

◆ propagatorWithMaterial

◆ propagatorWithoutMaterial

◆ rpcGeom

◆ rpcGeometryESToken_

◆ simMuonLabel

◆ simMuonToken

◆ simVertexLabel

◆ simVertexToken

◆ theEstimator

◆ theMaterialEffects

◆ theService