#include <OutInConversionSeedFinder.h>

Inheritance diagram for OutInConversionSeedFinder:

Public Member Functions
virtual void	makeSeeds (const edm::Handle< edm::View< reco::CaloCluster > > &allBc) const

virtual void	makeSeeds (const reco::CaloClusterPtr &aBC) const

	OutInConversionSeedFinder (const edm::ParameterSet &config)

virtual	~OutInConversionSeedFinder ()

Public Member Functions inherited from ConversionSeedFinder
	ConversionSeedFinder ()

	ConversionSeedFinder (const edm::ParameterSet &config)

const MeasurementTracker *	getMeasurementTracker () const

std::vector< const DetLayer * >	layerList () const

TrajectorySeedCollection	seeds ()

virtual void	setCandidate (float e, GlobalPoint pos) const

void	setEvent (const edm::Event &e)

void	setEventSetup (const edm::EventSetup &es)
	Initialize EventSetup objects at each event. More...

void	setMeasurementTracker (const MeasurementTracker *tracker) const

virtual	~ConversionSeedFinder ()

Private Types
typedef FreeTrajectoryState	FTS

typedef TrajectoryStateOnSurface	TSOS

Private Member Functions
void	completeSeed (const TrajectoryMeasurement &m1, FreeTrajectoryState &fts, const Propagator *, int layer) const

void	createSeed (const TrajectoryMeasurement &m1, const TrajectoryMeasurement &m2) const

FreeTrajectoryState	createSeedFTS (const TrajectoryMeasurement &m1, const TrajectoryMeasurement &m2) const

void	fillClusterSeeds (const reco::CaloClusterPtr &bc) const

GlobalPoint	fixPointRadius (const TrajectoryMeasurement &) const

MeasurementEstimator *	makeEstimator (const DetLayer *, float dphi) const

std::pair< FreeTrajectoryState, bool >	makeTrackState (int charge) const

void	startSeed (const FreeTrajectoryState &) const

Private Attributes
edm::ParameterSet	conf_

int	maxNumberOfOutInSeedsPerBC_

int	nSeedsPerBC_

float	the2ndHitdphi_

float	the2ndHitdzConst_

float	the2ndHitdznSigma_

std::vector < TrajectoryMeasurement >	theFirstMeasurements_

Additional Inherited Members
Protected Member Functions inherited from ConversionSeedFinder
void	findLayers () const

void	findLayers (const FreeTrajectoryState &fts) const

void	printLayer (int i) const

FreeTrajectoryState	trackStateFromClusters (int aCharge, const GlobalPoint &gpOrigine, PropagationDirection dir, float scaleFactor) const

Protected Attributes inherited from ConversionSeedFinder
edm::ParameterSet	conf_

PropagationDirection	dir_

float	theBCEnergy_

GlobalPoint	theBCPosition_

reco::BeamSpot	theBeamSpot_

edm::ESHandle < GeometricSearchTracker >	theGeomSearchTracker_

std::vector< const DetLayer * >	theLayerList_

const MeasurementTracker *	theMeasurementTracker_

edm::ESHandle< MagneticField >	theMF_

const Propagator *	thePropagatorAlongMomentum_

const Propagator *	thePropagatorOppositeToMomentum_

reco::CaloCluster *	theSC_

float	theSCenergy_

GlobalPoint	theSCPosition_

TrajectorySeedCollection	theSeeds_

const TrackingGeometry *	theTrackerGeom_

KFUpdator	theUpdator_

Detailed Description

Id:: OutInConversionSeedFinder.h,v 1.7 2008/02/15 16:47:15 nancy Exp

Date:: 2008/02/15 16:47:15

Revision:: 1.7

Author: Nancy Marinelli, U. of Notre Dame, US

Definition at line 34 of file OutInConversionSeedFinder.h.

Member Typedef Documentation

typedef FreeTrajectoryState OutInConversionSeedFinder::FTS

private

Definition at line 39 of file OutInConversionSeedFinder.h.

typedef TrajectoryStateOnSurface OutInConversionSeedFinder::TSOS

private

Definition at line 40 of file OutInConversionSeedFinder.h.

Constructor & Destructor Documentation

OutInConversionSeedFinder::OutInConversionSeedFinder ( const edm::ParameterSet & config )

Definition at line 30 of file OutInConversionSeedFinder.cc.

References LogDebug, maxNumberOfOutInSeedsPerBC_, the2ndHitdphi_, the2ndHitdzConst_, and the2ndHitdznSigma_.

                                                                                  : ConversionSeedFinder( conf ), conf_(conf)  
 {
 
   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder CTOR " << "\n";      
 
   //the2ndHitdphi_ = 0.01; 
   the2ndHitdphi_ = 0.03; 
   the2ndHitdzConst_ = 5.;
   the2ndHitdznSigma_ = 2.;
   maxNumberOfOutInSeedsPerBC_=50;    
    
     
 }

