TICL/plugins/GeneralInterpretationAlgo.cc

0001 #include "RecoHGCal/TICL/interface/TICLInterpretationAlgoBase.h"
0002 #include "RecoHGCal/TICL/plugins/GeneralInterpretationAlgo.h"
0003 #include "RecoParticleFlow/PFProducer/interface/PFMuonAlgo.h"
0004 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0005
0006 using namespace ticl;
0007
0008 using Vector = ticl::Trackster::Vector;
0009
0010 GeneralInterpretationAlgo::~GeneralInterpretationAlgo() {}
0011
0012 GeneralInterpretationAlgo::GeneralInterpretationAlgo(const edm::ParameterSet &conf, edm::ConsumesCollector cc)
0013     : TICLInterpretationAlgoBase(conf, cc),
0014       del_tk_ts_layer1_(conf.getParameter<double>("delta_tk_ts_layer1")),
0015       del_tk_ts_int_(conf.getParameter<double>("delta_tk_ts_interface")),
0016       timing_quality_threshold_(conf.getParameter<double>("timing_quality_threshold")) {}
0017
0018 void GeneralInterpretationAlgo::initialize(const HGCalDDDConstants *hgcons,
0019                                            const hgcal::RecHitTools rhtools,
0020                                            const edm::ESHandle<MagneticField> bfieldH,
0021                                            const edm::ESHandle<Propagator> propH) {
0022   hgcons_ = hgcons;
0023   rhtools_ = rhtools;
0024   buildLayers();
0025
0026   bfield_ = bfieldH;
0027   propagator_ = propH;
0028 }
0029
0030 void GeneralInterpretationAlgo::buildLayers() {
0031   // build disks at HGCal front & EM-Had interface for track propagation
0032
0033   float zVal = hgcons_->waferZ(1, true);
0034   std::pair<float, float> rMinMax = hgcons_->rangeR(zVal, true);
0035
0036   float zVal_interface = rhtools_.getPositionLayer(rhtools_.lastLayerEE()).z();
0037   std::pair<float, float> rMinMax_interface = hgcons_->rangeR(zVal_interface, true);
0038
0039   for (int iSide = 0; iSide < 2; ++iSide) {
0040     float zSide = (iSide == 0) ? (-1. * zVal) : zVal;
0041     firstDisk_[iSide] =
0042         std::make_unique<GeomDet>(Disk::build(Disk::PositionType(0, 0, zSide),
0043                                               Disk::RotationType(),
0044                                               SimpleDiskBounds(rMinMax.first, rMinMax.second, zSide - 0.5, zSide + 0.5))
0045                                       .get());
0046
0047     zSide = (iSide == 0) ? (-1. * zVal_interface) : zVal_interface;
0048     interfaceDisk_[iSide] = std::make_unique<GeomDet>(
0049         Disk::build(Disk::PositionType(0, 0, zSide),
0050                     Disk::RotationType(),
0051                     SimpleDiskBounds(rMinMax_interface.first, rMinMax_interface.second, zSide - 0.5, zSide + 0.5))
0052             .get());
0053   }
0054 }
0055 Vector GeneralInterpretationAlgo::propagateTrackster(const Trackster &t,
0056                                                      const unsigned idx,
0057                                                      float zVal,
0058                                                      std::array<TICLLayerTile, 2> &tracksterTiles) {
0059   // needs only the positive Z co-ordinate of the surface to propagate to
0060   // the correct sign is calculated inside according to the barycenter of trackster
0061   Vector const &baryc = t.barycenter();
0062   Vector directnv = t.eigenvectors(0);
0063
0064   // barycenter as direction for tracksters w/ poor PCA
0065   // propagation still done to get the cartesian coords
0066   // which are anyway converted to eta, phi in linking
0067   // -> can be simplified later
0068
0069   //FP: disable PCA propagation for the moment and fallback to barycenter position
0070   // if (t.eigenvalues()[0] / t.eigenvalues()[1] < 20)
0071   directnv = baryc.unit();
0072   zVal *= (baryc.Z() > 0) ? 1 : -1;
0073   float par = (zVal - baryc.Z()) / directnv.Z();
0074   float xOnSurface = par * directnv.X() + baryc.X();
0075   float yOnSurface = par * directnv.Y() + baryc.Y();
0076   Vector tPoint(xOnSurface, yOnSurface, zVal);
0077   if (tPoint.Eta() > 0) {
0078     tracksterTiles[1].fill(tPoint.Eta(), tPoint.Phi(), idx);
0079   } else if (tPoint.Eta() < 0) {
0080     tracksterTiles[0].fill(tPoint.Eta(), tPoint.Phi(), idx);
0081   }
0082
0083   return tPoint;
0084 }
0085
0086 void GeneralInterpretationAlgo::findTrackstersInWindow(const MultiVectorManager<Trackster> &tracksters,
0087                                                        const std::vector<std::pair<Vector, unsigned>> &seedingCollection,
0088                                                        const std::array<TICLLayerTile, 2> &tracksterTiles,
0089                                                        const std::vector<Vector> &tracksterPropPoints,
0090                                                        const float delta,
0091                                                        unsigned trackstersSize,
0092                                                        std::vector<std::vector<unsigned>> &resultCollection,
0093                                                        bool useMask = false) {
0094   // Finds tracksters in tracksterTiles within an eta-phi window
0095   // (given by delta) of the objects (track/trackster) in the seedingCollection.
0096   // Element i in resultCollection is the vector of trackster
0097   // indices found close to the i-th object in the seedingCollection.
0098   // If specified, Tracksters are masked once found as close to an object.
0099   std::vector<int> mask(trackstersSize, 0);
0100   const float delta2 = delta * delta;
0101
0102   for (auto &i : seedingCollection) {
0103     float seed_eta = i.first.Eta();
0104     float seed_phi = i.first.Phi();
0105     unsigned seedId = i.second;
0106     auto sideZ = seed_eta > 0;  //forward or backward region
0107     const TICLLayerTile &tile = tracksterTiles[sideZ];
0108     float eta_min = std::max(std::fabs(seed_eta) - delta, (float)TileConstants::minEta);
0109     float eta_max = std::min(std::fabs(seed_eta) + delta, (float)TileConstants::maxEta);
0110
0111     // get range of bins touched by delta
0112     std::array<int, 4> search_box = tile.searchBoxEtaPhi(eta_min, eta_max, seed_phi - delta, seed_phi + delta);
0113
0114     std::vector<unsigned> in_delta;
0115     // std::vector<float> distances2;
0116     std::vector<float> energies;
0117     for (int eta_i = search_box[0]; eta_i <= search_box[1]; ++eta_i) {
0118       for (int phi_i = search_box[2]; phi_i <= search_box[3]; ++phi_i) {
0119         const auto &in_tile = tile[tile.globalBin(eta_i, (phi_i % TileConstants::nPhiBins))];
0120         for (const unsigned &t_i : in_tile) {
0121           // calculate actual distances of tracksters to the seed for a more accurate cut
0122           auto sep2 = (tracksterPropPoints[t_i].Eta() - seed_eta) * (tracksterPropPoints[t_i].Eta() - seed_eta) +
0123                       (tracksterPropPoints[t_i].Phi() - seed_phi) * (tracksterPropPoints[t_i].Phi() - seed_phi);
0124           if (sep2 < delta2) {
0125             in_delta.push_back(t_i);
0126             // distances2.push_back(sep2);
0127             energies.push_back(tracksters[t_i].raw_energy());
0128           }
0129         }
0130       }
0131     }
0132
0133     // sort tracksters found in ascending order of their distances from the seed
0134     std::vector<unsigned> indices(in_delta.size());
0135     std::iota(indices.begin(), indices.end(), 0);
0136     std::sort(indices.begin(), indices.end(), [&](unsigned i, unsigned j) { return energies[i] > energies[j]; });
0137
0138     // push back sorted tracksters in the result collection
0139     for (const unsigned &index : indices) {
0140       const auto &t_i = in_delta[index];
0141       if (!mask[t_i]) {
0142         resultCollection[seedId].push_back(t_i);
0143         if (useMask)
0144           mask[t_i] = 1;
0145       }
0146     }
0147
0148   }  // seeding collection loop
0149 }
0150
0151 bool GeneralInterpretationAlgo::timeAndEnergyCompatible(float &total_raw_energy,
0152                                                         const reco::Track &track,
0153                                                         const Trackster &trackster,
0154                                                         const float &tkT,
0155                                                         const float &tkTErr,
0156                                                         const float &tkQual,
0157                                                         const float &tkBeta,
0158                                                         const GlobalPoint &tkMtdPos,
0159                                                         bool useMTDTiming) {
0160   float threshold = std::min(0.2 * trackster.raw_energy(), 10.0);
0161   bool energyCompatible = (total_raw_energy + trackster.raw_energy() < track.p() + threshold);
0162
0163   if (!useMTDTiming)
0164     return energyCompatible;
0165
0166   // compatible if trackster time is within 3sigma of
0167   // track time; compatible if either: no time assigned
0168   // to trackster or track
0169
0170   float tsT = trackster.time();
0171   float tsTErr = trackster.timeError();
0172
0173   bool timeCompatible = false;
0174   if (tsT == -99. or tkTErr == -1 or tkQual < timing_quality_threshold_)
0175     timeCompatible = true;
0176   else {
0177     const auto &barycenter = trackster.barycenter();
0178
0179     const auto deltaSoverV = std::sqrt((barycenter.x() - tkMtdPos.x()) * (barycenter.x() - tkMtdPos.x()) +
0180                                        (barycenter.y() - tkMtdPos.y()) * (barycenter.y() - tkMtdPos.y()) +
0181                                        (barycenter.z() - tkMtdPos.z()) * (barycenter.z() - tkMtdPos.z())) /
0182                              (tkBeta * 29.9792458);
0183
0184     const auto deltaT = tsT - tkT;
0185
0186     //  timeCompatible = (std::abs(deltaSoverV - deltaT) < maxDeltaT_ * sqrt(tsTErr * tsTErr + tkTErr * tkTErr));
0187     // use sqrt(2) * error on the track for the total error, because the time of the trackster is too small
0188     timeCompatible = std::abs(deltaSoverV - deltaT) < maxDeltaT_ * std::sqrt(tsTErr * tsTErr + tkTErr * tkTErr);
0189   }
0190
0191   if (TICLInterpretationAlgoBase::algo_verbosity_ > VerbosityLevel::Advanced) {
0192     if (!(energyCompatible))
0193       LogDebug("GeneralInterpretationAlgo")
0194           << "energy incompatible : track p " << track.p() << " trackster energy " << trackster.raw_energy() << "\n"
0195           << "                      total_raw_energy " << total_raw_energy << " greater than track p + threshold "
0196           << track.p() + threshold << "\n";
0197     if (!(timeCompatible))
0198       LogDebug("GeneralInterpretationAlgo") << "time incompatible : track time " << tkT << " +/- " << tkTErr
0199                                             << " trackster time " << tsT << " +/- " << tsTErr << "\n";
0200   }
0201
0202   return energyCompatible && timeCompatible;
0203 }
0204
0205 void GeneralInterpretationAlgo::makeCandidates(const Inputs &input,
0206                                                edm::Handle<MtdHostCollection> inputTiming_h,
0207                                                std::vector<Trackster> &resultTracksters,
0208                                                std::vector<int> &resultCandidate) {
0209   bool useMTDTiming = inputTiming_h.isValid();
0210   const auto tkH = input.tracksHandle;
0211   const auto maskTracks = input.maskedTracks;
0212   const auto &tracks = *tkH;
0213   const auto &tracksters = input.tracksters;
0214
0215   auto bFieldProd = bfield_.product();
0216   const Propagator &prop = (*propagator_);
0217
0218   // propagated point collections
0219   // elements in the propagated points collecions are used
0220   // to look for potential linkages in the appropriate tiles
0221   std::vector<std::pair<Vector, unsigned>> trackPColl;     // propagated track points and index of track in collection
0222   std::vector<std::pair<Vector, unsigned>> tkPropIntColl;  // tracks propagated to lastLayerEE
0223
0224   trackPColl.reserve(tracks.size());
0225   tkPropIntColl.reserve(tracks.size());
0226
0227   std::array<TICLLayerTile, 2> tracksterPropTiles = {};  // all Tracksters, propagated to layer 1
0228   std::array<TICLLayerTile, 2> tsPropIntTiles = {};      // all Tracksters, propagated to lastLayerEE
0229
0230   if (TICLInterpretationAlgoBase::algo_verbosity_ > VerbosityLevel::Advanced)
0231     LogDebug("GeneralInterpretationAlgo") << "------- Geometric Linking ------- \n";
0232
0233   // Propagate tracks
0234   std::vector<unsigned> candidateTrackIds;
0235   candidateTrackIds.reserve(tracks.size());
0236   for (unsigned i = 0; i < tracks.size(); ++i) {
0237     if (!maskTracks[i])
0238       continue;
0239     candidateTrackIds.push_back(i);
0240   }
0241
0242   std::sort(candidateTrackIds.begin(), candidateTrackIds.end(), [&](unsigned i, unsigned j) {
0243     return tracks[i].p() > tracks[j].p();
0244   });
0245
0246   for (auto const i : candidateTrackIds) {
0247     const auto &tk = tracks[i];
0248     int iSide = int(tk.eta() > 0);
0249     const auto &fts = trajectoryStateTransform::outerFreeState((tk), bFieldProd);
0250     // to the HGCal front
0251     const auto &tsos = prop.propagate(fts, firstDisk_[iSide]->surface());
0252     if (tsos.isValid()) {
0253       Vector trackP(tsos.globalPosition().x(), tsos.globalPosition().y(), tsos.globalPosition().z());
0254       trackPColl.emplace_back(trackP, i);
0255     }
0256     // to lastLayerEE
0257     const auto &tsos_int = prop.propagate(fts, interfaceDisk_[iSide]->surface());
0258     if (tsos_int.isValid()) {
0259       Vector trackP(tsos_int.globalPosition().x(), tsos_int.globalPosition().y(), tsos_int.globalPosition().z());
0260       tkPropIntColl.emplace_back(trackP, i);
0261     }
0262   }  // Tracks
0263   tkPropIntColl.shrink_to_fit();
0264   trackPColl.shrink_to_fit();
0265   candidateTrackIds.shrink_to_fit();
0266
0267   // Propagate tracksters
0268
0269   // Record postions of all tracksters propagated to layer 1 and lastLayerEE,
0270   // to be used later for distance calculation in the link finding stage
0271   // indexed by trackster index in event collection
0272   std::vector<Vector> tsAllProp;
0273   std::vector<Vector> tsAllPropInt;
0274   tsAllProp.reserve(tracksters.size());
0275   tsAllPropInt.reserve(tracksters.size());
0276   // Propagate tracksters
0277
0278   for (unsigned i = 0; i < tracksters.size(); ++i) {
0279     const auto &t = tracksters[i];
0280     if (TICLInterpretationAlgoBase::algo_verbosity_ > VerbosityLevel::Advanced)
0281       LogDebug("GeneralInterpretationAlgo")
0282           << "trackster " << i << " - eta " << t.barycenter().eta() << " phi " << t.barycenter().phi() << " time "
0283           << t.time() << " energy " << t.raw_energy() << "\n";
0284
0285     // to HGCal front
0286     float zVal = hgcons_->waferZ(1, true);
0287     auto tsP = propagateTrackster(t, i, zVal, tracksterPropTiles);
0288     tsAllProp.emplace_back(tsP);
0289
0290     // to lastLayerEE
0291     zVal = rhtools_.getPositionLayer(rhtools_.lastLayerEE()).z();
0292     tsP = propagateTrackster(t, i, zVal, tsPropIntTiles);
0293     tsAllPropInt.emplace_back(tsP);
0294
0295   }  // TS
0296
0297   // step 1: tracks -> all tracksters, at firstLayerEE
0298   std::vector<std::vector<unsigned>> tsNearTk(tracks.size());
0299   findTrackstersInWindow(
0300       tracksters, trackPColl, tracksterPropTiles, tsAllProp, del_tk_ts_layer1_, tracksters.size(), tsNearTk);
0301
0302   // step 2: tracks -> all tracksters, at lastLayerEE
0303   std::vector<std::vector<unsigned>> tsNearTkAtInt(tracks.size());
0304   findTrackstersInWindow(
0305       tracksters, tkPropIntColl, tsPropIntTiles, tsAllPropInt, del_tk_ts_int_, tracksters.size(), tsNearTkAtInt);
0306
0307   std::vector<unsigned int> chargedHadronsFromTk;
0308   std::vector<std::vector<unsigned int>> trackstersInTrackIndices;
0309   trackstersInTrackIndices.resize(tracks.size());
0310
0311   std::vector<bool> chargedMask(tracksters.size(), true);
0312   for (unsigned &i : candidateTrackIds) {
0313     if (tsNearTk[i].empty() && tsNearTkAtInt[i].empty()) {  // nothing linked to track, make charged hadrons
0314       continue;
0315     }
0316
0317     std::vector<unsigned int> chargedCandidate;
0318     float total_raw_energy = 0.f;
0319
0320     float track_time = 0.f;
0321     float track_timeErr = 0.f;
0322     float track_quality = 0.f;
0323     float track_beta = 0.f;
0324     GlobalPoint track_MtdPos{0.f, 0.f, 0.f};
0325     if (useMTDTiming) {
0326       auto const &inputTimingView = (*inputTiming_h).const_view();
0327       track_time = inputTimingView.time()[i];
0328       track_timeErr = inputTimingView.timeErr()[i];
0329       track_quality = inputTimingView.MVAquality()[i];
0330       track_beta = inputTimingView.beta()[i];
0331       track_MtdPos = {
0332           inputTimingView.posInMTD_x()[i], inputTimingView.posInMTD_y()[i], inputTimingView.posInMTD_z()[i]};
0333     }
0334
0335     for (auto const tsIdx : tsNearTk[i]) {
0336       if (chargedMask[tsIdx] && timeAndEnergyCompatible(total_raw_energy,
0337                                                         tracks[i],
0338                                                         tracksters[tsIdx],
0339                                                         track_time,
0340                                                         track_timeErr,
0341                                                         track_quality,
0342                                                         track_beta,
0343                                                         track_MtdPos,
0344                                                         useMTDTiming)) {
0345         chargedCandidate.push_back(tsIdx);
0346         chargedMask[tsIdx] = false;
0347         total_raw_energy += tracksters[tsIdx].raw_energy();
0348       }
0349     }
0350     for (const unsigned tsIdx : tsNearTkAtInt[i]) {  // do the same for tk -> ts links at the interface
0351       if (chargedMask[tsIdx] && timeAndEnergyCompatible(total_raw_energy,
0352                                                         tracks[i],
0353                                                         tracksters[tsIdx],
0354                                                         track_time,
0355                                                         track_timeErr,
0356                                                         track_quality,
0357                                                         track_beta,
0358                                                         track_MtdPos,
0359                                                         useMTDTiming)) {
0360         chargedCandidate.push_back(tsIdx);
0361         chargedMask[tsIdx] = false;
0362         total_raw_energy += tracksters[tsIdx].raw_energy();
0363       }
0364     }
0365     trackstersInTrackIndices[i] = chargedCandidate;
0366   }
0367
0368   for (size_t iTrack = 0; iTrack < trackstersInTrackIndices.size(); iTrack++) {
0369     if (!trackstersInTrackIndices[iTrack].empty()) {
0370       if (trackstersInTrackIndices[iTrack].size() == 1) {
0371         auto tracksterId = trackstersInTrackIndices[iTrack][0];
0372         resultCandidate[iTrack] = resultTracksters.size();
0373         resultTracksters.push_back(input.tracksters[tracksterId]);
0374       } else {
0375         // in this case mergeTracksters() clears the pid probabilities and the regressed energy is not set
0376         // TODO: fix probabilities when CNN will be splitted
0377         Trackster outTrackster;
0378         float regr_en = 0.f;
0379         bool isHadron = false;
0380         for (auto const tracksterId : trackstersInTrackIndices[iTrack]) {
0381           //maskTracksters[tracksterId] = 0;
0382           outTrackster.mergeTracksters(input.tracksters[tracksterId]);
0383           regr_en += input.tracksters[tracksterId].regressed_energy();
0384           if (input.tracksters[tracksterId].isHadronic())
0385             isHadron = true;
0386         }
0387         resultCandidate[iTrack] = resultTracksters.size();
0388         resultTracksters.push_back(outTrackster);
0389         resultTracksters.back().setRegressedEnergy(regr_en);
0390         // since a track has been linked it can only be electron or charged hadron
0391         if (isHadron)
0392           resultTracksters.back().setIdProbability(ticl::Trackster::ParticleType::charged_hadron, 1.f);
0393         else
0394           resultTracksters.back().setIdProbability(ticl::Trackster::ParticleType::electron, 1.f);
0395       }
0396     }
0397   }
0398
0399   for (size_t iTrackster = 0; iTrackster < input.tracksters.size(); iTrackster++) {
0400     if (chargedMask[iTrackster]) {
0401       resultTracksters.push_back(input.tracksters[iTrackster]);
0402     }
0403   }
0404 };
0405
0406 void GeneralInterpretationAlgo::fillPSetDescription(edm::ParameterSetDescription &desc) {
0407   desc.add<std::string>("cutTk",
0408                         "1.48 < abs(eta) < 3.0 && pt > 1. && quality(\"highPurity\") && "
0409                         "hitPattern().numberOfLostHits(\"MISSING_OUTER_HITS\") < 5");
0410   desc.add<double>("delta_tk_ts_layer1", 0.02);
0411   desc.add<double>("delta_tk_ts_interface", 0.03);
0412   desc.add<double>("timing_quality_threshold", 0.5);
0413   TICLInterpretationAlgoBase::fillPSetDescription(desc);
0414 }