MuonSimHitProducer Class Reference

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
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
std::unique_ptr< 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
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 104 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.

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 555 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, dqmMemoryStats::float, TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), PV3DBase< T, PVType, FrameType >::mag2(), magfield, rawparticle::makeMuon(), RawParticle::momentum(), amptDefaultParameters_cff::mu, fastSimProducer_cff::multipleScattering, 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 127 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.reset(propagatorWithMaterial->clone());
  SteppingHelixPropagator* SHpropagator =
      dynamic_cast<SteppingHelixPropagator*>(propagatorWithoutMaterial.get());  // 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 162 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(), phase1PixelTopology::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, 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_FULL_cff::xAxis, y, HLT_FULL_cff::yAxis, z, PV3DBase< T, PVType, FrameType >::z(), and HLT_FULL_cff::zAxis.

readParameters()

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

private

Definition at line 520 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 = std::make_unique<MaterialEffects>(matEff);
}

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

Referenced by MuonSimHitProducer().

Member Data Documentation

cscGeom

const CSCGeometry* MuonSimHitProducer::cscGeom

private

Definition at line 61 of file MuonSimHitProducer.cc.

Referenced by beginRun(), and produce().

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().

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().

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().

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().