TrackingTruthAccumulator.cc

Go to the documentation of this file.

 
#include "SimGeneral/TrackingAnalysis/plugins/TrackingTruthAccumulator.h"
 
#include <memory>
 
#include "FWCore/Framework/interface/ConsumesCollector.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h"
#include "SimDataFormats/TrackingHit/interface/PSimHit.h"
#include "SimGeneral/MixingModule/interface/DigiAccumulatorMixModFactory.h"
#include "SimGeneral/MixingModule/interface/PileUpEventPrincipal.h"
#include "SimGeneral/TrackingAnalysis/interface/EncodedTruthId.h"
 
#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
#include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Utilities/interface/isFinite.h"
#include "Geometry/Records/interface/TrackerTopologyRcd.h"
 
// Turn on integrity checking
//#define DO_DEBUG_TESTING
 
//---------------------------------------------------------------------------------
//---------------------------------------------------------------------------------
//---------------------------------------------------------------------------------
//------ ------
//------  Declarations for utility classes and functions in the unnamed ------
//------  namespace. The definitions are right at the bottom of the file. ------
//------ ------
//---------------------------------------------------------------------------------
//---------------------------------------------------------------------------------
//---------------------------------------------------------------------------------
 
namespace {
  struct DecayChainTrack {
    int simTrackIndex;
    struct DecayChainVertex *pParentVertex;
    // Some things have multiple decay vertices. Not sure how that works but it
    // seems to be mostly electrons and some photons.
    std::vector<struct DecayChainVertex *> daughterVertices;
    DecayChainTrack *pMergedBremSource;
    DecayChainTrack() : simTrackIndex(-1), pParentVertex(nullptr), pMergedBremSource(nullptr) {}
    DecayChainTrack(int newSimTrackIndex)
        : simTrackIndex(newSimTrackIndex), pParentVertex(nullptr), pMergedBremSource() {}
  };
 
  struct DecayChainVertex {
    int simVertexIndex;
    DecayChainTrack *pParentTrack;
    std::vector<DecayChainTrack *> daughterTracks;
    DecayChainVertex *pMergedBremSource;
    DecayChainVertex() : simVertexIndex(-1), pParentTrack(nullptr), pMergedBremSource(nullptr) {}
    DecayChainVertex(int newIndex) : simVertexIndex(newIndex), pParentTrack(nullptr), pMergedBremSource(nullptr) {}
  };
 
  struct DecayChain {
  public:
    DecayChain(const std::vector<SimTrack> &trackCollection, const std::vector<SimVertex> &vertexCollection);
    const size_t decayTracksSize;
    const size_t decayVerticesSize;
 
#if defined(DO_DEBUG_TESTING)
 
    void integrityCheck();
#endif
    const SimTrack &getSimTrack(const ::DecayChainTrack *pDecayTrack) const {
      return simTrackCollection_.at(pDecayTrack->simTrackIndex);
    }
    const SimVertex &getSimVertex(const ::DecayChainVertex *pDecayVertex) const {
      return simVertexCollection_.at(pDecayVertex->simVertexIndex);
    }
 
  private:
    void findBrem(const std::vector<SimTrack> &trackCollection, const std::vector<SimVertex> &vertexCollection);
    std::unique_ptr<::DecayChainTrack[]> decayTracks_;
    std::unique_ptr<::DecayChainVertex[]> decayVertices_;
 
    std::vector<::DecayChainVertex *> rootVertices_;
 
    // Keep a record of the constructor parameters
    const std::vector<SimTrack> &simTrackCollection_;
    const std::vector<SimVertex> &simVertexCollection_;
 
  public:
    const std::unique_ptr<::DecayChainTrack[]> &decayTracks;  
    const std::unique_ptr<::DecayChainVertex[]> &decayVertices;  
    const std::vector<::DecayChainVertex *> &rootVertices;       
  };
 
  class TrackingParticleFactory {
  public:
    TrackingParticleFactory(const ::DecayChain &decayChain,
                            const edm::Handle<std::vector<reco::GenParticle>> &hGenParticles,
                            const edm::Handle<edm::HepMCProduct> &hepMCproduct,
                            const edm::Handle<std::vector<int>> &hHepMCGenParticleIndices,
                            const std::vector<const PSimHit *> &simHits,
                            double volumeRadius,
                            double volumeZ,
                            double vertexDistanceCut,
                            bool allowDifferentProcessTypes);
    TrackingParticle createTrackingParticle(const DecayChainTrack *pTrack, const TrackerTopology *tTopo) const;
    TrackingVertex createTrackingVertex(const DecayChainVertex *pVertex) const;
    bool vectorIsInsideVolume(const math::XYZTLorentzVectorD &vector) const;
 
  private:
    const ::DecayChain &decayChain_;
    const edm::Handle<std::vector<reco::GenParticle>> &hGenParticles_;
    const edm::Handle<edm::HepMCProduct> &hepMCproduct_;
    std::vector<int> genParticleIndices_;
    const std::vector<const PSimHit *> &simHits_;
    const double volumeRadius_;
    const double volumeZ_;
    const double vertexDistanceCut2_;                        // distance based on which HepMC::GenVertexs
                                                             // are added to SimVertexs
    std::multimap<unsigned int, size_t> trackIdToHitIndex_;  
    bool allowDifferentProcessTypeForDifferentDetectors_;    
  };
 
  class OutputCollectionWrapper {
  public:
    OutputCollectionWrapper(const DecayChain &decayChain,
                            TrackingTruthAccumulator::OutputCollections &outputCollections);
    TrackingParticle *addTrackingParticle(const ::DecayChainTrack *pDecayTrack,
                                          const TrackingParticle &trackingParticle);
    TrackingVertex *addTrackingVertex(const ::DecayChainVertex *pDecayVertex, const TrackingVertex &trackingVertex);
    TrackingParticle *getTrackingParticle(const ::DecayChainTrack *pDecayTrack);
    void setProxy(const ::DecayChainTrack *pOriginalTrack, const ::DecayChainTrack *pProxyTrack);
    void setProxy(const ::DecayChainVertex *pOriginalVertex, const ::DecayChainVertex *pProxyVertex);
    TrackingVertex *getTrackingVertex(const ::DecayChainVertex *pDecayVertex);
    TrackingParticleRef getRef(const ::DecayChainTrack *pDecayTrack);
    TrackingVertexRef getRef(const ::DecayChainVertex *pDecayVertex);
 
  private:
    void associateToExistingObjects(const ::DecayChainVertex *pChainVertex);
    void associateToExistingObjects(const ::DecayChainTrack *pChainTrack);
    TrackingTruthAccumulator::OutputCollections &output_;
    std::vector<int> trackingParticleIndices_;
    std::vector<int> trackingVertexIndices_;
  };
 
  TrackingParticle *addTrackAndParentVertex(::DecayChainTrack *pDecayTrack,
                                            const TrackingParticle &trackingParticle,
                                            ::OutputCollectionWrapper *pOutput);
 
  void addTrack(::DecayChainTrack *pDecayChainTrack,
                const TrackingParticleSelector *pSelector,
                ::OutputCollectionWrapper *pUnmergedOutput,
                ::OutputCollectionWrapper *pMergedOutput,
                const ::TrackingParticleFactory &objectFactory,
                bool addAncestors,
                const TrackerTopology *tTopo);
 
}  // namespace
 
//---------------------------------------------------------------------------------
//---------------------------------------------------------------------------------
//---------------------------------------------------------------------------------
//------ ------
//------  Definitions for the TrackingTruthAccumulator methods ------
//------  are below. ------
//------ ------
//---------------------------------------------------------------------------------
//---------------------------------------------------------------------------------
//---------------------------------------------------------------------------------
 
