File indexing completed on 2024-04-06 12:31:07
0001
0002
0003
0004
0005 #include <cmath>
0006 #include <cstdlib>
0007 #include <iostream>
0008
0009 #include "HepPDT/ParticleID.hh"
0010
0011 #include "SimTracker/TrackHistory/interface/VertexClassifier.h"
0012
0013 #define update(a, b) \
0014 do { \
0015 (a) = (a) || (b); \
0016 } while (0)
0017
0018 VertexClassifier::VertexClassifier(edm::ParameterSet const &config, edm::ConsumesCollector collector)
0019 : VertexCategories(),
0020 tracer_(config, collector),
0021 hepMCLabel_(config.getUntrackedParameter<edm::InputTag>("hepMC")),
0022 particleDataTableToken_(collector.esConsumes()) {
0023 collector.consumes<edm::HepMCProduct>(hepMCLabel_);
0024
0025 tracer_.depth(-2);
0026
0027
0028 longLivedDecayLength_ = config.getUntrackedParameter<double>("longLivedDecayLength");
0029
0030
0031 vertexClusteringDistance_ = config.getUntrackedParameter<double>("vertexClusteringDistance");
0032 }
0033
0034 void VertexClassifier::newEvent(edm::Event const &event, edm::EventSetup const &setup) {
0035
0036 tracer_.newEvent(event, setup);
0037
0038
0039 event.getByLabel(hepMCLabel_, mcInformation_);
0040
0041
0042 particleDataTable_ = setup.getHandle(particleDataTableToken_);
0043
0044
0045 genPrimaryVertices();
0046 }
0047
0048 VertexClassifier const &VertexClassifier::evaluate(reco::VertexBaseRef const &vertex) {
0049
0050 reset();
0051
0052
0053 if (tracer_.evaluate(vertex)) {
0054
0055 simulationInformation();
0056
0057
0058 processesAtGenerator();
0059
0060
0061 processesAtSimulation();
0062
0063
0064 vertexInformation();
0065
0066
0067 unknownVertex();
0068 } else
0069 flags_[Fake] = true;
0070
0071 return *this;
0072 }
0073
0074 VertexClassifier const &VertexClassifier::evaluate(TrackingVertexRef const &vertex) {
0075
0076 reset();
0077
0078
0079 tracer_.evaluate(vertex);
0080
0081
0082 if (tracer_.recoVertex().isNonnull())
0083 flags_[Reconstructed] = true;
0084 else
0085 flags_[Reconstructed] = false;
0086
0087
0088 simulationInformation();
0089
0090
0091 processesAtGenerator();
0092
0093
0094 processesAtSimulation();
0095
0096
0097 vertexInformation();
0098
0099
0100 unknownVertex();
0101
0102 return *this;
0103 }
0104
0105 void VertexClassifier::simulationInformation() {
0106
0107 EncodedEventId eventId = tracer_.simVertex()->eventId();
0108
0109 flags_[SignalEvent] = !eventId.bunchCrossing() && !eventId.event();
0110 }
0111
0112 void VertexClassifier::processesAtGenerator() {
0113
0114 VertexHistory::GenVertexTrail const &genVertexTrail = tracer_.genVertexTrail();
0115
0116
0117 for (VertexHistory::GenVertexTrail::const_iterator ivertex = genVertexTrail.begin(); ivertex != genVertexTrail.end();
0118 ++ivertex) {
0119
0120
0121 HepMC::GenVertex *vertex = const_cast<HepMC::GenVertex *>(*ivertex);
0122
0123
0124 for (HepMC::GenVertex::particle_iterator iparent = vertex->particles_begin(HepMC::parents);
0125 iparent != vertex->particles_end(HepMC::parents);
0126 ++iparent) {
0127
0128 int pdgid = std::abs((*iparent)->pdg_id());
0129
0130 HepPDT::ParticleID particleID(pdgid);
0131
0132
0133 if (particleID.isValid()) {
0134
0135 ParticleData const *particleData = particleDataTable_->particle(particleID);
0136
0137 if (particleData) {
0138
0139 if (particleData->lifetime() > longLivedDecayLength_) {
0140
0141 update(flags_[BWeakDecay], particleID.hasBottom());
0142 update(flags_[CWeakDecay], particleID.hasCharm());
0143 update(flags_[LongLivedDecay], true);
0144 }
0145
0146 update(flags_[TauDecay], pdgid == 15);
0147 update(flags_[KsDecay], pdgid == 310);
0148 update(flags_[LambdaDecay], pdgid == 3122);
0149 update(flags_[JpsiDecay], pdgid == 443);
0150 update(flags_[XiDecay], pdgid == 3312);
0151 update(flags_[OmegaDecay], pdgid == 3334);
0152 update(flags_[SigmaPlusDecay], pdgid == 3222);
0153 update(flags_[SigmaMinusDecay], pdgid == 3112);
0154 }
0155 }
0156 }
0157 }
0158 }
0159
0160 void VertexClassifier::processesAtSimulation() {
0161 VertexHistory::SimVertexTrail const &simVertexTrail = tracer_.simVertexTrail();
0162
0163 for (VertexHistory::SimVertexTrail::const_iterator ivertex = simVertexTrail.begin(); ivertex != simVertexTrail.end();
0164 ++ivertex) {
0165
0166 int pdgid = 0;
0167
0168
0169 bool flag = false;
0170 TrackingVertex::tp_iterator itd, its;
0171
0172 for (its = (*ivertex)->sourceTracks_begin(); its != (*ivertex)->sourceTracks_end(); ++its) {
0173 for (itd = (*ivertex)->daughterTracks_begin(); itd != (*ivertex)->daughterTracks_end(); ++itd)
0174 if (itd != its) {
0175 flag = true;
0176 break;
0177 }
0178 if (flag)
0179 break;
0180 }
0181
0182 if (its != (*ivertex)->sourceTracks_end())
0183 pdgid = std::abs((*its)->pdgId());
0184 else
0185 pdgid = 0;
0186
0187
0188
0189 unsigned int processG4 = 0;
0190
0191 if ((*ivertex)->nG4Vertices() > 0) {
0192 processG4 = (*(*ivertex)->g4Vertices_begin()).processType();
0193 }
0194
0195 unsigned int process = g4toCMSProcMap_.processId(processG4);
0196
0197
0198 update(flags_[KnownProcess], process != CMS::Undefined && process != CMS::Unknown && process != CMS::Primary);
0199
0200 update(flags_[UndefinedProcess], process == CMS::Undefined);
0201 update(flags_[UnknownProcess], process == CMS::Unknown);
0202 update(flags_[PrimaryProcess], process == CMS::Primary);
0203 update(flags_[HadronicProcess], process == CMS::Hadronic);
0204 update(flags_[DecayProcess], process == CMS::Decay);
0205 update(flags_[ComptonProcess], process == CMS::Compton);
0206 update(flags_[AnnihilationProcess], process == CMS::Annihilation);
0207 update(flags_[EIoniProcess], process == CMS::EIoni);
0208 update(flags_[HIoniProcess], process == CMS::HIoni);
0209 update(flags_[MuIoniProcess], process == CMS::MuIoni);
0210 update(flags_[PhotonProcess], process == CMS::Photon);
0211 update(flags_[MuPairProdProcess], process == CMS::MuPairProd);
0212 update(flags_[ConversionsProcess], process == CMS::Conversions);
0213 update(flags_[EBremProcess], process == CMS::EBrem);
0214 update(flags_[SynchrotronRadiationProcess], process == CMS::SynchrotronRadiation);
0215 update(flags_[MuBremProcess], process == CMS::MuBrem);
0216 update(flags_[MuNuclProcess], process == CMS::MuNucl);
0217
0218
0219 for (TrackingVertex::tp_iterator iparticle = (*ivertex)->daughterTracks_begin();
0220 iparticle != (*ivertex)->daughterTracks_end();
0221 ++iparticle) {
0222 if ((*iparticle)->numberOfTrackerLayers()) {
0223
0224 if (process == CMS::Decay) {
0225
0226 HepPDT::ParticleID particleID(pdgid);
0227
0228 if (particleID.isValid()) {
0229
0230 ParticleData const *particleData = particleDataTable_->particle(particleID);
0231
0232 if (particleData) {
0233
0234 if (particleDataTable_->particle(particleID)->lifetime() > longLivedDecayLength_) {
0235
0236 update(flags_[BWeakDecay], particleID.hasBottom());
0237 update(flags_[CWeakDecay], particleID.hasCharm());
0238 update(flags_[LongLivedDecay], true);
0239 }
0240
0241 update(flags_[TauDecay], pdgid == 15);
0242 update(flags_[KsDecay], pdgid == 310);
0243 update(flags_[LambdaDecay], pdgid == 3122);
0244 update(flags_[JpsiDecay], pdgid == 443);
0245 update(flags_[XiDecay], pdgid == 3312);
0246 update(flags_[OmegaDecay], pdgid == 3334);
0247 update(flags_[SigmaPlusDecay], pdgid == 3222);
0248 update(flags_[SigmaMinusDecay], pdgid == 3112);
0249 }
0250 }
0251 }
0252 }
0253 }
0254 }
0255 }
0256
0257 void VertexClassifier::vertexInformation() {
0258
0259 typedef std::multimap<double, HepMC::ThreeVector> Clusters;
0260 typedef std::pair<double, HepMC::ThreeVector> ClusterPair;
0261
0262 Clusters clusters;
0263
0264
0265 GeneratedPrimaryVertex const &genpv = genpvs_.back();
0266
0267
0268 const VertexHistory::GenVertexTrail &genVertexTrail = tracer_.genVertexTrail();
0269
0270
0271 double const mm = 0.1;
0272
0273
0274 for (VertexHistory::GenVertexTrail::const_iterator ivertex = genVertexTrail.begin(); ivertex != genVertexTrail.end();
0275 ++ivertex) {
0276
0277 if (*ivertex) {
0278
0279 HepMC::ThreeVector p = (*ivertex)->point3d();
0280 double distance =
0281 sqrt(pow(p.x() * mm - genpv.x, 2) + pow(p.y() * mm - genpv.y, 2) + pow(p.z() * mm - genpv.z, 2));
0282
0283
0284 if (clusters.empty()) {
0285 clusters.insert(ClusterPair(distance, HepMC::ThreeVector(p.x() * mm, p.y() * mm, p.z() * mm)));
0286 continue;
0287 }
0288
0289
0290
0291 Clusters::const_iterator icluster = clusters.lower_bound(distance - vertexClusteringDistance_);
0292
0293 if (icluster == clusters.upper_bound(distance + vertexClusteringDistance_)) {
0294 clusters.insert(ClusterPair(distance, HepMC::ThreeVector(p.x() * mm, p.y() * mm, p.z() * mm)));
0295 continue;
0296 }
0297
0298 bool cluster = false;
0299
0300
0301
0302 for (; icluster != clusters.upper_bound(distance + vertexClusteringDistance_); ++icluster) {
0303 double difference = sqrt(pow(p.x() * mm - icluster->second.x(), 2) + pow(p.y() * mm - icluster->second.y(), 2) +
0304 pow(p.z() * mm - icluster->second.z(), 2));
0305
0306 if (difference < vertexClusteringDistance_) {
0307 cluster = true;
0308 break;
0309 }
0310 }
0311
0312 if (!cluster)
0313 clusters.insert(ClusterPair(distance, HepMC::ThreeVector(p.x() * mm, p.y() * mm, p.z() * mm)));
0314 }
0315 }
0316
0317 const VertexHistory::SimVertexTrail &simVertexTrail = tracer_.simVertexTrail();
0318
0319
0320 for (VertexHistory::SimVertexTrail::const_reverse_iterator ivertex = simVertexTrail.rbegin();
0321 ivertex != simVertexTrail.rend();
0322 ++ivertex) {
0323
0324 TrackingVertex::LorentzVector p = (*ivertex)->position();
0325
0326 double distance = sqrt(pow(p.x() - genpv.x, 2) + pow(p.y() - genpv.y, 2) + pow(p.z() - genpv.z, 2));
0327
0328
0329 if (clusters.empty()) {
0330 clusters.insert(ClusterPair(distance, HepMC::ThreeVector(p.x(), p.y(), p.z())));
0331 continue;
0332 }
0333
0334
0335
0336 Clusters::const_iterator icluster = clusters.lower_bound(distance - vertexClusteringDistance_);
0337
0338 if (icluster == clusters.upper_bound(distance + vertexClusteringDistance_)) {
0339 clusters.insert(ClusterPair(distance, HepMC::ThreeVector(p.x(), p.y(), p.z())));
0340 continue;
0341 }
0342
0343 bool cluster = false;
0344
0345
0346 for (; icluster != clusters.upper_bound(distance + vertexClusteringDistance_); ++icluster) {
0347 double difference = sqrt(pow(p.x() - icluster->second.x(), 2) + pow(p.y() - icluster->second.y(), 2) +
0348 pow(p.z() - icluster->second.z(), 2));
0349
0350 if (difference < vertexClusteringDistance_) {
0351 cluster = true;
0352 break;
0353 }
0354 }
0355
0356 if (!cluster)
0357 clusters.insert(ClusterPair(distance, HepMC::ThreeVector(p.x(), p.y(), p.z())));
0358 }
0359
0360 if (clusters.size() == 1)
0361 flags_[PrimaryVertex] = true;
0362 else if (clusters.size() == 2)
0363 flags_[SecondaryVertex] = true;
0364 else
0365 flags_[TertiaryVertex] = true;
0366 }
0367
0368 bool VertexClassifier::isFinalstateParticle(const HepMC::GenParticle *p) {
0369 return !p->end_vertex() && p->status() == 1;
0370 }
0371
0372 bool VertexClassifier::isCharged(const HepMC::GenParticle *p) {
0373 const ParticleData *part = particleDataTable_->particle(p->pdg_id());
0374 if (part)
0375 return part->charge() != 0;
0376 else {
0377
0378 return particleDataTable_->particle(-p->pdg_id()) != nullptr;
0379 }
0380 }
0381
0382 void VertexClassifier::genPrimaryVertices() {
0383 genpvs_.clear();
0384
0385 const HepMC::GenEvent *event = mcInformation_->GetEvent();
0386
0387 if (event) {
0388
0389 for (HepMC::GenEvent::vertex_const_iterator ivertex = event->vertices_begin(); ivertex != event->vertices_end();
0390 ++ivertex) {
0391 bool hasParentVertex = false;
0392
0393
0394
0395 for (HepMC::GenVertex::particle_iterator iparent = (*ivertex)->particles_begin(HepMC::parents);
0396 iparent != (*ivertex)->particles_end(HepMC::parents);
0397 ++iparent)
0398 if ((*iparent)->production_vertex()) {
0399 hasParentVertex = true;
0400 break;
0401 }
0402
0403
0404 if (hasParentVertex)
0405 continue;
0406
0407
0408 HepMC::FourVector pos = (*ivertex)->position();
0409
0410 double const mm = 0.1;
0411
0412 GeneratedPrimaryVertex pv(pos.x() * mm, pos.y() * mm, pos.z() * mm);
0413
0414 std::vector<GeneratedPrimaryVertex>::iterator ientry = genpvs_.begin();
0415
0416
0417 for (; ientry != genpvs_.end(); ++ientry) {
0418 double distance = sqrt(pow(pv.x - ientry->x, 2) + pow(pv.y - ientry->y, 2) + pow(pv.z - ientry->z, 2));
0419 if (distance < vertexClusteringDistance_)
0420 break;
0421 }
0422
0423
0424 if (ientry == genpvs_.end())
0425 ientry = genpvs_.insert(ientry, pv);
0426
0427
0428 ientry->genVertex.push_back((*ivertex)->barcode());
0429
0430
0431 for (HepMC::GenVertex::particle_iterator idecendants = (*ivertex)->particles_begin(HepMC::descendants);
0432 idecendants != (*ivertex)->particles_end(HepMC::descendants);
0433 ++idecendants) {
0434 if (isFinalstateParticle(*idecendants))
0435 if (find(ientry->finalstateParticles.begin(), ientry->finalstateParticles.end(), (*idecendants)->barcode()) ==
0436 ientry->finalstateParticles.end()) {
0437 ientry->finalstateParticles.push_back((*idecendants)->barcode());
0438 HepMC::FourVector m = (*idecendants)->momentum();
0439
0440 ientry->ptot.setPx(ientry->ptot.px() + m.px());
0441 ientry->ptot.setPy(ientry->ptot.py() + m.py());
0442 ientry->ptot.setPz(ientry->ptot.pz() + m.pz());
0443 ientry->ptot.setE(ientry->ptot.e() + m.e());
0444 ientry->ptsq += m.perp() * m.perp();
0445
0446 if (m.perp() > 0.8 && std::abs(m.pseudoRapidity()) < 2.5 && isCharged(*idecendants))
0447 ientry->nGenTrk++;
0448 }
0449 }
0450 }
0451 }
0452
0453 std::sort(genpvs_.begin(), genpvs_.end());
0454 }