#include <RoadSearchCloudMakerAlgorithm.h>

Public Member Functions
double	CheckXYIntersection (LocalPoint &ip1, LocalPoint &op1, LocalPoint &ip2, LocalPoint &op2)

double	CheckZPhiIntersection (double iPhi1, double iZ1, double oPhi1, double oZ1, double iPhi2, double iZ2, double oPhi2, double oZ2)

RoadSearchCloudCollection	Clean (RoadSearchCloudCollection *rawColl)

SiStripMatchedRecHit2D *	CorrectMatchedHit (const TrackingRecHit originalRH, const GluedGeomDet gluedDet, const TrackerGeometry tracker, const SiStripRecHitMatcher theHitMatcher, double k0, double phi0)

unsigned int	FillPixRecHitsIntoCloud (DetId id, const SiPixelRecHitCollection inputRecHits, double d0, double phi0, double k0, Roads::type roadType, double ringPhi, const TrackerGeometry tracker, RoadSearchCloud &cloud)

unsigned int	FillRecHitsIntoCloud (DetId id, const SiStripRecHit2DCollection inputRecHits, double d0, double phi0, double k0, Roads::type roadType, double ringPhi, const TrackerGeometry tracker, RoadSearchCloud &cloud)

unsigned int	FillRecHitsIntoCloudGeneral (DetId id, double d0, double phi0, double k0, double phi1, double k1, Roads::type roadType, double ringPhi, const TrackerGeometry tracker, const SiStripRecHitMatcher theHitMatcher, RoadSearchCloud &cloud)

bool	isBarrelSensor (DetId id)

bool	isSingleLayer (DetId id)

void	makecircle (double x1_cs, double y1_cs, double x2_cs, double y2_cs, double x3_cs, double y3_cs)

double	map_phi (double phi)

double	map_phi2 (double phi)

double	phiFromExtrapolation (double d0, double phi0, double k0, double ringRadius, Roads::type roadType)

double	phiMax (Roads::type roadType, double phi0, double k0)

	RoadSearchCloudMakerAlgorithm (const edm::ParameterSet &conf)

void	run (edm::Handle< RoadSearchSeedCollection > input, const SiStripRecHit2DCollection rphiRecHits, const SiStripRecHit2DCollection stereoRecHits, const SiStripMatchedRecHit2DCollection matchedRecHits, const SiPixelRecHitCollection pixRecHits, const edm::EventSetup &es, RoadSearchCloudCollection &output)
	Runs the algorithm. More...

double	ZPhiDeltaPhi (double phi1, double phi2, double phiExpect)

	~RoadSearchCloudMakerAlgorithm ()

Private Attributes
edm::ParameterSet	conf_

double	d0h

bool	doCleaning_

unsigned int	increaseMaxNumberOfConsecutiveMissedLayersPerCloud

unsigned int	increaseMaxNumberOfMissedLayersPerCloud

unsigned int	maxDetHitsInCloudPerDetId

double	maxFractionOfConsecutiveMissedLayersPerCloud

double	maxFractionOfMissedLayersPerCloud

unsigned int	maxRecHitsInCloud_

double	mergingFraction_

double	minFractionOfUsedLayersPerCloud

bool	NoFieldCosmic

double	omegah

std::ostringstream	output_

double	phi0h

DetHitAccess	recHitVectorClass

std::string	roadsLabel_

double	rphicsq

int	rphinhits

double	scalefactorRoadSeedWindow_

double	theMinimumHalfRoad

const SiPixelRecHitCollection	thePixRecHits

double	theRPhiRoadSize

double	theZPhiRoadSize

bool	UsePixels

Static Private Attributes
static double	epsilon = 0.000000001

Detailed Description

Definition at line 82 of file RoadSearchCloudMakerAlgorithm.h.

Constructor & Destructor Documentation

RoadSearchCloudMakerAlgorithm::RoadSearchCloudMakerAlgorithm ( const edm::ParameterSet & conf )

Definition at line 103 of file RoadSearchCloudMakerAlgorithm.cc.

References conf_, doCleaning_, edm::ParameterSet::getParameter(), increaseMaxNumberOfConsecutiveMissedLayersPerCloud, increaseMaxNumberOfMissedLayersPerCloud, maxDetHitsInCloudPerDetId, maxFractionOfConsecutiveMissedLayersPerCloud, maxFractionOfMissedLayersPerCloud, maxRecHitsInCloud_, mergingFraction_, minFractionOfUsedLayersPerCloud, NoFieldCosmic, recHitVectorClass, roadsLabel_, scalefactorRoadSeedWindow_, DetHitAccess::setMode(), DetHitAccess::standard, theMinimumHalfRoad, theRPhiRoadSize, theZPhiRoadSize, DetHitAccess::use_rphiRecHits(), DetHitAccess::use_stereoRecHits(), and UsePixels.

                                                                                         : conf_(conf) { 
   recHitVectorClass.setMode(DetHitAccess::standard);    
   recHitVectorClass.use_rphiRecHits(conf_.getParameter<bool>("UseRphiRecHits"));
   recHitVectorClass.use_stereoRecHits(conf_.getParameter<bool>("UseStereoRecHits"));
   
   
   theRPhiRoadSize =  conf_.getParameter<double>("RPhiRoadSize");
   theZPhiRoadSize =  conf_.getParameter<double>("ZPhiRoadSize");
   UsePixels = conf_.getParameter<bool>("UsePixelsinRS");
   NoFieldCosmic = conf_.getParameter<bool>("StraightLineNoBeamSpotCloud");
   theMinimumHalfRoad = conf_.getParameter<double>("MinimumHalfRoad");
   
   maxDetHitsInCloudPerDetId = conf_.getParameter<unsigned int>("MaxDetHitsInCloudPerDetId");
   minFractionOfUsedLayersPerCloud = conf_.getParameter<double>("MinimalFractionOfUsedLayersPerCloud");
   maxFractionOfMissedLayersPerCloud = conf_.getParameter<double>("MaximalFractionOfMissedLayersPerCloud");
   maxFractionOfConsecutiveMissedLayersPerCloud = conf_.getParameter<double>("MaximalFractionOfConsecutiveMissedLayersPerCloud");
   increaseMaxNumberOfConsecutiveMissedLayersPerCloud = conf_.getParameter<unsigned int>("IncreaseMaxNumberOfConsecutiveMissedLayersPerCloud");
   increaseMaxNumberOfMissedLayersPerCloud = conf_.getParameter<unsigned int>("IncreaseMaxNumberOfMissedLayersPerCloud");
   
   doCleaning_ = conf.getParameter<bool>("DoCloudCleaning");
   mergingFraction_ = conf.getParameter<double>("MergingFraction");
   maxRecHitsInCloud_ = (unsigned int)conf.getParameter<int>("MaxRecHitsInCloud");
   scalefactorRoadSeedWindow_ = conf.getParameter<double>("scalefactorRoadSeedWindow");
   
   roadsLabel_ = conf.getParameter<std::string>("RoadsLabel");
 
 }

RoadSearchCloudMakerAlgorithm::~RoadSearchCloudMakerAlgorithm ( )

Definition at line 131 of file RoadSearchCloudMakerAlgorithm.cc.

                                                               {
 
 }

Member Function Documentation

double RoadSearchCloudMakerAlgorithm::CheckXYIntersection	(	LocalPoint &	ip1,
		LocalPoint &	op1,
		LocalPoint &	ip2,
		LocalPoint &	op2
	)

Definition at line 1007 of file RoadSearchCloudMakerAlgorithm.cc.

References inter, min, PV3DBase< T, PVType, FrameType >::x(), and PV3DBase< T, PVType, FrameType >::y().

Referenced by FillRecHitsIntoCloudGeneral().

                                                                    {
 
   double deltaX = -999.;
   // just get the x coord of intersection of two line segments
   // check if intersection lies inside segments
   double det12 = inner1.x()*outer1.y() - inner1.y()*outer1.x();
   double det34 = inner2.x()*outer2.y() - inner2.y()*outer2.x();
 
   double xinter = (det12*(inner2.x()-outer2.x()) - det34*(inner1.x()-outer1.x()))/
     ((inner1.x()-outer1.x())*(inner2.y()-outer2.y()) - 
      (inner2.x()-outer2.x())*(inner1.y()-outer1.y()));
 
   bool inter = true;
   if (inner1.x() < outer1.x()){
     if ((xinter<inner1.x()) || (xinter>outer1.x())) inter = false;
   }
   else{
     if ((xinter>inner1.x()) || (xinter<outer1.x())) inter = false;
   }
 
   if (inner2.x() < outer2.x()){
     if ((xinter<inner2.x()) || (xinter>outer2.x())) inter = false;
   }
   else{
     if ((xinter>inner2.x()) || (xinter<outer2.x())) inter = false;
   }
 
   if (inter){
     deltaX = 0;
   }
   else{
     deltaX = min(fabs(inner1.x()-inner2.x()),fabs(outer1.x()-outer2.x()));
   }
   return deltaX;
 
 }

