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