TrackingTruthAccumulator::TrackingTruthAccumulator(const edm::ParameterSet &config,
                                                   edm::ProducesCollector producesCollector,
                                                   edm::ConsumesCollector &iC)
    : messageCategory_("TrackingTruthAccumulator"),
      volumeRadius_(config.getParameter<double>("volumeRadius")),
      volumeZ_(config.getParameter<double>("volumeZ")),
      vertexDistanceCut_(config.getParameter<double>("vertexDistanceCut")),
      ignoreTracksOutsideVolume_(config.getParameter<bool>("ignoreTracksOutsideVolume")),
      maximumPreviousBunchCrossing_(config.getParameter<unsigned int>("maximumPreviousBunchCrossing")),
      maximumSubsequentBunchCrossing_(config.getParameter<unsigned int>("maximumSubsequentBunchCrossing")),
      createUnmergedCollection_(config.getParameter<bool>("createUnmergedCollection")),
      createMergedCollection_(config.getParameter<bool>("createMergedBremsstrahlung")),
      createInitialVertexCollection_(config.getParameter<bool>("createInitialVertexCollection")),
      addAncestors_(config.getParameter<bool>("alwaysAddAncestors")),
      removeDeadModules_(config.getParameter<bool>("removeDeadModules")),
      simTrackLabel_(config.getParameter<edm::InputTag>("simTrackCollection")),
      simVertexLabel_(config.getParameter<edm::InputTag>("simVertexCollection")),
      collectionTags_(),
      genParticleLabel_(config.getParameter<edm::InputTag>("genParticleCollection")),
      hepMCproductLabel_(config.getParameter<edm::InputTag>("HepMCProductLabel")),
      allowDifferentProcessTypeForDifferentDetectors_(config.getParameter<bool>("allowDifferentSimHitProcesses")) {
  //
  // Make sure at least one of the merged and unmerged collections have been set
  // to be created.
  //
  if (!createUnmergedCollection_ && !createMergedCollection_)
    edm::LogError(messageCategory_) << "Both \"createUnmergedCollection\" and "
                                       "\"createMergedBremsstrahlung\" have been"
                                    << "set to false, which means no collections will be created";
 
  // Initialize selection for building TrackingParticles
  //
  if (config.exists("select")) {
    edm::ParameterSet param = config.getParameter<edm::ParameterSet>("select");
    selector_ = TrackingParticleSelector(param.getParameter<double>("ptMinTP"),
                                         param.getParameter<double>("ptMaxTP"),
                                         param.getParameter<double>("minRapidityTP"),
                                         param.getParameter<double>("maxRapidityTP"),
                                         param.getParameter<double>("tipTP"),
                                         param.getParameter<double>("lipTP"),
                                         param.getParameter<int>("minHitTP"),
                                         param.getParameter<bool>("signalOnlyTP"),
                                         param.getParameter<bool>("intimeOnlyTP"),
                                         param.getParameter<bool>("chargedOnlyTP"),
                                         param.getParameter<bool>("stableOnlyTP"),
                                         param.getParameter<std::vector<int>>("pdgIdTP"));
    selectorFlag_ = true;
 
    // Also set these two variables, which are used to drop out early if the
    // SimTrack doesn't conform. The selector_ requires a full TrackingParticle
    // object, but these two variables can veto things early.
    chargedOnly_ = param.getParameter<bool>("chargedOnlyTP");
    signalOnly_ = param.getParameter<bool>("signalOnlyTP");
  } else {
    selectorFlag_ = false;
    chargedOnly_ = false;
    signalOnly_ = false;
  }
 
  //
  // Need to state what collections are going to be added to the event. This
  // depends on which of the merged and unmerged collections have been
  // configured to be created.
  //
  if (createUnmergedCollection_) {
    producesCollector.produces<TrackingVertexCollection>();
    producesCollector.produces<TrackingParticleCollection>();
  }
 
  if (createMergedCollection_) {
    producesCollector.produces<TrackingParticleCollection>("MergedTrackTruth");
    producesCollector.produces<TrackingVertexCollection>("MergedTrackTruth");
  }
 
  if (createInitialVertexCollection_) {
    producesCollector.produces<TrackingVertexCollection>("InitialVertices");
  }
 
  iC.consumes<std::vector<SimTrack>>(simTrackLabel_);
  iC.consumes<std::vector<SimVertex>>(simVertexLabel_);
  iC.consumes<std::vector<reco::GenParticle>>(genParticleLabel_);
  iC.consumes<std::vector<int>>(genParticleLabel_);
  iC.consumes<std::vector<int>>(hepMCproductLabel_);
 
  // Fill the collection tags
  const edm::ParameterSet &simHitCollectionConfig = config.getParameterSet("simHitCollections");
  std::vector<std::string> parameterNames = simHitCollectionConfig.getParameterNames();
 
  for (const auto &parameterName : parameterNames) {
    std::vector<edm::InputTag> tags = simHitCollectionConfig.getParameter<std::vector<edm::InputTag>>(parameterName);
    collectionTags_.insert(collectionTags_.end(), tags.begin(), tags.end());
  }
 
  for (const auto &collectionTag : collectionTags_) {
    iC.consumes<std::vector<PSimHit>>(collectionTag);
  }
}
 
void TrackingTruthAccumulator::initializeEvent(edm::Event const &event, edm::EventSetup const &setup) {
  if (createUnmergedCollection_) {
    unmergedOutput_.pTrackingParticles = std::make_unique<TrackingParticleCollection>();
    unmergedOutput_.pTrackingVertices = std::make_unique<TrackingVertexCollection>();
    unmergedOutput_.refTrackingParticles =
        const_cast<edm::Event &>(event).getRefBeforePut<TrackingParticleCollection>();
    unmergedOutput_.refTrackingVertexes = const_cast<edm::Event &>(event).getRefBeforePut<TrackingVertexCollection>();
  }
 
  if (createMergedCollection_) {
    mergedOutput_.pTrackingParticles = std::make_unique<TrackingParticleCollection>();
    mergedOutput_.pTrackingVertices = std::make_unique<TrackingVertexCollection>();
    mergedOutput_.refTrackingParticles =
        const_cast<edm::Event &>(event).getRefBeforePut<TrackingParticleCollection>("MergedTrackTruth");
    mergedOutput_.refTrackingVertexes =
        const_cast<edm::Event &>(event).getRefBeforePut<TrackingVertexCollection>("MergedTrackTruth");
  }
 
  if (createInitialVertexCollection_) {
    pInitialVertices_ = std::make_unique<TrackingVertexCollection>();
  }
}
 
 
void TrackingTruthAccumulator::accumulate(edm::Event const &event, edm::EventSetup const &setup) {
  // Call the templated version that does the same for both signal and pileup
  // events
 
  edm::Handle<edm::HepMCProduct> hepmc;
  event.getByLabel(hepMCproductLabel_, hepmc);
 
  accumulateEvent(event, setup, hepmc);
}
 
void TrackingTruthAccumulator::accumulate(PileUpEventPrincipal const &event,
                                          edm::EventSetup const &setup,
                                          edm::StreamID const &) {
  // If this bunch crossing is outside the user configured limit, don't do
  // anything.
  if (event.bunchCrossing() >= -static_cast<int>(maximumPreviousBunchCrossing_) &&
      event.bunchCrossing() <= static_cast<int>(maximumSubsequentBunchCrossing_)) {
    // edm::LogInfo(messageCategory_) << "Analysing pileup event for bunch
    // crossing " << event.bunchCrossing();
 
    // simply create empty handle as we do not have a HepMCProduct in PU anyway
    edm::Handle<edm::HepMCProduct> hepmc;
    accumulateEvent(event, setup, hepmc);
  } else
    edm::LogInfo(messageCategory_) << "Skipping pileup event for bunch crossing " << event.bunchCrossing();
}
 
