CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Public Member Functions | Static Public Member Functions
MuonMETAlgo Class Reference

#include <MuonMETAlgo.h>

Public Member Functions

void GetMuDepDeltas (const reco::Muon *inputMuon, TrackDetMatchInfo &info, bool useTrackAssociatorPositions, bool useRecHits, bool useHO, double towerEtThreshold, double &deltax, double &deltay, double Bfield)
 
reco::CaloMET makeMET (const reco::CaloMET &fMet, double fSumEt, const std::vector< CorrMETData > &fCorrections, const reco::MET::LorentzVector &)
 
reco::MET makeMET (const reco::MET &, double fSumEt, const std::vector< CorrMETData > &fCorrections, const reco::MET::LorentzVector &fP4)
 
 MuonMETAlgo ()
 
template<class T >
void MuonMETAlgo_run (const edm::View< reco::Muon > &inputMuons, const edm::ValueMap< reco::MuonMETCorrectionData > &vm_muCorrData, const edm::View< T > &v_uncorMET, std::vector< T > *v_corMET)
 
virtual void run (const edm::View< reco::Muon > &inputMuons, const edm::ValueMap< reco::MuonMETCorrectionData > &vm_muCorrData, const edm::View< reco::MET > &uncorMET, reco::METCollection *corMET)
 
virtual void run (const edm::View< reco::Muon > &inputMuons, const edm::ValueMap< reco::MuonMETCorrectionData > &vm_muCorrData, const edm::View< reco::CaloMET > &uncorMET, reco::CaloMETCollection *corMET)
 
virtual ~MuonMETAlgo ()
 

Static Public Member Functions

static void correctMETforMuon (double &deltax, double &deltay, double bfield, int muonCharge, const math::XYZTLorentzVector muonP4, const math::XYZPoint muonVertex, MuonMETInfo &)
 

Detailed Description

Correct MET for muons in the events.

Version
1st Version August 30, 2007

Definition at line 32 of file MuonMETAlgo.h.

Constructor & Destructor Documentation

MuonMETAlgo::MuonMETAlgo ( )

Definition at line 433 of file MuonMETAlgo.cc.

433 {}
MuonMETAlgo::~MuonMETAlgo ( )
virtual

Definition at line 437 of file MuonMETAlgo.cc.

437 {}

Member Function Documentation

void MuonMETAlgo::correctMETforMuon ( double &  deltax,
double &  deltay,
double  bfield,
int  muonCharge,
const math::XYZTLorentzVector  muonP4,
const math::XYZPoint  muonVertex,
MuonMETInfo muonMETInfo 
)
static

Definition at line 194 of file MuonMETAlgo.cc.

References abs, funct::cos(), ExpressReco_HICollisions_FallBack::e, MuonMETInfo::ecalE, MuonMETInfo::ecalPos, MuonMETInfo::hcalE, MuonMETInfo::hcalPos, MuonMETInfo::hoE, MuonMETInfo::hoPos, funct::log(), Pi, funct::pow(), funct::sin(), mathSSE::sqrt(), cond::rpcobtemp::temp, MuonMETInfo::useAverage, MuonMETInfo::useHO, and MuonMETInfo::useTkAssociatorPositions.

