#include <HICMeasurementEstimator.h>

Inheritance diagram for HICMeasurementEstimator:

Public Member Functions
virtual void	chooseCuts (int &i)
virtual HICMeasurementEstimator *	clone () const
virtual std::pair< bool, double >	estimate (const TrajectoryStateOnSurface &, const TransientTrackingRecHit &) const
bool	estimate (const TrajectoryStateOnSurface &ts, const BoundPlane &plane) const
template<unsigned int D>
std::pair< bool, double >	estimate (const TrajectoryStateOnSurface &, const TransientTrackingRecHit &) const
virtual int	getDetectorCode (const DetLayer *a)
const MagneticField *	getField ()
cms::HICConst *	getHICConst ()
	HICMeasurementEstimator (const GeometricSearchTracker theTracker0, const MagneticField mf, double maxChi2, double nSigma=3.)
virtual std::vector< double >	setCuts (Trajectory &traj, const DetLayer *theCurrentLayer)
void	setHICConst (cms::HICConst *hh)
virtual void	setHICDetMap ()
virtual void	setLastLayer (DetLayer *&ll)
virtual void	setMult (int aMult=1)
virtual void	setSign (int &i)
Private Attributes
std::vector< BarrelDetLayer * >	bl
const MagneticField *	field
std::vector< ForwardDetLayer * >	fneg
std::vector< ForwardDetLayer * >	fpos
std::map< const DetLayer *, int >	theBackward
std::map< const DetLayer *, int >	theBarrel
double	theChi2Cut
int	theCutType
const DetLayer *	theFirstLayer
std::map< const DetLayer *, int >	theForward
cms::HICConst *	theHICConst
const DetLayer *	theLastLayer
const DetLayer *	theLayer
int	theLowMult
double	theMaxChi2
double	theNewCut
double	theNewCutB
int	theNSigma
double	thePhiBound
double	thePhiBoundMean
double	thePhiCut
double	thePhiCutMean
double	thePhiWin
double	thePhiWinB
double	thePhiWinMean
double	thePhiWinMeanB
int	theSign
const GeometricSearchTracker *	theTracker
int	theTrajectorySize
double	theZBound
double	theZBoundMean
double	theZCut
double	theZCutMean
double	theZWin
double	theZWinB
double	theZWinMean
double	theZWinMeanB

Detailed Description

A Chi2 Measurement Estimator. Computhes the Chi^2 of a TrajectoryState with a RecHit or a BoundPlane. The TrajectoryState must have errors. Works for any RecHit dimension. Ported from ORCA.

Date:: 2009/08/15 07:31:43

Revision:: 1.2

Author:: todorov, cerati

Definition at line 22 of file HICMeasurementEstimator.h.

Constructor & Destructor Documentation

HICMeasurementEstimator::HICMeasurementEstimator	(	const GeometricSearchTracker *	theTracker0,
		const MagneticField *	mf,
		double	maxChi2,
		double	nSigma = `3.`
	)		`[inline, explicit]`

Construct with cuts on chi2 and nSigma. The cut on Chi2 is used to define the acceptance of RecHits. The errors of the trajectory state are multiplied by nSigma to define acceptance of BoundPlane and maximalLocalDisplacement.

Definition at line 30 of file HICMeasurementEstimator.h.

References GeometricSearchTracker::barrelLayers(), bl, field, fneg, fpos, GeometricSearchTracker::negForwardLayers(), GeometricSearchTracker::posForwardLayers(), setHICDetMap(), and theTracker.

Referenced by clone().

                                                               :
  Chi2MeasurementEstimatorBase(maxChi2,nSigma)
    {
      theTracker = theTracker0;
      bl = theTracker->barrelLayers();
      fpos = theTracker->posForwardLayers();
      fneg = theTracker->negForwardLayers();
      field = mf;
      setHICDetMap();
    }

Member Function Documentation

void HICMeasurementEstimator::chooseCuts ( int & i ) [virtual]

Definition at line 278 of file HICMeasurementEstimator.cc.

References i.

Referenced by HICTrajectoryBuilder::findCompatibleMeasurements(), and HICTrajectoryBuilder::updateTrajectory().

