00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
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
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
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
00161 const int nphibin = 24;
00162 const int netabin = 24;
00163 RoadSearchCloudCollection CloudArray[nphibin][netabin];
00164
00165
00166 edm::ESHandle<Roads> roads;
00167 es.get<RoadMapRecord>().get(roadsLabel_, roads);
00168
00169
00170 const RoadSearchSeedCollection* inputSeeds = input.product();
00171
00172
00173 recHitVectorClass.setCollections(rphiRecHits,stereoRecHits,matchedRecHits,pixRecHits);
00174 recHitVectorClass.setMode(DetHitAccess::standard);
00175
00176
00177 edm::ESHandle<TrackerGeometry> tracker;
00178 es.get<TrackerDigiGeometryRecord>().get(tracker);
00179
00180
00181 SiStripRecHitMatcher* theHitMatcher = new SiStripRecHitMatcher(3.0);
00182
00183 edm::LogInfo("RoadSearch") << "Found " << inputSeeds->size() << " input seeds.";
00184
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
00197
00198
00199
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
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
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
00239 if ( roadType == Roads::RPhi ) {
00240 double dr = outerSeedHitGlobalPosition.perp() - innerSeedHitGlobalPosition.perp();
00241 const double dr_min = 1;
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
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;
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
00276 k0 = map_phi2(outerSeedHitGlobalPosition.phi() - innerSeedHitGlobalPosition.phi()) / dz;
00277 phi0 = map_phi(innerSeedHitGlobalPosition.phi() - k0 * innerSeedHitGlobalPosition.z());
00278
00279
00280 makecircle(innerSeedHitGlobalPosition.x(),innerSeedHitGlobalPosition.y(),
00281 outerSeedHitGlobalPosition.x(),outerSeedHitGlobalPosition.y(),
00282 0.0,0.0);
00283 phi1 = phi0h;
00284 k1 = omegah;
00285 }
00286 }
00287
00288
00289 if ( (phi0 != -99.) && (k0 != -99999999.99) ) {
00290 const Roads::RoadSet *roadSet = seed->getSet();
00291
00292
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
00304
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
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
00327
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
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
00376 totalLayers = roadSet->size() - layerCounter + 1;
00377
00378
00379
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
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
00412 if ( consecutiveMissedLayers > maxNumberOfConsecutiveMissedLayersPerCloud ) {
00413 LogDebug("RoadSearch") << "BREAK: seed # "<<seed-inputSeeds->begin()
00414 << " More than " << maxNumberOfConsecutiveMissedLayersPerCloud << " missed consecutive layers!";
00415 break;
00416 }
00417
00418
00419 if ( missedLayers > maxNumberOfMissedLayersPerCloud ) {
00420 LogDebug("RoadSearch") << "BREAK: seed # "<<seed-inputSeeds->begin()
00421 << " More than " << maxNumberOfMissedLayersPerCloud << " missed layers!";
00422 break;
00423 }
00424
00425
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
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
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
00524
00525
00526
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
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 {
00556 if ( isBarrelSensor(hitId) ) {
00557
00558
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
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
00614
00615 if ((usedRecHits < maxDetHitsInCloudPerDetId) && (cloud.size() < maxRecHitsInCloud_)) {
00616 cloud.addHit(recHit);
00617 ++usedRecHits;
00618 }
00619 }
00620
00621
00622 }
00623 }
00624 } else {
00625
00626
00627
00628 if (double_ring_layer && isSingleLayer(hitId)) {
00629
00630
00631
00632
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
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 {
00660 if ( isBarrelSensor(hitId) ) {
00661
00662
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
00670
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
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
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
00714
00715
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
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
00748 if ( std::abs(dx) < phiMax(roadType,phi0,k0) ) {
00749 cloud.addHit(recHit);
00750 ++usedRecHits;
00751 }
00752 } else {
00753
00754
00755
00756
00757 LocalPoint hit = recHit->localPosition();
00758
00759
00760 double hitphi = map_phi(tracker->idToDetUnit(id)->surface().toGlobal(hit).phi());
00761
00762
00763 GlobalPoint ghit = tracker->idToDet(recHit->geographicalId())->surface().toGlobal(recHit->localPosition());
00764
00765
00766 double hitRadius = sqrt(ghit.x()*ghit.x()+ghit.y()*ghit.y());
00767
00768
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 }
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
00808
00809
00810
00811
00812
00813
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
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
00843
00844
00845 LocalPoint hit = recHit->localPosition();
00846
00847
00848 double hitphi = map_phi(tracker->idToDetUnit(id)->surface().toGlobal(hit).phi());
00849
00850
00851 GlobalPoint ghit = tracker->idToDet(recHit->geographicalId())->surface().toGlobal(recHit->localPosition());
00852
00853
00854 double hitRadius = sqrt(ghit.x()*ghit.x()+ghit.y()*ghit.y());
00855
00856
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
00968
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
01014
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
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
01059
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){
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 {
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
01132
01133
01134 if ( inputCollection->empty() ){
01135 LogDebug("RoadSearch") << "Found " << output.size() << " clean clouds.";
01136 return output;
01137 }
01138
01139
01140
01141
01142
01143 if ( 1==inputCollection->size() ){
01144 output = *inputCollection;
01145
01146
01147
01148 LogDebug("RoadSearch") << "Found " << output.size() << " clean clouds.";
01149 return output;
01150 }
01151
01152
01153
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
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
01168
01169
01170
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
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
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 }
01228
01229 }
01230
01231 output.push_back(lone_cloud);
01232
01233 }
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 }