196  {
197 
198  double mu_p = muonP4.P();
199  double mu_pt = muonP4.Pt();
200  double mu_phi = muonP4.Phi();
201  double mu_eta = muonP4.Eta();
202  double mu_vz = muonVertex.z()/100.;
203  double mu_pz = muonP4.Pz();
204 
205  double ecalPhi, ecalTheta;
206  double hcalPhi, hcalTheta;
207  double hoPhi, hoTheta;
208 
209 
210  //should always be false for FWLite
211  //unless you want to supply co-ordinates at
212  //the calorimeter sub-detectors yourself
213  if(muonMETInfo.useTkAssociatorPositions) {
214  ecalPhi = muonMETInfo.ecalPos.Phi();
215  ecalTheta = muonMETInfo.ecalPos.Theta();
216  hcalPhi = muonMETInfo.hcalPos.Phi();
217  hcalTheta = muonMETInfo.hcalPos.Theta();
218  hoPhi = muonMETInfo.hoPos.Phi();
219  hoTheta = muonMETInfo.hoPos.Theta();
220  } else {
221 
222  /*
223  use the analytical solution for the
224  intersection of a helix with a cylinder
225  to find the positions of the muon
226  at the various calo surfaces
227  */
228 
229  //radii of subdetectors in meters
230  double rEcal = 1.290;
231  double rHcal = 1.9;
232  double rHo = 3.82;
233  if(abs(mu_eta) > 0.3) rHo = 4.07;
234  //distance from the center of detector to face of Ecal
235  double zFaceEcal = 3.209;
236  if(mu_eta < 0 ) zFaceEcal = -1*zFaceEcal;
237  //distance from the center of detector to face of Hcal
238  double zFaceHcal = 3.88;
239  if(mu_eta < 0 ) zFaceHcal = -1*zFaceHcal;
240 
241  //now we have to get Phi
242  //bending radius of the muon (units are meters)
243  double bendr = mu_pt*1000/(300*bfield);
244 
245  double tb_ecal = TMath::ACos(1-rEcal*rEcal/(2*bendr*bendr)); //helix time interval parameter
246  double tb_hcal = TMath::ACos(1-rHcal*rHcal/(2*bendr*bendr)); //helix time interval parameter
247  double tb_ho = TMath::ACos(1-rHo*rHo/(2*bendr*bendr)); //helix time interval parameter
248  double xEcal,yEcal,zEcal;
249  double xHcal,yHcal,zHcal;
250  double xHo, yHo,zHo;
251  //Ecal
252  //in the barrel and if not a looper
253  if(fabs(mu_pz*bendr*tb_ecal/mu_pt+mu_vz) < fabs(zFaceEcal) && rEcal < 2*bendr) {
254  xEcal = bendr*(TMath::Sin(tb_ecal+mu_phi)-TMath::Sin(mu_phi));
255  yEcal = bendr*(-TMath::Cos(tb_ecal+mu_phi)+TMath::Cos(mu_phi));
256  zEcal = bendr*tb_ecal*mu_pz/mu_pt + mu_vz;
257  } else { //endcap
258  if(mu_pz > 0) {
259  double te_ecal = (fabs(zFaceEcal) - mu_vz)*mu_pt/(bendr*mu_pz);
260  xEcal = bendr*(TMath::Sin(te_ecal+mu_phi) - TMath::Sin(mu_phi));
261  yEcal = bendr*(-TMath::Cos(te_ecal+mu_phi) + TMath::Cos(mu_phi));
262  zEcal = fabs(zFaceEcal);
263  } else {
264  double te_ecal = -(fabs(zFaceEcal) + mu_vz)*mu_pt/(bendr*mu_pz);
265  xEcal = bendr*(TMath::Sin(te_ecal+mu_phi) - TMath::Sin(mu_phi));
266  yEcal = bendr*(-TMath::Cos(te_ecal+mu_phi) + TMath::Cos(mu_phi));
267  zEcal = -fabs(zFaceEcal);
268  }
269  }
270 
271  //Hcal
272  if(fabs(mu_pz*bendr*tb_hcal/mu_pt+mu_vz) < fabs(zFaceHcal) && rEcal < 2*bendr) { //in the barrel
273  xHcal = bendr*(TMath::Sin(tb_hcal+mu_phi)-TMath::Sin(mu_phi));
274  yHcal = bendr*(-TMath::Cos(tb_hcal+mu_phi)+TMath::Cos(mu_phi));
275  zHcal = bendr*tb_hcal*mu_pz/mu_pt + mu_vz;
276  } else { //endcap
277  if(mu_pz > 0) {
278  double te_hcal = (fabs(zFaceHcal) - mu_vz)*mu_pt/(bendr*mu_pz);
279  xHcal = bendr*(TMath::Sin(te_hcal+mu_phi) - TMath::Sin(mu_phi));
280  yHcal = bendr*(-TMath::Cos(te_hcal+mu_phi) + TMath::Cos(mu_phi));
281  zHcal = fabs(zFaceHcal);
282  } else {
283  double te_hcal = -(fabs(zFaceHcal) + mu_vz)*mu_pt/(bendr*mu_pz);
284  xHcal = bendr*(TMath::Sin(te_hcal+mu_phi) - TMath::Sin(mu_phi));
285  yHcal = bendr*(-TMath::Cos(te_hcal+mu_phi) + TMath::Cos(mu_phi));
286  zHcal = -fabs(zFaceHcal);
287  }
288  }
289 
290  //Ho - just the barrel
291  xHo = bendr*(TMath::Sin(tb_ho+mu_phi)-TMath::Sin(mu_phi));
292  yHo = bendr*(-TMath::Cos(tb_ho+mu_phi)+TMath::Cos(mu_phi));
293  zHo = bendr*tb_ho*mu_pz/mu_pt + mu_vz;
294 
295  ecalTheta = TMath::ACos(zEcal/sqrt(pow(xEcal,2) + pow(yEcal,2)+pow(zEcal,2)));
296  ecalPhi = atan2(yEcal,xEcal);
297  hcalTheta = TMath::ACos(zHcal/sqrt(pow(xHcal,2) + pow(yHcal,2)+pow(zHcal,2)));
298  hcalPhi = atan2(yHcal,xHcal);
299  hoTheta = TMath::ACos(zHo/sqrt(pow(xHo,2) + pow(yHo,2)+pow(zHo,2)));
300  hoPhi = atan2(yHo,xHo);
301 
302  //2d radius in x-y plane
303  double r2dEcal = sqrt(pow(xEcal,2)+pow(yEcal,2));
304  double r2dHcal = sqrt(pow(xHcal,2)+pow(yHcal,2));
305  double r2dHo = sqrt(pow(xHo,2) +pow(yHo,2));
306 
307  /*
308  the above prescription is for right handed helicies only
309  Positively charged muons trace a left handed helix
310  so we correct for that
311  */
312  if(muonCharge > 0) {
313 
314  //Ecal
315  double dphi = mu_phi - ecalPhi;
316  if(fabs(dphi) > TMath::Pi())
317  dphi = 2*TMath::Pi() - fabs(dphi);
318  ecalPhi = mu_phi - fabs(dphi);
319  if(fabs(ecalPhi) > TMath::Pi()) {
320  double temp = 2*TMath::Pi() - fabs(ecalPhi);
321  ecalPhi = -1*temp*ecalPhi/fabs(ecalPhi);
322  }
323  xEcal = r2dEcal*TMath::Cos(ecalPhi);
324  yEcal = r2dEcal*TMath::Sin(ecalPhi);
325 
326  //Hcal
327  dphi = mu_phi - hcalPhi;
328  if(fabs(dphi) > TMath::Pi())
329  dphi = 2*TMath::Pi() - fabs(dphi);
330  hcalPhi = mu_phi - fabs(dphi);
331  if(fabs(hcalPhi) > TMath::Pi()) {
332  double temp = 2*TMath::Pi() - fabs(hcalPhi);
333  hcalPhi = -1*temp*hcalPhi/fabs(hcalPhi);
334  }
335  xHcal = r2dHcal*TMath::Cos(hcalPhi);
336  yHcal = r2dHcal*TMath::Sin(hcalPhi);
337 
338 
339  //Ho
340  dphi = mu_phi - hoPhi;
341  if(fabs(dphi) > TMath::Pi())
342  dphi = 2*TMath::Pi() - fabs(dphi);
343  hoPhi = mu_phi - fabs(dphi);
344  if(fabs(hoPhi) > TMath::Pi()) {
345  double temp = 2*TMath::Pi() - fabs(hoPhi);
346  hoPhi = -1*temp*hoPhi/fabs(hoPhi);
347  }
348  xHo = r2dHo*TMath::Cos(hoPhi);
349  yHo = r2dHo*TMath::Sin(hoPhi);
350 
351  }
352  }
353 
354  //for isolated muons
355  if(!muonMETInfo.useAverage) {
356 
357  double mu_Ex = muonMETInfo.ecalE*sin(ecalTheta)*cos(ecalPhi)
358  + muonMETInfo.hcalE*sin(hcalTheta)*cos(hcalPhi)
359  + muonMETInfo.hoE*sin(hoTheta)*cos(hoPhi);
360  double mu_Ey = muonMETInfo.ecalE*sin(ecalTheta)*sin(ecalPhi)
361  + muonMETInfo.hcalE*sin(hcalTheta)*sin(hcalPhi)
362  + muonMETInfo.hoE*sin(hoTheta)*sin(hoPhi);
363 
364  deltax += mu_Ex;
365  deltay += mu_Ey;
366 
367  } else { //non-isolated muons - derive the correction
368 
369  //dE/dx in matter for iron:
370  //-(11.4 + 0.96*fabs(log(p0*2.8)) + 0.033*p0*(1.0 - pow(p0, -0.33)) )*1e-3
371  //from http://cmslxr.fnal.gov/lxr/source/TrackPropagation/SteppingHelixPropagator/src/SteppingHelixPropagator.ccyes,
372  //line ~1100
373  //normalisation is at 50 GeV
374  double dEdx_normalization = -(11.4 + 0.96*fabs(log(50*2.8)) + 0.033*50*(1.0 - pow(50, -0.33)) )*1e-3;
375  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;
376 
377  double temp = 0.0;
378 
379  if(muonMETInfo.useHO) {
380  //for the Towers, with HO
381  if(fabs(mu_eta) < 0.2)
382  temp = 2.75*(1-0.00003*mu_p);
383  if(fabs(mu_eta) > 0.2 && fabs(mu_eta) < 1.0)
384  temp = (2.38+0.0144*fabs(mu_eta))*(1-0.0003*mu_p);
385  if(fabs(mu_eta) > 1.0 && fabs(mu_eta) < 1.3)
386  temp = 7.413-5.12*fabs(mu_eta);
387  if(fabs(mu_eta) > 1.3)
388  temp = 2.084-0.743*fabs(mu_eta);
389  } else {
390  if(fabs(mu_eta) < 1.0)
391  temp = 2.33*(1-0.0004*mu_p);
392  if(fabs(mu_eta) > 1.0 && fabs(mu_eta) < 1.3)
393  temp = (7.413-5.12*fabs(mu_eta))*(1-0.0003*mu_p);
394  if(fabs(mu_eta) > 1.3)
395  temp = 2.084-0.743*fabs(mu_eta);
396  }
397 
398  double dep = temp*dEdx_normalization/dEdx_numerator;
399  if(dep < 0.5) dep = 0;
400  //use the average phi of the 3 subdetectors
401  if(fabs(mu_eta) < 1.3) {
402  deltax += dep*cos((ecalPhi+hcalPhi+hoPhi)/3);
403  deltay += dep*sin((ecalPhi+hcalPhi+hoPhi)/3);
404  } else {
405  deltax += dep*cos( (ecalPhi+hcalPhi)/2);
406  deltay += dep*cos( (ecalPhi+hcalPhi)/2);
407  }
408  }
409 
410 
411 }
const double Pi
math::XYZPoint hcalPos
Definition: MuonMETInfo.h:12
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
#define abs(x)
Definition: mlp_lapack.h:159
float hcalE
Definition: MuonMETInfo.h:9
T sqrt(T t)
Definition: SSEVec.h:28
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
math::XYZPoint hoPos
Definition: MuonMETInfo.h:13
bool useTkAssociatorPositions
Definition: MuonMETInfo.h:16
float ecalE
Definition: MuonMETInfo.h:8
Log< T >::type log(const T &t)
Definition: Log.h:22
bool useAverage
Definition: MuonMETInfo.h:14
math::XYZPoint ecalPos
Definition: MuonMETInfo.h:11
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
void MuonMETAlgo::GetMuDepDeltas ( const reco::Muon inputMuon,
TrackDetMatchInfo info,
bool  useTrackAssociatorPositions,
bool  useRecHits,
bool  useHO,
double  towerEtThreshold,
double &  deltax,
double &  deltay,
double  Bfield 
)

