ConversionSeedGenerators/plugins/HitPairGeneratorFromLayerPairForPhotonConversion.cc

0001 #include "HitPairGeneratorFromLayerPairForPhotonConversion.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 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/Framework/interface/EventSetup.h"
0015
0016 using namespace GeomDetEnumerators;
0017 using namespace std;
0018 #include "RecoTracker/TkMSParametrization/interface/PixelRecoUtilities.h"
0019
0020 // #define mydebug_Seed
0021
0022 typedef PixelRecoRange<float> Range;
0023
0024 namespace {
0025   template <class T>
0026   inline T sqr(T t) {
0027     return t * t;
0028   }
0029 }  // namespace
0030
0031 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
0032 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0033 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0034 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
0035 #include "FWCore/Framework/interface/ESHandle.h"
0036 #include "DataFormats/Math/interface/deltaPhi.h"
0037
0038 HitPairGeneratorFromLayerPairForPhotonConversion::HitPairGeneratorFromLayerPairForPhotonConversion(
0039     edm::ConsumesCollector iC,
0040     unsigned int inner,
0041     unsigned int outer,
0042     LayerCacheType* layerCache,
0043     unsigned int nSize,
0044     unsigned int max)
0045     : theLayerCache(*layerCache),
0046       theOuterLayer(outer),
0047       theInnerLayer(inner),
0048       theMaxElement(max),
0049       theFieldToken(iC.esConsumes()) {}
0050
0051 void HitPairGeneratorFromLayerPairForPhotonConversion::hitPairs(const ConversionRegion& convRegion,
0052                                                                 const TrackingRegion& region,
0053                                                                 OrderedHitPairs& result,
0054                                                                 const Layers& layers,
0055                                                                 const edm::Event& event,
0056                                                                 const edm::EventSetup& es) {
0057   auto oldSize = result.size();
0058   auto maxNum = result.size() + theMaxElement;
0059
0060 #ifdef mydebug_Seed
0061   ss.str("");
0062 #endif
0063
0064   typedef OrderedHitPair::InnerRecHit InnerHit;
0065   typedef OrderedHitPair::OuterRecHit OuterHit;
0066   typedef RecHitsSortedInPhi::Hit Hit;
0067
0068   Layer innerLayerObj = layers[theInnerLayer];
0069   Layer outerLayerObj = layers[theOuterLayer];
0070
0071 #ifdef mydebug_Seed
0072   ss << "In " << innerLayerObj.name() << " Out " << outerLayerObj.name() << std::endl;
0073 #endif
0074
0075   if (!checkBoundaries(*innerLayerObj.detLayer(), convRegion, 40., 60.))
0076     return;  //FIXME, the maxSearchR(Z) are not optimized
0077   if (!checkBoundaries(*outerLayerObj.detLayer(), convRegion, 50., 60.))
0078     return;  //FIXME, the maxSearchR(Z) are not optimized
0079
0080   /*get hit sorted in phi for each layer: NB: doesn't apply any region cut*/
0081   const RecHitsSortedInPhi& innerHitsMap = theLayerCache(innerLayerObj, region);
0082   if (innerHitsMap.empty())
0083     return;
0084
0085   const RecHitsSortedInPhi& outerHitsMap = theLayerCache(outerLayerObj, region);
0086   if (outerHitsMap.empty())
0087     return;
0088   /*----------------*/
0089
0090   /*This object will check the compatibility of the his in phi among the two layers. */
0091   //InnerDeltaPhi deltaPhi(innerLayerObj.detLayer(), region, es);
0092
0093   static const float nSigmaRZ = std::sqrt(12.f);
0094   //  static const float nSigmaPhi = 3.f;
0095   vector<RecHitsSortedInPhi::Hit> innerHits, outerHits;
0096   float outerPhimin, outerPhimax;
0097   float innerPhimin, innerPhimax;
0098
0099   /*Getting only the Hits in the outer layer that are compatible with the conversion region*/
0100   const auto& field = es.getData(theFieldToken);
0101   if (!getPhiRange(outerPhimin, outerPhimax, *outerLayerObj.detLayer(), convRegion, field))
0102     return;
0103   outerHitsMap.hits(outerPhimin, outerPhimax, outerHits);
0104
0105 #ifdef mydebug_Seed
0106   ss << "\tophimin, ophimax " << outerPhimin << " " << outerPhimax << std::endl;
0107 #endif
0108
0109   /* loop on outer hits*/
0110   for (vector<RecHitsSortedInPhi::Hit>::const_iterator oh = outerHits.begin(); oh != outerHits.end(); ++oh) {
0111     RecHitsSortedInPhi::Hit ohit = (*oh);
0112 #ifdef mydebug_Seed
0113     GlobalPoint oPos = ohit->globalPosition();
0114
0115     ss << "\toPos " << oPos << " r " << oPos.perp() << " phi " << oPos.phi() << " cotTheta " << oPos.z() / oPos.perp()
0116        << std::endl;
0117 #endif
0118
0119     /*Check the compatibility of the ohit with the eta of the seeding track*/
0120     if (checkRZCompatibilityWithSeedTrack(ohit, *outerLayerObj.detLayer(), convRegion))
0121       continue;
0122
0123     /*
0124     //Do I need this? it uses a compatibility that probably I wouldn't
0125     //Removing for the time being
0126
0127     PixelRecoRange<float> phiRange = deltaPhi( oPos.perp(), oPos.phi(), oPos.z(), nSigmaPhi*(ohit->errorGlobalRPhi()));
0128     if (phiRange.empty()) continue;
0129     */
0130
0131     std::unique_ptr<const HitRZCompatibility> checkRZ = region.checkRZ(innerLayerObj.detLayer(), ohit);
0132     if (!checkRZ) {
0133 #ifdef mydebug_Seed
0134       ss << "*******\nNo valid checkRZ\n*******" << std::endl;
0135 #endif
0136       continue;
0137     }
0138
0139     /*Get only the inner hits compatible with the conversion region*/
0140     innerHits.clear();
0141     if (!getPhiRange(innerPhimin, innerPhimax, *innerLayerObj.detLayer(), convRegion, field))
0142       continue;
0143     innerHitsMap.hits(innerPhimin, innerPhimax, innerHits);
0144
0145 #ifdef mydebug_Seed
0146     ss << "\tiphimin, iphimax " << innerPhimin << " " << innerPhimax << std::endl;
0147 #endif
0148
0149     /*Loop on inner hits*/
0150     for (vector<RecHitsSortedInPhi::Hit>::const_iterator ih = innerHits.begin(), ieh = innerHits.end(); ih < ieh;
0151          ++ih) {
0152       GlobalPoint innPos = (*ih)->globalPosition();
0153
0154 #ifdef mydebug_Seed
0155       ss << "\tinnPos " << innPos << " r " << innPos.perp() << " phi " << innPos.phi() << " cotTheta "
0156          << innPos.z() / innPos.perp() << std::endl;
0157 #endif
0158
0159       /*Check the compatibility of the ohit with the eta of the seeding track*/
0160       if (checkRZCompatibilityWithSeedTrack(*ih, *innerLayerObj.detLayer(), convRegion))
0161         continue;
0162
0163       float r_reduced = std::sqrt(sqr(innPos.x() - region.origin().x()) + sqr(innPos.y() - region.origin().y()));
0164       Range allowed;
0165       Range hitRZ;
0166       if (innerLayerObj.detLayer()->location() == barrel) {
0167         allowed = checkRZ->range(r_reduced);
0168         float zErr = nSigmaRZ * (*ih)->errorGlobalZ();
0169         hitRZ = Range(innPos.z() - zErr, innPos.z() + zErr);
0170       } else {
0171         allowed = checkRZ->range(innPos.z());
0172         float rErr = nSigmaRZ * (*ih)->errorGlobalR();
0173         hitRZ = Range(r_reduced - rErr, r_reduced + rErr);
0174       }
0175       Range crossRange = allowed.intersection(hitRZ);
0176
0177 #ifdef mydebug_Seed
0178       ss << "\n\t\t allowed Range " << allowed.min() << " \t, " << allowed.max() << "\n\t\t hitRz   Range "
0179          << hitRZ.min() << " \t, " << hitRZ.max() << "\n\t\t Cross   Range " << crossRange.min() << " \t, "
0180          << crossRange.max() << "\n\t\t the seed track has origin " << convRegion.convPoint() << " \t cotTheta "
0181          << convRegion.cotTheta() << std::endl;
0182 #endif
0183
0184       if (!crossRange.empty()) {
0185 #ifdef mydebug_Seed
0186         ss << "\n\t\t !!!!ACCEPTED!!! \n\n";
0187 #endif
0188         if (theMaxElement != 0 && result.size() >= maxNum) {
0189           result.resize(oldSize);
0190 #ifdef mydebug_Seed
0191           edm::LogError("TooManySeeds") << "number of pairs exceed maximum, no pairs produced";
0192           std::cout << ss.str();
0193 #endif
0194           return;
0195         }
0196         result.push_back(OrderedHitPair(*ih, ohit));
0197       }
0198     }
0199   }
0200
0201 #ifdef mydebug_Seed
0202   std::cout << ss.str();
0203 #endif
0204 }
0205
0206 float HitPairGeneratorFromLayerPairForPhotonConversion::getLayerRadius(const DetLayer& layer) {
0207   if (layer.location() == GeomDetEnumerators::barrel) {
0208     const BarrelDetLayer& bl = static_cast<const BarrelDetLayer&>(layer);
0209     float rLayer = bl.specificSurface().radius();
0210
0211     // the maximal delta phi will be for the innermost hits
0212     float theThickness = layer.surface().bounds().thickness();
0213     return rLayer + 0.5f * theThickness;
0214   }
0215
0216   //Fixme
0217   return 0;
0218 }
0219
0220 float HitPairGeneratorFromLayerPairForPhotonConversion::getLayerZ(const DetLayer& layer) {
0221   if (layer.location() == GeomDetEnumerators::endcap) {
0222     float layerZ = layer.position().z();
0223     float theThickness = layer.surface().bounds().thickness();
0224     float layerZmax = layerZ > 0 ? layerZ + 0.5f * theThickness : layerZ - 0.5f * theThickness;
0225     return layerZmax;
0226   } else {
0227     //Fixme
0228     return 0;
0229   }
0230 }
0231 bool HitPairGeneratorFromLayerPairForPhotonConversion::checkBoundaries(const DetLayer& layer,
0232                                                                        const ConversionRegion& convRegion,
0233                                                                        float maxSearchR,
0234                                                                        float maxSearchZ) {
0235   if (layer.location() == GeomDetEnumerators::barrel) {
0236     float minZEndCap = 130;
0237     if (fabs(convRegion.convPoint().z()) > minZEndCap) {
0238 #ifdef mydebug_Seed
0239       ss << "\tthe conversion seems to be in the endcap. Zconv " << convRegion.convPoint().z() << std::endl;
0240       std::cout << ss.str();
0241 #endif
0242       return false;
0243     }
0244
0245     float R = getLayerRadius(layer);
0246
0247     if (convRegion.convPoint().perp() > R) {
0248 #ifdef mydebug_Seed
0249       ss << "\tthis layer is before the conversion : R layer " << R << " [ Rconv " << convRegion.convPoint().perp()
0250          << " Zconv " << convRegion.convPoint().z() << std::endl;
0251       std::cout << ss.str();
0252 #endif
0253       return false;
0254     }
0255
0256     if (R - convRegion.convPoint().perp() > maxSearchR) {
0257 #ifdef mydebug_Seed
0258       ss << "\tthis layer is far from the conversion more than cut " << maxSearchR << " cm. R layer " << R
0259          << " [ Rconv " << convRegion.convPoint().perp() << " Zconv " << convRegion.convPoint().z() << std::endl;
0260       std::cout << ss.str();
0261 #endif
0262       return false;
0263     }
0264
0265   } else if (layer.location() == GeomDetEnumerators::endcap) {
0266     float Z = getLayerZ(layer);
0267     if ((convRegion.convPoint().z() > 0 && convRegion.convPoint().z() > Z) ||
0268         (convRegion.convPoint().z() < 0 && convRegion.convPoint().z() < Z)) {
0269 #ifdef mydebug_Seed
0270       ss << "\tthis layer is before the conversion : Z layer " << Z << " [ Rconv " << convRegion.convPoint().perp()
0271          << " Zconv " << convRegion.convPoint().z() << std::endl;
0272       std::cout << ss.str();
0273 #endif
0274       return false;
0275     }
0276
0277     if (fabs(Z - convRegion.convPoint().z()) > maxSearchZ) {
0278 #ifdef mydebug_Seed
0279       ss << "\tthis layer is far from the conversion more than cut " << maxSearchZ << " cm. Z layer " << Z
0280          << " [ Rconv " << convRegion.convPoint().perp() << " Zconv " << convRegion.convPoint().z() << std::endl;
0281       std::cout << ss.str();
0282 #endif
0283       return false;
0284     }
0285   }
0286   return true;
0287 }
0288
0289 bool HitPairGeneratorFromLayerPairForPhotonConversion::getPhiRange(float& Phimin,
0290                                                                    float& Phimax,
0291                                                                    const DetLayer& layer,
0292                                                                    const ConversionRegion& convRegion,
0293                                                                    const MagneticField& field) {
0294   if (layer.location() == GeomDetEnumerators::barrel) {
0295     return getPhiRange(Phimin, Phimax, getLayerRadius(layer), convRegion, field);
0296   } else if (layer.location() == GeomDetEnumerators::endcap) {
0297     float Z = getLayerZ(layer);
0298     float R = Z / convRegion.cotTheta();
0299     return getPhiRange(Phimin, Phimax, R, convRegion, field);  //FIXME
0300   }
0301   return false;
0302 }
0303
0304 bool HitPairGeneratorFromLayerPairForPhotonConversion::getPhiRange(
0305     float& Phimin, float& Phimax, const float& layerR, const ConversionRegion& convRegion, const MagneticField& field) {
0306   Phimin = reco::deltaPhi(convRegion.convPoint().phi(), 0.);
0307
0308   float dphi;
0309   float ptmin = 0.1;
0310   float DeltaL = layerR - convRegion.convPoint().perp();
0311
0312   if (DeltaL < 0) {
0313     Phimin = 0;
0314     Phimax = 0;
0315     return false;
0316   }
0317
0318   float theRCurvatureMin = PixelRecoUtilities::bendingRadius(ptmin, field);
0319
0320   if (theRCurvatureMin < DeltaL)
0321     dphi = atan(DeltaL / layerR);
0322   else
0323     dphi = atan(theRCurvatureMin / layerR * (1 - sqrt(1 - sqr(DeltaL / theRCurvatureMin))));
0324
0325   if (convRegion.charge() > 0) {
0326     Phimax = Phimin;
0327     Phimin = Phimax - dphi;
0328   } else {
0329     Phimax = Phimin + dphi;
0330   }
0331
0332   //std::cout << dphi << " " << Phimin << " " << Phimax << " " << layerR << " " << DeltaL  << " " << convRegion.convPoint().phi() << " " << convRegion.convPoint().perp()<< std::endl;
0333   return true;
0334 }
0335
0336 bool HitPairGeneratorFromLayerPairForPhotonConversion::checkRZCompatibilityWithSeedTrack(
0337     const RecHitsSortedInPhi::Hit& hit, const DetLayer& layer, const ConversionRegion& convRegion) {
0338   static const float nSigmaRZ = std::sqrt(12.f);
0339   Range hitCotTheta;
0340
0341   double sigmaCotTheta = convRegion.errTheta() *
0342                          (1 + convRegion.cotTheta() * convRegion.cotTheta());  //Error Propagation from sigma theta.
0343   Range allowedCotTheta(convRegion.cotTheta() - nSigmaRZ * sigmaCotTheta,
0344                         convRegion.cotTheta() + nSigmaRZ * sigmaCotTheta);
0345
0346   double dz = hit->globalPosition().z() - convRegion.pvtxPoint().z();
0347   double r_reduced = std::sqrt(sqr(hit->globalPosition().x() - convRegion.pvtxPoint().x()) +
0348                                sqr(hit->globalPosition().y() - convRegion.pvtxPoint().y()));
0349
0350   if (layer.location() == GeomDetEnumerators::barrel) {
0351     float zErr = nSigmaRZ * hit->errorGlobalZ();
0352     hitCotTheta = Range(getCot(dz - zErr, r_reduced), getCot(dz + zErr, r_reduced));
0353   } else {
0354     float rErr = nSigmaRZ * hit->errorGlobalR();
0355     if (dz > 0)
0356       hitCotTheta = Range(getCot(dz, r_reduced + rErr), getCot(dz, r_reduced - rErr));
0357     else
0358       hitCotTheta = Range(getCot(dz, r_reduced - rErr), getCot(dz, r_reduced + rErr));
0359   }
0360
0361   Range crossRange = allowedCotTheta.intersection(hitCotTheta);
0362
0363 #ifdef mydebug_Seed
0364   ss << "\n\t\t cotTheta allowed Range " << allowedCotTheta.min() << " \t, " << allowedCotTheta.max()
0365      << "\n\t\t hitCotTheta   Range " << hitCotTheta.min() << " \t, " << hitCotTheta.max() << "\n\t\t Cross   Range "
0366      << crossRange.min() << " \t, " << crossRange.max() << "\n\t\t the seed track has origin " << convRegion.convPoint()
0367      << " \t cotTheta " << convRegion.cotTheta() << std::endl;
0368 #endif
0369
0370   return crossRange.empty();
0371 }
0372
0373 double HitPairGeneratorFromLayerPairForPhotonConversion::getCot(double dz, double dr) {
0374   if (std::abs(dr) > 1.e-4f)
0375     return dz / dr;
0376   else if (dz > 0)
0377     return 99999.f;
0378   else
0379     return -99999.f;
0380 }