OutInConversionSeedFinder::~OutInConversionSeedFinder ( )

virtual

Definition at line 47 of file OutInConversionSeedFinder.cc.

References LogDebug.

                                                       {
   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder DTOR " << "\n";
  
 }

Member Function Documentation

void OutInConversionSeedFinder::completeSeed	(	const TrajectoryMeasurement &	m1,
		FreeTrajectoryState &	fts,
		const Propagator *	propagator,
		int	layer
	)		const

private

Definition at line 381 of file OutInConversionSeedFinder.cc.

References GeomDetEnumerators::barrel, createSeed(), ConversionSeedFinder::getMeasurementTracker(), i, DetLayer::location(), LogDebug, LayerMeasurements::measurements(), Propagator::propagationDirection(), TrajectoryMeasurement::recHit(), mathSSE::sqrt(), GeometricSearchDet::surface(), the2ndHitdphi_, the2ndHitdzConst_, the2ndHitdznSigma_, and ConversionSeedFinder::theLayerList_.

Referenced by startSeed().

                                                                          {
 
   //std::cout <<  "OutInConversionSeedFinder::completeSeed ilayer " << ilayer << "\n";
 
   MeasurementEstimator * newEstimator=0;
   const DetLayer * layer = theLayerList_[ilayer];
   
 
   if ( layer->location() == GeomDetEnumerators::barrel ) {
     // z error for 2nd hit is  2 sigma quadded with 5 cm
     LogDebug("OutInConversionSeedFinder") << " Barrel OutInConversionSeedFinder::completeSeed " << the2ndHitdznSigma_ << " " << the2ndHitdzConst_ << " " << the2ndHitdphi_ << "\n";
     float dz = sqrt(the2ndHitdznSigma_*the2ndHitdznSigma_*m1.recHit()->globalPositionError().czz() 
             + the2ndHitdzConst_*the2ndHitdzConst_);
     newEstimator =
       new ConversionBarrelEstimator(-the2ndHitdphi_, the2ndHitdphi_, -dz, dz);
   }
   else {
     LogDebug("OutInConversionSeedFinder") << " EndCap OutInConversionSeedFinder::completeSeed " << the2ndHitdznSigma_ << " " << the2ndHitdzConst_ << " " << the2ndHitdphi_ << "\n";
     // z error for 2nd hit is 2sigma quadded with 5 cm
     //float m1dr = m1.recHit().globalPositionError().rerr(m1.recHit().globalPosition());
     float m1dr = sqrt(m1.recHit()->localPositionError().yy());
     float dr = sqrt(the2ndHitdznSigma_*the2ndHitdznSigma_*m1dr*m1dr 
             + the2ndHitdzConst_*the2ndHitdznSigma_);
     // LogDebug("OutInConversionSeedFinder") << "second hit forward dr " << dr << " this hit " << m1dr << endl;
     newEstimator =
       new ConversionForwardEstimator(-the2ndHitdphi_, the2ndHitdphi_, dr);
   }
 
   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::completeSeed  ilayer " << ilayer <<  "\n"; 
   
   // Get the measurements consistent with the FTS and the Estimator
   TSOS tsos(fts, layer->surface() );
   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::completeSeed propagationDirection  " << int(propagator->propagationDirection() ) << "\n";               
   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::completeSeed pointer to estimator " << newEstimator << "\n";
 
   LayerMeasurements theLayerMeasurements_(this->getMeasurementTracker() );
   std::vector<TrajectoryMeasurement> measurements = theLayerMeasurements_.measurements( *layer, tsos, *propagator, *newEstimator);
   //std::cout << "OutInConversionSeedFinder::completeSeed Found " << measurements.size() << " second hits " << "\n";
   delete newEstimator;
 
   for(unsigned int i = 0; i < measurements.size(); ++i) {
     if( measurements[i].recHit()->isValid()  ) {
       createSeed(m1, measurements[i]);
     }
   }
 
 
 
   //LogDebug("OutInConversionSeedFinder") << "COMPLETED " << theSeeds_.size() << " SEEDS " << "\n";
 }

void OutInConversionSeedFinder::createSeed	(	const TrajectoryMeasurement &	m1,
		const TrajectoryMeasurement &	m2
	)		const

private

Definition at line 436 of file OutInConversionSeedFinder.cc.

