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  double horesolscale=1.0;
31 
32  //track extrapolation
33  const reco::PFTrajectoryPoint& atVertex
35  const reco::PFTrajectoryPoint& atECAL
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  switch (cluster.layer()) {
52  case PFLayer::ECAL_BARREL: barrel = true;
54 #ifdef PFLOW_DEBUG
55  if( debug )
56  std::cout << "Fetching Ecal Resolution Maps"
57  << std::endl;
58 #endif
59  // did not reach ecal, cannot be associated with a cluster.
60  if( ! atECAL.isValid() ) return -1.;
61 
62  tracketa = atECAL.positionREP().Eta();
63  trackphi = atECAL.positionREP().Phi();
64  track_X = atECAL.position().X();
65  track_Y = atECAL.position().Y();
66  track_Z = atECAL.position().Z();
67 
68 /* distance
69  = std::sqrt( (track_X-clusterX)*(track_X-clusterX)
70  +(track_Y-clusterY)*(track_Y-clusterY)
71  +(track_Z-clusterZ)*(track_Z-clusterZ)
72  );
73 */
74  break;
75 
76  case PFLayer::HCAL_BARREL1: barrel = true;
78 #ifdef PFLOW_DEBUG
79  if( debug )
80  std::cout << "Fetching Hcal Resolution Maps"
81  << std::endl;
82 #endif
83  if( isBrem ) {
84  return -1.;
85  } else {
86  hcal=true;
87  const reco::PFTrajectoryPoint& atHCAL
89  const reco::PFTrajectoryPoint& atHCALExit
91  // did not reach hcal, cannot be associated with a cluster.
92  if( ! atHCAL.isValid() ) return -1.;
93 
94  // The link is computed between 0 and ~1 interaction length in HCAL
95  dHEta = atHCALExit.positionREP().Eta()-atHCAL.positionREP().Eta();
96  dHPhi = atHCALExit.positionREP().Phi()-atHCAL.positionREP().Phi();
97  if ( dHPhi > M_PI ) dHPhi = dHPhi - 2.*M_PI;
98  else if ( dHPhi < -M_PI ) dHPhi = dHPhi + 2.*M_PI;
99  tracketa = atHCAL.positionREP().Eta() + 0.1*dHEta;
100  trackphi = atHCAL.positionREP().Phi() + 0.1*dHPhi;
101  track_X = atHCAL.position().X();
102  track_Y = atHCAL.position().Y();
103  track_Z = atHCAL.position().Z();
104 /* distance
105  = -std::sqrt( (track_X-clusterX)*(track_X-clusterX)
106  +(track_Y-clusterY)*(track_Y-clusterY)
107  +(track_Z-clusterZ)*(track_Z-clusterZ)
108  );
109 */
110  }
111  break;
112 
113  case PFLayer::HCAL_BARREL2: barrel = true;
114 #ifdef PFLOW_DEBUG
115  if( debug )
116  std::cout << "Fetching HO Resolution Maps"
117  << std::endl;
118 #endif
119  if( isBrem ) {
120  return -1.;
121  } else {
122  hcal=true;
123  horesolscale=1.15;
124  const reco::PFTrajectoryPoint& atHO
126  // did not reach ho, cannot be associated with a cluster.
127  if( ! atHO.isValid() ) return -1.;
128 
129  //
130  tracketa = atHO.positionREP().Eta();
131  trackphi = atHO.positionREP().Phi();
132  track_X = atHO.position().X();
133  track_Y = atHO.position().Y();
134  track_Z = atHO.position().Z();
135 
136  // Is this check really useful ?
137  if ( fabs(track_Z) > 700.25 ) return -1.;
138 
139 
140 /* distance
141  = std::sqrt( (track_X-clusterX)*(track_X-clusterX)
142  +(track_Y-clusterY)*(track_Y-clusterY)
143  +(track_Z-clusterZ)*(track_Z-clusterZ)
144  );
145 */
146 #ifdef PFLOW_DEBUG
147 /* if( debug ) {
148  std::cout <<"dist "<<distance<<" "<<cluster.energy()<<" "<<cluster.layer()<<" "
149  <<track_X<<" "<<clusterX<<" "
150  <<track_Y<<" "<<clusterY<<" "
151  <<track_Z<<" "<<clusterZ<<std::endl;
152  }
153 */
154 #endif
155  }
156  break;
157 
158  case PFLayer::PS1:
159  case PFLayer::PS2:
160  //Note Alex: Nothing implemented for the
161  //PreShower (No resolution maps yet)
162 #ifdef PFLOW_DEBUG
163  if( debug )
164  std::cout << "No link by rechit possible for pre-shower yet!"
165  << std::endl;
166 #endif
167  return -1.;
168  default:
169  return -1.;
170  }
171 
172 
173  // Check that, if the cluster is in the endcap,
174  // 0) the track indeed points to the endcap at vertex (DISABLED)
175  // 1) the track extrapolation is in the endcap too !
176  // 2) the track is in the same end-cap !
177  // PJ - 10-May-09
178  if ( !barrel ) {
179  // if ( fabs(trackEta) < 1.0 ) return -1;
180  if ( !hcal && fabs(track_Z) < 300. ) return -1.;
181  if ( track_Z * clusterZ < 0. ) return -1.;
182  }
183  // Check that, if the cluster is in the barrel,
184  // 1) the track is in the barrel too !
185  if ( barrel ) {
186  if ( !hcal && fabs(track_Z) > 300. ) return -1.;
187  }
188 
189  double dist = LinkByRecHit::computeDist( clustereta, clusterphi,
190  tracketa, trackphi);
191 
192 #ifdef PFLOW_DEBUG
193  if(debug) std::cout<<"test link by rechit "<< dist <<" "<<std::endl;
194  if(debug){
195  std::cout<<" clustereta " << clustereta
196  <<" clusterphi " << clusterphi
197  <<" tracketa " << tracketa
198  <<" trackphi " << trackphi << std::endl;
199  }
200 #endif
201 
202  //Testing if Track can be linked by rechit to a cluster.
203  //A cluster can be linked to a track if the extrapolated position
204  //of the track to the ECAL ShowerMax/HCAL entrance falls within
205  //the boundaries of any cell that belongs to this cluster.
206 
207  const std::vector< reco::PFRecHitFraction >&
208  fracs = cluster.recHitFractions();
209 
210  bool linkedbyrechit = false;
211  //loop rechits
212  for(unsigned int rhit = 0; rhit < fracs.size(); ++rhit){
213 
214  const reco::PFRecHitRef& rh = fracs[rhit].recHitRef();
215  double fraction = fracs[rhit].fraction();
216  if(fraction < 1E-4) continue;
217  if(rh.isNull()) continue;
218 
219  //getting rechit center position
220  const reco::PFRecHit& rechit_cluster = *rh;
221  const math::XYZPoint& posxyz
222  = rechit_cluster.position();
223  const reco::PFRecHit::REPPoint& posrep
224  = rechit_cluster.positionREP();
225 
226  //getting rechit corners
227  const std::vector< math::XYZPoint >&
228  cornersxyz = rechit_cluster.getCornersXYZ();
229  const std::vector<reco::PFRecHit::REPPoint>& corners =
230  rechit_cluster.getCornersREP();
231  assert(corners.size() == 4);
232 
233  if( barrel || hcal ){ // barrel case matching in eta/phi
234  // (and HCAL endcap too!)
235 
236  //rechit size determination
237  // blown up by 50% (HCAL) to 100% (ECAL) to include cracks & gaps
238  // also blown up to account for multiple scattering at low pt.
239  double rhsizeEta
240  = fabs(corners[0].Eta() - corners[2].Eta());
241  double rhsizePhi
242  = fabs(corners[0].Phi() - corners[2].Phi());
243  if ( rhsizePhi > M_PI ) rhsizePhi = 2.*M_PI - rhsizePhi;
244  if ( hcal ) {
245  rhsizeEta = rhsizeEta * horesolscale * (1.50 + 0.5/fracs.size()) + 0.2*fabs(dHEta);
246  rhsizePhi = rhsizePhi * horesolscale * (1.50 + 0.5/fracs.size()) + 0.2*fabs(dHPhi);
247 
248  } else {
249  rhsizeEta *= 2.00 + 1.0/fracs.size()/std::min(1.,trackPt/2.);
250  rhsizePhi *= 2.00 + 1.0/fracs.size()/std::min(1.,trackPt/2.);
251  }
252 
253 #ifdef PFLOW_DEBUG
254  if( debug ) {
255  std::cout << rhit << " Hcal RecHit="
256  << posrep.Eta() << " "
257  << posrep.Phi() << " "
258  << rechit_cluster.energy()
259  << std::endl;
260  for ( unsigned jc=0; jc<4; ++jc )
261  std::cout<<"corners "<<jc<<" "<<corners[jc].Eta()
262  <<" "<<corners[jc].Phi()<<std::endl;
263 
264  std::cout << "RecHit SizeEta=" << rhsizeEta
265  << " SizePhi=" << rhsizePhi << std::endl;
266  }
267 #endif
268 
269  //distance track-rechit center
270  // const math::XYZPoint& posxyz
271  // = rechit_cluster.position();
272  double deta = fabs(posrep.Eta() - tracketa);
273  double dphi = fabs(posrep.Phi() - trackphi);
274  if ( dphi > M_PI ) dphi = 2.*M_PI - dphi;
275 
276 #ifdef PFLOW_DEBUG
277  if( debug ){
278  std::cout << "distance="
279  << deta << " "
280  << dphi << " ";
281  if(deta < (rhsizeEta/2.) && dphi < (rhsizePhi/2.))
282  std::cout << " link here !" << std::endl;
283  else std::cout << std::endl;
284  }
285 #endif
286 
287  if(deta < (rhsizeEta/2.) && dphi < (rhsizePhi/2.)){
288  linkedbyrechit = true;
289  break;
290  }
291  }
292  else { //ECAL & PS endcap case, matching in X,Y
293 
294 #ifdef PFLOW_DEBUG
295  if( debug ){
296  const math::XYZPoint& posxyz
297  = rechit_cluster.position();
298 
299  std::cout << "RH " << posxyz.X()
300  << " " << posxyz.Y()
301  << std::endl;
302 
303  std::cout << "TRACK " << track_X
304  << " " << track_Y
305  << std::endl;
306  }
307 #endif
308 
309  double x[5];
310  double y[5];
311 
312  for ( unsigned jc=0; jc<4; ++jc ) {
313  math::XYZPoint cornerposxyz = cornersxyz[jc];
314  x[jc] = cornerposxyz.X() + (cornerposxyz.X()-posxyz.X())
315  * (1.00+0.50/fracs.size()/std::min(1.,trackPt/2.));
316  y[jc] = cornerposxyz.Y() + (cornerposxyz.Y()-posxyz.Y())
317  * (1.00+0.50/fracs.size()/std::min(1.,trackPt/2.));
318 
319 #ifdef PFLOW_DEBUG
320  if( debug ){
321  std::cout<<"corners "<<jc
322  << " " << cornerposxyz.X()
323  << " " << cornerposxyz.Y()
324  << std::endl;
325  }
326 #endif
327  }//loop corners
328 
329  //need to close the polygon in order to
330  //use the TMath::IsInside fonction from root lib
331  x[4] = x[0];
332  y[4] = y[0];
333 
334  //Check if the extrapolation point of the track falls
335  //within the rechit boundaries
336  bool isinside = TMath::IsInside(track_X,
337  track_Y,
338  5,x,y);
339 
340  if( isinside ){
341  linkedbyrechit = true;
342  break;
343  }
344  }//
345 
346  }//loop rechits
347 
348  if( linkedbyrechit ) {
349 #ifdef PFLOW_DEBUG
350  if( debug )
351  std::cout << "Track and Cluster LINKED BY RECHIT" << std::endl;
352 #endif
353  /*
354  //if ( distance > 40. || distance < -100. )
355  double clusterr = std::sqrt(clusterX*clusterX+clusterY*clusterY);
356  double trackr = std::sqrt(track_X*track_X+track_Y*track_Y);
357  if ( distance > 40. )
358  std::cout << "Distance = " << distance
359  << ", Barrel/Hcal/Brem ? " << barrel << " " << hcal << " " << isBrem << std::endl
360  << " Cluster " << clusterr << " " << clusterZ << " " << clusterphi << " " << clustereta << std::endl
361  << " Track " << trackr << " " << track_Z << " " << trackphi << " " << tracketa << std::endl;
362  if ( !barrel && fabs(trackEta) < 1.0 ) {
363  double clusterr = std::sqrt(clusterX*clusterX+clusterY*clusterY);
364  double trackr = std::sqrt(track_X*track_X+track_Y*track_Y);
365  std::cout << "TrackEta/Pt = " << trackEta << " " << trackPt << ", distance = " << distance << std::endl
366  << ", Barrel/Hcal/Brem ? " << barrel << " " << hcal << " " << isBrem << std::endl
367  << " Cluster " << clusterr << " " << clusterZ << " " << clusterphi << " " << clustereta << std::endl
368  << " Track " << trackr << " " << track_Z << " " << trackphi << " " << tracketa << " " << trackEta << " " << trackPt << std::endl;
369  }
370  */
371  return dist;
372  } else {
373  return -1.;
374  }
375 
376 }
377 
378 
379 
380 double
382  const reco::PFCluster& clusterPS,
383  bool debug) {
384 
385 // 0.19 <-> strip_pitch
386 // 6.1 <-> strip_length
387  static double resPSpitch = 0.19/sqrt(12.);
388  static double resPSlength = 6.1/sqrt(12.);
389 
390  // Check that clusterECAL is in ECAL endcap and that clusterPS is a preshower cluster
391  if ( clusterECAL.layer() != PFLayer::ECAL_ENDCAP ||
392  ( clusterPS.layer() != PFLayer::PS1 &&
393  clusterPS.layer() != PFLayer::PS2 ) ) return -1.;
394 
395 #ifdef PFLOW_DEBUG
396  if( debug )
397  std::cout<<"entering test link by rechit function for ECAL and PS"<<std::endl;
398 #endif
399 
400  //ECAL cluster position
401  double zECAL = clusterECAL.position().Z();
402  double xECAL = clusterECAL.position().X();
403  double yECAL = clusterECAL.position().Y();
404 
405  // PS cluster position, extrapolated to ECAL
406  double zPS = clusterPS.position().Z();
407  double xPS = clusterPS.position().X(); //* zECAL/zPS;
408  double yPS = clusterPS.position().Y(); //* zECAL/zPS;
409 // MDN jan09 : check that zEcal and zPs have the same sign
410  if (zECAL*zPS <0.) return -1.;
411  double deltaX = 0.;
412  double deltaY = 0.;
413  double sqr12 = std::sqrt(12.);
414  switch (clusterPS.layer()) {
415  case PFLayer::PS1:
416  // vertical strips, measure x with pitch precision
417  deltaX = resPSpitch * sqr12;
418  deltaY = resPSlength * sqr12;
419  break;
420  case PFLayer::PS2:
421  // horizontal strips, measure y with pitch precision
422  deltaY = resPSpitch * sqr12;
423  deltaX = resPSlength * sqr12;
424  break;
425  default:
426  break;
427  }
428 
429  // Get the rechits
430  const std::vector< reco::PFRecHitFraction >& fracs = clusterECAL.recHitFractions();
431  bool linkedbyrechit = false;
432  //loop rechits
433  for(unsigned int rhit = 0; rhit < fracs.size(); ++rhit){
434 
435  const reco::PFRecHitRef& rh = fracs[rhit].recHitRef();
436  double fraction = fracs[rhit].fraction();
437  if(fraction < 1E-4) continue;
438  if(rh.isNull()) continue;
439 
440  //getting rechit center position
441  const reco::PFRecHit& rechit_cluster = *rh;
442 
443  //getting rechit corners
444  const std::vector< math::XYZPoint >& corners = rechit_cluster.getCornersXYZ();
445  assert(corners.size() == 4);
446 
447  const math::XYZPoint& posxyz = rechit_cluster.position() * zPS/zECAL;
448 #ifdef PFLOW_DEBUG
449  if( debug ){
450  std::cout << "Ecal rechit " << posxyz.X() << " " << posxyz.Y() << std::endl;
451  std::cout << "PS cluster " << xPS << " " << yPS << std::endl;
452  }
453 #endif
454 
455  double x[5];
456  double y[5];
457  for ( unsigned jc=0; jc<4; ++jc ) {
458  // corner position projected onto the preshower
459  math::XYZPoint cornerpos = corners[jc] * zPS/zECAL;
460  // Inflate the size by the size of the PS strips, and by 5% to include ECAL cracks.
461  x[jc] = cornerpos.X() + (cornerpos.X()-posxyz.X()) * (0.05 +1.0/fabs((cornerpos.X()-posxyz.X()))*deltaX/2.);
462  y[jc] = cornerpos.Y() + (cornerpos.Y()-posxyz.Y()) * (0.05 +1.0/fabs((cornerpos.Y()-posxyz.Y()))*deltaY/2.);
463 
464 #ifdef PFLOW_DEBUG
465  if( debug ){
466  std::cout<<"corners "<<jc
467  << " " << cornerpos.X() << " " << x[jc]
468  << " " << cornerpos.Y() << " " << y[jc]
469  << std::endl;
470  }
471 #endif
472  }//loop corners
473 
474  //need to close the polygon in order to
475  //use the TMath::IsInside fonction from root lib
476  x[4] = x[0];
477  y[4] = y[0];
478 
479  //Check if the extrapolation point of the track falls
480  //within the rechit boundaries
481  bool isinside = TMath::IsInside(xPS,yPS,5,x,y);
482 
483  if( isinside ){
484  linkedbyrechit = true;
485  break;
486  }
487 
488  }//loop rechits
489 
490  if( linkedbyrechit ) {
491  if( debug ) std::cout << "Cluster PS and Cluster ECAL LINKED BY RECHIT" << std::endl;
492  double dist = computeDist( xECAL/1000.,yECAL/1000.,
493  xPS/1000. ,yPS/1000,
494  false);
495  return dist;
496  } else {
497  return -1.;
498  }
499 
500 }
501 
502 
503 
504 double
506  const reco::PFCluster& clusterHFHAD,
507  bool debug) {
508 
509  math::XYZPoint posxyzEM = clusterHFEM.position();
510  math::XYZPoint posxyzHAD = clusterHFHAD.position();
511 
512  double dX = posxyzEM.X()-posxyzHAD.X();
513  double dY = posxyzEM.Y()-posxyzHAD.Y();
514  double sameZ = posxyzEM.Z()*posxyzHAD.Z();
515 
516  if(sameZ<0) return -1.;
517 
518  double dist2 = dX*dX + dY*dY;
519 
520  if( dist2 < 0.1 ) {
521  // less than one mm
522  double dist = sqrt( dist2 );
523  return dist;;
524  }
525  else
526  return -1.;
527 
528 }
529 
530 double
531 LinkByRecHit::computeDist( double eta1, double phi1,
532  double eta2, double phi2,
533  bool etaPhi ) {
534 
535  double phicor = etaPhi ? normalizedPhi(phi1 - phi2) : phi1 - phi2;
536 
537  // double chi2 =
538  // (eta1 - eta2)*(eta1 - eta2) / ( reta1*reta1+ reta2*reta2 ) +
539  // phicor*phicor / ( rphi1*rphi1+ rphi2*rphi2 );
540 
541  double dist = std::sqrt( (eta1 - eta2)*(eta1 - eta2)
542  + phicor*phicor);
543 
544  return dist;
545 
546 }
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
static double computeDist(double eta1, double phi1, double eta2, double phi2, bool etaPhi=true)
computes a chisquare
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
const std::vector< math::XYZPoint > & getCornersXYZ() const
rechit corners
Definition: PFRecHit.h:144
bool isNull() const
Checks for null.
Definition: Ref.h:247
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:46
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:121
A PFTrack holds several trajectory points, which basically contain the position and momentum of a tra...
x
Definition: VDTMath.h:216
#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