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