CMS 3D CMS Logo

Public Member Functions | Private Member Functions | Private Attributes

SiStripElectronAlgo Class Reference

#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 SiStripElectronAlgooperator= (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 SiStripRecHit2Dinnerhit_neg_
const SiStripRecHit2Dinnerhit_pos_
double intercept_neg_
double intercept_pos_
const MagneticFieldmagneticField_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< SiStripRecHit2DoutputMatchedHits_neg_
std::vector< SiStripRecHit2DoutputRphiHits_neg_
std::vector< SiStripRecHit2DoutputRphiHits_pos_
std::vector< SiStripRecHit2DoutputStereoHits_neg_
std::vector< SiStripRecHit2DoutputStereoHits_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 SiStripRecHit2DCollectionrphiHits_p_
std::map< const
SiStripRecHit2D *, unsigned
int > 
rphiKey_
double slope_neg_
double slope_pos_
const edm::Handle
< SiStripRecHit2DCollection > * 
stereoHits_hp_
const SiStripRecHit2DCollectionstereoHits_p_
std::map< const
SiStripRecHit2D *, unsigned
int > 
stereoKey_
const TrackerGeometrytracker_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 edmNew::DetSet< T >::begin(), edmNew::DetSetVector< T >::begin(), edmNew::DetSet< T >::detId(), edmNew::DetSetVector< T >::end(), edmNew::DetSet< T >::end(), funct::false, ExpressReco_HICollisions_FallBack::id, TrackerGeometry::idToDetUnit(), TIBDetId::layer(), TOBDetId::layer(), LogDebug, maxHitsOnDetId_, rphiHits_p_, edmNew::DetSet< T >::size(), stereoHits_p_, sistripsummary::TEC, sistripsummary::TIB, sistripsummary::TID, sistripsummary::TOB, tracker_p_, funct::true, TIDDetId::wheel(), and TECDetId::wheel().

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 = 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 435 of file SiStripElectronAlgo.cc.

References edmNew::DetSet< T >::begin(), edmNew::DetSetVector< T >::begin(), edmNew::DetSetVector< T >::end(), edmNew::DetSet< T >::end(), LogDebug, matchedHits_p_, maxHitsOnDetId_, sistripsummary::TIB, and sistripsummary::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 
)

Definition at line 150 of file SiStripElectronAlgo.cc.

References alongMomentum, chi2_neg_, chi2_pos_, correct_pT_neg_, correct_pT_pos_, ExpressReco_HICollisions_FallBack::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(), 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(), evf::utils::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 
)

Definition at line 98 of file SiStripElectronAlgo.cc.

References cmsDriverOptions::counter, edmNew::DetSetVector< T >::data(), hitUsed_, LogDebug, magneticField_p_, matchedHits_p_, matchedHitUsed_, matchedKey_, edm::ESHandle< T >::product(), edm::Handle< T >::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 
) [private]

Definition at line 479 of file SiStripElectronAlgo.cc.

References abs, ExpressReco_HICollisions_FallBack::chi2, chi2_neg_, chi2_pos_, coarseHitSelection(), coarseMatchedHitSelection(), correct_pT_neg_, correct_pT_pos_, funct::cos(), generateEDF::done, relval_parameters_module::energy, funct::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, position_neg_, position_pos_, pZ_neg_, pZ_pos_, csvReporter::r, redchi2_neg_, redchi2_pos_, funct::sin(), slope, slope_neg_, slope_pos_, python::multivaluedict::sort(), mathSSE::sqrt(), GeomDetEnumerators::TIB, GeomDetEnumerators::TOB, GeomDet::toGlobal(), tracker_p_, funct::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  && 
               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 " ;

  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
  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 [inline, private]

Definition at line 110 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;
      }

Member Data Documentation

Definition at line 178 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 158 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 180 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

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(), prepareEvent(), and projectPhiBand().

Definition at line 171 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 149 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 177 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 157 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 131 of file SiStripElectronAlgo.h.

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

Definition at line 135 of file SiStripElectronAlgo.h.

Definition at line 130 of file SiStripElectronAlgo.h.

Referenced by coarseMatchedHitSelection(), and prepareEvent().

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

Definition at line 143 of file SiStripElectronAlgo.h.

Referenced by prepareEvent(), and projectPhiBand().

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

Definition at line 140 of file SiStripElectronAlgo.h.

Referenced by prepareEvent().

unsigned int SiStripElectronAlgo::maxHitsOnDetId_ [private]

Definition at line 119 of file SiStripElectronAlgo.h.

Referenced by coarseHitSelection(), and coarseMatchedHitSelection().

Definition at line 122 of file SiStripElectronAlgo.h.

Referenced by projectPhiBand().

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

Definition at line 170 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 148 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 179 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 159 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 184 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 165 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 185 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 166 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 163 of file SiStripElectronAlgo.h.

Definition at line 183 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 164 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 120 of file SiStripElectronAlgo.h.

Referenced by projectPhiBand().

Definition at line 172 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 150 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 153 of file SiStripElectronAlgo.h.

Definition at line 173 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 151 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 174 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 152 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 121 of file SiStripElectronAlgo.h.

Referenced by projectPhiBand().

Definition at line 175 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 155 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 169 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

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

Definition at line 168 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 146 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 133 of file SiStripElectronAlgo.h.

Referenced by prepareEvent().

Definition at line 128 of file SiStripElectronAlgo.h.

Referenced by coarseHitSelection(), and prepareEvent().

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

Definition at line 138 of file SiStripElectronAlgo.h.

Referenced by prepareEvent().

Definition at line 176 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 156 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 134 of file SiStripElectronAlgo.h.

Referenced by prepareEvent().

Definition at line 129 of file SiStripElectronAlgo.h.

Referenced by coarseHitSelection(), and prepareEvent().

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

Definition at line 139 of file SiStripElectronAlgo.h.

Referenced by prepareEvent().

Definition at line 182 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().

Definition at line 162 of file SiStripElectronAlgo.h.

Referenced by findElectron(), and projectPhiBand().