CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
MuonSimHitProducer.cc
Go to the documentation of this file.
1 //
2 // Package: MuonSimHitProducer
3 // Class: MuonSimHitProducer
4 //
13 //
14 // Original Author: Martijn Mulders/Matthew Jones
15 // Created: Wed Jul 30 11:37:24 CET 2007
16 // Working: Fri Nov 9 09:39:33 CST 2007
17 //
18 // $Id: MuonSimHitProducer.cc,v 1.33 2010/05/20 15:03:51 aperrott Exp $
19 //
20 //
21 
22 // CMSSW headers
26 
27 // Fast Simulation headers
34 
35 // SimTrack
39 
40 // STL headers
41 #include <vector>
42 #include <iostream>
43 
44 // RecoMuon headers
49 
50 // Tracking Tools
53 
54 // Data Formats
58 
59 
61 // Geometry, Magnetic Field
67 
69 
70 // #include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
71 
72 //for debug only
73 //#define FAMOS_DEBUG
74 
75 //
76 // constructors and destructor
77 //
79  theEstimator(iConfig.getParameter<double>("Chi2EstimatorCut")),
80  propagatorWithoutMaterial(0)
81  {
82 
83  //
84  // Initialize the random number generator service
85  //
87  if ( ! rng.isAvailable() ) {
88  throw cms::Exception("Configuration") <<
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.";
93  }
94 
95  random = new RandomEngine(&(*rng));
96 
97  // Read relevant parameters
99  iConfig.getParameter<edm::ParameterSet>("TRACKS"),
100  iConfig.getParameter<edm::ParameterSet>("MaterialEffectsForMuons"));
101 
102  //
103  // register your products ... need to declare at least one possible product...
104  //
105  produces<edm::PSimHitContainer>("MuonCSCHits");
106  produces<edm::PSimHitContainer>("MuonDTHits");
107  produces<edm::PSimHitContainer>("MuonRPCHits");
108 
109  edm::ParameterSet serviceParameters =
110  iConfig.getParameter<edm::ParameterSet>("ServiceParameters");
111  theService = new MuonServiceProxy(serviceParameters);
112 }
113 
114 // ---- method called once each job just before starting event loop ----
115 void
117 
118  //services
120  edm::ESHandle<DTGeometry> dtGeometry;
121  edm::ESHandle<CSCGeometry> cscGeometry;
122  edm::ESHandle<RPCGeometry> rpcGeometry;
124 
125  es.get<IdealMagneticFieldRecord>().get(magField);
126  es.get<MuonGeometryRecord>().get("MisAligned",dtGeometry);
127  es.get<MuonGeometryRecord>().get("MisAligned",cscGeometry);
128  es.get<MuonGeometryRecord>().get(rpcGeometry);
129 
130  magfield = &(*magField);
131  dtGeom = &(*dtGeometry);
132  cscGeom = &(*cscGeometry);
133  rpcGeom = &(*rpcGeometry);
134 
135  theService->update(es);
136 
137  // A few propagators
138  propagatorWithMaterial = &(*(theService->propagator("SteppingHelixPropagatorAny")));
140  SteppingHelixPropagator* SHpropagator = dynamic_cast<SteppingHelixPropagator*> (propagatorWithoutMaterial); // Beuark!
141  SHpropagator->setMaterialMode(true); // switches OFF material effects;
142 
143 }
144 
146 {
147  // do anything here that needs to be done at destruction time
148  // (e.g. close files, deallocate resources etc.)
149 
150  if ( random ) {
151  delete random;
152  }
153 
156 }
157 
158 
159 //
160 // member functions
161 //
162 
163 // ------------ method called to produce the data ------------
164 
165 void
167  // using namespace edm;
168  // using namespace std;
169 
170  MuonPatternRecoDumper dumper;
171 
174  std::vector<PSimHit> theCSCHits;
175  std::vector<PSimHit> theDTHits;
176  std::vector<PSimHit> theRPCHits;
177 
180  iEvent.getByLabel(theSimModuleLabel_,simVertices);
181 
182  for ( unsigned int itrk=0; itrk<simMuons->size(); itrk++ ) {
183  const SimTrack &mySimTrack = (*simMuons)[itrk];
184  math::XYZTLorentzVector mySimP4(mySimTrack.momentum().x(),
185  mySimTrack.momentum().y(),
186  mySimTrack.momentum().z(),
187  mySimTrack.momentum().t());
188 
189  // Decaying hadrons are now in the list, and so are their muon daughter
190  // Ignore the hadrons here.
191  int pid = mySimTrack.type();
192  if ( abs(pid) != 13 ) continue;
193 
194  double t0 = 0;
195  GlobalPoint initialPosition;
196  int ivert = mySimTrack.vertIndex();
197  if ( ivert >= 0 ) {
198  t0 = (*simVertices)[ivert].position().t();
199  GlobalPoint xyzzy((*simVertices)[ivert].position().x(),
200  (*simVertices)[ivert].position().y(),
201  (*simVertices)[ivert].position().z());
202  initialPosition = xyzzy;
203  }
204 //
205 // Presumably t0 has dimensions of cm if not mm?
206 // Convert to ns for internal calculations.
207 // I wonder where we should get c from?
208 //
209  double tof = t0/29.98;
210 
211 #ifdef FAMOS_DEBUG
212  std::cout << " ===> MuonSimHitProducer::reconstruct() found SIMTRACK - pid = "
213  << pid ;
214  std::cout << " : pT = " << mySimP4.Pt()
215  << ", eta = " << mySimP4.Eta()
216  << ", phi = " << mySimP4.Phi() << std::endl;
217 #endif
218 
219 //
220 // Produce muons sim hits starting from undecayed simulated muons
221 //
222 
223  GlobalPoint startingPosition(mySimTrack.trackerSurfacePosition().x(),
224  mySimTrack.trackerSurfacePosition().y(),
225  mySimTrack.trackerSurfacePosition().z());
226  GlobalVector startingMomentum(mySimTrack.trackerSurfaceMomentum().x(),
227  mySimTrack.trackerSurfaceMomentum().y(),
228  mySimTrack.trackerSurfaceMomentum().z());
229 //
230 // Crap... there's no time-of-flight to the trackerSurfacePosition()...
231 // So, this will be wrong when the curvature can't be neglected, but that
232 // will be rather seldom... May as well ignore the mass too.
233 //
234  GlobalVector dtracker = startingPosition-initialPosition;
235  tof += dtracker.mag()/29.98;
236 
237 #ifdef FAMOS_DEBUG
238  std::cout << " the Muon START position " << startingPosition << std::endl;
239  std::cout << " the Muon START momentum " << startingMomentum << std::endl;
240 #endif
241 
242 //
243 // Some magic to define a TrajectoryStateOnSurface
244 //
245  PlaneBuilder pb;
246  GlobalVector zAxis = startingMomentum.unit();
247  GlobalVector yAxis(zAxis.y(),-zAxis.x(),0);
248  GlobalVector xAxis = yAxis.cross(zAxis);
249  Surface::RotationType rot = Surface::RotationType(xAxis,yAxis,zAxis);
250  PlaneBuilder::ReturnType startingPlane = pb.plane(startingPosition,rot);
251  GlobalTrajectoryParameters gtp(startingPosition,
252  startingMomentum,
253  (int)mySimTrack.charge(),
254  magfield);
255  TrajectoryStateOnSurface startingState(gtp,*startingPlane);
256 
257  std::vector<const DetLayer *> navLayers;
258  if ( fabs(startingState.globalMomentum().eta()) > 4.5 ) {
259  navLayers = navigation.compatibleEndcapLayers(*(startingState.freeState()),
260  alongMomentum);
261  }
262  else {
263  navLayers = navigation.compatibleLayers(*(startingState.freeState()),
264  alongMomentum);
265  }
266  /*
267  edm::ESHandle<Propagator> propagator =
268  theService->propagator("SteppingHelixPropagatorAny");
269  */
270 
271  if ( navLayers.empty() ) continue;
272 
273 #ifdef FAMOS_DEBUG
274  std::cout << "Found " << navLayers.size()
275  << " compatible DetLayers..." << std::endl;
276 #endif
277 
278  TrajectoryStateOnSurface propagatedState = startingState;
279  for ( unsigned int ilayer=0; ilayer<navLayers.size(); ilayer++ ) {
280 
281 #ifdef FAMOS_DEBUG
282  std::cout << "Propagating to layer " << ilayer << " "
283  << dumper.dumpLayer(navLayers[ilayer])
284  << std::endl;
285 #endif
286 
287  std::vector<DetWithState> comps =
288  navLayers[ilayer]->compatibleDets(propagatedState,*propagatorWithMaterial,theEstimator);
289  if ( comps.empty() ) continue;
290 
291 #ifdef FAMOS_DEBUG
292  std::cout << "Propagating " << propagatedState << std::endl;
293 #endif
294 
295  // Starting momentum
296  double pi = propagatedState.globalMomentum().mag();
297 
298  // Propagate with material effects (dE/dx average only)
299  SteppingHelixStateInfo shsStart(*(propagatedState.freeTrajectoryState()));
300  const SteppingHelixStateInfo& shsDest =
301  ((const SteppingHelixPropagator*)propagatorWithMaterial)->propagate(shsStart,navLayers[ilayer]->surface());
302  std::pair<TrajectoryStateOnSurface,double> next(shsDest.getStateOnSurface(navLayers[ilayer]->surface()),
303  shsDest.path());
304  // No need to continue if there is no valid propagation available.
305  // This happens rarely (~0.1% of ttbar events)
306  if ( !next.first.isValid() ) continue;
307  // This is the estimate of the number of radiation lengths traversed,
308  // together with the total path length
309  double radPath = shsDest.radPath();
310  double pathLength = next.second;
311 
312  // Now propagate without dE/dx (average)
313  // [To add the dE/dx fluctuations to the actual dE/dx]
314  std::pair<TrajectoryStateOnSurface,double> nextNoMaterial =
315  propagatorWithoutMaterial->propagateWithPath(propagatedState,navLayers[ilayer]->surface());
316 
317  // Update the propagated state
318  propagatedState = next.first;
319  double pf = propagatedState.globalMomentum().mag();
320 
321  // Insert dE/dx fluctuations and multiple scattering
322  // Skip this step if nextNoMaterial.first is not valid
323  // This happens rarely (~0.02% of ttbar events)
324  if ( theMaterialEffects && nextNoMaterial.first.isValid() ) applyMaterialEffects(propagatedState,nextNoMaterial.first,radPath);
325  // Check that the 'shaken' propagatedState is still valid, otherwise continue
326  if ( !propagatedState.isValid() ) continue;
327  // (No evidence that this ever happens)
328 //
329 // Consider this... 1 GeV muon has a velocity that is only 0.5% slower than c...
330 // We probably can safely ignore the mass for anything that makes it out to the
331 // muon chambers.
332 //
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;
337 
338 #ifdef FAMOS_DEBUG
339  std::cout << "Propagated to next surface... path length = " << pathLength
340  << " cm, dTOF = " << dtof << " ns" << std::endl;
341 #endif
342 
343  tof += dtof;
344 
345  for ( unsigned int icomp=0; icomp<comps.size(); icomp++ )
346  {
347  const GeomDet *gd = comps[icomp].first;
348  if ( gd->subDetector() == GeomDetEnumerators::DT ) {
350  const DTChamber *chamber = dtGeom->chamber(id);
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++ ) {
355  DTLayerId lid = layer[ilayer]->id();
356 #ifdef FAMOS_DEBUG
357  std::cout << " Extrapolated to DT ("
358  << lid.wheel() << ","
359  << lid.station() << ","
360  << lid.sector() << ","
361  << lid.superlayer() << ","
362  << lid.layer() << ")" << std::endl;
363 #endif
364 
365  const GeomDetUnit *det = dtGeom->idToDetUnit(lid);
366 
367  HelixArbitraryPlaneCrossing crossing(propagatedState.globalPosition().basicVector(),
368  propagatedState.globalMomentum().basicVector(),
369  propagatedState.transverseCurvature(),
370  anyDirection);
371  std::pair<bool,double> path = crossing.pathLength(det->surface());
372  if ( ! path.first ) continue;
373  LocalPoint lpos = det->toLocal(GlobalPoint(crossing.position(path.second)));
374  if ( ! det->surface().bounds().inside(lpos) ) continue;
375  const DTTopology& dtTopo = layer[ilayer]->specificTopology();
376  int wire = dtTopo.channel(lpos);
377  if (wire < dtTopo.firstChannel() || wire > dtTopo.lastChannel()) continue;
378  // no drift cell here (on the chamber edge or just outside)
379  // this hit would otherwise be discarded downstream in the digitizer
380 
381  DTWireId wid(lid,wire);
382  double thickness = det->surface().bounds().thickness();
383  LocalVector lmom = det->toLocal(GlobalVector(crossing.direction(path.second)));
384  lmom = lmom.unit()*propagatedState.localMomentum().mag();
385 
386  // Factor that takes into account the (rec)hits lost because of delta's, etc.:
387  // (Not fully satisfactory patch, but it seems to work...)
388  double pmu = lmom.mag();
389  double theDTHitIneff = pmu>0? exp(kDT*log(pmu)+fDT):0.;
390  if (random->flatShoot()<theDTHitIneff) continue;
391 
392  double eloss = 0;
393  double pz = fabs(lmom.z());
394  LocalPoint entry = lpos - 0.5*thickness*lmom/pz;
395  LocalPoint exit = lpos + 0.5*thickness*lmom/pz;
396  double dtof = path.second*rbeta;
397  int trkid = mySimTrack.trackId();
398  unsigned int id = wid.rawId();
399  short unsigned int processType = 2;
400  PSimHit hit(entry,exit,lmom.mag(),
401  tof+dtof,eloss,pid,id,trkid,lmom.theta(),lmom.phi(),processType);
402  theDTHits.push_back(hit);
403 
404  }
405  }
406  }
407  else if ( gd->subDetector() == GeomDetEnumerators::CSC ) {
408  CSCDetId id(gd->geographicalId());
409  const CSCChamber *chamber = cscGeom->chamber(id);
410  std::vector<const CSCLayer *> layers = chamber->layers();
411  for ( unsigned int ilayer=0; ilayer<layers.size(); ilayer++ ) {
412  CSCDetId lid = layers[ilayer]->id();
413 
414 #ifdef FAMOS_DEBUG
415  std::cout << " Extrapolated to CSC ("
416  << lid.endcap() << ","
417  << lid.ring() << ","
418  << lid.station() << ","
419  << lid.layer() << ")" << std::endl;
420 #endif
421 
422  const GeomDetUnit *det = cscGeom->idToDetUnit(lid);
423  HelixArbitraryPlaneCrossing crossing(propagatedState.globalPosition().basicVector(),
424  propagatedState.globalMomentum().basicVector(),
425  propagatedState.transverseCurvature(),
426  anyDirection);
427  std::pair<bool,double> path = crossing.pathLength(det->surface());
428  if ( ! path.first ) continue;
429  LocalPoint lpos = det->toLocal(GlobalPoint(crossing.position(path.second)));
430  // For CSCs the Bounds are for chamber frames not gas regions
431  // if ( ! det->surface().bounds().inside(lpos) ) continue;
432  // New function knows where the 'active' volume is:
433  const CSCLayerGeometry* laygeom = layers[ilayer]->geometry();
434  if ( ! laygeom->inside( lpos ) ) continue;
435  //double thickness = laygeom->thickness(); gives number which is about 20 times too big
436  double thickness = det->surface().bounds().thickness(); // this one works much better...
437  LocalVector lmom = det->toLocal(GlobalVector(crossing.direction(path.second)));
438  lmom = lmom.unit()*propagatedState.localMomentum().mag();
439 
440  // Factor that takes into account the (rec)hits lost because of delta's, etc.:
441  // (Not fully satisfactory patch, but it seems to work...)
442  double pmu = lmom.mag();
443  double theCSCHitIneff = pmu>0? exp(kCSC*log(pmu)+fCSC):0.;
444  // Take into account the different geometry in ME11:
445  if (id.station()==1 && id.ring()==1) theCSCHitIneff = theCSCHitIneff*0.442;
446  if (random->flatShoot()<theCSCHitIneff) continue;
447 
448  double eloss = 0;
449  double pz = fabs(lmom.z());
450  LocalPoint entry = lpos - 0.5*thickness*lmom/pz;
451  LocalPoint exit = lpos + 0.5*thickness*lmom/pz;
452  double dtof = path.second*rbeta;
453  int trkid = mySimTrack.trackId();
454  unsigned int id = lid.rawId();
455  short unsigned int processType = 2;
456  PSimHit hit(entry,exit,lmom.mag(),
457  tof+dtof,eloss,pid,id,trkid,lmom.theta(),lmom.phi(),processType);
458  theCSCHits.push_back(hit);
459  }
460  }
461  else if ( gd->subDetector() == GeomDetEnumerators::RPCBarrel ||
463  RPCDetId id(gd->geographicalId());
464  const RPCChamber *chamber = rpcGeom->chamber(id);
465  std::vector<const RPCRoll *> roll = chamber->rolls();
466  for ( unsigned int iroll=0; iroll<roll.size(); iroll++ ) {
467  RPCDetId rid = roll[iroll]->id();
468 
469 #ifdef FAMOS_DEBUG
470  std::cout << " Extrapolated to RPC ("
471  << rid.ring() << ","
472  << rid.station() << ","
473  << rid.sector() << ","
474  << rid.subsector() << ","
475  << rid.layer() << ","
476  << rid.roll() << ")" << std::endl;
477 #endif
478 
479  const GeomDetUnit *det = rpcGeom->idToDetUnit(rid);
480  HelixArbitraryPlaneCrossing crossing(propagatedState.globalPosition().basicVector(),
481  propagatedState.globalMomentum().basicVector(),
482  propagatedState.transverseCurvature(),
483  anyDirection);
484  std::pair<bool,double> path = crossing.pathLength(det->surface());
485  if ( ! path.first ) continue;
486  LocalPoint lpos = det->toLocal(GlobalPoint(crossing.position(path.second)));
487  if ( ! det->surface().bounds().inside(lpos) ) continue;
488  double thickness = det->surface().bounds().thickness();
489  LocalVector lmom = det->toLocal(GlobalVector(crossing.direction(path.second)));
490  lmom = lmom.unit()*propagatedState.localMomentum().mag();
491  double eloss = 0;
492  double pz = fabs(lmom.z());
493  LocalPoint entry = lpos - 0.5*thickness*lmom/pz;
494  LocalPoint exit = lpos + 0.5*thickness*lmom/pz;
495  double dtof = path.second*rbeta;
496  int trkid = mySimTrack.trackId();
497  unsigned int id = rid.rawId();
498  short unsigned int processType = 2;
499  PSimHit hit(entry,exit,lmom.mag(),
500  tof+dtof,eloss,pid,id,trkid,lmom.theta(),lmom.phi(),processType);
501  theRPCHits.push_back(hit);
502  }
503  }
504  else {
505  std::cout << "Extrapolated to unknown subdetector '" << gd->subDetector() << "'..." << std::endl;
506  }
507  }
508  }
509  }
510 
511  std::auto_ptr<edm::PSimHitContainer> pcsc(new edm::PSimHitContainer);
512  int n = 0;
513  for ( std::vector<PSimHit>::const_iterator i = theCSCHits.begin();
514  i != theCSCHits.end(); i++ ) {
515  pcsc->push_back(*i);
516  n += 1;
517  }
518  iEvent.put(pcsc,"MuonCSCHits");
519 
520  std::auto_ptr<edm::PSimHitContainer> pdt(new edm::PSimHitContainer);
521  n = 0;
522  for ( std::vector<PSimHit>::const_iterator i = theDTHits.begin();
523  i != theDTHits.end(); i++ ) {
524  pdt->push_back(*i);
525  n += 1;
526  }
527  iEvent.put(pdt,"MuonDTHits");
528 
529  std::auto_ptr<edm::PSimHitContainer> prpc(new edm::PSimHitContainer);
530  n = 0;
531  for ( std::vector<PSimHit>::const_iterator i = theRPCHits.begin();
532  i != theRPCHits.end(); i++ ) {
533  prpc->push_back(*i);
534  n += 1;
535  }
536  iEvent.put(prpc,"MuonRPCHits");
537 
538 }
539 
540 
541 // ------------ method called once each job just after ending the event loop ------------
542 void
544 {
545 }
546 
547 
548 void
550  const edm::ParameterSet& fastTracks,
551  const edm::ParameterSet& matEff) {
552  // Muons
553  theSimModuleLabel_ = fastMuons.getParameter<std::string>("simModuleLabel");
554  theSimModuleProcess_ = fastMuons.getParameter<std::string>("simModuleProcess");
555  theTrkModuleLabel_ = fastMuons.getParameter<std::string>("trackModuleLabel");
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];
562 
563  // Tracks
564  fullPattern_ = fastTracks.getUntrackedParameter<bool>("FullPatternRecognition");
565 
566  // The following should be on LogInfo
567  // std::cout << " MUON SIM HITS: FastSimulation parameters " << std::endl;
568  // std::cout << " ============================================== " << std::endl;
569  // if ( fullPattern_ )
570  // std::cout << " The FULL pattern recognition option is turned ON" << std::endl;
571  // else
572  // std::cout << " The FAST tracking option is turned ON" << std::endl;
573 
574  // Material Effects
575  theMaterialEffects = 0;
576  if ( matEff.getParameter<bool>("PairProduction") ||
577  matEff.getParameter<bool>("Bremsstrahlung") ||
578  matEff.getParameter<bool>("EnergyLoss") ||
579  matEff.getParameter<bool>("MultipleScattering") )
580  theMaterialEffects = new MaterialEffects(matEff,random);
581 
582 }
583 
584 void
587  double radPath) {
588 
589  // The energy loss simulator
591 
592  // The multiple scattering simulator
594 
595  // Initialize the Particle position, momentum and energy
596  const Surface& nextSurface = tsos.surface();
597  GlobalPoint gPos = energyLoss ? tsos.globalPosition() : tsosWithdEdx.globalPosition();
598  GlobalVector gMom = energyLoss ? tsos.globalMomentum() : tsosWithdEdx.globalMomentum();
599  double mu = 0.1056583692;
600  double en = std::sqrt(gMom.mag2()+mu*mu);
601 
602  // And now create the Particle
603  XYZTLorentzVector position(gPos.x(),gPos.y(),gPos.z(),0.);
604  XYZTLorentzVector momentum(gMom.x(),gMom.y(),gMom.z(),en);
605  float charge = (float)(tsos.charge());
606  ParticlePropagator theMuon(momentum,position,charge,0);
607  theMuon.setID(-(int)charge*13);
608 
609  // Recompute the energy loss to get the fluctuations
610  if ( energyLoss ) {
611  // Difference between with and without dE/dx (average only)
612  // (for corrections once fluctuations are applied)
613  GlobalPoint gPosWithdEdx = tsosWithdEdx.globalPosition();
614  GlobalVector gMomWithdEdx = tsosWithdEdx.globalMomentum();
615  double enWithdEdx = std::sqrt(gMomWithdEdx.mag2()+mu*mu);
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);
622 
623  // Energy loss in iron (+ fluctuations)
624  energyLoss->updateState(theMuon,radPath);
625 
626  // Correcting factors to account for slight differences in material descriptions
627  // (Material description is more accurate in the stepping helix propagator)
628  radPath *= -deltaMom.E()/energyLoss->mostLikelyLoss();
629  double fac = energyLoss->deltaMom().E()/energyLoss->mostLikelyLoss();
630 
631  // Particle momentum & position after energy loss + fluctuation
632  XYZTLorentzVector theNewMomentum = theMuon.momentum() + energyLoss->deltaMom() + fac * deltaMom;
633  XYZTLorentzVector theNewPosition = theMuon.vertex() + fac * deltaPos;
634  fac = (theNewMomentum.E()*theNewMomentum.E()-mu*mu)/theNewMomentum.Vect().Mag2();
635  fac = fac>0.? std::sqrt(fac) : 1E-9;
636  theMuon.SetXYZT(theNewMomentum.Px()*fac,theNewMomentum.Py()*fac,
637  theNewMomentum.Pz()*fac,theNewMomentum.E());
638  theMuon.setVertex(theNewPosition);
639 
640  }
641 
642  // Does the actual mutliple scattering
643  if ( multipleScattering ) {
644  // Pass the vector normal to the "next" surface
645  GlobalVector normal = nextSurface.tangentPlane(tsos.globalPosition())->normalVector();
646  multipleScattering->setNormalVector(normal);
647  // Compute the amount of multiple scattering after a given path length
648  multipleScattering->updateState(theMuon,radPath);
649  }
650 
651  // Fill the propagated state
652  GlobalPoint propagatedPosition(theMuon.X(),theMuon.Y(),theMuon.Z());
653  GlobalVector propagatedMomentum(theMuon.Px(),theMuon.Py(),theMuon.Pz());
654  GlobalTrajectoryParameters propagatedGtp(propagatedPosition,propagatedMomentum,(int)charge,magfield);
655  tsosWithdEdx = TrajectoryStateOnSurface(propagatedGtp,nextSurface);
656 
657 }
658 
659 
660 
661 //define this as a plug-in
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.
Definition: DTTopology.h:82
const math::XYZVectorD & trackerSurfacePosition() const
Definition: SimTrack.h:36
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
const MagneticField * magfield
int i
Definition: DBlmapReader.cc:9
T mag2() const
Definition: PV3DBase.h:60
const std::vector< const CSCLayer * > & layers() const
Return all layers.
Definition: CSCChamber.h:57
virtual void beginRun(edm::Run &run, const edm::EventSetup &es)
const DTGeometry * dtGeom
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.
Definition: RPCGeometry.cc:47
void setNormalVector(const GlobalVector &normal)
Sets the vector normal to the surface traversed.
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
FreeTrajectoryState * freeTrajectoryState(bool withErrors=true) const
int channel(const LocalPoint &p) const
Definition: DTTopology.cc:68
MaterialEffects * theMaterialEffects
Geom::Phi< T > phi() const
Definition: PV3DBase.h:63
MultipleScatteringSimulator * multipleScatteringSimulator() const
Return the Multiple Scattering engine.
ReturnType plane(const PositionType &pos, const RotationType &rot) const
Definition: PlaneBuilder.h:22
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
T y() const
Definition: PV3DBase.h:57
virtual SubDetector subDetector() const =0
Which subdetector.
virtual Propagator * clone() const =0
#define abs(x)
Definition: mlp_lapack.h:159
GlobalPoint globalPosition() const
LocalPoint toLocal(const GlobalPoint &gp) const
Conversion to the R.F. of the GeomDet.
Definition: GeomDet.h:64
Exp< T >::type exp(const T &t)
Definition: Exp.h:22
int layer() const
Return the layer number.
Definition: DTLayerId.h:55
const int firstChannel() const
Returns the wire number of the first wire.
Definition: DTTopology.h:80
float charge() const
charge
Definition: CoreSimTrack.cc:3
int layer() const
Definition: CSCDetId.h:63
double charge(const std::vector< uint8_t > &Ampls)
static int position[TOTALCHAMBERS][3]
Definition: ReadPGInfo.cc:509
int path() const
Definition: HLTadd.h:3
Geom::Theta< T > theta() const
Definition: PV3DBase.h:69
const RPCGeometry * rpcGeom
uint32_t rawId() const
get the raw id
Definition: DetId.h:45
LocalVector localMomentum() const
int endcap() const
Definition: CSCDetId.h:95
virtual float thickness() const =0
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:30
const RPCChamber * chamber(RPCDetId id) const
Definition: RPCGeometry.cc:66
virtual std::pair< TrajectoryStateOnSurface, double > propagateWithPath(const FreeTrajectoryState &, const Surface &) const
Definition: Propagator.cc:74
int iEvent
Definition: GenABIO.cc:243
const XYZTLorentzVector & momentum() const
the momentum fourvector
Definition: RawParticle.h:285
T mag() const
Definition: PV3DBase.h:61
int roll() const
Definition: RPCDetId.h:122
Definition: DDAxes.h:10
std::string theSimModuleProcess_
virtual const GeomDetUnit * idToDetUnit(DetId) const
Return the pointer to the GeomDetUnit corresponding to a given DetId.
Definition: CSCGeometry.cc:87
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:84
const std::vector< const DTSuperLayer * > & superLayers() const
Return the superlayers in the chamber.
Definition: DTChamber.cc:62
int ring() const
Definition: RPCDetId.h:74
const Propagator * propagatorWithMaterial
TrajectoryStateOnSurface getStateOnSurface(const Surface &surf, bool returnTangentPlane=false) const
T sqrt(T t)
Definition: SSEVec.h:28
Propagator * propagatorWithoutMaterial
std::pair< std::string, MonitorElement * > entry
Definition: ME_MAP.h:8
Vector3DBase< typename PreciseFloatType< T, U >::Type, FrameTag > cross(const Vector3DBase< U, FrameTag > &v) const
Definition: Vector3DBase.h:119
T z() const
Definition: PV3DBase.h:58
double Y() const
y of vertex
Definition: RawParticle.h:274
void setID(const int id)
Definition: RawParticle.cc:101
const XYZTLorentzVector & deltaMom() const
Returns the actual energy lost.
bool isAvailable() const
Definition: Service.h:47
double Z() const
z of vertex
Definition: RawParticle.h:275
const std::vector< const RPCRoll * > & rolls() const
Return the Rolls.
Definition: RPCChamber.cc:70
const CSCGeometry * cscGeom
DetId geographicalId() const
The label of this GeomDet.
Definition: GeomDet.h:74
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:359
int vertIndex() const
index of the vertex in the Event container (-1 if no vertex)
Definition: SimTrack.h:29
const XYZTLorentzVector & vertex() const
the vertex fourvector
Definition: RawParticle.h:284
int superlayer() const
Return the superlayer number (deprecated method name)
int ring() const
Definition: CSCDetId.h:77
Vector3DBase unit() const
Definition: Vector3DBase.h:57
int layer() const
Definition: RPCDetId.h:110
const Bounds & bounds() const
Definition: BoundSurface.h:89
virtual ReferenceCountingPointer< TangentPlane > tangentPlane(const GlobalPoint &) const =0
unsigned int trackId() const
Definition: CoreSimTrack.h:49
void updateState(ParticlePropagator &myTrack, double radlen)
Compute the material effect (calls the sub class)
Log< T >::type log(const T &t)
Definition: Log.h:22
const math::XYZTLorentzVectorD & trackerSurfaceMomentum() const
Definition: SimTrack.h:38
const T & get() const
Definition: EventSetup.h:55
const CSCChamber * chamber(CSCDetId id) const
Return the chamber corresponding to given DetId.
Definition: CSCGeometry.cc:117
const RandomEngine * random
const DTChamber * chamber(DTChamberId id) const
Return a DTChamber given its id.
Definition: DTGeometry.cc:105
double flatShoot(double xmin=0.0, double xmax=1.0) const
Definition: RandomEngine.h:30
int sector() const
Sector id: the group of chambers at same phi (and increasing r)
Definition: RPCDetId.h:104
double X() const
x of vertex
Definition: RawParticle.h:273
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...
Definition: RPCDetId.h:116
int type() const
particle type (HEP PDT convension)
Definition: CoreSimTrack.h:40
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_
int sector() const
Definition: DTChamberId.h:63
const math::XYZTLorentzVectorD & momentum() const
particle info...
Definition: CoreSimTrack.h:36
TkRotation< float > RotationType
const Surface & surface() const
virtual const GeomDetUnit * idToDetUnit(DetId) const
Return the pointer to the GeomDetUnit corresponding to a given DetId.
Definition: DTGeometry.cc:76
bool inside(const Local3DPoint &, const LocalError &, float scale=1.f) const
int station() const
Definition: CSCDetId.h:88
tuple cout
Definition: gather_cfg.py:41
std::vector< PSimHit > PSimHitContainer
double pi
int station() const
Return the station number.
Definition: DTChamberId.h:53
EnergyLossSimulator * energyLossSimulator() const
Return the Energy Loss engine.
int wheel() const
Return the wheel number.
Definition: DTChamberId.h:47
T x() const
Definition: PV3DBase.h:56
const BasicVectorType & basicVector() const
Definition: PV3DBase.h:54
virtual const BoundPlane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:37
void setMaterialMode(bool noMaterial)
Switch for material effects mode: no material effects if true.
void setVertex(const XYZTLorentzVector &vtx)
set the vertex
Definition: RawParticle.h:287
edm::ESHandle< Propagator > propagator(std::string propagatorName) const
get the propagator
Definition: Run.h:31
Global3DVector GlobalVector
Definition: GlobalVector.h:10
math::XYZTLorentzVector XYZTLorentzVector
Definition: RawParticle.h:15
Chi2MeasurementEstimator theEstimator
int station() const
Definition: RPCDetId.h:98