CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_4_5_patch3/src/RecoEgamma/EgammaElectronAlgos/src/EgAmbiguityTools.cc

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 //#include "RecoEgamma/EgammaElectronAlgos/interface/ElectronHcalHelper.h"
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   // retreive first valid hit
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   //get the Hit Pattern for the gsfTracks
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()))  //count only valid Hits
00078       continue;
00079     //if (gsfHitCounter1>1) continue; // test only the first hit of the track 1
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())) //count only valid Hits!
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 //        if (comp.equals(&(**elHitsIt1),&(**gsfHitsIt2))) {
00102         //std::cout << "found shared hit " << gsfHit2 << std::endl;
00103         shared++;
00104       }
00105     }//gsfHits2 iterator
00106   }//gsfHits1 iterator
00107 
00108   //std::cout << "[sharedHits] number of shared hits " << shared << std::endl;
00109   return shared;
00110 
00111 }
00112 
00113 int sharedDets(const GsfTrackRef& gsfTrackRef1, const
00114  GsfTrackRef& gsfTrackRef2 ) {
00115 
00116   //get the Hit Pattern for the gsfTracks
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()))  //count only valid Hits
00126       continue;
00127     //if (gsfHitCounter1>1) continue; // to test only the first hit of the track 1
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())) //count only valid Hits!
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     }//gsfHits2 iterator
00150   }//gsfHits1 iterator
00151 
00152   //std::cout << "[sharedHits] number of shared dets " << shared << std::endl;
00153   //return shared/min(gsfTrackRef1->numberOfValidHits(),gsfTrackRef2->numberOfValidHits());
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       // here we have common Xtal id
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   //std::cout << "[sharedEnergy] shared energy /min(energy1,energy2) " << fractionShared << std::endl;
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 }