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

Inheritance diagram for MuonSimHitProducer:

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

	~MuonSimHitProducer ()

Public Member Functions inherited from edm::EDProducer
	EDProducer ()

virtual	~EDProducer ()

Public Member Functions inherited from edm::ProducerBase
	ProducerBase ()

void	registerProducts (ProducerBase , ProductRegistry , ModuleDescription const &)

boost::function< void(const BranchDescription &)>	registrationCallback () const
	used by the fwk to register list of products More...

virtual	~ProducerBase ()

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

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

virtual void	endJob ()

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

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 RandomEngine *	random

const RPCGeometry *	rpcGeom

Chi2MeasurementEstimator	theEstimator

MaterialEffects *	theMaterialEffects

MuonServiceProxy *	theService

std::string	theSimModuleLabel_

std::string	theSimModuleProcess_

std::string	theTrkModuleLabel_

Additional Inherited Members
Public Types inherited from edm::EDProducer
typedef EDProducer	ModuleType

typedef WorkerT< EDProducer >	WorkerType

Public Types inherited from edm::ProducerBase
typedef ProductRegistryHelper::TypeLabelList	TypeLabelList

Static Public Member Functions inherited from edm::EDProducer
static const std::string &	baseType ()

static void	fillDescriptions (ConfigurationDescriptions &descriptions)

static void	prevalidate (ConfigurationDescriptions &descriptions)

Protected Member Functions inherited from edm::EDProducer
CurrentProcessingContext const *	currentContext () const

Protected Member Functions inherited from edm::ProducerBase
template<class TProducer , class TMethod >
void	callWhenNewProductsRegistered (TProducer *iProd, TMethod iMethod)

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 57 of file MuonSimHitProducer.h.

Constructor & Destructor Documentation

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

explicit

Definition at line 80 of file MuonSimHitProducer.cc.

References edm::hlt::Exception, edm::ParameterSet::getParameter(), edm::Service< T >::isAvailable(), MuonServiceProxy_cff::MuonServiceProxy, random, readParameters(), and theService.

                                                                     :
   theEstimator(iConfig.getParameter<double>("Chi2EstimatorCut")),
   propagatorWithoutMaterial(0)
  {
 
   //
   //  Initialize the random number generator service
   //
   edm::Service<edm::RandomNumberGenerator> rng;
   if ( ! rng.isAvailable() ) {
     throw cms::Exception("Configuration") <<
       "MuonSimHitProducer requires the RandomGeneratorService \n"
       "which is not present in the configuration file. \n"
       "You must add the service in the configuration file\n"
       "or remove the module that requires it.";
   }
 
   random = new RandomEngine(&(*rng));
 
   // 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);
 }

MuonSimHitProducer::~MuonSimHitProducer ( )

Definition at line 147 of file MuonSimHitProducer.cc.

References propagatorWithoutMaterial, random, and theMaterialEffects.

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

Member Function Documentation

void MuonSimHitProducer::applyMaterialEffects	(	TrajectoryStateOnSurface &	tsosWithdEdx,
		TrajectoryStateOnSurface &	tsos,
		double	radPath
	)

private

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

Definition at line 588 of file MuonSimHitProducer.cc.

References DeDxDiscriminatorTools::charge(), TrajectoryStateOnSurface::charge(), EnergyLossSimulator::deltaMom(), MaterialEffects::energyLossSimulator(), TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), PV3DBase< T, PVType, FrameType >::mag2(), magfield, RawParticle::momentum(), EnergyLossSimulator::mostLikelyLoss(), RPCpg::mu, MaterialEffects::multipleScatteringSimulator(), MaterialEffects::muonBremsstrahlungSimulator(), position, RawParticle::setID(), MaterialEffectsSimulator::setNormalVector(), RawParticle::setVertex(), mathSSE::sqrt(), TrajectoryStateOnSurface::surface(), Surface::tangentPlane(), theMaterialEffects, MaterialEffectsSimulator::updateState(), 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(momentum,position,charge,0);
   theMuon.setID(-(int)charge*13);
 
   // 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);
 
     // 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.momentum() + energyLoss->deltaMom() + fac * deltaMom;
     XYZTLorentzVector theNewPosition = theMuon.vertex() + fac * deltaPos;
     fac  = (theNewMomentum.E()*theNewMomentum.E()-mu*mu)/theNewMomentum.Vect().Mag2();
     fac  = fac>0.? std::sqrt(fac) : 1E-9;
     theMuon.SetXYZT(theNewMomentum.Px()*fac,theNewMomentum.Py()*fac,
                     theNewMomentum.Pz()*fac,theNewMomentum.E());    
     theMuon.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);
   }
   
   // Muon Bremsstrahlung
   if ( bremsstrahlung ) {
     // Compute the amount of Muon Bremsstrahlung after given path length
     bremsstrahlung->updateState(theMuon,radPath);
   }
 
 
   // Fill the propagated state
   GlobalPoint propagatedPosition(theMuon.X(),theMuon.Y(),theMuon.Z());
   GlobalVector propagatedMomentum(theMuon.Px(),theMuon.Py(),theMuon.Pz());
   GlobalTrajectoryParameters propagatedGtp(propagatedPosition,propagatedMomentum,(int)charge,magfield);
   tsosWithdEdx = TrajectoryStateOnSurface(propagatedGtp,nextSurface);
 
 }

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

privatevirtual

Reimplemented from edm::EDProducer.

