|
|
|
#include <FastSimulation/MuonSimHitProducer/src/MuonSimHitProducer.cc>
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.
| MuonSimHitProducer::MuonSimHitProducer | ( | const edm::ParameterSet & | iConfig | ) | [explicit] |
Definition at line 79 of file MuonSimHitProducer.cc.
References Exception, edm::ParameterSet::getParameter(), insideOut, edm::Service< T >::isAvailable(), MuonServiceProxy_cff::MuonServiceProxy, random, readParameters(), theService, and theUpdator.
00079 { 00080 00081 // 00082 // Initialize the random number generator service 00083 // 00084 edm::Service<edm::RandomNumberGenerator> rng; 00085 if ( ! rng.isAvailable() ) { 00086 throw cms::Exception("Configuration") << 00087 "MuonSimHitProducer requires the RandomGeneratorService \n" 00088 "which is not present in the configuration file. \n" 00089 "You must add the service in the configuration file\n" 00090 "or remove the module that requires it."; 00091 } 00092 00093 random = new RandomEngine(&(*rng)); 00094 00095 // Read relevant parameters 00096 readParameters(iConfig.getParameter<edm::ParameterSet>("MUONS"), 00097 iConfig.getParameter<edm::ParameterSet>("TRACKS"), 00098 iConfig.getParameter<edm::ParameterSet>("MaterialEffectsForMuons")); 00099 00100 // 00101 // register your products ... need to declare at least one possible product... 00102 // 00103 produces<edm::PSimHitContainer>("MuonCSCHits"); 00104 produces<edm::PSimHitContainer>("MuonDTHits"); 00105 produces<edm::PSimHitContainer>("MuonRPCHits"); 00106 00107 edm::ParameterSet serviceParameters = 00108 iConfig.getParameter<edm::ParameterSet>("ServiceParameters"); 00109 theService = new MuonServiceProxy(serviceParameters); 00110 edm::ParameterSet updatorParameters = 00111 iConfig.getParameter<edm::ParameterSet>("MuonTrajectoryUpdatorParameters"); 00112 theUpdator = new MuonTrajectoryUpdator(updatorParameters,insideOut); 00113 }
| MuonSimHitProducer::~MuonSimHitProducer | ( | ) |
Definition at line 146 of file MuonSimHitProducer.cc.
References random, and theMaterialEffects.
00147 { 00148 // do anything here that needs to be done at destruction time 00149 // (e.g. close files, deallocate resources etc.) 00150 00151 if ( random ) { 00152 delete random; 00153 } 00154 00155 if ( theMaterialEffects ) delete theMaterialEffects; 00156 }
| void MuonSimHitProducer::applyMaterialEffects | ( | TrajectoryStateOnSurface & | tsosWithdEdx, | |
| TrajectoryStateOnSurface & | tsos, | |||
| double | radPath | |||
| ) | [private] |
Simulate material effects in iron (dE/dx, multiple scattering).
Definition at line 567 of file MuonSimHitProducer.cc.
References TrajectoryStateOnSurface::charge(), EnergyLossSimulator::deltaMom(), MaterialEffects::energyLossSimulator(), TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), PV3DBase< T, PVType, FrameType >::mag2(), magfield, RawParticle::momentum(), EnergyLossSimulator::mostLikelyLoss(), MaterialEffects::multipleScatteringSimulator(), RawParticle::setID(), MaterialEffectsSimulator::setNormalVector(), RawParticle::setVertex(), funct::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().
00569 { 00570 00571 // The energy loss simulator 00572 EnergyLossSimulator* energyLoss = theMaterialEffects->energyLossSimulator(); 00573 00574 // The multiple scattering simulator 00575 MultipleScatteringSimulator* multipleScattering = theMaterialEffects->multipleScatteringSimulator(); 00576 00577 // Initialize the Particle position, momentum and energy 00578 const Surface& nextSurface = tsos.surface(); 00579 GlobalPoint gPos = energyLoss ? tsos.globalPosition() : tsosWithdEdx.globalPosition(); 00580 GlobalVector gMom = energyLoss ? tsos.globalMomentum() : tsosWithdEdx.globalMomentum(); 00581 double mu = 0.1056583692; 00582 double en = std::sqrt(gMom.mag2()+mu*mu); 00583 00584 // And now create the Particle 00585 XYZTLorentzVector position(gPos.x(),gPos.y(),gPos.z(),0.); 00586 XYZTLorentzVector momentum(gMom.x(),gMom.y(),gMom.z(),en); 00587 float charge = (float)(tsos.charge()); 00588 ParticlePropagator theMuon(momentum,position,charge,0); 00589 theMuon.setID(-(int)charge*13); 00590 00591 // Recompute the energy loss to get the fluctuations 00592 if ( energyLoss ) { 00593 // Difference between with and without dE/dx (average only) 00594 // (for corrections once fluctuations are applied) 00595 GlobalPoint gPosWithdEdx = tsosWithdEdx.globalPosition(); 00596 GlobalVector gMomWithdEdx = tsosWithdEdx.globalMomentum(); 00597 double enWithdEdx = std::sqrt(gMomWithdEdx.mag2()+mu*mu); 00598 XYZTLorentzVector 00599 deltaPos(gPosWithdEdx.x()-gPos.x(),gPosWithdEdx.y()-gPos.y(), 00600 gPosWithdEdx.z()-gPos.z(),0.); 00601 XYZTLorentzVector 00602 deltaMom(gMomWithdEdx.x()-gMom.x(),gMomWithdEdx.y()-gMom.y(), 00603 gMomWithdEdx.z()-gMom.z(), enWithdEdx -en); 00604 00605 // Energy loss in iron (+ fluctuations) 00606 energyLoss->updateState(theMuon,radPath); 00607 00608 // Correcting factors to account for slight differences in material descriptions 00609 // (Material description is more accurate in the stepping helix propagator) 00610 radPath *= -deltaMom.E()/energyLoss->mostLikelyLoss(); 00611 double fac = energyLoss->deltaMom().E()/energyLoss->mostLikelyLoss(); 00612 00613 // Particle momentum & position after energy loss + fluctuation 00614 XYZTLorentzVector theNewMomentum = theMuon.momentum() + energyLoss->deltaMom() + fac * deltaMom; 00615 XYZTLorentzVector theNewPosition = theMuon.vertex() + fac * deltaPos; 00616 fac = std::sqrt((theNewMomentum.E()*theNewMomentum.E()-mu*mu)/theNewMomentum.Vect().Mag2()); 00617 theMuon.SetXYZT(theNewMomentum.Px()*fac,theNewMomentum.Py()*fac, 00618 theNewMomentum.Pz()*fac,theNewMomentum.E()); 00619 theMuon.setVertex(theNewPosition); 00620 00621 } 00622 00623 // Does the actual mutliple scattering 00624 if ( multipleScattering ) { 00625 // Pass the vector normal to the "next" surface 00626 GlobalVector normal = nextSurface.tangentPlane(tsos.globalPosition())->normalVector(); 00627 multipleScattering->setNormalVector(normal); 00628 // Compute the amount of multiple scattering after a given path length 00629 multipleScattering->updateState(theMuon,radPath); 00630 } 00631 00632 // Fill the propagated state 00633 GlobalPoint propagatedPosition(theMuon.X(),theMuon.Y(),theMuon.Z()); 00634 GlobalVector propagatedMomentum(theMuon.Px(),theMuon.Py(),theMuon.Pz()); 00635 GlobalTrajectoryParameters propagatedGtp(propagatedPosition,propagatedMomentum,(int)charge,magfield); 00636 tsosWithdEdx = TrajectoryStateOnSurface(propagatedGtp,nextSurface); 00637 00638 }
| void MuonSimHitProducer::beginRun | ( | edm::Run & | run, | |
| const edm::EventSetup & | es | |||
| ) | [private, virtual] |
Reimplemented from edm::EDProducer.
Definition at line 117 of file MuonSimHitProducer.cc.
References Propagator::clone(), cscGeom, dtGeom, edm::EventSetup::get(), magfield, propagatorWithMaterial, propagatorWithoutMaterial, rpcGeom, SteppingHelixPropagator::setMaterialMode(), and theService.
00117 { 00118 00119 //services 00120 edm::ESHandle<MagneticField> magField; 00121 edm::ESHandle<DTGeometry> dtGeometry; 00122 edm::ESHandle<CSCGeometry> cscGeometry; 00123 edm::ESHandle<RPCGeometry> rpcGeometry; 00124 edm::ESHandle<Propagator> propagator; 00125 00126 es.get<IdealMagneticFieldRecord>().get(magField); 00127 es.get<MuonGeometryRecord>().get(dtGeometry); 00128 es.get<MuonGeometryRecord>().get(cscGeometry); 00129 es.get<MuonGeometryRecord>().get(rpcGeometry); 00130 00131 magfield = &(*magField); 00132 dtGeom = &(*dtGeometry); 00133 cscGeom = &(*cscGeometry); 00134 rpcGeom = &(*rpcGeometry); 00135 00136 theService->update(es); 00137 00138 // A few propagators 00139 propagatorWithMaterial = &(*(theService->propagator("SteppingHelixPropagatorAny"))); 00140 propagatorWithoutMaterial = propagatorWithMaterial->clone(); 00141 SteppingHelixPropagator* SHpropagator = dynamic_cast<SteppingHelixPropagator*> (propagatorWithoutMaterial); // Beuark! 00142 SHpropagator->setMaterialMode(true); // switches OFF material effects; 00143 00144 }
| void MuonSimHitProducer::produce | ( | edm::Event & | iEvent, | |
| const edm::EventSetup & | iSetup | |||
| ) | [private, virtual] |
Implements edm::EDProducer.
Definition at line 166 of file MuonSimHitProducer.cc.
References funct::abs(), alongMomentum, anyDirection, applyMaterialEffects(), PV3DBase< T, PVType, FrameType >::basicVector(), BoundSurface::bounds(), RPCGeometry::chamber(), DTGeometry::chamber(), CSCGeometry::chamber(), DTTopology::channel(), CoreSimTrack::charge(), GenMuonPlsPt100GeV_cfg::cout, GeomDetEnumerators::CSC, cscGeom, GeomDetEnumerators::DT, dtGeom, MuonPatternRecoDumper::dumpLayer(), CSCDetId::endcap(), lat::endl(), MuonTrajectoryUpdator::estimator(), cmsRelvalreport::exit, DTTopology::firstChannel(), TrajectoryStateOnSurface::freeTrajectoryState(), GeomDet::geographicalId(), edm::Event::getByLabel(), SteppingHelixStateInfo::getStateOnSurface(), TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), i, id, CSCGeometry::idToDetUnit(), RPCGeometry::idToDetUnit(), DTGeometry::idToDetUnit(), Bounds::inside(), TrajectoryStateOnSurface::isValid(), DTTopology::lastChannel(), DTLayerId::layer(), CSCDetId::layer(), RPCDetId::layer(), CSCChamber::layers(), TrajectoryStateOnSurface::localMomentum(), m, PV3DBase< T, PVType, FrameType >::mag(), magfield, CoreSimTrack::momentum(), n, SteppingHelixStateInfo::path(), path(), PV3DBase< T, PVType, FrameType >::phi(), pi, PlaneBuilder::plane(), Propagator::propagateWithPath(), propagatorWithMaterial, propagatorWithoutMaterial, edm::Event::put(), SteppingHelixStateInfo::radPath(), DetId::rawId(), RPCDetId::ring(), CSCDetId::ring(), RPCDetId::roll(), RPCChamber::rolls(), rot, GeomDetEnumerators::RPCBarrel, GeomDetEnumerators::RPCEndcap, rpcGeom, DTChamberId::sector(), RPCDetId::sector(), funct::sqrt(), RPCDetId::station(), DTChamberId::station(), CSCDetId::station(), GeomDet::subDetector(), RPCDetId::subsector(), DTSuperLayerId::superlayer(), DTChamber::superLayers(), GeomDet::surface(), theMaterialEffects, theService, theSimModuleLabel_, theSimModuleProcess_, PV3DBase< T, PVType, FrameType >::theta(), theUpdator, Bounds::thickness(), GeomDet::toLocal(), SimTrack::trackerSurfaceMomentum(), SimTrack::trackerSurfacePosition(), CoreSimTrack::trackId(), TrajectoryStateOnSurface::transverseCurvature(), CoreSimTrack::type(), Vector3DBase< T, FrameTag >::unit(), SimTrack::vertIndex(), DTChamberId::wheel(), PV3DBase< T, PVType, FrameType >::x(), x, PV3DBase< T, PVType, FrameType >::y(), y, PV3DBase< T, PVType, FrameType >::z(), and z.
00166 { 00167 // using namespace edm; 00168 // using namespace std; 00169 00170 MuonPatternRecoDumper dumper; 00171 00172 edm::Handle<std::vector<SimTrack> > simMuons; 00173 edm::Handle<std::vector<SimVertex> > simVertices; 00174 std::vector<PSimHit> theCSCHits; 00175 std::vector<PSimHit> theDTHits; 00176 std::vector<PSimHit> theRPCHits; 00177 00178 DirectMuonNavigation navigation(theService->detLayerGeometry()); 00179 iEvent.getByLabel(theSimModuleLabel_,theSimModuleProcess_,simMuons); 00180 iEvent.getByLabel(theSimModuleLabel_,simVertices); 00181 00182 for ( unsigned int itrk=0; itrk<simMuons->size(); itrk++ ) { 00183 const SimTrack &mySimTrack = (*simMuons)[itrk]; 00184 math::XYZTLorentzVector mySimP4(mySimTrack.momentum().x(), 00185 mySimTrack.momentum().y(), 00186 mySimTrack.momentum().z(), 00187 mySimTrack.momentum().t()); 00188 00189 // Decaying hadrons are now in the list, and so are their muon daughter 00190 // Ignore the hadrons here. 00191 int pid = mySimTrack.type(); 00192 if ( abs(pid) != 13 ) continue; 00193 00194 double t0 = 0; 00195 GlobalPoint initialPosition; 00196 int ivert = mySimTrack.vertIndex(); 00197 if ( ivert >= 0 ) { 00198 t0 = (*simVertices)[ivert].position().t(); 00199 GlobalPoint xyzzy((*simVertices)[ivert].position().x(), 00200 (*simVertices)[ivert].position().y(), 00201 (*simVertices)[ivert].position().z()); 00202 initialPosition = xyzzy; 00203 } 00204 // 00205 // Presumably t0 has dimensions of cm if not mm? 00206 // Convert to ns for internal calculations. 00207 // I wonder where we should get c from? 00208 // 00209 double tof = t0/29.98; 00210 00211 #ifdef FAMOS_DEBUG 00212 std::cout << " ===> MuonSimHitProducer::reconstruct() found SIMTRACK - pid = " 00213 << pid ; 00214 std::cout << " : pT = " << mySimP4.Pt() 00215 << ", eta = " << mySimP4.Eta() 00216 << ", phi = " << mySimP4.Phi() << std::endl; 00217 #endif 00218 00219 // 00220 // Produce muons sim hits starting from undecayed simulated muons 00221 // 00222 00223 GlobalPoint startingPosition(mySimTrack.trackerSurfacePosition().x(), 00224 mySimTrack.trackerSurfacePosition().y(), 00225 mySimTrack.trackerSurfacePosition().z()); 00226 GlobalVector startingMomentum(mySimTrack.trackerSurfaceMomentum().x(), 00227 mySimTrack.trackerSurfaceMomentum().y(), 00228 mySimTrack.trackerSurfaceMomentum().z()); 00229 // 00230 // Crap... there's no time-of-flight to the trackerSurfacePosition()... 00231 // So, this will be wrong when the curvature can't be neglected, but that 00232 // will be rather seldom... May as well ignore the mass too. 00233 // 00234 GlobalVector dtracker = startingPosition-initialPosition; 00235 tof += dtracker.mag()/29.98; 00236 00237 #ifdef FAMOS_DEBUG 00238 std::cout << " the Muon START position " << startingPosition << std::endl; 00239 std::cout << " the Muon START momentum " << startingMomentum << std::endl; 00240 #endif 00241 00242 // 00243 // Some magic to define a TrajectoryStateOnSurface 00244 // 00245 PlaneBuilder pb; 00246 GlobalVector zAxis = startingMomentum.unit(); 00247 GlobalVector yAxis(zAxis.y(),-zAxis.x(),0); 00248 GlobalVector xAxis = yAxis.cross(zAxis); 00249 Surface::RotationType rot = Surface::RotationType(xAxis,yAxis,zAxis); 00250 PlaneBuilder::ReturnType startingPlane = pb.plane(startingPosition,rot); 00251 GlobalTrajectoryParameters gtp(startingPosition, 00252 startingMomentum, 00253 (int)mySimTrack.charge(), 00254 magfield); 00255 TrajectoryStateOnSurface startingState(gtp,*startingPlane); 00256 00257 std::vector<const DetLayer *> navLayers; 00258 if ( fabs(startingState.globalMomentum().eta()) > 4.5 ) { 00259 navLayers = navigation.compatibleEndcapLayers(*(startingState.freeState()), 00260 alongMomentum); 00261 } 00262 else { 00263 navLayers = navigation.compatibleLayers(*(startingState.freeState()), 00264 alongMomentum); 00265 } 00266 /* 00267 edm::ESHandle<Propagator> propagator = 00268 theService->propagator("SteppingHelixPropagatorAny"); 00269 */ 00270 00271 if ( navLayers.empty() ) continue; 00272 00273 #ifdef FAMOS_DEBUG 00274 std::cout << "Found " << navLayers.size() 00275 << " compatible DetLayers..." << std::endl; 00276 #endif 00277 00278 TrajectoryStateOnSurface propagatedState = startingState; 00279 for ( unsigned int ilayer=0; ilayer<navLayers.size(); ilayer++ ) { 00280 00281 #ifdef FAMOS_DEBUG 00282 std::cout << "Propagating to layer " << ilayer << " " 00283 << dumper.dumpLayer(navLayers[ilayer]) 00284 << std::endl; 00285 #endif 00286 00287 std::vector<DetWithState> comps = 00288 navLayers[ilayer]->compatibleDets(propagatedState,*propagatorWithMaterial,*(theUpdator->estimator())); 00289 if ( comps.empty() ) continue; 00290 00291 #ifdef FAMOS_DEBUG 00292 std::cout << "Propagating " << propagatedState << std::endl; 00293 #endif 00294 00295 // Strating momentum 00296 double pi = propagatedState.globalMomentum().mag(); 00297 00298 // Propgate with material effects (dE/dx average only) 00299 SteppingHelixStateInfo shsStart(*(propagatedState.freeTrajectoryState())); 00300 const SteppingHelixStateInfo& shsDest = 00301 ((const SteppingHelixPropagator*)propagatorWithMaterial)->propagate(shsStart,navLayers[ilayer]->surface()); 00302 std::pair<TrajectoryStateOnSurface,double> next(shsDest.getStateOnSurface(navLayers[ilayer]->surface()), 00303 shsDest.path()); 00304 // No need to continue if there is no valid propagation available. 00305 // This happens rarely (~0.1% of ttbar events) 00306 if ( !next.first.isValid() ) continue; 00307 // This is the estimate of the number of radiation lengths traversed, 00308 // together with the total path length 00309 double radPath = shsDest.radPath(); 00310 double pathLength = next.second; 00311 00312 // Now propagate without dE/dx (average) 00313 // [To add the dE/dx fluctuations to the actual dE/dx] 00314 std::pair<TrajectoryStateOnSurface,double> nextNoMaterial = 00315 propagatorWithoutMaterial->propagateWithPath(propagatedState,navLayers[ilayer]->surface()); 00316 00317 // Update the propagated state 00318 propagatedState = next.first; 00319 double pf = propagatedState.globalMomentum().mag(); 00320 00321 // Insert dE/dx fluctuations and multiple scattering 00322 // Skip this step if nextNoMaterial.first is not valid 00323 // This happens rarely (~0.02% of ttbar events) 00324 if ( theMaterialEffects && nextNoMaterial.first.isValid() ) applyMaterialEffects(propagatedState,nextNoMaterial.first,radPath); 00325 // Check that the 'shaken' propagatedState is still valid, otherwise continue 00326 if ( !propagatedState.isValid() ) continue; 00327 // (No evidence that this ever happens) 00328 // 00329 // Consider this... 1 GeV muon has a velocity that is only 0.5% slower than c... 00330 // We probably can safely ignore the mass for anything that makes it out to the 00331 // muon chambers. 00332 // 00333 double pavg = 0.5*(pi+pf); 00334 double m = mySimP4.M(); 00335 double rbeta = sqrt(1+m*m/(pavg*pavg))/29.98; 00336 double dtof = pathLength*rbeta; 00337 00338 #ifdef FAMOS_DEBUG 00339 std::cout << "Propagated to next surface... path length = " << pathLength 00340 << " cm, dTOF = " << dtof << " ns" << std::endl; 00341 #endif 00342 00343 tof += dtof; 00344 00345 for ( unsigned int icomp=0; icomp<comps.size(); icomp++ ) 00346 { 00347 const GeomDet *gd = comps[icomp].first; 00348 if ( gd->subDetector() == GeomDetEnumerators::DT ) { 00349 DTChamberId id(gd->geographicalId()); 00350 const DTChamber *chamber = dtGeom->chamber(id); 00351 std::vector<const DTSuperLayer *> superlayer = chamber->superLayers(); 00352 for ( unsigned int isl=0; isl<superlayer.size(); isl++ ) { 00353 std::vector<const DTLayer *> layer = superlayer[isl]->layers(); 00354 for ( unsigned int ilayer=0; ilayer<layer.size(); ilayer++ ) { 00355 DTLayerId lid = layer[ilayer]->id(); 00356 #ifdef FAMOS_DEBUG 00357 std::cout << " Extrapolated to DT (" 00358 << lid.wheel() << "," 00359 << lid.station() << "," 00360 << lid.sector() << "," 00361 << lid.superlayer() << "," 00362 << lid.layer() << ")" << std::endl; 00363 #endif 00364 00365 const GeomDetUnit *det = dtGeom->idToDetUnit(lid); 00366 00367 HelixArbitraryPlaneCrossing crossing(propagatedState.globalPosition().basicVector(), 00368 propagatedState.globalMomentum().basicVector(), 00369 propagatedState.transverseCurvature(), 00370 anyDirection); 00371 std::pair<bool,double> path = crossing.pathLength(det->surface()); 00372 if ( ! path.first ) continue; 00373 LocalPoint lpos = det->toLocal(GlobalPoint(crossing.position(path.second))); 00374 if ( ! det->surface().bounds().inside(lpos) ) continue; 00375 const DTTopology& dtTopo = layer[ilayer]->specificTopology(); 00376 int wire = dtTopo.channel(lpos); 00377 if (wire < dtTopo.firstChannel() || wire > dtTopo.lastChannel()) continue; 00378 // no drift cell here (on the chamber edge or just outside) 00379 // this hit would otherwise be discarded downstream in the digitizer 00380 00381 DTWireId wid(lid,wire); 00382 double thickness = det->surface().bounds().thickness(); 00383 LocalVector lmom = det->toLocal(GlobalVector(crossing.direction(path.second))); 00384 lmom = lmom.unit()*propagatedState.localMomentum().mag(); 00385 double eloss = 0; 00386 double pz = fabs(lmom.z()); 00387 LocalPoint entry = lpos - 0.5*thickness*lmom/pz; 00388 LocalPoint exit = lpos + 0.5*thickness*lmom/pz; 00389 double dtof = path.second*rbeta; 00390 int trkid = mySimTrack.trackId(); 00391 unsigned int id = wid.rawId(); 00392 short unsigned int processType = 2; 00393 PSimHit hit(entry,exit,lmom.mag(), 00394 tof+dtof,eloss,pid,id,trkid,lmom.theta(),lmom.phi(),processType); 00395 theDTHits.push_back(hit); 00396 00397 } 00398 } 00399 } 00400 else if ( gd->subDetector() == GeomDetEnumerators::CSC ) { 00401 CSCDetId id(gd->geographicalId()); 00402 const CSCChamber *chamber = cscGeom->chamber(id); 00403 std::vector<const CSCLayer *> layer = chamber->layers(); 00404 for ( unsigned int ilayer=0; ilayer<layer.size(); ilayer++ ) { 00405 CSCDetId lid = layer[ilayer]->id(); 00406 00407 #ifdef FAMOS_DEBUG 00408 std::cout << " Extrapolated to CSC (" 00409 << lid.endcap() << "," 00410 << lid.ring() << "," 00411 << lid.station() << "," 00412 << lid.layer() << ")" << std::endl; 00413 #endif 00414 00415 const GeomDetUnit *det = cscGeom->idToDetUnit(lid); 00416 HelixArbitraryPlaneCrossing crossing(propagatedState.globalPosition().basicVector(), 00417 propagatedState.globalMomentum().basicVector(), 00418 propagatedState.transverseCurvature(), 00419 anyDirection); 00420 std::pair<bool,double> path = crossing.pathLength(det->surface()); 00421 if ( ! path.first ) continue; 00422 LocalPoint lpos = det->toLocal(GlobalPoint(crossing.position(path.second))); 00423 if ( ! det->surface().bounds().inside(lpos) ) continue; 00424 double thickness = det->surface().bounds().thickness(); 00425 LocalVector lmom = det->toLocal(GlobalVector(crossing.direction(path.second))); 00426 lmom = lmom.unit()*propagatedState.localMomentum().mag(); 00427 double eloss = 0; 00428 double pz = fabs(lmom.z()); 00429 LocalPoint entry = lpos - 0.5*thickness*lmom/pz; 00430 LocalPoint exit = lpos + 0.5*thickness*lmom/pz; 00431 double dtof = path.second*rbeta; 00432 int trkid = mySimTrack.trackId(); 00433 unsigned int id = lid.rawId(); 00434 short unsigned int processType = 2; 00435 PSimHit hit(entry,exit,lmom.mag(), 00436 tof+dtof,eloss,pid,id,trkid,lmom.theta(),lmom.phi(),processType); 00437 theCSCHits.push_back(hit); 00438 } 00439 } 00440 else if ( gd->subDetector() == GeomDetEnumerators::RPCBarrel || 00441 gd->subDetector() == GeomDetEnumerators::RPCEndcap ) { 00442 RPCDetId id(gd->geographicalId()); 00443 const RPCChamber *chamber = rpcGeom->chamber(id); 00444 std::vector<const RPCRoll *> roll = chamber->rolls(); 00445 for ( unsigned int iroll=0; iroll<roll.size(); iroll++ ) { 00446 RPCDetId rid = roll[iroll]->id(); 00447 00448 #ifdef FAMOS_DEBUG 00449 std::cout << " Extrapolated to RPC (" 00450 << rid.ring() << "," 00451 << rid.station() << "," 00452 << rid.sector() << "," 00453 << rid.subsector() << "," 00454 << rid.layer() << "," 00455 << rid.roll() << ")" << std::endl; 00456 #endif 00457 00458 const GeomDetUnit *det = rpcGeom->idToDetUnit(rid); 00459 HelixArbitraryPlaneCrossing crossing(propagatedState.globalPosition().basicVector(), 00460 propagatedState.globalMomentum().basicVector(), 00461 propagatedState.transverseCurvature(), 00462 anyDirection); 00463 std::pair<bool,double> path = crossing.pathLength(det->surface()); 00464 if ( ! path.first ) continue; 00465 LocalPoint lpos = det->toLocal(GlobalPoint(crossing.position(path.second))); 00466 if ( ! det->surface().bounds().inside(lpos) ) continue; 00467 double thickness = det->surface().bounds().thickness(); 00468 LocalVector lmom = det->toLocal(GlobalVector(crossing.direction(path.second))); 00469 lmom = lmom.unit()*propagatedState.localMomentum().mag(); 00470 double eloss = 0; 00471 double pz = fabs(lmom.z()); 00472 LocalPoint entry = lpos - 0.5*thickness*lmom/pz; 00473 LocalPoint exit = lpos + 0.5*thickness*lmom/pz; 00474 double dtof = path.second*rbeta; 00475 int trkid = mySimTrack.trackId(); 00476 unsigned int id = rid.rawId(); 00477 short unsigned int processType = 2; 00478 PSimHit hit(entry,exit,lmom.mag(), 00479 tof+dtof,eloss,pid,id,trkid,lmom.theta(),lmom.phi(),processType); 00480 theRPCHits.push_back(hit); 00481 } 00482 } 00483 else { 00484 std::cout << "Extrapolated to unknown subdetector '" << gd->subDetector() << "'..." << std::endl; 00485 } 00486 } 00487 } 00488 } 00489 00490 std::auto_ptr<edm::PSimHitContainer> pcsc(new edm::PSimHitContainer); 00491 int n = 0; 00492 for ( std::vector<PSimHit>::const_iterator i = theCSCHits.begin(); 00493 i != theCSCHits.end(); i++ ) { 00494 pcsc->push_back(*i); 00495 n += 1; 00496 } 00497 iEvent.put(pcsc,"MuonCSCHits"); 00498 00499 std::auto_ptr<edm::PSimHitContainer> pdt(new edm::PSimHitContainer); 00500 n = 0; 00501 for ( std::vector<PSimHit>::const_iterator i = theDTHits.begin(); 00502 i != theDTHits.end(); i++ ) { 00503 pdt->push_back(*i); 00504 n += 1; 00505 } 00506 iEvent.put(pdt,"MuonDTHits"); 00507 00508 std::auto_ptr<edm::PSimHitContainer> prpc(new edm::PSimHitContainer); 00509 n = 0; 00510 for ( std::vector<PSimHit>::const_iterator i = theRPCHits.begin(); 00511 i != theRPCHits.end(); i++ ) { 00512 prpc->push_back(*i); 00513 n += 1; 00514 } 00515 iEvent.put(prpc,"MuonRPCHits"); 00516 00517 }
| void MuonSimHitProducer::readParameters | ( | const edm::ParameterSet & | fastMuons, | |
| const edm::ParameterSet & | fastTracks, | |||
| const edm::ParameterSet & | matEff | |||
| ) | [private] |
Definition at line 528 of file MuonSimHitProducer.cc.
References fullPattern_, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), maxEta_, minEta_, random, theMaterialEffects, theSimModuleLabel_, theSimModuleProcess_, and theTrkModuleLabel_.
Referenced by MuonSimHitProducer().
00530 { 00531 // Muons 00532 theSimModuleLabel_ = fastMuons.getParameter<std::string>("simModuleLabel"); 00533 theSimModuleProcess_ = fastMuons.getParameter<std::string>("simModuleProcess"); 00534 theTrkModuleLabel_ = fastMuons.getParameter<std::string>("trackModuleLabel"); 00535 minEta_ = fastMuons.getParameter<double>("MinEta"); 00536 maxEta_ = fastMuons.getParameter<double>("MaxEta"); 00537 if (minEta_ > maxEta_) { 00538 double tempEta_ = maxEta_ ; 00539 maxEta_ = minEta_ ; 00540 minEta_ = tempEta_ ; 00541 } 00542 00543 // Tracks 00544 fullPattern_ = fastTracks.getUntrackedParameter<bool>("FullPatternRecognition"); 00545 00546 // The following should be on LogInfo 00547 // std::cout << " MUON SIM HITS: FastSimulation parameters " << std::endl; 00548 // std::cout << " ============================================== " << std::endl; 00549 // std::cout << " Sim Hits produced for muons in the pseudorapidity range : " 00550 // << minEta_ << " -> " << maxEta_ << std::endl; 00551 //if ( fullPattern_ ) 00552 //std::cout << " The FULL pattern recognition option is turned ON" << std::endl; 00553 //else 00554 //std::cout << " The FAST tracking option is turned ON" << std::endl; 00555 00556 // Material Effects 00557 theMaterialEffects = 0; 00558 if ( matEff.getParameter<bool>("PairProduction") || 00559 matEff.getParameter<bool>("Bremsstrahlung") || 00560 matEff.getParameter<bool>("EnergyLoss") || 00561 matEff.getParameter<bool>("MultipleScattering") ) 00562 theMaterialEffects = new MaterialEffects(matEff,random); 00563 00564 }
const CSCGeometry* MuonSimHitProducer::cscGeom [private] |
bool MuonSimHitProducer::doGL_ [private] |
Definition at line 93 of file MuonSimHitProducer.h.
bool MuonSimHitProducer::doL1_ [private] |
Definition at line 93 of file MuonSimHitProducer.h.
bool MuonSimHitProducer::doL3_ [private] |
Definition at line 93 of file MuonSimHitProducer.h.
const DTGeometry* MuonSimHitProducer::dtGeom [private] |
bool MuonSimHitProducer::fullPattern_ [private] |
const MagneticField* MuonSimHitProducer::magfield [private] |
Definition at line 69 of file MuonSimHitProducer.h.
Referenced by applyMaterialEffects(), beginRun(), and produce().
double MuonSimHitProducer::maxEta_ [private] |
double MuonSimHitProducer::minEta_ [private] |
const Propagator* MuonSimHitProducer::propagatorWithMaterial [private] |
const RandomEngine* MuonSimHitProducer::random [private] |
Definition at line 65 of file MuonSimHitProducer.h.
Referenced by MuonSimHitProducer(), readParameters(), and ~MuonSimHitProducer().
const RPCGeometry* MuonSimHitProducer::rpcGeom [private] |
Definition at line 76 of file MuonSimHitProducer.h.
Referenced by applyMaterialEffects(), produce(), readParameters(), and ~MuonSimHitProducer().
MuonServiceProxy* MuonSimHitProducer::theService [private] |
Definition at line 66 of file MuonSimHitProducer.h.
Referenced by beginRun(), MuonSimHitProducer(), and produce().
std::string MuonSimHitProducer::theSimModuleLabel_ [private] |
std::string MuonSimHitProducer::theSimModuleProcess_ [private] |
std::string MuonSimHitProducer::theTrkModuleLabel_ [private] |
Definition at line 67 of file MuonSimHitProducer.h.
Referenced by MuonSimHitProducer(), and produce().
1.5.4