Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:27:26 

0001 #include "RecoParticleFlow/PFClusterTools/interface/LinkByRecHit.h"
0002 #include "DataFormats/Math/interface/normalizedPhi.h"
0003 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
0004 #include "TMath.h"
0005 
0006 using BVector2D = Basic2DVector<double>;
0007 using Vector2D = Basic2DVector<double>::MathVector;
0008 namespace {
0009   const Vector2D one2D = BVector2D(1.0, 1.0).v;
0010   const Vector2D fivepercent2D = BVector2D(0.05, 0.05).v;
0011   // rechit with fraction below this value will be ignored in LinkByRecHit
0012   const float minPFRecHitFrac = 1E-4;
0013 }  // namespace
0014 
0015 // to enable debugs
0016 //#define PFLOW_DEBUG
0017 
0018 double LinkByRecHit::testTrackAndClusterByRecHit(const reco::PFRecTrack& track,
0019                                                  const reco::PFCluster& cluster,
0020                                                  bool isBrem,
0021                                                  bool debug) {
0022 #ifdef PFLOW_DEBUG
0023   if (debug)
0024     std::cout << "entering test link by rechit function" << std::endl;
0025 #endif
0026 
0027   //cluster position
0028   auto clustereta = cluster.positionREP().Eta();
0029   auto clusterphi = cluster.positionREP().Phi();
0030   auto clusterZ = cluster.position().Z();
0031 
0032   bool barrel = false;
0033   bool hcal = false;
0034   // double distance = 999999.9;
0035   double horesolscale = 1.0;
0036 
0037   //track extrapolation
0038   const reco::PFTrajectoryPoint& atVertex = track.extrapolatedPoint(reco::PFTrajectoryPoint::ClosestApproach);
0039   const reco::PFTrajectoryPoint& atECAL = track.extrapolatedPoint(reco::PFTrajectoryPoint::ECALShowerMax);
0040 
0041   //track at calo's
0042   double tracketa = 999999.9;
0043   double trackphi = 999999.9;
0044   double track_X = 999999.9;
0045   double track_Y = 999999.9;
0046   double track_Z = 999999.9;
0047   double dHEta = 0.;
0048   double dHPhi = 0.;
0049 
0050   // Quantities at vertex
0051   double trackPt = isBrem ? 999. : sqrt(atVertex.momentum().Vect().Perp2());
0052   // double trackEta = isBrem ? 999. : atVertex.momentum().Vect().Eta();
0053 
0054   switch (cluster.layer()) {
0055     case PFLayer::ECAL_BARREL:
0056       barrel = true;
0057       [[fallthrough]];
0058     case PFLayer::ECAL_ENDCAP:
0059 #ifdef PFLOW_DEBUG
0060       if (debug)
0061         std::cout << "Fetching Ecal Resolution Maps" << std::endl;
0062 #endif
0063       // did not reach ecal, cannot be associated with a cluster.
0064       if (!atECAL.isValid())
0065         return -1.;
0066 
0067       tracketa = atECAL.positionREP().Eta();
0068       trackphi = atECAL.positionREP().Phi();
0069       track_X = atECAL.position().X();
0070       track_Y = atECAL.position().Y();
0071       track_Z = atECAL.position().Z();
0072 
0073       /*    distance 
0074       = std::sqrt( (track_X-clusterX)*(track_X-clusterX)
0075           +(track_Y-clusterY)*(track_Y-clusterY)
0076           +(track_Z-clusterZ)*(track_Z-clusterZ)
0077            );
0078 */
0079       break;
0080 
0081     case PFLayer::HCAL_BARREL1:
0082       barrel = true;
0083       [[fallthrough]];
0084     case PFLayer::HCAL_ENDCAP:
0085 #ifdef PFLOW_DEBUG
0086       if (debug)
0087         std::cout << "Fetching Hcal Resolution Maps" << std::endl;
0088 #endif
0089       if (isBrem) {
0090         return -1.;
0091       } else {
0092         hcal = true;
0093         const reco::PFTrajectoryPoint& atHCAL = track.extrapolatedPoint(reco::PFTrajectoryPoint::HCALEntrance);
0094         const reco::PFTrajectoryPoint& atHCALExit = track.extrapolatedPoint(reco::PFTrajectoryPoint::HCALExit);
0095         // did not reach hcal, cannot be associated with a cluster.
0096         if (!atHCAL.isValid())
0097           return -1.;
0098 
0099         // The link is computed between 0 and ~1 interaction length in HCAL
0100         dHEta = atHCALExit.positionREP().Eta() - atHCAL.positionREP().Eta();
0101         dHPhi = atHCALExit.positionREP().Phi() - atHCAL.positionREP().Phi();
0102         if (dHPhi > M_PI)
0103           dHPhi = dHPhi - 2. * M_PI;
0104         else if (dHPhi < -M_PI)
0105           dHPhi = dHPhi + 2. * M_PI;
0106         tracketa = atHCAL.positionREP().Eta() + 0.1 * dHEta;
0107         trackphi = atHCAL.positionREP().Phi() + 0.1 * dHPhi;
0108         track_X = atHCAL.position().X();
0109         track_Y = atHCAL.position().Y();
0110         track_Z = atHCAL.position().Z();
0111         /*      distance 
0112     = -std::sqrt( (track_X-clusterX)*(track_X-clusterX)
0113              +(track_Y-clusterY)*(track_Y-clusterY)
0114              +(track_Z-clusterZ)*(track_Z-clusterZ)
0115              );
0116 */
0117       }
0118       break;
0119 
0120     case PFLayer::HCAL_BARREL2:
0121       barrel = true;
0122 #ifdef PFLOW_DEBUG
0123       if (debug)
0124         std::cout << "Fetching HO Resolution Maps" << std::endl;
0125 #endif
0126       if (isBrem) {
0127         return -1.;
0128       } else {
0129         hcal = true;
0130         horesolscale = 1.15;
0131         const reco::PFTrajectoryPoint& atHO = track.extrapolatedPoint(reco::PFTrajectoryPoint::HOLayer);
0132         // did not reach ho, cannot be associated with a cluster.
0133         if (!atHO.isValid())
0134           return -1.;
0135 
0136         //
0137         tracketa = atHO.positionREP().Eta();
0138         trackphi = atHO.positionREP().Phi();
0139         track_X = atHO.position().X();
0140         track_Y = atHO.position().Y();
0141         track_Z = atHO.position().Z();
0142 
0143         // Is this check really useful ?
0144         if (fabs(track_Z) > 700.25)
0145           return -1.;
0146 
0147 /*      distance 
0148     = std::sqrt( (track_X-clusterX)*(track_X-clusterX)
0149               +(track_Y-clusterY)*(track_Y-clusterY)
0150               +(track_Z-clusterZ)*(track_Z-clusterZ)
0151               );
0152 */
0153 #ifdef PFLOW_DEBUG
0154 /*      if( debug ) {
0155     std::cout <<"dist "<<distance<<" "<<cluster.energy()<<" "<<cluster.layer()<<" "
0156           <<track_X<<" "<<clusterX<<" "
0157           <<track_Y<<" "<<clusterY<<" "
0158           <<track_Z<<" "<<clusterZ<<std::endl;
0159       }
0160 */
0161 #endif
0162       }
0163       break;
0164 
0165     case PFLayer::PS1:
0166       [[fallthrough]];
0167     case PFLayer::PS2:
0168       //Note Alex: Nothing implemented for the
0169       //PreShower (No resolution maps yet)
0170 #ifdef PFLOW_DEBUG
0171       if (debug)
0172         std::cout << "No link by rechit possible for pre-shower yet!" << std::endl;
0173 #endif
0174       return -1.;
0175     default:
0176       return -1.;
0177   }
0178 
0179   // Check that, if the cluster is in the endcap,
0180   // 0) the track indeed points to the endcap at vertex (DISABLED)
0181   // 1) the track extrapolation is in the endcap too !
0182   // 2) the track is in the same end-cap !
0183   // PJ - 10-May-09
0184   if (!barrel) {
0185     // if ( fabs(trackEta) < 1.0 ) return -1;
0186     if (!hcal && fabs(track_Z) < 300.)
0187       return -1.;
0188     if (track_Z * clusterZ < 0.)
0189       return -1.;
0190   }
0191   // Check that, if the cluster is in the barrel,
0192   // 1) the track is in the barrel too !
0193   if (barrel) {
0194     if (!hcal && fabs(track_Z) > 300.)
0195       return -1.;
0196   }
0197 
0198   double dist = LinkByRecHit::computeDist(clustereta, clusterphi, tracketa, trackphi);
0199 
0200 #ifdef PFLOW_DEBUG
0201   if (debug)
0202     std::cout << "test link by rechit " << dist << " " << std::endl;
0203   if (debug) {
0204     std::cout << " clustereta " << clustereta << " clusterphi " << clusterphi << " tracketa " << tracketa
0205               << " trackphi " << trackphi << std::endl;
0206   }
0207 #endif
0208 
0209   //Testing if Track can be linked by rechit to a cluster.
0210   //A cluster can be linked to a track if the extrapolated position
0211   //of the track to the ECAL ShowerMax/HCAL entrance falls within
0212   //the boundaries of any cell that belongs to this cluster.
0213 
0214   const std::vector<reco::PFRecHitFraction>& fracs = cluster.recHitFractions();
0215 
0216   bool linkedbyrechit = false;
0217   //loop rechits
0218   for (unsigned int rhit = 0; rhit < fracs.size(); ++rhit) {
0219     const reco::PFRecHitRef& rh = fracs[rhit].recHitRef();
0220     double fraction = fracs[rhit].fraction();
0221     if (fraction < minPFRecHitFrac)
0222       continue;
0223     if (rh.isNull())
0224       continue;
0225 
0226     //getting rechit center position
0227     const auto& rechit_cluster = *rh;
0228     const auto& posxyz = rechit_cluster.position();
0229     const auto& posrep = rechit_cluster.positionREP();
0230 
0231     //getting rechit corners
0232     const auto& cornersxyz = rechit_cluster.getCornersXYZ();
0233     const auto& corners = rechit_cluster.getCornersREP();
0234 
0235     if (barrel || hcal) {  // barrel case matching in eta/phi
0236                            // (and HCAL endcap too!)
0237 
0238       //rechit size determination
0239       // blown up by 50% (HCAL) to 100% (ECAL) to include cracks & gaps
0240       // also blown up to account for multiple scattering at low pt.
0241       double rhsizeEta = std::abs(corners[3].eta() - corners[1].eta());
0242       double rhsizePhi = std::abs(corners[3].phi() - corners[1].phi());
0243       if (rhsizePhi > M_PI)
0244         rhsizePhi = 2. * M_PI - rhsizePhi;
0245       if (hcal) {
0246         const double mult = horesolscale * (1.50 + 0.5 / fracs.size());
0247         rhsizeEta = rhsizeEta * mult + 0.2 * std::abs(dHEta);
0248         rhsizePhi = rhsizePhi * mult + 0.2 * fabs(dHPhi);
0249 
0250       } else {
0251         const double mult = 2.00 + 1.0 / (fracs.size() * std::min(1., 0.5 * trackPt));
0252         rhsizeEta *= mult;
0253         rhsizePhi *= mult;
0254       }
0255 
0256 #ifdef PFLOW_DEBUG
0257       if (debug) {
0258         std::cout << rhit << " Hcal RecHit=" << posrep.Eta() << " " << posrep.Phi() << " " << rechit_cluster.energy()
0259                   << std::endl;
0260         for (unsigned jc = 0; jc < 4; ++jc)
0261           std::cout << "corners " << jc << " " << corners[jc].eta() << " " << corners[jc].phi() << std::endl;
0262 
0263         std::cout << "RecHit SizeEta=" << rhsizeEta << " SizePhi=" << rhsizePhi << std::endl;
0264       }
0265 #endif
0266 
0267       //distance track-rechit center
0268       // const math::XYZPoint& posxyz
0269       // = rechit_cluster.position();
0270       double deta = fabs(posrep.eta() - tracketa);
0271       double dphi = fabs(posrep.phi() - trackphi);
0272       if (dphi > M_PI)
0273         dphi = 2. * M_PI - dphi;
0274 
0275 #ifdef PFLOW_DEBUG
0276       if (debug) {
0277         std::cout << "distance=" << deta << " " << dphi << " ";
0278         if (deta < (0.5 * rhsizeEta) && dphi < (0.5 * rhsizePhi))
0279           std::cout << " link here !" << std::endl;
0280         else
0281           std::cout << std::endl;
0282       }
0283 #endif
0284 
0285       if (deta < (0.5 * rhsizeEta) && dphi < (0.5 * rhsizePhi)) {
0286         linkedbyrechit = true;
0287         break;
0288       }
0289     } else {  //ECAL & PS endcap case, matching in X,Y
0290 
0291 #ifdef PFLOW_DEBUG
0292       if (debug) {
0293         const auto& posxyz = rechit_cluster.position();
0294 
0295         std::cout << "RH " << posxyz.x() << " " << posxyz.y() << std::endl;
0296 
0297         std::cout << "TRACK " << track_X << " " << track_Y << std::endl;
0298       }
0299 #endif
0300 
0301       double x[5];
0302       double y[5];
0303 
0304       for (unsigned jc = 0; jc < 4; ++jc) {
0305         const auto& cornerposxyz = cornersxyz[jc];
0306         const double mult = (1.00 + 0.50 / (fracs.size() * std::min(1., 0.5 * trackPt)));
0307         x[3 - jc] = cornerposxyz.x() + (cornerposxyz.x() - posxyz.x()) * mult;
0308         y[3 - jc] = cornerposxyz.y() + (cornerposxyz.y() - posxyz.y()) * mult;
0309 
0310 #ifdef PFLOW_DEBUG
0311         if (debug) {
0312           std::cout << "corners " << jc << " " << cornerposxyz.x() << " " << cornerposxyz.y() << std::endl;
0313         }
0314 #endif
0315       }  //loop corners
0316 
0317       //need to close the polygon in order to
0318       //use the TMath::IsInside fonction from root lib
0319       x[4] = x[0];
0320       y[4] = y[0];
0321 
0322       //Check if the extrapolation point of the track falls
0323       //within the rechit boundaries
0324       bool isinside = TMath::IsInside(track_X, track_Y, 5, x, y);
0325 
0326       if (isinside) {
0327         linkedbyrechit = true;
0328         break;
0329       }
0330     }  //
0331 
0332   }  //loop rechits
0333 
0334   if (linkedbyrechit) {
0335 #ifdef PFLOW_DEBUG
0336     if (debug)
0337       std::cout << "Track and Cluster LINKED BY RECHIT" << std::endl;
0338 #endif
0339     /*    
0340     //if ( distance > 40. || distance < -100. ) 
0341     double clusterr = std::sqrt(clusterX*clusterX+clusterY*clusterY);
0342     double trackr = std::sqrt(track_X*track_X+track_Y*track_Y);
0343     if ( distance > 40. ) 
0344     std::cout << "Distance = " << distance 
0345     << ", Barrel/Hcal/Brem ? " << barrel << " " << hcal << " " << isBrem << std::endl
0346     << " Cluster " << clusterr << " " << clusterZ << " " << clusterphi << " " << clustereta << std::endl
0347     << " Track   " << trackr << " " << track_Z << " " << trackphi << " " << tracketa << std::endl;
0348     if ( !barrel && fabs(trackEta) < 1.0 ) { 
0349       double clusterr = std::sqrt(clusterX*clusterX+clusterY*clusterY);
0350       double trackr = std::sqrt(track_X*track_X+track_Y*track_Y);
0351       std::cout << "TrackEta/Pt = " << trackEta << " " << trackPt << ", distance = " << distance << std::endl 
0352         << ", Barrel/Hcal/Brem ? " << barrel << " " << hcal << " " << isBrem << std::endl
0353         << " Cluster " << clusterr << " " << clusterZ << " " << clusterphi << " " << clustereta << std::endl
0354         << " Track   " << trackr << " " << track_Z << " " << trackphi << " " << tracketa << " " << trackEta << " " << trackPt << std::endl;
0355     } 
0356     */
0357     return dist;
0358   } else {
0359     return -1.;
0360   }
0361 }
0362 
0363 double LinkByRecHit::testECALAndPSByRecHit(const reco::PFCluster& clusterECAL,
0364                                            const reco::PFCluster& clusterPS,
0365                                            bool debug) {
0366   // 0.19 <-> strip_pitch
0367   // 6.1  <-> strip_length
0368   static const double resPSpitch = 0.19 / sqrt(12.);
0369   static const double resPSlength = 6.1 / sqrt(12.);
0370 
0371   // Check that clusterECAL is in ECAL endcap and that clusterPS is a preshower cluster
0372   if (clusterECAL.layer() != PFLayer::ECAL_ENDCAP ||
0373       (clusterPS.layer() != PFLayer::PS1 && clusterPS.layer() != PFLayer::PS2))
0374     return -1.;
0375 
0376 #ifdef PFLOW_DEBUG
0377   if (debug)
0378     std::cout << "entering test link by rechit function for ECAL and PS" << std::endl;
0379 #endif
0380 
0381   //ECAL cluster position
0382   double zECAL = clusterECAL.position().Z();
0383   double xECAL = clusterECAL.position().X();
0384   double yECAL = clusterECAL.position().Y();
0385 
0386   // PS cluster position, extrapolated to ECAL
0387   double zPS = clusterPS.position().Z();
0388   double xPS = clusterPS.position().X();  //* zECAL/zPS;
0389   double yPS = clusterPS.position().Y();  //* zECAL/zPS;
0390                                           // MDN jan09 : check that zEcal and zPs have the same sign
0391   if (zECAL * zPS < 0.)
0392     return -1.;
0393   double deltaX = 0.;
0394   double deltaY = 0.;
0395   double sqr12 = std::sqrt(12.);
0396   switch (clusterPS.layer()) {
0397     case PFLayer::PS1:
0398       // vertical strips, measure x with pitch precision
0399       deltaX = resPSpitch * sqr12;
0400       deltaY = resPSlength * sqr12;
0401       break;
0402     case PFLayer::PS2:
0403       // horizontal strips, measure y with pitch precision
0404       deltaY = resPSpitch * sqr12;
0405       deltaX = resPSlength * sqr12;
0406       break;
0407     default:
0408       break;
0409   }
0410 
0411   auto deltaXY = BVector2D(deltaX, deltaY).v * 0.5;
0412   // Get the rechits
0413   auto zCorr = zPS / zECAL;
0414   const std::vector<reco::PFRecHitFraction>& fracs = clusterECAL.recHitFractions();
0415   bool linkedbyrechit = false;
0416   //loop rechits
0417   for (unsigned int rhit = 0; rhit < fracs.size(); ++rhit) {
0418     const auto& rh = fracs[rhit].recHitRef();
0419     double fraction = fracs[rhit].fraction();
0420     if (fraction < minPFRecHitFrac)
0421       continue;
0422     if (rh.isNull())
0423       continue;
0424 
0425     //getting rechit center position
0426     const reco::PFRecHit& rechit_cluster = *rh;
0427 
0428     //getting rechit corners
0429     auto const& corners = rechit_cluster.getCornersXYZ();
0430 
0431     auto posxy = BVector2D(rechit_cluster.position().xy()).v * zCorr;
0432 #ifdef PFLOW_DEBUG
0433     if (debug) {
0434       std::cout << "Ecal rechit " << posxy.x() << " " << posxy.y() << std::endl;
0435       std::cout << "PS cluster  " << xPS << " " << yPS << std::endl;
0436     }
0437 #endif
0438 
0439     double x[5];
0440     double y[5];
0441     for (unsigned jc = 0; jc < 4; ++jc) {
0442       // corner position projected onto the preshower
0443       Vector2D cornerpos = BVector2D(corners[jc].basicVector().xy()).v * zCorr;
0444       auto dist = (cornerpos - posxy);
0445       auto adist =
0446           BVector2D(std::abs(dist[0]), std::abs(dist[1])).v;  // all this beacuse icc does not support vector extension
0447       // Inflate the size by the size of the PS strips, and by 5% to include ECAL cracks.
0448       auto xy = cornerpos + (dist * (fivepercent2D + one2D / adist) * deltaXY);
0449       /*
0450       Vector2D xy(
0451           cornerpos.x() + (cornerpos.x()-posxy.x()) * (0.05 +1.0/std::abs((cornerpos.x()-posxy.x()))*deltaXY.x()),
0452           cornerpos.y() + (cornerpos.y()-posxy.y()) * (0.05 +1.0/std::abs((cornerpos.y()-posxy.y()))*deltaXY.y())
0453           );
0454       */
0455       x[3 - jc] = xy[0];
0456       y[3 - jc] = xy[1];
0457 
0458 #ifdef PFLOW_DEBUG
0459       if (debug) {
0460         std::cout << "corners " << jc << " " << cornerpos.x() << " " << x[3 - jc] << " " << cornerpos.y() << " "
0461                   << y[3 - jc] << std::endl;
0462       }
0463 #endif
0464     }  //loop corners
0465 
0466     //need to close the polygon in order to
0467     //use the TMath::IsInside fonction from root lib
0468     x[4] = x[0];
0469     y[4] = y[0];
0470 
0471     //Check if the extrapolation point of the track falls
0472     //within the rechit boundaries
0473     bool isinside = TMath::IsInside(xPS, yPS, 5, x, y);
0474 
0475     if (isinside) {
0476       linkedbyrechit = true;
0477       break;
0478     }
0479 
0480   }  //loop rechits
0481 
0482   if (linkedbyrechit) {
0483 #ifdef PFLOW_DEBUG
0484     if (debug)
0485       std::cout << "Cluster PS and Cluster ECAL LINKED BY RECHIT" << std::endl;
0486 #endif
0487     constexpr double scale = 1. / 1000.;
0488     double dist = computeDist(xECAL * scale, yECAL * scale, xPS * scale, yPS * scale, false);
0489     return dist;
0490   } else {
0491     return -1.;
0492   }
0493 }
0494 
0495 double LinkByRecHit::testHFEMAndHFHADByRecHit(const reco::PFCluster& clusterHFEM,
0496                                               const reco::PFCluster& clusterHFHAD,
0497                                               bool debug) {
0498   const auto& posxyzEM = clusterHFEM.position();
0499   const auto& posxyzHAD = clusterHFHAD.position();
0500 
0501   double sameZ = posxyzEM.Z() * posxyzHAD.Z();
0502   if (sameZ < 0)
0503     return -1.;
0504 
0505   double dX = posxyzEM.X() - posxyzHAD.X();
0506   double dY = posxyzEM.Y() - posxyzHAD.Y();
0507 
0508   double dist2 = dX * dX + dY * dY;
0509 
0510   if (dist2 < 0.1) {
0511     // less than one mm
0512     double dist = sqrt(dist2);
0513     return dist;
0514     ;
0515   } else
0516     return -1.;
0517 }
0518 
0519 double LinkByRecHit::computeDist(double eta1, double phi1, double eta2, double phi2, bool etaPhi) {
0520   auto phicor = etaPhi ? normalizedPhi(phi1 - phi2) : phi1 - phi2;
0521   auto etadiff = eta1 - eta2;
0522 
0523   // double chi2 =
0524   //  (eta1 - eta2)*(eta1 - eta2) / ( reta1*reta1+ reta2*reta2 ) +
0525   //  phicor*phicor / ( rphi1*rphi1+ rphi2*rphi2 );
0526 
0527   return std::sqrt(etadiff * etadiff + phicor * phicor);
0528 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!