void TrackingTruthAccumulator::finalizeEvent(edm::Event &event, edm::EventSetup const &setup) {
  if (createUnmergedCollection_) {
    edm::LogInfo("TrackingTruthAccumulator")
        << "Adding " << unmergedOutput_.pTrackingParticles->size() << " TrackingParticles and "
        << unmergedOutput_.pTrackingVertices->size() << " TrackingVertexs to the event.";
 
    event.put(std::move(unmergedOutput_.pTrackingParticles));
    event.put(std::move(unmergedOutput_.pTrackingVertices));
  }
 
  if (createMergedCollection_) {
    edm::LogInfo("TrackingTruthAccumulator")
        << "Adding " << mergedOutput_.pTrackingParticles->size() << " merged TrackingParticles and "
        << mergedOutput_.pTrackingVertices->size() << " merged TrackingVertexs to the event.";
 
    event.put(std::move(mergedOutput_.pTrackingParticles), "MergedTrackTruth");
    event.put(std::move(mergedOutput_.pTrackingVertices), "MergedTrackTruth");
  }
 
  if (createInitialVertexCollection_) {
    edm::LogInfo("TrackingTruthAccumulator")
        << "Adding " << pInitialVertices_->size() << " initial TrackingVertexs to the event.";
 
    event.put(std::move(pInitialVertices_), "InitialVertices");
  }
}
 
template <class T>
void TrackingTruthAccumulator::accumulateEvent(const T &event,
                                               const edm::EventSetup &setup,
                                               const edm::Handle<edm::HepMCProduct> &hepMCproduct) {
  //
  // Get the collections
  //
  edm::Handle<std::vector<SimTrack>> hSimTracks;
  edm::Handle<std::vector<SimVertex>> hSimVertices;
  edm::Handle<std::vector<reco::GenParticle>> hGenParticles;
  edm::Handle<std::vector<int>> hGenParticleIndices;
 
  event.getByLabel(simTrackLabel_, hSimTracks);
  event.getByLabel(simVertexLabel_, hSimVertices);
 
  try {
    event.getByLabel(genParticleLabel_, hGenParticles);
    event.getByLabel(genParticleLabel_, hGenParticleIndices);
  } catch (cms::Exception &exception) {
    //
    // The Monte Carlo is not always available, e.g. for pileup events. The
    // information is only used if it's available, but for some reason the
    // PileUpEventPrincipal wrapper throws an exception here rather than waiting
    // to see if the handle is used (as is the case for edm::Event). So I just
    // want to catch this exception and use the normal handle checking later on.
    //
  }
 
  // Retrieve tracker topology from geometry
  edm::ESHandle<TrackerTopology> tTopoHandle;
  setup.get<TrackerTopologyRcd>().get(tTopoHandle);
  const TrackerTopology *const tTopo = tTopoHandle.product();
 
  // Run through the collections and work out the decay chain of each
  // track/vertex. The information in SimTrack and SimVertex only allows
  // traversing upwards, but this will allow traversal in both directions. This
  // is required for things like grouping electrons that bremsstrahlung as one
  // TrackingParticle if "mergedBremsstrahlung" is set in the config file.
  DecayChain decayChain(*hSimTracks, *hSimVertices);
 
  // I only want to create these collections if they're actually required
  std::unique_ptr<::OutputCollectionWrapper> pUnmergedCollectionWrapper;
  std::unique_ptr<::OutputCollectionWrapper> pMergedCollectionWrapper;
  if (createUnmergedCollection_)
    pUnmergedCollectionWrapper = std::make_unique<::OutputCollectionWrapper>(decayChain, unmergedOutput_);
  if (createMergedCollection_)
    pMergedCollectionWrapper = std::make_unique<::OutputCollectionWrapper>(decayChain, mergedOutput_);
 
  std::vector<const PSimHit *> simHitPointers;
  fillSimHits(simHitPointers, event, setup);
  TrackingParticleFactory objectFactory(decayChain,
                                        hGenParticles,
                                        hepMCproduct,
                                        hGenParticleIndices,
                                        simHitPointers,
                                        volumeRadius_,
                                        volumeZ_,
                                        vertexDistanceCut_,
                                        allowDifferentProcessTypeForDifferentDetectors_);
 
#if defined(DO_DEBUG_TESTING)
  // While I'm testing, perform some checks.
  // TODO - drop this call once I'm happy it works in all situations.
  decayChain.integrityCheck();
#endif
 
  TrackingParticleSelector *pSelector = nullptr;
  if (selectorFlag_)
    pSelector = &selector_;
 
  // Run over all of the SimTracks, but because I'm interested in the decay
  // hierarchy do it through the DecayChainTrack objects. These are looped over
  // in sequence here but they have the hierarchy information for the functions
  // called to traverse the decay chain.
 
  for (size_t index = 0; index < decayChain.decayTracksSize; ++index) {
    ::DecayChainTrack *pDecayTrack = &decayChain.decayTracks[index];
    const SimTrack &simTrack = hSimTracks->at(pDecayTrack->simTrackIndex);
 
    // Perform some quick checks to see if we can drop out early. Note that
    // these are a subset of the cuts in the selector_ so the created
    // TrackingParticle could still fail. The selector_ requires the full
    // TrackingParticle to be made however, which can be computationally
    // expensive.
    if (chargedOnly_ && simTrack.charge() == 0)
      continue;
    if (signalOnly_ && (simTrack.eventId().bunchCrossing() != 0 || simTrack.eventId().event() != 0))
      continue;
 
    // Also perform a check to see if the production vertex is inside the
    // tracker volume (if required).
    if (ignoreTracksOutsideVolume_) {
      const SimVertex &simVertex = hSimVertices->at(pDecayTrack->pParentVertex->simVertexIndex);
      if (!objectFactory.vectorIsInsideVolume(simVertex.position()))
        continue;
    }
 
    // This function creates the TrackinParticle and adds it to the collection
    // if it passes the selection criteria specified in the configuration. If
    // the config specifies adding ancestors, the function is called recursively
    // to do that.
    ::addTrack(pDecayTrack,
               pSelector,
               pUnmergedCollectionWrapper.get(),
               pMergedCollectionWrapper.get(),
               objectFactory,
               addAncestors_,
               tTopo);
  }
 
  // If configured to create a collection of initial vertices, add them from
  // this bunch crossing. No selection is applied on this collection, but it
  // also has no links to the TrackingParticle decay products. There are a lot
  // of "initial vertices", I'm not entirely sure what they all are (nuclear
  // interactions in the detector maybe?), but the one for the main event is the
  // one with vertexId==0.
  if (createInitialVertexCollection_) {
    // Pretty sure the one with vertexId==0 is always the first one, but doesn't
    // hurt to check
    for (const auto &pRootVertex : decayChain.rootVertices) {
      const SimVertex &vertex = hSimVertices->at(decayChain.rootVertices[0]->simVertexIndex);
      if (vertex.vertexId() != 0)
        continue;
 
      pInitialVertices_->push_back(objectFactory.createTrackingVertex(pRootVertex));
      break;
    }
  }
}
 
