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