Definition at line 102 of file MuonMETAlgo.cc.

References reco::Muon::calEnergy(), reco::LeafCandidate::charge(), TrackDetMatchInfo::crossedTowers, MuonMETInfo::ecalE, MuonMETInfo::ecalPos, reco::MuonIsolation::emEt, reco::MuonEnergy::emS9, reco::Muon::globalTrack(), reco::MuonIsolation::hadEt, reco::MuonEnergy::hadS9, MuonMETInfo::hcalE, MuonMETInfo::hcalPos, MuonMETInfo::hoE, MuonMETInfo::hoPos, reco::MuonEnergy::hoS9, reco::Muon::innerTrack(), reco::Muon::isGlobalMuon(), reco::Muon::isIsolationValid(), reco::Muon::isolationR03(), reco::Muon::isTrackerMuon(), reco::Muon::outerTrack(), reco::MuonIsolation::sumPt, TrackDetMatchInfo::trkGlobPosAtEcal, TrackDetMatchInfo::trkGlobPosAtHcal, TrackDetMatchInfo::trkGlobPosAtHO, MuonMETInfo::useAverage, MuonMETInfo::useHO, ExpressReco_HICollisions_FallBack::useHO, MuonMETInfo::useTkAssociatorPositions, and reco::LeafCandidate::vertex().

