CMS 3D CMS Logo
PFCandConnector.cc
Go to the documentation of this file.
1 #include "RecoParticleFlow/PFProducer/interface/PFCandConnector.h"
2 #include "DataFormats/ParticleFlowReco/interface/PFDisplacedVertex.h"
3 #include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
4 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
5 #include "FWCore/MessageLogger/interface/MessageLogger.h"
6 
7 using namespace reco;
8 using namespace std;
9 
10 const double PFCandConnector::pion_mass2 = 0.0194;
11 
12 const reco::PFCandidate::Flags PFCandConnector::fT_TO_DISP_ = PFCandidate::T_TO_DISP;
13 const reco::PFCandidate::Flags PFCandConnector::fT_FROM_DISP_ = PFCandidate::T_FROM_DISP;
14 
15 void PFCandConnector::setParameters(bool bCorrect,
16  bool bCalibPrimary,
17  double dptRel_PrimaryTrack,
18  double dptRel_MergedTrack,
19  double ptErrorSecondary,
20  const std::vector<double>& nuclCalibFactors) {
21  bCorrect_ = bCorrect;
22  bCalibPrimary_ = bCalibPrimary;
23  dptRel_PrimaryTrack_ = dptRel_PrimaryTrack;
24  dptRel_MergedTrack_ = dptRel_MergedTrack;
25  ptErrorSecondary_ = ptErrorSecondary;
26 
27  if (nuclCalibFactors.size() == 5) {
28  fConst_[0] = nuclCalibFactors[0];
29  fConst_[1] = nuclCalibFactors[1];
30 
31  fNorm_[0] = nuclCalibFactors[2];
32  fNorm_[1] = nuclCalibFactors[3];
33 
34  fExp_[0] = nuclCalibFactors[4];
35  } else {
36  edm::LogWarning("PFCandConnector")
37  << "Wrong calibration factors for nuclear interactions. The calibration procedure would not be applyed."
38  << std::endl;
39  bCalibPrimary_ = false;
40  }
41 
42  std::string sCorrect = bCorrect_ ? "On" : "Off";
43  edm::LogInfo("PFCandConnector") << " ====================== The PFCandConnector is switched " << sCorrect.c_str()
44  << " ==================== " << std::endl;
45  std::string sCalibPrimary = bCalibPrimary_ ? "used for calibration" : "not used for calibration";
46  if (bCorrect_)
47  edm::LogInfo("PFCandConnector") << "Primary Tracks are " << sCalibPrimary.c_str() << std::endl;
48  if (bCorrect_ && bCalibPrimary_)
49  edm::LogInfo("PFCandConnector") << "Under the condition that the precision on the Primary track is better than "
50  << dptRel_PrimaryTrack_ << " % " << std::endl;
51  if (bCorrect_ && bCalibPrimary_)
52  edm::LogInfo("PFCandConnector") << " and on merged tracks better than " << dptRel_MergedTrack_ << " % "
53  << std::endl;
54  if (bCorrect_ && bCalibPrimary_)
55  edm::LogInfo("PFCandConnector") << " and secondary tracks in some cases more precise than "
56  << ptErrorSecondary_ << " GeV" << std::endl;
57  if (bCorrect_ && bCalibPrimary_)
58  edm::LogInfo("PFCandConnector") << "factor = (" << fConst_[0] << " + " << fConst_[1] << "*cFrac) - (" << fNorm_[0]
59  << " - " << fNorm_[1] << "cFrac)*exp( " << -1 * fExp_[0] << "*pT )" << std::endl;
60  edm::LogInfo("PFCandConnector") << " =========================================================== " << std::endl;
61 }
62 
63 reco::PFCandidateCollection PFCandConnector::connect(PFCandidateCollection& pfCand) const {
65  PFCandidateCollection pfC{};
67  std::vector<bool> bMask;
68  bMask.resize(pfCand.size(), false);
69 
70  // loop on primary
71  if (bCorrect_) {
72  LogTrace("PFCandConnector|connect") << "pfCand.size()=" << pfCand.size() << "bCalibPrimary_=" << bCalibPrimary_;
73 
74  for (unsigned int ce1 = 0; ce1 < pfCand.size(); ++ce1) {
75  if (isPrimaryNucl(pfCand.at(ce1))) {
76  LogTrace("PFCandConnector|connect")
77  << "" << endl
78  << "Nuclear Interaction w Primary Candidate " << ce1 << " " << pfCand.at(ce1) << endl
79  << " based on the Track " << pfCand.at(ce1).trackRef().key()
80  << " w pT = " << pfCand.at(ce1).trackRef()->pt() << " #pm "
81  << pfCand.at(ce1).trackRef()->ptError() / pfCand.at(ce1).trackRef()->pt() * 100 << " %"
82  << " ECAL = " << pfCand.at(ce1).ecalEnergy() << " HCAL = " << pfCand.at(ce1).hcalEnergy() << endl;
83 
84 #ifdef EDM_ML_DEBUG
85  (pfCand.at(ce1)).displacedVertexRef(fT_TO_DISP_)->Dump();
86 #endif
87 
88  analyseNuclearWPrim(pfCand, bMask, ce1);
89 
90 #ifdef EDM_ML_DEBUG
91  LogTrace("PFCandConnector|connect")
92  << "After Connection the candidate " << ce1 << " is " << pfCand.at(ce1) << endl
93  << endl;
94 
95  PFCandidate::ElementsInBlocks elementsInBlocks = pfCand.at(ce1).elementsInBlocks();
96  for (unsigned blockElem = 0; blockElem < elementsInBlocks.size(); blockElem++) {
97  if (blockElem == 0)
98  LogTrace("PFCandConnector|connect") << *(elementsInBlocks[blockElem].first) << endl;
99  LogTrace("PFCandConnector|connect") << " position " << elementsInBlocks[blockElem].second;
100  }
101 #endif
102  }
103  }
104 
105  for (unsigned int ce1 = 0; ce1 < pfCand.size(); ++ce1) {
106  if (!bMask[ce1] && isSecondaryNucl(pfCand.at(ce1))) {
107  LogTrace("PFCandConnector|connect")
108  << "" << endl
109  << "Nuclear Interaction w no Primary Candidate " << ce1 << " " << pfCand.at(ce1) << endl
110  << " based on the Track " << pfCand.at(ce1).trackRef().key()
111  << " w pT = " << pfCand.at(ce1).trackRef()->pt() << " #pm " << pfCand.at(ce1).trackRef()->ptError() << " %"
112  << " ECAL = " << pfCand.at(ce1).ecalEnergy() << " HCAL = " << pfCand.at(ce1).hcalEnergy()
113  << " dE(Trk-CALO) = "
114  << pfCand.at(ce1).trackRef()->p() - pfCand.at(ce1).ecalEnergy() - pfCand.at(ce1).hcalEnergy()
115  << " Nmissing hits = "
116  << pfCand.at(ce1).trackRef()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_OUTER_HITS) << endl;
117 
118 #ifdef EDM_ML_DEBUG
119  (pfCand.at(ce1)).displacedVertexRef(fT_FROM_DISP_)->Dump();
120 #endif
121 
122  analyseNuclearWSec(pfCand, bMask, ce1);
123 
124 #ifdef EDM_ML_DEBUG
125  LogTrace("PFCandConnector|connect") << "After Connection the candidate " << ce1 << " is " << pfCand.at(ce1)
126  << " and elements connected to it are: " << endl;
127 
128  PFCandidate::ElementsInBlocks elementsInBlocks = pfCand.at(ce1).elementsInBlocks();
129  for (unsigned blockElem = 0; blockElem < elementsInBlocks.size(); blockElem++) {
130  if (blockElem == 0)
131  LogTrace("PFCandConnector|connect") << *(elementsInBlocks[blockElem].first) << endl;
132  LogTrace("PFCandConnector|connect") << " position " << elementsInBlocks[blockElem].second;
133  }
134 #endif
135  }
136  }
137  }
138 
139  for (unsigned int ce1 = 0; ce1 < pfCand.size(); ++ce1)
140  if (!bMask[ce1])
141  pfC.push_back(pfCand.at(ce1));
142 
143  LogTrace("PFCandConnector|connect") << "end of function";
144 
145  return pfC;
146 }
147 
148 void PFCandConnector::analyseNuclearWPrim(PFCandidateCollection& pfCand,
149  std::vector<bool>& bMask,
150  unsigned int ce1) const {
151  PFDisplacedVertexRef ref1, ref2, ref1_bis;
152 
153  PFCandidate primaryCand = pfCand.at(ce1);
154 
155  // ------- look for the little friends -------- //
156 
157  const math::XYZTLorentzVectorD& momentumPrim = primaryCand.p4();
158 
159  math::XYZTLorentzVectorD momentumSec;
160 
161  momentumSec = momentumPrim / momentumPrim.E() * (primaryCand.ecalEnergy() + primaryCand.hcalEnergy());
162 
163  map<double, math::XYZTLorentzVectorD> candidatesWithTrackExcess;
164  map<double, math::XYZTLorentzVectorD> candidatesWithoutCalo;
165 
166  ref1 = primaryCand.displacedVertexRef(fT_TO_DISP_);
167 
168  for (unsigned int ce2 = 0; ce2 < pfCand.size(); ++ce2) {
169  if (ce2 != ce1 && isSecondaryNucl(pfCand.at(ce2))) {
170  ref2 = (pfCand.at(ce2)).displacedVertexRef(fT_FROM_DISP_);
171 
172  if (ref1 == ref2) {
173  LogTrace("PFCandConnector|analyseNuclearWPrim")
174  << "\t here is a Secondary Candidate " << ce2 << " " << pfCand.at(ce2) << endl
175  << "\t based on the Track " << pfCand.at(ce2).trackRef().key()
176  << " w p = " << pfCand.at(ce2).trackRef()->p() << " w pT = " << pfCand.at(ce2).trackRef()->pt() << " #pm "
177  << pfCand.at(ce2).trackRef()->ptError() << " %"
178  << " ECAL = " << pfCand.at(ce2).ecalEnergy() << " HCAL = " << pfCand.at(ce2).hcalEnergy()
179  << " dE(Trk-CALO) = "
180  << pfCand.at(ce2).trackRef()->p() - pfCand.at(ce2).ecalEnergy() - pfCand.at(ce2).hcalEnergy()
181  << " Nmissing hits = "
182  << pfCand.at(ce2).trackRef()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_OUTER_HITS) << endl;
183 
184  if (isPrimaryNucl(pfCand.at(ce2))) {
185  LogTrace("PFCandConnector|analyseNuclearWPrim") << "\t\t but it is also a Primary Candidate " << ce2 << endl;
186 
187  ref1_bis = (pfCand.at(ce2)).displacedVertexRef(fT_TO_DISP_);
188  if (ref1_bis.isNonnull())
189  analyseNuclearWPrim(pfCand, bMask, ce2);
190  }
191 
192  // Take now the parameters of the secondary track that are relevant and use them to construct the NI candidate
193 
194  PFCandidate::ElementsInBlocks elementsInBlocks = pfCand.at(ce2).elementsInBlocks();
195  PFCandidate::ElementsInBlocks elementsAlreadyInBlocks = pfCand.at(ce1).elementsInBlocks();
196  for (unsigned blockElem = 0; blockElem < elementsInBlocks.size(); blockElem++) {
197  bool isAlreadyHere = false;
198  for (unsigned alreadyBlock = 0; alreadyBlock < elementsAlreadyInBlocks.size(); alreadyBlock++) {
199  if (elementsAlreadyInBlocks[alreadyBlock].second == elementsInBlocks[blockElem].second)
200  isAlreadyHere = true;
201  }
202  if (!isAlreadyHere)
203  pfCand.at(ce1).addElementInBlock(elementsInBlocks[blockElem].first, elementsInBlocks[blockElem].second);
204  }
205 
206  double caloEn = pfCand.at(ce2).ecalEnergy() + pfCand.at(ce2).hcalEnergy();
207  double deltaEn = pfCand.at(ce2).p4().E() - caloEn;
208  int nMissOuterHits =
209  pfCand.at(ce2).trackRef()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_OUTER_HITS);
210 
211  // Check if the difference Track Calo is not too large and if we can trust the track, ie it doesn't miss too much hits.
212  if (deltaEn > 1 && nMissOuterHits > 1) {
213  math::XYZTLorentzVectorD momentumToAdd = pfCand.at(ce2).p4() * caloEn / pfCand.at(ce2).p4().E();
214  momentumSec += momentumToAdd;
215  LogTrace("PFCandConnector|analyseNuclearWPrim")
216  << "The difference track-calo s really large and the track miss at least 2 hits. A secondary NI may "
217  "have happened. Let's trust the calo energy"
218  << endl
219  << "add " << momentumToAdd << endl;
220 
221  } else {
222  // Check if the difference Track Calo is not too large and if we can trust the track, ie it doesn't miss too much hits.
223  if (caloEn > 0.01 && deltaEn > 1 && nMissOuterHits > 0) {
224  math::XYZTLorentzVectorD momentumExcess = pfCand.at(ce2).p4() * deltaEn / pfCand.at(ce2).p4().E();
225  candidatesWithTrackExcess[pfCand.at(ce2).trackRef()->pt() / pfCand.at(ce2).trackRef()->ptError()] =
226  momentumExcess;
227  } else if (caloEn < 0.01)
228  candidatesWithoutCalo[pfCand.at(ce2).trackRef()->pt() / pfCand.at(ce2).trackRef()->ptError()] =
229  pfCand.at(ce2).p4();
230  momentumSec += (pfCand.at(ce2)).p4();
231  }
232 
233  bMask[ce2] = true;
234  }
235  }
236  }
237 
238  // We have more primary energy than secondary: reject all secondary tracks which have no calo energy attached.
239 
240  if (momentumPrim.E() < momentumSec.E()) {
241  LogTrace("PFCandConnector|analyseNuclearWPrim")
242  << "Size of 0 calo Energy secondary candidates" << candidatesWithoutCalo.size() << endl;
243  for (map<double, math::XYZTLorentzVectorD>::iterator iter = candidatesWithoutCalo.begin();
244  iter != candidatesWithoutCalo.end() && momentumPrim.E() < momentumSec.E();
245  iter++)
246  if (momentumSec.E() > iter->second.E() + 0.1) {
247  momentumSec -= iter->second;
248 
249  LogTrace("PFCandConnector|analyseNuclearWPrim")
250  << "\t Remove a SecondaryCandidate with 0 calo energy " << iter->second << endl;
251  LogTrace("PFCandConnector|analyseNuclearWPrim")
252  << "momentumPrim.E() = " << momentumPrim.E() << " and momentumSec.E() = " << momentumSec.E() << endl;
253  }
254  }
255 
256  if (momentumPrim.E() < momentumSec.E()) {
257  LogTrace("PFCandConnector|analyseNuclearWPrim")
258  << "0 Calo Energy rejected but still not sufficient. Size of not enough calo Energy secondary candidates"
259  << candidatesWithTrackExcess.size() << endl;
260  for (map<double, math::XYZTLorentzVectorD>::iterator iter = candidatesWithTrackExcess.begin();
261  iter != candidatesWithTrackExcess.end() && momentumPrim.E() < momentumSec.E();
262  iter++)
263  if (momentumSec.E() > iter->second.E() + 0.1)
264  momentumSec -= iter->second;
265  }
266 
267  double dpt = pfCand.at(ce1).trackRef()->ptError() / pfCand.at(ce1).trackRef()->pt() * 100;
268 
269  if (momentumSec.E() < 0.1) {
270  bMask[ce1] = true;
271  return;
272  }
273 
274  // Rescale the secondary candidates to account for the loss of energy, but only if we can trust the primary track:
275  // if it has more energy than secondaries and is precise enough and secondary exist and was not eaten or rejected during the PFAlgo step.
276 
277  if (((ref1->isTherePrimaryTracks() && dpt < dptRel_PrimaryTrack_) ||
278  (ref1->isThereMergedTracks() && dpt < dptRel_MergedTrack_)) &&
279  momentumPrim.E() > momentumSec.E() && momentumSec.E() > 0.1) {
280  if (bCalibPrimary_) {
281  double factor = rescaleFactor(momentumPrim.Pt(), momentumSec.E() / momentumPrim.E());
282  LogTrace("PFCandConnector|analyseNuclearWPrim") << "factor = " << factor << endl;
283  if (factor * momentumPrim.Pt() < momentumSec.Pt())
284  momentumSec = momentumPrim;
285  else
286  momentumSec += (1 - factor) * momentumPrim;
287  }
288 
289  double px = momentumPrim.Px() * momentumSec.P() / momentumPrim.P();
290  double py = momentumPrim.Py() * momentumSec.P() / momentumPrim.P();
291  double pz = momentumPrim.Pz() * momentumSec.P() / momentumPrim.P();
292  double E = sqrt(px * px + py * py + pz * pz + pion_mass2);
293  math::XYZTLorentzVectorD momentum(px, py, pz, E);
294  pfCand.at(ce1).setP4(momentum);
295 
296  return;
297 
298  } else {
299  math::XYZVector primDir = ref1->primaryDirection();
300 
301  if (primDir.Mag2() < 0.1) {
302  // It might be 0 but this situation should never happend. Throw a warning if it happens.
303  edm::LogWarning("PFCandConnector") << "A Nuclear Interaction do not have primary direction" << std::endl;
304  pfCand.at(ce1).setP4(momentumSec);
305  return;
306  } else {
307  // rescale the primary direction to the optimal momentum. But take care of the factthat it shall not be completly 0 to avoid a warning if Jet Area.
308  double momentumS = momentumSec.P();
309  if (momentumS < 1e-4)
310  momentumS = 1e-4;
311  double px = momentumS * primDir.x();
312  double py = momentumS * primDir.y();
313  double pz = momentumS * primDir.z();
314  double E = sqrt(px * px + py * py + pz * pz + pion_mass2);
315 
316  math::XYZTLorentzVectorD momentum(px, py, pz, E);
317  pfCand.at(ce1).setP4(momentum);
318  return;
319  }
320  }
321 }
322 
323 void PFCandConnector::analyseNuclearWSec(PFCandidateCollection& pfCand,
324  std::vector<bool>& bMask,
325  unsigned int ce1) const {
326  PFDisplacedVertexRef ref1, ref2;
327 
328  // Check if the track excess was not too large and track may miss some outer hits. This may point to a secondary NI.
329 
330  double caloEn = pfCand.at(ce1).ecalEnergy() + pfCand.at(ce1).hcalEnergy();
331  double deltaEn = pfCand.at(ce1).p4().E() - caloEn;
332  int nMissOuterHits = pfCand.at(ce1).trackRef()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_OUTER_HITS);
333 
334  ref1 = pfCand.at(ce1).displacedVertexRef(fT_FROM_DISP_);
335 
336  // ------- check if an electron or a muon vas spotted as incoming track -------- //
337  // ------- this mean probably that the NI was fake thus we do not correct it -------- /
338 
339  if (ref1->isTherePrimaryTracks() || ref1->isThereMergedTracks()) {
340  std::vector<reco::Track> refittedTracks = ref1->refittedTracks();
341  for (unsigned it = 0; it < refittedTracks.size(); it++) {
342  reco::TrackBaseRef primaryBaseRef = ref1->originalTrack(refittedTracks[it]);
343  if (ref1->isIncomingTrack(primaryBaseRef))
344  LogTrace("PFCandConnector|analyseNuclearWSec")
345  << "There is a Primary track ref with pt = " << primaryBaseRef->pt() << endl;
346 
347  for (unsigned int ce = 0; ce < pfCand.size(); ++ce) {
348  // cout << "PFCand Id = " << (pfCand.at(ce)).particleId() << endl;
349  if ((pfCand.at(ce)).particleId() == reco::PFCandidate::e ||
350  (pfCand.at(ce)).particleId() == reco::PFCandidate::mu) {
351  LogTrace("PFCandConnector|analyseNuclearWSec")
352  << " It is an electron and it has a ref to a track " << (pfCand.at(ce)).trackRef().isNonnull() << endl;
353 
354  if ((pfCand.at(ce)).trackRef().isNonnull()) {
355  reco::TrackRef tRef = (pfCand.at(ce)).trackRef();
356  reco::TrackBaseRef bRef(tRef);
357  LogTrace("PFCandConnector|analyseNuclearWSec")
358  << "With Track Ref pt = " << (pfCand.at(ce)).trackRef()->pt() << endl;
359 
360  if (bRef == primaryBaseRef) {
361  if ((pfCand.at(ce)).particleId() == reco::PFCandidate::e)
362  LogTrace("PFCandConnector|analyseNuclearWSec")
363  << "It is a NI from electron. NI Discarded. Just release the candidate." << endl;
364  if ((pfCand.at(ce)).particleId() == reco::PFCandidate::mu)
365  LogTrace("PFCandConnector|analyseNuclearWSec")
366  << "It is a NI from muon. NI Discarded. Just release the candidate" << endl;
367 
368  // release the track but take care of not overcounting bad tracks. In fact those tracks was protected against destruction in
369  // PFAlgo. Now we treat them as if they was treated in PFAlgo
370 
371  if (caloEn < 0.1 && pfCand.at(ce1).trackRef()->ptError() > ptErrorSecondary_) {
372  edm::LogInfo("PFCandConnector|analyseNuclearWSec")
373  << "discarded track since no calo energy and ill measured" << endl;
374  bMask[ce1] = true;
375  }
376  if (caloEn > 0.1 && deltaEn > ptErrorSecondary_ &&
377  pfCand.at(ce1).trackRef()->ptError() > ptErrorSecondary_) {
378  edm::LogInfo("PFCandConnector|analyseNuclearWSec")
379  << "rescaled momentum of the track since no calo energy and ill measured" << endl;
380 
381  double factor = caloEn / pfCand.at(ce1).p4().E();
382  pfCand.at(ce1).rescaleMomentum(factor);
383  }
384 
385  return;
386  }
387  }
388  }
389  }
390  }
391  }
392 
393  PFCandidate secondaryCand = pfCand.at(ce1);
394 
395  math::XYZTLorentzVectorD momentumSec = secondaryCand.p4();
396 
397  if (deltaEn > ptErrorSecondary_ && nMissOuterHits > 1) {
398  math::XYZTLorentzVectorD momentumToAdd = pfCand.at(ce1).p4() * caloEn / pfCand.at(ce1).p4().E();
399  momentumSec = momentumToAdd;
400  LogTrace("PFCandConnector|analyseNuclearWSec")
401  << "The difference track-calo s really large and the track miss at least 2 hits. A secondary NI may have "
402  "happened. Let's trust the calo energy"
403  << endl
404  << "add " << momentumToAdd << endl;
405  }
406 
407  // ------- look for the little friends -------- //
408  for (unsigned int ce2 = ce1 + 1; ce2 < pfCand.size(); ++ce2) {
409  if (isSecondaryNucl(pfCand.at(ce2))) {
410  ref2 = (pfCand.at(ce2)).displacedVertexRef(fT_FROM_DISP_);
411 
412  if (ref1 == ref2) {
413  LogTrace("PFCandConnector|analyseNuclearWSec")
414  << "\t here is a Secondary Candidate " << ce2 << " " << pfCand.at(ce2) << endl
415  << "\t based on the Track " << pfCand.at(ce2).trackRef().key()
416  << " w pT = " << pfCand.at(ce2).trackRef()->pt() << " #pm " << pfCand.at(ce2).trackRef()->ptError() << " %"
417  << " ECAL = " << pfCand.at(ce2).ecalEnergy() << " HCAL = " << pfCand.at(ce2).hcalEnergy()
418  << " dE(Trk-CALO) = "
419  << pfCand.at(ce2).trackRef()->p() - pfCand.at(ce2).ecalEnergy() - pfCand.at(ce2).hcalEnergy()
420  << " Nmissing hits = "
421  << pfCand.at(ce2).trackRef()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_OUTER_HITS) << endl;
422 
423  // Take now the parameters of the secondary track that are relevant and use them to construct the NI candidate
424  PFCandidate::ElementsInBlocks elementsInBlocks = pfCand.at(ce2).elementsInBlocks();
425  PFCandidate::ElementsInBlocks elementsAlreadyInBlocks = pfCand.at(ce1).elementsInBlocks();
426  for (unsigned blockElem = 0; blockElem < elementsInBlocks.size(); blockElem++) {
427  bool isAlreadyHere = false;
428  for (unsigned alreadyBlock = 0; alreadyBlock < elementsAlreadyInBlocks.size(); alreadyBlock++) {
429  if (elementsAlreadyInBlocks[alreadyBlock].second == elementsInBlocks[blockElem].second)
430  isAlreadyHere = true;
431  }
432  if (!isAlreadyHere)
433  pfCand.at(ce1).addElementInBlock(elementsInBlocks[blockElem].first, elementsInBlocks[blockElem].second);
434  }
435 
436  double caloEn = pfCand.at(ce2).ecalEnergy() + pfCand.at(ce2).hcalEnergy();
437  double deltaEn = pfCand.at(ce2).p4().E() - caloEn;
438  int nMissOuterHits =
439  pfCand.at(ce2).trackRef()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_OUTER_HITS);
440  if (deltaEn > ptErrorSecondary_ && nMissOuterHits > 1) {
441  math::XYZTLorentzVectorD momentumToAdd = pfCand.at(ce2).p4() * caloEn / pfCand.at(ce2).p4().E();
442  momentumSec += momentumToAdd;
443  LogTrace("PFCandConnector|analyseNuclearWSec")
444  << "The difference track-calo s really large and the track miss at least 2 hits. A secondary NI may "
445  "have happened. Let's trust the calo energy"
446  << endl
447  << "add " << momentumToAdd << endl;
448  } else {
449  momentumSec += (pfCand.at(ce2)).p4();
450  }
451 
452  bMask[ce2] = true;
453  }
454  }
455  }
456 
457  math::XYZVector primDir = ref1->primaryDirection();
458 
459  if (primDir.Mag2() < 0.1) {
460  // It might be 0 but this situation should never happend. Throw a warning if it happens.
461  pfCand.at(ce1).setP4(momentumSec);
462  edm::LogWarning("PFCandConnector") << "A Nuclear Interaction do not have primary direction" << std::endl;
463  return;
464  } else {
465  // rescale the primary direction to the optimal momentum. But take care of the factthat it shall not be completly 0 to avoid a warning if Jet Area.
466  double momentumS = momentumSec.P();
467  if (momentumS < 1e-4)
468  momentumS = 1e-4;
469  double px = momentumS * primDir.x();
470  double py = momentumS * primDir.y();
471  double pz = momentumS * primDir.z();
472  double E = sqrt(px * px + py * py + pz * pz + pion_mass2);
473 
474  math::XYZTLorentzVectorD momentum(px, py, pz, E);
475 
476  pfCand.at(ce1).setP4(momentum);
477  return;
478  }
479 }
480 
481 bool PFCandConnector::isSecondaryNucl(const PFCandidate& pf) const {
482  PFDisplacedVertexRef ref1;
483  // nuclear
484  if (pf.flag(fT_FROM_DISP_)) {
485  ref1 = pf.displacedVertexRef(fT_FROM_DISP_);
486  // ref1->Dump();
487  if (!ref1.isNonnull())
488  return false;
489  else if (ref1->isNucl() || ref1->isNucl_Loose() || ref1->isNucl_Kink())
490  return true;
491  }
492 
493  return false;
494 }
495 
496 bool PFCandConnector::isPrimaryNucl(const PFCandidate& pf) const {
497  PFDisplacedVertexRef ref1;
498 
499  // nuclear
500  if (pf.flag(fT_TO_DISP_)) {
501  ref1 = pf.displacedVertexRef(fT_TO_DISP_);
502  //ref1->Dump();
503 
504  if (!ref1.isNonnull())
505  return false;
506  else if (ref1->isNucl() || ref1->isNucl_Loose() || ref1->isNucl_Kink())
507  return true;
508  }
509 
510  return false;
511 }
512 
513 double PFCandConnector::rescaleFactor(const double pt, const double cFrac) const {
514  /*
515  LOG NORMAL FIT
516  FCN=35.8181 FROM MIGRAD STATUS=CONVERGED 257 CALLS 258 TOTAL
517  EDM=8.85763e-09 STRATEGY= 1 ERROR MATRIX ACCURATE
518  EXT PARAMETER STEP FIRST
519  NO. NAME VALUE ERROR SIZE DERIVATIVE
520  1 p0 7.99434e-01 2.77264e-02 6.59108e-06 9.80247e-03
521  2 p1 1.51303e-01 2.89981e-02 1.16775e-05 6.99035e-03
522  3 p2 -5.03829e-01 2.87929e-02 1.90070e-05 1.37015e-03
523  4 p3 4.54043e-01 5.00908e-02 3.17625e-05 3.86622e-03
524  5 p4 -4.61736e-02 8.07940e-03 3.25775e-06 -1.37247e-02
525  */
526 
527  /*
528  FCN=34.4051 FROM MIGRAD STATUS=CONVERGED 221 CALLS 222 TOTAL
529  EDM=1.02201e-09 STRATEGY= 1 ERROR MATRIX UNCERTAINTY 2.3 per cent
530 
531  fConst
532  1 p0 7.99518e-01 2.23519e-02 1.41523e-06 4.05975e-04
533  2 p1 1.44619e-01 2.39398e-02 -7.68117e-07 -2.55775e-03
534 
535  fNorm
536  3 p2 -5.16571e-01 3.12362e-02 5.74932e-07 3.42292e-03
537  4 p3 4.69055e-01 5.09665e-02 1.94353e-07 1.69031e-03
538 
539  fExp
540  5 p4 -5.18044e-02 8.13458e-03 4.29815e-07 -1.07624e-02
541  */
542 
543  double fConst, fNorm, fExp;
544 
545  fConst = fConst_[0] + fConst_[1] * cFrac;
546  fNorm = fNorm_[0] - fNorm_[1] * cFrac;
547  fExp = fExp_[0];
548 
549  double factor = fConst - fNorm * exp(-fExp * pt);
550 
551  return factor;
552 }
553 
554 void PFCandConnector::fillPSetDescription(edm::ParameterSetDescription& iDesc) {
555  iDesc.add<bool>("bCorrect", true);
556  iDesc.add<bool>("bCalibPrimary", true);
557  iDesc.add<double>("dptRel_PrimaryTrack", 10.0);
558  iDesc.add<double>("dptRel_MergedTrack", 5.0);
559  iDesc.add<double>("ptErrorSecondary", 1.0);
560  iDesc.add<std::vector<double>>("nuclCalibFactors", {0.8, 0.15, 0.5, 0.5, 0.05});
561 }
reco::HitPattern::MISSING_OUTER_HITS
Definition: HitPattern.h:155
edm::ParameterSetDescription::add
ParameterDescriptionBase * add(U const &iLabel, T const &value)
Definition: ParameterSetDescription.h:95
PFCandConnector.h
MessageLogger.h
HLT_2018_cff.nuclCalibFactors
nuclCalibFactors
Definition: HLT_2018_cff.py:12317
reco::PFCandidate::e
Definition: PFCandidate.h:46
DiDispStaMuonMonitor_cfi.pt
pt
Definition: DiDispStaMuonMonitor_cfi.py:39
multPhiCorr_741_25nsDY_cfi.py
py
Definition: multPhiCorr_741_25nsDY_cfi.py:12
PFCandConnector::analyseNuclearWSec
void analyseNuclearWSec(reco::PFCandidateCollection &, std::vector< bool > &, unsigned int) const
Analyse nuclear interactions where a secondary track is present.
Definition: PFCandConnector.cc:323
math::XYZTLorentzVectorD
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > XYZTLorentzVectorD
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:14
edm::LogInfo
Definition: MessageLogger.h:254
edm::ParameterSetDescription
Definition: ParameterSetDescription.h:52
edm::second
U second(std::pair< T, U > const &p)
Definition: ParameterSet.cc:215
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:45
dqmdumpme.first
first
Definition: dqmdumpme.py:55
edm::Ref< PFDisplacedVertexCollection >
HLT_2018_cff.bCorrect
bCorrect
Definition: HLT_2018_cff.py:12318
reco::TrackBase::pt
double pt() const
track transverse momentum
Definition: TrackBase.h:608
PFCandConnector::isPrimaryNucl
bool isPrimaryNucl(const reco::PFCandidate &pf) const
Definition: PFCandConnector.cc:496
reco::PFCandidate::displacedVertexRef
reco::PFDisplacedVertexRef displacedVertexRef(Flags type) const
Definition: PFCandidate.cc:481
PFCandConnector::analyseNuclearWPrim
void analyseNuclearWPrim(reco::PFCandidateCollection &, std::vector< bool > &, unsigned int) const
Analyse nuclear interactions where a primary or merged track is present.
Definition: PFCandConnector.cc:148
reco::PFCandidate::mu
Definition: PFCandidate.h:47
HLT_2018_cff.bCalibPrimary
bCalibPrimary
Definition: HLT_2018_cff.py:12319
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
DQMScaleToClient_cfi.factor
factor
Definition: DQMScaleToClient_cfi.py:8
reco::PFCandidate::ElementsInBlocks
std::vector< ElementInBlock > ElementsInBlocks
Definition: PFCandidate.h:378
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
edm::LogWarning
Definition: MessageLogger.h:141
PFCandConnector::setParameters
void setParameters(const edm::ParameterSet &iCfgCandConnector)
Definition: PFCandConnector.h:30
PFDisplacedVertex.h
edm::Ref::isNonnull
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:238
reco::PFCandidate::T_FROM_DISP
Definition: PFCandidate.h:68
math::XYZVector
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:31
reco::LeafCandidate::p4
const LorentzVector & p4() const final
four-momentum Lorentz vector
Definition: LeafCandidate.h:114
p4
double p4[4]
Definition: TauolaWrapper.h:92
GsfTrack.h
reco::PFCandidate::ecalEnergy
double ecalEnergy() const
return corrected Ecal energy
Definition: PFCandidate.h:220
tkAlTrackRefitSequence_cff.refittedTracks
refittedTracks
Definition: tkAlTrackRefitSequence_cff.py:6
PFBlock.h
multPhiCorr_741_25nsDY_cfi.px
px
Definition: multPhiCorr_741_25nsDY_cfi.py:10
std
Definition: JetResolutionObject.h:76
reco::PFCandidate::hcalEnergy
double hcalEnergy() const
return corrected Hcal energy
Definition: PFCandidate.h:232
packedPFCandidateRefMixer_cfi.pf
pf
Definition: packedPFCandidateRefMixer_cfi.py:4
reco::PFCandidate::T_TO_DISP
Definition: PFCandidate.h:67
reco::PFCandidateCollection
std::vector< reco::PFCandidate > PFCandidateCollection
collection of PFCandidates
Definition: PFCandidateFwd.h:12
PFCandConnector::isSecondaryNucl
bool isSecondaryNucl(const reco::PFCandidate &pf) const
Definition: PFCandConnector.cc:481
PFCandConnector::fT_TO_DISP_
static const reco::PFCandidate::Flags fT_TO_DISP_
Definition: PFCandConnector.h:115
edm::RefToBase< reco::Track >
PFCandConnector::fT_FROM_DISP_
static const reco::PFCandidate::Flags fT_FROM_DISP_
Definition: PFCandConnector.h:116
reco::PFCandidate
Particle reconstructed by the particle flow algorithm.
Definition: PFCandidate.h:40
LogTrace
#define LogTrace(id)
Definition: MessageLogger.h:671
PFCandConnector::connect
reco::PFCandidateCollection connect(reco::PFCandidateCollection &pfCand) const
Definition: PFCandConnector.cc:63
JetChargeProducer_cfi.exp
exp
Definition: JetChargeProducer_cfi.py:6
PFCandConnector::rescaleFactor
double rescaleFactor(const double pt, const double cFrac) const
Return a calibration factor for a reconstructed nuclear interaction.
Definition: PFCandConnector.cc:513
reco::PFCandidate::Flags
Flags
Definition: PFCandidate.h:54
PFCandConnector::fillPSetDescription
static void fillPSetDescription(edm::ParameterSetDescription &iDesc)
Definition: PFCandConnector.cc:554
PFCandConnector::pion_mass2
static const double pion_mass2
Useful constants.
Definition: PFCandConnector.h:114
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
Generated for CMSSW Reference Manual by doxygen 1.8.16