---

SiStripElectronAlgo Class Reference

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

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

List of all members.

Public Member Functions
bool	findElectron (reco::SiStripElectronCollection &electronOut, TrackCandidateCollection &trackCandidateOut, const reco::SuperClusterRef &superclusterIn)
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, bool stereo, bool endcap)
void	coarseMatchedHitSelection (std::vector< const SiStripMatchedRecHit2D * > &coarseMatchedHitPointersOut)
const SiStripElectronAlgo &	operator= (const SiStripElectronAlgo &)
bool	projectPhiBand (float chargeHypothesis, const reco::SuperClusterRef &superclusterIn)
	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_
edm::RefVector < SiStripRecHit2DCollection >	outputMatchedHits_neg_
edm::RefVector < SiStripRecHit2DCollection >	outputRphiHits_neg_
edm::RefVector < SiStripRecHit2DCollection >	outputRphiHits_pos_
edm::RefVector < SiStripRecHit2DCollection >	outputStereoHits_neg_
edm::RefVector < SiStripRecHit2DCollection >	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 58 of file SiStripElectronAlgo.cc.

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

SiStripElectronAlgo::~SiStripElectronAlgo

(

)

[virtual]

Definition at line 78 of file SiStripElectronAlgo.cc.

{
}

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,
		bool	stereo,
		bool	endcap
	)			[private]

Definition at line 357 of file SiStripElectronAlgo.cc.

