00001
00002 #include <math.h>
00003 #include <cstdlib>
00004 #include <iostream>
00005
00006 #include "HepPDT/ParticleID.hh"
00007
00008 #include "SimTracker/TrackHistory/interface/TrackClassifier.h"
00009
00010 #define update(a, b) do { (a) = (a) | (b); } while(0)
00011
00012 TrackClassifier::TrackClassifier(edm::ParameterSet const & config) : TrackCategories(),
00013 hepMCLabel_( config.getUntrackedParameter<edm::InputTag>("hepMC") ),
00014 beamSpotLabel_( config.getUntrackedParameter<edm::InputTag>("beamSpot") ),
00015 tracer_(config),
00016 quality_(config)
00017 {
00018
00019 tracer_.depth(-2);
00020
00021
00022 badPull_ = config.getUntrackedParameter<double>("badPull");
00023
00024
00025 longLivedDecayLength_ = config.getUntrackedParameter<double>("longLivedDecayLength");
00026
00027
00028 float vertexClusteringDistance = config.getUntrackedParameter<double>("vertexClusteringDistance");
00029 vertexClusteringSqDistance_ = vertexClusteringDistance * vertexClusteringDistance;
00030
00031
00032 numberOfInnerLayers_ = config.getUntrackedParameter<unsigned int>("numberOfInnerLayers");
00033
00034
00035 minTrackerSimHits_ = config.getUntrackedParameter<unsigned int>("minTrackerSimHits");
00036 }
00037
00038
00039 void TrackClassifier::newEvent ( edm::Event const & event, edm::EventSetup const & setup )
00040 {
00041
00042 tracer_.newEvent(event, setup);
00043
00044
00045 quality_.newEvent(event, setup);
00046
00047
00048 event.getByLabel(hepMCLabel_, mcInformation_);
00049
00050
00051 setup.get<IdealMagneticFieldRecord>().get(magneticField_);
00052
00053
00054 setup.getData(particleDataTable_);
00055
00056
00057 event.getByLabel(beamSpotLabel_, beamSpot_);
00058
00059
00060 setup.get<TransientTrackRecord>().get("TransientTrackBuilder", transientTrackBuilder_);
00061
00062
00063 genPrimaryVertices();
00064
00065
00066 edm::ESHandle<TrackerTopology> tTopoHand;
00067 setup.get<IdealGeometryRecord>().get(tTopoHand);
00068 tTopo_=tTopoHand.product();
00069 }
00070
00071
00072 TrackClassifier const & TrackClassifier::evaluate (reco::TrackBaseRef const & track)
00073 {
00074
00075 reset();
00076
00077
00078 if ( tracer_.evaluate(track) )
00079 {
00080
00081 reconstructionInformation(track);
00082
00083
00084 simulationInformation();
00085
00086
00087 qualityInformation(track);
00088
00089
00090 hadronFlavor();
00091
00092
00093 processesAtGenerator();
00094
00095
00096 processesAtSimulation();
00097
00098
00099 vertexInformation();
00100
00101
00102 unknownTrack();
00103 }
00104 else
00105 flags_[Fake] = true;
00106
00107 return *this;
00108 }
00109
00110
00111 TrackClassifier const & TrackClassifier::evaluate (TrackingParticleRef const & track)
00112 {
00113
00114 reset();
00115
00116
00117 tracer_.evaluate(track);
00118
00119
00120 const reco::TrackBaseRef & recotrack = tracer_.recoTrack();
00121
00122
00123 if ( recotrack.isNonnull() )
00124 {
00125 flags_[Reconstructed] = true;
00126
00127 reconstructionInformation(recotrack);
00128
00129 qualityInformation(recotrack);
00130 }
00131 else
00132 flags_[Reconstructed] = false;
00133
00134
00135 simulationInformation();
00136
00137
00138 hadronFlavor();
00139
00140
00141 processesAtGenerator();
00142
00143
00144 processesAtSimulation();
00145
00146
00147 vertexInformation();
00148
00149
00150 unknownTrack();
00151
00152 return *this;
00153 }
00154
00155
00156 void TrackClassifier::reconstructionInformation(reco::TrackBaseRef const & track)
00157 {
00158 TrackingParticleRef tpr = tracer_.simParticle();
00159
00160
00161
00162 const SimTrack * assocTrack = &(*tpr->g4Track_begin());
00163
00164 FreeTrajectoryState ftsAtProduction(
00165 GlobalPoint(
00166 tpr->vertex().x(),
00167 tpr->vertex().y(),
00168 tpr->vertex().z()
00169 ),
00170 GlobalVector(
00171 assocTrack->momentum().x(),
00172 assocTrack->momentum().y(),
00173 assocTrack->momentum().z()
00174 ),
00175 TrackCharge(track->charge()),
00176 magneticField_.product()
00177 );
00178
00179 try
00180 {
00181 TSCPBuilderNoMaterial tscpBuilder;
00182 TrajectoryStateClosestToPoint tsAtClosestApproach = tscpBuilder(
00183 ftsAtProduction,
00184 GlobalPoint(beamSpot_->x0(), beamSpot_->y0(), beamSpot_->z0())
00185 );
00186
00187 GlobalVector v = tsAtClosestApproach.theState().position()
00188 - GlobalPoint(beamSpot_->x0(), beamSpot_->y0(), beamSpot_->z0());
00189 GlobalVector p = tsAtClosestApproach.theState().momentum();
00190
00191
00192 double dxySim = -v.x()*sin(p.phi()) + v.y()*cos(p.phi());
00193
00194
00195 double dzSim = v.z() - (v.x()*p.x() + v.y()*p.y())*p.z()/p.perp2();
00196
00197
00198 double dxyPull = std::abs(
00199 track->dxy( reco::TrackBase::Point(beamSpot_->x0(), beamSpot_->y0(), beamSpot_->z0()) ) - dxySim
00200 ) / track->dxyError();
00201
00202
00203 double dzPull = std::abs(
00204 track->dz( reco::TrackBase::Point(beamSpot_->x0(), beamSpot_->y0(), beamSpot_->z0()) ) - dzSim
00205 ) / track->dzError();
00206
00207
00208 flags_[Bad] = (dxyPull > badPull_ || dzPull > badPull_);
00209
00210 }
00211 catch (cms::Exception exception)
00212 {
00213 flags_[Bad] = true;
00214 }
00215 }
00216
00217
00218 void TrackClassifier::simulationInformation()
00219 {
00220
00221 EncodedEventId eventId = tracer_.simParticle()->eventId();
00222
00223 flags_[SignalEvent] = !eventId.bunchCrossing() && !eventId.event();
00224
00225 flags_[Muon] = (abs(tracer_.simParticle()->pdgId()) == 13);
00226
00227 flags_[TrackerSimHits] = tracer_.simParticle()->numberOfTrackerLayers() >= (int)minTrackerSimHits_;
00228 }
00229
00230
00231 void TrackClassifier::qualityInformation(reco::TrackBaseRef const & track)
00232 {
00233
00234 quality_.evaluate(tracer_.simParticleTrail(), track, tTopo_);
00235
00236 unsigned int maxLayers = std::min(numberOfInnerLayers_, quality_.numberOfLayers());
00237
00238
00239 for (unsigned int i = 0; i < maxLayers; i++)
00240 {
00241 const TrackQuality::Layer &layer = quality_.layer(i);
00242
00243
00244 for (unsigned int j = 0; j < layer.hits.size(); j++)
00245 {
00246 const TrackQuality::Layer::Hit &hit = layer.hits[j];
00247
00248
00249 if (hit.state == TrackQuality::Layer::Noise ||
00250 hit.state == TrackQuality::Layer::Misassoc)
00251 flags_[BadInnerHits] = true;
00252 else if (hit.state == TrackQuality::Layer::Shared)
00253 flags_[SharedInnerHits] = true;
00254 }
00255 }
00256 }
00257
00258
00259 void TrackClassifier::hadronFlavor()
00260 {
00261
00262 const HepMC::GenParticle * particle = tracer_.genParticle();
00263
00264
00265 if (particle)
00266 {
00267 HepPDT::ParticleID pid(particle->pdg_id());
00268 flags_[Bottom] = pid.hasBottom();
00269 flags_[Charm] = pid.hasCharm();
00270 flags_[Light] = !pid.hasCharm() && !pid.hasBottom();
00271 }
00272 }
00273
00274
00275 void TrackClassifier::processesAtGenerator()
00276 {
00277
00278 int pdgid = 0;
00279
00280
00281 TrackHistory::GenParticleTrail const & genParticleTrail = tracer_.genParticleTrail();
00282
00283
00284 for (TrackHistory::GenParticleTrail::const_iterator iparticle = genParticleTrail.begin(); iparticle != genParticleTrail.end(); ++iparticle)
00285 {
00286
00287 HepMC::GenVertex * productionVertex = (*iparticle)->production_vertex();
00288
00289
00290
00291
00292
00293 if (productionVertex)
00294 {
00295
00296 if ( productionVertex->particles_in_size() == 1 )
00297 {
00298
00299 pdgid = std::abs((*productionVertex->particles_in_const_begin())->pdg_id());
00300
00301 HepPDT::ParticleID particleID(pdgid);
00302
00303
00304 if (particleID.isValid())
00305 {
00306
00307 ParticleData const * particleData = particleDataTable_->particle(particleID);
00308
00309 if (particleData)
00310 {
00311
00312 if ( particleData->lifetime() > longLivedDecayLength_ )
00313 update(flags_[LongLivedDecay], true);
00314
00315 update(flags_[BWeakDecay], particleID.hasBottom());
00316 update(flags_[CWeakDecay], particleID.hasCharm());
00317
00318 int daughterId = abs((*iparticle)->pdg_id());
00319 update(flags_[FromBWeakDecayMuon], particleID.hasBottom() && daughterId == 13);
00320 update(flags_[FromCWeakDecayMuon], particleID.hasCharm() && daughterId == 13);
00321 }
00322
00323 update(flags_[ChargePionDecay], pdgid == 211);
00324 update(flags_[ChargeKaonDecay], pdgid == 321);
00325 update(flags_[TauDecay], pdgid == 15);
00326 update(flags_[KsDecay], pdgid == 310);
00327 update(flags_[LambdaDecay], pdgid == 3122);
00328 update(flags_[JpsiDecay], pdgid == 443);
00329 update(flags_[XiDecay], pdgid == 3312);
00330 update(flags_[SigmaPlusDecay], pdgid == 3222);
00331 update(flags_[SigmaMinusDecay], pdgid == 3112);
00332 }
00333 }
00334 }
00335 }
00336
00337 update(flags_[FromChargePionMuon], flags_[Muon] && flags_[ChargePionDecay]);
00338 update(flags_[FromChargeKaonMuon], flags_[Muon] && flags_[ChargeKaonDecay]);
00339 update(flags_[DecayOnFlightMuon], (flags_[FromChargePionMuon] || flags_[FromChargeKaonMuon]));
00340 }
00341
00342
00343 void TrackClassifier::processesAtSimulation()
00344 {
00345 TrackHistory::SimParticleTrail const & simParticleTrail = tracer_.simParticleTrail();
00346
00347
00348 for (
00349 TrackHistory::SimParticleTrail::const_iterator iparticle = simParticleTrail.begin();
00350 iparticle != simParticleTrail.end();
00351 ++iparticle
00352 )
00353 {
00354
00355 int pdgid = 0;
00356
00357
00358 TrackingVertexRef const & parentVertex = (*iparticle)->parentVertex();
00359
00360
00361 if ( parentVertex.isNonnull() )
00362 {
00363
00364 bool flag = false;
00365 TrackingVertex::tp_iterator itd, its;
00366
00367 for (its = parentVertex->sourceTracks_begin(); its != parentVertex->sourceTracks_end(); ++its)
00368 {
00369 for (itd = parentVertex->daughterTracks_begin(); itd != parentVertex->daughterTracks_end(); ++itd)
00370 if (itd != its)
00371 {
00372 flag = true;
00373 break;
00374 }
00375 if (flag)
00376 break;
00377 }
00378
00379
00380 if ( its != parentVertex->sourceTracks_end() )
00381 pdgid = std::abs((*its)->pdgId());
00382 else
00383 pdgid = 0;
00384 }
00385
00386
00387 if ((*iparticle)->numberOfHits()) continue;
00388
00389
00390 #warning "This file has been modified just to get it to compile without any regard as to whether it still functions as intended"
00391 #ifdef REMOVED_JUST_TO_GET_IT_TO_COMPILE__THIS_CODE_NEEDS_TO_BE_CHECKED
00392 unsigned short process = (*iparticle)->pSimHit_begin()->processType();
00393 #else
00394 unsigned short process = 0;
00395 #endif
00396
00397
00398 update(
00399 flags_[KnownProcess],
00400 process != G4::Undefined &&
00401 process != G4::Unknown &&
00402 process != G4::Primary
00403 );
00404
00405 update(flags_[UndefinedProcess], process == G4::Undefined);
00406 update(flags_[UnknownProcess], process == G4::Unknown);
00407 update(flags_[PrimaryProcess], process == G4::Primary);
00408 update(flags_[HadronicProcess], process == G4::Hadronic);
00409 update(flags_[DecayProcess], process == G4::Decay);
00410 update(flags_[ComptonProcess], process == G4::Compton);
00411 update(flags_[AnnihilationProcess], process == G4::Annihilation);
00412 update(flags_[EIoniProcess], process == G4::EIoni);
00413 update(flags_[HIoniProcess], process == G4::HIoni);
00414 update(flags_[MuIoniProcess], process == G4::MuIoni);
00415 update(flags_[PhotonProcess], process == G4::Photon);
00416 update(flags_[MuPairProdProcess], process == G4::MuPairProd);
00417 update(flags_[ConversionsProcess], process == G4::Conversions);
00418 update(flags_[EBremProcess], process == G4::EBrem);
00419 update(flags_[SynchrotronRadiationProcess], process == G4::SynchrotronRadiation);
00420 update(flags_[MuBremProcess], process == G4::MuBrem);
00421 update(flags_[MuNuclProcess], process == G4::MuNucl);
00422
00423
00424 HepPDT::ParticleID particleID(pdgid);
00425
00426
00427 if (particleID.isValid())
00428 {
00429
00430 ParticleData const * particleData = particleDataTable_->particle(particleID);
00431
00432 if (process == G4::Decay)
00433 {
00434
00435 if (particleData)
00436 {
00437
00438 if ( particleDataTable_->particle(particleID)->lifetime() > longLivedDecayLength_ )
00439 update(flags_[LongLivedDecay], true);
00440
00441
00442 update(flags_[BWeakDecay], particleID.hasBottom());
00443 update(flags_[CWeakDecay], particleID.hasCharm());
00444
00445
00446 int daughtId = abs((*iparticle)->pdgId());
00447 update(flags_[FromBWeakDecayMuon], particleID.hasBottom() && daughtId == 13);
00448 update(flags_[FromCWeakDecayMuon], particleID.hasCharm() && daughtId == 13);
00449 }
00450
00451 update(flags_[ChargePionDecay], pdgid == 211);
00452 update(flags_[ChargeKaonDecay], pdgid == 321);
00453 update(flags_[TauDecay], pdgid == 15);
00454 update(flags_[KsDecay], pdgid == 310);
00455 update(flags_[LambdaDecay], pdgid == 3122);
00456 update(flags_[JpsiDecay], pdgid == 443);
00457 update(flags_[XiDecay], pdgid == 3312);
00458 update(flags_[OmegaDecay], pdgid == 3334);
00459 update(flags_[SigmaPlusDecay], pdgid == 3222);
00460 update(flags_[SigmaMinusDecay], pdgid == 3112);
00461 }
00462 }
00463 }
00464
00465 update(flags_[FromChargePionMuon], flags_[Muon] && flags_[ChargePionDecay]);
00466 update(flags_[FromChargeKaonMuon], flags_[Muon] && flags_[ChargeKaonDecay]);
00467 update(flags_[DecayOnFlightMuon], flags_[FromChargePionMuon] || flags_[FromChargeKaonMuon]);
00468 }
00469
00470
00471 void TrackClassifier::vertexInformation()
00472 {
00473
00474 GeneratedPrimaryVertex const & genpv = genpvs_.back();
00475
00476
00477 TrackHistory::GenParticleTrail & genParticleTrail = const_cast<TrackHistory::GenParticleTrail &> (tracer_.genParticleTrail());
00478
00479
00480 int counter = 0;
00481
00482
00483 double const mm = 0.1;
00484
00485 double oldX = genpv.x;
00486 double oldY = genpv.y;
00487 double oldZ = genpv.z;
00488
00489
00490 for (
00491 TrackHistory::GenParticleTrail::reverse_iterator iparticle = genParticleTrail.rbegin();
00492 iparticle != genParticleTrail.rend();
00493 ++iparticle
00494 )
00495 {
00496
00497 HepMC::GenVertex * parent = (*iparticle)->production_vertex();
00498 if (parent)
00499 {
00500 HepMC::ThreeVector p = parent->point3d();
00501
00502 double distance2 = pow(p.x() * mm - genpv.x, 2) + pow(p.y() * mm - genpv.y, 2) + pow(p.z() * mm - genpv.z, 2);
00503 double difference2 = pow(p.x() * mm - oldX, 2) + pow(p.y() * mm - oldY, 2) + pow(p.z() * mm - oldZ, 2);
00504
00505
00506
00507
00508 if ( difference2 > vertexClusteringSqDistance_ )
00509 {
00510 if ( distance2 > vertexClusteringSqDistance_ ) counter++;
00511 oldX = p.x() * mm;
00512 oldY = p.y() * mm;
00513 oldZ = p.z() * mm;
00514 }
00515 }
00516 }
00517
00518 TrackHistory::SimParticleTrail & simParticleTrail = const_cast<TrackHistory::SimParticleTrail &> (tracer_.simParticleTrail());
00519
00520
00521 for (
00522 TrackHistory::SimParticleTrail::reverse_iterator iparticle = simParticleTrail.rbegin();
00523 iparticle != simParticleTrail.rend();
00524 ++iparticle
00525 )
00526 {
00527
00528 TrackingParticle::Point p = (*iparticle)->vertex();
00529
00530 double distance2 = pow(p.x() - genpv.x, 2) + pow(p.y() - genpv.y, 2) + pow(p.z() - genpv.z, 2);
00531 double difference2 = pow(p.x() - oldX, 2) + pow(p.y() - oldY, 2) + pow(p.z() - oldZ, 2);
00532
00533
00534
00535
00536 if ( difference2 > vertexClusteringSqDistance_ )
00537 {
00538 if ( distance2 > vertexClusteringSqDistance_ ) counter++;
00539 oldX = p.x();
00540 oldY = p.y();
00541 oldZ = p.z();
00542 }
00543 }
00544
00545 if ( !counter )
00546 flags_[PrimaryVertex] = true;
00547 else if ( counter == 1 )
00548 flags_[SecondaryVertex] = true;
00549 else
00550 flags_[TertiaryVertex] = true;
00551 }
00552
00553
00554 bool TrackClassifier::isFinalstateParticle(const HepMC::GenParticle * p)
00555 {
00556 return !p->end_vertex() && p->status() == 1;
00557 }
00558
00559
00560 bool TrackClassifier::isCharged(const HepMC::GenParticle * p)
00561 {
00562 const ParticleData * part = particleDataTable_->particle( p->pdg_id() );
00563 if (part)
00564 return part->charge()!=0;
00565 else
00566 {
00567
00568 return particleDataTable_->particle( -p->pdg_id() ) != 0;
00569 }
00570 }
00571
00572
00573 void TrackClassifier::genPrimaryVertices()
00574 {
00575 genpvs_.clear();
00576
00577 const HepMC::GenEvent * event = mcInformation_->GetEvent();
00578
00579 if (event)
00580 {
00581 int idx = 0;
00582
00583
00584 for ( HepMC::GenEvent::vertex_const_iterator ivertex = event->vertices_begin(); ivertex != event->vertices_end(); ++ivertex )
00585 {
00586 bool hasParentVertex = false;
00587
00588
00589 for (
00590 HepMC::GenVertex::particle_iterator iparent = (*ivertex)->particles_begin(HepMC::parents);
00591 iparent != (*ivertex)->particles_end(HepMC::parents);
00592 ++iparent
00593 )
00594 if ( (*iparent)->production_vertex() )
00595 {
00596 hasParentVertex = true;
00597 break;
00598 }
00599
00600
00601 if (hasParentVertex) continue;
00602
00603
00604 HepMC::FourVector pos = (*ivertex)->position();
00605
00606 double const mm = 0.1;
00607
00608 GeneratedPrimaryVertex pv(pos.x()*mm, pos.y()*mm, pos.z()*mm);
00609
00610 std::vector<GeneratedPrimaryVertex>::iterator ientry = genpvs_.begin();
00611
00612
00613 for (; ientry != genpvs_.end(); ++ientry)
00614 {
00615 double distance2 = pow(pv.x - ientry->x, 2) + pow(pv.y - ientry->y, 2) + pow(pv.z - ientry->z, 2);
00616 if ( distance2 < vertexClusteringSqDistance_ )
00617 break;
00618 }
00619
00620
00621 if (ientry == genpvs_.end())
00622 ientry = genpvs_.insert(ientry,pv);
00623
00624
00625 ientry->genVertex.push_back((*ivertex)->barcode());
00626
00627
00628 for (
00629 HepMC::GenVertex::particle_iterator idecendants = (*ivertex)->particles_begin(HepMC::descendants);
00630 idecendants != (*ivertex)->particles_end(HepMC::descendants);
00631 ++idecendants
00632 )
00633 {
00634 if (isFinalstateParticle(*idecendants))
00635 if ( find(ientry->finalstateParticles.begin(), ientry->finalstateParticles.end(), (*idecendants)->barcode()) == ientry->finalstateParticles.end() )
00636 {
00637 ientry->finalstateParticles.push_back((*idecendants)->barcode());
00638 HepMC::FourVector m = (*idecendants)->momentum();
00639
00640 ientry->ptot.setPx(ientry->ptot.px() + m.px());
00641 ientry->ptot.setPy(ientry->ptot.py() + m.py());
00642 ientry->ptot.setPz(ientry->ptot.pz() + m.pz());
00643 ientry->ptot.setE(ientry->ptot.e() + m.e());
00644 ientry->ptsq += m.perp() * m.perp();
00645
00646 if ( m.perp() > 0.8 && std::abs(m.pseudoRapidity()) < 2.5 && isCharged(*idecendants) ) ientry->nGenTrk++;
00647 }
00648 }
00649 idx++;
00650 }
00651 }
00652
00653 std::sort(genpvs_.begin(), genpvs_.end());
00654 }