File indexing completed on 2025-01-09 23:33:59
0001
0002 #include "RecoMuon/MuonSeedGenerator/test/MCSeedGenerator/MCMuonSeedGenerator.h"
0003
0004
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
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
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
0044 ParameterSet serviceParameters = parameterSet.getParameter<ParameterSet>("ServiceParameters");
0045
0046
0047 theService = new MuonServiceProxy(serviceParameters, consumesCollector());
0048
0049
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
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
0076 theService->update(setup);
0077
0078
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
0139
0140
0141
0142
0143
0144
0145
0146
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
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;
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
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 }