MtdTruthAccumulator Class Reference

Inheritance diagram for MtdTruthAccumulator:

Classes
struct	calo_particles
struct	OutputCollections

Public Member Functions
	MtdTruthAccumulator (const edm::ParameterSet &config, edm::ProducesCollector, edm::ConsumesCollector &iC)
Public Member Functions inherited from DigiAccumulatorMixMod
virtual void	beginLuminosityBlock (edm::LuminosityBlock const &lumi, edm::EventSetup const &setup)
virtual void	beginRun (edm::Run const &run, edm::EventSetup const &setup)
	DigiAccumulatorMixMod ()
	DigiAccumulatorMixMod (DigiAccumulatorMixMod const &)=delete
virtual void	endLuminosityBlock (edm::LuminosityBlock const &lumi, edm::EventSetup const &setup)
virtual void	endRun (edm::Run const &run, edm::EventSetup const &setup)
virtual void	finalizeBunchCrossing (edm::Event &event, edm::EventSetup const &setup, int bunchCrossing)
virtual PileupMixingContent *	getEventPileupInfo ()
virtual void	initializeBunchCrossing (edm::Event const &event, edm::EventSetup const &setup, int bunchCrossing)
DigiAccumulatorMixMod const &	operator= (DigiAccumulatorMixMod const &)=delete
virtual void	StorePileupInformation (std::vector< int > &numInteractionList, std::vector< int > &bunchCrossingList, std::vector< float > &TrueInteractionList, std::vector< edm::EventID > &eventList, int bunchSpace)
virtual	~DigiAccumulatorMixMod ()

Private Member Functions
void	accumulate (const edm::Event &event, const edm::EventSetup &setup) override
void	accumulate (const PileUpEventPrincipal &event, const edm::EventSetup &setup, edm::StreamID const &) override
template<class T >
void	accumulateEvent (const T &event, const edm::EventSetup &setup, const edm::Handle< edm::HepMCProduct > &hepMCproduct)
	Both forms of accumulate() delegate to this templated method. More...
template<class T >
void	fillSimHits (std::vector< std::pair< uint64_t, const PSimHit *>> &returnValue, std::unordered_map< int, std::map< uint64_t, std::tuple< float, float, LocalPoint >>> &simTrackDetIdMap, const T &event, const edm::EventSetup &setup)
	Fills the supplied vector with pointers to the SimHits, checking for bad modules if required. More...
void	finalizeEvent (edm::Event &event, const edm::EventSetup &setup) override
void	initializeEvent (const edm::Event &event, const edm::EventSetup &setup) override

Private Attributes
const unsigned int	bunchSpacing_
std::vector< edm::InputTag >	collectionTags_
edm::InputTag	genParticleLabel_
const MTDGeometry *	geom = nullptr
const edm::ESGetToken< MTDGeometry, MTDDigiGeometryRecord >	geomToken_
mtd::MTDGeomUtil	geomTools_
edm::InputTag	hepMCproductLabel_
	Needed to add HepMC::GenVertex to SimVertex. More...
edm::Handle< std::vector< SimTrack > >	hSimTracks
edm::Handle< std::vector< SimVertex > >	hSimVertices
bool	isEtl_
calo_particles	m_caloParticles
std::unordered_map< uint64_t, float >	m_detIdToTotalSimEnergy
std::unordered_multimap< Barcode_t, Index_t >	m_simHitBarcodeToIndex
const unsigned int	maximumPreviousBunchCrossing_
const unsigned int	maximumSubsequentBunchCrossing_
const double	maxPseudoRapidity_
const std::string	messageCategory_
const double	minEnergy_
const edm::ESGetToken< MTDTopology, MTDTopologyRcd >	mtdtopoToken_
OutputCollections	output_
const bool	premixStage1_
const edm::InputTag	simTrackLabel_
const edm::InputTag	simVertexLabel_
const MTDTopology *	topology = nullptr

Detailed Description

Definition at line 67 of file MtdTruthAccumulator.cc.

Constructor & Destructor Documentation

MtdTruthAccumulator()

MtdTruthAccumulator::MtdTruthAccumulator ( const edm::ParameterSet &  config, edm::ProducesCollector  producesCollector, edm::ConsumesCollector &  iC  )

explicit

Definition at line 254 of file MtdTruthAccumulator.cc.

References collectionTags_, edm::ConsumesCollector::consumes(), genParticleLabel_, hepMCproductLabel_, isEtl_, premixStage1_, edm::ProducesCollector::produces(), simTrackLabel_, simVertexLabel_, and triggerMatcherToHLTDebug_cfi::tags.

    : messageCategory_("MtdTruthAccumulator"),
      maximumPreviousBunchCrossing_(config.getParameter<unsigned int>("maximumPreviousBunchCrossing")),
      maximumSubsequentBunchCrossing_(config.getParameter<unsigned int>("maximumSubsequentBunchCrossing")),
      bunchSpacing_(config.getParameter<unsigned int>("bunchspace")),
      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")),
      geomToken_(iC.esConsumes<MTDGeometry, MTDDigiGeometryRecord>()),
      mtdtopoToken_(iC.esConsumes<MTDTopology, MTDTopologyRcd>()),
      minEnergy_(config.getParameter<double>("MinEnergy")),
      maxPseudoRapidity_(config.getParameter<double>("MaxPseudoRapidity")),
      premixStage1_(config.getParameter<bool>("premixStage1")) {
  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 (auto const &parameterName : parameterNames) {
    std::vector<edm::InputTag> tags = simHitCollectionConfig.getParameter<std::vector<edm::InputTag>>(parameterName);
    collectionTags_.insert(collectionTags_.end(), tags.begin(), tags.end());
  }

  for (auto const &collectionTag : collectionTags_) {
    iC.consumes<std::vector<PSimHit>>(collectionTag);
    isEtl_ = (collectionTag.instance().find("FastTimerHitsEndcap") != std::string::npos);
  }

  producesCollector.produces<MtdSimClusterCollection>("MergedMtdTruth");
  producesCollector.produces<MtdSimLayerClusterCollection>("MergedMtdTruthLC");
  producesCollector.produces<MtdSimTracksterCollection>("MergedMtdTruthST");
  producesCollector.produces<MtdCaloParticleCollection>("MergedMtdTruth");
  if (premixStage1_) {
    producesCollector.produces<std::vector<std::pair<uint64_t, float>>>("MergedMtdTruth");
  }
}

