TrackClassifier.cc

Go to the documentation of this file.

#include <cmath>
#include <cstdlib>
#include <iostream>

#include "HepPDT/ParticleID.hh"

#include "SimTracker/TrackHistory/interface/TrackClassifier.h"

#define update(a, b) \
  do {               \
    (a) = (a) | (b); \
  } while (0)

TrackClassifier::TrackClassifier(edm::ParameterSet const &config, edm::ConsumesCollector &&collector)
    : TrackCategories(),
      hepMCLabel_(config.getUntrackedParameter<edm::InputTag>("hepMC")),
      beamSpotLabel_(config.getUntrackedParameter<edm::InputTag>("beamSpot")),
      tracer_(config, std::move(collector)),
      quality_(config, collector) {
  collector.consumes<edm::HepMCProduct>(hepMCLabel_);
  collector.consumes<reco::BeamSpot>(beamSpotLabel_);

  // Set the history depth after hadronization
  tracer_.depth(-2);

  // Set the maximum d0pull for the bad category
  badPull_ = config.getUntrackedParameter<double>("badPull");

  // Set the minimum decay length for detecting long decays
  longLivedDecayLength_ = config.getUntrackedParameter<double>("longLivedDecayLength");

  // Set the distance for clustering vertices
  float vertexClusteringDistance = config.getUntrackedParameter<double>("vertexClusteringDistance");
  vertexClusteringSqDistance_ = vertexClusteringDistance * vertexClusteringDistance;

  // Set the number of innermost layers to check for bad hits
  numberOfInnerLayers_ = config.getUntrackedParameter<unsigned int>("numberOfInnerLayers");

  // Set the minimum number of simhits in the tracker
  minTrackerSimHits_ = config.getUntrackedParameter<unsigned int>("minTrackerSimHits");
}

void TrackClassifier::newEvent(edm::Event const &event, edm::EventSetup const &setup) {
  // Get the new event information for the tracer
  tracer_.newEvent(event, setup);

  // Get the new event information for the track quality analyser
  quality_.newEvent(event, setup);

  // Get hepmc of the event
  event.getByLabel(hepMCLabel_, mcInformation_);

  // Magnetic field
  setup.get<IdealMagneticFieldRecord>().get(magneticField_);

  // Get the partivle data table
  setup.getData(particleDataTable_);

  // get the beam spot
  event.getByLabel(beamSpotLabel_, beamSpot_);

  // Transient track builder
  setup.get<TransientTrackRecord>().get("TransientTrackBuilder", transientTrackBuilder_);

  // Create the list of primary vertices associated to the event
  genPrimaryVertices();

  // Retrieve tracker topology from geometry
  edm::ESHandle<TrackerTopology> tTopoHand;
  setup.get<TrackerTopologyRcd>().get(tTopoHand);
  tTopo_ = tTopoHand.product();
}

TrackClassifier const &TrackClassifier::evaluate(reco::TrackBaseRef const &track) {
  // Initializing the category vector
  reset();

  // Associate and evaluate the track history (check for fakes)
  if (tracer_.evaluate(track)) {
    // Classify all the tracks by their association and reconstruction
    // information
    reconstructionInformation(track);

    // Get all the information related to the simulation details
    simulationInformation();

    // Analyse the track reconstruction quality
    qualityInformation(track);

    // Get hadron flavor of the initial hadron
    hadronFlavor();

    // Get all the information related to decay process
    processesAtGenerator();

    // Get information about conversion and other interactions
    processesAtSimulation();

    // Get geometrical information about the vertices
    vertexInformation();

    // Check for unkown classification
    unknownTrack();
  } else
    flags_[Fake] = true;

  return *this;
}

