ATLAS Offline Software
Loading...
Searching...
No Matches
WprimeWZFlatPtMass6.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
3*/
4#include "Pythia8_i/UserHooksFactory.h"
5#include "Pythia8/PhaseSpace.h"
6#include "UserHooksUtils.h"
7#include "UserSetting.h"
8
9namespace Pythia8{
10 class WprimeWZFlatPtMass6;
11}
12
13namespace Pythia8 {
14
15 class WprimeWZFlatPtMass6 : public UserHooks {
16
17 public:
18
19 // Constructor.
20 WprimeWZFlatPtMass6(){}
21
22 // Destructor.
23 ~WprimeWZFlatPtMass6(){}
24
25 // Allow process cross section to be modified...
26 virtual bool canModifySigma() {return true;}
27
28 // ...which gives access to the event at the trial level, before selection.
29 virtual double multiplySigmaBy(const SigmaProcess* sigmaProcessPtr,
30 const PhaseSpace* phaseSpacePtr,
31 bool /* inEvent */) {
32 // All events should be 2 -> 1, but kill them if not.
33 if (sigmaProcessPtr->nFinal() != 1) return 0.;
34
35 // Weight cross section with BW propagator, i.e. to remove it.
36 // (inEvent = false for initialization).
37 // No inEvent criteria, want weight both for cross section
38 // and MC generation.
39
40 int idRes = sigmaProcessPtr->resonanceA();
41 double mRes = particleDataPtr->m0(idRes);
42 double wRes = particleDataPtr->mWidth(idRes);
43 double m2Res = mRes*mRes;
44 double gamMRat = wRes/mRes;
45 double sHat = phaseSpacePtr->sHat();
46 double weightBW = m2Res*m2Res + sHat*sHat*(1 + gamMRat*gamMRat) - 2.*sHat*m2Res;
47 double rH = sqrt(sHat);
48
49 return (m_flatpT(settingsPtr)) ? weightBW * pTWeight(rH) : weightBW * breitWignerDenom(rH/settingsPtr->parm("Beams:eCM"));
50 }
51 //bool canVetoResonanceDecays() { return true; }
52 bool canVetoProcessLevel() { return true; }
53
54 //bool doVetoResonanceDecays(Event& process) {
55 bool doVetoProcessLevel(Event& process) {
56
57 const int mode = m_flatMass(settingsPtr);
58 //subEvent(process,true);
59 //omitResonanceDecays(process,true);
60 //process = workEvent;
61 //process.list();
62
63 double apT = 1.0;
64 //std::cout<<"before changes"<<std::endl;
65 //process.list();
66 bool isVetoed = false;
67 if(mode == 1){
68 for (int i = 1; i < process.size(); ++i) {
69 // Select vector bosons
70 Particle& v = process[i];
71
72 if (v.idAbs() != 34) continue;
73
74 // Find W/Z daughter particles
75 Particle& d_W = process[v.daughter1()];
76 Particle& d_Z = process[v.daughter2()];
77 if(d_W.idAbs() !=24 ){
78 d_W = process[v.daughter2()];
79 d_Z = process[v.daughter1()];
80 }
81 Vec4 pv_W = d_W.p();
82 Vec4 pv_Z = d_Z.p();
83 Vec4 pv_W_orig = d_W.p();
84 Vec4 pv_Z_orig = d_Z.p();
85 double pTW = sqrt(pow(pv_W.px(),2) + pow(pv_W.py(),2));
86 double mW = sqrt(pow(pv_W.e(),2) - pow(pv_W.px(),2) - pow(pv_W.py(),2) - pow(pv_W.pz(),2));
87 double mZ = sqrt(pow(pv_Z.e(),2) - pow(pv_Z.px(),2) - pow(pv_Z.py(),2) - pow(pv_Z.pz(),2));
88 double mWflat = getFlatmW(mW);
89
90 if((pow(mWflat,2) - pow(mW,2))/pow(pTW,2) < 1.0){
91 apT = sqrt(1 - (pow(mWflat,2) - pow(mW,2))/pow(pTW,2));
92 }else{
93 return true;
94 }
95
97 pv_W.px(apT*pv_W.px());
98 pv_W.py(apT*pv_W.py());
99 d_W.p(pv_W);
100 d_W.m(mWflat);
101
103 pv_Z.px(apT*pv_Z.px());
104 pv_Z.py(apT*pv_Z.py());
105 d_Z.p(pv_Z);
106 d_Z.m(sqrt(pow(pv_Z.e(),2) - pow(pv_Z.px(),2) - pow(pv_Z.py(),2) - pow(pv_Z.pz(),2)));
107 double mZref = d_Z.m();
108
109 if (mZref > 600.0) {
110 isVetoed = true;
111 }
112
113 if (mWflat > 700.0){
114 isVetoed = true;
115 }
116
117 //rescale components of the W and Z boson daughters.
118 Particle& d1_W = process[d_W.daughter1()];
119 Particle& d2_W = process[d_W.daughter2()];
120 RescaleDaughters(d1_W,d2_W,pv_W,pv_W_orig,mW,mWflat);
121
122 Particle& d1_Z = process[d_Z.daughter1()];
123 Particle& d2_Z = process[d_Z.daughter2()];
124 RescaleDaughters(d1_Z,d2_Z,pv_Z,pv_Z_orig,mZ,mZref);
125
126 }//loop over particles
127 //std::cout<<"after changes"<<std::endl;
128 //process.list();
129
130 }//mode == 1
131 return isVetoed;
132 }
133
134 private:
135
136 double getFlatmW(double mW){
137
138 double mMin = 10.;
139 double mMax = 1000.;
140
141 double mWt = ((mMax - mMin)/log(mMax/mMin))*log(mW/mMin)+mMin;
142 return mWt;
143 }
144
145 void RescaleDaughters(Particle& d1_W,Particle& d2_W,Vec4 pv_W,Vec4 pv_W_orig,double mW,double mWflat){
146
147 Vec4 pd1_W = d1_W.p();
148 Vec4 pd2_W = d2_W.p();
150 pd1_W.bstback(pv_W_orig);
151 pd2_W.bstback(pv_W_orig);
152
153 double m1_W_orig = sqrt(pow(pd1_W.e(),2) - pow(pd1_W.px(),2) - pow(pd1_W.py(),2) - pow(pd1_W.pz(),2));
154 double m2_W_orig = sqrt(pow(pd2_W.e(),2) - pow(pd2_W.px(),2) - pow(pd2_W.py(),2) - pow(pd2_W.pz(),2));
155 double mWchild = sqrt(pow(pd1_W.e()+pd2_W.e(),2)-pow(pd1_W.px()+pd2_W.px(),2)-pow(pd1_W.py()+pd2_W.py(),2)-pow(pd1_W.pz()+pd2_W.pz(),2));
156
157 double r_W_orig = 1 - (pow(m1_W_orig,2) + pow(m2_W_orig,2))/pow(mW,2);
158 double r_W_flat = 1 - (pow(m1_W_orig,2) + pow(m2_W_orig,2))/pow(mWflat,2);
159 double mR_W = (mWflat*r_W_flat)/(mWchild*r_W_orig);
160
161 //rescale daughters to ensure invariant mass equal to the new W boson mass
162 pd1_W.px(mR_W*pd1_W.px());
163 pd1_W.py(mR_W*pd1_W.py());
164 pd1_W.pz(mR_W*pd1_W.pz());
165 pd1_W.e(sqrt(pow(pd1_W.px(),2)+pow(pd1_W.py(),2)+pow(pd1_W.pz(),2)+pow(m1_W_orig,2)));
166
167 pd2_W.px(mR_W*pd2_W.px());
168 pd2_W.py(mR_W*pd2_W.py());
169 pd2_W.pz(mR_W*pd2_W.pz());
170 pd2_W.e(sqrt(pow(pd2_W.px(),2)+pow(pd2_W.py(),2)+pow(pd2_W.pz(),2)+pow(m2_W_orig,2)));
171
173 pd1_W.bst(pv_W); d1_W.p(pd1_W);
174 pd2_W.bst(pv_W); d2_W.p(pd2_W);
175 double m1_W = sqrt(pow(pd1_W.e(),2) - pow(pd1_W.px(),2) - pow(pd1_W.py(),2) - pow(pd1_W.pz(),2));
176 double m2_W = sqrt(pow(pd2_W.e(),2) - pow(pd2_W.px(),2) - pow(pd2_W.py(),2) - pow(pd2_W.pz(),2));
177
178 //rescale mass of W boson daughters
179 d1_W.m(m1_W);
180 d2_W.m(m2_W);
181 }
182 double breitWignerDenom(double mFrac){
183
184 if(mFrac < 0.025) return breitWignerDenom(0.025);
185 if(mFrac > 0.6) return breitWignerDenom(0.6);
186
187 if(mFrac < 0.0425) return 1e-12/(-1.293+1.098e+2*mFrac-2.800e+3*mFrac*mFrac+2.345e+4*mFrac*mFrac*mFrac);
188 if(mFrac < 0.073) return 1.248e-12*(exp(1.158+18.34*mFrac));
189
190 return 5.733e-10*pow(mFrac,-3.798-0.6555*log(mFrac))/pow(1.427-mFrac,30.017);
191 }
192
193 double pTWeight(double rH){
194
195 double pe0 = 9.705/2000.;
196 double pe1 = -1.27668e-03;
197
198 double weightHighpT =1./(exp(pe0+pe1*rH));
199
200 double p0 = 0.00405295;
201 double p1 = -1.15389e-06;
202 double p2 = -8.83305e-10;
203 double p3 = 1.02983e-12;
204 double p4 = -3.64486e-16;
205 double p5 = 6.05783e-20;
206 double p6 = -4.74988e-24;
207 double p7 = 1.40627e-28;
208 double weightFinal = (p0+(p1*rH)+(p2*pow(rH,2))+(p3*pow(rH,3))+(p4*pow(rH,4))+(p5*pow(rH,5))+(p6*pow(rH,6))+(p7*pow(rH,7)));
209
210 if(rH < 400.) weightFinal *= 0.5;
211
212 return weightHighpT * weightFinal;
213 }
214
215 // This switch says whether to use the old style flattening of the Breit Wigner, or additionally flatten the PT spectrum.
216 // Off by default, for consistency with old production jobs
217 Pythia8_UserHooks::UserSetting<int> m_flatpT = Pythia8_UserHooks::UserSetting<int>("WprimeWZFlatPtMass6:FlatPt", 0);
218 Pythia8_UserHooks::UserSetting<int> m_flatMass = Pythia8_UserHooks::UserSetting<int>("WprimeWZFlatPtMass6:FlatMass",0);
219
220 };
221
222 Pythia8_UserHooks::UserHooksFactory::Creator<Pythia8::WprimeWZFlatPtMass6> WprimeWZFlatPtMass6Creator("WprimeWZFlatPtMass6");
223
224} // end namespace Pythia8
225
226
pow
constexpr int pow(int base, int exp) noexcept
Definition ap_fixedTest.cxx:15
Pythia8::WprimeWZFlatPtMass6::m_flatpT
Pythia8_UserHooks::UserSetting< int > m_flatpT
Definition WprimeWZFlatPtMass6.cxx:217
Pythia8::WprimeWZFlatPtMass6::RescaleDaughters
void RescaleDaughters(Particle &d1_W, Particle &d2_W, Vec4 pv_W, Vec4 pv_W_orig, double mW, double mWflat)
Definition WprimeWZFlatPtMass6.cxx:145
Pythia8::WprimeWZFlatPtMass6::m_flatMass
Pythia8_UserHooks::UserSetting< int > m_flatMass
Definition WprimeWZFlatPtMass6.cxx:218
Pythia8::WprimeWZFlatPtMass6::doVetoProcessLevel
bool doVetoProcessLevel(Event &process)
Definition WprimeWZFlatPtMass6.cxx:55
Pythia8::WprimeWZFlatPtMass6::multiplySigmaBy
virtual double multiplySigmaBy(const SigmaProcess *sigmaProcessPtr, const PhaseSpace *phaseSpacePtr, bool)
Definition WprimeWZFlatPtMass6.cxx:29
Pythia8_UserHooks::UserHooksFactory::Creator
Definition UserHooksFactory.h:54
Pythia8_UserHooks::UserSetting
Definition UserSetting.h:17
process
const std::string process
Definition fbtTestBasics.cxx:76
Pythia8
Author: James Monk (jmonk@cern.ch)
Definition IPythia8Custom.h:10
Pythia8::WprimeWZFlatPtMass6Creator
Pythia8_UserHooks::UserHooksFactory::Creator< Pythia8::WprimeWZFlatPtMass6 > WprimeWZFlatPtMass6Creator("WprimeWZFlatPtMass6")
Generated on for ATLAS Offline Software by doxygen 1.14.0