CMS 3D CMS Logo

/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/RecoTracker/RoadSearchCloudMaker/src/RoadSearchCloudMakerAlgorithm.cc

Go to the documentation of this file.
00001 //
00002 // Package:         RecoTracker/RoadSearchCloudMaker
00003 // Class:           RoadSearchCloudMakerAlgorithm
00004 // 
00005 // Description:     
00006 //                  Road categories determined by outer Seed RecHit
00007 //                      RPhi: outer Seed RecHit in the Barrel
00008 //                      ZPhi: outer Seed RecHit in the Disks
00009 //                  use inner and outer Seed RecHit and BeamSpot to calculate extrapolation
00010 //                      RPhi: phi = phi0 + asin(k r)
00011 //                      ZPhi: phi = phi0 + C z
00012 //                  Loop over RoadSet, access Rings of Road
00013 //                      get average radius of Ring
00014 //                      use extrapolation to calculate phi_ref at average Ring radius
00015 //                      determine window in phi for DetId lookup in the Ring
00016 //                              phi_ref ± phi_window
00017 //                                      PARAMETER: phi_window = pi/24
00018 //                      loop over DetId's in phi window in Ring
00019 //                              two cases (problem of phi = 0/2pi):
00020 //                                      lower window border < upper window border
00021 //                                      upper window border < lower window border
00022 //                  Loop over RecHits of DetId
00023 //                      check compatibility of RecHit with extrapolation (delta phi cut)
00024 //                      single layer sensor treatment
00025 //                              RPhi:
00026 //                                      stereo and barrel single layer sensor
00027 //                                              calculate delta-phi
00028 //                                      disk single layer sensor (r coordinate not well defined)
00029 //                                              calculate delta phi between global positions of maximal strip extension and reference
00030 //                              ZPhi:
00031 //                                      stereo sensor
00032 //                                              calculate delta-phi
00033 //                                      barrel single layer sensor (z coordinate not well defined)
00034 //                                              calculate delta phi between global positions of maximal strip extension and reference
00035 //                                      disk single layer sensor (tilted strips relative to local coordinate system of sensor
00036 //                                              calculate delta phi between global positions of maximal strip extension and reference
00037 //                  Check delta phi cut
00038 //                      cut value can be calculated based on extrapolation and Seed (pT dependent delta phi cut)
00039 //                      currently: constant delta phi cut (PARAMETER)
00040 //                              fill RecHit into Cloud
00041 //                                      do not fill more than 32 RecHits per DetID into cloud (PARAMETER)
00042 //                  first stage of Cloud cleaning cuts:
00043 //                      number of layers with hits in cloud (PARAMETER)
00044 //                      number of layers with no hits in cloud (PARAMETER)
00045 //                      number of consecutive layers with no hits in cloud (PARAMETER)
00046 //
00047 // Original Author: Oliver Gutsche, gutsche@fnal.gov
00048 // Created:         Sat Jan 14 22:00:00 UTC 2006
00049 //
00050 // $Author: gpetrucc $
00051 // $Date: 2010/12/13 17:33:48 $
00052 // $Revision: 1.57 $
00053 //
00054 
00055 #include <vector>
00056 #include <iostream>
00057 #include <cmath>
00058 
00059 #include "RecoTracker/RoadSearchCloudMaker/interface/RoadSearchCloudMakerAlgorithm.h"
00060 
00061 #include "FWCore/ServiceRegistry/interface/Service.h"
00062 
00063 #include "DataFormats/Common/interface/Handle.h"
00064 #include "FWCore/Framework/interface/ESHandle.h"
00065 #include "FWCore/Framework/interface/EventSetup.h"
00066 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00067 
00068 #include "DataFormats/Common/interface/Ref.h"
00069 #include "DataFormats/RoadSearchSeed/interface/RoadSearchSeedCollection.h"
00070 #include "DataFormats/RoadSearchCloud/interface/RoadSearchCloud.h"
00071 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
00072 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
00073 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
00074 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
00075 #include "DataFormats/SiStripDetId/interface/TIBDetId.h"
00076 #include "DataFormats/SiStripDetId/interface/TOBDetId.h"
00077 #include "DataFormats/SiStripDetId/interface/TIDDetId.h"
00078 #include "DataFormats/SiStripDetId/interface/TECDetId.h"
00079 #include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
00080 #include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
00081 
00082 #include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
00083 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00084 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00085 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
00086 #include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h"
00087 #include "Geometry/CommonTopologies/interface/RectangularStripTopology.h"
00088 
00089 #include "DataFormats/SiPixelDigi/interface/PixelDigi.h"
00090 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
00091 
00092 #include "TrackingTools/RoadSearchHitAccess/interface/RoadSearchDetIdHelper.h"
00093 
00094 #include "RecoTracker/RoadMapRecord/interface/RoadMapRecord.h"
00095 
00096 #include "RecoLocalTracker/SiStripRecHitConverter/interface/SiStripRecHitMatcher.h"
00097 #include "Geometry/TrackerGeometryBuilder/interface/GluedGeomDet.h"
00098 
00099 using namespace std;
00100 
00101 double RoadSearchCloudMakerAlgorithm::epsilon      =   0.000000001;
00102 
00103 RoadSearchCloudMakerAlgorithm::RoadSearchCloudMakerAlgorithm(const edm::ParameterSet& conf) : conf_(conf) { 
00104   recHitVectorClass.setMode(DetHitAccess::standard);    
00105   recHitVectorClass.use_rphiRecHits(conf_.getParameter<bool>("UseRphiRecHits"));
00106   recHitVectorClass.use_stereoRecHits(conf_.getParameter<bool>("UseStereoRecHits"));
00107   
00108   
00109   theRPhiRoadSize =  conf_.getParameter<double>("RPhiRoadSize");
00110   theZPhiRoadSize =  conf_.getParameter<double>("ZPhiRoadSize");
00111   UsePixels = conf_.getParameter<bool>("UsePixelsinRS");
00112   NoFieldCosmic = conf_.getParameter<bool>("StraightLineNoBeamSpotCloud");
00113   theMinimumHalfRoad = conf_.getParameter<double>("MinimumHalfRoad");
00114   
00115   maxDetHitsInCloudPerDetId = conf_.getParameter<unsigned int>("MaxDetHitsInCloudPerDetId");
00116   minFractionOfUsedLayersPerCloud = conf_.getParameter<double>("MinimalFractionOfUsedLayersPerCloud");
00117   maxFractionOfMissedLayersPerCloud = conf_.getParameter<double>("MaximalFractionOfMissedLayersPerCloud");
00118   maxFractionOfConsecutiveMissedLayersPerCloud = conf_.getParameter<double>("MaximalFractionOfConsecutiveMissedLayersPerCloud");
00119   increaseMaxNumberOfConsecutiveMissedLayersPerCloud = conf_.getParameter<unsigned int>("IncreaseMaxNumberOfConsecutiveMissedLayersPerCloud");
00120   increaseMaxNumberOfMissedLayersPerCloud = conf_.getParameter<unsigned int>("IncreaseMaxNumberOfMissedLayersPerCloud");
00121   
00122   doCleaning_ = conf.getParameter<bool>("DoCloudCleaning");
00123   mergingFraction_ = conf.getParameter<double>("MergingFraction");
00124   maxRecHitsInCloud_ = (unsigned int)conf.getParameter<int>("MaxRecHitsInCloud");
00125   scalefactorRoadSeedWindow_ = conf.getParameter<double>("scalefactorRoadSeedWindow");
00126   
00127   roadsLabel_ = conf.getParameter<std::string>("RoadsLabel");
00128 
00129 }
00130 
00131 RoadSearchCloudMakerAlgorithm::~RoadSearchCloudMakerAlgorithm() {
00132 
00133 }
00134 
00135 double RoadSearchCloudMakerAlgorithm::map_phi(double phi) {
00136   // map phi to [0,2pi]
00137   double result = phi;
00138   if ( result < -1.0*Geom::twoPi()) result = result + Geom::twoPi();
00139   if ( result < 0)                 result = Geom::twoPi() + result;
00140   if ( result > Geom::twoPi())     result = result - Geom::twoPi();
00141   return result;
00142 }
00143 
00144 double RoadSearchCloudMakerAlgorithm::map_phi2(double phi) {
00145   // map phi to [-pi,pi]
00146   double result = phi;
00147   if ( result < 1.0*Geom::pi() ) result = result + Geom::twoPi();
00148   if ( result >= Geom::pi())  result = result - Geom::twoPi();
00149   return result;
00150 }
00151 
00152 void RoadSearchCloudMakerAlgorithm::run(edm::Handle<RoadSearchSeedCollection> input,
00153                                         const SiStripRecHit2DCollection* rphiRecHits,
00154                                         const SiStripRecHit2DCollection* stereoRecHits,
00155                                         const SiStripMatchedRecHit2DCollection* matchedRecHits,
00156                                         const SiPixelRecHitCollection *pixRecHits,
00157                                         const edm::EventSetup& es,
00158                                         RoadSearchCloudCollection &output)
00159 {  
00160   // intermediate arrays for storing clouds for cleaning
00161   const int nphibin = 24;
00162   const int netabin = 24;
00163   RoadSearchCloudCollection CloudArray[nphibin][netabin];
00164   
00165   // get roads
00166   edm::ESHandle<Roads> roads;
00167   es.get<RoadMapRecord>().get(roadsLabel_, roads);
00168   
00169   // get RoadSearchSeed collection
00170   const RoadSearchSeedCollection* inputSeeds = input.product();
00171   
00172   // set collections for general hit access method
00173   recHitVectorClass.setCollections(rphiRecHits,stereoRecHits,matchedRecHits,pixRecHits);
00174   recHitVectorClass.setMode(DetHitAccess::standard);
00175   
00176   // get tracker geometry
00177   edm::ESHandle<TrackerGeometry> tracker;
00178   es.get<TrackerDigiGeometryRecord>().get(tracker);
00179   
00180   // get hit matcher
00181   SiStripRecHitMatcher* theHitMatcher = new SiStripRecHitMatcher(3.0);
00182 
00183   edm::LogInfo("RoadSearch") << "Found " << inputSeeds->size() << " input seeds.";   
00184   // loop over seeds
00185   for ( RoadSearchSeedCollection::const_iterator seed = inputSeeds->begin(); seed != inputSeeds->end(); ++seed) {
00186     
00187     const Roads::RoadSeed *roadSeed = seed->getSeed();
00188       
00189     if ( roadSeed == 0 ) {
00190       edm::LogWarning("RoadSearch") << "RoadSeed could not be resolved from RoadSearchSeed hits, discard seed!";
00191     } else {
00192         
00193       Roads::type roadType = roads->getRoadType(roadSeed);
00194       if (NoFieldCosmic) roadType = Roads::RPhi;
00195         
00196       // fixme: from here on, calculate with 1st and 3rd seed hit (inner and outer of initial circle)
00197       // fixme: adapt to new seed structure
00198         
00199       // get global positions of the hits, calculate before Road lookup to be used
00200       const TrackingRecHit* innerSeedRingHit = (*(seed->begin()));
00201       const TrackingRecHit* outerSeedRingHit = (*(seed->end() - 1));
00202         
00203       GlobalPoint innerSeedHitGlobalPosition = tracker->idToDet(innerSeedRingHit->geographicalId())->surface().toGlobal(innerSeedRingHit->localPosition());
00204       GlobalPoint outerSeedHitGlobalPosition = tracker->idToDet(outerSeedRingHit->geographicalId())->surface().toGlobal(outerSeedRingHit->localPosition());
00205  
00206       LogDebug("RoadSearch") << "Seed # " <<seed-inputSeeds->begin() << " inner hit (x/y/z): "
00207                              << innerSeedHitGlobalPosition.x() << " / "
00208                              << innerSeedHitGlobalPosition.y() << " / "
00209                              << innerSeedHitGlobalPosition.z();
00210       LogDebug("RoadSearch") << "Seed # " <<seed-inputSeeds->begin() << " outer hit (x/y/z): "
00211                              << outerSeedHitGlobalPosition.x() << " / "
00212                              << outerSeedHitGlobalPosition.y() << " / "
00213                              << outerSeedHitGlobalPosition.z();
00214       
00215       LogDebug("RoadSearch") << "Seed # " <<seed-inputSeeds->begin() << " inner hit (r/phi): "
00216                              << innerSeedHitGlobalPosition.perp() << " / "
00217                              << innerSeedHitGlobalPosition.phi();
00218       LogDebug("RoadSearch") << "Seed # " <<seed-inputSeeds->begin() << " outer hit (r/phi): "
00219                              << outerSeedHitGlobalPosition.perp() << " / "
00220                              << outerSeedHitGlobalPosition.phi();
00221 
00222 
00223       // extrapolation parameters, phio: [0,2pi]
00224       double d0 = 0.0;
00225       double phi0 = -99.;
00226       double k0   = -99999999.99;
00227       
00228       double phi1 = -99.;
00229       double k1   = -99999999.99;
00230       // get bins of eta and phi of outer seed hit;
00231       
00232       double outer_phi = map_phi(outerSeedHitGlobalPosition.phi());
00233       double outer_eta = outerSeedHitGlobalPosition.eta();
00234       
00235       int phibin = (int)(nphibin*(outer_phi/(2*Geom::pi())));
00236       int etabin = (int)(netabin*(outer_eta+3.0)/6.0);
00237         
00238       // calculate phi0 and k0 dependent on RoadType
00239       if ( roadType == Roads::RPhi ) {
00240         double dr = outerSeedHitGlobalPosition.perp() - innerSeedHitGlobalPosition.perp();
00241         const double dr_min = 1; // cm
00242         if ( dr < dr_min ) {
00243           edm::LogWarning("RoadSearch") << "RPhi road: seed Hits distance smaller than 1 cm, do not consider this seed.";
00244         } else {
00245           // calculate r-phi extrapolation: phi = phi0 + asin(k0 r)
00246           double det = innerSeedHitGlobalPosition.x() * outerSeedHitGlobalPosition.y() - innerSeedHitGlobalPosition.y() * outerSeedHitGlobalPosition.x();
00247           if ( det == 0 ) {
00248             edm::LogWarning("RoadSearch") << "RPhi road: 'det' == 0, do not consider this seed.";
00249           } else {
00250             double x0=0.0; double y0=0.0;
00251             double innerx=innerSeedHitGlobalPosition.x();
00252             double innery=innerSeedHitGlobalPosition.y();
00253             double outerx=outerSeedHitGlobalPosition.x();
00254             double outery=outerSeedHitGlobalPosition.y();
00255 
00256             if (NoFieldCosmic){
00257               phi0=atan2(outery-innery,outerx-innerx);
00258               double alpha=atan2(innery,innerx);
00259               double d1=sqrt(innerx*innerx+innery*innery);
00260               d0=d1*sin(alpha-phi0); x0=-d0*sin(phi0); y0=d0*cos(phi0); k0=0.0;
00261             }else{
00262               makecircle(innerx,innery,outerx,outery,x0,y0);
00263               phi0 = phi0h;
00264               k0 = omegah;
00265             }              
00266             LogDebug("RoadSearch") << "Seed # " <<seed-inputSeeds->begin() << " trajectory parameters: d0 = "<< d0 << " phi0 = " << phi0;
00267           }
00268         }
00269       } else {
00270         double dz = outerSeedHitGlobalPosition.z() - innerSeedHitGlobalPosition.z();
00271         const double dz_min = 1.e-6; // cm;
00272         if ( std::abs(dz) < dz_min ) {
00273           edm::LogWarning("RoadSearch") << "ZPhi road: seed Hits are less than .01 microns away in z, do not consider this seed.";
00274         } else {
00275           // calculate z-phi extrapolation: phi = phi0 + k0 z
00276           k0 = map_phi2(outerSeedHitGlobalPosition.phi() - innerSeedHitGlobalPosition.phi()) / dz;
00277           phi0 =  map_phi(innerSeedHitGlobalPosition.phi() - k0 * innerSeedHitGlobalPosition.z());
00278           
00279           // get approx pt for use in correcting matched hits
00280           makecircle(innerSeedHitGlobalPosition.x(),innerSeedHitGlobalPosition.y(),
00281                      outerSeedHitGlobalPosition.x(),outerSeedHitGlobalPosition.y(),
00282                      0.0,0.0); // x0,y0 = 0.0 for now
00283           phi1 = phi0h;
00284           k1 = omegah;
00285         }
00286       }
00287         
00288       // continue if valid extrapolation parameters have been found
00289       if ( (phi0 != -99.) && (k0 != -99999999.99) ) {
00290         const Roads::RoadSet *roadSet = seed->getSet();
00291 
00292         // create cloud
00293         RoadSearchCloud cloud;
00294           
00295         bool firstHitFound = false;
00296         unsigned int layerCounter = 0;
00297         unsigned int usedLayers = 0;
00298         unsigned int missedLayers = 0;
00299         unsigned int consecutiveMissedLayers = 0;
00300 
00301         unsigned int totalLayers = roadSet->size();
00302 
00303         // caluclate minNumberOfUsedLayersPerCloud, maxNumberOfMissedLayersPerCloud and maxNumberOfConsecutiveMissedLayersPerCloud 
00304         // by rounding to integer minFractionOfUsedLayersPerCloud. maxFractionOfMissedLayersPerCloud and maxFractionOfConsecutiveMissedLayersPerCloud
00305         unsigned int minNumberOfUsedLayersPerCloud = static_cast<unsigned int>(totalLayers * minFractionOfUsedLayersPerCloud + 0.5);
00306         if (minNumberOfUsedLayersPerCloud < 3) minNumberOfUsedLayersPerCloud = 3;
00307         unsigned int maxNumberOfMissedLayersPerCloud = static_cast<unsigned int>(totalLayers * maxFractionOfMissedLayersPerCloud + 0.5);
00308         unsigned int maxNumberOfConsecutiveMissedLayersPerCloud = static_cast<unsigned int>(totalLayers * maxFractionOfConsecutiveMissedLayersPerCloud + 0.5);
00309 
00310         // increase consecutive layer cuts between 0.9 and 1.5
00311         if (std::abs(outer_eta) > 0.9 && std::abs(outer_eta) < 1.5) {
00312           maxNumberOfConsecutiveMissedLayersPerCloud += increaseMaxNumberOfConsecutiveMissedLayersPerCloud;
00313           maxNumberOfMissedLayersPerCloud += increaseMaxNumberOfMissedLayersPerCloud;
00314         }
00315           
00316         for ( Roads::RoadSet::const_iterator roadSetVector = roadSet->begin();
00317               roadSetVector != roadSet->end();
00318               ++roadSetVector ) {
00319             
00320           ++layerCounter;
00321           unsigned int usedHitsInThisLayer = 0;
00322           bool intersectsLayer = false;
00323             
00324           for ( std::vector<const Ring*>::const_iterator ring = roadSetVector->begin(); ring != roadSetVector->end(); ++ring ) {
00325               
00326             // calculate phi-range for lookup of DetId's in Rings of RoadSet
00327             // calculate phi at radius of Ring using extrapolation, Ring radius average of Ring.rmin() and Ring.rmax()
00328             double ringRadius = (*ring)->getrmin() + ((*ring)->getrmax()-(*ring)->getrmin())/2;
00329             double ringZ      = (*ring)->getzmin() + ((*ring)->getzmax()-(*ring)->getzmin())/2;
00330             double ringPhi = 0.0;
00331             if ( roadType == Roads::RPhi ) {
00332               ringPhi = phiFromExtrapolation(d0,phi0,k0,ringRadius,roadType);
00333             } else {
00334               ringPhi = phiFromExtrapolation(d0,phi0,k0,ringZ,roadType);
00335             }
00336             if (ringPhi == -99) continue;
00337             intersectsLayer = true;
00338 
00339             LogDebug("RoadSearch") << "Seed # " <<seed-inputSeeds->begin() << " testing ring at R = " << ringRadius
00340                                    << " Z = " << ringZ << " ringPhi = " << ringPhi;
00341               
00342             int nDetIds = (*ring)->getNumDetIds();
00343             double theHalfRoad = theMinimumHalfRoad*(2.0*Geom::pi())/((double)nDetIds);
00344             // calculate range in phi around ringPhi
00345             double upperPhiRangeBorder = map_phi2(ringPhi + theHalfRoad);
00346             double lowerPhiRangeBorder = map_phi2(ringPhi - theHalfRoad);
00347 
00348             if ( lowerPhiRangeBorder <= upperPhiRangeBorder ) {
00349                 
00350               for ( Ring::const_iterator detid = (*ring)->lower_bound(lowerPhiRangeBorder); detid != (*ring)->upper_bound(upperPhiRangeBorder); ++detid) {
00351                 usedHitsInThisLayer += FillRecHitsIntoCloudGeneral(detid->second,d0,phi0,k0,phi1,k1,roadType,ringPhi,
00352                                                                    tracker.product(),theHitMatcher,cloud);
00353               }
00354                 
00355             } else {
00356               for ( Ring::const_iterator detid = (*ring)->lower_bound(lowerPhiRangeBorder); detid != (*ring)->end(); ++detid) {
00357                 usedHitsInThisLayer += FillRecHitsIntoCloudGeneral(detid->second,d0,phi0,k0,phi1,k1,roadType,ringPhi,
00358                                                                    tracker.product(),theHitMatcher,cloud);
00359               }
00360                 
00361               for ( Ring::const_iterator detid = (*ring)->begin(); detid != (*ring)->upper_bound(upperPhiRangeBorder); ++detid) {
00362                 usedHitsInThisLayer += FillRecHitsIntoCloudGeneral(detid->second,d0,phi0,k0,phi1,k1,roadType,ringPhi,
00363                                                                    tracker.product(),theHitMatcher,cloud);
00364               }
00365             }
00366           LogDebug("RoadSearch") << "Seed # " <<seed-inputSeeds->begin() << " now has " << usedHitsInThisLayer << "  hits in ring at R = " << ringRadius
00367                                  << " Z = " << ringZ << " ringPhi = " << ringPhi;
00368           }
00369             
00370           if ( !firstHitFound ) {
00371             if ( usedHitsInThisLayer > 0 ) {
00372 
00373               firstHitFound = true;
00374 
00375               // reset totalLayers according to first layer with hit
00376               totalLayers = roadSet->size() - layerCounter + 1;
00377 
00378               // re-caluclate minNumberOfUsedLayersPerCloud, maxNumberOfMissedLayersPerCloud and maxNumberOfConsecutiveMissedLayersPerCloud 
00379               // by rounding to integer minFractionOfUsedLayersPerCloud. maxFractionOfMissedLayersPerCloud and maxFractionOfConsecutiveMissedLayersPerCloud
00380               minNumberOfUsedLayersPerCloud = static_cast<unsigned int>(totalLayers * minFractionOfUsedLayersPerCloud + 0.5);
00381               if (minNumberOfUsedLayersPerCloud < 3) minNumberOfUsedLayersPerCloud = 3;
00382               maxNumberOfMissedLayersPerCloud = static_cast<unsigned int>(totalLayers * maxFractionOfMissedLayersPerCloud + 0.5);
00383               maxNumberOfConsecutiveMissedLayersPerCloud = static_cast<unsigned int>(totalLayers * maxFractionOfConsecutiveMissedLayersPerCloud + 0.5);
00384 
00385               // increase consecutive layer cuts between 0.9 and 1.5
00386               if (std::abs(outer_eta) > 0.9 && std::abs(outer_eta) < 1.5) {
00387                 maxNumberOfConsecutiveMissedLayersPerCloud += increaseMaxNumberOfConsecutiveMissedLayersPerCloud;
00388                 maxNumberOfMissedLayersPerCloud += increaseMaxNumberOfMissedLayersPerCloud;
00389               }
00390 
00391               ++usedLayers;
00392               consecutiveMissedLayers = 0;
00393 
00394             }
00395           } else {
00396             if (intersectsLayer){
00397               if ( usedHitsInThisLayer > 0 ) {
00398                 ++usedLayers;
00399                 consecutiveMissedLayers = 0;
00400               } else {
00401                 ++ missedLayers;
00402                 ++consecutiveMissedLayers;
00403               }
00404             }
00405             LogDebug("RoadSearch") << "Seed # "<<seed-inputSeeds->begin() << " Layer info: " 
00406                                    << " totalLayers: " << totalLayers 
00407                                    << " usedLayers: " << usedLayers 
00408                                    << " missedLayers: " << missedLayers
00409                                    << " consecutiveMissedLayers: " << consecutiveMissedLayers;
00410 
00411             // break condition, hole larger than maxNumberOfConsecutiveMissedLayersPerCloud
00412             if ( consecutiveMissedLayers > maxNumberOfConsecutiveMissedLayersPerCloud ) {
00413               LogDebug("RoadSearch") << "BREAK: seed # "<<seed-inputSeeds->begin() 
00414                                      << " More than " << maxNumberOfConsecutiveMissedLayersPerCloud << " missed consecutive layers!";
00415               break;
00416             }
00417 
00418             // break condition, already  missed too many layers
00419             if ( missedLayers > maxNumberOfMissedLayersPerCloud ) {
00420               LogDebug("RoadSearch") << "BREAK: seed # "<<seed-inputSeeds->begin() 
00421                                      << " More than " << maxNumberOfMissedLayersPerCloud << " missed layers!";
00422               break;
00423             }
00424 
00425             // break condition, cannot satisfy minimal number of used layers
00426             if ( totalLayers-missedLayers < minNumberOfUsedLayersPerCloud ) {
00427               LogDebug("RoadSearch") << "BREAK: seed # "<<seed-inputSeeds->begin() 
00428                                      << " Cannot satisfy at least " << minNumberOfUsedLayersPerCloud << " used layers!";
00429               break;
00430             }
00431           }       
00432             
00433         }
00434 
00435         if ( consecutiveMissedLayers <= maxNumberOfConsecutiveMissedLayersPerCloud ) {
00436           if ( usedLayers >= minNumberOfUsedLayersPerCloud ) {
00437             if ( missedLayers <= maxNumberOfMissedLayersPerCloud ) {
00438 
00439               CloudArray[phibin][etabin].push_back(cloud);
00440 
00441               if ( roadType == Roads::RPhi ){ 
00442                 LogDebug("RoadSearch") << "This r-phi seed # "<<seed-inputSeeds->begin()
00443                                        <<" yields a cloud with " <<cloud.size() <<" hits on " << usedLayers << " layers out of " << totalLayers;
00444               } else {
00445                 LogDebug("RoadSearch") << "This z-phi seed # "<<seed-inputSeeds->begin()
00446                                        <<" yields a cloud with " <<cloud.size() <<" hits on " << usedLayers << " layers out of " << totalLayers;
00447               }
00448             } else {
00449               LogDebug("RoadSearch") << "Missed layers: " << missedLayers << " More than " << maxNumberOfMissedLayersPerCloud << " missed layers!";
00450               if ( roadType == Roads::RPhi ){ 
00451                 LogDebug("RoadSearch") << "This r-phi seed # "<<seed-inputSeeds->begin() <<" yields no clouds";
00452               } else {
00453                 LogDebug("RoadSearch") << "This z-phi seed # "<<seed-inputSeeds->begin() <<" yields no clouds";
00454               }
00455             }
00456           }
00457           else {
00458             LogDebug("RoadSearch") << "Seed # "<<seed-inputSeeds->begin() <<" fails: used layers = " << usedLayers << " < " << minNumberOfUsedLayersPerCloud;
00459           }
00460         }
00461         else {
00462           LogDebug("RoadSearch") << "Seed # "<<seed-inputSeeds->begin() <<" fails: consecutive missed layers = " << consecutiveMissedLayers << " > " << maxNumberOfConsecutiveMissedLayersPerCloud;
00463         }
00464       }
00465     }
00466   }
00467   
00468   // Loop for initial cleaning
00469   for (int iphi=0; iphi<nphibin; ++iphi){
00470     for (int ieta=0; ieta<netabin; ++ieta){
00471       if (!CloudArray[iphi][ieta].empty()) {
00472         if (doCleaning_){
00473           RoadSearchCloudCollection temp = Clean(&CloudArray[iphi][ieta]);
00474           for ( RoadSearchCloudCollection::const_iterator ic = temp.begin(); ic!=temp.end(); ++ic)
00475             output.push_back(*ic);
00476         }
00477         else 
00478           for ( RoadSearchCloudCollection::const_iterator ic = CloudArray[iphi][ieta].begin(); 
00479                 ic!=CloudArray[iphi][ieta].end(); ++ic)
00480             output.push_back(*ic);
00481       }
00482     }
00483   }
00484 
00485   delete theHitMatcher;
00486   edm::LogInfo("RoadSearch") << "Found " << output.size() << " clouds."; 
00487   for ( RoadSearchCloudCollection::const_iterator ic = output.begin(); ic!=output.end(); ++ic)
00488     edm::LogInfo("RoadSearch") << "    Cloud " << ic-output.begin()<< " has " << ic->size() << " hits."; 
00489   
00490 }
00491 
00492 unsigned int RoadSearchCloudMakerAlgorithm::FillRecHitsIntoCloudGeneral(DetId id, double d0, double phi0, double k0, double phi1, double k1, 
00493                                                                         Roads::type roadType, double ringPhi,
00494                                                                         const TrackerGeometry *tracker, const SiStripRecHitMatcher* theHitMatcher,
00495                                                                         RoadSearchCloud &cloud) {
00496   
00497   unsigned int usedRecHits = 0;
00498   
00499   bool double_ring_layer = !isSingleLayer(id);
00500   
00501   std::vector<TrackingRecHit*> recHitVector = recHitVectorClass.getHitVector(&id);
00502   
00503   for ( std::vector<TrackingRecHit*>::const_iterator recHitIterator = recHitVector.begin(); recHitIterator != recHitVector.end(); ++recHitIterator) {
00504     
00505     if (      (unsigned int)id.subdetId() == StripSubdetector::TIB 
00506               || (unsigned int)id.subdetId() == StripSubdetector::TOB 
00507               || (unsigned int)id.subdetId() == StripSubdetector::TID 
00508               || (unsigned int)id.subdetId() == StripSubdetector::TEC ) {
00509       
00510       const SiStripRecHit2D *recHit = (SiStripRecHit2D*)(*recHitIterator);
00511       DetId hitId = recHit->geographicalId();
00512 
00513       GlobalPoint ghit = tracker->idToDet(recHit->geographicalId())->surface().toGlobal(recHit->localPosition());
00514       //LogDebug("RoadSearch") << "    Testing hit at (x/y/z): " << ghit.x() << " / " << ghit.y() << " / " << ghit.z();
00515       LogDebug("RoadSearch") << "    Testing hit at (x/y/z): " 
00516                              << tracker->idToDet(recHit->geographicalId())->surface().toGlobal(recHit->localPosition()).x() << " / " 
00517                              << tracker->idToDet(recHit->geographicalId())->surface().toGlobal(recHit->localPosition()).y() << " / " 
00518                              << tracker->idToDet(recHit->geographicalId())->surface().toGlobal(recHit->localPosition()).z();
00519 
00520       if ( roadType == Roads::RPhi ) {
00521         if (double_ring_layer && isSingleLayer(hitId)) {
00522           //
00523           //  This is where the barrel stereoRecHits end up for Roads::RPhi
00524           //
00525           
00526           // Adjust matched hit for track angle
00527 
00528           const GluedGeomDet *theGluedDet = dynamic_cast<const GluedGeomDet*>(tracker->idToDet(hitId));
00529 
00530           SiStripMatchedRecHit2D* theCorrectedHit = CorrectMatchedHit(*recHitIterator,theGluedDet,
00531                                                                         tracker, theHitMatcher,
00532                                                                         k0, phi0);
00533           if (theCorrectedHit != 0){
00534 
00535             GlobalPoint ghit = tracker->idToDet(theCorrectedHit->geographicalId())->surface().toGlobal(theCorrectedHit->localPosition());
00536             double hitRadius = sqrt(ghit.x()*ghit.x()+ghit.y()*ghit.y());
00537             double hitphi = map_phi(ghit.phi());
00538             double phi = phiFromExtrapolation(d0,phi0,k0,hitRadius,roadType);
00539             
00540             float dp = (hitphi-phi);
00541             float dx = hitRadius*tan(dp);
00542             
00543             LogDebug("RoadSearch") << "   Hit phi = " << hitphi << " expected phi = " << phi
00544                                        <<"  dx = " << dx << " for dxMax = " << phiMax(roadType,phi0,k0);
00545 
00546             // switch cut to dx instead of dphi
00547             if ( std::abs(dx) < phiMax(roadType,phi0,k0) ) {
00548               if ((usedRecHits < maxDetHitsInCloudPerDetId) && (cloud.size() < maxRecHitsInCloud_)) {
00549                 cloud.addHit(recHit);
00550                 ++usedRecHits;
00551               }
00552             }
00553             delete theCorrectedHit;
00554           }
00555         } else { // Single layer hits here
00556           if ( isBarrelSensor(hitId) ) {
00557             //
00558             //  This is where the barrel rphiRecHits end up for Roads::RPhi
00559             //
00560             GlobalPoint ghit = tracker->idToDet(recHit->geographicalId())->surface().toGlobal(recHit->localPosition());
00561             double hitRadius = sqrt(ghit.x()*ghit.x()+ghit.y()*ghit.y());
00562             double hitphi = map_phi(ghit.phi());
00563             double phi = phiFromExtrapolation(d0,phi0,k0,hitRadius,roadType);
00564               
00565             float dp = (hitphi-phi);
00566             float dx = hitRadius*tan(dp);
00567             LogDebug("RoadSearch") << "   Hit phi = " << hitphi << " expected phi = " << phi
00568                                        <<"  dx = " << dx << " for dxMax = " << phiMax(roadType,phi0,k0);
00569             // switch cut to dx instead of dphi
00570             if ( std::abs(dx) < phiMax(roadType,phi0,k0) ) {
00571               if ((usedRecHits < maxDetHitsInCloudPerDetId) && (cloud.size() < maxRecHitsInCloud_)) {
00572                 cloud.addHit(recHit);
00573                 ++usedRecHits;
00574               }
00575             }
00576           } 
00577           else {
00578 
00579             LocalPoint hit = recHit->localPosition();
00580             const StripTopology *topology = dynamic_cast<const StripTopology*>(&(tracker->idToDetUnit(hitId)->topology()));
00581             double stripAngle = topology->stripAngle(topology->strip(hit));
00582             double stripLength = topology->localStripLength(hit);
00583 
00584             LocalPoint innerHitLocal(hit.x()+stripLength/2*std::sin(stripAngle),hit.y()-stripLength/2*std::cos(stripAngle),0);
00585             LocalPoint outerHitLocal(hit.x()-stripLength/2*std::sin(stripAngle),hit.y()+stripLength/2*std::cos(stripAngle),0);
00586 
00587             double innerRadius = tracker->idToDetUnit(hitId)->surface().toGlobal(innerHitLocal).perp(); 
00588             double outerRadius = tracker->idToDetUnit(hitId)->surface().toGlobal(outerHitLocal).perp();
00589             double innerExtrapolatedPhi = phiFromExtrapolation(d0,phi0,k0,innerRadius,roadType);
00590             double outerExtrapolatedPhi = phiFromExtrapolation(d0,phi0,k0,outerRadius,roadType);
00591 
00592             GlobalPoint innerHitGlobal =tracker->idToDetUnit(hitId)->surface().toGlobal(innerHitLocal); 
00593             GlobalPoint outerHitGlobal =tracker->idToDetUnit(hitId)->surface().toGlobal(outerHitLocal); 
00594 
00595             GlobalPoint innerRoadGlobal(GlobalPoint::Cylindrical(innerRadius,innerExtrapolatedPhi,
00596                                                                  tracker->idToDetUnit(hitId)->surface().toGlobal(hit).z()));
00597             GlobalPoint outerRoadGlobal(GlobalPoint::Cylindrical(outerRadius,outerExtrapolatedPhi,
00598                                                                  tracker->idToDetUnit(hitId)->surface().toGlobal(hit).z()));
00599 
00600             LocalPoint innerRoadLocal = tracker->idToDetUnit(hitId)->surface().toLocal(innerRoadGlobal);
00601             LocalPoint outerRoadLocal = tracker->idToDetUnit(hitId)->surface().toLocal(outerRoadGlobal);
00602 
00603             double dxinter = CheckXYIntersection(innerHitLocal, outerHitLocal, 
00604                                                innerRoadLocal, outerRoadLocal);
00605 
00606             LogDebug("RoadSearch") << " Hit phi inner = " << innerHitGlobal.phi() << " and outer = " << outerHitGlobal.phi()
00607                                    << " expected inner phi = " << innerExtrapolatedPhi
00608                                    << " and outer phi = "      << outerExtrapolatedPhi
00609                                    <<"  dx = " << dxinter << " for dxMax = " << phiMax(roadType,phi0,k0);
00610 
00611             if ( fabs(dxinter) < phiMax(roadType,phi0,k0)) {
00612               //
00613               //  This is where the disk rphiRecHits end up for Roads::ZPhi
00614               //
00615               if ((usedRecHits < maxDetHitsInCloudPerDetId) && (cloud.size() < maxRecHitsInCloud_)) {
00616                 cloud.addHit(recHit);
00617                 ++usedRecHits;
00618               }
00619             }
00620             //else
00621               //std::cout<< " ===>>> HIT FAILS !!! " << std::endl;
00622           }
00623         } 
00624       } else {
00625         //
00626         // roadType == Roads::ZPhi
00627         //
00628         if (double_ring_layer && isSingleLayer(hitId)) {
00629 
00630           // Adjust matched hit for track angle
00631 
00632           //const SiStripMatchedRecHit2D *theRH = dynamic_cast<SiStripMatchedRecHit2D*>(*recHitIterator);
00633           const GluedGeomDet *theGluedDet = dynamic_cast<const GluedGeomDet*>(tracker->idToDet(hitId));
00634 
00635           SiStripMatchedRecHit2D* theCorrectedHit = CorrectMatchedHit(*recHitIterator,theGluedDet,
00636                                                                         tracker, theHitMatcher,
00637                                                                         k1, phi1);
00638           if (theCorrectedHit != 0){
00639 
00640             GlobalPoint ghit = tracker->idToDet(theCorrectedHit->geographicalId())->surface().toGlobal(theCorrectedHit->localPosition());
00641             double hitphi = map_phi(ghit.phi());
00642             double hitZ = ghit.z();
00643             double phi = phiFromExtrapolation(d0,phi0,k0,hitZ,roadType);
00644             
00645             float dp = hitphi-phi;
00646             float dx = hitZ*tan(dp);
00647               
00648             //
00649             //  This is where the disk stereoRecHits end up for Roads::ZPhi
00650             //
00651             if ( std::abs(dx) < phiMax(roadType,phi0,k1)) {
00652               if ((usedRecHits < maxDetHitsInCloudPerDetId) && (cloud.size() < maxRecHitsInCloud_)) {
00653                 cloud.addHit(recHit);
00654                 ++usedRecHits;
00655               }
00656             }
00657             delete theCorrectedHit;
00658           }
00659         } else { // Single layer hits here
00660           if ( isBarrelSensor(hitId) ) {
00661             //
00662             //  This is where the barrel (???) rphiRecHits end up for Roads::ZPhi
00663             //
00664             LocalPoint hit = recHit->localPosition();
00665             const StripTopology *topology = dynamic_cast<const StripTopology*>(&(tracker->idToDetUnit(hitId)->topology()));
00666             double stripAngle = topology->stripAngle(topology->strip(hit));
00667             double stripLength = topology->localStripLength(hit);
00668 
00669             //if (stripAngle!=0) std::cout<<"HEY, WE FOUND A HIT ON A STEREO MODULE!!!" << std::endl;
00670             // new method
00671             LocalPoint innerHitLocal(hit.x()+stripLength/2*std::sin(stripAngle),hit.y()-stripLength/2*std::cos(stripAngle),0);
00672             LocalPoint outerHitLocal(hit.x()-stripLength/2*std::sin(stripAngle),hit.y()+stripLength/2*std::cos(stripAngle),0);
00673             double innerHitPhi = map_phi(tracker->idToDetUnit(hitId)->surface().toGlobal(innerHitLocal).phi()); 
00674             double outerHitPhi = map_phi(tracker->idToDetUnit(hitId)->surface().toGlobal(outerHitLocal).phi());
00675             double innerHitZ = tracker->idToDetUnit(hitId)->surface().toGlobal(innerHitLocal).z(); 
00676             double outerHitZ = tracker->idToDetUnit(hitId)->surface().toGlobal(outerHitLocal).z();
00677             double innerExtrapolatedPhi = phiFromExtrapolation(d0,phi0,k0,innerHitZ,roadType);
00678             double outerExtrapolatedPhi = phiFromExtrapolation(d0,phi0,k0,outerHitZ,roadType);
00679 
00680             double midPointZ = 0.5*(innerHitZ+outerHitZ);
00681 
00682             double dPhiInter = CheckZPhiIntersection(innerHitPhi, innerHitZ, outerHitPhi, outerHitZ, 
00683                                                    innerExtrapolatedPhi, innerHitZ,
00684                                                    outerExtrapolatedPhi, outerHitZ);
00685 
00686             double dX = midPointZ*tan(dPhiInter);
00687 
00688             if (std::abs(dX) < 1.5*phiMax(roadType,phi0,k1)) {
00689               if ((usedRecHits < maxDetHitsInCloudPerDetId) && (cloud.size() < maxRecHitsInCloud_)) {
00690                 cloud.addHit(recHit);
00691                 ++usedRecHits;
00692               }
00693             }
00694 
00695           } else {
00696 
00697             //
00698             //  This is where the disk rphiRecHits end up for Roads::ZPhi
00699             //
00700             LocalPoint hit = recHit->localPosition();
00701             const StripTopology *topology = dynamic_cast<const StripTopology*>(&(tracker->idToDetUnit(hitId)->topology()));
00702             double stripAngle = topology->stripAngle(topology->strip(hit));
00703             double stripLength = topology->localStripLength(hit);
00704             // new method
00705             double hitZ = tracker->idToDetUnit(hitId)->surface().toGlobal(hit).z(); 
00706             double extrapolatedPhi = phiFromExtrapolation(d0,phi0,k0,hitZ,roadType);
00707 
00708             LocalPoint innerHitLocal(hit.x()+stripLength/2*std::sin(stripAngle),hit.y()-stripLength/2*std::cos(stripAngle),0);
00709             LocalPoint outerHitLocal(hit.x()-stripLength/2*std::sin(stripAngle),hit.y()+stripLength/2*std::cos(stripAngle),0);
00710 
00711             double innerHitPhi = map_phi(tracker->idToDetUnit(hitId)->surface().toGlobal(innerHitLocal).phi()); 
00712             double outerHitPhi = map_phi(tracker->idToDetUnit(hitId)->surface().toGlobal(outerHitLocal).phi());
00713             //double innerZ = tracker->idToDetUnit(hitId)->surface().toGlobal(innerHitLocal).z(); 
00714             //double outerZ = tracker->idToDetUnit(hitId)->surface().toGlobal(outerHitLocal).z();
00715             //if (innerZ != outerZ)  std::cout<<"HEY!!! innerZ = " << innerZ << " != outerZ = " << outerZ << std::endl;
00716 
00717             double deltaPhi = ZPhiDeltaPhi(innerHitPhi,outerHitPhi,extrapolatedPhi);
00718             double deltaX   =  hitZ*tan(deltaPhi);
00719             if (std::abs(deltaX) < phiMax(roadType,phi0,k1)){
00720               if ((usedRecHits < maxDetHitsInCloudPerDetId) && (cloud.size() < maxRecHitsInCloud_)) {
00721                 cloud.addHit(recHit);
00722                 ++usedRecHits;
00723               }
00724             }
00725           }
00726         } 
00727       }
00728     } else if ( (unsigned int)id.subdetId() == PixelSubdetector::PixelBarrel 
00729                 || (unsigned int)id.subdetId() == PixelSubdetector::PixelEndcap) {
00730       if ( UsePixels ) {
00731         
00732         const SiPixelRecHit *recHit = (SiPixelRecHit*)(*recHitIterator);
00733         
00734         if ( roadType == Roads::RPhi ) {
00735           
00736           if ( isBarrelSensor(id) ) {
00737             // Barrel Pixel, RoadType RPHI
00738             
00739             GlobalPoint ghit = tracker->idToDet(recHit->geographicalId())->surface().toGlobal(recHit->localPosition());
00740             double hitRadius = sqrt(ghit.x()*ghit.x()+ghit.y()*ghit.y());
00741             double hitphi = map_phi(ghit.phi());
00742             double phi = phiFromExtrapolation(d0,phi0,k0,hitRadius,roadType);
00743             
00744             float dp = hitphi-phi;
00745             float dx = hitRadius*tan(dp);
00746             
00747             // switch cut to dx instead of dphi
00748             if ( std::abs(dx) < phiMax(roadType,phi0,k0) ) {
00749               cloud.addHit(recHit);
00750               ++usedRecHits;
00751             }
00752           } else {
00753               
00754             // Forward Pixel,roadtype RPHI
00755               
00756             // Get Local Hit Position of the Pixel Hit
00757             LocalPoint hit = recHit->localPosition();
00758               
00759             // Get Phi of hit position 
00760             double hitphi = map_phi(tracker->idToDetUnit(id)->surface().toGlobal(hit).phi());
00761               
00762             // Get Global Hit position
00763             GlobalPoint ghit = tracker->idToDet(recHit->geographicalId())->surface().toGlobal(recHit->localPosition());
00764               
00765             // Get Hit Radis
00766             double hitRadius = sqrt(ghit.x()*ghit.x()+ghit.y()*ghit.y());
00767               
00768             // Get Phi from extrapolation
00769             double phi = phiFromExtrapolation(d0,phi0,k0,hitRadius,roadType);
00770               
00771             if ( std::abs(hitphi-phi) < phiMax(roadType,phi0,k0) ) {
00772               cloud.addHit(recHit);
00773               ++usedRecHits;
00774             }   
00775           }
00776         } else {
00777             
00778           GlobalPoint ghit = tracker->idToDet(recHit->geographicalId())->surface().toGlobal(recHit->localPosition());            
00779           double  phi = phiFromExtrapolation(d0,phi0,k0,ghit.z(),roadType);            
00780           double hitphi = map_phi(ghit.phi());
00781           double dphi = map_phi2(hitphi-phi);
00782           float dx = ghit.z()*tan(dphi);
00783             
00784           if ( std::abs(dx) < 0.25 ) {
00785             cloud.addHit(recHit);
00786             ++usedRecHits;
00787           }
00788         }
00789       }
00790     } else {
00791       edm::LogError("RoadSearch") << "recHitVector from general hit access function contains unknown detector id: " << (unsigned int)id.subdetId() << " rawId: " << id.rawId();
00792     }
00793     
00794   } //for loop over all recHits
00795   
00796   
00797   return usedRecHits;
00798 }
00799 
00800 unsigned int RoadSearchCloudMakerAlgorithm::FillPixRecHitsIntoCloud(DetId id, const SiPixelRecHitCollection *inputRecHits, 
00801                                                                     double d0, double phi0, double k0, Roads::type roadType, double ringPhi,
00802                                                                     const TrackerGeometry *tracker, RoadSearchCloud &cloud) {
00803   
00804   
00805   unsigned int usedRecHits = 0;
00806   
00807   // Get Geometry
00808   //const PixelTopology *topology = dynamic_cast<const PixelTopology*>(&(tracker->idToDetUnit(id)->topology()));
00809   
00810   
00811   // retrieve vector<SiPixelRecHit> for id
00812   // loop over SiPixelRecHit
00813   // check if compatible with cloud, fill into cloud
00814   
00815   SiPixelRecHitCollection::const_iterator recHitMatch = inputRecHits->find(id);
00816   if (recHitMatch == inputRecHits->end()) return usedRecHits;
00817 
00818   const SiPixelRecHitCollection::DetSet recHitRange = *recHitMatch;
00819   
00820   for ( SiPixelRecHitCollection::DetSet::const_iterator recHitIterator = recHitRange.begin(); 
00821         recHitIterator != recHitRange.end(); ++recHitIterator) {
00822     
00823     const SiPixelRecHit * recHit = &(*recHitIterator);
00824     
00825     if ( roadType == Roads::RPhi ) {
00826       
00827       if ( isBarrelSensor(id) ) {
00828         // Barrel Pixel, RoadType RPHI
00829         
00830         GlobalPoint ghit = tracker->idToDet(recHit->geographicalId())->surface().toGlobal(recHit->localPosition());
00831         double hitRadius = sqrt(ghit.x()*ghit.x()+ghit.y()*ghit.y());
00832         double hitphi = map_phi(ghit.phi());
00833         double phi = phiFromExtrapolation(d0,phi0,k0,hitRadius,roadType);
00834         
00835         if ( std::abs(hitphi-phi) < phiMax(roadType,phi0,k0) ) {
00836           cloud.addHit(recHit);
00837           ++usedRecHits;
00838         }
00839       } 
00840       else {
00841         
00842         // Forward Pixel,roadtype RPHI
00843         
00844         // Get Local Hit Position of the Pixel Hit
00845         LocalPoint hit = recHit->localPosition();
00846         
00847         // Get Phi of hit position 
00848         double hitphi = map_phi(tracker->idToDetUnit(id)->surface().toGlobal(hit).phi());
00849         
00850         // Get Global Hit position
00851         GlobalPoint ghit = tracker->idToDet(recHit->geographicalId())->surface().toGlobal(recHit->localPosition());
00852         
00853         // Get Hit Radis
00854         double hitRadius = sqrt(ghit.x()*ghit.x()+ghit.y()*ghit.y());
00855         
00856         // Get Phi from extrapolation
00857         double phi = phiFromExtrapolation(d0,phi0,k0,hitRadius,roadType);
00858         
00859         if ( std::abs(hitphi-phi) < phiMax(roadType,phi0,k0) ) {
00860           cloud.addHit(recHit);
00861           ++usedRecHits;
00862         }       
00863       }
00864     } 
00865     
00866     else {
00867       
00868       GlobalPoint ghit = tracker->idToDet(recHit->geographicalId())->surface().toGlobal(recHit->localPosition());
00869       
00870       double  phi = phiFromExtrapolation(d0,phi0,k0,ghit.z(),roadType);
00871       if ( (phi - phiMax(roadType,phi0,k0)) < ringPhi && (phi + phiMax(roadType,phi0,k0))>ringPhi ) {
00872         cloud.addHit(recHit);
00873         ++usedRecHits;
00874       }
00875     }
00876     
00877   }
00878   
00879   return usedRecHits;
00880 }
00881 
00882 bool RoadSearchCloudMakerAlgorithm::isSingleLayer(DetId id) {
00883   
00884   if ( (unsigned int)id.subdetId() == StripSubdetector::TIB ) {
00885     TIBDetId tibid(id.rawId()); 
00886     if ( !tibid.glued() ) {
00887       return true;
00888     }
00889   } else if ( (unsigned int)id.subdetId() == StripSubdetector::TOB ) {
00890     TOBDetId tobid(id.rawId()); 
00891     if ( !tobid.glued() ) {
00892       return true;
00893     }
00894   } else if ( (unsigned int)id.subdetId() == StripSubdetector::TID ) {
00895     TIDDetId tidid(id.rawId()); 
00896     if ( !tidid.glued() ) {
00897       return true;
00898     }
00899   } else if ( (unsigned int)id.subdetId() == StripSubdetector::TEC ) {
00900     TECDetId tecid(id.rawId()); 
00901     if ( !tecid.glued() ) {
00902       return true;
00903     }
00904   } else {
00905     return false;
00906   }
00907   
00908   return false;
00909 }
00910 
00911 bool RoadSearchCloudMakerAlgorithm::isBarrelSensor(DetId id) {
00912   
00913   if ( (unsigned int)id.subdetId() == StripSubdetector::TIB ) {
00914     return true;
00915   } else if ( (unsigned int)id.subdetId() == StripSubdetector::TOB ) {
00916     return true;
00917   } else if ( (unsigned int)id.subdetId() == PixelSubdetector::PixelBarrel ) {
00918     return true;
00919   } else {
00920     return false;
00921   }
00922   
00923 }
00924 
00925 double RoadSearchCloudMakerAlgorithm::phiFromExtrapolation(double d0, double phi0, double k0, double ringRadius, Roads::type roadType) {
00926   
00927   double ringPhi = -99.;
00928   if ( roadType == Roads::RPhi ) {
00929     double omega=k0, rl=ringRadius;
00930     double sp0=sin(phi0); double cp0=cos(phi0);  
00931     if (fabs(omega)>0.000005){
00932       double xc=-sp0*(d0+1.0/omega);    
00933       double yc=cp0*(d0+1.0/omega);
00934       double rh=fabs(1.0/omega);
00935       double bbb=fabs(d0+1.0/omega);
00936       double sss=0.5*(rl+rh+bbb);
00937       double ddd=sqrt((sss-bbb)*(sss-rh)/(sss*(sss-rl)));
00938       double phil1=2.0*atan(ddd);
00939       double phit=phi0+phil1; if (omega<0.0)phit=phi0-phil1;
00940       double xh=xc+sin(phit)/omega;
00941       double yh=yc-cos(phit)/omega;
00942       double phih=atan2(yh,xh);
00943       ringPhi = map_phi(phih);
00944     }
00945     else {
00946       double cee = rl*rl - d0*d0 -0.25*omega*omega - omega*d0;
00947       if (cee<0.0){return ringPhi;}
00948       double l = sqrt(cee);
00949       double xh=-sp0*d0+l*cp0-0.5*l*l*omega*sp0;
00950       double yh= cp0*d0+l*sp0+0.5*l*l*omega*cp0;
00951       double phih=atan2(yh,xh);
00952       ringPhi = map_phi(phih);     
00953     }
00954   } 
00955   else {
00956     ringPhi = map_phi(phi0 + k0 * ringRadius);
00957   }
00958   
00959   return ringPhi;
00960 }
00961 
00962 double RoadSearchCloudMakerAlgorithm::phiMax(Roads::type roadType,
00963                                              double phi0, double k0) {
00964   
00965   double dphi;
00966   if ( roadType == Roads::RPhi ) {
00967     // switch cut to dx instead of dphi
00968     // Still call this dphi, but will now be dx
00969     dphi = theRPhiRoadSize + 0.15*82.0*fabs(k0);
00970   }
00971   else if ( roadType == Roads::ZPhi ) {
00972     dphi = theZPhiRoadSize + 0.4*82.0*fabs(k0);
00973   }
00974   else {
00975     edm::LogWarning("RoadSearch") << "Bad roadType: "<< roadType;
00976     dphi = theZPhiRoadSize;
00977   }
00978   return dphi;
00979   
00980 }
00981 
00982 void RoadSearchCloudMakerAlgorithm::makecircle(double x1, double y1, 
00983                                                double x2,double y2, double x3, double y3){
00984   double x1t=x1-x3; double y1t=y1-y3; double r1s=x1t*x1t+y1t*y1t;
00985   double x2t=x2-x3; double y2t=y2-y3; double r2s=x2t*x2t+y2t*y2t;
00986   double rho=x1t*y2t-x2t*y1t;
00987   double xc, yc, rc, fac;
00988   if (fabs(rho)<RoadSearchCloudMakerAlgorithm::epsilon){
00989     rc=1.0/(RoadSearchCloudMakerAlgorithm::epsilon);
00990     fac=sqrt(x1t*x1t+y1t*y1t);
00991     xc=x2+y1t*rc/fac;
00992     yc=y2-x1t*rc/fac;
00993   }else{
00994     fac=0.5/rho;
00995     xc=fac*(r1s*y2t-r2s*y1t);
00996     yc=fac*(r2s*x1t-r1s*x2t); 
00997     rc=sqrt(xc*xc+yc*yc); xc+=x3; yc+=y3;
00998   }
00999   double s3=0.0;
01000   double f1=x1*yc-y1*xc; double f2=x2*yc-y2*xc; 
01001   double f3=x3*yc-y3*xc;
01002   if ((f1<0.0)&&(f2<0.0)&&(f3<=0.0))s3=1.0;
01003   if ((f1>0.0)&&(f2>0.0)&&(f3>=0.0))s3=-1.0;
01004   d0h=-s3*(sqrt(xc*xc+yc*yc)-rc);
01005   phi0h=atan2(yc,xc)+s3*Geom::halfPi();
01006   omegah=-s3/rc;
01007 }
01008 
01009 double RoadSearchCloudMakerAlgorithm::CheckXYIntersection(LocalPoint& inner1, LocalPoint& outer1,
01010                                                         LocalPoint& inner2, LocalPoint& outer2){
01011 
01012   double deltaX = -999.;
01013   // just get the x coord of intersection of two line segments
01014   // check if intersection lies inside segments
01015   double det12 = inner1.x()*outer1.y() - inner1.y()*outer1.x();
01016   double det34 = inner2.x()*outer2.y() - inner2.y()*outer2.x();
01017 
01018   double xinter = (det12*(inner2.x()-outer2.x()) - det34*(inner1.x()-outer1.x()))/
01019     ((inner1.x()-outer1.x())*(inner2.y()-outer2.y()) - 
01020      (inner2.x()-outer2.x())*(inner1.y()-outer1.y()));
01021 
01022   bool inter = true;
01023   if (inner1.x() < outer1.x()){
01024     if ((xinter<inner1.x()) || (xinter>outer1.x())) inter = false;
01025   }
01026   else{
01027     if ((xinter>inner1.x()) || (xinter<outer1.x())) inter = false;
01028   }
01029 
01030   if (inner2.x() < outer2.x()){
01031     if ((xinter<inner2.x()) || (xinter>outer2.x())) inter = false;
01032   }
01033   else{
01034     if ((xinter>inner2.x()) || (xinter<outer2.x())) inter = false;
01035   }
01036 
01037   if (inter){
01038     deltaX = 0;
01039   }
01040   else{
01041     deltaX = min(fabs(inner1.x()-inner2.x()),fabs(outer1.x()-outer2.x()));
01042   }
01043   return deltaX;
01044 
01045 }
01046 
01047 double RoadSearchCloudMakerAlgorithm::CheckZPhiIntersection(double iPhi1, double iZ1, double oPhi1, double oZ1,
01048                                                             double iPhi2, double iZ2, double oPhi2, double oZ2){
01049   
01050   // Have to make sure all are in the same hemisphere
01051   if ((iPhi1 > Geom::pi() || oPhi1 > Geom::pi() || iPhi2 > Geom::pi() || oPhi2 > Geom::pi()) &&
01052       (iPhi1 < Geom::pi() || oPhi1 < Geom::pi() || iPhi2 < Geom::pi() || oPhi2 < Geom::pi())){
01053     iPhi1 = map_phi2(iPhi1);  oPhi1 = map_phi2(oPhi1);
01054     iPhi2 = map_phi2(iPhi2);  oPhi2 = map_phi2(oPhi2);
01055   }
01056 
01057   double deltaPhi = -999.;
01058   // just get the x coord of intersection of two line segments
01059   // check if intersection lies inside segments
01060   double det12 = iZ1*oPhi1 - iPhi1*oZ1;
01061   double det34 = iZ2*oPhi2 - iPhi2*oZ2;
01062 
01063   double xinter = (det12*(iZ2-oZ2) - det34*(iZ1-oZ1))/
01064     ((iZ1-oZ1)*(iPhi2-oPhi2) - 
01065      (iZ2-oZ2)*(iPhi1-oPhi1));
01066 
01067   bool inter = true;
01068   if (iZ1 < oZ1){
01069     if ((xinter<iZ1) || (xinter>oZ1)) inter = false;
01070   }
01071   else{
01072     if ((xinter>iZ1) || (xinter<oZ1)) inter = false;
01073   }
01074 
01075   if (iZ2 < oZ2){
01076     if ((xinter<iZ2) || (xinter>oZ2)) inter = false;
01077   }
01078   else{
01079     if ((xinter>iZ2) || (xinter<oZ2)) inter = false;
01080   }
01081 
01082   if (inter){
01083     deltaPhi = 0;
01084   }
01085   else{
01086     deltaPhi = min(fabs(iPhi2-iPhi1),fabs(oPhi2-oPhi1));
01087   }
01088   return deltaPhi;
01089 
01090 }
01091 
01092 
01093 double RoadSearchCloudMakerAlgorithm::ZPhiDeltaPhi(double hitPhi1, double hitPhi2, double predictedPhi){
01094 
01095   double deltaPhi = -999.;
01096 
01097   double dPhiHits = map_phi2(hitPhi1-hitPhi2);
01098   double dPhi1 = map_phi2(hitPhi1-predictedPhi);
01099   double dPhi2 = map_phi2(hitPhi2-predictedPhi);
01100 
01101   if (dPhiHits >= 0){  // hitPhi1 >= hitPhi2
01102     if ( (dPhi1>=0.0) && (dPhi2 <= 0.0))
01103       deltaPhi = 0.0;
01104     else{
01105       if (std::abs(dPhi1)<std::abs(dPhi2))
01106         deltaPhi = dPhi1;
01107       else
01108         deltaPhi = dPhi2;
01109     }
01110   }
01111   else { // hitPhi1 < hitPhi2
01112     if ( (dPhi1<=0.0) && (dPhi2 >= 0.0))
01113       deltaPhi = 0.0;
01114     else{
01115       if (std::abs(dPhi1)<std::abs(dPhi2))
01116         deltaPhi = dPhi1;
01117       else
01118         deltaPhi = dPhi2;
01119     }
01120   }
01121 
01122   return deltaPhi;
01123 
01124 }
01125 
01126 RoadSearchCloudCollection RoadSearchCloudMakerAlgorithm::Clean(RoadSearchCloudCollection* inputCollection){
01127   
01128   RoadSearchCloudCollection output;
01129   
01130   //
01131   //  no raw clouds - nothing to try merging
01132   //
01133 
01134   if ( inputCollection->empty() ){
01135     LogDebug("RoadSearch") << "Found " << output.size() << " clean clouds.";
01136     return output;  
01137   }
01138   
01139   //
01140   //  1 raw cloud - nothing to try merging, but one cloud to duplicate
01141   //
01142   
01143   if ( 1==inputCollection->size() ){
01144     output = *inputCollection;
01145 //     RoadSearchCloud *temp = inputCollection->begin()->clone();
01146 //     output.push_back(*temp);
01147 //     delete temp;
01148     LogDebug("RoadSearch") << "Found " << output.size() << " clean clouds.";
01149     return output;
01150   }  
01151   
01152   //
01153   //  got > 1 raw cloud - something to try merging
01154   //
01155   std::vector<bool> already_gone(inputCollection->size());
01156   for (unsigned int i=0; i<inputCollection->size(); ++i) {
01157     already_gone[i] = false; 
01158   }
01159   
01160   int raw_cloud_ctr=0;
01161   // loop over clouds
01162   for ( RoadSearchCloudCollection::const_iterator raw_cloud = inputCollection->begin(); raw_cloud != inputCollection->end(); ++raw_cloud) {
01163     ++raw_cloud_ctr;
01164     
01165     if (already_gone[raw_cloud_ctr-1])continue;
01166     
01167     // produce output cloud where other clouds are merged in
01168     // create temp pointer for clone which will be deleted afterwards
01169 //     RoadSearchCloud *temp_lone_cloud = raw_cloud->clone();
01170 //     RoadSearchCloud lone_cloud = *temp_lone_cloud;
01171     RoadSearchCloud lone_cloud = *raw_cloud;
01172 
01173     int second_cloud_ctr=raw_cloud_ctr;
01174     for ( RoadSearchCloudCollection::const_iterator second_cloud = raw_cloud+1; second_cloud != inputCollection->end(); ++second_cloud) {
01175       second_cloud_ctr++;
01176       
01177       std::vector<const TrackingRecHit*> unshared_hits;
01178       
01179       if ( already_gone[second_cloud_ctr-1] )continue;
01180       
01181       for ( RoadSearchCloud::RecHitVector::const_iterator second_cloud_hit = second_cloud->begin_hits();
01182             second_cloud_hit != second_cloud->end_hits();
01183             ++ second_cloud_hit ) {
01184         bool is_shared = false;
01185         for ( RoadSearchCloud::RecHitVector::const_iterator lone_cloud_hit = lone_cloud.begin_hits();
01186               lone_cloud_hit != lone_cloud.end_hits();
01187               ++ lone_cloud_hit ) {
01188           
01189           if ((*lone_cloud_hit)->geographicalId() == (*second_cloud_hit)->geographicalId())
01190             if ((*lone_cloud_hit)->localPosition().x() == (*second_cloud_hit)->localPosition().x())
01191               if ((*lone_cloud_hit)->localPosition().y() == (*second_cloud_hit)->localPosition().y())
01192                 {is_shared=true; break;}
01193         }
01194         if (!is_shared)  unshared_hits.push_back(*second_cloud_hit);
01195           
01196         if ( ((float(unshared_hits.size())/float(lone_cloud.size())) > 
01197               ((float(second_cloud->size())/float(lone_cloud.size()))-mergingFraction_)) &&
01198              ((float(unshared_hits.size())/float(second_cloud->size())) > (1-mergingFraction_))){
01199           // You'll never merge these clouds..... Could quit now!
01200           break;
01201         }
01202           
01203         if (lone_cloud.size()+unshared_hits.size() > maxRecHitsInCloud_) {
01204           break;
01205         }
01206           
01207       }
01208         
01209       double f_lone_shared=double(second_cloud->size()-unshared_hits.size())/double(lone_cloud.size());
01210       double f_second_shared=double(second_cloud->size()-unshared_hits.size())/double(second_cloud->size());
01211       
01212       if ( ( (static_cast<unsigned int>(f_lone_shared*1E9) > static_cast<unsigned int>(mergingFraction_*1E9))||(static_cast<unsigned int>(f_second_shared*1E9) > static_cast<unsigned int>(mergingFraction_*1E9)) ) 
01213            && (lone_cloud.size()+unshared_hits.size() <= maxRecHitsInCloud_) ){
01214         
01215         LogDebug("RoadSearch") << " Merge CloudA: " << raw_cloud_ctr << " with  CloudB: " << second_cloud_ctr 
01216                                << " Shared fractions are " << f_lone_shared << " and " << f_second_shared;
01217           
01218         //
01219         //  got a cloud to merge
01220         //
01221         for (unsigned int k=0; k<unshared_hits.size(); ++k) {
01222           lone_cloud.addHit(unshared_hits[k]);
01223         }
01224           
01225         already_gone[second_cloud_ctr-1]=true;
01226         
01227       }//end got a cloud to merge
01228       
01229     }//interate over all second clouds
01230       
01231     output.push_back(lone_cloud);
01232     
01233   }//iterate over all raw clouds
01234   
01235   LogDebug("RoadSearch") << "Found " << output.size() << " clean clouds.";
01236   
01237   return output;
01238 }
01239 
01240 SiStripMatchedRecHit2D* RoadSearchCloudMakerAlgorithm::CorrectMatchedHit(const TrackingRecHit *originalHit,
01241                                                                          const GluedGeomDet* gluedDet,
01242                                                                          const TrackerGeometry *tracker,
01243                                                                          const SiStripRecHitMatcher* theHitMatcher,
01244                                                                          double k0, double phi0) {
01245  
01246           const SiStripMatchedRecHit2D *theRH = dynamic_cast<const SiStripMatchedRecHit2D*>(originalHit);
01247           if (theRH == 0) {
01248             std::cout<<" Could not cast original hit" << std::endl;
01249           }
01250           if (theRH != 0){
01251             const GeomDet *recHitGeomDet = tracker->idToDet(theRH->geographicalId());
01252             const GluedGeomDet *theGluedDet = dynamic_cast<const GluedGeomDet*>(recHitGeomDet);
01253              
01254             const GeomDetUnit* theMonoDet = theGluedDet->monoDet();
01255             const SiStripRecHit2D* theMonoHit   = theRH->monoHit();
01256             GlobalPoint monoRHpos = (theMonoDet->surface()).toGlobal(theMonoHit->localPosition());
01257              
01258             GlobalPoint gcenterofstrip=(theMonoDet->surface()).toGlobal(theMonoHit->localPosition());
01259              
01260             float gtrackangle_xy = map_phi2(phi0 + 2.0*asin(0.5*gcenterofstrip.perp()*k0));
01261             float rzangle = atan2(gcenterofstrip.perp(),gcenterofstrip.z());
01262  
01263             GlobalVector gtrackangle2(cos(gtrackangle_xy)*sin(rzangle),
01264                                       sin(gtrackangle_xy)*sin(rzangle),
01265                                       cos(rzangle));
01266             LocalVector trackdirection2=((tracker->idToDet(theRH->geographicalId()))->surface()).toLocal(gtrackangle2);
01267             GlobalVector gdir = theMonoDet->surface().toGlobal(trackdirection2);
01268  
01269             SiStripMatchedRecHit2D* theCorrectedHit = theHitMatcher->match(theRH,theGluedDet,trackdirection2);
01270             if (theCorrectedHit!=0) return theCorrectedHit;
01271           }
01272  
01273           return 0;
01274 }