{

      theChi2Cut = theNewCut;
//      cout<<" Choose Chi2Cut "<<theChi2Cut<<endl;
      
   if( i == 1 )
   {
      thePhiBound = thePhiWin;
      theZBound = theZWin;
      thePhiBoundMean = thePhiWinMean;
      theZBoundMean = theZWinMean;
//      cout<<" HICMeasurementEstimator::chooseCuts  "<<i<<" "<<thePhiBound<<" "<<theZBound<<endl;
   } 
      if( i == 2 )
   {
      thePhiBound = thePhiWinB;
      theZBound = theZWinB;
      thePhiBoundMean = thePhiWinMean;
      theZBoundMean = theZWinMean;
//      cout<<" HICMeasurementEstimator::chooseCuts  "<<i<<" "<<thePhiBound<<" "<<theZBound<<endl;
   }  
      if( i == 3 )
   {
      thePhiBound = thePhiCut;
      theZBound = theZCut;
      thePhiBoundMean = thePhiCutMean;
      theZBoundMean = theZCutMean;
//      cout<<" HICMeasurementEstimator::chooseCuts  "<<i<<" "<<thePhiBound<<" "<<theZBound<<endl;
   }
     
   theCutType = i;
}

virtual HICMeasurementEstimator* HICMeasurementEstimator::clone ( void ) const [inline, virtual]

Implements MeasurementEstimator.

Definition at line 71 of file HICMeasurementEstimator.h.

References HICMeasurementEstimator().

                                                 {
    return new HICMeasurementEstimator(*this);
  }

std::pair< bool, double > HICMeasurementEstimator::estimate	(	const TrajectoryStateOnSurface &	ts,
		const TransientTrackingRecHit &	hit
	)		const `[virtual]`

Returns pair( true, value) if the TrajectoryStateOnSurface is compatible with the RecHit, and pair( false, value) if it is not compatible. The TrajectoryStateOnSurface must be on the same Surface as the RecHit. For an estimator where there is no value computed, e.g. fixed window estimator, only the first(bool) part is of interest.

Implements Chi2MeasurementEstimatorBase.

Definition at line 12 of file HICMeasurementEstimator.cc.

References gather_cfg::cout, TrackingRecHit::dimension(), Exception, and TrajectoryStateOnSurface::isValid().

Referenced by HICTrajectoryBuilder::findCompatibleMeasurements(), and HICTrajectoryBuilder::updateTrajectory().

                                                                                 {
    std::pair<bool,double> flag(false,0.);
    if(!tsos.isValid()) {
#ifdef DEBUG
          std::cout<<" HICMeasurementEstimator::estimate::trajectory is not valid "<<std::endl;
#endif
    return flag;
    }

    switch (aRecHit.dimension()) {
        case 1: return estimate<1>(tsos,aRecHit);
        case 2: return estimate<2>(tsos,aRecHit);
        case 3: return estimate<3>(tsos,aRecHit);
        case 4: return estimate<4>(tsos,aRecHit);
        case 5: return estimate<5>(tsos,aRecHit);
    }
    throw cms::Exception("RecHit of invalid size (not 1,2,3,4,5)");
}

template<unsigned int D>

std::pair<bool,double> HICMeasurementEstimator::estimate	(	const TrajectoryStateOnSurface &	ts,
		const TransientTrackingRecHit &	hit
	)		const `[virtual]`

Returns pair( true, value) if the TrajectoryStateOnSurface is compatible with the RecHit, and pair( false, value) if it is not compatible. The TrajectoryStateOnSurface must be on the same Surface as the RecHit. For an estimator where there is no value computed, e.g. fixed window estimator, only the first(bool) part is of interest.

Implements Chi2MeasurementEstimatorBase.

bool HICMeasurementEstimator::estimate	(	const TrajectoryStateOnSurface &	ts,
		const BoundPlane &	plane
	)		const `[virtual]`

Returns true if the TrajectoryStateOnSurface is compatible with the BoundPlane, false otherwise. The TrajectoryStateOnSurface must be on the plane.

Reimplemented from Chi2MeasurementEstimatorBase.

Definition at line 98 of file HICMeasurementEstimator.cc.

References BoundSurface::bounds(), gather_cfg::cout, TrajectoryStateOnSurface::globalPosition(), Bounds::length(), PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::phi(), pi, GloballyPositioned< T >::position(), Bounds::width(), and PV3DBase< T, PVType, FrameType >::z().

