File indexing completed on 2025-01-18 10:06:34
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{}, *beamBPtr{};
0948 Info* infoPtr{};
0949 PartonSystems* partonSystemsPtr{};
0950 Rndm* rndmPtr{};
0951 Settings* settingsPtr{};
0952 Logger* loggerPtr{};
0953 VinciaCommon* vinComPtr{};
0954 AlphaEM* al{};
0955
0956
0957 vector<EWAntennaFF> antVecFinal;
0958 vector<EWAntennaII> antVecInitial;
0959 vector<EWAntennaFFres> antVecRes;
0960
0961
0962 EWAntenna* antTrial{};
0963 double q2Trial;
0964 bool lastWasSplitSav, lastWasDecSav, lastWasInitialSav, lastWasBelowCut;
0965 int ISav, KSav;
0966
0967
0968 unordered_map<pair<int, int>, vector<EWBranching> >
0969 *brMapFinal{}, *brMapInitial{}, *brMapResonance{};
0970
0971
0972 unordered_map<pair<int, int>, vector<pair<int, int> > >
0973 *cluMapFinal{}, *cluMapInitial{};
0974
0975
0976 AmpCalculator* ampCalcPtr{};
0977
0978
0979 bool isInit, doVetoHardEmissions;
0980 int verbose;
0981 double vetoHardEmissionsDeltaR2;
0982
0983 };
0984
0985
0986
0987
0988
0989 class VinciaEW : public VinciaModule {
0990
0991 public:
0992
0993
0994 VinciaEW(): isLoaded{false} {;}
0995
0996
0997 void initPtr(Info* infoPtrIn, VinciaCommon* vinComPtrIn) override {
0998 infoPtr = infoPtrIn;
0999 particleDataPtr = infoPtr->particleDataPtr;
1000 loggerPtr = infoPtr->loggerPtr;
1001 partonSystemsPtr = infoPtr->partonSystemsPtr;
1002 rndmPtr = infoPtr->rndmPtr;
1003 settingsPtr = infoPtr->settingsPtr;
1004 vinComPtr = vinComPtrIn;
1005 ampCalc.initPtr(infoPtr, &al, &vinComPtr->alphaStrong);
1006 isInitPtr = true;}
1007
1008
1009 void init(BeamParticle* beamAPtrIn = 0, BeamParticle* beamBPtrIn = 0)
1010 override;
1011
1012
1013 bool polarise(vector<Particle> &state) override {
1014 if (isLoaded) return ampCalc.polarise(state);
1015 else return false;}
1016
1017
1018
1019 bool prepare(int iSysIn, Event &event, bool isBelowHadIn=false) override;
1020
1021
1022 void update(Event &event, int iSysIn) override {
1023 if (verbose >= VinciaConstants::DEBUG)
1024 printOut(__METHOD_NAME__, "begin", VinciaConstants::DASHLEN);
1025 if (iSysIn != ewSystem.system()) return;
1026 else ewSystem.buildSystem(event);
1027 if (verbose >= VinciaConstants::DEBUG)
1028 printOut(__METHOD_NAME__, "end", VinciaConstants::DASHLEN);}
1029
1030
1031 void setVerbose(int verboseIn) override {
1032 verbose = verboseIn; ampCalc.setVerbose(verboseIn);}
1033
1034
1035 void load() override;
1036
1037
1038 double q2Next(Event&, double q2Start,double q2End) override;
1039
1040
1041 bool lastIsSplitting() override { return ewSystem.lastWasSplitToFermions(); }
1042 bool lastIsInitial() override { return ewSystem.lastWasInitial(); }
1043 bool lastIsResonanceDecay() override {
1044 return ewSystem.lastWasResonanceDecay(); }
1045
1046
1047 bool acceptTrial(Event& event) override {
1048 if (verbose >= VinciaConstants::DEBUG)
1049 printOut(__METHOD_NAME__, "begin", VinciaConstants::DASHLEN);
1050 bool success = false;
1051 if (ewSystem.hasTrial()) success = ewSystem.acceptTrial(event);
1052 else loggerPtr->ERROR_MSG("trial doesn't exist!");
1053 if (verbose >= VinciaConstants::DEBUG)
1054 printOut(__METHOD_NAME__, "end", VinciaConstants::DASHLEN);
1055 return success;}
1056
1057
1058 void updateEvent(Event& event) override {
1059 if (verbose >= VinciaConstants::DEBUG)
1060 printOut(__METHOD_NAME__, "begin", VinciaConstants::DASHLEN);
1061 if (ewSystem.hasTrial()) ewSystem.updateEvent(event);
1062 else loggerPtr->ERROR_MSG("trial doesn't exist!");
1063 if (verbose >=VinciaConstants::DEBUG) {
1064 printOut(__METHOD_NAME__,"Event after update:"); event.list();
1065 printOut(__METHOD_NAME__, "end", VinciaConstants::DASHLEN);}}
1066
1067
1068 void updatePartonSystems(Event& event) override {
1069 if (verbose >= VinciaConstants::DEBUG)
1070 printOut(__METHOD_NAME__, "begin", VinciaConstants::DASHLEN);
1071 if (ewSystem.hasTrial()) ewSystem.updatePartonSystems(event);
1072 else loggerPtr->ERROR_MSG("trial doesn't exist!");
1073 if (verbose >= VinciaConstants::DEBUG)
1074 printOut(__METHOD_NAME__, "end", VinciaConstants::DASHLEN);}
1075
1076
1077 void clear(int) override {
1078 ewSystem.clearAntennae();
1079 ampCalc.eventWeight(1);
1080 }
1081
1082
1083 unsigned int nBranchers() override {return ewSystem.nBranchers();}
1084 unsigned int nResDec() override {return ewSystem.nResDec();}
1085
1086
1087 void printData();
1088
1089
1090 void printBranchings();
1091
1092
1093 double q2min() override {return q2minSav;}
1094 double q2minColoured() override {return q2minSav;}
1095 double eventWeight() {return ampCalc.eventWeight();}
1096
1097
1098 void q2min(double q2minSavIn) {q2minSav = q2minSavIn;}
1099
1100
1101 int sysWin() override {return ewSystem.system();}
1102
1103
1104 bool isResonance(int pid) {return ewData.isRes(pid);}
1105
1106
1107
1108
1109 unordered_map<pair<int, int>, vector<pair<int, int> > >
1110 cluMapFinal, cluMapInitial;
1111
1112
1113
1114 unordered_map<pair<int, int>, vector<EWBranching> >
1115 brMapFinal, brMapInitial, brMapResonance;
1116
1117
1118 EWParticleData ewData;
1119
1120
1121 AmpCalculator ampCalc;
1122
1123 private:
1124
1125
1126
1127 bool readFile(string file);
1128 bool readLine(string line);
1129 bool attributeValue(string line, string attribute, string& val);
1130 template <class T> bool attributeValue(
1131 string line, string attribute, T& val) {
1132 string valString;
1133 if (!attributeValue(line, attribute, valString)) return false;
1134 istringstream valStream(valString);
1135 if ( !(valStream >> val) ) {
1136 loggerPtr->ERROR_MSG("failed to store attribute "
1137 + attribute + " " + valString);
1138 return false;
1139 } else return true;
1140 }
1141 bool addBranching(string line, unordered_map< pair<int, int>,
1142 vector<EWBranching> > & branchings, unordered_map< pair<int, int>,
1143 vector<pair<int, int> > > & clusterings, double headroom, bool decay);
1144 bool addParticle(int idIn, int polIn, bool isRes);
1145
1146
1147
1148
1149 double q2minSav;
1150
1151
1152 AlphaEM al;
1153
1154
1155 EWSystem ewSystem;
1156
1157
1158 double q2Trial;
1159
1160
1161 bool isLoaded, doEW, doFFbranchings, doIIbranchings, doRFbranchings,
1162 doBosonInterference;
1163
1164
1165 int nFlavZeroMass;
1166
1167
1168 double headroomFinal, headroomInitial;
1169
1170 };
1171
1172
1173
1174
1175
1176 class VinciaEWVetoHook : public UserHooks {
1177
1178 public:
1179
1180
1181 VinciaEWVetoHook() = default;
1182
1183
1184 bool canVetoISREmission() override {return mayVeto;}
1185 bool canVetoFSREmission() override {return mayVeto;}
1186
1187
1188 bool doVetoISREmission(int sizeOld, const Event &event, int iSys) override;
1189 bool doVetoFSREmission(int sizeOld, const Event &event, int iSys,
1190 bool inResonance = false) override;
1191
1192
1193 void init(shared_ptr<VinciaEW> ewShowerPtrIn);
1194
1195 private:
1196
1197
1198 bool doVetoEmission(int sizeOld, const Event &event, int iSys);
1199
1200 double findQCDScale(int sizeOld, const Event& event, int iSys);
1201
1202 double findEWScale(int sizeOld, const Event& event, int iSys);
1203
1204 double ktMeasure(const Event& event, int indexi, int indexj, double mI2);
1205
1206 double findktQCD(const Event& event, int indexi, int indexj);
1207
1208 double findktEW(const Event& event, int indexi, int indexj);
1209
1210 bool setLastFSREmission(int sizeOld, const Event& event);
1211
1212 bool setLastISREmission(int sizeOld, const Event& event);
1213
1214
1215
1216
1217 int verbose;
1218
1219
1220 bool mayVeto{true};
1221
1222
1223 double deltaR;
1224
1225
1226 double q2EW;
1227
1228
1229 bool lastIsQCD;
1230 double lastkT2;
1231
1232
1233 shared_ptr<VinciaEW> ewShowerPtr{};
1234
1235 };
1236
1237
1238
1239
1240 }
1241
1242 #endif