File indexing completed on 2024-09-07 04:37:30
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
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
0060
0061 Vector const &baryc = t.barycenter();
0062 Vector directnv = t.eigenvectors(0);
0063
0064
0065
0066
0067
0068
0069
0070
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
0095
0096
0097
0098
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;
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
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
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
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
0127 energies.push_back(tracksters[t_i].raw_energy());
0128 }
0129 }
0130 }
0131 }
0132
0133
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
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 }
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
0167
0168
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
0187
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
0219
0220
0221 std::vector<std::pair<Vector, unsigned>> trackPColl;
0222 std::vector<std::pair<Vector, unsigned>> tkPropIntColl;
0223
0224 trackPColl.reserve(tracks.size());
0225 tkPropIntColl.reserve(tracks.size());
0226
0227 std::array<TICLLayerTile, 2> tracksterPropTiles = {};
0228 std::array<TICLLayerTile, 2> tsPropIntTiles = {};
0229
0230 if (TICLInterpretationAlgoBase::algo_verbosity_ > VerbosityLevel::Advanced)
0231 LogDebug("GeneralInterpretationAlgo") << "------- Geometric Linking ------- \n";
0232
0233
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
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
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 }
0263 tkPropIntColl.shrink_to_fit();
0264 trackPColl.shrink_to_fit();
0265 candidateTrackIds.shrink_to_fit();
0266
0267
0268
0269
0270
0271
0272 std::vector<Vector> tsAllProp;
0273 std::vector<Vector> tsAllPropInt;
0274 tsAllProp.reserve(tracksters.size());
0275 tsAllPropInt.reserve(tracksters.size());
0276
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
0286 float zVal = hgcons_->waferZ(1, true);
0287 auto tsP = propagateTrackster(t, i, zVal, tracksterPropTiles);
0288 tsAllProp.emplace_back(tsP);
0289
0290
0291 zVal = rhtools_.getPositionLayer(rhtools_.lastLayerEE()).z();
0292 tsP = propagateTrackster(t, i, zVal, tsPropIntTiles);
0293 tsAllPropInt.emplace_back(tsP);
0294
0295 }
0296
0297
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
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()) {
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]) {
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
0376
0377 Trackster outTrackster;
0378 float regr_en = 0.f;
0379 bool isHadron = false;
0380 for (auto const tracksterId : trackstersInTrackIndices[iTrack]) {
0381
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
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 }