TrackClassifier const &TrackClassifier::evaluate(TrackingParticleRef const &track) {
  // Initializing the category vector
  reset();

  // Trace the history for the given TP
  tracer_.evaluate(track);

  // Collect the associated reco track
  const reco::TrackBaseRef &recotrack = tracer_.recoTrack();

  // If there is a reco truck then evaluate the simulated history
  if (recotrack.isNonnull()) {
    flags_[Reconstructed] = true;
    // Classify all the tracks by their association and reconstruction
    // information
    reconstructionInformation(recotrack);
    // Analyse the track reconstruction quality
    qualityInformation(recotrack);
  } else
    flags_[Reconstructed] = false;

  // Get all the information related to the simulation details
  simulationInformation();

  // Get hadron flavor of the initial hadron
  hadronFlavor();

  // Get all the information related to decay process
  processesAtGenerator();

  // Get information about conversion and other interactions
  processesAtSimulation();

  // Get geometrical information about the vertices
  vertexInformation();

  // Check for unkown classification
  unknownTrack();

  return *this;
}

void TrackClassifier::reconstructionInformation(reco::TrackBaseRef const &track) {
  TrackingParticleRef tpr = tracer_.simParticle();

  // Compute tracking particle parameters at point of closest approach to the
  // beamline

  const SimTrack *assocTrack = &(*tpr->g4Track_begin());

  FreeTrajectoryState ftsAtProduction(
      GlobalPoint(tpr->vertex().x(), tpr->vertex().y(), tpr->vertex().z()),
      GlobalVector(assocTrack->momentum().x(), assocTrack->momentum().y(), assocTrack->momentum().z()),
      TrackCharge(track->charge()),
      magneticField_.product());

  try {
    TSCPBuilderNoMaterial tscpBuilder;
    TrajectoryStateClosestToPoint tsAtClosestApproach =
        tscpBuilder(ftsAtProduction, GlobalPoint(beamSpot_->x0(), beamSpot_->y0(), beamSpot_->z0()));

    GlobalVector v =
        tsAtClosestApproach.theState().position() - GlobalPoint(beamSpot_->x0(), beamSpot_->y0(), beamSpot_->z0());
    GlobalVector p = tsAtClosestApproach.theState().momentum();

    // Simulated dxy
    double dxySim = -v.x() * sin(p.phi()) + v.y() * cos(p.phi());

    // Simulated dz
    double dzSim = v.z() - (v.x() * p.x() + v.y() * p.y()) * p.z() / p.perp2();

    // Calculate the dxy pull
    double dxyPull =
        std::abs(track->dxy(reco::TrackBase::Point(beamSpot_->x0(), beamSpot_->y0(), beamSpot_->z0())) - dxySim) /
        track->dxyError();

    // Calculate the dx pull
    double dzPull =
        std::abs(track->dz(reco::TrackBase::Point(beamSpot_->x0(), beamSpot_->y0(), beamSpot_->z0())) - dzSim) /
        track->dzError();

    // Return true if d0Pull > badD0Pull sigmas
    flags_[Bad] = (dxyPull > badPull_ || dzPull > badPull_);

  } catch (cms::Exception const &) {
    flags_[Bad] = true;
  }
}

void TrackClassifier::simulationInformation() {
  // Get the event id for the initial TP.
  EncodedEventId eventId = tracer_.simParticle()->eventId();
  // Check for signal events
  flags_[SignalEvent] = !eventId.bunchCrossing() && !eventId.event();
  // Check for muons
  flags_[Muon] = (abs(tracer_.simParticle()->pdgId()) == 13);
  // Check for the number of psimhit in tracker
  flags_[TrackerSimHits] = tracer_.simParticle()->numberOfTrackerLayers() >= (int)minTrackerSimHits_;
}

