CMSSW Reference Manual

00001 // File: MuonMETAlgo.cc
00002 // Description:  see MuonMETAlgo.h
00003 // Author: M. Schmitt, R. Cavanaugh, The University of Florida
00004 // Creation Date:  MHS May 31, 2005 Initial version.
00005 //
00006 //--------------------------------------------
00007 #include <math.h>
00008 #include <vector>
00009 #include "JetMETCorrections/Type1MET/interface/MuonMETAlgo.h"
00010 #include "DataFormats/METReco/interface/CaloMET.h"
00011 #include "DataFormats/METReco/interface/SpecificCaloMETData.h"
00012 #include "DataFormats/JetReco/interface/CaloJet.h"
00013 #include "DataFormats/TrackReco/interface/Track.h"
00014 #include "DataFormats/CaloTowers/interface/CaloTower.h"
00015 #include "FWCore/Framework/interface/ESHandle.h"
00016 #include "Geometry/Records/interface/IdealGeometryRecord.h"
00017 #include "MagneticField/Engine/interface/MagneticField.h"
00018 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00019 #include "DataFormats/Math/interface/Point3D.h"
00020 #include "TMath.h"
00021 
00022 
00023 
00024 using namespace std;
00025 using namespace reco;
00026 
00027 typedef math::XYZTLorentzVector LorentzVector;
00028 typedef math::XYZPoint Point;
00029 
00030 CaloMET MuonMETAlgo::makeMET (const CaloMET& fMet, 
00031                               double fSumEt, 
00032                               const vector<CorrMETData>& fCorrections, 
00033                               const MET::LorentzVector& fP4) {
00034   return CaloMET (fMet.getSpecific (), fSumEt, fCorrections, fP4, fMet.vertex ());
00035 }
00036   
00037   
00038   
00039 MET MuonMETAlgo::makeMET (const MET& fMet, 
00040                           double fSumEt, 
00041                           const vector<CorrMETData>& fCorrections, 
00042                           const MET::LorentzVector& fP4) {
00043   return MET (fSumEt, fCorrections, fP4, fMet.vertex ()); 
00044 }
00045 
00046 template <class T> void MuonMETAlgo::MuonMETAlgo_run(const edm::Event& iEvent,
00047                                                      const edm::EventSetup& iSetup, 
00048                                                      const edm::View<T>& v_uncorMET,
00049                                                      const edm::View<Muon>& inputMuons,
00050                                                      TrackDetectorAssociator& trackAssociator,
00051                                                      TrackAssociatorParameters& trackAssociatorParameters,
00052                                                      vector<T>* v_corMET,
00053                                                      bool useTrackAssociatorPositions,
00054                                                      bool useRecHits,
00055                                                      bool useHO,
00056                                                      double towerEtThreshold) {
00057   
00058   
00059   edm::ESHandle<MagneticField> magneticField;
00060   iSetup.get<IdealMagneticFieldRecord>().get(magneticField);
00061     
00062 //get the B-field at the origin
00063   double Bfield = magneticField->inTesla(GlobalPoint(0.,0.,0.)).z();
00064   T uncorMETObj = v_uncorMET.front();
00065     
00066   double uncorMETX = uncorMETObj.px();
00067   double uncorMETY = uncorMETObj.py();
00068     
00069   double corMETX = uncorMETX;
00070   double corMETY = uncorMETY;
00071   
00072   CorrMETData delta;
00073   double sumMuEt=0;
00074   double sumMuDepEt = 0;
00075   for(edm::View<Muon>::const_iterator mus_it = inputMuons.begin();
00076       mus_it != inputMuons.end(); mus_it++) {
00077 
00078     bool useAverage = false;
00079     //decide whether or not we want to correct on average based 
00080     //on isolation information from the muon
00081     double sumPt   = mus_it->isIsolationValid()? mus_it->isolationR03().sumPt       : 0.0;
00082     double sumEtEcal = mus_it->isIsolationValid() ? mus_it->isolationR03().emEt     : 0.0;
00083     double sumEtHcal    = mus_it->isIsolationValid() ? mus_it->isolationR03().hadEt : 0.0;
00084     
00085     if(sumPt > 3 || sumEtEcal + sumEtHcal > 5) useAverage = true;
00086 
00087     //get the energy using TrackAssociator if
00088     //the muon turns out to be isolated
00089     MuonMETInfo muMETInfo;
00090     muMETInfo.useAverage = useAverage;
00091     muMETInfo.useTkAssociatorPositions = useTrackAssociatorPositions;
00092     muMETInfo.useHO = useHO;
00093 
00094     TrackRef mu_track = mus_it->globalTrack();
00095     TrackDetMatchInfo info = 
00096       trackAssociator.associate(iEvent, iSetup,
00097                                 trackAssociator.getFreeTrajectoryState(iSetup, *mu_track),
00098                                 trackAssociatorParameters);
00099     if(useTrackAssociatorPositions) {
00100       muMETInfo.ecalPos  = info.trkGlobPosAtEcal;
00101       muMETInfo.hcalPos  = info.trkGlobPosAtHcal;
00102       muMETInfo.hoPos    = info.trkGlobPosAtHO;
00103     }
00104     
00105     if(!useAverage) {
00106 
00107       if(useRecHits) {
00108         muMETInfo.ecalE = mus_it->calEnergy().emS9;
00109         muMETInfo.hcalE = mus_it->calEnergy().hadS9;
00110         if(useHO) //muMETInfo.hoE is 0 by default
00111           muMETInfo.hoE   = mus_it->calEnergy().hoS9;
00112       } else {// use Towers (this is the default)
00113         //only include towers whose Et > 0.5 since 
00114         //by default the MET only includes towers with Et > 0.5
00115         vector<const CaloTower*> towers = info.crossedTowers;
00116         for(vector<const CaloTower*>::const_iterator it = towers.begin();
00117             it != towers.end(); it++) {
00118           if( (*it)->et() <  towerEtThreshold) continue;
00119           muMETInfo.ecalE =+ (*it)->emEt();
00120           muMETInfo.hcalE =+ (*it)->hadEt();
00121           if(useHO)
00122             muMETInfo.hoE =+ (*it)->outerEt();
00123         }
00124       }//use Towers
00125     }
00126   
00127     
00128     sumMuEt += mus_it->et();
00129     double metxBeforeCorr = corMETX;
00130     double metyBeforeCorr = corMETY;
00131     
00132     
00133     //The tracker has better resolution for pt < 200 GeV
00134     math::XYZTLorentzVector mup4;
00135     if(mus_it->globalTrack()->pt() < 200) {
00136       mup4 = LorentzVector(mus_it->innerTrack()->px(), mus_it->innerTrack()->py(),
00137                            mus_it->innerTrack()->pz(), mus_it->innerTrack()->p());
00138     } else {
00139       mup4 = LorentzVector(mus_it->globalTrack()->px(), mus_it->globalTrack()->py(),
00140                            mus_it->globalTrack()->pz(), mus_it->globalTrack()->p());
00141     }   
00142     
00143     //call function that does the work 
00144     correctMETforMuon(corMETX, corMETY, Bfield, mus_it->charge(),
00145                       mup4, mus_it->vertex(),
00146                       muMETInfo);
00147     double sumMuDepEx = corMETX - metxBeforeCorr + mus_it->px();
00148     double sumMuDepEy = corMETY - metyBeforeCorr + mus_it->py();
00149     sumMuDepEt += sqrt(pow(sumMuDepEx,2)+pow(sumMuDepEy,2));
00150     
00151   }
00152   
00153   double corMET   = sqrt(corMETX*corMETX   +   corMETY*corMETY    );
00154   delta.mex   = corMETX - uncorMETX;
00155   delta.mey   = corMETY - uncorMETY;
00156   //remove muon's contribution to sumEt and add in the 
00157   //muon Et from track
00158   delta.sumet = sumMuEt - sumMuDepEt;
00159     
00160   MET::LorentzVector correctedMET4vector( corMETX, corMETY, 0., corMET);
00161   //----------------- get previous corrections and push into new corrections
00162   std::vector<CorrMETData> corrections = uncorMETObj.mEtCorr();
00163   corrections.push_back( delta );
00164   
00165   T result = makeMET(uncorMETObj, uncorMETObj.sumEt()+delta.sumet, corrections, correctedMET4vector);
00166   v_corMET->push_back(result);
00167   
00168 }
00169    
00170 //----------------------------------------------------------------------------
00171 
00172 void MuonMETAlgo::correctMETforMuon(double& metx, double& mety, double bfield, int muonCharge,
00173                                     math::XYZTLorentzVector muonP4, math::XYZPoint muonVertex,
00174                                     MuonMETInfo& muonMETInfo) {
00175   
00176   double mu_p     = muonP4.P();
00177   double mu_pt    = muonP4.Pt();
00178   double mu_phi   = muonP4.Phi();
00179   double mu_eta   = muonP4.Eta();
00180   double mu_vz    = muonVertex.z()/100.;
00181   double mu_pz    = muonP4.Pz();
00182   
00183   
00184   
00185   //add in the muon's pt
00186   metx -= mu_pt*cos(mu_phi);
00187   mety -= mu_pt*sin(mu_phi);
00188 
00189   double ecalPhi, ecalTheta;
00190   double hcalPhi, hcalTheta;
00191   double hoPhi, hoTheta;
00192   
00193 
00194   //should always be false for FWLite
00195   //unless you want to supply co-ordinates at 
00196   //the calorimeter sub-detectors yourself
00197   if(muonMETInfo.useTkAssociatorPositions) {
00198     ecalPhi   = muonMETInfo.ecalPos.Phi();
00199     ecalTheta = muonMETInfo.ecalPos.Theta();
00200     hcalPhi   = muonMETInfo.hcalPos.Phi();
00201     hcalTheta = muonMETInfo.hcalPos.Theta();
00202     hoPhi     = muonMETInfo.hoPos.Phi();
00203     hoTheta   = muonMETInfo.hoPos.Theta();
00204   } else {
00205     
00206     /*
00207       use the analytical solution for the
00208       intersection of a helix with a cylinder
00209       to find the positions of the muon
00210       at the various calo surfaces
00211     */
00212     
00213     //radii of subdetectors in meters
00214     double rEcal = 1.290;
00215     double rHcal = 1.9;
00216     double rHo   = 3.82;
00217     if(abs(mu_eta) > 0.3) rHo = 4.07;
00218     //distance from the center of detector to face of Ecal
00219     double zFaceEcal = 3.209;
00220     if(mu_eta < 0 ) zFaceEcal = -1*zFaceEcal;
00221     //distance from the center of detector to face of Hcal
00222     double zFaceHcal = 3.88;
00223     if(mu_eta < 0 ) zFaceHcal = -1*zFaceHcal;    
00224     
00225     //now we have to get Phi
00226     //bending radius of the muon (units are meters)
00227     double bendr = mu_pt*1000/(300*bfield);
00228 
00229     double tb_ecal = TMath::ACos(1-rEcal*rEcal/(2*bendr*bendr)); //helix time interval parameter
00230     double tb_hcal = TMath::ACos(1-rHcal*rHcal/(2*bendr*bendr)); //helix time interval parameter
00231     double tb_ho   = TMath::ACos(1-rHo*rHo/(2*bendr*bendr));     //helix time interval parameter
00232     double xEcal,yEcal,zEcal;
00233     double xHcal,yHcal,zHcal; 
00234     double xHo, yHo,zHo;
00235     //Ecal
00236     //in the barrel and if not a looper
00237     if(fabs(mu_pz*bendr*tb_ecal/mu_pt+mu_vz) < fabs(zFaceEcal) && rEcal < 2*bendr) { 
00238       xEcal = bendr*(TMath::Sin(tb_ecal+mu_phi)-TMath::Sin(mu_phi));
00239       yEcal = bendr*(-TMath::Cos(tb_ecal+mu_phi)+TMath::Cos(mu_phi));
00240       zEcal = bendr*tb_ecal*mu_pz/mu_pt + mu_vz;
00241     } else { //endcap 
00242       if(mu_pz > 0) {
00243         double te_ecal = (fabs(zFaceEcal) - mu_vz)*mu_pt/(bendr*mu_pz);
00244         xEcal = bendr*(TMath::Sin(te_ecal+mu_phi) - TMath::Sin(mu_phi));
00245         yEcal = bendr*(-TMath::Cos(te_ecal+mu_phi) + TMath::Cos(mu_phi));
00246         zEcal = fabs(zFaceEcal);
00247       } else {
00248         double te_ecal = -(fabs(zFaceEcal) + mu_vz)*mu_pt/(bendr*mu_pz);
00249         xEcal = bendr*(TMath::Sin(te_ecal+mu_phi) - TMath::Sin(mu_phi));
00250         yEcal = bendr*(-TMath::Cos(te_ecal+mu_phi) + TMath::Cos(mu_phi));
00251         zEcal = -fabs(zFaceEcal);
00252       }
00253     }
00254 
00255     //Hcal
00256     if(fabs(mu_pz*bendr*tb_hcal/mu_pt+mu_vz) < fabs(zFaceHcal) && rEcal < 2*bendr) { //in the barrel
00257       xHcal = bendr*(TMath::Sin(tb_hcal+mu_phi)-TMath::Sin(mu_phi));
00258       yHcal = bendr*(-TMath::Cos(tb_hcal+mu_phi)+TMath::Cos(mu_phi));
00259       zHcal = bendr*tb_hcal*mu_pz/mu_pt + mu_vz;
00260     } else { //endcap 
00261       if(mu_pz > 0) {
00262         double te_hcal = (fabs(zFaceHcal) - mu_vz)*mu_pt/(bendr*mu_pz);
00263         xHcal = bendr*(TMath::Sin(te_hcal+mu_phi) - TMath::Sin(mu_phi));
00264         yHcal = bendr*(-TMath::Cos(te_hcal+mu_phi) + TMath::Cos(mu_phi));
00265         zHcal = fabs(zFaceHcal);
00266       } else {
00267         double te_hcal = -(fabs(zFaceHcal) + mu_vz)*mu_pt/(bendr*mu_pz);
00268         xHcal = bendr*(TMath::Sin(te_hcal+mu_phi) - TMath::Sin(mu_phi));
00269         yHcal = bendr*(-TMath::Cos(te_hcal+mu_phi) + TMath::Cos(mu_phi));
00270         zHcal = -fabs(zFaceHcal);
00271       }
00272     }
00273     
00274     //Ho - just the barrel
00275     xHo = bendr*(TMath::Sin(tb_ho+mu_phi)-TMath::Sin(mu_phi));
00276     yHo = bendr*(-TMath::Cos(tb_ho+mu_phi)+TMath::Cos(mu_phi));
00277     zHo = bendr*tb_ho*mu_pz/mu_pt + mu_vz;  
00278     
00279     ecalTheta = TMath::ACos(zEcal/sqrt(pow(xEcal,2) + pow(yEcal,2)+pow(zEcal,2)));
00280     ecalPhi   = atan2(yEcal,xEcal);
00281     hcalTheta = TMath::ACos(zHcal/sqrt(pow(xHcal,2) + pow(yHcal,2)+pow(zHcal,2)));
00282     hcalPhi   = atan2(yHcal,xHcal);
00283     hoTheta   = TMath::ACos(zHo/sqrt(pow(xHo,2) + pow(yHo,2)+pow(zHo,2)));
00284     hoPhi     = atan2(yHo,xHo);
00285 
00286     //2d radius in x-y plane
00287     double r2dEcal = sqrt(pow(xEcal,2)+pow(yEcal,2));
00288     double r2dHcal = sqrt(pow(xHcal,2)+pow(yHcal,2));
00289     double r2dHo   = sqrt(pow(xHo,2)  +pow(yHo,2));
00290     
00291     /*
00292       the above prescription is for right handed helicies only
00293       Positively charged muons trace a left handed helix
00294       so we correct for that 
00295     */
00296     if(muonCharge > 0) {
00297          
00298       //Ecal
00299       double dphi = mu_phi - ecalPhi;
00300       if(fabs(dphi) > TMath::Pi()) 
00301         dphi = 2*TMath::Pi() - fabs(dphi);
00302       ecalPhi = mu_phi - fabs(dphi);
00303       if(fabs(ecalPhi) > TMath::Pi()) {
00304         double temp = 2*TMath::Pi() - fabs(ecalPhi);
00305         ecalPhi = -1*temp*ecalPhi/fabs(ecalPhi);
00306       }
00307       xEcal = r2dEcal*TMath::Cos(ecalPhi);
00308       yEcal = r2dEcal*TMath::Sin(ecalPhi);
00309       
00310       //Hcal
00311       dphi = mu_phi - hcalPhi;
00312       if(fabs(dphi) > TMath::Pi()) 
00313         dphi = 2*TMath::Pi() - fabs(dphi);
00314       hcalPhi = mu_phi - fabs(dphi);
00315       if(fabs(hcalPhi) > TMath::Pi()) {
00316         double temp = 2*TMath::Pi() - fabs(hcalPhi);
00317         hcalPhi = -1*temp*hcalPhi/fabs(hcalPhi);
00318       }
00319       xHcal = r2dHcal*TMath::Cos(hcalPhi);
00320       yHcal = r2dHcal*TMath::Sin(hcalPhi);
00321          
00322       //Ho
00323       dphi = mu_phi - hoPhi;
00324       if(fabs(dphi) > TMath::Pi())
00325         dphi = 2*TMath::Pi() - fabs(dphi);
00326       hoPhi = mu_phi - fabs(dphi);
00327       if(fabs(hoPhi) > TMath::Pi()) {
00328         double temp = 2*TMath::Pi() - fabs(hoPhi);
00329         hoPhi = -1*temp*hoPhi/fabs(hoPhi);
00330       }
00331       xHo = r2dHo*TMath::Cos(hoPhi);
00332       yHo = r2dHo*TMath::Sin(hoPhi);
00333       
00334     }
00335   }
00336   
00337   //for isolated muons
00338   if(!muonMETInfo.useAverage) {
00339     
00340     double mu_Ex =  muonMETInfo.ecalE*sin(ecalTheta)*cos(ecalPhi)
00341       + muonMETInfo.hcalE*sin(hcalTheta)*cos(hcalPhi)
00342       + muonMETInfo.hoE*sin(hoTheta)*cos(hoPhi);
00343     double mu_Ey =  muonMETInfo.ecalE*sin(ecalTheta)*sin(ecalPhi)
00344       + muonMETInfo.hcalE*sin(hcalTheta)*sin(hcalPhi)
00345       + muonMETInfo.hoE*sin(hoTheta)*sin(hoPhi);
00346 
00347     metx += mu_Ex;
00348     mety += mu_Ey;
00349     
00350   } else { //non-isolated muons - derive the correction
00351     
00352     //dE/dx in matter for iron:
00353     //-(11.4 + 0.96*fabs(log(p0*2.8)) + 0.033*p0*(1.0 - pow(p0, -0.33)) )*1e-3
00354     //from https://cmslxr.fnal.gov/lxr/source/TrackPropagation/SteppingHelixPropagator/src/SteppingHelixPropagator.ccyes,
00355     //line ~1100
00356     //normalisation is at 50 GeV
00357     double dEdx_normalization = -(11.4 + 0.96*fabs(log(50*2.8)) + 0.033*50*(1.0 - pow(50, -0.33)) )*1e-3;
00358     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;
00359     
00360     double temp = 0.0;
00361     
00362     if(muonMETInfo.useHO) {
00363       //for the Towers, with HO
00364       if(fabs(mu_eta) < 0.2)
00365         temp = 2.75*(1-0.00003*mu_p);
00366       if(fabs(mu_eta) > 0.2 && fabs(mu_eta) < 1.0)
00367         temp = (2.38+0.0144*fabs(mu_eta))*(1-0.0003*mu_p);
00368       if(fabs(mu_eta) > 1.0 && fabs(mu_eta) < 1.3)
00369         temp = 7.413-5.12*fabs(mu_eta);
00370       if(fabs(mu_eta) > 1.3)
00371         temp = 2.084-0.743*fabs(mu_eta);
00372     } else {
00373       if(fabs(mu_eta) < 1.0)
00374         temp = 2.33*(1-0.0004*mu_p);
00375       if(fabs(mu_eta) > 1.0 && fabs(mu_eta) < 1.3)
00376         temp = (7.413-5.12*fabs(mu_eta))*(1-0.0003*mu_p);
00377       if(fabs(mu_eta) > 1.3)
00378         temp = 2.084-0.743*fabs(mu_eta);
00379     }
00380 
00381     double dep = temp*dEdx_normalization/dEdx_numerator;
00382     if(dep < 0.5) dep = 0;
00383     //use the average phi of the 3 subdetectors
00384     if(fabs(mu_eta) < 1.3) {
00385       metx += dep*cos((ecalPhi+hcalPhi+hoPhi)/3);
00386       mety += dep*sin((ecalPhi+hcalPhi+hoPhi)/3);
00387     } else {
00388       metx += dep*cos( (ecalPhi+hcalPhi)/2);
00389       mety += dep*cos( (ecalPhi+hcalPhi)/2);
00390     }
00391   }
00392 
00393   
00394 }
00395 //----------------------------------------------------------------------------
00396 void MuonMETAlgo::run(const edm::Event& iEvent,
00397                       const edm::EventSetup& iSetup,
00398                       const edm::View<reco::MET>& uncorMET, 
00399                       const edm::View<reco::Muon>& Muons, 
00400                       TrackDetectorAssociator& trackAssociator,
00401                       TrackAssociatorParameters& trackAssociatorParameters,
00402                       METCollection *corMET,
00403                       bool useTrackAssociatorPositions,
00404                       bool useRecHits,
00405                       bool useHO,
00406                       double towerEtThreshold) {
00407                      
00408   
00409   return MuonMETAlgo_run(iEvent, iSetup, uncorMET, Muons, 
00410                          trackAssociator, trackAssociatorParameters,
00411                          corMET, useTrackAssociatorPositions,
00412                          useRecHits, useHO, towerEtThreshold);
00413 }
00414 
00415 
00416 //----------------------------------------------------------------------------
00417 void MuonMETAlgo::run(const edm::Event& iEvent,
00418                       const edm::EventSetup& iSetup,
00419                       const edm::View<reco::CaloMET>& uncorMET, 
00420                       const edm::View<reco::Muon>& Muons,
00421                       TrackDetectorAssociator& trackAssociator,
00422                       TrackAssociatorParameters& trackAssociatorParameters,
00423                       CaloMETCollection *corMET,
00424                       bool useTrackAssociatorPositions,
00425                       bool useRecHits,
00426                       bool useHO,
00427                       double towerEtThreshold) {
00428                       
00429   
00430   return MuonMETAlgo_run(iEvent, iSetup, uncorMET, Muons,
00431                          trackAssociator, trackAssociatorParameters,
00432                          corMET, useTrackAssociatorPositions,
00433                          useRecHits, useHO, towerEtThreshold);
00434 }
00435 
00436 
00437 //----------------------------------------------------------------------------
00438 MuonMETAlgo::MuonMETAlgo() {}
00439 //----------------------------------------------------------------------------
00440   
00441 //----------------------------------------------------------------------------
00442 MuonMETAlgo::~MuonMETAlgo() {}
00443 //----------------------------------------------------------------------------
00444