EgammaElectronAlgos/src/EgAmbiguityTools.cc

0001 #include "RecoEgamma/EgammaElectronAlgos/interface/EgAmbiguityTools.h"
0002 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0003 #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
0004 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0005 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0006 #include "TrackingTools/GsfTools/interface/MultiTrajectoryStateTransform.h"
0007
0008 #include <algorithm>
0009
0010 using namespace edm;
0011 using namespace std;
0012 using namespace reco;
0013
0014 namespace egamma {
0015
0016   bool isBetterElectron(reco::GsfElectron const& e1, reco::GsfElectron const& e2) {
0017     return (std::abs(e1.eSuperClusterOverP() - 1) < std::abs(e2.eSuperClusterOverP() - 1));
0018   }
0019
0020   bool isInnermostElectron(reco::GsfElectron const& e1, reco::GsfElectron const& e2) {
0021     // retreive first valid hit

0022     int gsfHitCounter1 = 0;
0023     for (auto const& elHit : e1.gsfTrack()->recHits()) {
0024       if (elHit->isValid())
0025         break;
0026       gsfHitCounter1++;
0027     }
0028
0029     int gsfHitCounter2 = 0;
0030     for (auto const& elHit : e2.gsfTrack()->recHits()) {
0031       if (elHit->isValid())
0032         break;
0033       gsfHitCounter2++;
0034     }
0035
0036     uint32_t gsfHit1 = e1.gsfTrack()->hitPattern().getHitPattern(HitPattern::TRACK_HITS, gsfHitCounter1);
0037     uint32_t gsfHit2 = e2.gsfTrack()->hitPattern().getHitPattern(HitPattern::TRACK_HITS, gsfHitCounter2);
0038
0039     if (HitPattern::getSubStructure(gsfHit1) != HitPattern::getSubStructure(gsfHit2)) {
0040       return (HitPattern::getSubStructure(gsfHit1) < HitPattern::getSubStructure(gsfHit2));
0041     } else if (HitPattern::getLayer(gsfHit1) != HitPattern::getLayer(gsfHit2)) {
0042       return (HitPattern::getLayer(gsfHit1) < HitPattern::getLayer(gsfHit2));
0043     } else {
0044       return isBetterElectron(e1, e2);
0045     }
0046   }
0047
0048   int sharedHits(const GsfTrackRef& gsfTrackRef1, const GsfTrackRef& gsfTrackRef2) {
0049     //get the Hit Pattern for the gsfTracks

0050     HitPattern const& gsfHitPattern1 = gsfTrackRef1->hitPattern();
0051     HitPattern const& gsfHitPattern2 = gsfTrackRef2->hitPattern();
0052
0053     unsigned int shared = 0;
0054
0055     int gsfHitCounter1 = 0;
0056     for (trackingRecHit_iterator elHitsIt1 = gsfTrackRef1->recHitsBegin(); elHitsIt1 != gsfTrackRef1->recHitsEnd();
0057          elHitsIt1++, gsfHitCounter1++) {
0058       if (!(*elHitsIt1)->isValid()) {
0059         //count only valid Hits

0060         continue;
0061       }
0062       //if (gsfHitCounter1>1) continue; // test only the first hit of the track 1

0063       uint32_t gsfHit = gsfHitPattern1.getHitPattern(HitPattern::TRACK_HITS, gsfHitCounter1);
0064       if (!(HitPattern::pixelHitFilter(gsfHit) || HitPattern::stripTIBHitFilter(gsfHit) ||
0065             HitPattern::stripTOBHitFilter(gsfHit) || HitPattern::stripTECHitFilter(gsfHit) ||
0066             HitPattern::stripTIDHitFilter(gsfHit))) {
0067         continue;
0068       }
0069
0070       int gsfHitsCounter2 = 0;
0071       for (trackingRecHit_iterator gsfHitsIt2 = gsfTrackRef2->recHitsBegin(); gsfHitsIt2 != gsfTrackRef2->recHitsEnd();
0072            gsfHitsIt2++, gsfHitsCounter2++) {
0073         if (!(**gsfHitsIt2).isValid()) {
0074           //count only valid Hits

0075           continue;
0076         }
0077         uint32_t gsfHit2 = gsfHitPattern2.getHitPattern(HitPattern::TRACK_HITS, gsfHitsCounter2);
0078         if (!(HitPattern::pixelHitFilter(gsfHit2) || HitPattern::stripTIBHitFilter(gsfHit2) ||
0079               HitPattern::stripTOBHitFilter(gsfHit2) || HitPattern::stripTECHitFilter(gsfHit2) ||
0080               HitPattern::stripTIDHitFilter(gsfHit2))) {
0081           continue;
0082         }
0083         if ((*elHitsIt1)->sharesInput(&(**gsfHitsIt2), TrackingRecHit::some)) {
0084           //if (comp.equals(&(**elHitsIt1),&(**gsfHitsIt2))) {

0085           ////std::cout << "found shared hit " << gsfHit2 << std::endl;

0086           shared++;
0087         }
0088       }  //gsfHits2 iterator

0089     }  //gsfHits1 iterator

0090
0091     //std::cout << "[sharedHits] number of shared hits " << shared << std::endl;

0092     return shared;
0093   }
0094
0095   int sharedDets(const GsfTrackRef& gsfTrackRef1, const GsfTrackRef& gsfTrackRef2) {
0096     //get the Hit Pattern for the gsfTracks

0097     const HitPattern& gsfHitPattern1 = gsfTrackRef1->hitPattern();
0098     const HitPattern& gsfHitPattern2 = gsfTrackRef2->hitPattern();
0099
0100     unsigned int shared = 0;
0101
0102     int gsfHitCounter1 = 0;
0103     for (trackingRecHit_iterator elHitsIt1 = gsfTrackRef1->recHitsBegin(); elHitsIt1 != gsfTrackRef1->recHitsEnd();
0104          elHitsIt1++, gsfHitCounter1++) {
0105       if (!((**elHitsIt1).isValid())) {
0106         //count only valid Hits

0107         continue;
0108       }
0109       //if (gsfHitCounter1>1) continue; // to test only the first hit of the track 1

0110       uint32_t gsfHit = gsfHitPattern1.getHitPattern(HitPattern::TRACK_HITS, gsfHitCounter1);
0111       if (!(HitPattern::pixelHitFilter(gsfHit) || HitPattern::stripTIBHitFilter(gsfHit) ||
0112             HitPattern::stripTOBHitFilter(gsfHit) || HitPattern::stripTECHitFilter(gsfHit) ||
0113             HitPattern::stripTIDHitFilter(gsfHit))) {
0114         continue;
0115       }
0116
0117       int gsfHitsCounter2 = 0;
0118       for (trackingRecHit_iterator gsfHitsIt2 = gsfTrackRef2->recHitsBegin(); gsfHitsIt2 != gsfTrackRef2->recHitsEnd();
0119            gsfHitsIt2++, gsfHitsCounter2++) {
0120         if (!((**gsfHitsIt2).isValid())) {
0121           //count only valid Hits!

0122           continue;
0123         }
0124
0125         uint32_t gsfHit2 = gsfHitPattern2.getHitPattern(HitPattern::TRACK_HITS, gsfHitsCounter2);
0126         if (!(HitPattern::pixelHitFilter(gsfHit2) || HitPattern::stripTIBHitFilter(gsfHit2) ||
0127               HitPattern::stripTOBHitFilter(gsfHit2) || HitPattern::stripTECHitFilter(gsfHit2) ||
0128               HitPattern::stripTIDHitFilter(gsfHit2)))
0129           continue;
0130         if ((**elHitsIt1).geographicalId() == (**gsfHitsIt2).geographicalId())
0131           shared++;
0132       }  //gsfHits2 iterator

0133     }  //gsfHits1 iterator

0134
0135     //std::cout << "[sharedHits] number of shared dets " << shared << std::endl;

0136     //return shared/min(gsfTrackRef1->numberOfValidHits(),gsfTrackRef2->numberOfValidHits());

0137     return shared;
0138   }
0139
0140   float sharedEnergy(CaloCluster const& clu1,
0141                      CaloCluster const& clu2,
0142                      EcalRecHitCollection const& barrelRecHits,
0143                      EcalRecHitCollection const& endcapRecHits) {
0144     double fractionShared = 0;
0145
0146     for (auto const& h1 : clu1.hitsAndFractions()) {
0147       for (auto const& h2 : clu2.hitsAndFractions()) {
0148         if (h1.first != h2.first)
0149           continue;
0150
0151         // here we have common Xtal id

0152         EcalRecHitCollection::const_iterator itt;
0153         if (h1.first.subdetId() == EcalBarrel) {
0154           if ((itt = barrelRecHits.find(h1.first)) != barrelRecHits.end())
0155             fractionShared += itt->energy();
0156         } else if (h1.first.subdetId() == EcalEndcap) {
0157           if ((itt = endcapRecHits.find(h1.first)) != endcapRecHits.end())
0158             fractionShared += itt->energy();
0159         }
0160       }
0161     }
0162
0163     //std::cout << "[sharedEnergy] shared energy /min(energy1,energy2) " << fractionShared << std::endl;

0164     return fractionShared;
0165   }
0166
0167   float sharedEnergy(SuperClusterRef const& sc1,
0168                      SuperClusterRef const& sc2,
0169                      EcalRecHitCollection const& barrelRecHits,
0170                      EcalRecHitCollection const& endcapRecHits) {
0171     double energyShared = 0;
0172     for (CaloCluster_iterator icl1 = sc1->clustersBegin(); icl1 != sc1->clustersEnd(); icl1++) {
0173       for (CaloCluster_iterator icl2 = sc2->clustersBegin(); icl2 != sc2->clustersEnd(); icl2++) {
0174         energyShared += sharedEnergy(**icl1, **icl2, barrelRecHits, endcapRecHits);
0175       }
0176     }
0177     return energyShared;
0178   }
0179
0180 }  // namespace egamma