Member Function Documentation

accumulate() [1/2]

void MtdTruthAccumulator::accumulate ( const edm::Event &  event, const edm::EventSetup &  setup  )

overrideprivatevirtual

Create handle to edm::HepMCProduct here because event.getByLabel with edm::HepMCProduct only works for edm::Event but not for PileUpEventPrincipal; PileUpEventPrincipal::getByLabel tries to call T::value_type and T::iterator (where T is the type of the object one wants to get a handle to) which is only implemented for container-like objects like std::vector but not for edm::HepMCProduct!

Implements DigiAccumulatorMixMod.

Definition at line 326 of file MtdTruthAccumulator.cc.

References accumulateEvent(), hepMCproductLabel_, messageCategory_, and singleTopDQM_cfi::setup.

                                                                                      {
  edm::Handle<edm::HepMCProduct> hepmc;
  event.getByLabel(hepMCproductLabel_, hepmc);

  edm::LogInfo(messageCategory_) << " MtdTruthAccumulator::accumulate (signal)";
  accumulateEvent(event, setup, hepmc);
}

accumulate() [2/2]

void MtdTruthAccumulator::accumulate ( const PileUpEventPrincipal &  event, const edm::EventSetup &  setup, edm::StreamID const &    )

overrideprivatevirtual

Implements DigiAccumulatorMixMod.

Definition at line 334 of file MtdTruthAccumulator.cc.

References accumulateEvent(), maximumPreviousBunchCrossing_, maximumSubsequentBunchCrossing_, messageCategory_, and singleTopDQM_cfi::setup.

                                                          {
  if (event.bunchCrossing() >= -static_cast<int>(maximumPreviousBunchCrossing_) &&
      event.bunchCrossing() <= static_cast<int>(maximumSubsequentBunchCrossing_)) {
    // simply create empty handle as we do not have a HepMCProduct in PU anyway
    edm::Handle<edm::HepMCProduct> hepmc;
    edm::LogInfo(messageCategory_) << " MtdTruthAccumulator::accumulate (pileup) bunchCrossing="
                                   << event.bunchCrossing();
    accumulateEvent(event, setup, hepmc);
  } else {
    edm::LogInfo(messageCategory_) << "Skipping pileup event for bunch crossing " << event.bunchCrossing();
  }
}

accumulateEvent()

template<class T >

void MtdTruthAccumulator::accumulateEvent ( const T &  event, const edm::EventSetup &  setup, const edm::Handle< edm::HepMCProduct > &  hepMCproduct  )

private

Both forms of accumulate() delegate to this templated method.

Build the main decay graph and assign the SimTrack to each edge. The graph built here will only contain the particles that have a decay vertex associated to them. In order to recover also the particles that will not decay, we need to keep track of the SimTrack used here and add, a-posteriori, the ones not used, associating a ghost vertex (starting from the highest simulated vertex number), in order to build the edge and identify them immediately as stable (i.e. not decayed).

To take into account the multi-bremsstrahlung effects in which a single particle is emitting photons in different vertices keeping the same track index, we also collapsed those vertices into 1 unique vertex. The other approach of fully representing the decay chain keeping the same track index would have the problem of over-counting the contributions of that track, especially in terms of hits.

The 2 auxiliary vectors are structured as follow:

1. used_sim_tracks is a vector that has the same size as the overall number of simulated tracks. The associated integer is the vertexId of the decaying vertex for that track.
2. collapsed_vertices is a vector that has the same size as the overall number of simulated vertices. The vector's index is the vertexId itself, the associated value is the vertexId of the vertex on which this should collapse.

Definition at line 593 of file MtdTruthAccumulator.cc.

References funct::abs(), bunchSpacing_, gather_cfg::cout, EdgeProperty::cumulative_simHits, DEBUG, PA_ZEESkim_cff::decay, fillSimHits(), genParticleLabel_, hSimTracks, hSimVertices, mps_fire::i, heavyIonCSV_trainingSettings::idx, IfLogDebug, m_caloParticles, m_simHitBarcodeToIndex, maxPseudoRapidity_, messageCategory_, minEnergy_, CoreSimTrack::momentum(), SimTrack::noGenpart(), hltrates_dqm_sourceclient-live_cfg::offset, output_, put(), edm::second(), singleTopDQM_cfi::setup, EdgeProperty::simTrack, simTrackLabel_, simVertexLabel_, findQualityFiles::size, submitPVValidationJobs::t, DiMuonV_cfg::tracks, findQualityFiles::v, AlignmentTracksFromVertexSelector_cfi::vertices, and runTauDisplay::vis.

