Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 10:00:57 

0001 // -*- C++ -*-
0002 //
0003 // Package:    METAlgorithms
0004 // Class:      MuonMETAlgo
0005 //
0006 // Original Authors:  Michael Schmitt, Richard Cavanaugh The University of Florida
0007 //          Created:  August 30, 2007
0008 //
0009 
0010 //____________________________________________________________________________||
0011 #include <cmath>
0012 #include <vector>
0013 #include "RecoMET/METAlgorithms/interface/MuonMETAlgo.h"
0014 #include "DataFormats/METReco/interface/CaloMET.h"
0015 #include "DataFormats/METReco/interface/SpecificCaloMETData.h"
0016 #include "DataFormats/TrackReco/interface/Track.h"
0017 #include "DataFormats/CaloTowers/interface/CaloTower.h"
0018 
0019 #include "TMath.h"
0020 
0021 using namespace std;
0022 using namespace reco;
0023 
0024 typedef math::XYZTLorentzVector LorentzVector;
0025 typedef math::XYZPoint Point;
0026 
0027 //____________________________________________________________________________||
0028 CaloMET MuonMETAlgo::makeMET(const CaloMET& fMet,
0029                              double fSumEt,
0030                              const std::vector<CorrMETData>& fCorrections,
0031                              const MET::LorentzVector& fP4) {
0032   return CaloMET(fMet.getSpecific(), fSumEt, fCorrections, fP4, fMet.vertex());
0033 }
0034 
0035 //____________________________________________________________________________||
0036 MET MuonMETAlgo::makeMET(const MET& fMet,
0037                          double fSumEt,
0038                          const std::vector<CorrMETData>& fCorrections,
0039                          const MET::LorentzVector& fP4) {
0040   return MET(fSumEt, fCorrections, fP4, fMet.vertex());
0041 }
0042 
0043 //____________________________________________________________________________||
0044 template <class T>
0045 void MuonMETAlgo::MuonMETAlgo_run(const edm::View<reco::Muon>& inputMuons,
0046                                   const edm::ValueMap<reco::MuonMETCorrectionData>& vm_muCorrData,
0047                                   const edm::View<T>& v_uncorMET,
0048                                   std::vector<T>* v_corMET) {
0049   const T& uncorMETObj = v_uncorMET.front();
0050 
0051   double corMETX = uncorMETObj.px();
0052   double corMETY = uncorMETObj.py();
0053 
0054   double sumMuPx = 0.;
0055   double sumMuPy = 0.;
0056   double sumMuPt = 0.;
0057   double sumMuDepEx = 0.;
0058   double sumMuDepEy = 0.;
0059   double sumMuDepEt = 0.;
0060 
0061   for (unsigned int i = 0; i < inputMuons.size(); ++i) {
0062     reco::MuonMETCorrectionData muCorrData = vm_muCorrData[inputMuons.refAt(i)];
0063 
0064     if (muCorrData.type() == 0)
0065       continue;
0066 
0067     float deltax = muCorrData.corrX();
0068     float deltay = muCorrData.corrY();
0069 
0070     const reco::Muon* mu = &inputMuons[i];
0071     const LorentzVector& mup4 = mu->p4();
0072 
0073     sumMuPx += mup4.px();
0074     sumMuPy += mup4.py();
0075     sumMuPt += mup4.pt();
0076     sumMuDepEx += deltax;
0077     sumMuDepEy += deltay;
0078     sumMuDepEt += sqrt(deltax * deltax + deltay * deltay);
0079     corMETX = corMETX - mup4.px() + deltax;
0080     corMETY = corMETY - mup4.py() + deltay;
0081   }
0082 
0083   CorrMETData delta;
0084   delta.mex = sumMuDepEx - sumMuPx;
0085   delta.mey = sumMuDepEy - sumMuPy;
0086   delta.sumet = sumMuPt - sumMuDepEt;
0087   MET::LorentzVector correctedMET4vector(corMETX, corMETY, 0., sqrt(corMETX * corMETX + corMETY * corMETY));
0088   std::vector<CorrMETData> corrections = uncorMETObj.mEtCorr();
0089   corrections.push_back(delta);
0090 
0091   T result = makeMET(uncorMETObj, uncorMETObj.sumEt() + delta.sumet, corrections, correctedMET4vector);
0092   v_corMET->push_back(result);
0093 }
0094 
0095 //____________________________________________________________________________||
0096 void MuonMETAlgo::GetMuDepDeltas(const reco::Muon* inputMuon,
0097                                  TrackDetMatchInfo& info,
0098                                  bool useTrackAssociatorPositions,
0099                                  bool useRecHits,
0100                                  bool useHO,
0101                                  double towerEtThreshold,
0102                                  double& deltax,
0103                                  double& deltay,
0104                                  double Bfield) {
0105   bool useAverage = false;
0106   //decide whether or not we want to correct on average based
0107   //on isolation information from the muon
0108   double sumPt = inputMuon->isIsolationValid() ? inputMuon->isolationR03().sumPt : 0.0;
0109   double sumEtEcal = inputMuon->isIsolationValid() ? inputMuon->isolationR03().emEt : 0.0;
0110   double sumEtHcal = inputMuon->isIsolationValid() ? inputMuon->isolationR03().hadEt : 0.0;
0111 
0112   if (sumPt > 3 || sumEtEcal + sumEtHcal > 5)
0113     useAverage = true;
0114 
0115   //get the energy using TrackAssociator if
0116   //the muon turns out to be isolated
0117   MuonMETInfo muMETInfo;
0118   muMETInfo.useAverage = useAverage;
0119   muMETInfo.useTkAssociatorPositions = useTrackAssociatorPositions;
0120   muMETInfo.useHO = useHO;
0121 
0122   TrackRef mu_track;
0123   if (inputMuon->isGlobalMuon()) {
0124     mu_track = inputMuon->globalTrack();
0125   } else if (inputMuon->isTrackerMuon() || inputMuon->isRPCMuon() || inputMuon->isGEMMuon() || inputMuon->isME0Muon()) {
0126     mu_track = inputMuon->innerTrack();
0127   } else
0128     mu_track = inputMuon->outerTrack();
0129 
0130   if (useTrackAssociatorPositions) {
0131     muMETInfo.ecalPos = info.trkGlobPosAtEcal;
0132     muMETInfo.hcalPos = info.trkGlobPosAtHcal;
0133     muMETInfo.hoPos = info.trkGlobPosAtHO;
0134   }
0135 
0136   if (!useAverage) {
0137     if (useRecHits) {
0138       muMETInfo.ecalE = inputMuon->calEnergy().emS9;
0139       muMETInfo.hcalE = inputMuon->calEnergy().hadS9;
0140       if (useHO)  //muMETInfo.hoE is 0 by default
0141         muMETInfo.hoE = inputMuon->calEnergy().hoS9;
0142     } else {  // use Towers (this is the default)
0143       //only include towers whose Et > 0.5 since
0144       //by default the MET only includes towers with Et > 0.5
0145       std::vector<const CaloTower*> towers = info.crossedTowers;
0146       for (vector<const CaloTower*>::const_iterator it = towers.begin(); it != towers.end(); it++) {
0147         if ((*it)->et() < towerEtThreshold)
0148           continue;
0149         muMETInfo.ecalE += (*it)->emEnergy();
0150         muMETInfo.hcalE += (*it)->hadEnergy();
0151         if (useHO)
0152           muMETInfo.hoE += (*it)->outerEnergy();
0153       }
0154     }  //use Towers
0155   }
0156 
0157   //This needs to be fixed!!!!!
0158   //The tracker has better resolution for pt < 200 GeV
0159   math::XYZTLorentzVector mup4;
0160   if (inputMuon->isGlobalMuon()) {
0161     if (inputMuon->globalTrack()->pt() < 200) {
0162       mup4 = LorentzVector(inputMuon->innerTrack()->px(),
0163                            inputMuon->innerTrack()->py(),
0164                            inputMuon->innerTrack()->pz(),
0165                            inputMuon->innerTrack()->p());
0166     } else {
0167       mup4 = LorentzVector(inputMuon->globalTrack()->px(),
0168                            inputMuon->globalTrack()->py(),
0169                            inputMuon->globalTrack()->pz(),
0170                            inputMuon->globalTrack()->p());
0171     }
0172   } else if (inputMuon->isTrackerMuon() || inputMuon->isRPCMuon() || inputMuon->isGEMMuon() || inputMuon->isME0Muon()) {
0173     mup4 = LorentzVector(inputMuon->innerTrack()->px(),
0174                          inputMuon->innerTrack()->py(),
0175                          inputMuon->innerTrack()->pz(),
0176                          inputMuon->innerTrack()->p());
0177   } else
0178     mup4 = LorentzVector(inputMuon->outerTrack()->px(),
0179                          inputMuon->outerTrack()->py(),
0180                          inputMuon->outerTrack()->pz(),
0181                          inputMuon->outerTrack()->p());
0182 
0183   //call function that does the work
0184   correctMETforMuon(deltax, deltay, Bfield, inputMuon->charge(), mup4, inputMuon->vertex(), muMETInfo);
0185 }
0186 
0187 //____________________________________________________________________________||
0188 void MuonMETAlgo::correctMETforMuon(double& deltax,
0189                                     double& deltay,
0190                                     double bfield,
0191                                     int muonCharge,
0192                                     const math::XYZTLorentzVector& muonP4,
0193                                     const math::XYZPoint& muonVertex,
0194                                     MuonMETInfo& muonMETInfo) {
0195   double mu_p = muonP4.P();
0196   double mu_pt = muonP4.Pt();
0197   double mu_phi = muonP4.Phi();
0198   double mu_eta = muonP4.Eta();
0199   double mu_vz = muonVertex.z() / 100.;
0200   double mu_pz = muonP4.Pz();
0201 
0202   double ecalPhi, ecalTheta;
0203   double hcalPhi, hcalTheta;
0204   double hoPhi, hoTheta;
0205 
0206   //should always be false for FWLite
0207   //unless you want to supply co-ordinates at
0208   //the calorimeter sub-detectors yourself
0209   if (muonMETInfo.useTkAssociatorPositions) {
0210     ecalPhi = muonMETInfo.ecalPos.Phi();
0211     ecalTheta = muonMETInfo.ecalPos.Theta();
0212     hcalPhi = muonMETInfo.hcalPos.Phi();
0213     hcalTheta = muonMETInfo.hcalPos.Theta();
0214     hoPhi = muonMETInfo.hoPos.Phi();
0215     hoTheta = muonMETInfo.hoPos.Theta();
0216   } else {
0217     /*
0218       use the analytical solution for the
0219       intersection of a helix with a cylinder
0220       to find the positions of the muon
0221       at the various calo surfaces
0222     */
0223 
0224     //radii of subdetectors in meters
0225     double rEcal = 1.290;
0226     double rHcal = 1.9;
0227     double rHo = 3.82;
0228     if (abs(mu_eta) > 0.3)
0229       rHo = 4.07;
0230     //distance from the center of detector to face of Ecal
0231     double zFaceEcal = 3.209;
0232     if (mu_eta < 0)
0233       zFaceEcal = -1 * zFaceEcal;
0234     //distance from the center of detector to face of Hcal
0235     double zFaceHcal = 3.88;
0236     if (mu_eta < 0)
0237       zFaceHcal = -1 * zFaceHcal;
0238 
0239     //now we have to get Phi
0240     //bending radius of the muon (units are meters)
0241     double bendr = mu_pt * 1000 / (300 * bfield);
0242 
0243     double tb_ecal = TMath::ACos(1 - rEcal * rEcal / (2 * bendr * bendr));  //helix time interval parameter
0244     double tb_hcal = TMath::ACos(1 - rHcal * rHcal / (2 * bendr * bendr));  //helix time interval parameter
0245     double tb_ho = TMath::ACos(1 - rHo * rHo / (2 * bendr * bendr));        //helix time interval parameter
0246     double xEcal, yEcal, zEcal;
0247     double xHcal, yHcal, zHcal;
0248     double xHo, yHo, zHo;
0249     //Ecal
0250     //in the barrel and if not a looper
0251     if (fabs(mu_pz * bendr * tb_ecal / mu_pt + mu_vz) < fabs(zFaceEcal) && rEcal < 2 * bendr) {
0252       xEcal = bendr * (TMath::Sin(tb_ecal + mu_phi) - TMath::Sin(mu_phi));
0253       yEcal = bendr * (-TMath::Cos(tb_ecal + mu_phi) + TMath::Cos(mu_phi));
0254       zEcal = bendr * tb_ecal * mu_pz / mu_pt + mu_vz;
0255     } else {  //endcap
0256       if (mu_pz > 0) {
0257         double te_ecal = (fabs(zFaceEcal) - mu_vz) * mu_pt / (bendr * mu_pz);
0258         xEcal = bendr * (TMath::Sin(te_ecal + mu_phi) - TMath::Sin(mu_phi));
0259         yEcal = bendr * (-TMath::Cos(te_ecal + mu_phi) + TMath::Cos(mu_phi));
0260         zEcal = fabs(zFaceEcal);
0261       } else {
0262         double te_ecal = -(fabs(zFaceEcal) + mu_vz) * mu_pt / (bendr * mu_pz);
0263         xEcal = bendr * (TMath::Sin(te_ecal + mu_phi) - TMath::Sin(mu_phi));
0264         yEcal = bendr * (-TMath::Cos(te_ecal + mu_phi) + TMath::Cos(mu_phi));
0265         zEcal = -fabs(zFaceEcal);
0266       }
0267     }
0268 
0269     //Hcal
0270     if (fabs(mu_pz * bendr * tb_hcal / mu_pt + mu_vz) < fabs(zFaceHcal) && rEcal < 2 * bendr) {  //in the barrel
0271       xHcal = bendr * (TMath::Sin(tb_hcal + mu_phi) - TMath::Sin(mu_phi));
0272       yHcal = bendr * (-TMath::Cos(tb_hcal + mu_phi) + TMath::Cos(mu_phi));
0273       zHcal = bendr * tb_hcal * mu_pz / mu_pt + mu_vz;
0274     } else {  //endcap
0275       if (mu_pz > 0) {
0276         double te_hcal = (fabs(zFaceHcal) - mu_vz) * mu_pt / (bendr * mu_pz);
0277         xHcal = bendr * (TMath::Sin(te_hcal + mu_phi) - TMath::Sin(mu_phi));
0278         yHcal = bendr * (-TMath::Cos(te_hcal + mu_phi) + TMath::Cos(mu_phi));
0279         zHcal = fabs(zFaceHcal);
0280       } else {
0281         double te_hcal = -(fabs(zFaceHcal) + mu_vz) * mu_pt / (bendr * mu_pz);
0282         xHcal = bendr * (TMath::Sin(te_hcal + mu_phi) - TMath::Sin(mu_phi));
0283         yHcal = bendr * (-TMath::Cos(te_hcal + mu_phi) + TMath::Cos(mu_phi));
0284         zHcal = -fabs(zFaceHcal);
0285       }
0286     }
0287 
0288     //Ho - just the barrel
0289     xHo = bendr * (TMath::Sin(tb_ho + mu_phi) - TMath::Sin(mu_phi));
0290     yHo = bendr * (-TMath::Cos(tb_ho + mu_phi) + TMath::Cos(mu_phi));
0291     zHo = bendr * tb_ho * mu_pz / mu_pt + mu_vz;
0292 
0293     ecalTheta = TMath::ACos(zEcal / sqrt(pow(xEcal, 2) + pow(yEcal, 2) + pow(zEcal, 2)));
0294     ecalPhi = atan2(yEcal, xEcal);
0295     hcalTheta = TMath::ACos(zHcal / sqrt(pow(xHcal, 2) + pow(yHcal, 2) + pow(zHcal, 2)));
0296     hcalPhi = atan2(yHcal, xHcal);
0297     hoTheta = TMath::ACos(zHo / sqrt(pow(xHo, 2) + pow(yHo, 2) + pow(zHo, 2)));
0298     hoPhi = atan2(yHo, xHo);
0299 
0300     /*
0301       the above prescription is for right handed helicies only
0302       Positively charged muons trace a left handed helix
0303       so we correct for that 
0304     */
0305     if (muonCharge > 0) {
0306       //Ecal
0307       double dphi = mu_phi - ecalPhi;
0308       if (fabs(dphi) > TMath::Pi())
0309         dphi = 2 * TMath::Pi() - fabs(dphi);
0310       ecalPhi = mu_phi - fabs(dphi);
0311       if (fabs(ecalPhi) > TMath::Pi()) {
0312         double temp = 2 * TMath::Pi() - fabs(ecalPhi);
0313         ecalPhi = -1 * temp * ecalPhi / fabs(ecalPhi);
0314       }
0315 
0316       //Hcal
0317       dphi = mu_phi - hcalPhi;
0318       if (fabs(dphi) > TMath::Pi())
0319         dphi = 2 * TMath::Pi() - fabs(dphi);
0320       hcalPhi = mu_phi - fabs(dphi);
0321       if (fabs(hcalPhi) > TMath::Pi()) {
0322         double temp = 2 * TMath::Pi() - fabs(hcalPhi);
0323         hcalPhi = -1 * temp * hcalPhi / fabs(hcalPhi);
0324       }
0325 
0326       //Ho
0327       dphi = mu_phi - hoPhi;
0328       if (fabs(dphi) > TMath::Pi())
0329         dphi = 2 * TMath::Pi() - fabs(dphi);
0330       hoPhi = mu_phi - fabs(dphi);
0331       if (fabs(hoPhi) > TMath::Pi()) {
0332         double temp = 2 * TMath::Pi() - fabs(hoPhi);
0333         hoPhi = -1 * temp * hoPhi / fabs(hoPhi);
0334       }
0335     }
0336   }
0337 
0338   //for isolated muons
0339   if (!muonMETInfo.useAverage) {
0340     double mu_Ex = muonMETInfo.ecalE * sin(ecalTheta) * cos(ecalPhi) +
0341                    muonMETInfo.hcalE * sin(hcalTheta) * cos(hcalPhi) + muonMETInfo.hoE * sin(hoTheta) * cos(hoPhi);
0342     double mu_Ey = muonMETInfo.ecalE * sin(ecalTheta) * sin(ecalPhi) +
0343                    muonMETInfo.hcalE * sin(hcalTheta) * sin(hcalPhi) + muonMETInfo.hoE * sin(hoTheta) * sin(hoPhi);
0344 
0345     deltax += mu_Ex;
0346     deltay += mu_Ey;
0347 
0348   } else {  //non-isolated muons - derive the correction
0349 
0350     //dE/dx in matter for iron:
0351     //-(11.4 + 0.96*fabs(log(p0*2.8)) + 0.033*p0*(1.0 - pow(p0, -0.33)) )*1e-3
0352     //from http://cmslxr.fnal.gov/lxr/source/TrackPropagation/SteppingHelixPropagator/src/SteppingHelixPropagator.ccyes,
0353     //line ~1100
0354     //normalisation is at 50 GeV
0355     double dEdx_normalization = -(11.4 + 0.96 * fabs(log(50 * 2.8)) + 0.033 * 50 * (1.0 - pow(50, -0.33))) * 1e-3;
0356     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;
0357 
0358     double temp = 0.0;
0359 
0360     if (muonMETInfo.useHO) {
0361       //for the Towers, with HO
0362       if (fabs(mu_eta) < 0.2)
0363         temp = 2.75 * (1 - 0.00003 * mu_p);
0364       if (fabs(mu_eta) > 0.2 && fabs(mu_eta) < 1.0)
0365         temp = (2.38 + 0.0144 * fabs(mu_eta)) * (1 - 0.0003 * mu_p);
0366       if (fabs(mu_eta) > 1.0 && fabs(mu_eta) < 1.3)
0367         temp = 7.413 - 5.12 * fabs(mu_eta);
0368       if (fabs(mu_eta) > 1.3)
0369         temp = 2.084 - 0.743 * fabs(mu_eta);
0370     } else {
0371       if (fabs(mu_eta) < 1.0)
0372         temp = 2.33 * (1 - 0.0004 * mu_p);
0373       if (fabs(mu_eta) > 1.0 && fabs(mu_eta) < 1.3)
0374         temp = (7.413 - 5.12 * fabs(mu_eta)) * (1 - 0.0003 * mu_p);
0375       if (fabs(mu_eta) > 1.3)
0376         temp = 2.084 - 0.743 * fabs(mu_eta);
0377     }
0378 
0379     double dep = temp * dEdx_normalization / dEdx_numerator;
0380     if (dep < 0.5)
0381       dep = 0;
0382     //use the average phi of the 3 subdetectors
0383     if (fabs(mu_eta) < 1.3) {
0384       deltax += dep * cos((ecalPhi + hcalPhi + hoPhi) / 3);
0385       deltay += dep * sin((ecalPhi + hcalPhi + hoPhi) / 3);
0386     } else {
0387       deltax += dep * cos((ecalPhi + hcalPhi) / 2);
0388       deltay += dep * cos((ecalPhi + hcalPhi) / 2);
0389     }
0390   }
0391 }
0392 
0393 //____________________________________________________________________________||
0394 void MuonMETAlgo::run(const edm::View<reco::Muon>& inputMuons,
0395                       const edm::ValueMap<reco::MuonMETCorrectionData>& vm_muCorrData,
0396                       const edm::View<reco::MET>& uncorMET,
0397                       METCollection* corMET) {
0398   MuonMETAlgo_run(inputMuons, vm_muCorrData, uncorMET, corMET);
0399 }
0400 
0401 //____________________________________________________________________________||
0402 void MuonMETAlgo::run(const edm::View<reco::Muon>& inputMuons,
0403                       const edm::ValueMap<reco::MuonMETCorrectionData>& vm_muCorrData,
0404                       const edm::View<reco::CaloMET>& uncorMET,
0405                       CaloMETCollection* corMET) {
0406   MuonMETAlgo_run(inputMuons, vm_muCorrData, uncorMET, corMET);
0407 }
0408 
0409 //____________________________________________________________________________||
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!