ATLAS Offline Software
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 
9 namespace Pythia8{
10  class WprimeWZFlatPtMass6;
11 }
12 
13 namespace 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 
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
Pythia8::WprimeWZFlatPtMass6
Definition: WprimeWZFlatPtMass6.cxx:15
Pythia8::WprimeWZFlatPtMass6::getFlatmW
double getFlatmW(double mW)
Definition: WprimeWZFlatPtMass6.cxx:136
Pythia8::WprimeWZFlatPtMass6::canVetoProcessLevel
bool canVetoProcessLevel()
Definition: WprimeWZFlatPtMass6.cxx:52
Event
Definition: trigbs_orderedMerge.cxx:42
TRTCalib_cfilter.p1
p1
Definition: TRTCalib_cfilter.py:130
Pythia8::WprimeWZFlatPtMass6::canModifySigma
virtual bool canModifySigma()
Definition: WprimeWZFlatPtMass6.cxx:26
Pythia8::WprimeWZFlatPtMass6::m_flatpT
Pythia8_UserHooks::UserSetting< int > m_flatpT
Definition: WprimeWZFlatPtMass6.cxx:217
Pythia8_UserHooks::UserHooksFactory::Creator
Definition: UserHooksFactory.h:54
drawFromPickle.exp
exp
Definition: drawFromPickle.py:36
SUSY_SimplifiedModel_PostInclude.process
string process
Definition: SUSY_SimplifiedModel_PostInclude.py:43
xAOD::Particle
Particle_v1 Particle
Define the latest version of the particle class.
Definition: Event/xAOD/xAODParticleEvent/xAODParticleEvent/Particle.h:17
UserHooksFactory.h
Pythia8
Author: James Monk (jmonk@cern.ch)
Definition: IPythia8Custom.h:10
Pythia8_UserHooks::UserSetting< int >
TRTCalib_cfilter.p2
p2
Definition: TRTCalib_cfilter.py:131
Pythia8::WprimeWZFlatPtMass6::WprimeWZFlatPtMass6
WprimeWZFlatPtMass6()
Definition: WprimeWZFlatPtMass6.cxx:20
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
lumiFormat.i
int i
Definition: lumiFormat.py:85
Preparation.mode
mode
Definition: Preparation.py:107
UserSetting.h
Pythia8::WprimeWZFlatPtMass6::doVetoProcessLevel
bool doVetoProcessLevel(Event &process)
Definition: WprimeWZFlatPtMass6.cxx:55
Pythia8::WprimeWZFlatPtMass6::m_flatMass
Pythia8_UserHooks::UserSetting< int > m_flatMass
Definition: WprimeWZFlatPtMass6.cxx:218
Pythia8::WprimeWZFlatPtMass6::multiplySigmaBy
virtual double multiplySigmaBy(const SigmaProcess *sigmaProcessPtr, const PhaseSpace *phaseSpacePtr, bool)
Definition: WprimeWZFlatPtMass6.cxx:29
Pythia8::WprimeWZFlatPtMass6::breitWignerDenom
double breitWignerDenom(double mFrac)
Definition: WprimeWZFlatPtMass6.cxx:182
Pythia8::WprimeWZFlatPtMass6::~WprimeWZFlatPtMass6
~WprimeWZFlatPtMass6()
Definition: WprimeWZFlatPtMass6.cxx:23
Pythia8::WprimeWZFlatPtMass6Creator
Pythia8_UserHooks::UserHooksFactory::Creator< Pythia8::WprimeWZFlatPtMass6 > WprimeWZFlatPtMass6Creator("WprimeWZFlatPtMass6")
python.PyAthena.v
v
Definition: PyAthena.py:154
UserHooksUtils.h
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
Pythia8::WprimeWZFlatPtMass6::pTWeight
double pTWeight(double rH)
Definition: WprimeWZFlatPtMass6.cxx:193
TRTCalib_cfilter.p3
p3
Definition: TRTCalib_cfilter.py:132
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
TRTCalib_cfilter.p0
p0
Definition: TRTCalib_cfilter.py:129
Generated on Thu Oct 30 2025 21:22:18 for ATLAS Offline Software by   doxygen 1.8.18