Referenced by accumulate().

                                                                                          {
  edm::Handle<std::vector<reco::GenParticle>> hGenParticles;
  edm::Handle<std::vector<int>> hGenParticleIndices;

  event.getByLabel(simTrackLabel_, hSimTracks);
  event.getByLabel(simVertexLabel_, hSimVertices);

  event.getByLabel(genParticleLabel_, hGenParticles);
  event.getByLabel(genParticleLabel_, hGenParticleIndices);

  std::vector<std::pair<uint64_t, const PSimHit *>> simHitPointers;
  std::unordered_map<int, std::map<uint64_t, std::tuple<float, float, LocalPoint>>> simTrackDetIdMap;
  fillSimHits(simHitPointers, simTrackDetIdMap, event, setup);

  // Clear maps from previous event fill them for this one
  m_simHitBarcodeToIndex.clear();
  for (unsigned int i = 0; i < simHitPointers.size(); ++i) {
    m_simHitBarcodeToIndex.emplace(simHitPointers[i].second->trackId(), i);
  }

  auto const &tracks = *hSimTracks;
  auto const &vertices = *hSimVertices;
  std::unordered_map<int, int> trackid_to_track_index;
  std::unordered_map<uint32_t, float> vertex_time_map;
  DecayChain decay;

  for (uint32_t i = 0; i < vertices.size(); i++) {
    vertex_time_map[i] = vertices[i].position().t() * 1e9 + event.bunchCrossing() * static_cast<int>(bunchSpacing_);
  }

  IfLogDebug(DEBUG, messageCategory_) << " TRACKS" << std::endl;
  int idx = 0;
  for (auto const &t : tracks) {
    IfLogDebug(DEBUG, messageCategory_) << " " << idx << "\t" << t.trackId() << "\t" << t << std::endl;
    trackid_to_track_index[t.trackId()] = idx;
    idx++;
  }

  idx = 0;
  std::vector<int> used_sim_tracks(tracks.size(), 0);
  std::vector<int> collapsed_vertices(vertices.size(), 0);
  IfLogDebug(DEBUG, messageCategory_) << " VERTICES" << std::endl;
  for (auto const &v : vertices) {
    IfLogDebug(DEBUG, messageCategory_) << " " << idx++ << "\t" << v << std::endl;
    if (v.parentIndex() != -1) {
      auto const trk_idx = trackid_to_track_index[v.parentIndex()];
      auto origin_vtx = tracks[trk_idx].vertIndex();
      if (used_sim_tracks[trk_idx]) {
        // collapse the vertex into the original first vertex we saw associated
        // to this track. Omit adding the edge in order to avoid double
        // counting of the very same particles  and its associated hits.
        collapsed_vertices[v.vertexId()] = used_sim_tracks[trk_idx];
        continue;
      }
      // Perform the actual vertex collapsing, if needed.
      if (collapsed_vertices[origin_vtx])
        origin_vtx = collapsed_vertices[origin_vtx];
      add_edge(
          origin_vtx, v.vertexId(), EdgeProperty(&tracks[trk_idx], simTrackDetIdMap[v.parentIndex()].size(), 0), decay);
      used_sim_tracks[trk_idx] = v.vertexId();
    }
  }
  // Build the motherParticle property to each vertex
  auto const &vertexMothersProp = get(vertex_name, decay);
  // Now recover the particles that did not decay. Append them with an index
  // bigger than the size of the generated vertices.
  int offset = vertices.size();
  for (size_t i = 0; i < tracks.size(); ++i) {
    if (!used_sim_tracks[i]) {
      auto origin_vtx = tracks[i].vertIndex();
      // Perform the actual vertex collapsing, if needed.
      if (collapsed_vertices[origin_vtx])
        origin_vtx = collapsed_vertices[origin_vtx];
      add_edge(origin_vtx, offset, EdgeProperty(&tracks[i], simTrackDetIdMap[tracks[i].trackId()].size(), 0), decay);
      // The properties for "fake" vertices associated to stable particles have
      // to be set inside this loop, since they do not belong to the vertices
      // collection and would be skipped by that loop (coming next)
      put(vertexMothersProp, offset, VertexProperty(&tracks[i], 0));
      offset++;
    }
  }
  for (auto const &v : vertices) {
    if (v.parentIndex() != -1) {
      // Skip collapsed_vertices
      if (collapsed_vertices[v.vertexId()])
        continue;
      put(vertexMothersProp, v.vertexId(), VertexProperty(&tracks[trackid_to_track_index[v.parentIndex()]], 0));
    }
  }
  SimHitsAccumulator_dfs_visitor vis;
  depth_first_search(decay, visitor(vis));
  CaloParticle_dfs_visitor caloParticleCreator(
      output_,
      m_caloParticles,
      m_simHitBarcodeToIndex,
      simTrackDetIdMap,
      vertex_time_map,
      [&](EdgeProperty &edge_property) -> bool {
        // Apply selection on SimTracks in order to promote them to be
        // CaloParticles. The function returns TRUE if the particle satisfies
        // the selection, FALSE otherwise. Therefore the correct logic to select
        // the particle is to ask for TRUE as return value.
        return (edge_property.cumulative_simHits != 0 and !edge_property.simTrack->noGenpart() and
                edge_property.simTrack->momentum().E() > minEnergy_ and
                std::abs(edge_property.simTrack->momentum().Eta()) < maxPseudoRapidity_);
      });
  depth_first_search(decay, visitor(caloParticleCreator));

#if DEBUG
  boost::write_graphviz(std::cout,
                        decay,
                        make_label_writer(make_transform_value_property_map(&graphviz_vertex, get(vertex_name, decay))),
                        make_label_writer(make_transform_value_property_map(&graphviz_edge, get(edge_weight, decay))));
#endif
}

fillSimHits()

