Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-09-21 02:12:46 

0001 #include <cmath>
0002 #include <string>
0003 #include "RecoHGCal/TICL/plugins/LinkingAlgoByDirectionGeometric.h"
0004 
0005 #include "DataFormats/GeometrySurface/interface/BoundDisk.h"
0006 #include "DataFormats/HGCalReco/interface/Common.h"
0007 
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 
0010 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0011 
0012 #include "RecoParticleFlow/PFProducer/interface/PFMuonAlgo.h"
0013 
0014 using namespace ticl;
0015 
0016 LinkingAlgoByDirectionGeometric::LinkingAlgoByDirectionGeometric(const edm::ParameterSet &conf)
0017     : LinkingAlgoBase(conf),
0018       del_tk_ts_layer1_(conf.getParameter<double>("delta_tk_ts_layer1")),
0019       del_tk_ts_int_(conf.getParameter<double>("delta_tk_ts_interface")),
0020       del_ts_em_had_(conf.getParameter<double>("delta_ts_em_had")),
0021       del_ts_had_had_(conf.getParameter<double>("delta_ts_had_had")),
0022       timing_quality_threshold_(conf.getParameter<double>("track_time_quality_threshold")),
0023       cutTk_(conf.getParameter<std::string>("cutTk")) {}
0024 
0025 LinkingAlgoByDirectionGeometric::~LinkingAlgoByDirectionGeometric() {}
0026 
0027 void LinkingAlgoByDirectionGeometric::initialize(const HGCalDDDConstants *hgcons,
0028                                                  const hgcal::RecHitTools rhtools,
0029                                                  const edm::ESHandle<MagneticField> bfieldH,
0030                                                  const edm::ESHandle<Propagator> propH) {
0031   hgcons_ = hgcons;
0032   rhtools_ = rhtools;
0033   buildLayers();
0034 
0035   bfield_ = bfieldH;
0036   propagator_ = propH;
0037 }
0038 
0039 math::XYZVector LinkingAlgoByDirectionGeometric::propagateTrackster(const Trackster &t,
0040                                                                     const unsigned idx,
0041                                                                     float zVal,
0042                                                                     std::array<TICLLayerTile, 2> &tracksterTiles) {
0043   // needs only the positive Z co-ordinate of the surface to propagate to
0044   // the correct sign is calculated inside according to the barycenter of trackster
0045   Vector const &baryc = t.barycenter();
0046   Vector directnv = t.eigenvectors(0);
0047 
0048   // barycenter as direction for tracksters w/ poor PCA
0049   // propagation still done to get the cartesian coords
0050   // which are anyway converted to eta, phi in linking
0051   // -> can be simplified later
0052 
0053   //FP: disable PCA propagation for the moment and fallback to barycenter position
0054   // if (t.eigenvalues()[0] / t.eigenvalues()[1] < 20)
0055   directnv = baryc.unit();
0056 
0057   zVal *= (baryc.Z() > 0) ? 1 : -1;
0058 
0059   float par = (zVal - baryc.Z()) / directnv.Z();
0060   float xOnSurface = par * directnv.X() + baryc.X();
0061   float yOnSurface = par * directnv.Y() + baryc.Y();
0062   Vector tPoint(xOnSurface, yOnSurface, zVal);
0063   if (tPoint.Eta() > 0)
0064     tracksterTiles[1].fill(tPoint.Eta(), tPoint.Phi(), idx);
0065 
0066   else if (tPoint.Eta() < 0)
0067     tracksterTiles[0].fill(tPoint.Eta(), tPoint.Phi(), idx);
0068 
0069   return tPoint;
0070 }
0071 
0072 void LinkingAlgoByDirectionGeometric::findTrackstersInWindow(
0073     const std::vector<std::pair<Vector, unsigned>> &seedingCollection,
0074     const std::array<TICLLayerTile, 2> &tracksterTiles,
0075     const std::vector<Vector> &tracksterPropPoints,
0076     float delta,
0077     unsigned trackstersSize,
0078     std::vector<std::vector<unsigned>> &resultCollection,
0079     bool useMask = false) {
0080   // Finds tracksters in tracksterTiles within an eta-phi window
0081   // (given by delta) of the objects (track/trackster) in the seedingCollection.
0082   // Element i in resultCollection is the vector of trackster
0083   // indices found close to the i-th object in the seedingCollection.
0084   // If specified, Tracksters are masked once found as close to an object.
0085   std::vector<int> mask(trackstersSize, 0);
0086   float delta2 = delta * delta;
0087 
0088   for (auto &i : seedingCollection) {
0089     float seed_eta = i.first.Eta();
0090     float seed_phi = i.first.Phi();
0091     unsigned seedId = i.second;
0092     auto sideZ = seed_eta > 0;  //forward or backward region
0093     const TICLLayerTile &tile = tracksterTiles[sideZ];
0094     float eta_min = std::max(abs(seed_eta) - delta, (float)TileConstants::minEta);
0095     float eta_max = std::min(abs(seed_eta) + delta, (float)TileConstants::maxEta);
0096 
0097     // get range of bins touched by delta
0098     std::array<int, 4> search_box = tile.searchBoxEtaPhi(eta_min, eta_max, seed_phi - delta, seed_phi + delta);
0099 
0100     std::vector<unsigned> in_delta;
0101     std::vector<float> distances2;
0102     for (int eta_i = search_box[0]; eta_i <= search_box[1]; ++eta_i) {
0103       for (int phi_i = search_box[2]; phi_i <= search_box[3]; ++phi_i) {
0104         const auto &in_tile = tile[tile.globalBin(eta_i, (phi_i % TileConstants::nPhiBins))];
0105         for (const unsigned &t_i : in_tile) {
0106           // calculate actual distances of tracksters to the seed for a more accurate cut
0107           auto sep2 = (tracksterPropPoints[t_i].Eta() - seed_eta) * (tracksterPropPoints[t_i].Eta() - seed_eta) +
0108                       (tracksterPropPoints[t_i].Phi() - seed_phi) * (tracksterPropPoints[t_i].Phi() - seed_phi);
0109           if (sep2 < delta2) {
0110             in_delta.push_back(t_i);
0111             distances2.push_back(sep2);
0112           }
0113         }
0114       }
0115     }
0116 
0117     // sort tracksters found in ascending order of their distances from the seed
0118     std::vector<unsigned> indices(in_delta.size());
0119     std::iota(indices.begin(), indices.end(), 0);
0120     std::sort(indices.begin(), indices.end(), [&](unsigned i, unsigned j) { return distances2[i] < distances2[j]; });
0121 
0122     // push back sorted tracksters in the result collection
0123     for (const unsigned &index : indices) {
0124       const auto &t_i = in_delta[index];
0125       if (!mask[t_i]) {
0126         resultCollection[seedId].push_back(t_i);
0127         if (useMask)
0128           mask[t_i] = 1;
0129       }
0130     }
0131 
0132   }  // seeding collection loop
0133 }
0134 
0135 bool LinkingAlgoByDirectionGeometric::timeAndEnergyCompatible(float &total_raw_energy,
0136                                                               const reco::Track &track,
0137                                                               const Trackster &trackster,
0138                                                               const float &tkT,
0139                                                               const float &tkTErr,
0140                                                               const float &tkTimeQual) {
0141   float threshold = std::min(0.2 * trackster.raw_energy(), 10.0);
0142 
0143   bool energyCompatible = (total_raw_energy + trackster.raw_energy() < track.p() + threshold);
0144   // compatible if trackster time is within 3sigma of
0145   // track time; compatible if either: no time assigned
0146   // to trackster or track time quality is below threshold
0147   float tsT = trackster.time();
0148   float tsTErr = trackster.timeError();
0149 
0150   bool timeCompatible = false;
0151 
0152   if (tsT == -99. or tkTimeQual < timing_quality_threshold_)
0153     timeCompatible = true;
0154   else {
0155     timeCompatible = (std::abs(tsT - tkT) < maxDeltaT_ * sqrt(tsTErr * tsTErr + tkTErr * tkTErr));
0156   }
0157 
0158   if (LinkingAlgoBase::algo_verbosity_ > VerbosityLevel::Advanced) {
0159     if (!(energyCompatible))
0160       LogDebug("LinkingAlgoByDirectionGeometric")
0161           << "energy incompatible : track p " << track.p() << " trackster energy " << trackster.raw_energy() << "\n";
0162     if (!(timeCompatible))
0163       LogDebug("LinkingAlgoByDirectionGeometric") << "time incompatible : track time " << tkT << " +/- " << tkTErr
0164                                                   << " trackster time " << tsT << " +/- " << tsTErr << "\n";
0165   }
0166   return energyCompatible && timeCompatible;
0167 }
0168 
0169 void LinkingAlgoByDirectionGeometric::recordTrackster(const unsigned ts,  //trackster index
0170                                                       const std::vector<Trackster> &tracksters,
0171                                                       const edm::Handle<std::vector<Trackster>> tsH,
0172                                                       std::vector<unsigned> &ts_mask,
0173                                                       float &energy_in_candidate,
0174                                                       TICLCandidate &candidate) {
0175   if (ts_mask[ts])
0176     return;
0177   candidate.addTrackster(edm::Ptr<Trackster>(tsH, ts));
0178   ts_mask[ts] = 1;
0179   energy_in_candidate += tracksters[ts].raw_energy();
0180 }
0181 
0182 void LinkingAlgoByDirectionGeometric::dumpLinksFound(std::vector<std::vector<unsigned>> &resultCollection,
0183                                                      const char *label) const {
0184 #ifdef EDM_ML_DEBUG
0185   if (!(LinkingAlgoBase::algo_verbosity_ > VerbosityLevel::Advanced))
0186     return;
0187 
0188   LogDebug("LinkingAlgoByDirectionGeometric") << "All links found - " << label << "\n";
0189   LogDebug("LinkingAlgoByDirectionGeometric") << "(seed can either be a track or trackster depending on the step)\n";
0190   for (unsigned i = 0; i < resultCollection.size(); ++i) {
0191     LogDebug("LinkingAlgoByDirectionGeometric") << "seed " << i << " - tracksters : ";
0192     const auto &links = resultCollection[i];
0193     for (unsigned j = 0; j < links.size(); ++j) {
0194       LogDebug("LinkingAlgoByDirectionGeometric") << j;
0195     }
0196     LogDebug("LinkingAlgoByDirectionGeometric") << "\n";
0197   }
0198 #endif  // EDM_ML_DEBUG
0199 }
0200 
0201 void LinkingAlgoByDirectionGeometric::buildLayers() {
0202   // build disks at HGCal front & EM-Had interface for track propagation
0203 
0204   float zVal = hgcons_->waferZ(1, true);
0205   std::pair<float, float> rMinMax = hgcons_->rangeR(zVal, true);
0206 
0207   float zVal_interface = rhtools_.getPositionLayer(rhtools_.lastLayerEE()).z();
0208   std::pair<float, float> rMinMax_interface = hgcons_->rangeR(zVal_interface, true);
0209 
0210   for (int iSide = 0; iSide < 2; ++iSide) {
0211     float zSide = (iSide == 0) ? (-1. * zVal) : zVal;
0212     firstDisk_[iSide] =
0213         std::make_unique<GeomDet>(Disk::build(Disk::PositionType(0, 0, zSide),
0214                                               Disk::RotationType(),
0215                                               SimpleDiskBounds(rMinMax.first, rMinMax.second, zSide - 0.5, zSide + 0.5))
0216                                       .get());
0217 
0218     zSide = (iSide == 0) ? (-1. * zVal_interface) : zVal_interface;
0219     interfaceDisk_[iSide] = std::make_unique<GeomDet>(
0220         Disk::build(Disk::PositionType(0, 0, zSide),
0221                     Disk::RotationType(),
0222                     SimpleDiskBounds(rMinMax_interface.first, rMinMax_interface.second, zSide - 0.5, zSide + 0.5))
0223             .get());
0224   }
0225 }
0226 
0227 void LinkingAlgoByDirectionGeometric::linkTracksters(const edm::Handle<std::vector<reco::Track>> tkH,
0228                                                      const edm::ValueMap<float> &tkTime,
0229                                                      const edm::ValueMap<float> &tkTimeErr,
0230                                                      const edm::ValueMap<float> &tkTimeQual,
0231                                                      const std::vector<reco::Muon> &muons,
0232                                                      const edm::Handle<std::vector<Trackster>> tsH,
0233                                                      std::vector<TICLCandidate> &resultLinked,
0234                                                      std::vector<TICLCandidate> &chargedHadronsFromTk) {
0235   const auto &tracks = *tkH;
0236   const auto &tracksters = *tsH;
0237 
0238   auto bFieldProd = bfield_.product();
0239   const Propagator &prop = (*propagator_);
0240 
0241   // propagated point collections
0242   // elements in the propagated points collecions are used
0243   // to look for potential linkages in the appropriate tiles
0244   std::vector<std::pair<Vector, unsigned>> trackPColl;     // propagated track points and index of track in collection
0245   std::vector<std::pair<Vector, unsigned>> tkPropIntColl;  // tracks propagated to lastLayerEE
0246   std::vector<std::pair<Vector, unsigned>> tsPropIntColl;  // Tracksters in CE-E, propagated to lastLayerEE
0247   std::vector<std::pair<Vector, unsigned>> tsHadPropIntColl;  // Tracksters in CE-H, propagated to lastLayerEE
0248   trackPColl.reserve(tracks.size());
0249   tkPropIntColl.reserve(tracks.size());
0250   tsPropIntColl.reserve(tracksters.size());
0251   tsHadPropIntColl.reserve(tracksters.size());
0252   // tiles, element 0 is bw, 1 is fw
0253   std::array<TICLLayerTile, 2> tracksterPropTiles = {};  // all Tracksters, propagated to layer 1
0254   std::array<TICLLayerTile, 2> tsPropIntTiles = {};      // all Tracksters, propagated to lastLayerEE
0255   std::array<TICLLayerTile, 2> tsHadPropIntTiles = {};   // Tracksters in CE-H, propagated to lastLayerEE
0256 
0257   // linking : trackster is hadronic if its barycenter is in CE-H
0258   auto isHadron = [&](const Trackster &t) -> bool {
0259     auto boundary_z = rhtools_.getPositionLayer(rhtools_.lastLayerEE()).z();
0260     return (std::abs(t.barycenter().Z()) > boundary_z);
0261   };
0262 
0263   if (LinkingAlgoBase::algo_verbosity_ > VerbosityLevel::Advanced)
0264     LogDebug("LinkingAlgoByDirectionGeometric") << "------- Geometric Linking ------- \n";
0265 
0266   // Propagate tracks
0267   std::vector<unsigned> candidateTrackIds;
0268   candidateTrackIds.reserve(tracks.size());
0269   for (unsigned i = 0; i < tracks.size(); ++i) {
0270     const auto &tk = tracks[i];
0271     reco::TrackRef trackref = reco::TrackRef(tkH, i);
0272 
0273     // veto tracks associated to muons
0274     int muId = PFMuonAlgo::muAssocToTrack(trackref, muons);
0275 
0276     if (LinkingAlgoBase::algo_verbosity_ > VerbosityLevel::Advanced)
0277       LogDebug("LinkingAlgoByDirectionGeometric")
0278           << "track " << i << " - eta " << tk.eta() << " phi " << tk.phi() << " time " << tkTime[reco::TrackRef(tkH, i)]
0279           << " time qual " << tkTimeQual[reco::TrackRef(tkH, i)] << "  muid " << muId << "\n";
0280 
0281     if (!cutTk_((tk)) or muId != -1)
0282       continue;
0283 
0284     // record tracks that can be used to make a ticlcandidate
0285     candidateTrackIds.push_back(i);
0286 
0287     // don't consider tracks below 2 GeV for linking
0288     if (std::sqrt(tk.p() * tk.p() + ticl::mpion2) < tkEnergyCut_)
0289       continue;
0290 
0291     int iSide = int(tk.eta() > 0);
0292     const auto &fts = trajectoryStateTransform::outerFreeState((tk), bFieldProd);
0293     // to the HGCal front
0294     const auto &tsos = prop.propagate(fts, firstDisk_[iSide]->surface());
0295     if (tsos.isValid()) {
0296       Vector trackP(tsos.globalPosition().x(), tsos.globalPosition().y(), tsos.globalPosition().z());
0297       trackPColl.emplace_back(trackP, i);
0298     }
0299     // to lastLayerEE
0300     const auto &tsos_int = prop.propagate(fts, interfaceDisk_[iSide]->surface());
0301     if (tsos_int.isValid()) {
0302       Vector trackP(tsos_int.globalPosition().x(), tsos_int.globalPosition().y(), tsos_int.globalPosition().z());
0303       tkPropIntColl.emplace_back(trackP, i);
0304     }
0305   }  // Tracks
0306   tkPropIntColl.shrink_to_fit();
0307   trackPColl.shrink_to_fit();
0308   candidateTrackIds.shrink_to_fit();
0309 
0310   // Propagate tracksters
0311 
0312   // Record postions of all tracksters propagated to layer 1 and lastLayerEE,
0313   // to be used later for distance calculation in the link finding stage
0314   // indexed by trackster index in event collection
0315   std::vector<Vector> tsAllProp;
0316   std::vector<Vector> tsAllPropInt;
0317   tsAllProp.reserve(tracksters.size());
0318   tsAllPropInt.reserve(tracksters.size());
0319 
0320   for (unsigned i = 0; i < tracksters.size(); ++i) {
0321     const auto &t = tracksters[i];
0322     if (LinkingAlgoBase::algo_verbosity_ > VerbosityLevel::Advanced)
0323       LogDebug("LinkingAlgoByDirectionGeometric")
0324           << "trackster " << i << " - eta " << t.barycenter().eta() << " phi " << t.barycenter().phi() << " time "
0325           << t.time() << " energy " << t.raw_energy() << "\n";
0326 
0327     // to HGCal front
0328     float zVal = hgcons_->waferZ(1, true);
0329     auto tsP = propagateTrackster(t, i, zVal, tracksterPropTiles);
0330     tsAllProp.emplace_back(tsP);
0331 
0332     // to lastLayerEE
0333     zVal = rhtools_.getPositionLayer(rhtools_.lastLayerEE()).z();
0334     tsP = propagateTrackster(t, i, zVal, tsPropIntTiles);
0335     tsAllPropInt.emplace_back(tsP);
0336 
0337     if (!isHadron(t))  // EM tracksters
0338       tsPropIntColl.emplace_back(tsP, i);
0339     else {  // HAD
0340       tsHadPropIntTiles[(t.barycenter().Z() > 0) ? 1 : 0].fill(tsP.Eta(), tsP.Phi(), i);
0341       tsHadPropIntColl.emplace_back(tsP, i);
0342     }
0343   }  // TS
0344   tsPropIntColl.shrink_to_fit();
0345   tsHadPropIntColl.shrink_to_fit();
0346 
0347   // Track - Trackster link finding
0348   // step 3: tracks -> all tracksters, at layer 1
0349 
0350   std::vector<std::vector<unsigned>> tsNearTk(tracks.size());
0351   findTrackstersInWindow(trackPColl, tracksterPropTiles, tsAllProp, del_tk_ts_layer1_, tracksters.size(), tsNearTk);
0352 
0353   // step 4: tracks -> all tracksters, at lastLayerEE
0354 
0355   std::vector<std::vector<unsigned>> tsNearTkAtInt(tracks.size());
0356   findTrackstersInWindow(tkPropIntColl, tsPropIntTiles, tsAllPropInt, del_tk_ts_int_, tracksters.size(), tsNearTkAtInt);
0357 
0358   // Trackster - Trackster link finding
0359   // step 2: tracksters EM -> HAD, at lastLayerEE
0360 
0361   std::vector<std::vector<unsigned>> tsNearAtInt(tracksters.size());
0362   findTrackstersInWindow(
0363       tsPropIntColl, tsHadPropIntTiles, tsAllPropInt, del_ts_em_had_, tracksters.size(), tsNearAtInt);
0364 
0365   // step 1: tracksters HAD -> HAD, at lastLayerEE
0366 
0367   std::vector<std::vector<unsigned>> tsHadNearAtInt(tracksters.size());
0368   findTrackstersInWindow(
0369       tsHadPropIntColl, tsHadPropIntTiles, tsAllPropInt, del_ts_had_had_, tracksters.size(), tsHadNearAtInt);
0370 
0371 #ifdef EDM_ML_DEBUG
0372   dumpLinksFound(tsNearTk, "track -> tracksters at layer 1");
0373   dumpLinksFound(tsNearTkAtInt, "track -> tracksters at lastLayerEE");
0374   dumpLinksFound(tsNearAtInt, "EM -> HAD tracksters at lastLayerEE");
0375   dumpLinksFound(tsHadNearAtInt, "HAD -> HAD tracksters at lastLayerEE");
0376 #endif  //EDM_ML_DEBUG
0377 
0378   // make final collections
0379 
0380   std::vector<TICLCandidate> chargedCandidates;
0381   std::vector<unsigned int> chargedMask(tracksters.size(), 0);
0382   for (unsigned &i : candidateTrackIds) {
0383     if (tsNearTk[i].empty() && tsNearTkAtInt[i].empty()) {  // nothing linked to track, make charged hadrons
0384       TICLCandidate chargedHad;
0385       chargedHad.setTrackPtr(edm::Ptr<reco::Track>(tkH, i));
0386       chargedHadronsFromTk.push_back(chargedHad);
0387       continue;
0388     }
0389 
0390     TICLCandidate chargedCandidate;
0391     float total_raw_energy = 0.;
0392 
0393     auto tkRef = reco::TrackRef(tkH, i);
0394     auto track_time = tkTime[tkRef];
0395     auto track_timeErr = tkTimeErr[tkRef];
0396     auto track_timeQual = tkTimeQual[tkRef];
0397 
0398     for (const unsigned ts3_idx : tsNearTk[i]) {  // tk -> ts
0399       if (timeAndEnergyCompatible(
0400               total_raw_energy, tracks[i], tracksters[ts3_idx], track_time, track_timeErr, track_timeQual)) {
0401         recordTrackster(ts3_idx, tracksters, tsH, chargedMask, total_raw_energy, chargedCandidate);
0402       }
0403       for (const unsigned ts2_idx : tsNearAtInt[ts3_idx]) {  // ts_EM -> ts_HAD
0404         if (timeAndEnergyCompatible(
0405                 total_raw_energy, tracks[i], tracksters[ts2_idx], track_time, track_timeErr, track_timeQual)) {
0406           recordTrackster(ts2_idx, tracksters, tsH, chargedMask, total_raw_energy, chargedCandidate);
0407         }
0408         for (const unsigned ts1_idx : tsHadNearAtInt[ts2_idx]) {  // ts_HAD -> ts_HAD
0409           if (timeAndEnergyCompatible(
0410                   total_raw_energy, tracks[i], tracksters[ts1_idx], track_time, track_timeErr, track_timeQual)) {
0411             recordTrackster(ts1_idx, tracksters, tsH, chargedMask, total_raw_energy, chargedCandidate);
0412           }
0413         }
0414       }
0415       for (const unsigned ts1_idx : tsHadNearAtInt[ts3_idx]) {  // ts_HAD -> ts_HAD
0416         if (timeAndEnergyCompatible(
0417                 total_raw_energy, tracks[i], tracksters[ts1_idx], track_time, track_timeErr, track_timeQual)) {
0418           recordTrackster(ts1_idx, tracksters, tsH, chargedMask, total_raw_energy, chargedCandidate);
0419         }
0420       }
0421     }
0422     for (const unsigned ts4_idx : tsNearTkAtInt[i]) {  // do the same for tk -> ts links at the interface
0423       if (timeAndEnergyCompatible(
0424               total_raw_energy, tracks[i], tracksters[ts4_idx], track_time, track_timeErr, track_timeQual)) {
0425         recordTrackster(ts4_idx, tracksters, tsH, chargedMask, total_raw_energy, chargedCandidate);
0426       }
0427       for (const unsigned ts2_idx : tsNearAtInt[ts4_idx]) {
0428         if (timeAndEnergyCompatible(
0429                 total_raw_energy, tracks[i], tracksters[ts2_idx], track_time, track_timeErr, track_timeQual)) {
0430           recordTrackster(ts2_idx, tracksters, tsH, chargedMask, total_raw_energy, chargedCandidate);
0431         }
0432         for (const unsigned ts1_idx : tsHadNearAtInt[ts2_idx]) {
0433           if (timeAndEnergyCompatible(
0434                   total_raw_energy, tracks[i], tracksters[ts1_idx], track_time, track_timeErr, track_timeQual)) {
0435             recordTrackster(ts1_idx, tracksters, tsH, chargedMask, total_raw_energy, chargedCandidate);
0436           }
0437         }
0438       }
0439       for (const unsigned ts1_idx : tsHadNearAtInt[ts4_idx]) {
0440         if (timeAndEnergyCompatible(
0441                 total_raw_energy, tracks[i], tracksters[ts1_idx], track_time, track_timeErr, track_timeQual)) {
0442           recordTrackster(ts1_idx, tracksters, tsH, chargedMask, total_raw_energy, chargedCandidate);
0443         }
0444       }
0445     }
0446 
0447     // do not create a candidate if no tracksters were added to candidate
0448     // can happen if all the tracksters linked to that track were already masked
0449     if (!chargedCandidate.tracksters().empty()) {
0450       chargedCandidate.setTrackPtr(edm::Ptr<reco::Track>(tkH, i));
0451       chargedCandidates.push_back(chargedCandidate);
0452     } else {  // create charged hadron
0453       TICLCandidate chargedHad;
0454       chargedHad.setTrackPtr(edm::Ptr<reco::Track>(tkH, i));
0455       chargedHadronsFromTk.push_back(chargedHad);
0456     }
0457   }
0458 
0459   std::vector<TICLCandidate> neutralCandidates;
0460   std::vector<int> neutralMask(tracksters.size(), 0);
0461   for (unsigned i = 0; i < tracksters.size(); ++i) {
0462     if (chargedMask[i])
0463       continue;
0464 
0465     TICLCandidate neutralCandidate;
0466     if (tsNearAtInt[i].empty() && tsHadNearAtInt[i].empty() && !neutralMask[i]) {  // nothing linked to this ts
0467       neutralCandidate.addTrackster(edm::Ptr<Trackster>(tsH, i));
0468       neutralMask[i] = 1;
0469       neutralCandidates.push_back(neutralCandidate);
0470       continue;
0471     }
0472     if (!neutralMask[i]) {
0473       neutralCandidate.addTrackster(edm::Ptr<Trackster>(tsH, i));
0474       neutralMask[i] = 1;
0475     }
0476     for (const unsigned ts2_idx : tsNearAtInt[i]) {
0477       if (chargedMask[ts2_idx])
0478         continue;
0479       if (!neutralMask[ts2_idx]) {
0480         neutralCandidate.addTrackster(edm::Ptr<Trackster>(tsH, ts2_idx));
0481         neutralMask[ts2_idx] = 1;
0482       }
0483       for (const unsigned ts1_idx : tsHadNearAtInt[ts2_idx]) {
0484         if (chargedMask[ts1_idx])
0485           continue;
0486         if (!neutralMask[ts1_idx]) {
0487           neutralCandidate.addTrackster(edm::Ptr<Trackster>(tsH, ts1_idx));
0488           neutralMask[ts1_idx] = 1;
0489         }
0490       }
0491     }
0492     for (const unsigned ts1_idx : tsHadNearAtInt[i]) {
0493       if (chargedMask[ts1_idx])
0494         continue;
0495       if (!neutralMask[ts1_idx]) {
0496         neutralCandidate.addTrackster(edm::Ptr<Trackster>(tsH, ts1_idx));
0497         neutralMask[ts1_idx] = 1;
0498       }
0499     }
0500     // filter empty candidates
0501     if (!neutralCandidate.tracksters().empty()) {
0502       neutralCandidates.push_back(neutralCandidate);
0503     }
0504   }
0505 
0506   resultLinked.insert(std::end(resultLinked), std::begin(neutralCandidates), std::end(neutralCandidates));
0507   resultLinked.insert(std::end(resultLinked), std::begin(chargedCandidates), std::end(chargedCandidates));
0508 
0509 }  // linkTracksters
0510 
0511 void LinkingAlgoByDirectionGeometric::fillPSetDescription(edm::ParameterSetDescription &desc) {
0512   desc.add<std::string>("cutTk",
0513                         "1.48 < abs(eta) < 3.0 && pt > 1. && quality(\"highPurity\") && "
0514                         "hitPattern().numberOfLostHits(\"MISSING_OUTER_HITS\") < 5");
0515   desc.add<double>("delta_tk_ts_layer1", 0.02);
0516   desc.add<double>("delta_tk_ts_interface", 0.03);
0517   desc.add<double>("delta_ts_em_had", 0.03);
0518   desc.add<double>("delta_ts_had_had", 0.03);
0519   desc.add<double>("track_time_quality_threshold", 0.5);
0520   LinkingAlgoBase::fillPSetDescription(desc);
0521 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!