References createSeedFTS(), TrajectoryMeasurement::estimate(), PV3DBase< T, PVType, FrameType >::eta(), TrajectoryStateOnSurface::freeTrajectoryState(), TrajectoryStateOnSurface::globalPosition(), TrajectoryStateOnSurface::isValid(), TrajectoryMeasurement::layer(), LogDebug, maxNumberOfOutInSeedsPerBC_, nSeedsPerBC_, oppositeToMomentum, PTrajectoryStateOnDet::parameters(), PV3DBase< T, PVType, FrameType >::perp(), TrajectoryStateTransform::persistentState(), PV3DBase< T, PVType, FrameType >::phi(), LocalTrajectoryParameters::position(), Propagator::propagate(), edm::OwnVector< T, P >::push_back(), TrajectoryMeasurement::recHit(), ConversionSeedFinder::thePropagatorOppositeToMomentum_, ConversionSeedFinder::theSeeds_, ConversionSeedFinder::theUpdator_, and KFUpdator::update().

Referenced by completeSeed().

                                                                {
 
   //std::cout  << "OutInConversionSeedFinder::createSeed  from hit1 " << m1.recHit()->globalPosition() << " r1 " << m1.recHit()->globalPosition().perp() << " and hit2 " << m2.recHit()->globalPosition() << " r2 " << m2.recHit()->globalPosition().perp() << "\n";
   
   FreeTrajectoryState fts = createSeedFTS(m1, m2);
 
   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::createSeed First point errors " <<m1.recHit()->parametersError() << "\n";
   LogDebug("OutInConversionSeedFinder") << "original cluster FTS " << fts <<"\n";
 
   //std::cout  << "OutInConversionSeedFinder::createSeed propagation dir " << int( thePropagatorOppositeToMomentum_->propagationDirection() ) << "\n"; 
   TrajectoryStateOnSurface state1 = thePropagatorOppositeToMomentum_->propagate(fts,  m1.recHit()->det()->surface());
 
   // LogDebug("OutInConversionSeedFinder") << "hit surface " << h1.det().surface().position() << endl;
   // LogDebug("OutInConversionSeedFinder") << "prop to " << typeid(h1.det().surface()).name() << endl;
   // LogDebug("OutInConversionSeedFinder") << "prop to first hit " << state1 << endl; 
   // LogDebug("OutInConversionSeedFinder") << "update to " << h1.globalPosition() << endl;
 
   if ( state1.isValid() ) {
     TrajectoryStateOnSurface updatedState1 = theUpdator_.update(state1,  *m1.recHit() );
 
     if ( updatedState1.isValid() ) {
       TrajectoryStateOnSurface state2 = thePropagatorOppositeToMomentum_->propagate(*updatedState1.freeTrajectoryState(),  m2.recHit()->det()->surface());
 
       if ( state2.isValid() ) {
 
     TrajectoryStateOnSurface updatedState2 = theUpdator_.update(state2, *m2.recHit() );
     TrajectoryMeasurement meas1(state1, updatedState1,  m1.recHit()  , m1.estimate(), m1.layer());
     TrajectoryMeasurement meas2(state2, updatedState2,  m2.recHit()  , m2.estimate(), m2.layer());
     
         edm::OwnVector<TrackingRecHit> myHits;
     myHits.push_back(meas1.recHit()->hit()->clone());
     myHits.push_back(meas2.recHit()->hit()->clone());
 
         if ( nSeedsPerBC_ >= maxNumberOfOutInSeedsPerBC_ ) return;
    
     TrajectoryStateTransform tsTransform;
     PTrajectoryStateOnDet* ptsod= tsTransform.persistentState(state2, meas2.recHit()->hit()->geographicalId().rawId()  );
 
     LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::createSeed new seed  from state " << state2.globalPosition()  <<  "\n";
     LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::createSeed new seed  ptsod " <<  ptsod->parameters().position() << " R " << ptsod->parameters().position().perp() << " phi " << ptsod->parameters().position().phi() << " eta " << ptsod->parameters().position().eta() << "\n" ;
     
     theSeeds_.push_back(TrajectorySeed( *ptsod, myHits, oppositeToMomentum ));
         nSeedsPerBC_++;
  
         delete ptsod;  
         
       }
     }
   }
 
 }

FreeTrajectoryState OutInConversionSeedFinder::createSeedFTS	(	const TrajectoryMeasurement &	m1,
		const TrajectoryMeasurement &	m2
	)		const

private

Definition at line 504 of file OutInConversionSeedFinder.cc.

References alpha, DeDxDiscriminatorTools::charge(), TrajectoryStateOnSurface::charge(), funct::cos(), fixPointRadius(), m, PV3DBase< T, PVType, FrameType >::perp(), phi, PV3DBase< T, PVType, FrameType >::phi(), TrajectoryMeasurement::predictedState(), ExpressReco_HICollisions_FallBack::pt, funct::sin(), ConversionSeedFinder::theMF_, ConversionSeedFinder::theSCenergy_, ConversionSeedFinder::theSCPosition_, and PV3DBase< T, PVType, FrameType >::theta().