template<class T >

void MtdTruthAccumulator::fillSimHits ( std::vector< std::pair< uint64_t, const PSimHit *>> &  returnValue, std::unordered_map< int, std::map< uint64_t, std::tuple< float, float, LocalPoint >>> &  simTrackDetIdMap, const T &  event, const edm::EventSetup &  setup  )

private

Fills the supplied vector with pointers to the SimHits, checking for bad modules if required.

Definition at line 738 of file MtdTruthAccumulator.cc.

References collectionTags_, angle_units::operators::convertMmToCm(), geant_units::operators::convertUnitsTo(), DEBUG, geomTools_, mtd::MTDGeomUtil::globalPosition(), l1ctLayer2EG_cff::id, IfLogDebug, createfilelist::int, mtd::MTDGeomUtil::layer(), m_detIdToTotalSimEnergy, messageCategory_, mtd::MTDGeomUtil::module(), muonClassificationByHits_cfi::pixel, mtd::MTDGeomUtil::pixelInModule(), position, rpcPointValidation_cfi::simHit, x, PV3DBase< T, PVType, FrameType >::x(), y, PV3DBase< T, PVType, FrameType >::y(), z, and PV3DBase< T, PVType, FrameType >::z().

Referenced by accumulateEvent().

                                {
  using namespace geant_units::operators;
  using namespace angle_units::operators;
  for (auto const &collectionTag : collectionTags_) {
    edm::Handle<std::vector<PSimHit>> hSimHits;
    event.getByLabel(collectionTag, hSimHits);

    for (auto const &simHit : *hSimHits) {
      DetId id(0);

      // --- Use only hits compatible with the in-time bunch-crossing
      if (simHit.tof() < 0 || simHit.tof() > 25.)
        continue;

      id = simHit.detUnitId();

      if (id == DetId(0)) {
        edm::LogWarning(messageCategory_) << "Invalid DetId for the current simHit!";
        continue;
      }

      if (simHit.trackId() == 0) {
        continue;
      }

      returnValue.emplace_back(id, &simHit);

      // get an unique id: for BTL the detId is unique (one for each crystal), for ETL the detId is not enough
      // also row and column are needed. An unique number is created from detId, row, col
      // Get row and column
      const auto &position = simHit.localPosition();

      LocalPoint simscaled(convertMmToCm(position.x()), convertMmToCm(position.y()), convertMmToCm(position.z()));
      std::pair<uint8_t, uint8_t> pixel = geomTools_.pixelInModule(id, simscaled);
      // create the unique id
      uint64_t uniqueId = static_cast<uint64_t>(id.rawId()) << 32;
      uniqueId |= pixel.first << 16;
      uniqueId |= pixel.second;

      std::get<0>(simTrackDetIdMap[simHit.trackId()][uniqueId]) += simHit.energyLoss();
      m_detIdToTotalSimEnergy[uniqueId] += simHit.energyLoss();
      // --- Get the time of the first SIM hit in the cell
      if (std::get<1>(simTrackDetIdMap[simHit.trackId()][uniqueId]) == 0. ||
          simHit.tof() < std::get<1>(simTrackDetIdMap[simHit.trackId()][uniqueId])) {
        std::get<1>(simTrackDetIdMap[simHit.trackId()][uniqueId]) = simHit.tof();
      }

      float xSim = std::get<2>(simTrackDetIdMap[simHit.trackId()][uniqueId]).x() + simscaled.x() * simHit.energyLoss();
      float ySim = std::get<2>(simTrackDetIdMap[simHit.trackId()][uniqueId]).y() + simscaled.y() * simHit.energyLoss();
      float zSim = std::get<2>(simTrackDetIdMap[simHit.trackId()][uniqueId]).z() + simscaled.z() * simHit.energyLoss();
      LocalPoint posSim(xSim, ySim, zSim);
      std::get<2>(simTrackDetIdMap[simHit.trackId()][uniqueId]) = posSim;

#ifdef PRINT_DEBUG
      IfLogDebug(DEBUG, messageCategory_)
          << "hitId " << id.rawId() << " from track " << simHit.trackId() << " in layer " << geomTools_.layer(id)
          << ", module " << geomTools_.module(id) << ", pixel ( " << (int)geomTools_.pixelInModule(id, simscaled).first
          << ", " << (int)geomTools_.pixelInModule(id, simscaled).second << " )\n global pos(cm) "
          << geomTools_.globalPosition(id, simscaled) << ", time(ns) " << simHit.tof() << ", energy(MeV) "
          << convertUnitsTo(0.001_MeV, simHit.energyLoss()) << std::endl;
#endif
    }  // end of loop over simHits

  }  // End of loop over InputTags

  for (auto &tkIt : simTrackDetIdMap) {
    for (auto &uIt : tkIt.second) {
      float accEnergy = std::get<0>(uIt.second);
      float xSim = std::get<2>(uIt.second).x() / accEnergy;
      float ySim = std::get<2>(uIt.second).y() / accEnergy;
      float zSim = std::get<2>(uIt.second).z() / accEnergy;
      LocalPoint posSim(xSim, ySim, zSim);
      std::get<2>(uIt.second) = posSim;
    }
  }
}

finalizeEvent()

void MtdTruthAccumulator::finalizeEvent ( edm::Event &  event, const edm::EventSetup &  setup  )

overrideprivatevirtual

Implements DigiAccumulatorMixMod.

Definition at line 349 of file MtdTruthAccumulator.cc.

