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 
00252   theSetup= &setup;
00253   NavigationSetter theSetter(*theNavigationSchool);
00254 
00255   // get initial TrajectorySeeds if necessary
00256   //  if (fromTrackerSeeds_) e.getByLabel(initialSeeds_, theInitialSeedColl);
00257 
00258   // get the beamspot from the Event:
00259   //e.getByType(theBeamSpot);
00260   e.getByLabel(beamSpotTag_,theBeamSpot);
00261 
00262   // if required get the vertices
00263   if (useRecoVertex_) e.getByLabel(verticesTag_,theVertices);
00264 
00265   if (!fromTrackerSeeds_)
00266    { theMeasurementTracker->update(e) ; }
00267 
00268   for  (unsigned int i=0;i<sclRefs.size();++i) {
00269     // Find the seeds
00270     recHits_.clear();
00271 
00272     LogDebug ("run") << "new cluster, calling seedsFromThisCluster";
00273     seedsFromThisCluster(sclRefs[i],hoe1s[i],hoe2s[i],out);
00274   }
00275 
00276   LogDebug ("run") << ": For event "<<e.id();
00277   LogDebug ("run") <<"Nr of superclusters after filter: "<<sclRefs.size()
00278    <<", no. of ElectronSeeds found  = " << out.size();
00279 }
00280 
00281 void ElectronSeedGenerator::seedsFromThisCluster
00282 ( edm::Ref<reco::SuperClusterCollection> seedCluster,
00283   float hoe1, float hoe2,
00284   reco::ElectronSeedCollection & out )
00285 {
00286   float clusterEnergy = seedCluster->energy() ;
00287   GlobalPoint clusterPos
00288     ( seedCluster->position().x(),
00289       seedCluster->position().y(),
00290       seedCluster->position().z() ) ;
00291   reco::ElectronSeed::CaloClusterRef caloCluster(seedCluster) ;
00292 
00293   //LogDebug("") << "[ElectronSeedGenerator::seedsFromThisCluster] new supercluster with energy: " << clusterEnergy ;
00294   //LogDebug("") << "[ElectronSeedGenerator::seedsFromThisCluster] and position: " << clusterPos ;
00295 
00296   if (dynamicphiroad_)
00297    {
00298     float clusterEnergyT = clusterEnergy / cosh( EleRelPoint(clusterPos,theBeamSpot->position()).eta() ) ;
00299 
00300     float deltaPhi1 ;
00301     if (clusterEnergyT < lowPtThreshold_)
00302      { deltaPhi1= deltaPhi1Low_ ; }
00303     else if (clusterEnergyT > highPtThreshold_)
00304      { deltaPhi1= deltaPhi1High_ ; }
00305     else
00306      { deltaPhi1 = deltaPhi1Coef1_ + deltaPhi1Coef2_/clusterEnergyT ; }
00307 
00308     float ephimin1 = -deltaPhi1*sizeWindowENeg_ ;
00309     float ephimax1 =  deltaPhi1*(1.-sizeWindowENeg_);
00310     float pphimin1 = -deltaPhi1*(1.-sizeWindowENeg_);
00311     float pphimax1 =  deltaPhi1*sizeWindowENeg_;
00312 
00313     float phimin2B  = -deltaPhi2B_/2. ;
00314     float phimax2B  =  deltaPhi2B_/2. ;
00315     float phimin2F  = -deltaPhi2F_/2. ;
00316     float phimax2F  =  deltaPhi2F_/2. ;
00317 
00318 
00319     myMatchEle->set1stLayer(ephimin1,ephimax1);
00320     myMatchPos->set1stLayer(pphimin1,pphimax1);
00321     myMatchEle->set2ndLayer(phimin2B,phimax2B, phimin2F,phimax2F);
00322     myMatchPos->set2ndLayer(phimin2B,phimax2B, phimin2F,phimax2F);
00323    }
00324 
00325   PropagationDirection dir = alongMomentum ;
00326 
00327   if (!useRecoVertex_) // here use the beam spot position
00328    {
00329     double sigmaZ=theBeamSpot->sigmaZ();
00330     double sigmaZ0Error=theBeamSpot->sigmaZ0Error();
00331     double sq=sqrt(sigmaZ*sigmaZ+sigmaZ0Error*sigmaZ0Error);
00332     double myZmin1=theBeamSpot->position().z()-nSigmasDeltaZ1_*sq;
00333     double myZmax1=theBeamSpot->position().z()+nSigmasDeltaZ1_*sq;
00334 
00335     GlobalPoint vertexPos ;
00336     ele_convert(theBeamSpot->position(),vertexPos) ;
00337 
00338     myMatchEle->set1stLayerZRange(myZmin1,myZmax1);
00339     myMatchPos->set1stLayerZRange(myZmin1,myZmax1);
00340 
00341     if (!fromTrackerSeeds_)
00342      {
00343       // try electron
00344       std::vector<std::pair<RecHitWithDist,ConstRecHitPointer> > elePixelHits
00345        = myMatchEle->compatibleHits(clusterPos,vertexPos,clusterEnergy,-1.) ;
00346       GlobalPoint eleVertex(theBeamSpot->position().x(),theBeamSpot->position().y(),myMatchEle->getVertex()) ;
00347       seedsFromRecHits(elePixelHits,dir,eleVertex,caloCluster,out,false) ;
00348       // try positron
00349       std::vector<std::pair<RecHitWithDist,ConstRecHitPointer> > posPixelHits
00350        = myMatchPos->compatibleHits(clusterPos,vertexPos,clusterEnergy,1.) ;
00351       GlobalPoint posVertex(theBeamSpot->position().x(),theBeamSpot->position().y(),myMatchPos->getVertex()) ;
00352       seedsFromRecHits(posPixelHits,dir,posVertex,caloCluster,out,true) ;
00353      }
00354     else
00355      {
00356       // try electron
00357       std::vector<SeedWithInfo> elePixelSeeds
00358        = myMatchEle->compatibleSeeds(theInitialSeedColl,clusterPos,vertexPos,clusterEnergy,-1.) ;
00359       seedsFromTrajectorySeeds(elePixelSeeds,caloCluster,hoe1,hoe2,out,false) ;
00360       // try positron
00361       std::vector<SeedWithInfo> posPixelSeeds
00362        = myMatchPos->compatibleSeeds(theInitialSeedColl,clusterPos,vertexPos,clusterEnergy,1.) ;
00363       seedsFromTrajectorySeeds(posPixelSeeds,caloCluster,hoe1,hoe2,out,true) ;
00364      }
00365 
00366    }
00367   else // here we use the reco vertices
00368    {
00369 
00370     myMatchEle->setUseRecoVertex(true) ; //Hit matchers need to know that the vertex is known
00371     myMatchPos->setUseRecoVertex(true) ;
00372 
00373     const  std::vector<reco::Vertex> * vtxCollection = theVertices.product() ;
00374     std::vector<reco::Vertex>::const_iterator vtxIter ;
00375     for (vtxIter = vtxCollection->begin(); vtxIter != vtxCollection->end() ; vtxIter++)
00376      {
00377 
00378       GlobalPoint vertexPos(vtxIter->position().x(),vtxIter->position().y(),vtxIter->position().z());
00379       double myZmin1, myZmax1 ;
00380       if (vertexPos.z()==theBeamSpot->position().z())
00381        { // in case vetex not found
00382         double sigmaZ=theBeamSpot->sigmaZ();
00383         double sigmaZ0Error=theBeamSpot->sigmaZ0Error();
00384         double sq=sqrt(sigmaZ*sigmaZ+sigmaZ0Error*sigmaZ0Error);
00385         myZmin1=theBeamSpot->position().z()-nSigmasDeltaZ1_*sq;
00386         myZmax1=theBeamSpot->position().z()+nSigmasDeltaZ1_*sq;
00387        }
00388       else
00389        { // a vertex has been recoed
00390         myZmin1=vtxIter->position().z()-deltaZ1WithVertex_;
00391         myZmax1=vtxIter->position().z()+deltaZ1WithVertex_;
00392        }
00393 
00394       myMatchEle->set1stLayerZRange(myZmin1,myZmax1);
00395       myMatchPos->set1stLayerZRange(myZmin1,myZmax1);
00396 
00397       if (!fromTrackerSeeds_)
00398        {
00399         // try electron
00400         std::vector<std::pair<RecHitWithDist,ConstRecHitPointer> > elePixelHits
00401          = myMatchEle->compatibleHits(clusterPos,vertexPos,clusterEnergy,-1.) ;
00402         seedsFromRecHits(elePixelHits,dir,vertexPos,caloCluster,out,false) ;
00403         // try positron
00404               std::vector<std::pair<RecHitWithDist,ConstRecHitPointer> > posPixelHits
00405                = myMatchPos->compatibleHits(clusterPos,vertexPos,clusterEnergy,1.) ;
00406         seedsFromRecHits(posPixelHits,dir,vertexPos,caloCluster,out,true) ;
00407        }
00408       else
00409        {
00410         // try electron
00411         std::vector<SeedWithInfo> elePixelSeeds
00412          = myMatchEle->compatibleSeeds(theInitialSeedColl,clusterPos,vertexPos,clusterEnergy,-1.) ;
00413         seedsFromTrajectorySeeds(elePixelSeeds,caloCluster,hoe1,hoe2,out,false) ;
00414         // try positron
00415         std::vector<SeedWithInfo> posPixelSeeds
00416          = myMatchPos->compatibleSeeds(theInitialSeedColl,clusterPos,vertexPos,clusterEnergy,1.) ;
00417         seedsFromTrajectorySeeds(posPixelSeeds,caloCluster,hoe1,hoe2,out,true) ;
00418        }
00419      }
00420    }
00421 
00422   return ;
00423  }
00424 
00425 void ElectronSeedGenerator::seedsFromRecHits
00426  ( std::vector<std::pair<RecHitWithDist,ConstRecHitPointer> > & pixelHits,
00427    PropagationDirection & dir,
00428    const GlobalPoint & vertexPos, const reco::ElectronSeed::CaloClusterRef & cluster,
00429    reco::ElectronSeedCollection & out,
00430    bool positron )
00431  {
00432   if (!pixelHits.empty())
00433    { LogDebug("ElectronSeedGenerator") << "Compatible "<<(positron?"positron":"electron")<<" hits found." ; }
00434 
00435   std::vector<std::pair<RecHitWithDist,ConstRecHitPointer> >::iterator v ;
00436   for ( v = pixelHits.begin() ; v != pixelHits.end() ; v++ )
00437    {
00438     if (!positron)
00439      { (*v).first.invert() ; }
00440     if (!prepareElTrackSeed((*v).first.recHit(),(*v).second,vertexPos))
00441      { continue ; }
00442     reco::ElectronSeed seed(pts_,recHits_,dir) ;
00443     seed.setCaloCluster(cluster) ;
00444     addSeed(seed,0,positron,out) ;
00445    }
00446  }
00447 
00448 void ElectronSeedGenerator::seedsFromTrajectorySeeds
00449  ( const std::vector<SeedWithInfo> & pixelSeeds,
00450    const reco::ElectronSeed::CaloClusterRef & cluster,
00451    float hoe1, float hoe2,
00452    reco::ElectronSeedCollection & out,
00453    bool positron )
00454  {
00455   if (!pixelSeeds.empty())
00456    { LogDebug("ElectronSeedGenerator") << "Compatible "<<(positron?"positron":"electron")<<" seeds found." ; }
00457 
00458   std::vector<SeedWithInfo>::const_iterator s;
00459   for ( s = pixelSeeds.begin() ; s != pixelSeeds.end() ; s++ )
00460    {
00461     reco::ElectronSeed seed(s->seed()) ;
00462     seed.setCaloCluster(cluster,s->hitsMask(),s->subDet2(),s->subDet1(),hoe1,hoe2) ;
00463     addSeed(seed,&*s,positron,out) ;
00464    }
00465  }
00466 
00467 void ElectronSeedGenerator::addSeed
00468  ( reco::ElectronSeed & seed,
00469    const SeedWithInfo * info,
00470    bool positron,
00471    reco::ElectronSeedCollection & out )
00472  {
00473   if (!info)
00474    { out.push_back(seed) ; return ; }
00475 
00476   if (positron)
00477    { seed.setPosAttributes(info->dRz2(),info->dPhi2(),info->dRz1(),info->dPhi1()) ; }
00478   else
00479    { seed.setNegAttributes(info->dRz2(),info->dPhi2(),info->dRz1(),info->dPhi1()) ; }
00480   reco::ElectronSeedCollection::iterator resItr ;
00481   for ( resItr=out.begin() ; resItr!=out.end() ; ++resItr )
00482    {
00483     if ( (seed.caloCluster()==resItr->caloCluster()) &&
00484          (seed.hitsMask()==resItr->hitsMask()) &&
00485          equivalent(seed,*resItr) )
00486      {
00487       if (positron)
00488        {
00489         if ( resItr->dRz2Pos()==std::numeric_limits<float>::infinity() &&
00490              resItr->dRz2()!=std::numeric_limits<float>::infinity() )
00491          {
00492           resItr->setPosAttributes(info->dRz2(),info->dPhi2(),info->dRz1(),info->dPhi1()) ;
00493           seed.setNegAttributes(resItr->dRz2(),resItr->dPhi2(),resItr->dRz1(),resItr->dPhi1()) ;
00494           break ;
00495          }
00496         else
00497          {
00498           if ( resItr->dRz2Pos()!=std::numeric_limits<float>::infinity() )
00499            {
00500             if ( resItr->dRz2Pos()!=seed.dRz2Pos() )
00501              {
00502               edm::LogWarning("ElectronSeedGenerator|BadValue")
00503                <<"this similar old seed already has another dRz2Pos"
00504                <<"\nold seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)resItr->hitsMask()<<"/"<<resItr->dRz2()<<"/"<<resItr->dPhi2()<<"/"<<resItr->dRz2Pos()<<"/"<<resItr->dPhi2Pos()
00505                <<"\nnew seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)seed.hitsMask()<<"/"<<seed.dRz2()<<"/"<<seed.dPhi2()<<"/"<<seed.dRz2Pos()<<"/"<<seed.dPhi2Pos() ;
00506              }
00507 //            else
00508 //             {
00509 //              edm::LogWarning("ElectronSeedGenerator|UnexpectedValue")
00510 //               <<"this old seed already knows its dRz2Pos, we suspect duplicates in input trajectry seeds"
00511 //               <<"\nold seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)resItr->hitsMask()<<"/"<<resItr->dRz2()<<"/"<<resItr->dPhi2()<<"/"<<resItr->dRz2Pos()<<"/"<<resItr->dPhi2Pos()
00512 //               <<"\nnew seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)seed.hitsMask()<<"/"<<seed.dRz2()<<"/"<<seed.dPhi2()<<"/"<<seed.dRz2Pos()<<"/"<<seed.dPhi2Pos() ;
00513 //             }
00514             }
00515 //          if (resItr->dRz2()==std::numeric_limits<float>::infinity())
00516 //           {
00517 //            edm::LogWarning("ElectronSeedGenerator|BadValue")
00518 //             <<"this old seed has no dRz2, we suspect duplicates in input trajectry seeds"
00519 //             <<"\nold seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)resItr->hitsMask()<<"/"<<resItr->dRz2()<<"/"<<resItr->dPhi2()<<"/"<<resItr->dRz2Pos()<<"/"<<resItr->dPhi2Pos()
00520 //             <<"\nnew seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)seed.hitsMask()<<"/"<<seed.dRz2()<<"/"<<seed.dPhi2()<<"/"<<seed.dRz2Pos()<<"/"<<seed.dPhi2Pos() ;
00521 //           }
00522          }
00523        }
00524       else
00525        {
00526         if ( resItr->dRz2()==std::numeric_limits<float>::infinity()
00527           && resItr->dRz2Pos()!=std::numeric_limits<float>::infinity() )
00528          {
00529           resItr->setNegAttributes(info->dRz2(),info->dPhi2(),info->dRz1(),info->dPhi1()) ;
00530           seed.setPosAttributes(resItr->dRz2Pos(),resItr->dPhi2Pos(),resItr->dRz1Pos(),resItr->dPhi1Pos()) ;
00531           break ;
00532          }
00533         else
00534          {
00535           if ( resItr->dRz2()!=std::numeric_limits<float>::infinity() )
00536            {
00537             if (resItr->dRz2()!=seed.dRz2())
00538              {
00539               edm::LogWarning("ElectronSeedGenerator|BadValue")
00540                <<"this old seed already has another dRz2"
00541                <<"\nold seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)resItr->hitsMask()<<"/"<<resItr->dRz2()<<"/"<<resItr->dPhi2()<<"/"<<resItr->dRz2Pos()<<"/"<<resItr->dPhi2Pos()
00542                <<"\nnew seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)seed.hitsMask()<<"/"<<seed.dRz2()<<"/"<<seed.dPhi2()<<"/"<<seed.dRz2Pos()<<"/"<<seed.dPhi2Pos() ;
00543              }
00544     //        else
00545     //         {
00546     //          edm::LogWarning("ElectronSeedGenerator|UnexpectedValue")
00547     //           <<"this old seed already knows its dRz2, we suspect duplicates in input trajectry seeds"
00548     //           <<"\nold seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)resItr->hitsMask()<<"/"<<resItr->dRz2()<<"/"<<resItr->dPhi2()<<"/"<<resItr->dRz2Pos()<<"/"<<resItr->dPhi2Pos()
00549     //           <<"\nnew seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)seed.hitsMask()<<"/"<<seed.dRz2()<<"/"<<seed.dPhi2()<<"/"<<seed.dRz2Pos()<<"/"<<seed.dPhi2Pos() ;
00550     //          seed.setPosAttributes(resItr->dRz2Pos(),resItr->dPhi2Pos(),resItr->dRz1Pos(),resItr->dPhi1Pos()) ;
00551     //         }
00552            }
00553 //          if (resItr->dRz2Pos()==std::numeric_limits<float>::infinity())
00554 //           {
00555 //            edm::LogWarning("ElectronSeedGenerator|BadValue")
00556 //             <<"this old seed has no dRz2Pos"
00557 //             <<"\nold seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)resItr->hitsMask()<<"/"<<resItr->dRz2()<<"/"<<resItr->dPhi2()<<"/"<<resItr->dRz2Pos()<<"/"<<resItr->dPhi2Pos()
00558 //             <<"\nnew seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)seed.hitsMask()<<"/"<<seed.dRz2()<<"/"<<seed.dPhi2()<<"/"<<seed.dRz2Pos()<<"/"<<seed.dPhi2Pos() ;
00559 //           }
00560          }
00561        }
00562      }
00563    }
00564 
00565   out.push_back(seed) ;
00566  }
00567 
00568 bool ElectronSeedGenerator::prepareElTrackSeed
00569  ( ConstRecHitPointer innerhit,
00570    ConstRecHitPointer outerhit,
00571    const GlobalPoint& vertexPos )
00572  {
00573 
00574   // debug prints
00575   LogDebug("") <<"[ElectronSeedGenerator::prepareElTrackSeed] inner PixelHit   x,y,z "<<innerhit->globalPosition();
00576   LogDebug("") <<"[ElectronSeedGenerator::prepareElTrackSeed] outer PixelHit   x,y,z "<<outerhit->globalPosition();
00577 
00578   recHits_.clear();
00579 
00580   SiPixelRecHit *pixhit=0;
00581   SiStripMatchedRecHit2D *striphit=0;
00582   const SiPixelRecHit* constpixhit = dynamic_cast <const SiPixelRecHit*> (innerhit->hit());
00583   if (constpixhit) {
00584     pixhit=new SiPixelRecHit(*constpixhit);
00585     recHits_.push_back(pixhit);
00586   } else  return false;
00587   constpixhit =  dynamic_cast <const SiPixelRecHit *> (outerhit->hit());
00588   if (constpixhit) {
00589     pixhit=new SiPixelRecHit(*constpixhit);
00590     recHits_.push_back(pixhit);
00591   } else {
00592     const SiStripMatchedRecHit2D * conststriphit=dynamic_cast <const SiStripMatchedRecHit2D *> (outerhit->hit());
00593     if (conststriphit) {
00594       striphit = new SiStripMatchedRecHit2D(*conststriphit);
00595       recHits_.push_back(striphit);
00596     } else return false;
00597   }
00598 
00599   typedef TrajectoryStateOnSurface     TSOS;
00600 
00601   // make a spiral
00602   FastHelix helix(outerhit->globalPosition(),innerhit->globalPosition(),vertexPos,*theSetup);
00603   if ( !helix.isValid()) {
00604     return false;
00605   }
00606   FreeTrajectoryState fts = helix.stateAtVertex();
00607   TSOS propagatedState = thePropagator->propagate(fts,innerhit->det()->surface()) ;
00608   if (!propagatedState.isValid())
00609     return false;
00610   TSOS updatedState = theUpdator->update(propagatedState, *innerhit);
00611 
00612   TSOS propagatedState_out = thePropagator->propagate(updatedState,outerhit->det()->surface()) ;
00613   if (!propagatedState_out.isValid())
00614     return false;
00615   TSOS updatedState_out = theUpdator->update(propagatedState_out, *outerhit);
00616   // debug prints
00617   LogDebug("") <<"[ElectronSeedGenerator::prepareElTrackSeed] final TSOS, position: "<<updatedState_out.globalPosition()<<" momentum: "<<updatedState_out.globalMomentum();
00618   LogDebug("") <<"[ElectronSeedGenerator::prepareElTrackSeed] final TSOS Pt: "<<updatedState_out.globalMomentum().perp();
00619   pts_ =  trajectoryStateTransform::persistentState(updatedState_out, outerhit->geographicalId().rawId());
00620 
00621   return true;
00622  }
/data/refman/pasoursint/CMSSW_6_1_1/src/RecoEgamma/EgammaElectronAlgos/src/ElectronSeedGenerator.cc