References lat::endl(), false, edm::RangeMap< ID, C, P >::get(), id, edm::RangeMap< ID, C, P >::ids(), TrackerGeometry::idToDetUnit(), LogDebug, maxHitsOnDetId_, rphiHits_p_, stereoHits_p_, GeomDetEnumerators::TEC, theLayer, GeomDetEnumerators::TIB, GeomDetEnumerators::TID, GeomDetEnumerators::TOB, tracker_p_, and true.

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
  const std::vector<DetId> ids = (stereo ? stereoHits_p_->ids() : rphiHits_p_->ids());
  for (std::vector<DetId>::const_iterator id = ids.begin();  id != ids.end();  ++id) {

    // Get the hits on this detector id
    SiStripRecHit2DCollection::range hits = (stereo ? stereoHits_p_->get(*id) : rphiHits_p_->get(*id));

    // Count the number of hits on this detector id
    unsigned int numberOfHits = 0;
    for (SiStripRecHit2DCollection::const_iterator hit = hits.first;  hit != hits.second;  ++hit) {
      numberOfHits++;
      if (numberOfHits > maxHitsOnDetId_) { break; }
    }
      
    // 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::const_iterator hit = hits.first;  
           hit != hits.second;  ++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 = TIBDetId(*id).layer(); 
          if(TIBDetId(*id).stereo()==1) { isStereoDet = true ; }
        } else if
          (tracker_p_->idToDetUnit(hit->geographicalId())->type().subDetector() == GeomDetEnumerators::TOB) { 
          theDet = "TOB" ;
          theLayer = TOBDetId(*id).layer(); 
          if(TOBDetId(*id).stereo()==1) { isStereoDet = true ; }
        }else if
          (tracker_p_->idToDetUnit(hit->geographicalId())->type().subDetector() == GeomDetEnumerators::TID) { 
          theDet = "TID" ;
          theLayer = TIDDetId(*id).wheel();  // or ring  ?
          if(TIDDetId(*id).stereo()==1) { isStereoDet = true ; }
        }else if
          (tracker_p_->idToDetUnit(hit->geographicalId())->type().subDetector() == GeomDetEnumerators::TEC) { 
          theDet = "TEC" ;
          theLayer = TECDetId(*id).wheel();  // or ring or petal ?
          if(TECDetId(*id).stereo()==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 438 of file SiStripElectronAlgo.cc.

References lat::endl(), edm::RangeMap< ID, C, P >::get(), id, edm::RangeMap< ID, C, P >::ids(), LogDebug, matchedHits_p_, maxHitsOnDetId_, StripSubdetector::TIB, and StripSubdetector::TOB.

Referenced by projectPhiBand().

{
  
  // Loop over the detector ids
  const std::vector<DetId> ids = matchedHits_p_->ids() ;
  for (std::vector<DetId>::const_iterator id = ids.begin();  id != ids.end();  ++id) {
    
    // Get the hits on this detector id
    SiStripMatchedRecHit2DCollection::range hits = matchedHits_p_->get(*id) ;
    
    // Count the number of hits on this detector id
    unsigned int numberOfHits = 0;
    for (SiStripMatchedRecHit2DCollection::const_iterator hit = hits.first;  hit != hits.second;  ++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::const_iterator hit = hits.first;  hit != hits.second;  ++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
	)

Definition at line 150 of file SiStripElectronAlgo.cc.

References alongMomentum, chi2_neg_, chi2_pos_, correct_pT_neg_, correct_pT_pos_, HLT_VtxMuL3::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_, muonGeometry::perp(), 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(), state, 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);
  bool positronSuccess = projectPhiBand( 1., superclusterIn);

  // electron hypothesis did better than electron
  if ((electronSuccess  &&  !positronSuccess)  ||
      (electronSuccess  &&  positronSuccess  &&  redchi2_neg_ <= redchi2_pos_)) {
      
    // Initial uncertainty for tracking
    AlgebraicSymMatrix errors(5,1);  // makes identity 5x5 matrix, indexed from (1,1) to (5,5)
    errors(1,1) = 3.;      // uncertainty**2 in 1/momentum
    errors(2,2) = 0.01;    // uncertainty**2 in lambda (lambda == pi/2 - polar angle theta)
    errors(3,3) = 0.0001;  // uncertainty**2 in phi
    errors(4,4) = 0.01;    // uncertainty**2 in x_transverse (where x is in cm)
    errors(5,5) = 0.01;    // uncertainty**2 in y_transverse (where y is in cm)

    // 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"; 
    
    std::sort(outputHits_neg_.begin(), outputHits_neg_.end(),
              CompareHits(TrackingRecHitLessFromGlobalPosition(((TrackingGeometry*)(tracker_p_)), alongMomentum)));
    

    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();

    
    //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());
    
    TrajectoryStateTransform transformer;
    PTrajectoryStateOnDet* PTraj = transformer.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_));
    
    delete PTraj; // ShR: fix memory leak reported per perfrmance task force
    
    return true;
  }
  
  // positron hypothesis did better than electron
  if ((!electronSuccess  &&  positronSuccess)  ||
      (electronSuccess  &&  positronSuccess  &&  redchi2_neg_ > redchi2_pos_)) {
      
    // Initial uncertainty for tracking
    AlgebraicSymMatrix errors(5,1);  // makes identity 5x5 matrix, indexed from (1,1) to (5,5)
    errors(1,1) = 3.;      // uncertainty**2 in 1/momentum
    errors(2,2) = 0.01;    // uncertainty**2 in lambda (lambda == pi/2 - polar angle theta)
    errors(3,3) = 0.0001;  // uncertainty**2 in phi
    errors(4,4) = 0.01;    // uncertainty**2 in x_transverse (where x is in cm)
    errors(5,5) = 0.01;    // uncertainty**2 in y_transverse (where y is in cm)

    // 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"; 
    
    std::sort(outputHits_pos_.begin(), outputHits_pos_.end(),
              CompareHits(TrackingRecHitLessFromGlobalPosition(((TrackingGeometry*)(tracker_p_)), alongMomentum)));
    
    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();

    //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());
    
    TrajectoryStateTransform transformer;
    PTrajectoryStateOnDet* PTraj = transformer.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_));
    
    delete PTraj; // JED: fix of 2nd memory leak reported per perfrmance task force
    
    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
	)

Referenced by SiStripElectronProducer::produce().

bool SiStripElectronAlgo::projectPhiBand	(	float	chargeHypothesis,
		const reco::SuperClusterRef &	superclusterIn
	)			[private]

Definition at line 482 of file SiStripElectronAlgo.cc.

