File indexing completed on 2024-04-06 12:26:40
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010 #include "RecoMET/METAlgorithms/interface/CaloSpecificAlgo.h"
0011 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0012 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0013 #include "DataFormats/CaloTowers/interface/CaloTower.h"
0014
0015 #include <iostream>
0016
0017 using namespace reco;
0018
0019
0020
0021
0022
0023
0024
0025
0026 reco::CaloMET CaloSpecificAlgo::addInfo(edm::Handle<edm::View<Candidate> > towers,
0027 const CommonMETData &umet,
0028 bool noHF,
0029 double globalThreshold) {
0030 SpecificCaloMETData specific;
0031 CommonMETData met = umet;
0032 double totalEt = 0.0;
0033 double totalEm = 0.0;
0034
0035 double sumEtInpHF = 0.0;
0036 double sumEtInmHF = 0.0;
0037 double MExInpHF = 0.0;
0038 double MEyInpHF = 0.0;
0039 double MExInmHF = 0.0;
0040 double MEyInmHF = 0.0;
0041
0042 for (edm::View<Candidate>::const_iterator towerCand = towers->begin(); towerCand != towers->end(); ++towerCand) {
0043 const CaloTower *calotower = dynamic_cast<const CaloTower *>(&(*towerCand));
0044 if (!calotower)
0045 continue;
0046 if (calotower->et() < globalThreshold)
0047 continue;
0048 update_totalEt_totalEm(totalEt, totalEm, calotower, noHF);
0049 update_MaxTowerEm_MaxTowerHad(specific.MaxEtInEmTowers, specific.MaxEtInHadTowers, calotower, noHF);
0050 update_EmEtInEB_EmEtInEE(specific.EmEtInEB, specific.EmEtInEE, calotower);
0051 update_HadEtInHB_HadEtInHE_HadEtInHO_HadEtInHF_EmEtInHF(specific.HadEtInHB,
0052 specific.HadEtInHE,
0053 specific.HadEtInHO,
0054 specific.HadEtInHF,
0055 specific.EmEtInHF,
0056 calotower,
0057 noHF);
0058 update_sumEtInpHF_MExInpHF_MEyInpHF_sumEtInmHF_MExInmHF_MEyInmHF(
0059 sumEtInpHF, MExInpHF, MEyInpHF, sumEtInmHF, MExInmHF, MEyInmHF, calotower);
0060 }
0061
0062 double totalHad = totalEt - totalEm;
0063
0064 if (noHF)
0065 remove_HF_from_MET(met, sumEtInpHF, MExInpHF, MEyInpHF, sumEtInmHF, MExInmHF, MEyInmHF);
0066
0067 if (!noHF)
0068 add_MET_in_HF(specific, sumEtInpHF, MExInpHF, MEyInpHF, sumEtInmHF, MExInmHF, MEyInmHF);
0069
0070 specific.EtFractionHadronic = (totalEt == 0.0) ? 0.0 : totalHad / totalEt;
0071 specific.EtFractionEm = (totalEt == 0.0) ? 0.0 : totalEm / totalEt;
0072
0073 const LorentzVector p4(met.mex, met.mey, 0.0, met.met);
0074 const Point vtx(0.0, 0.0, 0.0);
0075 CaloMET caloMET(specific, met.sumet, p4, vtx);
0076 return caloMET;
0077 }
0078
0079
0080 void CaloSpecificAlgo::update_totalEt_totalEm(double &totalEt, double &totalEm, const CaloTower *calotower, bool noHF) {
0081 if (noHF) {
0082 DetId detIdHcal = find_DetId_of_HCAL_cell_in_constituent_of(calotower);
0083 if (!detIdHcal.null()) {
0084 HcalSubdetector subdet = HcalDetId(detIdHcal).subdet();
0085 if (subdet == HcalForward)
0086 return;
0087 }
0088 }
0089
0090 totalEt += calotower->et();
0091 totalEm += calotower->emEt();
0092 }
0093
0094
0095 void CaloSpecificAlgo::update_MaxTowerEm_MaxTowerHad(float &MaxTowerEm,
0096 float &MaxTowerHad,
0097 const CaloTower *calotower,
0098 bool noHF) {
0099 DetId detIdHcal = find_DetId_of_HCAL_cell_in_constituent_of(calotower);
0100 DetId detIdEcal = find_DetId_of_ECAL_cell_in_constituent_of(calotower);
0101
0102 if (!detIdHcal.null()) {
0103 HcalSubdetector subdet = HcalDetId(detIdHcal).subdet();
0104 if (subdet == HcalBarrel || subdet == HcalOuter || subdet == HcalEndcap || (!noHF && subdet == HcalForward)) {
0105 if (calotower->hadEt() > MaxTowerHad)
0106 MaxTowerHad = calotower->hadEt();
0107 if (calotower->emEt() > MaxTowerEm)
0108 MaxTowerEm = calotower->emEt();
0109 }
0110 }
0111
0112 if (!detIdEcal.null()) {
0113 if (calotower->emEt() > MaxTowerEm)
0114 MaxTowerEm = calotower->emEt();
0115 }
0116 }
0117
0118
0119 void CaloSpecificAlgo::update_EmEtInEB_EmEtInEE(float &EmEtInEB, float &EmEtInEE, const CaloTower *calotower) {
0120 DetId detIdEcal = find_DetId_of_ECAL_cell_in_constituent_of(calotower);
0121 if (detIdEcal.null())
0122 return;
0123
0124 EcalSubdetector subdet = EcalSubdetector(detIdEcal.subdetId());
0125 if (subdet == EcalBarrel) {
0126 EmEtInEB += calotower->emEt();
0127 } else if (subdet == EcalEndcap) {
0128 EmEtInEE += calotower->emEt();
0129 }
0130 }
0131
0132
0133 void CaloSpecificAlgo::update_HadEtInHB_HadEtInHE_HadEtInHO_HadEtInHF_EmEtInHF(float &HadEtInHB,
0134 float &HadEtInHE,
0135 float &HadEtInHO,
0136 float &HadEtInHF,
0137 float &EmEtInHF,
0138 const CaloTower *calotower,
0139 bool noHF) {
0140 DetId detIdHcal = find_DetId_of_HCAL_cell_in_constituent_of(calotower);
0141 if (detIdHcal.null())
0142 return;
0143
0144 HcalSubdetector subdet = HcalDetId(detIdHcal).subdet();
0145 if (subdet == HcalBarrel || subdet == HcalOuter) {
0146 HadEtInHB += calotower->hadEt();
0147 HadEtInHO += calotower->outerEt();
0148 }
0149
0150 if (subdet == HcalEndcap) {
0151 HadEtInHE += calotower->hadEt();
0152 }
0153
0154 if (subdet == HcalForward && !noHF) {
0155 HadEtInHF += calotower->hadEt();
0156 EmEtInHF += calotower->emEt();
0157 }
0158 }
0159
0160
0161 void CaloSpecificAlgo::update_sumEtInpHF_MExInpHF_MEyInpHF_sumEtInmHF_MExInmHF_MEyInmHF(double &sumEtInpHF,
0162 double &MExInpHF,
0163 double &MEyInpHF,
0164 double &sumEtInmHF,
0165 double &MExInmHF,
0166 double &MEyInmHF,
0167 const CaloTower *calotower) {
0168 DetId detIdHcal = find_DetId_of_HCAL_cell_in_constituent_of(calotower);
0169 if (detIdHcal.null())
0170 return;
0171
0172 HcalSubdetector subdet = HcalDetId(detIdHcal).subdet();
0173 if (!(subdet == HcalForward))
0174 return;
0175
0176 if (calotower->eta() >= 0) {
0177 sumEtInpHF += calotower->et();
0178 MExInpHF -= (calotower->et() * cos(calotower->phi()));
0179 MEyInpHF -= (calotower->et() * sin(calotower->phi()));
0180 } else {
0181 sumEtInmHF += calotower->et();
0182 MExInmHF -= (calotower->et() * cos(calotower->phi()));
0183 MEyInmHF -= (calotower->et() * sin(calotower->phi()));
0184 }
0185 }
0186
0187
0188 void CaloSpecificAlgo::remove_HF_from_MET(CommonMETData &met,
0189 double sumEtInpHF,
0190 double MExInpHF,
0191 double MEyInpHF,
0192 double sumEtInmHF,
0193 double MExInmHF,
0194 double MEyInmHF) {
0195 met.mex -= (MExInmHF + MExInpHF);
0196 met.mey -= (MEyInmHF + MEyInpHF);
0197 met.sumet -= (sumEtInpHF + sumEtInmHF);
0198 met.met = sqrt(met.mex * met.mex + met.mey * met.mey);
0199 }
0200
0201
0202 void CaloSpecificAlgo::add_MET_in_HF(SpecificCaloMETData &specific,
0203 double sumEtInpHF,
0204 double MExInpHF,
0205 double MEyInpHF,
0206 double sumEtInmHF,
0207 double MExInmHF,
0208 double MEyInmHF) {
0209 LorentzVector METpHF(MExInpHF, MEyInpHF, 0, sqrt(MExInpHF * MExInpHF + MEyInpHF * MEyInpHF));
0210 LorentzVector METmHF(MExInmHF, MEyInmHF, 0, sqrt(MExInmHF * MExInmHF + MEyInmHF * MEyInmHF));
0211 specific.CaloMETInpHF = METpHF.pt();
0212 specific.CaloMETInmHF = METmHF.pt();
0213 specific.CaloMETPhiInpHF = METpHF.Phi();
0214 specific.CaloMETPhiInmHF = METmHF.Phi();
0215 specific.CaloSETInpHF = sumEtInpHF;
0216 specific.CaloSETInmHF = sumEtInmHF;
0217 }
0218
0219
0220 DetId CaloSpecificAlgo::find_DetId_of_HCAL_cell_in_constituent_of(const CaloTower *calotower) {
0221 DetId ret;
0222 for (int cell = calotower->constituentsSize() - 1; cell >= 0; --cell) {
0223 DetId id = calotower->constituent(cell);
0224 if (id.det() == DetId::Hcal) {
0225 ret = id;
0226 break;
0227 }
0228 }
0229 return ret;
0230 }
0231
0232
0233 DetId CaloSpecificAlgo::find_DetId_of_ECAL_cell_in_constituent_of(const CaloTower *calotower) {
0234 DetId ret;
0235 for (int cell = calotower->constituentsSize() - 1; cell >= 0; --cell) {
0236 DetId id = calotower->constituent(cell);
0237 if (id.det() == DetId::Ecal) {
0238 ret = id;
0239 break;
0240 }
0241 }
0242 return ret;
0243 }
0244
0245
