#include <RecoEgamma/EgammaElectronAlgos/interface/SiStripElectronAlgo.h>

Public Member Functions
bool	findElectron (reco::SiStripElectronCollection &electronOut, TrackCandidateCollection &trackCandidateOut, const reco::SuperClusterRef &superclusterIn, const TrackerTopology *tTopo)

void	prepareEvent (const edm::ESHandle< TrackerGeometry > &tracker, const edm::Handle< SiStripRecHit2DCollection > &rphiHits, const edm::Handle< SiStripRecHit2DCollection > &stereoHits, const edm::Handle< SiStripMatchedRecHit2DCollection > &matchedHits, const edm::ESHandle< MagneticField > &magneticField)

	SiStripElectronAlgo (unsigned int maxHitsOnDetId, double originUncertainty, double phiBandWidth, double maxNormResid, unsigned int minHits, double maxReducedChi2)

virtual	~SiStripElectronAlgo ()

Private Member Functions
void	coarseBarrelMonoHitSelection (std::vector< const SiStripRecHit2D * > &monoHitPointersOut)

void	coarseEndcapMonoHitSelection (std::vector< const SiStripRecHit2D * > &monoHitPointersOut)

void	coarseHitSelection (std::vector< const SiStripRecHit2D * > &hitPointersOut, const TrackerTopology *tTopo, bool stereo, bool endcap)

void	coarseMatchedHitSelection (std::vector< const SiStripMatchedRecHit2D * > &coarseMatchedHitPointersOut)

const SiStripElectronAlgo &	operator= (const SiStripElectronAlgo &)

bool	projectPhiBand (float chargeHypothesis, const reco::SuperClusterRef &superclusterIn, const TrackerTopology *tTopo)

	SiStripElectronAlgo (const SiStripElectronAlgo &)

double	unwrapPhi (double phi) const

Private Attributes
double	chi2_neg_

double	chi2_pos_

double	correct_pT_neg_

double	correct_pT_pos_

std::map< const TrackingRecHit *, bool >	hitUsed_

const SiStripRecHit2D *	innerhit_neg_

const SiStripRecHit2D *	innerhit_pos_

double	intercept_neg_

double	intercept_pos_

const MagneticField *	magneticField_p_

const edm::Handle < SiStripMatchedRecHit2DCollection > *	matchedHits_hp_

const SiStripMatchedRecHit2DCollection *	matchedHits_p_

std::map< const TrackingRecHit *, bool >	matchedHitUsed_

std::map< const SiStripMatchedRecHit2D *, unsigned int >	matchedKey_

unsigned int	maxHitsOnDetId_

double	maxNormResid_

double	maxReducedChi2_

unsigned int	minHits_

GlobalVector	momentum_neg_

GlobalVector	momentum_pos_

int	ndof_neg_

int	ndof_pos_

unsigned	numberOfBarrelRphiHits_neg_

unsigned int	numberOfBarrelRphiHits_pos_

unsigned	numberOfEndcapZphiHits_neg_

unsigned int	numberOfEndcapZphiHits_pos_

unsigned int	numberOfMatchedHits_pos_

unsigned	numberOfStereoHits_neg_

unsigned int	numberOfStereoHits_pos_

double	originUncertainty_

std::vector< const TrackingRecHit * >	outputHits_neg_

std::vector< const TrackingRecHit * >	outputHits_pos_

std::vector< SiStripRecHit2D >	outputMatchedHits_neg_

std::vector< SiStripRecHit2D >	outputRphiHits_neg_

std::vector< SiStripRecHit2D >	outputRphiHits_pos_

std::vector< SiStripRecHit2D >	outputStereoHits_neg_

std::vector< SiStripRecHit2D >	outputStereoHits_pos_

double	phiBandWidth_

double	phiVsRSlope_neg_

double	phiVsRSlope_pos_

GlobalPoint	position_neg_

GlobalPoint	position_pos_

double	pZ_neg_

double	pZ_pos_

double	redchi2_neg_

double	redchi2_pos_

const edm::Handle < SiStripRecHit2DCollection > *	rphiHits_hp_

const SiStripRecHit2DCollection *	rphiHits_p_

std::map< const SiStripRecHit2D *, unsigned int >	rphiKey_

double	slope_neg_

double	slope_pos_

const edm::Handle < SiStripRecHit2DCollection > *	stereoHits_hp_

const SiStripRecHit2DCollection *	stereoHits_p_

std::map< const SiStripRecHit2D *, unsigned int >	stereoKey_

const TrackerGeometry *	tracker_p_

double	zVsRSlope_neg_

double	zVsRSlope_pos_

Detailed Description

Description: <one line="" class="" summary>="">

Usage: <usage>

Definition at line 58 of file SiStripElectronAlgo.h.

Constructor & Destructor Documentation

SiStripElectronAlgo::SiStripElectronAlgo	(	unsigned int	maxHitsOnDetId,
		double	originUncertainty,
		double	phiBandWidth,
		double	maxNormResid,
		unsigned int	minHits,
		double	maxReducedChi2
	)

Definition at line 52 of file SiStripElectronAlgo.cc.

   : maxHitsOnDetId_(maxHitsOnDetId)
   , originUncertainty_(originUncertainty)
   , phiBandWidth_(phiBandWidth)
   , maxNormResid_(maxNormResid)
   , minHits_(minHits)
   , maxReducedChi2_(maxReducedChi2)
 {
 }

SiStripElectronAlgo::~SiStripElectronAlgo ( )

virtual

Definition at line 72 of file SiStripElectronAlgo.cc.

73 {

74 }

SiStripElectronAlgo::SiStripElectronAlgo ( const SiStripElectronAlgo & )

private

Member Function Documentation

void SiStripElectronAlgo::coarseBarrelMonoHitSelection ( std::vector< const SiStripRecHit2D * > & monoHitPointersOut )

private

void SiStripElectronAlgo::coarseEndcapMonoHitSelection ( std::vector< const SiStripRecHit2D * > & monoHitPointersOut )

private

void SiStripElectronAlgo::coarseHitSelection	(	std::vector< const SiStripRecHit2D * > &	hitPointersOut,
		const TrackerTopology *	tTopo,
		bool	stereo,
		bool	endcap
	)

private

Definition at line 354 of file SiStripElectronAlgo.cc.

