Reference Manual

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    EgammaElectronAlgos
00004 // Class:      ElectronSeedGenerator.
00005 //
00013 //
00014 // Original Author:  Ursula Berthon, Claude Charlot
00015 //         Created:  Mon Mar 27 13:22:06 CEST 2006
00016 //
00017 
00018 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronSeedGenerator.h"
00019 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronUtilities.h"
00020 
00021 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
00022 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
00023 #include "RecoTracker/TkSeedGenerator/interface/FastHelix.h"
00024 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
00025 #include "RecoTracker/Record/interface/NavigationSchoolRecord.h"
00026 
00027 //#include "DataFormats/EgammaReco/interface/EcalCluster.h"
00028 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
00029 
00030 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
00031 //#include "DataFormats/BeamSpot/interface/BeamSpot.h"
00032 //#include "DataFormats/Common/interface/Handle.h"
00033 
00034 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00035 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00036 
00037 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
00038 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
00039 #include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"
00040 
00041 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00042 
00043 #include <vector>
00044 #include <utility>
00045 
00046 ElectronSeedGenerator::ElectronSeedGenerator(const edm::ParameterSet &pset)
00047  : dynamicphiroad_(pset.getParameter<bool>("dynamicPhiRoad")),
00048    fromTrackerSeeds_(pset.getParameter<bool>("fromTrackerSeeds")),
00049    useRecoVertex_(false),
00050    verticesTag_("offlinePrimaryVerticesWithBS"),
00051    beamSpotTag_("offlineBeamSpot"),
00052    lowPtThreshold_(pset.getParameter<double>("LowPtThreshold")),
00053    highPtThreshold_(pset.getParameter<double>("HighPtThreshold")),
00054    nSigmasDeltaZ1_(pset.getParameter<double>("nSigmasDeltaZ1")),
00055    deltaZ1WithVertex_(0.5),
00056    sizeWindowENeg_(pset.getParameter<double>("SizeWindowENeg")),
00057    deltaPhi1Low_(pset.getParameter<double>("DeltaPhi1Low")),
00058    deltaPhi1High_(pset.getParameter<double>("DeltaPhi1High")),
00059    deltaPhi1Coef1_(0.), deltaPhi1Coef2_(0.),
00060    myMatchEle(0), myMatchPos(0),
00061    thePropagator(0),
00062    theMeasurementTracker(0),
00063    theSetup(0), 
00064    cacheIDMagField_(0),/*cacheIDGeom_(0),*/cacheIDNavSchool_(0),cacheIDCkfComp_(0),cacheIDTrkGeom_(0)
00065  {
00066   // so that deltaPhi1 = deltaPhi1Coef1_ + deltaPhi1Coef2_/clusterEnergyT
00067   if (dynamicphiroad_)
00068    {
00069     deltaPhi1Coef2_ = (deltaPhi1Low_-deltaPhi1High_)/(1./lowPtThreshold_-1./highPtThreshold_) ;
00070     deltaPhi1Coef1_ = deltaPhi1Low_ - deltaPhi1Coef2_/lowPtThreshold_ ;
00071    }
00072 
00073   // use of a theMeasurementTrackerName
00074   if (pset.exists("measurementTrackerName"))
00075    { theMeasurementTrackerName = pset.getParameter<std::string>("measurementTrackerName") ; }
00076 
00077   // use of reco vertex
00078   if (pset.exists("useRecoVertex"))
00079    { useRecoVertex_ = pset.getParameter<bool>("useRecoVertex") ; }
00080   if (pset.exists("vertices"))
00081    { verticesTag_ = pset.getParameter<edm::InputTag>("vertices") ; }
00082   if (pset.exists("deltaZ1WithVertex"))
00083    { deltaZ1WithVertex_ = pset.getParameter<double>("deltaZ1WithVertex") ; }
00084 
00085   // new beamSpot tag
00086   if (pset.exists("beamSpot"))
00087    { beamSpotTag_ = pset.getParameter<edm::InputTag>("beamSpot") ; }
00088 
00089   // new B/F configurables
00090   if (pset.exists("DeltaPhi2"))
00091    { deltaPhi2B_ = deltaPhi2F_ = pset.getParameter<double>("DeltaPhi2") ; }
00092   else
00093    {
00094     deltaPhi2B_ = pset.getParameter<double>("DeltaPhi2B") ;
00095     deltaPhi2F_ = pset.getParameter<double>("DeltaPhi2F") ;
00096    }
00097   if (pset.exists("PhiMin2"))
00098    { phiMin2B_ = phiMin2F_ = pset.getParameter<double>("PhiMin2") ; }
00099   else
00100    {
00101     phiMin2B_ = pset.getParameter<double>("PhiMin2B") ;
00102     phiMin2F_ = pset.getParameter<double>("PhiMin2F") ;
00103    }
00104   if (pset.exists("PhiMax2"))
00105    { phiMax2B_ = phiMax2F_ = pset.getParameter<double>("PhiMax2") ; }
00106   else
00107    {
00108     phiMax2B_ = pset.getParameter<double>("PhiMax2B") ;
00109     phiMax2F_ = pset.getParameter<double>("PhiMax2F") ;
00110    }
00111 
00112   // Instantiate the pixel hit matchers
00113   myMatchEle = new PixelHitMatcher
00114    ( pset.getParameter<double>("ePhiMin1"),
00115                  pset.getParameter<double>("ePhiMax1"),
00116                  phiMin2B_, phiMax2B_, phiMin2F_, phiMax2F_,
00117      pset.getParameter<double>("z2MinB"),
00118      pset.getParameter<double>("z2MaxB"),
00119      pset.getParameter<double>("r2MinF"),
00120      pset.getParameter<double>("r2MaxF"),
00121      pset.getParameter<double>("rMinI"),
00122      pset.getParameter<double>("rMaxI"),
00123      pset.getParameter<bool>("searchInTIDTEC") ) ;
00124 
00125   myMatchPos = new PixelHitMatcher
00126    ( pset.getParameter<double>("pPhiMin1"),
00127                  pset.getParameter<double>("pPhiMax1"),
00128                  phiMin2B_, phiMax2B_, phiMin2F_, phiMax2F_,
00129      pset.getParameter<double>("z2MinB"),
00130      pset.getParameter<double>("z2MaxB"),
00131      pset.getParameter<double>("r2MinF"),
00132      pset.getParameter<double>("r2MaxF"),
00133      pset.getParameter<double>("rMinI"),
00134      pset.getParameter<double>("rMaxI"),
00135      pset.getParameter<bool>("searchInTIDTEC") ) ;
00136 
00137   theUpdator = new KFUpdator() ;
00138  }
00139 
00140 ElectronSeedGenerator::~ElectronSeedGenerator()
00141  {
00142   delete myMatchEle ;
00143   delete myMatchPos ;
00144   delete thePropagator ;
00145   delete theUpdator ;
00146  }
00147 
00148 void ElectronSeedGenerator::setupES(const edm::EventSetup& setup) {
00149 
00150   // get records if necessary (called once per event)
00151   bool tochange=false;
00152 
00153   if (cacheIDMagField_!=setup.get<IdealMagneticFieldRecord>().cacheIdentifier()) {
00154     setup.get<IdealMagneticFieldRecord>().get(theMagField);
00155     cacheIDMagField_=setup.get<IdealMagneticFieldRecord>().cacheIdentifier();
00156     if (thePropagator) delete thePropagator;
00157     thePropagator = new PropagatorWithMaterial(alongMomentum,.000511,&(*theMagField));
00158     tochange=true;
00159   }
00160 
00161   if (!fromTrackerSeeds_ && cacheIDCkfComp_!=setup.get<CkfComponentsRecord>().cacheIdentifier()) {
00162     edm::ESHandle<MeasurementTracker> measurementTrackerHandle;
00163     setup.get<CkfComponentsRecord>().get(theMeasurementTrackerName,measurementTrackerHandle);
00164     cacheIDCkfComp_=setup.get<CkfComponentsRecord>().cacheIdentifier();
00165     theMeasurementTracker = measurementTrackerHandle.product();
00166     tochange=true;
00167   }
00168 
00169   //edm::ESHandle<TrackerGeometry> trackerGeometryHandle;
00170   if (cacheIDTrkGeom_!=setup.get<TrackerDigiGeometryRecord>().cacheIdentifier()) {
00171     cacheIDTrkGeom_=setup.get<TrackerDigiGeometryRecord>().cacheIdentifier();
00172     setup.get<TrackerDigiGeometryRecord>().get(theTrackerGeometry);
00173     tochange=true; //FIXME
00174   }
00175 
00176   if (tochange) {
00177     myMatchEle->setES(&(*theMagField),theMeasurementTracker,theTrackerGeometry.product());
00178     myMatchPos->setES(&(*theMagField),theMeasurementTracker,theTrackerGeometry.product());
00179   }
00180 
00181   if (cacheIDNavSchool_!=setup.get<NavigationSchoolRecord>().cacheIdentifier()) {
00182     edm::ESHandle<NavigationSchool> nav;
00183     setup.get<NavigationSchoolRecord>().get("SimpleNavigationSchool", nav);
00184     cacheIDNavSchool_=setup.get<NavigationSchoolRecord>().cacheIdentifier();
00185     theNavigationSchool = nav.product();
00186   }
00187 
00188 //  if (cacheIDGeom_!=setup.get<TrackerRecoGeometryRecord>().cacheIdentifier()) {
00189 //    setup.get<TrackerRecoGeometryRecord>().get( theGeomSearchTracker );
00190 //    cacheIDGeom_=setup.get<TrackerRecoGeometryRecord>().cacheIdentifier();
00191 //  }
00192 
00193 }
00194 
00195 void display_seed( const std::string & title1, const std::string & title2, const reco::ElectronSeed & seed, edm::ESHandle<TrackerGeometry> trackerGeometry )
00196  {
00197   const PTrajectoryStateOnDet & startingState = seed.startingState() ;
00198   const LocalTrajectoryParameters & parameters = startingState.parameters() ;
00199   std::cout<<title1
00200     <<" ("<<seed.subDet2()<<"/"<<seed.dRz2()<<"/"<<seed.dPhi2()<<")"
00201     <<" ("<<seed.direction()<<"/"<<startingState.detId()<<"/"<<startingState.surfaceSide()<<"/"<<parameters.charge()<<"/"<<parameters.position()<<"/"<<parameters.momentum()<<")"
00202     <<std::endl ;
00203  }
00204 
00205 bool equivalent( const TrajectorySeed & s1, const TrajectorySeed & s2 )
00206  {
00207   if (s1.nHits()!=s2.nHits()) return false ;
00208 
00209   unsigned int nHits ;
00210   TrajectorySeed::range r1 = s1.recHits(), r2 = s2.recHits() ;
00211   TrajectorySeed::const_iterator i1, i2 ;
00212   for ( i1=r1.first, i2=r2.first, nHits=0 ; i1!=r1.second ; ++i1, ++i2, ++nHits )
00213    {
00214     if ( !i1->isValid() || !i2->isValid() ) return false ;
00215     if ( i1->geographicalId()!=i2->geographicalId() ) return false ;
00216     if ( ! ( i1->localPosition()==i2->localPosition() ) ) return false ;
00217    }
00218 
00219   return true ;
00220  }
00221 
00222 void  ElectronSeedGenerator::run
00223  ( edm::Event & e, const edm::EventSetup & setup,
00224    const reco::SuperClusterRefVector & sclRefs, const std::vector<float> & hoe1s, const std::vector<float> & hoe2s,
00225    TrajectorySeedCollection * seeds, reco::ElectronSeedCollection & out )
00226  {
00227   theInitialSeedColl = seeds ;
00228 //  bool duplicateTrajectorySeeds =false ;
00229 //  unsigned int i,j ;
00230 //  for (i=0;i<seeds->size();++i)
00231 //    for (j=i+1;j<seeds->size();++j)
00232 //     {
00233 //      const TrajectorySeed & s1 =(*seeds)[i] ;
00234 //      const TrajectorySeed & s2 =(*seeds)[j] ;
00235 //      if ( equivalent(s1,s2) )
00236 //       {
00237 //        const PTrajectoryStateOnDet & ss1 = s1.startingState() ;
00238 //        const LocalTrajectoryParameters & p1 = ss1.parameters() ;
00239 //        const PTrajectoryStateOnDet & ss2 = s2.startingState() ;
00240 //        const LocalTrajectoryParameters & p2 = ss2.parameters() ;
00241 //        duplicateTrajectorySeeds = true ;
00242 //        std::cout<<"Same hits for "
00243 //          <<"\n  s["<<i<<"] ("<<s1.direction()<<"/"<<ss1.detId()<<"/"<<ss1.surfaceSide()<<"/"<<p1.charge()<<"/"<<p1.position()<<"/"<<p1.momentum()<<")"
00244 //          <<"\n  s["<<j<<"] ("<<s2.direction()<<"/"<<ss2.detId()<<"/"<<ss2.surfaceSide()<<"/"<<p2.charge()<<"/"<<p2.position()<<"/"<<p2.momentum()<<")"
00245 //          <<std::endl ;
00246 //       }
00247 //     }
00248 //  if (duplicateTrajectorySeeds)
00249 //   { edm::LogWarning("ElectronSeedGenerator|DuplicateTrajectorySeeds")<<"We see several identical trajectory seeds." ; }
00250 
00251   //Retrieve tracker topology from geometry
00252   edm::ESHandle<TrackerTopology> tTopoHand;
00253   setup.get<IdealGeometryRecord>().get(tTopoHand);
00254   const TrackerTopology *tTopo=tTopoHand.product();
00255 
00256   theSetup= &setup;
00257   NavigationSetter theSetter(*theNavigationSchool);
00258 
00259   // get initial TrajectorySeeds if necessary
00260   //  if (fromTrackerSeeds_) e.getByLabel(initialSeeds_, theInitialSeedColl);
00261 
00262   // get the beamspot from the Event:
00263   //e.getByType(theBeamSpot);
00264   e.getByLabel(beamSpotTag_,theBeamSpot);
00265 
00266   // if required get the vertices
00267   if (useRecoVertex_) e.getByLabel(verticesTag_,theVertices);
00268 
00269   if (!fromTrackerSeeds_)
00270    { theMeasurementTracker->update(e) ; }
00271 
00272   for  (unsigned int i=0;i<sclRefs.size();++i) {
00273     // Find the seeds
00274     recHits_.clear();
00275 
00276     LogDebug ("run") << "new cluster, calling seedsFromThisCluster";
00277     seedsFromThisCluster(sclRefs[i],hoe1s[i],hoe2s[i],out,tTopo);
00278   }
00279 
00280   LogDebug ("run") << ": For event "<<e.id();
00281   LogDebug ("run") <<"Nr of superclusters after filter: "<<sclRefs.size()
00282    <<", no. of ElectronSeeds found  = " << out.size();
00283 }
00284 
00285 void ElectronSeedGenerator::seedsFromThisCluster
00286 ( edm::Ref<reco::SuperClusterCollection> seedCluster,
00287   float hoe1, float hoe2,
00288   reco::ElectronSeedCollection & out, const TrackerTopology *tTopo )
00289 {
00290   float clusterEnergy = seedCluster->energy() ;
00291   GlobalPoint clusterPos
00292     ( seedCluster->position().x(),
00293       seedCluster->position().y(),
00294       seedCluster->position().z() ) ;
00295   reco::ElectronSeed::CaloClusterRef caloCluster(seedCluster) ;
00296 
00297   //LogDebug("") << "[ElectronSeedGenerator::seedsFromThisCluster] new supercluster with energy: " << clusterEnergy ;
00298   //LogDebug("") << "[ElectronSeedGenerator::seedsFromThisCluster] and position: " << clusterPos ;
00299 
00300   if (dynamicphiroad_)
00301    {
00302     float clusterEnergyT = clusterEnergy / cosh( EleRelPoint(clusterPos,theBeamSpot->position()).eta() ) ;
00303 
00304     float deltaPhi1 ;
00305     if (clusterEnergyT < lowPtThreshold_)
00306      { deltaPhi1= deltaPhi1Low_ ; }
00307     else if (clusterEnergyT > highPtThreshold_)
00308      { deltaPhi1= deltaPhi1High_ ; }
00309     else
00310      { deltaPhi1 = deltaPhi1Coef1_ + deltaPhi1Coef2_/clusterEnergyT ; }
00311 
00312     float ephimin1 = -deltaPhi1*sizeWindowENeg_ ;
00313     float ephimax1 =  deltaPhi1*(1.-sizeWindowENeg_);
00314     float pphimin1 = -deltaPhi1*(1.-sizeWindowENeg_);
00315     float pphimax1 =  deltaPhi1*sizeWindowENeg_;
00316 
00317     float phimin2B  = -deltaPhi2B_/2. ;
00318     float phimax2B  =  deltaPhi2B_/2. ;
00319     float phimin2F  = -deltaPhi2F_/2. ;
00320     float phimax2F  =  deltaPhi2F_/2. ;
00321 
00322 
00323     myMatchEle->set1stLayer(ephimin1,ephimax1);
00324     myMatchPos->set1stLayer(pphimin1,pphimax1);
00325     myMatchEle->set2ndLayer(phimin2B,phimax2B, phimin2F,phimax2F);
00326     myMatchPos->set2ndLayer(phimin2B,phimax2B, phimin2F,phimax2F);
00327    }
00328 
00329   PropagationDirection dir = alongMomentum ;
00330 
00331   if (!useRecoVertex_) // here use the beam spot position
00332    {
00333     double sigmaZ=theBeamSpot->sigmaZ();
00334     double sigmaZ0Error=theBeamSpot->sigmaZ0Error();
00335     double sq=sqrt(sigmaZ*sigmaZ+sigmaZ0Error*sigmaZ0Error);
00336     double myZmin1=theBeamSpot->position().z()-nSigmasDeltaZ1_*sq;
00337     double myZmax1=theBeamSpot->position().z()+nSigmasDeltaZ1_*sq;
00338 
00339     GlobalPoint vertexPos ;
00340     ele_convert(theBeamSpot->position(),vertexPos) ;
00341 
00342     myMatchEle->set1stLayerZRange(myZmin1,myZmax1);
00343     myMatchPos->set1stLayerZRange(myZmin1,myZmax1);
00344 
00345     if (!fromTrackerSeeds_)
00346      {
00347       // try electron
00348       std::vector<std::pair<RecHitWithDist,ConstRecHitPointer> > elePixelHits
00349        = myMatchEle->compatibleHits(clusterPos,vertexPos,
00350                                     clusterEnergy,-1., tTopo) ;
00351       GlobalPoint eleVertex(theBeamSpot->position().x(),theBeamSpot->position().y(),myMatchEle->getVertex()) ;
00352       seedsFromRecHits(elePixelHits,dir,eleVertex,caloCluster,out,false) ;
00353       // try positron
00354       std::vector<std::pair<RecHitWithDist,ConstRecHitPointer> > posPixelHits
00355         = myMatchPos->compatibleHits(clusterPos,vertexPos,clusterEnergy,1.,tTopo) ;
00356       GlobalPoint posVertex(theBeamSpot->position().x(),theBeamSpot->position().y(),myMatchPos->getVertex()) ;
00357       seedsFromRecHits(posPixelHits,dir,posVertex,caloCluster,out,true) ;
00358      }
00359     else
00360      {
00361       // try electron
00362       std::vector<SeedWithInfo> elePixelSeeds
00363        = myMatchEle->compatibleSeeds(theInitialSeedColl,clusterPos,vertexPos,clusterEnergy,-1.) ;
00364       seedsFromTrajectorySeeds(elePixelSeeds,caloCluster,hoe1,hoe2,out,false) ;
00365       // try positron
00366       std::vector<SeedWithInfo> posPixelSeeds
00367        = myMatchPos->compatibleSeeds(theInitialSeedColl,clusterPos,vertexPos,clusterEnergy,1.) ;
00368       seedsFromTrajectorySeeds(posPixelSeeds,caloCluster,hoe1,hoe2,out,true) ;
00369      }
00370 
00371    }
00372   else // here we use the reco vertices
00373    {
00374 
00375     myMatchEle->setUseRecoVertex(true) ; //Hit matchers need to know that the vertex is known
00376     myMatchPos->setUseRecoVertex(true) ;
00377 
00378     const  std::vector<reco::Vertex> * vtxCollection = theVertices.product() ;
00379     std::vector<reco::Vertex>::const_iterator vtxIter ;
00380     for (vtxIter = vtxCollection->begin(); vtxIter != vtxCollection->end() ; vtxIter++)
00381      {
00382 
00383       GlobalPoint vertexPos(vtxIter->position().x(),vtxIter->position().y(),vtxIter->position().z());
00384       double myZmin1, myZmax1 ;
00385       if (vertexPos.z()==theBeamSpot->position().z())
00386        { // in case vetex not found
00387         double sigmaZ=theBeamSpot->sigmaZ();
00388         double sigmaZ0Error=theBeamSpot->sigmaZ0Error();
00389         double sq=sqrt(sigmaZ*sigmaZ+sigmaZ0Error*sigmaZ0Error);
00390         myZmin1=theBeamSpot->position().z()-nSigmasDeltaZ1_*sq;
00391         myZmax1=theBeamSpot->position().z()+nSigmasDeltaZ1_*sq;
00392        }
00393       else
00394        { // a vertex has been recoed
00395         myZmin1=vtxIter->position().z()-deltaZ1WithVertex_;
00396         myZmax1=vtxIter->position().z()+deltaZ1WithVertex_;
00397        }
00398 
00399       myMatchEle->set1stLayerZRange(myZmin1,myZmax1);
00400       myMatchPos->set1stLayerZRange(myZmin1,myZmax1);
00401 
00402       if (!fromTrackerSeeds_)
00403        {
00404         // try electron
00405         std::vector<std::pair<RecHitWithDist,ConstRecHitPointer> > elePixelHits
00406           = myMatchEle->compatibleHits(clusterPos,vertexPos,clusterEnergy,-1.,tTopo) ;
00407         seedsFromRecHits(elePixelHits,dir,vertexPos,caloCluster,out,false) ;
00408         // try positron
00409               std::vector<std::pair<RecHitWithDist,ConstRecHitPointer> > posPixelHits
00410                 = myMatchPos->compatibleHits(clusterPos,vertexPos,clusterEnergy,1.,tTopo) ;
00411         seedsFromRecHits(posPixelHits,dir,vertexPos,caloCluster,out,true) ;
00412        }
00413       else
00414        {
00415         // try electron
00416         std::vector<SeedWithInfo> elePixelSeeds
00417          = myMatchEle->compatibleSeeds(theInitialSeedColl,clusterPos,vertexPos,clusterEnergy,-1.) ;
00418         seedsFromTrajectorySeeds(elePixelSeeds,caloCluster,hoe1,hoe2,out,false) ;
00419         // try positron
00420         std::vector<SeedWithInfo> posPixelSeeds
00421          = myMatchPos->compatibleSeeds(theInitialSeedColl,clusterPos,vertexPos,clusterEnergy,1.) ;
00422         seedsFromTrajectorySeeds(posPixelSeeds,caloCluster,hoe1,hoe2,out,true) ;
00423        }
00424      }
00425    }
00426 
00427   return ;
00428  }
00429 
00430 void ElectronSeedGenerator::seedsFromRecHits
00431  ( std::vector<std::pair<RecHitWithDist,ConstRecHitPointer> > & pixelHits,
00432    PropagationDirection & dir,
00433    const GlobalPoint & vertexPos, const reco::ElectronSeed::CaloClusterRef & cluster,
00434    reco::ElectronSeedCollection & out,
00435    bool positron )
00436  {
00437   if (!pixelHits.empty())
00438    { LogDebug("ElectronSeedGenerator") << "Compatible "<<(positron?"positron":"electron")<<" hits found." ; }
00439 
00440   std::vector<std::pair<RecHitWithDist,ConstRecHitPointer> >::iterator v ;
00441   for ( v = pixelHits.begin() ; v != pixelHits.end() ; v++ )
00442    {
00443     if (!positron)
00444      { (*v).first.invert() ; }
00445     if (!prepareElTrackSeed((*v).first.recHit(),(*v).second,vertexPos))
00446      { continue ; }
00447     reco::ElectronSeed seed(pts_,recHits_,dir) ;
00448     seed.setCaloCluster(cluster) ;
00449     addSeed(seed,0,positron,out) ;
00450    }
00451  }
00452 
00453 void ElectronSeedGenerator::seedsFromTrajectorySeeds
00454  ( const std::vector<SeedWithInfo> & pixelSeeds,
00455    const reco::ElectronSeed::CaloClusterRef & cluster,
00456    float hoe1, float hoe2,
00457    reco::ElectronSeedCollection & out,
00458    bool positron )
00459  {
00460   if (!pixelSeeds.empty())
00461    { LogDebug("ElectronSeedGenerator") << "Compatible "<<(positron?"positron":"electron")<<" seeds found." ; }
00462 
00463   std::vector<SeedWithInfo>::const_iterator s;
00464   for ( s = pixelSeeds.begin() ; s != pixelSeeds.end() ; s++ )
00465    {
00466     reco::ElectronSeed seed(s->seed()) ;
00467     seed.setCaloCluster(cluster,s->hitsMask(),s->subDet2(),s->subDet1(),hoe1,hoe2) ;
00468     addSeed(seed,&*s,positron,out) ;
00469    }
00470  }
00471 
00472 void ElectronSeedGenerator::addSeed
00473  ( reco::ElectronSeed & seed,
00474    const SeedWithInfo * info,
00475    bool positron,
00476    reco::ElectronSeedCollection & out )
00477  {
00478   if (!info)
00479    { out.push_back(seed) ; return ; }
00480 
00481   if (positron)
00482    { seed.setPosAttributes(info->dRz2(),info->dPhi2(),info->dRz1(),info->dPhi1()) ; }
00483   else
00484    { seed.setNegAttributes(info->dRz2(),info->dPhi2(),info->dRz1(),info->dPhi1()) ; }
00485   reco::ElectronSeedCollection::iterator resItr ;
00486   for ( resItr=out.begin() ; resItr!=out.end() ; ++resItr )
00487    {
00488     if ( (seed.caloCluster()==resItr->caloCluster()) &&
00489          (seed.hitsMask()==resItr->hitsMask()) &&
00490          equivalent(seed,*resItr) )
00491      {
00492       if (positron)
00493        {
00494         if ( resItr->dRz2Pos()==std::numeric_limits<float>::infinity() &&
00495              resItr->dRz2()!=std::numeric_limits<float>::infinity() )
00496          {
00497           resItr->setPosAttributes(info->dRz2(),info->dPhi2(),info->dRz1(),info->dPhi1()) ;
00498           seed.setNegAttributes(resItr->dRz2(),resItr->dPhi2(),resItr->dRz1(),resItr->dPhi1()) ;
00499           break ;
00500          }
00501         else
00502          {
00503           if ( resItr->dRz2Pos()!=std::numeric_limits<float>::infinity() )
00504            {
00505             if ( resItr->dRz2Pos()!=seed.dRz2Pos() )
00506              {
00507               edm::LogWarning("ElectronSeedGenerator|BadValue")
00508                <<"this similar old seed already has another dRz2Pos"
00509                <<"\nold seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)resItr->hitsMask()<<"/"<<resItr->dRz2()<<"/"<<resItr->dPhi2()<<"/"<<resItr->dRz2Pos()<<"/"<<resItr->dPhi2Pos()
00510                <<"\nnew seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)seed.hitsMask()<<"/"<<seed.dRz2()<<"/"<<seed.dPhi2()<<"/"<<seed.dRz2Pos()<<"/"<<seed.dPhi2Pos() ;
00511              }
00512 //            else
00513 //             {
00514 //              edm::LogWarning("ElectronSeedGenerator|UnexpectedValue")
00515 //               <<"this old seed already knows its dRz2Pos, we suspect duplicates in input trajectry seeds"
00516 //               <<"\nold seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)resItr->hitsMask()<<"/"<<resItr->dRz2()<<"/"<<resItr->dPhi2()<<"/"<<resItr->dRz2Pos()<<"/"<<resItr->dPhi2Pos()
00517 //               <<"\nnew seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)seed.hitsMask()<<"/"<<seed.dRz2()<<"/"<<seed.dPhi2()<<"/"<<seed.dRz2Pos()<<"/"<<seed.dPhi2Pos() ;
00518 //             }
00519             }
00520 //          if (resItr->dRz2()==std::numeric_limits<float>::infinity())
00521 //           {
00522 //            edm::LogWarning("ElectronSeedGenerator|BadValue")
00523 //             <<"this old seed has no dRz2, we suspect duplicates in input trajectry seeds"
00524 //             <<"\nold seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)resItr->hitsMask()<<"/"<<resItr->dRz2()<<"/"<<resItr->dPhi2()<<"/"<<resItr->dRz2Pos()<<"/"<<resItr->dPhi2Pos()
00525 //             <<"\nnew seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)seed.hitsMask()<<"/"<<seed.dRz2()<<"/"<<seed.dPhi2()<<"/"<<seed.dRz2Pos()<<"/"<<seed.dPhi2Pos() ;
00526 //           }
00527          }
00528        }
00529       else
00530        {
00531         if ( resItr->dRz2()==std::numeric_limits<float>::infinity()
00532           && resItr->dRz2Pos()!=std::numeric_limits<float>::infinity() )
00533          {
00534           resItr->setNegAttributes(info->dRz2(),info->dPhi2(),info->dRz1(),info->dPhi1()) ;
00535           seed.setPosAttributes(resItr->dRz2Pos(),resItr->dPhi2Pos(),resItr->dRz1Pos(),resItr->dPhi1Pos()) ;
00536           break ;
00537          }
00538         else
00539          {
00540           if ( resItr->dRz2()!=std::numeric_limits<float>::infinity() )
00541            {
00542             if (resItr->dRz2()!=seed.dRz2())
00543              {
00544               edm::LogWarning("ElectronSeedGenerator|BadValue")
00545                <<"this old seed already has another dRz2"
00546                <<"\nold seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)resItr->hitsMask()<<"/"<<resItr->dRz2()<<"/"<<resItr->dPhi2()<<"/"<<resItr->dRz2Pos()<<"/"<<resItr->dPhi2Pos()
00547                <<"\nnew seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)seed.hitsMask()<<"/"<<seed.dRz2()<<"/"<<seed.dPhi2()<<"/"<<seed.dRz2Pos()<<"/"<<seed.dPhi2Pos() ;
00548              }
00549     //        else
00550     //         {
00551     //          edm::LogWarning("ElectronSeedGenerator|UnexpectedValue")
00552     //           <<"this old seed already knows its dRz2, we suspect duplicates in input trajectry seeds"
00553     //           <<"\nold seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)resItr->hitsMask()<<"/"<<resItr->dRz2()<<"/"<<resItr->dPhi2()<<"/"<<resItr->dRz2Pos()<<"/"<<resItr->dPhi2Pos()
00554     //           <<"\nnew seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)seed.hitsMask()<<"/"<<seed.dRz2()<<"/"<<seed.dPhi2()<<"/"<<seed.dRz2Pos()<<"/"<<seed.dPhi2Pos() ;
00555     //          seed.setPosAttributes(resItr->dRz2Pos(),resItr->dPhi2Pos(),resItr->dRz1Pos(),resItr->dPhi1Pos()) ;
00556     //         }
00557            }
00558 //          if (resItr->dRz2Pos()==std::numeric_limits<float>::infinity())
00559 //           {
00560 //            edm::LogWarning("ElectronSeedGenerator|BadValue")
00561 //             <<"this old seed has no dRz2Pos"
00562 //             <<"\nold seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)resItr->hitsMask()<<"/"<<resItr->dRz2()<<"/"<<resItr->dPhi2()<<"/"<<resItr->dRz2Pos()<<"/"<<resItr->dPhi2Pos()
00563 //             <<"\nnew seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)seed.hitsMask()<<"/"<<seed.dRz2()<<"/"<<seed.dPhi2()<<"/"<<seed.dRz2Pos()<<"/"<<seed.dPhi2Pos() ;
00564 //           }
00565          }
00566        }
00567      }
00568    }
00569 
00570   out.push_back(seed) ;
00571  }
00572 
00573 bool ElectronSeedGenerator::prepareElTrackSeed
00574  ( ConstRecHitPointer innerhit,
00575    ConstRecHitPointer outerhit,
00576    const GlobalPoint& vertexPos )
00577  {
00578 
00579   // debug prints
00580   LogDebug("") <<"[ElectronSeedGenerator::prepareElTrackSeed] inner PixelHit   x,y,z "<<innerhit->globalPosition();
00581   LogDebug("") <<"[ElectronSeedGenerator::prepareElTrackSeed] outer PixelHit   x,y,z "<<outerhit->globalPosition();
00582 
00583   recHits_.clear();
00584 
00585   SiPixelRecHit *pixhit=0;
00586   SiStripMatchedRecHit2D *striphit=0;
00587   const SiPixelRecHit* constpixhit = dynamic_cast <const SiPixelRecHit*> (innerhit->hit());
00588   if (constpixhit) {
00589     pixhit=new SiPixelRecHit(*constpixhit);
00590     recHits_.push_back(pixhit);
00591   } else  return false;
00592   constpixhit =  dynamic_cast <const SiPixelRecHit *> (outerhit->hit());
00593   if (constpixhit) {
00594     pixhit=new SiPixelRecHit(*constpixhit);
00595     recHits_.push_back(pixhit);
00596   } else {
00597     const SiStripMatchedRecHit2D * conststriphit=dynamic_cast <const SiStripMatchedRecHit2D *> (outerhit->hit());
00598     if (conststriphit) {
00599       striphit = new SiStripMatchedRecHit2D(*conststriphit);
00600       recHits_.push_back(striphit);
00601     } else return false;
00602   }
00603 
00604   typedef TrajectoryStateOnSurface     TSOS;
00605 
00606   // FIXME to be optimized outside the loop
00607   edm::ESHandle<MagneticField> bfield;
00608   theSetup->get<IdealMagneticFieldRecord>().get(bfield);
00609   float nomField = bfield->nominalValue();
00610 
00611   // make a spiral
00612   FastHelix helix(outerhit->globalPosition(),innerhit->globalPosition(),vertexPos,nomField,&*bfield);
00613   if ( !helix.isValid()) {
00614     return false;
00615   }
00616   FreeTrajectoryState fts(helix.stateAtVertex());
00617   TSOS propagatedState = thePropagator->propagate(fts,innerhit->det()->surface()) ;
00618   if (!propagatedState.isValid())
00619     return false;
00620   TSOS updatedState = theUpdator->update(propagatedState, *innerhit);
00621 
00622   TSOS propagatedState_out = thePropagator->propagate(updatedState,outerhit->det()->surface()) ;
00623   if (!propagatedState_out.isValid())
00624     return false;
00625   TSOS updatedState_out = theUpdator->update(propagatedState_out, *outerhit);
00626   // debug prints
00627   LogDebug("") <<"[ElectronSeedGenerator::prepareElTrackSeed] final TSOS, position: "<<updatedState_out.globalPosition()<<" momentum: "<<updatedState_out.globalMomentum();
00628   LogDebug("") <<"[ElectronSeedGenerator::prepareElTrackSeed] final TSOS Pt: "<<updatedState_out.globalMomentum().perp();
00629   pts_ =  trajectoryStateTransform::persistentState(updatedState_out, outerhit->geographicalId().rawId());
00630 
00631   return true;
00632  }
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/RecoEgamma/EgammaElectronAlgos/src/ElectronSeedGenerator.cc
