CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
RecoTauElectronRejectionPlugin.cc
Go to the documentation of this file.
1 /*
2  * ===========================================================================
3  *
4  * Filename: RecoTauElectronRejectionPlugin.cc
5  *
6  * Description: Add electron rejection information to PFTau
7  *
8  * Authors: Chi Nhan Nguyen, Simone Gennai, Evan Friis
9  *
10  * ===========================================================================
11  */
12 
19 #include <Math/VectorUtil.h>
20 #include <algorithm>
21 
22 namespace reco {
23  namespace tau {
24 
26  public:
29  void operator()(PFTau&) const override;
30 
31  private:
39  };
40 
44  // Load parameters
45  ElecPreIDLeadTkMatch_maxDR_ = pset.getParameter<double>("ElecPreIDLeadTkMatch_maxDR");
46  EcalStripSumE_minClusEnergy_ = pset.getParameter<double>("EcalStripSumE_minClusEnergy");
47  EcalStripSumE_deltaEta_ = pset.getParameter<double>("EcalStripSumE_deltaEta");
48  EcalStripSumE_deltaPhiOverQ_minValue_ = pset.getParameter<double>("EcalStripSumE_deltaPhiOverQ_minValue");
49  EcalStripSumE_deltaPhiOverQ_maxValue_ = pset.getParameter<double>("EcalStripSumE_deltaPhiOverQ_maxValue");
50  maximumForElectrionPreIDOutput_ = pset.getParameter<double>("maximumForElectrionPreIDOutput");
51  DataType_ = pset.getParameter<std::string>("DataType");
52  }
53 
54  namespace {
55  bool checkPos(std::vector<math::XYZPoint> CalPos, math::XYZPoint CandPos) {
56  bool flag = false;
57  for (unsigned int i = 0; i < CalPos.size(); i++) {
58  if (CalPos[i] == CandPos) {
59  flag = true;
60  break;
61  }
62  }
63  return flag;
64  }
65  } // namespace
66 
68  // copy pasted from PFRecoTauAlgorithm...
69  double myECALenergy = 0.;
70  double myHCALenergy = 0.;
71  double myHCALenergy3x3 = 0.;
72  double myMaximumHCALPFClusterE = 0.;
73  double myMaximumHCALPFClusterEt = 0.;
74  double myStripClusterE = 0.;
75  double myEmfrac = -1.;
76  double myElectronPreIDOutput = -1111.;
77  bool myElecPreid = false;
78  reco::TrackRef myElecTrk;
79 
80  typedef std::pair<reco::PFBlockRef, unsigned> ElementInBlock;
81  typedef std::vector<ElementInBlock> ElementsInBlocks;
82 
83  PFCandidatePtr myleadPFChargedCand = tau.leadPFChargedHadrCand();
84  // Build list of PFCands in tau
85  std::vector<PFCandidatePtr> myPFCands;
86  myPFCands.reserve(tau.isolationPFCands().size() + tau.signalPFCands().size());
87 
88  std::copy(tau.isolationPFCands().begin(), tau.isolationPFCands().end(), std::back_inserter(myPFCands));
89  std::copy(tau.signalPFCands().begin(), tau.signalPFCands().end(), std::back_inserter(myPFCands));
90 
91  //Use the electron rejection only in case there is a charged leading pion
92  if (myleadPFChargedCand.isNonnull()) {
93  myElectronPreIDOutput = myleadPFChargedCand->mva_e_pi();
94 
95  math::XYZPointF myElecTrkEcalPos = myleadPFChargedCand->positionAtECALEntrance();
96  myElecTrk = myleadPFChargedCand->trackRef(); //Electron candidate
97 
98  if (myElecTrk.isNonnull()) {
99  //FROM AOD
100  if (DataType_ == "AOD") {
101  // Corrected Cluster energies
102  for (int i = 0; i < (int)myPFCands.size(); i++) {
103  myHCALenergy += myPFCands[i]->hcalEnergy();
104  myECALenergy += myPFCands[i]->ecalEnergy();
105 
106  math::XYZPointF candPos;
107  if (myPFCands[i]->particleId() == 1 || myPFCands[i]->particleId() == 2) //if charged hadron or electron
108  candPos = myPFCands[i]->positionAtECALEntrance();
109  else
110  candPos = math::XYZPointF(myPFCands[i]->px(), myPFCands[i]->py(), myPFCands[i]->pz());
111 
112  double deltaR = ROOT::Math::VectorUtil::DeltaR(myElecTrkEcalPos, candPos);
113  double deltaPhi = ROOT::Math::VectorUtil::DeltaPhi(myElecTrkEcalPos, candPos);
114  double deltaEta = std::abs(myElecTrkEcalPos.eta() - candPos.eta());
115  double deltaPhiOverQ = deltaPhi / (double)myElecTrk->charge();
116 
118  deltaPhiOverQ > EcalStripSumE_deltaPhiOverQ_minValue_ &&
119  deltaPhiOverQ < EcalStripSumE_deltaPhiOverQ_maxValue_) {
120  myStripClusterE += myPFCands[i]->ecalEnergy();
121  }
122  if (deltaR < 0.184) {
123  myHCALenergy3x3 += myPFCands[i]->hcalEnergy();
124  }
125  if (myPFCands[i]->hcalEnergy() > myMaximumHCALPFClusterE) {
126  myMaximumHCALPFClusterE = myPFCands[i]->hcalEnergy();
127  }
128  if ((myPFCands[i]->hcalEnergy() * fabs(sin(candPos.Theta()))) > myMaximumHCALPFClusterEt) {
129  myMaximumHCALPFClusterEt = (myPFCands[i]->hcalEnergy() * fabs(sin(candPos.Theta())));
130  }
131  }
132 
133  } else if (DataType_ == "RECO") { //From RECO
134  // Against double counting of clusters
135  std::vector<math::XYZPoint> hcalPosV;
136  hcalPosV.clear();
137  std::vector<math::XYZPoint> ecalPosV;
138  ecalPosV.clear();
139  for (int i = 0; i < (int)myPFCands.size(); i++) {
140  const ElementsInBlocks& elts = myPFCands[i]->elementsInBlocks();
141  for (ElementsInBlocks::const_iterator it = elts.begin(); it != elts.end(); ++it) {
142  const reco::PFBlock& block = *(it->first);
143  unsigned indexOfElementInBlock = it->second;
145  assert(indexOfElementInBlock < elements.size());
146 
147  const reco::PFBlockElement& element = elements[indexOfElementInBlock];
148 
149  if (element.type() == reco::PFBlockElement::HCAL) {
150  math::XYZPoint clusPos = element.clusterRef()->position();
151  double en = (double)element.clusterRef()->energy();
152  double et = (double)element.clusterRef()->energy() * fabs(sin(clusPos.Theta()));
153  if (en > myMaximumHCALPFClusterE) {
154  myMaximumHCALPFClusterE = en;
155  }
156  if (et > myMaximumHCALPFClusterEt) {
157  myMaximumHCALPFClusterEt = et;
158  }
159  if (!checkPos(hcalPosV, clusPos)) {
160  hcalPosV.push_back(clusPos);
161  myHCALenergy += en;
162  double deltaR = ROOT::Math::VectorUtil::DeltaR(myElecTrkEcalPos, clusPos);
163  if (deltaR < 0.184) {
164  myHCALenergy3x3 += en;
165  }
166  }
167  } else if (element.type() == reco::PFBlockElement::ECAL) {
168  double en = (double)element.clusterRef()->energy();
169  math::XYZPoint clusPos = element.clusterRef()->position();
170  if (!checkPos(ecalPosV, clusPos)) {
171  ecalPosV.push_back(clusPos);
172  myECALenergy += en;
173  double deltaPhi = ROOT::Math::VectorUtil::DeltaPhi(myElecTrkEcalPos, clusPos);
174  double deltaEta = std::abs(myElecTrkEcalPos.eta() - clusPos.eta());
175  double deltaPhiOverQ = deltaPhi / (double)myElecTrk->charge();
177  deltaPhiOverQ > EcalStripSumE_deltaPhiOverQ_minValue_ &&
178  deltaPhiOverQ < EcalStripSumE_deltaPhiOverQ_maxValue_) {
179  myStripClusterE += en;
180  }
181  }
182  }
183  } //end elements in blocks
184  } //end loop over PFcands
185  } //end RECO case
186  } // end check for null electrk
187  } // end check for null pfChargedHadrCand
188 
189  if ((myHCALenergy + myECALenergy) > 0.)
190  myEmfrac = myECALenergy / (myHCALenergy + myECALenergy);
191  tau.setemFraction((float)myEmfrac);
192 
193  // scale the appropriate quantities by the momentum of the electron if it exists
194  if (myElecTrk.isNonnull()) {
195  float myElectronMomentum = (float)myElecTrk->p();
196  if (myElectronMomentum > 0.) {
197  myHCALenergy /= myElectronMomentum;
198  myMaximumHCALPFClusterE /= myElectronMomentum;
199  myHCALenergy3x3 /= myElectronMomentum;
200  myStripClusterE /= myElectronMomentum;
201  }
202  }
203  tau.sethcalTotOverPLead((float)myHCALenergy);
204  tau.sethcalMaxOverPLead((float)myMaximumHCALPFClusterE);
205  tau.sethcal3x3OverPLead((float)myHCALenergy3x3);
206  tau.setecalStripSumEOverPLead((float)myStripClusterE);
207  tau.setmaximumHCALPFClusterEt(myMaximumHCALPFClusterEt);
208  tau.setelectronPreIDOutput(myElectronPreIDOutput);
209  if (myElecTrk.isNonnull())
210  tau.setelectronPreIDTrack(myElecTrk);
211  if (myElectronPreIDOutput > maximumForElectrionPreIDOutput_)
212  myElecPreid = true;
213  tau.setelectronPreIDDecision(myElecPreid);
214 
215  // These need to be filled!
216  //tau.setbremsRecoveryEOverPLead(my...);
217 
218  /* End elecron rejection */
219  }
220  } // namespace tau
221 } // namespace reco
225  "RecoTauElectronRejectionPlugin");
constexpr double deltaPhi(double phi1, double phi2)
Definition: deltaPhi.h:26
Abstract base class for a PFBlock element (track, cluster...)
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:238
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< float > > XYZPointF
point in space with cartesian internal representation
Definition: Point3D.h:10
assert(be >=bs)
static const double deltaEta
Definition: CaloConstants.h:8
std::vector< ElementInBlock > ElementsInBlocks
bool isNonnull() const
Checks for non-null.
Definition: Ptr.h:146
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
constexpr auto deltaR(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
Definition: deltaR.h:30
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
std::pair< reco::PFBlockRef, unsigned > ElementInBlock
RecoTauElectronRejectionPlugin(const edm::ParameterSet &pset, edm::ConsumesCollector &&iC)
fixed size matrix
#define DEFINE_EDM_PLUGIN(factory, type, name)
def move(src, dest)
Definition: eostools.py:511
Block of elements.
Definition: PFBlock.h:26