CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/RecoEgamma/EgammaPhotonAlgos/src/OutInConversionSeedFinder.cc

Go to the documentation of this file.
00001 #include "RecoEgamma/EgammaPhotonAlgos/interface/OutInConversionSeedFinder.h"
00002 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionBarrelEstimator.h"
00003 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionForwardEstimator.h"
00004 //
00005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00006 // Field
00007 #include "MagneticField/Engine/interface/MagneticField.h"
00008 //
00009 #include "DataFormats/CLHEP/interface/AlgebraicObjects.h"
00010 // Geometry
00011 #include "Geometry/Records/interface/IdealGeometryRecord.h"
00012 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00013 //
00014 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
00015 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00016 #include "TrackingTools/DetLayers/interface/DetLayer.h"
00017 #include "TrackingTools/MeasurementDet/interface/LayerMeasurements.h"
00018 //
00019 #include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h" 
00020 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
00021 
00022 //
00023 #include "CLHEP/Units/GlobalPhysicalConstants.h"
00024 #include "CLHEP/Geometry/Point3D.h"
00025 #include "CLHEP/Geometry/Vector3D.h" 
00026  #include "CLHEP/Geometry/Transform3D.h"  
00027 #include <cfloat>
00028 
00029 
00030 OutInConversionSeedFinder::OutInConversionSeedFinder( const edm::ParameterSet& conf ): ConversionSeedFinder( conf ), conf_(conf)  
00031 {
00032 
00033   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder CTOR " << "\n";      
00034 
00035   //the2ndHitdphi_ = 0.01; 
00036   the2ndHitdphi_ = 0.03; 
00037   the2ndHitdzConst_ = 5.;
00038   the2ndHitdznSigma_ = 2.;
00039   maxNumberOfOutInSeedsPerBC_=50;    
00040    
00041     
00042 }
00043  
00044 
00045 
00046 
00047 OutInConversionSeedFinder::~OutInConversionSeedFinder() {
00048   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder DTOR " << "\n";
00049  
00050 }
00051 
00052 
00053 
00054 void OutInConversionSeedFinder::makeSeeds( const edm::Handle<edm::View<reco::CaloCluster> > &  allBC )  const  {
00055 
00056   //  std::cout  << "  OutInConversionSeedFinder::makeSeeds() " << "\n";
00057   theSeeds_.clear();
00058 
00059 
00060   // debug  
00061   //  std::cout << "  OutInConversionSeedFinder::makeSeeds() SC position " << theSCPosition_.x() << " " << theSCPosition_.y() << " " << theSCPosition_.z() << "\n";
00062   // std::cout << " SC eta  " <<  theSCPosition_.eta() <<  " SC phi  " <<  theSCPosition_.phi() << "\n";  
00063   // std::cout << "  OutInConversionSeedFinder::makeSeeds() SC energy " << theSCenergy_  << "\n";  
00064   //
00065 
00066   findLayers();
00067 
00068   
00069   //  std::cout  << " Check Calo cluster collection size " << allBC->size() << "\n";
00070   
00071   float  theSCPhi=theSCPosition_.phi();
00072   float  theSCEta=theSCPosition_.eta();
00073 
00074   
00075 
00076   //  Loop over the Calo Clusters  in the event looking for seeds 
00077   reco::CaloClusterCollection::const_iterator bcItr;
00078   LogDebug("OutInConversionSeedFinder") << "  OutInConversionSeedFinder::makeSeeds() All BC in the event " << "\n";
00079   for (unsigned i = 0; i < allBC->size(); ++i ) {
00080   
00081     //for(bcItr = allBC.begin(); bcItr != allBC.end(); bcItr++) {
00082     nSeedsPerBC_=0;
00083 
00084     theBCEnergy_=allBC->ptrAt(i)->energy();
00085     theBCPosition_ = GlobalPoint(allBC->ptrAt(i)->position().x(), allBC->ptrAt(i)->position().y(), allBC->ptrAt(i)->position().z() ) ;
00086     float theBcEta=  theBCPosition_.eta();
00087     float theBcPhi=  theBCPosition_.phi();
00088     //    float  dPhi= theBcPhi-theSCPhi;
00089 
00090     //    std::cout << "  OutInConversionSeedFinder::makeSeeds() BC eta  " << theBcEta << " phi " <<  theBcPhi << " BC energy " << theBCEnergy_ << " dPhi " << fabs(theBcPhi-theSCPhi) << " dEta " <<  fabs(theBcEta-theSCEta) << "\n";
00091     
00092     if ( theBCEnergy_ < 1.5 ) continue;
00093 
00094     LogDebug("OutInConversionSeedFinder") << "  OutInConversionSeedFinder::makeSeeds() Passing the >=1.5 GeV cut  BC eta  " << theBcEta << " phi " <<  theBcPhi << " BC energy " << theBCEnergy_ << "\n";
00095 
00096     if (  fabs(theBcEta-theSCEta) < 0.015  && fabs(theBcPhi-theSCPhi) < 0.3 ) { 
00097       LogDebug("OutInConversionSeedFinder") << "  OutInConversionSeedFinder::makeSeeds() in et and phi range passed to the analysis " << "\n";
00098       fillClusterSeeds( allBC->ptrAt(i)  );
00099     }
00100     
00101 
00102   }
00103 
00104 
00105   //  std::cout << "Built vector of seeds of size  " << theSeeds_.size() <<  "\n" ;
00106   
00108   /*
00109   int nSeed=0;
00110   for ( std::vector<TrajectorySeed>::const_iterator iSeed= theSeeds_.begin(); iSeed != theSeeds_.end(); ++iSeed) {
00111     nSeed++;
00112     PTrajectoryStateOnDet  ptsod=iSeed->startingState();
00113     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" ;
00114     
00115     
00116     DetId tmpId = DetId( iSeed->startingState().detId());
00117     const GeomDet* tmpDet  = this->getMeasurementTracker()->geomTracker()->idToDet( tmpId );
00118     GlobalVector gv = tmpDet->surface().toGlobal( iSeed->startingState().parameters().momentum() );
00119     
00120     LogDebug("OutInConversionSeedFinder") << "seed perp,phi,eta : " 
00121                                           << gv.perp() << " , " 
00122                                           << gv.phi() << " , " 
00123                                           << gv.eta() << "\n" ; ;
00124     
00125 
00126 
00127 
00128     TrajectorySeed::range hitRange = iSeed->recHits();
00129     for (TrajectorySeed::const_iterator ihit = hitRange.first; ihit != hitRange.second; ihit++) {
00130    
00131       if ( ihit->isValid() ) {
00132 
00133         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" ;
00134 
00135       }
00136     }
00137   } 
00138   
00139   */
00140 
00141   
00142   
00143 }
00144 
00145 
00146 
00147 void OutInConversionSeedFinder::makeSeeds( const reco::CaloClusterPtr&  aBC )  const  {
00148 
00149   theSeeds_.clear();
00150 
00151   findLayers();
00152 
00153   float  theSCPhi=theSCPosition_.phi();
00154   float  theSCEta=theSCPosition_.eta();
00155 
00156   nSeedsPerBC_=0;
00157 
00158   theBCEnergy_=aBC->energy();
00159   theBCPosition_ = GlobalPoint(aBC->position().x(), aBC->position().y(), aBC->position().z() ) ;
00160   float theBcEta=  theBCPosition_.eta();
00161   float theBcPhi=  theBCPosition_.phi();
00162   //  float  dPhi= theBcPhi-theSCPhi;
00163 
00164   if ( theBCEnergy_ < 1.5 ) return;
00165 
00166   if (  fabs(theBcEta-theSCEta) < 0.015  && fabs(theBcPhi-theSCPhi) < 0.25 ) {
00167     fillClusterSeeds( aBC);
00168   }
00169 
00170 }
00171 
00172 
00173  
00174 
00175 
00176 void OutInConversionSeedFinder::fillClusterSeeds(const reco::CaloClusterPtr& bc) const {
00177 
00178   
00179   theFirstMeasurements_.clear();
00180   FreeTrajectoryState fts;  
00181 
00183   if ( makeTrackState(-1).second ) {
00184     fts = makeTrackState(-1).first;
00185     startSeed(fts);
00186   }
00188 
00189   if ( makeTrackState(1).second ) {
00190     fts = makeTrackState(1).first;
00191     startSeed(fts);
00192   }
00193 
00194 }
00195 
00196 
00197 
00198 std::pair<FreeTrajectoryState,bool>  OutInConversionSeedFinder::makeTrackState(int  charge) const {
00199 
00200   std::pair<FreeTrajectoryState,bool> result;
00201   result.second=false;
00202  
00203 
00204   //std::cout << "  OutInConversionSeedFinder:makeTrackState " << "\n";
00205 
00206 
00207   //  Old GlobalPoint gpOrigine(theBCPosition_.x()*0.3, theBCPosition_.y()*0.3, theBCPosition_.z()*0.3) ;
00208  //  GlobalPoint gpOrigine(0.,0.,0.);
00209 
00210   GlobalPoint  gpOrigine(theBeamSpot_.position().x(),theBeamSpot_.position().y(),theBeamSpot_.position().z()); 
00211   GlobalVector gvBcRadius = theBCPosition_ - gpOrigine ;
00212   HepGeom::Point3D<double>  radiusBc(gvBcRadius.x(),gvBcRadius.y(),gvBcRadius.z()) ;
00213   HepGeom::Point3D<double>  momentumWithoutCurvature = radiusBc.unit() * theBCEnergy_ ;
00214 
00215   // compute momentum direction at calo
00216   double curvature = theMF_->inTesla(theBCPosition_).z() * c_light * 1.e-3 / momentumWithoutCurvature.perp() ;
00217   curvature /= 100. ; // in cm-1 !!
00218 
00219   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState gpOrigine " << gpOrigine.x() << " " <<  gpOrigine.y() << " " <<  gpOrigine.z() << " momentumWithoutCurvature " << momentumWithoutCurvature.mag() << " curvature " << curvature << "\n";
00220 
00221   // define rotation angle
00222   float R = theBCPosition_.perp();
00223   float r = gpOrigine.perp();
00224   float rho = 1./curvature;
00225   // from the formula for the intersection of two circles
00226   // turns out to be about 2/3 of the deflection of the old formula
00227   float d = sqrt(r*r+rho*rho);
00228   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));
00229   //float u = rho + rho/d/d*(R*R-rho*rho) ;
00230   if ( u <=R )   result.second=true;
00231 
00232   double sinAlpha = 0.5*u/R;
00233   if ( sinAlpha>(1.-10*DBL_EPSILON) )  sinAlpha = 1.-10*DBL_EPSILON;
00234   else if ( sinAlpha<-(1.-10*DBL_EPSILON) )  sinAlpha = -(1.-10*DBL_EPSILON);
00235   
00236   double newdphi = charge * asin( sinAlpha) ;
00237 
00238   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState charge " << charge << " R " << R << " u/R " << u/R << " asin(0.5*u/R) " << asin(sinAlpha) << "\n";
00239 
00240   HepGeom::Transform3D rotation =  HepGeom::Rotate3D(newdphi, HepGeom::Vector3D<double> (0., 0. ,1.));
00241 
00242 
00243   HepGeom::Point3D<double>  momentumInTracker = momentumWithoutCurvature.transform(rotation) ;
00244   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState  R " << R << " r " << r << " rho " << rho  << " d " << d  << " u " << u << " newdphi " << newdphi << " momentumInTracker " <<  momentumInTracker << "\n";
00245 
00246   HepGeom::Point3D<double>  hepStartingPoint(gpOrigine.x(), gpOrigine.y(), gpOrigine.z()) ;
00247 
00248   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState hepStartingPoint " << hepStartingPoint << "\n";
00249 
00250   hepStartingPoint.transform(rotation);
00251 
00252   GlobalPoint startingPoint(hepStartingPoint.x(), hepStartingPoint.y(), hepStartingPoint.z());
00253 
00254   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeTrackState startingPoint " << startingPoint << " calo position " << theBCPosition_ << "\n";
00255   GlobalVector gvTracker(momentumInTracker.x(), momentumInTracker.y(), momentumInTracker.z());
00256   GlobalTrajectoryParameters gtp(startingPoint, gvTracker, charge, &(*theMF_) );
00257   // error matrix
00258   AlgebraicSymMatrix55 m = AlgebraicMatrixID();
00259   m(0,0) = 0.1; m(1,1) = 0.1 ; m(2,2) = 0.1 ;
00260   m(3,3) = 0.1 ; m(4,4) = 0.1;
00261   
00262   //std::cout << "OutInConversionSeedFinder::makeTrackState " <<  FreeTrajectoryState(gtp, CurvilinearTrajectoryError(m) ) << "\n";
00263    
00264   result.first= FreeTrajectoryState(gtp, CurvilinearTrajectoryError(m) ) ;
00265   return result;
00266 
00267 }
00268 
00269 
00270 void OutInConversionSeedFinder::startSeed(const FreeTrajectoryState & fts) const {
00271 
00272 
00273   //  std::cout << "OutInConversionSeedFinder::startSeed layer list " << this->layerList().size() <<  "\n";
00274   //std::cout << "OutInConversionSeedFinder::startSeed  fts " << fts <<  "\n";  
00275 
00276   std::vector<const DetLayer*> myLayers=layerList();
00277   if ( myLayers.size() > 3 ) {
00278     
00279     for(unsigned int ilayer = myLayers.size()-1; ilayer >= myLayers.size()-2; --ilayer) {
00280       const DetLayer * layer = myLayers[ilayer];
00281       
00282       
00283       // allow the z of the hit to be within a straight line from a vertex
00284       // of +-15 cm to the cluster
00285       //      float dphi = 0.015;
00286       float dphi = 0.030;
00287 
00288       MeasurementEstimator * newEstimator = makeEstimator(layer, dphi);
00289 
00290      
00291       //std::cout << "OutInSeedFinder::startSeed propagationDirection  " << int(thePropagatorAlongMomentum_->propagationDirection() ) << "\n";       
00292       
00293       TSOS tsos(fts, layer->surface() );
00294       
00295       LogDebug("OutInConversionSeedFinder") << "OutInSeedFinder::startSeed  after  TSOS tsos(fts, layer->surface() ) " << "\n";
00296       
00297       LayerMeasurements theLayerMeasurements_(this->getMeasurementTracker() );
00298       theFirstMeasurements_ = theLayerMeasurements_.measurements( *layer, tsos, *thePropagatorAlongMomentum_, *newEstimator);
00299       
00300       //std::cout << "OutInSeedFinder::startSeed  after  theFirstMeasurements_   " << theFirstMeasurements_.size() <<  "\n";
00301       
00302       if(theFirstMeasurements_.size() > 1) // always a dummy returned, too
00303         LogDebug("OutInConversionSeedFinder") <<  " Found " << theFirstMeasurements_.size()-1 << " 1st hits in seed" << "\n";
00304       
00305       delete newEstimator;
00306       
00307       LogDebug("OutInConversionSeedFinder") << "OutInSeedFinder::startSeed Layer " << ilayer << " theFirstMeasurements_.size " << theFirstMeasurements_.size() << "\n";
00308       
00309       for(unsigned int i = 0; i < theFirstMeasurements_.size(); ++i) {
00310         TrajectoryMeasurement m1 = theFirstMeasurements_[i];
00311         if(m1.recHit()->isValid()) {
00312           
00313           // update the fts to start from this point.  much better than starting from
00314           // extrapolated point along the line
00315           GlobalPoint hitPoint = m1.recHit()->globalPosition();
00316           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";
00317 
00318           
00319           FreeTrajectoryState newfts = trackStateFromClusters(fts.charge(), hitPoint, alongMomentum, 0.8);
00320           LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::startSeed  newfts " << newfts << "\n";
00321           LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::startSeed propagationDirection  after switching " << int(thePropagatorOppositeToMomentum_->propagationDirection() ) << "\n";        
00322           //  std::cout << "OutInConversionSeedFinder::startSeed propagationDirection  after switching " << int(thePropagatorOppositeToMomentum_->propagationDirection() ) << "\n";        
00323 
00324        
00325           completeSeed(m1, newfts, thePropagatorOppositeToMomentum_, ilayer-1);
00326           // skip a layer, if you haven't already skipped the first layer
00327           if(ilayer == myLayers.size()-1) {
00328             completeSeed(m1, newfts, thePropagatorOppositeToMomentum_, ilayer-2);
00329           }
00330         }
00331       }
00332       
00333     } // loop over layers
00334   }
00335 
00336 
00337 
00338   
00339 }
00340 
00341 
00342 
00343 MeasurementEstimator * OutInConversionSeedFinder::makeEstimator(const DetLayer * layer, float dphi) const {
00344  
00345   //std::cout  << "OutInConversionSeedFinder::makeEstimator  " << "\n";
00346 
00347   MeasurementEstimator * newEstimator=0;
00348 
00349   if (layer->location() == GeomDetEnumerators::barrel ) {
00350     
00351     const BarrelDetLayer * barrelLayer = dynamic_cast<const BarrelDetLayer*>(layer);
00352     LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeEstimator Barrel  r = " << barrelLayer->specificSurface().radius() << " " << "\n";        
00353     float r = barrelLayer->specificSurface().radius();
00354     float zrange = 15.* (1.-r/theBCPosition_.perp());
00355     newEstimator = new ConversionBarrelEstimator(-dphi, dphi, -zrange, zrange);
00356   }
00357 
00358 
00359 
00360   if (layer->location() == GeomDetEnumerators::endcap ) {   
00361    
00362     const ForwardDetLayer * forwardLayer = dynamic_cast<const ForwardDetLayer*>(layer);
00363     LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::makeEstimator Endcap r = " << forwardLayer->specificSurface().innerRadius() << " R " << forwardLayer->specificSurface().outerRadius()  <<  " Z " << forwardLayer->specificSurface().position().z() << "\n";  
00364     
00365     float zc = fabs(theBCPosition_.z());
00366     float z =  fabs(forwardLayer->surface().position().z());
00367     
00368     float rrange = 15. * theBCPosition_.perp() * (zc - z) / (zc*zc - 15.*zc);
00369     newEstimator = new ConversionForwardEstimator(-dphi, dphi, rrange);
00370   }
00371 
00372 
00373 
00374 
00375   return newEstimator;
00376 }
00377 
00378 
00379 
00380 
00381 void OutInConversionSeedFinder::completeSeed(const TrajectoryMeasurement & m1, 
00382                                              FreeTrajectoryState & fts, 
00383                                              const Propagator* propagator, int ilayer) const {
00384 
00385   //std::cout <<  "OutInConversionSeedFinder::completeSeed ilayer " << ilayer << "\n";
00386 
00387   MeasurementEstimator * newEstimator=0;
00388   const DetLayer * layer = theLayerList_[ilayer];
00389   
00390 
00391   if ( layer->location() == GeomDetEnumerators::barrel ) {
00392     // z error for 2nd hit is  2 sigma quadded with 5 cm
00393     LogDebug("OutInConversionSeedFinder") << " Barrel OutInConversionSeedFinder::completeSeed " << the2ndHitdznSigma_ << " " << the2ndHitdzConst_ << " " << the2ndHitdphi_ << "\n";
00394     float dz = sqrt(the2ndHitdznSigma_*the2ndHitdznSigma_*m1.recHit()->globalPositionError().czz() 
00395                     + the2ndHitdzConst_*the2ndHitdzConst_);
00396     newEstimator =
00397       new ConversionBarrelEstimator(-the2ndHitdphi_, the2ndHitdphi_, -dz, dz);
00398   }
00399   else {
00400     LogDebug("OutInConversionSeedFinder") << " EndCap OutInConversionSeedFinder::completeSeed " << the2ndHitdznSigma_ << " " << the2ndHitdzConst_ << " " << the2ndHitdphi_ << "\n";
00401     // z error for 2nd hit is 2sigma quadded with 5 cm
00402     //float m1dr = m1.recHit().globalPositionError().rerr(m1.recHit().globalPosition());
00403     float m1dr = sqrt(m1.recHit()->localPositionError().yy());
00404     float dr = sqrt(the2ndHitdznSigma_*the2ndHitdznSigma_*m1dr*m1dr 
00405                     + the2ndHitdzConst_*the2ndHitdznSigma_);
00406     // LogDebug("OutInConversionSeedFinder") << "second hit forward dr " << dr << " this hit " << m1dr << endl;
00407     newEstimator =
00408       new ConversionForwardEstimator(-the2ndHitdphi_, the2ndHitdphi_, dr);
00409   }
00410 
00411   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::completeSeed  ilayer " << ilayer <<  "\n"; 
00412   
00413   // Get the measurements consistent with the FTS and the Estimator
00414   TSOS tsos(fts, layer->surface() );
00415   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::completeSeed propagationDirection  " << int(propagator->propagationDirection() ) << "\n";               
00416   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::completeSeed pointer to estimator " << newEstimator << "\n";
00417 
00418   LayerMeasurements theLayerMeasurements_(this->getMeasurementTracker() );
00419   std::vector<TrajectoryMeasurement> measurements = theLayerMeasurements_.measurements( *layer, tsos, *propagator, *newEstimator);
00420   //std::cout << "OutInConversionSeedFinder::completeSeed Found " << measurements.size() << " second hits " << "\n";
00421   delete newEstimator;
00422 
00423   for(unsigned int i = 0; i < measurements.size(); ++i) {
00424     if( measurements[i].recHit()->isValid()  ) {
00425       createSeed(m1, measurements[i]);
00426     }
00427   }
00428 
00429 
00430 
00431   //LogDebug("OutInConversionSeedFinder") << "COMPLETED " << theSeeds_.size() << " SEEDS " << "\n";
00432 }
00433 
00434 
00435 
00436 void OutInConversionSeedFinder::createSeed(const TrajectoryMeasurement & m1, 
00437                                            const TrajectoryMeasurement & m2) const {
00438 
00439   //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";
00440   
00441 
00442   FreeTrajectoryState fts = createSeedFTS(m1, m2);
00443 
00444 
00445   LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::createSeed First point errors " <<m1.recHit()->parametersError() << "\n";
00446   LogDebug("OutInConversionSeedFinder") << "original cluster FTS " << fts <<"\n";
00447 
00448 
00449 
00450   //std::cout  << "OutInConversionSeedFinder::createSeed propagation dir " << int( thePropagatorOppositeToMomentum_->propagationDirection() ) << "\n"; 
00451   TrajectoryStateOnSurface state1 = thePropagatorOppositeToMomentum_->propagate(fts,  m1.recHit()->det()->surface());
00452 
00453   // LogDebug("OutInConversionSeedFinder") << "hit surface " << h1.det().surface().position() << endl;
00454   // LogDebug("OutInConversionSeedFinder") << "prop to " << typeid(h1.det().surface()).name() << endl;
00455   // LogDebug("OutInConversionSeedFinder") << "prop to first hit " << state1 << endl; 
00456   // LogDebug("OutInConversionSeedFinder") << "update to " << h1.globalPosition() << endl;
00457 
00458   if ( state1.isValid() ) {
00459     TrajectoryStateOnSurface updatedState1 = theUpdator_.update(state1,  *m1.recHit() );
00460 
00461     if ( updatedState1.isValid() ) {
00462       TrajectoryStateOnSurface state2 = thePropagatorOppositeToMomentum_->propagate(*updatedState1.freeTrajectoryState(),  m2.recHit()->det()->surface());
00463 
00464       if ( state2.isValid() ) {
00465 
00466         TrajectoryStateOnSurface updatedState2 = theUpdator_.update(state2, *m2.recHit() );
00467         TrajectoryMeasurement meas1(state1, updatedState1,  m1.recHit()  , m1.estimate(), m1.layer());
00468         TrajectoryMeasurement meas2(state2, updatedState2,  m2.recHit()  , m2.estimate(), m2.layer());
00469         
00470         edm::OwnVector<TrackingRecHit> myHits;
00471         myHits.push_back(meas1.recHit()->hit()->clone());
00472         myHits.push_back(meas2.recHit()->hit()->clone());
00473 
00474         if ( nSeedsPerBC_ >= maxNumberOfOutInSeedsPerBC_ ) return;
00475    
00476         TrajectoryStateTransform tsTransform;
00477         PTrajectoryStateOnDet* ptsod= tsTransform.persistentState(state2, meas2.recHit()->hit()->geographicalId().rawId()  );
00478 
00479         LogDebug("OutInConversionSeedFinder") << "OutInConversionSeedFinder::createSeed new seed  from state " << state2.globalPosition()  <<  "\n";
00480         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" ;
00481         
00482 
00483 
00484         theSeeds_.push_back(TrajectorySeed( *ptsod, myHits, oppositeToMomentum ));
00485         nSeedsPerBC_++;
00486  
00487         delete ptsod;  
00488         
00489 
00490       }
00491     }
00492   }
00493 
00494 
00495 
00496 
00497 
00498 }
00499 
00500 
00501 
00502 
00503 
00504 FreeTrajectoryState OutInConversionSeedFinder::createSeedFTS(const TrajectoryMeasurement & m1,
00505                                                              const TrajectoryMeasurement & m2) const {
00506 
00507   //std::cout  << "OutInConversionSeedFinder::createSeedFTS " << "\n";
00508 
00509   GlobalPoint xmeas = fixPointRadius(m1);
00510   GlobalPoint xvert = fixPointRadius(m2);
00511 
00512 
00513   float pt = theSCenergy_ * sin(theSCPosition_.theta());
00514   float pz = theSCenergy_ * cos(theSCPosition_.theta());
00515 
00516 
00517 
00518   // doesn't work at all for endcap, where r is badly measured
00519   //float dphidr = (p1.phi()-p2.phi())/(p1.perp()-p2.perp());
00520   //int charge = (dphidr > 0.) ? -1 : 1;
00521   int charge = m2.predictedState().charge();
00522 
00523   double BInTesla = theMF_->inTesla(xmeas).z();
00524   GlobalVector xdiff = xmeas -xvert;
00525  
00526   double phi= xdiff.phi();
00527   double pxOld = pt*cos(phi);
00528   double pyOld = pt*sin(phi);
00529   double RadCurv = 100*pt/(BInTesla*0.29979);
00530   double alpha = asin(0.5*xdiff.perp()/RadCurv);
00531   float ca = cos(charge*alpha);
00532   float sa = sin(charge*alpha);
00533   double pxNew =   ca*pxOld + sa*pyOld;
00534   double pyNew =  -sa*pxOld + ca*pyOld;
00535   GlobalVector pNew(pxNew, pyNew, pz);
00536 
00537   GlobalTrajectoryParameters gp(xmeas, pNew, charge, &(*theMF_) );
00538 
00539   AlgebraicSymMatrix55 m = AlgebraicMatrixID();
00540   m(0,0) = 0.05; m(1,1) = 0.02 ; m(2,2) = 0.007 ;
00541   m(3,3) = 10. ; m(4,4) = 10. ;
00542   return FreeTrajectoryState(gp, CurvilinearTrajectoryError(m));
00543 
00544 
00545 }
00546 
00547 
00548 
00549 
00550 GlobalPoint OutInConversionSeedFinder::fixPointRadius(const  TrajectoryMeasurement& m1) const {
00551   GlobalPoint p1 = m1.recHit()->globalPosition();
00552   GlobalPoint p2;
00553   if(m1.layer()->location() == GeomDetEnumerators::barrel) {
00554     p2 = p1;
00555   } else {
00556     float z = p1.z();
00557     float phi = p1.phi();
00558     float theta = theSCPosition_.theta();
00559     float r = p1.z() * tan(theta);
00560     p2 = GlobalPoint(r*cos(phi), r*sin(phi), z);
00561     // LogDebug("OutInConversionSeedFinder") << "fixing point radius " << p2 << " from " << p1 << endl;
00562   }
00563   return p2;
00564 }
00565 
00566 
00567 
00568 
00569