Referenced by projectPhiBand().

 {
   // This function is not time-efficient.  If you want to improve the
   // speed of the algorithm, you'll probably want to change this
   // function.  There may be a more efficienct way to extract hits,
   // and it would definitely help to put a geographical cut on the
   // DetIds.  (How does one determine the global position of a given
   // DetId?  Is tracker_p_->idToDet(id)->surface().toGlobal(LocalPosition(0,0,0))
   // expensive?)
 
   // Loop over the detector ids
   SiStripRecHit2DCollection::const_iterator itdet = (stereo ? stereoHits_p_->begin() : rphiHits_p_->begin());
   SiStripRecHit2DCollection::const_iterator eddet = (stereo ? stereoHits_p_->end()   : rphiHits_p_->end()  );
   for (; itdet != eddet; ++itdet) {
     // Get the hits on this detector id
     SiStripRecHit2DCollection::DetSet hits = *itdet;
     DetId id(hits.detId());
 
     // Count the number of hits on this detector id
     unsigned int numberOfHits = hits.size();
       
     // Only take the hits if there aren't too many
     // (Would it be better to loop only once, fill a temporary list,
     // and copy that if numberOfHits <= maxHitsOnDetId_?)
     if (numberOfHits <= maxHitsOnDetId_) {
       for (SiStripRecHit2DCollection::DetSet::const_iterator hit = hits.begin();  
            hit != hits.end();  ++hit) {
         // check that hit is valid first !
     if(!(*hit).isValid()) {
       LogDebug("") << " InValid hit skipped in coarseHitSelection " << std::endl ;
       continue ;
     }
         std::string theDet = "null";
         int theLayer = -999;
         bool isStereoDet = false ;
         if(tracker_p_->idToDetUnit(hit->geographicalId())->type().subDetector() == GeomDetEnumerators::TIB) { 
           theDet = "TIB" ;
           theLayer = tTopo->tibLayer(id); 
           if(tTopo->tibStereo(id)==1) { isStereoDet = true ; }
         } else if
           (tracker_p_->idToDetUnit(hit->geographicalId())->type().subDetector() == GeomDetEnumerators::TOB) { 
           theDet = "TOB" ;
           theLayer = tTopo->tobLayer(id); 
           if(tTopo->tobStereo(id)==1) { isStereoDet = true ; }
         }else if
           (tracker_p_->idToDetUnit(hit->geographicalId())->type().subDetector() == GeomDetEnumerators::TID) { 
           theDet = "TID" ;
           theLayer = tTopo->tidWheel(id);  // or ring  ?
           if(tTopo->tidStereo(id)==1) { isStereoDet = true ; }
         }else if
           (tracker_p_->idToDetUnit(hit->geographicalId())->type().subDetector() == GeomDetEnumerators::TEC) { 
           theDet = "TEC" ;
           theLayer = tTopo->tecWheel(id);  // or ring or petal ?
           if(tTopo->tecStereo(id)==1) { isStereoDet = true ; }
         } else {
           LogDebug("") << " UHOH BIG PROBLEM - Unrecognized SI Layer" ;
           LogDebug("") << " Det "<< theDet << " Lay " << theLayer ;
           assert(1!=1) ;
         }
 
     if ((endcap  &&  stereo && (theDet=="TID" || theDet== "TEC") && isStereoDet ) ||
             (endcap  &&  !stereo && (theDet=="TID" || theDet== "TEC") && !isStereoDet )  ||
             (!endcap  && stereo && (theDet=="TIB" || theDet=="TOB") &&  isStereoDet )    ||
             (!endcap  &&  !stereo && (theDet=="TIB" || theDet=="TOB" )&& !isStereoDet )
             ) {  
               
 
       hitPointersOut.push_back(&(*hit));
 
     } // end if this is the right subdetector
       } // end loop over hits
     } // end if this detector id doesn't have too many hits on it
   } // end loop over detector ids
 }

void SiStripElectronAlgo::coarseMatchedHitSelection ( std::vector< const SiStripMatchedRecHit2D * > & coarseMatchedHitPointersOut )

private

Definition at line 433 of file SiStripElectronAlgo.cc.

References edmNew::DetSet< T >::begin(), edmNew::DetSetVector< T >::begin(), edmNew::DetSet< T >::end(), edmNew::DetSetVector< T >::end(), LogDebug, matchedHits_p_, maxHitsOnDetId_, StripSubdetector::TIB, and StripSubdetector::TOB.

Referenced by projectPhiBand().

 {
   
   // Loop over the detector ids
   SiStripMatchedRecHit2DCollection::const_iterator itdet = matchedHits_p_->begin(), eddet = matchedHits_p_->end();
   for (; itdet != eddet; ++itdet) {
     
     // Get the hits on this detector id
     SiStripMatchedRecHit2DCollection::DetSet hits = *itdet ;
     
     // Count the number of hits on this detector id
     unsigned int numberOfHits = 0;
     for (SiStripMatchedRecHit2DCollection::DetSet::const_iterator hit = hits.begin();  hit != hits.end();  ++hit) {
       if ( !((hit->geographicalId()).subdetId() == StripSubdetector::TIB) &&
            !( (hit->geographicalId()).subdetId() == StripSubdetector::TOB )) { break;}
       numberOfHits++;
       if (numberOfHits > maxHitsOnDetId_) { break; }
     }
     
     // Only take the hits if there aren't too many
     if (numberOfHits <= maxHitsOnDetId_) {
       for (SiStripMatchedRecHit2DCollection::DetSet::const_iterator hit = hits.begin();  hit != hits.end();  ++hit) {
     if(!(*hit).isValid()) {
       LogDebug("") << " InValid hit skipped in coarseMatchedHitSelection " << std::endl ;
       continue ;
     }
         if ( !((hit->geographicalId()).subdetId() == StripSubdetector::TIB) &&
              !( (hit->geographicalId()).subdetId() == StripSubdetector::TOB )) { break;}
         
         coarseMatchedHitPointersOut.push_back(&(*hit));
       } // end loop over hits
       
     } // end if this detector id doesn't have too many hits on it
   } // end loop over detector ids
   
 
 }// end of matchedHitSelection

bool SiStripElectronAlgo::findElectron	(	reco::SiStripElectronCollection &	electronOut,
		TrackCandidateCollection &	trackCandidateOut,
		const reco::SuperClusterRef &	superclusterIn,
		const TrackerTopology *	tTopo
	)

Definition at line 144 of file SiStripElectronAlgo.cc.

References alongMomentum, assert(), chi2_neg_, chi2_pos_, correct_pT_neg_, correct_pT_pos_, benchmark_cfg::errors, TrackingRecHit::geographicalId(), hitUsed_, TrackerGeometry::idToDet(), innerhit_neg_, innerhit_pos_, intercept_neg_, intercept_pos_, LogDebug, magneticField_p_, momentum_neg_, momentum_pos_, ndof_neg_, ndof_pos_, numberOfBarrelRphiHits_neg_, numberOfBarrelRphiHits_pos_, numberOfEndcapZphiHits_neg_, numberOfEndcapZphiHits_pos_, numberOfStereoHits_neg_, numberOfStereoHits_pos_, outputHits_neg_, outputHits_pos_, outputRphiHits_neg_, outputRphiHits_pos_, outputStereoHits_neg_, outputStereoHits_pos_, perp(), trajectoryStateTransform::persistentState(), phi, phiVsRSlope_neg_, phiVsRSlope_pos_, position_neg_, position_pos_, projectPhiBand(), edm::OwnVector< T, P >::push_back(), pZ_neg_, pZ_pos_, DetId::rawId(), redchi2_neg_, redchi2_pos_, slope_neg_, slope_pos_, python.multivaluedict::sort(), tracker_p_, z, zVsRSlope_neg_, and zVsRSlope_pos_.