double RoadSearchCloudMakerAlgorithm::CheckZPhiIntersection	(	double	iPhi1,
		double	iZ1,
		double	oPhi1,
		double	oZ1,
		double	iPhi2,
		double	iZ2,
		double	oPhi2,
		double	oZ2
	)

Definition at line 1045 of file RoadSearchCloudMakerAlgorithm.cc.

References SiPixelRawToDigiRegional_cfi::deltaPhi, inter, map_phi2(), min, and Geom::pi().

Referenced by FillRecHitsIntoCloudGeneral().

                                                                                    {
   
   // Have to make sure all are in the same hemisphere
   if ((iPhi1 > Geom::pi() || oPhi1 > Geom::pi() || iPhi2 > Geom::pi() || oPhi2 > Geom::pi()) &&
       (iPhi1 < Geom::pi() || oPhi1 < Geom::pi() || iPhi2 < Geom::pi() || oPhi2 < Geom::pi())){
     iPhi1 = map_phi2(iPhi1);  oPhi1 = map_phi2(oPhi1);
     iPhi2 = map_phi2(iPhi2);  oPhi2 = map_phi2(oPhi2);
   }
 
   double deltaPhi = -999.;
   // just get the x coord of intersection of two line segments
   // check if intersection lies inside segments
   double det12 = iZ1*oPhi1 - iPhi1*oZ1;
   double det34 = iZ2*oPhi2 - iPhi2*oZ2;
 
   double xinter = (det12*(iZ2-oZ2) - det34*(iZ1-oZ1))/
     ((iZ1-oZ1)*(iPhi2-oPhi2) - 
      (iZ2-oZ2)*(iPhi1-oPhi1));
 
   bool inter = true;
   if (iZ1 < oZ1){
     if ((xinter<iZ1) || (xinter>oZ1)) inter = false;
   }
   else{
     if ((xinter>iZ1) || (xinter<oZ1)) inter = false;
   }
 
   if (iZ2 < oZ2){
     if ((xinter<iZ2) || (xinter>oZ2)) inter = false;
   }
   else{
     if ((xinter>iZ2) || (xinter<oZ2)) inter = false;
   }
 
   if (inter){
     deltaPhi = 0;
   }
   else{
     deltaPhi = min(fabs(iPhi2-iPhi1),fabs(oPhi2-oPhi1));
   }
   return deltaPhi;
 
 }

RoadSearchCloudCollection RoadSearchCloudMakerAlgorithm::Clean ( RoadSearchCloudCollection * rawColl )

Definition at line 1124 of file RoadSearchCloudMakerAlgorithm.cc.

References RoadSearchCloud::addHit(), RoadSearchCloud::begin_hits(), RoadSearchCloud::end_hits(), i, gen::k, LogDebug, maxRecHitsInCloud_, mergingFraction_, convertSQLitetoXML_cfg::output, and RoadSearchCloud::size().

Referenced by run().

                                                                                                         {
   
   RoadSearchCloudCollection output;
   
   //
   //  no raw clouds - nothing to try merging
   //
 
   if ( inputCollection->empty() ){
     LogDebug("RoadSearch") << "Found " << output.size() << " clean clouds.";
     return output;  
   }
   
   //
   //  1 raw cloud - nothing to try merging, but one cloud to duplicate
   //
   
   if ( 1==inputCollection->size() ){
     output = *inputCollection;
 //     RoadSearchCloud *temp = inputCollection->begin()->clone();
 //     output.push_back(*temp);
 //     delete temp;
     LogDebug("RoadSearch") << "Found " << output.size() << " clean clouds.";
     return output;
   }  
   
   //
   //  got > 1 raw cloud - something to try merging
   //
   std::vector<bool> already_gone(inputCollection->size());
   for (unsigned int i=0; i<inputCollection->size(); ++i) {
     already_gone[i] = false; 
   }
   
   int raw_cloud_ctr=0;
   // loop over clouds
   for ( RoadSearchCloudCollection::const_iterator raw_cloud = inputCollection->begin(); raw_cloud != inputCollection->end(); ++raw_cloud) {
     ++raw_cloud_ctr;
     
     if (already_gone[raw_cloud_ctr-1])continue;
     
     // produce output cloud where other clouds are merged in
     // create temp pointer for clone which will be deleted afterwards
 //     RoadSearchCloud *temp_lone_cloud = raw_cloud->clone();
 //     RoadSearchCloud lone_cloud = *temp_lone_cloud;
     RoadSearchCloud lone_cloud = *raw_cloud;
 
     int second_cloud_ctr=raw_cloud_ctr;
     for ( RoadSearchCloudCollection::const_iterator second_cloud = raw_cloud+1; second_cloud != inputCollection->end(); ++second_cloud) {
       second_cloud_ctr++;
       
       std::vector<const TrackingRecHit*> unshared_hits;
       
       if ( already_gone[second_cloud_ctr-1] )continue;
       
       for ( RoadSearchCloud::RecHitVector::const_iterator second_cloud_hit = second_cloud->begin_hits();
             second_cloud_hit != second_cloud->end_hits();
             ++ second_cloud_hit ) {
         bool is_shared = false;
         for ( RoadSearchCloud::RecHitVector::const_iterator lone_cloud_hit = lone_cloud.begin_hits();
               lone_cloud_hit != lone_cloud.end_hits();
               ++ lone_cloud_hit ) {
           
           if ((*lone_cloud_hit)->geographicalId() == (*second_cloud_hit)->geographicalId())
             if ((*lone_cloud_hit)->localPosition().x() == (*second_cloud_hit)->localPosition().x())
               if ((*lone_cloud_hit)->localPosition().y() == (*second_cloud_hit)->localPosition().y())
         {is_shared=true; break;}
         }
     if (!is_shared)  unshared_hits.push_back(*second_cloud_hit);
           
     if ( ((float(unshared_hits.size())/float(lone_cloud.size())) > 
           ((float(second_cloud->size())/float(lone_cloud.size()))-mergingFraction_)) &&
          ((float(unshared_hits.size())/float(second_cloud->size())) > (1-mergingFraction_))){
       // You'll never merge these clouds..... Could quit now!
       break;
     }
           
     if (lone_cloud.size()+unshared_hits.size() > maxRecHitsInCloud_) {
       break;
     }
           
       }
         
       double f_lone_shared=double(second_cloud->size()-unshared_hits.size())/double(lone_cloud.size());
       double f_second_shared=double(second_cloud->size()-unshared_hits.size())/double(second_cloud->size());
       
       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)) ) 
        && (lone_cloud.size()+unshared_hits.size() <= maxRecHitsInCloud_) ){
     
     LogDebug("RoadSearch") << " Merge CloudA: " << raw_cloud_ctr << " with  CloudB: " << second_cloud_ctr 
                    << " Shared fractions are " << f_lone_shared << " and " << f_second_shared;
           
     //
     //  got a cloud to merge
     //
     for (unsigned int k=0; k<unshared_hits.size(); ++k) {
       lone_cloud.addHit(unshared_hits[k]);
     }
           
     already_gone[second_cloud_ctr-1]=true;
         
       }//end got a cloud to merge
       
     }//interate over all second clouds
       
     output.push_back(lone_cloud);
     
   }//iterate over all raw clouds
   
   LogDebug("RoadSearch") << "Found " << output.size() << " clean clouds.";
   
   return output;
 }

SiStripMatchedRecHit2D * RoadSearchCloudMakerAlgorithm::CorrectMatchedHit	(	const TrackingRecHit *	originalRH,
		const GluedGeomDet *	gluedDet,
		const TrackerGeometry *	tracker,
		const SiStripRecHitMatcher *	theHitMatcher,
		double	k0,
		double	phi0
	)

Definition at line 1238 of file RoadSearchCloudMakerAlgorithm.cc.

