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
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;
00066 if(!checkBoundaries(*theOuterLayer.detLayer(),convRegion,50.,60.)) return;
00067
00068
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
00077
00078
00079 static const float nSigmaRZ = std::sqrt(12.f);
00080
00081 vector<RecHitsSortedInPhi::Hit> innerHits, outerHits;
00082 float outerPhimin, outerPhimax;
00083 float innerPhimin, innerPhimax;
00084
00085
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
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
00103 if(checkRZCompatibilityWithSeedTrack(ohit,*theOuterLayer.detLayer(),convRegion)) continue;
00104
00105
00106
00107
00108
00109
00110
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
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
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
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
00200 float theThickness = layer.surface().bounds().thickness();
00201 return rLayer + 0.5f*theThickness;
00202 }
00203
00204
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
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);
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
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());
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