Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Timing propagation in HGCAL RECO #28740

Merged
merged 22 commits into from Feb 4, 2020
Merged
Show file tree
Hide file tree
Changes from 17 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
2 changes: 2 additions & 0 deletions DataFormats/Common/src/classes_def.xml
Expand Up @@ -146,11 +146,13 @@
<class name="edm::ValueMap<bool>" />
<class name="edm::ValueMap<float>" />
<class name="edm::ValueMap<double>" />
<class name="edm::ValueMap<std::pair<float, float>>" />
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should this go in DataFormats/StdDictionaries?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I had missed that this was a ValueMap. It definitely should stay here.

<class name="edm::Wrapper<edm::ValueMap<int> >" splitLevel="0"/>
<class name="edm::Wrapper<edm::ValueMap<bool> >" splitLevel="0"/>
<class name="edm::Wrapper<edm::ValueMap<unsigned int> >" splitLevel="0"/>
<class name="edm::Wrapper<edm::ValueMap<float> >" splitLevel="0"/>
<class name="edm::Wrapper<edm::ValueMap<double> >" splitLevel="0"/>
<class name="edm::Wrapper<edm::ValueMap<std::pair<float, float>> >" splitLevel="0"/>
<class name="std::vector<edm::EventAuxiliary>"/>
<class name="edm::Wrapper<std::vector<edm::EventAuxiliary> >"/>
<class name="edm::Wrapper<std::vector<edm::ErrorSummaryEntry> >"/>
Expand Down
3 changes: 2 additions & 1 deletion DataFormats/HGCRecHit/interface/HGCRecHit.h
Expand Up @@ -89,7 +89,7 @@ class HGCRecHit : public CaloRecHit {
void setOutOfTimeEnergy(float energy);
void setSignalOverSigmaNoise(float sOverNoise);

void setTimeError(uint8_t timeErrBits);
void setTimeError(float timeErr);

/// set the flags (from Flags or ESFlags)
void setFlag(int flag) { flagBits_ |= (0x1 << flag); }
Expand All @@ -105,6 +105,7 @@ class HGCRecHit : public CaloRecHit {
/// store rechit condition (see Flags enum) in a bit-wise way
uint32_t flagBits_;
uint8_t signalOverSigmaNoise_;
float timeError_;
};

std::ostream& operator<<(std::ostream& s, const HGCRecHit& hit);
Expand Down
21 changes: 4 additions & 17 deletions DataFormats/HGCRecHit/src/HGCRecHit.cc
Expand Up @@ -70,25 +70,12 @@ void HGCRecHit::setSignalOverSigmaNoise(float sOverNoise) {

float HGCRecHit::signalOverSigmaNoise() const { return (float)signalOverSigmaNoise_ * 0.125f; }

void HGCRecHit::setTimeError(uint8_t timeErrBits) {
// take the bits and put them in the right spot
setAux((~0xFF & aux()) | timeErrBits);
void HGCRecHit::setTimeError(float timeErr) {
//expected resolution on single cell for that given S/N
timeError_ = timeErr;
}

float HGCRecHit::timeError() const {
uint32_t timeErrorBits = 0xFF & aux();
// all bits off --> time reco bailed out (return negative value)
if ((0xFF & timeErrorBits) == 0x00)
return -1;
// all bits on --> time error over 5 ns (return large value)
if ((0xFF & timeErrorBits) == 0xFF)
return 10000;

float LSB = 1.26008;
uint8_t exponent = timeErrorBits >> 5;
uint8_t significand = timeErrorBits & ~(0x7 << 5);
return pow(2., exponent) * significand * LSB / 1000.f;
}
float HGCRecHit::timeError() const { return timeError_; }

bool HGCRecHit::isTimeValid() const {
if (timeError() <= 0)
Expand Down
3 changes: 2 additions & 1 deletion DataFormats/HGCRecHit/src/classes_def.xml
Expand Up @@ -9,7 +9,8 @@
<class name="edm::RefVector<edm::SortedCollection<HGCUncalibratedRecHit,edm::StrictWeakOrdering<HGCUncalibratedRecHit> >,HGCUncalibratedRecHit,edm::refhelper::FindUsingAdvance<edm::SortedCollection<HGCUncalibratedRecHit,edm::StrictWeakOrdering<HGCUncalibratedRecHit> >,HGCUncalibratedRecHit> >"/>
<class name="edm::RefProd<edm::SortedCollection<HGCUncalibratedRecHit,edm::StrictWeakOrdering<HGCUncalibratedRecHit> > >"/>

<class name="HGCRecHit" ClassVersion="13">
<class name="HGCRecHit" ClassVersion="14">
<version ClassVersion="14" checksum="1119576013"/>
<version ClassVersion="13" checksum="311955435"/>
<version ClassVersion="12" checksum="3568613480"/>
<version ClassVersion="11" checksum="3873633606"/>
Expand Down
4 changes: 4 additions & 0 deletions DataFormats/HGCalReco/interface/Trackster.h
Expand Up @@ -29,6 +29,10 @@ namespace ticl {
edm::ProductID seedID;
int seedIndex;

// -99, -1 if not available. ns units otherwise
float time;
float timeError;

// regressed energy
float regressed_energy;

Expand Down
3 changes: 2 additions & 1 deletion DataFormats/HGCalReco/src/classes_def.xml
@@ -1,6 +1,7 @@

<lcgdict>
<class name="ticl::Trackster" ClassVersion="5">
<class name="ticl::Trackster" ClassVersion="6">
<version ClassVersion="6" checksum="3003722924"/>
<version ClassVersion="5" checksum="1679804166"/>
<version ClassVersion="4" checksum="2963413313"/>
<version ClassVersion="3" checksum="3878166595"/>
Expand Down
1 change: 1 addition & 0 deletions RecoHGCal/TICL/BuildFile.xml
Expand Up @@ -2,6 +2,7 @@
<use name="DataFormats/HGCalReco"/>
<use name="DataFormats/VertexReco"/>
<use name="FWCore/PluginManager"/>
<use name="RecoLocalCalo/HGCalRecProducers"/>
<export>
<lib name="1"/>
</export>
1 change: 0 additions & 1 deletion RecoHGCal/TICL/plugins/BuildFile.xml
Expand Up @@ -19,7 +19,6 @@
<use name="Geometry/Records"/>
<use name="PhysicsTools/TensorFlow"/>
<use name="RecoHGCal/TICL"/>

<library file="*.cc" name="RecoHGCalTICLPlugins">
<flags EDM_PLUGIN="1"/>
</library>
15 changes: 9 additions & 6 deletions RecoHGCal/TICL/plugins/HGCGraph.cc
Expand Up @@ -13,7 +13,7 @@ void HGCGraph::makeAndConnectDoublets(const TICLLayerTiles &histo,
int nPhiBins,
const std::vector<reco::CaloCluster> &layerClusters,
const std::vector<float> &mask,
const edm::ValueMap<float> &layerClustersTime,
const edm::ValueMap<std::pair<float, float>> &layerClustersTime,
int deltaIEta,
int deltaIPhi,
float minCosTheta,
Expand Down Expand Up @@ -125,12 +125,15 @@ void HGCGraph::makeAndConnectDoublets(const TICLLayerTiles &histo,

bool HGCGraph::areTimeCompatible(int innerIdx,
int outerIdx,
const edm::ValueMap<float> &layerClustersTime,
const edm::ValueMap<std::pair<float, float>> &layerClustersTime,
float maxDeltaTime) {
float timeIn = layerClustersTime.get(innerIdx);
float timeOut = layerClustersTime.get(outerIdx);
float timeIn = layerClustersTime.get(innerIdx).first;
float timeInE = layerClustersTime.get(innerIdx).second;
float timeOut = layerClustersTime.get(outerIdx).first;
float timeOutE = layerClustersTime.get(outerIdx).second;

return (timeIn == -99 || timeOut == -99 || std::abs(timeIn - timeOut) < maxDeltaTime);
return (timeIn == -99 || timeOut == -99 ||
std::abs(timeIn - timeOut) < maxDeltaTime * sqrt(timeInE * timeInE + timeOutE * timeOutE));
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would it make sense to also check that the Time Error is meaningful, i.e. it is >0?
I'm not sure if we can have a reasonable time but not a meaningful error.

Copy link
Contributor Author

@amartelli amartelli Jan 17, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No, should not be the case:

  • the error for recHits is sampled from a positive function [between 0.02 and 5.5 ns] for hits with time (and set to -1 for recHits in the scintillator and those without time)
  • only hits with time error >= 0 are contributing to layer clusters
  • all time and timeError estimates in cascade are from a weighted mean, so starting with positive input errors only a positive error on the mean is allowed
  • I don't see any way the timeError can be < 0, except when time = -99 (-99, -1 are the default values returned for no time). But this is already taken into account

If apart from this timeError is <= 0, it is for a bug somewhere else in the chain.

Can I suggest to leave this as it is? (to help spotting the bug in case?)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sure, thanks for checking!

}

//also return a vector of seedIndex for the reconstructed tracksters
Expand All @@ -141,7 +144,7 @@ void HGCGraph::findNtuplets(std::vector<HGCDoublet::HGCntuplet> &foundNtuplets,
unsigned int maxOutInHops) {
HGCDoublet::HGCntuplet tmpNtuplet;
tmpNtuplet.reserve(minClustersPerNtuplet);
std::vector<std::pair<unsigned int, unsigned int> > outInToVisit;
std::vector<std::pair<unsigned int, unsigned int>> outInToVisit;
for (auto rootDoublet : theRootDoublets_) {
tmpNtuplet.clear();
outInToVisit.clear();
Expand Down
7 changes: 5 additions & 2 deletions RecoHGCal/TICL/plugins/HGCGraph.h
Expand Up @@ -19,7 +19,7 @@ class HGCGraph {
int nPhiBins,
const std::vector<reco::CaloCluster> &layerClusters,
const std::vector<float> &mask,
const edm::ValueMap<float> &layerClustersTime,
const edm::ValueMap<std::pair<float, float>> &layerClustersTime,
int deltaIEta,
int deltaIPhi,
float minCosThetai,
Expand All @@ -28,7 +28,10 @@ class HGCGraph {
int maxNumberOfLayers,
float maxDeltaTime);

bool areTimeCompatible(int innerIdx, int outerIdx, const edm::ValueMap<float> &layerClustersTime, float maxDeltaTime);
bool areTimeCompatible(int innerIdx,
int outerIdx,
const edm::ValueMap<std::pair<float, float>> &layerClustersTime,
float maxDeltaTime);

std::vector<HGCDoublet> &getAllDoublets() { return allDoublets_; }
void findNtuplets(std::vector<HGCDoublet::HGCntuplet> &foundNtuplets,
Expand Down
Expand Up @@ -90,6 +90,7 @@ void MultiClustersFromTrackstersProducer::produce(edm::Event& evt, const edm::Ev
temp.setCorrectedEnergy(total_weight);
temp.setPosition(math::XYZPoint(baricenter[0], baricenter[1], baricenter[2]));
temp.setAlgoId(reco::CaloCluster::hgcal_em);
temp.setTime(trackster.time, trackster.timeError);
multiclusters->push_back(temp);
});

Expand Down
4 changes: 2 additions & 2 deletions RecoHGCal/TICL/plugins/PatternRecognitionAlgoBase.h
Expand Up @@ -32,15 +32,15 @@ namespace ticl {
const edm::EventSetup& es;
const std::vector<reco::CaloCluster>& layerClusters;
const std::vector<float>& mask;
const edm::ValueMap<float>& layerClustersTime;
const edm::ValueMap<std::pair<float, float>>& layerClustersTime;
const TICLLayerTiles& tiles;
const std::vector<TICLSeedingRegion>& regions;

Inputs(const edm::Event& eV,
const edm::EventSetup& eS,
const std::vector<reco::CaloCluster>& lC,
const std::vector<float>& mS,
const edm::ValueMap<float>& lT,
const edm::ValueMap<std::pair<float, float>>& lT,
const TICLLayerTiles& tL,
const std::vector<TICLSeedingRegion>& rG)
: ev(eV), es(eS), layerClusters(lC), mask(mS), layerClustersTime(lT), tiles(tL), regions(rG) {}
Expand Down
36 changes: 34 additions & 2 deletions RecoHGCal/TICL/plugins/PatternRecognitionbyCA.cc
Expand Up @@ -8,6 +8,7 @@
#include "FWCore/Utilities/interface/Exception.h"
#include "PatternRecognitionbyCA.h"
#include "HGCGraph.h"
#include "RecoLocalCalo/HGCalRecProducers/interface/ComputeClusterTime.h"

using namespace ticl;

Expand Down Expand Up @@ -70,13 +71,36 @@ void PatternRecognitionbyCA::makeTracksters(const PatternRecognitionAlgoBase::In
//#ifdef FP_DEBUG
const auto &doublets = theGraph_->getAllDoublets();
int tracksterId = 0;

for (auto const &ntuplet : foundNtuplets) {
std::set<unsigned int> effective_cluster_idx;
std::pair<std::set<unsigned int>::iterator, bool> retVal;

std::vector<float> times;
std::vector<float> timeErrors;

for (auto const &doublet : ntuplet) {
auto innerCluster = doublets[doublet].innerClusterId();
auto outerCluster = doublets[doublet].outerClusterId();
effective_cluster_idx.insert(innerCluster);
effective_cluster_idx.insert(outerCluster);

retVal = effective_cluster_idx.insert(innerCluster);
if (retVal.second) {
float time = input.layerClustersTime.get(innerCluster).first;
if (time > -99) {
times.push_back(time);
timeErrors.push_back(1. / pow(input.layerClustersTime.get(innerCluster).second, 2));
}
}

retVal = effective_cluster_idx.insert(outerCluster);
if (retVal.second) {
float time = input.layerClustersTime.get(outerCluster).first;
if (time > -99) {
times.push_back(time);
timeErrors.push_back(1. / pow(input.layerClustersTime.get(outerCluster).second, 2));
}
}

if (algo_verbosity_ > Advanced) {
LogDebug("HGCPatterRecoByCA") << "New doublet " << doublet << " for trackster: " << result.size() << " InnerCl "
<< innerCluster << " " << input.layerClusters[innerCluster].x() << " "
Expand All @@ -99,6 +123,14 @@ void PatternRecognitionbyCA::makeTracksters(const PatternRecognitionAlgoBase::In
//if a seeding region does not lead to any trackster
tmp.seedID = input.regions[0].collectionID;
tmp.seedIndex = seedIndices[tracksterId];

std::pair<float, float> timeTrackster(-99., -1.);
if (times.size() >= 3) {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Any specific reason why you need at least 3 measurements? Maybe it would make sense to "hide" this into the timeEstimator code? And make the number 3 a const there?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

3 was tuned looking at the full shower performance using recHits (TDR time).

It is kept and propagated in the chain, as reasonable to maximize efficiency and still allow some clustering of the information (truncation and weighted mean)
I don't exclude that this could be tuned in the future and it might also be different for layerClusters and tracksters.

What if I replace this by a configurable parameter, leaving it in the calling code?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ciao @amartelli maybe we could do both? If you think the parameter could be different for different use cases, then it makes sense to leave it to the caller and make it parametrizable. Yet, maybe you could include it as a parameter into the fixSizeHighestDensity method so that the caller does not need to add additional checks on his side?
You could default it if you do not feel like chasing all places in CMSSW where this method is already called (if there are few places, it's better to add it explicitly).

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree, it's reasonable, I'll change it. thanks

hgcalsimclustertime::ComputeClusterTime timeEstimator;
timeTrackster = timeEstimator.fixSizeHighestDensity(times, timeErrors);
}
tmp.time = timeTrackster.first;
tmp.timeError = timeTrackster.second;
std::copy(std::begin(effective_cluster_idx), std::end(effective_cluster_idx), std::back_inserter(tmp.vertices));
result.push_back(tmp);
tracksterId++;
Expand Down
9 changes: 5 additions & 4 deletions RecoHGCal/TICL/plugins/TrackstersProducer.cc
Expand Up @@ -45,7 +45,7 @@ class TrackstersProducer : public edm::stream::EDProducer<edm::GlobalCache<Track
edm::EDGetTokenT<std::vector<reco::CaloCluster>> clusters_token_;
edm::EDGetTokenT<std::vector<float>> filtered_layerclusters_mask_token_;
edm::EDGetTokenT<std::vector<float>> original_layerclusters_mask_token_;
edm::EDGetTokenT<edm::ValueMap<float>> clustersTime_token_;
edm::EDGetTokenT<edm::ValueMap<std::pair<float, float>>> clustersTime_token_;
edm::EDGetTokenT<TICLLayerTiles> layer_clusters_tiles_token_;
edm::EDGetTokenT<std::vector<TICLSeedingRegion>> seeding_regions_token_;

Expand Down Expand Up @@ -77,7 +77,8 @@ TrackstersProducer::TrackstersProducer(const edm::ParameterSet& ps, const Tracks
clusters_token_ = consumes<std::vector<reco::CaloCluster>>(ps.getParameter<edm::InputTag>("layer_clusters"));
filtered_layerclusters_mask_token_ = consumes<std::vector<float>>(ps.getParameter<edm::InputTag>("filtered_mask"));
original_layerclusters_mask_token_ = consumes<std::vector<float>>(ps.getParameter<edm::InputTag>("original_mask"));
clustersTime_token_ = consumes<edm::ValueMap<float>>(ps.getParameter<edm::InputTag>("time_layerclusters"));
clustersTime_token_ =
consumes<edm::ValueMap<std::pair<float, float>>>(ps.getParameter<edm::InputTag>("time_layerclusters"));
layer_clusters_tiles_token_ = consumes<TICLLayerTiles>(ps.getParameter<edm::InputTag>("layer_clusters_tiles"));
seeding_regions_token_ = consumes<std::vector<TICLSeedingRegion>>(ps.getParameter<edm::InputTag>("seeding_regions"));

Expand All @@ -99,7 +100,7 @@ void TrackstersProducer::fillDescriptions(edm::ConfigurationDescriptions& descri
desc.add<double>("min_cos_pointing", -1.);
desc.add<int>("missing_layers", 0);
desc.add<int>("min_clusters_per_ntuplet", 10);
desc.add<double>("max_delta_time", 0.09);
desc.add<double>("max_delta_time", 3.); //nsigma
desc.add<bool>("out_in_dfs", true);
desc.add<int>("max_out_in_hops", 10);
desc.add<std::string>("eid_graph_path", "RecoHGCal/TICL/data/tf_models/energy_id_v0.pb");
Expand All @@ -119,7 +120,7 @@ void TrackstersProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
edm::Handle<std::vector<reco::CaloCluster>> cluster_h;
edm::Handle<std::vector<float>> filtered_layerclusters_mask_h;
edm::Handle<std::vector<float>> original_layerclusters_mask_h;
edm::Handle<edm::ValueMap<float>> time_clusters_h;
edm::Handle<edm::ValueMap<std::pair<float, float>>> time_clusters_h;
edm::Handle<TICLLayerTiles> layer_clusters_tiles_h;
edm::Handle<std::vector<TICLSeedingRegion>> seeding_regions_h;

Expand Down
72 changes: 72 additions & 0 deletions RecoLocalCalo/HGCalRecProducers/interface/ComputeClusterTime.h
@@ -0,0 +1,72 @@
#ifndef RecoLocalCalo_HGCalRecProducers_ComputeClusterTime_h
#define RecoLocalCalo_HGCalRecProducers_ComputeClusterTime_h

// user include files
#include <algorithm>
#include <cmath>
#include <numeric>
#include <vector>
#include "TF1.h"
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

remove unneeded include?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes, thanks


// functions to select the hits to compute the time of a given cluster
// start with the only hits with timing information
// average among the hits contained in the chosen time interval
// weighted average wrt resolution or preferred function

// N.B. time is corrected wrt vtx-calorimeter distance
// with straight line and light speed hypothesis
// for charged tracks or heavy particles (longer track length or beta < 1)
// need to correct the offset at analysis level

namespace hgcalsimclustertime {

/*
std::vector<size_t> decrease_sorted_indices(const std::vector<float>& v){
// initialize original index locations
std::vector<size_t> idx(v.size());
for (size_t i = 0; i != idx.size(); ++i) idx[i] = i;
// sort indices based on comparing values in v (decreasing order)
std::sort(idx.begin(), idx.end(),
[&v](size_t i1, size_t i2) {return v[i1] < v[i2];} );
return idx;
};
*/

class ComputeClusterTime {
public:
ComputeClusterTime(float Xmix, float Xmax, float Cterm, float Aterm);
ComputeClusterTime();

void setParameters(float Xmix, float Xmax, float Cterm, float Aterm);

//time resolution parametrization
float timeResolution(float xVal);

float getTimeError(std::string type, float xVal);

//time-interval based on that ~210ps wide and with the highest number of hits
//apply weights if provided => weighted mean
//return also error on the mean
std::pair<float, float> fixSizeHighestDensity(std::vector<float>& time,
std::vector<float> weight = std::vector<float>(),
float deltaT = 0.210, /*time window in ns*/
float timeWidthBy = 0.5);

//same as before but provide weights
//configure the reweight through type
/* std::pair<float,float> fixSizeHighestDensityEnergyResWeig(std::vector<float>& t, */
/* std::vector<float>& w, */
/* std::string& type, */
/* float deltaT = 0.210, //time window in ns */
/* float timeWidthBy = 0.5) {}; */

private:
float _Xmin;
float _Xmax;
float _Cterm;
float _Aterm;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

leading underscores are frowned upon; please use a trailing (or prepend something, like "m") to denote the member data

};

} // namespace hgcalsimclustertime

#endif