void TrackClassifier::qualityInformation(reco::TrackBaseRef const &track) {
  // run the hit-by-hit reconstruction quality analysis
  quality_.evaluate(tracer_.simParticleTrail(), track, tTopo_);

  unsigned int maxLayers = std::min(numberOfInnerLayers_, quality_.numberOfLayers());

  // check the innermost layers for bad hits
  for (unsigned int i = 0; i < maxLayers; i++) {
    const TrackQuality::Layer &layer = quality_.layer(i);

    // check all hits in that layer
    for (unsigned int j = 0; j < layer.hits.size(); j++) {
      const TrackQuality::Layer::Hit &hit = layer.hits[j];

      // In those cases the bad hit was used by track reconstruction
      if (hit.state == TrackQuality::Layer::Noise || hit.state == TrackQuality::Layer::Misassoc)
        flags_[BadInnerHits] = true;
      else if (hit.state == TrackQuality::Layer::Shared)
        flags_[SharedInnerHits] = true;
    }
  }
}

void TrackClassifier::hadronFlavor() {
  // Get the initial hadron from the recoGenParticleTrail
  const reco::GenParticle *particle = tracer_.recoGenParticle();

  // Check for the initial hadron
  if (particle) {
    HepPDT::ParticleID pid(particle->pdgId());
    flags_[Bottom] = pid.hasBottom();
    flags_[Charm] = pid.hasCharm();
    flags_[Light] = !pid.hasCharm() && !pid.hasBottom();
  }
}

void TrackClassifier::processesAtGenerator() {
  // pdgid of the "in" particle to the production vertex
  int pdgid = 0;

  // Get the generated particles from track history (reco::GenParticle in the
  // recoGenParticleTrail)
  TrackHistory::RecoGenParticleTrail const &recoGenParticleTrail = tracer_.recoGenParticleTrail();

  // Loop over the generated particles (reco::GenParticle in the
  // recoGenParticleTrail)
  for (TrackHistory::RecoGenParticleTrail::const_iterator iparticle = recoGenParticleTrail.begin();
       iparticle != recoGenParticleTrail.end();
       ++iparticle) {
    pdgid = std::abs((*iparticle)->pdgId());
    // Get particle type
    HepPDT::ParticleID particleID(pdgid);

    // Check if the particle type is valid one
    if (particleID.isValid()) {
      // Get particle data
      ParticleData const *particleData = particleDataTable_->particle(particleID);
      // Check if the particle exist in the table
      if (particleData) {
        // Check if their life time is bigger than longLivedDecayLength_
        if (particleData->lifetime() > longLivedDecayLength_)
          update(flags_[LongLivedDecay], true);
        // Check for B and C weak decays
        update(flags_[BWeakDecay], particleID.hasBottom());
        update(flags_[CWeakDecay], particleID.hasCharm());
        // Check for B and C pure leptonic decay
        std::set<int> daughterIds;
        size_t ndau = (*iparticle)->numberOfDaughters();
        for (size_t i = 0; i < ndau; ++i) {
          daughterIds.insert((*iparticle)->daughter(i)->pdgId());
        }
        update(flags_[FromBWeakDecayMuon], particleID.hasBottom() && (daughterIds.find(13) != daughterIds.end()));
        update(flags_[FromCWeakDecayMuon], particleID.hasCharm() && (daughterIds.find(13) != daughterIds.end()));
      }
      // Check Tau, Ks and Lambda decay
      update(flags_[ChargePionDecay], pdgid == 211);
      update(flags_[ChargeKaonDecay], pdgid == 321);
      update(flags_[TauDecay], pdgid == 15);
      update(flags_[KsDecay], pdgid == 310);
      update(flags_[LambdaDecay], pdgid == 3122);
      update(flags_[JpsiDecay], pdgid == 443);
      update(flags_[XiDecay], pdgid == 3312);
      update(flags_[SigmaPlusDecay], pdgid == 3222);
      update(flags_[SigmaMinusDecay], pdgid == 3112);
    }
  }
  // Decays in flight
  update(flags_[FromChargePionMuon], flags_[Muon] && flags_[ChargePionDecay]);
  update(flags_[FromChargeKaonMuon], flags_[Muon] && flags_[ChargeKaonDecay]);
  update(flags_[DecayOnFlightMuon], (flags_[FromChargePionMuon] || flags_[FromChargeKaonMuon]));
}

