diff --git a/DataFormats/EgammaCandidates/interface/Photon.h b/DataFormats/EgammaCandidates/interface/Photon.h index f4b3bc4e2bc1c..ddae30524c9b7 100644 --- a/DataFormats/EgammaCandidates/interface/Photon.h +++ b/DataFormats/EgammaCandidates/interface/Photon.h @@ -581,11 +581,9 @@ namespace reco { ///MVA based beam halo tagger - trained for EE and for pT > 200 GeV float haloTaggerMVAVal() const { return haloTaggerMVAVal_; } - + ///set the haloTaggerMVAVal here - void setHaloTaggerMVAVal(float x) { - haloTaggerMVAVal_ = x; - } + void setHaloTaggerMVAVal(float x) { haloTaggerMVAVal_ = x; } private: /// check overlap with another candidate diff --git a/DataFormats/EgammaCandidates/src/Photon.cc b/DataFormats/EgammaCandidates/src/Photon.cc index 27892904b9144..7633e35a966f5 100644 --- a/DataFormats/EgammaCandidates/src/Photon.cc +++ b/DataFormats/EgammaCandidates/src/Photon.cc @@ -4,7 +4,11 @@ using namespace reco; Photon::Photon(const LorentzVector& p4, const Point& caloPos, const PhotonCoreRef& core, const Point& vtx) - : RecoCandidate(0, p4, vtx, 22), caloPosition_(caloPos), photonCore_(core), pixelSeed_(false),haloTaggerMVAVal_(-999) {} + : RecoCandidate(0, p4, vtx, 22), + caloPosition_(caloPos), + photonCore_(core), + pixelSeed_(false), + haloTaggerMVAVal_(-999) {} Photon::Photon(const Photon& rhs) : RecoCandidate(rhs), @@ -19,8 +23,8 @@ Photon::Photon(const Photon& rhs) saturationInfo_(rhs.saturationInfo_), eCorrections_(rhs.eCorrections_), mipVariableBlock_(rhs.mipVariableBlock_), - pfIsolation_(rhs.pfIsolation_), - haloTaggerMVAVal_(rhs.haloTaggerMVAVal_) {} + pfIsolation_(rhs.pfIsolation_), + haloTaggerMVAVal_(rhs.haloTaggerMVAVal_) {} Photon::~Photon() {} diff --git a/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc b/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc index 30713677ace12..b85643ac35874 100644 --- a/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc +++ b/RecoEgamma/EgammaPhotonProducers/src/GEDPhotonProducer.cc @@ -418,8 +418,8 @@ GEDPhotonProducer::GEDPhotonProducer(const edm::ParameterSet& config, const Cach if (recoStep_.isFinal() && runMVABasedHaloTagger_) { edm::ParameterSet mvaBasedHaloVariableSet = config.getParameter("mvaBasedHaloVariableSet"); - photonMVABasedHaloTagger_ = std::make_unique(mvaBasedHaloVariableSet, consumesCollector()); - + photonMVABasedHaloTagger_ = + std::make_unique(mvaBasedHaloVariableSet, consumesCollector()); } ///Get the set for PF cluster isolation calculator @@ -1068,7 +1068,7 @@ void GEDPhotonProducer::fillPhotonCollection(edm::Event& evt, reco::Photon newCandidate(*phoRef); iSC++; - if (runMVABasedHaloTagger_) { ///sets values only for EE, for EB it always returns 1 + if (runMVABasedHaloTagger_) { ///sets values only for EE, for EB it always returns 1 float BHmva = photonMVABasedHaloTagger_->calculateMVA(&newCandidate, evt, es); newCandidate.setHaloTaggerMVAVal(BHmva); } diff --git a/RecoEgamma/EgammaPhotonProducers/src/PhotonProducer.cc b/RecoEgamma/EgammaPhotonProducers/src/PhotonProducer.cc index 4db6d3e9b0074..e733bc925bad5 100644 --- a/RecoEgamma/EgammaPhotonProducers/src/PhotonProducer.cc +++ b/RecoEgamma/EgammaPhotonProducers/src/PhotonProducer.cc @@ -238,9 +238,9 @@ PhotonProducer::PhotonProducer(const edm::ParameterSet& config) if (runMVABasedHaloTagger_) { edm::ParameterSet mvaBasedHaloVariableSet = config.getParameter("mvaBasedHaloVariableSet"); - photonMVABasedHaloTagger_ = std::make_unique(mvaBasedHaloVariableSet, consumesCollector()); + photonMVABasedHaloTagger_ = + std::make_unique(mvaBasedHaloVariableSet, consumesCollector()); } - // Register the product produces(PhotonCollection_); @@ -498,7 +498,7 @@ void PhotonProducer::fillPhotonCollection(edm::Event& evt, newCandidate.setMIPVariables(mipVar); } - if (runMVABasedHaloTagger_) { ///sets values only for EE, for EB it always returns 1 + if (runMVABasedHaloTagger_) { ///sets values only for EE, for EB it always returns 1 double BHmva = photonMVABasedHaloTagger_->calculateMVA(&newCandidate, evt, es); newCandidate.setHaloTaggerMVAVal(BHmva); } diff --git a/RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h b/RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h index f3498f526fc8b..b25d0e2f3e35f 100644 --- a/RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h +++ b/RecoEgamma/PhotonIdentification/interface/PhotonMVABasedHaloTagger.h @@ -20,7 +20,7 @@ #include "CondFormats/GBRForest/interface/GBRForest.h" #include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h" #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaHcalIsolation.h" - #include "CondFormats/EcalObjects/interface/EcalPFRecHitThresholds.h" +#include "CondFormats/EcalObjects/interface/EcalPFRecHitThresholds.h" #include "CondFormats/DataRecord/interface/EcalPFRecHitThresholdsRcd.h" #include @@ -29,64 +29,52 @@ class PhotonMVABasedHaloTagger { PhotonMVABasedHaloTagger(const edm::ParameterSet& conf, edm::ConsumesCollector&& iC); virtual ~PhotonMVABasedHaloTagger(){}; - double calculateMVA(const reco::Photon* pho, - const edm::Event&iEvent, - const edm::EventSetup& es - ); - - + double calculateMVA(const reco::Photon* pho, const edm::Event& iEvent, const edm::EventSetup& es); + void calphoClusCoordinECAL(const CaloGeometry* geo, - const reco::Photon*, - const EcalPFRecHitThresholds *thresholds, - edm::Handle ecalrechitCollHandle - ); - + const reco::Photon*, + const EcalPFRecHitThresholds* thresholds, + edm::Handle ecalrechitCollHandle); void calmatchedHBHECoordForBothHypothesis(const CaloGeometry* geo, - const reco::Photon*, - edm::Handle hbheHitHandle - ); - - + const reco::Photon*, + edm::Handle hbheHitHandle); + void calmatchedESCoordForBothHypothesis(const CaloGeometry* geo, - const reco::Photon*, - edm::Handle esrechitCollHandle - ); - - double calAngleBetweenEEAndSubDet(int nhits, - double subdetClusX, - double subdetClusY, - double subdetClusZ); - + const reco::Photon*, + edm::Handle esrechitCollHandle); + + double calAngleBetweenEEAndSubDet(int nhits, double subdetClusX, double subdetClusY, double subdetClusZ); private: int hcalClusNhits_samedPhi_, hcalClusNhits_samedR_; - int ecalClusNhits_, preshowerNhits_samedPhi_, preshowerNhits_samedR_; - double hcalClusX_samedPhi_, hcalClusY_samedPhi_, hcalClusZ_samedPhi_, hcalClusX_samedR_, hcalClusY_samedR_, hcalClusZ_samedR_; + int ecalClusNhits_, preshowerNhits_samedPhi_, preshowerNhits_samedR_; + double hcalClusX_samedPhi_, hcalClusY_samedPhi_, hcalClusZ_samedPhi_, hcalClusX_samedR_, hcalClusY_samedR_, + hcalClusZ_samedR_; double hcalClusE_samedPhi_, hcalClusE_samedR_; - + double ecalClusX_, ecalClusY_, ecalClusZ_; - double preshowerX_samedPhi_, preshowerY_samedPhi_, preshowerZ_samedPhi_, preshowerX_samedR_, preshowerY_samedR_, preshowerZ_samedR_; + double preshowerX_samedPhi_, preshowerY_samedPhi_, preshowerZ_samedPhi_, preshowerX_samedR_, preshowerY_samedR_, + preshowerZ_samedR_; double ecalClusE_, preshowerE_samedPhi_, preshowerE_samedR_; double noiseThrES_; EgammaHcalIsolation::arrayHB recHitEThresholdHB_; EgammaHcalIsolation::arrayHE recHitEThresholdHE_; - std::string trainingFileName_; + std::string trainingFileName_; const edm::ESGetToken geometryToken_; const edm::ESGetToken ecalPFRechitThresholdsToken_; const EcalClusterLazyTools::ESGetTokens ecalClusterToolsESGetTokens_; - - edm::ESHandle pG_; - edm::EDGetTokenT rhoLabel_; + + edm::ESHandle pG_; + edm::EDGetTokenT rhoLabel_; edm::EDGetTokenT EBecalCollection_; edm::EDGetTokenT EEecalCollection_; edm::EDGetTokenT ESCollection_; edm::EDGetTokenT HBHERecHitsCollection_; std::unique_ptr gbr_; - }; #endif // PhotonMVABasedHaloTagger_H diff --git a/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc b/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc index 9154c8a7efc9c..29fdc57e3c23f 100644 --- a/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc +++ b/RecoEgamma/PhotonIdentification/src/PhotonMVABasedHaloTagger.cc @@ -15,15 +15,13 @@ #include #include "CommonTools/MVAUtils/interface/GBRForestTools.h" - PhotonMVABasedHaloTagger::PhotonMVABasedHaloTagger(const edm::ParameterSet& conf, edm::ConsumesCollector&& iC) - : geometryToken_(iC.esConsumes()), - ecalPFRechitThresholdsToken_(iC.esConsumes()), - ecalClusterToolsESGetTokens_(std::move(iC)) -{ + : geometryToken_(iC.esConsumes()), + ecalPFRechitThresholdsToken_(iC.esConsumes()), + ecalClusterToolsESGetTokens_(std::move(iC)) { trainingFileName_ = conf.getParameter("trainingFileName"); - rhoLabel_ = iC.consumes(conf.getParameter("rhoLabel")); - + rhoLabel_ = iC.consumes(conf.getParameter("rhoLabel")); + EBecalCollection_ = iC.consumes(conf.getParameter("barrelEcalRecHitCollection")); EEecalCollection_ = iC.consumes(conf.getParameter("endcapEcalRecHitCollection")); ESCollection_ = iC.consumes(conf.getParameter("esRecHitCollection")); @@ -31,34 +29,34 @@ PhotonMVABasedHaloTagger::PhotonMVABasedHaloTagger(const edm::ParameterSet& conf recHitEThresholdHB_ = conf.getParameter("recHitEThresholdHB"); recHitEThresholdHE_ = conf.getParameter("recHitEThresholdHE"); - noiseThrES_ = conf.getParameter("noiseThrES"); + noiseThrES_ = conf.getParameter("noiseThrES"); gbr_ = createGBRForest(trainingFileName_); - } -double PhotonMVABasedHaloTagger::calculateMVA(const reco::Photon* pho, const edm::Event& iEvent, const edm::EventSetup& es){ - - +double PhotonMVABasedHaloTagger::calculateMVA(const reco::Photon* pho, + const edm::Event& iEvent, + const edm::EventSetup& es) { bool isEB = pho->isEB(); - if(isEB) return 1.0; /// this MVA is useful and trained only for the EE photons. For EB, there are a lot of other useful handles which can reject beam halo efficiently - + if (isEB) + return 1.0; /// this MVA is useful and trained only for the EE photons. For EB, there are a lot of other useful handles which can reject beam halo efficiently + //rho handle edm::Handle rhoHandle; iEvent.getByToken(rhoLabel_, rhoHandle); - double rho_ = *(rhoHandle.product()); + double rho_ = *(rhoHandle.product()); // Get all the RecHits edm::Handle barrelHitHandle; iEvent.getByToken(EBecalCollection_, barrelHitHandle); - + edm::Handle endcapHitHandle; iEvent.getByToken(EEecalCollection_, endcapHitHandle); edm::Handle esHitHandle; iEvent.getByToken(ESCollection_, esHitHandle); - + edm::Handle hbheHitHandle; iEvent.getByToken(HBHERecHitsCollection_, hbheHitHandle); @@ -68,33 +66,42 @@ double PhotonMVABasedHaloTagger::calculateMVA(const reco::Photon* pho, const edm ///ECAL PF rechit thresholds auto const& thresholds = es.getData(ecalPFRechitThresholdsToken_); - - noZS::EcalClusterLazyTools lazyToolnoZS(iEvent, ecalClusterToolsESGetTokens_.get(es), EBecalCollection_, EEecalCollection_); + noZS::EcalClusterLazyTools lazyToolnoZS( + iEvent, ecalClusterToolsESGetTokens_.get(es), EBecalCollection_, EEecalCollection_); ///calculate the energy weighted X, Y and Z position of the photon cluster - if(isEB) calphoClusCoordinECAL(geo, pho, &thresholds, barrelHitHandle); - else calphoClusCoordinECAL(geo, pho, &thresholds, endcapHitHandle); - + if (isEB) + calphoClusCoordinECAL(geo, pho, &thresholds, barrelHitHandle); + else + calphoClusCoordinECAL(geo, pho, &thresholds, endcapHitHandle); + ///calculate the HBHE cluster position hypothesis calmatchedHBHECoordForBothHypothesis(geo, pho, hbheHitHandle); calmatchedESCoordForBothHypothesis(geo, pho, esHitHandle); ///this function works for EE only. Above ones work for EB as well in case later one wants to put a similar function for EB without returning 1 - double angle_EE_HE_samedPhi = calAngleBetweenEEAndSubDet(hcalClusNhits_samedPhi_, hcalClusX_samedPhi_, hcalClusY_samedPhi_, hcalClusZ_samedPhi_); //essentially caculates the angle and energy variables in the two hypothesis between EE and HE + double angle_EE_HE_samedPhi = calAngleBetweenEEAndSubDet( + hcalClusNhits_samedPhi_, + hcalClusX_samedPhi_, + hcalClusY_samedPhi_, + hcalClusZ_samedPhi_); //essentially caculates the angle and energy variables in the two hypothesis between EE and HE - double angle_EE_HE_samedR = calAngleBetweenEEAndSubDet(hcalClusNhits_samedR_, hcalClusX_samedR_, hcalClusY_samedR_, hcalClusZ_samedR_); + double angle_EE_HE_samedR = + calAngleBetweenEEAndSubDet(hcalClusNhits_samedR_, hcalClusX_samedR_, hcalClusY_samedR_, hcalClusZ_samedR_); - double angle_EE_ES_samedPhi = calAngleBetweenEEAndSubDet(preshowerNhits_samedPhi_, preshowerX_samedPhi_, preshowerY_samedPhi_, preshowerZ_samedPhi_); + double angle_EE_ES_samedPhi = calAngleBetweenEEAndSubDet( + preshowerNhits_samedPhi_, preshowerX_samedPhi_, preshowerY_samedPhi_, preshowerZ_samedPhi_); + + double angle_EE_ES_samedR = + calAngleBetweenEEAndSubDet(preshowerNhits_samedR_, preshowerX_samedR_, preshowerY_samedR_, preshowerZ_samedR_); - double angle_EE_ES_samedR = calAngleBetweenEEAndSubDet(preshowerNhits_samedR_, preshowerX_samedR_, preshowerY_samedR_, preshowerZ_samedR_); - ////set all the above calculated variables as input to the MVA - const auto& vCov = lazyToolnoZS.localCovariances( *(pho->superCluster()->seed()) ); + const auto& vCov = lazyToolnoZS.localCovariances(*(pho->superCluster()->seed())); double spp = (isnan(vCov[2]) ? 0. : sqrt(vCov[2])); - + ///https://cmssdt.cern.ch/lxr/source/RecoEgamma/ElectronIdentification/src/ElectronMVAEstimator.cc float vars[15]; @@ -112,221 +119,208 @@ double PhotonMVABasedHaloTagger::calculateMVA(const reco::Photon* pho, const edm vars[10] = angle_EE_HE_samedR; vars[11] = angle_EE_ES_samedPhi; vars[12] = angle_EE_HE_samedPhi; - vars[13] = (pho->superCluster()->preshowerEnergyPlane1() + pho->superCluster()->preshowerEnergyPlane2())/pho->superCluster()->rawEnergy(); + vars[13] = (pho->superCluster()->preshowerEnergyPlane1() + pho->superCluster()->preshowerEnergyPlane2()) / + pho->superCluster()->rawEnergy(); vars[14] = rho_; double BHmva = gbr_->GetAdaBoostClassifier(vars); return BHmva; - } - void PhotonMVABasedHaloTagger::calphoClusCoordinECAL(const CaloGeometry* geo, - const reco::Photon* pho, - const EcalPFRecHitThresholds *thresholds, - edm::Handle ecalrechitCollHandle - ){ - - ecalClusX_ = 0; ecalClusY_ = 0; ecalClusZ_ = 0; + const reco::Photon* pho, + const EcalPFRecHitThresholds* thresholds, + edm::Handle ecalrechitCollHandle) { + ecalClusX_ = 0; + ecalClusY_ = 0; + ecalClusZ_ = 0; ecalClusNhits_ = 0; ecalClusE_ = 0; double phoSCEta = pho->superCluster()->eta(); double phoSCPhi = pho->superCluster()->phi(); - + const EcalRecHitCollection* ecalRecHits = ecalrechitCollHandle.product(); - - for(EcalRecHitCollection::const_iterator ecalrechit = ecalRecHits->begin(); ecalrechit != ecalRecHits->end(); ecalrechit++){ - + + for (EcalRecHitCollection::const_iterator ecalrechit = ecalRecHits->begin(); ecalrechit != ecalRecHits->end(); + ecalrechit++) { auto const det = ecalrechit->id(); double rhE = ecalrechit->energy(); - const GlobalPoint &rechitPoint = geo->getPosition(det); - + const GlobalPoint& rechitPoint = geo->getPosition(det); + double rhEta = rechitPoint.eta(); double rhPhi = rechitPoint.phi(); double rhX = rechitPoint.x(); double rhY = rechitPoint.y(); double rhZ = rechitPoint.z(); - + if ((thresholds == nullptr)) { throw cms::Exception("EmptyPFRechHitThresCollection") - << "In PhotonMVABasedHaloTagger::calphoClusCoordinECAL, EcalPFRecHitThresholds cannot be = nulptr"; + << "In PhotonMVABasedHaloTagger::calphoClusCoordinECAL, EcalPFRecHitThresholds cannot be = nulptr"; } - + float rhThres = 0.0; if (thresholds != nullptr) { - rhThres = (*thresholds)[det]; + rhThres = (*thresholds)[det]; } - + if (rhE <= rhThres) continue; - - if(phoSCEta * rhEta < 0) continue; + if (phoSCEta * rhEta < 0) + continue; //double rho = sqrt(pow(rhX,2) + pow(rhY,2)); double dPhi = deltaPhi(rhPhi, phoSCPhi); double dEta = fabs(rhEta - phoSCEta); - - double dR = sqrt( pow(dPhi,2) + pow(dEta,2) ); - - if(dR < 0.2){ - + double dR = sqrt(pow(dPhi, 2) + pow(dEta, 2)); + + if (dR < 0.2) { ecalClusX_ += rhX * rhE; ecalClusY_ += rhY * rhE; ecalClusZ_ += rhZ * rhE; ecalClusE_ += rhE; ecalClusNhits_++; - } - }//for(int ih=0; ih 0){ //should always be > 0 for an EM cluster - ecalClusX_ = ecalClusX_/ecalClusE_; - ecalClusY_ = ecalClusY_/ecalClusE_; - ecalClusZ_ = ecalClusZ_/ecalClusE_; - }//if(ecalClusNhits_>0) - + } //for(int ih=0; ih 0) { //should always be > 0 for an EM cluster + ecalClusX_ = ecalClusX_ / ecalClusE_; + ecalClusY_ = ecalClusY_ / ecalClusE_; + ecalClusZ_ = ecalClusZ_ / ecalClusE_; + } //if(ecalClusNhits_>0) } void PhotonMVABasedHaloTagger::calmatchedHBHECoordForBothHypothesis( - const CaloGeometry* geo, - const reco::Photon* pho, - edm::Handle hcalhitsCollHBHE - ){ - - - hcalClusX_samedPhi_ = 0; hcalClusY_samedPhi_ = 0; hcalClusZ_samedPhi_ = 0; - + const CaloGeometry* geo, const reco::Photon* pho, edm::Handle hcalhitsCollHBHE) { + hcalClusX_samedPhi_ = 0; + hcalClusY_samedPhi_ = 0; + hcalClusZ_samedPhi_ = 0; + hcalClusNhits_samedPhi_ = 0; hcalClusE_samedPhi_ = 0; - - hcalClusX_samedR_ = 0; hcalClusY_samedR_ = 0; hcalClusZ_samedR_ = 0; + + hcalClusX_samedR_ = 0; + hcalClusY_samedR_ = 0; + hcalClusZ_samedR_ = 0; hcalClusNhits_samedR_ = 0; hcalClusE_samedR_ = 0; - - const HBHERecHitCollection* HBHERecHits = hcalhitsCollHBHE.product(); - + + const HBHERecHitCollection* HBHERecHits = hcalhitsCollHBHE.product(); + double phoSCEta = pho->superCluster()->eta(); double phoSCPhi = pho->superCluster()->phi(); - // Loop over HBHERecHit's - for (HBHERecHitCollection::const_iterator hbherechit = HBHERecHits->begin(); hbherechit != HBHERecHits->end(); hbherechit++) { - + for (HBHERecHitCollection::const_iterator hbherechit = HBHERecHits->begin(); hbherechit != HBHERecHits->end(); + hbherechit++) { HcalDetId det = hbherechit->id(); - const GlobalPoint & rechitPoint = geo->getPosition(det); - + const GlobalPoint& rechitPoint = geo->getPosition(det); + double rhEta = rechitPoint.eta(); double rhPhi = rechitPoint.phi(); double rhX = rechitPoint.x(); double rhY = rechitPoint.y(); double rhZ = rechitPoint.z(); double rhE = hbherechit->energy(); - + int depth = det.depth(); - + if ((det.subdet() == HcalBarrel and (depth < 1 or depth > int(recHitEThresholdHB_.size()))) or - (det.subdet() == HcalEndcap and (depth < 1 or depth > int(recHitEThresholdHE_.size())))) + (det.subdet() == HcalEndcap and (depth < 1 or depth > int(recHitEThresholdHE_.size())))) edm::LogWarning("PhotonMVABasedHaloTagger") - << " hit in subdet " << det.subdet() << " has an unaccounted for depth of " << depth << "!!"; - - const bool goodHBe = det.subdet() == HcalBarrel and rhE > recHitEThresholdHB_[depth-1]; - const bool goodHEe = det.subdet() == HcalEndcap and rhE > recHitEThresholdHE_[depth-1]; + << " hit in subdet " << det.subdet() << " has an unaccounted for depth of " << depth << "!!"; + + const bool goodHBe = det.subdet() == HcalBarrel and rhE > recHitEThresholdHB_[depth - 1]; + const bool goodHEe = det.subdet() == HcalEndcap and rhE > recHitEThresholdHE_[depth - 1]; if (!(goodHBe or goodHEe)) continue; - - if(phoSCEta * rhEta < 0) continue; ///Should be on the same side of Z - - - double rho = sqrt(pow(rhX,2) + pow(rhY,2)); + if (phoSCEta * rhEta < 0) + continue; ///Should be on the same side of Z + + double rho = sqrt(pow(rhX, 2) + pow(rhY, 2)); double dEta = fabs(phoSCEta - rhEta); double dPhi = deltaPhi(phoSCPhi, rhPhi); - - + bool isRHBehindECAL = false; - - double dRho = sqrt( pow(rhX-ecalClusX_,2) + pow(rhY-ecalClusY_,2) ) ; - - - if(rho>=31 && rho<=172 && dRho <= 26 && fabs(dPhi)<0.15){ ///only valid for the EE; this is 26 cm; hit within 3x3 of HCAL centered at the EECAL xtal + + double dRho = sqrt(pow(rhX - ecalClusX_, 2) + pow(rhY - ecalClusY_, 2)); + + if (rho >= 31 && rho <= 172 && dRho <= 26 && + fabs(dPhi) < 0.15) { ///only valid for the EE; this is 26 cm; hit within 3x3 of HCAL centered at the EECAL xtal hcalClusX_samedPhi_ += rhX * rhE; hcalClusY_samedPhi_ += rhY * rhE; hcalClusZ_samedPhi_ += rhZ * rhE; hcalClusE_samedPhi_ += rhE; hcalClusNhits_samedPhi_++; isRHBehindECAL = true; - - }//if(rho>=31 && rho<=172) - - double dR = sqrt( pow(dEta,2) + pow(dPhi,2) ); - - if(dR<0.15 && !isRHBehindECAL){ ///dont use hits which are just behind the ECAL in the same phi region + + } //if(rho>=31 && rho<=172) + + double dR = sqrt(pow(dEta, 2) + pow(dPhi, 2)); + + if (dR < 0.15 && !isRHBehindECAL) { ///dont use hits which are just behind the ECAL in the same phi region hcalClusX_samedR_ += rhX * rhE; hcalClusY_samedR_ += rhY * rhE; hcalClusZ_samedR_ += rhZ * rhE; hcalClusE_samedR_ += rhE; hcalClusNhits_samedR_++; - } - - - }//for(int ih=0; ih0){ - hcalClusX_samedPhi_ = hcalClusX_samedPhi_/hcalClusE_samedPhi_; - hcalClusY_samedPhi_ = hcalClusY_samedPhi_/hcalClusE_samedPhi_; - hcalClusZ_samedPhi_ = hcalClusZ_samedPhi_/hcalClusE_samedPhi_; - }//if(hcalClusNhits_samedPhi_>0) - - if(hcalClusNhits_samedR_>0){ - hcalClusX_samedR_ = hcalClusX_samedR_/hcalClusE_samedR_; - hcalClusY_samedR_ = hcalClusY_samedR_/hcalClusE_samedR_; - hcalClusZ_samedR_ = hcalClusZ_samedR_/hcalClusE_samedR_; - }//if(hcalClusNhits_samedR_>0) + + } //for(int ih=0; ih 0) { + hcalClusX_samedPhi_ = hcalClusX_samedPhi_ / hcalClusE_samedPhi_; + hcalClusY_samedPhi_ = hcalClusY_samedPhi_ / hcalClusE_samedPhi_; + hcalClusZ_samedPhi_ = hcalClusZ_samedPhi_ / hcalClusE_samedPhi_; + } //if(hcalClusNhits_samedPhi_>0) + + if (hcalClusNhits_samedR_ > 0) { + hcalClusX_samedR_ = hcalClusX_samedR_ / hcalClusE_samedR_; + hcalClusY_samedR_ = hcalClusY_samedR_ / hcalClusE_samedR_; + hcalClusZ_samedR_ = hcalClusZ_samedR_ / hcalClusE_samedR_; + } //if(hcalClusNhits_samedR_>0) } -void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis(const CaloGeometry* geo, - const reco::Photon *pho, - edm::Handle esrechitCollHandle - ){ - - preshowerX_samedPhi_ = 0; preshowerY_samedPhi_ = 0; preshowerZ_samedPhi_ = 0; +void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis( + const CaloGeometry* geo, const reco::Photon* pho, edm::Handle esrechitCollHandle) { + preshowerX_samedPhi_ = 0; + preshowerY_samedPhi_ = 0; + preshowerZ_samedPhi_ = 0; preshowerNhits_samedPhi_ = 0; preshowerE_samedPhi_ = 0; - preshowerX_samedR_ = 0; preshowerY_samedR_ = 0; preshowerZ_samedR_ = 0; + preshowerX_samedR_ = 0; + preshowerY_samedR_ = 0; + preshowerZ_samedR_ = 0; preshowerNhits_samedR_ = 0; preshowerE_samedR_ = 0; const EcalRecHitCollection* ESRecHits = esrechitCollHandle.product(); - + double phoSCEta = pho->superCluster()->eta(); double phoSCPhi = pho->superCluster()->phi(); - + double tmpDiffdRho = 999; double matchX_samephi = -999; double matchY_samephi = -999; bool foundESRH_samephi = false; - + double tmpDiffdRho_samedR = 999; double matchX_samedR = -999; double matchY_samedR = -999; bool foundESRH_samedR = false; ///get theta and phi of the coordinates of photon - double tan_theta = 1./sinh(phoSCEta); - + double tan_theta = 1. / sinh(phoSCEta); + double cos_phi = cos(phoSCPhi); double sin_phi = sin(phoSCPhi); - for( EcalRecHitCollection::const_iterator esrechit = ESRecHits->begin(); esrechit != ESRecHits->end(); esrechit++ ){ - - - const GlobalPoint & rechitPoint = geo->getPosition(esrechit->id()); + for (EcalRecHitCollection::const_iterator esrechit = ESRecHits->begin(); esrechit != ESRecHits->end(); esrechit++) { + const GlobalPoint& rechitPoint = geo->getPosition(esrechit->id()); double rhEta = rechitPoint.eta(); //double rhPhi = rechitPoint.phi(); @@ -334,58 +328,55 @@ void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis(const CaloGeom double rhY = rechitPoint.y(); double rhZ = rechitPoint.z(); double rhE = esrechit->energy(); - - - if(phoSCEta * rhEta < 0) continue; + if (phoSCEta * rhEta < 0) + continue; + + if (rhE < noiseThrES_) + continue; - if(rhE < noiseThrES_) continue; - /*double rho = sqrt(pow(rhX,2) + pow(rhY,2)); double dEta = fabs(phoSCEta - rhEta); double dPhi = deltaPhi(phoSCPhi, rhPhi); */ - ////try to include RH according to the strips, 11 in X and 11 in Y - /////////First calculate RH nearest in phi and eta to that of the photon SC - + ////try to include RH according to the strips, 11 in X and 11 in Y + /////////First calculate RH nearest in phi and eta to that of the photon SC + //////same phi ----> the X and Y should be similar - double dRho = sqrt( pow(rhX-ecalClusX_,2) + pow(rhY-ecalClusY_,2) ) ; + double dRho = sqrt(pow(rhX - ecalClusX_, 2) + pow(rhY - ecalClusY_, 2)); + + if (dRho < tmpDiffdRho && + dRho < 2.2) { ////i.e. hit is required to be within the ----> seems better match with the data compared to 2.47 - if(dRho < tmpDiffdRho && dRho < 2.2){ ////i.e. hit is required to be within the ----> seems better match with the data compared to 2.47 - tmpDiffdRho = dRho; matchX_samephi = rhX; matchY_samephi = rhY; foundESRH_samephi = true; } - - ////////same eta - ///calculate the expected x and y at the position of hte rechit + ////////same eta + ///calculate the expected x and y at the position of hte rechit double exp_ESRho = rhZ * tan_theta; double exp_ESX = cos_phi * exp_ESRho; double exp_ESY = sin_phi * exp_ESRho; - double dRho_samedR = sqrt( pow(rhX-exp_ESX, 2) + pow(rhY-exp_ESY,2) ); - - - if(dRho_samedR < tmpDiffdRho_samedR){ + double dRho_samedR = sqrt(pow(rhX - exp_ESX, 2) + pow(rhY - exp_ESY, 2)); + if (dRho_samedR < tmpDiffdRho_samedR) { tmpDiffdRho_samedR = dRho_samedR; matchX_samedR = rhX; matchY_samedR = rhY; foundESRH_samedR = true; } - }/// for( EcalRecHitCollection::const_iterator esrechit = ESRecHits->begin(); esrechit != ESRecHits->end(); esrechit++ ){ - - ////Now calculate the sum in +/- 5 strips in X and y around the matched RH - //+/5 strips mean = 5*~2mm = +/-10 mm = 1 cm + } /// for( EcalRecHitCollection::const_iterator esrechit = ESRecHits->begin(); esrechit != ESRecHits->end(); esrechit++ ){ - for( EcalRecHitCollection::const_iterator esrechit = ESRecHits->begin(); esrechit != ESRecHits->end(); esrechit++ ){ + ////Now calculate the sum in +/- 5 strips in X and y around the matched RH + //+/5 strips mean = 5*~2mm = +/-10 mm = 1 cm - const GlobalPoint & rechitPoint = geo->getPosition(esrechit->id()); + for (EcalRecHitCollection::const_iterator esrechit = ESRecHits->begin(); esrechit != ESRecHits->end(); esrechit++) { + const GlobalPoint& rechitPoint = geo->getPosition(esrechit->id()); double rhEta = rechitPoint.eta(); double rhPhi = rechitPoint.phi(); @@ -393,76 +384,69 @@ void PhotonMVABasedHaloTagger::calmatchedESCoordForBothHypothesis(const CaloGeom double rhY = rechitPoint.y(); double rhZ = rechitPoint.z(); double rhE = esrechit->energy(); - - - if(phoSCEta * rhEta < 0) continue; - if(rhE < noiseThrES_) continue; + if (phoSCEta * rhEta < 0) + continue; + if (rhE < noiseThrES_) + continue; bool isRHBehindECAL = false; - + ///same phi double dX_samephi = fabs(matchX_samephi - rhX); double dY_samephi = fabs(matchY_samephi - rhY); - if( (dX_samephi < 1 && dY_samephi < 1) && foundESRH_samephi){ - + if ((dX_samephi < 1 && dY_samephi < 1) && foundESRH_samephi) { preshowerX_samedPhi_ += rhX * rhE; preshowerY_samedPhi_ += rhY * rhE; preshowerZ_samedPhi_ += rhZ * rhE; - preshowerE_samedPhi_ += rhE; + preshowerE_samedPhi_ += rhE; preshowerNhits_samedPhi_++; isRHBehindECAL = true; - } - ///same dR + ///same dR double dX_samedR = fabs(matchX_samedR - rhX); double dY_samedR = fabs(matchY_samedR - rhY); - if(!isRHBehindECAL && foundESRH_samedR && (dX_samedR < 1 && dY_samedR < 1) ){ - preshowerX_samedR_ += rhX * rhE; - preshowerY_samedR_ += rhY * rhE; - preshowerZ_samedR_ += rhZ * rhE; + if (!isRHBehindECAL && foundESRH_samedR && (dX_samedR < 1 && dY_samedR < 1)) { + preshowerX_samedR_ += rhX * rhE; + preshowerY_samedR_ += rhY * rhE; + preshowerZ_samedR_ += rhZ * rhE; preshowerE_samedR_ += rhE; preshowerNhits_samedR_++; } - - - }///for(int ih=0; ih0){ - preshowerX_samedPhi_ = preshowerX_samedPhi_/preshowerE_samedPhi_; - preshowerY_samedPhi_ = preshowerY_samedPhi_/preshowerE_samedPhi_; - preshowerZ_samedPhi_ = preshowerZ_samedPhi_/preshowerE_samedPhi_; - }//if(preshowerNhits_samedPhi_>0) - - - if(preshowerNhits_samedR_>0){ - preshowerX_samedR_ = preshowerX_samedR_/preshowerE_samedR_; - preshowerY_samedR_ = preshowerY_samedR_/preshowerE_samedR_; - preshowerZ_samedR_ = preshowerZ_samedR_/preshowerE_samedR_; - }//if(preshowerNhits_samedR_>0) -} + } ///for(int ih=0; ih 0) { + preshowerX_samedPhi_ = preshowerX_samedPhi_ / preshowerE_samedPhi_; + preshowerY_samedPhi_ = preshowerY_samedPhi_ / preshowerE_samedPhi_; + preshowerZ_samedPhi_ = preshowerZ_samedPhi_ / preshowerE_samedPhi_; + } //if(preshowerNhits_samedPhi_>0) + + if (preshowerNhits_samedR_ > 0) { + preshowerX_samedR_ = preshowerX_samedR_ / preshowerE_samedR_; + preshowerY_samedR_ = preshowerY_samedR_ / preshowerE_samedR_; + preshowerZ_samedR_ = preshowerZ_samedR_ / preshowerE_samedR_; + } //if(preshowerNhits_samedR_>0) +} -double PhotonMVABasedHaloTagger::calAngleBetweenEEAndSubDet(int nhits, double subdetClusX, double subdetClusY, double subdetClusZ){ - +double PhotonMVABasedHaloTagger::calAngleBetweenEEAndSubDet(int nhits, + double subdetClusX, + double subdetClusY, + double subdetClusZ) { ////get the angle of the line joining the ECAL cluster and the subdetector wrt Z axis for any hypothesis - + double angle = -999; - - if(ecalClusNhits_>0 && nhits>0 ){ - - double dR = sqrt( pow(subdetClusX-ecalClusX_,2) + pow(subdetClusY-ecalClusY_,2) + pow(subdetClusZ-ecalClusZ_,2) ); - - double cosTheta = fabs(subdetClusZ-ecalClusZ_)/dR; - + + if (ecalClusNhits_ > 0 && nhits > 0) { + double dR = + sqrt(pow(subdetClusX - ecalClusX_, 2) + pow(subdetClusY - ecalClusY_, 2) + pow(subdetClusZ - ecalClusZ_, 2)); + + double cosTheta = fabs(subdetClusZ - ecalClusZ_) / dR; + angle = acos(cosTheta); } - + return angle; - } - - -