{

//  cout<<" start estimate plane "<<endl;
  double pi = 4.*atan(1.);
  double twopi = 2.*pi;
  float theZError = plane.bounds().length() + 4.;
  float thePhiError = 2.*plane.bounds().width()/plane.position().perp(); 
// Change 02.07.08
//  float thePhiError = 4.*plane.bounds().width()/plane.position().perp();

#ifdef DEBUG  
  cout<<" ======================================================================================== ";
  cout<<" Estimate detector::   tsos      :     detector   :   Error "<<endl;
  cout<<" R                 "<<ts.globalPosition().perp()<<" "<<plane.position().perp()<<" "<<theZError<<endl;
  cout<<" Phi               "<<ts.globalPosition().phi()<<" "<<plane.position().phi()<<" "<<thePhiError<<endl;
  cout<<" Z                 "<<ts.globalPosition().z()<<" "<<plane.position().z()<<" "<<theZError<<endl;
#endif

  bool flag = false;
  if(fabs(ts.globalPosition().perp()-plane.position().perp())<theZError){
   if(fabs(ts.globalPosition().z()-plane.position().z())<theZError){
   float phi1 = ts.globalPosition().phi();
   float phi2 = plane.position().phi();
   if(phi1<0.) phi1 = twopi+phi1;
   if(phi2<0.) phi2 = twopi+phi2;
   float dfi = fabs(phi1-phi2);
   if(dfi>pi) dfi = twopi-dfi;
      if(dfi<thePhiError) flag = true;
   }
  }
#ifdef DEBUG
    cout<<" Estimate = "<<flag<<endl;
#endif

  return flag;
  
}

int HICMeasurementEstimator::getDetectorCode ( const DetLayer * a ) [virtual]

Definition at line 312 of file HICMeasurementEstimator.cc.

References Reference_intrackfit_cff::barrel, Reference_intrackfit_cff::endcap, DetLayer::location(), GloballyPositioned< T >::position(), edm::second(), and GeometricSearchDet::surface().

{
     int layer = 0;
     if( a->location() == GeomDetEnumerators::barrel ) 
     {
           layer = (*theBarrel.find(a)).second;
     }     
     if( a->location() == GeomDetEnumerators::endcap ) 
     {
           if( a->surface().position().z() > 0. ) { layer = 100+(*theForward.find(a)).second;}
           if( a->surface().position().z() < 0. ) { layer = -100-(*theBackward.find(a)).second;}           
     }
     return layer;
}

const MagneticField* HICMeasurementEstimator::getField ( ) [inline]

Definition at line 69 of file HICMeasurementEstimator.h.

References field.

Referenced by HICTrajectoryBuilder::findCompatibleMeasurements(), and HICTrajectoryBuilder::updateTrajectory().

{return field;}

cms::HICConst* HICMeasurementEstimator::getHICConst ( ) [inline]

Definition at line 55 of file HICMeasurementEstimator.h.

References theHICConst.

Referenced by HICTrajectoryBuilder::trajectories().

{return theHICConst;}

vector< double > HICMeasurementEstimator::setCuts	(	Trajectory &	traj,
		const DetLayer *	theCurrentLayer
	)		`[virtual]`

Definition at line 138 of file HICMeasurementEstimator.cc.

References a, b, Reference_intrackfit_cff::barrel, gather_cfg::cout, Trajectory::data(), GeomDetEnumerators::endcap, Reference_intrackfit_cff::endcap, edm::first(), DetLayer::location(), Trajectory::measurements(), GeomDetEnumerators::PixelBarrel, GeomDetEnumerators::PixelEndcap, GloballyPositioned< T >::position(), edm::second(), DetLayer::subDetector(), and GeometricSearchDet::surface().

Referenced by HICTrajectoryBuilder::findCompatibleMeasurements(), and HICTrajectoryBuilder::updateTrajectory().