References funct::cos(), gather_cfg::cout, TrackingRecHit::geographicalId(), TrackerGeometry::idToDet(), reco::ParticleMasses::k0, BaseTrackerRecHit::localPosition(), map_phi2(), SiStripRecHitMatcher::match(), GluedGeomDet::monoDet(), SiStripMatchedRecHit2D::monoHit(), PV3DBase< T, PVType, FrameType >::perp(), funct::sin(), GeomDet::surface(), toLocal(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by FillRecHitsIntoCloudGeneral().

                                                                                                  {
           // VI January 2012
           // this is not supported anymore w/o cpe
     
           const SiStripMatchedRecHit2D *theRH = dynamic_cast<const SiStripMatchedRecHit2D*>(originalHit);
           if (theRH == 0) {
             std::cout<<" Could not cast original hit" << std::endl;
           }
           if (theRH != 0){
             const GeomDet *recHitGeomDet = tracker->idToDet(theRH->geographicalId());
             const GluedGeomDet *theGluedDet = dynamic_cast<const GluedGeomDet*>(recHitGeomDet);
              
             const GeomDetUnit* theMonoDet = theGluedDet->monoDet();
             const SiStripRecHit2D theMonoHit   = theRH->monoHit();
             //GlobalPoint monoRHpos = (theMonoDet->surface()).toGlobal(theMonoHit.localPosition());
              
             GlobalPoint gcenterofstrip=(theMonoDet->surface()).toGlobal(theMonoHit.localPosition());
              
             float gtrackangle_xy = map_phi2(phi0 + 2.0*asin(0.5*gcenterofstrip.perp()*k0));
             float rzangle = atan2(gcenterofstrip.perp(),gcenterofstrip.z());
  
             GlobalVector gtrackangle2(cos(gtrackangle_xy)*sin(rzangle),
                                       sin(gtrackangle_xy)*sin(rzangle),
                                       cos(rzangle));
             LocalVector trackdirection2=((tracker->idToDet(theRH->geographicalId()))->surface()).toLocal(gtrackangle2);
             //GlobalVector gdir = theMonoDet->surface().toGlobal(trackdirection2);
  
             SiStripMatchedRecHit2D* theCorrectedHit = theHitMatcher->match(theRH,theGluedDet,trackdirection2);
             if (theCorrectedHit!=0) return theCorrectedHit;
           }
  
           return 0;
 }

unsigned int RoadSearchCloudMakerAlgorithm::FillPixRecHitsIntoCloud	(	DetId	id,
		const SiPixelRecHitCollection *	inputRecHits,
		double	d0,
		double	phi0,
		double	k0,
		Roads::type	roadType,
		double	ringPhi,
		const TrackerGeometry *	tracker,
		RoadSearchCloud &	cloud
	)

Definition at line 798 of file RoadSearchCloudMakerAlgorithm.cc.

References abs, RoadSearchCloud::addHit(), edmNew::DetSet< T >::begin(), edmNew::DetSet< T >::end(), edmNew::DetSetVector< T >::end(), edmNew::DetSetVector< T >::find(), TrackingRecHit::geographicalId(), TrackerGeometry::idToDet(), TrackerGeometry::idToDetUnit(), isBarrelSensor(), BaseTrackerRecHit::localPosition(), map_phi(), phi, PV3DBase< T, PVType, FrameType >::phi(), phiFromExtrapolation(), phiMax(), Roads::RPhi, mathSSE::sqrt(), GeomDet::surface(), GeomDet::toGlobal(), Surface::toGlobal(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

                                                                                                                             {
   
   
   unsigned int usedRecHits = 0;
   
   // Get Geometry
   //const PixelTopology *topology = dynamic_cast<const PixelTopology*>(&(tracker->idToDetUnit(id)->topology()));
   
   
   // retrieve vector<SiPixelRecHit> for id
   // loop over SiPixelRecHit
   // check if compatible with cloud, fill into cloud
   
   SiPixelRecHitCollection::const_iterator recHitMatch = inputRecHits->find(id);
   if (recHitMatch == inputRecHits->end()) return usedRecHits;
 
   const SiPixelRecHitCollection::DetSet recHitRange = *recHitMatch;
   
   for ( SiPixelRecHitCollection::DetSet::const_iterator recHitIterator = recHitRange.begin(); 
         recHitIterator != recHitRange.end(); ++recHitIterator) {
     
     const SiPixelRecHit * recHit = &(*recHitIterator);
     
     if ( roadType == Roads::RPhi ) {
       
       if ( isBarrelSensor(id) ) {
         // Barrel Pixel, RoadType RPHI
         
         GlobalPoint ghit = tracker->idToDet(recHit->geographicalId())->surface().toGlobal(recHit->localPosition());
         double hitRadius = sqrt(ghit.x()*ghit.x()+ghit.y()*ghit.y());
         double hitphi = map_phi(ghit.phi());
         double phi = phiFromExtrapolation(d0,phi0,k0,hitRadius,roadType);
         
         if ( std::abs(hitphi-phi) < phiMax(roadType,phi0,k0) ) {
           cloud.addHit(recHit);
           ++usedRecHits;
         }
       } 
       else {
         
         // Forward Pixel,roadtype RPHI
         
         // Get Local Hit Position of the Pixel Hit
         LocalPoint hit = recHit->localPosition();
         
         // Get Phi of hit position 
         double hitphi = map_phi(tracker->idToDetUnit(id)->surface().toGlobal(hit).phi());
         
         // Get Global Hit position
         GlobalPoint ghit = tracker->idToDet(recHit->geographicalId())->surface().toGlobal(recHit->localPosition());
         
         // Get Hit Radis
         double hitRadius = sqrt(ghit.x()*ghit.x()+ghit.y()*ghit.y());
         
         // Get Phi from extrapolation
         double phi = phiFromExtrapolation(d0,phi0,k0,hitRadius,roadType);
         
         if ( std::abs(hitphi-phi) < phiMax(roadType,phi0,k0) ) {
           cloud.addHit(recHit);
           ++usedRecHits;
         }   
       }
     } 
     
     else {
       
       GlobalPoint ghit = tracker->idToDet(recHit->geographicalId())->surface().toGlobal(recHit->localPosition());
       
       double  phi = phiFromExtrapolation(d0,phi0,k0,ghit.z(),roadType);
       if ( (phi - phiMax(roadType,phi0,k0)) < ringPhi && (phi + phiMax(roadType,phi0,k0))>ringPhi ) {
         cloud.addHit(recHit);
         ++usedRecHits;
       }
     }
     
   }
   
   return usedRecHits;
 }

unsigned int RoadSearchCloudMakerAlgorithm::FillRecHitsIntoCloud	(	DetId	id,
		const SiStripRecHit2DCollection *	inputRecHits,
		double	d0,
		double	phi0,
		double	k0,
		Roads::type	roadType,
		double	ringPhi,
		const TrackerGeometry *	tracker,
		RoadSearchCloud &	cloud
	)

unsigned int RoadSearchCloudMakerAlgorithm::FillRecHitsIntoCloudGeneral	(	DetId	id,
		double	d0,
		double	phi0,
		double	k0,
		double	phi1,
		double	k1,
		Roads::type	roadType,
		double	ringPhi,
		const TrackerGeometry *	tracker,
		const SiStripRecHitMatcher *	theHitMatcher,
		RoadSearchCloud &	cloud
	)

Definition at line 492 of file RoadSearchCloudMakerAlgorithm.cc.

References abs, RoadSearchCloud::addHit(), CheckXYIntersection(), CheckZPhiIntersection(), CorrectMatchedHit(), funct::cos(), SiPixelRawToDigiRegional_cfi::deltaPhi, TrackingRecHit::geographicalId(), DetHitAccess::getHitVector(), TrackerGeometry::idToDet(), TrackerGeometry::idToDetUnit(), isBarrelSensor(), isSingleLayer(), BaseTrackerRecHit::localPosition(), StripTopology::localStripLength(), LogDebug, map_phi(), map_phi2(), maxDetHitsInCloudPerDetId, maxRecHitsInCloud_, PV3DBase< T, PVType, FrameType >::perp(), phi, PV3DBase< T, PVType, FrameType >::phi(), phiFromExtrapolation(), phiMax(), PixelSubdetector::PixelBarrel, PixelSubdetector::PixelEndcap, recHitVectorClass, Roads::RPhi, funct::sin(), RoadSearchCloud::size(), mathSSE::sqrt(), StripTopology::strip(), StripTopology::stripAngle(), GeomDet::surface(), funct::tan(), StripSubdetector::TEC, StripSubdetector::TIB, StripSubdetector::TID, StripSubdetector::TOB, GeomDet::toGlobal(), Surface::toGlobal(), GloballyPositioned< T >::toLocal(), GeomDetUnit::topology(), UsePixels, PV3DBase< T, PVType, FrameType >::x(), vdt::x, detailsBasic3DVector::y, PV3DBase< T, PVType, FrameType >::y(), detailsBasic3DVector::z, PV3DBase< T, PVType, FrameType >::z(), and ZPhiDeltaPhi().

Referenced by run().

                                                                                                 {
   
   unsigned int usedRecHits = 0;
   
   bool double_ring_layer = !isSingleLayer(id);
   
   std::vector<TrackingRecHit*> recHitVector = recHitVectorClass.getHitVector(&id);
   
   for ( std::vector<TrackingRecHit*>::const_iterator recHitIterator = recHitVector.begin(); recHitIterator != recHitVector.end(); ++recHitIterator) {
     
     if (      (unsigned int)id.subdetId() == StripSubdetector::TIB 
               || (unsigned int)id.subdetId() == StripSubdetector::TOB 
               || (unsigned int)id.subdetId() == StripSubdetector::TID 
               || (unsigned int)id.subdetId() == StripSubdetector::TEC ) {
       
       const SiStripRecHit2D *recHit = (SiStripRecHit2D*)(*recHitIterator);
       DetId hitId = recHit->geographicalId();
 
       LogDebug("RoadSearch") << "    Testing hit at (x/y/z): " 
                  << tracker->idToDet(recHit->geographicalId())->surface().toGlobal(recHit->localPosition()).x() << " / " 
                  << tracker->idToDet(recHit->geographicalId())->surface().toGlobal(recHit->localPosition()).y() << " / " 
                  << tracker->idToDet(recHit->geographicalId())->surface().toGlobal(recHit->localPosition()).z();
 
       if ( roadType == Roads::RPhi ) {
         if (double_ring_layer && isSingleLayer(hitId)) {
           //
           //  This is where the barrel stereoRecHits end up for Roads::RPhi
           //
           
           // Adjust matched hit for track angle
 
           const GluedGeomDet *theGluedDet = dynamic_cast<const GluedGeomDet*>(tracker->idToDet(hitId));
 
           SiStripMatchedRecHit2D* theCorrectedHit = CorrectMatchedHit(*recHitIterator,theGluedDet,
                                                                         tracker, theHitMatcher,
                                                                         k0, phi0);
           if (theCorrectedHit != 0){
 
             GlobalPoint ghit = tracker->idToDet(theCorrectedHit->geographicalId())->surface().toGlobal(theCorrectedHit->localPosition());
             double hitRadius = sqrt(ghit.x()*ghit.x()+ghit.y()*ghit.y());
             double hitphi = map_phi(ghit.phi());
             double phi = phiFromExtrapolation(d0,phi0,k0,hitRadius,roadType);
             
             float dp = (hitphi-phi);
             float dx = hitRadius*tan(dp);
             
         LogDebug("RoadSearch") << "   Hit phi = " << hitphi << " expected phi = " << phi
                        <<"  dx = " << dx << " for dxMax = " << phiMax(roadType,phi0,k0);
 
             // switch cut to dx instead of dphi
             if ( std::abs(dx) < phiMax(roadType,phi0,k0) ) {
               if ((usedRecHits < maxDetHitsInCloudPerDetId) && (cloud.size() < maxRecHitsInCloud_)) {
                 cloud.addHit(recHit);
                 ++usedRecHits;
               }
             }
             delete theCorrectedHit;
       }
     } else { // Single layer hits here
       if ( isBarrelSensor(hitId) ) {
         //
         //  This is where the barrel rphiRecHits end up for Roads::RPhi
         //
         GlobalPoint ghit = tracker->idToDet(recHit->geographicalId())->surface().toGlobal(recHit->localPosition());
         double hitRadius = sqrt(ghit.x()*ghit.x()+ghit.y()*ghit.y());
         double hitphi = map_phi(ghit.phi());
         double phi = phiFromExtrapolation(d0,phi0,k0,hitRadius,roadType);
               
         float dp = (hitphi-phi);
         float dx = hitRadius*tan(dp);
         LogDebug("RoadSearch") << "   Hit phi = " << hitphi << " expected phi = " << phi
                        <<"  dx = " << dx << " for dxMax = " << phiMax(roadType,phi0,k0);
         // switch cut to dx instead of dphi
         if ( std::abs(dx) < phiMax(roadType,phi0,k0) ) {
               if ((usedRecHits < maxDetHitsInCloudPerDetId) && (cloud.size() < maxRecHitsInCloud_)) {
         cloud.addHit(recHit);
         ++usedRecHits;
           }
         }
       } 
       else {
 
         LocalPoint hit = recHit->localPosition();
         const StripTopology *topology = dynamic_cast<const StripTopology*>(&(tracker->idToDetUnit(hitId)->topology()));
         double stripAngle = topology->stripAngle(topology->strip(hit));
         double stripLength = topology->localStripLength(hit);
 
         LocalPoint innerHitLocal(hit.x()+stripLength/2*std::sin(stripAngle),hit.y()-stripLength/2*std::cos(stripAngle),0);
         LocalPoint outerHitLocal(hit.x()-stripLength/2*std::sin(stripAngle),hit.y()+stripLength/2*std::cos(stripAngle),0);
 
         double innerRadius = tracker->idToDetUnit(hitId)->surface().toGlobal(innerHitLocal).perp(); 
         double outerRadius = tracker->idToDetUnit(hitId)->surface().toGlobal(outerHitLocal).perp();
         double innerExtrapolatedPhi = phiFromExtrapolation(d0,phi0,k0,innerRadius,roadType);
         double outerExtrapolatedPhi = phiFromExtrapolation(d0,phi0,k0,outerRadius,roadType);
 
         GlobalPoint innerHitGlobal =tracker->idToDetUnit(hitId)->surface().toGlobal(innerHitLocal); 
         GlobalPoint outerHitGlobal =tracker->idToDetUnit(hitId)->surface().toGlobal(outerHitLocal); 
 
         GlobalPoint innerRoadGlobal(GlobalPoint::Cylindrical(innerRadius,innerExtrapolatedPhi,
                                  tracker->idToDetUnit(hitId)->surface().toGlobal(hit).z()));
         GlobalPoint outerRoadGlobal(GlobalPoint::Cylindrical(outerRadius,outerExtrapolatedPhi,
                                  tracker->idToDetUnit(hitId)->surface().toGlobal(hit).z()));
 
         LocalPoint innerRoadLocal = tracker->idToDetUnit(hitId)->surface().toLocal(innerRoadGlobal);
         LocalPoint outerRoadLocal = tracker->idToDetUnit(hitId)->surface().toLocal(outerRoadGlobal);
 
         double dxinter = CheckXYIntersection(innerHitLocal, outerHitLocal, 
                            innerRoadLocal, outerRoadLocal);
 
         LogDebug("RoadSearch") << " Hit phi inner = " << innerHitGlobal.phi() << " and outer = " << outerHitGlobal.phi()
                    << " expected inner phi = " << innerExtrapolatedPhi
                    << " and outer phi = "      << outerExtrapolatedPhi
                    <<"  dx = " << dxinter << " for dxMax = " << phiMax(roadType,phi0,k0);
 
         if ( fabs(dxinter) < phiMax(roadType,phi0,k0)) {
           //
           //  This is where the disk rphiRecHits end up for Roads::ZPhi
           //
           if ((usedRecHits < maxDetHitsInCloudPerDetId) && (cloud.size() < maxRecHitsInCloud_)) {
         cloud.addHit(recHit);
         ++usedRecHits;
           }
         }
         //else
           //std::cout<< " ===>>> HIT FAILS !!! " << std::endl;
       }
     } 
       } else {
     //
     // roadType == Roads::ZPhi
     //
     if (double_ring_layer && isSingleLayer(hitId)) {
 
           // Adjust matched hit for track angle
 
           //const SiStripMatchedRecHit2D *theRH = dynamic_cast<SiStripMatchedRecHit2D*>(*recHitIterator);
           const GluedGeomDet *theGluedDet = dynamic_cast<const GluedGeomDet*>(tracker->idToDet(hitId));
 
           SiStripMatchedRecHit2D* theCorrectedHit = CorrectMatchedHit(*recHitIterator,theGluedDet,
                                                                         tracker, theHitMatcher,
                                                                         k1, phi1);
           if (theCorrectedHit != 0){
 
             GlobalPoint ghit = tracker->idToDet(theCorrectedHit->geographicalId())->surface().toGlobal(theCorrectedHit->localPosition());
             double hitphi = map_phi(ghit.phi());
             double hitZ = ghit.z();
         double phi = phiFromExtrapolation(d0,phi0,k0,hitZ,roadType);
             
         float dp = hitphi-phi;
         float dx = hitZ*tan(dp);
               
         //
         //  This is where the disk stereoRecHits end up for Roads::ZPhi
         //
         if ( std::abs(dx) < phiMax(roadType,phi0,k1)) {
           if ((usedRecHits < maxDetHitsInCloudPerDetId) && (cloud.size() < maxRecHitsInCloud_)) {
         cloud.addHit(recHit);
         ++usedRecHits;
           }
         }
         delete theCorrectedHit;
       }
     } else { // Single layer hits here
       if ( isBarrelSensor(hitId) ) {
         //
         //  This is where the barrel (???) rphiRecHits end up for Roads::ZPhi
         //
         LocalPoint hit = recHit->localPosition();
         const StripTopology *topology = dynamic_cast<const StripTopology*>(&(tracker->idToDetUnit(hitId)->topology()));
         double stripAngle = topology->stripAngle(topology->strip(hit));
         double stripLength = topology->localStripLength(hit);
 
         //if (stripAngle!=0) std::cout<<"HEY, WE FOUND A HIT ON A STEREO MODULE!!!" << std::endl;
         // new method
         LocalPoint innerHitLocal(hit.x()+stripLength/2*std::sin(stripAngle),hit.y()-stripLength/2*std::cos(stripAngle),0);
         LocalPoint outerHitLocal(hit.x()-stripLength/2*std::sin(stripAngle),hit.y()+stripLength/2*std::cos(stripAngle),0);
         double innerHitPhi = map_phi(tracker->idToDetUnit(hitId)->surface().toGlobal(innerHitLocal).phi()); 
         double outerHitPhi = map_phi(tracker->idToDetUnit(hitId)->surface().toGlobal(outerHitLocal).phi());
         double innerHitZ = tracker->idToDetUnit(hitId)->surface().toGlobal(innerHitLocal).z(); 
         double outerHitZ = tracker->idToDetUnit(hitId)->surface().toGlobal(outerHitLocal).z();
         double innerExtrapolatedPhi = phiFromExtrapolation(d0,phi0,k0,innerHitZ,roadType);
         double outerExtrapolatedPhi = phiFromExtrapolation(d0,phi0,k0,outerHitZ,roadType);
 
         double midPointZ = 0.5*(innerHitZ+outerHitZ);
 
         double dPhiInter = CheckZPhiIntersection(innerHitPhi, innerHitZ, outerHitPhi, outerHitZ, 
                            innerExtrapolatedPhi, innerHitZ,
                            outerExtrapolatedPhi, outerHitZ);
 
         double dX = midPointZ*tan(dPhiInter);
 
         if (std::abs(dX) < 1.5*phiMax(roadType,phi0,k1)) {
           if ((usedRecHits < maxDetHitsInCloudPerDetId) && (cloud.size() < maxRecHitsInCloud_)) {
         cloud.addHit(recHit);
         ++usedRecHits;
           }
         }
 
       } else {
 
         //
         //  This is where the disk rphiRecHits end up for Roads::ZPhi
         //
         LocalPoint hit = recHit->localPosition();
         const StripTopology *topology = dynamic_cast<const StripTopology*>(&(tracker->idToDetUnit(hitId)->topology()));
         double stripAngle = topology->stripAngle(topology->strip(hit));
         double stripLength = topology->localStripLength(hit);
         // new method
         double hitZ = tracker->idToDetUnit(hitId)->surface().toGlobal(hit).z(); 
         double extrapolatedPhi = phiFromExtrapolation(d0,phi0,k0,hitZ,roadType);
 
         LocalPoint innerHitLocal(hit.x()+stripLength/2*std::sin(stripAngle),hit.y()-stripLength/2*std::cos(stripAngle),0);
         LocalPoint outerHitLocal(hit.x()-stripLength/2*std::sin(stripAngle),hit.y()+stripLength/2*std::cos(stripAngle),0);
 
         double innerHitPhi = map_phi(tracker->idToDetUnit(hitId)->surface().toGlobal(innerHitLocal).phi()); 
         double outerHitPhi = map_phi(tracker->idToDetUnit(hitId)->surface().toGlobal(outerHitLocal).phi());
         //double innerZ = tracker->idToDetUnit(hitId)->surface().toGlobal(innerHitLocal).z(); 
         //double outerZ = tracker->idToDetUnit(hitId)->surface().toGlobal(outerHitLocal).z();
         //if (innerZ != outerZ)  std::cout<<"HEY!!! innerZ = " << innerZ << " != outerZ = " << outerZ << std::endl;
 
         double deltaPhi = ZPhiDeltaPhi(innerHitPhi,outerHitPhi,extrapolatedPhi);
         double deltaX   =  hitZ*tan(deltaPhi);
         if (std::abs(deltaX) < phiMax(roadType,phi0,k1)){
           if ((usedRecHits < maxDetHitsInCloudPerDetId) && (cloud.size() < maxRecHitsInCloud_)) {
         cloud.addHit(recHit);
         ++usedRecHits;
           }
         }
       }
     } 
       }
     } else if ( (unsigned int)id.subdetId() == PixelSubdetector::PixelBarrel 
                 || (unsigned int)id.subdetId() == PixelSubdetector::PixelEndcap) {
       if ( UsePixels ) {
         
         const SiPixelRecHit *recHit = (SiPixelRecHit*)(*recHitIterator);
         
         if ( roadType == Roads::RPhi ) {
           
           if ( isBarrelSensor(id) ) {
             // Barrel Pixel, RoadType RPHI
             
             GlobalPoint ghit = tracker->idToDet(recHit->geographicalId())->surface().toGlobal(recHit->localPosition());
             double hitRadius = sqrt(ghit.x()*ghit.x()+ghit.y()*ghit.y());
             double hitphi = map_phi(ghit.phi());
             double phi = phiFromExtrapolation(d0,phi0,k0,hitRadius,roadType);
             
             float dp = hitphi-phi;
             float dx = hitRadius*tan(dp);
             
             // switch cut to dx instead of dphi
             if ( std::abs(dx) < phiMax(roadType,phi0,k0) ) {
               cloud.addHit(recHit);
               ++usedRecHits;
             }
       } else {
               
         // Forward Pixel,roadtype RPHI
               
         // Get Local Hit Position of the Pixel Hit
         LocalPoint hit = recHit->localPosition();
               
         // Get Phi of hit position 
         double hitphi = map_phi(tracker->idToDetUnit(id)->surface().toGlobal(hit).phi());
               
         // Get Global Hit position
         GlobalPoint ghit = tracker->idToDet(recHit->geographicalId())->surface().toGlobal(recHit->localPosition());
               
         // Get Hit Radis
         double hitRadius = sqrt(ghit.x()*ghit.x()+ghit.y()*ghit.y());
               
         // Get Phi from extrapolation
         double phi = phiFromExtrapolation(d0,phi0,k0,hitRadius,roadType);
               
         if ( std::abs(hitphi-phi) < phiMax(roadType,phi0,k0) ) {
           cloud.addHit(recHit);
           ++usedRecHits;
         }   
       }
     } else {
             
       GlobalPoint ghit = tracker->idToDet(recHit->geographicalId())->surface().toGlobal(recHit->localPosition());            
       double  phi = phiFromExtrapolation(d0,phi0,k0,ghit.z(),roadType);            
       double hitphi = map_phi(ghit.phi());
       double dphi = map_phi2(hitphi-phi);
       float dx = ghit.z()*tan(dphi);
             
       if ( std::abs(dx) < 0.25 ) {
         cloud.addHit(recHit);
         ++usedRecHits;
       }
     }
       }
     } else {
       edm::LogError("RoadSearch") << "recHitVector from general hit access function contains unknown detector id: " << (unsigned int)id.subdetId() << " rawId: " << id.rawId();
     }
     
   } //for loop over all recHits
   
   
   return usedRecHits;
 }

bool RoadSearchCloudMakerAlgorithm::isBarrelSensor ( DetId id )

Definition at line 909 of file RoadSearchCloudMakerAlgorithm.cc.

References PixelSubdetector::PixelBarrel, StripSubdetector::TIB, and StripSubdetector::TOB.

Referenced by FillPixRecHitsIntoCloud(), and FillRecHitsIntoCloudGeneral().

                                                            {
   
   if ( (unsigned int)id.subdetId() == StripSubdetector::TIB ) {
     return true;
   } else if ( (unsigned int)id.subdetId() == StripSubdetector::TOB ) {
     return true;
   } else if ( (unsigned int)id.subdetId() == PixelSubdetector::PixelBarrel ) {
     return true;
   } else {
     return false;
   }
   
 }

bool RoadSearchCloudMakerAlgorithm::isSingleLayer ( DetId id )

Definition at line 880 of file RoadSearchCloudMakerAlgorithm.cc.

References SiStripDetId::glued(), StripSubdetector::TEC, StripSubdetector::TIB, StripSubdetector::TID, and StripSubdetector::TOB.

Referenced by FillRecHitsIntoCloudGeneral().

                                                           {
   
   if ( (unsigned int)id.subdetId() == StripSubdetector::TIB ) {
     TIBDetId tibid(id.rawId()); 
     if ( !tibid.glued() ) {
       return true;
     }
   } else if ( (unsigned int)id.subdetId() == StripSubdetector::TOB ) {
     TOBDetId tobid(id.rawId()); 
     if ( !tobid.glued() ) {
       return true;
     }
   } else if ( (unsigned int)id.subdetId() == StripSubdetector::TID ) {
     TIDDetId tidid(id.rawId()); 
     if ( !tidid.glued() ) {
       return true;
     }
   } else if ( (unsigned int)id.subdetId() == StripSubdetector::TEC ) {
     TECDetId tecid(id.rawId()); 
     if ( !tecid.glued() ) {
       return true;
     }
   } else {
     return false;
   }
   
   return false;
 }

void RoadSearchCloudMakerAlgorithm::makecircle	(	double	x1_cs,
		double	y1_cs,
		double	x2_cs,
		double	y2_cs,
		double	x3_cs,
		double	y3_cs
	)

Definition at line 980 of file RoadSearchCloudMakerAlgorithm.cc.

References d0h, epsilon, python.connectstrParser::f1, python.connectstrParser::f2, python.connectstrParser::f3, Geom::halfPi(), omegah, phi0h, rho, and mathSSE::sqrt().

Referenced by run().

                                                                                          {
   double x1t=x1-x3; double y1t=y1-y3; double r1s=x1t*x1t+y1t*y1t;
   double x2t=x2-x3; double y2t=y2-y3; double r2s=x2t*x2t+y2t*y2t;
   double rho=x1t*y2t-x2t*y1t;
   double xc, yc, rc, fac;
   if (fabs(rho)<RoadSearchCloudMakerAlgorithm::epsilon){
     rc=1.0/(RoadSearchCloudMakerAlgorithm::epsilon);
     fac=sqrt(x1t*x1t+y1t*y1t);
     xc=x2+y1t*rc/fac;
     yc=y2-x1t*rc/fac;
   }else{
     fac=0.5/rho;
     xc=fac*(r1s*y2t-r2s*y1t);
     yc=fac*(r2s*x1t-r1s*x2t); 
     rc=sqrt(xc*xc+yc*yc); xc+=x3; yc+=y3;
   }
   double s3=0.0;
   double f1=x1*yc-y1*xc; double f2=x2*yc-y2*xc; 
   double f3=x3*yc-y3*xc;
   if ((f1<0.0)&&(f2<0.0)&&(f3<=0.0))s3=1.0;
   if ((f1>0.0)&&(f2>0.0)&&(f3>=0.0))s3=-1.0;
   d0h=-s3*(sqrt(xc*xc+yc*yc)-rc);
   phi0h=atan2(yc,xc)+s3*Geom::halfPi();
   omegah=-s3/rc;
 }

double RoadSearchCloudMakerAlgorithm::map_phi ( double phi )

Definition at line 135 of file RoadSearchCloudMakerAlgorithm.cc.

References phi, query::result, and Geom::twoPi().

Referenced by FillPixRecHitsIntoCloud(), FillRecHitsIntoCloudGeneral(), phiFromExtrapolation(), and run().

                                                         {
   // map phi to [0,2pi]
   double result = phi;
   if ( result < -1.0*Geom::twoPi()) result = result + Geom::twoPi();
   if ( result < 0)                 result = Geom::twoPi() + result;
   if ( result > Geom::twoPi())     result = result - Geom::twoPi();
   return result;
 }

double RoadSearchCloudMakerAlgorithm::map_phi2 ( double phi )

Definition at line 144 of file RoadSearchCloudMakerAlgorithm.cc.

References phi, Geom::pi(), query::result, and Geom::twoPi().

Referenced by CheckZPhiIntersection(), CorrectMatchedHit(), FillRecHitsIntoCloudGeneral(), run(), and ZPhiDeltaPhi().

                                                          {
   // map phi to [-pi,pi]
   double result = phi;
   if ( result < 1.0*Geom::pi() ) result = result + Geom::twoPi();
   if ( result >= Geom::pi())  result = result - Geom::twoPi();
   return result;
 }

double RoadSearchCloudMakerAlgorithm::phiFromExtrapolation	(	double	d0,
		double	phi0,
		double	k0,
		double	ringRadius,
		Roads::type	roadType
	)

Definition at line 923 of file RoadSearchCloudMakerAlgorithm.cc.

References funct::cos(), reco::ParticleMasses::k0, prof2calltree::l, map_phi(), Roads::RPhi, funct::sin(), and mathSSE::sqrt().

Referenced by FillPixRecHitsIntoCloud(), FillRecHitsIntoCloudGeneral(), and run().

                                                                                                                                    {
   
   double ringPhi = -99.;
   if ( roadType == Roads::RPhi ) {
     double omega=k0, rl=ringRadius;
     double sp0=sin(phi0); double cp0=cos(phi0);  
     if (fabs(omega)>0.000005){
       double xc=-sp0*(d0+1.0/omega);    
       double yc=cp0*(d0+1.0/omega);
       double rh=fabs(1.0/omega);
       double bbb=fabs(d0+1.0/omega);
       double sss=0.5*(rl+rh+bbb);
       double ddd=sqrt((sss-bbb)*(sss-rh)/(sss*(sss-rl)));
       double phil1=2.0*atan(ddd);
       double phit=phi0+phil1; if (omega<0.0)phit=phi0-phil1;
       double xh=xc+sin(phit)/omega;
       double yh=yc-cos(phit)/omega;
       double phih=atan2(yh,xh);
       ringPhi = map_phi(phih);
     }
     else {
       double cee = rl*rl - d0*d0 -0.25*omega*omega - omega*d0;
       if (cee<0.0){return ringPhi;}
       double l = sqrt(cee);
       double xh=-sp0*d0+l*cp0-0.5*l*l*omega*sp0;
       double yh= cp0*d0+l*sp0+0.5*l*l*omega*cp0;
       double phih=atan2(yh,xh);
       ringPhi = map_phi(phih);     
     }
   } 
   else {
     ringPhi = map_phi(phi0 + k0 * ringRadius);
   }
   
   return ringPhi;
 }

double RoadSearchCloudMakerAlgorithm::phiMax	(	Roads::type	roadType,
		double	phi0,
		double	k0
	)

Definition at line 960 of file RoadSearchCloudMakerAlgorithm.cc.

References Roads::RPhi, theRPhiRoadSize, theZPhiRoadSize, and Roads::ZPhi.

Referenced by FillPixRecHitsIntoCloud(), and FillRecHitsIntoCloudGeneral().

                                                                      {
   
   double dphi;
   if ( roadType == Roads::RPhi ) {
     // switch cut to dx instead of dphi
     // Still call this dphi, but will now be dx
     dphi = theRPhiRoadSize + 0.15*82.0*fabs(k0);
   }
   else if ( roadType == Roads::ZPhi ) {
     dphi = theZPhiRoadSize + 0.4*82.0*fabs(k0);
   }
   else {
     edm::LogWarning("RoadSearch") << "Bad roadType: "<< roadType;
     dphi = theZPhiRoadSize;
   }
   return dphi;
   
 }

void RoadSearchCloudMakerAlgorithm::run	(	edm::Handle< RoadSearchSeedCollection >	input,
		const SiStripRecHit2DCollection *	rphiRecHits,
		const SiStripRecHit2DCollection *	stereoRecHits,
		const SiStripMatchedRecHit2DCollection *	matchedRecHits,
		const SiPixelRecHitCollection *	pixRecHits,
		const edm::EventSetup &	es,
		RoadSearchCloudCollection &	output
	)

Runs the algorithm.

Definition at line 152 of file RoadSearchCloudMakerAlgorithm.cc.

References abs, alpha, begin, Clean(), funct::cos(), cond::rpcobgas::detid, doCleaning_, relativeConstraints::empty, PV3DBase< T, PVType, FrameType >::eta(), FillRecHitsIntoCloudGeneral(), TrackingRecHit::geographicalId(), edm::EventSetup::get(), increaseMaxNumberOfConsecutiveMissedLayersPerCloud, increaseMaxNumberOfMissedLayersPerCloud, reco::ParticleMasses::k0, TrackingRecHit::localPosition(), LogDebug, makecircle(), map_phi(), map_phi2(), maxFractionOfConsecutiveMissedLayersPerCloud, maxFractionOfMissedLayersPerCloud, minFractionOfUsedLayersPerCloud, netabin, NoFieldCosmic, omegah, PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::phi(), phi0h, phiFromExtrapolation(), Geom::pi(), edm::Handle< T >::product(), edm::ESHandle< class >::product(), recHitVectorClass, relativeConstraints::ring, RoadMapESSource_cfi::roads, roadsLabel_, Roads::RPhi, DetHitAccess::setCollections(), DetHitAccess::setMode(), funct::sin(), RoadSearchCloud::size(), mathSSE::sqrt(), DetHitAccess::standard, groupFilesInBlocks::temp, theMinimumHalfRoad, patCandidatesForDimuonsSequences_cff::tracker, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by RoadSearchCloudMaker::produce().

 {  
   // intermediate arrays for storing clouds for cleaning
   const int nphibin = 24;
   const int netabin = 24;
   RoadSearchCloudCollection CloudArray[nphibin][netabin];
   
   // get roads
   edm::ESHandle<Roads> roads;
   es.get<RoadMapRecord>().get(roadsLabel_, roads);
   
   // get RoadSearchSeed collection
   const RoadSearchSeedCollection* inputSeeds = input.product();
   
   // set collections for general hit access method
   recHitVectorClass.setCollections(rphiRecHits,stereoRecHits,matchedRecHits,pixRecHits);
   recHitVectorClass.setMode(DetHitAccess::standard);
   
   // get tracker geometry
   edm::ESHandle<TrackerGeometry> tracker;
   es.get<TrackerDigiGeometryRecord>().get(tracker);
   
   // get hit matcher
   SiStripRecHitMatcher* theHitMatcher = new SiStripRecHitMatcher(3.0);
 
   edm::LogInfo("RoadSearch") << "Found " << inputSeeds->size() << " input seeds.";   
   // loop over seeds
   for ( RoadSearchSeedCollection::const_iterator seed = inputSeeds->begin(); seed != inputSeeds->end(); ++seed) {
     
     const Roads::RoadSeed *roadSeed = seed->getSeed();
       
     if ( roadSeed == 0 ) {
       edm::LogWarning("RoadSearch") << "RoadSeed could not be resolved from RoadSearchSeed hits, discard seed!";
     } else {
         
       Roads::type roadType = roads->getRoadType(roadSeed);
       if (NoFieldCosmic) roadType = Roads::RPhi;
         
       // fixme: from here on, calculate with 1st and 3rd seed hit (inner and outer of initial circle)
       // fixme: adapt to new seed structure
         
       // get global positions of the hits, calculate before Road lookup to be used
       const TrackingRecHit* innerSeedRingHit = (*(seed->begin()));
       const TrackingRecHit* outerSeedRingHit = (*(seed->end() - 1));
         
       GlobalPoint innerSeedHitGlobalPosition = tracker->idToDet(innerSeedRingHit->geographicalId())->surface().toGlobal(innerSeedRingHit->localPosition());
       GlobalPoint outerSeedHitGlobalPosition = tracker->idToDet(outerSeedRingHit->geographicalId())->surface().toGlobal(outerSeedRingHit->localPosition());
  
       LogDebug("RoadSearch") << "Seed # " <<seed-inputSeeds->begin() << " inner hit (x/y/z): "
                  << innerSeedHitGlobalPosition.x() << " / "
                  << innerSeedHitGlobalPosition.y() << " / "
                  << innerSeedHitGlobalPosition.z();
       LogDebug("RoadSearch") << "Seed # " <<seed-inputSeeds->begin() << " outer hit (x/y/z): "
                  << outerSeedHitGlobalPosition.x() << " / "
                  << outerSeedHitGlobalPosition.y() << " / "
                  << outerSeedHitGlobalPosition.z();
       
       LogDebug("RoadSearch") << "Seed # " <<seed-inputSeeds->begin() << " inner hit (r/phi): "
                  << innerSeedHitGlobalPosition.perp() << " / "
                  << innerSeedHitGlobalPosition.phi();
       LogDebug("RoadSearch") << "Seed # " <<seed-inputSeeds->begin() << " outer hit (r/phi): "
                  << outerSeedHitGlobalPosition.perp() << " / "
                  << outerSeedHitGlobalPosition.phi();
 
 
       // extrapolation parameters, phio: [0,2pi]
       double d0 = 0.0;
       double phi0 = -99.;
       double k0   = -99999999.99;
       
       double phi1 = -99.;
       double k1   = -99999999.99;
       // get bins of eta and phi of outer seed hit;
       
       double outer_phi = map_phi(outerSeedHitGlobalPosition.phi());
       double outer_eta = outerSeedHitGlobalPosition.eta();
       
       int phibin = (int)(nphibin*(outer_phi/(2*Geom::pi())));
       int etabin = (int)(netabin*(outer_eta+3.0)/6.0);
         
       // calculate phi0 and k0 dependent on RoadType
       if ( roadType == Roads::RPhi ) {
     double dr = outerSeedHitGlobalPosition.perp() - innerSeedHitGlobalPosition.perp();
     const double dr_min = 1; // cm
     if ( dr < dr_min ) {
       edm::LogWarning("RoadSearch") << "RPhi road: seed Hits distance smaller than 1 cm, do not consider this seed.";
     } else {
       // calculate r-phi extrapolation: phi = phi0 + asin(k0 r)
       double det = innerSeedHitGlobalPosition.x() * outerSeedHitGlobalPosition.y() - innerSeedHitGlobalPosition.y() * outerSeedHitGlobalPosition.x();
       if ( det == 0 ) {
         edm::LogWarning("RoadSearch") << "RPhi road: 'det' == 0, do not consider this seed.";
       } else {
         double x0=0.0; double y0=0.0;
         double innerx=innerSeedHitGlobalPosition.x();
         double innery=innerSeedHitGlobalPosition.y();
         double outerx=outerSeedHitGlobalPosition.x();
         double outery=outerSeedHitGlobalPosition.y();
 
         if (NoFieldCosmic){
           phi0=atan2(outery-innery,outerx-innerx);
           double alpha=atan2(innery,innerx);
           double d1=sqrt(innerx*innerx+innery*innery);
           d0=d1*sin(alpha-phi0); x0=-d0*sin(phi0); y0=d0*cos(phi0); k0=0.0;
         }else{
           makecircle(innerx,innery,outerx,outery,x0,y0);
           phi0 = phi0h;
           k0 = omegah;
         }              
         LogDebug("RoadSearch") << "Seed # " <<seed-inputSeeds->begin() << " trajectory parameters: d0 = "<< d0 << " phi0 = " << phi0;
       }
     }
       } else {
     double dz = outerSeedHitGlobalPosition.z() - innerSeedHitGlobalPosition.z();
     const double dz_min = 1.e-6; // cm;
     if ( std::abs(dz) < dz_min ) {
       edm::LogWarning("RoadSearch") << "ZPhi road: seed Hits are less than .01 microns away in z, do not consider this seed.";
     } else {
       // calculate z-phi extrapolation: phi = phi0 + k0 z
       k0 = map_phi2(outerSeedHitGlobalPosition.phi() - innerSeedHitGlobalPosition.phi()) / dz;
       phi0 =  map_phi(innerSeedHitGlobalPosition.phi() - k0 * innerSeedHitGlobalPosition.z());
       
       // get approx pt for use in correcting matched hits
       makecircle(innerSeedHitGlobalPosition.x(),innerSeedHitGlobalPosition.y(),
              outerSeedHitGlobalPosition.x(),outerSeedHitGlobalPosition.y(),
              0.0,0.0); // x0,y0 = 0.0 for now
       phi1 = phi0h;
       k1 = omegah;
     }
       }
         
       // continue if valid extrapolation parameters have been found
       if ( (phi0 != -99.) && (k0 != -99999999.99) ) {
     const Roads::RoadSet *roadSet = seed->getSet();
 
     // create cloud
     RoadSearchCloud cloud;
           
     bool firstHitFound = false;
     unsigned int layerCounter = 0;
     unsigned int usedLayers = 0;
     unsigned int missedLayers = 0;
     unsigned int consecutiveMissedLayers = 0;
 
     unsigned int totalLayers = roadSet->size();
 
     // caluclate minNumberOfUsedLayersPerCloud, maxNumberOfMissedLayersPerCloud and maxNumberOfConsecutiveMissedLayersPerCloud 
     // by rounding to integer minFractionOfUsedLayersPerCloud. maxFractionOfMissedLayersPerCloud and maxFractionOfConsecutiveMissedLayersPerCloud
     unsigned int minNumberOfUsedLayersPerCloud = static_cast<unsigned int>(totalLayers * minFractionOfUsedLayersPerCloud + 0.5);
     if (minNumberOfUsedLayersPerCloud < 3) minNumberOfUsedLayersPerCloud = 3;
     unsigned int maxNumberOfMissedLayersPerCloud = static_cast<unsigned int>(totalLayers * maxFractionOfMissedLayersPerCloud + 0.5);
     unsigned int maxNumberOfConsecutiveMissedLayersPerCloud = static_cast<unsigned int>(totalLayers * maxFractionOfConsecutiveMissedLayersPerCloud + 0.5);
 
     // increase consecutive layer cuts between 0.9 and 1.5
     if (std::abs(outer_eta) > 0.9 && std::abs(outer_eta) < 1.5) {
       maxNumberOfConsecutiveMissedLayersPerCloud += increaseMaxNumberOfConsecutiveMissedLayersPerCloud;
       maxNumberOfMissedLayersPerCloud += increaseMaxNumberOfMissedLayersPerCloud;
     }
           
     for ( Roads::RoadSet::const_iterator roadSetVector = roadSet->begin();
           roadSetVector != roadSet->end();
           ++roadSetVector ) {
             
       ++layerCounter;
       unsigned int usedHitsInThisLayer = 0;
       bool intersectsLayer = false;
             
       for ( std::vector<const Ring*>::const_iterator ring = roadSetVector->begin(); ring != roadSetVector->end(); ++ring ) {
               
         // calculate phi-range for lookup of DetId's in Rings of RoadSet
         // calculate phi at radius of Ring using extrapolation, Ring radius average of Ring.rmin() and Ring.rmax()
         double ringRadius = (*ring)->getrmin() + ((*ring)->getrmax()-(*ring)->getrmin())/2;
         double ringZ      = (*ring)->getzmin() + ((*ring)->getzmax()-(*ring)->getzmin())/2;
         double ringPhi = 0.0;
         if ( roadType == Roads::RPhi ) {
           ringPhi = phiFromExtrapolation(d0,phi0,k0,ringRadius,roadType);
         } else {
           ringPhi = phiFromExtrapolation(d0,phi0,k0,ringZ,roadType);
         }
         if (ringPhi == -99) continue;
         intersectsLayer = true;
 
         LogDebug("RoadSearch") << "Seed # " <<seed-inputSeeds->begin() << " testing ring at R = " << ringRadius
                    << " Z = " << ringZ << " ringPhi = " << ringPhi;
               
         int nDetIds = (*ring)->getNumDetIds();
         double theHalfRoad = theMinimumHalfRoad*(2.0*Geom::pi())/((double)nDetIds);
         // calculate range in phi around ringPhi
         double upperPhiRangeBorder = map_phi2(ringPhi + theHalfRoad);
         double lowerPhiRangeBorder = map_phi2(ringPhi - theHalfRoad);
 
         if ( lowerPhiRangeBorder <= upperPhiRangeBorder ) {
                 
           for ( Ring::const_iterator detid = (*ring)->lower_bound(lowerPhiRangeBorder); detid != (*ring)->upper_bound(upperPhiRangeBorder); ++detid) {
         usedHitsInThisLayer += FillRecHitsIntoCloudGeneral(detid->second,d0,phi0,k0,phi1,k1,roadType,ringPhi,
                                    tracker.product(),theHitMatcher,cloud);
           }
                 
         } else {
           for ( Ring::const_iterator detid = (*ring)->lower_bound(lowerPhiRangeBorder); detid != (*ring)->end(); ++detid) {
         usedHitsInThisLayer += FillRecHitsIntoCloudGeneral(detid->second,d0,phi0,k0,phi1,k1,roadType,ringPhi,
                                    tracker.product(),theHitMatcher,cloud);
           }
                 
           for ( Ring::const_iterator detid = (*ring)->begin(); detid != (*ring)->upper_bound(upperPhiRangeBorder); ++detid) {
         usedHitsInThisLayer += FillRecHitsIntoCloudGeneral(detid->second,d0,phi0,k0,phi1,k1,roadType,ringPhi,
                                    tracker.product(),theHitMatcher,cloud);
           }
         }
       LogDebug("RoadSearch") << "Seed # " <<seed-inputSeeds->begin() << " now has " << usedHitsInThisLayer << "  hits in ring at R = " << ringRadius
                  << " Z = " << ringZ << " ringPhi = " << ringPhi;
       }
             
       if ( !firstHitFound ) {
         if ( usedHitsInThisLayer > 0 ) {
 
           firstHitFound = true;
 
           // reset totalLayers according to first layer with hit
           totalLayers = roadSet->size() - layerCounter + 1;
 
           // re-caluclate minNumberOfUsedLayersPerCloud, maxNumberOfMissedLayersPerCloud and maxNumberOfConsecutiveMissedLayersPerCloud 
           // by rounding to integer minFractionOfUsedLayersPerCloud. maxFractionOfMissedLayersPerCloud and maxFractionOfConsecutiveMissedLayersPerCloud
           minNumberOfUsedLayersPerCloud = static_cast<unsigned int>(totalLayers * minFractionOfUsedLayersPerCloud + 0.5);
           if (minNumberOfUsedLayersPerCloud < 3) minNumberOfUsedLayersPerCloud = 3;
           maxNumberOfMissedLayersPerCloud = static_cast<unsigned int>(totalLayers * maxFractionOfMissedLayersPerCloud + 0.5);
           maxNumberOfConsecutiveMissedLayersPerCloud = static_cast<unsigned int>(totalLayers * maxFractionOfConsecutiveMissedLayersPerCloud + 0.5);
 
           // increase consecutive layer cuts between 0.9 and 1.5
           if (std::abs(outer_eta) > 0.9 && std::abs(outer_eta) < 1.5) {
         maxNumberOfConsecutiveMissedLayersPerCloud += increaseMaxNumberOfConsecutiveMissedLayersPerCloud;
         maxNumberOfMissedLayersPerCloud += increaseMaxNumberOfMissedLayersPerCloud;
           }
 
           ++usedLayers;
           consecutiveMissedLayers = 0;
 
         }
       } else {
         if (intersectsLayer){
           if ( usedHitsInThisLayer > 0 ) {
         ++usedLayers;
         consecutiveMissedLayers = 0;
           } else {
         ++ missedLayers;
         ++consecutiveMissedLayers;
           }
         }
         LogDebug("RoadSearch") << "Seed # "<<seed-inputSeeds->begin() << " Layer info: " 
                    << " totalLayers: " << totalLayers 
                    << " usedLayers: " << usedLayers 
                    << " missedLayers: " << missedLayers
                    << " consecutiveMissedLayers: " << consecutiveMissedLayers;
 
         // break condition, hole larger than maxNumberOfConsecutiveMissedLayersPerCloud
         if ( consecutiveMissedLayers > maxNumberOfConsecutiveMissedLayersPerCloud ) {
           LogDebug("RoadSearch") << "BREAK: seed # "<<seed-inputSeeds->begin() 
                      << " More than " << maxNumberOfConsecutiveMissedLayersPerCloud << " missed consecutive layers!";
           break;
         }
 
         // break condition, already  missed too many layers
         if ( missedLayers > maxNumberOfMissedLayersPerCloud ) {
           LogDebug("RoadSearch") << "BREAK: seed # "<<seed-inputSeeds->begin() 
                      << " More than " << maxNumberOfMissedLayersPerCloud << " missed layers!";
           break;
         }
 
         // break condition, cannot satisfy minimal number of used layers
         if ( totalLayers-missedLayers < minNumberOfUsedLayersPerCloud ) {
           LogDebug("RoadSearch") << "BREAK: seed # "<<seed-inputSeeds->begin() 
                      << " Cannot satisfy at least " << minNumberOfUsedLayersPerCloud << " used layers!";
           break;
         }
       }   
             
     }
 
     if ( consecutiveMissedLayers <= maxNumberOfConsecutiveMissedLayersPerCloud ) {
       if ( usedLayers >= minNumberOfUsedLayersPerCloud ) {
         if ( missedLayers <= maxNumberOfMissedLayersPerCloud ) {
 
           CloudArray[phibin][etabin].push_back(cloud);
 
           if ( roadType == Roads::RPhi ){ 
         LogDebug("RoadSearch") << "This r-phi seed # "<<seed-inputSeeds->begin()
                        <<" yields a cloud with " <<cloud.size() <<" hits on " << usedLayers << " layers out of " << totalLayers;
           } else {
         LogDebug("RoadSearch") << "This z-phi seed # "<<seed-inputSeeds->begin()
                        <<" yields a cloud with " <<cloud.size() <<" hits on " << usedLayers << " layers out of " << totalLayers;
           }
         } else {
           LogDebug("RoadSearch") << "Missed layers: " << missedLayers << " More than " << maxNumberOfMissedLayersPerCloud << " missed layers!";
           if ( roadType == Roads::RPhi ){ 
         LogDebug("RoadSearch") << "This r-phi seed # "<<seed-inputSeeds->begin() <<" yields no clouds";
           } else {
         LogDebug("RoadSearch") << "This z-phi seed # "<<seed-inputSeeds->begin() <<" yields no clouds";
           }
         }
       }
       else {
         LogDebug("RoadSearch") << "Seed # "<<seed-inputSeeds->begin() <<" fails: used layers = " << usedLayers << " < " << minNumberOfUsedLayersPerCloud;
       }
     }
     else {
       LogDebug("RoadSearch") << "Seed # "<<seed-inputSeeds->begin() <<" fails: consecutive missed layers = " << consecutiveMissedLayers << " > " << maxNumberOfConsecutiveMissedLayersPerCloud;
     }
       }
     }
   }
   
   // Loop for initial cleaning
   for (int iphi=0; iphi<nphibin; ++iphi){
     for (int ieta=0; ieta<netabin; ++ieta){
       if (!CloudArray[iphi][ieta].empty()) {
         if (doCleaning_){
           RoadSearchCloudCollection temp = Clean(&CloudArray[iphi][ieta]);
           for ( RoadSearchCloudCollection::const_iterator ic = temp.begin(); ic!=temp.end(); ++ic)
             output.push_back(*ic);
         }
         else 
           for ( RoadSearchCloudCollection::const_iterator ic = CloudArray[iphi][ieta].begin(); 
                 ic!=CloudArray[iphi][ieta].end(); ++ic)
             output.push_back(*ic);
       }
     }
   }
 
   delete theHitMatcher;
   edm::LogInfo("RoadSearch") << "Found " << output.size() << " clouds."; 
   for ( RoadSearchCloudCollection::const_iterator ic = output.begin(); ic!=output.end(); ++ic)
     edm::LogInfo("RoadSearch") << "    Cloud " << ic-output.begin()<< " has " << ic->size() << " hits."; 
   
 }

double RoadSearchCloudMakerAlgorithm::ZPhiDeltaPhi	(	double	phi1,
		double	phi2,
		double	phiExpect
	)

Definition at line 1091 of file RoadSearchCloudMakerAlgorithm.cc.

References abs, SiPixelRawToDigiRegional_cfi::deltaPhi, and map_phi2().

Referenced by FillRecHitsIntoCloudGeneral().

                                                                                                      {
 
   double deltaPhi = -999.;
 
   double dPhiHits = map_phi2(hitPhi1-hitPhi2);
   double dPhi1 = map_phi2(hitPhi1-predictedPhi);
   double dPhi2 = map_phi2(hitPhi2-predictedPhi);
 
   if (dPhiHits >= 0){  // hitPhi1 >= hitPhi2
     if ( (dPhi1>=0.0) && (dPhi2 <= 0.0))
       deltaPhi = 0.0;
     else{
       if (std::abs(dPhi1)<std::abs(dPhi2))
     deltaPhi = dPhi1;
       else
     deltaPhi = dPhi2;
     }
   }
   else { // hitPhi1 < hitPhi2
     if ( (dPhi1<=0.0) && (dPhi2 >= 0.0))
       deltaPhi = 0.0;
     else{
       if (std::abs(dPhi1)<std::abs(dPhi2))
     deltaPhi = dPhi1;
       else
     deltaPhi = dPhi2;
     }
   }
 
   return deltaPhi;
 
 }