Referenced by createSeed().

                                                                          {
 
   //std::cout  << "OutInConversionSeedFinder::createSeedFTS " << "\n";
 
   GlobalPoint xmeas = fixPointRadius(m1);
   GlobalPoint xvert = fixPointRadius(m2);
 
 
   float pt = theSCenergy_ * sin(theSCPosition_.theta());
   float pz = theSCenergy_ * cos(theSCPosition_.theta());
 
 
 
   // doesn't work at all for endcap, where r is badly measured
   //float dphidr = (p1.phi()-p2.phi())/(p1.perp()-p2.perp());
   //int charge = (dphidr > 0.) ? -1 : 1;
   int charge = m2.predictedState().charge();
 
   double BInTesla = theMF_->inTesla(xmeas).z();
   GlobalVector xdiff = xmeas -xvert;
  
   double phi= xdiff.phi();
   double pxOld = pt*cos(phi);
   double pyOld = pt*sin(phi);
   double RadCurv = 100*pt/(BInTesla*0.29979);
   double alpha = asin(0.5*xdiff.perp()/RadCurv);
   float ca = cos(charge*alpha);
   float sa = sin(charge*alpha);
   double pxNew =   ca*pxOld + sa*pyOld;
   double pyNew =  -sa*pxOld + ca*pyOld;
   GlobalVector pNew(pxNew, pyNew, pz);
 
   GlobalTrajectoryParameters gp(xmeas, pNew, charge, &(*theMF_) );
 
   AlgebraicSymMatrix55 m = AlgebraicMatrixID();
   m(0,0) = 0.05; m(1,1) = 0.02 ; m(2,2) = 0.007 ;
   m(3,3) = 10. ; m(4,4) = 10. ;
   return FreeTrajectoryState(gp, CurvilinearTrajectoryError(m));
 
 
 }

void OutInConversionSeedFinder::fillClusterSeeds ( const reco::CaloClusterPtr & bc ) const

private

negative charge state

positive charge state

Definition at line 176 of file OutInConversionSeedFinder.cc.

References makeTrackState(), edm::second(), startSeed(), and theFirstMeasurements_.

Referenced by makeSeeds().

                                                                                    {
 
   
   theFirstMeasurements_.clear();
   FreeTrajectoryState fts;  
 
   if ( makeTrackState(-1).second ) {
     fts = makeTrackState(-1).first;
     startSeed(fts);
   }
 
   if ( makeTrackState(1).second ) {
     fts = makeTrackState(1).first;
     startSeed(fts);
   }
 
 }

GlobalPoint OutInConversionSeedFinder::fixPointRadius ( const TrajectoryMeasurement & m1 ) const

private

Definition at line 550 of file OutInConversionSeedFinder.cc.

References GeomDetEnumerators::barrel, funct::cos(), TrajectoryMeasurement::layer(), DetLayer::location(), p1, p2, phi, PV3DBase< T, PVType, FrameType >::phi(), csvReporter::r, TrajectoryMeasurement::recHit(), funct::sin(), funct::tan(), ConversionSeedFinder::theSCPosition_, PV3DBase< T, PVType, FrameType >::theta(), theta(), z, and PV3DBase< T, PVType, FrameType >::z().

Referenced by createSeedFTS().

                                                                                             {
   GlobalPoint p1 = m1.recHit()->globalPosition();
   GlobalPoint p2;
   if(m1.layer()->location() == GeomDetEnumerators::barrel) {
     p2 = p1;
   } else {
     float z = p1.z();
     float phi = p1.phi();
     float theta = theSCPosition_.theta();
     float r = p1.z() * tan(theta);
     p2 = GlobalPoint(r*cos(phi), r*sin(phi), z);
     // LogDebug("OutInConversionSeedFinder") << "fixing point radius " << p2 << " from " << p1 << endl;
   }
   return p2;
 }

MeasurementEstimator * OutInConversionSeedFinder::makeEstimator	(	const DetLayer *	layer,
		float	dphi
	)		const

private

Definition at line 343 of file OutInConversionSeedFinder.cc.

References GeomDetEnumerators::barrel, GeomDetEnumerators::endcap, BoundDisk::innerRadius(), DetLayer::location(), LogDebug, BoundDisk::outerRadius(), PV3DBase< T, PVType, FrameType >::perp(), GloballyPositioned< T >::position(), csvReporter::r, Cylinder::radius(), ForwardDetLayer::specificSurface(), BarrelDetLayer::specificSurface(), ForwardDetLayer::surface(), ConversionSeedFinder::theBCPosition_, z, and PV3DBase< T, PVType, FrameType >::z().