template <class T>
void TrackingTruthAccumulator::fillSimHits(std::vector<const PSimHit *> &returnValue,
                                           const T &event,
                                           const edm::EventSetup &setup) {
  // loop over the collections
  for (const auto &collectionTag : collectionTags_) {
    edm::Handle<std::vector<PSimHit>> hSimHits;
    event.getByLabel(collectionTag, hSimHits);
 
    // TODO - implement removing the dead modules
    for (const auto &simHit : *hSimHits) {
      returnValue.push_back(&simHit);
    }
 
  }  // end of loop over InputTags
 
  // sort the SimHits according to their time of flight,
  // necessary for looping over them "in order" in
  // TrackingParticleFactory::createTrackingParticle()
  std::sort(returnValue.begin(), returnValue.end(), [](const PSimHit *a, const PSimHit *b) {
    const auto atof =
        edm::isFinite(a->timeOfFlight()) ? a->timeOfFlight() : std::numeric_limits<decltype(a->timeOfFlight())>::max();
    const auto btof =
        edm::isFinite(b->timeOfFlight()) ? b->timeOfFlight() : std::numeric_limits<decltype(b->timeOfFlight())>::max();
    return atof < btof;
  });
}
 
//---------------------------------------------------------------------------------
//---------------------------------------------------------------------------------
//---------------------------------------------------------------------------------
//------ ------
//------  End of the definitions for the TrackingTruthAccumulator methods.
//------
//------  Definitions for the functions and classes declared in the unnamed
//------
//------  namespace are below. Documentation for those is above, where ------
//------  they're declared. ------
//------ ------
//---------------------------------------------------------------------------------
//---------------------------------------------------------------------------------
//---------------------------------------------------------------------------------
 
namespace  // Unnamed namespace for things only used in this file
{
  //---------------------------------------------------------------------------------
  //---------------------------------------------------------------------------------
  //----   TrackingParticleFactory methods
  //----------------------------------------
  //---------------------------------------------------------------------------------
  //---------------------------------------------------------------------------------
 
  ::TrackingParticleFactory::TrackingParticleFactory(const ::DecayChain &decayChain,
                                                     const edm::Handle<std::vector<reco::GenParticle>> &hGenParticles,
                                                     const edm::Handle<edm::HepMCProduct> &hepMCproduct,
                                                     const edm::Handle<std::vector<int>> &hHepMCGenParticleIndices,
                                                     const std::vector<const PSimHit *> &simHits,
                                                     double volumeRadius,
                                                     double volumeZ,
                                                     double vertexDistanceCut,
                                                     bool allowDifferentProcessTypes)
      : decayChain_(decayChain),
        hGenParticles_(hGenParticles),
        hepMCproduct_(hepMCproduct),
        simHits_(simHits),
        volumeRadius_(volumeRadius),
        volumeZ_(volumeZ),
        vertexDistanceCut2_(vertexDistanceCut * vertexDistanceCut),
        allowDifferentProcessTypeForDifferentDetectors_(allowDifferentProcessTypes) {
    // Need to create a multimap to get from a SimTrackId to all of the hits in
    // it. The SimTrackId is an unsigned int.
    for (size_t index = 0; index < simHits_.size(); ++index) {
      trackIdToHitIndex_.insert(std::make_pair(simHits_[index]->trackId(), index));
    }
 
    if (hHepMCGenParticleIndices.isValid())  // Monte Carlo might not be available
                                             // for the pileup events
    {
      genParticleIndices_.resize(hHepMCGenParticleIndices->size() + 1);
 
      // What I need is the reverse mapping of this vector. The sizes are already
      // equivalent because I set the size in the initialiser list.
      for (size_t recoGenParticleIndex = 0; recoGenParticleIndex < hHepMCGenParticleIndices->size();
           ++recoGenParticleIndex) {
        size_t hepMCGenParticleIndex = (*hHepMCGenParticleIndices)[recoGenParticleIndex];
 
        // They should be the same size, give or take a fencepost error, so this
        // should never happen - but just in case
        if (genParticleIndices_.size() <= hepMCGenParticleIndex)
          genParticleIndices_.resize(hepMCGenParticleIndex + 1);
 
        genParticleIndices_[hepMCGenParticleIndex] = recoGenParticleIndex;
      }
    }
  }
 
  TrackingParticle TrackingParticleFactory::createTrackingParticle(const ::DecayChainTrack *pChainTrack,
                                                                   const TrackerTopology *tTopo) const {
    typedef math::XYZTLorentzVectorD LorentzVector;
    typedef math::XYZPoint Vector;
 
    const SimTrack &simTrack = decayChain_.getSimTrack(pChainTrack);
    const SimVertex &parentSimVertex = decayChain_.getSimVertex(pChainTrack->pParentVertex);
 
    LorentzVector position(0, 0, 0, 0);
    if (!simTrack.noVertex())
      position = parentSimVertex.position();
 
    int pdgId = simTrack.type();
 
    TrackingParticle returnValue;
    // N.B. The sim track is added a few lines below, the parent and decay
    // vertices are added in another function later on.
 
    //
    // If there is some valid Monte Carlo for this track, take some information
    // from that. Only do so if it is from the signal event however. Not sure why
    // but that's what the old code did.
    //
    if (simTrack.eventId().event() == 0 &&
        simTrack.eventId().bunchCrossing() == 0)  // if this is a track in the signal event
    {
      int hepMCGenParticleIndex = simTrack.genpartIndex();
      if (hepMCGenParticleIndex >= 0 && hepMCGenParticleIndex < static_cast<int>(genParticleIndices_.size()) &&
          hGenParticles_.isValid()) {
        int recoGenParticleIndex = genParticleIndices_[hepMCGenParticleIndex];
        reco::GenParticleRef generatorParticleRef(hGenParticles_, recoGenParticleIndex);
        pdgId = generatorParticleRef->pdgId();
        returnValue.addGenParticle(generatorParticleRef);
      }
    }
 
    returnValue.addG4Track(simTrack);
 
    // I need to run over the sim hits and count up the different types of hits.
    // To be honest I don't fully understand this next section of code, I copied
    // it from the old TrackingTruthProducer to provide the same results. The
    // different types of hits that I need to count are:
    size_t matchedHits = 0;          // As far as I can tell, this is the number of tracker layers with hits
                                     // on them, i.e. taking account of overlaps.
    size_t numberOfHits = 0;         // The total number of hits not taking account of overlaps.
    size_t numberOfTrackerHits = 0;  // The number of tracker hits not taking account of overlaps.
 
    bool init = true;
    int processType = 0;
    int particleType = 0;
    int oldLayer = 0;
    int newLayer = 0;
    DetId oldDetector;
    DetId newDetector;
 
    // Loop over the SimHits associated to this SimTrack
    // in the order defined by time of flight, which is
    // probably the best quantity available for going
    // through the hits "in order" (ok, most important is
    // to get the first hit right because processType and
    // particleType are taken from it)
    for (auto iHitIndex = trackIdToHitIndex_.lower_bound(simTrack.trackId()),
              end = trackIdToHitIndex_.upper_bound(simTrack.trackId());
         iHitIndex != end;
         ++iHitIndex) {
      const auto &pSimHit = simHits_[iHitIndex->second];
 
      // Skip hits with particle type different from SimTrack pdgId
      if (pSimHit->particleType() != pdgId)
        continue;
 
      // Initial condition for consistent simhit selection
      if (init) {
        processType = pSimHit->processType();
        particleType = pSimHit->particleType();
        newDetector = DetId(pSimHit->detUnitId());
        init = false;
      }
 
      oldDetector = newDetector;
      newDetector = DetId(pSimHit->detUnitId());
 
      // Fast sim seems to have different processType between tracker and muons,
      // so if this flag has been set allow the processType to change if the
      // detector changes.
      if (allowDifferentProcessTypeForDifferentDetectors_ && newDetector.det() != oldDetector.det())
        processType = pSimHit->processType();
 
      // Check for delta and interaction products discards
      if (processType == pSimHit->processType() && particleType == pSimHit->particleType()) {
        ++numberOfHits;
        oldLayer = newLayer;
        newLayer = 0;
        if (newDetector.det() == DetId::Tracker) {
          ++numberOfTrackerHits;
 
          newLayer = tTopo->layer(newDetector);
 
          // Count hits using layers for glued detectors
          if ((oldLayer != newLayer || (oldLayer == newLayer && oldDetector.subdetId() != newDetector.subdetId())))
            ++matchedHits;
        }
      }
    }  // end of loop over the sim hits for this sim track
 
    returnValue.setNumberOfTrackerLayers(matchedHits);
    returnValue.setNumberOfHits(numberOfHits);
    returnValue.setNumberOfTrackerHits(numberOfTrackerHits);
 
    return returnValue;
  }
 
