first release of TICLv5

CommonTools/RecoAlgos/interface/MultiVectorManager.h

@@ -0,0 +1,79 @@
+// Author: Felice Pantaleo (CERN), 2023, felice.pantaleo@cern.ch
+#ifndef MultiVectorManager_h
+#define MultiVectorManager_h
+
+#include <vector>
+#include <cassert>
+#include <algorithm>
+#include <iterator>
+
+template <typename T>
+class MultiVectorManager {
+public:
+  void addVector(const std::vector<T>& vec) {
+    vectors.emplace_back(vec.begin(), vec.end());
+    offsets.push_back(totalSize);
+    totalSize += vec.size();
+  }
+
+  T& operator[](size_t globalIndex) {
+    return const_cast<T&>(static_cast<const MultiVectorManager*>(this)->operator[](globalIndex));
+  }
+
+  const T& operator[](size_t globalIndex) const {
+    assert(globalIndex < totalSize && "Global index out of range");
+
+    auto it = std::upper_bound(offsets.begin(), offsets.end(), globalIndex);
+    size_t vectorIndex = std::distance(offsets.begin(), it) - 1;
+    size_t localIndex = globalIndex - offsets[vectorIndex];
+
+    return vectors[vectorIndex][localIndex];
+  }
+
+  size_t getGlobalIndex(size_t vectorIndex, size_t localIndex) const {
+    assert(vectorIndex < vectors.size() && "Vector index out of range");
+
+    const auto& vec = vectors[vectorIndex];
+    assert(localIndex < vec.size() && "Local index out of range");
+
+    return offsets[vectorIndex] + localIndex;
+  }
+
+  size_t size() const { return totalSize; }
+
+  std::pair<size_t, size_t> getVectorAndLocalIndex(size_t globalIndex) const {
+    assert(globalIndex < totalSize && "Global index out of range");
+
+    auto it = std::upper_bound(offsets.begin(), offsets.end(), globalIndex);
+    size_t vectorIndex = std::distance(offsets.begin(), it) - 1;
+    size_t localIndex = globalIndex - offsets[vectorIndex];
+
+    return {vectorIndex, localIndex};
+  }
+
+  class Iterator {
+  public:
+    Iterator(const MultiVectorManager& manager, size_t index) : manager(manager), currentIndex(index) {}
+
+    bool operator!=(const Iterator& other) const { return currentIndex != other.currentIndex; }
+
+    T& operator*() const { return const_cast<T&>(manager[currentIndex]); }
+
+    void operator++() { ++currentIndex; }
+
+  private:
+    const MultiVectorManager& manager;
+    size_t currentIndex;
+  };
+
+  Iterator begin() const { return Iterator(*this, 0); }
+
+  Iterator end() const { return Iterator(*this, totalSize); }
+
+private:
+  std::vector<std::vector<T>> vectors;
+  std::vector<size_t> offsets;
+  size_t totalSize = 0;
+};
+
+#endif

Configuration/EventContent/python/EventContent_cff.py

@@ -568,6 +568,8 @@ def SwapKeepAndDrop(l):
     outputCommands = FEVTEventContent.outputCommands + RecoLocalFastTimeFEVT.outputCommands)
 phase2_timing_layer.toModify(FEVTEventContent, 
     outputCommands = FEVTEventContent.outputCommands + RecoMTDFEVT.outputCommands)
+from Configuration.ProcessModifiers.ticl_v5_cff import ticl_v5
+ticl_v5.toModify(FEVTEventContent, outputCommands=FEVTEventContent.outputCommands+TICLv5_FEVT.outputCommands)
 
 FEVTHLTALLEventContent = cms.PSet(
     outputCommands = cms.untracked.vstring('drop *'),
@@ -651,6 +653,8 @@ def SwapKeepAndDrop(l):
                                   'keep *_hltSiStripClusters2ApproxClusters_*_*',
                                   'keep DetIds_hltSiStripRawToDigi_*_*'
                               ])
+
+ticl_v5.toModify(FEVTDEBUGEventContent, outputCommands=FEVTDEBUGEventContent.outputCommands+TICLv5_FEVT.outputCommands)
 #
 #
 # FEVTDEBUGHLT Data Tier definition

Configuration/ProcessModifiers/python/ticl_v5_cff.py

@@ -0,0 +1,5 @@
+import FWCore.ParameterSet.Config as cms
+
+# This modifier is for running TICL v5.
+
+ticl_v5 =  cms.Modifier()