Referenced by startSeed().

                                                                                                         {
  
   //std::cout  << "OutInConversionSeedFinder::makeEstimator  " << "\n";
 
   MeasurementEstimator * newEstimator=0;
 
   if (layer->location() == GeomDetEnumerators::barrel ) {
     
     const BarrelDetLayer * barrelLayer = dynamic_cast<const BarrelDetLayer*>(layer);
     LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeEstimator Barrel  r = " << barrelLayer->specificSurface().radius() << " " << "\n";        
     float r = barrelLayer->specificSurface().radius();
     float zrange = 15.* (1.-r/theBCPosition_.perp());
     newEstimator = new ConversionBarrelEstimator(-dphi, dphi, -zrange, zrange);
   }
 
   if (layer->location() == GeomDetEnumerators::endcap ) {   
    
     const ForwardDetLayer * forwardLayer = dynamic_cast<const ForwardDetLayer*>(layer);
     LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeEstimator Endcap r = " << forwardLayer->specificSurface().innerRadius() << " R " << forwardLayer->specificSurface().outerRadius()  <<  " Z " << forwardLayer->specificSurface().position().z() << "\n";  
     
     float zc = fabs(theBCPosition_.z());
     float z =  fabs(forwardLayer->surface().position().z());
     
     float rrange = 15. * theBCPosition_.perp() * (zc - z) / (zc*zc - 15.*zc);
     newEstimator = new ConversionForwardEstimator(-dphi, dphi, rrange);
   }
 
   return newEstimator;
 }

void OutInConversionSeedFinder::makeSeeds ( const edm::Handle< edm::View< reco::CaloCluster > > & allBc ) const

virtual

Implements ConversionSeedFinder.

Definition at line 54 of file OutInConversionSeedFinder.cc.

References PV3DBase< T, PVType, FrameType >::eta(), fillClusterSeeds(), ConversionSeedFinder::findLayers(), i, LogDebug, nSeedsPerBC_, PV3DBase< T, PVType, FrameType >::phi(), ConversionSeedFinder::theBCEnergy_, ConversionSeedFinder::theBCPosition_, ConversionSeedFinder::theSCPosition_, and ConversionSeedFinder::theSeeds_.

