Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:00:50 

0001 //
0002 // Package:    MuonSimHitProducer
0003 // Class:      MuonSimHitProducer
0004 //
0005 /**\class MuonSimHitProducer FastSimulation/MuonSimHitProducer/src/MuonSimHitProducer.cc
0006 
0007  Description:
0008     Fast simulation producer of Muon Sim Hits (to be used for realistic Muon reconstruction)
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 
0013 */
0014 //
0015 // Original Author:  Martijn Mulders/Matthew Jones
0016 //         Created:  Wed Jul 30 11:37:24 CET 2007
0017 //         Working:  Fri Nov  9 09:39:33 CST 2007
0018 //
0019 // $Id: MuonSimHitProducer.cc,v 1.36 2011/10/07 08:25:42 aperrott Exp $
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   // Parameters to emulate the muonSimHit association inefficiency due to delta's
0073   double kDT;
0074   double fDT;
0075   double kCSC;
0076   double fCSC;
0077 
0078   /// Simulate material effects in iron (dE/dx, multiple scattering)
0079   void applyMaterialEffects(TrajectoryStateOnSurface& tsosWithdEdx,
0080                             TrajectoryStateOnSurface& tsos,
0081                             double radPath,
0082                             RandomEngineAndDistribution const*,
0083                             HepPDT::ParticleDataTable const&);
0084 
0085   // ----------- parameters ----------------------------
0086   bool fullPattern_;
0087   bool doL1_, doL3_, doGL_;
0088 
0089   // tags
0090   edm::InputTag simMuonLabel;
0091   edm::InputTag simVertexLabel;
0092 
0093   // tokens
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 //for debug only
0105 //#define FAMOS_DEBUG
0106 
0107 //
0108 // constructors and destructor
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   // Read relevant parameters
0119   readParameters(iConfig.getParameter<edm::ParameterSet>("MUONS"),
0120                  iConfig.getParameter<edm::ParameterSet>("TRACKS"),
0121                  iConfig.getParameter<edm::ParameterSet>("MaterialEffectsForMuons"));
0122 
0123   //
0124   //  register your products ... need to declare at least one possible product...
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   // consumes
0134   simMuonToken = consumes<std::vector<SimTrack> >(simMuonLabel);
0135   simVertexToken = consumes<std::vector<SimVertex> >(simVertexLabel);
0136 }
0137 
0138 // ---- method called once each job just before starting event loop ----
0139 void MuonSimHitProducer::beginRun(edm::Run const& run, const edm::EventSetup& es) {
0140   //services
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   // A few propagators
0150   propagatorWithMaterial = &(*(theService->propagator("SteppingHelixPropagatorAny")));
0151   propagatorWithoutMaterial.reset(propagatorWithMaterial->clone());
0152   SteppingHelixPropagator* SHpropagator =
0153       dynamic_cast<SteppingHelixPropagator*>(propagatorWithoutMaterial.get());  // Beuark!
0154   SHpropagator->setMaterialMode(true);                                          // switches OFF material effects;
0155 }
0156 
0157 //
0158 // member functions
0159 //
0160 
0161 // ------------ method called to produce the data  ------------
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     // Decaying hadrons are now in the list, and so are their muon daughter
0186     // Ignore the hadrons here.
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     //  Presumably t0 has dimensions of cm if not mm?
0202     //  Convert to ns for internal calculations.
0203     //  I wonder where we should get c from?
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     //  Produce muons sim hits starting from undecayed simulated muons
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     //  Crap... there's no time-of-flight to the trackerSurfacePosition()...
0224     //  So, this will be wrong when the curvature can't be neglected, but that
0225     //  will be rather seldom...  May as well ignore the mass too.
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     //  Some magic to define a TrajectoryStateOnSurface
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     edm::ESHandle<Propagator> propagator =
0255       theService->propagator("SteppingHelixPropagatorAny");
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       // Starting momentum
0281       double pi = propagatedState.globalMomentum().mag();
0282 
0283       // Propagate with material effects (dE/dx average only)
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       // No need to continue if there is no valid propagation available.
0291       // This happens rarely (~0.1% of ttbar events)
0292       if (!next.first.isValid())
0293         continue;
0294       // This is the estimate of the number of radiation lengths traversed,
0295       // together with the total path length
0296       double radPath = shsDest.radPath();
0297       double pathLength = next.second;
0298 
0299       // Now propagate without dE/dx (average)
0300       // [To add the dE/dx fluctuations to the actual dE/dx]
0301       std::pair<TrajectoryStateOnSurface, double> nextNoMaterial =
0302           propagatorWithoutMaterial->propagateWithPath(propagatedState, navLayers[ilayer]->surface());
0303 
0304       // Update the propagated state
0305       propagatedState = next.first;
0306       double pf = propagatedState.globalMomentum().mag();
0307 
0308       // Insert dE/dx fluctuations and multiple scattering
0309       // Skip this step if nextNoMaterial.first is not valid
0310       // This happens rarely (~0.02% of ttbar events)
0311       if (theMaterialEffects && nextNoMaterial.first.isValid())
0312         applyMaterialEffects(propagatedState, nextNoMaterial.first, radPath, &random, *pdg);
0313       // Check that the 'shaken' propagatedState is still valid, otherwise continue
0314       if (!propagatedState.isValid())
0315         continue;
0316       // (No evidence that this ever happens)
0317       //
0318       //  Consider this... 1 GeV muon has a velocity that is only 0.5% slower than c...
0319       //  We probably can safely ignore the mass for anything that makes it out to the
0320       //  muon chambers.
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               // no drift cell here (on the chamber edge or just outside)
0366               // this hit would otherwise be discarded downstream in the digitizer
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               // Factor that takes into account the (rec)hits lost because of delta's, etc.:
0374               // (Not fully satisfactory patch, but it seems to work...)
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             // For CSCs the Bounds are for chamber frames not gas regions
0415             //      if ( ! det->surface().bounds().inside(lpos) ) continue;
0416             // New function knows where the 'active' volume is:
0417             const CSCLayerGeometry* laygeom = layers[ilayer]->geometry();
0418             if (!laygeom->inside(lpos))
0419               continue;
0420             //double thickness = laygeom->thickness(); gives number which is about 20 times too big
0421             double thickness = det->surface().bounds().thickness();  // this one works much better...
0422             LocalVector lmom = det->toLocal(GlobalVector(crossing.direction(path.second)));
0423             lmom = lmom.unit() * propagatedState.localMomentum().mag();
0424 
0425             // Factor that takes into account the (rec)hits lost because of delta's, etc.:
0426             // (Not fully satisfactory patch, but it seems to work...)
0427             double pmu = lmom.mag();
0428             double theCSCHitIneff = pmu > 0 ? exp(kCSC * log(pmu) + fCSC) : 0.;
0429             // Take into account the different geometry in ME11:
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   // Muons
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   // Tracks
0529   fullPattern_ = fastTracks.getUntrackedParameter<bool>("FullPatternRecognition");
0530 
0531   // The following should be on LogInfo
0532   //  std::cout << " MUON SIM HITS: FastSimulation parameters " << std::endl;
0533   //  std::cout << " ============================================== " << std::endl;
0534   //  if ( fullPattern_ )
0535   //    std::cout << " The FULL pattern recognition option is turned ON" << std::endl;
0536   //  else
0537   //    std::cout << " The FAST tracking option is turned ON" << std::endl;
0538 
0539   // Material Effects
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   // The energy loss simulator
0553   EnergyLossSimulator* energyLoss = theMaterialEffects->energyLossSimulator();
0554 
0555   // The multiple scattering simulator
0556   MultipleScatteringSimulator* multipleScattering = theMaterialEffects->multipleScatteringSimulator();
0557 
0558   // The Muon Bremsstrahlung simulator
0559   MuonBremsstrahlungSimulator* bremsstrahlung = theMaterialEffects->muonBremsstrahlungSimulator();
0560 
0561   // Initialize the Particle position, momentum and energy
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   // And now create the Particle
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   // Recompute the energy loss to get the fluctuations
0575   if (energyLoss) {
0576     // Difference between with and without dE/dx (average only)
0577     // (for corrections once fluctuations are applied)
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     // Energy loss in iron (+ fluctuations)
0587     energyLoss->updateState(theMuon, radPath, random);
0588 
0589     // Correcting factors to account for slight differences in material descriptions
0590     // (Material description is more accurate in the stepping helix propagator)
0591     radPath *= -deltaMom.E() / energyLoss->mostLikelyLoss();
0592     double fac = energyLoss->deltaMom().E() / energyLoss->mostLikelyLoss();
0593 
0594     // Particle momentum & position after energy loss + fluctuation
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   // Does the actual mutliple scattering
0605   if (multipleScattering) {
0606     // Pass the vector normal to the "next" surface
0607     GlobalVector normal = nextSurface.tangentPlane(tsos.globalPosition())->normalVector();
0608     multipleScattering->setNormalVector(normal);
0609     // Compute the amount of multiple scattering after a given path length
0610     multipleScattering->updateState(theMuon, radPath, random);
0611   }
0612 
0613   // Muon Bremsstrahlung
0614   if (bremsstrahlung) {
0615     // Compute the amount of Muon Bremsstrahlung after given path length
0616     bremsstrahlung->updateState(theMuon, radPath, random);
0617   }
0618 
0619   // Fill the propagated state
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 //define this as a plug-in
0627 #include "FWCore/Framework/interface/MakerMacros.h"
0628 DEFINE_FWK_MODULE(MuonSimHitProducer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!