Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-09-10 02:59:11

0001 /** @file
0002  * @brief Definitions of the TrackingTruthAccumulator methods.
0003  *
0004  * There are a lot of utility classes and functions in this file. This makes it
0005  * quite long but I didn't want to confuse the directory structure with lots of
0006  * extra files. My reasoning was that lots of people look at the directory
0007  * contents but only the really interested ones will look in this particular
0008  * file, and the utility stuff isn't used elsewhere.
0009  *
0010  * There are three main sections to this file, from top to bottom:
0011  * 1 - Declarations of extra utility stuff not in the header file. All in the
0012  * unnamed namespace. 2 - Definitions of the TrackingTruthAccumulator methods.
0013  * 3 - Definitions of the utility classes and functions declared in (1).
0014  *
0015  * @author Mark Grimes (mark.grimes@bristol.ac.uk)
0016  * @date circa Oct/2012 to Feb/2013
0017  *
0018  * Changelog:
0019  * 05/May/2015 Mark Grimes - Added functionality to add a collection of just the
0020  * initial vertices for FastTimer studies.
0021  *
0022  * 17/Jul/2014 Dominik Nowatschin (dominik.nowatschin@cern.ch) - added SimVertex
0023  * and a ref to HepMC::Genvertex to TrackingVertex in
0024  * TrackingParticleFactory::createTrackingVertex; handle to edm::HepMCProduct is
0025  * created directly in TrackingTruthAccumulator::accumulate and not in
0026  * accumulateEvent as edm::Event and PileUpEventPrincipal have different
0027  * getByLabel() functions
0028  *
0029  * 07/Feb/2013 Mark Grimes - Reorganised and added a bit more documentation.
0030  * Still not enough though. 12/Mar/2012 Mark Grimes - Updated TrackingParticle
0031  * creation to fit in with Subir Sarkar's re-implementation of TrackingParticle.
0032  */
0033 #include "SimGeneral/TrackingAnalysis/plugins/TrackingTruthAccumulator.h"
0034 
0035 #include <memory>
0036 
0037 #include "FWCore/Framework/interface/ConsumesCollector.h"
0038 #include "FWCore/Framework/interface/Event.h"
0039 #include "FWCore/Framework/interface/EventSetup.h"
0040 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0041 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0042 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0043 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0044 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h"
0045 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0046 #include "SimGeneral/MixingModule/interface/DigiAccumulatorMixModFactory.h"
0047 #include "SimGeneral/MixingModule/interface/PileUpEventPrincipal.h"
0048 #include "SimGeneral/TrackingAnalysis/interface/EncodedTruthId.h"
0049 
0050 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0051 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0052 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0053 #include "FWCore/Utilities/interface/isFinite.h"
0054 
0055 // Turn on integrity checking
0056 //#define DO_DEBUG_TESTING
0057 
0058 //---------------------------------------------------------------------------------
0059 //---------------------------------------------------------------------------------
0060 //---------------------------------------------------------------------------------
0061 //------ ------
0062 //------  Declarations for utility classes and functions in the unnamed ------
0063 //------  namespace. The definitions are right at the bottom of the file. ------
0064 //------ ------
0065 //---------------------------------------------------------------------------------
0066 //---------------------------------------------------------------------------------
0067 //---------------------------------------------------------------------------------
0068 
0069 namespace {
0070   /** @brief Class to represent tracks in the decay chain.
0071  * @author Mark Grimes (mark.grimes@bristol.ac.uk)
0072  * @date 30/Oct/2012
0073  */
0074   struct DecayChainTrack {
0075     int simTrackIndex;
0076     struct DecayChainVertex *pParentVertex;
0077     // Some things have multiple decay vertices. Not sure how that works but it
0078     // seems to be mostly electrons and some photons.
0079     std::vector<struct DecayChainVertex *> daughterVertices;
0080     DecayChainTrack *pMergedBremSource;
0081     DecayChainTrack() : simTrackIndex(-1), pParentVertex(nullptr), pMergedBremSource(nullptr) {}
0082     DecayChainTrack(int newSimTrackIndex)
0083         : simTrackIndex(newSimTrackIndex), pParentVertex(nullptr), pMergedBremSource() {}
0084   };
0085 
0086   /** @brief Class to represent a vertex in the decay chain, and it's relationship
0087  * with parents and daughters.
0088  * @author Mark Grimes (mark.grimes@bristol.ac.uk)
0089  * @date 30/Oct/2012
0090  */
0091   struct DecayChainVertex {
0092     int simVertexIndex;
0093     DecayChainTrack *pParentTrack;
0094     std::vector<DecayChainTrack *> daughterTracks;
0095     DecayChainVertex *pMergedBremSource;
0096     DecayChainVertex() : simVertexIndex(-1), pParentTrack(nullptr), pMergedBremSource(nullptr) {}
0097     DecayChainVertex(int newIndex) : simVertexIndex(newIndex), pParentTrack(nullptr), pMergedBremSource(nullptr) {}
0098   };
0099 
0100   /** @brief Intermediary class. Mainly here to handle memory safely.
0101  *
0102  * Reconstructs the parent and child information from SimVertex and SimTracks to
0103  * create the decay chain. SimVertex and SimTrack only have information about
0104  * their parent, so to get information about the children the whole collection
0105  * has to be analysed. This function does that and returns a tree of
0106  * DecayChainTrack and DecayChainVertex classes that represent the full decay
0107  * chain.
0108  *
0109  * @author Mark Grimes (mark.grimes@bristol.ac.uk)
0110  * @date 31/Oct/2012
0111  */
0112   struct DecayChain {
0113   public:
0114     DecayChain(const std::vector<SimTrack> &trackCollection, const std::vector<SimVertex> &vertexCollection);
0115     const size_t decayTracksSize;
0116     const size_t decayVerticesSize;
0117 
0118 #if defined(DO_DEBUG_TESTING)
0119     /** @brief Testing check. Won't actually be called when the code becomes
0120    * production.
0121    *
0122    * Checks that there are no dangling objects not associated in the decay
0123    * chain.
0124    */
0125     void integrityCheck();
0126 #endif
0127     const SimTrack &getSimTrack(const ::DecayChainTrack *pDecayTrack) const {
0128       return simTrackCollection_.at(pDecayTrack->simTrackIndex);
0129     }
0130     const SimVertex &getSimVertex(const ::DecayChainVertex *pDecayVertex) const {
0131       return simVertexCollection_.at(pDecayVertex->simVertexIndex);
0132     }
0133 
0134   private:
0135     void findBrem(const std::vector<SimTrack> &trackCollection, const std::vector<SimVertex> &vertexCollection);
0136     std::unique_ptr<::DecayChainTrack[]> decayTracks_;
0137     std::unique_ptr<::DecayChainVertex[]> decayVertices_;
0138 
0139     /// The vertices at the top of the decay chain. These are a subset of the ones
0140     /// in the decayVertices member. Don't delete them.
0141     std::vector<::DecayChainVertex *> rootVertices_;
0142 
0143     // Keep a record of the constructor parameters
0144     const std::vector<SimTrack> &simTrackCollection_;
0145     const std::vector<SimVertex> &simVertexCollection_;
0146 
0147   public:
0148     const std::unique_ptr<::DecayChainTrack[]> &decayTracks;  ///< Reference maps to decayTracks_ for easy external const
0149                                                               ///< access.
0150     const std::unique_ptr<::DecayChainVertex[]> &decayVertices;  ///< Reference maps to decayVertices_ for easy external
0151                                                                  ///< const access.
0152     const std::vector<::DecayChainVertex *> &rootVertices;       ///< Reference maps to rootVertices_ for easy external
0153                                                                  ///< const access.
0154   };
0155 
0156   /** @brief Class to create TrackingParticle and TrackingVertex objects.
0157  * @author Mark Grimes (mark.grimes@bristol.ac.uk)
0158  * @date 12/Nov/2012
0159  */
0160   class TrackingParticleFactory {
0161   public:
0162     TrackingParticleFactory(const ::DecayChain &decayChain,
0163                             const edm::Handle<std::vector<reco::GenParticle>> &hGenParticles,
0164                             const edm::Handle<edm::HepMCProduct> &hepMCproduct,
0165                             const edm::Handle<std::vector<int>> &hHepMCGenParticleIndices,
0166                             const std::vector<const PSimHit *> &simHits,
0167                             double volumeRadius,
0168                             double volumeZ,
0169                             double vertexDistanceCut,
0170                             bool allowDifferentProcessTypes);
0171     TrackingParticle createTrackingParticle(const DecayChainTrack *pTrack, const TrackerTopology *tTopo) const;
0172     TrackingVertex createTrackingVertex(const DecayChainVertex *pVertex) const;
0173     bool vectorIsInsideVolume(const math::XYZTLorentzVectorD &vector) const;
0174 
0175   private:
0176     const ::DecayChain &decayChain_;
0177     const edm::Handle<std::vector<reco::GenParticle>> &hGenParticles_;
0178     const edm::Handle<edm::HepMCProduct> &hepMCproduct_;
0179     std::vector<int> genParticleIndices_;
0180     const std::vector<const PSimHit *> &simHits_;
0181     const double volumeRadius_;
0182     const double volumeZ_;
0183     const double vertexDistanceCut2_;                        // distance based on which HepMC::GenVertexs
0184                                                              // are added to SimVertexs
0185     std::multimap<unsigned int, size_t> trackIdToHitIndex_;  ///< A multimap linking SimTrack::trackId() to the hit
0186                                                              ///< index in pSimHits_
0187     bool allowDifferentProcessTypeForDifferentDetectors_;    ///< See the comment for
0188                                                              ///< the same member in
0189                                                              ///< TrackingTruthAccumulator
0190   };
0191 
0192   /** @brief Handles adding objects to the output collections correctly, and
0193  * checks to see if a particular object already exists.
0194  * @author Mark Grimes (mark.grimes@bristol.ac.uk)
0195  * @date 09/Nov/2012
0196  */
0197   class OutputCollectionWrapper {
0198   public:
0199     OutputCollectionWrapper(const DecayChain &decayChain,
0200                             TrackingTruthAccumulator::OutputCollections &outputCollections);
0201     TrackingParticle *addTrackingParticle(const ::DecayChainTrack *pDecayTrack,
0202                                           const TrackingParticle &trackingParticle);
0203     TrackingVertex *addTrackingVertex(const ::DecayChainVertex *pDecayVertex, const TrackingVertex &trackingVertex);
0204     TrackingParticle *getTrackingParticle(const ::DecayChainTrack *pDecayTrack);
0205     /// @brief After this call, any call to getTrackingParticle or getRef with
0206     /// pOriginalTrack will return the TrackingParticle for pProxyTrack instead.
0207     void setProxy(const ::DecayChainTrack *pOriginalTrack, const ::DecayChainTrack *pProxyTrack);
0208     void setProxy(const ::DecayChainVertex *pOriginalVertex, const ::DecayChainVertex *pProxyVertex);
0209     TrackingVertex *getTrackingVertex(const ::DecayChainVertex *pDecayVertex);
0210     TrackingParticleRef getRef(const ::DecayChainTrack *pDecayTrack);
0211     TrackingVertexRef getRef(const ::DecayChainVertex *pDecayVertex);
0212 
0213   private:
0214     void associateToExistingObjects(const ::DecayChainVertex *pChainVertex);
0215     void associateToExistingObjects(const ::DecayChainTrack *pChainTrack);
0216     TrackingTruthAccumulator::OutputCollections &output_;
0217     std::vector<int> trackingParticleIndices_;
0218     std::vector<int> trackingVertexIndices_;
0219   };
0220 
0221   /** @brief Adds the supplied TrackingParticle and its parent TrackingVertex to
0222  * the output collection. Checks to make sure they don't already exist first.
0223  * @author Mark Grimes (mark.grimes@bristol.ac.uk)
0224  * @date 12/Nov/2012
0225  */
0226   TrackingParticle *addTrackAndParentVertex(::DecayChainTrack *pDecayTrack,
0227                                             const TrackingParticle &trackingParticle,
0228                                             ::OutputCollectionWrapper *pOutput);
0229 
0230   /** @brief Creates a TrackingParticle for the supplied DecayChainTrack and adds
0231  * it to the output if it passes selection. Gets called recursively if adding
0232  * parents.
0233  * @author Mark Grimes (mark.grimes@bristol.ac.uk)
0234  * @date 05/Nov/2012
0235  */
0236   void addTrack(::DecayChainTrack *pDecayChainTrack,
0237                 const TrackingParticleSelector *pSelector,
0238                 ::OutputCollectionWrapper *pUnmergedOutput,
0239                 ::OutputCollectionWrapper *pMergedOutput,
0240                 const ::TrackingParticleFactory &objectFactory,
0241                 bool addAncestors,
0242                 const TrackerTopology *tTopo);
0243 
0244 }  // namespace
0245 
0246 //---------------------------------------------------------------------------------
0247 //---------------------------------------------------------------------------------
0248 //---------------------------------------------------------------------------------
0249 //------ ------
0250 //------  Definitions for the TrackingTruthAccumulator methods ------
0251 //------  are below. ------
0252 //------ ------
0253 //---------------------------------------------------------------------------------
0254 //---------------------------------------------------------------------------------
0255 //---------------------------------------------------------------------------------
0256 
0257 TrackingTruthAccumulator::TrackingTruthAccumulator(const edm::ParameterSet &config,
0258                                                    edm::ProducesCollector producesCollector,
0259                                                    edm::ConsumesCollector &iC)
0260     : messageCategory_("TrackingTruthAccumulator"),
0261       volumeRadius_(config.getParameter<double>("volumeRadius")),
0262       volumeZ_(config.getParameter<double>("volumeZ")),
0263       vertexDistanceCut_(config.getParameter<double>("vertexDistanceCut")),
0264       ignoreTracksOutsideVolume_(config.getParameter<bool>("ignoreTracksOutsideVolume")),
0265       maximumPreviousBunchCrossing_(config.getParameter<unsigned int>("maximumPreviousBunchCrossing")),
0266       maximumSubsequentBunchCrossing_(config.getParameter<unsigned int>("maximumSubsequentBunchCrossing")),
0267       createUnmergedCollection_(config.getParameter<bool>("createUnmergedCollection")),
0268       createMergedCollection_(config.getParameter<bool>("createMergedBremsstrahlung")),
0269       createInitialVertexCollection_(config.getParameter<bool>("createInitialVertexCollection")),
0270       addAncestors_(config.getParameter<bool>("alwaysAddAncestors")),
0271       removeDeadModules_(config.getParameter<bool>("removeDeadModules")),
0272       simTrackLabel_(config.getParameter<edm::InputTag>("simTrackCollection")),
0273       simVertexLabel_(config.getParameter<edm::InputTag>("simVertexCollection")),
0274       collectionTags_(),
0275       genParticleLabel_(config.getParameter<edm::InputTag>("genParticleCollection")),
0276       hepMCproductLabel_(config.getParameter<edm::InputTag>("HepMCProductLabel")),
0277       tTopoToken_(iC.esConsumes()),
0278       allowDifferentProcessTypeForDifferentDetectors_(config.getParameter<bool>("allowDifferentSimHitProcesses")) {
0279   //
0280   // Make sure at least one of the merged and unmerged collections have been set
0281   // to be created.
0282   //
0283   if (!createUnmergedCollection_ && !createMergedCollection_)
0284     edm::LogError(messageCategory_) << "Both \"createUnmergedCollection\" and "
0285                                        "\"createMergedBremsstrahlung\" have been"
0286                                     << "set to false, which means no collections will be created";
0287 
0288   // Initialize selection for building TrackingParticles
0289   //
0290   if (config.exists("select")) {
0291     edm::ParameterSet param = config.getParameter<edm::ParameterSet>("select");
0292     selector_ = TrackingParticleSelector(param.getParameter<double>("ptMinTP"),
0293                                          param.getParameter<double>("ptMaxTP"),
0294                                          param.getParameter<double>("minRapidityTP"),
0295                                          param.getParameter<double>("maxRapidityTP"),
0296                                          param.getParameter<double>("tipTP"),
0297                                          param.getParameter<double>("lipTP"),
0298                                          param.getParameter<int>("minHitTP"),
0299                                          param.getParameter<bool>("signalOnlyTP"),
0300                                          param.getParameter<bool>("intimeOnlyTP"),
0301                                          param.getParameter<bool>("chargedOnlyTP"),
0302                                          param.getParameter<bool>("stableOnlyTP"),
0303                                          param.getParameter<std::vector<int>>("pdgIdTP"));
0304     selectorFlag_ = true;
0305 
0306     // Also set these two variables, which are used to drop out early if the
0307     // SimTrack doesn't conform. The selector_ requires a full TrackingParticle
0308     // object, but these two variables can veto things early.
0309     chargedOnly_ = param.getParameter<bool>("chargedOnlyTP");
0310     signalOnly_ = param.getParameter<bool>("signalOnlyTP");
0311   } else {
0312     selectorFlag_ = false;
0313     chargedOnly_ = false;
0314     signalOnly_ = false;
0315   }
0316 
0317   //
0318   // Need to state what collections are going to be added to the event. This
0319   // depends on which of the merged and unmerged collections have been
0320   // configured to be created.
0321   //
0322   if (createUnmergedCollection_) {
0323     producesCollector.produces<TrackingVertexCollection>();
0324     producesCollector.produces<TrackingParticleCollection>();
0325   }
0326 
0327   if (createMergedCollection_) {
0328     producesCollector.produces<TrackingParticleCollection>("MergedTrackTruth");
0329     producesCollector.produces<TrackingVertexCollection>("MergedTrackTruth");
0330   }
0331 
0332   if (createInitialVertexCollection_) {
0333     producesCollector.produces<TrackingVertexCollection>("InitialVertices");
0334   }
0335 
0336   iC.consumes<std::vector<SimTrack>>(simTrackLabel_);
0337   iC.consumes<std::vector<SimVertex>>(simVertexLabel_);
0338   iC.consumes<std::vector<reco::GenParticle>>(genParticleLabel_);
0339   iC.consumes<std::vector<int>>(genParticleLabel_);
0340   iC.consumes<std::vector<int>>(hepMCproductLabel_);
0341 
0342   // Fill the collection tags
0343   const edm::ParameterSet &simHitCollectionConfig = config.getParameterSet("simHitCollections");
0344   std::vector<std::string> parameterNames = simHitCollectionConfig.getParameterNames();
0345 
0346   for (const auto &parameterName : parameterNames) {
0347     std::vector<edm::InputTag> tags = simHitCollectionConfig.getParameter<std::vector<edm::InputTag>>(parameterName);
0348     collectionTags_.insert(collectionTags_.end(), tags.begin(), tags.end());
0349   }
0350 
0351   for (const auto &collectionTag : collectionTags_) {
0352     iC.consumes<std::vector<PSimHit>>(collectionTag);
0353   }
0354 }
0355 
0356 void TrackingTruthAccumulator::initializeEvent(edm::Event const &event, edm::EventSetup const &setup) {
0357   if (createUnmergedCollection_) {
0358     unmergedOutput_.pTrackingParticles = std::make_unique<TrackingParticleCollection>();
0359     unmergedOutput_.pTrackingVertices = std::make_unique<TrackingVertexCollection>();
0360     unmergedOutput_.refTrackingParticles =
0361         const_cast<edm::Event &>(event).getRefBeforePut<TrackingParticleCollection>();
0362     unmergedOutput_.refTrackingVertexes = const_cast<edm::Event &>(event).getRefBeforePut<TrackingVertexCollection>();
0363   }
0364 
0365   if (createMergedCollection_) {
0366     mergedOutput_.pTrackingParticles = std::make_unique<TrackingParticleCollection>();
0367     mergedOutput_.pTrackingVertices = std::make_unique<TrackingVertexCollection>();
0368     mergedOutput_.refTrackingParticles =
0369         const_cast<edm::Event &>(event).getRefBeforePut<TrackingParticleCollection>("MergedTrackTruth");
0370     mergedOutput_.refTrackingVertexes =
0371         const_cast<edm::Event &>(event).getRefBeforePut<TrackingVertexCollection>("MergedTrackTruth");
0372   }
0373 
0374   if (createInitialVertexCollection_) {
0375     pInitialVertices_ = std::make_unique<TrackingVertexCollection>();
0376   }
0377 }
0378 
0379 /// create handle to edm::HepMCProduct here because event.getByLabel with
0380 /// edm::HepMCProduct only works for edm::Event but not for
0381 /// PileUpEventPrincipal; PileUpEventPrincipal::getByLabel tries to call
0382 /// T::value_type and T::iterator (where T is the type of the object one wants
0383 /// to get a handle to) which is only implemented for container-like objects
0384 /// like std::vector but not for edm::HepMCProduct!
0385 
0386 void TrackingTruthAccumulator::accumulate(edm::Event const &event, edm::EventSetup const &setup) {
0387   // Call the templated version that does the same for both signal and pileup
0388   // events
0389 
0390   edm::Handle<edm::HepMCProduct> hepmc;
0391   event.getByLabel(hepMCproductLabel_, hepmc);
0392 
0393   accumulateEvent(event, setup, hepmc);
0394 }
0395 
0396 void TrackingTruthAccumulator::accumulate(PileUpEventPrincipal const &event,
0397                                           edm::EventSetup const &setup,
0398                                           edm::StreamID const &) {
0399   // If this bunch crossing is outside the user configured limit, don't do
0400   // anything.
0401   if (event.bunchCrossing() >= -static_cast<int>(maximumPreviousBunchCrossing_) &&
0402       event.bunchCrossing() <= static_cast<int>(maximumSubsequentBunchCrossing_)) {
0403     // edm::LogInfo(messageCategory_) << "Analysing pileup event for bunch
0404     // crossing " << event.bunchCrossing();
0405 
0406     // simply create empty handle as we do not have a HepMCProduct in PU anyway
0407     edm::Handle<edm::HepMCProduct> hepmc;
0408     accumulateEvent(event, setup, hepmc);
0409   } else
0410     edm::LogInfo(messageCategory_) << "Skipping pileup event for bunch crossing " << event.bunchCrossing();
0411 }
0412 
0413 void TrackingTruthAccumulator::finalizeEvent(edm::Event &event, edm::EventSetup const &setup) {
0414   if (createUnmergedCollection_) {
0415     edm::LogInfo("TrackingTruthAccumulator")
0416         << "Adding " << unmergedOutput_.pTrackingParticles->size() << " TrackingParticles and "
0417         << unmergedOutput_.pTrackingVertices->size() << " TrackingVertexs to the event.";
0418 
0419     event.put(std::move(unmergedOutput_.pTrackingParticles));
0420     event.put(std::move(unmergedOutput_.pTrackingVertices));
0421   }
0422 
0423   if (createMergedCollection_) {
0424     edm::LogInfo("TrackingTruthAccumulator")
0425         << "Adding " << mergedOutput_.pTrackingParticles->size() << " merged TrackingParticles and "
0426         << mergedOutput_.pTrackingVertices->size() << " merged TrackingVertexs to the event.";
0427 
0428     event.put(std::move(mergedOutput_.pTrackingParticles), "MergedTrackTruth");
0429     event.put(std::move(mergedOutput_.pTrackingVertices), "MergedTrackTruth");
0430   }
0431 
0432   if (createInitialVertexCollection_) {
0433     edm::LogInfo("TrackingTruthAccumulator")
0434         << "Adding " << pInitialVertices_->size() << " initial TrackingVertexs to the event.";
0435 
0436     event.put(std::move(pInitialVertices_), "InitialVertices");
0437   }
0438 }
0439 
0440 template <class T>
0441 void TrackingTruthAccumulator::accumulateEvent(const T &event,
0442                                                const edm::EventSetup &setup,
0443                                                const edm::Handle<edm::HepMCProduct> &hepMCproduct) {
0444   //
0445   // Get the collections
0446   //
0447   edm::Handle<std::vector<SimTrack>> hSimTracks;
0448   edm::Handle<std::vector<SimVertex>> hSimVertices;
0449   edm::Handle<std::vector<reco::GenParticle>> hGenParticles;
0450   edm::Handle<std::vector<int>> hGenParticleIndices;
0451 
0452   event.getByLabel(simTrackLabel_, hSimTracks);
0453   event.getByLabel(simVertexLabel_, hSimVertices);
0454 
0455   try {
0456     event.getByLabel(genParticleLabel_, hGenParticles);
0457     event.getByLabel(genParticleLabel_, hGenParticleIndices);
0458   } catch (cms::Exception &exception) {
0459     //
0460     // The Monte Carlo is not always available, e.g. for pileup events. The
0461     // information is only used if it's available, but for some reason the
0462     // PileUpEventPrincipal wrapper throws an exception here rather than waiting
0463     // to see if the handle is used (as is the case for edm::Event). So I just
0464     // want to catch this exception and use the normal handle checking later on.
0465     //
0466   }
0467 
0468   // Retrieve tracker topology from geometry
0469   const TrackerTopology *const tTopo = &setup.getData(tTopoToken_);
0470 
0471   // Run through the collections and work out the decay chain of each
0472   // track/vertex. The information in SimTrack and SimVertex only allows
0473   // traversing upwards, but this will allow traversal in both directions. This
0474   // is required for things like grouping electrons that bremsstrahlung as one
0475   // TrackingParticle if "mergedBremsstrahlung" is set in the config file.
0476   DecayChain decayChain(*hSimTracks, *hSimVertices);
0477 
0478   // I only want to create these collections if they're actually required
0479   std::unique_ptr<::OutputCollectionWrapper> pUnmergedCollectionWrapper;
0480   std::unique_ptr<::OutputCollectionWrapper> pMergedCollectionWrapper;
0481   if (createUnmergedCollection_)
0482     pUnmergedCollectionWrapper = std::make_unique<::OutputCollectionWrapper>(decayChain, unmergedOutput_);
0483   if (createMergedCollection_)
0484     pMergedCollectionWrapper = std::make_unique<::OutputCollectionWrapper>(decayChain, mergedOutput_);
0485 
0486   std::vector<const PSimHit *> simHitPointers;
0487   fillSimHits(simHitPointers, event, setup);
0488   TrackingParticleFactory objectFactory(decayChain,
0489                                         hGenParticles,
0490                                         hepMCproduct,
0491                                         hGenParticleIndices,
0492                                         simHitPointers,
0493                                         volumeRadius_,
0494                                         volumeZ_,
0495                                         vertexDistanceCut_,
0496                                         allowDifferentProcessTypeForDifferentDetectors_);
0497 
0498 #if defined(DO_DEBUG_TESTING)
0499   // While I'm testing, perform some checks.
0500   // TODO - drop this call once I'm happy it works in all situations.
0501   decayChain.integrityCheck();
0502 #endif
0503 
0504   TrackingParticleSelector *pSelector = nullptr;
0505   if (selectorFlag_)
0506     pSelector = &selector_;
0507 
0508   // Run over all of the SimTracks, but because I'm interested in the decay
0509   // hierarchy do it through the DecayChainTrack objects. These are looped over
0510   // in sequence here but they have the hierarchy information for the functions
0511   // called to traverse the decay chain.
0512 
0513   for (size_t index = 0; index < decayChain.decayTracksSize; ++index) {
0514     ::DecayChainTrack *pDecayTrack = &decayChain.decayTracks[index];
0515     const SimTrack &simTrack = hSimTracks->at(pDecayTrack->simTrackIndex);
0516 
0517     // Perform some quick checks to see if we can drop out early. Note that
0518     // these are a subset of the cuts in the selector_ so the created
0519     // TrackingParticle could still fail. The selector_ requires the full
0520     // TrackingParticle to be made however, which can be computationally
0521     // expensive.
0522     if (chargedOnly_ && simTrack.charge() == 0)
0523       continue;
0524     if (signalOnly_ && (simTrack.eventId().bunchCrossing() != 0 || simTrack.eventId().event() != 0))
0525       continue;
0526 
0527     // Also perform a check to see if the production vertex is inside the
0528     // tracker volume (if required).
0529     if (ignoreTracksOutsideVolume_) {
0530       const SimVertex &simVertex = hSimVertices->at(pDecayTrack->pParentVertex->simVertexIndex);
0531       if (!objectFactory.vectorIsInsideVolume(simVertex.position()))
0532         continue;
0533     }
0534 
0535     // This function creates the TrackinParticle and adds it to the collection
0536     // if it passes the selection criteria specified in the configuration. If
0537     // the config specifies adding ancestors, the function is called recursively
0538     // to do that.
0539     ::addTrack(pDecayTrack,
0540                pSelector,
0541                pUnmergedCollectionWrapper.get(),
0542                pMergedCollectionWrapper.get(),
0543                objectFactory,
0544                addAncestors_,
0545                tTopo);
0546   }
0547 
0548   // If configured to create a collection of initial vertices, add them from
0549   // this bunch crossing. No selection is applied on this collection, but it
0550   // also has no links to the TrackingParticle decay products. There are a lot
0551   // of "initial vertices", I'm not entirely sure what they all are (nuclear
0552   // interactions in the detector maybe?), but the one for the main event is the
0553   // one with vertexId==0.
0554   if (createInitialVertexCollection_) {
0555     // Pretty sure the one with vertexId==0 is always the first one, but doesn't
0556     // hurt to check
0557     for (const auto &pRootVertex : decayChain.rootVertices) {
0558       const SimVertex &vertex = hSimVertices->at(decayChain.rootVertices[0]->simVertexIndex);
0559       if (vertex.vertexId() != 0)
0560         continue;
0561 
0562       pInitialVertices_->push_back(objectFactory.createTrackingVertex(pRootVertex));
0563       break;
0564     }
0565   }
0566 }
0567 
0568 template <class T>
0569 void TrackingTruthAccumulator::fillSimHits(std::vector<const PSimHit *> &returnValue,
0570                                            const T &event,
0571                                            const edm::EventSetup &setup) {
0572   // loop over the collections
0573   for (const auto &collectionTag : collectionTags_) {
0574     edm::Handle<std::vector<PSimHit>> hSimHits;
0575     event.getByLabel(collectionTag, hSimHits);
0576 
0577     // TODO - implement removing the dead modules
0578     for (const auto &simHit : *hSimHits) {
0579       returnValue.push_back(&simHit);
0580     }
0581 
0582   }  // end of loop over InputTags
0583 
0584   // sort the SimHits according to their time of flight,
0585   // necessary for looping over them "in order" in
0586   // TrackingParticleFactory::createTrackingParticle()
0587   std::sort(returnValue.begin(), returnValue.end(), [](const PSimHit *a, const PSimHit *b) {
0588     const auto atof =
0589         edm::isFinite(a->timeOfFlight()) ? a->timeOfFlight() : std::numeric_limits<decltype(a->timeOfFlight())>::max();
0590     const auto btof =
0591         edm::isFinite(b->timeOfFlight()) ? b->timeOfFlight() : std::numeric_limits<decltype(b->timeOfFlight())>::max();
0592     return atof < btof;
0593   });
0594 }
0595 
0596 //---------------------------------------------------------------------------------
0597 //---------------------------------------------------------------------------------
0598 //---------------------------------------------------------------------------------
0599 //------ ------
0600 //------  End of the definitions for the TrackingTruthAccumulator methods.
0601 //------
0602 //------  Definitions for the functions and classes declared in the unnamed
0603 //------
0604 //------  namespace are below. Documentation for those is above, where ------
0605 //------  they're declared. ------
0606 //------ ------
0607 //---------------------------------------------------------------------------------
0608 //---------------------------------------------------------------------------------
0609 //---------------------------------------------------------------------------------
0610 
0611 namespace  // Unnamed namespace for things only used in this file
0612 {
0613   //---------------------------------------------------------------------------------
0614   //---------------------------------------------------------------------------------
0615   //----   TrackingParticleFactory methods
0616   //----------------------------------------
0617   //---------------------------------------------------------------------------------
0618   //---------------------------------------------------------------------------------
0619 
0620   ::TrackingParticleFactory::TrackingParticleFactory(const ::DecayChain &decayChain,
0621                                                      const edm::Handle<std::vector<reco::GenParticle>> &hGenParticles,
0622                                                      const edm::Handle<edm::HepMCProduct> &hepMCproduct,
0623                                                      const edm::Handle<std::vector<int>> &hHepMCGenParticleIndices,
0624                                                      const std::vector<const PSimHit *> &simHits,
0625                                                      double volumeRadius,
0626                                                      double volumeZ,
0627                                                      double vertexDistanceCut,
0628                                                      bool allowDifferentProcessTypes)
0629       : decayChain_(decayChain),
0630         hGenParticles_(hGenParticles),
0631         hepMCproduct_(hepMCproduct),
0632         simHits_(simHits),
0633         volumeRadius_(volumeRadius),
0634         volumeZ_(volumeZ),
0635         vertexDistanceCut2_(vertexDistanceCut * vertexDistanceCut),
0636         allowDifferentProcessTypeForDifferentDetectors_(allowDifferentProcessTypes) {
0637     // Need to create a multimap to get from a SimTrackId to all of the hits in
0638     // it. The SimTrackId is an unsigned int.
0639     for (size_t index = 0; index < simHits_.size(); ++index) {
0640       trackIdToHitIndex_.insert(std::make_pair(simHits_[index]->trackId(), index));
0641     }
0642 
0643     if (hHepMCGenParticleIndices.isValid())  // Monte Carlo might not be available
0644                                              // for the pileup events
0645     {
0646       genParticleIndices_.resize(hHepMCGenParticleIndices->size() + 1);
0647 
0648       // What I need is the reverse mapping of this vector. The sizes are already
0649       // equivalent because I set the size in the initialiser list.
0650       for (size_t recoGenParticleIndex = 0; recoGenParticleIndex < hHepMCGenParticleIndices->size();
0651            ++recoGenParticleIndex) {
0652         size_t hepMCGenParticleIndex = (*hHepMCGenParticleIndices)[recoGenParticleIndex];
0653 
0654         // They should be the same size, give or take a fencepost error, so this
0655         // should never happen - but just in case
0656         if (genParticleIndices_.size() <= hepMCGenParticleIndex)
0657           genParticleIndices_.resize(hepMCGenParticleIndex + 1);
0658 
0659         genParticleIndices_[hepMCGenParticleIndex] = recoGenParticleIndex;
0660       }
0661     }
0662   }
0663 
0664   TrackingParticle TrackingParticleFactory::createTrackingParticle(const ::DecayChainTrack *pChainTrack,
0665                                                                    const TrackerTopology *tTopo) const {
0666     typedef math::XYZTLorentzVectorD LorentzVector;
0667     typedef math::XYZPoint Vector;
0668 
0669     const SimTrack &simTrack = decayChain_.getSimTrack(pChainTrack);
0670     const SimVertex &parentSimVertex = decayChain_.getSimVertex(pChainTrack->pParentVertex);
0671 
0672     LorentzVector position(0, 0, 0, 0);
0673     if (!simTrack.noVertex())
0674       position = parentSimVertex.position();
0675 
0676     int pdgId = simTrack.type();
0677 
0678     TrackingParticle returnValue;
0679     // N.B. The sim track is added a few lines below, the parent and decay
0680     // vertices are added in another function later on.
0681 
0682     //
0683     // If there is some valid Monte Carlo for this track, take some information
0684     // from that. Only do so if it is from the signal event however. Not sure why
0685     // but that's what the old code did.
0686     //
0687     if (simTrack.eventId().event() == 0 &&
0688         simTrack.eventId().bunchCrossing() == 0)  // if this is a track in the signal event
0689     {
0690       int hepMCGenParticleIndex = simTrack.genpartIndex();
0691       if (hepMCGenParticleIndex >= 0 && hepMCGenParticleIndex < static_cast<int>(genParticleIndices_.size()) &&
0692           hGenParticles_.isValid()) {
0693         int recoGenParticleIndex = genParticleIndices_[hepMCGenParticleIndex];
0694         reco::GenParticleRef generatorParticleRef(hGenParticles_, recoGenParticleIndex);
0695         pdgId = generatorParticleRef->pdgId();
0696         returnValue.addGenParticle(generatorParticleRef);
0697       }
0698     }
0699 
0700     returnValue.addG4Track(simTrack);
0701 
0702     // I need to run over the sim hits and count up the different types of hits.
0703     // To be honest I don't fully understand this next section of code, I copied
0704     // it from the old TrackingTruthProducer to provide the same results. The
0705     // different types of hits that I need to count are:
0706     size_t matchedHits = 0;          // As far as I can tell, this is the number of tracker layers with hits
0707                                      // on them, i.e. taking account of overlaps.
0708     size_t numberOfHits = 0;         // The total number of hits not taking account of overlaps.
0709     size_t numberOfTrackerHits = 0;  // The number of tracker hits not taking account of overlaps.
0710 
0711     bool init = true;
0712     int processType = 0;
0713     int particleType = 0;
0714     int oldLayer = 0;
0715     int newLayer = 0;
0716     DetId oldDetector;
0717     DetId newDetector;
0718 
0719     // Loop over the SimHits associated to this SimTrack
0720     // in the order defined by time of flight, which is
0721     // probably the best quantity available for going
0722     // through the hits "in order" (ok, most important is
0723     // to get the first hit right because processType and
0724     // particleType are taken from it)
0725     for (auto iHitIndex = trackIdToHitIndex_.lower_bound(simTrack.trackId()),
0726               end = trackIdToHitIndex_.upper_bound(simTrack.trackId());
0727          iHitIndex != end;
0728          ++iHitIndex) {
0729       const auto &pSimHit = simHits_[iHitIndex->second];
0730 
0731       // Skip hits with particle type different from SimTrack pdgId
0732       if (pSimHit->particleType() != pdgId)
0733         continue;
0734 
0735       // Initial condition for consistent simhit selection
0736       if (init) {
0737         processType = pSimHit->processType();
0738         particleType = pSimHit->particleType();
0739         newDetector = DetId(pSimHit->detUnitId());
0740         init = false;
0741       }
0742 
0743       oldDetector = newDetector;
0744       newDetector = DetId(pSimHit->detUnitId());
0745 
0746       // Fast sim seems to have different processType between tracker and muons,
0747       // so if this flag has been set allow the processType to change if the
0748       // detector changes.
0749       if (allowDifferentProcessTypeForDifferentDetectors_ && newDetector.det() != oldDetector.det())
0750         processType = pSimHit->processType();
0751 
0752       // Check for delta and interaction products discards
0753       if (processType == pSimHit->processType() && particleType == pSimHit->particleType()) {
0754         ++numberOfHits;
0755         oldLayer = newLayer;
0756         newLayer = 0;
0757         if (newDetector.det() == DetId::Tracker) {
0758           ++numberOfTrackerHits;
0759 
0760           newLayer = tTopo->layer(newDetector);
0761 
0762           // Count hits using layers for glued detectors
0763           if ((oldLayer != newLayer || (oldLayer == newLayer && oldDetector.subdetId() != newDetector.subdetId())))
0764             ++matchedHits;
0765         }
0766       }
0767     }  // end of loop over the sim hits for this sim track
0768 
0769     returnValue.setNumberOfTrackerLayers(matchedHits);
0770     returnValue.setNumberOfHits(numberOfHits);
0771     returnValue.setNumberOfTrackerHits(numberOfTrackerHits);
0772 
0773     return returnValue;
0774   }
0775 
0776   TrackingVertex TrackingParticleFactory::createTrackingVertex(const ::DecayChainVertex *pChainVertex) const {
0777     typedef math::XYZTLorentzVectorD LorentzVector;
0778     typedef math::XYZPoint Vector;
0779     typedef edm::Ref<edm::HepMCProduct, HepMC::GenVertex> GenVertexRef;
0780 
0781     const SimVertex &simVertex = decayChain_.getSimVertex(pChainVertex);
0782 
0783     bool isInVolume = this->vectorIsInsideVolume(simVertex.position());
0784 
0785     TrackingVertex returnValue(simVertex.position(), isInVolume, simVertex.eventId());
0786 
0787     // add the SimVertex to the TrackingVertex
0788     returnValue.addG4Vertex(simVertex);
0789 
0790     // also add refs to nearby HepMC::GenVertexs; the way this is done (i.e. based
0791     // on the position) is transcribed over from the old TrackingTruthProducer
0792     if (simVertex.eventId().event() == 0 && simVertex.eventId().bunchCrossing() == 0 &&
0793         hepMCproduct_.isValid())  // if this is a track in the signal event
0794     {
0795       const HepMC::GenEvent *genEvent = hepMCproduct_->GetEvent();
0796 
0797       if (genEvent != nullptr) {
0798         Vector tvPosition(returnValue.position().x(), returnValue.position().y(), returnValue.position().z());
0799 
0800         for (HepMC::GenEvent::vertex_const_iterator iGenVertex = genEvent->vertices_begin();
0801              iGenVertex != genEvent->vertices_end();
0802              ++iGenVertex) {
0803           HepMC::ThreeVector rawPosition = (*iGenVertex)->position();
0804 
0805           Vector genPosition(rawPosition.x() * 0.1, rawPosition.y() * 0.1, rawPosition.z() * 0.1);
0806 
0807           auto distance2 = (tvPosition - genPosition).mag2();
0808 
0809           if (distance2 < vertexDistanceCut2_)
0810             returnValue.addGenVertex(GenVertexRef(hepMCproduct_, (*iGenVertex)->barcode()));
0811         }
0812       }
0813     }
0814 
0815     return returnValue;
0816   }
0817 
0818   bool ::TrackingParticleFactory::vectorIsInsideVolume(const math::XYZTLorentzVectorD &vector) const {
0819     return (vector.Pt() < volumeRadius_ && std::abs(vector.z()) < volumeZ_);
0820   }
0821 
0822   //---------------------------------------------------------------------------------
0823   //---------------------------------------------------------------------------------
0824   //----   DecayChain methods
0825   //-----------------------------------------------------
0826   //---------------------------------------------------------------------------------
0827   //---------------------------------------------------------------------------------
0828 
0829   ::DecayChain::DecayChain(const std::vector<SimTrack> &trackCollection, const std::vector<SimVertex> &vertexCollection)
0830       : decayTracksSize(trackCollection.size()),
0831         decayVerticesSize(vertexCollection.size()),
0832         decayTracks_(new DecayChainTrack[decayTracksSize]),
0833         decayVertices_(new DecayChainVertex[decayVerticesSize]),
0834         simTrackCollection_(trackCollection),
0835         simVertexCollection_(vertexCollection),
0836         decayTracks(decayTracks_),  // These const references map onto the actual
0837                                     // members for public const access
0838         decayVertices(decayVertices_),
0839         rootVertices(rootVertices_) {
0840     // I need some maps to be able to get object pointers from the track/vertex ID
0841     std::map<int, ::DecayChainTrack *> trackIdToDecayTrack;
0842     std::map<int, ::DecayChainVertex *> vertexIdToDecayVertex;
0843 
0844     // First create a DecayChainTrack for every SimTrack and make a note of the
0845     // trackIds in the map. Also add a pointer to the daughter list of the parent
0846     // DecayChainVertex, which might include creating the vertex object if it
0847     // doesn't already exist.
0848     size_t decayVertexIndex = 0;  // The index of the next free slot in the DecayChainVertex array.
0849     for (size_t index = 0; index < trackCollection.size(); ++index) {
0850       ::DecayChainTrack *pDecayTrack = &decayTracks_[index];  // Get a pointer for ease of use
0851 
0852       // This is the array index of the SimTrack that this DecayChainTrack
0853       // corresponds to. At the moment this is a 1 to 1 relationship with the
0854       // DecayChainTrack index, but they won't necessarily be accessed through the
0855       // array later so it's still required to store it.
0856       pDecayTrack->simTrackIndex = index;
0857 
0858       trackIdToDecayTrack[trackCollection[index].trackId()] = pDecayTrack;
0859 
0860       int parentVertexIndex = trackCollection[index].vertIndex();
0861       if (parentVertexIndex >= 0) {
0862         // Get the DecayChainVertex corresponding to this SimVertex, or initialise
0863         // it if it hasn't been done already.
0864         ::DecayChainVertex *&pParentVertex = vertexIdToDecayVertex[parentVertexIndex];
0865         if (pParentVertex == nullptr) {
0866           // Note that I'm using a reference, so changing pParentVertex will
0867           // change the entry in the map too.
0868           pParentVertex = &decayVertices_[decayVertexIndex];
0869           ++decayVertexIndex;
0870           pParentVertex->simVertexIndex = parentVertexIndex;
0871         }
0872         pParentVertex->daughterTracks.push_back(pDecayTrack);
0873         pDecayTrack->pParentVertex = pParentVertex;
0874       } else
0875         throw std::runtime_error(
0876             "TrackingTruthAccumulator: Found a track with "
0877             "an invalid parent vertex index.");
0878     }
0879 
0880     // This assert was originally in to check the internal consistency of the
0881     // decay chain. Fast sim pileup seems to have a load of vertices with no
0882     // tracks pointing to them though, so fast sim fails this assert if pileup is
0883     // added. I don't think the problem is vital however, so if the assert is
0884     // taken out these extra vertices are ignored.
0885     // assert( vertexIdToDecayVertex.size()==vertexCollection.size() &&
0886     // vertexCollection.size()==decayVertexIndex );
0887 
0888     // I still need to set DecayChainTrack::daughterVertices and
0889     // DecayChainVertex::pParentTrack. The information to do this comes from
0890     // SimVertex::parentIndex. I couldn't do this before because I need all of the
0891     // DecayChainTracks initialised.
0892     for (auto &decayVertexMapPair : vertexIdToDecayVertex) {
0893       ::DecayChainVertex *pDecayVertex = decayVertexMapPair.second;
0894       int parentTrackIndex = vertexCollection[pDecayVertex->simVertexIndex].parentIndex();
0895       if (parentTrackIndex != -1) {
0896         std::map<int, ::DecayChainTrack *>::iterator iParentTrackMapPair = trackIdToDecayTrack.find(parentTrackIndex);
0897         if (iParentTrackMapPair == trackIdToDecayTrack.end()) {
0898           std::stringstream errorStream;
0899           errorStream << "TrackingTruthAccumulator: Something has gone wrong "
0900                          "with the indexing. Parent track index is "
0901                       << parentTrackIndex << ".";
0902           throw std::runtime_error(errorStream.str());
0903         }
0904 
0905         ::DecayChainTrack *pParentTrackHierarchy = iParentTrackMapPair->second;
0906 
0907         pParentTrackHierarchy->daughterVertices.push_back(pDecayVertex);
0908         pDecayVertex->pParentTrack = pParentTrackHierarchy;
0909       } else
0910         rootVertices_.push_back(pDecayVertex);  // Has no parent so is at the top of the decay chain.
0911     }  // end of loop over the vertexIdToDecayVertex map
0912 
0913     findBrem(trackCollection, vertexCollection);
0914 
0915   }  // end of ::DecayChain constructor
0916 
0917 #if defined(DO_DEBUG_TESTING)
0918   // Function documentation is with the declaration above. This function is only
0919   // used for testing.
0920   void ::DecayChain::integrityCheck() {
0921     //
0922     // Start off checking the tracks
0923     //
0924     for (size_t index = 0; index < decayTracksSize; ++index) {
0925       const auto &decayTrack = decayTracks[index];
0926       //
0927       // Tracks should always have a production vertex
0928       //
0929       if (decayTrack.pParentVertex == NULL)
0930         throw std::runtime_error(
0931             "TrackingTruthAccumulator.cc integrityCheck(): "
0932             "Found DecayChainTrack with no parent vertex.");
0933 
0934       //
0935       // Make sure the parent has this as a child. Also make sure it's only listed
0936       // once.
0937       //
0938       size_t numberOfTimesListed = 0;
0939       for (const auto pSiblingTrack : decayTrack.pParentVertex->daughterTracks) {
0940         if (pSiblingTrack == &decayTrack)
0941           ++numberOfTimesListed;
0942       }
0943       if (numberOfTimesListed != 1)
0944         throw std::runtime_error(
0945             "TrackingTruthAccumulator.cc integrityCheck(): "
0946             "Found DecayChainTrack whose parent does not "
0947             "have it listed once and only once as a child.");
0948 
0949       //
0950       // Make sure all of the children have this listed as a parent.
0951       //
0952       for (const auto pDaughterVertex : decayTrack.daughterVertices) {
0953         if (pDaughterVertex->pParentTrack != &decayTrack)
0954           throw std::runtime_error(
0955               "TrackingTruthAccumulator.cc integrityCheck(): Found "
0956               "DecayChainTrack whose child does not have it listed as the "
0957               "parent.");
0958       }
0959 
0960       //
0961       // Follow the chain up to the root vertex, and make sure that is listed in
0962       // rootVertices
0963       //
0964       const DecayChainVertex *pAncestorVertex = decayTrack.pParentVertex;
0965       while (pAncestorVertex->pParentTrack != NULL) {
0966         if (pAncestorVertex->pParentTrack->pParentVertex == NULL)
0967           throw std::runtime_error(
0968               "TrackingTruthAccumulator.cc integrityCheck(): Found "
0969               "DecayChainTrack with no parent vertex higher in the decay chain.");
0970         pAncestorVertex = pAncestorVertex->pParentTrack->pParentVertex;
0971       }
0972       if (std::find(rootVertices.begin(), rootVertices.end(), pAncestorVertex) == rootVertices.end()) {
0973         throw std::runtime_error(
0974             "TrackingTruthAccumulator.cc integrityCheck(): Found DecayChainTrack "
0975             "whose root vertex is not recorded anywhere.");
0976       }
0977     }  // end of loop over decayTracks
0978 
0979     //
0980     // Now check the vertices
0981     //
0982     for (size_t index = 0; index < decayVerticesSize; ++index) {
0983       const auto &decayVertex = decayVertices[index];
0984 
0985       //
0986       // Make sure this, or a vertex higher in the chain, is in the list of root
0987       // vertices.
0988       //
0989       const DecayChainVertex *pAncestorVertex = &decayVertex;
0990       while (pAncestorVertex->pParentTrack != NULL) {
0991         if (pAncestorVertex->pParentTrack->pParentVertex == NULL)
0992           throw std::runtime_error(
0993               "TrackingTruthAccumulator.cc integrityCheck(): Found "
0994               "DecayChainTrack with no parent vertex higher in the vertex decay "
0995               "chain.");
0996         pAncestorVertex = pAncestorVertex->pParentTrack->pParentVertex;
0997       }
0998       if (std::find(rootVertices.begin(), rootVertices.end(), pAncestorVertex) == rootVertices.end()) {
0999         throw std::runtime_error(
1000             "TrackingTruthAccumulator.cc integrityCheck(): Found DecayChainTrack "
1001             "whose root vertex is not recorded anywhere.");
1002       }
1003 
1004       //
1005       // Make sure the parent has this listed in its daughters once and only once.
1006       //
1007       if (decayVertex.pParentTrack != NULL) {
1008         size_t numberOfTimesListed = 0;
1009         for (const auto pSibling : decayVertex.pParentTrack->daughterVertices) {
1010           if (pSibling == &decayVertex)
1011             ++numberOfTimesListed;
1012         }
1013         if (numberOfTimesListed != 1)
1014           throw std::runtime_error(
1015               "TrackingTruthAccumulator.cc integrityCheck(): Found "
1016               "DecayChainVertex whose parent does not have it listed once and "
1017               "only once as a child.");
1018       }
1019 
1020       //
1021       // Make sure all of the children have this listed as the parent
1022       //
1023       for (const auto pDaughter : decayVertex.daughterTracks) {
1024         if (pDaughter->pParentVertex != &decayVertex)
1025           throw std::runtime_error(
1026               "TrackingTruthAccumulator.cc integrityCheck(): Found "
1027               "DecayChainVertex whose child does not have it listed as the "
1028               "parent.");
1029       }
1030     }  // end of loop over decay vertices
1031 
1032     std::cout << "TrackingTruthAccumulator.cc integrityCheck() completed successfully" << std::endl;
1033   }  // end of ::DecayChain::integrityCheck()
1034 #endif
1035 
1036   void ::DecayChain::findBrem(const std::vector<SimTrack> &trackCollection,
1037                               const std::vector<SimVertex> &vertexCollection) {
1038     for (size_t index = 0; index < decayVerticesSize; ++index) {
1039       auto &vertex = decayVertices_[index];
1040 
1041       // Make sure parent is an electron
1042       if (vertex.pParentTrack == nullptr)
1043         continue;
1044       int parentTrackPDG = trackCollection[vertex.pParentTrack->simTrackIndex].type();
1045       if (std::abs(parentTrackPDG) != 11)
1046         continue;
1047 
1048       size_t numberOfElectrons = 0;
1049       size_t numberOfNonElectronsOrPhotons = 0;
1050       for (auto &pDaughterTrack : vertex.daughterTracks) {
1051         const auto &simTrack = trackCollection[pDaughterTrack->simTrackIndex];
1052         if (simTrack.type() == 11 || simTrack.type() == -11)
1053           ++numberOfElectrons;
1054         else if (simTrack.type() != 22)
1055           ++numberOfNonElectronsOrPhotons;
1056       }
1057       if (numberOfElectrons == 1 && numberOfNonElectronsOrPhotons == 0) {
1058         // This is a valid brem, run through and tell all daughters to use the
1059         // parent as the brem
1060         for (auto &pDaughterTrack : vertex.daughterTracks)
1061           pDaughterTrack->pMergedBremSource = vertex.pParentTrack;
1062         vertex.pMergedBremSource = vertex.pParentTrack->pParentVertex;
1063       }
1064     }
1065   }  // end of ::DecayChain::findBrem()
1066 
1067   //---------------------------------------------------------------------------------
1068   //---------------------------------------------------------------------------------
1069   //----   OutputCollectionWrapper methods
1070   //----------------------------------------
1071   //---------------------------------------------------------------------------------
1072   //---------------------------------------------------------------------------------
1073 
1074   ::OutputCollectionWrapper::OutputCollectionWrapper(const ::DecayChain &decayChain,
1075                                                      TrackingTruthAccumulator::OutputCollections &outputCollections)
1076       : output_(outputCollections),
1077         trackingParticleIndices_(decayChain.decayTracksSize, -1),
1078         trackingVertexIndices_(decayChain.decayVerticesSize, -1) {
1079     // No operation besides the initialiser list
1080   }
1081 
1082   TrackingParticle * ::OutputCollectionWrapper::addTrackingParticle(const ::DecayChainTrack *pDecayTrack,
1083                                                                     const TrackingParticle &trackingParticle) {
1084     if (trackingParticleIndices_[pDecayTrack->simTrackIndex] != -1)
1085       throw std::runtime_error(
1086           "OutputCollectionWrapper::addTrackingParticle - "
1087           "trying to add a particle twice");
1088 
1089     trackingParticleIndices_[pDecayTrack->simTrackIndex] = output_.pTrackingParticles->size();
1090     output_.pTrackingParticles->push_back(trackingParticle);
1091 
1092     // Clear any associations in case there were any beforehand
1093     output_.pTrackingParticles->back().clearDecayVertices();
1094     output_.pTrackingParticles->back().clearParentVertex();
1095 
1096     // Associate to the objects that are already in the output collections
1097     associateToExistingObjects(pDecayTrack);
1098 
1099     return &output_.pTrackingParticles->back();
1100   }
1101 
1102   TrackingVertex * ::OutputCollectionWrapper::addTrackingVertex(const ::DecayChainVertex *pDecayVertex,
1103                                                                 const TrackingVertex &trackingVertex) {
1104     if (trackingVertexIndices_[pDecayVertex->simVertexIndex] != -1)
1105       throw std::runtime_error(
1106           "OutputCollectionWrapper::addTrackingVertex - "
1107           "trying to add a vertex twice");
1108 
1109     trackingVertexIndices_[pDecayVertex->simVertexIndex] = output_.pTrackingVertices->size();
1110     output_.pTrackingVertices->push_back(trackingVertex);
1111 
1112     // Associate the new vertex to any parents or children that are already in the
1113     // output collections
1114     associateToExistingObjects(pDecayVertex);
1115 
1116     return &output_.pTrackingVertices->back();
1117   }
1118 
1119   TrackingParticle * ::OutputCollectionWrapper::getTrackingParticle(const ::DecayChainTrack *pDecayTrack) {
1120     const int index = trackingParticleIndices_[pDecayTrack->simTrackIndex];
1121     if (index == -1)
1122       return nullptr;
1123     else
1124       return &(*output_.pTrackingParticles)[index];
1125   }
1126 
1127   TrackingVertex * ::OutputCollectionWrapper::getTrackingVertex(const ::DecayChainVertex *pDecayVertex) {
1128     const int index = trackingVertexIndices_[pDecayVertex->simVertexIndex];
1129     if (index == -1)
1130       return nullptr;
1131     else
1132       return &(*output_.pTrackingVertices)[index];
1133   }
1134 
1135   TrackingParticleRef OutputCollectionWrapper::getRef(const ::DecayChainTrack *pDecayTrack) {
1136     const int index = trackingParticleIndices_[pDecayTrack->simTrackIndex];
1137     if (index == -1)
1138       throw std::runtime_error(
1139           "OutputCollectionWrapper::getRefTrackingParticle - ref requested for a "
1140           "non existent TrackingParticle");
1141     else
1142       return TrackingParticleRef(output_.refTrackingParticles, index);
1143   }
1144 
1145   TrackingVertexRef OutputCollectionWrapper::getRef(const ::DecayChainVertex *pDecayVertex) {
1146     const int index = trackingVertexIndices_[pDecayVertex->simVertexIndex];
1147     if (index == -1)
1148       throw std::runtime_error(
1149           "OutputCollectionWrapper::getRefTrackingParticle - ref requested for a "
1150           "non existent TrackingVertex");
1151     else
1152       return TrackingVertexRef(output_.refTrackingVertexes, index);
1153   }
1154 
1155   void ::OutputCollectionWrapper::setProxy(const ::DecayChainTrack *pOriginalTrack,
1156                                            const ::DecayChainTrack *pProxyTrack) {
1157     int &index = trackingParticleIndices_[pOriginalTrack->simTrackIndex];
1158     if (index != -1)
1159       throw std::runtime_error(
1160           "OutputCollectionWrapper::setProxy() was called for a TrackingParticle "
1161           "that has already been created");
1162     // Note that index is a reference so this call changes the value in the vector
1163     // too
1164     index = trackingParticleIndices_[pProxyTrack->simTrackIndex];
1165   }
1166 
1167   void ::OutputCollectionWrapper::setProxy(const ::DecayChainVertex *pOriginalVertex,
1168                                            const ::DecayChainVertex *pProxyVertex) {
1169     int &index = trackingVertexIndices_[pOriginalVertex->simVertexIndex];
1170     const int newIndex = trackingVertexIndices_[pProxyVertex->simVertexIndex];
1171 
1172     if (index != -1 && index != newIndex)
1173       throw std::runtime_error(
1174           "OutputCollectionWrapper::setProxy() was called for a TrackingVertex "
1175           "that has already been created");
1176 
1177     // Note that index is a reference so this call changes the value in the vector
1178     // too
1179     index = newIndex;
1180   }
1181 
1182   void ::OutputCollectionWrapper::associateToExistingObjects(const ::DecayChainVertex *pChainVertex) {
1183     // First make sure the DecayChainVertex supplied has been turned into a
1184     // TrackingVertex
1185     TrackingVertex *pTrackingVertex = getTrackingVertex(pChainVertex);
1186     if (pTrackingVertex == nullptr)
1187       throw std::runtime_error("associateToExistingObjects was passed a non existent TrackingVertex");
1188 
1189     //
1190     // Associate to the parent track (if there is one)
1191     //
1192     ::DecayChainTrack *pParentChainTrack = pChainVertex->pParentTrack;
1193     if (pParentChainTrack != nullptr)  // Make sure there is a parent track first
1194     {
1195       // There is a parent track, but it might not have been turned into a
1196       // TrackingParticle yet
1197       TrackingParticle *pParentTrackingParticle = getTrackingParticle(pParentChainTrack);
1198       if (pParentTrackingParticle != nullptr) {
1199         pParentTrackingParticle->addDecayVertex(getRef(pChainVertex));
1200         pTrackingVertex->addParentTrack(getRef(pParentChainTrack));
1201       }
1202       // Don't worry about the "else" case - the parent track might not
1203       // necessarily get added to the output collection at all. A check is done on
1204       // daughter vertices when tracks are added, so the association will be
1205       // picked up then.
1206     }
1207 
1208     //
1209     // A parent TrackingVertex is always associated to a daughter TrackingParticle
1210     // when the TrackingParticle is created. Hence there is no need to check the
1211     // list of daughter tracks.
1212     //
1213   }
1214 
1215   void ::OutputCollectionWrapper::associateToExistingObjects(const ::DecayChainTrack *pChainTrack) {
1216     //
1217     // First make sure this DecayChainTrack has been turned into a
1218     // TrackingParticle
1219     //
1220     TrackingParticle *pTrackingParticle = getTrackingParticle(pChainTrack);
1221     if (pTrackingParticle == nullptr)
1222       throw std::runtime_error(
1223           "associateToExistingObjects was passed a non "
1224           "existent TrackingParticle");
1225 
1226     // Get the parent vertex. This should always already have been turned into a
1227     // TrackingVertex, and there should always be a parent DecayChainVertex.
1228     ::DecayChainVertex *pParentChainVertex = pChainTrack->pParentVertex;
1229     TrackingVertex *pParentTrackingVertex = getTrackingVertex(pParentChainVertex);
1230 
1231     //
1232     // First associate to the parent vertex.
1233     //
1234     pTrackingParticle->setParentVertex(getRef(pParentChainVertex));
1235     pParentTrackingVertex->addDaughterTrack(getRef(pChainTrack));
1236 
1237     // If there are any daughter vertices that have already been made into a
1238     // TrackingVertex make sure they know about each other. If the daughter
1239     // vertices haven't been made into TrackingVertexs yet, I won't do anything
1240     // and create the association when the vertex is made.
1241     for (auto pDaughterChainVertex : pChainTrack->daughterVertices) {
1242       TrackingVertex *pDaughterTrackingVertex = getTrackingVertex(pDaughterChainVertex);
1243       if (pDaughterTrackingVertex != nullptr) {
1244         pTrackingParticle->addDecayVertex(getRef(pDaughterChainVertex));
1245         pDaughterTrackingVertex->addParentTrack(getRef(pChainTrack));
1246       }
1247     }
1248   }
1249 
1250   TrackingParticle *addTrackAndParentVertex(::DecayChainTrack *pDecayTrack,
1251                                             const TrackingParticle &trackingParticle,
1252                                             ::OutputCollectionWrapper *pOutput) {
1253     // See if this TrackingParticle has already been created (could be if the
1254     // DecayChainTracks are looped over in a funny order). If it has then there's
1255     // no need to do anything.
1256     TrackingParticle *pTrackingParticle = pOutput->getTrackingParticle(pDecayTrack);
1257     if (pTrackingParticle == nullptr) {
1258       // Need to make sure the production vertex has been created first
1259       if (pOutput->getTrackingVertex(pDecayTrack->pParentVertex) == nullptr) {
1260         // TrackingVertex doesn't exist in the output collection yet. However,
1261         // it's already been created in the addTrack() function and a temporary
1262         // reference to it set in the TrackingParticle. I'll use that reference to
1263         // create a copy in the output collection. When the TrackingParticle is
1264         // added to the output collection a few lines below the temporary
1265         // reference to the parent vertex will be cleared, and the correct one
1266         // referring to the output collection will be set.
1267         pOutput->addTrackingVertex(pDecayTrack->pParentVertex, *trackingParticle.parentVertex());
1268       }
1269 
1270       pTrackingParticle = pOutput->addTrackingParticle(pDecayTrack, trackingParticle);
1271     }
1272 
1273     return pTrackingParticle;
1274   }
1275 
1276   void addTrack(::DecayChainTrack *pDecayChainTrack,
1277                 const TrackingParticleSelector *pSelector,
1278                 ::OutputCollectionWrapper *pUnmergedOutput,
1279                 ::OutputCollectionWrapper *pMergedOutput,
1280                 const ::TrackingParticleFactory &objectFactory,
1281                 bool addAncestors,
1282                 const TrackerTopology *tTopo) {
1283     if (pDecayChainTrack == nullptr)
1284       return;  // This is required for when the addAncestors_ recursive call
1285                // reaches the top of the chain
1286 
1287     // Check to see if this particle has already been processed while traversing
1288     // up the parents of another split in the decay chain. The check in the line
1289     // above only stops when the top of the chain is reached, whereas this will
1290     // stop when a previously traversed split is reached.
1291     {  // block to limit the scope of temporary variables
1292       bool alreadyProcessed = true;
1293       if (pUnmergedOutput != nullptr) {
1294         if (pUnmergedOutput->getTrackingParticle(pDecayChainTrack) == nullptr)
1295           alreadyProcessed = false;
1296       }
1297       if (pMergedOutput != nullptr) {
1298         if (pMergedOutput->getTrackingParticle(pDecayChainTrack) == nullptr)
1299           alreadyProcessed = false;
1300       }
1301       if (alreadyProcessed)
1302         return;
1303     }
1304 
1305     // Create a TrackingParticle.
1306     TrackingParticle newTrackingParticle = objectFactory.createTrackingParticle(pDecayChainTrack, tTopo);
1307 
1308     // The selector checks the impact parameters from the vertex, so I need to
1309     // have a valid reference to the parent vertex in the TrackingParticle before
1310     // that can be called. TrackingParticle needs an edm::Ref for the parent
1311     // TrackingVertex though. I still don't know if this is going to be added to
1312     // the collection so I can't take it from there, so I need to create a
1313     // temporary one. When the addTrackAndParentVertex() is called (assuming it
1314     // passes selection) it will use the temporary reference to create a copy of
1315     // the parent vertex, put that in the output collection, and then set the
1316     // reference in the TrackingParticle properly.
1317     TrackingVertexCollection dummyCollection;  // Only needed to create an edm::Ref
1318     dummyCollection.push_back(objectFactory.createTrackingVertex(pDecayChainTrack->pParentVertex));
1319     TrackingVertexRef temporaryRef(&dummyCollection, 0);
1320     newTrackingParticle.setParentVertex(temporaryRef);
1321 
1322     // If a selector has been supplied apply it on the new TrackingParticle and
1323     // return if it fails.
1324     if (pSelector) {
1325       if (!(*pSelector)(newTrackingParticle))
1326         return;  // Return if the TrackingParticle fails selection
1327     }
1328 
1329     // Add the ancestors first (if required) so that the collection has some kind
1330     // of chronological order. I don't know how important that is but other code
1331     // might assume chronological order. If adding ancestors, no selection is
1332     // applied. Note that I've already checked that all DecayChainTracks have a
1333     // pParentVertex.
1334     if (addAncestors)
1335       addTrack(pDecayChainTrack->pParentVertex->pParentTrack,
1336                nullptr,
1337                pUnmergedOutput,
1338                pMergedOutput,
1339                objectFactory,
1340                addAncestors,
1341                tTopo);
1342 
1343     // If creation of the unmerged collection has been turned off in the config
1344     // this pointer will be null.
1345     if (pUnmergedOutput != nullptr)
1346       addTrackAndParentVertex(pDecayChainTrack, newTrackingParticle, pUnmergedOutput);
1347 
1348     // If creation of the merged collection has been turned off in the config this
1349     // pointer will be null.
1350     if (pMergedOutput != nullptr) {
1351       ::DecayChainTrack *pBremParentChainTrack = pDecayChainTrack;
1352       while (pBremParentChainTrack->pMergedBremSource != nullptr)
1353         pBremParentChainTrack = pBremParentChainTrack->pMergedBremSource;
1354 
1355       if (pBremParentChainTrack != pDecayChainTrack) {
1356         TrackingParticle *pBremParentTrackingParticle =
1357             addTrackAndParentVertex(pBremParentChainTrack, newTrackingParticle, pMergedOutput);
1358         // The original particle in the bremsstrahlung decay chain has been
1359         // created (or retrieved if it already existed), now copy in the
1360         // extra information.
1361         // TODO - copy extra information.
1362 
1363         if (std::abs(newTrackingParticle.pdgId()) == 22) {
1364           // Photons are added as separate TrackingParticles, but with the
1365           // production vertex changed to be the production vertex of the original
1366           // electron.
1367 
1368           // Set up a proxy, so that all requests for the parent TrackingVertex
1369           // get redirected to the brem parent's TrackingVertex.
1370           pMergedOutput->setProxy(pDecayChainTrack->pParentVertex, pBremParentChainTrack->pParentVertex);
1371 
1372           // Now that pMergedOutput thinks the parent vertex is the brem parent's
1373           // vertex I can call this and it will set the TrackingParticle parent
1374           // vertex correctly to the brem parent vertex.
1375           addTrackAndParentVertex(pDecayChainTrack, newTrackingParticle, pMergedOutput);
1376         } else if (std::abs(newTrackingParticle.pdgId()) == 11) {
1377           // Electrons have their G4 tracks and SimHits copied to the parent
1378           // TrackingParticle.
1379           for (const auto &trackSegment : newTrackingParticle.g4Tracks()) {
1380             pBremParentTrackingParticle->addG4Track(trackSegment);
1381           }
1382 
1383           // Also copy the generator particle references
1384           for (const auto &genParticleRef : newTrackingParticle.genParticles()) {
1385             pBremParentTrackingParticle->addGenParticle(genParticleRef);
1386           }
1387 
1388           pBremParentTrackingParticle->setNumberOfHits(pBremParentTrackingParticle->numberOfHits() +
1389                                                        newTrackingParticle.numberOfHits());
1390           pBremParentTrackingParticle->setNumberOfTrackerHits(pBremParentTrackingParticle->numberOfTrackerHits() +
1391                                                               newTrackingParticle.numberOfTrackerHits());
1392           pBremParentTrackingParticle->setNumberOfTrackerLayers(pBremParentTrackingParticle->numberOfTrackerLayers() +
1393                                                                 newTrackingParticle.numberOfTrackerLayers());
1394 
1395           // Set a proxy in the output collection wrapper so that any attempt to
1396           // get objects for this DecayChainTrack again get redirected to the brem
1397           // parent.
1398           pMergedOutput->setProxy(pDecayChainTrack, pBremParentChainTrack);
1399         }
1400       } else {
1401         // This is not the result of bremsstrahlung so add to the collection as
1402         // normal.
1403         addTrackAndParentVertex(pDecayChainTrack, newTrackingParticle, pMergedOutput);
1404       }
1405     }  // end of "if( pMergedOutput!=NULL )", i.e. end of "if bremsstrahlung
1406     // merging is turned on"
1407 
1408   }  // end of addTrack function
1409 
1410 }  // namespace
1411 
1412 // Register with the framework
1413 DEFINE_DIGI_ACCUMULATOR(TrackingTruthAccumulator);