Referenced by SoftConversionTrackCandidateProducer::buildCollections(), and ConversionTrackCandidateProducer::buildCollections().

                                                                                                        {
 
   //  std::cout  << "  OutInConversionSeedFinder::makeSeeds() " << "\n";
   theSeeds_.clear();
 
 
   // debug  
   //  std::cout << "  OutInConversionSeedFinder::makeSeeds() SC position " << theSCPosition_.x() << " " << theSCPosition_.y() << " " << theSCPosition_.z() << "\n";
   // std::cout << " SC eta  " <<  theSCPosition_.eta() <<  " SC phi  " <<  theSCPosition_.phi() << "\n";  
   // std::cout << "  OutInConversionSeedFinder::makeSeeds() SC energy " << theSCenergy_  << "\n";  
   //
 
   findLayers();
 
   
   //  std::cout  << " Check Calo cluster collection size " << allBC->size() << "\n";
   
   float  theSCPhi=theSCPosition_.phi();
   float  theSCEta=theSCPosition_.eta();
 
   
 
   //  Loop over the Calo Clusters  in the event looking for seeds 
   reco::CaloClusterCollection::const_iterator bcItr;
   LogDebug("OutInConversionSeedFinder") << "  OutInConversionSeedFinder::makeSeeds() All BC in the event " << "\n";
   for (unsigned i = 0; i < allBC->size(); ++i ) {
   
     //for(bcItr = allBC.begin(); bcItr != allBC.end(); bcItr++) {
     nSeedsPerBC_=0;
 
     theBCEnergy_=allBC->ptrAt(i)->energy();
     theBCPosition_ = GlobalPoint(allBC->ptrAt(i)->position().x(), allBC->ptrAt(i)->position().y(), allBC->ptrAt(i)->position().z() ) ;
     float theBcEta=  theBCPosition_.eta();
     float theBcPhi=  theBCPosition_.phi();
     //    float  dPhi= theBcPhi-theSCPhi;
 
     //    std::cout << "  OutInConversionSeedFinder::makeSeeds() BC eta  " << theBcEta << " phi " <<  theBcPhi << " BC energy " << theBCEnergy_ << " dPhi " << fabs(theBcPhi-theSCPhi) << " dEta " <<  fabs(theBcEta-theSCEta) << "\n";
     
     if ( theBCEnergy_ < 1.5 ) continue;
 
     LogDebug("OutInConversionSeedFinder") << "  OutInConversionSeedFinder::makeSeeds() Passing the >=1.5 GeV cut  BC eta  " << theBcEta << " phi " <<  theBcPhi << " BC energy " << theBCEnergy_ << "\n";
 
     if (  fabs(theBcEta-theSCEta) < 0.015  && fabs(theBcPhi-theSCPhi) < 0.3 ) { 
       LogDebug("OutInConversionSeedFinder") << "  OutInConversionSeedFinder::makeSeeds() in et and phi range passed to the analysis " << "\n";
       fillClusterSeeds( allBC->ptrAt(i)  );
     }
     
 
   }
 
 
   //  std::cout << "Built vector of seeds of size  " << theSeeds_.size() <<  "\n" ;
   
   /*
   int nSeed=0;
   for ( std::vector<TrajectorySeed>::const_iterator iSeed= theSeeds_.begin(); iSeed != theSeeds_.end(); ++iSeed) {
     nSeed++;
     PTrajectoryStateOnDet  ptsod=iSeed->startingState();
     LogDebug("OutInConversionSeedFinder") << nSeed << ")  Direction " << iSeed->direction() << " Num of hits " << iSeed->nHits() <<  " starting state position " << ptsod.parameters().position() << " R " << ptsod.parameters().position().perp() << " phi " << ptsod.parameters().position().phi() << " eta " << ptsod.parameters().position().eta() << "\n" ;
     
     
     DetId tmpId = DetId( iSeed->startingState().detId());
     const GeomDet* tmpDet  = this->getMeasurementTracker()->geomTracker()->idToDet( tmpId );
     GlobalVector gv = tmpDet->surface().toGlobal( iSeed->startingState().parameters().momentum() );
     
     LogDebug("OutInConversionSeedFinder") << "seed perp,phi,eta : " 
                       << gv.perp() << " , " 
                       << gv.phi() << " , " 
                       << gv.eta() << "\n" ; ;
     
 
 
 
     TrajectorySeed::range hitRange = iSeed->recHits();
     for (TrajectorySeed::const_iterator ihit = hitRange.first; ihit != hitRange.second; ihit++) {
    
       if ( ihit->isValid() ) {
 
     LogDebug("OutInConversionSeedFinder") << " Valid hit global position " << this->getMeasurementTracker()->geomTracker()->idToDet((ihit)->geographicalId())->surface().toGlobal((ihit)->localPosition()) << " R " << this->getMeasurementTracker()->geomTracker()->idToDet((ihit)->geographicalId())->surface().toGlobal((ihit)->localPosition()).perp() << " phi " << this->getMeasurementTracker()->geomTracker()->idToDet((ihit)->geographicalId())->surface().toGlobal((ihit)->localPosition()).phi() << " eta " << this->getMeasurementTracker()->geomTracker()->idToDet((ihit)->geographicalId())->surface().toGlobal((ihit)->localPosition()).eta() <<    "\n" ;
 
       }
     }
   } 
   
   */
 
   
   
 }

void OutInConversionSeedFinder::makeSeeds ( const reco::CaloClusterPtr & aBC ) const

virtual

Definition at line 147 of file OutInConversionSeedFinder.cc.

References PV3DBase< T, PVType, FrameType >::eta(), fillClusterSeeds(), ConversionSeedFinder::findLayers(), nSeedsPerBC_, PV3DBase< T, PVType, FrameType >::phi(), ConversionSeedFinder::theBCEnergy_, ConversionSeedFinder::theBCPosition_, ConversionSeedFinder::theSCPosition_, and ConversionSeedFinder::theSeeds_.

                                                                                   {
 
   theSeeds_.clear();
 
   findLayers();
 
   float  theSCPhi=theSCPosition_.phi();
   float  theSCEta=theSCPosition_.eta();
 
   nSeedsPerBC_=0;
 
   theBCEnergy_=aBC->energy();
   theBCPosition_ = GlobalPoint(aBC->position().x(), aBC->position().y(), aBC->position().z() ) ;
   float theBcEta=  theBCPosition_.eta();
   float theBcPhi=  theBCPosition_.phi();
   //  float  dPhi= theBcPhi-theSCPhi;
 
   if ( theBCEnergy_ < 1.5 ) return;
 
   if (  fabs(theBcEta-theSCEta) < 0.015  && fabs(theBcPhi-theSCPhi) < 0.25 ) {
     fillClusterSeeds( aBC);
   }
 
 }

std::pair< FreeTrajectoryState, bool > OutInConversionSeedFinder::makeTrackState ( int charge ) const

private