Referenced by cms::MuonMETValueMapProducer::produce().

109  {
110 
111 
112  bool useAverage = false;
113  //decide whether or not we want to correct on average based
114  //on isolation information from the muon
115  double sumPt = inputMuon->isIsolationValid()? inputMuon->isolationR03().sumPt : 0.0;
116  double sumEtEcal = inputMuon->isIsolationValid() ? inputMuon->isolationR03().emEt : 0.0;
117  double sumEtHcal = inputMuon->isIsolationValid() ? inputMuon->isolationR03().hadEt : 0.0;
118 
119  if(sumPt > 3 || sumEtEcal + sumEtHcal > 5) useAverage = true;
120 
121  //get the energy using TrackAssociator if
122  //the muon turns out to be isolated
123  MuonMETInfo muMETInfo;
124  muMETInfo.useAverage = useAverage;
125  muMETInfo.useTkAssociatorPositions = useTrackAssociatorPositions;
126  muMETInfo.useHO = useHO;
127 
128 
129  TrackRef mu_track;
130  if(inputMuon->isGlobalMuon()) {
131  mu_track = inputMuon->globalTrack();
132  } else if(inputMuon->isTrackerMuon()) {
133  mu_track = inputMuon->innerTrack();
134  } else
135  mu_track = inputMuon->outerTrack();
136 
137  if(useTrackAssociatorPositions) {
138  muMETInfo.ecalPos = info.trkGlobPosAtEcal;
139  muMETInfo.hcalPos = info.trkGlobPosAtHcal;
140  muMETInfo.hoPos = info.trkGlobPosAtHO;
141  }
142 
143  if(!useAverage) {
144 
145  if(useRecHits) {
146  muMETInfo.ecalE = inputMuon->calEnergy().emS9;
147  muMETInfo.hcalE = inputMuon->calEnergy().hadS9;
148  if(useHO) //muMETInfo.hoE is 0 by default
149  muMETInfo.hoE = inputMuon->calEnergy().hoS9;
150  } else {// use Towers (this is the default)
151  //only include towers whose Et > 0.5 since
152  //by default the MET only includes towers with Et > 0.5
153  std::vector<const CaloTower*> towers = info.crossedTowers;
154  for(vector<const CaloTower*>::const_iterator it = towers.begin();
155  it != towers.end(); it++) {
156  if((*it)->et() < towerEtThreshold) continue;
157  muMETInfo.ecalE += (*it)->emEnergy();
158  muMETInfo.hcalE += (*it)->hadEnergy();
159  if(useHO)
160  muMETInfo.hoE +=(*it)->outerEnergy();
161  }
162  }//use Towers
163  }
164 
165 
166 
167  //This needs to be fixed!!!!!
168  //The tracker has better resolution for pt < 200 GeV
170  if(inputMuon->isGlobalMuon()) {
171  if(inputMuon->globalTrack()->pt() < 200) {
172  mup4 = LorentzVector(inputMuon->innerTrack()->px(), inputMuon->innerTrack()->py(),
173  inputMuon->innerTrack()->pz(), inputMuon->innerTrack()->p());
174  } else {
175  mup4 = LorentzVector(inputMuon->globalTrack()->px(), inputMuon->globalTrack()->py(),
176  inputMuon->globalTrack()->pz(), inputMuon->globalTrack()->p());
177  }
178  } else if(inputMuon->isTrackerMuon()) {
179  mup4 = LorentzVector(inputMuon->innerTrack()->px(), inputMuon->innerTrack()->py(),
180  inputMuon->innerTrack()->pz(), inputMuon->innerTrack()->p());
181  } else
182  mup4 = LorentzVector(inputMuon->outerTrack()->px(), inputMuon->outerTrack()->py(),
183  inputMuon->outerTrack()->pz(), inputMuon->outerTrack()->p());
184 
185 
186  //call function that does the work
187  correctMETforMuon(deltax, deltay, Bfield, inputMuon->charge(),
188  mup4, inputMuon->vertex(),
189  muMETInfo);
190 }
float hadEt
hcal sum-Et
Definition: MuonIsolation.h:9
math::XYZPoint trkGlobPosAtHO
math::XYZPoint hcalPos
Definition: MuonMETInfo.h:12
float sumPt
sum-pt of tracks
Definition: MuonIsolation.h:7
std::vector< const CaloTower * > crossedTowers
virtual TrackRef innerTrack() const
Definition: Muon.h:38
virtual const Point & vertex() const
vertex position
bool isTrackerMuon() const
Definition: Muon.h:149
bool isGlobalMuon() const
Definition: Muon.h:148
float emS9
energy deposited in 3x3 ECAL crystal shape around central crystal
Definition: MuonEnergy.h:18
float hcalE
Definition: MuonMETInfo.h:9
math::XYZPoint trkGlobPosAtHcal
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:30
float emEt
ecal sum-Et
Definition: MuonIsolation.h:8
virtual int charge() const
electric charge
float hoS9
energy deposited in 3x3 HO tower shape around central tower
Definition: MuonEnergy.h:34
math::XYZPoint hoPos
Definition: MuonMETInfo.h:13
bool useTkAssociatorPositions
Definition: MuonMETInfo.h:16
float ecalE
Definition: MuonMETInfo.h:8
virtual TrackRef outerTrack() const
reference to Track reconstructed in the muon detector only
Definition: Muon.h:41
MuonEnergy calEnergy() const
get energy deposition information
Definition: Muon.h:62
bool useAverage
Definition: MuonMETInfo.h:14
math::XYZPoint ecalPos
Definition: MuonMETInfo.h:11
math::XYZPoint trkGlobPosAtEcal
Track position at different parts of the calorimeter.
bool isIsolationValid() const
Definition: Muon.h:112
static void correctMETforMuon(double &deltax, double &deltay, double bfield, int muonCharge, const math::XYZTLorentzVector muonP4, const math::XYZPoint muonVertex, MuonMETInfo &)
Definition: MuonMETAlgo.cc:194
float hadS9
energy deposited in 3x3 HCAL tower shape around central tower
Definition: MuonEnergy.h:28
const MuonIsolation & isolationR03() const
Definition: Muon.h:109
math::PtEtaPhiELorentzVectorF LorentzVector
virtual TrackRef globalTrack() const
reference to Track reconstructed in both tracked and muon detector
Definition: Muon.h:44
reco::CaloMET MuonMETAlgo::makeMET ( const reco::CaloMET fMet,
double  fSumEt,
const std::vector< CorrMETData > &  fCorrections,
const reco::MET::LorentzVector  
)
reco::MET MuonMETAlgo::makeMET ( const reco::MET ,
double  fSumEt,
const std::vector< CorrMETData > &  fCorrections,
const reco::MET::LorentzVector fP4 
)
template<class T >
void MuonMETAlgo::MuonMETAlgo_run ( const edm::View< reco::Muon > &  inputMuons,
const edm::ValueMap< reco::MuonMETCorrectionData > &  vm_muCorrData,
const edm::View< T > &  v_uncorMET,
std::vector< T > *  v_corMET 
)