void TrackClassifier::processesAtSimulation() {
  TrackHistory::SimParticleTrail const &simParticleTrail = tracer_.simParticleTrail();

  // Loop over the simulated particles
  for (TrackHistory::SimParticleTrail::const_iterator iparticle = simParticleTrail.begin();
       iparticle != simParticleTrail.end();
       ++iparticle) {
    // pdgid of the real source parent vertex
    int pdgid = 0;

    // Get a reference to the TP's parent vertex
    TrackingVertexRef const &parentVertex = (*iparticle)->parentVertex();

    // Look for the original source track
    if (parentVertex.isNonnull()) {
      // select the original source in case of combined vertices
      bool flag = false;
      TrackingVertex::tp_iterator itd, its;

      for (its = parentVertex->sourceTracks_begin(); its != parentVertex->sourceTracks_end(); ++its) {
        for (itd = parentVertex->daughterTracks_begin(); itd != parentVertex->daughterTracks_end(); ++itd)
          if (itd != its) {
            flag = true;
            break;
          }
        if (flag)
          break;
      }

      // Collect the pdgid of the original source track
      if (its != parentVertex->sourceTracks_end())
        pdgid = std::abs((*its)->pdgId());
      else
        pdgid = 0;
    }

    unsigned int processG4 = 0;

    // Check existence of SimVerteces assigned
    if (parentVertex->nG4Vertices() > 0) {
      processG4 = (*(parentVertex->g4Vertices_begin())).processType();
    }

    unsigned int process = g4toCMSProcMap_.processId(processG4);

    // Flagging all the different processes
    update(flags_[KnownProcess], process != CMS::Undefined && process != CMS::Unknown && process != CMS::Primary);

    update(flags_[UndefinedProcess], process == CMS::Undefined);
    update(flags_[UnknownProcess], process == CMS::Unknown);
    update(flags_[PrimaryProcess], process == CMS::Primary);
    update(flags_[HadronicProcess], process == CMS::Hadronic);
    update(flags_[DecayProcess], process == CMS::Decay);
    update(flags_[ComptonProcess], process == CMS::Compton);
    update(flags_[AnnihilationProcess], process == CMS::Annihilation);
    update(flags_[EIoniProcess], process == CMS::EIoni);
    update(flags_[HIoniProcess], process == CMS::HIoni);
    update(flags_[MuIoniProcess], process == CMS::MuIoni);
    update(flags_[PhotonProcess], process == CMS::Photon);
    update(flags_[MuPairProdProcess], process == CMS::MuPairProd);
    update(flags_[ConversionsProcess], process == CMS::Conversions);
    update(flags_[EBremProcess], process == CMS::EBrem);
    update(flags_[SynchrotronRadiationProcess], process == CMS::SynchrotronRadiation);
    update(flags_[MuBremProcess], process == CMS::MuBrem);
    update(flags_[MuNuclProcess], process == CMS::MuNucl);

    // Get particle type
    HepPDT::ParticleID particleID(pdgid);

    // Check if the particle type is valid one
    if (particleID.isValid()) {
      // Get particle data
      ParticleData const *particleData = particleDataTable_->particle(particleID);
      // Special treatment for decays
      if (process == CMS::Decay) {
        // Check if the particle exist in the table
        if (particleData) {
          // Check if their life time is bigger than 1e-14
          if (particleDataTable_->particle(particleID)->lifetime() > longLivedDecayLength_)
            update(flags_[LongLivedDecay], true);

          // Check for B and C weak decays
          update(flags_[BWeakDecay], particleID.hasBottom());
          update(flags_[CWeakDecay], particleID.hasCharm());

          // Check for B or C pure leptonic decays
          int daughtId = abs((*iparticle)->pdgId());
          update(flags_[FromBWeakDecayMuon], particleID.hasBottom() && daughtId == 13);
          update(flags_[FromCWeakDecayMuon], particleID.hasCharm() && daughtId == 13);
        }
        // Check decays
        update(flags_[ChargePionDecay], pdgid == 211);
        update(flags_[ChargeKaonDecay], pdgid == 321);
        update(flags_[TauDecay], pdgid == 15);
        update(flags_[KsDecay], pdgid == 310);
        update(flags_[LambdaDecay], pdgid == 3122);
        update(flags_[JpsiDecay], pdgid == 443);
        update(flags_[XiDecay], pdgid == 3312);
        update(flags_[OmegaDecay], pdgid == 3334);
        update(flags_[SigmaPlusDecay], pdgid == 3222);
        update(flags_[SigmaMinusDecay], pdgid == 3112);
      }
    }
  }
  // Decays in flight
  update(flags_[FromChargePionMuon], flags_[Muon] && flags_[ChargePionDecay]);
  update(flags_[FromChargeKaonMuon], flags_[Muon] && flags_[ChargeKaonDecay]);
  update(flags_[DecayOnFlightMuon], flags_[FromChargePionMuon] || flags_[FromChargeKaonMuon]);
}