References chi2_neg_, chi2_pos_, coarseHitSelection(), coarseMatchedHitSelection(), correct_pT_neg_, correct_pT_pos_, funct::cos(), lat::endl(), relval_parameters_module::energy, false, TrackingRecHit::geographicalId(), hitUsed_, i, TrackerGeometry::idToDet(), TrackerGeometry::idToDetUnit(), innerhit_neg_, innerhit_pos_, intercept_neg_, intercept_pos_, MagneticField::inTesla(), j, BaseSiTrackerRecHit2DLocalPos::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_, PV3DBase< T, PVType, FrameType >::phi(), phi, phiBandWidth_, phiVsRSlope_neg_, phiVsRSlope_pos_, position_neg_, position_pos_, edm::RefVector< C, T, F >::push_back(), pZ_neg_, pZ_pos_, r, redchi2_neg_, redchi2_pos_, rphiHits_hp_, rphiKey_, funct::sin(), slope, slope_neg_, slope_pos_, python::multivaluedict::sort(), funct::sqrt(), stereoHits_hp_, stereoKey_, GeomDetEnumerators::TIB, GeomDetEnumerators::TOB, GeomDet::toGlobal(), tracker_p_, true, unwrapPhi(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), PV3DBase< T, PVType, FrameType >::z(), 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,     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, 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);
 
      const SiStripRecHit2D::ClusterRef & stereocluster=(*hit)->cluster();
    
      bool thisHitIsMatched = false ;

      if (!hitUsed_[*hit]) {
 
        unsigned int matcheduselist_size = matcheduselist.size();
        for (unsigned int i = 0;  i < matcheduselist_size;  i++) {
          if (matcheduselist[i]) {
            const SiStripRecHit2D::ClusterRef &  mystereocluster = matchedhitlist[i]->stereoHit()->cluster();
            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);

      const SiStripRecHit2D::ClusterRef &  monocluster=(*hit)->cluster();
    

      if (!hitUsed_[*hit]) {
        if( (tracker_p_->idToDetUnit((*hit)->geographicalId())->type().subDetector() == GeomDetEnumerators::TIB  && 
             (TIBDetId((*hit)->geographicalId()).layer()==1 || TIBDetId((*hit)->geographicalId()).layer()==2)) || 
            (tracker_p_->idToDetUnit((*hit)->geographicalId())->type().subDetector() == GeomDetEnumerators::TOB  
             && (TOBDetId((*hit)->geographicalId()).layer()==1 || TOBDetId((*hit)->geographicalId()).layer()==2)) ) {
          bool thisHitIsMatched = false ;
          unsigned int matcheduselist_size = matcheduselist.size();
          for (unsigned int i = 0;  i < matcheduselist_size;  i++) {
            if (matcheduselist[i]) {
              const SiStripRecHit2D::ClusterRef &  mymonocluster = matchedhitlist[i]->monoHit()->cluster();
              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  && 
               fabs(z - zFit) < 12.)  ||
              (tracker_p_->idToDetUnit((*hit)->geographicalId())->type().subDetector() == GeomDetEnumerators::TOB  && 
               fabs(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 " ;

  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";
  // 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() ;
      debugstr5<< " DetId " << detIDUsed << "\n";
      detIdList.push_back(detIDUsed) ;
    }
  }
  debugstr5 << " There are " << detIdList.size() << " hits on detectors \n";
  // 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());

  debugstr5 << " There are " << detIdList.size() << " unique detectors \n";

  //now we have a list of unique detectors used



  debugstr5 << " Looping over detectors to choose best hit " << " \n";
  // 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()) {
                debugstr5 << " Found 2 hits on same Si Detector " 
                          << ((hit2)->geographicalId()).rawId() << "\n";
                // 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;
                debugstr5 << " Chi1 " << chi1 << " Chi2 " << chi2 <<"\n";
                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] ) { 
      const SiStripRecHit2D* hit = hitlist[i];
      double localchi2 = (philist[i]-(rlist[i]-scr)*phiVsRSlope)*(philist[i]-(rlist[i]-scr)*phiVsRSlope)*w2list[i] ;
      if(localchi2 > chi2HitMax ) {
        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" ;
        uselist[i]=0 ;
      }
    }
  }

  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" ;


  // 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 ) {
    debugstr5 << " Too few hits "<< nHitsLeft 
              << " left - return false " << " \n";
    return false ;
  }
  LogDebug("") << debugstr5.str();

  
  // 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 ;  
        
        debugstr4 << " totalWeight " << totalWeight 
                  << " totalWeightRatio " << totalWeightRatio
                  << " localWeightRatio "<< localWeightRatio
                  << " phiBar " << phiBar
                  << " rBar " << rBar 
                  << " r2Bar " << r2Bar
                  << " phirBar " << phirBar 
                  << " \n ";


      } // 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_) {
      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";
      uselist[biggest_index] = false;
    }
    else {
      done = true;
    }
  } // end loop over trial fits

  debugstr4 <<" Fit " 
            << " Intercept  " << intercept
            << " slope " << slope 
            << " chi2 " << chi2
            << " \n" ;

  LogDebug("") <<   debugstr4.str() ;


  // 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
  edm::RefVector<SiStripRecHit2DCollection> outputRphiHits;
  edm::RefVector<SiStripRecHit2DCollection> outputStereoHits;

  typedef edm::Ref<SiStripRecHit2DCollection> 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(SiStripRecHit2DRef(*stereoHits_hp_, stereoKey_[hit]));
      }
      else if (typelist[i] == 1) {
        numberOfBarrelRphiHits++;

        // Copy this hit for the TrajectorySeed
        outputHits.push_back(hit);
        outputRphiHits.push_back(SiStripRecHit2DRef(*rphiHits_hp_, rphiKey_[hit]));
      }
      else if (typelist[i] == 2) {
        numberOfEndcapZphiHits++;

        // Copy this hit for the TrajectorySeed
        outputHits.push_back(hit);
        outputRphiHits.push_back(SiStripRecHit2DRef(*rphiHits_hp_, rphiKey_[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 [inline, private]

Definition at line 110 of file SiStripElectronAlgo.h.

Referenced by projectPhiBand().

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

---

Member Data Documentation

double SiStripElectronAlgo::chi2_neg_ [private]

Definition at line 178 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

double SiStripElectronAlgo::chi2_pos_ [private]

Definition at line 158 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

double SiStripElectronAlgo::correct_pT_neg_ [private]

Definition at line 180 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

double SiStripElectronAlgo::correct_pT_pos_ [private]

Definition at line 160 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

std::map<const TrackingRecHit*, bool> SiStripElectronAlgo::hitUsed_ [private]

Definition at line 142 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

const SiStripRecHit2D* SiStripElectronAlgo::innerhit_neg_ [private]

Definition at line 171 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

const SiStripRecHit2D* SiStripElectronAlgo::innerhit_pos_ [private]

Definition at line 149 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

double SiStripElectronAlgo::intercept_neg_ [private]

Definition at line 177 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

double SiStripElectronAlgo::intercept_pos_ [private]

Definition at line 157 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

const MagneticField* SiStripElectronAlgo::magneticField_p_ [private]

Definition at line 131 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

const edm::Handle<SiStripMatchedRecHit2DCollection>* SiStripElectronAlgo::matchedHits_hp_ [private]

Definition at line 135 of file SiStripElectronAlgo.h.

const SiStripMatchedRecHit2DCollection* SiStripElectronAlgo::matchedHits_p_ [private]

Definition at line 130 of file SiStripElectronAlgo.h.

Referenced by coarseMatchedHitSelection().

std::map<const TrackingRecHit*, bool> SiStripElectronAlgo::matchedHitUsed_ [private]

Definition at line 143 of file SiStripElectronAlgo.h.

Referenced by projectPhiBand().

std::map<const SiStripMatchedRecHit2D*, unsigned int> SiStripElectronAlgo::matchedKey_ [private]

Definition at line 140 of file SiStripElectronAlgo.h.

unsigned int SiStripElectronAlgo::maxHitsOnDetId_ [private]

Definition at line 119 of file SiStripElectronAlgo.h.

Referenced by coarseHitSelection(), and coarseMatchedHitSelection().

double SiStripElectronAlgo::maxNormResid_ [private]

Definition at line 122 of file SiStripElectronAlgo.h.

Referenced by projectPhiBand().

double SiStripElectronAlgo::maxReducedChi2_ [private]

Definition at line 124 of file SiStripElectronAlgo.h.

Referenced by projectPhiBand().

unsigned int SiStripElectronAlgo::minHits_ [private]

Definition at line 123 of file SiStripElectronAlgo.h.

Referenced by projectPhiBand().

GlobalVector SiStripElectronAlgo::momentum_neg_ [private]

Definition at line 170 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

GlobalVector SiStripElectronAlgo::momentum_pos_ [private]

Definition at line 148 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

int SiStripElectronAlgo::ndof_neg_ [private]

Definition at line 179 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

int SiStripElectronAlgo::ndof_pos_ [private]

Definition at line 159 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

unsigned SiStripElectronAlgo::numberOfBarrelRphiHits_neg_ [private]

Definition at line 184 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

unsigned int SiStripElectronAlgo::numberOfBarrelRphiHits_pos_ [private]

Definition at line 165 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

unsigned SiStripElectronAlgo::numberOfEndcapZphiHits_neg_ [private]

Definition at line 185 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

unsigned int SiStripElectronAlgo::numberOfEndcapZphiHits_pos_ [private]

Definition at line 166 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

unsigned int SiStripElectronAlgo::numberOfMatchedHits_pos_ [private]

Definition at line 163 of file SiStripElectronAlgo.h.

unsigned SiStripElectronAlgo::numberOfStereoHits_neg_ [private]

Definition at line 183 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

unsigned int SiStripElectronAlgo::numberOfStereoHits_pos_ [private]

Definition at line 164 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

double SiStripElectronAlgo::originUncertainty_ [private]

Definition at line 120 of file SiStripElectronAlgo.h.

Referenced by projectPhiBand().

std::vector<const TrackingRecHit*> SiStripElectronAlgo::outputHits_neg_ [private]

Definition at line 172 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

std::vector<const TrackingRecHit*> SiStripElectronAlgo::outputHits_pos_ [private]

Definition at line 150 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

edm::RefVector<SiStripRecHit2DCollection> SiStripElectronAlgo::outputMatchedHits_neg_ [private]

Definition at line 153 of file SiStripElectronAlgo.h.

edm::RefVector<SiStripRecHit2DCollection> SiStripElectronAlgo::outputRphiHits_neg_ [private]

Definition at line 173 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

edm::RefVector<SiStripRecHit2DCollection> SiStripElectronAlgo::outputRphiHits_pos_ [private]

Definition at line 151 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

edm::RefVector<SiStripRecHit2DCollection> SiStripElectronAlgo::outputStereoHits_neg_ [private]

Definition at line 174 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

edm::RefVector<SiStripRecHit2DCollection> SiStripElectronAlgo::outputStereoHits_pos_ [private]

Definition at line 152 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

double SiStripElectronAlgo::phiBandWidth_ [private]

Definition at line 121 of file SiStripElectronAlgo.h.

Referenced by projectPhiBand().

double SiStripElectronAlgo::phiVsRSlope_neg_ [private]

Definition at line 175 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

double SiStripElectronAlgo::phiVsRSlope_pos_ [private]

Definition at line 155 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

GlobalPoint SiStripElectronAlgo::position_neg_ [private]

Definition at line 169 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

GlobalPoint SiStripElectronAlgo::position_pos_ [private]

Definition at line 147 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

double SiStripElectronAlgo::pZ_neg_ [private]

Definition at line 181 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

double SiStripElectronAlgo::pZ_pos_ [private]

Definition at line 161 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

double SiStripElectronAlgo::redchi2_neg_ [private]

Definition at line 168 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

double SiStripElectronAlgo::redchi2_pos_ [private]

Definition at line 146 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

const edm::Handle<SiStripRecHit2DCollection>* SiStripElectronAlgo::rphiHits_hp_ [private]

Definition at line 133 of file SiStripElectronAlgo.h.

Referenced by projectPhiBand().

const SiStripRecHit2DCollection* SiStripElectronAlgo::rphiHits_p_ [private]

Definition at line 128 of file SiStripElectronAlgo.h.

Referenced by coarseHitSelection().

std::map<const SiStripRecHit2D*, unsigned int> SiStripElectronAlgo::rphiKey_ [private]

Definition at line 138 of file SiStripElectronAlgo.h.

Referenced by projectPhiBand().

double SiStripElectronAlgo::slope_neg_ [private]

Definition at line 176 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

double SiStripElectronAlgo::slope_pos_ [private]

Definition at line 156 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

const edm::Handle<SiStripRecHit2DCollection>* SiStripElectronAlgo::stereoHits_hp_ [private]

Definition at line 134 of file SiStripElectronAlgo.h.

Referenced by projectPhiBand().

const SiStripRecHit2DCollection* SiStripElectronAlgo::stereoHits_p_ [private]

Definition at line 129 of file SiStripElectronAlgo.h.

Referenced by coarseHitSelection().

std::map<const SiStripRecHit2D*, unsigned int> SiStripElectronAlgo::stereoKey_ [private]

Definition at line 139 of file SiStripElectronAlgo.h.

Referenced by projectPhiBand().

const TrackerGeometry* SiStripElectronAlgo::tracker_p_ [private]

Definition at line 127 of file SiStripElectronAlgo.h.

Referenced by coarseHitSelection(), findElectron(), and projectPhiBand().

double SiStripElectronAlgo::zVsRSlope_neg_ [private]

Definition at line 182 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

double SiStripElectronAlgo::zVsRSlope_pos_ [private]

Definition at line 162 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

---

The documentation for this class was generated from the following files:

• RecoEgamma/EgammaElectronAlgos/interface/SiStripElectronAlgo.h
• RecoEgamma/EgammaElectronAlgos/src/SiStripElectronAlgo.cc

---