File indexing completed on 2024-04-06 12:19:25
0001 #if !defined(__CINT__) && !defined(__MAKECINT__)
0002 #include "DataFormats/FWLite/interface/Handle.h"
0003 #include "DataFormats/FWLite/interface/Event.h"
0004 #include "DataFormats/METReco/interface/CaloMET.h"
0005 #include "DataFormats/METReco/interface/METCollection.h"
0006 #include "DataFormats/METReco/interface/MET.h"
0007 #include "DataFormats/METReco/interface/CaloMETCollection.h"
0008 #include "DataFormats/METReco/interface/SpecificCaloMETData.h"
0009 #include "DataFormats/JetReco/interface/CaloJet.h"
0010 #include "DataFormats/TrackReco/interface/Track.h"
0011 #include "DataFormats/MuonReco/interface/Muon.h"
0012 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0013 #endif
0014
0015 namespace reco {
0016 class Muon;
0017 class CaloMET;
0018 };
0019 #include "DataFormats/Math/interface/Point3D.h"
0020 #include "DataFormats/Math/interface/LorentzVector.h"
0021 #include "RecoMET/METAlgorithms/interface/MuonMETInfo.h"
0022
0023 #include <vector>
0024 #include "TMath.h"
0025 #include "TFile.h"
0026
0027
0028
0029 void correctMETforMuon(double& metx, double& mety, double bfield, int muonCharge,
0030 math::XYZTLorentzVector muonP4, math::XYZPoint muonVertex,
0031 MuonMETInfo& muonMETInfo) {
0032
0033
0034
0035 double mu_p = muonP4.P();
0036 double mu_pt = muonP4.Pt();
0037 double mu_phi = muonP4.Phi();
0038 double mu_eta = muonP4.Eta();
0039 double mu_vz = muonVertex.z()/100.;
0040 double mu_pz = muonP4.Pz();
0041
0042
0043
0044 metx -= mu_pt*cos(mu_phi);
0045 mety -= mu_pt*sin(mu_phi);
0046
0047 double ecalPhi, ecalTheta;
0048 double hcalPhi, hcalTheta;
0049 double hoPhi, hoTheta;
0050
0051
0052
0053
0054
0055 if(muonMETInfo.useTkAssociatorPositions) {
0056 ecalPhi = muonMETInfo.ecalPos.Phi();
0057 ecalTheta = muonMETInfo.ecalPos.Theta();
0058 hcalPhi = muonMETInfo.hcalPos.Phi();
0059 hcalTheta = muonMETInfo.hcalPos.Theta();
0060 hoPhi = muonMETInfo.hoPos.Phi();
0061 hoTheta = muonMETInfo.hoPos.Theta();
0062 } else {
0063
0064
0065
0066
0067
0068
0069
0070 double rEcal = 1.290;
0071 double rHcal = 1.9;
0072 double rHo = 3.82;
0073 if(abs(mu_eta) > 0.3) rHo = 4.07;
0074
0075 double zFaceEcal = 3.209;
0076 if(mu_eta < 0 ) zFaceEcal = -1*zFaceEcal;
0077
0078 double zFaceHcal = 3.88;
0079 if(mu_eta < 0 ) zFaceHcal = -1*zFaceHcal;
0080
0081
0082
0083 double bendr = mu_pt*1000/(300*bfield);
0084
0085 double tb_ecal = TMath::ACos(1-rEcal*rEcal/(2*bendr*bendr));
0086 double tb_hcal = TMath::ACos(1-rHcal*rHcal/(2*bendr*bendr));
0087 double tb_ho = TMath::ACos(1-rHo*rHo/(2*bendr*bendr));
0088 double xEcal,yEcal,zEcal;
0089 double xHcal,yHcal,zHcal;
0090 double xHo, yHo,zHo;
0091
0092
0093 if(fabs(mu_pz*bendr*tb_ecal/mu_pt+mu_vz) < fabs(zFaceEcal) && rEcal < 2*bendr) {
0094 xEcal = bendr*(TMath::Sin(tb_ecal+mu_phi)-TMath::Sin(mu_phi));
0095 yEcal = bendr*(-TMath::Cos(tb_ecal+mu_phi)+TMath::Cos(mu_phi));
0096 zEcal = bendr*tb_ecal*mu_pz/mu_pt + mu_vz;
0097 } else {
0098 if(mu_pz > 0) {
0099 double te_ecal = (fabs(zFaceEcal) - mu_vz)*mu_pt/(bendr*mu_pz);
0100 xEcal = bendr*(TMath::Sin(te_ecal+mu_phi) - TMath::Sin(mu_phi));
0101 yEcal = bendr*(-TMath::Cos(te_ecal+mu_phi) + TMath::Cos(mu_phi));
0102 zEcal = fabs(zFaceEcal);
0103 } else {
0104 double te_ecal = -(fabs(zFaceEcal) + mu_vz)*mu_pt/(bendr*mu_pz);
0105 xEcal = bendr*(TMath::Sin(te_ecal+mu_phi) - TMath::Sin(mu_phi));
0106 yEcal = bendr*(-TMath::Cos(te_ecal+mu_phi) + TMath::Cos(mu_phi));
0107 zEcal = -fabs(zFaceEcal);
0108 }
0109 }
0110
0111
0112 if(fabs(mu_pz*bendr*tb_hcal/mu_pt+mu_vz) < fabs(zFaceHcal) && rEcal < 2*bendr) {
0113 xHcal = bendr*(TMath::Sin(tb_hcal+mu_phi)-TMath::Sin(mu_phi));
0114 yHcal = bendr*(-TMath::Cos(tb_hcal+mu_phi)+TMath::Cos(mu_phi));
0115 zHcal = bendr*tb_hcal*mu_pz/mu_pt + mu_vz;
0116 } else {
0117 if(mu_pz > 0) {
0118 double te_hcal = (fabs(zFaceHcal) - mu_vz)*mu_pt/(bendr*mu_pz);
0119 xHcal = bendr*(TMath::Sin(te_hcal+mu_phi) - TMath::Sin(mu_phi));
0120 yHcal = bendr*(-TMath::Cos(te_hcal+mu_phi) + TMath::Cos(mu_phi));
0121 zHcal = fabs(zFaceHcal);
0122 } else {
0123 double te_hcal = -(fabs(zFaceHcal) + mu_vz)*mu_pt/(bendr*mu_pz);
0124 xHcal = bendr*(TMath::Sin(te_hcal+mu_phi) - TMath::Sin(mu_phi));
0125 yHcal = bendr*(-TMath::Cos(te_hcal+mu_phi) + TMath::Cos(mu_phi));
0126 zHcal = -fabs(zFaceHcal);
0127 }
0128 }
0129
0130
0131 xHo = bendr*(TMath::Sin(tb_ho+mu_phi)-TMath::Sin(mu_phi));
0132 yHo = bendr*(-TMath::Cos(tb_ho+mu_phi)+TMath::Cos(mu_phi));
0133 zHo = bendr*tb_ho*mu_pz/mu_pt + mu_vz;
0134
0135 ecalTheta = TMath::ACos(zEcal/sqrt(pow(xEcal,2) + pow(yEcal,2)+pow(zEcal,2)));
0136 ecalPhi = atan2(yEcal,xEcal);
0137 hcalTheta = TMath::ACos(zHcal/sqrt(pow(xHcal,2) + pow(yHcal,2)+pow(zHcal,2)));
0138 hcalPhi = atan2(yHcal,xHcal);
0139 hoTheta = TMath::ACos(zHo/sqrt(pow(xHo,2) + pow(yHo,2)+pow(zHo,2)));
0140 hoPhi = atan2(yHo,xHo);
0141
0142
0143 double r2dEcal = sqrt(pow(xEcal,2)+pow(yEcal,2));
0144 double r2dHcal = sqrt(pow(xHcal,2)+pow(yHcal,2));
0145 double r2dHo = sqrt(pow(xHo,2) +pow(yHo,2));
0146
0147
0148
0149
0150
0151
0152 if(muonCharge > 0) {
0153
0154
0155 double dphi = mu_phi - ecalPhi;
0156 if(fabs(dphi) > TMath::Pi())
0157 dphi = 2*TMath::Pi() - fabs(dphi);
0158 ecalPhi = mu_phi - fabs(dphi);
0159 if(fabs(ecalPhi) > TMath::Pi()) {
0160 double temp = 2*TMath::Pi() - fabs(ecalPhi);
0161 ecalPhi = -1*temp*ecalPhi/fabs(ecalPhi);
0162 }
0163 xEcal = r2dEcal*TMath::Cos(ecalPhi);
0164 yEcal = r2dEcal*TMath::Sin(ecalPhi);
0165
0166
0167 dphi = mu_phi - hcalPhi;
0168 if(fabs(dphi) > TMath::Pi())
0169 dphi = 2*TMath::Pi() - fabs(dphi);
0170 hcalPhi = mu_phi - fabs(dphi);
0171 if(fabs(hcalPhi) > TMath::Pi()) {
0172 double temp = 2*TMath::Pi() - fabs(hcalPhi);
0173 hcalPhi = -1*temp*hcalPhi/fabs(hcalPhi);
0174 }
0175 xHcal = r2dHcal*TMath::Cos(hcalPhi);
0176 yHcal = r2dHcal*TMath::Sin(hcalPhi);
0177
0178
0179 dphi = mu_phi - hoPhi;
0180 if(fabs(dphi) > TMath::Pi())
0181 dphi = 2*TMath::Pi() - fabs(dphi);
0182 hoPhi = mu_phi - fabs(dphi);
0183 if(fabs(hoPhi) > TMath::Pi()) {
0184 double temp = 2*TMath::Pi() - fabs(hoPhi);
0185 hoPhi = -1*temp*hoPhi/fabs(hoPhi);
0186 }
0187 xHo = r2dHo*TMath::Cos(hoPhi);
0188 yHo = r2dHo*TMath::Sin(hoPhi);
0189
0190 }
0191 }
0192
0193
0194 if(!muonMETInfo.useHO) muonMETInfo.hoE = 0;
0195
0196
0197 if(!muonMETInfo.useAverage) {
0198
0199 double mu_Ex = muonMETInfo.ecalE*sin(ecalTheta)*cos(ecalPhi)
0200 + muonMETInfo.hcalE*sin(hcalTheta)*cos(hcalPhi)
0201 + muonMETInfo.hoE*sin(hoTheta)*cos(hoPhi);
0202 double mu_Ey = muonMETInfo.ecalE*sin(ecalTheta)*sin(ecalPhi)
0203 + muonMETInfo.hcalE*sin(hcalTheta)*sin(hcalPhi)
0204 + muonMETInfo.hoE*sin(hoTheta)*sin(hoPhi);
0205
0206 metx += mu_Ex;
0207 mety += mu_Ey;
0208
0209 } else {
0210
0211
0212
0213
0214
0215
0216 double dEdx_normalization = -(11.4 + 0.96*fabs(log(50*2.8)) + 0.033*50*(1.0 - pow(50, -0.33)) )*1e-3;
0217 double dEdx_numerator = -(11.4 + 0.96*fabs(log(mu_p*2.8)) + 0.033*mu_p*(1.0 - pow(mu_p, -0.33)) )*1e-3;
0218
0219 double temp = 0.0;
0220
0221 if(muonMETInfo.useHO) {
0222
0223 if(fabs(mu_eta) < 0.2)
0224 temp = 2.75*(1-0.00003*mu_p);
0225 if(fabs(mu_eta) > 0.2 && fabs(mu_eta) < 1.0)
0226 temp = (2.38+0.0144*fabs(mu_eta))*(1-0.0003*mu_p);
0227 if(fabs(mu_eta) > 1.0 && fabs(mu_eta) < 1.3)
0228 temp = 7.413-5.12*fabs(mu_eta);
0229 if(fabs(mu_eta) > 1.3)
0230 temp = 2.084-0.743*fabs(mu_eta);
0231 } else {
0232 if(fabs(mu_eta) < 1.0)
0233 temp = 2.33*(1-0.0004*mu_p);
0234 if(fabs(mu_eta) > 1.0 && fabs(mu_eta) < 1.3)
0235 temp = (7.413-5.12*fabs(mu_eta))*(1-0.0003*mu_p);
0236 if(fabs(mu_eta) > 1.3)
0237 temp = 2.084-0.743*fabs(mu_eta);
0238 }
0239
0240 double dep = temp*dEdx_normalization/dEdx_numerator;
0241 if(dep < 0.5) dep = 0.0;
0242
0243 if(fabs(mu_eta) < 1.3) {
0244 metx += dep*cos((ecalPhi+hcalPhi+hoPhi)/3);
0245 mety += dep*sin((ecalPhi+hcalPhi+hoPhi)/3);
0246 } else {
0247 metx += dep*cos( (ecalPhi+hcalPhi)/2);
0248 mety += dep*cos( (ecalPhi+hcalPhi)/2);
0249 }
0250 }
0251
0252
0253 }
0254
0255
0256
0257 void correctMETforMuon(double& metx, double&mety, double bfield,
0258 const reco::Muon* muon, MuonMETInfo& muonMETInfo){
0259
0260
0261 math::XYZTLorentzVector mup4;
0262 if(muon->combinedMuon()->pt() < 200) {
0263 mup4 = math::XYZTLorentzVector(muon->innerTrack()->px(), muon->innerTrack()->py(),
0264 muon->innerTrack()->pz(), muon->innerTrack()->p());
0265 } else {
0266 mup4 = math::XYZTLorentzVector(muon->combinedMuon()->px(), muon->combinedMuon()->py(),
0267 muon->combinedMuon()->pz(), muon->combinedMuon()->p());
0268 }
0269
0270 correctMETforMuon(metx, mety, bfield, muon->charge(),
0271 mup4, muon->vertex(),
0272 muonMETInfo);
0273 }
0274
0275
0276
0277 void print_data() {
0278 TFile file("0AC27F01-D260-DD11-8207-0018F3D096EE.root");
0279
0280 fwlite::Event ev(&file);
0281
0282 for( ev.toBegin();
0283 ! ev.atEnd();
0284 ++ev) {
0285 fwlite::Handle<std::vector<reco::Muon> > muonCollection;
0286 muonCollection.getByLabel(ev,"muons");
0287
0288 fwlite::Handle<std::vector<reco::CaloMET> > metCollection;
0289 metCollection.getByLabel(ev,"met");
0290
0291
0292
0293
0294
0295
0296 reco::CaloMETCollection::const_iterator metIt = metCollection->begin();
0297 double metx = metIt->px();
0298 double mety = metIt->py();
0299
0300
0301 for(reco::MuonCollection::const_iterator it = muonCollection->begin();
0302 it != muonCollection->end(); it++) {
0303
0304 if(it->isGlobalMuon() ==0 ) continue;
0305
0306
0307 std::cout << "Before Correcting for the muon of pt: " << it->pt()
0308 << " the MET is: " << sqrt(pow(metx,2)+pow(mety,2)) << endl;
0309
0310 MuonMETInfo muonMETinfo;
0311 muonMETinfo.ecalE = it->isEnergyValid() ? it->calEnergy().emS9 : 0.0;
0312 muonMETinfo.hcalE = it->isEnergyValid() ? it->calEnergy().hadS9: 0.0;
0313 muonMETinfo.hoE = it->isEnergyValid() ? it->calEnergy().hoS9 : 0.0;
0314
0315
0316
0317 muonMETinfo.useTkAssociatorPositions = false;
0318
0319
0320
0321
0322 correctMETforMuon(metx,mety,3.8112,&(*it), muonMETinfo);
0323
0324 std::cout << "After correcting for the muon, the MET is: "
0325 << sqrt(pow(metx,2)+pow(mety,2)) << endl;
0326
0327 }
0328
0329 }
0330 }