VertexFinder/src/InputData.cc

0001 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0002 #include "FWCore/Framework/interface/Event.h"
0003 #include "FWCore/Framework/interface/EventSetup.h"
0004 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0005 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0006 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0007 #include "SimDataFormats/EncodedEventId/interface/EncodedEventId.h"
0008 #include "SimDataFormats/Track/interface/SimTrack.h"
0009 #include "SimDataFormats/Associations/interface/TTClusterAssociationMap.h"
0010 #include "SimDataFormats/Associations/interface/TTStubAssociationMap.h"
0011 #include "L1Trigger/VertexFinder/interface/AnalysisSettings.h"
0012 #include "L1Trigger/VertexFinder/interface/InputData.h"
0013
0014 #include <map>
0015
0016 using namespace std;
0017
0018 namespace l1tVertexFinder {
0019
0020   InputData::InputData() {}
0021
0022   InputData::InputData(const edm::Event& iEvent,
0023                        const edm::EventSetup& iSetup,
0024                        const AnalysisSettings& settings,
0025                        const edm::EDGetTokenT<edm::HepMCProduct> hepMCToken,
0026                        const edm::EDGetTokenT<edm::View<reco::GenParticle>> genParticlesToken,
0027                        const edm::EDGetTokenT<edm::View<TrackingParticle>> tpToken,
0028                        const edm::EDGetTokenT<edm::ValueMap<l1tVertexFinder::TP>> tpValueMapToken,
0029                        const edm::EDGetTokenT<DetSetVec> stubToken,
0030                        edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> tTopoToken,
0031                        edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tGeomToken) {
0032     // Get the tracker geometry info needed to unpack the stub info.
0033     const TrackerTopology& tTopo = iSetup.getData(tTopoToken);
0034     const TrackerGeometry& tGeom = iSetup.getData(tGeomToken);
0035
0036     // Get stub info, by looping over modules and then stubs inside each module.
0037     // Also get the association map from stubs to tracking particles.
0038     edm::Handle<DetSetVec> ttStubHandle;
0039     iEvent.getByToken(stubToken, ttStubHandle);
0040
0041     std::set<DetId> lStubDetIds;
0042     for (DetSetVec::const_iterator p_module = ttStubHandle->begin(); p_module != ttStubHandle->end(); p_module++) {
0043       for (DetSet::const_iterator p_ttstub = p_module->begin(); p_ttstub != p_module->end(); p_ttstub++) {
0044         lStubDetIds.insert(p_ttstub->getDetId());
0045       }
0046     }
0047
0048     std::map<DetId, DetId> stubGeoDetIdMap;
0049     for (auto gd = tGeom.dets().begin(); gd != tGeom.dets().end(); gd++) {
0050       DetId detid = (*gd)->geographicalId();
0051       if (detid.subdetId() != StripSubdetector::TOB && detid.subdetId() != StripSubdetector::TID)
0052         continue;  // only run on OT
0053       if (!tTopo.isLower(detid))
0054         continue;                             // loop on the stacks: choose the lower arbitrarily
0055       DetId stackDetid = tTopo.stack(detid);  // Stub module detid
0056
0057       if (lStubDetIds.count(stackDetid) > 0) {
0058         assert(stubGeoDetIdMap.count(stackDetid) == 0);
0059         stubGeoDetIdMap[stackDetid] = detid;
0060       }
0061     }
0062     assert(lStubDetIds.size() == stubGeoDetIdMap.size());
0063
0064     // Get TrackingParticle info
0065     edm::Handle<edm::View<TrackingParticle>> tpHandle;
0066     edm::Handle<edm::ValueMap<TP>> tpValueMapHandle;
0067     iEvent.getByToken(tpToken, tpHandle);
0068     iEvent.getByToken(tpValueMapToken, tpValueMapHandle);
0069     edm::ValueMap<TP> tpValueMap = *tpValueMapHandle;
0070
0071     for (unsigned int i = 0; i < tpHandle->size(); i++) {
0072       if (tpValueMap[tpHandle->refAt(i)].use()) {
0073         tpPtrToRefMap_[tpHandle->ptrAt(i)] = tpHandle->refAt(i);
0074       }
0075     }
0076
0077     // Find the various vertices
0078     genPt_ = 0.;
0079     genPt_PU_ = 0.;
0080     for (const auto& [edmPtr, edmRef] : tpPtrToRefMap_) {
0081       TP tp = tpValueMap[edmRef];
0082       if (tp.physicsCollision()) {
0083         genPt_ += tp->pt();
0084       } else {
0085         genPt_PU_ += tp->pt();
0086       }
0087       if (settings.debug() > 2) {
0088         edm::LogInfo("InputData") << "InputData::genPt in the event " << genPt_;
0089       }
0090
0091       if (tp.useForAlgEff()) {
0092         vertex_.insert(tp);
0093       } else if (tp.useForVertexReco()) {
0094         bool found = false;
0095         for (unsigned int i = 0; i < vertices_.size(); ++i) {
0096           if (tp->vz() == vertices_[i].vz()) {
0097             vertices_[i].insert(tp);
0098             found = true;
0099             break;
0100           }
0101         }
0102         if (!found) {
0103           Vertex vertex(tp->vz());
0104           vertex.insert(tp);
0105           vertices_.push_back(vertex);
0106         }
0107       }
0108     }
0109
0110     for (const Vertex& vertex : vertices_) {
0111       if (vertex.numTracks() >= settings.vx_minTracks())
0112         recoVertices_.push_back(vertex);
0113     }
0114
0115     if (settings.debug() > 2)
0116       edm::LogInfo("InputData") << "InputData::" << vertices_.size() << " pileup vertices in the event, "
0117                                 << recoVertices_.size() << " reconstructable";
0118
0119     vertex_.computeParameters();
0120     if (settings.debug() > 2)
0121       edm::LogInfo("InputData") << "InputData::Vertex " << vertex_.z0() << " containing " << vertex_.numTracks()
0122                                 << " total pT " << vertex_.pT();
0123
0124     for (unsigned int i = 0; i < vertices_.size(); ++i) {
0125       vertices_[i].computeParameters();
0126     }
0127
0128     for (unsigned int i = 0; i < recoVertices_.size(); ++i) {
0129       recoVertices_[i].computeParameters();
0130     }
0131
0132     std::sort(vertices_.begin(), vertices_.end(), SortVertexByPt());
0133     std::sort(recoVertices_.begin(), recoVertices_.end(), SortVertexByPt());
0134
0135     // Form the HepMC and GenParticle based vertices
0136     edm::Handle<edm::HepMCProduct> HepMCEvt;
0137     iEvent.getByToken(hepMCToken, HepMCEvt);
0138
0139     edm::Handle<edm::View<reco::GenParticle>> GenParticleHandle;
0140     iEvent.getByToken(genParticlesToken, GenParticleHandle);
0141
0142     if (!HepMCEvt.isValid() && !GenParticleHandle.isValid()) {
0143       throw cms::Exception("Neither the edm::HepMCProduct nor the generator particles are available.");
0144     }
0145     if (HepMCEvt.isValid()) {
0146       const HepMC::GenEvent* MCEvt = HepMCEvt->GetEvent();
0147       for (HepMC::GenEvent::vertex_const_iterator ivertex = MCEvt->vertices_begin(); ivertex != MCEvt->vertices_end();
0148            ++ivertex) {
0149         bool hasParentVertex = false;
0150
0151         // Loop over the parents looking to see if they are coming from a production vertex
0152         for (HepMC::GenVertex::particle_iterator iparent = (*ivertex)->particles_begin(HepMC::parents);
0153              iparent != (*ivertex)->particles_end(HepMC::parents);
0154              ++iparent)
0155           if ((*iparent)->production_vertex()) {
0156             hasParentVertex = true;
0157             break;
0158           }
0159
0160         // Reject those vertices with parent vertices
0161         if (hasParentVertex)
0162           continue;
0163         // Get the position of the vertex
0164         HepMC::FourVector pos = (*ivertex)->position();
0165         const double mm = 0.1;  // [mm] --> [cm]
0166         hepMCVertex_ = Vertex(pos.z() * mm);
0167         break;  // there should be one single primary vertex
0168       }  // end loop over gen vertices
0169     }
0170     if (GenParticleHandle.isValid()) {
0171       for (const auto& genpart : *GenParticleHandle) {
0172         if ((genpart.status() != 1) || (genpart.numberOfMothers() == 0))  // not stable or one of the incoming hadrons
0173           continue;
0174         genVertex_ = Vertex(genpart.vz());
0175         break;
0176       }
0177     }
0178     if ((hepMCVertex_.vz() == 0.0) && (genVertex_.vz() == 0.0)) {
0179       throw cms::Exception("Neither the HepMC vertex nor the generator particle vertex were found.");
0180     }
0181
0182   }  // end InputData::InputData
0183
0184   InputData::~InputData() {}
0185
0186 }  // end namespace l1tVertexFinder