[1] | 1 | // main31.cc is a part of the PYTHIA event generator. |
---|
| 2 | // Copyright (C) 2012 Richard Corke, Torbjorn Sjostrand. |
---|
| 3 | // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details. |
---|
| 4 | // Please respect the MCnet Guidelines, see GUIDELINES for details. |
---|
| 5 | |
---|
| 6 | #include "Pythia.h" |
---|
| 7 | using namespace Pythia8; |
---|
| 8 | |
---|
| 9 | //========================================================================== |
---|
| 10 | |
---|
| 11 | // Use userhooks to veto PYTHIA emissions above the POWHEG scale. |
---|
| 12 | |
---|
| 13 | class PowhegHooks : public UserHooks { |
---|
| 14 | |
---|
| 15 | public: |
---|
| 16 | |
---|
| 17 | // Constructor and destructor. |
---|
| 18 | PowhegHooks(int nFinalIn, int vetoModeIn, int vetoCountIn, |
---|
| 19 | int pThardModeIn, int pTemtModeIn, int emittedModeIn, |
---|
| 20 | int pTdefModeIn, int MPIvetoModeIn) : nFinal(nFinalIn), |
---|
| 21 | vetoMode(vetoModeIn), vetoCount(vetoCountIn), |
---|
| 22 | pThardMode(pThardModeIn), pTemtMode(pTemtModeIn), |
---|
| 23 | emittedMode(emittedModeIn), pTdefMode(pTdefModeIn), |
---|
| 24 | MPIvetoMode(MPIvetoModeIn) {}; |
---|
| 25 | ~PowhegHooks() {} |
---|
| 26 | |
---|
| 27 | //-------------------------------------------------------------------------- |
---|
| 28 | |
---|
| 29 | // Routines to calculate the pT (according to pTdefMode) in a splitting: |
---|
| 30 | // ISR: i (radiator after) -> j (emitted after) k (radiator before) |
---|
| 31 | // FSR: i (radiator before) -> j (emitted after) k (radiator after) |
---|
| 32 | // For the Pythia pT definition, a recoiler (after) must be specified. |
---|
| 33 | |
---|
| 34 | // Compute the Pythia pT separation. Based on pTLund function in History.cc |
---|
| 35 | double pTpythia(const Event &e, int RadAfterBranch, int EmtAfterBranch, |
---|
| 36 | int RecAfterBranch, bool FSR) { |
---|
| 37 | |
---|
| 38 | // Convenient shorthands for later |
---|
| 39 | Vec4 radVec = e[RadAfterBranch].p(); |
---|
| 40 | Vec4 emtVec = e[EmtAfterBranch].p(); |
---|
| 41 | Vec4 recVec = e[RecAfterBranch].p(); |
---|
| 42 | int radID = e[RadAfterBranch].id(); |
---|
| 43 | |
---|
| 44 | // Calculate virtuality of splitting |
---|
| 45 | double sign = (FSR) ? 1. : -1.; |
---|
| 46 | Vec4 Q(radVec + sign * emtVec); |
---|
| 47 | double Qsq = sign * Q.m2Calc(); |
---|
| 48 | |
---|
| 49 | // Mass term of radiator |
---|
| 50 | double m2Rad = (abs(radID) >= 4 && abs(radID) < 7) ? |
---|
| 51 | pow2(particleDataPtr->m0(radID)) : 0.; |
---|
| 52 | |
---|
| 53 | // z values for FSR and ISR |
---|
| 54 | double z, pTnow; |
---|
| 55 | if (FSR) { |
---|
| 56 | // Construct 2 -> 3 variables |
---|
| 57 | Vec4 sum = radVec + recVec + emtVec; |
---|
| 58 | double m2Dip = sum.m2Calc(); |
---|
| 59 | double x1 = 2. * (sum * radVec) / m2Dip; |
---|
| 60 | double x3 = 2. * (sum * emtVec) / m2Dip; |
---|
| 61 | z = x1 / (x1 + x3); |
---|
| 62 | pTnow = z * (1. - z); |
---|
| 63 | |
---|
| 64 | } else { |
---|
| 65 | // Construct dipoles before/after splitting |
---|
| 66 | Vec4 qBR(radVec - emtVec + recVec); |
---|
| 67 | Vec4 qAR(radVec + recVec); |
---|
| 68 | z = qBR.m2Calc() / qAR.m2Calc(); |
---|
| 69 | pTnow = (1. - z); |
---|
| 70 | } |
---|
| 71 | |
---|
| 72 | // Virtuality with correct sign |
---|
| 73 | pTnow *= (Qsq - sign * m2Rad); |
---|
| 74 | |
---|
| 75 | // Can get negative pT for massive splittings |
---|
| 76 | if (pTnow < 0.) { |
---|
| 77 | cout << "Warning: pTpythia was negative" << endl; |
---|
| 78 | return -1.; |
---|
| 79 | } |
---|
| 80 | |
---|
| 81 | #ifdef DBGOUTPUT |
---|
| 82 | cout << "pTpythia: rad = " << RadAfterBranch << ", emt = " |
---|
| 83 | << EmtAfterBranch << ", rec = " << RecAfterBranch |
---|
| 84 | << ", pTnow = " << sqrt(pTnow) << endl; |
---|
| 85 | #endif |
---|
| 86 | |
---|
| 87 | // Return pT |
---|
| 88 | return sqrt(pTnow); |
---|
| 89 | } |
---|
| 90 | |
---|
| 91 | // Compute the POWHEG pT separation between i and j |
---|
| 92 | double pTpowheg(const Event &e, int i, int j, bool FSR) { |
---|
| 93 | |
---|
| 94 | // pT value for FSR and ISR |
---|
| 95 | double pTnow = 0.; |
---|
| 96 | if (FSR) { |
---|
| 97 | // POWHEG d_ij (in CM frame). Note that the incoming beams have not |
---|
| 98 | // been updated in the parton systems pointer yet (i.e. prior to any |
---|
| 99 | // potential recoil). |
---|
| 100 | int iInA = partonSystemsPtr->getInA(0); |
---|
| 101 | int iInB = partonSystemsPtr->getInB(0); |
---|
| 102 | double betaZ = - ( e[iInA].pz() + e[iInB].pz() ) / |
---|
| 103 | ( e[iInA].e() + e[iInB].e() ); |
---|
| 104 | Vec4 iVecBst(e[i].p()), jVecBst(e[j].p()); |
---|
| 105 | iVecBst.bst(0., 0., betaZ); |
---|
| 106 | jVecBst.bst(0., 0., betaZ); |
---|
| 107 | pTnow = sqrt( (iVecBst + jVecBst).m2Calc() * |
---|
| 108 | iVecBst.e() * jVecBst.e() / |
---|
| 109 | pow2(iVecBst.e() + jVecBst.e()) ); |
---|
| 110 | |
---|
| 111 | } else { |
---|
| 112 | // POWHEG pT_ISR is just kinematic pT |
---|
| 113 | pTnow = e[j].pT(); |
---|
| 114 | } |
---|
| 115 | |
---|
| 116 | // Check result |
---|
| 117 | if (pTnow < 0.) { |
---|
| 118 | cout << "Warning: pTpowheg was negative" << endl; |
---|
| 119 | return -1.; |
---|
| 120 | } |
---|
| 121 | |
---|
| 122 | #ifdef DBGOUTPUT |
---|
| 123 | cout << "pTpowheg: i = " << i << ", j = " << j |
---|
| 124 | << ", pTnow = " << pTnow << endl; |
---|
| 125 | #endif |
---|
| 126 | |
---|
| 127 | return pTnow; |
---|
| 128 | } |
---|
| 129 | |
---|
| 130 | // Calculate pT for a splitting based on pTdefMode. |
---|
| 131 | // If j is -1, all final-state partons are tried. |
---|
| 132 | // If i, k, r and xSR are -1, then all incoming and outgoing |
---|
| 133 | // partons are tried. |
---|
| 134 | // xSR set to 0 means ISR, while xSR set to 1 means FSR |
---|
| 135 | double pTcalc(const Event &e, int i, int j, int k, int r, int xSRin) { |
---|
| 136 | |
---|
| 137 | // Loop over ISR and FSR if necessary |
---|
| 138 | double pTemt = -1., pTnow; |
---|
| 139 | int xSR1 = (xSRin == -1) ? 0 : xSRin; |
---|
| 140 | int xSR2 = (xSRin == -1) ? 2 : xSRin + 1; |
---|
| 141 | for (int xSR = xSR1; xSR < xSR2; xSR++) { |
---|
| 142 | // FSR flag |
---|
| 143 | bool FSR = (xSR == 0) ? false : true; |
---|
| 144 | |
---|
| 145 | // If all necessary arguments have been given, then directly calculate. |
---|
| 146 | // POWHEG ISR and FSR, need i and j. |
---|
| 147 | if ((pTdefMode == 0 || pTdefMode == 1) && i > 0 && j > 0) { |
---|
| 148 | pTemt = pTpowheg(e, i, j, (pTdefMode == 0) ? false : FSR); |
---|
| 149 | |
---|
| 150 | // Pythia ISR, need i, j and r. |
---|
| 151 | } else if (!FSR && pTdefMode == 2 && i > 0 && j > 0 && r > 0) { |
---|
| 152 | pTemt = pTpythia(e, i, j, r, FSR); |
---|
| 153 | |
---|
| 154 | // Pythia FSR, need k, j and r. |
---|
| 155 | } else if (FSR && pTdefMode == 2 && j > 0 && k > 0 && r > 0) { |
---|
| 156 | pTemt = pTpythia(e, k, j, r, FSR); |
---|
| 157 | |
---|
| 158 | // Otherwise need to try all possible combintations. |
---|
| 159 | } else { |
---|
| 160 | // Start by finding incoming legs to the hard system after |
---|
| 161 | // branching (radiator after branching, i for ISR). |
---|
| 162 | // Use partonSystemsPtr to find incoming just prior to the |
---|
| 163 | // branching and track mothers. |
---|
| 164 | int iInA = partonSystemsPtr->getInA(0); |
---|
| 165 | int iInB = partonSystemsPtr->getInB(0); |
---|
| 166 | while (e[iInA].mother1() != 1) { iInA = e[iInA].mother1(); } |
---|
| 167 | while (e[iInB].mother1() != 2) { iInB = e[iInB].mother1(); } |
---|
| 168 | |
---|
| 169 | // If we do not have j, then try all final-state partons |
---|
| 170 | int jNow = (j > 0) ? j : 0; |
---|
| 171 | int jMax = (j > 0) ? j + 1 : e.size(); |
---|
| 172 | for (; jNow < jMax; jNow++) { |
---|
| 173 | |
---|
| 174 | // Final-state and coloured jNow only |
---|
| 175 | if (!e[jNow].isFinal() || e[jNow].colType() == 0) continue; |
---|
| 176 | |
---|
| 177 | // POWHEG |
---|
| 178 | if (pTdefMode == 0 || pTdefMode == 1) { |
---|
| 179 | |
---|
| 180 | // ISR - only done once as just kinematical pT |
---|
| 181 | if (!FSR) { |
---|
| 182 | pTnow = pTpowheg(e, iInA, jNow, (pTdefMode == 0) ? false : FSR); |
---|
| 183 | if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow); |
---|
| 184 | |
---|
| 185 | // FSR - try all outgoing partons from system before branching |
---|
| 186 | // as i. Note that for the hard system, there is no |
---|
| 187 | // "before branching" information. |
---|
| 188 | } else { |
---|
| 189 | |
---|
| 190 | int outSize = partonSystemsPtr->sizeOut(0); |
---|
| 191 | for (int iMem = 0; iMem < outSize; iMem++) { |
---|
| 192 | int iNow = partonSystemsPtr->getOut(0, iMem); |
---|
| 193 | |
---|
| 194 | // Coloured only, i != jNow and no carbon copies |
---|
| 195 | if (iNow == jNow || e[iNow].colType() == 0) continue; |
---|
| 196 | if (jNow == e[iNow].daughter1() |
---|
| 197 | && jNow == e[iNow].daughter2()) continue; |
---|
| 198 | |
---|
| 199 | pTnow = pTpowheg(e, iNow, jNow, (pTdefMode == 0) |
---|
| 200 | ? false : FSR); |
---|
| 201 | if (pTnow > 0.) pTemt = (pTemt < 0) |
---|
| 202 | ? pTnow : min(pTemt, pTnow); |
---|
| 203 | } // for (iMem) |
---|
| 204 | |
---|
| 205 | } // if (!FSR) |
---|
| 206 | |
---|
| 207 | // Pythia |
---|
| 208 | } else if (pTdefMode == 2) { |
---|
| 209 | |
---|
| 210 | // ISR - other incoming as recoiler |
---|
| 211 | if (!FSR) { |
---|
| 212 | pTnow = pTpythia(e, iInA, jNow, iInB, FSR); |
---|
| 213 | if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow); |
---|
| 214 | pTnow = pTpythia(e, iInB, jNow, iInA, FSR); |
---|
| 215 | if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow); |
---|
| 216 | |
---|
| 217 | // FSR - try all final-state coloured partons as radiator |
---|
| 218 | // after emission (k). |
---|
| 219 | } else { |
---|
| 220 | for (int kNow = 0; kNow < e.size(); kNow++) { |
---|
| 221 | if (kNow == jNow || !e[kNow].isFinal() || |
---|
| 222 | e[kNow].colType() == 0) continue; |
---|
| 223 | |
---|
| 224 | // For this kNow, need to have a recoiler. |
---|
| 225 | // Try two incoming. |
---|
| 226 | pTnow = pTpythia(e, kNow, jNow, iInA, FSR); |
---|
| 227 | if (pTnow > 0.) pTemt = (pTemt < 0) |
---|
| 228 | ? pTnow : min(pTemt, pTnow); |
---|
| 229 | pTnow = pTpythia(e, kNow, jNow, iInB, FSR); |
---|
| 230 | if (pTnow > 0.) pTemt = (pTemt < 0) |
---|
| 231 | ? pTnow : min(pTemt, pTnow); |
---|
| 232 | |
---|
| 233 | // Try all other outgoing. |
---|
| 234 | for (int rNow = 0; rNow < e.size(); rNow++) { |
---|
| 235 | if (rNow == kNow || rNow == jNow || |
---|
| 236 | !e[rNow].isFinal() || e[rNow].colType() == 0) continue; |
---|
| 237 | pTnow = pTpythia(e, kNow, jNow, rNow, FSR); |
---|
| 238 | if (pTnow > 0.) pTemt = (pTemt < 0) |
---|
| 239 | ? pTnow : min(pTemt, pTnow); |
---|
| 240 | } // for (rNow) |
---|
| 241 | |
---|
| 242 | } // for (kNow) |
---|
| 243 | } // if (!FSR) |
---|
| 244 | } // if (pTdefMode) |
---|
| 245 | } // for (j) |
---|
| 246 | } |
---|
| 247 | } // for (xSR) |
---|
| 248 | |
---|
| 249 | #ifdef DBGOUTPUT |
---|
| 250 | cout << "pTcalc: i = " << i << ", j = " << j << ", k = " << k |
---|
| 251 | << ", r = " << r << ", xSR = " << xSRin |
---|
| 252 | << ", pTemt = " << pTemt << endl; |
---|
| 253 | #endif |
---|
| 254 | |
---|
| 255 | return pTemt; |
---|
| 256 | } |
---|
| 257 | |
---|
| 258 | //-------------------------------------------------------------------------- |
---|
| 259 | |
---|
| 260 | // Extraction of pThard based on the incoming event. |
---|
| 261 | // Assume that all the final-state particles are in a continuous block |
---|
| 262 | // at the end of the event and the final entry is the POWHEG emission. |
---|
| 263 | // If there is no POWHEG emission, then pThard is set to Qfac. |
---|
| 264 | |
---|
| 265 | bool canVetoMPIStep() { return true; } |
---|
| 266 | int numberVetoMPIStep() { return 1; } |
---|
| 267 | bool doVetoMPIStep(int nMPI, const Event &e) { |
---|
| 268 | // Extra check on nMPI |
---|
| 269 | if (nMPI > 1) return false; |
---|
| 270 | |
---|
| 271 | // Find if there is a POWHEG emission. Go backwards through the |
---|
| 272 | // event record until there is a non-final particle. Also sum pT and |
---|
| 273 | // find pT_1 for possible MPI vetoing |
---|
| 274 | int count = 0; |
---|
| 275 | double pT1 = 0., pTsum = 0.; |
---|
| 276 | for (int i = e.size() - 1; i > 0; i--) { |
---|
| 277 | if (e[i].isFinal()) { |
---|
| 278 | count++; |
---|
| 279 | pT1 = e[i].pT(); |
---|
| 280 | pTsum += e[i].pT(); |
---|
| 281 | } else break; |
---|
| 282 | } |
---|
| 283 | // Extra check that we have the correct final state |
---|
| 284 | if (count != nFinal && count != nFinal + 1) { |
---|
| 285 | cout << "Error: wrong number of final state particles in event" << endl; |
---|
| 286 | exit(1); |
---|
| 287 | } |
---|
| 288 | // Flag if POWHEG radiation present and index |
---|
| 289 | bool isEmt = (count == nFinal) ? false : true; |
---|
| 290 | int iEmt = (isEmt) ? e.size() - 1 : -1; |
---|
| 291 | |
---|
| 292 | // If there is no radiation or if pThardMode is 0 then set pThard to Qfac. |
---|
| 293 | if (!isEmt || pThardMode == 0) { |
---|
| 294 | pThard = infoPtr->QFac(); |
---|
| 295 | |
---|
| 296 | // If pThardMode is 1 then the pT of the POWHEG emission is checked against |
---|
| 297 | // all other incoming and outgoing partons, with the minimal value taken |
---|
| 298 | } else if (pThardMode == 1) { |
---|
| 299 | pThard = pTcalc(e, -1, iEmt, -1, -1, -1); |
---|
| 300 | |
---|
| 301 | // If pThardMode is 2, then the pT of all final-state partons is checked |
---|
| 302 | // against all other incoming and outgoing partons, with the minimal value |
---|
| 303 | // taken |
---|
| 304 | } else if (pThardMode == 2) { |
---|
| 305 | pThard = pTcalc(e, -1, -1, -1, -1, -1); |
---|
| 306 | |
---|
| 307 | } |
---|
| 308 | |
---|
| 309 | // Find MPI veto pT if necessary |
---|
| 310 | if (MPIvetoMode == 1) { |
---|
| 311 | pTMPI = (isEmt) ? pTsum / 2. : pT1; |
---|
| 312 | } |
---|
| 313 | |
---|
| 314 | #ifdef DBGOUTPUT |
---|
| 315 | cout << "doVetoMPIStep: Qfac = " << infoPtr->QFac() |
---|
| 316 | << ", pThard = " << pThard << endl << endl; |
---|
| 317 | #endif |
---|
| 318 | |
---|
| 319 | // Initialise other variables |
---|
| 320 | accepted = false; |
---|
| 321 | nAcceptSeq = nISRveto = nFSRveto = 0; |
---|
| 322 | |
---|
| 323 | // Do not veto the event |
---|
| 324 | return false; |
---|
| 325 | } |
---|
| 326 | |
---|
| 327 | //-------------------------------------------------------------------------- |
---|
| 328 | |
---|
| 329 | // ISR veto |
---|
| 330 | |
---|
| 331 | bool canVetoISREmission() { return (vetoMode == 0) ? false : true; } |
---|
| 332 | bool doVetoISREmission(int, const Event &e, int iSys) { |
---|
| 333 | // Must be radiation from the hard system |
---|
| 334 | if (iSys != 0) return false; |
---|
| 335 | |
---|
| 336 | // If we already have accepted 'vetoCount' emissions in a row, do nothing |
---|
| 337 | if (vetoMode == 1 && nAcceptSeq >= vetoCount) return false; |
---|
| 338 | |
---|
| 339 | // Pythia radiator after, emitted and recoiler after. |
---|
| 340 | int iRadAft = -1, iEmt = -1, iRecAft = -1; |
---|
| 341 | for (int i = e.size() - 1; i > 0; i--) { |
---|
| 342 | if (iRadAft == -1 && e[i].status() == -41) iRadAft = i; |
---|
| 343 | else if (iEmt == -1 && e[i].status() == 43) iEmt = i; |
---|
| 344 | else if (iRecAft == -1 && e[i].status() == -42) iRecAft = i; |
---|
| 345 | if (iRadAft != -1 && iEmt != -1 && iRecAft != -1) break; |
---|
| 346 | } |
---|
| 347 | if (iRadAft == -1 || iEmt == -1 || iRecAft == -1) { |
---|
| 348 | e.list(); |
---|
| 349 | cout << "Error: couldn't find Pythia ISR emission" << endl; |
---|
| 350 | exit(1); |
---|
| 351 | } |
---|
| 352 | |
---|
| 353 | // pTemtMode == 0: pT of emitted w.r.t. radiator |
---|
| 354 | // pTemtMode == 1: min(pT of emitted w.r.t. all incoming/outgoing) |
---|
| 355 | // pTemtMode == 2: min(pT of all outgoing w.r.t. all incoming/outgoing) |
---|
| 356 | int xSR = (pTemtMode == 0) ? 0 : -1; |
---|
| 357 | int i = (pTemtMode == 0) ? iRadAft : -1; |
---|
| 358 | int j = (pTemtMode != 2) ? iEmt : -1; |
---|
| 359 | int k = -1; |
---|
| 360 | int r = (pTemtMode == 0) ? iRecAft : -1; |
---|
| 361 | double pTemt = pTcalc(e, i, j, k, r, xSR); |
---|
| 362 | |
---|
| 363 | #ifdef DBGOUTPUT |
---|
| 364 | cout << "doVetoISREmission: pTemt = " << pTemt << endl << endl; |
---|
| 365 | #endif |
---|
| 366 | |
---|
| 367 | // Veto if pTemt > pThard |
---|
| 368 | if (pTemt > pThard) { |
---|
| 369 | nAcceptSeq = 0; |
---|
| 370 | nISRveto++; |
---|
| 371 | return true; |
---|
| 372 | } |
---|
| 373 | |
---|
| 374 | // Else mark that an emission has been accepted and continue |
---|
| 375 | nAcceptSeq++; |
---|
| 376 | accepted = true; |
---|
| 377 | return false; |
---|
| 378 | } |
---|
| 379 | |
---|
| 380 | //-------------------------------------------------------------------------- |
---|
| 381 | |
---|
| 382 | // FSR veto |
---|
| 383 | |
---|
| 384 | bool canVetoFSREmission() { return (vetoMode == 0) ? false : true; } |
---|
| 385 | bool doVetoFSREmission(int, const Event &e, int iSys, bool) { |
---|
| 386 | // Must be radiation from the hard system |
---|
| 387 | if (iSys != 0) return false; |
---|
| 388 | |
---|
| 389 | // If we already have accepted 'vetoCount' emissions in a row, do nothing |
---|
| 390 | if (vetoMode == 1 && nAcceptSeq >= vetoCount) return false; |
---|
| 391 | |
---|
| 392 | // Pythia radiator (before and after), emitted and recoiler (after) |
---|
| 393 | int iRecAft = e.size() - 1; |
---|
| 394 | int iEmt = e.size() - 2; |
---|
| 395 | int iRadAft = e.size() - 3; |
---|
| 396 | int iRadBef = e[iEmt].mother1(); |
---|
| 397 | if ( (e[iRecAft].status() != 52 && e[iRecAft].status() != -53) || |
---|
| 398 | e[iEmt].status() != 51 || e[iRadAft].status() != 51) { |
---|
| 399 | e.list(); |
---|
| 400 | cout << "Error: couldn't find Pythia FSR emission" << endl; |
---|
| 401 | exit(1); |
---|
| 402 | } |
---|
| 403 | |
---|
| 404 | // Behaviour based on pTemtMode: |
---|
| 405 | // 0 - pT of emitted w.r.t. radiator before |
---|
| 406 | // 1 - min(pT of emitted w.r.t. all incoming/outgoing) |
---|
| 407 | // 2 - min(pT of all outgoing w.r.t. all incoming/outgoing) |
---|
| 408 | int xSR = (pTemtMode == 0) ? 1 : -1; |
---|
| 409 | int i = (pTemtMode == 0) ? iRadBef : -1; |
---|
| 410 | int k = (pTemtMode == 0) ? iRadAft : -1; |
---|
| 411 | int r = (pTemtMode == 0) ? iRecAft : -1; |
---|
| 412 | |
---|
| 413 | // When pTemtMode is 0 or 1, iEmt has been selected |
---|
| 414 | double pTemt; |
---|
| 415 | if (pTemtMode == 0 || pTemtMode == 1) { |
---|
| 416 | // Which parton is emitted, based on emittedMode: |
---|
| 417 | // 0 - Pythia definition of emitted |
---|
| 418 | // 1 - Pythia definition of radiated after emission |
---|
| 419 | // 2 - Random selection of emitted or radiated after emission |
---|
| 420 | // 3 - Try both emitted and radiated after emission |
---|
| 421 | int j = iRadAft; |
---|
| 422 | if (emittedMode == 0 || (emittedMode == 2 && rndmPtr->flat() < 0.5)) j++; |
---|
| 423 | |
---|
| 424 | for (int jLoop = 0; jLoop < 2; jLoop++) { |
---|
| 425 | if (jLoop == 0) pTemt = pTcalc(e, i, j, k, r, xSR); |
---|
| 426 | else if (jLoop == 1) pTemt = min(pTemt, pTcalc(e, i, j, k, r, xSR)); |
---|
| 427 | |
---|
| 428 | // For emittedMode == 3, have tried iRadAft, now try iEmt |
---|
| 429 | if (emittedMode != 3) break; |
---|
| 430 | if (k != -1) swap(j, k); else j = iEmt; |
---|
| 431 | } |
---|
| 432 | |
---|
| 433 | // If pTemtMode is 2, then try all final-state partons as emitted |
---|
| 434 | } else if (pTemtMode == 2) { |
---|
| 435 | pTemt = pTcalc(e, i, -1, k, r, xSR); |
---|
| 436 | |
---|
| 437 | } |
---|
| 438 | |
---|
| 439 | #ifdef DBGOUTPUT |
---|
| 440 | cout << "doVetoFSREmission: pTemt = " << pTemt << endl << endl; |
---|
| 441 | #endif |
---|
| 442 | |
---|
| 443 | // Veto if pTemt > pThard |
---|
| 444 | if (pTemt > pThard) { |
---|
| 445 | nAcceptSeq = 0; |
---|
| 446 | nFSRveto++; |
---|
| 447 | return true; |
---|
| 448 | } |
---|
| 449 | |
---|
| 450 | // Else mark that an emission has been accepted and continue |
---|
| 451 | nAcceptSeq++; |
---|
| 452 | accepted = true; |
---|
| 453 | return false; |
---|
| 454 | } |
---|
| 455 | |
---|
| 456 | //-------------------------------------------------------------------------- |
---|
| 457 | |
---|
| 458 | // MPI veto |
---|
| 459 | |
---|
| 460 | bool canVetoMPIEmission() { return (MPIvetoMode == 0) ? false : true; } |
---|
| 461 | bool doVetoMPIEmission(int, const Event &e) { |
---|
| 462 | if (MPIvetoMode == 1) { |
---|
| 463 | if (e[e.size() - 1].pT() > pTMPI) { |
---|
| 464 | #ifdef DBGOUTPUT |
---|
| 465 | cout << "doVetoMPIEmission: pTnow = " << e[e.size() - 1].pT() |
---|
| 466 | << ", pTMPI = " << pTMPI << endl << endl; |
---|
| 467 | #endif |
---|
| 468 | return true; |
---|
| 469 | } |
---|
| 470 | } |
---|
| 471 | return false; |
---|
| 472 | } |
---|
| 473 | |
---|
| 474 | //-------------------------------------------------------------------------- |
---|
| 475 | |
---|
| 476 | // Functions to return information |
---|
| 477 | |
---|
| 478 | int getNISRveto() { return nISRveto; } |
---|
| 479 | int getNFSRveto() { return nFSRveto; } |
---|
| 480 | |
---|
| 481 | private: |
---|
| 482 | int nFinal, vetoMode, vetoCount, pThardMode, pTemtMode, |
---|
| 483 | emittedMode, pTdefMode, MPIvetoMode; |
---|
| 484 | double pThard, pTMPI; |
---|
| 485 | bool accepted; |
---|
| 486 | // The number of accepted emissions (in a row) |
---|
| 487 | int nAcceptSeq; |
---|
| 488 | // Statistics on vetos |
---|
| 489 | unsigned long int nISRveto, nFSRveto; |
---|
| 490 | |
---|
| 491 | }; |
---|
| 492 | |
---|
| 493 | //========================================================================== |
---|
| 494 | |
---|
| 495 | int main(int, char **) { |
---|
| 496 | |
---|
| 497 | // Generator |
---|
| 498 | Pythia pythia; |
---|
| 499 | |
---|
| 500 | // Add further settings that can be set in the configuration file |
---|
| 501 | pythia.settings.addMode("POWHEG:nFinal", 2, true, false, 1, 0); |
---|
| 502 | pythia.settings.addMode("POWHEG:veto", 0, true, true, 0, 2); |
---|
| 503 | pythia.settings.addMode("POWHEG:vetoCount", 3, true, false, 0, 0); |
---|
| 504 | pythia.settings.addMode("POWHEG:pThard", 0, true, true, 0, 2); |
---|
| 505 | pythia.settings.addMode("POWHEG:pTemt", 0, true, true, 0, 2); |
---|
| 506 | pythia.settings.addMode("POWHEG:emitted", 0, true, true, 0, 3); |
---|
| 507 | pythia.settings.addMode("POWHEG:pTdef", 0, true, true, 0, 2); |
---|
| 508 | pythia.settings.addMode("POWHEG:MPIveto", 0, true, true, 0, 1); |
---|
| 509 | |
---|
| 510 | // Load configuration file |
---|
| 511 | pythia.readFile("main31.cmnd"); |
---|
| 512 | |
---|
| 513 | // Read in main settings |
---|
| 514 | int nEvent = pythia.settings.mode("Main:numberOfEvents"); |
---|
| 515 | int nError = pythia.settings.mode("Main:timesAllowErrors"); |
---|
| 516 | // Read in POWHEG settings |
---|
| 517 | int nFinal = pythia.settings.mode("POWHEG:nFinal"); |
---|
| 518 | int vetoMode = pythia.settings.mode("POWHEG:veto"); |
---|
| 519 | int vetoCount = pythia.settings.mode("POWHEG:vetoCount"); |
---|
| 520 | int pThardMode = pythia.settings.mode("POWHEG:pThard"); |
---|
| 521 | int pTemtMode = pythia.settings.mode("POWHEG:pTemt"); |
---|
| 522 | int emittedMode = pythia.settings.mode("POWHEG:emitted"); |
---|
| 523 | int pTdefMode = pythia.settings.mode("POWHEG:pTdef"); |
---|
| 524 | int MPIvetoMode = pythia.settings.mode("POWHEG:MPIveto"); |
---|
| 525 | bool loadHooks = (vetoMode > 0 || MPIvetoMode > 0); |
---|
| 526 | |
---|
| 527 | // Add in user hooks for shower vetoing |
---|
| 528 | PowhegHooks *powhegHooks = NULL; |
---|
| 529 | if (loadHooks) { |
---|
| 530 | |
---|
| 531 | // Set ISR and FSR to start at the kinematical limit |
---|
| 532 | if (vetoMode > 0) { |
---|
| 533 | pythia.readString("SpaceShower:pTmaxMatch = 2"); |
---|
| 534 | pythia.readString("TimeShower:pTmaxMatch = 2"); |
---|
| 535 | } |
---|
| 536 | |
---|
| 537 | // Set MPI to start at the kinematical limit |
---|
| 538 | if (MPIvetoMode > 0) { |
---|
| 539 | pythia.readString("MultipartonInteractions:pTmaxMatch = 2"); |
---|
| 540 | } |
---|
| 541 | |
---|
| 542 | powhegHooks = new PowhegHooks(nFinal, vetoMode, vetoCount, |
---|
| 543 | pThardMode, pTemtMode, emittedMode, pTdefMode, MPIvetoMode); |
---|
| 544 | pythia.setUserHooksPtr((UserHooks *) powhegHooks); |
---|
| 545 | } |
---|
| 546 | |
---|
| 547 | // Initialise and list settings |
---|
| 548 | pythia.init(); |
---|
| 549 | |
---|
| 550 | // Counters for number of ISR/FSR emissions vetoed |
---|
| 551 | unsigned long int nISRveto = 0, nFSRveto = 0; |
---|
| 552 | |
---|
| 553 | // Begin event loop; generate until nEvent events are processed |
---|
| 554 | // or end of LHEF file |
---|
| 555 | int iEvent = 0, iError = 0; |
---|
| 556 | while (true) { |
---|
| 557 | |
---|
| 558 | // Generate the next event |
---|
| 559 | if (!pythia.next()) { |
---|
| 560 | |
---|
| 561 | // If failure because reached end of file then exit event loop |
---|
| 562 | if (pythia.info.atEndOfFile()) break; |
---|
| 563 | |
---|
| 564 | // Otherwise count event failure and continue/exit as necessary |
---|
| 565 | cout << "Warning: event " << iEvent << " failed" << endl; |
---|
| 566 | if (++iError == nError) { |
---|
| 567 | cout << "Error: too many event failures.. exiting" << endl; |
---|
| 568 | break; |
---|
| 569 | } |
---|
| 570 | |
---|
| 571 | continue; |
---|
| 572 | } |
---|
| 573 | |
---|
| 574 | /* |
---|
| 575 | * Process dependent checks and analysis may be inserted here |
---|
| 576 | */ |
---|
| 577 | |
---|
| 578 | // Update ISR/FSR veto counters |
---|
| 579 | if (loadHooks) { |
---|
| 580 | nISRveto += powhegHooks->getNISRveto(); |
---|
| 581 | nFSRveto += powhegHooks->getNFSRveto(); |
---|
| 582 | } |
---|
| 583 | |
---|
| 584 | // If nEvent is set, check and exit loop if necessary |
---|
| 585 | ++iEvent; |
---|
| 586 | if (nEvent != 0 && iEvent == nEvent) break; |
---|
| 587 | |
---|
| 588 | } // End of event loop. |
---|
| 589 | |
---|
| 590 | // Statistics, histograms and veto information |
---|
| 591 | pythia.stat(); |
---|
| 592 | cout << "Number of ISR emissions vetoed: " << nISRveto << endl; |
---|
| 593 | cout << "Number of FSR emissions vetoed: " << nFSRveto << endl; |
---|
| 594 | cout << endl; |
---|
| 595 | |
---|
| 596 | // Done. |
---|
| 597 | return 0; |
---|
| 598 | } |
---|
| 599 | |
---|