Configuration/PyReleaseValidation/python/upgradeWorkflowComponents.py

@@ -689,15 +689,15 @@ def condition(self, fragment, stepList, key, hasHarvest):
 upgradeWFs['ticl_FastJet'].step3 = {'--procModifiers': 'fastJetTICL'}
 upgradeWFs['ticl_FastJet'].step4 = {'--procModifiers': 'fastJetTICL'}
 
-class UpgradeWorkflow_ticl_v3(UpgradeWorkflow):
+class UpgradeWorkflow_ticl_v5(UpgradeWorkflow):
     def setup_(self, step, stepName, stepDict, k, properties):
         if 'RecoGlobal' in step:
             stepDict[stepName][k] = merge([self.step3, stepDict[step][k]])
         if 'HARVESTGlobal' in step:
             stepDict[stepName][k] = merge([self.step4, stepDict[step][k]])
     def condition(self, fragment, stepList, key, hasHarvest):
         return (fragment=="TTbar_14TeV" or 'CloseByP' in fragment or 'Eta1p7_2p7' in fragment) and '2026' in key
-upgradeWFs['ticl_v3'] = UpgradeWorkflow_ticl_v3(
+upgradeWFs['ticl_v5'] = UpgradeWorkflow_ticl_v5(
     steps = [
         'RecoGlobal',
         'HARVESTGlobal'
@@ -706,12 +706,11 @@ def condition(self, fragment, stepList, key, hasHarvest):
         'RecoGlobal',
         'HARVESTGlobal'
     ],
-    suffix = '_ticl_v3',
+    suffix = '_ticl_v5',
     offset = 0.203,
 )
-upgradeWFs['ticl_v3'].step3 = {'--procModifiers': 'ticl_v3'}
-upgradeWFs['ticl_v3'].step4 = {'--procModifiers': 'ticl_v3'}
-
+upgradeWFs['ticl_v5'].step3 = {'--procModifiers': 'ticl_v5'}
+upgradeWFs['ticl_v5'].step4 = {'--procModifiers': 'ticl_v5'}
 
 # Track DNN workflows
 class UpgradeWorkflow_trackdnn(UpgradeWorkflow):

DataFormats/CaloRecHit/interface/CaloCluster.h

@@ -39,6 +39,7 @@ namespace reco {
       particleFlow = 5,
       hgcal_em = 6,
       hgcal_had = 7,
+      hgcal_scintillator = 8,
       hfnose = 9,
       undefined = 1000
     };

DataFormats/GeometryVector/src/classes.h

@@ -1,3 +1,5 @@
+#include "DataFormats/Common/interface/ValueMap.h"
+#include "DataFormats/Common/interface/Wrapper.h"
 #include "DataFormats/GeometryVector/interface/Phi.h"
 namespace {
   Geom::Phi<double, Geom::MinusPiToPi> dummy;

DataFormats/GeometryVector/src/classes_def.xml

@@ -88,6 +88,8 @@
    <version ClassVersion="11" checksum="122348193"/>
    <version ClassVersion="10" checksum="3782271882"/>
   </class>
+  <class name="edm::ValueMap<Point3DBase<float,GlobalTag> >"/>
+  <class name="edm::Wrapper<edm::ValueMap<Point3DBase<float,GlobalTag> > >"/>
   <class name="Vector2DBase<float,GlobalTag>" ClassVersion="11">
    <version ClassVersion="11" checksum="2993275990"/>
    <version ClassVersion="10" checksum="662089398"/>

DataFormats/HGCalReco/BuildFile.xml

@@ -1,5 +1,7 @@
 <use name="DataFormats/Math"/>
 <use name="DataFormats/Provenance"/>
+<use name="DataFormats/Portable"/>
+<use name="DataFormats/SoATemplate"/>
 <use name="DataFormats/Candidate"/>
 <use name="DataFormats/Common"/>
 <use name="DataFormats/TrackReco"/>

DataFormats/HGCalReco/interface/MtdHostCollection.h

@@ -0,0 +1,12 @@
+#ifndef DataFormats_HGCalReco_MtdHostCollection_h
+#define DataFormats_HGCalReco_MtdHostCollection_h
+
+#include "DataFormats/Portable/interface/PortableHostCollection.h"
+#include "DataFormats/HGCalReco/interface/MtdSoA.h"
+
+// MtdSoA in host memory
+using MtdHostCollection = PortableHostCollection<MtdSoA>;
+using MtdHostCollectionView = PortableHostCollection<MtdSoA>::View;
+using MtdHostCollectionConstView = PortableHostCollection<MtdSoA>::ConstView;
+
+#endif

DataFormats/HGCalReco/interface/MtdSoA.h

@@ -0,0 +1,27 @@
+#ifndef DataFormats_HGCalReco_MtdSoA_h
+#define DataFormats_HGCalReco_MtdSoA_h
+
+#include "DataFormats/SoATemplate/interface/SoALayout.h"
+
+GENERATE_SOA_LAYOUT(MtdSoALayout,
+                    SOA_COLUMN(int32_t, trackAsocMTD),
+                    SOA_COLUMN(float, time0),
+                    SOA_COLUMN(float, time0Err),
+                    SOA_COLUMN(float, time),
+                    SOA_COLUMN(float, timeErr),
+                    SOA_COLUMN(float, MVAquality),
+                    SOA_COLUMN(float, pathLength),
+                    SOA_COLUMN(float, beta),
+                    SOA_COLUMN(float, posInMTD_x),
+                    SOA_COLUMN(float, posInMTD_y),
+                    SOA_COLUMN(float, posInMTD_z),
+                    SOA_COLUMN(float, momentumWithMTD),
+                    SOA_COLUMN(float, probPi),
+                    SOA_COLUMN(float, probK),
+                    SOA_COLUMN(float, probP))
+
+using MtdSoA = MtdSoALayout<>;
+using MtdSoAView = MtdSoA::View;
+using MtdSoAConstView = MtdSoA::ConstView;
+
+#endif

DataFormats/HGCalReco/interface/TICLCandidate.h

@@ -7,6 +7,9 @@
 #include "DataFormats/Math/interface/Point3D.h"
 #include "DataFormats/TrackReco/interface/Track.h"
 #include "DataFormats/TrackReco/interface/TrackFwd.h"
+#include "RecoHGCal/TICL/interface/commons.h"
+
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
 
 // A TICLCandidate is a lightweight physics object made from one or multiple Tracksters.
 
@@ -15,23 +18,82 @@ class TICLCandidate : public reco::LeafCandidate {
   typedef ticl::Trackster::ParticleType ParticleType;
 
   TICLCandidate(Charge q, const LorentzVector& p4)
-      : LeafCandidate(q, p4), time_(0.f), timeError_(-1.f), rawEnergy_(0.f), idProbabilities_{} {}
+      : LeafCandidate(q, p4), idProbabilities_{}, time_(0.f), timeError_(-1.f), rawEnergy_(0.f) {}
 
-  TICLCandidate() : LeafCandidate(), time_(0.f), timeError_(-1.f), rawEnergy_(0.f), idProbabilities_{} {}
+  TICLCandidate() : LeafCandidate(), idProbabilities_{}, time_(0.f), timeError_(-1.f), rawEnergy_(0.f) {}
 
   TICLCandidate(const edm::Ptr<ticl::Trackster>& trackster)
       : LeafCandidate(),
+        tracksters_({trackster}),
+        idProbabilities_{},
         time_(trackster->time()),
         timeError_(trackster->timeError()),
-        rawEnergy_(0.f),
-        tracksters_({trackster}),
-        idProbabilities_{} {}
+        MTDtime_{0.f},
+        MTDtimeError_{-1.f},
+        rawEnergy_(0.f) {}
+
+  TICLCandidate(const edm::Ptr<reco::Track> trackPtr, const edm::Ptr<ticl::Trackster>& tracksterPtr)
+      : LeafCandidate(), tracksters_{}, trackPtr_(trackPtr), time_(0.f), timeError_(-1.f) {
+    if (trackPtr_.isNull() and tracksterPtr.isNull())
+      throw cms::Exception("NullPointerError")
+          << "TICLCandidate constructor: at least one between track and trackster must be valid";
+
+    if (tracksterPtr.isNonnull()) {
+      tracksters_.push_back(tracksterPtr);
+      auto const& trackster = tracksters_[0].get();
+      idProbabilities_ = trackster->id_probabilities();
+      if (trackPtr_.isNonnull()) {
+        auto pdgId = trackster->isHadronic() ? 211 : 11;
+        auto const& tk = trackPtr_.get();
+        setPdgId(pdgId * tk->charge());
+        setCharge(tk->charge());
+        rawEnergy_ = trackster->raw_energy();
+        auto const& regrE = trackster->regressed_energy();
+        math::XYZTLorentzVector p4(regrE * tk->momentum().unit().x(),
+                                   regrE * tk->momentum().unit().y(),
+                                   regrE * tk->momentum().unit().z(),
+                                   regrE);
+        setP4(p4);
+
+      } else {
+        auto pdgId = trackster->isHadronic() ? 130 : 22;
+        setPdgId(pdgId);
+        setCharge(0);
+        rawEnergy_ = trackster->raw_energy();
+        const float& regrE = trackster->regressed_energy();
+        math::XYZTLorentzVector p4(regrE * trackster->barycenter().unit().x(),
+                                   regrE * trackster->barycenter().unit().y(),
+                                   regrE * trackster->barycenter().unit().z(),
+                                   regrE);
+        setP4(p4);
+      }
+    } else {
+      //candidate from track only
+      auto const& tk = trackPtr_.get();
+      setPdgId(211 * tk->charge());
+      setCharge(tk->charge());
+      const float energy = std::sqrt(tk->p() * tk->p() + ticl::mpion2);
+      setRawEnergy(energy);
+      math::PtEtaPhiMLorentzVector p4Polar(tk->pt(), tk->eta(), tk->phi(), ticl::mpion);
+      setP4(p4Polar);
+    }
+  }
 
   inline float time() const { return time_; }
   inline float timeError() const { return timeError_; }
 
-  void setTime(float time) { time_ = time; };
-  void setTimeError(float timeError) { timeError_ = timeError; }
+  void setTime(float time, float timeError) {
+    time_ = time;
+    timeError_ = timeError;
+  };
+
+  inline float MTDtime() const { return MTDtime_; }
+  inline float MTDtimeError() const { return MTDtimeError_; }
+
+  void setMTDTime(float time, float timeError) {
+    MTDtime_ = time;
+    MTDtimeError_ = timeError;
+  };
 
   inline const edm::Ptr<reco::Track> trackPtr() const { return trackPtr_; }
   void setTrackPtr(const edm::Ptr<reco::Track>& trackPtr) { trackPtr_ = trackPtr; }
@@ -65,17 +127,17 @@ class TICLCandidate : public reco::LeafCandidate {
   inline void setIdProbability(ParticleType type, float value) { idProbabilities_[int(type)] = value; }
 
 private:
-  float time_;
-  float timeError_;
-  edm::Ptr<reco::Track> trackPtr_;
-
-  float rawEnergy_;
-
   // vector of Ptr so Tracksters can come from different collections
   // and there can be derived classes
   std::vector<edm::Ptr<ticl::Trackster> > tracksters_;
-
+  edm::Ptr<reco::Track> trackPtr_;
   // Since it contains multiple tracksters, duplicate the probability interface
   std::array<float, 8> idProbabilities_;
+
+  float time_;
+  float timeError_;
+  float MTDtime_;
+  float MTDtimeError_;
+  float rawEnergy_;
 };
 #endif

DataFormats/HGCalReco/interface/TICLLayerTile.h

@@ -12,9 +12,7 @@ class TICLLayerTileT {
 public:
   typedef T type;
 
-  void fill(double eta, double phi, unsigned int layerClusterId) {
-    tile_[globalBin(eta, phi)].push_back(layerClusterId);
-  }
+  void fill(float eta, float phi, unsigned int layerClusterId) { tile_[globalBin(eta, phi)].push_back(layerClusterId); }
 
   int etaBin(float eta) const {
     constexpr float etaRange = T::maxEta - T::minEta;
@@ -62,7 +60,7 @@ class TICLLayerTileT {
 
   int globalBin(int etaBin, int phiBin) const { return phiBin + etaBin * T::nPhiBins; }
 
-  int globalBin(double eta, double phi) const { return phiBin(phi) + etaBin(eta) * T::nPhiBins; }
+  int globalBin(float eta, float phi) const { return phiBin(phi) + etaBin(eta) * T::nPhiBins; }
 
   void clear() {
     auto nBins = T::nEtaBins * T::nPhiBins;
@@ -96,7 +94,7 @@ class TICLGenericTile {
   // numbering is not handled internally. It is the user's responsibility to
   // properly use and consistently access it here.
   const auto& operator[](int index) const { return tiles_[index]; }
-  void fill(int index, double eta, double phi, unsigned int objectId) { tiles_[index].fill(eta, phi, objectId); }
+  void fill(int index, float eta, float phi, unsigned int objectId) { tiles_[index].fill(eta, phi, objectId); }
 
 private:
   T tiles_;