Definition at line 48 of file MuonMETAlgo.cc.

References reco::MuonMETCorrectionData::corrX(), reco::MuonMETCorrectionData::corrY(), delta, CorrMETData::mex, CorrMETData::mey, reco::LeafCandidate::p4(), edm::View< T >::refAt(), query::result, edm::View< T >::size(), mathSSE::sqrt(), CorrMETData::sumet, and reco::MuonMETCorrectionData::type().

51  {
52  T uncorMETObj = v_uncorMET.front();
53 
54  double uncorMETX = uncorMETObj.px();
55  double uncorMETY = uncorMETObj.py();
56 
57  double corMETX = uncorMETX;
58  double corMETY = uncorMETY;
59 
61  double sumMuPx = 0.;
62  double sumMuPy = 0.;
63  double sumMuDepEx = 0.;
64  double sumMuDepEy = 0.;
65 
66  unsigned int nMuons = inputMuons.size();
67  for(unsigned int iMu = 0; iMu<nMuons; iMu++) {
68  const reco::Muon *mu = &inputMuons[iMu]; //new
69  reco::MuonMETCorrectionData muCorrData = (vm_muCorrData)[inputMuons.refAt(iMu)];
70  int flag = muCorrData.type();
71  float deltax = muCorrData.corrX();
72  float deltay = muCorrData.corrY();
73 
74 
75  LorentzVector mup4;
76  if (flag == 0) //this muon is not used to correct the MET
77  continue;
78 
79  //if we're here, then the muon should be used to correct the MET using the default fit
80  mup4 = mu->p4();
81 
82  sumMuPx += mup4.px();
83  sumMuPy += mup4.py();
84  sumMuDepEx += deltax;
85  sumMuDepEy += deltay;
86  corMETX = corMETX - mup4.px() + deltax;
87  corMETY = corMETY - mup4.py() + deltay;
88 
89  }
90  delta.mex = sumMuDepEx - sumMuPx;
91  delta.mey = sumMuDepEy - sumMuPy;
92  delta.sumet = sqrt(sumMuPx*sumMuPx + sumMuPy*sumMuPy) - sqrt(sumMuDepEx*sumMuDepEx + sumMuDepEy*sumMuDepEy);
93  MET::LorentzVector correctedMET4vector( corMETX, corMETY, 0., sqrt(corMETX*corMETX + corMETY*corMETY));
94  std::vector<CorrMETData> corrections = uncorMETObj.mEtCorr();
95  corrections.push_back(delta);
96 
97  T result = makeMET(uncorMETObj, uncorMETObj.sumEt()+delta.sumet, corrections, correctedMET4vector);
98  v_corMET->push_back(result);
99 
100 }
dbl * delta
Definition: mlp_gen.cc:36
long int flag
Definition: mlp_lapack.h:47
RefToBase< value_type > refAt(size_type i) const
math::XYZTLorentzVector LorentzVector
Definition: HLTMuonBPAG.h:55
T sqrt(T t)
Definition: SSEVec.h:28
tuple result
Definition: query.py:137
double sumet
Definition: CorrMETData.h:26
size_type size() const
Structure containing data common to all types of MET.
Definition: CorrMETData.h:20
double mey
Definition: CorrMETData.h:25
double mex
Definition: CorrMETData.h:24
virtual const LorentzVector & p4() const
four-momentum Lorentz vector
reco::CaloMET makeMET(const reco::CaloMET &fMet, double fSumEt, const std::vector< CorrMETData > &fCorrections, const reco::MET::LorentzVector &)
math::PtEtaPhiELorentzVectorF LorentzVector
virtual void MuonMETAlgo::run ( const edm::View< reco::Muon > &  inputMuons,
const edm::ValueMap< reco::MuonMETCorrectionData > &  vm_muCorrData,
const edm::View< reco::MET > &  uncorMET,
reco::METCollection corMET 
)
virtual

Referenced by cms::MuonMET::produce().

virtual void MuonMETAlgo::run ( const edm::View< reco::Muon > &  inputMuons,
const edm::ValueMap< reco::MuonMETCorrectionData > &  vm_muCorrData,
const edm::View< reco::CaloMET > &  uncorMET,
reco::CaloMETCollection corMET 
)
virtual