  TrackingVertex TrackingParticleFactory::createTrackingVertex(const ::DecayChainVertex *pChainVertex) const {
    typedef math::XYZTLorentzVectorD LorentzVector;
    typedef math::XYZPoint Vector;
    typedef edm::Ref<edm::HepMCProduct, HepMC::GenVertex> GenVertexRef;
 
    const SimVertex &simVertex = decayChain_.getSimVertex(pChainVertex);
 
    bool isInVolume = this->vectorIsInsideVolume(simVertex.position());
 
    TrackingVertex returnValue(simVertex.position(), isInVolume, simVertex.eventId());
 
    // add the SimVertex to the TrackingVertex
    returnValue.addG4Vertex(simVertex);
 
    // also add refs to nearby HepMC::GenVertexs; the way this is done (i.e. based
    // on the position) is transcribed over from the old TrackingTruthProducer
    if (simVertex.eventId().event() == 0 && simVertex.eventId().bunchCrossing() == 0 &&
        hepMCproduct_.isValid())  // if this is a track in the signal event
    {
      const HepMC::GenEvent *genEvent = hepMCproduct_->GetEvent();
 
      if (genEvent != nullptr) {
        Vector tvPosition(returnValue.position().x(), returnValue.position().y(), returnValue.position().z());
 
        for (HepMC::GenEvent::vertex_const_iterator iGenVertex = genEvent->vertices_begin();
             iGenVertex != genEvent->vertices_end();
             ++iGenVertex) {
          HepMC::ThreeVector rawPosition = (*iGenVertex)->position();
 
          Vector genPosition(rawPosition.x() * 0.1, rawPosition.y() * 0.1, rawPosition.z() * 0.1);
 
          auto distance2 = (tvPosition - genPosition).mag2();
 
          if (distance2 < vertexDistanceCut2_)
            returnValue.addGenVertex(GenVertexRef(hepMCproduct_, (*iGenVertex)->barcode()));
        }
      }
    }
 
    return returnValue;
  }
 
  bool ::TrackingParticleFactory::vectorIsInsideVolume(const math::XYZTLorentzVectorD &vector) const {
    return (vector.Pt() < volumeRadius_ && std::abs(vector.z()) < volumeZ_);
  }
 
  //---------------------------------------------------------------------------------
  //---------------------------------------------------------------------------------
  //----   DecayChain methods
  //-----------------------------------------------------
  //---------------------------------------------------------------------------------
  //---------------------------------------------------------------------------------
 
