Reference Manual

Go to the documentation of this file.
00001 #include "RecoEgamma/EgammaElectronAlgos/interface/EgAmbiguityTools.h"
00002 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
00003 
00004 #include "RecoEgamma/EgammaElectronAlgos/interface/GsfElectronAlgo.h"
00005 #include "RecoEgamma/EgammaElectronAlgos/interface/EgAmbiguityTools.h"
00006 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronClassification.h"
00007 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronMomentumCorrector.h"
00008 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronEnergyCorrector.h"
00009 //#include "RecoEgamma/EgammaElectronAlgos/interface/ElectronHcalHelper.h"
00010 
00011 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
00012 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
00013 #include "DataFormats/EgammaReco/interface/ElectronSeed.h"
00014 #include "DataFormats/EgammaReco/interface/ElectronSeedFwd.h"
00015 #include "DataFormats/TrackCandidate/interface/TrackCandidate.h"
00016 #include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
00017 #include "DataFormats/Math/interface/LorentzVector.h"
00018 #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
00019 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
00020 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
00021 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00022 #include "DataFormats/EcalDetId/interface/EEDetId.h"
00023 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
00024 
00025 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00026 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00027 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
00028 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00029 #include "Geometry/Records/interface/CaloTopologyRecord.h"
00030 
00031 #include "TrackingTools/PatternTools/interface/Trajectory.h"
00032 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
00033 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00034 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00035 #include "TrackingTools/GsfTools/interface/MultiTrajectoryStateTransform.h"
00036 #include "TrackingTools/GsfTools/interface/MultiTrajectoryStateMode.h"
00037 
00038 #include "TrackingTools/TransientTrackingRecHit/interface/RecHitComparatorByPosition.h"
00039 
00040 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
00041 
00042 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
00043 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
00044 
00045 #include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
00046 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00047 
00048 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00049 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
00050 
00051 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00052 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00053 
00054 #include "CLHEP/Units/GlobalPhysicalConstants.h"
00055 #include <TMath.h>
00056 #include <Math/VectorUtil.h>
00057 #include <Math/Point3D.h>
00058 #include <sstream>
00059 #include <algorithm>
00060 
00061 using namespace edm ;
00062 using namespace std ;
00063 using namespace reco ;
00064 
00065 namespace EgAmbiguityTools
00066  {
00067 
00068 bool isBetter( const reco::GsfElectron * e1, const reco::GsfElectron * e2 )
00069  { return (std::abs(e1->eSuperClusterOverP()-1)<std::abs(e2->eSuperClusterOverP()-1)) ; }
00070 
00071 bool isInnerMost::operator()
00072  ( const reco::GsfElectron * e1, const reco::GsfElectron * e2 )
00073  {
00074   reco::HitPattern gsfHitPattern1 = e1->gsfTrack()->hitPattern() ;
00075   reco::HitPattern gsfHitPattern2 = e2->gsfTrack()->hitPattern() ;
00076 
00077   // retreive first valid hit
00078   int gsfHitCounter1 = 0 ;
00079   trackingRecHit_iterator elHitsIt1 ;
00080   for
00081    ( elHitsIt1 = e1->gsfTrack()->recHitsBegin() ;
00082      elHitsIt1 != e1->gsfTrack()->recHitsEnd() ;
00083      elHitsIt1++, gsfHitCounter1++ )
00084    { if (((**elHitsIt1).isValid())) break ; }
00085 
00086   int gsfHitCounter2 = 0 ;
00087   trackingRecHit_iterator elHitsIt2 ;
00088   for
00089    ( elHitsIt2 = e2->gsfTrack()->recHitsBegin() ;
00090      elHitsIt2 != e2->gsfTrack()->recHitsEnd() ;
00091      elHitsIt2++, gsfHitCounter2++ )
00092    { if (((**elHitsIt2).isValid())) break ; }
00093 
00094   uint32_t gsfHit1 = gsfHitPattern1.getHitPattern(gsfHitCounter1) ;
00095   uint32_t gsfHit2 = gsfHitPattern2.getHitPattern(gsfHitCounter2) ;
00096   if (gsfHitPattern1.getSubStructure(gsfHit1)!=gsfHitPattern2.getSubStructure(gsfHit2))
00097    { return (gsfHitPattern1.getSubStructure(gsfHit1)<gsfHitPattern2.getSubStructure(gsfHit2)) ; }
00098   else if (gsfHitPattern1.getLayer(gsfHit1)!=gsfHitPattern2.getLayer(gsfHit2))
00099    { return (gsfHitPattern1.getLayer(gsfHit1)<gsfHitPattern2.getLayer(gsfHit2)) ; }
00100   else
00101    { return isBetter(e1,e2) ; }
00102  }
00103 
00104 int sharedHits( const GsfTrackRef & gsfTrackRef1, const GsfTrackRef & gsfTrackRef2 )
00105  {
00106   //get the Hit Pattern for the gsfTracks
00107   const HitPattern & gsfHitPattern1 = gsfTrackRef1->hitPattern() ;
00108   const HitPattern & gsfHitPattern2 = gsfTrackRef2->hitPattern() ;
00109 
00110   unsigned int shared = 0;
00111 
00112   int gsfHitCounter1 = 0;
00113   for(trackingRecHit_iterator elHitsIt1 = gsfTrackRef1->recHitsBegin();
00114       elHitsIt1 != gsfTrackRef1->recHitsEnd(); elHitsIt1++, gsfHitCounter1++) {
00115     if(!((**elHitsIt1).isValid()))  //count only valid Hits
00116       continue;
00117     //if (gsfHitCounter1>1) continue; // test only the first hit of the track 1
00118     uint32_t gsfHit = gsfHitPattern1.getHitPattern(gsfHitCounter1);
00119     if(!(gsfHitPattern1.pixelHitFilter(gsfHit)
00120          || gsfHitPattern1.stripTIBHitFilter(gsfHit)
00121          || gsfHitPattern1.stripTOBHitFilter(gsfHit)
00122          || gsfHitPattern1.stripTECHitFilter(gsfHit)
00123          || gsfHitPattern1.stripTIDHitFilter(gsfHit) )
00124     ) continue;
00125     int gsfHitsCounter2 = 0;
00126     for(trackingRecHit_iterator gsfHitsIt2 = gsfTrackRef2->recHitsBegin();
00127         gsfHitsIt2 != gsfTrackRef2->recHitsEnd(); gsfHitsIt2++, gsfHitsCounter2++) {
00128       if(!((**gsfHitsIt2).isValid())) //count only valid Hits!
00129         continue;
00130 
00131       uint32_t gsfHit2 = gsfHitPattern2.getHitPattern(gsfHitsCounter2);
00132       if(!(gsfHitPattern2.pixelHitFilter(gsfHit2)
00133          || gsfHitPattern2.stripTIBHitFilter(gsfHit2)
00134          || gsfHitPattern2.stripTOBHitFilter(gsfHit2)
00135          || gsfHitPattern2.stripTECHitFilter(gsfHit2)
00136          || gsfHitPattern2.stripTIDHitFilter(gsfHit2) )
00137       ) continue;
00138       if( (**elHitsIt1).sharesInput(&(**gsfHitsIt2), TrackingRecHit::some) ) {
00139 //        if (comp.equals(&(**elHitsIt1),&(**gsfHitsIt2))) {
00140         //std::cout << "found shared hit " << gsfHit2 << std::endl;
00141         shared++;
00142       }
00143     }//gsfHits2 iterator
00144   }//gsfHits1 iterator
00145 
00146   //std::cout << "[sharedHits] number of shared hits " << shared << std::endl;
00147   return shared;
00148 
00149 }
00150 
00151 int sharedDets(const GsfTrackRef& gsfTrackRef1, const
00152  GsfTrackRef& gsfTrackRef2 ) {
00153 
00154   //get the Hit Pattern for the gsfTracks
00155   const HitPattern& gsfHitPattern1 = gsfTrackRef1->hitPattern();
00156   const HitPattern& gsfHitPattern2 = gsfTrackRef2->hitPattern();
00157 
00158   unsigned int shared = 0;
00159 
00160   int gsfHitCounter1 = 0;
00161   for(trackingRecHit_iterator elHitsIt1 = gsfTrackRef1->recHitsBegin();
00162       elHitsIt1 != gsfTrackRef1->recHitsEnd(); elHitsIt1++, gsfHitCounter1++) {
00163     if(!((**elHitsIt1).isValid()))  //count only valid Hits
00164       continue;
00165     //if (gsfHitCounter1>1) continue; // to test only the first hit of the track 1
00166     uint32_t gsfHit = gsfHitPattern1.getHitPattern(gsfHitCounter1);
00167     if(!(gsfHitPattern1.pixelHitFilter(gsfHit)
00168          || gsfHitPattern1.stripTIBHitFilter(gsfHit)
00169          || gsfHitPattern1.stripTOBHitFilter(gsfHit)
00170          || gsfHitPattern1.stripTECHitFilter(gsfHit)
00171          || gsfHitPattern1.stripTIDHitFilter(gsfHit) )
00172     ) continue;
00173     int gsfHitsCounter2 = 0;
00174     for(trackingRecHit_iterator gsfHitsIt2 = gsfTrackRef2->recHitsBegin();
00175         gsfHitsIt2 != gsfTrackRef2->recHitsEnd(); gsfHitsIt2++, gsfHitsCounter2++) {
00176       if(!((**gsfHitsIt2).isValid())) //count only valid Hits!
00177         continue;
00178 
00179       uint32_t gsfHit2 = gsfHitPattern2.getHitPattern(gsfHitsCounter2);
00180       if(!(gsfHitPattern2.pixelHitFilter(gsfHit2)
00181            || gsfHitPattern2.stripTIBHitFilter(gsfHit2)
00182            || gsfHitPattern1.stripTOBHitFilter(gsfHit2)
00183            || gsfHitPattern2.stripTECHitFilter(gsfHit2)
00184            || gsfHitPattern2.stripTIDHitFilter(gsfHit2) )
00185       ) continue;
00186       if ((**elHitsIt1).geographicalId() == (**gsfHitsIt2).geographicalId()) shared++;
00187     }//gsfHits2 iterator
00188   }//gsfHits1 iterator
00189 
00190   //std::cout << "[sharedHits] number of shared dets " << shared << std::endl;
00191   //return shared/min(gsfTrackRef1->numberOfValidHits(),gsfTrackRef2->numberOfValidHits());
00192   return shared;
00193 
00194 }
00195 
00196 float sharedEnergy(const CaloCluster *clu1, const CaloCluster *clu2,
00197        edm::Handle<EcalRecHitCollection> & reducedEBRecHits,
00198        edm::Handle<EcalRecHitCollection> & reducedEERecHits ) {
00199 
00200   double fractionShared = 0;
00201 
00202   std::vector< std::pair<DetId, float> > v_id1 = clu1->hitsAndFractions();
00203   std::vector< std::pair<DetId, float> > v_id2 = clu2->hitsAndFractions();
00204   std::vector< std::pair<DetId, float> >::iterator ih1;
00205   std::vector< std::pair<DetId, float> >::iterator ih2;
00206 
00207   for(ih1 = v_id1.begin();ih1 != v_id1.end(); ih1++) {
00208 
00209     for(ih2 = v_id2.begin();ih2 != v_id2.end(); ih2++) {
00210 
00211       if ( (*ih1).first != (*ih2).first ) continue;
00212 
00213       // here we have common Xtal id
00214       EcalRecHitCollection::const_iterator itt;
00215       if ((*ih1).first.subdetId() == EcalBarrel) {
00216         if ((itt=reducedEBRecHits->find((*ih1).first))!=reducedEBRecHits->end())
00217         fractionShared += itt->energy();
00218       } else if ((*ih1).first.subdetId() == EcalEndcap) {
00219         if ((itt=reducedEERecHits->find((*ih1).first))!=reducedEERecHits->end())
00220         fractionShared += itt->energy();
00221       }
00222 
00223     }
00224   }
00225 
00226   //std::cout << "[sharedEnergy] shared energy /min(energy1,energy2) " << fractionShared << std::endl;
00227   return fractionShared;
00228 
00229 }
00230 
00231 float sharedEnergy(const SuperClusterRef& sc1, const SuperClusterRef& sc2,
00232        edm::Handle<EcalRecHitCollection> & reducedEBRecHits,
00233        edm::Handle<EcalRecHitCollection> & reducedEERecHits ) {
00234 
00235   double energyShared = 0;
00236   for(CaloCluster_iterator icl1=sc1->clustersBegin();icl1!=sc1->clustersEnd(); icl1++) {
00237     for(CaloCluster_iterator icl2=sc2->clustersBegin();icl2!=sc2->clustersEnd(); icl2++) {
00238       energyShared += sharedEnergy(&(**icl1),&(**icl2),reducedEBRecHits,reducedEERecHits );
00239     }
00240   }
00241   return energyShared;
00242 
00243 }
00244 
00245 }
/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/RecoEgamma/EgammaElectronAlgos/src/EgAmbiguityTools.cc