Definition at line 198 of file OutInConversionSeedFinder.cc.

References PixelRecoUtilities::curvature(), LogDebug, m, PV3DBase< T, PVType, FrameType >::perp(), reco::BeamSpot::position(), csvReporter::r, dttmaxenums::R, query::result, rho, mathSSE::sqrt(), ConversionSeedFinder::theBCEnergy_, ConversionSeedFinder::theBCPosition_, ConversionSeedFinder::theBeamSpot_, and ConversionSeedFinder::theMF_.

Referenced by fillClusterSeeds().

                                                                                               {
 
   std::pair<FreeTrajectoryState,bool> result;
   result.second=false;
  
 
   //std::cout << "  OutInConversionSeedFinder:makeTrackState " << "\n";
 
 
   //  Old GlobalPoint gpOrigine(theBCPosition_.x()*0.3, theBCPosition_.y()*0.3, theBCPosition_.z()*0.3) ;
  //  GlobalPoint gpOrigine(0.,0.,0.);
 
   GlobalPoint  gpOrigine(theBeamSpot_.position().x(),theBeamSpot_.position().y(),theBeamSpot_.position().z()); 
   GlobalVector gvBcRadius = theBCPosition_ - gpOrigine ;
   HepGeom::Point3D<double>  radiusBc(gvBcRadius.x(),gvBcRadius.y(),gvBcRadius.z()) ;
   HepGeom::Point3D<double>  momentumWithoutCurvature = radiusBc.unit() * theBCEnergy_ ;
 
   // compute momentum direction at calo
   double curvature = theMF_->inTesla(theBCPosition_).z() * c_light * 1.e-3 / momentumWithoutCurvature.perp() ;
   curvature /= 100. ; // in cm-1 !!
 
   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState gpOrigine " << gpOrigine.x() << " " <<  gpOrigine.y() << " " <<  gpOrigine.z() << " momentumWithoutCurvature " << momentumWithoutCurvature.mag() << " curvature " << curvature << "\n";
 
   // define rotation angle
   float R = theBCPosition_.perp();
   float r = gpOrigine.perp();
   float rho = 1./curvature;
   // from the formula for the intersection of two circles
   // turns out to be about 2/3 of the deflection of the old formula
   float d = sqrt(r*r+rho*rho);
   float u = rho + rho/d/d*(R*R-rho*rho) - r/d/d*sqrt((R*R-r*r+2*rho*R)*(R*R-r*r+2*rho*R));
   //float u = rho + rho/d/d*(R*R-rho*rho) ;
   if ( u <=R )   result.second=true;
 
   double sinAlpha = 0.5*u/R;
   if ( sinAlpha>(1.-10*DBL_EPSILON) )  sinAlpha = 1.-10*DBL_EPSILON;
   else if ( sinAlpha<-(1.-10*DBL_EPSILON) )  sinAlpha = -(1.-10*DBL_EPSILON);
   
   double newdphi = charge * asin( sinAlpha) ;
 
   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState charge " << charge << " R " << R << " u/R " << u/R << " asin(0.5*u/R) " << asin(sinAlpha) << "\n";
 
   HepGeom::Transform3D rotation =  HepGeom::Rotate3D(newdphi, HepGeom::Vector3D<double> (0., 0. ,1.));
 
 
   HepGeom::Point3D<double>  momentumInTracker = momentumWithoutCurvature.transform(rotation) ;
   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState  R " << R << " r " << r << " rho " << rho  << " d " << d  << " u " << u << " newdphi " << newdphi << " momentumInTracker " <<  momentumInTracker << "\n";
 
   HepGeom::Point3D<double>  hepStartingPoint(gpOrigine.x(), gpOrigine.y(), gpOrigine.z()) ;
 
   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState hepStartingPoint " << hepStartingPoint << "\n";
 
   hepStartingPoint.transform(rotation);
 
   GlobalPoint startingPoint(hepStartingPoint.x(), hepStartingPoint.y(), hepStartingPoint.z());
 
   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState startingPoint " << startingPoint << " calo position " << theBCPosition_ << "\n";
   GlobalVector gvTracker(momentumInTracker.x(), momentumInTracker.y(), momentumInTracker.z());
   GlobalTrajectoryParameters gtp(startingPoint, gvTracker, charge, &(*theMF_) );
   // error matrix
   AlgebraicSymMatrix55 m = AlgebraicMatrixID();
   m(0,0) = 0.1; m(1,1) = 0.1 ; m(2,2) = 0.1 ;
   m(3,3) = 0.1 ; m(4,4) = 0.1;
   
   //std::cout << "OutInConversionSeedFinder::makeTrackState " <<  FreeTrajectoryState(gtp, CurvilinearTrajectoryError(m) ) << "\n";
    
   result.first= FreeTrajectoryState(gtp, CurvilinearTrajectoryError(m) ) ;
   return result;
 
 }