{
     vector<double> theCuts;
     const DetLayer* a = traj.data().back().layer();
     const DetLayer* first = traj.data().front().layer();
//     const DetLayer* last = traj.data().front().layer();
     thePhiWinMean = 0.;
     theZWinMean = 0.;
     thePhiWin = 0.;
     theZWin = 0.;
     theNewCut = 11.; // change 5->10 03.07.2008 //10-11 23.06.09
     theNewCutB = 5.;
     
     thePhiWinMeanB = 0.002;
     theZWinMeanB = 0.;
     thePhiWinB = 0.008;
     theZWinB = 17.;
          
     theZCutMean = 0.;
     thePhiCutMean = 0.;
     thePhiCut = 0.;
     theZCut = 0.;
     
     theLayer = b;
     theLastLayer = a;

     
     theTrajectorySize = traj.data().size();
     
     
     if( theBarrel.size() == 0 || theForward.size() == 0 )
     {
#ifdef DEBUG
       cout<<" HICMeasurementEstimator::setCuts:: no datector map "<<endl;
#endif
        return theCuts;
     }
     
     if( a->location() == GeomDetEnumerators::barrel )
     {
       if( first->location() == GeomDetEnumerators::barrel ) 
       {
        thePhiWin = (*theHICConst).phiwinbar[(*theBarrel.find(first)).second][(*theBarrel.find(a)).second][(*theBarrel.find(b)).second];
        theZWin = (*theHICConst).zwinbar[(*theBarrel.find(first)).second][(*theBarrel.find(a)).second][(*theBarrel.find(b)).second];
        thePhiCut = (*theHICConst).phicutbar[(*theBarrel.find(first)).second][(*theBarrel.find(theLastLayer)).second][(*theBarrel.find(b)).second];
        theZCut = (*theHICConst).zcutbar[(*theBarrel.find(first)).second][(*theBarrel.find(theLastLayer)).second][(*theBarrel.find(b)).second];
//        cout<<" Barrel first -Barrel cuts::layers "<<(*theBarrel.find(first)).second<<" "<<(*theBarrel.find(a)).second<<" "<<(*theBarrel.find(b)).second<<endl; 
//        cout<<" Barrel first -Barrel cuts "<< thePhiWin<<" "<<theZWin <<" "<<thePhiCut <<" "<<theZCut<<endl;
       }
         else
         {
          if(first->surface().position().z() > 0. )
          {
          thePhiWin = (*theHICConst).phiwinfbb[(*theForward.find(first)).second][(*theBarrel.find(a)).second][(*theBarrel.find(b)).second];
          theZWin = (*theHICConst).zwinfbb[(*theForward.find(first)).second][(*theBarrel.find(a)).second][(*theBarrel.find(b)).second];
          thePhiCut = (*theHICConst).phicutfbb[(*theForward.find(first)).second][(*theBarrel.find(theLastLayer)).second][(*theBarrel.find(b)).second];
          theZCut = (*theHICConst).zcutfbb[(*theForward.find(first)).second][(*theBarrel.find(theLastLayer)).second][(*theBarrel.find(b)).second];
//          cout<<" Endcap first positive -Barrel cuts::layers "<<(*theForward.find(first)).second<<" "<<(*theBarrel.find(a)).second<<" "<<(*theBarrel.find(b)).second<<endl;       
//          cout<<" Endcap first positive -Barrel cuts "<< thePhiWin<<" "<<theZWin <<" "<<thePhiCut <<" "<<theZCut<<endl;
          } else {
          thePhiWin = (*theHICConst).phiwinfbb[(*theBackward.find(first)).second][(*theBarrel.find(a)).second][(*theBarrel.find(b)).second];
          theZWin = (*theHICConst).zwinfbb[(*theBackward.find(first)).second][(*theBarrel.find(a)).second][(*theBarrel.find(b)).second];
          thePhiCut = (*theHICConst).phicutfbb[(*theBackward.find(first)).second][(*theBarrel.find(theLastLayer)).second][(*theBarrel.find(b)).second];
          theZCut = (*theHICConst).zcutfbb[(*theBackward.find(first)).second][(*theBarrel.find(theLastLayer)).second][(*theBarrel.find(b)).second];
//          cout<<" Endcap first negative -Barrel cuts::layers "<<(*theBackward.find(first)).second<<" "<<(*theBarrel.find(a)).second<<" "<<(*theBarrel.find(b)).second<<endl;
//          cout<<" Endcap first negative -Barrel cuts "<< thePhiWin<<" "<<theZWin <<" "<<thePhiCut <<" "<<theZCut<<endl;
         }
         } //


        theCuts.push_back(thePhiWin); theCuts.push_back(theZWin);
        theCuts.push_back(thePhiCut); theCuts.push_back(theZCut);

        return theCuts;
     }
     if( a->location() == GeomDetEnumerators::endcap && b->location() == GeomDetEnumerators::endcap)
     {
        if( a->surface().position().z() > 0. )
        {
        thePhiWin = (*theHICConst).phiwinfrw[(*theForward.find(first)).second][(*theForward.find(a)).second][(*theForward.find(b)).second];
        theZWin = (*theHICConst).zwinfrw[(*theForward.find(first)).second][(*theForward.find(a)).second][(*theForward.find(b)).second];
        thePhiCut = (*theHICConst).phicutfrw[(*theForward.find(first)).second][(*theForward.find(theLastLayer)).second][(*theForward.find(b)).second];
        theZCut = (*theHICConst).zcutfrw[(*theForward.find(first)).second][(*theForward.find(theLastLayer)).second][(*theForward.find(b)).second];
        }
          else
          {
           thePhiWin = (*theHICConst).phiwinfrw[(*theBackward.find(first)).second][(*theBackward.find(a)).second][(*theBackward.find(b)).second];
           theZWin = (*theHICConst).zwinfrw[(*theBackward.find(first)).second][(*theBackward.find(a)).second][(*theBackward.find(b)).second];
           thePhiCut = (*theHICConst).phicutfrw[(*theBackward.find(first)).second][(*theBackward.find(theLastLayer)).second][(*theBackward.find(b)).second];
           theZCut = (*theHICConst).zcutfrw[(*theBackward.find(first)).second][(*theBackward.find(theLastLayer)).second][(*theBackward.find(b)).second];
          }
        
        if( theLowMult == 1 )
        {
        if( b->subDetector() == GeomDetEnumerators::PixelEndcap ) theNewCut = 20.;
        if( traj.measurements().size() == 1 ) theNewCut = 20.; 
        theNewCutB = 30.;
        }
        
        thePhiWinMeanB = 0.004;
        thePhiWinB = 0.05;


        theCuts.push_back(thePhiWin); theCuts.push_back(theZWin);
        theCuts.push_back(thePhiCut); theCuts.push_back(theZCut);
        
        return theCuts;
     }
     if( a->location() == GeomDetEnumerators::endcap && b->location() == GeomDetEnumerators::barrel )
     {
       
       if( a->surface().position().z() > 0. )
       {
        thePhiWin = (*theHICConst).phiwinbfrw[(*theForward.find(first)).second][(*theForward.find(a)).second][(*theBarrel.find(b)).second];
        theZWin = (*theHICConst).zwinbfrw[(*theForward.find(first)).second][(*theForward.find(a)).second][(*theBarrel.find(b)).second];
        thePhiCut = (*theHICConst).phicutbfrw[(*theForward.find(first)).second][(*theForward.find(theLastLayer)).second][(*theBarrel.find(b)).second];
        theZCut = (*theHICConst).zcutbfrw[(*theForward.find(first)).second][(*theForward.find(theLastLayer)).second][(*theBarrel.find(b)).second];
       }
          else
          {
            thePhiWin = (*theHICConst).phiwinbfrw[(*theBackward.find(first)).second][(*theBackward.find(a)).second][(*theBarrel.find(b)).second];
            theZWin = (*theHICConst).zwinbfrw[(*theBackward.find(first)).second][(*theBackward.find(a)).second][(*theBarrel.find(b)).second];
            thePhiCut = (*theHICConst).phicutbfrw[(*theBackward.find(first)).second][(*theBackward.find(theLastLayer)).second][(*theBarrel.find(b)).second];
            theZCut = (*theHICConst).zcutbfrw[(*theBackward.find(first)).second][(*theBackward.find(theLastLayer)).second][(*theBarrel.find(b)).second];
          }     
        
        if( b->subDetector() ==  GeomDetEnumerators::PixelBarrel) theNewCut = 20.;

        thePhiWinMeanB = 0.004;
        thePhiWinB = 0.016;

        theCuts.push_back(thePhiWin); theCuts.push_back(theZWin);
        theCuts.push_back(thePhiCut); theCuts.push_back(theZCut);
        
        return theCuts;
     }
//     cout<<" HICMeasurementEstimator::setCuts::Error: unknown detector layer "<<endl;
     return theCuts;
}

