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
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
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
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()))
00116 continue;
00117
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()))
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
00140
00141 shared++;
00142 }
00143 }
00144 }
00145
00146
00147 return shared;
00148
00149 }
00150
00151 int sharedDets(const GsfTrackRef& gsfTrackRef1, const
00152 GsfTrackRef& gsfTrackRef2 ) {
00153
00154
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()))
00164 continue;
00165
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()))
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 }
00188 }
00189
00190
00191
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
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
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 }