TopTools/src/MEzCalculator.cc

0001 #include "TopQuarkAnalysis/TopTools/interface/MEzCalculator.h"
0002 #include "TMath.h"
0003
0004 /// constructor
0005 MEzCalculator::MEzCalculator() {
0006   isComplex_ = false;
0007   isMuon_ = true;
0008 }
0009
0010 /// destructor
0011 MEzCalculator::~MEzCalculator() {}
0012
0013 /// member functions
0014 double MEzCalculator::Calculate(int type) {
0015   if (type < 0 || type > 3)
0016     throw cms::Exception("UnimplementedFeature") << "Type " << type << " not supported in MEzCalculator.\n";
0017
0018   double M_W = 80.4;
0019   double M_mu = 0.10566;
0020   double M_e = 0.511e-3;
0021   double M_lepton = M_mu;
0022   if (!isMuon_)
0023     M_lepton = M_e;
0024
0025   double emu = lepton_.energy();
0026   double pxmu = lepton_.px();
0027   double pymu = lepton_.py();
0028   double pzmu = lepton_.pz();
0029   double pxnu = MET_.px();
0030   double pynu = MET_.py();
0031   double pznu = 0.;
0032
0033   // use pznu = - B/2*A +/- sqrt(B*B-4*A*C)/(2*A)
0034
0035   double a = M_W * M_W - M_lepton * M_lepton + 2.0 * (pxmu * pxnu + pymu * pynu);
0036   double A = 4.0 * (emu * emu - pzmu * pzmu);
0037   double B = -4.0 * a * pzmu;
0038   double C = 4.0 * emu * emu * (pxnu * pxnu + pynu * pynu) - a * a;
0039
0040   double tmproot = B * B - 4.0 * A * C;
0041
0042   if (tmproot < 0) {
0043     isComplex_ = true;
0044     pznu = -B / (2 * A);  // take real part of complex roots
0045   } else {
0046     isComplex_ = false;
0047     double tmpsol1 = (-B + TMath::Sqrt(tmproot)) / (2.0 * A);
0048     double tmpsol2 = (-B - TMath::Sqrt(tmproot)) / (2.0 * A);
0049
0050     if (type == 0) {
0051       // two real roots, pick the one closest to pz of muon
0052       if (TMath::Abs(tmpsol2 - pzmu) < TMath::Abs(tmpsol1 - pzmu)) {
0053         pznu = tmpsol2;
0054       } else
0055         pznu = tmpsol1;
0056       // if pznu is > 300 pick the most central root
0057       if (pznu > 300.) {
0058         if (TMath::Abs(tmpsol1) < TMath::Abs(tmpsol2))
0059           pznu = tmpsol1;
0060         else
0061           pznu = tmpsol2;
0062       }
0063     }
0064     if (type == 1) {
0065       // two real roots, pick the one closest to pz of muon
0066       if (TMath::Abs(tmpsol2 - pzmu) < TMath::Abs(tmpsol1 - pzmu)) {
0067         pznu = tmpsol2;
0068       } else
0069         pznu = tmpsol1;
0070     }
0071     if (type == 2) {
0072       // pick the most central root.
0073       if (TMath::Abs(tmpsol1) < TMath::Abs(tmpsol2))
0074         pznu = tmpsol1;
0075       else
0076         pznu = tmpsol2;
0077     }
0078     if (type == 3) {
0079       // pick the largest value of the cosine
0080       TVector3 p3w, p3mu;
0081       p3w.SetXYZ(pxmu + pxnu, pymu + pynu, pzmu + tmpsol1);
0082       p3mu.SetXYZ(pxmu, pymu, pzmu);
0083
0084       double sinthcm1 = 2. * (p3mu.Perp(p3w)) / M_W;
0085       p3w.SetXYZ(pxmu + pxnu, pymu + pynu, pzmu + tmpsol2);
0086       double sinthcm2 = 2. * (p3mu.Perp(p3w)) / M_W;
0087
0088       double costhcm1 = TMath::Sqrt(1. - sinthcm1 * sinthcm1);
0089       double costhcm2 = TMath::Sqrt(1. - sinthcm2 * sinthcm2);
0090
0091       if (costhcm1 > costhcm2)
0092         pznu = tmpsol1;
0093       else
0094         pznu = tmpsol2;
0095     }
0096   }
0097
0098   //Particle neutrino;
0099   //neutrino.setP4( LorentzVector(pxnu, pynu, pznu, TMath::Sqrt(pxnu*pxnu + pynu*pynu + pznu*pznu ))) ;
0100
0101   return pznu;
0102 }