Heppy/src/FSRWeightAlgo.cc

0001 #include "PhysicsTools/Heppy/interface/FSRWeightAlgo.h"
0002
0003 #include "DataFormats/Candidate/interface/CompositeCandidate.h"
0004 #include "CommonTools/CandUtils/interface/AddFourMomenta.h"
0005 #include "CommonTools/CandUtils/interface/Booster.h"
0006 #include <Math/VectorUtil.h>
0007
0008 namespace heppy {
0009
0010   double FSRWeightAlgo::weight() const {
0011     double weight = 1.;
0012
0013     unsigned int gensize = genParticles_.size();
0014     for (unsigned int i = 0; i < gensize; ++i) {
0015       const reco::GenParticle& lepton = genParticles_[i];
0016       if (lepton.status() != 3)
0017         continue;
0018       int leptonId = lepton.pdgId();
0019       if (abs(leptonId) != 11 && abs(leptonId) != 13 && abs(leptonId) != 15)
0020         continue;
0021       if (lepton.numberOfMothers() != 1)
0022         continue;
0023       const reco::Candidate* boson = lepton.mother();
0024       int bosonId = abs(boson->pdgId());
0025       if (bosonId != 23 && bosonId != 24)
0026         continue;
0027       double bosonMass = boson->mass();
0028       double leptonMass = lepton.mass();
0029       double leptonEnergy = lepton.energy();
0030       double cosLeptonTheta = cos(lepton.theta());
0031       double sinLeptonTheta = sin(lepton.theta());
0032       double leptonPhi = lepton.phi();
0033
0034       int trueKey = i;
0035       if (lepton.numberOfDaughters() == 0) {
0036         continue;
0037       } else if (lepton.numberOfDaughters() == 1) {
0038         int otherleptonKey = lepton.daughterRef(0).key();
0039         const reco::GenParticle& otherlepton = genParticles_[otherleptonKey];
0040         if (otherlepton.pdgId() != leptonId)
0041           continue;
0042         if (otherlepton.numberOfDaughters() <= 1)
0043           continue;
0044         trueKey = otherleptonKey;
0045       }
0046
0047       const reco::GenParticle& trueLepton = genParticles_[trueKey];
0048       unsigned int nDaughters = trueLepton.numberOfDaughters();
0049
0050       for (unsigned int j = 0; j < nDaughters; ++j) {
0051         const reco::Candidate* photon = trueLepton.daughter(j);
0052         if (photon->pdgId() != 22)
0053           continue;
0054         double photonEnergy = photon->energy();
0055         double cosPhotonTheta = cos(photon->theta());
0056         double sinPhotonTheta = sin(photon->theta());
0057         double photonPhi = photon->phi();
0058         double costheta =
0059             sinLeptonTheta * sinPhotonTheta * cos(leptonPhi - photonPhi) + cosLeptonTheta * cosPhotonTheta;
0060         // Missing O(alpha) terms in soft-collinear approach
0061         // Only for W, from hep-ph/0303260
0062         if (bosonId == 24) {
0063           double betaLepton = sqrt(1 - pow(leptonMass / leptonEnergy, 2));
0064           double delta = -8 * photonEnergy * (1 - betaLepton * costheta) / pow(bosonMass, 3) /
0065                          (1 - pow(leptonMass / bosonMass, 2)) / (4 - pow(leptonMass / bosonMass, 2)) * leptonEnergy *
0066                          (pow(leptonMass, 2) / bosonMass + 2 * photonEnergy);
0067           weight *= (1 + delta);
0068         }
0069         // Missing NLO QED orders in QED parton shower approach
0070         // Change coupling scale from 0 to kT to estimate this effect
0071         weight *= alphaRatio(photonEnergy * sqrt(1 - pow(costheta, 2)));
0072       }
0073     }
0074
0075     return weight;
0076   }
0077
0078   double FSRWeightAlgo::alphaRatio(double pt) const {
0079     double pigaga = 0.;
0080
0081     // Leptonic contribution (just one loop, precise at < 0.3% level)
0082     const double alphapi = 1 / 137.036 / M_PI;
0083     const double mass_e = 0.0005;
0084     const double mass_mu = 0.106;
0085     const double mass_tau = 1.777;
0086     const double mass_Z = 91.2;
0087     if (pt > mass_e)
0088       pigaga += alphapi * (2 * log(pt / mass_e) / 3. - 5. / 9.);
0089     if (pt > mass_mu)
0090       pigaga += alphapi * (2 * log(pt / mass_mu) / 3. - 5. / 9.);
0091     if (pt > mass_tau)
0092       pigaga += alphapi * (2 * log(pt / mass_tau) / 3. - 5. / 9.);
0093
0094     // Hadronic vaccum contribution
0095     // Using simple effective parametrization from Physics Letters B 513 (2001) 46.
0096     // Top contribution neglected
0097     double A = 0.;
0098     double B = 0.;
0099     double C = 0.;
0100     if (pt < 0.7) {
0101       A = 0.0;
0102       B = 0.0023092;
0103       C = 3.9925370;
0104     } else if (pt < 2.0) {
0105       A = 0.0;
0106       B = 0.0022333;
0107       C = 4.2191779;
0108     } else if (pt < 4.0) {
0109       A = 0.0;
0110       B = 0.0024402;
0111       C = 3.2496684;
0112     } else if (pt < 10.0) {
0113       A = 0.0;
0114       B = 0.0027340;
0115       C = 2.0995092;
0116     } else if (pt < mass_Z) {
0117       A = 0.0010485;
0118       B = 0.0029431;
0119       C = 1.0;
0120     } else if (pt < 10000.) {
0121       A = 0.0012234;
0122       B = 0.0029237;
0123       C = 1.0;
0124     } else {
0125       A = 0.0016894;
0126       B = 0.0028984;
0127       C = 1.0;
0128     }
0129     pigaga += A + B * log(1. + C * pt * pt);
0130
0131     // Done
0132     return 1. / (1. - pigaga);
0133   }
0134
0135 }  // namespace heppy