  ::DecayChain::DecayChain(const std::vector<SimTrack> &trackCollection, const std::vector<SimVertex> &vertexCollection)
      : decayTracksSize(trackCollection.size()),
        decayVerticesSize(vertexCollection.size()),
        decayTracks_(new DecayChainTrack[decayTracksSize]),
        decayVertices_(new DecayChainVertex[decayVerticesSize]),
        simTrackCollection_(trackCollection),
        simVertexCollection_(vertexCollection),
        decayTracks(decayTracks_),  // These const references map onto the actual
                                    // members for public const access
        decayVertices(decayVertices_),
        rootVertices(rootVertices_) {
    // I need some maps to be able to get object pointers from the track/vertex ID
    std::map<int, ::DecayChainTrack *> trackIdToDecayTrack;
    std::map<int, ::DecayChainVertex *> vertexIdToDecayVertex;
 
    // First create a DecayChainTrack for every SimTrack and make a note of the
    // trackIds in the map. Also add a pointer to the daughter list of the parent
    // DecayChainVertex, which might include creating the vertex object if it
    // doesn't already exist.
    size_t decayVertexIndex = 0;  // The index of the next free slot in the DecayChainVertex array.
    for (size_t index = 0; index < trackCollection.size(); ++index) {
      ::DecayChainTrack *pDecayTrack = &decayTracks_[index];  // Get a pointer for ease of use
 
      // This is the array index of the SimTrack that this DecayChainTrack
      // corresponds to. At the moment this is a 1 to 1 relationship with the
      // DecayChainTrack index, but they won't necessarily be accessed through the
      // array later so it's still required to store it.
      pDecayTrack->simTrackIndex = index;
 
      trackIdToDecayTrack[trackCollection[index].trackId()] = pDecayTrack;
 
      int parentVertexIndex = trackCollection[index].vertIndex();
      if (parentVertexIndex >= 0) {
        // Get the DecayChainVertex corresponding to this SimVertex, or initialise
        // it if it hasn't been done already.
        ::DecayChainVertex *&pParentVertex = vertexIdToDecayVertex[parentVertexIndex];
        if (pParentVertex == nullptr) {
          // Note that I'm using a reference, so changing pParentVertex will
          // change the entry in the map too.
          pParentVertex = &decayVertices_[decayVertexIndex];
          ++decayVertexIndex;
          pParentVertex->simVertexIndex = parentVertexIndex;
        }
        pParentVertex->daughterTracks.push_back(pDecayTrack);
        pDecayTrack->pParentVertex = pParentVertex;
      } else
        throw std::runtime_error(
            "TrackingTruthAccumulator: Found a track with "
            "an invalid parent vertex index.");
    }
 
    // This assert was originally in to check the internal consistency of the
    // decay chain. Fast sim pileup seems to have a load of vertices with no
    // tracks pointing to them though, so fast sim fails this assert if pileup is
    // added. I don't think the problem is vital however, so if the assert is
    // taken out these extra vertices are ignored.
    // assert( vertexIdToDecayVertex.size()==vertexCollection.size() &&
    // vertexCollection.size()==decayVertexIndex );
 
    // I still need to set DecayChainTrack::daughterVertices and
    // DecayChainVertex::pParentTrack. The information to do this comes from
    // SimVertex::parentIndex. I couldn't do this before because I need all of the
    // DecayChainTracks initialised.
    for (auto &decayVertexMapPair : vertexIdToDecayVertex) {
      ::DecayChainVertex *pDecayVertex = decayVertexMapPair.second;
      int parentTrackIndex = vertexCollection[pDecayVertex->simVertexIndex].parentIndex();
      if (parentTrackIndex != -1) {
        std::map<int, ::DecayChainTrack *>::iterator iParentTrackMapPair = trackIdToDecayTrack.find(parentTrackIndex);
        if (iParentTrackMapPair == trackIdToDecayTrack.end()) {
          std::stringstream errorStream;
          errorStream << "TrackingTruthAccumulator: Something has gone wrong "
                         "with the indexing. Parent track index is "
                      << parentTrackIndex << ".";
          throw std::runtime_error(errorStream.str());
        }
 
        ::DecayChainTrack *pParentTrackHierarchy = iParentTrackMapPair->second;
 
        pParentTrackHierarchy->daughterVertices.push_back(pDecayVertex);
        pDecayVertex->pParentTrack = pParentTrackHierarchy;
      } else
        rootVertices_.push_back(pDecayVertex);  // Has no parent so is at the top of the decay chain.
    }                                           // end of loop over the vertexIdToDecayVertex map
 
    findBrem(trackCollection, vertexCollection);
 
  }  // end of ::DecayChain constructor
 
#if defined(DO_DEBUG_TESTING)
  // Function documentation is with the declaration above. This function is only
  // used for testing.
  void ::DecayChain::integrityCheck() {
    //
    // Start off checking the tracks
    //
    for (size_t index = 0; index < decayTracksSize; ++index) {
      const auto &decayTrack = decayTracks[index];
      //
      // Tracks should always have a production vertex
      //
      if (decayTrack.pParentVertex == NULL)
        throw std::runtime_error(
            "TrackingTruthAccumulator.cc integrityCheck(): "
            "Found DecayChainTrack with no parent vertex.");
 
      //
      // Make sure the parent has this as a child. Also make sure it's only listed
      // once.
      //
      size_t numberOfTimesListed = 0;
      for (const auto pSiblingTrack : decayTrack.pParentVertex->daughterTracks) {
        if (pSiblingTrack == &decayTrack)
          ++numberOfTimesListed;
      }
      if (numberOfTimesListed != 1)
        throw std::runtime_error(
            "TrackingTruthAccumulator.cc integrityCheck(): "
            "Found DecayChainTrack whose parent does not "
            "have it listed once and only once as a child.");
 
      //
      // Make sure all of the children have this listed as a parent.
      //
      for (const auto pDaughterVertex : decayTrack.daughterVertices) {
        if (pDaughterVertex->pParentTrack != &decayTrack)
          throw std::runtime_error(
              "TrackingTruthAccumulator.cc integrityCheck(): Found "
              "DecayChainTrack whose child does not have it listed as the "
              "parent.");
      }
 
      //
      // Follow the chain up to the root vertex, and make sure that is listed in
      // rootVertices
      //
      const DecayChainVertex *pAncestorVertex = decayTrack.pParentVertex;
      while (pAncestorVertex->pParentTrack != NULL) {
        if (pAncestorVertex->pParentTrack->pParentVertex == NULL)
          throw std::runtime_error(
              "TrackingTruthAccumulator.cc integrityCheck(): Found "
              "DecayChainTrack with no parent vertex higher in the decay chain.");
        pAncestorVertex = pAncestorVertex->pParentTrack->pParentVertex;
      }
      if (std::find(rootVertices.begin(), rootVertices.end(), pAncestorVertex) == rootVertices.end()) {
        throw std::runtime_error(
            "TrackingTruthAccumulator.cc integrityCheck(): Found DecayChainTrack "
            "whose root vertex is not recorded anywhere.");
      }
    }  // end of loop over decayTracks
 
    //
    // Now check the vertices
    //
    for (size_t index = 0; index < decayVerticesSize; ++index) {
      const auto &decayVertex = decayVertices[index];
 
      //
      // Make sure this, or a vertex higher in the chain, is in the list of root
      // vertices.
      //
      const DecayChainVertex *pAncestorVertex = &decayVertex;
      while (pAncestorVertex->pParentTrack != NULL) {
        if (pAncestorVertex->pParentTrack->pParentVertex == NULL)
          throw std::runtime_error(
              "TrackingTruthAccumulator.cc integrityCheck(): Found "
              "DecayChainTrack with no parent vertex higher in the vertex decay "
              "chain.");
        pAncestorVertex = pAncestorVertex->pParentTrack->pParentVertex;
      }
      if (std::find(rootVertices.begin(), rootVertices.end(), pAncestorVertex) == rootVertices.end()) {
        throw std::runtime_error(
            "TrackingTruthAccumulator.cc integrityCheck(): Found DecayChainTrack "
            "whose root vertex is not recorded anywhere.");
      }
 
      //
      // Make sure the parent has this listed in its daughters once and only once.
      //
      if (decayVertex.pParentTrack != NULL) {
        size_t numberOfTimesListed = 0;
        for (const auto pSibling : decayVertex.pParentTrack->daughterVertices) {
          if (pSibling == &decayVertex)
            ++numberOfTimesListed;
        }
        if (numberOfTimesListed != 1)
          throw std::runtime_error(
              "TrackingTruthAccumulator.cc integrityCheck(): Found "
              "DecayChainVertex whose parent does not have it listed once and "
              "only once as a child.");
      }
 
      //
      // Make sure all of the children have this listed as the parent
      //
      for (const auto pDaughter : decayVertex.daughterTracks) {
        if (pDaughter->pParentVertex != &decayVertex)
          throw std::runtime_error(
              "TrackingTruthAccumulator.cc integrityCheck(): Found "
              "DecayChainVertex whose child does not have it listed as the "
              "parent.");
      }
    }  // end of loop over decay vertices
 
    std::cout << "TrackingTruthAccumulator.cc integrityCheck() completed successfully" << std::endl;
  }  // end of ::DecayChain::integrityCheck()
#endif
 
  void ::DecayChain::findBrem(const std::vector<SimTrack> &trackCollection,
                              const std::vector<SimVertex> &vertexCollection) {
    for (size_t index = 0; index < decayVerticesSize; ++index) {
      auto &vertex = decayVertices_[index];
 
      // Make sure parent is an electron
      if (vertex.pParentTrack == nullptr)
        continue;
      int parentTrackPDG = trackCollection[vertex.pParentTrack->simTrackIndex].type();
      if (std::abs(parentTrackPDG) != 11)
        continue;
 
      size_t numberOfElectrons = 0;
      size_t numberOfNonElectronsOrPhotons = 0;
      for (auto &pDaughterTrack : vertex.daughterTracks) {
        const auto &simTrack = trackCollection[pDaughterTrack->simTrackIndex];
        if (simTrack.type() == 11 || simTrack.type() == -11)
          ++numberOfElectrons;
        else if (simTrack.type() != 22)
          ++numberOfNonElectronsOrPhotons;
      }
      if (numberOfElectrons == 1 && numberOfNonElectronsOrPhotons == 0) {
        // This is a valid brem, run through and tell all daughters to use the
        // parent as the brem
        for (auto &pDaughterTrack : vertex.daughterTracks)
          pDaughterTrack->pMergedBremSource = vertex.pParentTrack;
        vertex.pMergedBremSource = vertex.pParentTrack->pParentVertex;
      }
    }
  }  // end of ::DecayChain::findBrem()
 
  //---------------------------------------------------------------------------------
  //---------------------------------------------------------------------------------
  //----   OutputCollectionWrapper methods
  //----------------------------------------
  //---------------------------------------------------------------------------------
  //---------------------------------------------------------------------------------
 
  ::OutputCollectionWrapper::OutputCollectionWrapper(const ::DecayChain &decayChain,
                                                     TrackingTruthAccumulator::OutputCollections &outputCollections)
      : output_(outputCollections),
        trackingParticleIndices_(decayChain.decayTracksSize, -1),
        trackingVertexIndices_(decayChain.decayVerticesSize, -1) {
    // No operation besides the initialiser list
  }
 