Definition at line 118 of file MuonSimHitProducer.cc.

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

                                                                    {
 
   //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);
 
   theService->update(es);
 
   // A few propagators
   propagatorWithMaterial = &(*(theService->propagator("SteppingHelixPropagatorAny")));
   propagatorWithoutMaterial = propagatorWithMaterial->clone();
   SteppingHelixPropagator* SHpropagator = dynamic_cast<SteppingHelixPropagator*> (propagatorWithoutMaterial); // Beuark!
   SHpropagator->setMaterialMode(true); // switches OFF material effects;
 
 }

void MuonSimHitProducer::endJob ( void )

privatevirtual

Reimplemented from edm::EDProducer.

Definition at line 545 of file MuonSimHitProducer.cc.

546 {

547 }

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

privatevirtual

Implements edm::EDProducer.

Definition at line 168 of file MuonSimHitProducer.cc.

References abs, alongMomentum, anyDirection, applyMaterialEffects(), PV3DBase< T, PVType, FrameType >::basicVector(), BoundSurface::bounds(), RPCGeometry::chamber(), CSCGeometry::chamber(), DTGeometry::chamber(), DTTopology::channel(), CoreSimTrack::charge(), gather_cfg::cout, Vector3DBase< T, FrameTag >::cross(), GeomDetEnumerators::CSC, cscGeom, MuonServiceProxy::detLayerGeometry(), GeomDetEnumerators::DT, dtGeom, MuonPatternRecoDumper::dumpLayer(), CSCDetId::endcap(), cmsRelvalreport::exit, create_public_lumi_plots::exp, fCSC, fDT, DTTopology::firstChannel(), RandomEngine::flatShoot(), TrajectoryStateOnSurface::freeTrajectoryState(), GeomDet::geographicalId(), edm::Event::getByLabel(), SteppingHelixStateInfo::getStateOnSurface(), TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), i, errorMatrix2Lands_multiChannel::id, RPCGeometry::idToDetUnit(), DTGeometry::idToDetUnit(), CSCGeometry::idToDetUnit(), Bounds::inside(), CSCLayerGeometry::inside(), TrajectoryStateOnSurface::isValid(), kCSC, kDT, DTTopology::lastChannel(), DTLayerId::layer(), CSCDetId::layer(), RPCDetId::layer(), CSCChamber::layers(), TrajectoryStateOnSurface::localMomentum(), create_public_lumi_plots::log, m, PV3DBase< T, PVType, FrameType >::mag(), magfield, CoreSimTrack::momentum(), n, GetRecoTauVFromDQM_MC_cff::next, scaleCards::path, SteppingHelixStateInfo::path(), PV3DBase< T, PVType, FrameType >::phi(), pi, evf::utils::pid, PlaneBuilder::plane(), position, Propagator::propagateWithPath(), propagatorWithMaterial, propagatorWithoutMaterial, edm::Event::put(), SteppingHelixStateInfo::radPath(), random, DetId::rawId(), relativeConstraints::ring, RPCDetId::ring(), CSCDetId::ring(), RPCDetId::roll(), RPCChamber::rolls(), makeMuonMisalignmentScenario::rot, GeomDetEnumerators::RPCBarrel, GeomDetEnumerators::RPCEndcap, rpcGeom, DTChamberId::sector(), RPCDetId::sector(), mathSSE::sqrt(), DTChamberId::station(), relativeConstraints::station, CSCDetId::station(), RPCDetId::station(), GeomDet::subDetector(), RPCDetId::subsector(), DTSuperLayerId::superlayer(), DTChamber::superLayers(), GeomDet::surface(), theEstimator, theMaterialEffects, theService, theSimModuleLabel_, theSimModuleProcess_, PV3DBase< T, PVType, FrameType >::theta(), Bounds::thickness(), GeomDet::toLocal(), SimTrack::trackerSurfaceMomentum(), SimTrack::trackerSurfacePosition(), CoreSimTrack::trackId(), TrajectoryStateOnSurface::transverseCurvature(), CoreSimTrack::type(), Vector3DBase< T, FrameTag >::unit(), SimTrack::vertIndex(), DTChamberId::wheel(), x, PV3DBase< T, PVType, FrameType >::x(), detailsBasic3DVector::y, PV3DBase< T, PVType, FrameType >::y(), detailsBasic3DVector::z, and PV3DBase< T, PVType, FrameType >::z().

Referenced by JSONExport.JsonExport::export(), HTMLExport.HTMLExport::export(), and HTMLExport.HTMLExportStatic::export().

                                                                         {
   // using namespace edm;
   // using namespace std;
 
   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.getByLabel(theSimModuleLabel_,theSimModuleProcess_,simMuons);
   iEvent.getByLabel(theSimModuleLabel_,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 ) 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()));
       const SteppingHelixStateInfo& shsDest = 
     ((const SteppingHelixPropagator*)propagatorWithMaterial)->propagate(shsStart,navLayers[ilayer]->surface());
       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);
       // 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::auto_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(pcsc,"MuonCSCHits");
 
   std::auto_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(pdt,"MuonDTHits");
 
   std::auto_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(prpc,"MuonRPCHits");
 
 }

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

private

Definition at line 551 of file MuonSimHitProducer.cc.

References fCSC, fDT, fullPattern_, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), kCSC, kDT, random, theMaterialEffects, theSimModuleLabel_, theSimModuleProcess_, and theTrkModuleLabel_.

Referenced by MuonSimHitProducer().

                                                   {
   // Muons
   theSimModuleLabel_   = fastMuons.getParameter<std::string>("simModuleLabel");
   theSimModuleProcess_ = fastMuons.getParameter<std::string>("simModuleProcess");
   theTrkModuleLabel_   = fastMuons.getParameter<std::string>("trackModuleLabel");
   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 = 0;
   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,random);
 
 }