test/MCSeedGenerator/MCMuonSeedGenerator.cc

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