test
CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
KDTreeLinkerTrackHcal.cc
Go to the documentation of this file.
2 
4 #include "TMath.h"
5 
6 // the text name is different so that we can easily
7 // construct it when calling the factory
10  "KDTreeTrackAndHCALLinker");
11 
12 
15 {
17  setPhiOffset(0.32);
18 }
19 
21 {
22  clear();
23 }
24 
25 void
27 {
28  if( track->trackRefPF()->extrapolatedPoint( reco::PFTrajectoryPoint::HCALEntrance ).isValid() ) {
29  targetSet_.insert(track);
30  }
31 }
32 
33 
34 void
36 {
37  reco::PFClusterRef clusterref = hcalCluster->clusterRef();
38 
39  // This test is more or less done in PFBlockAlgo.h. In others cases, it should be switch on.
40  // if (!((clusterref->layer() == PFLayer::HCAL_ENDCAP) ||
41  // (clusterref->layer() == PFLayer::HCAL_BARREL1)))
42  // return;
43 
44  const std::vector<reco::PFRecHitFraction> &fraction = clusterref->recHitFractions();
45 
46  // We create a list of hcalCluster
47  fieldClusterSet_.insert(hcalCluster);
48  for(size_t rhit = 0; rhit < fraction.size(); ++rhit) {
49  const reco::PFRecHitRef& rh = fraction[rhit].recHitRef();
50  double fract = fraction[rhit].fraction();
51 
52  if ((rh.isNull()) || (fract < 1E-4))
53  continue;
54 
55  const reco::PFRecHit& rechit = *rh;
56 
57  // We save the links rechit to HcalClusters
58  rechit2ClusterLinks_[&rechit].insert(hcalCluster);
59 
60  // We create a liste of rechits
61  rechitsSet_.insert(&rechit);
62  }
63 }
64 
65 void
67 {
68  // List of pseudo-rechits that will be used to create the KDTree
69  std::vector<KDTreeNodeInfo> eltList;
70 
71  // Filling of this list
72  for(RecHitSet::const_iterator it = rechitsSet_.begin();
73  it != rechitsSet_.end(); it++) {
74 
75  const reco::PFRecHit::REPPoint &posrep = (*it)->positionREP();
76 
77  KDTreeNodeInfo rh1 (*it, posrep.eta(), posrep.phi());
78  eltList.push_back(rh1);
79 
80  // Here we solve the problem of phi circular set by duplicating some rechits
81  // too close to -Pi (or to Pi) and adding (substracting) to them 2 * Pi.
82  if (rh1.dim2 > (M_PI - getPhiOffset())) {
83  double phi = rh1.dim2 - 2 * M_PI;
84  KDTreeNodeInfo rh2(*it, posrep.eta(), phi);
85  eltList.push_back(rh2);
86  }
87 
88  if (rh1.dim2 < (M_PI * -1.0 + getPhiOffset())) {
89  double phi = rh1.dim2 + 2 * M_PI;
90  KDTreeNodeInfo rh3(*it, posrep.eta(), phi);
91  eltList.push_back(rh3);
92  }
93  }
94 
95  // Here we define the upper/lower bounds of the 2D space (eta/phi).
96  double phimin = -1.0 * M_PI - getPhiOffset();
97  double phimax = M_PI + getPhiOffset();
98 
99  // etamin-etamax, phimin-phimax
100  KDTreeBox region(-3.0, 3.0, phimin, phimax);
101 
102  // We may now build the KDTree
103  tree_.build(eltList, region);
104 }
105 
106 void
108 {
109  // Must of the code has been taken from LinkByRecHit.cc
110 
111  // We iterate over the tracks.
112  for(BlockEltSet::iterator it = targetSet_.begin();
113  it != targetSet_.end(); it++) {
114 
115  reco::PFRecTrackRef trackref = (*it)->trackRefPF();
116 
117  const reco::PFTrajectoryPoint& atHCAL =
118  trackref->extrapolatedPoint(reco::PFTrajectoryPoint::HCALEntrance);
119  const reco::PFTrajectoryPoint& atHCALExit =
120  trackref->extrapolatedPoint(reco::PFTrajectoryPoint::HCALExit);
121 
122  // The track didn't reach hcal
123  if( ! atHCAL.isValid()) continue;
124 
125  double dHeta = atHCALExit.positionREP().eta() - atHCAL.positionREP().eta();
126  double dHphi = atHCALExit.positionREP().phi() - atHCAL.positionREP().phi();
127  if ( dHphi > M_PI ) dHphi = dHphi - 2. * M_PI;
128  else if ( dHphi < -M_PI ) dHphi = dHphi + 2. * M_PI;
129 
130  double tracketa = atHCAL.positionREP().eta() + 0.1 * dHeta;
131  double trackphi = atHCAL.positionREP().phi() + 0.1 * dHphi;
132 
133  if (trackphi > M_PI) trackphi -= 2 * M_PI;
134  else if (trackphi < -M_PI) trackphi += 2 * M_PI;
135 
136  // Estimate the maximal envelope in phi/eta that will be used to find rechit candidates.
137  // Same envelope for cap et barrel rechits.
138  double inflation = 1.;
139  double rangeeta = (getCristalPhiEtaMaxSize() * (1.5 + 0.5) + 0.2 * fabs(dHeta)) * inflation;
140  double rangephi = (getCristalPhiEtaMaxSize() * (1.5 + 0.5) + 0.2 * fabs(dHphi)) * inflation;
141 
142  // We search for all candidate recHits, ie all recHits contained in the maximal size envelope.
143  std::vector<KDTreeNodeInfo> recHits;
144  KDTreeBox trackBox(tracketa - rangeeta, tracketa + rangeeta,
145  trackphi - rangephi, trackphi + rangephi);
146  tree_.search(trackBox, recHits);
147 
148  // Here we check all rechit candidates using the non-approximated method.
149  for(std::vector<KDTreeNodeInfo>::const_iterator rhit = recHits.begin();
150  rhit != recHits.end(); ++rhit) {
151 
152  const auto &rhrep = rhit->ptr->positionREP();
153  const auto & corners = rhit->ptr->getCornersREP();
154 
155  double rhsizeeta = fabs(corners[3].eta() - corners[1].eta());
156  double rhsizephi = fabs(corners[3].phi() - corners[1].phi());
157  if ( rhsizephi > M_PI ) rhsizephi = 2.*M_PI - rhsizephi;
158 
159  double deta = fabs(rhrep.eta() - tracketa);
160  double dphi = fabs(rhrep.phi() - trackphi);
161  if ( dphi > M_PI ) dphi = 2.*M_PI - dphi;
162 
163  // Find all clusters associated to given rechit
164  RecHit2BlockEltMap::iterator ret = rechit2ClusterLinks_.find(rhit->ptr);
165 
166  for(BlockEltSet::iterator clusterIt = ret->second.begin();
167  clusterIt != ret->second.end(); clusterIt++) {
168 
169  const reco::PFClusterRef clusterref = (*clusterIt)->clusterRef();
170  int fracsNbr = clusterref->recHitFractions().size();
171 
172  double _rhsizeeta = rhsizeeta * (1.5 + 0.5 / fracsNbr) + 0.2 * fabs(dHeta);
173  double _rhsizephi = rhsizephi * (1.5 + 0.5 / fracsNbr) + 0.2 * fabs(dHphi);
174 
175  // Check if the track and the cluster are linked
176  if(deta < (_rhsizeeta / 2.) && dphi < (_rhsizephi / 2.))
177  cluster2TargetLinks_[*clusterIt].insert(*it);
178  }
179  }
180  }
181 }
182 
183 void
185 {
186  //TODO YG : Check if cluster positionREP() is valid ?
187 
188  // Here we save in each HCAL cluster the list of phi/eta values of linked clusters.
189  for (BlockElt2BlockEltMap::iterator it = cluster2TargetLinks_.begin();
190  it != cluster2TargetLinks_.end(); ++it) {
191  reco::PFMultiLinksTC multitracks(true);
192 
193  for (BlockEltSet::iterator jt = it->second.begin();
194  jt != it->second.end(); ++jt) {
195 
196  reco::PFRecTrackRef trackref = (*jt)->trackRefPF();
197  const reco::PFTrajectoryPoint& atHCAL =
198  trackref->extrapolatedPoint(reco::PFTrajectoryPoint::HCALEntrance);
199  double tracketa = atHCAL.positionREP().eta();
200  double trackphi = atHCAL.positionREP().phi();
201 
202  multitracks.linkedClusters.push_back(std::make_pair(trackphi, tracketa));
203  }
204 
205  it->first->setMultilinks(multitracks);
206  }
207 
208  // We set the multilinks flag of the track to true. It will allow us to
209  // use in an optimized way our algo results in the recursive linking algo.
210  for (BlockEltSet::iterator it = fieldClusterSet_.begin();
211  it != fieldClusterSet_.end(); ++it)
212  (*it)->setIsValidMultilinks(true);
213 
214 }
215 
216 void
218 {
219  targetSet_.clear();
220  fieldClusterSet_.clear();
221 
222  rechitsSet_.clear();
223 
224  rechit2ClusterLinks_.clear();
225  cluster2TargetLinks_.clear();
226 
227  tree_.clear();
228 }
float phi() const
momentum azimuthal angle
Definition: PtEtaPhiMass.h:53
const REPPoint & positionREP() const
trajectory position in (rho, eta, phi) base
void insertFieldClusterElt(reco::PFBlockElement *hcalCluster)
Abstract base class for a PFBlock element (track, cluster...)
tuple ret
prodAgent to be discontinued
void build(std::vector< KDTreeNodeInfoT< DATA, DIM > > &eltList, const KDTreeBoxT< DIM > &region)
virtual const PFClusterRef & clusterRef() const
void search(const KDTreeBoxT< DIM > &searchBox, std::vector< KDTreeNodeInfoT< DATA, DIM > > &resRecHitList)
Particle flow rechit (rechit + geometry and topology information). See clustering algorithm in PFClus...
Definition: PFRecHit.h:31
RecHit2BlockEltMap rechit2ClusterLinks_
void setPhiOffset(double phiOffset)
PFMultilinksType linkedClusters
BlockElt2BlockEltMap cluster2TargetLinks_
const Fraction< n, m >::type & fract()
Definition: Fraction.h:38
virtual const PFRecTrackRef & trackRefPF() const
bool isNull() const
Checks for null.
Definition: Ref.h:249
void insertTargetElt(reco::PFBlockElement *track)
#define M_PI
void setCristalPhiEtaMaxSize(float size)
bool isValid() const
is this point valid ?
float getCristalPhiEtaMaxSize() const
#define DEFINE_EDM_PLUGIN(factory, type, name)
A PFTrack holds several trajectory points, which basically contain the position and momentum of a tra...
float eta() const
momentum pseudorapidity
Definition: PtEtaPhiMass.h:51
float getPhiOffset() const