File indexing completed on 2025-09-15 09:06:52
0001
0002
0003
0004
0005
0006
0007
0008
0009 #ifndef Pythia8_VinciaEW_H
0010 #define Pythia8_VinciaEW_H
0011
0012
0013 #include "Pythia8/Event.h"
0014 #include "Pythia8/StandardModel.h"
0015 #include "Pythia8/PartonSystems.h"
0016 #include "Pythia8/PythiaComplex.h"
0017 #include "Pythia8/BeamParticle.h"
0018 #include "Pythia8/UserHooks.h"
0019
0020
0021 #include "Pythia8/VinciaCommon.h"
0022 #include "Pythia8/VinciaQED.h"
0023
0024 namespace Pythia8 {
0025
0026
0027
0028
0029
0030 class EWParticle {
0031
0032 public:
0033
0034
0035 EWParticle() = default;
0036 EWParticle(double massIn, double widthIn, bool isResIn) :
0037 mass(massIn), width(widthIn), isRes(isResIn) {};
0038
0039
0040 double mass{0.}, width{0.};
0041 bool isRes{false};
0042
0043 };
0044
0045
0046
0047
0048
0049 class EWParticleData {
0050
0051 public:
0052
0053
0054 bool find(int id, int pol){
0055 return data.find(make_pair(id, pol)) != data.end();}
0056
0057
0058 void add(int id, int pol, double massIn, double widthIn, bool isResIn) {
0059 data[make_pair(id, pol)] = EWParticle(massIn, widthIn, isResIn);}
0060
0061
0062 double mass(int id, int pol) {
0063 return find(id, pol) ? data[make_pair(id, pol)].mass : 0.;}
0064
0065
0066 double mass(int id) {
0067
0068 if (find(id, 1)) return data[make_pair(id, 1)].mass;
0069 return find(id, 0) ? data[make_pair(id, 0)].mass : 0.;
0070 }
0071
0072
0073 double width(int id, int pol) {
0074 return find(id, pol) ? data[make_pair(id, pol)].width : 0.;}
0075
0076
0077 bool isRes(int id, int pol) {
0078 return find(id, pol) && data[make_pair(id, pol)].isRes;}
0079
0080
0081 bool isRes(int id) {
0082
0083 if (find(id, 1)) return data[make_pair(id, 1)].isRes;
0084 else if (find(id, 0)) return data[make_pair(id, 0)].isRes;
0085 else return false;
0086 }
0087
0088
0089 EWParticle* get(int id, int pol) { return &data.at(make_pair(id, pol)); }
0090 unordered_map<pair<int, int>, EWParticle>::iterator begin() {
0091 return data.begin();}
0092 unordered_map<pair<int, int>, EWParticle>::iterator end() {
0093 return data.end();}
0094
0095
0096 unordered_map<pair<int,int>, EWParticle> data;
0097
0098 };
0099
0100
0101
0102
0103
0104 class AntWrapper {
0105
0106 public:
0107
0108
0109 AntWrapper(double valIn, int poliIn, int poljIn):
0110 val(valIn), poli(poliIn), polj(poljIn) {}
0111
0112
0113 void print() {cout << "(" << poli << ", " << polj << ") " << val;}
0114
0115
0116 double val;
0117 int poli, polj;
0118
0119 };
0120
0121
0122
0123
0124
0125 class AmpWrapper {
0126
0127 public:
0128
0129
0130 AmpWrapper(complex valIn, int poliIn, int poljIn):
0131 val(valIn), poli(poliIn), polj(poljIn) {}
0132
0133
0134 AntWrapper norm() {return AntWrapper(std::norm(val), poli, polj);}
0135
0136
0137 AmpWrapper& operator+=(complex c) {this->val += c; return *this;}
0138 AmpWrapper& operator*=(complex c) {this->val *= c; return *this;}
0139 void print() {cout << "(" << poli << ", " << polj << ") " << val;}
0140
0141
0142 complex val;
0143 int poli, polj;
0144
0145 };
0146
0147
0148
0149
0150
0151 class AmpCalculator {
0152
0153 public:
0154
0155
0156 void initPtr(Info* infoPtrIn, AlphaEM* alphaPtrIn,
0157 AlphaStrong* alphaSPtrIn) {
0158 infoPtr = infoPtrIn;
0159 partonSystemsPtr = infoPtr->partonSystemsPtr;
0160 rndmPtr = infoPtr->rndmPtr;
0161 settingsPtr = infoPtr->settingsPtr;
0162 loggerPtr = infoPtr->loggerPtr;
0163 alphaPtr = alphaPtrIn;
0164 alphaSPtr = alphaSPtrIn;
0165 isInitPtr = true;
0166 }
0167
0168
0169 void init(EWParticleData* dataIn, unordered_map< pair<int, int>,
0170 vector<pair<int, int> > >* cluMapFinalIn,
0171 unordered_map< pair<int, int>, vector<pair<int, int> > >*
0172 cluMapInitialIn);
0173
0174
0175 void setVerbose(int verboseIn) {verbose = verboseIn;}
0176
0177
0178 complex spinProd(int pol, const Vec4& ka, const Vec4& kb);
0179 complex spinProd(int pol, const Vec4& ka, const Vec4& pa, const Vec4& kb);
0180 complex spinProd(int pol, const Vec4& ka, const Vec4& pa, const Vec4& pb,
0181 const Vec4& kb);
0182 complex spinProd(int pol, const Vec4& ka, const Vec4& pa, const Vec4& pb,
0183 const Vec4& pc, const Vec4& kb);
0184 complex spinProd(int pol, const Vec4& ka, const Vec4& pa, const Vec4& pb,
0185 const Vec4& pc, const Vec4& pd, const Vec4& kb);
0186 Vec4 spinProdFlat(string method, const Vec4& ka, const Vec4& pa);
0187
0188
0189 void initCoup(bool va, int id1, int id2, int pol, bool m);
0190
0191
0192 void initFSRAmp(bool va, int id1, int id2, int pol,
0193 const Vec4& pi, const Vec4 &pj, const double& mMot, const double& widthQ2);
0194
0195
0196 bool zdenFSRAmp(const string& method, const Vec4& pi, const Vec4& pj,
0197 bool check);
0198
0199
0200 void initISRAmp(bool va, int id1, int id2, int pol,
0201 const Vec4& pa, const Vec4 &pj, double& mA);
0202
0203
0204 bool zdenISRAmp(const string& method, const Vec4& pa, const Vec4& pj,
0205 bool check);
0206
0207
0208
0209
0210
0211 complex ftofvFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0212 int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
0213 complex ftofhFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0214 int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
0215 complex fbartofbarvFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0216 int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
0217 complex fbartofbarhFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0218 int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
0219 complex vTtoffbarFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0220 int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
0221 complex vTtovhFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0222 int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
0223 complex vTtovvFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0224 int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
0225 complex vLtoffbarFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0226 int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
0227 complex vLtovhFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0228 int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
0229 complex vLtovvFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0230 int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
0231 complex htoffbarFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0232 int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
0233 complex htovvFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0234 int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
0235 complex htohhFSRAmp(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0236 int idj, double mMot, double widthQ2, int polMot, int poli, int polj);
0237
0238
0239 complex ftofvISRAmp(const Vec4& pa, const Vec4& pj, int idA, int ida,
0240 int idj, double mA, int polA, int pola, int polj);
0241 complex ftofhISRAmp(const Vec4& pa, const Vec4& pj, int idA, int ida,
0242 int idj, double mA, int polA, int pola, int polj);
0243 complex fbartofbarvISRAmp(const Vec4& pa, const Vec4& pj, int idA, int ida,
0244 int idj, double mA, int polA, int pola, int polj);
0245 complex fbartofbarhISRAmp(const Vec4& pa, const Vec4& pj, int idA, int ida,
0246 int idj, double mA, int polA, int pola, int polj);
0247
0248
0249 complex branchAmpFSR(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0250 int idj, double mMot, double widthQ2, int polMot, int poli=9, int polj=9);
0251 complex branchAmpISR(const Vec4& pa, const Vec4& pj, int idA, int ida,
0252 int idj, double mA, int polA, int pola=9, int polj=9);
0253
0254
0255 double branchKernelFF(const Vec4& pi, const Vec4& pj, int idMot, int idi,
0256 int idj, double mMot, double widthQ2, int polMot, int poli, int polj) {
0257 return norm(branchAmpFSR(pi, pj, idMot, idi, idj, mMot, widthQ2,
0258 polMot, poli, polj));}
0259
0260
0261 vector<AntWrapper> branchKernelFF(Vec4 pi, Vec4 pj, int idMot, int idi,
0262 int idj, double mMot, double widthQ2, int polMot);
0263
0264
0265 double branchKernelII(Vec4 pa, Vec4 pj, int idA, int ida,
0266 int idj, double mA, int polA, int pola, int polj) {
0267 return norm(branchAmpISR(pa, pj, idA, ida, idj, mA, polA, pola, polj));}
0268
0269
0270 vector<AntWrapper> branchKernelII(Vec4 pa, Vec4 pj, int idA, int ida,
0271 int idj, double mA, int polA);
0272
0273
0274 void initFFAnt(bool va, int id1, int id2, int pol, const double& Q2,
0275 const double& widthQ2, const double& xi, const double& xj,
0276 const double& mMot, const double& miIn, const double& mjIn);
0277
0278
0279 void hmsgFFAnt(int polMot, int poli, int polj) {
0280 stringstream ss;
0281 ss << "helicity combination was not found:\n "
0282 << "polMot = " << polMot << " poli = " << poli << " polj = " << polj;
0283 loggerPtr->ERROR_MSG(ss.str());}
0284
0285
0286 void initIIAnt(int id1, int id2, int pol, const double& Q2,
0287 const double& xA, const double& xj,
0288 const double& mA, const double& maIn, const double& mjIn);
0289
0290
0291 void hmsgIIAnt(int polA, int pola, int polj) {
0292 stringstream ss;
0293 ss << "helicity combination was not found:\n "
0294 << "polA = " << polA << " pola = " << pola << " polj = " << polj;
0295 loggerPtr->ERROR_MSG(ss.str());}
0296
0297
0298
0299
0300 double ftofvFFAnt(double Q2, double widthQ2, double xi, double xj,
0301 int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
0302 int polMot, int poli, int polj);
0303 double ftofhFFAnt(double Q2, double widthQ2, double xi, double xj,
0304 int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
0305 int polMot, int poli, int polj);
0306 double fbartofbarvFFAnt(double Q2, double widthQ2, double xi, double xj,
0307 int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
0308 int polMot, int poli, int polj);
0309 double fbartofbarhFFAnt(double Q2, double widthQ2, double xi, double xj,
0310 int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
0311 int polMot, int poli, int polj);
0312 double vtoffbarFFAnt(double Q2, double widthQ2, double xi, double xj,
0313 int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
0314 int polMot, int poli, int polj);
0315 double vtovhFFAnt(double Q2, double widthQ2, double xi, double xj,
0316 int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
0317 int polMot, int poli, int polj);
0318 double vtovvFFAnt(double Q2, double widthQ2, double xi, double xj,
0319 int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
0320 int polMot, int poli, int polj);
0321 double htoffbarFFAnt(double Q2, double widthQ2, double xi, double xj,
0322 int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
0323 int polMot, int poli, int polj);
0324 double htovvFFAnt(double Q2, double widthQ2, double xi, double xj,
0325 int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
0326 int polMot, int poli, int polj);
0327 double htohhFFAnt(double Q2, double widthQ2, double xi, double xj,
0328 int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
0329 int polMot, int poli, int polj);
0330
0331
0332
0333
0334 double ftofvIIAnt(double Q2, double xA, double xj,
0335 int idA, int ida, int idj, double mA, double maIn, double mjIn,
0336 int polA, int pola, int polj);
0337 double fbartofbarvIIAnt(double Q2, double xA, double xj,
0338 int idA, int ida, int idj, double mA, double maIn, double mjIn,
0339 int polA, int pola, int polj);
0340
0341
0342 double antFuncFF(double Q2, double widthQ2, double xi, double xj,
0343 int idMot, int idi, int idj, double mMot, double miIn, double mjIn,
0344 int polMot, int poli, int polj);
0345
0346
0347 vector<AntWrapper> antFuncFF(double Q2, double widthQ2, double xi,
0348 double xj, int idMot, int idi, int idj, double mMot, double miIn,
0349 double mjIn, int polMot);
0350
0351
0352 double antFuncII(double Q2, double xA, double xj,
0353 int idA, int ida, int idj, double mA, double maIn, double mjIn,
0354 int polA, int pola, int polj);
0355
0356
0357 vector<AntWrapper> antFuncII(double Q2, double xA, double xj,
0358 int idA, int ida, int idj, double mA, double maIn, double mjIn, int polA);
0359
0360
0361 void initFSRSplit(bool va, int id1, int id2, int pol,
0362 const double& mMot, const double& miIn, const double& mjIn) {
0363 mi = miIn; mj = mjIn; mMot2 = pow2(mMot); mi2 = pow2(mi); mj2 = pow2(mj);
0364 initCoup(va, id1, id2, pol, true);}
0365
0366
0367 bool zdenFSRSplit(const string& method, const double& Q2, const double& z,
0368 bool check);
0369
0370
0371 void hmsgFSRSplit(int polMot, int poli, int polj) {
0372 stringstream ss;
0373 ss << "helicity combination was not found:\n "
0374 << "polMot = " << polMot << " poli = " << poli << " polj = " << polj;
0375 loggerPtr->ERROR_MSG(ss.str());}
0376
0377
0378 void initISRSplit(bool va, int id1, int id2, int pol,
0379 const double& mA, const double& maIn, const double& mjIn) {
0380 ma = maIn; mj = mjIn; mA2 = pow2(mA); ma2 = pow2(ma); mj2 = pow2(mj);
0381 initCoup(va, id1, id2, pol, mA > VinciaConstants::NANO);}
0382
0383
0384 bool zdenISRSplit(const string& method, const double& Q2, const double& z,
0385 bool flip, bool check);
0386
0387
0388 void hmsgISRSplit(int polA, int pola, int polj) {
0389 stringstream ss;
0390 ss << "helicity combination was not found:\n "
0391 << "polA = " << polA << " pola = " << pola << " polj = " << polj;
0392 loggerPtr->ERROR_MSG(ss.str());}
0393
0394
0395 double ftofvFSRSplit(double Q2, double z, int idMot, int idi, int idj,
0396 double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
0397 double ftofhFSRSplit(double Q2, double z, int idMot, int idi, int idj,
0398 double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
0399 double fbartofbarvFSRSplit(double Q2, double z, int idMot, int idi, int idj,
0400 double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
0401 double fbartofbarhFSRSplit(double Q2, double z, int idMot, int idi, int idj,
0402 double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
0403 double vTtoffbarFSRSplit(double Q2, double z, int idMot, int idi, int idj,
0404 double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
0405 double vTtovhFSRSplit(double Q2, double z, int idMot, int idi, int idj,
0406 double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
0407 double vTtovvFSRSplit(double Q2, double z, int idMot, int idi, int idj,
0408 double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
0409 double vLtoffbarFSRSplit(double Q2, double z, int idMot, int idi, int idj,
0410 double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
0411 double vLtovhFSRSplit(double Q2, double z, int idMot, int idi, int idj,
0412 double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
0413 double vLtovvFSRSplit(double Q2, double z, int idMot, int idi, int idj,
0414 double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
0415 double htoffbarFSRSplit(double Q2, double z, int idMot, int idi, int idj,
0416 double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
0417 double htovvFSRSplit(double Q2, double z, int idMot, int idi, int idj,
0418 double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
0419 double htohhFSRSplit(double Q2, double z, int idMot, int idi, int idj,
0420 double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
0421
0422
0423 double ftofvISRSplit(double Q2, double z, int idA, int ida, int idj,
0424 double mA, double maIn, double mjIn, int polA, int pola, int polj);
0425 double ftofhISRSplit(double Q2, double z, int idA, int ida, int idj,
0426 double mA, double maIn, double mjIn, int polA, int pola, int polj);
0427 double fbartofbarvISRSplit(double Q2, double z, int idA, int ida, int idj,
0428 double mA, double maIn, double mjIn, int polA, int pola, int polj);
0429 double fbartofbarhISRSplit(double Q2, double z, int idA, int ida, int idj,
0430 double mA, double maIn, double mjIn, int polA, int pola, int polj);
0431
0432
0433 double splitFuncFSR(double Q2, double z, int idMot, int idi, int idj,
0434 double mMot, double miIn, double mjIn, int polMot, int poli, int polj);
0435 double splitFuncISR(double Q2, double z, int idA, int ida, int idj,
0436 double mA, double maIn, double mjIn, int polA, int pola, int polj);
0437
0438
0439 double getPartialWidth(int idMot, int idi, int idj, double mMot, int polMot);
0440
0441
0442 double getTotalWidth(int idMot, double mMot, int polMot);
0443
0444
0445 double getBreitWigner(int id, double m, int pol);
0446 double getBreitWignerOverestimate(int id, double m, int pol);
0447
0448
0449 double sampleMass(int id, int pol);
0450
0451
0452 void applyBosonInterferenceFactor(Event &event, int XYEv, Vec4 pi, Vec4 pj,
0453 int idi, int idj, int poli, int polj);
0454
0455
0456 bool polarise(vector<Particle> &state);
0457
0458
0459 double eventWeight() {return eventWeightSave;}
0460 void eventWeight(double eventWeightIn) {eventWeightSave = eventWeightIn;}
0461
0462
0463
0464
0465 EWParticleData* dataPtr{};
0466
0467
0468
0469 unordered_map<pair<int, int>, double> vMap, aMap, gMap, vCKM;
0470
0471 private:
0472
0473
0474
0475
0476 unordered_map<int, vector<double> > cBW;
0477
0478 unordered_map<pair<int,int>, double> nBW, np;
0479
0480
0481 double eventWeightSave{1};
0482
0483
0484 double mw, mw2, sw, sw2;
0485
0486
0487 int bwMatchMode;
0488
0489
0490 unordered_map<pair<int, int>, vector<pair<int, int> > >* cluMapFinal{};
0491 unordered_map<pair<int, int>, vector<pair<int, int> > >* cluMapInitial{};
0492
0493
0494 vector<int> fermionPols, vectorPols, scalarPols;
0495
0496
0497 double v, a, vPls, vMin, g;
0498
0499
0500 double mMot2, mi, mi2, mj, mj2, mA2, ma, ma2, isrQ2;
0501 complex M, fsrQ2;
0502
0503
0504 Vec4 kij, ki, kj, pij, kaj, ka, paj;
0505 double wij, wi, wj, wij2, wi2, wj2, waj, wa, waj2, wa2;
0506
0507
0508 double Q4, Q4gam, Q2til, ant;
0509
0510
0511 Info* infoPtr{};
0512 PartonSystems* partonSystemsPtr{};
0513 Rndm* rndmPtr{};
0514 Settings* settingsPtr{};
0515 Logger* loggerPtr{};
0516 AlphaEM* alphaPtr{};
0517 AlphaStrong* alphaSPtr{};
0518
0519
0520 bool isInit{false};
0521 bool isInitPtr{false};
0522 int verbose{0};
0523 };
0524
0525
0526
0527
0528
0529 class EWBranching {
0530
0531 public:
0532
0533
0534 EWBranching(int idMotIn, int idiIn, int idjIn, int polMotIn,
0535 double c0In = 0, double c1In = 0, double c2In = 0, double c3In = 0):
0536 idMot(idMotIn), idi(idiIn), idj(idjIn), polMot(polMotIn), c0(c0In),
0537 c1(c1In), c2(c2In), c3(c3In),
0538 isSplitToFermions(abs(idMot) > 20 && abs(idi) < 20 && abs(idj) < 20) {;}
0539
0540
0541 void print() {cout <<" (" << idMot << ", " << polMot << ") -> " << idi <<
0542 "," << idj << ": (" << c0 << ", " << c1 << ", " << c2 << ", " << c3 <<
0543 ") \n";}
0544
0545
0546 int idMot, idi, idj, polMot;
0547
0548 double c0, c1, c2, c3;
0549
0550 bool isSplitToFermions;
0551
0552 };
0553
0554
0555
0556
0557
0558 class EWAntenna {
0559
0560 public:
0561
0562
0563 EWAntenna(): iSys(-1), shat(0), doBosonInterference(false) {};
0564 virtual ~EWAntenna() = default;
0565
0566
0567 void print() {
0568 stringstream ss;
0569 ss << "Brancher = (" << iMot << ", " << polMot
0570 << "), Recoiler = " << iRec;
0571 printOut(__METHOD_NAME__, ss.str());
0572 for (int i = 0; i < (int)brVec.size(); i++) brVec[i].print();}
0573
0574
0575 void initPtr(Info* infoPtrIn, VinciaCommon* vinComPtrIn, AlphaEM* alphaPtrIn,
0576 AmpCalculator* ampCalcPtrIn) {
0577 infoPtr = infoPtrIn;
0578 rndmPtr = infoPtr->rndmPtr;
0579 loggerPtr = infoPtr->loggerPtr;
0580 partonSystemsPtr = infoPtr->partonSystemsPtr;
0581 vinComPtr = vinComPtrIn;
0582 alphaPtr = alphaPtrIn;
0583 ampCalcPtr = ampCalcPtrIn;
0584 };
0585
0586
0587 virtual bool init(Event &event, int iMotIn, int iRecIn, int iSysIn,
0588 vector<EWBranching> &branchings, Settings* settingsPtr) = 0;
0589
0590
0591 void setVerbose(int verboseIn) {verbose = verboseIn;}
0592
0593
0594 virtual double generateTrial(double q2Start, double q2End,
0595 double alphaIn) = 0;
0596
0597
0598 virtual bool acceptTrial(Event& event) = 0;
0599
0600
0601 virtual void updateEvent(Event& event) = 0;
0602 virtual void updatePartonSystems(Event& event);
0603
0604
0605 int getIndexMot() {return iMot;};
0606 int getIndexRec() {return iRec;};
0607
0608
0609 bool isSplitToFermions() {
0610 return brTrial != nullptr && brTrial->isSplitToFermions;}
0611 virtual bool isInitial() {return false;}
0612 virtual bool isResonanceDecay() {return false;}
0613
0614
0615 bool selectChannel(int idx, const double& cSum, const map<double, int>&
0616 cSumSoFar, int& idi, int& idj, double& mi2, double& mj2);
0617
0618 protected:
0619
0620
0621 int iMot, iRec, idMot, idRec, polMot;
0622
0623 Vec4 pMot, pRec;
0624
0625 double sAnt, mMot, mMot2, mRec, mRec2;
0626
0627 double alpha;
0628
0629 int iSys;
0630
0631
0632 vector<EWBranching> brVec;
0633
0634
0635 bool hasTrial;
0636 double q2Trial, sijTrial, sjkTrial;
0637 int poliTrial, poljTrial;
0638
0639
0640 vector<Vec4> pNew;
0641
0642
0643 double c0Sum, c1Sum, c2Sum, c3Sum;
0644 map<double, int> c0SumSoFar, c1SumSoFar, c2SumSoFar, c3SumSoFar;
0645
0646
0647 double q2Match;
0648
0649
0650 int jNew;
0651 unordered_map<int,int> iReplace;
0652 double shat;
0653
0654
0655 EWBranching* brTrial{};
0656 Info* infoPtr{};
0657 Rndm* rndmPtr{};
0658 Logger* loggerPtr{};
0659 PartonSystems* partonSystemsPtr{};
0660 VinciaCommon* vinComPtr{};
0661 AlphaEM* alphaPtr{};
0662 AmpCalculator* ampCalcPtr{};
0663
0664
0665 bool doBosonInterference;
0666
0667
0668 int verbose;
0669
0670 };
0671
0672
0673
0674
0675
0676 class EWAntennaFF : public EWAntenna {
0677 public:
0678
0679
0680 virtual bool init(Event &event, int iMotIn, int iRecIn, int iSysIn,
0681 vector<EWBranching> &branchings, Settings* settingsPtr) override;
0682 virtual double generateTrial(double q2Start, double q2End, double alphaIn)
0683 override;
0684 virtual bool acceptTrial(Event &event) override;
0685 virtual void updateEvent(Event &event) override;
0686
0687 private:
0688
0689
0690 double mAnt2, sqrtKallen;
0691
0692 int kMapFinal;
0693
0694 bool vetoResonanceProduction;
0695
0696 };
0697
0698
0699
0700
0701
0702 class EWAntennaFFres : public EWAntennaFF {
0703
0704 public:
0705
0706
0707 bool init(Event &event, int iMotIn, int iRecIn, int iSysIn,
0708 vector<EWBranching> &branchings, Settings* settingsPtr) override;
0709 bool isResonanceDecay() override {return true;}
0710 double generateTrial(double q2Start, double q2End, double alphaIn) override;
0711 bool acceptTrial(Event &event) override;
0712 void updateEvent(Event &event) override;
0713
0714
0715 bool genForceDecay(Event &event);
0716
0717 private:
0718
0719
0720 bool trialIsResDecay;
0721
0722 int bwMatchMode;
0723
0724 double q2Dec, q2EW;
0725
0726 bool doDecayOnly{false};
0727
0728 };
0729
0730
0731
0732
0733
0734 class EWAntennaII : public EWAntenna {
0735
0736 public:
0737
0738
0739 EWAntennaII(BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn):
0740 beamAPtr(beamAPtrIn), beamBPtr(beamBPtrIn), shh(0), xMot(0), xRec(0),
0741 vetoResonanceProduction(false), TINYPDFtrial(1e-10) {;}
0742
0743
0744 bool init(Event &event, int iMotIn, int iRecIn, int iSysIn,
0745 vector<EWBranching> &branchings,Settings* settingsPtr) override;
0746 bool isInitial() override {return true;}
0747 double generateTrial(double q2Start, double q2End, double alphaIn) override;
0748 bool acceptTrial(Event &event) override;
0749 void updateEvent(Event &event) override;
0750 void updatePartonSystems(Event &event) override;
0751
0752 private:
0753
0754
0755
0756 BeamParticle* beamAPtr{};
0757 BeamParticle* beamBPtr{};
0758
0759 double shh;
0760
0761 double xMot, xRec;
0762
0763 bool vetoResonanceProduction;
0764
0765 vector<Vec4> pRecVec;
0766 vector<int> iRecVec;
0767
0768
0769 double TINYPDFtrial;
0770
0771 };
0772
0773
0774
0775
0776
0777 class EWSystem {
0778
0779 public:
0780
0781
0782 EWSystem(): antTrial(nullptr) {clearLastTrial();}
0783 EWSystem(unordered_map<pair<int, int>, vector<EWBranching> >* brMapFinalIn,
0784 unordered_map<pair<int, int>, vector<EWBranching> >* brMapInitialIn,
0785 unordered_map<pair<int, int>, vector<EWBranching> >* brMapResonanceIn,
0786 unordered_map<pair<int, int>, vector<pair<int, int> > >* cluMapFinalIn,
0787 unordered_map<pair<int, int>, vector<pair<int, int> > >* cluMapInitialIn,
0788 AmpCalculator * ampCalcIn):
0789 shh(0), iSysSav(0), resDecOnlySav(false), q2Cut(0), q2Trial(0),
0790 lastWasSplitSav(false), lastWasDecSav(false), lastWasInitialSav(false),
0791 lastWasBelowCut(false), ISav(0), KSav(0), brMapFinal(brMapFinalIn),
0792 brMapInitial(brMapInitialIn), brMapResonance(brMapResonanceIn),
0793 cluMapFinal(cluMapFinalIn), cluMapInitial(cluMapInitialIn),
0794 ampCalcPtr(ampCalcIn), isInit(false), doVetoHardEmissions(false),
0795 verbose(false), vetoHardEmissionsDeltaR2(0) {clearLastTrial();}
0796
0797
0798 void initPtr(Info* infoPtrIn, VinciaCommon* vinComPtrIn, AlphaEM* alIn) {
0799 infoPtr = infoPtrIn;
0800 partonSystemsPtr = infoPtr->partonSystemsPtr;
0801 rndmPtr = infoPtr->rndmPtr;
0802 settingsPtr = infoPtr->settingsPtr;
0803 loggerPtr = infoPtr->loggerPtr;
0804 vinComPtr = vinComPtrIn;
0805 al = alIn;}
0806
0807
0808 void init(BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn) {
0809 beamAPtr = beamAPtrIn; beamBPtr = beamBPtrIn;
0810 doVetoHardEmissions = settingsPtr->flag("Vincia:EWoverlapVeto");
0811 vetoHardEmissionsDeltaR2 =
0812 pow2(settingsPtr->parm("Vincia:EWoverlapVetoDeltaR"));
0813 isInit = true;}
0814
0815
0816 void setVerbose(int verboseIn) {verbose = verboseIn;}
0817
0818
0819 bool prepare(Event &event, int iSysIn, double q2CutIn, bool resDecOnlyIn) {
0820 iSysSav = iSysIn; resDecOnlySav = resDecOnlyIn; q2Cut = q2CutIn;
0821 shh = infoPtr->s(); return buildSystem(event);}
0822
0823
0824 bool buildSystem(Event &event);
0825
0826
0827 int system() {return iSysSav;}
0828
0829
0830 double q2Next(double q2Start, double q2End);
0831
0832
0833 template <class T> void addAntenna(T ant, vector<T>& antVec,
0834 Event &event, int iMot, int iRec,
0835 unordered_map<pair<int, int>, vector<EWBranching> >* brMapPtr) {
0836 if (iMot == 0) return;
0837 int idA(event[iMot].id()), polA(event[iMot].pol());
0838 if (idA == 21) return;
0839 auto it = brMapPtr->find(make_pair(idA, polA));
0840 if (it != brMapPtr->end()) {
0841
0842 ant.setVerbose(verbose);
0843
0844 ant.initPtr(infoPtr, vinComPtr, al, ampCalcPtr);
0845
0846 if (ant.init(event, iMot, iRec, iSysSav, it->second, settingsPtr)) {
0847 antVec.push_back(std::move(ant));
0848 if (verbose >= VinciaConstants::DEBUG) {
0849 stringstream ss;
0850 ss << "Added EW antenna with iEv = "
0851 << iMot << " and iRec = "<< iRec<< " in system "<< iSysSav;
0852 printOut(__METHOD_NAME__, ss.str());}}}}
0853
0854
0855 template <class T> void generateTrial(vector<T> & antVec, double q2Start,
0856 double q2End, double alphaIn) {
0857 if (q2End > q2Start) return;
0858 for (int i = 0; i < (int)antVec.size(); i++) {
0859
0860 double q2New = antVec[i].generateTrial(q2Start,q2End,alphaIn);
0861
0862 if (q2New > q2Trial && q2New > q2End) {
0863
0864 q2Trial = q2New;
0865 antTrial = &(antVec[i]);
0866 lastWasDecSav = antTrial->isResonanceDecay();
0867 lastWasInitialSav = antTrial->isInitial();
0868
0869
0870 lastWasSplitSav = lastWasDecSav ? true : antTrial->isSplitToFermions();
0871 lastWasBelowCut = (q2Trial < q2Cut)? true : false;
0872 ISav = antTrial->getIndexMot(); KSav = antTrial->getIndexRec();}}}
0873
0874
0875 void generateTrial(Event &event, vector<EWAntenna> & antVec,double q2Start,
0876 double q2End, double alphaIn);
0877
0878
0879 bool acceptTrial(Event &event) {
0880 bool passed = antTrial->acceptTrial(event);
0881 if (verbose >= VinciaConstants::DEBUG)
0882 printOut(__METHOD_NAME__, passed ? "Passed veto" : "Vetoed branching");
0883 return passed;}
0884
0885
0886 void updateEvent(Event &event) {
0887 if (verbose >= VinciaConstants::DEBUG)
0888 printOut(__METHOD_NAME__, "begin", VinciaConstants::DASHLEN);
0889 if (antTrial != nullptr) antTrial->updateEvent(event);
0890 else loggerPtr->ERROR_MSG("trial doesn't exist!");
0891 if (verbose >= VinciaConstants::DEBUG)
0892 printOut(__METHOD_NAME__, "end", VinciaConstants::DASHLEN);}
0893
0894
0895 void updatePartonSystems(Event &event) {
0896 if (verbose >= VinciaConstants::DEBUG)
0897 printOut(__METHOD_NAME__, "begin", VinciaConstants::DASHLEN);
0898 if (antTrial!=nullptr) antTrial->updatePartonSystems(event);
0899 else loggerPtr->ERROR_MSG("trial doesn't exist!");
0900 if (verbose >= VinciaConstants::DEBUG)
0901 printOut(__METHOD_NAME__, "end", VinciaConstants::DASHLEN);}
0902
0903
0904 void printAntennae() {
0905 for (int i = 0; i < (int)antVecFinal.size(); i++) antVecFinal[i].print();
0906 for (int i = 0; i < (int)antVecRes.size(); i++) antVecRes[i].print();
0907 for (int i = 0; i < (int)antVecInitial.size(); i++)
0908 antVecInitial[i].print();}
0909
0910
0911 unsigned int nBranchers() {return antVecFinal.size() +
0912 antVecInitial.size() + antVecRes.size();}
0913 unsigned int nResDec() {return antVecRes.size();}
0914
0915
0916 bool hasTrial() {return antTrial != nullptr;}
0917
0918
0919 bool lastWasSplitToFermions() {return lastWasSplitSav;}
0920 bool lastWasInitial() {return lastWasInitialSav;}
0921 bool lastWasResonanceDecay() {return lastWasDecSav;}
0922
0923
0924 void clearAntennae() {antVecFinal.clear(); antVecInitial.clear();
0925 antVecRes.clear(); clearLastTrial();}
0926
0927
0928 void clearLastTrial() {
0929 q2Trial = 0; antTrial = nullptr; lastWasSplitSav = false;
0930 lastWasDecSav = false; lastWasInitialSav = false; lastWasBelowCut = false;
0931 ISav = 0; KSav = 0;}
0932
0933 private:
0934
0935
0936
0937 double shh;
0938
0939
0940 int iSysSav;
0941 bool resDecOnlySav;
0942
0943
0944 double q2Cut;
0945
0946
0947 BeamParticle* beamAPtr{};
0948 BeamParticle* beamBPtr{};
0949 Info* infoPtr{};
0950 PartonSystems* partonSystemsPtr{};
0951 Rndm* rndmPtr{};
0952 Settings* settingsPtr{};
0953 Logger* loggerPtr{};
0954 VinciaCommon* vinComPtr{};
0955 AlphaEM* al{};
0956
0957
0958 vector<EWAntennaFF> antVecFinal;
0959 vector<EWAntennaII> antVecInitial;
0960 vector<EWAntennaFFres> antVecRes;
0961
0962
0963 EWAntenna* antTrial{};
0964 double q2Trial;
0965 bool lastWasSplitSav, lastWasDecSav, lastWasInitialSav, lastWasBelowCut;
0966 int ISav, KSav;
0967
0968
0969 unordered_map<pair<int, int>, vector<EWBranching> >
0970 *brMapFinal{}, *brMapInitial{}, *brMapResonance{};
0971
0972
0973 unordered_map<pair<int, int>, vector<pair<int, int> > >
0974 *cluMapFinal{}, *cluMapInitial{};
0975
0976
0977 AmpCalculator* ampCalcPtr{};
0978
0979
0980 bool isInit, doVetoHardEmissions;
0981 int verbose;
0982 double vetoHardEmissionsDeltaR2;
0983
0984 };
0985
0986
0987
0988
0989
0990 class VinciaEW : public VinciaModule {
0991
0992 public:
0993
0994
0995 VinciaEW(): isLoaded{false} {;}
0996
0997
0998 void initPtr(Info* infoPtrIn, VinciaCommon* vinComPtrIn) override {
0999 infoPtr = infoPtrIn;
1000 particleDataPtr = infoPtr->particleDataPtr;
1001 loggerPtr = infoPtr->loggerPtr;
1002 partonSystemsPtr = infoPtr->partonSystemsPtr;
1003 rndmPtr = infoPtr->rndmPtr;
1004 settingsPtr = infoPtr->settingsPtr;
1005 vinComPtr = vinComPtrIn;
1006 ampCalc.initPtr(infoPtr, &al, &vinComPtr->alphaStrong);
1007 isInitPtr = true;}
1008
1009
1010 void init(BeamParticle* beamAPtrIn = nullptr,
1011 BeamParticle* beamBPtrIn = nullptr) override;
1012
1013
1014 bool polarise(vector<Particle> &state) override {
1015 if (isLoaded) return ampCalc.polarise(state);
1016 else return false;}
1017
1018
1019
1020 bool prepare(int iSysIn, Event &event, int scaleRegionIn = 0) override;
1021
1022
1023 void update(Event &event, int iSysIn) override {
1024 if (verbose >= VinciaConstants::DEBUG)
1025 printOut(__METHOD_NAME__, "begin", VinciaConstants::DASHLEN);
1026 if (iSysIn != ewSystem.system()) return;
1027 else ewSystem.buildSystem(event);
1028 if (verbose >= VinciaConstants::DEBUG)
1029 printOut(__METHOD_NAME__, "end", VinciaConstants::DASHLEN);}
1030
1031
1032 void setVerbose(int verboseIn) override {
1033 verbose = verboseIn; ampCalc.setVerbose(verboseIn);}
1034
1035
1036 void load() override;
1037
1038
1039 double q2Next(Event&, double q2Start,double q2End) override;
1040
1041
1042 bool lastIsSplitting() override { return ewSystem.lastWasSplitToFermions(); }
1043 bool lastIsInitial() override { return ewSystem.lastWasInitial(); }
1044 bool lastIsResonanceDecay() override {
1045 return ewSystem.lastWasResonanceDecay(); }
1046
1047
1048 bool acceptTrial(Event& event) override {
1049 if (verbose >= VinciaConstants::DEBUG)
1050 printOut(__METHOD_NAME__, "begin", VinciaConstants::DASHLEN);
1051 bool success = false;
1052 if (ewSystem.hasTrial()) success = ewSystem.acceptTrial(event);
1053 else loggerPtr->ERROR_MSG("trial doesn't exist!");
1054 if (verbose >= VinciaConstants::DEBUG)
1055 printOut(__METHOD_NAME__, "end", VinciaConstants::DASHLEN);
1056 return success;}
1057
1058
1059 void updateEvent(Event& event) override {
1060 if (verbose >= VinciaConstants::DEBUG)
1061 printOut(__METHOD_NAME__, "begin", VinciaConstants::DASHLEN);
1062 if (ewSystem.hasTrial()) ewSystem.updateEvent(event);
1063 else loggerPtr->ERROR_MSG("trial doesn't exist!");
1064 if (verbose >=VinciaConstants::DEBUG) {
1065 printOut(__METHOD_NAME__,"Event after update:"); event.list();
1066 printOut(__METHOD_NAME__, "end", VinciaConstants::DASHLEN);}}
1067
1068
1069 void updatePartonSystems(Event& event) override {
1070 if (verbose >= VinciaConstants::DEBUG)
1071 printOut(__METHOD_NAME__, "begin", VinciaConstants::DASHLEN);
1072 if (ewSystem.hasTrial()) ewSystem.updatePartonSystems(event);
1073 else loggerPtr->ERROR_MSG("trial doesn't exist!");
1074 if (verbose >= VinciaConstants::DEBUG)
1075 printOut(__METHOD_NAME__, "end", VinciaConstants::DASHLEN);}
1076
1077
1078 void clear(int) override {
1079 ewSystem.clearAntennae();
1080 ampCalc.eventWeight(1);
1081 }
1082
1083
1084 unsigned int nBranchers() override {return ewSystem.nBranchers();}
1085 unsigned int nResDec() override {return ewSystem.nResDec();}
1086
1087
1088 void printData();
1089
1090
1091 void printBranchings();
1092
1093
1094 double q2min() override {return q2minSav;}
1095 double q2minColoured() override {return q2minSav;}
1096 double eventWeight() {return ampCalc.eventWeight();}
1097
1098
1099 void q2min(double q2minSavIn) {q2minSav = q2minSavIn;}
1100
1101
1102 int sysWin() override {return ewSystem.system();}
1103
1104
1105 bool isResonance(int pid) {return ewData.isRes(pid);}
1106
1107
1108
1109
1110 unordered_map<pair<int, int>, vector<pair<int, int> > >
1111 cluMapFinal, cluMapInitial;
1112
1113
1114
1115 unordered_map<pair<int, int>, vector<EWBranching> >
1116 brMapFinal, brMapInitial, brMapResonance;
1117
1118
1119 EWParticleData ewData;
1120
1121
1122 AmpCalculator ampCalc;
1123
1124 private:
1125
1126
1127
1128 bool readFile(string file);
1129 bool readLine(string line);
1130 bool attributeValue(string line, string attribute, string& val);
1131 template <class T> bool attributeValue(
1132 string line, string attribute, T& val) {
1133 string valString;
1134 if (!attributeValue(line, attribute, valString)) return false;
1135 istringstream valStream(valString);
1136 if ( !(valStream >> val) ) {
1137 loggerPtr->ERROR_MSG("failed to store attribute "
1138 + attribute + " " + valString);
1139 return false;
1140 } else return true;
1141 }
1142 bool addBranching(string line, unordered_map< pair<int, int>,
1143 vector<EWBranching> > & branchings, unordered_map< pair<int, int>,
1144 vector<pair<int, int> > > & clusterings, double headroom, bool decay);
1145 bool addParticle(int idIn, int polIn, bool isRes);
1146
1147
1148
1149
1150 double q2minSav;
1151
1152
1153 AlphaEM al;
1154
1155
1156 EWSystem ewSystem;
1157
1158
1159 double q2Trial;
1160
1161
1162 bool isLoaded, doEW, doFFbranchings, doIIbranchings, doRFbranchings,
1163 doBosonInterference;
1164
1165
1166 int nFlavZeroMass;
1167
1168
1169 double headroomFinal, headroomInitial;
1170
1171 };
1172
1173
1174
1175
1176
1177 class VinciaEWVetoHook : public UserHooks {
1178
1179 public:
1180
1181
1182 VinciaEWVetoHook() = default;
1183
1184
1185 bool canVetoISREmission() override {return mayVeto;}
1186 bool canVetoFSREmission() override {return mayVeto;}
1187
1188
1189 bool doVetoISREmission(int sizeOld, const Event &event, int iSys) override;
1190 bool doVetoFSREmission(int sizeOld, const Event &event, int iSys,
1191 bool inResonance = false) override;
1192
1193
1194 void init(shared_ptr<VinciaEW> ewShowerPtrIn);
1195
1196 private:
1197
1198
1199 bool doVetoEmission(int sizeOld, const Event &event, int iSys);
1200
1201 double findQCDScale(int sizeOld, const Event& event, int iSys);
1202
1203 double findEWScale(int sizeOld, const Event& event, int iSys);
1204
1205 double ktMeasure(const Event& event, int indexi, int indexj, double mI2);
1206
1207 double findktQCD(const Event& event, int indexi, int indexj);
1208
1209 double findktEW(const Event& event, int indexi, int indexj);
1210
1211 bool setLastFSREmission(int sizeOld, const Event& event);
1212
1213 bool setLastISREmission(int sizeOld, const Event& event);
1214
1215
1216
1217
1218 int verbose;
1219
1220
1221 bool mayVeto{true};
1222
1223
1224 double deltaR;
1225
1226
1227 double q2EW;
1228
1229
1230 bool lastIsQCD;
1231 double lastkT2;
1232
1233
1234 shared_ptr<VinciaEW> ewShowerPtr{};
1235
1236 };
1237
1238
1239
1240
1241 }
1242
1243 #endif