RecoTau/plugins/RecoTauElectronRejectionPlugin.cc

0001 /*
0002  * ===========================================================================
0003  *
0004  *       Filename:  RecoTauElectronRejectionPlugin.cc
0005  *
0006  *    Description:  Add electron rejection information to PFTau
0007  *
0008  *         Authors:  Chi Nhan Nguyen, Simone Gennai, Evan Friis
0009  *
0010  * ===========================================================================
0011  */
0012
0013 #include "RecoTauTag/RecoTau/interface/RecoTauBuilderPlugins.h"
0014 #include "DataFormats/TrackReco/interface/Track.h"
0015 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0016 #include "DataFormats/ParticleFlowReco/interface/PFBlockElement.h"
0017 #include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
0018 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0019 #include <Math/VectorUtil.h>
0020 #include <algorithm>
0021
0022 namespace reco {
0023   namespace tau {
0024
0025     class RecoTauElectronRejectionPlugin : public RecoTauModifierPlugin {
0026     public:
0027       explicit RecoTauElectronRejectionPlugin(const edm::ParameterSet& pset, edm::ConsumesCollector&& iC);
0028       ~RecoTauElectronRejectionPlugin() override {}
0029       void operator()(PFTau&) const override;
0030
0031     private:
0032       double ElecPreIDLeadTkMatch_maxDR_;
0033       double EcalStripSumE_minClusEnergy_;
0034       double EcalStripSumE_deltaEta_;
0035       double EcalStripSumE_deltaPhiOverQ_minValue_;
0036       double EcalStripSumE_deltaPhiOverQ_maxValue_;
0037       double maximumForElectrionPreIDOutput_;
0038       std::string DataType_;
0039     };
0040
0041     RecoTauElectronRejectionPlugin::RecoTauElectronRejectionPlugin(const edm::ParameterSet& pset,
0042                                                                    edm::ConsumesCollector&& iC)
0043         : RecoTauModifierPlugin(pset, std::move(iC)) {
0044       // Load parameters
0045       ElecPreIDLeadTkMatch_maxDR_ = pset.getParameter<double>("ElecPreIDLeadTkMatch_maxDR");
0046       EcalStripSumE_minClusEnergy_ = pset.getParameter<double>("EcalStripSumE_minClusEnergy");
0047       EcalStripSumE_deltaEta_ = pset.getParameter<double>("EcalStripSumE_deltaEta");
0048       EcalStripSumE_deltaPhiOverQ_minValue_ = pset.getParameter<double>("EcalStripSumE_deltaPhiOverQ_minValue");
0049       EcalStripSumE_deltaPhiOverQ_maxValue_ = pset.getParameter<double>("EcalStripSumE_deltaPhiOverQ_maxValue");
0050       maximumForElectrionPreIDOutput_ = pset.getParameter<double>("maximumForElectrionPreIDOutput");
0051       DataType_ = pset.getParameter<std::string>("DataType");
0052     }
0053
0054     namespace {
0055       bool checkPos(std::vector<math::XYZPoint> CalPos, math::XYZPoint CandPos) {
0056         bool flag = false;
0057         for (unsigned int i = 0; i < CalPos.size(); i++) {
0058           if (CalPos[i] == CandPos) {
0059             flag = true;
0060             break;
0061           }
0062         }
0063         return flag;
0064       }
0065     }  // namespace
0066
0067     void RecoTauElectronRejectionPlugin::operator()(PFTau& tau) const {
0068       // copy pasted from PFRecoTauAlgorithm...
0069       double myECALenergy = 0.;
0070       double myHCALenergy = 0.;
0071       double myHCALenergy3x3 = 0.;
0072       double myMaximumHCALPFClusterE = 0.;
0073       double myMaximumHCALPFClusterEt = 0.;
0074       double myStripClusterE = 0.;
0075       double myEmfrac = -1.;
0076       double myElectronPreIDOutput = -1111.;
0077       bool myElecPreid = false;
0078       reco::TrackRef myElecTrk;
0079
0080       typedef std::pair<reco::PFBlockRef, unsigned> ElementInBlock;
0081       typedef std::vector<ElementInBlock> ElementsInBlocks;
0082
0083       PFCandidatePtr myleadPFChargedCand = tau.leadPFChargedHadrCand();
0084       // Build list of PFCands in tau
0085       std::vector<PFCandidatePtr> myPFCands;
0086       myPFCands.reserve(tau.isolationPFCands().size() + tau.signalPFCands().size());
0087
0088       std::copy(tau.isolationPFCands().begin(), tau.isolationPFCands().end(), std::back_inserter(myPFCands));
0089       std::copy(tau.signalPFCands().begin(), tau.signalPFCands().end(), std::back_inserter(myPFCands));
0090
0091       //Use the electron rejection only in case there is a charged leading pion
0092       if (myleadPFChargedCand.isNonnull()) {
0093         myElectronPreIDOutput = myleadPFChargedCand->mva_e_pi();
0094
0095         math::XYZPointF myElecTrkEcalPos = myleadPFChargedCand->positionAtECALEntrance();
0096         myElecTrk = myleadPFChargedCand->trackRef();  //Electron candidate
0097
0098         if (myElecTrk.isNonnull()) {
0099           //FROM AOD
0100           if (DataType_ == "AOD") {
0101             // Corrected Cluster energies
0102             for (int i = 0; i < (int)myPFCands.size(); i++) {
0103               myHCALenergy += myPFCands[i]->hcalEnergy();
0104               myECALenergy += myPFCands[i]->ecalEnergy();
0105
0106               math::XYZPointF candPos;
0107               if (myPFCands[i]->particleId() == 1 || myPFCands[i]->particleId() == 2)  //if charged hadron or electron
0108                 candPos = myPFCands[i]->positionAtECALEntrance();
0109               else
0110                 candPos = math::XYZPointF(myPFCands[i]->px(), myPFCands[i]->py(), myPFCands[i]->pz());
0111
0112               double deltaR = ROOT::Math::VectorUtil::DeltaR(myElecTrkEcalPos, candPos);
0113               double deltaPhi = ROOT::Math::VectorUtil::DeltaPhi(myElecTrkEcalPos, candPos);
0114               double deltaEta = std::abs(myElecTrkEcalPos.eta() - candPos.eta());
0115               double deltaPhiOverQ = deltaPhi / (double)myElecTrk->charge();
0116
0117               if (myPFCands[i]->ecalEnergy() >= EcalStripSumE_minClusEnergy_ && deltaEta < EcalStripSumE_deltaEta_ &&
0118                   deltaPhiOverQ > EcalStripSumE_deltaPhiOverQ_minValue_ &&
0119                   deltaPhiOverQ < EcalStripSumE_deltaPhiOverQ_maxValue_) {
0120                 myStripClusterE += myPFCands[i]->ecalEnergy();
0121               }
0122               if (deltaR < 0.184) {
0123                 myHCALenergy3x3 += myPFCands[i]->hcalEnergy();
0124               }
0125               if (myPFCands[i]->hcalEnergy() > myMaximumHCALPFClusterE) {
0126                 myMaximumHCALPFClusterE = myPFCands[i]->hcalEnergy();
0127               }
0128               if ((myPFCands[i]->hcalEnergy() * fabs(sin(candPos.Theta()))) > myMaximumHCALPFClusterEt) {
0129                 myMaximumHCALPFClusterEt = (myPFCands[i]->hcalEnergy() * fabs(sin(candPos.Theta())));
0130               }
0131             }
0132
0133           } else if (DataType_ == "RECO") {  //From RECO
0134             // Against double counting of clusters
0135             std::vector<math::XYZPoint> hcalPosV;
0136             hcalPosV.clear();
0137             std::vector<math::XYZPoint> ecalPosV;
0138             ecalPosV.clear();
0139             for (int i = 0; i < (int)myPFCands.size(); i++) {
0140               const ElementsInBlocks& elts = myPFCands[i]->elementsInBlocks();
0141               for (ElementsInBlocks::const_iterator it = elts.begin(); it != elts.end(); ++it) {
0142                 const reco::PFBlock& block = *(it->first);
0143                 unsigned indexOfElementInBlock = it->second;
0144                 const edm::OwnVector<reco::PFBlockElement>& elements = block.elements();
0145                 assert(indexOfElementInBlock < elements.size());
0146
0147                 const reco::PFBlockElement& element = elements[indexOfElementInBlock];
0148
0149                 if (element.type() == reco::PFBlockElement::HCAL) {
0150                   math::XYZPoint clusPos = element.clusterRef()->position();
0151                   double en = (double)element.clusterRef()->energy();
0152                   double et = (double)element.clusterRef()->energy() * fabs(sin(clusPos.Theta()));
0153                   if (en > myMaximumHCALPFClusterE) {
0154                     myMaximumHCALPFClusterE = en;
0155                   }
0156                   if (et > myMaximumHCALPFClusterEt) {
0157                     myMaximumHCALPFClusterEt = et;
0158                   }
0159                   if (!checkPos(hcalPosV, clusPos)) {
0160                     hcalPosV.push_back(clusPos);
0161                     myHCALenergy += en;
0162                     double deltaR = ROOT::Math::VectorUtil::DeltaR(myElecTrkEcalPos, clusPos);
0163                     if (deltaR < 0.184) {
0164                       myHCALenergy3x3 += en;
0165                     }
0166                   }
0167                 } else if (element.type() == reco::PFBlockElement::ECAL) {
0168                   double en = (double)element.clusterRef()->energy();
0169                   math::XYZPoint clusPos = element.clusterRef()->position();
0170                   if (!checkPos(ecalPosV, clusPos)) {
0171                     ecalPosV.push_back(clusPos);
0172                     myECALenergy += en;
0173                     double deltaPhi = ROOT::Math::VectorUtil::DeltaPhi(myElecTrkEcalPos, clusPos);
0174                     double deltaEta = std::abs(myElecTrkEcalPos.eta() - clusPos.eta());
0175                     double deltaPhiOverQ = deltaPhi / (double)myElecTrk->charge();
0176                     if (en >= EcalStripSumE_minClusEnergy_ && deltaEta < EcalStripSumE_deltaEta_ &&
0177                         deltaPhiOverQ > EcalStripSumE_deltaPhiOverQ_minValue_ &&
0178                         deltaPhiOverQ < EcalStripSumE_deltaPhiOverQ_maxValue_) {
0179                       myStripClusterE += en;
0180                     }
0181                   }
0182                 }
0183               }  //end elements in blocks
0184             }    //end loop over PFcands
0185           }      //end RECO case
0186         }        // end check for null electrk
0187       }          // end check for null pfChargedHadrCand
0188
0189       if ((myHCALenergy + myECALenergy) > 0.)
0190         myEmfrac = myECALenergy / (myHCALenergy + myECALenergy);
0191       tau.setemFraction((float)myEmfrac);
0192
0193       // scale the appropriate quantities by the momentum of the electron if it exists
0194       if (myElecTrk.isNonnull()) {
0195         float myElectronMomentum = (float)myElecTrk->p();
0196         if (myElectronMomentum > 0.) {
0197           myHCALenergy /= myElectronMomentum;
0198           myMaximumHCALPFClusterE /= myElectronMomentum;
0199           myHCALenergy3x3 /= myElectronMomentum;
0200           myStripClusterE /= myElectronMomentum;
0201         }
0202       }
0203       tau.sethcalTotOverPLead((float)myHCALenergy);
0204       tau.sethcalMaxOverPLead((float)myMaximumHCALPFClusterE);
0205       tau.sethcal3x3OverPLead((float)myHCALenergy3x3);
0206       tau.setecalStripSumEOverPLead((float)myStripClusterE);
0207       tau.setmaximumHCALPFClusterEt(myMaximumHCALPFClusterEt);
0208       tau.setelectronPreIDOutput(myElectronPreIDOutput);
0209       if (myElecTrk.isNonnull())
0210         tau.setelectronPreIDTrack(myElecTrk);
0211       if (myElectronPreIDOutput > maximumForElectrionPreIDOutput_)
0212         myElecPreid = true;
0213       tau.setelectronPreIDDecision(myElecPreid);
0214
0215       // These need to be filled!
0216       //tau.setbremsRecoveryEOverPLead(my...);
0217
0218       /* End elecron rejection */
0219     }
0220   }  // namespace tau
0221 }  // namespace reco
0222 #include "FWCore/Framework/interface/MakerMacros.h"
0223 DEFINE_EDM_PLUGIN(RecoTauModifierPluginFactory,
0224                   reco::tau::RecoTauElectronRejectionPlugin,
0225                   "RecoTauElectronRejectionPlugin");