Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2025-01-09 23:33:59

0001 // Class Header
0002 #include "RecoMuon/MuonSeedGenerator/test/MCSeedGenerator/MCMuonSeedGenerator.h"
0003 
0004 // Data Formats
0005 #include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
0006 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0007 #include "DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h"
0008 #include "DataFormats/Common/interface/Handle.h"
0009 
0010 // Framework
0011 #include "FWCore/Framework/interface/ConsumesCollector.h"
0012 #include "FWCore/Framework/interface/EventSetup.h"
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 
0017 #include "FWCore/Framework/interface/MakerMacros.h"
0018 
0019 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0020 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0021 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0022 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0023 
0024 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0025 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0026 
0027 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0028 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0029 
0030 DEFINE_FWK_MODULE(MCMuonSeedGenerator2);
0031 
0032 using namespace std;
0033 using namespace edm;
0034 
0035 // constructors
0036 MCMuonSeedGenerator2::MCMuonSeedGenerator2(const edm::ParameterSet& parameterSet) : theSeedType(FromTracks) {
0037   theCSCSimHitLabel = parameterSet.getParameter<InputTag>("CSCSimHit");
0038   theDTSimHitLabel = parameterSet.getParameter<InputTag>("DTSimHit");
0039   theRPCSimHitLabel = parameterSet.getParameter<InputTag>("RPCSimHit");
0040   theSimTrackLabel = parameterSet.getParameter<InputTag>("SimTrack");
0041   theSimVertexLabel = parameterSet.getParameter<InputTag>("SimVertex");
0042 
0043   // service parameters
0044   ParameterSet serviceParameters = parameterSet.getParameter<ParameterSet>("ServiceParameters");
0045 
0046   // the services
0047   theService = new MuonServiceProxy(serviceParameters, consumesCollector());
0048 
0049   // Error Scale
0050   theErrorScale = parameterSet.getParameter<double>("ErrorScale");
0051 
0052   string seedType = parameterSet.getParameter<string>("SeedType");
0053 
0054   if (seedType == "FromHits")
0055     theSeedType = FromHits;
0056   else if (seedType == "FromTracks")
0057     theSeedType = FromTracks;
0058   else
0059     LogError("Muon|RecoMuon|MCMuonSeedGenerator2") << "Wrong seed type!!! " << seedType;
0060 
0061   produces<TrajectorySeedCollection>();
0062 }
0063 
0064 // destructor
0065 MCMuonSeedGenerator2::~MCMuonSeedGenerator2() {
0066   if (theService)
0067     delete theService;
0068 }
0069 
0070 void MCMuonSeedGenerator2::produce(edm::Event& event, const edm::EventSetup& setup) {
0071   const std::string metname = "Muon|RecoMuon|MCMuonSeedGenerator2";
0072 
0073   auto output = std::make_unique<TrajectorySeedCollection>();
0074 
0075   // Update the services
0076   theService->update(setup);
0077 
0078   // Get the SimHit collection from the event
0079   Handle<PSimHitContainer> dtSimHits;
0080   event.getByLabel(theDTSimHitLabel.instance(), theDTSimHitLabel.label(), dtSimHits);
0081 
0082   Handle<PSimHitContainer> cscSimHits;
0083   event.getByLabel(theCSCSimHitLabel.instance(), theCSCSimHitLabel.label(), cscSimHits);
0084 
0085   Handle<PSimHitContainer> rpcSimHits;
0086   event.getByLabel(theRPCSimHitLabel.instance(), theRPCSimHitLabel.label(), rpcSimHits);
0087 
0088   Handle<SimTrackContainer> simTracks;
0089   event.getByLabel(theSimTrackLabel.label(), simTracks);
0090 
0091   Handle<SimVertexContainer> simVertices;
0092   event.getByLabel<SimVertexContainer>(theSimVertexLabel, simVertices);
0093 
0094   map<unsigned int, vector<const PSimHit*> > mapOfMuonSimHits;
0095 
0096   for (PSimHitContainer::const_iterator simhit = dtSimHits->begin(); simhit != dtSimHits->end(); ++simhit) {
0097     if (abs(simhit->particleType()) != 13)
0098       continue;
0099     mapOfMuonSimHits[simhit->trackId()].push_back(&*simhit);
0100   }
0101 
0102   for (PSimHitContainer::const_iterator simhit = cscSimHits->begin(); simhit != cscSimHits->end(); ++simhit) {
0103     if (abs(simhit->particleType()) != 13)
0104       continue;
0105     mapOfMuonSimHits[simhit->trackId()].push_back(&*simhit);
0106   }
0107 
0108   for (PSimHitContainer::const_iterator simhit = rpcSimHits->begin(); simhit != rpcSimHits->end(); ++simhit) {
0109     if (abs(simhit->particleType()) != 13)
0110       continue;
0111     mapOfMuonSimHits[simhit->trackId()].push_back(&*simhit);
0112   }
0113 
0114   for (SimTrackContainer::const_iterator simTrack = simTracks->begin(); simTrack != simTracks->end(); ++simTrack) {
0115     if (abs(simTrack->type()) != 13)
0116       continue;
0117 
0118     map<unsigned int, vector<const PSimHit*> >::const_iterator mapIterator = mapOfMuonSimHits.find(simTrack->trackId());
0119 
0120     if (mapIterator == mapOfMuonSimHits.end()) {
0121       LogTrace(metname) << "...Very strange, no sim hits associated to the sim track!"
0122                         << "\n"
0123                         << "SimTrack's eta: " << simTrack->momentum().eta();
0124       continue;
0125     }
0126 
0127     vector<const PSimHit*> muonSimHits = mapIterator->second;
0128 
0129     if (muonSimHits.size() < 1)
0130       continue;
0131 
0132     stable_sort(muonSimHits.begin(), muonSimHits.end(), RadiusComparatorInOut(theService->trackingGeometry()));
0133 
0134     const PSimHit* innerSimHit = muonSimHits.front();
0135 
0136     TrajectorySeed* seed = 0;
0137 
0138     // switch(theSeedType){
0139 
0140     //     case FromHits:{
0141     //       seed = createSeedFromHit(innerSimHit);
0142     //       break;}
0143 
0144     //     case FromTracks:{
0145     //       seed = createSeedFromTrack(*simTrack, (*simVertices)[simTrack->vertIndex()], DetId(innerSimHit->detUnitId()));
0146     //       break;}
0147     //     }
0148 
0149     if (theSeedType == FromTracks)
0150       seed = createSeedFromTrack(*simTrack, (*simVertices)[simTrack->vertIndex()], DetId(innerSimHit->detUnitId()));
0151     else if (theSeedType == FromHits)
0152       seed = createSeedFromHit(innerSimHit);
0153     else {
0154       LogError(metname) << "ERROR!!";
0155       continue;
0156     }
0157 
0158     if (seed)
0159       output->push_back(*seed);
0160   }
0161 
0162   event.put(std::move(output));
0163 }
0164 
0165 TrajectorySeed* MCMuonSeedGenerator2::createSeedFromHit(const PSimHit* innerSimHit) {
0166   const std::string metname = "Muon|RecoMuon|MCMuonSeedGenerator2";
0167   MuonPatternRecoDumper debug;
0168 
0169   LogTrace(metname) << "Seed from sim hit";
0170 
0171   const GeomDet* geomDet = theService->trackingGeometry()->idToDet(DetId(innerSimHit->detUnitId()));
0172   LogTrace(metname) << "Seed geom det: " << debug.dumpMuonId(geomDet->geographicalId());
0173 
0174   GlobalPoint glbPosition = geomDet->toGlobal(innerSimHit->localPosition());
0175   GlobalVector glbMomentum = geomDet->toGlobal(innerSimHit->momentumAtEntry());
0176 
0177   const GlobalTrajectoryParameters globalParameters(glbPosition,
0178                                                     glbMomentum,
0179                                                     -innerSimHit->particleType() / abs(innerSimHit->particleType()),
0180                                                     &*theService->magneticField());
0181 
0182   AlgebraicSymMatrix66 covarianceMatrix = AlgebraicMatrixID();
0183   covarianceMatrix *= theErrorScale;
0184 
0185   const CartesianTrajectoryError cartesianErrors(covarianceMatrix);
0186 
0187   TrajectoryStateOnSurface tsos(globalParameters, cartesianErrors, geomDet->surface());
0188 
0189   LogTrace(metname) << "State on: " << debug.dumpMuonId(DetId(innerSimHit->detUnitId()));
0190   LogTrace(metname) << debug.dumpTSOS(tsos);
0191 
0192   // convert the TSOS into a PTSOD
0193   PTrajectoryStateOnDet seedTSOS = trajectoryStateTransform::persistentState(tsos, geomDet->geographicalId().rawId());
0194 
0195   edm::OwnVector<TrackingRecHit> container;
0196   TrajectorySeed* seed = new TrajectorySeed(seedTSOS, container, alongMomentum);
0197 
0198   return seed;
0199 }
0200 
0201 TrajectorySeed* MCMuonSeedGenerator2::createSeedFromTrack(const SimTrack& simTrack,
0202                                                           const SimVertex& simVertex,
0203                                                           DetId detId) {
0204   const std::string metname = "Muon|RecoMuon|MCMuonSeedGenerator2";
0205   MuonPatternRecoDumper debug;
0206 
0207   LogTrace(metname) << "Seed from sim track";
0208 
0209   TrajectorySeed* seed = 0;
0210 
0211   const GeomDet* geomDet = theService->trackingGeometry()->idToDet(detId);
0212   LogTrace(metname) << "Seed geom det: " << debug.dumpMuonId(geomDet->geographicalId());
0213 
0214   GlobalVector simMomentum(simTrack.momentum().x(), simTrack.momentum().y(), simTrack.momentum().z());
0215 
0216   GlobalPoint simPosition(0., 0., 0.);
0217   if (simTrack.vertIndex() >= 0)
0218     simPosition = GlobalPoint(simVertex.position().x(), simVertex.position().y(), simVertex.position().z());
0219 
0220   AlgebraicSymMatrix66 matI = AlgebraicMatrixID();
0221   matI *= 100;  // 1e-20;
0222   CartesianTrajectoryError simCov(matI);
0223 
0224   GlobalTrajectoryParameters parameters(
0225       simPosition, simMomentum, int(simTrack.charge()), &*theService->magneticField());
0226   FreeTrajectoryState simFTS(parameters, simCov);
0227 
0228   LogTrace(metname) << "FTS from the Seed";
0229   LogTrace(metname) << debug.dumpFTS(simFTS);
0230 
0231   TrajectoryStateOnSurface simSeedTSOS =
0232       theService->propagator("SteppingHelixPropagatorAlong")->propagate(simFTS, geomDet->surface());
0233 
0234   if (!simSeedTSOS.isValid()) {
0235     LogTrace(metname) << "Propagation from IP to the hit layer failed";
0236     return seed;
0237   }
0238 
0239   // convert the TSOS into a PTSOD
0240   PTrajectoryStateOnDet const& seedTSOS =
0241       trajectoryStateTransform::persistentState(simSeedTSOS, geomDet->geographicalId().rawId());
0242 
0243   LogTrace(metname) << "State on: " << debug.dumpMuonId(detId);
0244   LogTrace(metname) << debug.dumpTSOS(simSeedTSOS);
0245 
0246   edm::OwnVector<TrackingRecHit> container;
0247   seed = new TrajectorySeed(seedTSOS, container, alongMomentum);
0248 
0249   return seed;
0250 }