|
|
|
#include <OutInConversionSeedFinder.h>
Definition at line 34 of file OutInConversionSeedFinder.h.
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.
| OutInConversionSeedFinder::OutInConversionSeedFinder | ( | const edm::ParameterSet & | config | ) |
Definition at line 30 of file OutInConversionSeedFinder.cc.
References conf_, edm::ParameterSet::getParameter(), LogDebug, maxNumberOfOutInSeedsPerBC_, the2ndHitdphi_, the2ndHitdzConst_, and the2ndHitdznSigma_.
: ConversionSeedFinder( conf ), conf_(conf) { LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder CTOR " << "\n"; maxNumberOfOutInSeedsPerBC_ = conf_.getParameter<int>("maxNumOfSeedsOutIn"); //the2ndHitdphi_ = 0.01; the2ndHitdphi_ = 0.03; the2ndHitdzConst_ = 5.; the2ndHitdznSigma_ = 2.; }
| OutInConversionSeedFinder::~OutInConversionSeedFinder | ( | ) | [virtual] |
Definition at line 48 of file OutInConversionSeedFinder.cc.
References LogDebug.
{
LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder DTOR " << "\n";
}
| void OutInConversionSeedFinder::completeSeed | ( | const TrajectoryMeasurement & | m1, |
| FreeTrajectoryState & | fts, | ||
| const Propagator * | propagator, | ||
| int | layer | ||
| ) | const [private] |
Definition at line 382 of file OutInConversionSeedFinder.cc.
References Reference_intrackfit_cff::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 437 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 505 of file OutInConversionSeedFinder.cc.
References alpha, TrajectoryStateOnSurface::charge(), DeDxDiscriminatorTools::charge(), funct::cos(), fixPointRadius(), m, PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::phi(), phi, TrajectoryMeasurement::predictedState(), 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 177 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);
}
theFirstMeasurements_.clear();
}
| GlobalPoint OutInConversionSeedFinder::fixPointRadius | ( | const TrajectoryMeasurement & | m1 | ) | const [private] |
Definition at line 551 of file OutInConversionSeedFinder.cc.
References Reference_intrackfit_cff::barrel, funct::cos(), TrajectoryMeasurement::layer(), DetLayer::location(), p1, p2, PV3DBase< T, PVType, FrameType >::phi(), phi, csvReporter::r, TrajectoryMeasurement::recHit(), funct::sin(), funct::tan(), ConversionSeedFinder::theSCPosition_, theta(), PV3DBase< T, PVType, FrameType >::theta(), PV3DBase< T, PVType, FrameType >::z(), and 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 344 of file OutInConversionSeedFinder.cc.
References Reference_intrackfit_cff::barrel, Reference_intrackfit_cff::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_, PV3DBase< T, PVType, FrameType >::z(), and 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 55 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 ConversionTrackCandidateProducer::buildCollections().
{
theSeeds_.clear();
// std::cout << " OutInConversionSeedFinder::makeSeeds() " << "\n";
// 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 148 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 199 of file OutInConversionSeedFinder.cc.
References PixelRecoUtilities::curvature(), LogDebug, m, PV3DBase< T, PVType, FrameType >::perp(), reco::BeamSpot::position(), dttmaxenums::R, csvReporter::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 271 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
}
}
Reimplemented from ConversionSeedFinder.
Definition at line 57 of file OutInConversionSeedFinder.h.
Referenced by OutInConversionSeedFinder().
int OutInConversionSeedFinder::maxNumberOfOutInSeedsPerBC_ [private] |
Definition at line 83 of file OutInConversionSeedFinder.h.
Referenced by createSeed(), and OutInConversionSeedFinder().
int OutInConversionSeedFinder::nSeedsPerBC_ [mutable, private] |
Definition at line 82 of file OutInConversionSeedFinder.h.
Referenced by createSeed(), and makeSeeds().
float OutInConversionSeedFinder::the2ndHitdphi_ [private] |
Definition at line 78 of file OutInConversionSeedFinder.h.
Referenced by completeSeed(), and OutInConversionSeedFinder().
float OutInConversionSeedFinder::the2ndHitdzConst_ [private] |
Definition at line 79 of file OutInConversionSeedFinder.h.
Referenced by completeSeed(), and OutInConversionSeedFinder().
float OutInConversionSeedFinder::the2ndHitdznSigma_ [private] |
Definition at line 80 of file OutInConversionSeedFinder.h.
Referenced by completeSeed(), and OutInConversionSeedFinder().
std::vector<TrajectoryMeasurement> OutInConversionSeedFinder::theFirstMeasurements_ [mutable, private] |
Definition at line 81 of file OutInConversionSeedFinder.h.
Referenced by fillClusterSeeds(), and startSeed().
1.7.3