void OutInConversionSeedFinder::startSeed ( const FreeTrajectoryState & fts ) const

private

Definition at line 270 of file OutInConversionSeedFinder.cc.

References alongMomentum, FreeTrajectoryState::charge(), completeSeed(), ConversionSeedFinder::getMeasurementTracker(), i, ConversionSeedFinder::layerList(), LogDebug, makeEstimator(), LayerMeasurements::measurements(), Propagator::propagationDirection(), TrajectoryMeasurement::recHit(), GeometricSearchDet::surface(), theFirstMeasurements_, ConversionSeedFinder::thePropagatorAlongMomentum_, ConversionSeedFinder::thePropagatorOppositeToMomentum_, and ConversionSeedFinder::trackStateFromClusters().

Referenced by fillClusterSeeds().

                                                                                {
 
 
   //  std::cout << "OutInConversionSeedFinder::startSeed layer list " << this->layerList().size() <<  "\n";
   //std::cout << "OutInConversionSeedFinder::startSeed  fts " << fts <<  "\n";  
 
   std::vector<const DetLayer*> myLayers=layerList();
   if ( myLayers.size() > 3 ) {
     
     for(unsigned int ilayer = myLayers.size()-1; ilayer >= myLayers.size()-2; --ilayer) {
       const DetLayer * layer = myLayers[ilayer];
       
       
       // allow the z of the hit to be within a straight line from a vertex
       // of +-15 cm to the cluster
       //      float dphi = 0.015;
       float dphi = 0.030;
 
       MeasurementEstimator * newEstimator = makeEstimator(layer, dphi);
 
      
       //std::cout << "OutInSeedFinder::startSeed propagationDirection  " << int(thePropagatorAlongMomentum_->propagationDirection() ) << "\n";       
       
       TSOS tsos(fts, layer->surface() );
       
       LogDebug("OutInConversionSeedFinder") << "OutInSeedFinder::startSeed  after  TSOS tsos(fts, layer->surface() ) " << "\n";
       
       LayerMeasurements theLayerMeasurements_(this->getMeasurementTracker() );
       theFirstMeasurements_ = theLayerMeasurements_.measurements( *layer, tsos, *thePropagatorAlongMomentum_, *newEstimator);
       
       //std::cout << "OutInSeedFinder::startSeed  after  theFirstMeasurements_   " << theFirstMeasurements_.size() <<  "\n";
       
       if(theFirstMeasurements_.size() > 1) // always a dummy returned, too
     LogDebug("OutInConversionSeedFinder") <<  " Found " << theFirstMeasurements_.size()-1 << " 1st hits in seed" << "\n";
       
       delete newEstimator;
       
       LogDebug("OutInConversionSeedFinder") << "OutInSeedFinder::startSeed Layer " << ilayer << " theFirstMeasurements_.size " << theFirstMeasurements_.size() << "\n";
       
       for(unsigned int i = 0; i < theFirstMeasurements_.size(); ++i) {
     TrajectoryMeasurement m1 = theFirstMeasurements_[i];
     if(m1.recHit()->isValid()) {
       
       // update the fts to start from this point.  much better than starting from
       // extrapolated point along the line
       GlobalPoint hitPoint = m1.recHit()->globalPosition();
       LogDebug("OutInConversionSeedFinder")  << " Valid hit at R  " <<   m1.recHit()->globalPosition().perp() << " Z " <<  m1.recHit()->globalPosition().z() << " eta " << m1.recHit()->globalPosition().eta() << " phi " << m1.recHit()->globalPosition().phi()  << " xyz " <<  m1.recHit()->globalPosition() << "\n";
 
       
       FreeTrajectoryState newfts = trackStateFromClusters(fts.charge(), hitPoint, alongMomentum, 0.8);
       LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::startSeed  newfts " << newfts << "\n";
       LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::startSeed propagationDirection  after switching " << int(thePropagatorOppositeToMomentum_->propagationDirection() ) << "\n";        
       //  std::cout << "OutInConversionSeedFinder::startSeed propagationDirection  after switching " << int(thePropagatorOppositeToMomentum_->propagationDirection() ) << "\n";        
 
        
       completeSeed(m1, newfts, thePropagatorOppositeToMomentum_, ilayer-1);
       // skip a layer, if you haven't already skipped the first layer
       if(ilayer == myLayers.size()-1) {
         completeSeed(m1, newfts, thePropagatorOppositeToMomentum_, ilayer-2);
       }
     }
       }
       
     } // loop over layers
   }
 
 
 
   
 }