CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
WeakEffectsWeightProducer.cc
Go to the documentation of this file.
4 
6 public:
9  virtual void beginJob() ;
10  virtual void produce(edm::Event &, const edm::EventSetup&);
11  virtual void endJob() ;
12 private:
14  double rhoParameter_;
15 
16  double alphaQED(double q2);
17  double sigma0_qqbarll(unsigned int quark_type, double Q, double rho);
18 
19 };
20 
28 
30 
33  genParticlesTag_(pset.getUntrackedParameter<edm::InputTag> ("GenParticlesTag", edm::InputTag("genParticles"))),
34  rhoParameter_(pset.getUntrackedParameter<double> ("RhoParameter", 1.004))
35 {
36  produces<double>();
37 }
38 
41 
44 }
45 
48 }
49 
52  if (iEvent.isRealData()) return;
53 
55  iEvent.getByLabel("genParticles", genParticles);
56  unsigned int gensize = genParticles->size();
57 
58  std::auto_ptr<double> weight (new double);
59 
60  // Set default weight
61  (*weight) = 1.;
62 
63  // Only DY implemented for the time being
64  for(unsigned int i = 0; i<gensize; ++i) {
65  const reco::GenParticle& part = (*genParticles)[i];
66  int status = part.status();
67  if (status!=3) break;
68  int id = part.pdgId();
69  if (id!=23) continue;
70  double Q = part.mass();
71  unsigned int nmothers = part.numberOfMothers();
72  if (nmothers<=0) continue;
73  size_t key = part.motherRef(0).key();
74  unsigned int quark_id = abs((*genParticles)[key].pdgId());
75  if (quark_id>0 && quark_id<6) {
76  (*weight) *= sigma0_qqbarll(quark_id, Q, rhoParameter_)
77  / sigma0_qqbarll(quark_id, Q, 1.0);
78  }
79  break;
80  }
81 
82  //printf(" \t >>>>> WeakEffectsWeightProducer: Final weight = %f\n", (*weight));
83  iEvent.put(weight);
84 }
85 
87  double pigaga = -0.010449239475366825 - 0.0023228196282246765*log(q2)- 0.0288 - 0.002980*(log(q2/8464.)+0.006307*(q2/8464.-1.));
88  return (1./137.0359895) / (1.+pigaga);
89 }
90 
91 double WeakEffectsWeightProducer::sigma0_qqbarll(unsigned int quark_id, double Q, double rho) {
92  double MZ = 91.188;
93  double GZ = 2.495;
94  double sin2eff = 0.232;
95 
96  double vl = -0.5 + 2.*sin2eff;
97  double al = -0.5;
98 
99  double qq = 0.;
100  double vq = 0.;
101  double aq = 0.;
102  double alphaW = 2.7e-3 * pow(log(Q*Q/80.4/80.4),2);
103  double alphaZ = 2.7e-3 * pow(log(Q*Q/MZ/MZ),2);
104  double sudakov_factor = 1.;
105  if (abs(quark_id)%2==1) {
106  qq = -1./3.;
107  vq = -0.5 - 2.*qq*sin2eff;
108  aq = -0.5;
109  sudakov_factor = 1 + (-2.139 + 0.864)*alphaW - 0.385*alphaZ;
110  } else {
111  qq = 2./3.;
112  vq = 0.5 - 2.*qq*sin2eff;
113  aq = 0.5;
114  sudakov_factor = 1 + (-3.423 + 1.807)*alphaW - 0.557*alphaZ;
115  }
116 
117  double alfarn = alphaQED(Q*Q);
118  double zcoupl = sqrt(2.) * 1.166389e-5 * MZ*MZ / 4. / M_PI;
119  double gll = zcoupl * MZ/3. * (vl*vl + al*al);
120  double gdd = zcoupl * MZ/3. * (vq*vq + aq*aq);
121  double denom = (Q*Q-MZ*MZ)*(Q*Q-MZ*MZ)+ pow(Q,4)*GZ*GZ/MZ/MZ;
122  double qed = M_PI * qq*qq * alfarn*alfarn / Q/Q;
123  double zint = rho * 2*M_PI * zcoupl * alfarn * qq * vq*vl * (Q*Q-MZ*MZ) / denom;
124  double zonly = rho * rho * 9.*M_PI * gll * gdd / MZ/MZ * Q*Q / denom;
125 
126  return (qed + zint + zonly) * sudakov_factor;
127 }
128 
int i
Definition: DBlmapReader.cc:9
virtual int pdgId() const
PDG identifier.
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
virtual int status() const
status word
Definition: DDAxes.h:10
#define abs(x)
Definition: mlp_lapack.h:159
virtual void produce(edm::Event &, const edm::EventSetup &)
#define MZ
double sigma0_qqbarll(unsigned int quark_type, double Q, double rho)
bool isRealData() const
Definition: EventBase.h:60
double q2[4]
Definition: TauolaWrapper.h:88
virtual double mass() const
mass
int iEvent
Definition: GenABIO.cc:243
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:84
virtual size_t numberOfMothers() const
number of mothers
T sqrt(T t)
Definition: SSEVec.h:28
tuple pset
Definition: CrabTask.py:85
daughters::value_type motherRef(size_type i=0) const
reference to mother at given position
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:359
#define M_PI
Definition: BFit3D.cc:3
Log< T >::type log(const T &t)
Definition: Log.h:22
part
Definition: HCALResponse.h:21
list key
Definition: combine.py:13
tuple status
Definition: ntuplemaker.py:245
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
WeakEffectsWeightProducer(const edm::ParameterSet &pset)