CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
LinkByRecHit.cc
Go to the documentation of this file.
4 #include "TMath.h"
5 
6 // to enable debugs
7 //#define PFLOW_DEBUG
8 
9 double
11  const reco::PFCluster& cluster,
12  bool isBrem,
13  bool debug) {
14 
15 #ifdef PFLOW_DEBUG
16  if( debug )
17  std::cout<<"entering test link by rechit function"<<std::endl;
18 #endif
19 
20  //cluster position
21  double clustereta = cluster.positionREP().Eta();
22  double clusterphi = cluster.positionREP().Phi();
23  double clusterX = cluster.position().X();
24  double clusterY = cluster.position().Y();
25  double clusterZ = cluster.position().Z();
26 
27  bool barrel = false;
28  bool hcal = false;
29  double distance = 999999.9;
30 
31  //track extrapolation
32  const reco::PFTrajectoryPoint& atVertex
34  const reco::PFTrajectoryPoint& atECAL
36 
37 
38  //track at calo's
39  double tracketa = 999999.9;
40  double trackphi = 999999.9;
41  double track_X = 999999.9;
42  double track_Y = 999999.9;
43  double track_Z = 999999.9;
44  double dHEta = 0.;
45  double dHPhi = 0.;
46 
47  // Quantities at vertex
48  double trackPt = isBrem ? 999. : sqrt(atVertex.momentum().Vect().Perp2());
49  // double trackEta = isBrem ? 999. : atVertex.momentum().Vect().Eta();
50 
51 
52  switch (cluster.layer()) {
53  case PFLayer::ECAL_BARREL: barrel = true;
55 #ifdef PFLOW_DEBUG
56  if( debug )
57  std::cout << "Fetching Ecal Resolution Maps"
58  << std::endl;
59 #endif
60  // did not reach ecal, cannot be associated with a cluster.
61  if( ! atECAL.isValid() ) return -1.;
62 
63  tracketa = atECAL.positionREP().Eta();
64  trackphi = atECAL.positionREP().Phi();
65  track_X = atECAL.position().X();
66  track_Y = atECAL.position().Y();
67  track_Z = atECAL.position().Z();
68 
69  distance
70  = std::sqrt( (track_X-clusterX)*(track_X-clusterX)
71  +(track_Y-clusterY)*(track_Y-clusterY)
72  +(track_Z-clusterZ)*(track_Z-clusterZ)
73  );
74 
75  break;
76 
77  case PFLayer::HCAL_BARREL1: barrel = true;
79 #ifdef PFLOW_DEBUG
80  if( debug )
81  std::cout << "Fetching Hcal Resolution Maps"
82  << std::endl;
83 #endif
84  if( isBrem ) {
85  return -1.;
86  } else {
87  hcal=true;
88  const reco::PFTrajectoryPoint& atHCAL
90  const reco::PFTrajectoryPoint& atHCALExit
92  // did not reach hcal, cannot be associated with a cluster.
93  if( ! atHCAL.isValid() ) return -1.;
94 
95  // The link is computed between 0 and ~1 interaction length in HCAL
96  dHEta = atHCALExit.positionREP().Eta()-atHCAL.positionREP().Eta();
97  dHPhi = atHCALExit.positionREP().Phi()-atHCAL.positionREP().Phi();
98  if ( dHPhi > M_PI ) dHPhi = dHPhi - 2.*M_PI;
99  else if ( dHPhi < -M_PI ) dHPhi = dHPhi + 2.*M_PI;
100  tracketa = atHCAL.positionREP().Eta() + 0.1*dHEta;
101  trackphi = atHCAL.positionREP().Phi() + 0.1*dHPhi;
102  track_X = atHCAL.position().X();
103  track_Y = atHCAL.position().Y();
104  track_Z = atHCAL.position().Z();
105  distance
106  = -std::sqrt( (track_X-clusterX)*(track_X-clusterX)
107  +(track_Y-clusterY)*(track_Y-clusterY)
108  +(track_Z-clusterZ)*(track_Z-clusterZ)
109  );
110 
111  }
112  break;
113  case PFLayer::PS1:
114  case PFLayer::PS2:
115  //Note Alex: Nothing implemented for the
116  //PreShower (No resolution maps yet)
117 #ifdef PFLOW_DEBUG
118  if( debug )
119  std::cout << "No link by rechit possible for pre-shower yet!"
120  << std::endl;
121 #endif
122  return -1.;
123  default:
124  return -1.;
125  }
126 
127 
128  // Check that, if the cluster is in the endcap,
129  // 0) the track indeed points to the endcap at vertex (DISABLED)
130  // 1) the track extrapolation is in the endcap too !
131  // 2) the track is in the same end-cap !
132  // PJ - 10-May-09
133  if ( !barrel ) {
134  // if ( fabs(trackEta) < 1.0 ) return -1;
135  if ( !hcal && fabs(track_Z) < 300. ) return -1.;
136  if ( track_Z * clusterZ < 0. ) return -1.;
137  }
138  // Check that, if the cluster is in the barrel,
139  // 1) the track is in the barrel too !
140  if ( barrel )
141  if ( !hcal && fabs(track_Z) > 300. ) return -1.;
142 
143  // Finally check that, if the track points to the central barrel (|eta| < 1),
144  // it cannot be linked to a cluster in Endcaps (avoid low pt loopers)
145 
146 
147  double dist = LinkByRecHit::computeDist( clustereta, clusterphi,
148  tracketa, trackphi);
149 
150 #ifdef PFLOW_DEBUG
151  if(debug) std::cout<<"test link by rechit "<< dist <<" "<<std::endl;
152  if(debug){
153  std::cout<<" clustereta " << clustereta
154  <<" clusterphi " << clusterphi
155  <<" tracketa " << tracketa
156  <<" trackphi " << trackphi << std::endl;
157  }
158 #endif
159 
160  //Testing if Track can be linked by rechit to a cluster.
161  //A cluster can be linked to a track if the extrapolated position
162  //of the track to the ECAL ShowerMax/HCAL entrance falls within
163  //the boundaries of any cell that belongs to this cluster.
164 
165  const std::vector< reco::PFRecHitFraction >&
166  fracs = cluster.recHitFractions();
167 
168  bool linkedbyrechit = false;
169  //loop rechits
170  for(unsigned int rhit = 0; rhit < fracs.size(); ++rhit){
171 
172  const reco::PFRecHitRef& rh = fracs[rhit].recHitRef();
173  double fraction = fracs[rhit].fraction();
174  if(fraction < 1E-4) continue;
175  if(rh.isNull()) continue;
176 
177  //getting rechit center position
178  const reco::PFRecHit& rechit_cluster = *rh;
179  const math::XYZPoint& posxyz
180  = rechit_cluster.position();
181  const reco::PFRecHit::REPPoint& posrep
182  = rechit_cluster.positionREP();
183 
184  //getting rechit corners
185  const std::vector< math::XYZPoint >&
186  cornersxyz = rechit_cluster.getCornersXYZ();
187  const std::vector<reco::PFRecHit::REPPoint>& corners =
188  rechit_cluster.getCornersREP();
189  assert(corners.size() == 4);
190 
191  if( barrel || hcal ){ // barrel case matching in eta/phi
192  // (and HCAL endcap too!)
193 
194  //rechit size determination
195  // blown up by 50% (HCAL) to 100% (ECAL) to include cracks & gaps
196  // also blown up to account for multiple scattering at low pt.
197  double rhsizeEta
198  = fabs(corners[0].Eta() - corners[2].Eta());
199  double rhsizePhi
200  = fabs(corners[0].Phi() - corners[2].Phi());
201  if ( rhsizePhi > M_PI ) rhsizePhi = 2.*M_PI - rhsizePhi;
202  if ( hcal ) {
203  rhsizeEta = rhsizeEta * (1.50 + 0.5/fracs.size()) + 0.2*fabs(dHEta);
204  rhsizePhi = rhsizePhi * (1.50 + 0.5/fracs.size()) + 0.2*fabs(dHPhi);
205 
206  } else {
207  rhsizeEta *= 2.00 + 1.0/fracs.size()/std::min(1.,trackPt/2.);
208  rhsizePhi *= 2.00 + 1.0/fracs.size()/std::min(1.,trackPt/2.);
209  }
210 
211 #ifdef PFLOW_DEBUG
212  if( debug ) {
213  std::cout << rhit << " Hcal RecHit="
214  << posrep.Eta() << " "
215  << posrep.Phi() << " "
216  << rechit_cluster.energy()
217  << std::endl;
218  for ( unsigned jc=0; jc<4; ++jc )
219  std::cout<<"corners "<<jc<<" "<<corners[jc].Eta()
220  <<" "<<corners[jc].Phi()<<std::endl;
221 
222  std::cout << "RecHit SizeEta=" << rhsizeEta
223  << " SizePhi=" << rhsizePhi << std::endl;
224  }
225 #endif
226 
227  //distance track-rechit center
228  // const math::XYZPoint& posxyz
229  // = rechit_cluster.position();
230  double deta = fabs(posrep.Eta() - tracketa);
231  double dphi = fabs(posrep.Phi() - trackphi);
232  if ( dphi > M_PI ) dphi = 2.*M_PI - dphi;
233 
234 #ifdef PFLOW_DEBUG
235  if( debug ){
236  std::cout << "distance="
237  << deta << " "
238  << dphi << " ";
239  if(deta < (rhsizeEta/2.) && dphi < (rhsizePhi/2.))
240  std::cout << " link here !" << std::endl;
241  else std::cout << std::endl;
242  }
243 #endif
244 
245  if(deta < (rhsizeEta/2.) && dphi < (rhsizePhi/2.)){
246  linkedbyrechit = true;
247  break;
248  }
249  }
250  else { //ECAL & PS endcap case, matching in X,Y
251 
252 #ifdef PFLOW_DEBUG
253  if( debug ){
254  const math::XYZPoint& posxyz
255  = rechit_cluster.position();
256 
257  std::cout << "RH " << posxyz.X()
258  << " " << posxyz.Y()
259  << std::endl;
260 
261  std::cout << "TRACK " << track_X
262  << " " << track_Y
263  << std::endl;
264  }
265 #endif
266 
267  double x[5];
268  double y[5];
269 
270  for ( unsigned jc=0; jc<4; ++jc ) {
271  math::XYZPoint cornerposxyz = cornersxyz[jc];
272  x[jc] = cornerposxyz.X() + (cornerposxyz.X()-posxyz.X())
273  * (1.00+0.50/fracs.size()/std::min(1.,trackPt/2.));
274  y[jc] = cornerposxyz.Y() + (cornerposxyz.Y()-posxyz.Y())
275  * (1.00+0.50/fracs.size()/std::min(1.,trackPt/2.));
276 
277 #ifdef PFLOW_DEBUG
278  if( debug ){
279  std::cout<<"corners "<<jc
280  << " " << cornerposxyz.X()
281  << " " << cornerposxyz.Y()
282  << std::endl;
283  }
284 #endif
285  }//loop corners
286 
287  //need to close the polygon in order to
288  //use the TMath::IsInside fonction from root lib
289  x[4] = x[0];
290  y[4] = y[0];
291 
292  //Check if the extrapolation point of the track falls
293  //within the rechit boundaries
294  bool isinside = TMath::IsInside(track_X,
295  track_Y,
296  5,x,y);
297 
298  if( isinside ){
299  linkedbyrechit = true;
300  break;
301  }
302  }//
303 
304  }//loop rechits
305 
306  if( linkedbyrechit ) {
307 #ifdef PFLOW_DEBUG
308  if( debug )
309  std::cout << "Track and Cluster LINKED BY RECHIT" << std::endl;
310 #endif
311  /*
312  //if ( distance > 40. || distance < -100. )
313  double clusterr = std::sqrt(clusterX*clusterX+clusterY*clusterY);
314  double trackr = std::sqrt(track_X*track_X+track_Y*track_Y);
315  if ( distance > 40. )
316  std::cout << "Distance = " << distance
317  << ", Barrel/Hcal/Brem ? " << barrel << " " << hcal << " " << isBrem << std::endl
318  << " Cluster " << clusterr << " " << clusterZ << " " << clusterphi << " " << clustereta << std::endl
319  << " Track " << trackr << " " << track_Z << " " << trackphi << " " << tracketa << std::endl;
320  if ( !barrel && fabs(trackEta) < 1.0 ) {
321  double clusterr = std::sqrt(clusterX*clusterX+clusterY*clusterY);
322  double trackr = std::sqrt(track_X*track_X+track_Y*track_Y);
323  std::cout << "TrackEta/Pt = " << trackEta << " " << trackPt << ", distance = " << distance << std::endl
324  << ", Barrel/Hcal/Brem ? " << barrel << " " << hcal << " " << isBrem << std::endl
325  << " Cluster " << clusterr << " " << clusterZ << " " << clusterphi << " " << clustereta << std::endl
326  << " Track " << trackr << " " << track_Z << " " << trackphi << " " << tracketa << " " << trackEta << " " << trackPt << std::endl;
327  }
328  */
329  return dist;
330  } else {
331  return -1.;
332  }
333 
334 }
335 
336 
337 
338 double
340  const reco::PFCluster& clusterPS,
341  bool debug) {
342 
343 // 0.19 <-> strip_pitch
344 // 6.1 <-> strip_length
345  static double resPSpitch = 0.19/sqrt(12.);
346  static double resPSlength = 6.1/sqrt(12.);
347 
348  // Check that clusterECAL is in ECAL endcap and that clusterPS is a preshower cluster
349  if ( clusterECAL.layer() != PFLayer::ECAL_ENDCAP ||
350  ( clusterPS.layer() != PFLayer::PS1 &&
351  clusterPS.layer() != PFLayer::PS2 ) ) return -1.;
352 
353 #ifdef PFLOW_DEBUG
354  if( debug )
355  std::cout<<"entering test link by rechit function for ECAL and PS"<<std::endl;
356 #endif
357 
358  //ECAL cluster position
359  double zECAL = clusterECAL.position().Z();
360  double xECAL = clusterECAL.position().X();
361  double yECAL = clusterECAL.position().Y();
362 
363  // PS cluster position, extrapolated to ECAL
364  double zPS = clusterPS.position().Z();
365  double xPS = clusterPS.position().X(); //* zECAL/zPS;
366  double yPS = clusterPS.position().Y(); //* zECAL/zPS;
367 // MDN jan09 : check that zEcal and zPs have the same sign
368  if (zECAL*zPS <0.) return -1.;
369  double deltaX = 0.;
370  double deltaY = 0.;
371  double sqr12 = std::sqrt(12.);
372  switch (clusterPS.layer()) {
373  case PFLayer::PS1:
374  // vertical strips, measure x with pitch precision
375  deltaX = resPSpitch * sqr12;
376  deltaY = resPSlength * sqr12;
377  break;
378  case PFLayer::PS2:
379  // horizontal strips, measure y with pitch precision
380  deltaY = resPSpitch * sqr12;
381  deltaX = resPSlength * sqr12;
382  break;
383  default:
384  break;
385  }
386 
387  // Get the rechits
388  const std::vector< reco::PFRecHitFraction >& fracs = clusterECAL.recHitFractions();
389  bool linkedbyrechit = false;
390  //loop rechits
391  for(unsigned int rhit = 0; rhit < fracs.size(); ++rhit){
392 
393  const reco::PFRecHitRef& rh = fracs[rhit].recHitRef();
394  double fraction = fracs[rhit].fraction();
395  if(fraction < 1E-4) continue;
396  if(rh.isNull()) continue;
397 
398  //getting rechit center position
399  const reco::PFRecHit& rechit_cluster = *rh;
400 
401  //getting rechit corners
402  const std::vector< math::XYZPoint >& corners = rechit_cluster.getCornersXYZ();
403  assert(corners.size() == 4);
404 
405  const math::XYZPoint& posxyz = rechit_cluster.position() * zPS/zECAL;
406 #ifdef PFLOW_DEBUG
407  if( debug ){
408  std::cout << "Ecal rechit " << posxyz.X() << " " << posxyz.Y() << std::endl;
409  std::cout << "PS cluster " << xPS << " " << yPS << std::endl;
410  }
411 #endif
412 
413  double x[5];
414  double y[5];
415  for ( unsigned jc=0; jc<4; ++jc ) {
416  // corner position projected onto the preshower
417  math::XYZPoint cornerpos = corners[jc] * zPS/zECAL;
418  // Inflate the size by the size of the PS strips, and by 5% to include ECAL cracks.
419  x[jc] = cornerpos.X() + (cornerpos.X()-posxyz.X()) * (0.05 +1.0/fabs((cornerpos.X()-posxyz.X()))*deltaX/2.);
420  y[jc] = cornerpos.Y() + (cornerpos.Y()-posxyz.Y()) * (0.05 +1.0/fabs((cornerpos.Y()-posxyz.Y()))*deltaY/2.);
421 
422 #ifdef PFLOW_DEBUG
423  if( debug ){
424  std::cout<<"corners "<<jc
425  << " " << cornerpos.X() << " " << x[jc]
426  << " " << cornerpos.Y() << " " << y[jc]
427  << std::endl;
428  }
429 #endif
430  }//loop corners
431 
432  //need to close the polygon in order to
433  //use the TMath::IsInside fonction from root lib
434  x[4] = x[0];
435  y[4] = y[0];
436 
437  //Check if the extrapolation point of the track falls
438  //within the rechit boundaries
439  bool isinside = TMath::IsInside(xPS,yPS,5,x,y);
440 
441  if( isinside ){
442  linkedbyrechit = true;
443  break;
444  }
445 
446  }//loop rechits
447 
448  if( linkedbyrechit ) {
449  if( debug ) std::cout << "Cluster PS and Cluster ECAL LINKED BY RECHIT" << std::endl;
450  double dist = computeDist( xECAL/1000.,yECAL/1000.,
451  xPS/1000. ,yPS/1000);
452  return dist;
453  } else {
454  return -1.;
455  }
456 
457 }
458 
459 
460 
461 double
463  const reco::PFCluster& clusterHFHAD,
464  bool debug) {
465 
466  math::XYZPoint posxyzEM = clusterHFEM.position();
467  math::XYZPoint posxyzHAD = clusterHFHAD.position();
468 
469  double dX = posxyzEM.X()-posxyzHAD.X();
470  double dY = posxyzEM.Y()-posxyzHAD.Y();
471  double sameZ = posxyzEM.Z()*posxyzHAD.Z();
472 
473  if(sameZ<0) return -1.;
474 
475  double dist2 = dX*dX + dY*dY;
476 
477  if( dist2 < 0.1 ) {
478  // less than one mm
479  double dist = sqrt( dist2 );
480  return dist;;
481  }
482  else
483  return -1.;
484 
485 }
486 
487 double
488 LinkByRecHit::computeDist( double eta1, double phi1,
489  double eta2, double phi2 ) {
490 
491  double phicor = normalizedPhi(phi1 - phi2);
492 
493  // double chi2 =
494  // (eta1 - eta2)*(eta1 - eta2) / ( reta1*reta1+ reta2*reta2 ) +
495  // phicor*phicor / ( rphi1*rphi1+ rphi2*rphi2 );
496 
497  double dist = std::sqrt( (eta1 - eta2)*(eta1 - eta2)
498  + phicor*phicor);
499 
500  return dist;
501 
502 }
PFLayer::Layer layer() const
cluster layer, see PFLayer.h in this directory
Definition: PFCluster.cc:81
const REPPoint & positionREP() const
trajectory position in (rho, eta, phi) base
const math::XYZPoint & position() const
cluster centroid position
Definition: CaloCluster.h:123
reconstructed track used as an input to particle flow
Definition: PFRecTrack.h:22
ROOT::Math::PositionVector3D< ROOT::Math::CylindricalEta3D< Double32_t > > REPPoint
Definition: PFRecHit.h:39
Particle flow cluster, see clustering algorithm in PFClusterAlgo.
Definition: PFCluster.h:42
static double testECALAndPSByRecHit(const reco::PFCluster &clusterECAL, const reco::PFCluster &clusterPS, bool debug=false)
const math::XYZPoint & position() const
cartesian position (x, y, z)
#define min(a, b)
Definition: mlp_lapack.h:161
static double computeDist(double eta1, double phi1, double eta2, double phi2)
computes a chisquare
const std::vector< math::XYZPoint > & getCornersXYZ() const
rechit corners
Definition: PFRecHit.h:144
bool isNull() const
Checks for null.
Definition: Ref.h:246
Particle flow rechit (rechit + geometry and topology information). See clustering algorithm in PFClus...
Definition: PFRecHit.h:31
Point of closest approach from beam axis (initial point in the case of PFSimParticle) ...
T sqrt(T t)
Definition: SSEVec.h:28
const REPPoint & positionREP() const
cluster position: rho, eta, phi
Definition: PFCluster.h:76
const reco::PFTrajectoryPoint & extrapolatedPoint(unsigned layerid) const
Definition: PFTrack.cc:76
const math::XYZPoint & position() const
is seed ?
Definition: PFRecHit.h:135
const REPPointVector & getCornersREP() const
rechit corners
Definition: PFRecHit.h:148
const math::XYZTLorentzVector & momentum() const
4-momenta quadrivector
#define M_PI
Definition: BFit3D.cc:3
bool isValid() const
is this point valid ?
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:13
static double testHFEMAndHFHADByRecHit(const reco::PFCluster &clusterHFEM, const reco::PFCluster &clusterHFHAD, bool debug=false)
test association between HFEM and HFHAD, by rechit
double energy() const
rechit energy
Definition: PFRecHit.h:105
const std::vector< reco::PFRecHitFraction > & recHitFractions() const
vector of rechit fractions
Definition: PFCluster.h:63
tuple cout
Definition: gather_cfg.py:41
A PFTrack holds several trajectory points, which basically contain the position and momentum of a tra...
Definition: DDAxes.h:10
#define debug
Definition: MEtoEDMFormat.h:34
double normalizedPhi(double phi)
Definition: normalizedPhi.cc:5
static double testTrackAndClusterByRecHit(const reco::PFRecTrack &track, const reco::PFCluster &cluster, bool isBrem=false, bool debug=false)
Definition: LinkByRecHit.cc:10
const REPPoint & positionREP() const
rechit cell centre rho, eta, phi. call calculatePositionREP before !
Definition: PFRecHit.cc:121