EgammaPhotonAlgos/src/OutInConversionSeedFinder.cc

0001 #include "RecoEgamma/EgammaPhotonAlgos/interface/OutInConversionSeedFinder.h"
0002 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionBarrelEstimator.h"
0003 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionForwardEstimator.h"
0004 #include "DataFormats/Math/interface/deltaPhi.h"
0005 //
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 // Field
0008 #include "MagneticField/Engine/interface/MagneticField.h"
0009 //
0010 #include "DataFormats/CLHEP/interface/AlgebraicObjects.h"
0011 // Geometry
0012 //
0013 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0014 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0015 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0016 #include "TrackingTools/MeasurementDet/interface/LayerMeasurements.h"
0017 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0018 //
0019
0020 //
0021 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0022 #include "CLHEP/Geometry/Point3D.h"
0023 #include "CLHEP/Geometry/Vector3D.h"
0024 #include "CLHEP/Geometry/Transform3D.h"
0025 #include <cfloat>
0026
0027 namespace {
0028   inline double ptFast(const double energy, const math::XYZPoint& position, const math::XYZPoint& origin) {
0029     const auto v = position - origin;
0030     return energy * std::sqrt(v.perp2() / v.mag2());
0031   }
0032 }  // namespace
0033
0034 OutInConversionSeedFinder::OutInConversionSeedFinder(const edm::ParameterSet& conf, edm::ConsumesCollector&& iC)
0035     : ConversionSeedFinder(conf, iC), conf_(conf) {
0036   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder CTOR "
0037                                         << "\n";
0038
0039   maxNumberOfOutInSeedsPerBC_ = conf_.getParameter<int>("maxNumOfSeedsOutIn");
0040   bcEtcut_ = conf_.getParameter<double>("bcEtCut");
0041   bcEcut_ = conf_.getParameter<double>("bcECut");
0042   useEtCut_ = conf_.getParameter<bool>("useEtCut");
0043   //the2ndHitdphi_ = 0.01;
0044   the2ndHitdphi_ = 0.03;
0045   the2ndHitdzConst_ = 5.;
0046   the2ndHitdznSigma_ = 2.;
0047 }
0048
0049 OutInConversionSeedFinder::~OutInConversionSeedFinder() {
0050   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder DTOR "
0051                                         << "\n";
0052 }
0053
0054 void OutInConversionSeedFinder::makeSeeds(const edm::Handle<edm::View<reco::CaloCluster> >& allBC) {
0055   theSeeds_.clear();
0056
0057   //  std::cout  << "  OutInConversionSeedFinder::makeSeeds() " << "\n";
0058
0059   // debug
0060   //  std::cout << "  OutInConversionSeedFinder::makeSeeds() SC position " << theSCPosition_.x() << " " << theSCPosition_.y() << " " << theSCPosition_.z() << "\n";
0061   // std::cout << " SC eta  " <<  theSCPosition_.eta() <<  " SC phi  " <<  theSCPosition_.phi() << "\n";
0062   // std::cout << "  OutInConversionSeedFinder::makeSeeds() SC energy " << theSCenergy_  << "\n";
0063   //
0064
0065   findLayers();
0066
0067   //  std::cout  << " Check Calo cluster collection size " << allBC->size() << "\n";
0068
0069   float theSCPhi = theSCPosition_.phi();
0070   float theSCEta = theSCPosition_.eta();
0071
0072   //  Loop over the Calo Clusters  in the event looking for seeds
0073   LogDebug("OutInConversionSeedFinder") << "  OutInConversionSeedFinder::makeSeeds() All BC in the event "
0074                                         << "\n";
0075   for (unsigned i = 0; i < allBC->size(); ++i) {
0076     nSeedsPerBC_ = 0;
0077
0078     const reco::CaloCluster& theBC = allBC->at(i);
0079     const math::XYZPoint& rawBCpos = theBC.position();
0080
0081     theBCPosition_ = GlobalPoint(rawBCpos.x(), rawBCpos.y(), rawBCpos.z());
0082     float theBcEta = theBCPosition_.eta();
0083     float theBcPhi = theBCPosition_.phi();
0084     //    float  dPhi= theBcPhi-theSCPhi;
0085     theBCEnergy_ = theBC.energy();
0086
0087     float EtOrECut = bcEcut_;
0088     if (useEtCut_) {
0089       theBCEnergy_ = ptFast(theBCEnergy_, rawBCpos, math::XYZPoint(0, 0, 0));
0090       EtOrECut = bcEtcut_;
0091     }
0092
0093     if (theBCEnergy_ < EtOrECut)
0094       continue;
0095     // std::cout << "  OutInConversionSeedFinder::makeSeeds() BC eta  " << theBcEta << " phi " <<  theBcPhi << " BC transverse energy " << theBCEnergy_ << " dPhi " << fabs(theBcPhi-theSCPhi) << " dEta " <<  fabs(theBcEta-theSCEta) << "\n";
0096
0097     LogDebug("OutInConversionSeedFinder")
0098         << "  OutInConversionSeedFinder::makeSeeds() Passing the >=1.5 GeV cut  BC eta  " << theBcEta << " phi "
0099         << theBcPhi << " BC energy " << theBCEnergy_ << "\n";
0100
0101     if (fabs(theBcEta - theSCEta) < 0.015 && reco::deltaPhi(theBcPhi, theSCPhi) < 0.3) {
0102       LogDebug("OutInConversionSeedFinder")
0103           << "  OutInConversionSeedFinder::makeSeeds() in et and phi range passed to the analysis "
0104           << "\n";
0105       fillClusterSeeds(allBC->ptrAt(i));
0106     }
0107   }
0108
0109   //  std::cout << "Built vector of seeds of size  " << theSeeds_.size() <<  "\n" ;
0110
0111   ///// This part is only for local debugging: will be trhown away when no longer needed
0112   /*
0113   int nSeed=0;
0114   for ( std::vector<TrajectorySeed>::const_iterator iSeed= theSeeds_.begin(); iSeed != theSeeds_.end(); ++iSeed) {
0115     nSeed++;
0116     PTrajectoryStateOnDet  ptsod=iSeed->startingState();
0117     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" ;
0118
0119
0120     DetId tmpId = DetId( iSeed->startingState().detId());
0121     const GeomDet* tmpDet  = this->getMeasurementTracker()->geomTracker()->idToDet( tmpId );
0122     GlobalVector gv = tmpDet->surface().toGlobal( iSeed->startingState().parameters().momentum() );
0123
0124     LogDebug("OutInConversionSeedFinder") << "seed perp,phi,eta : "
0125                       << gv.perp() << " , "
0126                       << gv.phi() << " , "
0127                       << gv.eta() << "\n" ; ;
0128
0129
0130
0131
0132     TrajectorySeed::range hitRange = iSeed->recHits();
0133     for (TrajectorySeed::const_iterator ihit = hitRange.first; ihit != hitRange.second; ihit++) {
0134
0135       if ( ihit->isValid() ) {
0136
0137     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" ;
0138
0139       }
0140     }
0141   }
0142
0143   */
0144 }
0145
0146 void OutInConversionSeedFinder::makeSeeds(const reco::CaloClusterPtr& aBC) {
0147   theSeeds_.clear();
0148
0149   findLayers();
0150
0151   float theSCPhi = theSCPosition_.phi();
0152   float theSCEta = theSCPosition_.eta();
0153
0154   nSeedsPerBC_ = 0;
0155
0156   // theBCEnergy_=aBC->energy();
0157   theBCEnergy_ = ptFast(aBC->energy(), aBC->position(), math::XYZPoint(0, 0, 0));
0158   theBCPosition_ = GlobalPoint(aBC->position().x(), aBC->position().y(), aBC->position().z());
0159   float theBcEta = theBCPosition_.eta();
0160   float theBcPhi = theBCPosition_.phi();
0161   //  float  dPhi= theBcPhi-theSCPhi;
0162
0163   if (theBCEnergy_ < bcEtcut_)
0164     return;
0165
0166   if (fabs(theBcEta - theSCEta) < 0.015 && fabs(theBcPhi - theSCPhi) < 0.25) {
0167     fillClusterSeeds(aBC);
0168   }
0169 }
0170
0171 void OutInConversionSeedFinder::fillClusterSeeds(const reco::CaloClusterPtr& bc) {
0172   theFirstMeasurements_.clear();
0173
0174   /// negative charge state
0175   auto const stateNeg = makeTrackState(-1);
0176   if (stateNeg.second) {
0177     startSeed(stateNeg.first);
0178   }
0179   /// positive charge state
0180
0181   auto const statePos = makeTrackState(1);
0182   if (statePos.second) {
0183     startSeed(statePos.first);
0184   }
0185   theFirstMeasurements_.clear();
0186 }
0187
0188 std::pair<FreeTrajectoryState, bool> OutInConversionSeedFinder::makeTrackState(int charge) const {
0189   std::pair<FreeTrajectoryState, bool> result;
0190   result.second = false;
0191
0192   //std::cout << "  OutInConversionSeedFinder:makeTrackState " << "\n";
0193
0194   //  Old GlobalPoint gpOrigine(theBCPosition_.x()*0.3, theBCPosition_.y()*0.3, theBCPosition_.z()*0.3) ;
0195   //  GlobalPoint gpOrigine(0.,0.,0.);
0196
0197   GlobalPoint gpOrigine(theBeamSpot_.position().x(), theBeamSpot_.position().y(), theBeamSpot_.position().z());
0198   GlobalVector gvBcRadius = theBCPosition_ - gpOrigine;
0199   HepGeom::Point3D<double> radiusBc(gvBcRadius.x(), gvBcRadius.y(), gvBcRadius.z());
0200   HepGeom::Point3D<double> momentumWithoutCurvature = radiusBc.unit() * theBCEnergy_;
0201
0202   // compute momentum direction at calo
0203   double curvature = theMF_->inTesla(theBCPosition_).z() * c_light * 1.e-3 / momentumWithoutCurvature.perp();
0204   curvature /= 100.;  // in cm-1 !!
0205   if (curvature == 0)
0206     return result;
0207
0208   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState gpOrigine " << gpOrigine.x()
0209                                         << " " << gpOrigine.y() << " " << gpOrigine.z() << " momentumWithoutCurvature "
0210                                         << momentumWithoutCurvature.mag() << " curvature " << curvature << "\n";
0211
0212   // define rotation angle
0213   float R = theBCPosition_.perp();
0214   float r = gpOrigine.perp();
0215   float rho = 1. / curvature;
0216   // from the formula for the intersection of two circles
0217   // turns out to be about 2/3 of the deflection of the old formula
0218   float d = sqrt(r * r + rho * rho);
0219   float u = rho + rho / d / d * (R * R - rho * rho) -
0220             r / d / d * sqrt((R * R - r * r + 2 * rho * R) * (R * R - r * r + 2 * rho * R));
0221   //float u = rho + rho/d/d*(R*R-rho*rho) ;
0222   if (u <= R)
0223     result.second = true;
0224   else
0225     return result;
0226
0227   double sinAlpha = 0.5 * u / R;
0228   if (sinAlpha > (1. - 10 * DBL_EPSILON))
0229     sinAlpha = 1. - 10 * DBL_EPSILON;
0230   else if (sinAlpha < -(1. - 10 * DBL_EPSILON))
0231     sinAlpha = -(1. - 10 * DBL_EPSILON);
0232
0233   double newdphi = charge * asin(sinAlpha);
0234
0235   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState charge " << charge << " R " << R
0236                                         << " u/R " << u / R << " asin(0.5*u/R) " << asin(sinAlpha) << "\n";
0237
0238   HepGeom::Transform3D rotation = HepGeom::Rotate3D(newdphi, HepGeom::Vector3D<double>(0., 0., 1.));
0239
0240   HepGeom::Point3D<double> momentumInTracker = momentumWithoutCurvature.transform(rotation);
0241   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState  R " << R << " r " << r << " rho "
0242                                         << rho << " d " << d << " u " << u << " newdphi " << newdphi
0243                                         << " momentumInTracker " << momentumInTracker << "\n";
0244
0245   HepGeom::Point3D<double> hepStartingPoint(gpOrigine.x(), gpOrigine.y(), gpOrigine.z());
0246
0247   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState hepStartingPoint "
0248                                         << hepStartingPoint << "\n";
0249
0250   hepStartingPoint.transform(rotation);
0251
0252   GlobalPoint startingPoint(hepStartingPoint.x(), hepStartingPoint.y(), hepStartingPoint.z());
0253
0254   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState startingPoint " << startingPoint
0255                                         << " calo position " << theBCPosition_ << "\n";
0256   GlobalVector gvTracker(momentumInTracker.x(), momentumInTracker.y(), momentumInTracker.z());
0257   GlobalTrajectoryParameters gtp(startingPoint, gvTracker, charge, &(*theMF_));
0258   // error matrix
0259   AlgebraicSymMatrix55 m = AlgebraicMatrixID();
0260   m(0, 0) = 0.1;
0261   m(1, 1) = 0.1;
0262   m(2, 2) = 0.1;
0263   m(3, 3) = 0.1;
0264   m(4, 4) = 0.1;
0265
0266   //  std::cout << "OutInConversionSeedFinder::makeTrackState " <<  FreeTrajectoryState(gtp, CurvilinearTrajectoryError(m) ) << "\n";
0267
0268   result.first = FreeTrajectoryState(gtp, CurvilinearTrajectoryError(m));
0269   return result;
0270 }
0271
0272 void OutInConversionSeedFinder::startSeed(const FreeTrajectoryState& fts) {
0273   //  std::cout << "OutInConversionSeedFinder::startSeed layer list " << this->layerList().size() <<  "\n";
0274   //std::cout << "OutInConversionSeedFinder::startSeed  fts " << fts <<  "\n";
0275
0276   std::vector<const DetLayer*> myLayers = layerList();
0277   if (myLayers.size() > 3) {
0278     for (unsigned int ilayer = myLayers.size() - 1; ilayer >= myLayers.size() - 2; --ilayer) {
0279       const DetLayer* layer = myLayers[ilayer];
0280
0281       // allow the z of the hit to be within a straight line from a vertex
0282       // of +-15 cm to the cluster
0283       //      float dphi = 0.015;
0284       float dphi = 0.030;
0285
0286       MeasurementEstimator* newEstimator = makeEstimator(layer, dphi);
0287
0288       //std::cout << "OutInSeedFinder::startSeed propagationDirection  " << int(thePropagatorAlongMomentum_->propagationDirection() ) << "\n";
0289
0290       TSOS tsos(fts, layer->surface());
0291
0292       LogDebug("OutInConversionSeedFinder") << "OutInSeedFinder::startSeed  after  TSOS tsos(fts, layer->surface() ) "
0293                                             << "\n";
0294
0295       LayerMeasurements theLayerMeasurements_(*this->getMeasurementTracker(), *theTrackerData_);
0296       theFirstMeasurements_ =
0297           theLayerMeasurements_.measurements(*layer, tsos, *thePropagatorAlongMomentum_, *newEstimator);
0298
0299       //std::cout << "OutInSeedFinder::startSeed  after  theFirstMeasurements_   " << theFirstMeasurements_.size() <<  "\n";
0300
0301       if (theFirstMeasurements_.size() > 1)  // always a dummy returned, too
0302         LogDebug("OutInConversionSeedFinder") << " Found " << theFirstMeasurements_.size() - 1 << " 1st hits in seed"
0303                                               << "\n";
0304
0305       delete newEstimator;
0306
0307       LogDebug("OutInConversionSeedFinder") << "OutInSeedFinder::startSeed Layer " << ilayer
0308                                             << " theFirstMeasurements_.size " << theFirstMeasurements_.size() << "\n";
0309
0310       for (unsigned int i = 0; i < theFirstMeasurements_.size(); ++i) {
0311         TrajectoryMeasurement m1 = theFirstMeasurements_[i];
0312         if (m1.recHit()->isValid()) {
0313           // update the fts to start from this point.  much better than starting from
0314           // extrapolated point along the line
0315           GlobalPoint hitPoint = m1.recHit()->globalPosition();
0316           LogDebug("OutInConversionSeedFinder")
0317               << " Valid hit at R  " << m1.recHit()->globalPosition().perp() << " Z "
0318               << m1.recHit()->globalPosition().z() << " eta " << m1.recHit()->globalPosition().eta() << " phi "
0319               << m1.recHit()->globalPosition().phi() << " xyz " << m1.recHit()->globalPosition() << "\n";
0320
0321           FreeTrajectoryState newfts = trackStateFromClusters(fts.charge(), hitPoint, alongMomentum, 0.8);
0322           //std::cout << "OutInConversionSeedFinder::startSeed  newfts " << newfts << "\n";
0323           LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::startSeed  newfts " << newfts << "\n";
0324           LogDebug("OutInConversionSeedFinder")
0325               << "OutInConversionSeedFinder::startSeed propagationDirection  after switching "
0326               << int(thePropagatorOppositeToMomentum_->propagationDirection()) << "\n";
0327           //  std::cout << "OutInConversionSeedFinder::startSeed propagationDirection  after switching " << int(thePropagatorOppositeToMomentum_->propagationDirection() ) << "\n";
0328
0329           completeSeed(m1, newfts, thePropagatorOppositeToMomentum_, ilayer - 1);
0330           // skip a layer, if you haven't already skipped the first layer
0331           if (ilayer == myLayers.size() - 1) {
0332             completeSeed(m1, newfts, thePropagatorOppositeToMomentum_, ilayer - 2);
0333           }
0334         }
0335       }
0336
0337     }  // loop over layers
0338   }
0339 }
0340
0341 MeasurementEstimator* OutInConversionSeedFinder::makeEstimator(const DetLayer* layer, float dphi) const {
0342   //std::cout  << "OutInConversionSeedFinder::makeEstimator  " << "\n";
0343
0344   MeasurementEstimator* newEstimator = nullptr;
0345
0346   if (layer->location() == GeomDetEnumerators::barrel) {
0347     const BarrelDetLayer* barrelLayer = dynamic_cast<const BarrelDetLayer*>(layer);
0348     LogDebug("OutInConversionSeedFinder")
0349         << "OutInConversionSeedFinder::makeEstimator Barrel  r = " << barrelLayer->specificSurface().radius() << " "
0350         << "\n";
0351     float r = barrelLayer->specificSurface().radius();
0352     float zrange = 15. * (1. - r / theBCPosition_.perp());
0353     newEstimator = new ConversionBarrelEstimator(-dphi, dphi, -zrange, zrange);
0354   }
0355
0356   if (layer->location() == GeomDetEnumerators::endcap) {
0357     const ForwardDetLayer* forwardLayer = dynamic_cast<const ForwardDetLayer*>(layer);
0358     LogDebug("OutInConversionSeedFinder")
0359         << "OutInConversionSeedFinder::makeEstimator Endcap r = " << forwardLayer->specificSurface().innerRadius()
0360         << " R " << forwardLayer->specificSurface().outerRadius() << " Z "
0361         << forwardLayer->specificSurface().position().z() << "\n";
0362
0363     float zc = fabs(theBCPosition_.z());
0364     float z = fabs(forwardLayer->surface().position().z());
0365
0366     float rrange = 15. * theBCPosition_.perp() * (zc - z) / (zc * zc - 15. * zc);
0367     newEstimator = new ConversionForwardEstimator(-dphi, dphi, rrange);
0368   }
0369
0370   return newEstimator;
0371 }
0372
0373 void OutInConversionSeedFinder::completeSeed(const TrajectoryMeasurement& m1,
0374                                              const FreeTrajectoryState& fts,
0375                                              const Propagator* propagator,
0376                                              int ilayer) {
0377   //std::cout <<  "OutInConversionSeedFinder::completeSeed ilayer " << ilayer << "\n";
0378
0379   MeasurementEstimator* newEstimator = nullptr;
0380   const DetLayer* layer = theLayerList_[ilayer];
0381
0382   if (layer->location() == GeomDetEnumerators::barrel) {
0383     // z error for 2nd hit is  2 sigma quadded with 5 cm
0384     LogDebug("OutInConversionSeedFinder") << " Barrel OutInConversionSeedFinder::completeSeed " << the2ndHitdznSigma_
0385                                           << " " << the2ndHitdzConst_ << " " << the2ndHitdphi_ << "\n";
0386     float dz = sqrt(the2ndHitdznSigma_ * the2ndHitdznSigma_ * m1.recHit()->globalPositionError().czz() +
0387                     the2ndHitdzConst_ * the2ndHitdzConst_);
0388     newEstimator = new ConversionBarrelEstimator(-the2ndHitdphi_, the2ndHitdphi_, -dz, dz);
0389   } else {
0390     LogDebug("OutInConversionSeedFinder") << " EndCap OutInConversionSeedFinder::completeSeed " << the2ndHitdznSigma_
0391                                           << " " << the2ndHitdzConst_ << " " << the2ndHitdphi_ << "\n";
0392     // z error for 2nd hit is 2sigma quadded with 5 cm
0393     //float m1dr = m1.recHit().globalPositionError().rerr(m1.recHit().globalPosition());
0394     float m1dr = sqrt(m1.recHit()->localPositionError().yy());
0395     float dr = sqrt(the2ndHitdznSigma_ * the2ndHitdznSigma_ * m1dr * m1dr + the2ndHitdzConst_ * the2ndHitdznSigma_);
0396     // LogDebug("OutInConversionSeedFinder") << "second hit forward dr " << dr << " this hit " << m1dr << endl;
0397     newEstimator = new ConversionForwardEstimator(-the2ndHitdphi_, the2ndHitdphi_, dr);
0398   }
0399
0400   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::completeSeed  ilayer " << ilayer << "\n";
0401
0402   // Get the measurements consistent with the FTS and the Estimator
0403   TSOS tsos(fts, layer->surface());
0404   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::completeSeed propagationDirection  "
0405                                         << int(propagator->propagationDirection()) << "\n";
0406   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::completeSeed pointer to estimator "
0407                                         << newEstimator << "\n";
0408
0409   LayerMeasurements theLayerMeasurements_(*this->getMeasurementTracker(), *theTrackerData_);
0410   std::vector<TrajectoryMeasurement> measurements =
0411       theLayerMeasurements_.measurements(*layer, tsos, *propagator, *newEstimator);
0412   //std::cout << "OutInConversionSeedFinder::completeSeed Found " << measurements.size() << " second hits " << "\n";
0413   delete newEstimator;
0414
0415   for (unsigned int i = 0; i < measurements.size(); ++i) {
0416     if (measurements[i].recHit()->isValid()) {
0417       createSeed(m1, measurements[i]);
0418     }
0419   }
0420
0421   //LogDebug("OutInConversionSeedFinder") << "COMPLETED " << theSeeds_.size() << " SEEDS " << "\n";
0422 }
0423
0424 void OutInConversionSeedFinder::createSeed(const TrajectoryMeasurement& m1, const TrajectoryMeasurement& m2) {
0425   //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";
0426
0427   FreeTrajectoryState fts = createSeedFTS(m1, m2);
0428
0429   //std::cout << "OutInConversionSeedFinder::createSeed First point errors " <<m1.recHit()->parametersError() << "\n";
0430   // std::cout << "original cluster FTS " << fts <<"\n";
0431   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::createSeed First point errors "
0432                                         << m1.recHit()->parametersError() << "\n";
0433   LogDebug("OutInConversionSeedFinder") << "original cluster FTS " << fts << "\n";
0434
0435   //std::cout  << "OutInConversionSeedFinder::createSeed propagation dir " << int( thePropagatorOppositeToMomentum_->propagationDirection() ) << "\n";
0436   TrajectoryStateOnSurface state1 = thePropagatorOppositeToMomentum_->propagate(fts, m1.recHit()->det()->surface());
0437
0438   // LogDebug("OutInConversionSeedFinder") << "hit surface " << h1.det().surface().position() << endl;
0439   // LogDebug("OutInConversionSeedFinder") << "prop to " << typeid(h1.det().surface()).name() << endl;
0440   // LogDebug("OutInConversionSeedFinder") << "prop to first hit " << state1 << endl;
0441   // LogDebug("OutInConversionSeedFinder") << "update to " << h1.globalPosition() << endl;
0442
0443   if (state1.isValid()) {
0444     TrajectoryStateOnSurface updatedState1 = theUpdator_.update(state1, *m1.recHit());
0445
0446     if (updatedState1.isValid()) {
0447       TrajectoryStateOnSurface state2 = thePropagatorOppositeToMomentum_->propagate(
0448           *updatedState1.freeTrajectoryState(), m2.recHit()->det()->surface());
0449
0450       if (state2.isValid()) {
0451         TrajectoryStateOnSurface updatedState2 = theUpdator_.update(state2, *m2.recHit());
0452         TrajectoryMeasurement meas1(state1, updatedState1, m1.recHit(), m1.estimate(), m1.layer());
0453         TrajectoryMeasurement meas2(state2, updatedState2, m2.recHit(), m2.estimate(), m2.layer());
0454
0455         edm::OwnVector<TrackingRecHit> myHits;
0456         myHits.push_back(meas1.recHit()->hit()->clone());
0457         myHits.push_back(meas2.recHit()->hit()->clone());
0458
0459         if (nSeedsPerBC_ >= maxNumberOfOutInSeedsPerBC_)
0460           return;
0461
0462         PTrajectoryStateOnDet ptsod =
0463             trajectoryStateTransform::persistentState(state2, meas2.recHit()->hit()->geographicalId().rawId());
0464
0465         LogDebug("OutInConversionSeedFinder")
0466             << "OutInConversionSeedFinder::createSeed new seed  from state " << state2.globalPosition() << "\n";
0467         LogDebug("OutInConversionSeedFinder")
0468             << "OutInConversionSeedFinder::createSeed new seed  ptsod " << ptsod.parameters().position() << " R "
0469             << ptsod.parameters().position().perp() << " phi " << ptsod.parameters().position().phi() << " eta "
0470             << ptsod.parameters().position().eta() << "\n";
0471
0472         theSeeds_.push_back(TrajectorySeed(ptsod, myHits, oppositeToMomentum));
0473         nSeedsPerBC_++;
0474       }
0475     }
0476   }
0477 }
0478
0479 FreeTrajectoryState OutInConversionSeedFinder::createSeedFTS(const TrajectoryMeasurement& m1,
0480                                                              const TrajectoryMeasurement& m2) const {
0481   GlobalPoint xmeas = fixPointRadius(m1);
0482   GlobalPoint xvert = fixPointRadius(m2);
0483
0484   float pt = theSCenergy_ * sin(theSCPosition_.theta());
0485   float pz = theSCenergy_ * cos(theSCPosition_.theta());
0486
0487   // doesn't work at all for endcap, where r is badly measured
0488   //float dphidr = (p1.phi()-p2.phi())/(p1.perp()-p2.perp());
0489   //int charge = (dphidr > 0.) ? -1 : 1;
0490   int charge = m2.predictedState().charge();
0491
0492   double BInTesla = theMF_->inTesla(xmeas).z();
0493   GlobalVector xdiff = xmeas - xvert;
0494
0495   double phi = xdiff.phi();
0496   double pxOld = pt * cos(phi);
0497   double pyOld = pt * sin(phi);
0498   double RadCurv = 100 * pt / (BInTesla * 0.29979);
0499   double alpha = asin(0.5 * xdiff.perp() / RadCurv);
0500   float ca = cos(charge * alpha);
0501   float sa = sin(charge * alpha);
0502   double pxNew = ca * pxOld + sa * pyOld;
0503   double pyNew = -sa * pxOld + ca * pyOld;
0504   GlobalVector pNew(pxNew, pyNew, pz);
0505
0506   GlobalTrajectoryParameters gp(xmeas, pNew, charge, &(*theMF_));
0507
0508   AlgebraicSymMatrix55 m = AlgebraicMatrixID();
0509   m(0, 0) = 0.05;
0510   m(1, 1) = 0.02;
0511   m(2, 2) = 0.007;
0512   m(3, 3) = 10.;
0513   m(4, 4) = 10.;
0514   //std::cout  << "OutInConversionSeedFinder::createSeedFTS " <<  FreeTrajectoryState(gp, CurvilinearTrajectoryError(m))  << "\n";
0515   return FreeTrajectoryState(gp, CurvilinearTrajectoryError(m));
0516 }
0517
0518 GlobalPoint OutInConversionSeedFinder::fixPointRadius(const TrajectoryMeasurement& m1) const {
0519   GlobalPoint p1 = m1.recHit()->globalPosition();
0520   GlobalPoint p2;
0521   if (m1.layer()->location() == GeomDetEnumerators::barrel) {
0522     p2 = p1;
0523   } else {
0524     float z = p1.z();
0525     float phi = p1.phi();
0526     float theta = theSCPosition_.theta();
0527     float r = p1.z() * tan(theta);
0528     p2 = GlobalPoint(r * cos(phi), r * sin(phi), z);
0529     // LogDebug("OutInConversionSeedFinder") << "fixing point radius " << p2 << " from " << p1 << endl;
0530   }
0531   return p2;
0532 }