References a, b, geant_units::operators::convertUnitsTo(), filterCSVwithJSON::copy, MTDTopologyMode::crysLayoutFromTopoMode(), DEBUG, hcalRecHitTable_cff::detId, hcalRecHitTable_cff::energy, dqmdumpme::first, alignBH_cfg::fixed, HLT_2024v10_cff::fraction, geom, geomTools_, MTDTopology::getMTDTopologyMode(), mps_fire::i, l1ctLayer2EG_cff::id, MTDGeometry::idToDet(), IfLogDebug, dqmdumpme::indices, createfilelist::int, mtd::MTDGeomUtil::isBTL(), mtd::MTDGeomUtil::isETL(), dqmiolumiharvest::j, mtd::MTDGeomUtil::layer(), m_caloParticles, m_detIdToTotalSimEnergy, m_simHitBarcodeToIndex, messageCategory_, mtd::MTDGeomUtil::module(), eostools::move(), or, output_, MtdTruthAccumulator::OutputCollections::pCaloParticles, RectangularMTDTopology::pixelToModuleLocalPoint(), MtdTruthAccumulator::OutputCollections::pMtdSimLayerClusters, MtdTruthAccumulator::OutputCollections::pMtdSimTracksters, mtd::MTDGeomUtil::position(), position, premixStage1_, MtdTruthAccumulator::OutputCollections::pSimClusters, MtdTruthAccumulator::calo_particles::sc_start_, MtdTruthAccumulator::calo_particles::sc_stop_, edm::second(), jetUpdater_cfi::sort, ProxyMTDTopology::specificTopology(), edm::swap(), MtdTruthAccumulator::calo_particles::swap(), GeomDet::topology(), and topology.

                                                                                   {
  using namespace geant_units::operators;

  edm::LogInfo(messageCategory_) << "Adding " << output_.pSimClusters->size() << " SimParticles and "
                                 << output_.pCaloParticles->size() << " CaloParticles to the event.";

  // We need to normalize the hits and energies into hits and fractions (since
  // we have looped over all pileup events)
  // For premixing stage1 we keep the energies, they will be normalized to
  // fractions in stage2

  if (premixStage1_) {
    auto totalEnergies = std::make_unique<std::vector<std::pair<uint64_t, float>>>();
    totalEnergies->reserve(m_detIdToTotalSimEnergy.size());
    std::copy(m_detIdToTotalSimEnergy.begin(), m_detIdToTotalSimEnergy.end(), std::back_inserter(*totalEnergies));
    std::sort(totalEnergies->begin(), totalEnergies->end());
    event.put(std::move(totalEnergies), "MergedMtdTruth");
  } else {
    for (auto &sc : *(output_.pSimClusters)) {
      auto hitsAndEnergies = sc.hits_and_fractions();
      sc.clearHitsAndFractions();
      sc.clearHitsEnergy();
      for (auto &hAndE : hitsAndEnergies) {
        const float totalenergy = m_detIdToTotalSimEnergy[hAndE.first];
        float fraction = 0.;
        if (totalenergy > 0)
          fraction = hAndE.second / totalenergy;
        else
          edm::LogWarning(messageCategory_)
              << "TotalSimEnergy for hit " << hAndE.first << " is 0! The fraction for this hit cannot be computed.";
        sc.addHitAndFraction(hAndE.first, fraction);
        sc.addHitEnergy(hAndE.second);
      }
    }
  }

#ifdef PRINT_DEBUG
  IfLogDebug(DEBUG, messageCategory_) << "SIMCLUSTERS LIST:" << std::endl;
  for (const auto &sc : *(output_.pSimClusters)) {
    IfLogDebug(DEBUG, messageCategory_) << std::fixed << std::setprecision(3) << "SimCluster from CP with:"
                                        << "\n  charge " << sc.charge() << "\n  pdgId  " << sc.pdgId() << "\n  energy "
                                        << sc.energy() << " GeV\n  eta    " << sc.eta() << "\n  phi    " << sc.phi()
                                        << "\n  number of cells = " << sc.hits_and_fractions().size() << std::endl;
    for (unsigned int i = 0; i < sc.hits_and_fractions().size(); ++i) {
      DetId id(sc.detIds_and_rows()[i].first);
      IfLogDebug(DEBUG, messageCategory_)
          << std::fixed << std::setprecision(3) << " hit " << id.rawId() << " on layer " << geomTools_.layer(id)
          << " module " << geomTools_.module(id) << " row " << (unsigned int)(sc.detIds_and_rows()[i].second).first
          << " col " << (unsigned int)(sc.detIds_and_rows()[i].second).second << " at time "
          << sc.hits_and_times()[i].second << " ns" << std::endl;
    }
    IfLogDebug(DEBUG, messageCategory_) << "--------------\n";
  }
  IfLogDebug(DEBUG, messageCategory_) << std::endl;
#endif

  // save the SimCluster orphan handle so we can fill the calo particles
  auto scHandle = event.put(std::move(output_.pSimClusters), "MergedMtdTruth");

  // reserve for the best case scenario: already splitted
  output_.pMtdSimLayerClusters->reserve(scHandle->size());
  output_.pMtdSimTracksters->reserve(scHandle->size());

  uint32_t SC_index = 0;
  uint32_t LC_index = 0;
  for (const auto &sc : *scHandle) {
    auto const &hAndF = sc.hits_and_fractions();
    auto const &hAndE = sc.hits_and_energies();
    auto const &hAndT = sc.hits_and_times();
    auto const &hAndP = sc.hits_and_positions();
    auto const &hAndR = sc.detIds_and_rows();
    // create a vector with the indices of the hits in the simCluster
    std::vector<int> indices(hAndF.size());
    std::iota(indices.begin(), indices.end(), 0);
    // sort the hits indices based on the unique indices created before
    std::sort(indices.begin(), indices.end(), [&](int a, int b) { return hAndF[a].first < hAndF[b].first; });

    // now split the sc: loop on the sorted indices and save the first hit in a
    // temporary simCluster. If the following hit is in the same module and row (column),
    // but next column (row) put it in the temporary simcluster as well, otherwise
    // put the temporary simcluster in the collection and start creating a new one
    std::vector<uint32_t> LC_indices;
    MtdSimLayerCluster tmpLC(sc.g4Tracks()[0]);
    int prev = indices[0];
    DetId prevId(hAndR[prev].first);

    float SimLCenergy = 0.;
    float SimLCx = 0., SimLCy = 0., SimLCz = 0.;

    auto push_back_hit = [&](const int &ind) {
      tmpLC.addHitAndFraction(hAndF[ind].first, hAndF[ind].second);
      tmpLC.addHitEnergy(hAndE[ind].second);
      tmpLC.addHitTime(hAndT[ind].second);
      tmpLC.addHitPosition(hAndP[ind].second);
    };

    auto update_clu_info = [&](const int &ind) {
      double energy = hAndE[ind].second;
      auto position = hAndP[ind].second;
      if (geomTools_.isBTL((DetId)hAndR[ind].first)) {
        BTLDetId detId{(DetId)hAndR[ind].first};
        DetId geoId = detId.geographicalId(MTDTopologyMode::crysLayoutFromTopoMode(topology->getMTDTopologyMode()));
        const MTDGeomDet *thedet = geom->idToDet(geoId);
        const ProxyMTDTopology &topoproxy = static_cast<const ProxyMTDTopology &>(thedet->topology());
        const RectangularMTDTopology &topo = static_cast<const RectangularMTDTopology &>(topoproxy.specificTopology());
        position =
            topo.pixelToModuleLocalPoint(hAndP[ind].second, (hAndR[ind].second).first, (hAndR[ind].second).second);
      }
      SimLCenergy += energy;
      SimLCx += position.x() * energy;
      SimLCy += position.y() * energy;
      SimLCz += position.z() * energy;
    };

    auto push_back_clu = [&](const uint32_t &SC_index, uint32_t &LC_index) {
      tmpLC.addCluEnergy(SimLCenergy);
      LocalPoint SimLCpos(SimLCx / SimLCenergy, SimLCy / SimLCenergy, SimLCz / SimLCenergy);
      tmpLC.addCluLocalPos(SimLCpos);
      SimLCenergy = 0.;
      SimLCx = 0.;
      SimLCy = 0.;
      SimLCz = 0.;
      tmpLC.addCluIndex(SC_index);
      tmpLC.computeClusterTime();
      tmpLC.setTrackIdOffset(sc.trackIdOffset());  // add trackIdoffset
      output_.pMtdSimLayerClusters->push_back(tmpLC);
      LC_indices.push_back(LC_index);
      LC_index++;
      tmpLC.clear();
    };

    // fill tmpLC with the first hit
    push_back_hit(prev);
    update_clu_info(prev);
    for (const auto &ind : indices) {
      if (ind == indices[0])
        continue;
      DetId id(hAndR[ind].first);
      if (geomTools_.isETL(id) != geomTools_.isETL(prevId) or geomTools_.layer(id) != geomTools_.layer(prevId) or
          geomTools_.module(id) != geomTools_.module(prevId) or
          ((hAndR[ind].second).first == (hAndR[prev].second).first and
           (hAndR[ind].second).second != (hAndR[prev].second).second + 1) or
          ((hAndR[ind].second).second == (hAndR[prev].second).second and
           (hAndR[ind].second).first != (hAndR[prev].second).first + 1)) {
        // the next hit is not adjacent to the previous one, put the current temporary cluster in the collection
        // and the hit will be put in an empty temporary cluster
        push_back_clu(SC_index, LC_index);
      }
      // add the hit to the temporary cluster
      push_back_hit(ind);
      update_clu_info(ind);
      prev = ind;
      DetId newId(hAndR[prev].first);
      prevId = newId;
    }
    // add the remaining temporary cluster to the collection
    push_back_clu(SC_index, LC_index);

    // now the simTrackster: find position and time of the first simHit
    // bc right now there is no method to ask the simTrack for pos/time
    // at MTD entrance
    float timeAtEntrance = 99.;
    uint32_t idAtEntrance = 0;
    for (uint32_t i = 0; i < (uint32_t)hAndT.size(); i++) {
      if (hAndT[i].second < timeAtEntrance) {
        timeAtEntrance = hAndT[i].second;
        idAtEntrance = i;
      }
    }

    // sort LCs in the SimTrackster by time
    auto &MtdSimLayerClusters = output_.pMtdSimLayerClusters;
    std::sort(LC_indices.begin(), LC_indices.end(), [&MtdSimLayerClusters](int i, int j) {
      return (*MtdSimLayerClusters)[i].simLCTime() < (*MtdSimLayerClusters)[j].simLCTime();
    });

    GlobalPoint posAtEntrance = geomTools_
                                    .position((DetId)hAndR[idAtEntrance].first,
                                              (hAndR[idAtEntrance].second).first,
                                              (hAndR[idAtEntrance].second).second)
                                    .second;
    output_.pMtdSimTracksters->emplace_back(sc, LC_indices, timeAtEntrance, posAtEntrance);
    SC_index++;
  }

#ifdef PRINT_DEBUG
  IfLogDebug(DEBUG, messageCategory_) << "SIMLAYERCLUSTERS LIST: \n";
  for (auto &sc : *output_.pMtdSimLayerClusters) {
    IfLogDebug(DEBUG, messageCategory_) << std::fixed << std::setprecision(3) << "SimLayerCluster with:"
                                        << "\n  CP charge " << sc.charge() << "\n  CP pdgId  " << sc.pdgId()
                                        << "\n  CP energy " << sc.energy() << " GeV\n  CP eta    " << sc.eta()
                                        << "\n  CP phi    " << sc.phi()
                                        << "\n  number of cells = " << sc.hits_and_fractions().size() << std::endl;
    for (unsigned int i = 0; i < sc.hits_and_fractions().size(); ++i) {
      DetId id(sc.detIds_and_rows()[i].first);
      IfLogDebug(DEBUG, messageCategory_)
          << std::fixed << std::setprecision(3) << " hit " << sc.detIds_and_rows()[i].first << " on layer "
          << geomTools_.layer(id) << " at time " << sc.hits_and_times()[i].second << " ns" << std::endl;
    }
    IfLogDebug(DEBUG, messageCategory_) << std::fixed << std::setprecision(3) << "  Cluster time " << sc.simLCTime()
                                        << " ns \n Cluster pos" << sc.simLCPos() << " cm\n"
                                        << std::fixed << std::setprecision(6) << " Cluster energy "
                                        << convertUnitsTo(0.001_MeV, sc.simLCEnergy()) << " MeV" << std::endl;
    IfLogDebug(DEBUG, messageCategory_) << "--------------\n";
  }
  IfLogDebug(DEBUG, messageCategory_) << std::endl;

  IfLogDebug(DEBUG, messageCategory_) << "SIMTRACKSTERS LIST: \n";
  for (auto &sc : *output_.pMtdSimTracksters) {
    IfLogDebug(DEBUG, messageCategory_) << std::fixed << std::setprecision(3) << "SimTrackster with:"
                                        << "\n  CP charge " << sc.charge() << "\n  CP pdgId  " << sc.pdgId()
                                        << "\n  CP energy " << sc.energy() << " GeV\n  CP eta    " << sc.eta()
                                        << "\n  CP phi    " << sc.phi()
                                        << "\n number of layer clusters = " << sc.numberOfClusters()
                                        << "\n time of first simhit " << sc.time() << " ns\n position of first simhit"
                                        << sc.position() << "cm" << std::endl;
    IfLogDebug(DEBUG, messageCategory_) << "  LCs indices: ";
    for (const auto &lc : sc.clusters())
      IfLogDebug(DEBUG, messageCategory_) << lc << ", ";
    IfLogDebug(DEBUG, messageCategory_) << "\n--------------\n";
  }
  IfLogDebug(DEBUG, messageCategory_) << std::endl;
#endif

  event.put(std::move(output_.pMtdSimLayerClusters), "MergedMtdTruthLC");
  event.put(std::move(output_.pMtdSimTracksters), "MergedMtdTruthST");

  // now fill the calo particles
  for (unsigned i = 0; i < output_.pCaloParticles->size(); ++i) {
    auto &cp = (*output_.pCaloParticles)[i];
    for (unsigned j = m_caloParticles.sc_start_[i]; j < m_caloParticles.sc_stop_[i]; ++j) {
      edm::Ref<MtdSimClusterCollection> ref(scHandle, j);
      cp.addSimCluster(ref);
    }
  }
  event.put(std::move(output_.pCaloParticles), "MergedMtdTruth");

  calo_particles().swap(m_caloParticles);

  std::unordered_map<uint64_t, float>().swap(m_detIdToTotalSimEnergy);
  std::unordered_multimap<Barcode_t, Index_t>().swap(m_simHitBarcodeToIndex);
}

