Reference Manual

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