79 theEstimator(iConfig.getParameter<double>(
"Chi2EstimatorCut")),
80 propagatorWithoutMaterial(0)
89 "MuonSimHitProducer requires the RandomGeneratorService \n"
90 "which is not present in the configuration file. \n"
91 "You must add the service in the configuration file\n"
92 "or remove the module that requires it.";
105 produces<edm::PSimHitContainer>(
"MuonCSCHits");
106 produces<edm::PSimHitContainer>(
"MuonDTHits");
107 produces<edm::PSimHitContainer>(
"MuonRPCHits");
174 std::vector<PSimHit> theCSCHits;
175 std::vector<PSimHit> theDTHits;
176 std::vector<PSimHit> theRPCHits;
182 for (
unsigned int itrk=0; itrk<simMuons->size(); itrk++ ) {
183 const SimTrack &mySimTrack = (*simMuons)[itrk];
192 if (
abs(pid) != 13 )
continue;
198 t0 = (*simVertices)[ivert].position().t();
202 initialPosition = xyzzy;
209 double tof = t0/29.98;
212 std::cout <<
" ===> MuonSimHitProducer::reconstruct() found SIMTRACK - pid = "
215 <<
", eta = " << mySimP4.Eta()
216 <<
", phi = " << mySimP4.Phi() << std::endl;
234 GlobalVector dtracker = startingPosition-initialPosition;
235 tof += dtracker.
mag()/29.98;
238 std::cout <<
" the Muon START position " << startingPosition << std::endl;
239 std::cout <<
" the Muon START momentum " << startingMomentum << std::endl;
257 std::vector<const DetLayer *> navLayers;
258 if ( fabs(startingState.globalMomentum().eta()) > 4.5 ) {
259 navLayers = navigation.compatibleEndcapLayers(*(startingState.freeState()),
263 navLayers = navigation.compatibleLayers(*(startingState.freeState()),
271 if ( navLayers.empty() )
continue;
274 std::cout <<
"Found " << navLayers.size()
275 <<
" compatible DetLayers..." << std::endl;
279 for (
unsigned int ilayer=0; ilayer<navLayers.size(); ilayer++ ) {
282 std::cout <<
"Propagating to layer " << ilayer <<
" "
287 std::vector<DetWithState> comps =
289 if ( comps.empty() )
continue;
292 std::cout <<
"Propagating " << propagatedState << std::endl;
302 std::pair<TrajectoryStateOnSurface,double> next(shsDest.
getStateOnSurface(navLayers[ilayer]->surface()),
306 if ( !next.first.isValid() )
continue;
309 double radPath = shsDest.
radPath();
310 double pathLength = next.second;
314 std::pair<TrajectoryStateOnSurface,double> nextNoMaterial =
318 propagatedState = next.first;
326 if ( !propagatedState.
isValid() )
continue;
333 double pavg = 0.5*(pi+pf);
334 double m = mySimP4.M();
335 double rbeta =
sqrt(1+m*m/(pavg*pavg))/29.98;
336 double dtof = pathLength*rbeta;
339 std::cout <<
"Propagated to next surface... path length = " << pathLength
340 <<
" cm, dTOF = " << dtof <<
" ns" << std::endl;
345 for (
unsigned int icomp=0; icomp<comps.size(); icomp++ )
347 const GeomDet *gd = comps[icomp].first;
351 std::vector<const DTSuperLayer *> superlayer = chamber->
superLayers();
352 for (
unsigned int isl=0; isl<superlayer.size(); isl++ ) {
353 std::vector<const DTLayer *> layer = superlayer[isl]->layers();
354 for (
unsigned int ilayer=0; ilayer<layer.size(); ilayer++ ) {
358 << lid.
wheel() <<
","
362 << lid.
layer() <<
")" << std::endl;
371 std::pair<bool,double>
path = crossing.pathLength(det->
surface());
372 if ( ! path.first )
continue;
375 const DTTopology& dtTopo = layer[ilayer]->specificTopology();
376 int wire = dtTopo.
channel(lpos);
388 double pmu = lmom.
mag();
393 double pz = fabs(lmom.
z());
396 double dtof = path.second*rbeta;
397 int trkid = mySimTrack.
trackId();
398 unsigned int id = wid.
rawId();
399 short unsigned int processType = 2;
401 tof+dtof,eloss,
pid,
id,trkid,lmom.
theta(),lmom.
phi(),processType);
402 theDTHits.push_back(
hit);
410 std::vector<const CSCLayer *> layers = chamber->
layers();
411 for (
unsigned int ilayer=0; ilayer<layers.size(); ilayer++ ) {
412 CSCDetId lid = layers[ilayer]->id();
419 << lid.
layer() <<
")" << std::endl;
427 std::pair<bool,double>
path = crossing.pathLength(det->
surface());
428 if ( ! path.first )
continue;
434 if ( ! laygeom->
inside( lpos ) )
continue;
442 double pmu = lmom.
mag();
445 if (
id.
station()==1 &&
id.
ring()==1) theCSCHitIneff = theCSCHitIneff*0.442;
449 double pz = fabs(lmom.
z());
452 double dtof = path.second*rbeta;
453 int trkid = mySimTrack.
trackId();
454 unsigned int id = lid.
rawId();
455 short unsigned int processType = 2;
457 tof+dtof,eloss,
pid,
id,trkid,lmom.
theta(),lmom.
phi(),processType);
458 theCSCHits.push_back(
hit);
465 std::vector<const RPCRoll *> roll = chamber->
rolls();
466 for (
unsigned int iroll=0; iroll<roll.size(); iroll++ ) {
475 << rid.
layer() <<
","
476 << rid.
roll() <<
")" << std::endl;
484 std::pair<bool,double>
path = crossing.pathLength(det->
surface());
485 if ( ! path.first )
continue;
492 double pz = fabs(lmom.
z());
495 double dtof = path.second*rbeta;
496 int trkid = mySimTrack.
trackId();
497 unsigned int id = rid.
rawId();
498 short unsigned int processType = 2;
500 tof+dtof,eloss,
pid,
id,trkid,lmom.
theta(),lmom.
phi(),processType);
501 theRPCHits.push_back(
hit);
505 std::cout <<
"Extrapolated to unknown subdetector '" << gd->
subDetector() <<
"'..." << std::endl;
513 for ( std::vector<PSimHit>::const_iterator
i = theCSCHits.begin();
514 i != theCSCHits.end();
i++ ) {
518 iEvent.
put(pcsc,
"MuonCSCHits");
522 for ( std::vector<PSimHit>::const_iterator
i = theDTHits.begin();
523 i != theDTHits.end();
i++ ) {
527 iEvent.
put(pdt,
"MuonDTHits");
531 for ( std::vector<PSimHit>::const_iterator
i = theRPCHits.begin();
532 i != theRPCHits.end();
i++ ) {
536 iEvent.
put(prpc,
"MuonRPCHits");
556 std::vector<double> simHitIneffDT = fastMuons.
getParameter<std::vector<double> >(
"simHitDTIneffParameters");
557 std::vector<double> simHitIneffCSC = fastMuons.
getParameter<std::vector<double> >(
"simHitCSCIneffParameters");
558 kDT = simHitIneffDT[0];
559 fDT = simHitIneffDT[1];
560 kCSC = simHitIneffCSC[0];
561 fCSC = simHitIneffCSC[1];
599 double mu = 0.1056583692;
607 theMuon.
setID(-(
int)charge*13);
617 deltaPos(gPosWithdEdx.
x()-gPos.
x(),gPosWithdEdx.
y()-gPos.
y(),
618 gPosWithdEdx.
z()-gPos.
z(),0.);
620 deltaMom(gMomWithdEdx.
x()-gMom.
x(),gMomWithdEdx.
y()-gMom.
y(),
621 gMomWithdEdx.
z()-gMom.
z(), enWithdEdx -en);
634 fac = (theNewMomentum.E()*theNewMomentum.E()-mu*mu)/theNewMomentum.Vect().Mag2();
636 theMuon.SetXYZT(theNewMomentum.Px()*fac,theNewMomentum.Py()*fac,
637 theNewMomentum.Pz()*fac,theNewMomentum.E());
643 if ( multipleScattering ) {
652 GlobalPoint propagatedPosition(theMuon.
X(),theMuon.
Y(),theMuon.
Z());
653 GlobalVector propagatedMomentum(theMuon.Px(),theMuon.Py(),theMuon.Pz());
void update(const edm::EventSetup &setup)
update the services each event
MuonSimHitProducer(const edm::ParameterSet &)
double mostLikelyLoss() const
Return most probable energy loss.
edm::ESHandle< MuonDetLayerGeometry > detLayerGeometry() const
get the detLayer geometry
const int lastChannel() const
Returns the wire number of the last wire.
const math::XYZVectorD & trackerSurfacePosition() const
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
const MagneticField * magfield
const std::vector< const CSCLayer * > & layers() const
Return all layers.
virtual void beginRun(edm::Run &run, const edm::EventSetup &es)
const DTGeometry * dtGeom
TrackCharge charge() const
std::string dumpLayer(const DetLayer *layer) const
virtual bool inside(const Local3DPoint &) const =0
Determine if the point is inside the bounds.
virtual const GeomDetUnit * idToDetUnit(DetId) const
Return the pointer to the GeomDetUnit corresponding to a given DetId.
void setNormalVector(const GlobalVector &normal)
Sets the vector normal to the surface traversed.
#define DEFINE_FWK_MODULE(type)
FreeTrajectoryState * freeTrajectoryState(bool withErrors=true) const
int channel(const LocalPoint &p) const
MaterialEffects * theMaterialEffects
Geom::Phi< T > phi() const
MultipleScatteringSimulator * multipleScatteringSimulator() const
Return the Multiple Scattering engine.
ReturnType plane(const PositionType &pos, const RotationType &rot) const
Global3DPoint GlobalPoint
virtual SubDetector subDetector() const =0
Which subdetector.
virtual Propagator * clone() const =0
GlobalPoint globalPosition() const
LocalPoint toLocal(const GlobalPoint &gp) const
Conversion to the R.F. of the GeomDet.
Exp< T >::type exp(const T &t)
int layer() const
Return the layer number.
const int firstChannel() const
Returns the wire number of the first wire.
float charge() const
charge
static int position[TOTALCHAMBERS][3]
Geom::Theta< T > theta() const
const RPCGeometry * rpcGeom
uint32_t rawId() const
get the raw id
LocalVector localMomentum() const
virtual float thickness() const =0
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
const RPCChamber * chamber(RPCDetId id) const
virtual std::pair< TrajectoryStateOnSurface, double > propagateWithPath(const FreeTrajectoryState &, const Surface &) const
const XYZTLorentzVector & momentum() const
the momentum fourvector
std::string theSimModuleProcess_
virtual const GeomDetUnit * idToDetUnit(DetId) const
Return the pointer to the GeomDetUnit corresponding to a given DetId.
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
const std::vector< const DTSuperLayer * > & superLayers() const
Return the superlayers in the chamber.
const Propagator * propagatorWithMaterial
TrajectoryStateOnSurface getStateOnSurface(const Surface &surf, bool returnTangentPlane=false) const
Propagator * propagatorWithoutMaterial
std::pair< std::string, MonitorElement * > entry
Vector3DBase< typename PreciseFloatType< T, U >::Type, FrameTag > cross(const Vector3DBase< U, FrameTag > &v) const
double Y() const
y of vertex
const XYZTLorentzVector & deltaMom() const
Returns the actual energy lost.
double Z() const
z of vertex
const std::vector< const RPCRoll * > & rolls() const
Return the Rolls.
const CSCGeometry * cscGeom
DetId geographicalId() const
The label of this GeomDet.
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
int vertIndex() const
index of the vertex in the Event container (-1 if no vertex)
const XYZTLorentzVector & vertex() const
the vertex fourvector
int superlayer() const
Return the superlayer number (deprecated method name)
Vector3DBase unit() const
const Bounds & bounds() const
virtual ReferenceCountingPointer< TangentPlane > tangentPlane(const GlobalPoint &) const =0
unsigned int trackId() const
void updateState(ParticlePropagator &myTrack, double radlen)
Compute the material effect (calls the sub class)
Log< T >::type log(const T &t)
const math::XYZTLorentzVectorD & trackerSurfaceMomentum() const
const CSCChamber * chamber(CSCDetId id) const
Return the chamber corresponding to given DetId.
const RandomEngine * random
const DTChamber * chamber(DTChamberId id) const
Return a DTChamber given its id.
double flatShoot(double xmin=0.0, double xmax=1.0) const
int sector() const
Sector id: the group of chambers at same phi (and increasing r)
double X() const
x of vertex
MuonServiceProxy * theService
void applyMaterialEffects(TrajectoryStateOnSurface &tsosWithdEdx, TrajectoryStateOnSurface &tsos, double radPath)
Simulate material effects in iron (dE/dx, multiple scattering)
int subsector() const
SubSector id : some sectors are divided along the phi direction in subsectors (from 1 to 4 in Barrel...
int type() const
particle type (HEP PDT convension)
std::string theSimModuleLabel_
void readParameters(const edm::ParameterSet &, const edm::ParameterSet &, const edm::ParameterSet &)
virtual void produce(edm::Event &, const edm::EventSetup &)
GlobalVector globalMomentum() const
std::string theTrkModuleLabel_
const math::XYZTLorentzVectorD & momentum() const
particle info...
TkRotation< float > RotationType
const Surface & surface() const
virtual const GeomDetUnit * idToDetUnit(DetId) const
Return the pointer to the GeomDetUnit corresponding to a given DetId.
bool inside(const Local3DPoint &, const LocalError &, float scale=1.f) const
std::vector< PSimHit > PSimHitContainer
int station() const
Return the station number.
EnergyLossSimulator * energyLossSimulator() const
Return the Energy Loss engine.
int wheel() const
Return the wheel number.
const BasicVectorType & basicVector() const
virtual const BoundPlane & surface() const
The nominal surface of the GeomDet.
void setMaterialMode(bool noMaterial)
Switch for material effects mode: no material effects if true.
void setVertex(const XYZTLorentzVector &vtx)
set the vertex
edm::ESHandle< Propagator > propagator(std::string propagatorName) const
get the propagator
Global3DVector GlobalVector
math::XYZTLorentzVector XYZTLorentzVector
Chi2MeasurementEstimator theEstimator
double transverseCurvature() const