Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-10-01 22:41:00 

0001 #include "HitQuadrupletGeneratorFromLayerPairForPhotonConversion.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 
0004 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0005 #include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
0006 #include "TrackingTools/DetLayers/interface/ForwardDetLayer.h"
0007 
0008 #include "RecoTracker/TkTrackingRegions/interface/HitRZCompatibility.h"
0009 #include "RecoTracker/TkTrackingRegions/interface/TrackingRegion.h"
0010 #include "RecoTracker/TkTrackingRegions/interface/TrackingRegionBase.h"
0011 #include "RecoTracker/TkHitPairs/interface/OrderedHitPairs.h"
0012 
0013 using namespace GeomDetEnumerators;
0014 using namespace std;
0015 #include "RecoTracker/TkMSParametrization/interface/PixelRecoUtilities.h"
0016 
0017 //#define mydebug_QSeed
0018 
0019 typedef PixelRecoRange<float> Range;
0020 template <class T>
0021 inline T sqr(T t) {
0022   return t * t;
0023 }
0024 
0025 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
0026 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0027 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0028 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
0029 #include "FWCore/Framework/interface/ESHandle.h"
0030 #include "DataFormats/Math/interface/deltaPhi.h"
0031 
0032 HitQuadrupletGeneratorFromLayerPairForPhotonConversion::HitQuadrupletGeneratorFromLayerPairForPhotonConversion(
0033     unsigned int inner, unsigned int outer, LayerCacheType *layerCache, unsigned int max)
0034     : theLayerCache(*layerCache), theOuterLayer(outer), theInnerLayer(inner), theMaxElement(max) {
0035   ss = new std::stringstream;
0036 }
0037 
0038 void HitQuadrupletGeneratorFromLayerPairForPhotonConversion::hitPairs(const TrackingRegion &region,
0039                                                                       OrderedHitPairs &result,
0040                                                                       const Layers &layers) {
0041   ss->str("");
0042 
0043   typedef OrderedHitPair::InnerRecHit InnerHit;
0044   typedef OrderedHitPair::OuterRecHit OuterHit;
0045   typedef RecHitsSortedInPhi::Hit Hit;
0046 
0047   Layer innerLayerObj = layers[theInnerLayer];
0048   Layer outerLayerObj = layers[theOuterLayer];
0049 
0050   size_t totCountP2 = 0, totCountP1 = 0, totCountM2 = 0, totCountM1 = 0, selCount = 0;
0051 #ifdef mydebug_QSeed
0052   (*ss) << "In " << innerLayerObj.name() << " Out " << theOuterLayer.name() << std::endl;
0053 #endif
0054 
0055   /*get hit sorted in phi for each layer: NB: doesn't apply any region cut*/
0056   const RecHitsSortedInPhi &innerHitsMap = theLayerCache(innerLayerObj, region);
0057   if (innerHitsMap.empty())
0058     return;
0059 
0060   const RecHitsSortedInPhi &outerHitsMap = theLayerCache(outerLayerObj, region);
0061   if (outerHitsMap.empty())
0062     return;
0063   /*----------------*/
0064 
0065   /*This object will check the compatibility of the his in phi among the two layers. */
0066   //InnerDeltaPhi deltaPhi(*innerLayerObj.detLayer(), region, es);
0067 
0068   vector<RecHitsSortedInPhi::Hit> innerHits;
0069   //  float outerPhimin, outerPhimax;
0070   float innerPhimin, innerPhimax;
0071   float maxDeltaPhi = 1.;  //sguazz
0072   //float maxDeltaPhi=1.;
0073 
0074   RecHitsSortedInPhi::Range outerHits = outerHitsMap.all();
0075 
0076   RecHitsSortedInPhi::HitIter nextoh;
0077   for (RecHitsSortedInPhi::HitIter oh = outerHits.first; oh != outerHits.second; ++oh) {
0078     RecHitsSortedInPhi::Hit ohit = (*oh).hit();
0079     GlobalPoint oPos = ohit->globalPosition();
0080 
0081     totCountP2++;
0082     std::unique_ptr<const HitRZCompatibility> checkRZ = region.checkRZ(outerLayerObj.detLayer(), ohit);
0083     for (nextoh = oh + 1; nextoh != outerHits.second; ++nextoh) {
0084       RecHitsSortedInPhi::Hit nohit = (*nextoh).hit();
0085       GlobalPoint noPos = nohit->globalPosition();
0086 
0087 #ifdef mydebug_QSeed
0088       (*ss) << "\toPos " << oPos << " r " << oPos.perp() << " phi " << oPos.phi() << " cotTheta "
0089             << oPos.z() / oPos.perp() << std::endl;
0090       (*ss) << "\tnoPos " << noPos << " r " << noPos.perp() << " phi " << noPos.phi() << " cotTheta "
0091             << noPos.z() / noPos.perp() << std::endl;
0092       (*ss) << "\tDeltaPhi " << reco::deltaPhi(noPos.barePhi(), oPos.barePhi()) << std::endl;
0093 #endif
0094 
0095       if (fabs(reco::deltaPhi(noPos.barePhi(), oPos.barePhi())) > maxDeltaPhi)
0096         break;
0097 
0098       totCountM2++;
0099 
0100       /*Check the compatibility of the ohit with the eta of the seeding track*/
0101       if (failCheckRZCompatibility(nohit, *outerLayerObj.detLayer(), checkRZ.get(), region))
0102         continue;
0103 
0104       /*  
0105     //Do I need this? it uses a compatibility that probably I wouldn't 
0106     //Removing for the time being
0107 
0108     PixelRecoRange<float> phiRange = deltaPhi( oPos.perp(), oPos.phi(), oPos.z(), nSigmaPhi*(ohit->errorGlobalRPhi()));    
0109     if (phiRange.empty()) continue;
0110     */
0111 
0112       /*Get only the inner hits compatible with the conversion region*/
0113       innerPhimin = ohit->globalPosition().phi();
0114       innerPhimax = nohit->globalPosition().phi();
0115       // checkPhiRange(innerPhimin,innerPhimax);
0116 
0117       innerHits.clear();
0118       innerHitsMap.hits(innerPhimin, innerPhimax, innerHits);
0119 
0120 #ifdef mydebug_QSeed
0121       (*ss) << "\tiphimin, iphimax " << innerPhimin << " " << innerPhimax << std::endl;
0122 #endif
0123 
0124       std::unique_ptr<const HitRZCompatibility> checkRZb = region.checkRZ(innerLayerObj.detLayer(), ohit);
0125       std::unique_ptr<const HitRZCompatibility> checkRZc = region.checkRZ(innerLayerObj.detLayer(), nohit);
0126 
0127       /*Loop on inner hits*/
0128       vector<RecHitsSortedInPhi::Hit>::const_iterator ieh = innerHits.end();
0129       for (vector<RecHitsSortedInPhi::Hit>::const_iterator ih = innerHits.begin(); ih < ieh; ++ih) {
0130         RecHitsSortedInPhi::Hit ihit = *ih;
0131 
0132 #ifdef mydebug_QSeed
0133         GlobalPoint innPos = (*ih)->globalPosition();
0134         (*ss) << "\toPos " << oPos << " r " << oPos.perp() << " phi " << oPos.phi() << " cotTheta "
0135               << oPos.z() / oPos.perp() << std::endl;
0136         (*ss) << "\tnoPos " << noPos << " r " << noPos.perp() << " phi " << noPos.phi() << " cotTheta "
0137               << noPos.z() / noPos.perp() << std::endl;
0138         (*ss) << "\tinnPos " << innPos << " r " << innPos.perp() << " phi " << innPos.phi() << " cotTheta "
0139               << innPos.z() / innPos.perp() << std::endl;
0140 #endif
0141 
0142         totCountP1++;
0143 
0144         /*Check the compatibility of the ihit with the two outer hits*/
0145         if (failCheckRZCompatibility(ihit, *innerLayerObj.detLayer(), checkRZb.get(), region) ||
0146             failCheckRZCompatibility(ihit, *innerLayerObj.detLayer(), checkRZc.get(), region))
0147           continue;
0148 
0149         for (vector<RecHitsSortedInPhi::Hit>::const_iterator nextih = ih + 1; nextih != ieh; ++nextih) {
0150           RecHitsSortedInPhi::Hit nihit = *nextih;
0151 
0152 #ifdef mydebug_QSeed
0153           GlobalPoint ninnPos = (*nextih)->globalPosition();
0154           (*ss) << "\toPos " << oPos << " r " << oPos.perp() << " phi " << oPos.phi() << " cotTheta "
0155                 << oPos.z() / oPos.perp() << std::endl;
0156           (*ss) << "\tnoPos " << noPos << " r " << noPos.perp() << " phi " << noPos.phi() << " cotTheta "
0157                 << noPos.z() / noPos.perp() << std::endl;
0158           (*ss) << "\tinnPos " << innPos << " r " << innPos.perp() << " phi " << innPos.phi() << " cotTheta "
0159                 << innPos.z() / innPos.perp() << std::endl;
0160           (*ss) << "\tninnPos " << ninnPos << " r " << ninnPos.perp() << " phi " << ninnPos.phi() << " cotTheta "
0161                 << ninnPos.z() / ninnPos.perp() << std::endl;
0162 #endif
0163 
0164           totCountM1++;
0165 
0166           /*Check the compatibility of the nihit with the two outer hits*/
0167           if (failCheckRZCompatibility(nihit, *innerLayerObj.detLayer(), checkRZb.get(), region) ||
0168               failCheckRZCompatibility(nihit, *innerLayerObj.detLayer(), checkRZc.get(), region))
0169             continue;
0170 
0171           /*Sguazz modifica qui*/
0172           if (failCheckSlopeTest(ohit, nohit, ihit, nihit, region))
0173             continue;
0174 
0175           if (theMaxElement != 0 && result.size() >= theMaxElement) {
0176             result.clear();
0177             edm::LogError("TooManyQuads") << "number of Quad combinations exceed maximum, no quads produced";
0178 #ifdef mydebug_QSeed
0179             (*ss) << "In " << innerLayerObj.name() << " Out " << outerLayerObj.name() << "\tP2 " << totCountP2
0180                   << "\tM2 " << totCountM2 << "\tP1 " << totCountP1 << "\tM1 " << totCountM1 << "\tsel " << selCount
0181                   << std::endl;
0182             std::cout << (*ss).str();
0183 #endif
0184             return;
0185           }
0186           selCount++;
0187           result.push_back(OrderedHitPair(ihit, ohit));
0188           result.push_back(OrderedHitPair(nihit, nohit));
0189           //#ifdef mydebug_QSeed
0190           //(*ss) << "sizeOfresul " << result.size() << std::endl;
0191           //#endif
0192         }
0193       }
0194     }
0195   }
0196 #ifdef mydebug_QSeed
0197   (*ss) << "In " << innerLayerObj.name() << " Out " << outerLayerObj.name() << "\tP2 " << totCountP2 << "\tM2 "
0198         << totCountM2 << "\tP1 " << totCountP1 << "\tM1 " << totCountM1 << "\tsel " << selCount << std::endl;
0199   std::cout << (*ss).str();
0200 #endif
0201 }
0202 
0203 bool HitQuadrupletGeneratorFromLayerPairForPhotonConversion::failCheckRZCompatibility(const RecHitsSortedInPhi::Hit &hit,
0204                                                                                       const DetLayer &layer,
0205                                                                                       const HitRZCompatibility *checkRZ,
0206                                                                                       const TrackingRegion &region) {
0207   if (!checkRZ) {
0208 #ifdef mydebug_QSeed
0209     (*ss) << "*******\nNo valid checkRZ\n*******" << std::endl;
0210 #endif
0211     return true;
0212   }
0213 
0214   static const float nSigmaRZ = std::sqrt(12.f);
0215   float r_reduced = std::sqrt(sqr(hit->globalPosition().x() - region.origin().x()) +
0216                               sqr(hit->globalPosition().y() - region.origin().y()));
0217   Range allowed;
0218   Range hitRZ;
0219   if (layer.location() == barrel) {
0220     allowed = checkRZ->range(r_reduced);
0221     float zErr = nSigmaRZ * hit->errorGlobalZ();
0222     hitRZ = Range(hit->globalPosition().z() - zErr, hit->globalPosition().z() + zErr);
0223   } else {
0224     allowed = checkRZ->range(hit->globalPosition().z());
0225     float rErr = nSigmaRZ * hit->errorGlobalR();
0226     hitRZ = Range(r_reduced - rErr, r_reduced + rErr);
0227   }
0228   Range crossRange = allowed.intersection(hitRZ);
0229 
0230 #ifdef mydebug_QSeed
0231   (*ss) << "\n\t\t allowed Range " << allowed.min() << " \t, " << allowed.max() << "\n\t\t hitRz   Range "
0232         << hitRZ.min() << " \t, " << hitRZ.max() << "\n\t\t Cross   Range " << crossRange.min() << " \t, "
0233         << crossRange.max() << std::endl;
0234   if (!crossRange.empty())
0235     (*ss) << "\n\t\t !!!!ACCEPTED!!! \n\n";
0236 #endif
0237 
0238   return crossRange.empty();
0239 }
0240 
0241 bool HitQuadrupletGeneratorFromLayerPairForPhotonConversion::failCheckSlopeTest(const RecHitsSortedInPhi::Hit &ohit,
0242                                                                                 const RecHitsSortedInPhi::Hit &nohit,
0243                                                                                 const RecHitsSortedInPhi::Hit &ihit,
0244                                                                                 const RecHitsSortedInPhi::Hit &nihit,
0245                                                                                 const TrackingRegion &region) {
0246   double r[5], z[5], ez[5];
0247   //  double pr[2], pz[2], e2pz[2], mr[2], mz[2], e2mz[2];
0248 
0249   //
0250   //Hits
0251   r[0] = ohit->globalPosition().perp();
0252   z[0] = ohit->globalPosition().z();
0253   ez[0] = getEffectiveErrorOnZ(ohit, region);
0254   //
0255   r[1] = nohit->globalPosition().perp();
0256   z[1] = nohit->globalPosition().z();
0257   ez[1] = getEffectiveErrorOnZ(nohit, region);
0258   //
0259   r[2] = nihit->globalPosition().perp();
0260   z[2] = nihit->globalPosition().z();
0261   ez[2] = getEffectiveErrorOnZ(nihit, region);
0262   //
0263   r[3] = ihit->globalPosition().perp();
0264   z[3] = ihit->globalPosition().z();
0265   ez[3] = getEffectiveErrorOnZ(ihit, region);
0266   //
0267   //R (r) ordering of the 4 hit arrays
0268   bubbleSortVsR(4, r, z, ez);
0269   //
0270   //Vertex
0271   r[4] = region.origin().perp();
0272   z[4] = region.origin().z();
0273   double vError = region.originZBound();
0274   if (vError > 15.)
0275     vError = 1.;
0276   ez[4] = 3. * vError;
0277 
0278   //
0279   //Sequence of checks
0280   //
0281   //Inner segment == vertex
0282   double rInn = r[4];
0283   double zInnMin = z[4] - ez[4];
0284   double zInnMax = z[4] + ez[4];
0285   //
0286   // Int == 2, Out == 3
0287   double rOut = r[3];
0288   double zOutMin = z[3] - ez[3];
0289   double zOutMax = z[3] + ez[3];
0290   double rInt = r[2];
0291   double zIntMin = z[2] - ez[2];
0292   double zIntMax = z[2] + ez[2];
0293   if (failCheckSegmentZCompatibility(rInn, zInnMin, zInnMax, rInt, zIntMin, zIntMax, rOut, zOutMin, zOutMax))
0294     return true;
0295   //
0296   // Int == 1, Out == 2 (with updated limits)
0297   rOut = rInt;
0298   zOutMin = zIntMin;
0299   zOutMax = zIntMax;
0300   rInt = r[1];
0301   zIntMin = z[1] - ez[1];
0302   zIntMax = z[1] + ez[1];
0303   if (failCheckSegmentZCompatibility(rInn, zInnMin, zInnMax, rInt, zIntMin, zIntMax, rOut, zOutMin, zOutMax))
0304     return true;
0305   //
0306   // Int == 0, Out == 1 (with updated limits)
0307   rOut = rInt;
0308   zOutMin = zIntMin;
0309   zOutMax = zIntMax;
0310   rInt = r[0];
0311   zIntMin = z[0] - ez[0];
0312   zIntMax = z[0] + ez[0];
0313   if (failCheckSegmentZCompatibility(rInn, zInnMin, zInnMax, rInt, zIntMin, zIntMax, rOut, zOutMin, zOutMax))
0314     return true;
0315 
0316   //
0317   // Test is ok!!!
0318   return false;
0319 }
0320 
0321 double HitQuadrupletGeneratorFromLayerPairForPhotonConversion::verySimpleFit(
0322     int size, double *ax, double *ay, double *e2y, double &p0, double &e2p0, double &p1) {
0323   //#include "verySimpleFit.icc"
0324   return 0;
0325 }
0326 
0327 double HitQuadrupletGeneratorFromLayerPairForPhotonConversion::getSqrEffectiveErrorOnZ(
0328     const RecHitsSortedInPhi::Hit &hit, const TrackingRegion &region) {
0329   //
0330   //Fit-wise the effective error on Z is approximately the sum in quadrature of the error on Z
0331   //and the error on R correctly projected by using hit-vertex direction
0332 
0333   double sqrProjFactor =
0334       sqr((hit->globalPosition().z() - region.origin().z()) / (hit->globalPosition().perp() - region.origin().perp()));
0335   return (hit->globalPositionError().czz() + sqrProjFactor * hit->globalPositionError().rerr(hit->globalPosition()));
0336 }
0337 
0338 double HitQuadrupletGeneratorFromLayerPairForPhotonConversion::getEffectiveErrorOnZ(const RecHitsSortedInPhi::Hit &hit,
0339                                                                                     const TrackingRegion &region) {
0340   //
0341   //Fit-wise the effective error on Z is approximately the sum in quadrature of the error on Z
0342   //and the error on R correctly projected by using hit-vertex direction
0343   double sqrProjFactor =
0344       sqr((hit->globalPosition().z() - region.origin().z()) / (hit->globalPosition().perp() - region.origin().perp()));
0345   double effErr =
0346       sqrt(hit->globalPositionError().czz() + sqrProjFactor * hit->globalPositionError().rerr(hit->globalPosition()));
0347   if (effErr > 2.) {
0348     effErr *= 1.8;  //Single Side error is strip length * sqrt( 12.) = 3.46
0349     //empirically found that the error on ss hits value it is already twice as large (two sigmas)!
0350     //Multiply by sqrt(12.)/2.=1.73 to have effErr equal to the strip lenght (1.8 to allow for some margin)
0351     //effErr*=2.5; //Used in some tests
0352   } else {
0353     effErr *= 2.;  //Tight //Double side... allowing for 2 sigma variation
0354     //effErr*=5.; //Loose //Double side... allowing for 2 sigma variation
0355   }
0356   return effErr;
0357 }
0358 
0359 void HitQuadrupletGeneratorFromLayerPairForPhotonConversion::bubbleSortVsR(int n, double *ar, double *az, double *aez) {
0360   bool swapped = true;
0361   int j = 0;
0362   double tmpr, tmpz, tmpez;
0363   while (swapped) {
0364     swapped = false;
0365     j++;
0366     for (int i = 0; i < n - j; i++) {
0367       if (ar[i] > ar[i + 1]) {
0368         tmpr = ar[i];
0369         ar[i] = ar[i + 1];
0370         ar[i + 1] = tmpr;
0371         tmpz = az[i];
0372         az[i] = az[i + 1];
0373         az[i + 1] = tmpz;
0374         tmpez = aez[i];
0375         aez[i] = aez[i + 1];
0376         aez[i + 1] = tmpez;
0377         swapped = true;
0378       }
0379     }
0380   }
0381 }
0382 
0383 bool HitQuadrupletGeneratorFromLayerPairForPhotonConversion::failCheckSegmentZCompatibility(double &rInn,
0384                                                                                             double &zInnMin,
0385                                                                                             double &zInnMax,
0386                                                                                             double &rInt,
0387                                                                                             double &zIntMin,
0388                                                                                             double &zIntMax,
0389                                                                                             double &rOut,
0390                                                                                             double &zOutMin,
0391                                                                                             double &zOutMax) {
0392   //
0393   // Check the compatibility in z of an INTermediate segment between an INNer segment and an OUTer segment;
0394   // when true is returned zIntMin and zIntMax are replaced with allowed range values
0395 
0396   //Left side
0397   double zLeft = getZAtR(rInn, zInnMin, rInt, rOut, zOutMin);
0398   if (zIntMax < zLeft)
0399     return true;
0400   //Right side
0401   double zRight = getZAtR(rInn, zInnMax, rInt, rOut, zOutMax);
0402   if (zIntMin > zRight)
0403     return true;
0404   if (zIntMin < zLeft && zIntMax < zRight) {
0405     zIntMax = zLeft;
0406     return false;
0407   }
0408   if (zIntMin > zLeft && zIntMax > zRight) {
0409     zIntMax = zRight;
0410     return false;
0411   }
0412 
0413   //Segment is fully contained
0414   return false;
0415 }
0416 
0417 double HitQuadrupletGeneratorFromLayerPairForPhotonConversion::getZAtR(
0418     double &rInn, double &zInn, double &r, double &rOut, double &zOut) {
0419   //    z - zInn      r - rInn                                  r - rInn
0420   //  ----------- = ----------- ==> z = zInn + (zOut - zInn) * -----------
0421   //  zOut - zInn   rOut - rInn                                rOut - rInn
0422 
0423   return zInn + (zOut - zInn) * (r - rInn) / (rOut - rInn);
0424 }
[ 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!