void TrackClassifier::vertexInformation() {
  // Get the main primary vertex from the list
  GeneratedPrimaryVertex const &genpv = genpvs_.back();

  // Get the generated history of the tracks
  const TrackHistory::GenParticleTrail &genParticleTrail = tracer_.genParticleTrail();

  // Vertex counter
  int counter = 0;

  // Unit transformation from mm to cm
  double const mm = 0.1;

  double oldX = genpv.x;
  double oldY = genpv.y;
  double oldZ = genpv.z;

  // Loop over the generated particles
  for (TrackHistory::GenParticleTrail::const_reverse_iterator iparticle = genParticleTrail.rbegin();
       iparticle != genParticleTrail.rend();
       ++iparticle) {
    // Look for those with production vertex
    HepMC::GenVertex *parent = (*iparticle)->production_vertex();
    if (parent) {
      HepMC::ThreeVector p = parent->point3d();

      double distance2 = pow(p.x() * mm - genpv.x, 2) + pow(p.y() * mm - genpv.y, 2) + pow(p.z() * mm - genpv.z, 2);
      double difference2 = pow(p.x() * mm - oldX, 2) + pow(p.y() * mm - oldY, 2) + pow(p.z() * mm - oldZ, 2);

      // std::cout << "Distance2 : " << distance2 << " (" << p.x() * mm << ","
      // << p.y() * mm << "," << p.z() * mm << ")" << std::endl; std::cout <<
      // "Difference2 : " << difference2 << std::endl;

      if (difference2 > vertexClusteringSqDistance_) {
        if (distance2 > vertexClusteringSqDistance_)
          counter++;
        oldX = p.x() * mm;
        oldY = p.y() * mm;
        oldZ = p.z() * mm;
      }
    }
  }

  const TrackHistory::SimParticleTrail &simParticleTrail = tracer_.simParticleTrail();

  // Loop over the generated particles
  for (TrackHistory::SimParticleTrail::const_reverse_iterator iparticle = simParticleTrail.rbegin();
       iparticle != simParticleTrail.rend();
       ++iparticle) {
    // Look for those with production vertex
    TrackingParticle::Point p = (*iparticle)->vertex();

    double distance2 = pow(p.x() - genpv.x, 2) + pow(p.y() - genpv.y, 2) + pow(p.z() - genpv.z, 2);
    double difference2 = pow(p.x() - oldX, 2) + pow(p.y() - oldY, 2) + pow(p.z() - oldZ, 2);

    // std::cout << "Distance2 : " << distance2 << " (" << p.x() << "," << p.y()
    // << "," << p.z() << ")" << std::endl; std::cout << "Difference2 : " <<
    // difference2 << std::endl;

    if (difference2 > vertexClusteringSqDistance_) {
      if (distance2 > vertexClusteringSqDistance_)
        counter++;
      oldX = p.x();
      oldY = p.y();
      oldZ = p.z();
    }
  }

  if (!counter)
    flags_[PrimaryVertex] = true;
  else if (counter == 1)
    flags_[SecondaryVertex] = true;
  else
    flags_[TertiaryVertex] = true;
}

