ConversionSeedGenerators/plugins/HitQuadrupletGeneratorFromLayerPairForPhotonConversion.cc

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