71 void assignTimeToCandidates(std::vector<TICLCandidate> &resultCandidates,
106 const float tkEnergyCut_ = 2.0f;
112 static constexpr float c_light_ = CLHEP::c_light * CLHEP::ns / CLHEP::cm;
123 useMTDTiming_(ps.getParameter<
bool>(
"useMTDTiming")),
124 useTimingAverage_(ps.getParameter<
bool>(
"useTimingAverage")),
125 timingQualityThreshold_(ps.getParameter<double>(
"timingQualityThreshold")),
128 detector_(ps.getParameter<
std::
string>(
"detector")),
129 propName_(ps.getParameter<
std::
string>(
"propagator")),
132 trackingGeometry_token_(
134 cutTk_(ps.getParameter<
std::
string>(
"cutTk")) {
136 for (
auto const &
tag : ps.
getParameter<std::vector<edm::InputTag>>(
"egamma_tracksters_collections")) {
141 for (
auto const &
tag : ps.
getParameter<std::vector<edm::InputTag>>(
"egamma_tracksterlinks_collections")) {
150 for (
auto const &
tag : ps.
getParameter<std::vector<edm::InputTag>>(
"general_tracksters_collections")) {
154 for (
auto const &
tag : ps.
getParameter<std::vector<edm::InputTag>>(
"general_tracksterlinks_collections")) {
162 for (
auto const &
tag : ps.
getParameter<std::vector<edm::InputTag>>(
"original_masks")) {
166 std::string detectorName_ = (
detector_ ==
"HFNose") ?
"HGCalHFNoseSensitive" :
"HGCalEESensitive";
168 esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(
edm::ESInputTag(
"", detectorName_));
173 produces<std::vector<TICLCandidate>>();
176 produces<std::vector<Trackster>>();
199 const edm::Handle<std::vector<reco::Muon>> &muons_h,
201 const float tkEnergyCut_,
202 std::vector<bool> &maskTracks) {
203 auto const &
tracks = *tkH;
204 for (
unsigned i = 0;
i <
tracks.size(); ++
i) {
212 if (!cutTk_((tk))
or (muId != -1 and
PFMuonAlgo::isMuon(muonref) and not(*muons_h)[muId].isTrackerMuon())) {
213 maskTracks[
i] =
false;
219 maskTracks[
i] =
false;
224 maskTracks[
i] =
true;
229 auto resultTracksters = std::make_unique<std::vector<Trackster>>();
230 auto resultTrackstersMerged = std::make_unique<std::vector<Trackster>>();
231 auto linkedResultTracksters = std::make_unique<std::vector<std::vector<unsigned int>>>();
240 const auto &
tracks = *tracks_h;
246 inputTimingView = (*inputTiming_h).const_view();
254 std::vector<float> original_global_mask(
layerClusters.size(), 1.f);
257 for (
unsigned int j = 0;
j < tmp_mask.size(); ++
j) {
258 original_global_mask[
j] *= tmp_mask[
j];
262 auto resultMask = std::make_unique<std::vector<float>>(original_global_mask);
269 generalTrackstersManager.addVector(*general_tracksters_h[
i]);
272 std::vector<edm::Handle<std::vector<std::vector<unsigned>>>> general_tracksterlinks_h(
274 std::vector<std::vector<unsigned>> generalTracksterLinksGlobalId;
277 for (
unsigned int j = 0;
j < general_tracksterlinks_h[
i]->size(); ++
j) {
278 generalTracksterLinksGlobalId.emplace_back();
279 auto &links_vector = generalTracksterLinksGlobalId.back();
280 links_vector.resize((*general_tracksterlinks_h[
i])[
j].
size());
281 for (
unsigned int k = 0;
k < links_vector.size(); ++
k) {
282 links_vector[
k] = generalTrackstersManager.getGlobalIndex(
i, (*general_tracksterlinks_h[
i])[
j][
k]);
287 std::vector<bool> maskTracks;
288 maskTracks.resize(
tracks.size());
295 generalTrackstersManager,
296 generalTracksterLinksGlobalId,
300 auto resultCandidates = std::make_unique<std::vector<TICLCandidate>>();
301 std::vector<int> trackstersInTrackIndices(
tracks.size(), -1);
312 std::vector<bool> maskTracksters(resultTracksters->size(),
true);
315 for (
size_t iTrack = 0; iTrack <
tracks.size(); iTrack++) {
316 if (maskTracks[iTrack]) {
317 auto const tracksterId = trackstersInTrackIndices[iTrack];
319 if (tracksterId != -1 and !maskTracksters.empty()) {
322 resultCandidates->push_back(chargedCandidate);
323 maskTracksters[tracksterId] =
false;
331 if (muonRef.
isNonnull() and muonRef->isGlobalMuon()) {
333 chargedCandidate.
setPdgId(13 * trackPtr.get()->charge());
335 resultCandidates->push_back(chargedCandidate);
341 for (
size_t iTrackster = 0; iTrackster < maskTracksters.size(); iTrackster++) {
342 if (maskTracksters[iTrackster]) {
346 resultCandidates->push_back(neutralCandidate);
360 float pathlength =
propagator->propagateWithPath(stateAtBeamspot, surf_inn.surface()).
second;
363 const auto &t_inn_out =
propagator->propagateWithPath(fts_inn, surf_out.surface());
365 if (t_inn_out.first.isValid()) {
366 pathlength += t_inn_out.second;
368 std::pair<float, float> rMinMax =
hgcons_->
rangeR(zVal,
true);
371 float zSide = (iSide == 0) ? (-1. * zVal) : zVal;
372 const auto &disk = std::make_unique<GeomDet>(
377 const auto &tsos =
propagator->propagateWithPath(fts_out, disk->surface());
379 if (tsos.first.isValid()) {
380 pathlength += tsos.second;
386 <<
"Not able to use the track to compute the path length. A straight line will be used instead.";
395 template <
typename F>
400 for (
auto &
cand : resultCandidates) {
403 float invTimeErr = 0.f;
404 float timeErr = -1.f;
406 for (
const auto &tr :
cand.tracksters()) {
407 if (tr->timeError() > 0) {
408 const auto invTimeESq =
pow(tr->timeError(), -2);
409 const auto x = tr->barycenter().X();
410 const auto y = tr->barycenter().Y();
411 const auto z = tr->barycenter().Z();
413 if (
cand.trackPtr().get() !=
nullptr) {
415 if (
useMTDTiming_ and inputTimingView.timeErr()[trackIndex] > 0) {
416 const auto xMtd = inputTimingView.posInMTD_x()[trackIndex];
417 const auto yMtd = inputTimingView.posInMTD_y()[trackIndex];
418 const auto zMtd = inputTimingView.posInMTD_z()[trackIndex];
420 beta = inputTimingView.beta()[trackIndex];
421 path =
std::sqrt((
x - xMtd) * (
x - xMtd) + (
y - yMtd) * (
y - yMtd) + (
z - zMtd) * (
z - zMtd)) +
422 inputTimingView.pathLength()[trackIndex];
424 float pathLength =
func(*(
cand.trackPtr().get()),
z);
431 invTimeErr += invTimeESq;
434 if (invTimeErr > 0) {
437 timeErr =
sqrt(1.
f / invTimeErr);
448 const auto timeHGC =
cand.time();
449 const auto timeEHGC =
cand.timeError();
450 const auto timeMTD = inputTimingView.time0()[trackIndex];
451 const auto timeEMTD = inputTimingView.time0Err()[trackIndex];
455 const auto invTimeESqHGC =
pow(timeEHGC, -2);
456 const auto invTimeESqMTD =
pow(timeEMTD, -2);
457 timeErr = 1.f / (invTimeESqHGC + invTimeESqMTD);
458 time = (timeHGC * invTimeESqHGC + timeMTD * invTimeESqMTD) * timeErr;
459 timeErr =
sqrt(timeErr);
460 }
else if (timeEMTD > 0) {
466 cand.setMTDTime(inputTimingView.time()[trackIndex], inputTimingView.timeErr()[trackIndex]);
476 desc.add<std::vector<edm::InputTag>>(
"egamma_tracksters_collections", {
edm::InputTag(
"ticlTracksterLinks")});
477 desc.add<std::vector<edm::InputTag>>(
"egamma_tracksterlinks_collections", {
edm::InputTag(
"ticlTracksterLinks")});
478 desc.add<std::vector<edm::InputTag>>(
"general_tracksters_collections", {
edm::InputTag(
"ticlTracksterLinks")});
479 desc.add<std::vector<edm::InputTag>>(
"general_tracksterlinks_collections", {
edm::InputTag(
"ticlTracksterLinks")});
480 desc.add<std::vector<edm::InputTag>>(
"original_masks",
481 {
edm::InputTag(
"hgcalMergeLayerClusters",
"InitialLayerClustersMask")});
489 desc.add<
bool>(
"useMTDTiming",
true);
490 desc.add<
bool>(
"useTimingAverage",
true);
491 desc.add<
double>(
"timingQualityThreshold", 0.5);
493 "1.48 < abs(eta) < 3.0 && pt > 1. && quality(\"highPurity\") && " 494 "hitPattern().numberOfLostHits(\"MISSING_OUTER_HITS\") < 5");
495 descriptions.
add(
"ticlCandidateProducer",
desc);
static DiskPointer build(Args &&... args)
static int muAssocToTrack(const reco::TrackRef &trackref, const reco::MuonCollection &muons)
TICLCandidateProducer(const edm::ParameterSet &ps)
T getParameter(std::string const &) const
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
bool get(ProductID const &oid, Handle< PROD > &result) const
void beginRun(edm::Run const &iEvent, edm::EventSetup const &es) override
const StringCutObjectSelector< reco::Track > cutTk_
static bool isMuon(const reco::PFBlockElement &elt)
std::once_flag initializeGeometry_
const edm::EDGetTokenT< edm::ValueMap< std::pair< float, float > > > clustersTime_token_
void assignTimeToCandidates(std::vector< TICLCandidate > &resultCandidates, edm::Handle< std::vector< reco::Track >> track_h, MtdHostCollection::ConstView &inputTimingView, F func) const
Global3DPoint GlobalPoint
ParameterDescriptionNode * addNode(ParameterDescriptionNode const &node)
bool isNonnull() const
Checks for non-null.
bool getByToken(EDGetToken token, Handle< PROD > &result) const
const edm::ESGetToken< MagneticField, IdealMagneticFieldRecord > bfield_token_
const edm::EDGetTokenT< std::vector< reco::Track > > tracks_token_
static std::string const input
const float timingQualityThreshold_
edm::ESHandle< Propagator > propagator_
U second(std::pair< T, U > const &p)
~TICLCandidateProducer() override
const edm::ESGetToken< Propagator, TrackingComponentsRecord > propagator_token_
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
const bool useTimingAverage_
edm::ESHandle< GlobalTrackingGeometry > trackingGeometry_
T const * product() const
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e< void, edm::EventID const &, edm::Timestamp const & > We also list in braces which AR_WATCH_USING_METHOD_ is used for those or
edm::Ref< MuonCollection > MuonRef
presistent reference to a Muon
#define DEFINE_FWK_MODULE(type)
edm::ESGetToken< HGCalDDDConstants, IdealGeometryRecord > hdc_token_
void produce(edm::Event &, const edm::EventSetup &) override
ESHandle< T > getHandle(const ESGetToken< T, R > &iToken) const
std::vector< edm::EDGetTokenT< std::vector< Trackster > > > general_tracksters_tokens_
std::pair< double, double > rangeR(double z, bool reco) const
hgcal::RecHitTools rhtools_
typename Layout::ConstView ConstView
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
std::vector< edm::EDGetTokenT< std::vector< Trackster > > > egamma_tracksters_tokens_
edm::Ref< TrackCollection > TrackRef
persistent reference to a Track
void add(std::string const &label, ParameterSetDescription const &psetDescription)
const std::string detector_
const edm::ESGetToken< GlobalTrackingGeometry, GlobalTrackingGeometryRecord > trackingGeometry_token_
const edm::EDGetTokenT< std::vector< reco::CaloCluster > > clusters_token_
const HGCalDDDConstants * hgcons_
const std::string propName_
static constexpr float c_light_
std::vector< edm::EDGetTokenT< std::vector< std::vector< unsigned > > > > general_tracksterlinks_tokens_
const edm::ESGetToken< CaloGeometry, CaloGeometryRecord > geometry_token_
void assignPCAtoTracksters(std::vector< Trackster > &, const std::vector< reco::CaloCluster > &, const edm::ValueMap< std::pair< float, float >> &, double, bool computeLocalTime=false, bool energyWeight=true)
edm::EDGetTokenT< MtdHostCollection > inputTimingToken_
std::vector< edm::EDGetTokenT< std::vector< float > > > original_masks_tokens_
std::unique_ptr< TICLInterpretationAlgoBase< reco::Track > > generalInterpretationAlgo_
Log< level::Warning, false > LogWarning
static uInt32 F(BLOWFISH_CTX *ctx, uInt32 x)
edm::ESHandle< MagneticField > bfield_
static constexpr float timeRes
void setPdgId(int pdgId) final
std::vector< edm::EDGetTokenT< std::vector< std::vector< unsigned > > > > egamma_tracksterlinks_tokens_
void filterTracks(edm::Handle< std::vector< reco::Track >> tkH, const edm::Handle< std::vector< reco::Muon >> &muons_h, const StringCutObjectSelector< reco::Track > cutTk_, const float tkEnergyCut_, std::vector< bool > &maskTracks)
Power< A, B >::type pow(const A &a, const B &b)
const edm::EDGetTokenT< std::vector< reco::Muon > > muons_token_
Global3DVector GlobalVector