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), pts_(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     delete pts_;
00446     pts_=0;
00447    }
00448  }
00449 
00450 void ElectronSeedGenerator::seedsFromTrajectorySeeds
00451  ( const std::vector<SeedWithInfo> & pixelSeeds,
00452    const reco::ElectronSeed::CaloClusterRef & cluster,
00453    float hoe1, float hoe2,
00454    reco::ElectronSeedCollection & out,
00455    bool positron )
00456  {
00457   if (!pixelSeeds.empty())
00458    { LogDebug("ElectronSeedGenerator") << "Compatible "<<(positron?"positron":"electron")<<" seeds found." ; }
00459 
00460   std::vector<SeedWithInfo>::const_iterator s;
00461   for ( s = pixelSeeds.begin() ; s != pixelSeeds.end() ; s++ )
00462    {
00463     reco::ElectronSeed seed(s->seed()) ;
00464     seed.setCaloCluster(cluster,s->hitsMask(),s->subDet2(),s->subDet1(),hoe1,hoe2) ;
00465     addSeed(seed,&*s,positron,out) ;
00466    }
00467  }
00468 
00469 void ElectronSeedGenerator::addSeed
00470  ( reco::ElectronSeed & seed,
00471    const SeedWithInfo * info,
00472    bool positron,
00473    reco::ElectronSeedCollection & out )
00474  {
00475   if (!info)
00476    { out.push_back(seed) ; return ; }
00477 
00478   if (positron)
00479    { seed.setPosAttributes(info->dRz2(),info->dPhi2(),info->dRz1(),info->dPhi1()) ; }
00480   else
00481    { seed.setNegAttributes(info->dRz2(),info->dPhi2(),info->dRz1(),info->dPhi1()) ; }
00482   reco::ElectronSeedCollection::iterator resItr ;
00483   for ( resItr=out.begin() ; resItr!=out.end() ; ++resItr )
00484    {
00485     if ( (seed.caloCluster()==resItr->caloCluster()) &&
00486          (seed.hitsMask()==resItr->hitsMask()) &&
00487          equivalent(seed,*resItr) )
00488      {
00489       if (positron)
00490        {
00491         if ( resItr->dRz2Pos()==std::numeric_limits<float>::infinity() &&
00492              resItr->dRz2()!=std::numeric_limits<float>::infinity() )
00493          {
00494           resItr->setPosAttributes(info->dRz2(),info->dPhi2(),info->dRz1(),info->dPhi1()) ;
00495           seed.setNegAttributes(resItr->dRz2(),resItr->dPhi2(),resItr->dRz1(),resItr->dPhi1()) ;
00496           break ;
00497          }
00498         else
00499          {
00500           if ( resItr->dRz2Pos()!=std::numeric_limits<float>::infinity() )
00501            {
00502             if ( resItr->dRz2Pos()!=seed.dRz2Pos() )
00503              {
00504               edm::LogWarning("ElectronSeedGenerator|BadValue")
00505                <<"this similar old seed already has another dRz2Pos"
00506                <<"\nold seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)resItr->hitsMask()<<"/"<<resItr->dRz2()<<"/"<<resItr->dPhi2()<<"/"<<resItr->dRz2Pos()<<"/"<<resItr->dPhi2Pos()
00507                <<"\nnew seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)seed.hitsMask()<<"/"<<seed.dRz2()<<"/"<<seed.dPhi2()<<"/"<<seed.dRz2Pos()<<"/"<<seed.dPhi2Pos() ;
00508              }
00509 //            else
00510 //             {
00511 //              edm::LogWarning("ElectronSeedGenerator|UnexpectedValue")
00512 //               <<"this old seed already knows its dRz2Pos, we suspect duplicates in input trajectry seeds"
00513 //               <<"\nold seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)resItr->hitsMask()<<"/"<<resItr->dRz2()<<"/"<<resItr->dPhi2()<<"/"<<resItr->dRz2Pos()<<"/"<<resItr->dPhi2Pos()
00514 //               <<"\nnew seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)seed.hitsMask()<<"/"<<seed.dRz2()<<"/"<<seed.dPhi2()<<"/"<<seed.dRz2Pos()<<"/"<<seed.dPhi2Pos() ;
00515 //             }
00516             }
00517 //          if (resItr->dRz2()==std::numeric_limits<float>::infinity())
00518 //           {
00519 //            edm::LogWarning("ElectronSeedGenerator|BadValue")
00520 //             <<"this old seed has no dRz2, 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        }
00526       else
00527        {
00528         if ( resItr->dRz2()==std::numeric_limits<float>::infinity()
00529           && resItr->dRz2Pos()!=std::numeric_limits<float>::infinity() )
00530          {
00531           resItr->setNegAttributes(info->dRz2(),info->dPhi2(),info->dRz1(),info->dPhi1()) ;
00532           seed.setPosAttributes(resItr->dRz2Pos(),resItr->dPhi2Pos(),resItr->dRz1Pos(),resItr->dPhi1Pos()) ;
00533           break ;
00534          }
00535         else
00536          {
00537           if ( resItr->dRz2()!=std::numeric_limits<float>::infinity() )
00538            {
00539             if (resItr->dRz2()!=seed.dRz2())
00540              {
00541               edm::LogWarning("ElectronSeedGenerator|BadValue")
00542                <<"this old seed already has another dRz2"
00543                <<"\nold seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)resItr->hitsMask()<<"/"<<resItr->dRz2()<<"/"<<resItr->dPhi2()<<"/"<<resItr->dRz2Pos()<<"/"<<resItr->dPhi2Pos()
00544                <<"\nnew seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)seed.hitsMask()<<"/"<<seed.dRz2()<<"/"<<seed.dPhi2()<<"/"<<seed.dRz2Pos()<<"/"<<seed.dPhi2Pos() ;
00545              }
00546     //        else
00547     //         {
00548     //          edm::LogWarning("ElectronSeedGenerator|UnexpectedValue")
00549     //           <<"this old seed already knows its dRz2, we suspect duplicates in input trajectry seeds"
00550     //           <<"\nold seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)resItr->hitsMask()<<"/"<<resItr->dRz2()<<"/"<<resItr->dPhi2()<<"/"<<resItr->dRz2Pos()<<"/"<<resItr->dPhi2Pos()
00551     //           <<"\nnew seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)seed.hitsMask()<<"/"<<seed.dRz2()<<"/"<<seed.dPhi2()<<"/"<<seed.dRz2Pos()<<"/"<<seed.dPhi2Pos() ;
00552     //          seed.setPosAttributes(resItr->dRz2Pos(),resItr->dPhi2Pos(),resItr->dRz1Pos(),resItr->dPhi1Pos()) ;
00553     //         }
00554            }
00555 //          if (resItr->dRz2Pos()==std::numeric_limits<float>::infinity())
00556 //           {
00557 //            edm::LogWarning("ElectronSeedGenerator|BadValue")
00558 //             <<"this old seed has no dRz2Pos"
00559 //             <<"\nold seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)resItr->hitsMask()<<"/"<<resItr->dRz2()<<"/"<<resItr->dPhi2()<<"/"<<resItr->dRz2Pos()<<"/"<<resItr->dPhi2Pos()
00560 //             <<"\nnew seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: "<<(unsigned int)seed.hitsMask()<<"/"<<seed.dRz2()<<"/"<<seed.dPhi2()<<"/"<<seed.dRz2Pos()<<"/"<<seed.dPhi2Pos() ;
00561 //           }
00562          }
00563        }
00564      }
00565    }
00566 
00567   out.push_back(seed) ;
00568  }
00569 
00570 bool ElectronSeedGenerator::prepareElTrackSeed
00571  ( ConstRecHitPointer innerhit,
00572    ConstRecHitPointer outerhit,
00573    const GlobalPoint& vertexPos )
00574  {
00575 
00576   // debug prints
00577   LogDebug("") <<"[ElectronSeedGenerator::prepareElTrackSeed] inner PixelHit   x,y,z "<<innerhit->globalPosition();
00578   LogDebug("") <<"[ElectronSeedGenerator::prepareElTrackSeed] outer PixelHit   x,y,z "<<outerhit->globalPosition();
00579 
00580   pts_=0;
00581   recHits_.clear();
00582 
00583   SiPixelRecHit *pixhit=0;
00584   SiStripMatchedRecHit2D *striphit=0;
00585   const SiPixelRecHit* constpixhit = dynamic_cast <const SiPixelRecHit*> (innerhit->hit());
00586   if (constpixhit) {
00587     pixhit=new SiPixelRecHit(*constpixhit);
00588     recHits_.push_back(pixhit);
00589   } else  return false;
00590   constpixhit =  dynamic_cast <const SiPixelRecHit *> (outerhit->hit());
00591   if (constpixhit) {
00592     pixhit=new SiPixelRecHit(*constpixhit);
00593     recHits_.push_back(pixhit);
00594   } else {
00595     const SiStripMatchedRecHit2D * conststriphit=dynamic_cast <const SiStripMatchedRecHit2D *> (outerhit->hit());
00596     if (conststriphit) {
00597       striphit = new SiStripMatchedRecHit2D(*conststriphit);
00598       recHits_.push_back(striphit);
00599     } else return false;
00600   }
00601 
00602   typedef TrajectoryStateOnSurface     TSOS;
00603 
00604   // make a spiral
00605   FastHelix helix(outerhit->globalPosition(),innerhit->globalPosition(),vertexPos,*theSetup);
00606   if ( !helix.isValid()) {
00607     return false;
00608   }
00609   FreeTrajectoryState fts = helix.stateAtVertex();
00610   TSOS propagatedState = thePropagator->propagate(fts,innerhit->det()->surface()) ;
00611   if (!propagatedState.isValid())
00612     return false;
00613   TSOS updatedState = theUpdator->update(propagatedState, *innerhit);
00614 
00615   TSOS propagatedState_out = thePropagator->propagate(updatedState,outerhit->det()->surface()) ;
00616   if (!propagatedState_out.isValid())
00617     return false;
00618   TSOS updatedState_out = theUpdator->update(propagatedState_out, *outerhit);
00619   // debug prints
00620   LogDebug("") <<"[ElectronSeedGenerator::prepareElTrackSeed] final TSOS, position: "<<updatedState_out.globalPosition()<<" momentum: "<<updatedState_out.globalMomentum();
00621   LogDebug("") <<"[ElectronSeedGenerator::prepareElTrackSeed] final TSOS Pt: "<<updatedState_out.globalMomentum().perp();
00622   pts_ =  transformer_.persistentState(updatedState_out, outerhit->geographicalId().rawId());
00623 
00624   return true;
00625  }
CMSSW_4_4_3_patch1/src/RecoEgamma/EgammaElectronAlgos/src/ElectronSeedGenerator.cc