void HICMeasurementEstimator::setHICConst ( cms::HICConst * hh ) [inline]

Definition at line 57 of file HICMeasurementEstimator.h.

References theHICConst.

Referenced by cms::HITrackVertexMaker::produceTracks().

{theHICConst = hh;}

void HICMeasurementEstimator::setHICDetMap ( ) [virtual]

Definition at line 327 of file HICMeasurementEstimator.cc.

References gather_cfg::cout.

Referenced by HICMeasurementEstimator().

{

#ifdef DEBUG
  std::cout<<" Set Detector Map... "<<std::endl;
#endif
  
  int ila=0;
  for ( std::vector<BarrelDetLayer*>::const_iterator ilayer = bl.begin(); ilayer != bl.end(); ilayer++)
  {
     theBarrel[(*ilayer)]=ila;
     ila++;
  }
//
// The same for forward part.
//

   int ilf1 = 0; 
   int ilf2 = 0;
   for ( vector<ForwardDetLayer*>::const_iterator ilayer = fpos.begin(); 
                                                    ilayer != fpos.end(); ilayer++)
   {
                theForward[(*ilayer)] = ilf1;
                ilf1++;
   }
   for ( vector<ForwardDetLayer*>::const_iterator ilayer = fneg.begin(); 
                                                    ilayer != fneg.end(); ilayer++)
   {
#ifdef DEBUG
            cout<<" HICDetectorMap::negative layers "<<(**ilayer).position().z()<<" "<<ilf2<<endl;
#endif
                   theBackward[(*ilayer)] = ilf2;
                   ilf2++;
   }

}