bool TrackClassifier::isFinalstateParticle(const HepMC::GenParticle *p) { return !p->end_vertex() && p->status() == 1; }

bool TrackClassifier::isCharged(const HepMC::GenParticle *p) {
  const ParticleData *part = particleDataTable_->particle(p->pdg_id());
  if (part)
    return part->charge() != 0;
  else {
    // the new/improved particle table doesn't know anti-particles
    return particleDataTable_->particle(-p->pdg_id()) != nullptr;
  }
}

void TrackClassifier::genPrimaryVertices() {
  genpvs_.clear();

  const HepMC::GenEvent *event = mcInformation_->GetEvent();

  if (event) {
    int idx = 0;

    // Loop over the different GenVertex
    for (HepMC::GenEvent::vertex_const_iterator ivertex = event->vertices_begin(); ivertex != event->vertices_end();
         ++ivertex) {
      bool hasParentVertex = false;

      // Loop over the parents looking to see if they are coming from a
      // production vertex
      for (HepMC::GenVertex::particle_iterator iparent = (*ivertex)->particles_begin(HepMC::parents);
           iparent != (*ivertex)->particles_end(HepMC::parents);
           ++iparent)
        if ((*iparent)->production_vertex()) {
          hasParentVertex = true;
          break;
        }

      // Reject those vertices with parent vertices
      if (hasParentVertex)
        continue;

      // Get the position of the vertex
      HepMC::FourVector pos = (*ivertex)->position();

      double const mm = 0.1;

      GeneratedPrimaryVertex pv(pos.x() * mm, pos.y() * mm, pos.z() * mm);

      std::vector<GeneratedPrimaryVertex>::iterator ientry = genpvs_.begin();

      // Search for a VERY close vertex in the list
      for (; ientry != genpvs_.end(); ++ientry) {
        double distance2 = pow(pv.x - ientry->x, 2) + pow(pv.y - ientry->y, 2) + pow(pv.z - ientry->z, 2);
        if (distance2 < vertexClusteringSqDistance_)
          break;
      }

      // Check if there is not a VERY close vertex and added to the list
      if (ientry == genpvs_.end())
        ientry = genpvs_.insert(ientry, pv);

      // Add the vertex barcodes to the new or existent vertices
      ientry->genVertex.push_back((*ivertex)->barcode());

      // Collect final state descendants
      for (HepMC::GenVertex::particle_iterator idecendants = (*ivertex)->particles_begin(HepMC::descendants);
           idecendants != (*ivertex)->particles_end(HepMC::descendants);
           ++idecendants) {
        if (isFinalstateParticle(*idecendants))
          if (find(ientry->finalstateParticles.begin(), ientry->finalstateParticles.end(), (*idecendants)->barcode()) ==
              ientry->finalstateParticles.end()) {
            ientry->finalstateParticles.push_back((*idecendants)->barcode());
            HepMC::FourVector m = (*idecendants)->momentum();

            ientry->ptot.setPx(ientry->ptot.px() + m.px());
            ientry->ptot.setPy(ientry->ptot.py() + m.py());
            ientry->ptot.setPz(ientry->ptot.pz() + m.pz());
            ientry->ptot.setE(ientry->ptot.e() + m.e());
            ientry->ptsq += m.perp() * m.perp();

            if (m.perp() > 0.8 && std::abs(m.pseudoRapidity()) < 2.5 && isCharged(*idecendants))
              ientry->nGenTrk++;
          }
      }
      idx++;
    }
  }

  std::sort(genpvs_.begin(), genpvs_.end());
}