File indexing completed on 2023-03-17 11:00:50
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023 #include "DataFormats/GeometrySurface/interface/PlaneBuilder.h"
0024 #include "DataFormats/GeometrySurface/interface/TangentPlane.h"
0025 #include "DataFormats/MuonDetId/interface/DTWireId.h"
0026 #include "FastSimulation/MaterialEffects/interface/EnergyLossSimulator.h"
0027 #include "FastSimulation/MaterialEffects/interface/MaterialEffects.h"
0028 #include "FastSimulation/MaterialEffects/interface/MultipleScatteringSimulator.h"
0029 #include "FastSimulation/MaterialEffects/interface/MuonBremsstrahlungSimulator.h"
0030 #include "FastSimulation/ParticlePropagator/interface/ParticlePropagator.h"
0031 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0032 #include "FWCore/Framework/interface/ConsumesCollector.h"
0033 #include "FWCore/Framework/interface/stream/EDProducer.h"
0034 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0035 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0036 #include "Geometry/RPCGeometry/interface/RPCGeometry.h"
0037 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0038 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0039 #include "RecoMuon/MeasurementDet/interface/MuonDetLayerMeasurements.h"
0040 #include "RecoMuon/Navigation/interface/DirectMuonNavigation.h"
0041 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0042 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0043 #include "SimDataFormats/Track/interface/SimTrack.h"
0044 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0045 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0046 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0047 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
0048 #include "TrackingTools/GeomPropagators/interface/HelixArbitraryPlaneCrossing.h"
0049 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0050
0051 class MuonSimHitProducer : public edm::stream::EDProducer<> {
0052 public:
0053 explicit MuonSimHitProducer(const edm::ParameterSet&);
0054
0055 private:
0056 MuonServiceProxy* theService;
0057 Chi2MeasurementEstimator theEstimator;
0058
0059 const MagneticField* magfield;
0060 const DTGeometry* dtGeom;
0061 const CSCGeometry* cscGeom;
0062 const RPCGeometry* rpcGeom;
0063 const Propagator* propagatorWithMaterial;
0064 std::unique_ptr<Propagator> propagatorWithoutMaterial;
0065
0066 std::unique_ptr<MaterialEffects> theMaterialEffects;
0067
0068 void beginRun(edm::Run const& run, const edm::EventSetup& es) override;
0069 void produce(edm::Event&, const edm::EventSetup&) override;
0070 void readParameters(const edm::ParameterSet&, const edm::ParameterSet&, const edm::ParameterSet&);
0071
0072
0073 double kDT;
0074 double fDT;
0075 double kCSC;
0076 double fCSC;
0077
0078
0079 void applyMaterialEffects(TrajectoryStateOnSurface& tsosWithdEdx,
0080 TrajectoryStateOnSurface& tsos,
0081 double radPath,
0082 RandomEngineAndDistribution const*,
0083 HepPDT::ParticleDataTable const&);
0084
0085
0086 bool fullPattern_;
0087 bool doL1_, doL3_, doGL_;
0088
0089
0090 edm::InputTag simMuonLabel;
0091 edm::InputTag simVertexLabel;
0092
0093
0094 edm::EDGetTokenT<std::vector<SimTrack> > simMuonToken;
0095 edm::EDGetTokenT<std::vector<SimVertex> > simVertexToken;
0096
0097 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magneticFieldESToken_;
0098 const edm::ESGetToken<DTGeometry, MuonGeometryRecord> dtGeometryESToken_;
0099 const edm::ESGetToken<CSCGeometry, MuonGeometryRecord> cscGeometryESToken_;
0100 const edm::ESGetToken<RPCGeometry, MuonGeometryRecord> rpcGeometryESToken_;
0101 const edm::ESGetToken<HepPDT::ParticleDataTable, edm::DefaultRecord> particleDataTableESToken_;
0102 };
0103
0104
0105
0106
0107
0108
0109
0110 MuonSimHitProducer::MuonSimHitProducer(const edm::ParameterSet& iConfig)
0111 : theEstimator(iConfig.getParameter<double>("Chi2EstimatorCut")),
0112 propagatorWithoutMaterial(nullptr),
0113 magneticFieldESToken_(esConsumes<edm::Transition::BeginRun>()),
0114 dtGeometryESToken_(esConsumes<edm::Transition::BeginRun>(edm::ESInputTag("", "MisAligned"))),
0115 cscGeometryESToken_(esConsumes<edm::Transition::BeginRun>(edm::ESInputTag("", "MisAligned"))),
0116 rpcGeometryESToken_(esConsumes<edm::Transition::BeginRun>()),
0117 particleDataTableESToken_(esConsumes()) {
0118
0119 readParameters(iConfig.getParameter<edm::ParameterSet>("MUONS"),
0120 iConfig.getParameter<edm::ParameterSet>("TRACKS"),
0121 iConfig.getParameter<edm::ParameterSet>("MaterialEffectsForMuons"));
0122
0123
0124
0125
0126 produces<edm::PSimHitContainer>("MuonCSCHits");
0127 produces<edm::PSimHitContainer>("MuonDTHits");
0128 produces<edm::PSimHitContainer>("MuonRPCHits");
0129
0130 edm::ParameterSet serviceParameters = iConfig.getParameter<edm::ParameterSet>("ServiceParameters");
0131 theService = new MuonServiceProxy(serviceParameters, consumesCollector(), MuonServiceProxy::UseEventSetupIn::Run);
0132
0133
0134 simMuonToken = consumes<std::vector<SimTrack> >(simMuonLabel);
0135 simVertexToken = consumes<std::vector<SimVertex> >(simVertexLabel);
0136 }
0137
0138
0139 void MuonSimHitProducer::beginRun(edm::Run const& run, const edm::EventSetup& es) {
0140
0141 magfield = &es.getData(magneticFieldESToken_);
0142 dtGeom = &es.getData(dtGeometryESToken_);
0143 cscGeom = &es.getData(cscGeometryESToken_);
0144 rpcGeom = &es.getData(rpcGeometryESToken_);
0145
0146 bool duringEvent = false;
0147 theService->update(es, duringEvent);
0148
0149
0150 propagatorWithMaterial = &(*(theService->propagator("SteppingHelixPropagatorAny")));
0151 propagatorWithoutMaterial.reset(propagatorWithMaterial->clone());
0152 SteppingHelixPropagator* SHpropagator =
0153 dynamic_cast<SteppingHelixPropagator*>(propagatorWithoutMaterial.get());
0154 SHpropagator->setMaterialMode(true);
0155 }
0156
0157
0158
0159
0160
0161
0162
0163 void MuonSimHitProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0164 auto pdg = iSetup.getHandle(particleDataTableESToken_);
0165
0166 RandomEngineAndDistribution random(iEvent.streamID());
0167
0168 MuonPatternRecoDumper dumper;
0169
0170 edm::Handle<std::vector<SimTrack> > simMuons;
0171 edm::Handle<std::vector<SimVertex> > simVertices;
0172 std::vector<PSimHit> theCSCHits;
0173 std::vector<PSimHit> theDTHits;
0174 std::vector<PSimHit> theRPCHits;
0175
0176 DirectMuonNavigation navigation(theService->detLayerGeometry());
0177 iEvent.getByToken(simMuonToken, simMuons);
0178 iEvent.getByToken(simVertexToken, simVertices);
0179
0180 for (unsigned int itrk = 0; itrk < simMuons->size(); itrk++) {
0181 const SimTrack& mySimTrack = (*simMuons)[itrk];
0182 math::XYZTLorentzVector mySimP4(
0183 mySimTrack.momentum().x(), mySimTrack.momentum().y(), mySimTrack.momentum().z(), mySimTrack.momentum().t());
0184
0185
0186
0187 int pid = mySimTrack.type();
0188 if (abs(pid) != 13 && abs(pid) != 1000024)
0189 continue;
0190 double t0 = 0;
0191 GlobalPoint initialPosition;
0192 int ivert = mySimTrack.vertIndex();
0193 if (ivert >= 0) {
0194 t0 = (*simVertices)[ivert].position().t();
0195 GlobalPoint xyzzy((*simVertices)[ivert].position().x(),
0196 (*simVertices)[ivert].position().y(),
0197 (*simVertices)[ivert].position().z());
0198 initialPosition = xyzzy;
0199 }
0200
0201
0202
0203
0204
0205 double tof = t0 / 29.98;
0206
0207 #ifdef FAMOS_DEBUG
0208 std::cout << " ===> MuonSimHitProducer::reconstruct() found SIMTRACK - pid = " << pid;
0209 std::cout << " : pT = " << mySimP4.Pt() << ", eta = " << mySimP4.Eta() << ", phi = " << mySimP4.Phi() << std::endl;
0210 #endif
0211
0212
0213
0214
0215
0216 GlobalPoint startingPosition(mySimTrack.trackerSurfacePosition().x(),
0217 mySimTrack.trackerSurfacePosition().y(),
0218 mySimTrack.trackerSurfacePosition().z());
0219 GlobalVector startingMomentum(mySimTrack.trackerSurfaceMomentum().x(),
0220 mySimTrack.trackerSurfaceMomentum().y(),
0221 mySimTrack.trackerSurfaceMomentum().z());
0222
0223
0224
0225
0226
0227 GlobalVector dtracker = startingPosition - initialPosition;
0228 tof += dtracker.mag() / 29.98;
0229
0230 #ifdef FAMOS_DEBUG
0231 std::cout << " the Muon START position " << startingPosition << std::endl;
0232 std::cout << " the Muon START momentum " << startingMomentum << std::endl;
0233 #endif
0234
0235
0236
0237
0238 PlaneBuilder pb;
0239 GlobalVector zAxis = startingMomentum.unit();
0240 GlobalVector yAxis(zAxis.y(), -zAxis.x(), 0);
0241 GlobalVector xAxis = yAxis.cross(zAxis);
0242 Surface::RotationType rot = Surface::RotationType(xAxis, yAxis, zAxis);
0243 PlaneBuilder::ReturnType startingPlane = pb.plane(startingPosition, rot);
0244 GlobalTrajectoryParameters gtp(startingPosition, startingMomentum, (int)mySimTrack.charge(), magfield);
0245 TrajectoryStateOnSurface startingState(gtp, *startingPlane);
0246
0247 std::vector<const DetLayer*> navLayers;
0248 if (fabs(startingState.globalMomentum().eta()) > 4.5) {
0249 navLayers = navigation.compatibleEndcapLayers(*(startingState.freeState()), alongMomentum);
0250 } else {
0251 navLayers = navigation.compatibleLayers(*(startingState.freeState()), alongMomentum);
0252 }
0253
0254
0255
0256
0257
0258 if (navLayers.empty())
0259 continue;
0260
0261 #ifdef FAMOS_DEBUG
0262 std::cout << "Found " << navLayers.size() << " compatible DetLayers..." << std::endl;
0263 #endif
0264
0265 TrajectoryStateOnSurface propagatedState = startingState;
0266 for (unsigned int ilayer = 0; ilayer < navLayers.size(); ilayer++) {
0267 #ifdef FAMOS_DEBUG
0268 std::cout << "Propagating to layer " << ilayer << " " << dumper.dumpLayer(navLayers[ilayer]) << std::endl;
0269 #endif
0270
0271 std::vector<DetWithState> comps =
0272 navLayers[ilayer]->compatibleDets(propagatedState, *propagatorWithMaterial, theEstimator);
0273 if (comps.empty())
0274 continue;
0275
0276 #ifdef FAMOS_DEBUG
0277 std::cout << "Propagating " << propagatedState << std::endl;
0278 #endif
0279
0280
0281 double pi = propagatedState.globalMomentum().mag();
0282
0283
0284 SteppingHelixStateInfo shsStart(*(propagatedState.freeTrajectoryState()));
0285 SteppingHelixStateInfo shsDest;
0286 ((const SteppingHelixPropagator*)propagatorWithMaterial)
0287 ->propagate(shsStart, navLayers[ilayer]->surface(), shsDest);
0288 std::pair<TrajectoryStateOnSurface, double> next(shsDest.getStateOnSurface(navLayers[ilayer]->surface()),
0289 shsDest.path());
0290
0291
0292 if (!next.first.isValid())
0293 continue;
0294
0295
0296 double radPath = shsDest.radPath();
0297 double pathLength = next.second;
0298
0299
0300
0301 std::pair<TrajectoryStateOnSurface, double> nextNoMaterial =
0302 propagatorWithoutMaterial->propagateWithPath(propagatedState, navLayers[ilayer]->surface());
0303
0304
0305 propagatedState = next.first;
0306 double pf = propagatedState.globalMomentum().mag();
0307
0308
0309
0310
0311 if (theMaterialEffects && nextNoMaterial.first.isValid())
0312 applyMaterialEffects(propagatedState, nextNoMaterial.first, radPath, &random, *pdg);
0313
0314 if (!propagatedState.isValid())
0315 continue;
0316
0317
0318
0319
0320
0321
0322 double pavg = 0.5 * (pi + pf);
0323 double m = mySimP4.M();
0324 double rbeta = sqrt(1 + m * m / (pavg * pavg)) / 29.98;
0325 double dtof = pathLength * rbeta;
0326
0327 #ifdef FAMOS_DEBUG
0328 std::cout << "Propagated to next surface... path length = " << pathLength << " cm, dTOF = " << dtof << " ns"
0329 << std::endl;
0330 #endif
0331
0332 tof += dtof;
0333
0334 for (unsigned int icomp = 0; icomp < comps.size(); icomp++) {
0335 const GeomDet* gd = comps[icomp].first;
0336 if (gd->subDetector() == GeomDetEnumerators::DT) {
0337 DTChamberId id(gd->geographicalId());
0338 const DTChamber* chamber = dtGeom->chamber(id);
0339 std::vector<const DTSuperLayer*> superlayer = chamber->superLayers();
0340 for (unsigned int isl = 0; isl < superlayer.size(); isl++) {
0341 std::vector<const DTLayer*> layer = superlayer[isl]->layers();
0342 for (unsigned int ilayer = 0; ilayer < layer.size(); ilayer++) {
0343 DTLayerId lid = layer[ilayer]->id();
0344 #ifdef FAMOS_DEBUG
0345 std::cout << " Extrapolated to DT (" << lid.wheel() << "," << lid.station() << "," << lid.sector()
0346 << "," << lid.superlayer() << "," << lid.layer() << ")" << std::endl;
0347 #endif
0348
0349 const GeomDetUnit* det = dtGeom->idToDetUnit(lid);
0350
0351 HelixArbitraryPlaneCrossing crossing(propagatedState.globalPosition().basicVector(),
0352 propagatedState.globalMomentum().basicVector(),
0353 propagatedState.transverseCurvature(),
0354 anyDirection);
0355 std::pair<bool, double> path = crossing.pathLength(det->surface());
0356 if (!path.first)
0357 continue;
0358 LocalPoint lpos = det->toLocal(GlobalPoint(crossing.position(path.second)));
0359 if (!det->surface().bounds().inside(lpos))
0360 continue;
0361 const DTTopology& dtTopo = layer[ilayer]->specificTopology();
0362 int wire = dtTopo.channel(lpos);
0363 if (wire - dtTopo.firstChannel() < 0 || wire - dtTopo.lastChannel() > 0)
0364 continue;
0365
0366
0367
0368 DTWireId wid(lid, wire);
0369 double thickness = det->surface().bounds().thickness();
0370 LocalVector lmom = det->toLocal(GlobalVector(crossing.direction(path.second)));
0371 lmom = lmom.unit() * propagatedState.localMomentum().mag();
0372
0373
0374
0375 double pmu = lmom.mag();
0376 double theDTHitIneff = pmu > 0 ? exp(kDT * log(pmu) + fDT) : 0.;
0377 if (random.flatShoot() < theDTHitIneff)
0378 continue;
0379
0380 double eloss = 0;
0381 double pz = fabs(lmom.z());
0382 LocalPoint entry = lpos - 0.5 * thickness * lmom / pz;
0383 LocalPoint exit = lpos + 0.5 * thickness * lmom / pz;
0384 double dtof = path.second * rbeta;
0385 int trkid = mySimTrack.trackId();
0386 unsigned int id = wid.rawId();
0387 short unsigned int processType = 2;
0388 PSimHit hit(
0389 entry, exit, lmom.mag(), tof + dtof, eloss, pid, id, trkid, lmom.theta(), lmom.phi(), processType);
0390 theDTHits.push_back(hit);
0391 }
0392 }
0393 } else if (gd->subDetector() == GeomDetEnumerators::CSC) {
0394 CSCDetId id(gd->geographicalId());
0395 const CSCChamber* chamber = cscGeom->chamber(id);
0396 std::vector<const CSCLayer*> layers = chamber->layers();
0397 for (unsigned int ilayer = 0; ilayer < layers.size(); ilayer++) {
0398 CSCDetId lid = layers[ilayer]->id();
0399
0400 #ifdef FAMOS_DEBUG
0401 std::cout << " Extrapolated to CSC (" << lid.endcap() << "," << lid.ring() << "," << lid.station() << ","
0402 << lid.layer() << ")" << std::endl;
0403 #endif
0404
0405 const GeomDetUnit* det = cscGeom->idToDetUnit(lid);
0406 HelixArbitraryPlaneCrossing crossing(propagatedState.globalPosition().basicVector(),
0407 propagatedState.globalMomentum().basicVector(),
0408 propagatedState.transverseCurvature(),
0409 anyDirection);
0410 std::pair<bool, double> path = crossing.pathLength(det->surface());
0411 if (!path.first)
0412 continue;
0413 LocalPoint lpos = det->toLocal(GlobalPoint(crossing.position(path.second)));
0414
0415
0416
0417 const CSCLayerGeometry* laygeom = layers[ilayer]->geometry();
0418 if (!laygeom->inside(lpos))
0419 continue;
0420
0421 double thickness = det->surface().bounds().thickness();
0422 LocalVector lmom = det->toLocal(GlobalVector(crossing.direction(path.second)));
0423 lmom = lmom.unit() * propagatedState.localMomentum().mag();
0424
0425
0426
0427 double pmu = lmom.mag();
0428 double theCSCHitIneff = pmu > 0 ? exp(kCSC * log(pmu) + fCSC) : 0.;
0429
0430 if (id.station() == 1 && id.ring() == 1)
0431 theCSCHitIneff = theCSCHitIneff * 0.442;
0432 if (random.flatShoot() < theCSCHitIneff)
0433 continue;
0434
0435 double eloss = 0;
0436 double pz = fabs(lmom.z());
0437 LocalPoint entry = lpos - 0.5 * thickness * lmom / pz;
0438 LocalPoint exit = lpos + 0.5 * thickness * lmom / pz;
0439 double dtof = path.second * rbeta;
0440 int trkid = mySimTrack.trackId();
0441 unsigned int id = lid.rawId();
0442 short unsigned int processType = 2;
0443 PSimHit hit(
0444 entry, exit, lmom.mag(), tof + dtof, eloss, pid, id, trkid, lmom.theta(), lmom.phi(), processType);
0445 theCSCHits.push_back(hit);
0446 }
0447 } else if (gd->subDetector() == GeomDetEnumerators::RPCBarrel ||
0448 gd->subDetector() == GeomDetEnumerators::RPCEndcap) {
0449 RPCDetId id(gd->geographicalId());
0450 const RPCChamber* chamber = rpcGeom->chamber(id);
0451 std::vector<const RPCRoll*> roll = chamber->rolls();
0452 for (unsigned int iroll = 0; iroll < roll.size(); iroll++) {
0453 RPCDetId rid = roll[iroll]->id();
0454
0455 #ifdef FAMOS_DEBUG
0456 std::cout << " Extrapolated to RPC (" << rid.ring() << "," << rid.station() << "," << rid.sector() << ","
0457 << rid.subsector() << "," << rid.layer() << "," << rid.roll() << ")" << std::endl;
0458 #endif
0459
0460 const GeomDetUnit* det = rpcGeom->idToDetUnit(rid);
0461 HelixArbitraryPlaneCrossing crossing(propagatedState.globalPosition().basicVector(),
0462 propagatedState.globalMomentum().basicVector(),
0463 propagatedState.transverseCurvature(),
0464 anyDirection);
0465 std::pair<bool, double> path = crossing.pathLength(det->surface());
0466 if (!path.first)
0467 continue;
0468 LocalPoint lpos = det->toLocal(GlobalPoint(crossing.position(path.second)));
0469 if (!det->surface().bounds().inside(lpos))
0470 continue;
0471 double thickness = det->surface().bounds().thickness();
0472 LocalVector lmom = det->toLocal(GlobalVector(crossing.direction(path.second)));
0473 lmom = lmom.unit() * propagatedState.localMomentum().mag();
0474 double eloss = 0;
0475 double pz = fabs(lmom.z());
0476 LocalPoint entry = lpos - 0.5 * thickness * lmom / pz;
0477 LocalPoint exit = lpos + 0.5 * thickness * lmom / pz;
0478 double dtof = path.second * rbeta;
0479 int trkid = mySimTrack.trackId();
0480 unsigned int id = rid.rawId();
0481 short unsigned int processType = 2;
0482 PSimHit hit(
0483 entry, exit, lmom.mag(), tof + dtof, eloss, pid, id, trkid, lmom.theta(), lmom.phi(), processType);
0484 theRPCHits.push_back(hit);
0485 }
0486 } else {
0487 std::cout << "Extrapolated to unknown subdetector '" << gd->subDetector() << "'..." << std::endl;
0488 }
0489 }
0490 }
0491 }
0492
0493 std::unique_ptr<edm::PSimHitContainer> pcsc(new edm::PSimHitContainer);
0494 for (std::vector<PSimHit>::const_iterator i = theCSCHits.begin(); i != theCSCHits.end(); i++) {
0495 pcsc->push_back(*i);
0496 }
0497 iEvent.put(std::move(pcsc), "MuonCSCHits");
0498
0499 std::unique_ptr<edm::PSimHitContainer> pdt(new edm::PSimHitContainer);
0500 for (std::vector<PSimHit>::const_iterator i = theDTHits.begin(); i != theDTHits.end(); i++) {
0501 pdt->push_back(*i);
0502 }
0503 iEvent.put(std::move(pdt), "MuonDTHits");
0504
0505 std::unique_ptr<edm::PSimHitContainer> prpc(new edm::PSimHitContainer);
0506 for (std::vector<PSimHit>::const_iterator i = theRPCHits.begin(); i != theRPCHits.end(); i++) {
0507 prpc->push_back(*i);
0508 }
0509 iEvent.put(std::move(prpc), "MuonRPCHits");
0510 }
0511
0512 void MuonSimHitProducer::readParameters(const edm::ParameterSet& fastMuons,
0513 const edm::ParameterSet& fastTracks,
0514 const edm::ParameterSet& matEff) {
0515
0516 std::string _simModuleLabel = fastMuons.getParameter<std::string>("simModuleLabel");
0517 std::string _simModuleProcess = fastMuons.getParameter<std::string>("simModuleProcess");
0518 simMuonLabel = edm::InputTag(_simModuleLabel, _simModuleProcess);
0519 simVertexLabel = edm::InputTag(_simModuleLabel);
0520
0521 std::vector<double> simHitIneffDT = fastMuons.getParameter<std::vector<double> >("simHitDTIneffParameters");
0522 std::vector<double> simHitIneffCSC = fastMuons.getParameter<std::vector<double> >("simHitCSCIneffParameters");
0523 kDT = simHitIneffDT[0];
0524 fDT = simHitIneffDT[1];
0525 kCSC = simHitIneffCSC[0];
0526 fCSC = simHitIneffCSC[1];
0527
0528
0529 fullPattern_ = fastTracks.getUntrackedParameter<bool>("FullPatternRecognition");
0530
0531
0532
0533
0534
0535
0536
0537
0538
0539
0540 theMaterialEffects = nullptr;
0541 if (matEff.getParameter<bool>("PairProduction") || matEff.getParameter<bool>("Bremsstrahlung") ||
0542 matEff.getParameter<bool>("MuonBremsstrahlung") || matEff.getParameter<bool>("EnergyLoss") ||
0543 matEff.getParameter<bool>("MultipleScattering"))
0544 theMaterialEffects = std::make_unique<MaterialEffects>(matEff);
0545 }
0546
0547 void MuonSimHitProducer::applyMaterialEffects(TrajectoryStateOnSurface& tsosWithdEdx,
0548 TrajectoryStateOnSurface& tsos,
0549 double radPath,
0550 RandomEngineAndDistribution const* random,
0551 HepPDT::ParticleDataTable const& table) {
0552
0553 EnergyLossSimulator* energyLoss = theMaterialEffects->energyLossSimulator();
0554
0555
0556 MultipleScatteringSimulator* multipleScattering = theMaterialEffects->multipleScatteringSimulator();
0557
0558
0559 MuonBremsstrahlungSimulator* bremsstrahlung = theMaterialEffects->muonBremsstrahlungSimulator();
0560
0561
0562 const Surface& nextSurface = tsos.surface();
0563 GlobalPoint gPos = energyLoss ? tsos.globalPosition() : tsosWithdEdx.globalPosition();
0564 GlobalVector gMom = energyLoss ? tsos.globalMomentum() : tsosWithdEdx.globalMomentum();
0565 double mu = 0.1056583692;
0566 double en = std::sqrt(gMom.mag2() + mu * mu);
0567
0568
0569 XYZTLorentzVector position(gPos.x(), gPos.y(), gPos.z(), 0.);
0570 XYZTLorentzVector momentum(gMom.x(), gMom.y(), gMom.z(), en);
0571 float charge = (float)(tsos.charge());
0572 ParticlePropagator theMuon(rawparticle::makeMuon(charge < 1., momentum, position), nullptr, nullptr, &table);
0573
0574
0575 if (energyLoss) {
0576
0577
0578 GlobalPoint gPosWithdEdx = tsosWithdEdx.globalPosition();
0579 GlobalVector gMomWithdEdx = tsosWithdEdx.globalMomentum();
0580 double enWithdEdx = std::sqrt(gMomWithdEdx.mag2() + mu * mu);
0581 XYZTLorentzVector deltaPos(
0582 gPosWithdEdx.x() - gPos.x(), gPosWithdEdx.y() - gPos.y(), gPosWithdEdx.z() - gPos.z(), 0.);
0583 XYZTLorentzVector deltaMom(
0584 gMomWithdEdx.x() - gMom.x(), gMomWithdEdx.y() - gMom.y(), gMomWithdEdx.z() - gMom.z(), enWithdEdx - en);
0585
0586
0587 energyLoss->updateState(theMuon, radPath, random);
0588
0589
0590
0591 radPath *= -deltaMom.E() / energyLoss->mostLikelyLoss();
0592 double fac = energyLoss->deltaMom().E() / energyLoss->mostLikelyLoss();
0593
0594
0595 XYZTLorentzVector theNewMomentum = theMuon.particle().momentum() + energyLoss->deltaMom() + fac * deltaMom;
0596 XYZTLorentzVector theNewPosition = theMuon.particle().vertex() + fac * deltaPos;
0597 fac = (theNewMomentum.E() * theNewMomentum.E() - mu * mu) / theNewMomentum.Vect().Mag2();
0598 fac = fac > 0. ? std::sqrt(fac) : 1E-9;
0599 theMuon.particle().setMomentum(
0600 theNewMomentum.Px() * fac, theNewMomentum.Py() * fac, theNewMomentum.Pz() * fac, theNewMomentum.E());
0601 theMuon.particle().setVertex(theNewPosition);
0602 }
0603
0604
0605 if (multipleScattering) {
0606
0607 GlobalVector normal = nextSurface.tangentPlane(tsos.globalPosition())->normalVector();
0608 multipleScattering->setNormalVector(normal);
0609
0610 multipleScattering->updateState(theMuon, radPath, random);
0611 }
0612
0613
0614 if (bremsstrahlung) {
0615
0616 bremsstrahlung->updateState(theMuon, radPath, random);
0617 }
0618
0619
0620 GlobalPoint propagatedPosition(theMuon.particle().X(), theMuon.particle().Y(), theMuon.particle().Z());
0621 GlobalVector propagatedMomentum(theMuon.particle().Px(), theMuon.particle().Py(), theMuon.particle().Pz());
0622 GlobalTrajectoryParameters propagatedGtp(propagatedPosition, propagatedMomentum, (int)charge, magfield);
0623 tsosWithdEdx = TrajectoryStateOnSurface(propagatedGtp, nextSurface);
0624 }
0625
0626
0627 #include "FWCore/Framework/interface/MakerMacros.h"
0628 DEFINE_FWK_MODULE(MuonSimHitProducer);