initializeEvent()

void MtdTruthAccumulator::initializeEvent ( const edm::Event &  event, const edm::EventSetup &  setup  )

overrideprivatevirtual

Implements DigiAccumulatorMixMod.

Definition at line 300 of file MtdTruthAccumulator.cc.

References geom, geomToken_, geomTools_, m_detIdToTotalSimEnergy, mtdtopoToken_, output_, MtdTruthAccumulator::OutputCollections::pCaloParticles, MtdTruthAccumulator::OutputCollections::pMtdSimLayerClusters, MtdTruthAccumulator::OutputCollections::pMtdSimTracksters, MtdTruthAccumulator::OutputCollections::pSimClusters, mtd::MTDGeomUtil::setGeometry(), mtd::MTDGeomUtil::setTopology(), singleTopDQM_cfi::setup, and topology.

                                                                                           {
  output_.pSimClusters = std::make_unique<MtdSimClusterCollection>();
  output_.pCaloParticles = std::make_unique<MtdCaloParticleCollection>();

  output_.pMtdSimLayerClusters = std::make_unique<MtdSimLayerClusterCollection>();
  output_.pMtdSimTracksters = std::make_unique<MtdSimTracksterCollection>();

  m_detIdToTotalSimEnergy.clear();

  auto geometryHandle = setup.getTransientHandle(geomToken_);
  geom = geometryHandle.product();

  auto topologyHandle = setup.getTransientHandle(mtdtopoToken_);
  topology = topologyHandle.product();

  geomTools_.setGeometry(geom);
  geomTools_.setTopology(topology);
}