Referenced by SiStripElectronProducer::produce().

 {
   // Try each of the two charge hypotheses, but only take one
   bool electronSuccess = projectPhiBand(-1., superclusterIn, tTopo);
   bool positronSuccess = projectPhiBand( 1., superclusterIn, tTopo);
 
   // electron hypothesis did better than electron
   if ((electronSuccess  &&  !positronSuccess)  ||
       (electronSuccess  &&  positronSuccess  &&  redchi2_neg_ <= redchi2_pos_)) {
       
     // Initial uncertainty for tracking
     AlgebraicSymMatrix55 errors;  // makes 5x5 matrix, indexed from (0,0) to (44)
     errors(0,0) = 3.;      // uncertainty**2 in 1/momentum
     errors(1,1) = 0.01;    // uncertainty**2 in lambda (lambda == pi/2 - polar angle theta)
     errors(2,2) = 0.0001;  // uncertainty**2 in phi
     errors(3,3) = 0.01;    // uncertainty**2 in x_transverse (where x is in cm)
     errors(4,4) = 0.01;    // uncertainty**2 in y_transverse (where y is in cm)
 
 #ifdef EDM_ML_DEBUG 
     // JED Debugging possible double hit sort problem
     std::ostringstream debugstr6;
     debugstr6 << " HERE BEFORE SORT electron case " << " \n" ;
     for (std::vector<const TrackingRecHit*>::iterator itHit=outputHits_neg_.begin(); 
      itHit != outputHits_neg_.end(); ++itHit) {
       debugstr6 <<" HIT "<<((*itHit)->geographicalId()).rawId()<<" \n"
         <<"    Local Position: "<<(*itHit)->localPosition()<<" \n"
         <<"    Global Rho:  "
         <<(((TrackingGeometry*)(tracker_p_))->idToDet((*itHit)->geographicalId())->surface().toGlobal((*itHit)->localPosition()).perp())
         <<" Phi "
         <<(((TrackingGeometry*)(tracker_p_))->idToDet((*itHit)->geographicalId())->surface().toGlobal((*itHit)->localPosition()).phi())
         << " Z "
         <<(((TrackingGeometry*)(tracker_p_))->idToDet((*itHit)->geographicalId())->surface().toGlobal((*itHit)->localPosition()).z())<< " \n";
     }
     // end of dump 
     
     
     
     // JED call dump 
     debugstr6 << " Calling sort alongMomentum " << " \n"; 
 #endif   
  
     std::sort(outputHits_neg_.begin(), outputHits_neg_.end(),
               CompareHits(TrackingRecHitLessFromGlobalPosition(((TrackingGeometry*)(tracker_p_)), alongMomentum)));
     
 #ifdef EDM_ML_DEBUG 
     debugstr6 << " Done with sort " << " \n";
     
     debugstr6 << " HERE AFTER SORT electron case " << " \n";
     for (std::vector<const TrackingRecHit*>::iterator itHit=outputHits_neg_.begin(); 
      itHit != outputHits_neg_.end(); ++itHit) {
       debugstr6 <<" HIT "<<((*itHit)->geographicalId()).rawId()<<" \n"
         <<"    Local Position: "<<(*itHit)->localPosition()<<" \n"
         <<"    Global Rho:  "
         <<(((TrackingGeometry*)(tracker_p_))->idToDet((*itHit)->geographicalId())->surface().toGlobal((*itHit)->localPosition()).perp())
         <<" Phi "
         <<(((TrackingGeometry*)(tracker_p_))->idToDet((*itHit)->geographicalId())->surface().toGlobal((*itHit)->localPosition()).phi())
         << " Z "
         <<(((TrackingGeometry*)(tracker_p_))->idToDet((*itHit)->geographicalId())->surface().toGlobal((*itHit)->localPosition()).z())<< " \n";; 
     }
     // end of dump 
 
     LogDebug("")<< debugstr6.str();
 #endif
     
     //create an OwnVector needed by the classes which will be stored to the Event
     edm::OwnVector<TrackingRecHit> hits;
     for(std::vector<const TrackingRecHit*>::iterator itHit =outputHits_neg_.begin();
     itHit != outputHits_neg_.end();
     ++itHit) {
       hits.push_back( (*itHit)->clone());
       if( !(hitUsed_.find(*itHit) != hitUsed_.end()) ) {
     LogDebug("") << " Assert failure " ;
     assert(hitUsed_.find(*itHit) != hitUsed_.end());
       }
       hitUsed_[*itHit] = true;
     }
     
     TrajectoryStateOnSurface state(
                    GlobalTrajectoryParameters(position_neg_, momentum_neg_, -1, magneticField_p_),
                    CurvilinearTrajectoryError(errors),
                    tracker_p_->idToDet(innerhit_neg_->geographicalId())->surface());
     
     
     PTrajectoryStateOnDet const & PTraj = trajectoryStateTransform::persistentState(state, innerhit_neg_->geographicalId().rawId());
     TrajectorySeed trajectorySeed(PTraj, hits, alongMomentum);
     trackCandidateOut.push_back(TrackCandidate(hits, trajectorySeed, PTraj));
     
     electronOut.push_back(reco::SiStripElectron(superclusterIn,
                         -1,
                         outputRphiHits_neg_,
                         outputStereoHits_neg_,
                         phiVsRSlope_neg_,
                         slope_neg_,
                         intercept_neg_ + superclusterIn->position().phi(),
                         chi2_neg_,
                         ndof_neg_,
                         correct_pT_neg_,
                         pZ_neg_,
                         zVsRSlope_neg_,
                         numberOfStereoHits_neg_,
                         numberOfBarrelRphiHits_neg_,
                         numberOfEndcapZphiHits_neg_));
     
     return true;
   }
   
   // positron hypothesis did better than electron
   if ((!electronSuccess  &&  positronSuccess)  ||
       (electronSuccess  &&  positronSuccess  &&  redchi2_neg_ > redchi2_pos_)) {
       
     // Initial uncertainty for tracking
     AlgebraicSymMatrix55 errors;  // same as abovr
     errors(0,0) = 3.;      // uncertainty**2 in 1/momentum
     errors(1,1) = 0.01;    // uncertainty**2 in lambda (lambda == pi/2 - polar angle theta)
     errors(2,2) = 0.0001;  // uncertainty**2 in phi
     errors(3,3) = 0.01;    // uncertainty**2 in x_transverse (where x is in cm)
     errors(4,4) = 0.01;    // uncertainty**2 in y_transverse (where y is in cm)
 
 #ifdef EDM_ML_DEBUG 
     // JED Debugging possible double hit sort problem
     std::ostringstream debugstr7;
     debugstr7 << " HERE BEFORE SORT Positron case " << " \n";
     for (std::vector<const TrackingRecHit*>::iterator itHit = outputHits_pos_.begin(); 
      itHit != outputHits_pos_.end(); ++itHit) {
       debugstr7 <<" HIT "<<((*itHit)->geographicalId()).rawId()<<" \n"
         <<"    Local Position: "<<(*itHit)->localPosition()<<" \n"
         <<"    Global Rho:  "
         <<(((TrackingGeometry*)(tracker_p_))->idToDet((*itHit)->geographicalId())->surface().toGlobal((*itHit)->localPosition()).perp())
         <<" Phi "
         <<(((TrackingGeometry*)(tracker_p_))->idToDet((*itHit)->geographicalId())->surface().toGlobal((*itHit)->localPosition()).phi())
         << " Z "
         <<(((TrackingGeometry*)(tracker_p_))->idToDet((*itHit)->geographicalId())->surface().toGlobal((*itHit)->localPosition()).z())<< " \n";
     }
     // end of dump 
     
     debugstr7 << " Calling sort alongMomentum " << " \n"; 
 #endif
     
     std::sort(outputHits_pos_.begin(), outputHits_pos_.end(),
           CompareHits(TrackingRecHitLessFromGlobalPosition(((TrackingGeometry*)(tracker_p_)), alongMomentum)));
     
 #ifdef EDM_ML_DEBUG 
     debugstr7 << " Done with sort " << " \n";
     
     // JED Debugging possible double hit sort problem
     debugstr7 << " HERE AFTER SORT Positron case " << " \n";
     for (std::vector<const TrackingRecHit*>::iterator itHit = outputHits_pos_.begin(); 
      itHit != outputHits_pos_.end(); ++itHit) {
       debugstr7 <<" HIT "<<((*itHit)->geographicalId()).rawId()<<" \n"
         <<"    Local Position: "<<(*itHit)->localPosition()<<" \n"
         <<"    Global Rho:  "
         <<(((TrackingGeometry*)(tracker_p_))->idToDet((*itHit)->geographicalId())->surface().toGlobal((*itHit)->localPosition()).perp())
         <<" Phi "
         <<(((TrackingGeometry*)(tracker_p_))->idToDet((*itHit)->geographicalId())->surface().toGlobal((*itHit)->localPosition()).phi())
         << " Z "
         <<(((TrackingGeometry*)(tracker_p_))->idToDet((*itHit)->geographicalId())->surface().toGlobal((*itHit)->localPosition()).z())<< " \n"; 
     }
     // end of dump 
     LogDebug("") << debugstr7.str();
 #endif
 
     //create an OwnVector needed by the classes which will be stored to the Event
     edm::OwnVector<TrackingRecHit> hits;
     for(std::vector<const TrackingRecHit*>::iterator itHit =outputHits_pos_.begin();
     itHit != outputHits_pos_.end();
     ++itHit) {
       hits.push_back( (*itHit)->clone());
       assert(hitUsed_.find(*itHit) != hitUsed_.end());
       hitUsed_[*itHit] = true;
     }
     
     TrajectoryStateOnSurface state(
                    GlobalTrajectoryParameters(position_pos_, momentum_pos_, 1, magneticField_p_),
                    CurvilinearTrajectoryError(errors),
                    tracker_p_->idToDet(innerhit_pos_->geographicalId())->surface());
     
     
     PTrajectoryStateOnDet const & PTraj = trajectoryStateTransform::persistentState(state, innerhit_pos_->geographicalId().rawId());
     TrajectorySeed trajectorySeed(PTraj, hits, alongMomentum);
     trackCandidateOut.push_back(TrackCandidate(hits, trajectorySeed, PTraj));
     
     electronOut.push_back(reco::SiStripElectron(superclusterIn,
                         1,
                         outputRphiHits_pos_,
                         outputStereoHits_pos_,
                         phiVsRSlope_pos_,
                         slope_pos_,
                         intercept_pos_ + superclusterIn->position().phi(),
                         chi2_pos_,
                         ndof_pos_,
                         correct_pT_pos_,
                         pZ_pos_,
                         zVsRSlope_pos_,
                         numberOfStereoHits_pos_,
                         numberOfBarrelRphiHits_pos_,
                         numberOfEndcapZphiHits_pos_));
     
     return true;
   }
   
   return false;
 }

