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/EgAmbiguityTools.h"
00005
00006
00007 #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
00008 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
00009 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
00010
00011
00012 #include "TrackingTools/GsfTools/interface/MultiTrajectoryStateTransform.h"
00013
00014
00015
00016
00017
00018
00019
00020 #include <sstream>
00021 #include <algorithm>
00022
00023 using namespace edm ;
00024 using namespace std ;
00025 using namespace reco ;
00026
00027 namespace EgAmbiguityTools
00028 {
00029
00030 bool isBetter( const reco::GsfElectron * e1, const reco::GsfElectron * e2 )
00031 { return (std::abs(e1->eSuperClusterOverP()-1)<std::abs(e2->eSuperClusterOverP()-1)) ; }
00032
00033 bool isInnerMost::operator()
00034 ( const reco::GsfElectron * e1, const reco::GsfElectron * e2 )
00035 {
00036 reco::HitPattern gsfHitPattern1 = e1->gsfTrack()->hitPattern() ;
00037 reco::HitPattern gsfHitPattern2 = e2->gsfTrack()->hitPattern() ;
00038
00039
00040 int gsfHitCounter1 = 0 ;
00041 trackingRecHit_iterator elHitsIt1 ;
00042 for
00043 ( elHitsIt1 = e1->gsfTrack()->recHitsBegin() ;
00044 elHitsIt1 != e1->gsfTrack()->recHitsEnd() ;
00045 elHitsIt1++, gsfHitCounter1++ )
00046 { if (((**elHitsIt1).isValid())) break ; }
00047
00048 int gsfHitCounter2 = 0 ;
00049 trackingRecHit_iterator elHitsIt2 ;
00050 for
00051 ( elHitsIt2 = e2->gsfTrack()->recHitsBegin() ;
00052 elHitsIt2 != e2->gsfTrack()->recHitsEnd() ;
00053 elHitsIt2++, gsfHitCounter2++ )
00054 { if (((**elHitsIt2).isValid())) break ; }
00055
00056 uint32_t gsfHit1 = gsfHitPattern1.getHitPattern(gsfHitCounter1) ;
00057 uint32_t gsfHit2 = gsfHitPattern2.getHitPattern(gsfHitCounter2) ;
00058 if (gsfHitPattern1.getSubStructure(gsfHit1)!=gsfHitPattern2.getSubStructure(gsfHit2))
00059 { return (gsfHitPattern1.getSubStructure(gsfHit1)<gsfHitPattern2.getSubStructure(gsfHit2)) ; }
00060 else if (gsfHitPattern1.getLayer(gsfHit1)!=gsfHitPattern2.getLayer(gsfHit2))
00061 { return (gsfHitPattern1.getLayer(gsfHit1)<gsfHitPattern2.getLayer(gsfHit2)) ; }
00062 else
00063 { return isBetter(e1,e2) ; }
00064 }
00065
00066 int sharedHits( const GsfTrackRef & gsfTrackRef1, const GsfTrackRef & gsfTrackRef2 )
00067 {
00068
00069 const HitPattern & gsfHitPattern1 = gsfTrackRef1->hitPattern() ;
00070 const HitPattern & gsfHitPattern2 = gsfTrackRef2->hitPattern() ;
00071
00072 unsigned int shared = 0;
00073
00074 int gsfHitCounter1 = 0;
00075 for(trackingRecHit_iterator elHitsIt1 = gsfTrackRef1->recHitsBegin();
00076 elHitsIt1 != gsfTrackRef1->recHitsEnd(); elHitsIt1++, gsfHitCounter1++) {
00077 if(!((**elHitsIt1).isValid()))
00078 continue;
00079
00080 uint32_t gsfHit = gsfHitPattern1.getHitPattern(gsfHitCounter1);
00081 if(!(gsfHitPattern1.pixelHitFilter(gsfHit)
00082 || gsfHitPattern1.stripTIBHitFilter(gsfHit)
00083 || gsfHitPattern1.stripTOBHitFilter(gsfHit)
00084 || gsfHitPattern1.stripTECHitFilter(gsfHit)
00085 || gsfHitPattern1.stripTIDHitFilter(gsfHit) )
00086 ) continue;
00087 int gsfHitsCounter2 = 0;
00088 for(trackingRecHit_iterator gsfHitsIt2 = gsfTrackRef2->recHitsBegin();
00089 gsfHitsIt2 != gsfTrackRef2->recHitsEnd(); gsfHitsIt2++, gsfHitsCounter2++) {
00090 if(!((**gsfHitsIt2).isValid()))
00091 continue;
00092
00093 uint32_t gsfHit2 = gsfHitPattern2.getHitPattern(gsfHitsCounter2);
00094 if(!(gsfHitPattern2.pixelHitFilter(gsfHit2)
00095 || gsfHitPattern2.stripTIBHitFilter(gsfHit2)
00096 || gsfHitPattern2.stripTOBHitFilter(gsfHit2)
00097 || gsfHitPattern2.stripTECHitFilter(gsfHit2)
00098 || gsfHitPattern2.stripTIDHitFilter(gsfHit2) )
00099 ) continue;
00100 if( (**elHitsIt1).sharesInput(&(**gsfHitsIt2), TrackingRecHit::some) ) {
00101
00102
00103 shared++;
00104 }
00105 }
00106 }
00107
00108
00109 return shared;
00110
00111 }
00112
00113 int sharedDets(const GsfTrackRef& gsfTrackRef1, const
00114 GsfTrackRef& gsfTrackRef2 ) {
00115
00116
00117 const HitPattern& gsfHitPattern1 = gsfTrackRef1->hitPattern();
00118 const HitPattern& gsfHitPattern2 = gsfTrackRef2->hitPattern();
00119
00120 unsigned int shared = 0;
00121
00122 int gsfHitCounter1 = 0;
00123 for(trackingRecHit_iterator elHitsIt1 = gsfTrackRef1->recHitsBegin();
00124 elHitsIt1 != gsfTrackRef1->recHitsEnd(); elHitsIt1++, gsfHitCounter1++) {
00125 if(!((**elHitsIt1).isValid()))
00126 continue;
00127
00128 uint32_t gsfHit = gsfHitPattern1.getHitPattern(gsfHitCounter1);
00129 if(!(gsfHitPattern1.pixelHitFilter(gsfHit)
00130 || gsfHitPattern1.stripTIBHitFilter(gsfHit)
00131 || gsfHitPattern1.stripTOBHitFilter(gsfHit)
00132 || gsfHitPattern1.stripTECHitFilter(gsfHit)
00133 || gsfHitPattern1.stripTIDHitFilter(gsfHit) )
00134 ) continue;
00135 int gsfHitsCounter2 = 0;
00136 for(trackingRecHit_iterator gsfHitsIt2 = gsfTrackRef2->recHitsBegin();
00137 gsfHitsIt2 != gsfTrackRef2->recHitsEnd(); gsfHitsIt2++, gsfHitsCounter2++) {
00138 if(!((**gsfHitsIt2).isValid()))
00139 continue;
00140
00141 uint32_t gsfHit2 = gsfHitPattern2.getHitPattern(gsfHitsCounter2);
00142 if(!(gsfHitPattern2.pixelHitFilter(gsfHit2)
00143 || gsfHitPattern2.stripTIBHitFilter(gsfHit2)
00144 || gsfHitPattern1.stripTOBHitFilter(gsfHit2)
00145 || gsfHitPattern2.stripTECHitFilter(gsfHit2)
00146 || gsfHitPattern2.stripTIDHitFilter(gsfHit2) )
00147 ) continue;
00148 if ((**elHitsIt1).geographicalId() == (**gsfHitsIt2).geographicalId()) shared++;
00149 }
00150 }
00151
00152
00153
00154 return shared;
00155
00156 }
00157
00158 float sharedEnergy(const CaloCluster *clu1, const CaloCluster *clu2,
00159 edm::Handle<EcalRecHitCollection> & reducedEBRecHits,
00160 edm::Handle<EcalRecHitCollection> & reducedEERecHits ) {
00161
00162 double fractionShared = 0;
00163
00164 std::vector< std::pair<DetId, float> > v_id1 = clu1->hitsAndFractions();
00165 std::vector< std::pair<DetId, float> > v_id2 = clu2->hitsAndFractions();
00166 std::vector< std::pair<DetId, float> >::iterator ih1;
00167 std::vector< std::pair<DetId, float> >::iterator ih2;
00168
00169 for(ih1 = v_id1.begin();ih1 != v_id1.end(); ih1++) {
00170
00171 for(ih2 = v_id2.begin();ih2 != v_id2.end(); ih2++) {
00172
00173 if ( (*ih1).first != (*ih2).first ) continue;
00174
00175
00176 EcalRecHitCollection::const_iterator itt;
00177 if ((*ih1).first.subdetId() == EcalBarrel) {
00178 if ((itt=reducedEBRecHits->find((*ih1).first))!=reducedEBRecHits->end())
00179 fractionShared += itt->energy();
00180 } else if ((*ih1).first.subdetId() == EcalEndcap) {
00181 if ((itt=reducedEERecHits->find((*ih1).first))!=reducedEERecHits->end())
00182 fractionShared += itt->energy();
00183 }
00184
00185 }
00186 }
00187
00188
00189 return fractionShared;
00190
00191 }
00192
00193 float sharedEnergy(const SuperClusterRef& sc1, const SuperClusterRef& sc2,
00194 edm::Handle<EcalRecHitCollection> & reducedEBRecHits,
00195 edm::Handle<EcalRecHitCollection> & reducedEERecHits ) {
00196
00197 double energyShared = 0;
00198 for(CaloCluster_iterator icl1=sc1->clustersBegin();icl1!=sc1->clustersEnd(); icl1++) {
00199 for(CaloCluster_iterator icl2=sc2->clustersBegin();icl2!=sc2->clustersEnd(); icl2++) {
00200 energyShared += sharedEnergy(&(**icl1),&(**icl2),reducedEBRecHits,reducedEERecHits );
00201 }
00202 }
00203 return energyShared;
00204
00205 }
00206
00207 }