File indexing completed on 2025-02-05 23:51:52
0001
0002 #include <cmath>
0003 #include <cstdlib>
0004 #include <iostream>
0005
0006 #include "HepPDT/ParticleID.hh"
0007
0008 #include "SimTracker/TrackHistory/interface/TrackClassifier.h"
0009
0010 #define update(a, b) \
0011 do { \
0012 (a) = (a) || (b); \
0013 } while (0)
0014
0015 TrackClassifier::TrackClassifier(edm::ParameterSet const &config, edm::ConsumesCollector &&collector)
0016 : TrackCategories(),
0017 hepMCLabel_(config.getUntrackedParameter<edm::InputTag>("hepMC")),
0018 beamSpotLabel_(config.getUntrackedParameter<edm::InputTag>("beamSpot")),
0019 tracer_(config, std::move(collector)),
0020 quality_(config, collector),
0021 magneticFieldToken_(collector.esConsumes<MagneticField, IdealMagneticFieldRecord>()),
0022 particleDataTableToken_(collector.esConsumes()),
0023 transientTrackBuilderToken_(collector.esConsumes<TransientTrackBuilder, TransientTrackRecord>()),
0024 tTopoHandToken_(collector.esConsumes<TrackerTopology, TrackerTopologyRcd>()) {
0025 collector.consumes<edm::HepMCProduct>(hepMCLabel_);
0026 collector.consumes<reco::BeamSpot>(beamSpotLabel_);
0027
0028
0029 tracer_.depth(-2);
0030
0031
0032 badPull_ = config.getUntrackedParameter<double>("badPull");
0033
0034
0035 longLivedDecayLength_ = config.getUntrackedParameter<double>("longLivedDecayLength");
0036
0037
0038 float vertexClusteringDistance = config.getUntrackedParameter<double>("vertexClusteringDistance");
0039 vertexClusteringSqDistance_ = vertexClusteringDistance * vertexClusteringDistance;
0040
0041
0042 numberOfInnerLayers_ = config.getUntrackedParameter<unsigned int>("numberOfInnerLayers");
0043
0044
0045 minTrackerSimHits_ = config.getUntrackedParameter<unsigned int>("minTrackerSimHits");
0046 }
0047
0048 void TrackClassifier::newEvent(edm::Event const &event, edm::EventSetup const &setup) {
0049
0050 tracer_.newEvent(event, setup);
0051
0052
0053 quality_.newEvent(event, setup);
0054
0055
0056 event.getByLabel(hepMCLabel_, mcInformation_);
0057
0058
0059 magneticField_ = setup.getHandle(magneticFieldToken_);
0060
0061
0062 particleDataTable_ = setup.getHandle(particleDataTableToken_);
0063
0064
0065 event.getByLabel(beamSpotLabel_, beamSpot_);
0066
0067
0068 transientTrackBuilder_ = setup.getHandle(transientTrackBuilderToken_);
0069
0070
0071 genPrimaryVertices();
0072
0073
0074 tTopo_ = &setup.getData(tTopoHandToken_);
0075 }
0076
0077 TrackClassifier const &TrackClassifier::evaluate(reco::TrackBaseRef const &track) {
0078
0079 reset();
0080
0081
0082 if (tracer_.evaluate(track)) {
0083
0084
0085 reconstructionInformation(track);
0086
0087
0088 simulationInformation();
0089
0090
0091 qualityInformation(track);
0092
0093
0094 hadronFlavor();
0095
0096
0097 processesAtGenerator();
0098
0099
0100 processesAtSimulation();
0101
0102
0103 vertexInformation();
0104
0105
0106 unknownTrack();
0107 } else
0108 flags_[Fake] = true;
0109
0110 return *this;
0111 }
0112
0113 void TrackClassifier::fillPSetDescription(edm::ParameterSetDescription &desc) {
0114 TrackHistory::fillPSetDescription(desc);
0115
0116 edm::ParameterSetDescription trackQualityDesc;
0117 TrackQuality::fillPSetDescription(trackQualityDesc);
0118 desc.add<edm::ParameterSetDescription>("hitAssociator", trackQualityDesc);
0119
0120 desc.addUntracked<edm::InputTag>("hepMC", edm::InputTag("generatorSmeared"));
0121 desc.addUntracked<edm::InputTag>("beamSpot", edm::InputTag("offlineBeamSpot"));
0122 desc.addUntracked<double>("badPull", 3.0);
0123 desc.addUntracked<double>("longLivedDecayLength", 1e-14);
0124 desc.addUntracked<double>("vertexClusteringDistance", 0.0001);
0125 desc.addUntracked<unsigned int>("numberOfInnerLayers", 2);
0126 desc.addUntracked<unsigned int>("minTrackerSimHits", 3);
0127 }
0128
0129 TrackClassifier const &TrackClassifier::evaluate(TrackingParticleRef const &track) {
0130
0131 reset();
0132
0133
0134 tracer_.evaluate(track);
0135
0136
0137 const reco::TrackBaseRef &recotrack = tracer_.recoTrack();
0138
0139
0140 if (recotrack.isNonnull()) {
0141 flags_[Reconstructed] = true;
0142
0143
0144 reconstructionInformation(recotrack);
0145
0146 qualityInformation(recotrack);
0147 } else
0148 flags_[Reconstructed] = false;
0149
0150
0151 simulationInformation();
0152
0153
0154 hadronFlavor();
0155
0156
0157 processesAtGenerator();
0158
0159
0160 processesAtSimulation();
0161
0162
0163 vertexInformation();
0164
0165
0166 unknownTrack();
0167
0168 return *this;
0169 }
0170
0171 void TrackClassifier::reconstructionInformation(reco::TrackBaseRef const &track) {
0172 TrackingParticleRef tpr = tracer_.simParticle();
0173
0174
0175
0176
0177 const SimTrack *assocTrack = &(*tpr->g4Track_begin());
0178
0179 FreeTrajectoryState ftsAtProduction(
0180 GlobalPoint(tpr->vertex().x(), tpr->vertex().y(), tpr->vertex().z()),
0181 GlobalVector(assocTrack->momentum().x(), assocTrack->momentum().y(), assocTrack->momentum().z()),
0182 TrackCharge(track->charge()),
0183 magneticField_.product());
0184
0185 try {
0186 TSCPBuilderNoMaterial tscpBuilder;
0187 TrajectoryStateClosestToPoint tsAtClosestApproach =
0188 tscpBuilder(ftsAtProduction, GlobalPoint(beamSpot_->x0(), beamSpot_->y0(), beamSpot_->z0()));
0189
0190 GlobalVector v =
0191 tsAtClosestApproach.theState().position() - GlobalPoint(beamSpot_->x0(), beamSpot_->y0(), beamSpot_->z0());
0192 GlobalVector p = tsAtClosestApproach.theState().momentum();
0193
0194
0195 double dxySim = -v.x() * sin(p.phi()) + v.y() * cos(p.phi());
0196
0197
0198 double dzSim = v.z() - (v.x() * p.x() + v.y() * p.y()) * p.z() / p.perp2();
0199
0200
0201 double dxyPull =
0202 std::abs(track->dxy(reco::TrackBase::Point(beamSpot_->x0(), beamSpot_->y0(), beamSpot_->z0())) - dxySim) /
0203 track->dxyError();
0204
0205
0206 double dzPull =
0207 std::abs(track->dz(reco::TrackBase::Point(beamSpot_->x0(), beamSpot_->y0(), beamSpot_->z0())) - dzSim) /
0208 track->dzError();
0209
0210
0211 flags_[Bad] = (dxyPull > badPull_ || dzPull > badPull_);
0212
0213 } catch (cms::Exception const &) {
0214 flags_[Bad] = true;
0215 }
0216 }
0217
0218 void TrackClassifier::simulationInformation() {
0219
0220 EncodedEventId eventId = tracer_.simParticle()->eventId();
0221
0222 flags_[SignalEvent] = !eventId.bunchCrossing() && !eventId.event();
0223
0224 flags_[Muon] = (abs(tracer_.simParticle()->pdgId()) == 13);
0225
0226 flags_[TrackerSimHits] = tracer_.simParticle()->numberOfTrackerLayers() >= (int)minTrackerSimHits_;
0227 }
0228
0229 void TrackClassifier::qualityInformation(reco::TrackBaseRef const &track) {
0230
0231 quality_.evaluate(tracer_.simParticleTrail(), track, tTopo_);
0232
0233 unsigned int maxLayers = std::min(numberOfInnerLayers_, quality_.numberOfLayers());
0234
0235
0236 for (unsigned int i = 0; i < maxLayers; i++) {
0237 const TrackQuality::Layer &layer = quality_.layer(i);
0238
0239
0240 for (unsigned int j = 0; j < layer.hits.size(); j++) {
0241 const TrackQuality::Layer::Hit &hit = layer.hits[j];
0242
0243
0244 if (hit.state == TrackQuality::Layer::Noise || hit.state == TrackQuality::Layer::Misassoc)
0245 flags_[BadInnerHits] = true;
0246 else if (hit.state == TrackQuality::Layer::Shared)
0247 flags_[SharedInnerHits] = true;
0248 }
0249 }
0250 }
0251
0252 void TrackClassifier::hadronFlavor() {
0253
0254 const reco::GenParticle *particle = tracer_.recoGenParticle();
0255
0256
0257 if (particle) {
0258 HepPDT::ParticleID pid(particle->pdgId());
0259 flags_[Bottom] = pid.hasBottom();
0260 flags_[Charm] = pid.hasCharm();
0261 flags_[Light] = !pid.hasCharm() && !pid.hasBottom();
0262 }
0263 }
0264
0265 void TrackClassifier::processesAtGenerator() {
0266
0267 int pdgid = 0;
0268
0269
0270
0271 TrackHistory::RecoGenParticleTrail const &recoGenParticleTrail = tracer_.recoGenParticleTrail();
0272
0273
0274
0275 for (TrackHistory::RecoGenParticleTrail::const_iterator iparticle = recoGenParticleTrail.begin();
0276 iparticle != recoGenParticleTrail.end();
0277 ++iparticle) {
0278 pdgid = std::abs((*iparticle)->pdgId());
0279
0280 HepPDT::ParticleID particleID(pdgid);
0281
0282
0283 if (particleID.isValid()) {
0284
0285 ParticleData const *particleData = particleDataTable_->particle(particleID);
0286
0287 if (particleData) {
0288
0289 if (particleData->lifetime() > longLivedDecayLength_)
0290 update(flags_[LongLivedDecay], true);
0291
0292 update(flags_[BWeakDecay], particleID.hasBottom());
0293 update(flags_[CWeakDecay], particleID.hasCharm());
0294
0295 std::set<int> daughterIds;
0296 size_t ndau = (*iparticle)->numberOfDaughters();
0297 for (size_t i = 0; i < ndau; ++i) {
0298 daughterIds.insert((*iparticle)->daughter(i)->pdgId());
0299 }
0300 update(flags_[FromBWeakDecayMuon], particleID.hasBottom() && (daughterIds.find(13) != daughterIds.end()));
0301 update(flags_[FromCWeakDecayMuon], particleID.hasCharm() && (daughterIds.find(13) != daughterIds.end()));
0302 }
0303
0304 update(flags_[ChargePionDecay], pdgid == 211);
0305 update(flags_[ChargeKaonDecay], pdgid == 321);
0306 update(flags_[TauDecay], pdgid == 15);
0307 update(flags_[KsDecay], pdgid == 310);
0308 update(flags_[LambdaDecay], pdgid == 3122);
0309 update(flags_[JpsiDecay], pdgid == 443);
0310 update(flags_[XiDecay], pdgid == 3312);
0311 update(flags_[SigmaPlusDecay], pdgid == 3222);
0312 update(flags_[SigmaMinusDecay], pdgid == 3112);
0313 }
0314 }
0315
0316 update(flags_[FromChargePionMuon], flags_[Muon] && flags_[ChargePionDecay]);
0317 update(flags_[FromChargeKaonMuon], flags_[Muon] && flags_[ChargeKaonDecay]);
0318 update(flags_[DecayOnFlightMuon], (flags_[FromChargePionMuon] || flags_[FromChargeKaonMuon]));
0319 }
0320
0321 void TrackClassifier::processesAtSimulation() {
0322 TrackHistory::SimParticleTrail const &simParticleTrail = tracer_.simParticleTrail();
0323
0324
0325 for (TrackHistory::SimParticleTrail::const_iterator iparticle = simParticleTrail.begin();
0326 iparticle != simParticleTrail.end();
0327 ++iparticle) {
0328
0329 int pdgid = 0;
0330
0331
0332 TrackingVertexRef const &parentVertex = (*iparticle)->parentVertex();
0333
0334
0335 if (parentVertex.isNonnull()) {
0336
0337 bool flag = false;
0338 TrackingVertex::tp_iterator itd, its;
0339
0340 for (its = parentVertex->sourceTracks_begin(); its != parentVertex->sourceTracks_end(); ++its) {
0341 for (itd = parentVertex->daughterTracks_begin(); itd != parentVertex->daughterTracks_end(); ++itd)
0342 if (itd != its) {
0343 flag = true;
0344 break;
0345 }
0346 if (flag)
0347 break;
0348 }
0349
0350
0351 if (its != parentVertex->sourceTracks_end())
0352 pdgid = std::abs((*its)->pdgId());
0353 else
0354 pdgid = 0;
0355 }
0356
0357 unsigned int processG4 = 0;
0358
0359
0360 if (parentVertex->nG4Vertices() > 0) {
0361 processG4 = (*(parentVertex->g4Vertices_begin())).processType();
0362 }
0363
0364 unsigned int process = g4toCMSProcMap_.processId(processG4);
0365
0366
0367 update(flags_[KnownProcess], process != CMS::Undefined && process != CMS::Unknown && process != CMS::Primary);
0368
0369 update(flags_[UndefinedProcess], process == CMS::Undefined);
0370 update(flags_[UnknownProcess], process == CMS::Unknown);
0371 update(flags_[PrimaryProcess], process == CMS::Primary);
0372 update(flags_[HadronicProcess], process == CMS::Hadronic);
0373 update(flags_[DecayProcess], process == CMS::Decay);
0374 update(flags_[ComptonProcess], process == CMS::Compton);
0375 update(flags_[AnnihilationProcess], process == CMS::Annihilation);
0376 update(flags_[EIoniProcess], process == CMS::EIoni);
0377 update(flags_[HIoniProcess], process == CMS::HIoni);
0378 update(flags_[MuIoniProcess], process == CMS::MuIoni);
0379 update(flags_[PhotonProcess], process == CMS::Photon);
0380 update(flags_[MuPairProdProcess], process == CMS::MuPairProd);
0381 update(flags_[ConversionsProcess], process == CMS::Conversions);
0382 update(flags_[EBremProcess], process == CMS::EBrem);
0383 update(flags_[SynchrotronRadiationProcess], process == CMS::SynchrotronRadiation);
0384 update(flags_[MuBremProcess], process == CMS::MuBrem);
0385 update(flags_[MuNuclProcess], process == CMS::MuNucl);
0386
0387
0388 HepPDT::ParticleID particleID(pdgid);
0389
0390
0391 if (particleID.isValid()) {
0392
0393 ParticleData const *particleData = particleDataTable_->particle(particleID);
0394
0395 if (process == CMS::Decay) {
0396
0397 if (particleData) {
0398
0399 if (particleDataTable_->particle(particleID)->lifetime() > longLivedDecayLength_)
0400 update(flags_[LongLivedDecay], true);
0401
0402
0403 update(flags_[BWeakDecay], particleID.hasBottom());
0404 update(flags_[CWeakDecay], particleID.hasCharm());
0405
0406
0407 int daughtId = abs((*iparticle)->pdgId());
0408 update(flags_[FromBWeakDecayMuon], particleID.hasBottom() && daughtId == 13);
0409 update(flags_[FromCWeakDecayMuon], particleID.hasCharm() && daughtId == 13);
0410 }
0411
0412 update(flags_[ChargePionDecay], pdgid == 211);
0413 update(flags_[ChargeKaonDecay], pdgid == 321);
0414 update(flags_[TauDecay], pdgid == 15);
0415 update(flags_[KsDecay], pdgid == 310);
0416 update(flags_[LambdaDecay], pdgid == 3122);
0417 update(flags_[JpsiDecay], pdgid == 443);
0418 update(flags_[XiDecay], pdgid == 3312);
0419 update(flags_[OmegaDecay], pdgid == 3334);
0420 update(flags_[SigmaPlusDecay], pdgid == 3222);
0421 update(flags_[SigmaMinusDecay], pdgid == 3112);
0422 }
0423 }
0424 }
0425
0426 update(flags_[FromChargePionMuon], flags_[Muon] && flags_[ChargePionDecay]);
0427 update(flags_[FromChargeKaonMuon], flags_[Muon] && flags_[ChargeKaonDecay]);
0428 update(flags_[DecayOnFlightMuon], flags_[FromChargePionMuon] || flags_[FromChargeKaonMuon]);
0429 }
0430
0431 void TrackClassifier::vertexInformation() {
0432
0433 GeneratedPrimaryVertex const &genpv = genpvs_.back();
0434
0435
0436 const TrackHistory::GenParticleTrail &genParticleTrail = tracer_.genParticleTrail();
0437
0438
0439 int counter = 0;
0440
0441
0442 double const mm = 0.1;
0443
0444 double oldX = genpv.x;
0445 double oldY = genpv.y;
0446 double oldZ = genpv.z;
0447
0448
0449 for (TrackHistory::GenParticleTrail::const_reverse_iterator iparticle = genParticleTrail.rbegin();
0450 iparticle != genParticleTrail.rend();
0451 ++iparticle) {
0452
0453 HepMC::GenVertex *parent = (*iparticle)->production_vertex();
0454 if (parent) {
0455 HepMC::ThreeVector p = parent->point3d();
0456
0457 double distance2 = pow(p.x() * mm - genpv.x, 2) + pow(p.y() * mm - genpv.y, 2) + pow(p.z() * mm - genpv.z, 2);
0458 double difference2 = pow(p.x() * mm - oldX, 2) + pow(p.y() * mm - oldY, 2) + pow(p.z() * mm - oldZ, 2);
0459
0460
0461
0462
0463
0464 if (difference2 > vertexClusteringSqDistance_) {
0465 if (distance2 > vertexClusteringSqDistance_)
0466 counter++;
0467 oldX = p.x() * mm;
0468 oldY = p.y() * mm;
0469 oldZ = p.z() * mm;
0470 }
0471 }
0472 }
0473
0474 const TrackHistory::SimParticleTrail &simParticleTrail = tracer_.simParticleTrail();
0475
0476
0477 for (TrackHistory::SimParticleTrail::const_reverse_iterator iparticle = simParticleTrail.rbegin();
0478 iparticle != simParticleTrail.rend();
0479 ++iparticle) {
0480
0481 TrackingParticle::Point p = (*iparticle)->vertex();
0482
0483 double distance2 = pow(p.x() - genpv.x, 2) + pow(p.y() - genpv.y, 2) + pow(p.z() - genpv.z, 2);
0484 double difference2 = pow(p.x() - oldX, 2) + pow(p.y() - oldY, 2) + pow(p.z() - oldZ, 2);
0485
0486
0487
0488
0489
0490 if (difference2 > vertexClusteringSqDistance_) {
0491 if (distance2 > vertexClusteringSqDistance_)
0492 counter++;
0493 oldX = p.x();
0494 oldY = p.y();
0495 oldZ = p.z();
0496 }
0497 }
0498
0499 if (!counter)
0500 flags_[PrimaryVertex] = true;
0501 else if (counter == 1)
0502 flags_[SecondaryVertex] = true;
0503 else
0504 flags_[TertiaryVertex] = true;
0505 }
0506
0507 bool TrackClassifier::isFinalstateParticle(const HepMC::GenParticle *p) { return !p->end_vertex() && p->status() == 1; }
0508
0509 bool TrackClassifier::isCharged(const HepMC::GenParticle *p) {
0510 const ParticleData *part = particleDataTable_->particle(p->pdg_id());
0511 if (part)
0512 return part->charge() != 0;
0513 else {
0514
0515 return particleDataTable_->particle(-p->pdg_id()) != nullptr;
0516 }
0517 }
0518
0519 void TrackClassifier::genPrimaryVertices() {
0520 genpvs_.clear();
0521
0522 const HepMC::GenEvent *event = mcInformation_->GetEvent();
0523
0524 if (event) {
0525
0526 for (HepMC::GenEvent::vertex_const_iterator ivertex = event->vertices_begin(); ivertex != event->vertices_end();
0527 ++ivertex) {
0528 bool hasParentVertex = false;
0529
0530
0531
0532 for (HepMC::GenVertex::particle_iterator iparent = (*ivertex)->particles_begin(HepMC::parents);
0533 iparent != (*ivertex)->particles_end(HepMC::parents);
0534 ++iparent)
0535 if ((*iparent)->production_vertex()) {
0536 hasParentVertex = true;
0537 break;
0538 }
0539
0540
0541 if (hasParentVertex)
0542 continue;
0543
0544
0545 HepMC::FourVector pos = (*ivertex)->position();
0546
0547 double const mm = 0.1;
0548
0549 GeneratedPrimaryVertex pv(pos.x() * mm, pos.y() * mm, pos.z() * mm);
0550
0551 std::vector<GeneratedPrimaryVertex>::iterator ientry = genpvs_.begin();
0552
0553
0554 for (; ientry != genpvs_.end(); ++ientry) {
0555 double distance2 = pow(pv.x - ientry->x, 2) + pow(pv.y - ientry->y, 2) + pow(pv.z - ientry->z, 2);
0556 if (distance2 < vertexClusteringSqDistance_)
0557 break;
0558 }
0559
0560
0561 if (ientry == genpvs_.end())
0562 ientry = genpvs_.insert(ientry, pv);
0563
0564
0565 ientry->genVertex.push_back((*ivertex)->barcode());
0566
0567
0568 for (HepMC::GenVertex::particle_iterator idecendants = (*ivertex)->particles_begin(HepMC::descendants);
0569 idecendants != (*ivertex)->particles_end(HepMC::descendants);
0570 ++idecendants) {
0571 if (isFinalstateParticle(*idecendants))
0572 if (find(ientry->finalstateParticles.begin(), ientry->finalstateParticles.end(), (*idecendants)->barcode()) ==
0573 ientry->finalstateParticles.end()) {
0574 ientry->finalstateParticles.push_back((*idecendants)->barcode());
0575 HepMC::FourVector m = (*idecendants)->momentum();
0576
0577 ientry->ptot.setPx(ientry->ptot.px() + m.px());
0578 ientry->ptot.setPy(ientry->ptot.py() + m.py());
0579 ientry->ptot.setPz(ientry->ptot.pz() + m.pz());
0580 ientry->ptot.setE(ientry->ptot.e() + m.e());
0581 ientry->ptsq += m.perp() * m.perp();
0582
0583 if (m.perp() > 0.8 && std::abs(m.pseudoRapidity()) < 2.5 && isCharged(*idecendants))
0584 ientry->nGenTrk++;
0585 }
0586 }
0587 }
0588 }
0589
0590 std::sort(genpvs_.begin(), genpvs_.end());
0591 }