const SiStripElectronAlgo& SiStripElectronAlgo::operator= ( const SiStripElectronAlgo & )

private

void SiStripElectronAlgo::prepareEvent	(	const edm::ESHandle< TrackerGeometry > &	tracker,
		const edm::Handle< SiStripRecHit2DCollection > &	rphiHits,
		const edm::Handle< SiStripRecHit2DCollection > &	stereoHits,
		const edm::Handle< SiStripMatchedRecHit2DCollection > &	matchedHits,
		const edm::ESHandle< MagneticField > &	magneticField
	)

Definition at line 92 of file SiStripElectronAlgo.cc.

References counter, edmNew::DetSetVector< T >::data(), hitUsed_, LogDebug, magneticField_p_, matchedHits_p_, matchedHitUsed_, matchedKey_, edm::Handle< T >::product(), edm::ESHandle< class >::product(), rphiHits_hp_, rphiHits_p_, rphiKey_, stereoHits_hp_, stereoHits_p_, stereoKey_, and tracker_p_.

Referenced by SiStripElectronProducer::produce().

 {
   LogDebug("") << " In prepareEvent " ; 
 
   tracker_p_ = tracker.product();
   rphiHits_p_ = rphiHits.product();
   stereoHits_p_ = stereoHits.product();
   matchedHits_p_ = matchedHits.product();
   magneticField_p_ = magneticField.product();
 
   rphiHits_hp_ = &rphiHits;
   stereoHits_hp_ = &stereoHits;
 
   // Keep a table that relates hit pointers to their index (key) in the collections
   rphiKey_.clear();
   stereoKey_.clear();
   // Keep track of which hits have been used already (so a hit is assigned to only one electron)
   hitUsed_.clear();
   matchedHitUsed_.clear();
 
   unsigned int counter = 0;
   for (SiStripRecHit2DCollection::DataContainer::const_iterator it = rphiHits_p_->data().begin();  it != rphiHits_p_->data().end();  ++it) {
     rphiKey_[&(*it)] = counter;
     hitUsed_[&(*it)] = false;
     counter++;
   }
 
   counter = 0;
   for (SiStripRecHit2DCollection::DataContainer::const_iterator it = stereoHits_p_->data().begin();  it != stereoHits_p_->data().end();  ++it) {
     stereoKey_[&(*it)] = counter;
     hitUsed_[&(*it)] = false;
     counter++;
   }
 
   counter = 0;
   for (SiStripMatchedRecHit2DCollection::DataContainer::const_iterator it = matchedHits_p_->data().begin();  it != matchedHits_p_->data().end();  ++it) {
     matchedKey_[&(*it)] = counter;
     matchedHitUsed_[&(*it)] = false;
     counter++;
   }
   
   LogDebug("") << " Leaving prepareEvent " ;
 }

bool SiStripElectronAlgo::projectPhiBand	(	float	chargeHypothesis,
		const reco::SuperClusterRef &	superclusterIn,
		const TrackerTopology *	tTopo
	)

private

Definition at line 477 of file SiStripElectronAlgo.cc.

References funct::abs(), assert(), chi2_neg_, chi2_pos_, coarseHitSelection(), coarseMatchedHitSelection(), correct_pT_neg_, correct_pT_pos_, funct::cos(), run_regression::done, relval_parameters_module::energy, edm::false, TrackingRecHit::geographicalId(), hitUsed_, i, TrackerGeometry::idToDet(), TrackerGeometry::idToDetUnit(), innerhit_neg_, innerhit_pos_, intercept_neg_, intercept_pos_, MagneticField::inTesla(), j, BaseTrackerRecHit::localPosition(), LogDebug, magneticField_p_, matchedHitUsed_, maxNormResid_, maxReducedChi2_, minHits_, momentum_neg_, momentum_pos_, ndof_neg_, ndof_pos_, numberOfBarrelRphiHits_neg_, numberOfBarrelRphiHits_pos_, numberOfEndcapZphiHits_neg_, numberOfEndcapZphiHits_pos_, numberOfStereoHits_neg_, numberOfStereoHits_pos_, originUncertainty_, outputHits_neg_, outputHits_pos_, outputRphiHits_neg_, outputRphiHits_pos_, outputStereoHits_neg_, outputStereoHits_pos_, phi, PV3DBase< T, PVType, FrameType >::phi(), phiBandWidth_, phiVsRSlope_neg_, phiVsRSlope_pos_, position, position_neg_, position_pos_, pZ_neg_, pZ_pos_, alignCSCRings::r, redchi2_neg_, redchi2_pos_, funct::sin(), slope, slope_neg_, slope_pos_, python.multivaluedict::sort(), mathSSE::sqrt(), GeomDetEnumerators::TIB, TrackerTopology::tibLayer(), GeomDetEnumerators::TOB, TrackerTopology::tobLayer(), GeomDet::toGlobal(), tracker_p_, funct::true, unwrapPhi(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), z, PV3DBase< T, PVType, FrameType >::z(), zVsRSlope_neg_, and zVsRSlope_pos_.

