Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:25:12 

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