  TrackingParticle * ::OutputCollectionWrapper::addTrackingParticle(const ::DecayChainTrack *pDecayTrack,
                                                                    const TrackingParticle &trackingParticle) {
    if (trackingParticleIndices_[pDecayTrack->simTrackIndex] != -1)
      throw std::runtime_error(
          "OutputCollectionWrapper::addTrackingParticle - "
          "trying to add a particle twice");
 
    trackingParticleIndices_[pDecayTrack->simTrackIndex] = output_.pTrackingParticles->size();
    output_.pTrackingParticles->push_back(trackingParticle);
 
    // Clear any associations in case there were any beforehand
    output_.pTrackingParticles->back().clearDecayVertices();
    output_.pTrackingParticles->back().clearParentVertex();
 
    // Associate to the objects that are already in the output collections
    associateToExistingObjects(pDecayTrack);
 
    return &output_.pTrackingParticles->back();
  }
 
  TrackingVertex * ::OutputCollectionWrapper::addTrackingVertex(const ::DecayChainVertex *pDecayVertex,
                                                                const TrackingVertex &trackingVertex) {
    if (trackingVertexIndices_[pDecayVertex->simVertexIndex] != -1)
      throw std::runtime_error(
          "OutputCollectionWrapper::addTrackingVertex - "
          "trying to add a vertex twice");
 
    trackingVertexIndices_[pDecayVertex->simVertexIndex] = output_.pTrackingVertices->size();
    output_.pTrackingVertices->push_back(trackingVertex);
 
    // Associate the new vertex to any parents or children that are already in the
    // output collections
    associateToExistingObjects(pDecayVertex);
 
    return &output_.pTrackingVertices->back();
  }
 
  TrackingParticle * ::OutputCollectionWrapper::getTrackingParticle(const ::DecayChainTrack *pDecayTrack) {
    const int index = trackingParticleIndices_[pDecayTrack->simTrackIndex];
    if (index == -1)
      return nullptr;
    else
      return &(*output_.pTrackingParticles)[index];
  }
 
  TrackingVertex * ::OutputCollectionWrapper::getTrackingVertex(const ::DecayChainVertex *pDecayVertex) {
    const int index = trackingVertexIndices_[pDecayVertex->simVertexIndex];
    if (index == -1)
      return nullptr;
    else
      return &(*output_.pTrackingVertices)[index];
  }
 
  TrackingParticleRef OutputCollectionWrapper::getRef(const ::DecayChainTrack *pDecayTrack) {
    const int index = trackingParticleIndices_[pDecayTrack->simTrackIndex];
    if (index == -1)
      throw std::runtime_error(
          "OutputCollectionWrapper::getRefTrackingParticle - ref requested for a "
          "non existent TrackingParticle");
    else
      return TrackingParticleRef(output_.refTrackingParticles, index);
  }
 
  TrackingVertexRef OutputCollectionWrapper::getRef(const ::DecayChainVertex *pDecayVertex) {
    const int index = trackingVertexIndices_[pDecayVertex->simVertexIndex];
    if (index == -1)
      throw std::runtime_error(
          "OutputCollectionWrapper::getRefTrackingParticle - ref requested for a "
          "non existent TrackingVertex");
    else
      return TrackingVertexRef(output_.refTrackingVertexes, index);
  }
 
  void ::OutputCollectionWrapper::setProxy(const ::DecayChainTrack *pOriginalTrack,
                                           const ::DecayChainTrack *pProxyTrack) {
    int &index = trackingParticleIndices_[pOriginalTrack->simTrackIndex];
    if (index != -1)
      throw std::runtime_error(
          "OutputCollectionWrapper::setProxy() was called for a TrackingParticle "
          "that has already been created");
    // Note that index is a reference so this call changes the value in the vector
    // too
    index = trackingParticleIndices_[pProxyTrack->simTrackIndex];
  }
 
  void ::OutputCollectionWrapper::setProxy(const ::DecayChainVertex *pOriginalVertex,
                                           const ::DecayChainVertex *pProxyVertex) {
    int &index = trackingVertexIndices_[pOriginalVertex->simVertexIndex];
    const int newIndex = trackingVertexIndices_[pProxyVertex->simVertexIndex];
 
    if (index != -1 && index != newIndex)
      throw std::runtime_error(
          "OutputCollectionWrapper::setProxy() was called for a TrackingVertex "
          "that has already been created");
 
    // Note that index is a reference so this call changes the value in the vector
    // too
    index = newIndex;
  }
 
  void ::OutputCollectionWrapper::associateToExistingObjects(const ::DecayChainVertex *pChainVertex) {
    // First make sure the DecayChainVertex supplied has been turned into a
    // TrackingVertex
    TrackingVertex *pTrackingVertex = getTrackingVertex(pChainVertex);
    if (pTrackingVertex == nullptr)
      throw std::runtime_error("associateToExistingObjects was passed a non existent TrackingVertex");
 
    //
    // Associate to the parent track (if there is one)
    //
    ::DecayChainTrack *pParentChainTrack = pChainVertex->pParentTrack;
    if (pParentChainTrack != nullptr)  // Make sure there is a parent track first
    {
      // There is a parent track, but it might not have been turned into a
      // TrackingParticle yet
      TrackingParticle *pParentTrackingParticle = getTrackingParticle(pParentChainTrack);
      if (pParentTrackingParticle != nullptr) {
        pParentTrackingParticle->addDecayVertex(getRef(pChainVertex));
        pTrackingVertex->addParentTrack(getRef(pParentChainTrack));
      }
      // Don't worry about the "else" case - the parent track might not
      // necessarily get added to the output collection at all. A check is done on
      // daughter vertices when tracks are added, so the association will be
      // picked up then.
    }
 
    //
    // A parent TrackingVertex is always associated to a daughter TrackingParticle
    // when the TrackingParticle is created. Hence there is no need to check the
    // list of daughter tracks.
    //
  }
 
  void ::OutputCollectionWrapper::associateToExistingObjects(const ::DecayChainTrack *pChainTrack) {
    //
    // First make sure this DecayChainTrack has been turned into a
    // TrackingParticle
    //
    TrackingParticle *pTrackingParticle = getTrackingParticle(pChainTrack);
    if (pTrackingParticle == nullptr)
      throw std::runtime_error(
          "associateToExistingObjects was passed a non "
          "existent TrackingParticle");
 
    // Get the parent vertex. This should always already have been turned into a
    // TrackingVertex, and there should always be a parent DecayChainVertex.
    ::DecayChainVertex *pParentChainVertex = pChainTrack->pParentVertex;
    TrackingVertex *pParentTrackingVertex = getTrackingVertex(pParentChainVertex);
 
    //
    // First associate to the parent vertex.
    //
    pTrackingParticle->setParentVertex(getRef(pParentChainVertex));
    pParentTrackingVertex->addDaughterTrack(getRef(pChainTrack));
 
    // If there are any daughter vertices that have already been made into a
    // TrackingVertex make sure they know about each other. If the daughter
    // vertices haven't been made into TrackingVertexs yet, I won't do anything
    // and create the association when the vertex is made.
    for (auto pDaughterChainVertex : pChainTrack->daughterVertices) {
      TrackingVertex *pDaughterTrackingVertex = getTrackingVertex(pDaughterChainVertex);
      if (pDaughterTrackingVertex != nullptr) {
        pTrackingParticle->addDecayVertex(getRef(pDaughterChainVertex));
        pDaughterTrackingVertex->addParentTrack(getRef(pChainTrack));
      }
    }
  }
 
