Go to the documentation of this file.00001 #include "RecoParticleFlow/PFClusterTools/interface/LinkByRecHit.h"
00002 #include "DataFormats/Math/interface/normalizedPhi.h"
00003 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
00004 #include "TMath.h"
00005
00006
00007
00008
00009 double
00010 LinkByRecHit::testTrackAndClusterByRecHit ( const reco::PFRecTrack& track,
00011 const reco::PFCluster& cluster,
00012 bool isBrem,
00013 bool debug) {
00014
00015 #ifdef PFLOW_DEBUG
00016 if( debug )
00017 std::cout<<"entering test link by rechit function"<<std::endl;
00018 #endif
00019
00020
00021 double clustereta = cluster.positionREP().Eta();
00022 double clusterphi = cluster.positionREP().Phi();
00023
00024
00025 double clusterZ = cluster.position().Z();
00026
00027 bool barrel = false;
00028 bool hcal = false;
00029
00030 double horesolscale=1.0;
00031
00032
00033 const reco::PFTrajectoryPoint& atVertex
00034 = track.extrapolatedPoint( reco::PFTrajectoryPoint::ClosestApproach );
00035 const reco::PFTrajectoryPoint& atECAL
00036 = track.extrapolatedPoint( reco::PFTrajectoryPoint::ECALShowerMax );
00037
00038
00039 double tracketa = 999999.9;
00040 double trackphi = 999999.9;
00041 double track_X = 999999.9;
00042 double track_Y = 999999.9;
00043 double track_Z = 999999.9;
00044 double dHEta = 0.;
00045 double dHPhi = 0.;
00046
00047
00048 double trackPt = isBrem ? 999. : sqrt(atVertex.momentum().Vect().Perp2());
00049
00050
00051 switch (cluster.layer()) {
00052 case PFLayer::ECAL_BARREL: barrel = true;
00053 case PFLayer::ECAL_ENDCAP:
00054 #ifdef PFLOW_DEBUG
00055 if( debug )
00056 std::cout << "Fetching Ecal Resolution Maps"
00057 << std::endl;
00058 #endif
00059
00060 if( ! atECAL.isValid() ) return -1.;
00061
00062 tracketa = atECAL.positionREP().Eta();
00063 trackphi = atECAL.positionREP().Phi();
00064 track_X = atECAL.position().X();
00065 track_Y = atECAL.position().Y();
00066 track_Z = atECAL.position().Z();
00067
00068
00069
00070
00071
00072
00073
00074 break;
00075
00076 case PFLayer::HCAL_BARREL1: barrel = true;
00077 case PFLayer::HCAL_ENDCAP:
00078 #ifdef PFLOW_DEBUG
00079 if( debug )
00080 std::cout << "Fetching Hcal Resolution Maps"
00081 << std::endl;
00082 #endif
00083 if( isBrem ) {
00084 return -1.;
00085 } else {
00086 hcal=true;
00087 const reco::PFTrajectoryPoint& atHCAL
00088 = track.extrapolatedPoint( reco::PFTrajectoryPoint::HCALEntrance );
00089 const reco::PFTrajectoryPoint& atHCALExit
00090 = track.extrapolatedPoint( reco::PFTrajectoryPoint::HCALExit );
00091
00092 if( ! atHCAL.isValid() ) return -1.;
00093
00094
00095 dHEta = atHCALExit.positionREP().Eta()-atHCAL.positionREP().Eta();
00096 dHPhi = atHCALExit.positionREP().Phi()-atHCAL.positionREP().Phi();
00097 if ( dHPhi > M_PI ) dHPhi = dHPhi - 2.*M_PI;
00098 else if ( dHPhi < -M_PI ) dHPhi = dHPhi + 2.*M_PI;
00099 tracketa = atHCAL.positionREP().Eta() + 0.1*dHEta;
00100 trackphi = atHCAL.positionREP().Phi() + 0.1*dHPhi;
00101 track_X = atHCAL.position().X();
00102 track_Y = atHCAL.position().Y();
00103 track_Z = atHCAL.position().Z();
00104
00105
00106
00107
00108
00109
00110 }
00111 break;
00112
00113 case PFLayer::HCAL_BARREL2: barrel = true;
00114 #ifdef PFLOW_DEBUG
00115 if( debug )
00116 std::cout << "Fetching HO Resolution Maps"
00117 << std::endl;
00118 #endif
00119 if( isBrem ) {
00120 return -1.;
00121 } else {
00122 hcal=true;
00123 horesolscale=1.15;
00124 const reco::PFTrajectoryPoint& atHO
00125 = track.extrapolatedPoint( reco::PFTrajectoryPoint::HOLayer );
00126
00127 if( ! atHO.isValid() ) return -1.;
00128
00129
00130 tracketa = atHO.positionREP().Eta();
00131 trackphi = atHO.positionREP().Phi();
00132 track_X = atHO.position().X();
00133 track_Y = atHO.position().Y();
00134 track_Z = atHO.position().Z();
00135
00136
00137 if ( fabs(track_Z) > 700.25 ) return -1.;
00138
00139
00140
00141
00142
00143
00144
00145
00146 #ifdef PFLOW_DEBUG
00147
00148
00149
00150
00151
00152
00153
00154 #endif
00155 }
00156 break;
00157
00158 case PFLayer::PS1:
00159 case PFLayer::PS2:
00160
00161
00162 #ifdef PFLOW_DEBUG
00163 if( debug )
00164 std::cout << "No link by rechit possible for pre-shower yet!"
00165 << std::endl;
00166 #endif
00167 return -1.;
00168 default:
00169 return -1.;
00170 }
00171
00172
00173
00174
00175
00176
00177
00178 if ( !barrel ) {
00179
00180 if ( !hcal && fabs(track_Z) < 300. ) return -1.;
00181 if ( track_Z * clusterZ < 0. ) return -1.;
00182 }
00183
00184
00185 if ( barrel ) {
00186 if ( !hcal && fabs(track_Z) > 300. ) return -1.;
00187 }
00188
00189 double dist = LinkByRecHit::computeDist( clustereta, clusterphi,
00190 tracketa, trackphi);
00191
00192 #ifdef PFLOW_DEBUG
00193 if(debug) std::cout<<"test link by rechit "<< dist <<" "<<std::endl;
00194 if(debug){
00195 std::cout<<" clustereta " << clustereta
00196 <<" clusterphi " << clusterphi
00197 <<" tracketa " << tracketa
00198 <<" trackphi " << trackphi << std::endl;
00199 }
00200 #endif
00201
00202
00203
00204
00205
00206
00207 const std::vector< reco::PFRecHitFraction >&
00208 fracs = cluster.recHitFractions();
00209
00210 bool linkedbyrechit = false;
00211
00212 for(unsigned int rhit = 0; rhit < fracs.size(); ++rhit){
00213
00214 const reco::PFRecHitRef& rh = fracs[rhit].recHitRef();
00215 double fraction = fracs[rhit].fraction();
00216 if(fraction < 1E-4) continue;
00217 if(rh.isNull()) continue;
00218
00219
00220 const reco::PFRecHit& rechit_cluster = *rh;
00221 const math::XYZPoint& posxyz
00222 = rechit_cluster.position();
00223 const reco::PFRecHit::REPPoint& posrep
00224 = rechit_cluster.positionREP();
00225
00226
00227 const std::vector< math::XYZPoint >&
00228 cornersxyz = rechit_cluster.getCornersXYZ();
00229 const std::vector<reco::PFRecHit::REPPoint>& corners =
00230 rechit_cluster.getCornersREP();
00231 assert(corners.size() == 4);
00232
00233 if( barrel || hcal ){
00234
00235
00236
00237
00238
00239 double rhsizeEta
00240 = fabs(corners[0].Eta() - corners[2].Eta());
00241 double rhsizePhi
00242 = fabs(corners[0].Phi() - corners[2].Phi());
00243 if ( rhsizePhi > M_PI ) rhsizePhi = 2.*M_PI - rhsizePhi;
00244 if ( hcal ) {
00245 rhsizeEta = rhsizeEta * horesolscale * (1.50 + 0.5/fracs.size()) + 0.2*fabs(dHEta);
00246 rhsizePhi = rhsizePhi * horesolscale * (1.50 + 0.5/fracs.size()) + 0.2*fabs(dHPhi);
00247
00248 } else {
00249 rhsizeEta *= 2.00 + 1.0/fracs.size()/std::min(1.,trackPt/2.);
00250 rhsizePhi *= 2.00 + 1.0/fracs.size()/std::min(1.,trackPt/2.);
00251 }
00252
00253 #ifdef PFLOW_DEBUG
00254 if( debug ) {
00255 std::cout << rhit << " Hcal RecHit="
00256 << posrep.Eta() << " "
00257 << posrep.Phi() << " "
00258 << rechit_cluster.energy()
00259 << std::endl;
00260 for ( unsigned jc=0; jc<4; ++jc )
00261 std::cout<<"corners "<<jc<<" "<<corners[jc].Eta()
00262 <<" "<<corners[jc].Phi()<<std::endl;
00263
00264 std::cout << "RecHit SizeEta=" << rhsizeEta
00265 << " SizePhi=" << rhsizePhi << std::endl;
00266 }
00267 #endif
00268
00269
00270
00271
00272 double deta = fabs(posrep.Eta() - tracketa);
00273 double dphi = fabs(posrep.Phi() - trackphi);
00274 if ( dphi > M_PI ) dphi = 2.*M_PI - dphi;
00275
00276 #ifdef PFLOW_DEBUG
00277 if( debug ){
00278 std::cout << "distance="
00279 << deta << " "
00280 << dphi << " ";
00281 if(deta < (rhsizeEta/2.) && dphi < (rhsizePhi/2.))
00282 std::cout << " link here !" << std::endl;
00283 else std::cout << std::endl;
00284 }
00285 #endif
00286
00287 if(deta < (rhsizeEta/2.) && dphi < (rhsizePhi/2.)){
00288 linkedbyrechit = true;
00289 break;
00290 }
00291 }
00292 else {
00293
00294 #ifdef PFLOW_DEBUG
00295 if( debug ){
00296 const math::XYZPoint& posxyz
00297 = rechit_cluster.position();
00298
00299 std::cout << "RH " << posxyz.X()
00300 << " " << posxyz.Y()
00301 << std::endl;
00302
00303 std::cout << "TRACK " << track_X
00304 << " " << track_Y
00305 << std::endl;
00306 }
00307 #endif
00308
00309 double x[5];
00310 double y[5];
00311
00312 for ( unsigned jc=0; jc<4; ++jc ) {
00313 math::XYZPoint cornerposxyz = cornersxyz[jc];
00314 x[jc] = cornerposxyz.X() + (cornerposxyz.X()-posxyz.X())
00315 * (1.00+0.50/fracs.size()/std::min(1.,trackPt/2.));
00316 y[jc] = cornerposxyz.Y() + (cornerposxyz.Y()-posxyz.Y())
00317 * (1.00+0.50/fracs.size()/std::min(1.,trackPt/2.));
00318
00319 #ifdef PFLOW_DEBUG
00320 if( debug ){
00321 std::cout<<"corners "<<jc
00322 << " " << cornerposxyz.X()
00323 << " " << cornerposxyz.Y()
00324 << std::endl;
00325 }
00326 #endif
00327 }
00328
00329
00330
00331 x[4] = x[0];
00332 y[4] = y[0];
00333
00334
00335
00336 bool isinside = TMath::IsInside(track_X,
00337 track_Y,
00338 5,x,y);
00339
00340 if( isinside ){
00341 linkedbyrechit = true;
00342 break;
00343 }
00344 }
00345
00346 }
00347
00348 if( linkedbyrechit ) {
00349 #ifdef PFLOW_DEBUG
00350 if( debug )
00351 std::cout << "Track and Cluster LINKED BY RECHIT" << std::endl;
00352 #endif
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371 return dist;
00372 } else {
00373 return -1.;
00374 }
00375
00376 }
00377
00378
00379
00380 double
00381 LinkByRecHit::testECALAndPSByRecHit( const reco::PFCluster& clusterECAL,
00382 const reco::PFCluster& clusterPS,
00383 bool debug) {
00384
00385
00386
00387 static double resPSpitch = 0.19/sqrt(12.);
00388 static double resPSlength = 6.1/sqrt(12.);
00389
00390
00391 if ( clusterECAL.layer() != PFLayer::ECAL_ENDCAP ||
00392 ( clusterPS.layer() != PFLayer::PS1 &&
00393 clusterPS.layer() != PFLayer::PS2 ) ) return -1.;
00394
00395 #ifdef PFLOW_DEBUG
00396 if( debug )
00397 std::cout<<"entering test link by rechit function for ECAL and PS"<<std::endl;
00398 #endif
00399
00400
00401 double zECAL = clusterECAL.position().Z();
00402 double xECAL = clusterECAL.position().X();
00403 double yECAL = clusterECAL.position().Y();
00404
00405
00406 double zPS = clusterPS.position().Z();
00407 double xPS = clusterPS.position().X();
00408 double yPS = clusterPS.position().Y();
00409
00410 if (zECAL*zPS <0.) return -1.;
00411 double deltaX = 0.;
00412 double deltaY = 0.;
00413 double sqr12 = std::sqrt(12.);
00414 switch (clusterPS.layer()) {
00415 case PFLayer::PS1:
00416
00417 deltaX = resPSpitch * sqr12;
00418 deltaY = resPSlength * sqr12;
00419 break;
00420 case PFLayer::PS2:
00421
00422 deltaY = resPSpitch * sqr12;
00423 deltaX = resPSlength * sqr12;
00424 break;
00425 default:
00426 break;
00427 }
00428
00429
00430 const std::vector< reco::PFRecHitFraction >& fracs = clusterECAL.recHitFractions();
00431 bool linkedbyrechit = false;
00432
00433 for(unsigned int rhit = 0; rhit < fracs.size(); ++rhit){
00434
00435 const reco::PFRecHitRef& rh = fracs[rhit].recHitRef();
00436 double fraction = fracs[rhit].fraction();
00437 if(fraction < 1E-4) continue;
00438 if(rh.isNull()) continue;
00439
00440
00441 const reco::PFRecHit& rechit_cluster = *rh;
00442
00443
00444 const std::vector< math::XYZPoint >& corners = rechit_cluster.getCornersXYZ();
00445 assert(corners.size() == 4);
00446
00447 const math::XYZPoint& posxyz = rechit_cluster.position() * zPS/zECAL;
00448 #ifdef PFLOW_DEBUG
00449 if( debug ){
00450 std::cout << "Ecal rechit " << posxyz.X() << " " << posxyz.Y() << std::endl;
00451 std::cout << "PS cluster " << xPS << " " << yPS << std::endl;
00452 }
00453 #endif
00454
00455 double x[5];
00456 double y[5];
00457 for ( unsigned jc=0; jc<4; ++jc ) {
00458
00459 math::XYZPoint cornerpos = corners[jc] * zPS/zECAL;
00460
00461 x[jc] = cornerpos.X() + (cornerpos.X()-posxyz.X()) * (0.05 +1.0/fabs((cornerpos.X()-posxyz.X()))*deltaX/2.);
00462 y[jc] = cornerpos.Y() + (cornerpos.Y()-posxyz.Y()) * (0.05 +1.0/fabs((cornerpos.Y()-posxyz.Y()))*deltaY/2.);
00463
00464 #ifdef PFLOW_DEBUG
00465 if( debug ){
00466 std::cout<<"corners "<<jc
00467 << " " << cornerpos.X() << " " << x[jc]
00468 << " " << cornerpos.Y() << " " << y[jc]
00469 << std::endl;
00470 }
00471 #endif
00472 }
00473
00474
00475
00476 x[4] = x[0];
00477 y[4] = y[0];
00478
00479
00480
00481 bool isinside = TMath::IsInside(xPS,yPS,5,x,y);
00482
00483 if( isinside ){
00484 linkedbyrechit = true;
00485 break;
00486 }
00487
00488 }
00489
00490 if( linkedbyrechit ) {
00491 if( debug ) std::cout << "Cluster PS and Cluster ECAL LINKED BY RECHIT" << std::endl;
00492 double dist = computeDist( xECAL/1000.,yECAL/1000.,
00493 xPS/1000. ,yPS/1000,
00494 false);
00495 return dist;
00496 } else {
00497 return -1.;
00498 }
00499
00500 }
00501
00502
00503
00504 double
00505 LinkByRecHit::testHFEMAndHFHADByRecHit(const reco::PFCluster& clusterHFEM,
00506 const reco::PFCluster& clusterHFHAD,
00507 bool debug) {
00508
00509 math::XYZPoint posxyzEM = clusterHFEM.position();
00510 math::XYZPoint posxyzHAD = clusterHFHAD.position();
00511
00512 double dX = posxyzEM.X()-posxyzHAD.X();
00513 double dY = posxyzEM.Y()-posxyzHAD.Y();
00514 double sameZ = posxyzEM.Z()*posxyzHAD.Z();
00515
00516 if(sameZ<0) return -1.;
00517
00518 double dist2 = dX*dX + dY*dY;
00519
00520 if( dist2 < 0.1 ) {
00521
00522 double dist = sqrt( dist2 );
00523 return dist;;
00524 }
00525 else
00526 return -1.;
00527
00528 }
00529
00530 double
00531 LinkByRecHit::computeDist( double eta1, double phi1,
00532 double eta2, double phi2,
00533 bool etaPhi ) {
00534
00535 double phicor = etaPhi ? normalizedPhi(phi1 - phi2) : phi1 - phi2;
00536
00537
00538
00539
00540
00541 double dist = std::sqrt( (eta1 - eta2)*(eta1 - eta2)
00542 + phicor*phicor);
00543
00544 return dist;
00545
00546 }