CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_2_9/src/RecoParticleFlow/PFClusterTools/src/LinkByRecHit.cc

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 // to enable debugs
00007 //#define PFLOW_DEBUG
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   //cluster position
00021   double clustereta  = cluster.positionREP().Eta();
00022   double clusterphi  = cluster.positionREP().Phi();
00023   //double clusterX    = cluster.position().X();
00024   //double clusterY    = cluster.position().Y();
00025   double clusterZ    = cluster.position().Z();
00026 
00027   bool barrel = false;
00028   bool hcal = false;
00029   // double distance = 999999.9;
00030   double horesolscale=1.0;
00031 
00032   //track extrapolation
00033   const reco::PFTrajectoryPoint& atVertex 
00034     = track.extrapolatedPoint( reco::PFTrajectoryPoint::ClosestApproach );
00035   const reco::PFTrajectoryPoint& atECAL 
00036     = track.extrapolatedPoint( reco::PFTrajectoryPoint::ECALShowerMax );
00037 
00038   //track at calo's
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   // Quantities at vertex
00048   double trackPt = isBrem ? 999. : sqrt(atVertex.momentum().Vect().Perp2());
00049   // double trackEta = isBrem ? 999. : atVertex.momentum().Vect().Eta();
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     // did not reach ecal, cannot be associated with a cluster.
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 /*    distance 
00069       = std::sqrt( (track_X-clusterX)*(track_X-clusterX)
00070                   +(track_Y-clusterY)*(track_Y-clusterY)
00071                   +(track_Z-clusterZ)*(track_Z-clusterZ)
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       // did not reach hcal, cannot be associated with a cluster.
00092       if( ! atHCAL.isValid() ) return -1.;   
00093       
00094       // The link is computed between 0 and ~1 interaction length in HCAL
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 /*      distance 
00105         = -std::sqrt( (track_X-clusterX)*(track_X-clusterX)
00106                      +(track_Y-clusterY)*(track_Y-clusterY)
00107                      +(track_Z-clusterZ)*(track_Z-clusterZ)
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       // did not reach ho, cannot be associated with a cluster.
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       // Is this check really useful ?
00137       if ( fabs(track_Z) > 700.25 ) return -1.;
00138 
00139 
00140 /*      distance 
00141         = std::sqrt( (track_X-clusterX)*(track_X-clusterX)
00142                       +(track_Y-clusterY)*(track_Y-clusterY)
00143                       +(track_Z-clusterZ)*(track_Z-clusterZ)
00144                       );
00145 */
00146 #ifdef PFLOW_DEBUG
00147 /*      if( debug ) {
00148         std::cout <<"dist "<<distance<<" "<<cluster.energy()<<" "<<cluster.layer()<<" "
00149                   <<track_X<<" "<<clusterX<<" "
00150                   <<track_Y<<" "<<clusterY<<" "
00151                   <<track_Z<<" "<<clusterZ<<std::endl;
00152       }
00153 */
00154 #endif
00155     }
00156     break;
00157 
00158   case PFLayer::PS1:
00159   case PFLayer::PS2:
00160     //Note Alex: Nothing implemented for the
00161     //PreShower (No resolution maps yet)
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   // Check that, if the cluster is in the endcap, 
00174   // 0) the track indeed points to the endcap at vertex (DISABLED)
00175   // 1) the track extrapolation is in the endcap too !
00176   // 2) the track is in the same end-cap !
00177   // PJ - 10-May-09
00178   if ( !barrel ) { 
00179     // if ( fabs(trackEta) < 1.0 ) return -1; 
00180     if ( !hcal && fabs(track_Z) < 300. ) return -1.;
00181     if ( track_Z * clusterZ < 0. ) return -1.;
00182   }
00183   // Check that, if the cluster is in the barrel, 
00184   // 1) the track is in the barrel too !
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   //Testing if Track can be linked by rechit to a cluster.
00203   //A cluster can be linked to a track if the extrapolated position 
00204   //of the track to the ECAL ShowerMax/HCAL entrance falls within 
00205   //the boundaries of any cell that belongs to this cluster.
00206 
00207   const std::vector< reco::PFRecHitFraction >& 
00208     fracs = cluster.recHitFractions();
00209   
00210   bool linkedbyrechit = false;
00211   //loop rechits
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     //getting rechit center position
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     //getting rechit corners
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 ){ // barrel case matching in eta/phi 
00234                           // (and HCAL endcap too!)
00235       
00236       //rechit size determination 
00237       // blown up by 50% (HCAL) to 100% (ECAL) to include cracks & gaps
00238       // also blown up to account for multiple scattering at low pt.
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       //distance track-rechit center
00270       // const math::XYZPoint& posxyz 
00271       // = rechit_cluster.position();
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 { //ECAL & PS endcap case, matching in X,Y
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       }//loop corners
00328       
00329       //need to close the polygon in order to
00330       //use the TMath::IsInside fonction from root lib
00331       x[4] = x[0];
00332       y[4] = y[0];
00333       
00334       //Check if the extrapolation point of the track falls 
00335       //within the rechit boundaries
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   }//loop rechits
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     //if ( distance > 40. || distance < -100. ) 
00355     double clusterr = std::sqrt(clusterX*clusterX+clusterY*clusterY);
00356     double trackr = std::sqrt(track_X*track_X+track_Y*track_Y);
00357     if ( distance > 40. ) 
00358     std::cout << "Distance = " << distance 
00359     << ", Barrel/Hcal/Brem ? " << barrel << " " << hcal << " " << isBrem << std::endl
00360     << " Cluster " << clusterr << " " << clusterZ << " " << clusterphi << " " << clustereta << std::endl
00361     << " Track   " << trackr << " " << track_Z << " " << trackphi << " " << tracketa << std::endl;
00362     if ( !barrel && fabs(trackEta) < 1.0 ) { 
00363       double clusterr = std::sqrt(clusterX*clusterX+clusterY*clusterY);
00364       double trackr = std::sqrt(track_X*track_X+track_Y*track_Y);
00365       std::cout << "TrackEta/Pt = " << trackEta << " " << trackPt << ", distance = " << distance << std::endl 
00366                 << ", Barrel/Hcal/Brem ? " << barrel << " " << hcal << " " << isBrem << std::endl
00367                 << " Cluster " << clusterr << " " << clusterZ << " " << clusterphi << " " << clustereta << std::endl
00368                 << " Track   " << trackr << " " << track_Z << " " << trackphi << " " << tracketa << " " << trackEta << " " << trackPt << std::endl;
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 // 0.19 <-> strip_pitch
00386 // 6.1  <-> strip_length
00387   static double resPSpitch = 0.19/sqrt(12.);
00388   static double resPSlength = 6.1/sqrt(12.);
00389 
00390   // Check that clusterECAL is in ECAL endcap and that clusterPS is a preshower cluster
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   //ECAL cluster position
00401   double zECAL  = clusterECAL.position().Z();
00402   double xECAL  = clusterECAL.position().X();
00403   double yECAL  = clusterECAL.position().Y();
00404 
00405   // PS cluster position, extrapolated to ECAL
00406   double zPS = clusterPS.position().Z();
00407   double xPS = clusterPS.position().X(); //* zECAL/zPS;
00408   double yPS = clusterPS.position().Y(); //* zECAL/zPS;
00409 // MDN jan09 : check that zEcal and zPs have the same sign
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     // vertical strips, measure x with pitch precision
00417     deltaX = resPSpitch * sqr12;
00418     deltaY = resPSlength * sqr12;
00419     break;
00420   case PFLayer::PS2:
00421     // horizontal strips, measure y with pitch precision
00422     deltaY = resPSpitch * sqr12;
00423     deltaX = resPSlength * sqr12;
00424     break;
00425   default:
00426     break;
00427   }
00428 
00429   // Get the rechits
00430   const std::vector< reco::PFRecHitFraction >&  fracs = clusterECAL.recHitFractions();
00431   bool linkedbyrechit = false;
00432   //loop rechits
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     //getting rechit center position
00441     const reco::PFRecHit& rechit_cluster = *rh;
00442     
00443     //getting rechit corners
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       // corner position projected onto the preshower
00459       math::XYZPoint cornerpos = corners[jc] * zPS/zECAL;
00460       // Inflate the size by the size of the PS strips, and by 5% to include ECAL cracks.
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     }//loop corners
00473     
00474     //need to close the polygon in order to
00475     //use the TMath::IsInside fonction from root lib
00476     x[4] = x[0];
00477     y[4] = y[0];
00478     
00479     //Check if the extrapolation point of the track falls 
00480     //within the rechit boundaries
00481     bool isinside = TMath::IsInside(xPS,yPS,5,x,y);
00482       
00483     if( isinside ){
00484       linkedbyrechit = true;
00485       break;
00486     }
00487 
00488   }//loop rechits
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     // less than one mm
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   // double chi2 =  
00538   //  (eta1 - eta2)*(eta1 - eta2) / ( reta1*reta1+ reta2*reta2 ) +
00539   //  phicor*phicor / ( rphi1*rphi1+ rphi2*rphi2 );
00540 
00541   double dist = std::sqrt( (eta1 - eta2)*(eta1 - eta2) 
00542                           + phicor*phicor);
00543 
00544   return dist;
00545 
00546 }