Member Data Documentation

bunchSpacing_

const unsigned int MtdTruthAccumulator::bunchSpacing_

private

Definition at line 109 of file MtdTruthAccumulator.cc.

Referenced by accumulateEvent().

collectionTags_

std::vector<edm::InputTag> MtdTruthAccumulator::collectionTags_

private

Definition at line 116 of file MtdTruthAccumulator.cc.

Referenced by fillSimHits(), and MtdTruthAccumulator().

genParticleLabel_

edm::InputTag MtdTruthAccumulator::genParticleLabel_

private

Definition at line 117 of file MtdTruthAccumulator.cc.

Referenced by accumulateEvent(), and MtdTruthAccumulator().

geom

const MTDGeometry* MtdTruthAccumulator::geom = nullptr

private

Definition at line 155 of file MtdTruthAccumulator.cc.

Referenced by finalizeEvent(), and initializeEvent().

geomToken_

const edm::ESGetToken<MTDGeometry, MTDDigiGeometryRecord> MtdTruthAccumulator::geomToken_

private

Definition at line 120 of file MtdTruthAccumulator.cc.

Referenced by initializeEvent().

geomTools_

mtd::MTDGeomUtil MtdTruthAccumulator::geomTools_

private

Definition at line 124 of file MtdTruthAccumulator.cc.