  TrackingParticle *addTrackAndParentVertex(::DecayChainTrack *pDecayTrack,
                                            const TrackingParticle &trackingParticle,
                                            ::OutputCollectionWrapper *pOutput) {
    // See if this TrackingParticle has already been created (could be if the
    // DecayChainTracks are looped over in a funny order). If it has then there's
    // no need to do anything.
    TrackingParticle *pTrackingParticle = pOutput->getTrackingParticle(pDecayTrack);
    if (pTrackingParticle == nullptr) {
      // Need to make sure the production vertex has been created first
      if (pOutput->getTrackingVertex(pDecayTrack->pParentVertex) == nullptr) {
        // TrackingVertex doesn't exist in the output collection yet. However,
        // it's already been created in the addTrack() function and a temporary
        // reference to it set in the TrackingParticle. I'll use that reference to
        // create a copy in the output collection. When the TrackingParticle is
        // added to the output collection a few lines below the temporary
        // reference to the parent vertex will be cleared, and the correct one
        // referring to the output collection will be set.
        pOutput->addTrackingVertex(pDecayTrack->pParentVertex, *trackingParticle.parentVertex());
      }
 
      pTrackingParticle = pOutput->addTrackingParticle(pDecayTrack, trackingParticle);
    }
 
    return pTrackingParticle;
  }
 
  void addTrack(::DecayChainTrack *pDecayChainTrack,
                const TrackingParticleSelector *pSelector,
                ::OutputCollectionWrapper *pUnmergedOutput,
                ::OutputCollectionWrapper *pMergedOutput,
                const ::TrackingParticleFactory &objectFactory,
                bool addAncestors,
                const TrackerTopology *tTopo) {
    if (pDecayChainTrack == nullptr)
      return;  // This is required for when the addAncestors_ recursive call
               // reaches the top of the chain
 
    // Check to see if this particle has already been processed while traversing
    // up the parents of another split in the decay chain. The check in the line
    // above only stops when the top of the chain is reached, whereas this will
    // stop when a previously traversed split is reached.
    {  // block to limit the scope of temporary variables
      bool alreadyProcessed = true;
      if (pUnmergedOutput != nullptr) {
        if (pUnmergedOutput->getTrackingParticle(pDecayChainTrack) == nullptr)
          alreadyProcessed = false;
      }
      if (pMergedOutput != nullptr) {
        if (pMergedOutput->getTrackingParticle(pDecayChainTrack) == nullptr)
          alreadyProcessed = false;
      }
      if (alreadyProcessed)
        return;
    }
 
    // Create a TrackingParticle.
    TrackingParticle newTrackingParticle = objectFactory.createTrackingParticle(pDecayChainTrack, tTopo);
 
    // The selector checks the impact parameters from the vertex, so I need to
    // have a valid reference to the parent vertex in the TrackingParticle before
    // that can be called. TrackingParticle needs an edm::Ref for the parent
    // TrackingVertex though. I still don't know if this is going to be added to
    // the collection so I can't take it from there, so I need to create a
    // temporary one. When the addTrackAndParentVertex() is called (assuming it
    // passes selection) it will use the temporary reference to create a copy of
    // the parent vertex, put that in the output collection, and then set the
    // reference in the TrackingParticle properly.
    TrackingVertexCollection dummyCollection;  // Only needed to create an edm::Ref
    dummyCollection.push_back(objectFactory.createTrackingVertex(pDecayChainTrack->pParentVertex));
    TrackingVertexRef temporaryRef(&dummyCollection, 0);
    newTrackingParticle.setParentVertex(temporaryRef);
 
    // If a selector has been supplied apply it on the new TrackingParticle and
    // return if it fails.
    if (pSelector) {
      if (!(*pSelector)(newTrackingParticle))
        return;  // Return if the TrackingParticle fails selection
    }
 
    // Add the ancestors first (if required) so that the collection has some kind
    // of chronological order. I don't know how important that is but other code
    // might assume chronological order. If adding ancestors, no selection is
    // applied. Note that I've already checked that all DecayChainTracks have a
    // pParentVertex.
    if (addAncestors)
      addTrack(pDecayChainTrack->pParentVertex->pParentTrack,
               nullptr,
               pUnmergedOutput,
               pMergedOutput,
               objectFactory,
               addAncestors,
               tTopo);
 
    // If creation of the unmerged collection has been turned off in the config
    // this pointer will be null.
    if (pUnmergedOutput != nullptr)
      addTrackAndParentVertex(pDecayChainTrack, newTrackingParticle, pUnmergedOutput);
 
    // If creation of the merged collection has been turned off in the config this
    // pointer will be null.
    if (pMergedOutput != nullptr) {
      ::DecayChainTrack *pBremParentChainTrack = pDecayChainTrack;
      while (pBremParentChainTrack->pMergedBremSource != nullptr)
        pBremParentChainTrack = pBremParentChainTrack->pMergedBremSource;
 
      if (pBremParentChainTrack != pDecayChainTrack) {
        TrackingParticle *pBremParentTrackingParticle =
            addTrackAndParentVertex(pBremParentChainTrack, newTrackingParticle, pMergedOutput);
        // The original particle in the bremsstrahlung decay chain has been
        // created (or retrieved if it already existed), now copy in the
        // extra information.
        // TODO - copy extra information.
 
        if (std::abs(newTrackingParticle.pdgId()) == 22) {
          // Photons are added as separate TrackingParticles, but with the
          // production vertex changed to be the production vertex of the original
          // electron.
 
          // Set up a proxy, so that all requests for the parent TrackingVertex
          // get redirected to the brem parent's TrackingVertex.
          pMergedOutput->setProxy(pDecayChainTrack->pParentVertex, pBremParentChainTrack->pParentVertex);
 
          // Now that pMergedOutput thinks the parent vertex is the brem parent's
          // vertex I can call this and it will set the TrackingParticle parent
          // vertex correctly to the brem parent vertex.
          addTrackAndParentVertex(pDecayChainTrack, newTrackingParticle, pMergedOutput);
        } else if (std::abs(newTrackingParticle.pdgId()) == 11) {
          // Electrons have their G4 tracks and SimHits copied to the parent
          // TrackingParticle.
          for (const auto &trackSegment : newTrackingParticle.g4Tracks()) {
            pBremParentTrackingParticle->addG4Track(trackSegment);
          }
 
          // Also copy the generator particle references
          for (const auto &genParticleRef : newTrackingParticle.genParticles()) {
            pBremParentTrackingParticle->addGenParticle(genParticleRef);
          }
 
          pBremParentTrackingParticle->setNumberOfHits(pBremParentTrackingParticle->numberOfHits() +
                                                       newTrackingParticle.numberOfHits());
          pBremParentTrackingParticle->setNumberOfTrackerHits(pBremParentTrackingParticle->numberOfTrackerHits() +
                                                              newTrackingParticle.numberOfTrackerHits());
          pBremParentTrackingParticle->setNumberOfTrackerLayers(pBremParentTrackingParticle->numberOfTrackerLayers() +
                                                                newTrackingParticle.numberOfTrackerLayers());
 
          // Set a proxy in the output collection wrapper so that any attempt to
          // get objects for this DecayChainTrack again get redirected to the brem
          // parent.
          pMergedOutput->setProxy(pDecayChainTrack, pBremParentChainTrack);
        }
      } else {
        // This is not the result of bremsstrahlung so add to the collection as
        // normal.
        addTrackAndParentVertex(pDecayChainTrack, newTrackingParticle, pMergedOutput);
      }
    }  // end of "if( pMergedOutput!=NULL )", i.e. end of "if bremsstrahlung
       // merging is turned on"
 
  }  // end of addTrack function
 
}  // namespace
 
// Register with the framework
DEFINE_DIGI_ACCUMULATOR(TrackingTruthAccumulator);