File indexing completed on 2023-03-17 11:15:48
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
0061
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
0070
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
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
0095
0096
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
0132 return 1. / (1. - pigaga);
0133 }
0134
0135 }