CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/RecoTracker/ConversionSeedGenerators/src/HitPairGeneratorFromLayerPairForPhotonConversion.cc

Go to the documentation of this file.
00001 #include "RecoTracker/ConversionSeedGenerators/interface/HitPairGeneratorFromLayerPairForPhotonConversion.h"
00002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00003 
00004 #include "TrackingTools/DetLayers/interface/DetLayer.h"
00005 #include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
00006 #include "TrackingTools/DetLayers/interface/ForwardDetLayer.h"
00007 
00008 #include "RecoTracker/TkTrackingRegions/interface/HitRZCompatibility.h"
00009 #include "RecoTracker/TkTrackingRegions/interface/TrackingRegion.h"
00010 #include "RecoTracker/TkTrackingRegions/interface/TrackingRegionBase.h"
00011 #include "RecoTracker/TkHitPairs/interface/OrderedHitPairs.h"
00012 #include "RecoTracker/TkHitPairs/src/InnerDeltaPhi.h"
00013 
00014 #include "FWCore/Framework/interface/Event.h"
00015 #include "FWCore/Framework/interface/EventSetup.h"
00016 #include "RecoTracker/TkSeedingLayers/interface/SeedingLayer.h"
00017 
00018 using namespace GeomDetEnumerators;
00019 using namespace ctfseeding;
00020 using namespace std;
00021 #include "RecoTracker/TkMSParametrization/interface/PixelRecoUtilities.h"
00022 
00023 //#define mydebug_Seed
00024 
00025 typedef PixelRecoRange<float> Range;
00026 template<class T> inline T sqr( T t) {return t*t;}
00027 
00028 
00029 
00030 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
00031 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
00032 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00033 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
00034 #include "FWCore/Framework/interface/ESHandle.h"
00035 #include "DataFormats/Math/interface/deltaPhi.h"
00036 
00037 HitPairGeneratorFromLayerPairForPhotonConversion::HitPairGeneratorFromLayerPairForPhotonConversion(
00038                                                              const Layer& inner, 
00039                                                              const Layer& outer, 
00040                                                              LayerCacheType* layerCache,
00041                                                              unsigned int nSize,
00042                                                              unsigned int max)
00043   : HitPairGenerator(nSize),
00044     theLayerCache(*layerCache), theOuterLayer(outer), theInnerLayer(inner)
00045 {
00046   theMaxElement=max;
00047   ss = new std::stringstream;
00048 }
00049 
00050 void HitPairGeneratorFromLayerPairForPhotonConversion::hitPairs(const ConversionRegion& convRegion,
00051                                                                 const TrackingRegion & region, OrderedHitPairs & result,
00052                                                                 const edm::Event& event, const edm::EventSetup& es)
00053 {
00054 
00055   ss->str("");
00056 
00057   typedef OrderedHitPair::InnerRecHit InnerHit;
00058   typedef OrderedHitPair::OuterRecHit OuterHit;
00059   typedef RecHitsSortedInPhi::Hit Hit;
00060 
00061 #ifdef mydebug_Seed
00062   (*ss) << "In " << theInnerLayer.name() << " Out " << theOuterLayer.name() << std::endl;
00063 #endif
00064 
00065   if(!checkBoundaries(*theInnerLayer.detLayer(),convRegion,40.,60.)) return; //FIXME, the maxSearchR(Z) are not optimized
00066   if(!checkBoundaries(*theOuterLayer.detLayer(),convRegion,50.,60.)) return; //FIXME, the maxSearchR(Z) are not optimized
00067 
00068   /*get hit sorted in phi for each layer: NB: doesn't apply any region cut*/
00069   const RecHitsSortedInPhi & innerHitsMap = theLayerCache(&theInnerLayer, region, event, es);
00070   if (innerHitsMap.empty()) return;
00071  
00072   const RecHitsSortedInPhi& outerHitsMap = theLayerCache(&theOuterLayer, region, event, es);
00073   if (outerHitsMap.empty()) return;
00074   /*----------------*/
00075 
00076   /*This object will check the compatibility of the his in phi among the two layers. */
00077   //InnerDeltaPhi deltaPhi(*theInnerLayer.detLayer(), region, es);
00078 
00079   static const float nSigmaRZ = std::sqrt(12.f);
00080   //  static const float nSigmaPhi = 3.f;
00081   vector<RecHitsSortedInPhi::Hit> innerHits, outerHits;
00082   float outerPhimin, outerPhimax;
00083   float innerPhimin, innerPhimax;
00084 
00085   /*Getting only the Hits in the outer layer that are compatible with the conversion region*/
00086   if(!getPhiRange(outerPhimin,outerPhimax,*theOuterLayer.detLayer(),convRegion,es)) return;
00087   outerHitsMap.hits( outerPhimin, outerPhimax, outerHits);
00088 
00089 #ifdef mydebug_Seed
00090   (*ss) << "\tophimin, ophimax " << outerPhimin << " " << outerPhimax << std::endl;
00091 #endif
00092 
00093   /* loop on outer hits*/
00094   for (vector<RecHitsSortedInPhi::Hit>::const_iterator oh = outerHits.begin(); oh!= outerHits.end(); ++oh) { 
00095     RecHitsSortedInPhi::Hit ohit = (*oh);
00096 #ifdef mydebug_Seed
00097     GlobalPoint oPos = ohit->globalPosition();  
00098     
00099     (*ss) << "\toPos " << oPos << " r " << oPos.perp() << " phi " << oPos.phi() <<  " cotTheta " << oPos.z()/oPos.perp() << std::endl;
00100 #endif
00101 
00102     /*Check the compatibility of the ohit with the eta of the seeding track*/
00103     if(checkRZCompatibilityWithSeedTrack(ohit,*theOuterLayer.detLayer(),convRegion)) continue;
00104 
00105     /*  
00106     //Do I need this? it uses a compatibility that probably I wouldn't 
00107     //Removing for the time being
00108 
00109     PixelRecoRange<float> phiRange = deltaPhi( oPos.perp(), oPos.phi(), oPos.z(), nSigmaPhi*(ohit->errorGlobalRPhi()));    
00110     if (phiRange.empty()) continue;
00111     */
00112 
00113     const HitRZCompatibility *checkRZ = region.checkRZ(theInnerLayer.detLayer(), ohit, es);
00114     if(!checkRZ) {
00115 #ifdef mydebug_Seed
00116       (*ss) << "*******\nNo valid checkRZ\n*******" << std::endl;
00117 #endif
00118       continue;
00119     }
00120     
00121     /*Get only the inner hits compatible with the conversion region*/
00122     innerHits.clear();
00123     if(!getPhiRange(innerPhimin,innerPhimax,*theInnerLayer.detLayer(),convRegion,es)) continue;
00124     innerHitsMap.hits(innerPhimin, innerPhimax, innerHits);
00125 
00126 #ifdef mydebug_Seed
00127     (*ss) << "\tiphimin, iphimax " << innerPhimin << " " << innerPhimax << std::endl;
00128 #endif    
00129 
00130     /*Loop on inner hits*/
00131     for ( vector<RecHitsSortedInPhi::Hit>::const_iterator ih=innerHits.begin(), ieh = innerHits.end(); ih < ieh; ++ih) {  
00132       GlobalPoint innPos = (*ih)->globalPosition();
00133       
00134 
00135 #ifdef mydebug_Seed
00136       (*ss) << "\tinnPos " << innPos <<  " r " << innPos.perp() << " phi " << innPos.phi() << " cotTheta " << innPos.z()/innPos.perp() <<  std::endl;
00137 #endif
00138 
00139       /*Check the compatibility of the ohit with the eta of the seeding track*/
00140       if(checkRZCompatibilityWithSeedTrack(*ih,*theInnerLayer.detLayer(),convRegion)) continue;
00141       
00142       float r_reduced = std::sqrt( sqr(innPos.x()-region.origin().x())+sqr(innPos.y()-region.origin().y()));
00143       Range allowed;
00144       Range hitRZ;
00145       if (theInnerLayer.detLayer()->location() == barrel) {
00146         allowed = checkRZ->range(r_reduced);
00147         float zErr = nSigmaRZ * (*ih)->errorGlobalZ();
00148         hitRZ = Range(innPos.z()-zErr, innPos.z()+zErr);
00149       } else {
00150         allowed = checkRZ->range(innPos.z());
00151         float rErr = nSigmaRZ * (*ih)->errorGlobalR();
00152         hitRZ = Range(r_reduced-rErr, r_reduced+rErr);
00153       }
00154       Range crossRange = allowed.intersection(hitRZ);
00155 
00156 #ifdef mydebug_Seed      
00157       (*ss) 
00158         << "\n\t\t allowed Range " << allowed.min() << " \t, " << allowed.max() 
00159         << "\n\t\t hitRz   Range " << hitRZ.min()   << " \t, " << hitRZ.max() 
00160         << "\n\t\t Cross   Range " << crossRange.min()   << " \t, " << crossRange.max() 
00161         << "\n\t\t the seed track has origin " << convRegion.convPoint() << " \t cotTheta " << convRegion.cotTheta() 
00162         << std::endl;
00163 #endif
00164       
00165       if (! crossRange.empty() ) {
00166 #ifdef mydebug_Seed      
00167       (*ss) 
00168         << "\n\t\t !!!!ACCEPTED!!! \n\n";
00169 #endif
00170         if (theMaxElement!=0 && result.size() >= theMaxElement){
00171           result.clear();
00172 #ifdef mydebug_Seed      
00173           edm::LogError("TooManySeeds")<<"number of pairs exceed maximum, no pairs produced";
00174 #endif
00175           delete checkRZ;
00176 
00177 #ifdef mydebug_Seed
00178           std::cout << (*ss).str();
00179 #endif
00180           return;
00181         }
00182         result.push_back( OrderedHitPair( *ih, ohit) );
00183       }
00184     }
00185     delete checkRZ;
00186   }
00187 
00188 #ifdef mydebug_Seed
00189   std::cout << (*ss).str();
00190 #endif
00191 }
00192 
00193 
00194 float HitPairGeneratorFromLayerPairForPhotonConversion::getLayerRadius(const DetLayer& layer){
00195   if (layer.location() == GeomDetEnumerators::barrel){ 
00196     const BarrelDetLayer& bl = static_cast<const BarrelDetLayer&>(layer); 
00197     float rLayer = bl.specificSurface().radius(); 
00198     
00199     // the maximal delta phi will be for the innermost hits
00200     float theThickness = layer.surface().bounds().thickness();
00201     return rLayer + 0.5f*theThickness;
00202   }
00203 
00204   //Fixme
00205   return 0;
00206 }
00207 
00208 float HitPairGeneratorFromLayerPairForPhotonConversion::getLayerZ(const DetLayer& layer){
00209   if (layer.location() == GeomDetEnumerators::endcap){ 
00210     float layerZ = layer.position().z();
00211     float theThickness = layer.surface().bounds().thickness();
00212     float layerZmax = layerZ > 0 ? layerZ+0.5f*theThickness: layerZ-0.5f*theThickness;
00213     return layerZmax;
00214   }else{
00215     //Fixme
00216     return 0;
00217   }
00218 }
00219 bool HitPairGeneratorFromLayerPairForPhotonConversion::checkBoundaries(const DetLayer& layer,const ConversionRegion& convRegion,float maxSearchR, float maxSearchZ){
00220 
00221   if(layer.location() == GeomDetEnumerators::barrel){
00222 
00223     float minZEndCap=130;
00224     if(fabs(convRegion.convPoint().z()) > minZEndCap){
00225 #ifdef mydebug_Seed
00226       (*ss) << "\tthe conversion seems to be in the endcap. Zconv " << convRegion.convPoint().z() << std::endl;
00227       std::cout << (*ss).str();
00228 #endif
00229       return false;
00230     }
00231       
00232     float R=getLayerRadius(layer);
00233     
00234     if(convRegion.convPoint().perp()>R){
00235 #ifdef mydebug_Seed
00236       (*ss) << "\tthis layer is before the conversion : R layer " << R << " [ Rconv " << convRegion.convPoint().perp() << " Zconv " << convRegion.convPoint().z()<< std::endl;
00237       std::cout << (*ss).str();
00238 #endif
00239       return false;
00240     }
00241     
00242     if(R - convRegion.convPoint().perp() > maxSearchR ){
00243 #ifdef mydebug_Seed
00244       (*ss) << "\tthis layer is far from the conversion more than cut " << maxSearchR << " cm. R layer " << R << " [ Rconv " << convRegion.convPoint().perp() << " Zconv " << convRegion.convPoint().z()<< std::endl;
00245       std::cout << (*ss).str();
00246 #endif
00247       return false;
00248     }
00249     
00250   }else if (layer.location() == GeomDetEnumerators::endcap){
00251     
00252     float Z=getLayerZ(layer);
00253     if(
00254        (convRegion.convPoint().z()>0 && convRegion.convPoint().z()>Z)
00255        ||
00256        (convRegion.convPoint().z()<0 && convRegion.convPoint().z()<Z)
00257        ) {
00258 #ifdef mydebug_Seed
00259       (*ss) << "\tthis layer is before the conversion : Z layer " << Z << " [ Rconv " << convRegion.convPoint().perp()<< " Zconv " << convRegion.convPoint().z() << std::endl;
00260       std::cout << (*ss).str();
00261 #endif
00262       return false;
00263     }    
00264 
00265      
00266     if(fabs(Z - convRegion.convPoint().z()) > maxSearchZ ){
00267 #ifdef mydebug_Seed
00268       (*ss) << "\tthis layer is far from the conversion more than cut " << maxSearchZ << " cm. Z layer " << Z << " [ Rconv " << convRegion.convPoint().perp()<< " Zconv " << convRegion.convPoint().z() << std::endl;
00269       std::cout << (*ss).str();
00270 #endif
00271       return false;
00272     }
00273     
00274   }
00275   return true;
00276 }
00277 
00278 bool HitPairGeneratorFromLayerPairForPhotonConversion::getPhiRange(float& Phimin, float& Phimax, const DetLayer& layer, const ConversionRegion &convRegion, const edm::EventSetup& es){
00279    if(layer.location() == GeomDetEnumerators::barrel){
00280      return getPhiRange(Phimin,Phimax,getLayerRadius(layer),convRegion,es);
00281    }else if (layer.location() == GeomDetEnumerators::endcap){
00282      float Z=getLayerZ(layer);
00283      float R=Z/convRegion.cotTheta();
00284      return getPhiRange(Phimin,Phimax,R,convRegion,es); //FIXME
00285    }
00286    return false;
00287 }
00288 
00289 bool HitPairGeneratorFromLayerPairForPhotonConversion::getPhiRange(float& Phimin, float& Phimax, const float& layerR, const ConversionRegion &convRegion, const edm::EventSetup& es){
00290   Phimin = reco::deltaPhi(convRegion.convPoint().phi(),0.);
00291   
00292   float dphi;
00293   float ptmin=0.1;
00294   float DeltaL=layerR-convRegion.convPoint().perp();
00295 
00296   if(DeltaL<0){
00297     Phimin=0;
00298     Phimax=0;
00299     return false;
00300   }
00301 
00302   float theRCurvatureMin = PixelRecoUtilities::bendingRadius(ptmin,es);
00303 
00304   if(theRCurvatureMin<DeltaL)
00305     dphi = atan(DeltaL/layerR);
00306   else
00307     dphi = atan(theRCurvatureMin/layerR * ( 1 - sqrt(1-sqr(DeltaL/theRCurvatureMin)) ) );
00308   
00309   if(convRegion.charge()>0){
00310     Phimax=Phimin;
00311     Phimin=Phimax-dphi;
00312   }else{
00313     Phimax=Phimin+dphi;
00314   }
00315 
00316   //std::cout << dphi << " " << Phimin << " " << Phimax << " " << layerR << " " << DeltaL  << " " << convRegion.convPoint().phi() << " " << convRegion.convPoint().perp()<< std::endl;
00317   return true;
00318 }  
00319 
00320 bool HitPairGeneratorFromLayerPairForPhotonConversion::
00321 checkRZCompatibilityWithSeedTrack(const RecHitsSortedInPhi::Hit& hit,const DetLayer& layer,const ConversionRegion& convRegion){
00322   static const float nSigmaRZ = std::sqrt(12.f);
00323   Range hitCotTheta;
00324   
00325   double sigmaCotTheta = convRegion.errTheta() * (1+convRegion.cotTheta()*convRegion.cotTheta()); //Error Propagation from sigma theta.
00326   Range allowedCotTheta(convRegion.cotTheta()-nSigmaRZ*sigmaCotTheta,convRegion.cotTheta()+nSigmaRZ*sigmaCotTheta);
00327   
00328   double dz        = hit->globalPosition().z()-convRegion.pvtxPoint().z();
00329   double r_reduced = std::sqrt( sqr(hit->globalPosition().x()-convRegion.pvtxPoint().x())+sqr(hit->globalPosition().y()-convRegion.pvtxPoint().y()));
00330 
00331   if (layer.location() == GeomDetEnumerators::barrel){ 
00332     float zErr = nSigmaRZ * hit->errorGlobalZ();
00333     hitCotTheta = Range(getCot(dz-zErr,r_reduced),getCot(dz+zErr,r_reduced));
00334   }else{
00335     float rErr = nSigmaRZ * hit->errorGlobalR();
00336     if(dz>0)  
00337       hitCotTheta = Range(getCot(dz,r_reduced+rErr),getCot(dz,r_reduced-rErr));
00338     else
00339       hitCotTheta = Range(getCot(dz,r_reduced-rErr), getCot(dz,r_reduced+rErr));
00340   }
00341 
00342   Range crossRange = allowedCotTheta.intersection(hitCotTheta);
00343 
00344 #ifdef mydebug_Seed      
00345       (*ss) 
00346         << "\n\t\t cotTheta allowed Range " << allowedCotTheta.min() << " \t, " << allowedCotTheta.max() 
00347         << "\n\t\t hitCotTheta   Range " << hitCotTheta.min()   << " \t, " << hitCotTheta.max() 
00348         << "\n\t\t Cross   Range " << crossRange.min()   << " \t, " << crossRange.max() 
00349         << "\n\t\t the seed track has origin " << convRegion.convPoint() << " \t cotTheta " << convRegion.cotTheta()  
00350         << std::endl;
00351 #endif
00352 
00353   return crossRange.empty();
00354 }
00355 
00356 
00357 double HitPairGeneratorFromLayerPairForPhotonConversion::
00358 getCot(double dz, double dr){
00359   if ( std::abs(dr) > 1.e-4f )    return dz/dr;
00360   else
00361     if(dz>0) return  99999.f; 
00362     else     return -99999.f; 
00363 }
00364