Referenced by fillSimHits(), finalizeEvent(), and initializeEvent().

hepMCproductLabel_

edm::InputTag MtdTruthAccumulator::hepMCproductLabel_

private

Needed to add HepMC::GenVertex to SimVertex.

Definition at line 119 of file MtdTruthAccumulator.cc.

Referenced by accumulate(), and MtdTruthAccumulator().

hSimTracks

edm::Handle<std::vector<SimTrack> > MtdTruthAccumulator::hSimTracks

private

Definition at line 113 of file MtdTruthAccumulator.cc.

Referenced by accumulateEvent().

hSimVertices

edm::Handle<std::vector<SimVertex> > MtdTruthAccumulator::hSimVertices

private

Definition at line 114 of file MtdTruthAccumulator.cc.

Referenced by accumulateEvent().

isEtl_

bool MtdTruthAccumulator::isEtl_

private

Definition at line 129 of file MtdTruthAccumulator.cc.

Referenced by MtdTruthAccumulator().

m_caloParticles

calo_particles MtdTruthAccumulator::m_caloParticles

private

Definition at line 158 of file MtdTruthAccumulator.cc.

Referenced by accumulateEvent(), and finalizeEvent().

m_detIdToTotalSimEnergy

std::unordered_map<uint64_t, float> MtdTruthAccumulator::m_detIdToTotalSimEnergy

private

Definition at line 93 of file MtdTruthAccumulator.cc.

Referenced by fillSimHits(), finalizeEvent(), and initializeEvent().

m_simHitBarcodeToIndex

std::unordered_multimap<Barcode_t, Index_t> MtdTruthAccumulator::m_simHitBarcodeToIndex

private

Definition at line 94 of file MtdTruthAccumulator.cc.

Referenced by accumulateEvent(), and finalizeEvent().

maximumPreviousBunchCrossing_

const unsigned int MtdTruthAccumulator::maximumPreviousBunchCrossing_

private

The maximum bunch crossing BEFORE the signal crossing to create TrackinParticles for. Use positive values. If set to zero no previous bunches are added and only in-time, signal and after bunches (defined by maximumSubsequentBunchCrossing_) are used.

Definition at line 101 of file MtdTruthAccumulator.cc.

Referenced by accumulate().

maximumSubsequentBunchCrossing_

const unsigned int MtdTruthAccumulator::maximumSubsequentBunchCrossing_

private

The maximum bunch crossing AFTER the signal crossing to create TrackinParticles for. E.g. if set to zero only uses the signal and in time pileup (and previous bunches defined by the maximumPreviousBunchCrossing_ parameter).

Definition at line 107 of file MtdTruthAccumulator.cc.

Referenced by accumulate().

maxPseudoRapidity_

const double MtdTruthAccumulator::maxPseudoRapidity_

private

Definition at line 126 of file MtdTruthAccumulator.cc.

Referenced by accumulateEvent().

messageCategory_

const std::string MtdTruthAccumulator::messageCategory_

private

Definition at line 91 of file MtdTruthAccumulator.cc.

Referenced by accumulate(), accumulateEvent(), fillSimHits(), and finalizeEvent().

minEnergy_

const double MtdTruthAccumulator::minEnergy_

private

Definition at line 126 of file MtdTruthAccumulator.cc.

Referenced by accumulateEvent().

mtdtopoToken_

const edm::ESGetToken<MTDTopology, MTDTopologyRcd> MtdTruthAccumulator::mtdtopoToken_

private

Definition at line 121 of file MtdTruthAccumulator.cc.

Referenced by initializeEvent().

output_

OutputCollections MtdTruthAccumulator::output_

private

Definition at line 157 of file MtdTruthAccumulator.cc.

Referenced by accumulateEvent(), finalizeEvent(), and initializeEvent().

premixStage1_

const bool MtdTruthAccumulator::premixStage1_

private

Definition at line 127 of file MtdTruthAccumulator.cc.

Referenced by finalizeEvent(), and MtdTruthAccumulator().

simTrackLabel_

const edm::InputTag MtdTruthAccumulator::simTrackLabel_

private

Definition at line 111 of file MtdTruthAccumulator.cc.

Referenced by accumulateEvent(), and MtdTruthAccumulator().

simVertexLabel_

const edm::InputTag MtdTruthAccumulator::simVertexLabel_

private

Definition at line 112 of file MtdTruthAccumulator.cc.

Referenced by accumulateEvent(), and MtdTruthAccumulator().

topology

const MTDTopology* MtdTruthAccumulator::topology = nullptr

private

Definition at line 156 of file MtdTruthAccumulator.cc.

Referenced by finalizeEvent(), and initializeEvent().