CMS 3D CMS Logo

TOFPIDProducer.cc
Go to the documentation of this file.
8 
15 #include "CLHEP/Units/GlobalPhysicalConstants.h"
16 #include "CLHEP/Units/GlobalSystemOfUnits.h"
17 
18 using namespace std;
19 using namespace edm;
20 
22 public:
24 
25  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
26 
27  template <class H, class T>
28  void fillValueMap(edm::Event& iEvent,
29  const edm::Handle<H>& handle,
30  const std::vector<T>& vec,
31  const std::string& name) const;
32 
33  void produce(edm::Event& ev, const edm::EventSetup& es) final;
34 
35 private:
36  static constexpr char t0Name[] = "t0";
37  static constexpr char sigmat0Name[] = "sigmat0";
38  static constexpr char t0safeName[] = "t0safe";
39  static constexpr char sigmat0safeName[] = "sigmat0safe";
40  static constexpr char probPiName[] = "probPi";
41  static constexpr char probKName[] = "probK";
42  static constexpr char probPName[] = "probP";
43 
52  double vtxMaxSigmaT_;
53  double maxDz_;
55  double minProbHeavy_;
56  double fixedT0Error_;
57 };
58 
60  : tracksToken_(consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("tracksSrc"))),
61  t0Token_(consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("t0Src"))),
62  tmtdToken_(consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tmtdSrc"))),
63  sigmat0Token_(consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("sigmat0Src"))),
64  sigmatmtdToken_(consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("sigmatmtdSrc"))),
65  pathLengthToken_(consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("pathLengthSrc"))),
66  pToken_(consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("pSrc"))),
67  vtxsToken_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vtxsSrc"))),
68  vtxMaxSigmaT_(iConfig.getParameter<double>("vtxMaxSigmaT")),
69  maxDz_(iConfig.getParameter<double>("maxDz")),
70  maxDtSignificance_(iConfig.getParameter<double>("maxDtSignificance")),
71  minProbHeavy_(iConfig.getParameter<double>("minProbHeavy")),
72  fixedT0Error_(iConfig.getParameter<double>("fixedT0Error")) {
73  produces<edm::ValueMap<float>>(t0Name);
74  produces<edm::ValueMap<float>>(sigmat0Name);
75  produces<edm::ValueMap<float>>(t0safeName);
76  produces<edm::ValueMap<float>>(sigmat0safeName);
77  produces<edm::ValueMap<float>>(probPiName);
78  produces<edm::ValueMap<float>>(probKName);
79  produces<edm::ValueMap<float>>(probPName);
80 }
81 
82 // Configuration descriptions
85  desc.add<edm::InputTag>("tracksSrc", edm::InputTag("generalTracks"))->setComment("Input tracks collection");
86  desc.add<edm::InputTag>("t0Src", edm::InputTag("trackExtenderWithMTD:generalTrackt0"))
87  ->setComment("Input ValueMap for track time at beamline");
88  desc.add<edm::InputTag>("tmtdSrc", edm::InputTag("trackExtenderWithMTD:generalTracktmtd"))
89  ->setComment("Input ValueMap for track time at MTD");
90  desc.add<edm::InputTag>("sigmat0Src", edm::InputTag("trackExtenderWithMTD:generalTracksigmat0"))
91  ->setComment("Input ValueMap for track time uncertainty at beamline");
92  desc.add<edm::InputTag>("sigmatmtdSrc", edm::InputTag("trackExtenderWithMTD:generalTracksigmatmtd"))
93  ->setComment("Input ValueMap for track time uncertainty at MTD");
94  desc.add<edm::InputTag>("pathLengthSrc", edm::InputTag("trackExtenderWithMTD:generalTrackPathLength"))
95  ->setComment("Input ValueMap for track path lengh from beamline to MTD");
96  desc.add<edm::InputTag>("pSrc", edm::InputTag("trackExtenderWithMTD:generalTrackp"))
97  ->setComment("Input ValueMap for track momentum magnitude (normally from refit with MTD hits)");
98  desc.add<edm::InputTag>("vtxsSrc", edm::InputTag("unsortedOfflinePrimaryVertices4DwithPID"))
99  ->setComment("Input primary vertex collection");
100  desc.add<double>("vtxMaxSigmaT", 0.025)
101  ->setComment("Maximum primary vertex time uncertainty for use in particle id [ns]");
102  desc.add<double>("maxDz", 0.1)
103  ->setComment("Maximum distance in z for track-primary vertex association for particle id [cm]");
104  desc.add<double>("maxDtSignificance", 5.0)
105  ->setComment(
106  "Maximum distance in time (normalized by uncertainty) for track-primary vertex association for particle id");
107  desc.add<double>("minProbHeavy", 0.75)
108  ->setComment("Minimum probability for a particle to be a kaon or proton before reassigning the timestamp");
109  desc.add<double>("fixedT0Error", 0.)->setComment("Use a fixed T0 uncertainty [ns]");
110 
111  descriptions.add("tofPIDProducer", desc);
112 }
113 
114 template <class H, class T>
116  const edm::Handle<H>& handle,
117  const std::vector<T>& vec,
118  const std::string& name) const {
119  auto out = std::make_unique<edm::ValueMap<T>>();
121  filler.insert(handle, vec.begin(), vec.end());
122  filler.fill();
123  iEvent.put(std::move(out), name);
124 }
125 
127  constexpr double m_k = 0.493677; //[GeV]
128  constexpr double m_p = 0.9382720813; //[GeV]
129  constexpr double c_cm_ns = CLHEP::c_light * CLHEP::ns / CLHEP::cm; //[cm/ns]
130  constexpr double c_inv = 1.0 / c_cm_ns;
131 
133  ev.getByToken(tracksToken_, tracksH);
134  const auto& tracks = *tracksH;
135 
137  ev.getByToken(t0Token_, t0H);
138  const auto& t0In = *t0H;
139 
141  ev.getByToken(tmtdToken_, tmtdH);
142  const auto& tmtdIn = *tmtdH;
143 
145  ev.getByToken(sigmat0Token_, sigmat0H);
146  const auto& sigmat0In = *sigmat0H;
147 
149  ev.getByToken(sigmatmtdToken_, sigmatmtdH);
150  const auto& sigmatmtdIn = *sigmatmtdH;
151 
153  ev.getByToken(pathLengthToken_, pathLengthH);
154  const auto& pathLengthIn = *pathLengthH;
155 
157  ev.getByToken(pToken_, pH);
158  const auto& pIn = *pH;
159 
161  ev.getByToken(vtxsToken_, vtxsH);
162  const auto& vtxs = *vtxsH;
163 
164  //output value maps (PID probabilities and recalculated time at beamline)
165  std::vector<float> t0OutRaw;
166  std::vector<float> sigmat0OutRaw;
167  std::vector<float> t0safeOutRaw;
168  std::vector<float> sigmat0safeOutRaw;
169  std::vector<float> probPiOutRaw;
170  std::vector<float> probKOutRaw;
171  std::vector<float> probPOutRaw;
172 
173  //Do work here
174  for (unsigned int itrack = 0; itrack < tracks.size(); ++itrack) {
175  const reco::Track& track = tracks[itrack];
176  const reco::TrackRef trackref(tracksH, itrack);
177  float t0 = t0In[trackref];
178  float t0safe = t0;
179  float sigmat0safe = sigmat0In[trackref];
180  float sigmatmtd = (sigmatmtdIn[trackref] > 0. && fixedT0Error_ > 0.) ? fixedT0Error_ : sigmatmtdIn[trackref];
181  float sigmat0 = sigmatmtd;
182 
183  float prob_pi = -1.;
184  float prob_k = -1.;
185  float prob_p = -1.;
186 
187  if (sigmat0 > 0.) {
188  double rsigmazsq = 1. / track.dzError() / track.dzError();
189  double rsigmat = 1. / sigmatmtd;
190 
191  //find associated vertex
192  int vtxidx = -1;
193  int vtxidxmindz = -1;
194  int vtxidxminchisq = -1;
195  double mindz = maxDz_;
196  double minchisq = std::numeric_limits<double>::max();
197  //first try based on association weights, but keep track of closest in z and z-t as well
198  for (unsigned int ivtx = 0; ivtx < vtxs.size(); ++ivtx) {
199  const reco::Vertex& vtx = vtxs[ivtx];
200  float w = vtx.trackWeight(trackref);
201  if (w > 0.5) {
202  vtxidx = ivtx;
203  break;
204  }
205  double dz = std::abs(track.dz(vtx.position()));
206  if (dz < mindz) {
207  mindz = dz;
208  vtxidxmindz = ivtx;
209  }
210  if (vtx.tError() > 0. && vtx.tError() < vtxMaxSigmaT_) {
211  double dt = std::abs(t0 - vtx.t());
212  double dtsig = dt * rsigmat;
213  double chisq = dz * dz * rsigmazsq + dtsig * dtsig;
214  if (dz < maxDz_ && dtsig < maxDtSignificance_ && chisq < minchisq) {
215  minchisq = chisq;
216  vtxidxminchisq = ivtx;
217  }
218  }
219  }
220 
221  //if no vertex found based on association weights, fall back to closest in z or z-t
222  if (vtxidx < 0) {
223  //if closest vertex in z does not have valid time information, just use it,
224  //otherwise use the closest vertex in z-t plane with timing info, with a fallback to the closest in z
225  if (vtxidxmindz >= 0 && !(vtxs[vtxidxmindz].tError() > 0. && vtxs[vtxidxmindz].tError() < vtxMaxSigmaT_)) {
226  vtxidx = vtxidxmindz;
227  } else if (vtxidxminchisq >= 0) {
228  vtxidx = vtxidxminchisq;
229  } else if (vtxidxmindz >= 0) {
230  vtxidx = vtxidxmindz;
231  }
232  }
233 
234  //testing mass hypotheses only possible if there is an associated vertex with time information
235  if (vtxidx >= 0 && vtxs[vtxidx].tError() > 0. && vtxs[vtxidx].tError() < vtxMaxSigmaT_) {
236  //compute chisq in z-t plane for nominal vertex and mass hypothesis (pion)
237  const reco::Vertex& vtxnom = vtxs[vtxidx];
238  double dznom = std::abs(track.dz(vtxnom.position()));
239  double dtnom = std::abs(t0 - vtxnom.t());
240  double dtsignom = dtnom * rsigmat;
241  double chisqnom = dznom * dznom * rsigmazsq + dtsignom * dtsignom;
242 
243  //recompute t0 for alternate mass hypotheses
244  double t0_best = t0;
245 
246  //reliable match, revert to raw mtd time uncertainty
247  if (dtsignom < maxDtSignificance_) {
248  sigmat0safe = sigmatmtd;
249  }
250 
251  double tmtd = tmtdIn[trackref];
252  double pathlength = pathLengthIn[trackref];
253  double magp = pIn[trackref];
254 
255  double gammasq_k = 1. + magp * magp / m_k / m_k;
256  double beta_k = std::sqrt(1. - 1. / gammasq_k);
257  double t0_k = tmtd - pathlength / beta_k * c_inv;
258 
259  double gammasq_p = 1. + magp * magp / m_p / m_p;
260  double beta_p = std::sqrt(1. - 1. / gammasq_p);
261  double t0_p = tmtd - pathlength / beta_p * c_inv;
262 
263  double chisqmin = chisqnom;
264 
265  double chisqmin_pi = chisqnom;
266  double chisqmin_k = std::numeric_limits<double>::max();
267  double chisqmin_p = std::numeric_limits<double>::max();
268  //loop through vertices and check for better matches
269  for (const reco::Vertex& vtx : vtxs) {
270  if (!(vtx.tError() > 0. && vtx.tError() < vtxMaxSigmaT_)) {
271  continue;
272  }
273 
274  double dz = std::abs(track.dz(vtx.position()));
275  if (dz >= maxDz_) {
276  continue;
277  }
278 
279  double chisqdz = dz * dz * rsigmazsq;
280 
281  double dt_k = std::abs(t0_k - vtx.t());
282  double dtsig_k = dt_k * rsigmat;
283  double chisq_k = chisqdz + dtsig_k * dtsig_k;
284 
285  if (dtsig_k < maxDtSignificance_ && chisq_k < chisqmin_k) {
286  chisqmin_k = chisq_k;
287  }
288 
289  double dt_p = std::abs(t0_p - vtx.t());
290  double dtsig_p = dt_p * rsigmat;
291  double chisq_p = chisqdz + dtsig_p * dtsig_p;
292 
293  if (dtsig_p < maxDtSignificance_ && chisq_p < chisqmin_p) {
294  chisqmin_p = chisq_p;
295  }
296 
297  if (dtsig_k < maxDtSignificance_ && chisq_k < chisqmin) {
298  chisqmin = chisq_k;
299  t0_best = t0_k;
300  t0safe = t0_k;
301  sigmat0safe = sigmatmtd;
302  }
303  if (dtsig_p < maxDtSignificance_ && chisq_p < chisqmin) {
304  chisqmin = chisq_p;
305  t0_best = t0_p;
306  t0safe = t0_p;
307  sigmat0safe = sigmatmtd;
308  }
309  }
310 
311  //compute PID probabilities
312  //*TODO* deal with heavier nucleons and/or BSM case here?
313  double rawprob_pi = exp(-0.5 * chisqmin_pi);
314  double rawprob_k = exp(-0.5 * chisqmin_k);
315  double rawprob_p = exp(-0.5 * chisqmin_p);
316 
317  double normprob = 1. / (rawprob_pi + rawprob_k + rawprob_p);
318 
319  prob_pi = rawprob_pi * normprob;
320  prob_k = rawprob_k * normprob;
321  prob_p = rawprob_p * normprob;
322 
323  double prob_heavy = 1. - prob_pi;
324 
325  if (prob_heavy > minProbHeavy_) {
326  t0 = t0_best;
327  }
328  }
329  }
330 
331  t0OutRaw.push_back(t0);
332  sigmat0OutRaw.push_back(sigmat0);
333  t0safeOutRaw.push_back(t0safe);
334  sigmat0safeOutRaw.push_back(sigmat0safe);
335  probPiOutRaw.push_back(prob_pi);
336  probKOutRaw.push_back(prob_k);
337  probPOutRaw.push_back(prob_p);
338  }
339 
340  fillValueMap(ev, tracksH, t0OutRaw, t0Name);
341  fillValueMap(ev, tracksH, sigmat0OutRaw, sigmat0Name);
342  fillValueMap(ev, tracksH, t0safeOutRaw, t0safeName);
343  fillValueMap(ev, tracksH, sigmat0safeOutRaw, sigmat0safeName);
344  fillValueMap(ev, tracksH, probPiOutRaw, probPiName);
345  fillValueMap(ev, tracksH, probKOutRaw, probKName);
346  fillValueMap(ev, tracksH, probPOutRaw, probPName);
347 }
348 
349 //define this as a plug-in
ConfigurationDescriptions.h
TOFPIDProducer::pToken_
edm::EDGetTokenT< edm::ValueMap< float > > pToken_
Definition: TOFPIDProducer.cc:50
TOFPIDProducer::probPName
static constexpr char probPName[]
Definition: TOFPIDProducer.cc:42
TOFPIDProducer::vtxMaxSigmaT_
double vtxMaxSigmaT_
Definition: TOFPIDProducer.cc:52
PDWG_EXOHSCP_cff.tracks
tracks
Definition: PDWG_EXOHSCP_cff.py:28
edm::ParameterSetDescription::add
ParameterDescriptionBase * add(U const &iLabel, T const &value)
Definition: ParameterSetDescription.h:95
dqmMemoryStats.float
float
Definition: dqmMemoryStats.py:127
ESHandle.h
TOFPIDProducer::fixedT0Error_
double fixedT0Error_
Definition: TOFPIDProducer.cc:56
patZpeak.handle
handle
Definition: patZpeak.py:23
edm::EDGetTokenT< reco::TrackCollection >
TOFPIDProducer::vtxsToken_
edm::EDGetTokenT< reco::VertexCollection > vtxsToken_
Definition: TOFPIDProducer.cc:51
edm
HLT enums.
Definition: AlignableModifier.h:19
TOFPIDProducer
Definition: TOFPIDProducer.cc:21
c_cm_ns
constexpr double c_cm_ns
Definition: Phase2TrackerDigitizerAlgorithm.h:55
reco::VertexCollection
std::vector< Vertex > VertexCollection
collection of Vertex objects
Definition: VertexFwd.h:9
edm::ParameterSetDescription
Definition: ParameterSetDescription.h:52
EDProducer.h
reco::Vertex::position
const Point & position() const
position
Definition: Vertex.h:114
TOFPIDProducer::sigmat0safeName
static constexpr char sigmat0safeName[]
Definition: TOFPIDProducer.cc:39
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:45
edm::Handle
Definition: AssociativeIterator.h:50
edm::Ref< TrackCollection >
dt
float dt
Definition: AMPTWrapper.h:136
MakerMacros.h
Track.h
TOFPIDProducer::minProbHeavy_
double minProbHeavy_
Definition: TOFPIDProducer.cc:55
TrackFwd.h
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
BeamSpot.h
edm::ConfigurationDescriptions::add
void add(std::string const &label, ParameterSetDescription const &psetDescription)
Definition: ConfigurationDescriptions.cc:57
TOFPIDProducer::tracksToken_
edm::EDGetTokenT< reco::TrackCollection > tracksToken_
Definition: TOFPIDProducer.cc:44
FrontierCondition_GT_autoExpress_cfi.t0
t0
Definition: FrontierCondition_GT_autoExpress_cfi.py:148
w
const double w
Definition: UKUtility.cc:23
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
fillDescriptions
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
reco::Track
Definition: Track.h:27
TOFPIDProducer::maxDtSignificance_
double maxDtSignificance_
Definition: TOFPIDProducer.cc:54
edm::ConfigurationDescriptions
Definition: ConfigurationDescriptions.h:28
badGlobalMuonTaggersAOD_cff.vtx
vtx
Definition: badGlobalMuonTaggersAOD_cff.py:5
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
Vertex.h
TOFPIDProducer::t0Token_
edm::EDGetTokenT< edm::ValueMap< float > > t0Token_
Definition: TOFPIDProducer.cc:45
TOFPIDProducer::t0Name
static constexpr char t0Name[]
Definition: TOFPIDProducer.cc:36
TOFPIDProducer::probKName
static constexpr char probKName[]
Definition: TOFPIDProducer.cc:41
TOFPIDProducer::tmtdToken_
edm::EDGetTokenT< edm::ValueMap< float > > tmtdToken_
Definition: TOFPIDProducer.cc:46
TOFPIDProducer::t0safeName
static constexpr char t0safeName[]
Definition: TOFPIDProducer.cc:38
HLT_2018_cff.InputTag
InputTag
Definition: HLT_2018_cff.py:79016
edm::ParameterSet
Definition: ParameterSet.h:36
Event.h
SiStripPI::max
Definition: SiStripPayloadInspectorHelper.h:169
trigObjTnPSource_cfi.filler
filler
Definition: trigObjTnPSource_cfi.py:21
iEvent
int iEvent
Definition: GenABIO.cc:224
TOFPIDProducer::fillDescriptions
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
Definition: TOFPIDProducer.cc:83
TOFPIDProducer::pathLengthToken_
edm::EDGetTokenT< edm::ValueMap< float > > pathLengthToken_
Definition: TOFPIDProducer.cc:49
TOFPIDProducer::sigmatmtdToken_
edm::EDGetTokenT< edm::ValueMap< float > > sigmatmtdToken_
Definition: TOFPIDProducer.cc:48
edm::stream::EDProducer
Definition: EDProducer.h:38
edm::EventSetup
Definition: EventSetup.h:57
ValueMap.h
VertexFwd.h
TOFPIDProducer::maxDz_
double maxDz_
Definition: TOFPIDProducer.cc:53
TOFPIDProducer::TOFPIDProducer
TOFPIDProducer(const ParameterSet &pset)
Definition: TOFPIDProducer.cc:59
TOFPIDProducer::sigmat0Token_
edm::EDGetTokenT< edm::ValueMap< float > > sigmat0Token_
Definition: TOFPIDProducer.cc:47
TOFPIDProducer::probPiName
static constexpr char probPiName[]
Definition: TOFPIDProducer.cc:40
eostools.move
def move(src, dest)
Definition: eostools.py:511
std
Definition: JetResolutionObject.h:76
c_inv
constexpr double c_inv
Definition: Phase2TrackerDigitizerAlgorithm.h:56
PVValHelper::dz
Definition: PVValidationHelpers.h:50
Frameworkfwd.h
TOFPIDProducer::fillValueMap
void fillValueMap(edm::Event &iEvent, const edm::Handle< H > &handle, const std::vector< T > &vec, const std::string &name) const
Definition: TOFPIDProducer.cc:115
reco::Vertex::t
double t() const
t coordinate
Definition: Vertex.h:122
edm::ValueMap
Definition: ValueMap.h:107
ev
bool ev
Definition: Hydjet2Hadronizer.cc:95
Skims_PA_cff.name
name
Definition: Skims_PA_cff.py:17
EventSetup.h
TOFPIDProducer::produce
void produce(edm::Event &ev, const edm::EventSetup &es) final
Definition: TOFPIDProducer.cc:126
MillePedeFileConverter_cfg.out
out
Definition: MillePedeFileConverter_cfg.py:31
HLT_2018_cff.track
track
Definition: HLT_2018_cff.py:10352
edm::helper::Filler
Definition: ValueMap.h:22
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
ParameterSet.h
edm::ParameterDescriptionNode::setComment
void setComment(std::string const &value)
Definition: ParameterDescriptionNode.cc:106
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
edm::Event
Definition: Event.h:73
edm::InputTag
Definition: InputTag.h:15
reco::TrackCollection
std::vector< Track > TrackCollection
collection of Tracks
Definition: TrackFwd.h:14
reco::Vertex
Definition: Vertex.h:35
muonDTDigis_cfi.pset
pset
Definition: muonDTDigis_cfi.py:27
TOFPIDProducer::sigmat0Name
static constexpr char sigmat0Name[]
Definition: TOFPIDProducer.cc:37