/ElectroWeakAnalysis/Utilities/src/FSRWeightProducer.cc

  1#include <memory>
  2
  3#include "FWCore/Framework/interface/Frameworkfwd.h"
  4#include "FWCore/Framework/interface/EDProducer.h"
  5
  6#include "FWCore/Framework/interface/Event.h"
  7#include "FWCore/Framework/interface/MakerMacros.h"
  8
  9#include "FWCore/ParameterSet/interface/ParameterSet.h"
 10
 11#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
 12#include "DataFormats/Candidate/interface/CompositeCandidate.h"
 13#include "CommonTools/CandUtils/interface/AddFourMomenta.h"
 14#include "CommonTools/CandUtils/interface/Booster.h"
 15#include <Math/VectorUtil.h>
 16
 17//
 18// class declaration
 19//
 20class FSRWeightProducer : public edm::EDProducer {
 21   public:
 22      explicit FSRWeightProducer(const edm::ParameterSet&);
 23      ~FSRWeightProducer();
 24
 25   private:
 26      virtual void beginJob() override ;
 27      virtual void produce(edm::Event&, const edm::EventSetup&) override;
 28      virtual void endJob() override ;
 29      double alphaRatio(double) ;
 30
 31      edm::InputTag genTag_;
 32};
 33
 34
 35#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
 36/////////////////////////////////////////////////////////////////////////////////////
 37FSRWeightProducer::FSRWeightProducer(const edm::ParameterSet& pset) {
 38      genTag_ = pset.getUntrackedParameter<edm::InputTag> ("GenTag", edm::InputTag("genParticles"));
 39
 40      produces<double>();
 41} 
 42
 43/////////////////////////////////////////////////////////////////////////////////////
 44FSRWeightProducer::~FSRWeightProducer(){}
 45
 46/////////////////////////////////////////////////////////////////////////////////////
 47void FSRWeightProducer::beginJob() {}
 48
 49/////////////////////////////////////////////////////////////////////////////////////
 50void FSRWeightProducer::endJob(){}
 51
 52/////////////////////////////////////////////////////////////////////////////////////
 53void FSRWeightProducer::produce(edm::Event& iEvent, const edm::EventSetup&) {
 54
 55      if (iEvent.isRealData()) return;
 56
 57      edm::Handle<reco::GenParticleCollection> genParticles;
 58      iEvent.getByLabel(genTag_, genParticles);
 59
 60      std::auto_ptr<double> weight (new double);
 61
 62      // Set a default weight to start with
 63      (*weight) = 1.;
 64
 65      unsigned int gensize = genParticles->size();
 66      for (unsigned int i = 0; i<gensize; ++i) {
 67            const reco::GenParticle& lepton = (*genParticles)[i];
 68            if (lepton.status()!=3) continue;
 69            int leptonId = lepton.pdgId();
 70            if (abs(leptonId)!=11 && abs(leptonId)!=13 && abs(leptonId)!=15) continue;
 71            if (lepton.numberOfMothers()!=1) continue;
 72            const reco::Candidate * boson = lepton.mother();
 73            int bosonId = abs(boson->pdgId());
 74            if (bosonId!=23  && bosonId!=24) continue;
 75            double bosonMass = boson->mass();
 76            double leptonMass = lepton.mass();
 77            double leptonEnergy = lepton.energy();
 78            double cosLeptonTheta = cos(lepton.theta());
 79            double sinLeptonTheta = sin(lepton.theta());
 80            double leptonPhi = lepton.phi();
 81
 82            int trueKey = i;
 83            if (lepton.numberOfDaughters()==0) { 
 84                  continue;
 85            } else if (lepton.numberOfDaughters()==1) { 
 86                  int otherleptonKey = lepton.daughterRef(0).key();
 87                  const reco::GenParticle& otherlepton = (*genParticles)[otherleptonKey];
 88                  if (otherlepton.pdgId()!=leptonId) continue;
 89                  if (otherlepton.numberOfDaughters()<=1) continue;
 90                  trueKey = otherleptonKey;
 91            }
 92
 93            const reco::GenParticle& trueLepton = (*genParticles)[trueKey];
 94            unsigned int nDaughters = trueLepton.numberOfDaughters();
 95
 96            for (unsigned int j = 0; j<nDaughters; ++j) {
 97                  const reco::Candidate * photon = trueLepton.daughter(j);
 98                  if (photon->pdgId()!=22) continue;
 99                  double photonEnergy = photon->energy();
100                  double cosPhotonTheta = cos(photon->theta());
101                  double sinPhotonTheta = sin(photon->theta());
102                  double photonPhi = photon->phi();
103                  double costheta = sinLeptonTheta*sinPhotonTheta*cos(leptonPhi-photonPhi)
104                                  + cosLeptonTheta*cosPhotonTheta;
105                  // Missing O(alpha) terms in soft-collinear approach
106                  // Only for W, from hep-ph/0303260
107                  if (bosonId==24) {
108                        double betaLepton = sqrt(1-pow(leptonMass/leptonEnergy,2));
109                        double delta = - 8*photonEnergy *(1-betaLepton*costheta)
110                          / pow(bosonMass,3) 
111                          / (1-pow(leptonMass/bosonMass,2))
112                          / (4-pow(leptonMass/bosonMass,2))
113                          * leptonEnergy * (pow(leptonMass,2)/bosonMass+2*photonEnergy);
114                        (*weight) *= (1 + delta);
115                  }
116                  // Missing NLO QED orders in QED parton shower approach
117                  // Change coupling scale from 0 to kT to estimate this effect
118                  (*weight) *= alphaRatio(photonEnergy*sqrt(1-pow(costheta,2)));
119            }
120      }
121
122
123      iEvent.put(weight);
124}
125
126double FSRWeightProducer::alphaRatio(double pt) {
127
128      double pigaga = 0.;
129
130      // Leptonic contribution (just one loop, precise at < 0.3% level)
131      const double alphapi = 1/137.036/M_PI;
132      const double mass_e = 0.0005;
133      const double mass_mu = 0.106;
134      const double mass_tau = 1.777;
135      const double mass_Z = 91.2;
136      if (pt>mass_e) pigaga += alphapi * (2*log(pt/mass_e)/3.-5./9.);
137      if (pt>mass_mu) pigaga += alphapi * (2*log(pt/mass_mu)/3.-5./9.);
138      if (pt>mass_tau) pigaga += alphapi * (2*log(pt/mass_tau)/3.-5./9.);
139
140      // Hadronic vaccum contribution
141      // Using simple effective parametrization from Physics Letters B 513 (2001) 46.
142      // Top contribution neglected
143      double A = 0.; 
144      double B = 0.; 
145      double C = 0.; 
146      if (pt<0.7) {
147            A = 0.0; B = 0.0023092; C = 3.9925370;
148      } else if (pt<2.0) {
149            A = 0.0; B = 0.0022333; C = 4.2191779;
150      } else if (pt<4.0) {
151            A = 0.0; B = 0.0024402; C = 3.2496684;
152      } else if (pt<10.0) {
153            A = 0.0; B = 0.0027340; C = 2.0995092;
154      } else if (pt<mass_Z) {
155            A = 0.0010485; B = 0.0029431; C = 1.0;
156      } else if (pt<10000.) {
157            A = 0.0012234; B = 0.0029237; C = 1.0;
158      } else {
159            A = 0.0016894; B = 0.0028984; C = 1.0;
160      }
161      pigaga += A + B*log(1.+C*pt*pt);
162
163      // Done
164      return 1./(1.-pigaga);
165}
166
167DEFINE_FWK_MODULE(FSRWeightProducer);