Referenced by findElectron().

 {
   // This algorithm projects a phi band into the tracker three times:
   // (a) for all stereo hits, (b) for barrel rphi hits, and (c) for
   // endcap zphi hits.  While accumulating stereo hits in step (a),
   // we fit r vs z to a line.  This resolves the ambiguity in z for
   // rphi hits and the ambiguity in r for zphi hits.  We can then cut
   // on the z of rphi hits (a little wider than one strip length),
   // and we can convert the z of zphi hits into r to apply the phi
   // band cut.  (We don't take advantage of the endcap strips'
   // segmentation in r.)
   // 
   // As we project a phi band into the tracker, we count hits within
   // that band and performs a linear fit for phi vs r.  The number of
   // hits and reduced chi^2 from the fit are used to select a good
   // candidate.
 
   // set limits on Si hits - smallest error that makes sense = 1 micron ?
   static const double rphiHitSmallestError = 0.0001 ;
   static const double stereoHitSmallestError = 0.0001 ;
   //not used since endcap is not implemented  
   // static const double zphiHitSmallestError = 0.0001 ;
 
 
   static const double stereoPitchAngle = 0.100 ; // stereo angle in rad
   static const double cos2SiPitchAngle = cos(stereoPitchAngle)*cos(stereoPitchAngle) ;
   static const double sin2SiPitchAngle = sin(stereoPitchAngle)*sin(stereoPitchAngle) ;
   // overall misalignment fudge to be added in quad to position errors.
   // this is a rough approx to values reported in tracking meet 5/16/2007
   static const double rphiErrFudge = 0.0200 ;
   static const double stereoErrFudge = 0.0200 ;
 
   // max chi2 of a hit on an SiDet relative to the prediction
   static const double chi2HitMax = 25.0 ;
 
   // Minimum number of hits to consider on a candidate
   static const unsigned int nHitsLeftMinimum = 3  ;
 
   // Create and fill vectors of pointers to hits
   std::vector<const SiStripRecHit2D*> stereoHits;
   std::vector<const SiStripRecHit2D*> rphiBarrelHits;
   std::vector<const SiStripRecHit2D*> zphiEndcapHits;
 
   //                                 stereo? endcap?
   coarseHitSelection(stereoHits, tTopo,    true,   false);
 
   // skip endcap stereo for now
   //  LogDebug("") << " Getting endcap stereo hits " ;
   // coarseHitSelection(stereoHits,     true,   true);
 
   std::ostringstream debugstr1;
   debugstr1 << " Getting barrel rphi hits " << " \n" ;
 
   coarseHitSelection(rphiBarrelHits, tTopo, false,  false);
 
   //  LogDebug("") << " Getting endcap zphi hits " ;
   //  coarseHitSelection(zphiEndcapHits, false,  true);
 
   debugstr1 << " Getting matched hits "  << " \n" ;
   std::vector<const SiStripMatchedRecHit2D*> matchedHits;
   coarseMatchedHitSelection(matchedHits);
 
 
   // Determine how to project from the supercluster into the tracker
   double energy = superclusterIn->energy();
   double pT = energy * superclusterIn->position().rho()/sqrt(superclusterIn->x()*superclusterIn->x() +
                                  superclusterIn->y()*superclusterIn->y() +
                                  superclusterIn->z()*superclusterIn->z());
   // cf Jackson p. 581-2, a little geometry
   double phiVsRSlope = -3.00e-3 * chargeHypothesis * magneticField_p_->inTesla(GlobalPoint(superclusterIn->x(), superclusterIn->y(), 0.)).z() / pT / 2.;
 
   // Shorthand for supercluster radius, z
   const double scr = superclusterIn->position().rho();
   const double scz = superclusterIn->position().z();
 
   // These are used to fit all hits to a line in phi(r)
   std::vector<bool> uselist;
   std::vector<double> rlist, philist, w2list;
   std::vector<int> typelist;  // stereo = 0, rphi barrel = 1, and zphi disk = 2 (only used in this function)
   std::vector<const SiStripRecHit2D*> hitlist;
 
   std::vector<bool> matcheduselist;
   std::vector<const SiStripMatchedRecHit2D*> matchedhitlist;
 
   // These are used to fit the stereo hits to a line in z(r), constrained to pass through supercluster
   double zSlopeFitNumer = 0.;
   double zSlopeFitDenom = 0.;
 
   
   debugstr1 << " There are a total of " << stereoHits.size()  << " stereoHits in this event " << " \n" 
         << " There are a total of " <<  rphiBarrelHits.size() << " rphiBarrelHits in this event " << " \n"
         << " There are a total of " <<  zphiEndcapHits.size() << " zphiEndcapHits in this event " << " \n\n";
 
 
   LogDebug("") << debugstr1.str() ;
   
 
 
   // Loop over all matched hits
   // make a list of good matched rechits.
   // in the stereo and rphi loops check to see if the hit is associated with a matchedhit
   LogDebug("") << " Loop over matched hits " << " \n";
 
   for (std::vector<const SiStripMatchedRecHit2D*>::const_iterator hit = matchedHits.begin() ;
        hit != matchedHits.end() ; ++ hit) {
     
     assert(matchedHitUsed_.find(*hit) != matchedHitUsed_.end());
 
     if (!matchedHitUsed_[*hit]) {
 
       // Calculate the 3-D position of this hit
       GlobalPoint position = tracker_p_->idToDet((*hit)->geographicalId())->surface().toGlobal((*hit)->localPosition());
       double r = sqrt(position.x()*position.x() + position.y()*position.y());
       double phi = unwrapPhi(position.phi() - superclusterIn->position().phi());  // phi is relative to supercluster
       double z = position.z();
 
       // Cut a triangle in the z-r plane using the supercluster's eta/dip angle information
       // and the fact that the electron originated *near* the origin
       if ((-originUncertainty_ + (scz + originUncertainty_)*(r/scr)) < z  &&  z < (originUncertainty_ + (scz - originUncertainty_)*(r/scr))) {
      
     // Cut a narrow band around the supercluster's projection in phi
     if (unwrapPhi((r-scr)*phiVsRSlope - phiBandWidth_) < phi  &&  phi < unwrapPhi((r-scr)*phiVsRSlope + phiBandWidth_)) {
      
 
       matcheduselist.push_back(true);
           matchedhitlist.push_back(*hit);
       
 
 
       // Use this hit to fit z(r)
       zSlopeFitNumer += -(scr - r) * (z - scz);
       zSlopeFitDenom += (scr - r) * (scr - r);
         
     } // end cut on phi band
       } // end cut on electron originating *near* the origin
     } // end assign disjoint sets of hits to electrons
   } // end loop over matched hits
 
   // Calculate the linear fit for z(r)
   double zVsRSlope;
   if (zSlopeFitDenom > 0.) {
     zVsRSlope = zSlopeFitNumer / zSlopeFitDenom;
   }
   else {
     // zVsRSlope assumes electron is from origin if there were no stereo hits
     zVsRSlope = scz/scr;
   }
 
   //  // Loop over all stereo hits
   LogDebug("") << " Loop over stereo hits" << " \n";
 
   // check if the stereo hit is matched to one of the matched hit
   unsigned int numberOfStereoHits = 0;
   for (std::vector<const SiStripRecHit2D*>::const_iterator hit = stereoHits.begin();  hit != stereoHits.end();  ++hit) {
     assert(hitUsed_.find(*hit) != hitUsed_.end());
 
     // Calculate the 3-D position of this hit
     GlobalPoint position = tracker_p_->idToDet((*hit)->geographicalId())->surface().toGlobal((*hit)->localPosition());
     double r_stereo = sqrt(position.x()*position.x() + position.y()*position.y());
     double phi_stereo = unwrapPhi(position.phi() - superclusterIn->position().phi());  // phi is relative to supercluster
     //    double z_stereo = position.z();
 
     // stereo is has a pitch of 100 mrad - so consider both components
     double r_stereo_err = sqrt((*hit)->localPositionError().xx()*cos2SiPitchAngle +
                                (*hit)->localPositionError().yy()*sin2SiPitchAngle ) ; 
 
 
 
     // make sure that the error for this hit is sensible, ie > 1 micron
     // otherwise skip this hit
     if(r_stereo_err > stereoHitSmallestError ) {
       r_stereo_err = sqrt(r_stereo_err*r_stereo_err+stereoErrFudge*stereoErrFudge);
  
       OmniClusterRef const & stereocluster=(*hit)->omniClusterRef();
     
       bool thisHitIsMatched = false ;
 
       if (!hitUsed_[*hit]) {
  
     unsigned int matcheduselist_size = matcheduselist.size();
     for (unsigned int i = 0;  i < matcheduselist_size;  i++) {
       if (matcheduselist[i]) {
             OmniClusterRef const &  mystereocluster = matchedhitlist[i]->stereoClusterRef();
         if( stereocluster == mystereocluster ) {
           thisHitIsMatched = true ;
           //    LogDebug("")<< "     This hit is matched " << tracker_p_->idToDet(matchedhitlist[i]->stereoHit()->geographicalId())->surface().toGlobal(matchedhitlist[i]->stereoHit()->localPosition()) << std::endl;
           //      break ;
         }
       } // check if matcheduselist okay 
     }// loop over matched hits 
       
     if(thisHitIsMatched) {
       // Use this hit to fit phi(r)
       uselist.push_back(true);
       rlist.push_back(r_stereo);
       philist.push_back(phi_stereo);
       w2list.push_back(1./(r_stereo_err/r_stereo)/(r_stereo_err/r_stereo));  // weight**2 == 1./uncertainty**2
       typelist.push_back(0);
       hitlist.push_back(*hit);
     } // thisHitIsMatched
       } //  if(!hitUsed)
 
     } //  end of check on hit position error size 
 
   } // end loop over stereo hits
   
   LogDebug("") << " There are " << uselist.size()  << " good hits after stereo loop "  ;
  
 
   // Loop over barrel rphi hits
   LogDebug("") << " Looping over barrel rphi hits " ;
   unsigned int rphiMatchedCounter = 0 ;
   unsigned int rphiUnMatchedCounter = 0 ;
   unsigned int numberOfBarrelRphiHits = 0;
   for (std::vector<const SiStripRecHit2D*>::const_iterator hit = rphiBarrelHits.begin();  hit != rphiBarrelHits.end();  ++hit) {
     assert(hitUsed_.find(*hit) != hitUsed_.end());
     // Calculate the 2.5-D position of this hit
     GlobalPoint position = tracker_p_->idToDet((*hit)->geographicalId())->surface().toGlobal((*hit)->localPosition());
     double r = sqrt(position.x()*position.x() + position.y()*position.y());
     double phi = unwrapPhi(position.phi() - superclusterIn->position().phi());  // phi is relative to supercluster
     double z = position.z();
     double r_mono_err = sqrt((*hit)->localPositionError().xx()) ;
 
     // only consider hits with errors that make sense
     if( r_mono_err > rphiHitSmallestError) {
       // inflate the reported error
       r_mono_err=sqrt(r_mono_err*r_mono_err+rphiErrFudge*rphiErrFudge);
 
       OmniClusterRef const & monocluster=(*hit)->omniClusterRef();
     
 
       if (!hitUsed_[*hit]) {
     if( (tracker_p_->idToDetUnit((*hit)->geographicalId())->type().subDetector() == GeomDetEnumerators::TIB  && 
          (tTopo->tibLayer((*hit)->geographicalId())==1 || tTopo->tibLayer((*hit)->geographicalId())==2)) || 
         (tracker_p_->idToDetUnit((*hit)->geographicalId())->type().subDetector() == GeomDetEnumerators::TOB  
          && (tTopo->tobLayer((*hit)->geographicalId())==1 || tTopo->tobLayer((*hit)->geographicalId())==2)) ) {
       bool thisHitIsMatched = false ;
       unsigned int matcheduselist_size = matcheduselist.size();
       for (unsigned int i = 0;  i < matcheduselist_size;  i++) {
         if (matcheduselist[i]) {
               OmniClusterRef const &  mymonocluster = matchedhitlist[i]->monoClusterRef();
           if( monocluster == mymonocluster ) {
         thisHitIsMatched = true ;
           } 
         } // check if matcheduselist okay 
       }// loop over matched hits 
         
         
       if( thisHitIsMatched ) {
         // Use this hit to fit phi(r)
         uselist.push_back(true);
         rlist.push_back(r);
         philist.push_back(phi);
         w2list.push_back(1./(r_mono_err/r)/(r_mono_err/r));  // weight**2 == 1./uncertainty**2
         typelist.push_back(1);
         hitlist.push_back(*hit);
         rphiMatchedCounter++;
       } // end of matched hit check
 
     } else {
 
 
       // The expected z position of this hit, according to the z(r) fit
       double zFit = zVsRSlope * (r - scr) + scz;
         
       // Cut on the Z of the strip
       // TIB strips are 11 cm long, TOB strips are 19 cm long (can I get these from a function?)
       if ((tracker_p_->idToDetUnit((*hit)->geographicalId())->type().subDetector() == GeomDetEnumerators::TIB  && 
            std::abs(z - zFit) < 12.)  ||
           (tracker_p_->idToDetUnit((*hit)->geographicalId())->type().subDetector() == GeomDetEnumerators::TOB  && 
            std::abs(z - zFit) < 20.)    ) {
           
         // Cut a narrow band around the supercluster's projection in phi
         if (unwrapPhi((r-scr)*phiVsRSlope - phiBandWidth_) < phi  &&  phi < unwrapPhi((r-scr)*phiVsRSlope + phiBandWidth_)) {
             
           // Use this hit to fit phi(r)
           uselist.push_back(true);
           rlist.push_back(r);
           philist.push_back(phi);
           w2list.push_back(1./(r_mono_err/r)/(r_mono_err/r));  // weight**2 == 1./uncertainty**2
           typelist.push_back(1);
           hitlist.push_back(*hit);
           rphiUnMatchedCounter++;
             
         } // end cut on phi band
       } // end cut on strip z
     } // loop over TIB/TOB layer 1,2 
       } // end assign disjoint sets of hits to electrons
     } // end of check on rphi hit position error size
   } // end loop over barrel rphi hits
 
   LogDebug("") << " There are " << rphiMatchedCounter <<" matched rphi hits"; 
   LogDebug("") << " There are " << rphiUnMatchedCounter <<" unmatched rphi hits";
   LogDebug("") << " There are " << uselist.size() << " good stereo+rphi hits " ;
 
 
 
 
 
 
 
   // Loop over endcap zphi hits
   LogDebug("") << " Looping over barrel zphi hits " ;
 
 
   unsigned int numberOfEndcapZphiHits = 0;
   for (std::vector<const SiStripRecHit2D*>::const_iterator hit = zphiEndcapHits.begin();  
        hit != zphiEndcapHits.end();  ++hit) {
     assert(hitUsed_.find(*hit) != hitUsed_.end());
     if (!hitUsed_[*hit]) {
       // Calculate the 2.5-D position of this hit
       GlobalPoint position = tracker_p_->idToDet((*hit)->geographicalId())->surface().toGlobal((*hit)->localPosition());
       double z = position.z();
       double phi = unwrapPhi(position.phi() - superclusterIn->position().phi());  // phi is relative to supercluster
       //      double r=sqrt(position.x()*position.x()+position.y()*position.y()) ;
 
       // The expected r position of this hit, according to the z(r) fit
       double rFit = (z - scz)/zVsRSlope + scr;
 
       // I don't know the r widths of the endcap strips, otherwise I
       // would apply a cut on r similar to the rphi z cut
 
       // Cut a narrow band around the supercluster's projection in phi,
       // and be sure the hit isn't in a reflected band (extrapolation of large z differences into negative r)
       if (rFit > 0.  &&
       unwrapPhi((rFit-scr)*phiVsRSlope - phiBandWidth_) < phi  &&  phi < unwrapPhi((rFit-scr)*phiVsRSlope + phiBandWidth_)) {
 
     // Use this hit to fit phi(r)
     uselist.push_back(true);
     rlist.push_back(rFit);
     philist.push_back(phi);
     w2list.push_back(1./(0.05/rFit)/(0.05/rFit));  // weight**2 == 1./uncertainty**2
     typelist.push_back(2);
     hitlist.push_back(*hit);
 
       } // end cut on phi band
     } // end assign disjoint sets of hits to electrons
   } // end loop over endcap zphi hits
  
   LogDebug("") << " There are " << uselist.size() << " good stereo+rphi+zphi hits " ;
 
 #ifdef EDM_ML_DEBUG 
   std::ostringstream debugstr5;
   debugstr5 << " List of hits before uniqification " << " \n" ;
   for (unsigned int i = 0;  i < uselist.size();  i++) {
     if ( uselist[i] ) {
       const SiStripRecHit2D* hit = hitlist[i];
       debugstr5 << " DetID " << ((hit)->geographicalId()).rawId()
         << " R " << rlist[i] 
         << " Phi " << philist[i]
         << " Weight " << w2list[i] 
         << " PhiPred " << (rlist[i]-scr)*phiVsRSlope  
         << " Chi2 " << (philist[i]-(rlist[i]-scr)*phiVsRSlope)*(philist[i]-(rlist[i]-scr)*phiVsRSlope)*w2list[i]
         << " \n" ;
     }
   }
   debugstr5 << " \n\n\n" ;
   
   debugstr5 << " Counting number of unique detectors " << " \n" ;
 
   debugstr5 << " These are all the detectors with hits " << " \n";
 #endif
   // Loop over hits, and find the best hit for each SiDetUnit - keep only those
   // get a list of DetIds in uselist
   std::vector<uint32_t> detIdList ;
 
   for (unsigned int i = 0;  i < uselist.size();  i++) {
     if (uselist[i]) {
       const SiStripRecHit2D* hit = hitlist[i];
       uint32_t detIDUsed = ((hit)->geographicalId()).rawId() ;
 #ifdef EDM_ML_DEBUG 
       debugstr5<< " DetId " << detIDUsed << "\n";
 #endif
       detIdList.push_back(detIDUsed) ;
     }
   }
 
 #ifdef EDM_ML_DEBUG 
   debugstr5 << " There are " << detIdList.size() << " hits on detectors \n";
 #endif
   // now sort and then uniq this list of detId
   std::sort(detIdList.begin(), detIdList.end());
   detIdList.erase(
           std::unique(detIdList.begin(), detIdList.end()), detIdList.end());
 #ifdef EDM_ML_DEBUG 
   debugstr5 << " There are " << detIdList.size() << " unique detectors \n";
 #endif
   //now we have a list of unique detectors used
 
 
 #ifdef EDM_ML_DEBUG 
   debugstr5 << " Looping over detectors to choose best hit " << " \n";
 #endif
   // Loop over detectors ID and hits to create list of hits on same DetId
   for (unsigned int idet = 0 ; idet < detIdList.size() ; idet++ ) {
     for (unsigned int i = 0;  i+1 < uselist.size();  i++) {
       if (uselist[i]) {
     // Get Chi2 of this hit relative to predicted hit
     const SiStripRecHit2D* hit1 = hitlist[i];
     double r_hit1 = rlist[i];
     double phi_hit1 = philist[i];
     double w2_hit1 = w2list[i] ;
     double phi_pred1 = (r_hit1-scr)*phiVsRSlope ; 
     double chi1 = (phi_hit1-phi_pred1)*(phi_hit1-phi_pred1)*w2_hit1;
     if(detIdList[idet]== ((hit1)->geographicalId()).rawId()) {
       for (unsigned int j = i+1;  j < uselist.size();  j++) {
         if (uselist[j]) {
           const SiStripRecHit2D* hit2 = hitlist[j];
           if(detIdList[idet]== ((hit2)->geographicalId()).rawId()) {
 #ifdef EDM_ML_DEBUG 
         debugstr5 << " Found 2 hits on same Si Detector " 
               << ((hit2)->geographicalId()).rawId() << "\n";
 #endif
         // Get Chi2 of this hit relative to predicted hit
         double r_hit2 = rlist[j];
         double phi_hit2 = philist[j];
         double w2_hit2 = w2list[j] ;
         double phi_pred2 = (r_hit2-scr)*phiVsRSlope ; 
         double chi2 = (phi_hit2-phi_pred2)*(phi_hit2-phi_pred2)*w2_hit2;
 #ifdef EDM_ML_DEBUG 
         debugstr5 << " Chi1 " << chi1 << " Chi2 " << chi2 <<"\n";
 #endif
         if( chi1< chi2 ){
           uselist[j] = 0;
         }else{
           uselist[i] = 0;
         }
 
           } // end of Det check
         } // end of j- usehit check
       } // end of j hit list loop
       
     } // end of detector check
       } // end of uselist check
     } // end of i-hit list loop
 
 
   } // end of DetId loop
 
 
   
   // now let's through out hits with a predicted chi > chi2HitMax 
   for ( unsigned int i = 0;  i < uselist.size();  i++ ) { 
     if ( uselist[i] ) { 
       double localchi2 = (philist[i]-(rlist[i]-scr)*phiVsRSlope)*(philist[i]-(rlist[i]-scr)*phiVsRSlope)*w2list[i] ;
       if(localchi2 > chi2HitMax ) {
 #ifdef EDM_ML_DEBUG 
         const SiStripRecHit2D* hit = hitlist[i];
     debugstr5 << " Throwing out "
           <<" DetID " << ((hit)->geographicalId()).rawId()
           << " R " << rlist[i] 
           << " Phi " << philist[i]
           << " Weight " << w2list[i] 
           << " PhiPred " << (rlist[i]-scr)*phiVsRSlope  
           << " Chi2 " << (philist[i]-(rlist[i]-scr)*phiVsRSlope)*(philist[i]-(rlist[i]-scr)*phiVsRSlope)*w2list[i]
           << " \n" ;
 #endif
     uselist[i]=0 ;
       }
     }
   }
 
 #ifdef EDM_ML_DEBUG 
   debugstr5 << " List of hits after uniqification " << " \n" ;
   for (unsigned int i = 0;  i < uselist.size();  i++) {
     if ( uselist[i] ) {
       const SiStripRecHit2D* hit = hitlist[i];
       debugstr5 << " DetID " << ((hit)->geographicalId()).rawId()
         << " R " << rlist[i] 
         << " Phi " << philist[i]
         << " Weight " << w2list[i] 
         << " PhiPred " << (rlist[i]-scr)*phiVsRSlope  
         << " Chi2 " << (philist[i]-(rlist[i]-scr)*phiVsRSlope)*(philist[i]-(rlist[i]-scr)*phiVsRSlope)*w2list[i]
         << " \n" ;
     }
   }
   debugstr5 << " \n\n\n" ;
 #endif
 
   // need to check that there are more than 2 hits left here!
   unsigned int nHitsLeft =0;
   for (unsigned int i = 0;  i < uselist.size();  i++) {
     if ( uselist[i] ) {
       nHitsLeft++;
     }
   }
   if(nHitsLeft < nHitsLeftMinimum ) {
 #ifdef EDM_ML_DEBUG 
     debugstr5 << " Too few hits "<< nHitsLeft 
           << " left - return false " << " \n";
 #endif
     return false ;
   }
 #ifdef EDM_ML_DEBUG 
   LogDebug("") << debugstr5.str();
 #endif
   
   // Calculate a linear phi(r) fit and drop hits until the biggest contributor to chi^2 is less than maxNormResid_
   bool done = false;
   double intercept = 0 ;
   double slope = 0 ;
   double chi2 = 0;
 
   std::ostringstream debugstr4;
   debugstr4 <<" Calc of phi(r) "<<" \n";
 
   while (!done) {
     // Use an iterative update of <psi>, <r> etc instead in an
     // attempt to minize divisions of  large numbers
     // The linear fit  psi = slope*r + intercept
     //             slope is (<psi*r>-<psi>*<r>) / (<r^2>-<r>^2>)
     //             intercept is <psi> - slope*<r>
     
     double phiBar   = 0.;
     double phiBarOld   = 0.;
     double rBar  = 0.;
     double rBarOld  = 0.;
     double r2Bar = 0.;
     double r2BarOld = 0.;
     double phirBar = 0.;
     double phirBarOld = 0.;
     double totalWeight = 0.;
     double totalWeightOld = 0.; 
     unsigned int uselist_size = uselist.size();
     for (unsigned int i = 0;  i < uselist_size;  i++) {
       if (uselist[i]) {
         double r = rlist[i];
         double phi = philist[i];
         double weight2 = w2list[i];
     
         totalWeight = totalWeightOld + weight2 ;
     
         //weight2 is 1/sigma^2. Typically sigma is 100micron pitch
         // over root(12) = 30 microns so weight2 is about 10^5-10^6
     
         double totalWeightRatio = totalWeightOld/totalWeight ;
         double localWeightRatio = weight2/totalWeight ;
     
         phiBar = phiBarOld*totalWeightRatio + phi*localWeightRatio ; 
         rBar = rBarOld*totalWeightRatio + r*localWeightRatio ; 
         r2Bar= r2BarOld*totalWeightRatio + r*r*localWeightRatio ; 
         phirBar = phirBarOld*totalWeightRatio + phi*r*localWeightRatio ;
     
     totalWeightOld = totalWeight ;
         phiBarOld = phiBar ;
         rBarOld = rBar ;
         r2BarOld = r2Bar ;
         phirBarOld = phirBar ;  
 #ifdef EDM_ML_DEBUG 
     debugstr4 << " totalWeight " << totalWeight 
                   << " totalWeightRatio " << totalWeightRatio
                   << " localWeightRatio "<< localWeightRatio
                   << " phiBar " << phiBar
                   << " rBar " << rBar 
                   << " r2Bar " << r2Bar
                   << " phirBar " << phirBar 
                   << " \n ";
 #endif
 
       } // end of use hit loop
     } // end of hit loop to calculate a linear fit
     slope = (phirBar-phiBar*rBar)/(r2Bar-rBar*rBar);
     intercept = phiBar-slope*rBar ;
 
     debugstr4 << " end of loop slope " << slope 
               << " intercept " << intercept << " \n";
 
     // Calculate chi^2
     chi2 = 0.;
     double biggest_normresid = -1.;
     unsigned int biggest_index = 0;
     for (unsigned int i = 0;  i < uselist_size;  i++) {
       if (uselist[i]) {
     double r = rlist[i];
     double phi = philist[i];
     double weight2 = w2list[i];
 
     double normresid = weight2 * (phi - intercept - slope*r)*(phi - intercept - slope*r);
     chi2 += normresid;
 
     if (normresid > biggest_normresid) {
       biggest_normresid = normresid;
       biggest_index = i;
     }
       }
     } // end loop over hits to calculate chi^2 and find its biggest contributer
 
     if (biggest_normresid > maxNormResid_) {
 #ifdef EDM_ML_DEBUG
       debugstr4 << "Dropping hit from fit due to Chi2 " << " \n" ;
       const SiStripRecHit2D* hit = hitlist[biggest_index];
       debugstr4 << " DetID " << ((hit)->geographicalId()).rawId()
         << " R " << rlist[biggest_index]
         << " Phi " << philist[biggest_index]
         << " Weight " << w2list[biggest_index]
         << " PhiPred " << (rlist[biggest_index]-scr)*phiVsRSlope  
         << " Chi2 " << (philist[biggest_index]-(rlist[biggest_index]-scr)*phiVsRSlope)*(philist[biggest_index]-(rlist[biggest_index]-scr)*phiVsRSlope)*w2list[biggest_index] 
         << " normresid " <<  biggest_normresid
         << " \n\n";
 #endif
       uselist[biggest_index] = false;
     }
     else {
       done = true;
     }
   } // end loop over trial fits
 #ifdef EDM_ML_DEBUG 
   debugstr4 <<" Fit " 
             << " Intercept  " << intercept
             << " slope " << slope 
             << " chi2 " << chi2
             << " \n" ;
 
   LogDebug("") <<   debugstr4.str() ;
 #endif
 
   // Now we have intercept, slope, and chi2; uselist to tell us which hits are used, and hitlist for the hits
 
   // Identify the innermost hit
   const SiStripRecHit2D* innerhit = (SiStripRecHit2D*)(0);
   double innerhitRadius = -1.;  // meaningless until innerhit is defined
 
   // Copy hits into an OwnVector, which we put in the TrackCandidate
   std::vector<const TrackingRecHit*> outputHits;
   // Reference rphi and stereo hits into RefVectors, which we put in the SiStripElectron
   std::vector<SiStripRecHit2D> outputRphiHits;
   std::vector<SiStripRecHit2D> outputStereoHits;
 
   typedef edm::Ref<SiStripRecHit2DCollection,SiStripRecHit2D> SiStripRecHit2DRef;
 
 
   for (unsigned int i = 0;  i < uselist.size();  i++) {
     if (uselist[i]) {
       double r = rlist[i];
       const SiStripRecHit2D* hit = hitlist[i];
 
       // Keep track of the innermost hit
       if (innerhit == (SiStripRecHit2D*)(0)  ||  r < innerhitRadius) {
     innerhit = hit;
     innerhitRadius = r;
       }
      
       if (typelist[i] == 0) {
     numberOfStereoHits++;
 
     // Copy this hit for the TrajectorySeed
     outputHits.push_back(hit);
     outputStereoHits.push_back(*hit);
       }
       else if (typelist[i] == 1) {
     numberOfBarrelRphiHits++;
 
     // Copy this hit for the TrajectorySeed
     outputHits.push_back(hit);
     outputRphiHits.push_back(*hit);
       }
       else if (typelist[i] == 2) {
     numberOfEndcapZphiHits++;
 
     // Copy this hit for the TrajectorySeed
     outputHits.push_back(hit);
     outputRphiHits.push_back(*hit);
       }
     }
   } // end loop over all hits, after having culled the ones with big residuals
 
   unsigned int totalNumberOfHits = numberOfStereoHits + numberOfBarrelRphiHits + numberOfEndcapZphiHits;
   double reducedChi2 = (totalNumberOfHits > 2 ? chi2 / (totalNumberOfHits - 2) : 1e10);
 
   // Select this candidate if it passes minHits_ and maxReducedChi2_ cuts
   if (totalNumberOfHits >= minHits_  &&  reducedChi2 <= maxReducedChi2_) {
     // GlobalTrajectoryParameters evaluated at the position of the innerhit
     GlobalPoint position = tracker_p_->idToDet(innerhit->geographicalId())->surface().toGlobal(innerhit->localPosition());
 
     // Use our phi(r) linear fit to correct pT (pT is inversely proportional to slope)
     // (By applying a correction instead of going back to first
     // principles, any correction to the phiVsRSlope definition
     // above will be automatically propagated here.)
     double correct_pT = pT * phiVsRSlope / slope;
 
     // Our phi(r) linear fit returns phi relative to the supercluster phi for a given radius
     double phi = intercept + slope*sqrt(position.x()*position.x() + position.y()*position.y()) + superclusterIn->position().phi();
 
     // Our z(r) linear fit gives us a better measurement of eta/dip angle
     double pZ = correct_pT * zVsRSlope;
 
     // Now we can construct a trajectory momentum out of linear fits to hits
     GlobalVector momentum = GlobalVector(correct_pT*cos(phi), correct_pT*sin(phi), pZ);
 
     if (chargeHypothesis > 0.) {
       redchi2_pos_ = chi2 / double(totalNumberOfHits - 2);
       position_pos_ = position;
       momentum_pos_ = momentum;
       innerhit_pos_ = innerhit;
       outputHits_pos_ = outputHits;
       outputRphiHits_pos_ = outputRphiHits;
       outputStereoHits_pos_ = outputStereoHits;
       phiVsRSlope_pos_ = phiVsRSlope;
       slope_pos_ = slope;
       intercept_pos_ = intercept;
       chi2_pos_ = chi2;
       ndof_pos_ = totalNumberOfHits - 2;
       correct_pT_pos_ = correct_pT;
       pZ_pos_ = pZ;
       zVsRSlope_pos_ = zVsRSlope;
       numberOfStereoHits_pos_ = numberOfStereoHits;
       numberOfBarrelRphiHits_pos_ = numberOfBarrelRphiHits;
       numberOfEndcapZphiHits_pos_ = numberOfEndcapZphiHits;
     }
     else {
       redchi2_neg_ = chi2 / double(totalNumberOfHits - 2);
       position_neg_ = position;
       momentum_neg_ = momentum;
       innerhit_neg_ = innerhit;
       outputHits_neg_ = outputHits;
       outputRphiHits_neg_ = outputRphiHits;
       outputStereoHits_neg_ = outputStereoHits;
       phiVsRSlope_neg_ = phiVsRSlope;
       slope_neg_ = slope;
       intercept_neg_ = intercept;
       chi2_neg_ = chi2;
       ndof_neg_ = totalNumberOfHits - 2;
       correct_pT_neg_ = correct_pT;
       pZ_neg_ = pZ;
       zVsRSlope_neg_ = zVsRSlope;
       numberOfStereoHits_neg_ = numberOfStereoHits;
       numberOfBarrelRphiHits_neg_ = numberOfBarrelRphiHits;
       numberOfEndcapZphiHits_neg_ = numberOfEndcapZphiHits;
     }
 
     LogDebug("") << " return from projectFindBand with True \n" << std::endl ;
     return true;
   } // end if this is a good electron candidate
 
   // Signal for a failed electron candidate
   LogDebug("") << " return from projectFindBand with False \n" << std::endl ;
   return false;
 }

double SiStripElectronAlgo::unwrapPhi ( double phi ) const

inlineprivate

Definition at line 114 of file SiStripElectronAlgo.h.

References M_PI, and phi.

Referenced by projectPhiBand().

                                          {
      while (phi > M_PI) { phi -= 2.*M_PI; }
      while (phi < -M_PI) { phi += 2.*M_PI; }
      return phi;
       }