[1] | 1 | // Pythia.cc is a part of the PYTHIA event generator. |
---|
| 2 | // Copyright (C) 2012 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 | // Function definitions (not found in the header) for the Pythia class. |
---|
| 7 | |
---|
| 8 | #include "Pythia.h" |
---|
| 9 | |
---|
| 10 | // Access time information. |
---|
| 11 | #include <ctime> |
---|
| 12 | |
---|
| 13 | // Allow string and character manipulation. |
---|
| 14 | #include <cctype> |
---|
| 15 | |
---|
| 16 | namespace Pythia8 { |
---|
| 17 | |
---|
| 18 | //========================================================================== |
---|
| 19 | |
---|
| 20 | // The Pythia class. |
---|
| 21 | |
---|
| 22 | //-------------------------------------------------------------------------- |
---|
| 23 | |
---|
| 24 | // The current Pythia (sub)version number, to agree with XML version. |
---|
| 25 | const double Pythia::VERSIONNUMBERCODE = 8.170; |
---|
| 26 | |
---|
| 27 | //-------------------------------------------------------------------------- |
---|
| 28 | |
---|
| 29 | // Constants: could be changed here if desired, but normally should not. |
---|
| 30 | // These are of technical nature, as described for each. |
---|
| 31 | |
---|
| 32 | // Maximum number of tries to produce parton level from given input. |
---|
| 33 | const int Pythia::NTRY = 10; |
---|
| 34 | |
---|
| 35 | // Negative integer to denote that no subrun has been set. |
---|
| 36 | const int Pythia::SUBRUNDEFAULT = -999; |
---|
| 37 | |
---|
| 38 | //-------------------------------------------------------------------------- |
---|
| 39 | |
---|
| 40 | // Constructor. |
---|
| 41 | |
---|
| 42 | Pythia::Pythia(string xmlDir) { |
---|
| 43 | |
---|
| 44 | // Initial values for pointers to PDF's. |
---|
| 45 | useNewPdfA = false; |
---|
| 46 | useNewPdfB = false; |
---|
| 47 | useNewPdfHard = false; |
---|
| 48 | useNewPdfPomA = false; |
---|
| 49 | useNewPdfPomB = false; |
---|
| 50 | pdfAPtr = 0; |
---|
| 51 | pdfBPtr = 0; |
---|
| 52 | pdfHardAPtr = 0; |
---|
| 53 | pdfHardBPtr = 0; |
---|
| 54 | pdfPomAPtr = 0; |
---|
| 55 | pdfPomBPtr = 0; |
---|
| 56 | |
---|
| 57 | // Initial values for pointers to Les Houches Event objects. |
---|
| 58 | doLHA = false; |
---|
| 59 | useNewLHA = false; |
---|
| 60 | lhaUpPtr = 0; |
---|
| 61 | |
---|
| 62 | //Initial value for couplings pointer |
---|
| 63 | couplingsPtr = &couplings; |
---|
| 64 | |
---|
| 65 | // Initial value for pointer to external decay handler. |
---|
| 66 | decayHandlePtr = 0; |
---|
| 67 | |
---|
| 68 | // Initial value for pointer to user hooks. |
---|
| 69 | userHooksPtr = 0; |
---|
| 70 | |
---|
| 71 | // Initial value for pointer to merging hooks. |
---|
| 72 | doMerging = false; |
---|
| 73 | hasMergingHooks = false; |
---|
| 74 | hasOwnMergingHooks = false; |
---|
| 75 | mergingHooksPtr = 0; |
---|
| 76 | |
---|
| 77 | // Initial value for pointer to beam shape. |
---|
| 78 | useNewBeamShape = false; |
---|
| 79 | beamShapePtr = 0; |
---|
| 80 | |
---|
| 81 | // Initial values for pointers to timelike and spacelike showers. |
---|
| 82 | useNewTimes = false; |
---|
| 83 | useNewSpace = false; |
---|
| 84 | timesDecPtr = 0; |
---|
| 85 | timesPtr = 0; |
---|
| 86 | spacePtr = 0; |
---|
| 87 | |
---|
| 88 | // Find path to data files, i.e. xmldoc directory location. |
---|
| 89 | // Environment variable takes precedence, else use constructor input. |
---|
| 90 | xmlPath = ""; |
---|
| 91 | const char* PYTHIA8DATA = "PYTHIA8DATA"; |
---|
| 92 | char* envPath = getenv(PYTHIA8DATA); |
---|
| 93 | if (envPath != 0 && *envPath != '\0') { |
---|
| 94 | int i = 0; |
---|
| 95 | while (*(envPath+i) != '\0') xmlPath += *(envPath+(i++)); |
---|
| 96 | } |
---|
| 97 | else xmlPath = xmlDir; |
---|
| 98 | if (xmlPath[ xmlPath.length() - 1 ] != '/') xmlPath += "/"; |
---|
| 99 | |
---|
| 100 | // Read in files with all flags, modes, parms and words. |
---|
| 101 | settings.initPtr( &info); |
---|
| 102 | string initFile = xmlPath + "Index.xml"; |
---|
| 103 | isConstructed = settings.init( initFile); |
---|
| 104 | if (!isConstructed) { |
---|
| 105 | info.errorMsg("Abort from Pythia::Pythia: settings unavailable"); |
---|
| 106 | return; |
---|
| 107 | } |
---|
| 108 | |
---|
| 109 | // Check that XML version number matches code version number. |
---|
| 110 | double versionNumberXML = parm("Pythia:versionNumber"); |
---|
| 111 | isConstructed = (abs(versionNumberXML - VERSIONNUMBERCODE) < 0.0005); |
---|
| 112 | if (!isConstructed) { |
---|
| 113 | ostringstream errCode; |
---|
| 114 | errCode << fixed << setprecision(3) << ": in code " << VERSIONNUMBERCODE |
---|
| 115 | << " but in XML " << versionNumberXML; |
---|
| 116 | info.errorMsg("Abort from Pythia::Pythia: unmatched version numbers", |
---|
| 117 | errCode.str()); |
---|
| 118 | return; |
---|
| 119 | } |
---|
| 120 | |
---|
| 121 | // Read in files with all particle data. |
---|
| 122 | particleData.initPtr( &info, &settings, &rndm, couplingsPtr); |
---|
| 123 | string dataFile = xmlPath + "ParticleData.xml"; |
---|
| 124 | isConstructed = particleData.init( dataFile); |
---|
| 125 | if (!isConstructed) { |
---|
| 126 | info.errorMsg("Abort from Pythia::Pythia: particle data unavailable"); |
---|
| 127 | return; |
---|
| 128 | } |
---|
| 129 | |
---|
| 130 | // Write the Pythia banner to output. |
---|
| 131 | banner(); |
---|
| 132 | |
---|
| 133 | // Not initialized until at the end of the init() call. |
---|
| 134 | isInit = false; |
---|
| 135 | info.addCounter(0); |
---|
| 136 | |
---|
| 137 | } |
---|
| 138 | |
---|
| 139 | //-------------------------------------------------------------------------- |
---|
| 140 | |
---|
| 141 | // Destructor. |
---|
| 142 | |
---|
| 143 | Pythia::~Pythia() { |
---|
| 144 | |
---|
| 145 | // Delete the PDF's created with new. |
---|
| 146 | if (useNewPdfHard && pdfHardAPtr != pdfAPtr) delete pdfHardAPtr; |
---|
| 147 | if (useNewPdfHard && pdfHardBPtr != pdfBPtr) delete pdfHardBPtr; |
---|
| 148 | if (useNewPdfA) delete pdfAPtr; |
---|
| 149 | if (useNewPdfB) delete pdfBPtr; |
---|
| 150 | if (useNewPdfPomA) delete pdfPomAPtr; |
---|
| 151 | if (useNewPdfPomB) delete pdfPomBPtr; |
---|
| 152 | |
---|
| 153 | // Delete the Les Houches object created with new. |
---|
| 154 | if (useNewLHA) delete lhaUpPtr; |
---|
| 155 | |
---|
| 156 | // Delete the MergingHooks object created with new. |
---|
| 157 | if (hasOwnMergingHooks) delete mergingHooksPtr; |
---|
| 158 | |
---|
| 159 | // Delete the BeamShape object created with new. |
---|
| 160 | if (useNewBeamShape) delete beamShapePtr; |
---|
| 161 | |
---|
| 162 | // Delete the timelike and spacelike showers created with new. |
---|
| 163 | if (useNewTimes) delete timesPtr; |
---|
| 164 | if (useNewSpace) delete spacePtr; |
---|
| 165 | |
---|
| 166 | } |
---|
| 167 | |
---|
| 168 | //-------------------------------------------------------------------------- |
---|
| 169 | |
---|
| 170 | // Read in one update for a setting or particle data from a single line. |
---|
| 171 | |
---|
| 172 | bool Pythia::readString(string line, bool warn) { |
---|
| 173 | |
---|
| 174 | // Check that constructor worked. |
---|
| 175 | if (!isConstructed) return false; |
---|
| 176 | |
---|
| 177 | // If empty line then done. |
---|
| 178 | if (line.find_first_not_of(" \n\t\v\b\r\f\a") == string::npos) return true; |
---|
| 179 | |
---|
| 180 | // If first character is not a letter/digit, then taken to be a comment. |
---|
| 181 | int firstChar = line.find_first_not_of(" \n\t\v\b\r\f\a"); |
---|
| 182 | if (!isalnum(line[firstChar])) return true; |
---|
| 183 | |
---|
| 184 | // Send on particle data to the ParticleData database. |
---|
| 185 | if (isdigit(line[firstChar])) |
---|
| 186 | return particleData.readString(line, warn); |
---|
| 187 | |
---|
| 188 | // Everything else sent on to Settings. |
---|
| 189 | return settings.readString(line, warn); |
---|
| 190 | |
---|
| 191 | } |
---|
| 192 | |
---|
| 193 | //-------------------------------------------------------------------------- |
---|
| 194 | |
---|
| 195 | // Read in updates for settings or particle data from user-defined file. |
---|
| 196 | |
---|
| 197 | bool Pythia::readFile(string fileName, bool warn, int subrun) { |
---|
| 198 | |
---|
| 199 | // Check that constructor worked. |
---|
| 200 | if (!isConstructed) return false; |
---|
| 201 | |
---|
| 202 | // Open file for reading. |
---|
| 203 | const char* cstring = fileName.c_str(); |
---|
| 204 | ifstream is(cstring); |
---|
| 205 | if (!is.good()) { |
---|
| 206 | info.errorMsg("Error in Pythia::readFile: did not find file", fileName); |
---|
| 207 | return false; |
---|
| 208 | } |
---|
| 209 | |
---|
| 210 | // Hand over real work to next method. |
---|
| 211 | return readFile( is, warn, subrun); |
---|
| 212 | |
---|
| 213 | } |
---|
| 214 | |
---|
| 215 | //-------------------------------------------------------------------------- |
---|
| 216 | |
---|
| 217 | // Read in updates for settings or particle data |
---|
| 218 | // from user-defined stream (or file). |
---|
| 219 | |
---|
| 220 | bool Pythia::readFile(istream& is, bool warn, int subrun) { |
---|
| 221 | |
---|
| 222 | // Check that constructor worked. |
---|
| 223 | if (!isConstructed) return false; |
---|
| 224 | |
---|
| 225 | // Read in one line at a time. |
---|
| 226 | string line; |
---|
| 227 | bool accepted = true; |
---|
| 228 | int subrunNow = SUBRUNDEFAULT; |
---|
| 229 | while ( getline(is, line) ) { |
---|
| 230 | |
---|
| 231 | // Check whether entered new subrun. |
---|
| 232 | int subrunLine = readSubrun( line, warn); |
---|
| 233 | if (subrunLine >= 0) subrunNow = subrunLine; |
---|
| 234 | |
---|
| 235 | // Process the line if in correct subrun. |
---|
| 236 | if ( (subrunNow == subrun || subrunNow == SUBRUNDEFAULT) |
---|
| 237 | && !readString( line, warn) ) accepted = false; |
---|
| 238 | |
---|
| 239 | // Reached end of input file. |
---|
| 240 | }; |
---|
| 241 | return accepted; |
---|
| 242 | |
---|
| 243 | } |
---|
| 244 | |
---|
| 245 | //-------------------------------------------------------------------------- |
---|
| 246 | |
---|
| 247 | // Routine to pass in pointers to PDF's. Usage optional. |
---|
| 248 | |
---|
| 249 | bool Pythia::setPDFPtr( PDF* pdfAPtrIn, PDF* pdfBPtrIn, PDF* pdfHardAPtrIn, |
---|
| 250 | PDF* pdfHardBPtrIn, PDF* pdfPomAPtrIn, PDF* pdfPomBPtrIn) { |
---|
| 251 | |
---|
| 252 | // Delete any PDF's created in a previous init call. |
---|
| 253 | if (useNewPdfHard && pdfHardAPtr != pdfAPtr) delete pdfHardAPtr; |
---|
| 254 | if (useNewPdfHard && pdfHardBPtr != pdfBPtr) delete pdfHardBPtr; |
---|
| 255 | if (useNewPdfA) delete pdfAPtr; |
---|
| 256 | if (useNewPdfB) delete pdfBPtr; |
---|
| 257 | if (useNewPdfPomA) delete pdfPomAPtr; |
---|
| 258 | if (useNewPdfPomB) delete pdfPomBPtr; |
---|
| 259 | |
---|
| 260 | // Reset pointers to be empty. |
---|
| 261 | useNewPdfA = false; |
---|
| 262 | useNewPdfB = false; |
---|
| 263 | useNewPdfHard = false; |
---|
| 264 | useNewPdfPomA = false; |
---|
| 265 | useNewPdfPomB = false; |
---|
| 266 | pdfAPtr = 0; |
---|
| 267 | pdfBPtr = 0; |
---|
| 268 | pdfHardAPtr = 0; |
---|
| 269 | pdfHardBPtr = 0; |
---|
| 270 | pdfPomAPtr = 0; |
---|
| 271 | pdfPomBPtr = 0; |
---|
| 272 | |
---|
| 273 | // Switch off external PDF's by zero as input. |
---|
| 274 | if (pdfAPtrIn == 0 && pdfBPtrIn == 0) return true; |
---|
| 275 | |
---|
| 276 | // The two PDF objects cannot be one and the same. |
---|
| 277 | if (pdfAPtrIn == pdfBPtrIn) return false; |
---|
| 278 | |
---|
| 279 | // Save pointers. |
---|
| 280 | pdfAPtr = pdfAPtrIn; |
---|
| 281 | pdfBPtr = pdfBPtrIn; |
---|
| 282 | |
---|
| 283 | // By default same pointers for hard-process PDF's. |
---|
| 284 | pdfHardAPtr = pdfAPtrIn; |
---|
| 285 | pdfHardBPtr = pdfBPtrIn; |
---|
| 286 | |
---|
| 287 | // Optionally allow separate pointers for hard process. |
---|
| 288 | if (pdfHardAPtrIn != 0 && pdfHardBPtrIn != 0) { |
---|
| 289 | if (pdfHardAPtrIn == pdfHardBPtrIn) return false; |
---|
| 290 | pdfHardAPtr = pdfHardAPtrIn; |
---|
| 291 | pdfHardBPtr = pdfHardBPtrIn; |
---|
| 292 | } |
---|
| 293 | |
---|
| 294 | // Optionally allow pointers for Pomerons in the proton. |
---|
| 295 | if (pdfPomAPtrIn != 0 && pdfPomBPtrIn != 0) { |
---|
| 296 | if (pdfPomAPtrIn == pdfPomBPtrIn) return false; |
---|
| 297 | pdfPomAPtr = pdfPomAPtrIn; |
---|
| 298 | pdfPomBPtr = pdfPomBPtrIn; |
---|
| 299 | } |
---|
| 300 | |
---|
| 301 | // Done. |
---|
| 302 | return true; |
---|
| 303 | } |
---|
| 304 | |
---|
| 305 | //-------------------------------------------------------------------------- |
---|
| 306 | |
---|
| 307 | // Routine to initialize with the variable values of the Beams kind. |
---|
| 308 | |
---|
| 309 | bool Pythia::init() { |
---|
| 310 | |
---|
| 311 | // Check that constructor worked. |
---|
| 312 | isInit = false; |
---|
| 313 | if (!isConstructed) { |
---|
| 314 | info.errorMsg("Abort from Pythia::init: constructur " |
---|
| 315 | "initialization failed"); |
---|
| 316 | return false; |
---|
| 317 | } |
---|
| 318 | |
---|
| 319 | // Begin initialization. Find which frame type to use. |
---|
| 320 | info.addCounter(1); |
---|
| 321 | frameType = mode("Beams:frameType"); |
---|
| 322 | |
---|
| 323 | // Initialization with internal processes: read in and set values. |
---|
| 324 | if (frameType < 4 ) { |
---|
| 325 | doLHA = false; |
---|
| 326 | boostType = frameType; |
---|
| 327 | idA = mode("Beams:idA"); |
---|
| 328 | idB = mode("Beams:idB"); |
---|
| 329 | eCM = parm("Beams:eCM"); |
---|
| 330 | eA = parm("Beams:eA"); |
---|
| 331 | eB = parm("Beams:eB"); |
---|
| 332 | pxA = parm("Beams:pxA"); |
---|
| 333 | pyA = parm("Beams:pyA"); |
---|
| 334 | pzA = parm("Beams:pzA"); |
---|
| 335 | pxB = parm("Beams:pxB"); |
---|
| 336 | pyB = parm("Beams:pyB"); |
---|
| 337 | pzB = parm("Beams:pzB"); |
---|
| 338 | |
---|
| 339 | // Initialization with a Les Houches Event File or an LHAup object. |
---|
| 340 | } else { |
---|
| 341 | doLHA = true; |
---|
| 342 | boostType = 2; |
---|
| 343 | string lhef = word("Beams:LHEF"); |
---|
| 344 | string lhefHeader = word("Beams:LHEFheader"); |
---|
| 345 | bool readHeaders = flag("Beams:readLHEFheaders"); |
---|
| 346 | bool skipInit = flag("Beams:newLHEFsameInit"); |
---|
| 347 | int nSkipAtInit = mode("Beams:nSkipLHEFatInit"); |
---|
| 348 | |
---|
| 349 | // For file input: renew file stream or (re)new Les Houches object. |
---|
| 350 | if (frameType == 4) { |
---|
| 351 | const char* cstring1 = lhef.c_str(); |
---|
| 352 | if (useNewLHA && skipInit) lhaUpPtr->newEventFile(cstring1); |
---|
| 353 | else { |
---|
| 354 | if (useNewLHA) delete lhaUpPtr; |
---|
| 355 | // Header is optional, so use NULL pointer to indicate no value. |
---|
| 356 | const char* cstring2 = (lhefHeader == "void") ? |
---|
| 357 | NULL : lhefHeader.c_str(); |
---|
| 358 | lhaUpPtr = new LHAupLHEF(cstring1, cstring2, readHeaders); |
---|
| 359 | useNewLHA = true; |
---|
| 360 | } |
---|
| 361 | |
---|
| 362 | // Check that file was properly opened. |
---|
| 363 | if (!lhaUpPtr->fileFound()) { |
---|
| 364 | info.errorMsg("Abort from Pythia::init: " |
---|
| 365 | "Les Houches Event File not found"); |
---|
| 366 | return false; |
---|
| 367 | } |
---|
| 368 | |
---|
| 369 | // For object input: at least check that not null pointer. |
---|
| 370 | } else { |
---|
| 371 | if (lhaUpPtr == 0) { |
---|
| 372 | info.errorMsg("Abort from Pythia::init: " |
---|
| 373 | "LHAup object not found"); |
---|
| 374 | return false; |
---|
| 375 | } |
---|
| 376 | |
---|
| 377 | // LHAup object generic abort using fileFound() routine |
---|
| 378 | if (!lhaUpPtr->fileFound()) { |
---|
| 379 | info.errorMsg("Abort from Pythia::init: " |
---|
| 380 | "LHAup initialisation error"); |
---|
| 381 | return false; |
---|
| 382 | } |
---|
| 383 | } |
---|
| 384 | |
---|
| 385 | // Send in pointer to info. Store or replace LHA pointer in other classes. |
---|
| 386 | lhaUpPtr->setPtr( &info); |
---|
| 387 | processLevel.setLHAPtr( lhaUpPtr); |
---|
| 388 | |
---|
| 389 | // If second time around, only with new file, then simplify. |
---|
| 390 | // Optionally skip ahead a number of events at beginning of file. |
---|
| 391 | if (skipInit) { |
---|
| 392 | isInit = true; |
---|
| 393 | info.addCounter(2); |
---|
| 394 | if (nSkipAtInit > 0) lhaUpPtr->skipEvent(nSkipAtInit); |
---|
| 395 | return true; |
---|
| 396 | } |
---|
| 397 | |
---|
| 398 | // Set up values related to merging hooks. |
---|
| 399 | if (frameType == 4) { |
---|
| 400 | // Set up values related to CKKW-L merging. |
---|
| 401 | doUserMerging = settings.flag("Merging:doUserMerging"); |
---|
| 402 | doMGMerging = settings.flag("Merging:doMGMerging"); |
---|
| 403 | doKTMerging = settings.flag("Merging:doKTMerging"); |
---|
| 404 | doPTLundMerging = settings.flag("Merging:doPTLundMerging"); |
---|
| 405 | doCutBasedMerging = settings.flag("Merging:doCutBasedMerging"); |
---|
| 406 | doMerging = doUserMerging || doMGMerging || doKTMerging |
---|
| 407 | || doPTLundMerging || doCutBasedMerging; |
---|
| 408 | // Set up MergingHooks object |
---|
| 409 | if (!doUserMerging){ |
---|
| 410 | if (hasOwnMergingHooks && mergingHooksPtr) delete mergingHooksPtr; |
---|
| 411 | mergingHooksPtr = new MergingHooks(); |
---|
| 412 | hasOwnMergingHooks = true; |
---|
| 413 | } |
---|
| 414 | hasMergingHooks = (mergingHooksPtr != 0); |
---|
| 415 | // Merging hooks required for merging. If no merging hooks pointer is |
---|
| 416 | // available, exit. |
---|
| 417 | if (!hasMergingHooks) { |
---|
| 418 | info.errorMsg("Abort from Pythia::init: " |
---|
| 419 | "no merging hooks object has been provided"); |
---|
| 420 | return false; |
---|
| 421 | } else mergingHooksPtr->setLHEInputFile( lhef); |
---|
| 422 | // Initialise counting of Les Houches Events significantly above the |
---|
| 423 | // merging scale. |
---|
| 424 | info.setCounter(41,0); |
---|
| 425 | } |
---|
| 426 | |
---|
| 427 | // Set LHAinit information (in some external program). |
---|
| 428 | if (settings.flag("ProcessLevel:all") && !lhaUpPtr->setInit()) { |
---|
| 429 | info.errorMsg("Abort from Pythia::init: " |
---|
| 430 | "Les Houches initialization failed"); |
---|
| 431 | return false; |
---|
| 432 | } |
---|
| 433 | |
---|
| 434 | // Extract beams from values set in an LHAinit object. |
---|
| 435 | idA = lhaUpPtr->idBeamA(); |
---|
| 436 | idB = lhaUpPtr->idBeamB(); |
---|
| 437 | eA = lhaUpPtr->eBeamA(); |
---|
| 438 | eB = lhaUpPtr->eBeamB(); |
---|
| 439 | // Optionally skip ahead a number of events at beginning of file. |
---|
| 440 | if (nSkipAtInit > 0) lhaUpPtr->skipEvent(nSkipAtInit); |
---|
| 441 | } |
---|
| 442 | |
---|
| 443 | // Set up values related to user hooks. |
---|
| 444 | hasUserHooks = (userHooksPtr != 0); |
---|
| 445 | doVetoProcess = false; |
---|
| 446 | doVetoPartons = false; |
---|
| 447 | if (hasUserHooks) { |
---|
| 448 | userHooksPtr->initPtr( &info, &settings, &particleData, &rndm, |
---|
| 449 | &beamA, &beamB, &beamPomA, &beamPomB, couplingsPtr, &partonSystems, |
---|
| 450 | &sigmaTot); |
---|
| 451 | if (!userHooksPtr->initAfterBeams()) { |
---|
| 452 | info.errorMsg("Abort from Pythia::init: could not initialise UserHooks"); |
---|
| 453 | return false; |
---|
| 454 | } |
---|
| 455 | doVetoProcess = userHooksPtr->canVetoProcessLevel(); |
---|
| 456 | doVetoPartons = userHooksPtr->canVetoPartonLevel(); |
---|
| 457 | } |
---|
| 458 | |
---|
| 459 | // Back to common initialization. Reset error counters. |
---|
| 460 | nErrEvent = 0; |
---|
| 461 | info.errorReset(); |
---|
| 462 | info.setTooLowPTmin(false); |
---|
| 463 | info.sigmaReset(); |
---|
| 464 | |
---|
| 465 | // Initialize data members extracted from database. |
---|
| 466 | doProcessLevel = settings.flag("ProcessLevel:all"); |
---|
| 467 | doPartonLevel = settings.flag("PartonLevel:all"); |
---|
| 468 | doHadronLevel = settings.flag("HadronLevel:all"); |
---|
| 469 | doDiffraction = settings.flag("SoftQCD:all") |
---|
| 470 | || settings.flag("SoftQCD:singleDiffractive") |
---|
| 471 | || settings.flag("SoftQCD:doubleDiffractive") |
---|
| 472 | || settings.flag("SoftQCD:centralDiffractive"); |
---|
| 473 | doResDec = settings.flag("Standalone:allowResDec"); |
---|
| 474 | doFSRinRes = doPartonLevel && settings.flag("PartonLevel:FSR") |
---|
| 475 | && settings.flag("PartonLevel:FSRinResonances"); |
---|
| 476 | decayRHadrons = settings.flag("RHadrons:allowDecay"); |
---|
| 477 | doMomentumSpread = settings.flag("Beams:allowMomentumSpread"); |
---|
| 478 | doVertexSpread = settings.flag("Beams:allowVertexSpread"); |
---|
| 479 | abortIfVeto = settings.flag("Check:abortIfVeto"); |
---|
| 480 | checkEvent = settings.flag("Check:event"); |
---|
| 481 | checkHistory = settings.flag("Check:history"); |
---|
| 482 | nErrList = settings.mode("Check:nErrList"); |
---|
| 483 | epTolErr = settings.parm("Check:epTolErr"); |
---|
| 484 | epTolWarn = settings.parm("Check:epTolWarn"); |
---|
| 485 | |
---|
| 486 | |
---|
| 487 | // Initialise merging hooks. |
---|
| 488 | if (hasMergingHooks && doMerging ) |
---|
| 489 | mergingHooksPtr->init( settings, &info, &particleData ); |
---|
| 490 | |
---|
| 491 | // Initialize the random number generator. |
---|
| 492 | if ( settings.flag("Random:setSeed") ) |
---|
| 493 | rndm.init( settings.mode("Random:seed") ); |
---|
| 494 | |
---|
| 495 | // Check that combinations of settings are allowed; change if not. |
---|
| 496 | checkSettings(); |
---|
| 497 | |
---|
| 498 | // Initialize couplings (needed to initialize resonances). |
---|
| 499 | // Check if SUSY couplings need to be read in |
---|
| 500 | if( !initSLHA()) info.errorMsg("Error in Pythia::init: " |
---|
| 501 | "Could not read SLHA file"); |
---|
| 502 | if (couplingsPtr->isSUSY) { |
---|
| 503 | // Initialize the SM and SUSY. |
---|
| 504 | coupSUSY.init( settings, &rndm); |
---|
| 505 | coupSUSY.initSUSY(&slha, &settings, &particleData); |
---|
| 506 | couplingsPtr = (Couplings *) &coupSUSY; |
---|
| 507 | } else { |
---|
| 508 | // Initialize the SM couplings. |
---|
| 509 | couplingsPtr->init( settings, &rndm); |
---|
| 510 | } |
---|
| 511 | |
---|
| 512 | // Reset couplingsPtr to the correct place. |
---|
| 513 | particleData.initPtr( &info, &settings, &rndm, couplingsPtr); |
---|
| 514 | |
---|
| 515 | // Set headers to distinguish the two event listing kinds. |
---|
| 516 | int startColTag = settings.mode("Event:startColTag"); |
---|
| 517 | process.init("(hard process)", &particleData, startColTag); |
---|
| 518 | event.init("(complete event)", &particleData, startColTag); |
---|
| 519 | |
---|
| 520 | // Final setup stage of particle data, notably resonance widths. |
---|
| 521 | particleData.initWidths( resonancePtrs); |
---|
| 522 | |
---|
| 523 | // Set up R-hadrons particle data, where relevant. |
---|
| 524 | rHadrons.init( &info, settings, &particleData, &rndm); |
---|
| 525 | |
---|
| 526 | // Set up objects for timelike and spacelike showers. |
---|
| 527 | if (timesDecPtr == 0 || timesPtr == 0) { |
---|
| 528 | TimeShower* timesNow = new TimeShower(); |
---|
| 529 | if (timesDecPtr == 0) timesDecPtr = timesNow; |
---|
| 530 | if (timesPtr == 0) timesPtr = timesNow; |
---|
| 531 | useNewTimes = true; |
---|
| 532 | } |
---|
| 533 | if (spacePtr == 0) { |
---|
| 534 | spacePtr = new SpaceShower(); |
---|
| 535 | useNewSpace = true; |
---|
| 536 | } |
---|
| 537 | |
---|
| 538 | // Initialize showers, especially for simple showers in decays. |
---|
| 539 | timesPtr->initPtr( &info, &settings, &particleData, &rndm, couplingsPtr, |
---|
| 540 | &partonSystems, userHooksPtr); |
---|
| 541 | timesDecPtr->initPtr( &info, &settings, &particleData, &rndm, couplingsPtr, |
---|
| 542 | &partonSystems, userHooksPtr); |
---|
| 543 | spacePtr->initPtr( &info, &settings, &particleData, &rndm, |
---|
| 544 | &partonSystems, userHooksPtr); |
---|
| 545 | timesDecPtr->init( 0, 0); |
---|
| 546 | |
---|
| 547 | // Set up values related to beam shape. |
---|
| 548 | if (beamShapePtr == 0) { |
---|
| 549 | beamShapePtr = new BeamShape(); |
---|
| 550 | useNewBeamShape = true; |
---|
| 551 | } |
---|
| 552 | beamShapePtr->init( settings, &rndm); |
---|
| 553 | |
---|
| 554 | // Check that beams and beam combination can be handled. |
---|
| 555 | if (!checkBeams()) { |
---|
| 556 | info.errorMsg("Abort from Pythia::init: " |
---|
| 557 | "checkBeams initialization failed"); |
---|
| 558 | return false; |
---|
| 559 | } |
---|
| 560 | |
---|
| 561 | // Do not set up beam kinematics when no process level. |
---|
| 562 | if (!doProcessLevel) boostType = 1; |
---|
| 563 | else { |
---|
| 564 | |
---|
| 565 | // Set up beam kinematics. |
---|
| 566 | if (!initKinematics()) { |
---|
| 567 | info.errorMsg("Abort from Pythia::init: " |
---|
| 568 | "kinematics initialization failed"); |
---|
| 569 | return false; |
---|
| 570 | } |
---|
| 571 | |
---|
| 572 | // Set up pointers to PDFs. |
---|
| 573 | if (!initPDFs()) { |
---|
| 574 | info.errorMsg("Abort from Pythia::init: PDF initialization failed"); |
---|
| 575 | return false; |
---|
| 576 | } |
---|
| 577 | |
---|
| 578 | // Set up the two beams and the common remnant system. |
---|
| 579 | StringFlav* flavSelPtr = hadronLevel.getStringFlavPtr(); |
---|
| 580 | beamA.init( idA, pzAcm, eA, mA, &info, settings, &particleData, &rndm, |
---|
| 581 | pdfAPtr, pdfHardAPtr, isUnresolvedA, flavSelPtr); |
---|
| 582 | beamB.init( idB, pzBcm, eB, mB, &info, settings, &particleData, &rndm, |
---|
| 583 | pdfBPtr, pdfHardBPtr, isUnresolvedB, flavSelPtr); |
---|
| 584 | |
---|
| 585 | // Optionally set up new alternative beams for these Pomerons. |
---|
| 586 | if ( doDiffraction) { |
---|
| 587 | beamPomA.init( 990, 0.5 * eCM, 0.5 * eCM, 0., &info, settings, |
---|
| 588 | &particleData, &rndm, pdfPomAPtr, pdfPomAPtr, false, flavSelPtr); |
---|
| 589 | beamPomB.init( 990, -0.5 * eCM, 0.5 * eCM, 0., &info, settings, |
---|
| 590 | &particleData, &rndm, pdfPomBPtr, pdfPomBPtr, false, flavSelPtr); |
---|
| 591 | } |
---|
| 592 | } |
---|
| 593 | |
---|
| 594 | // Send info/pointers to process level for initialization. |
---|
| 595 | if ( doProcessLevel && !processLevel.init( &info, settings, &particleData, |
---|
| 596 | &rndm, &beamA, &beamB, couplingsPtr, &sigmaTot, doLHA, &slha, userHooksPtr, |
---|
| 597 | sigmaPtrs) ) { |
---|
| 598 | info.errorMsg("Abort from Pythia::init: " |
---|
| 599 | "processLevel initialization failed"); |
---|
| 600 | return false; |
---|
| 601 | } |
---|
| 602 | |
---|
| 603 | // Optionally only initialize resonance decays. |
---|
| 604 | if ( !doProcessLevel && doResDec) processLevel.initDecays( &info, |
---|
| 605 | &particleData, &rndm, lhaUpPtr); |
---|
| 606 | |
---|
| 607 | // Send info/pointers to parton level for initialization. |
---|
| 608 | if ( doPartonLevel && doProcessLevel && !partonLevel.init( &info, settings, |
---|
| 609 | &particleData, &rndm, &beamA, &beamB, &beamPomA, &beamPomB, couplingsPtr, |
---|
| 610 | &partonSystems, &sigmaTot, timesDecPtr, timesPtr, spacePtr, &rHadrons, |
---|
| 611 | userHooksPtr, mergingHooksPtr, false) ) { |
---|
| 612 | info.errorMsg("Abort from Pythia::init: " |
---|
| 613 | "partonLevel initialization failed" ); |
---|
| 614 | return false; |
---|
| 615 | } |
---|
| 616 | |
---|
| 617 | // Optionally only initialize final-state showers in resonance decays. |
---|
| 618 | if ( (!doProcessLevel || !doPartonLevel) && doFSRinRes) partonLevel.init( |
---|
| 619 | &info, settings, &particleData, &rndm, 0, 0, 0, 0, couplingsPtr, |
---|
| 620 | &partonSystems, 0, timesDecPtr, 0, 0, &rHadrons, 0, 0, false); |
---|
| 621 | |
---|
| 622 | // Send info/pointers to parton level for trial shower initialization. |
---|
| 623 | if ( doMerging |
---|
| 624 | && !trialPartonLevel.init( &info, settings, &particleData, |
---|
| 625 | &rndm, &beamA, &beamB, &beamPomA, &beamPomB, couplingsPtr, |
---|
| 626 | &partonSystems, &sigmaTot, timesDecPtr, timesPtr, spacePtr, &rHadrons, |
---|
| 627 | NULL, mergingHooksPtr, true) ) { |
---|
| 628 | info.errorMsg("Abort from Pythia::init: " |
---|
| 629 | "trialPartonLevel initialization failed"); |
---|
| 630 | return false; |
---|
| 631 | } |
---|
| 632 | |
---|
| 633 | // Send info/pointers to hadron level for initialization. |
---|
| 634 | // Note: forceHadronLevel() can come, so we must always initialize. |
---|
| 635 | if ( !hadronLevel.init( &info, settings, &particleData, &rndm, |
---|
| 636 | couplingsPtr, timesDecPtr, &rHadrons, decayHandlePtr, |
---|
| 637 | handledParticles) ) { |
---|
| 638 | info.errorMsg("Abort from Pythia::init: " |
---|
| 639 | "hadronLevel initialization failed"); |
---|
| 640 | return false; |
---|
| 641 | } |
---|
| 642 | |
---|
| 643 | // Optionally check particle data table for inconsistencies. |
---|
| 644 | if ( settings.flag("Check:particleData") ) |
---|
| 645 | particleData.checkTable( settings.mode("Check:levelParticleData") ); |
---|
| 646 | |
---|
| 647 | // Optionally show settings and particle data, changed or all. |
---|
| 648 | bool showCS = settings.flag("Init:showChangedSettings"); |
---|
| 649 | bool showAS = settings.flag("Init:showAllSettings"); |
---|
| 650 | bool showCPD = settings.flag("Init:showChangedParticleData"); |
---|
| 651 | bool showCRD = settings.flag("Init:showChangedResonanceData"); |
---|
| 652 | bool showAPD = settings.flag("Init:showAllParticleData"); |
---|
| 653 | int show1PD = settings.mode("Init:showOneParticleData"); |
---|
| 654 | bool showPro = settings.flag("Init:showProcesses"); |
---|
| 655 | if (showCS) settings.listChanged(); |
---|
| 656 | if (showAS) settings.listAll(); |
---|
| 657 | if (show1PD > 0) particleData.list(show1PD); |
---|
| 658 | if (showCPD) particleData.listChanged(showCRD); |
---|
| 659 | if (showAPD) particleData.listAll(); |
---|
| 660 | |
---|
| 661 | // Listing options for the next() routine. |
---|
| 662 | nCount = settings.mode("Next:numberCount"); |
---|
| 663 | nShowLHA = settings.mode("Next:numberShowLHA"); |
---|
| 664 | nShowInfo = settings.mode("Next:numberShowInfo"); |
---|
| 665 | nShowProc = settings.mode("Next:numberShowProcess"); |
---|
| 666 | nShowEvt = settings.mode("Next:numberShowEvent"); |
---|
| 667 | showSaV = settings.flag("Next:showScaleAndVertex"); |
---|
| 668 | showMaD = settings.flag("Next:showMothersAndDaughters"); |
---|
| 669 | |
---|
| 670 | // Succeeded. |
---|
| 671 | isInit = true; |
---|
| 672 | info.addCounter(2); |
---|
| 673 | if (useNewLHA && showPro) lhaUpPtr->listInit(); |
---|
| 674 | return true; |
---|
| 675 | |
---|
| 676 | } |
---|
| 677 | |
---|
| 678 | //-------------------------------------------------------------------------- |
---|
| 679 | |
---|
| 680 | // Routine to initialize with CM energy. |
---|
| 681 | |
---|
| 682 | bool Pythia::init( int idAin, int idBin, double eCMin) { |
---|
| 683 | |
---|
| 684 | // Set arguments in Settings database. |
---|
| 685 | settings.mode("Beams:idA", idAin); |
---|
| 686 | settings.mode("Beams:idB", idBin); |
---|
| 687 | settings.mode("Beams:frameType", 1); |
---|
| 688 | settings.parm("Beams:eCM", eCMin); |
---|
| 689 | |
---|
| 690 | // Send on to common initialization. |
---|
| 691 | return init(); |
---|
| 692 | |
---|
| 693 | } |
---|
| 694 | |
---|
| 695 | //-------------------------------------------------------------------------- |
---|
| 696 | |
---|
| 697 | // Routine to initialize with two collinear beams, energies specified. |
---|
| 698 | |
---|
| 699 | bool Pythia::init( int idAin, int idBin, double eAin, double eBin) { |
---|
| 700 | |
---|
| 701 | // Set arguments in Settings database. |
---|
| 702 | settings.mode("Beams:idA", idAin); |
---|
| 703 | settings.mode("Beams:idB", idBin); |
---|
| 704 | settings.mode("Beams:frameType", 2); |
---|
| 705 | settings.parm("Beams:eA", eAin); |
---|
| 706 | settings.parm("Beams:eB", eBin); |
---|
| 707 | |
---|
| 708 | // Send on to common initialization. |
---|
| 709 | return init(); |
---|
| 710 | |
---|
| 711 | } |
---|
| 712 | |
---|
| 713 | //-------------------------------------------------------------------------- |
---|
| 714 | |
---|
| 715 | // Routine to initialize with two beams specified by three-momenta. |
---|
| 716 | |
---|
| 717 | bool Pythia::init( int idAin, int idBin, double pxAin, double pyAin, |
---|
| 718 | double pzAin, double pxBin, double pyBin, double pzBin) { |
---|
| 719 | |
---|
| 720 | // Set arguments in Settings database. |
---|
| 721 | settings.mode("Beams:idA", idAin); |
---|
| 722 | settings.mode("Beams:idB", idBin); |
---|
| 723 | settings.mode("Beams:frameType", 3); |
---|
| 724 | settings.parm("Beams:pxA", pxAin); |
---|
| 725 | settings.parm("Beams:pyA", pyAin); |
---|
| 726 | settings.parm("Beams:pzA", pzAin); |
---|
| 727 | settings.parm("Beams:pxB", pxBin); |
---|
| 728 | settings.parm("Beams:pyB", pyBin); |
---|
| 729 | settings.parm("Beams:pzB", pzBin); |
---|
| 730 | |
---|
| 731 | // Send on to common initialization. |
---|
| 732 | return init(); |
---|
| 733 | |
---|
| 734 | } |
---|
| 735 | |
---|
| 736 | //-------------------------------------------------------------------------- |
---|
| 737 | |
---|
| 738 | // Routine to initialize when all info is given in a Les Houches Event File. |
---|
| 739 | |
---|
| 740 | bool Pythia::init( string LesHouchesEventFile, bool skipInit) { |
---|
| 741 | |
---|
| 742 | // Set arguments in Settings database. |
---|
| 743 | settings.mode("Beams:frameType", 4); |
---|
| 744 | settings.word("Beams:LHEF", LesHouchesEventFile); |
---|
| 745 | settings.flag("Beams:newLHEFsameInit", skipInit); |
---|
| 746 | |
---|
| 747 | // Send on to common initialization. |
---|
| 748 | return init(); |
---|
| 749 | |
---|
| 750 | } |
---|
| 751 | |
---|
| 752 | //-------------------------------------------------------------------------- |
---|
| 753 | |
---|
| 754 | // Routine to initialize when beam info is given in an LHAup object. |
---|
| 755 | |
---|
| 756 | bool Pythia::init( LHAup* lhaUpPtrIn) { |
---|
| 757 | |
---|
| 758 | // Store LHAup object and set arguments in Settings database. |
---|
| 759 | setLHAupPtr( lhaUpPtrIn); |
---|
| 760 | settings.mode("Beams:frameType", 5); |
---|
| 761 | |
---|
| 762 | // Send on to common initialization. |
---|
| 763 | return init(); |
---|
| 764 | |
---|
| 765 | } |
---|
| 766 | |
---|
| 767 | //-------------------------------------------------------------------------- |
---|
| 768 | |
---|
| 769 | // Check that combinations of settings are allowed; change if not. |
---|
| 770 | |
---|
| 771 | void Pythia::checkSettings() { |
---|
| 772 | |
---|
| 773 | // Double rescattering not allowed if ISR or FSR. |
---|
| 774 | if ((settings.flag("PartonLevel:ISR") || settings.flag("PartonLevel:FSR")) |
---|
| 775 | && settings.flag("MultipartonInteractions:allowDoubleRescatter")) { |
---|
| 776 | info.errorMsg("Warning in Pythia::checkSettings: " |
---|
| 777 | "double rescattering switched off since showering is on"); |
---|
| 778 | settings.flag("MultipartonInteractions:allowDoubleRescatter", false); |
---|
| 779 | } |
---|
| 780 | |
---|
| 781 | } |
---|
| 782 | |
---|
| 783 | //-------------------------------------------------------------------------- |
---|
| 784 | |
---|
| 785 | // Check that beams and beam combination can be handled. Set up unresolved. |
---|
| 786 | |
---|
| 787 | bool Pythia::checkBeams() { |
---|
| 788 | |
---|
| 789 | // Absolute flavours. If not to do process level then no check needed. |
---|
| 790 | int idAabs = abs(idA); |
---|
| 791 | int idBabs = abs(idB); |
---|
| 792 | if (!doProcessLevel) return true; |
---|
| 793 | |
---|
| 794 | // Neutrino beams always unresolved, charged lepton ones conditionally. |
---|
| 795 | bool isLeptonA = (idAabs > 10 && idAabs < 17); |
---|
| 796 | bool isLeptonB = (idBabs > 10 && idBabs < 17); |
---|
| 797 | bool isUnresLep = !settings.flag("PDF:lepton"); |
---|
| 798 | isUnresolvedA = isLeptonA && (idAabs%2 == 0 || isUnresLep); |
---|
| 799 | isUnresolvedB = isLeptonB && (idBabs%2 == 0 || isUnresLep); |
---|
| 800 | |
---|
| 801 | // Lepton-lepton collisions OK (including neutrinos) if both (un)resolved. |
---|
| 802 | if (isLeptonA && isLeptonB && isUnresolvedA == isUnresolvedB) return true; |
---|
| 803 | |
---|
| 804 | // MBR model only implemented for pp/ppbar/pbarp collisions. |
---|
| 805 | int PomFlux = settings.mode("Diffraction:PomFlux"); |
---|
| 806 | if (PomFlux == 5) { |
---|
| 807 | bool ispp = (idAabs == 2212 && idBabs == 2212); |
---|
| 808 | bool ispbarpbar = (idA == -2212 && idB == -2212); |
---|
| 809 | if (ispp && !ispbarpbar) return true; |
---|
| 810 | info.errorMsg("Error in Pythia::init: cannot handle this beam combination" |
---|
| 811 | " with PomFlux == 5"); |
---|
| 812 | return false; |
---|
| 813 | } |
---|
| 814 | |
---|
| 815 | // Hadron-hadron collisions OK, with Pomeron counted as hadron. |
---|
| 816 | bool isHadronA = (idAabs == 2212) || (idA == 111) || (idAabs == 211) |
---|
| 817 | || (idA == 990); |
---|
| 818 | bool isHadronB = (idBabs == 2212) || (idA == 111)|| (idBabs == 211) |
---|
| 819 | || (idB == 990); |
---|
| 820 | if (isHadronA && isHadronB) return true; |
---|
| 821 | |
---|
| 822 | // If no case above then failed. |
---|
| 823 | info.errorMsg("Error in Pythia::init: cannot handle this beam combination"); |
---|
| 824 | return false; |
---|
| 825 | |
---|
| 826 | } |
---|
| 827 | |
---|
| 828 | //-------------------------------------------------------------------------- |
---|
| 829 | |
---|
| 830 | // Calculate kinematics at initialization. Store beam four-momenta. |
---|
| 831 | |
---|
| 832 | bool Pythia::initKinematics() { |
---|
| 833 | |
---|
| 834 | // Find masses. Initial guess that we are in CM frame. |
---|
| 835 | mA = particleData.m0(idA); |
---|
| 836 | mB = particleData.m0(idB); |
---|
| 837 | betaZ = 0.; |
---|
| 838 | gammaZ = 1.; |
---|
| 839 | |
---|
| 840 | // Collinear beams not in CM frame: find CM energy. |
---|
| 841 | if (boostType == 2) { |
---|
| 842 | eA = max(eA, mA); |
---|
| 843 | eB = max(eB, mB); |
---|
| 844 | pzA = sqrt(eA*eA - mA*mA); |
---|
| 845 | pzB = -sqrt(eB*eB - mB*mB); |
---|
| 846 | pAinit = Vec4( 0., 0., pzA, eA); |
---|
| 847 | pBinit = Vec4( 0., 0., pzB, eB); |
---|
| 848 | eCM = sqrt( pow2(eA + eB) - pow2(pzA + pzB) ); |
---|
| 849 | |
---|
| 850 | // Find boost to rest frame. |
---|
| 851 | betaZ = (pzA + pzB) / (eA + eB); |
---|
| 852 | gammaZ = (eA + eB) / eCM; |
---|
| 853 | if (abs(betaZ) < 1e-10) boostType = 1; |
---|
| 854 | } |
---|
| 855 | |
---|
| 856 | // Completely general beam directions: find CM energy. |
---|
| 857 | else if (boostType == 3) { |
---|
| 858 | eA = sqrt( pxA*pxA + pyA*pyA + pzA*pzA + mA*mA); |
---|
| 859 | eB = sqrt( pxB*pxB + pyB*pyB + pzB*pzB + mB*mB); |
---|
| 860 | pAinit = Vec4( pxA, pyA, pzA, eA); |
---|
| 861 | pBinit = Vec4( pxB, pyB, pzB, eB); |
---|
| 862 | eCM = (pAinit + pBinit).mCalc(); |
---|
| 863 | |
---|
| 864 | // Find boost+rotation needed to move from/to CM frame. |
---|
| 865 | MfromCM.reset(); |
---|
| 866 | MfromCM.fromCMframe( pAinit, pBinit); |
---|
| 867 | MtoCM = MfromCM; |
---|
| 868 | MtoCM.invert(); |
---|
| 869 | } |
---|
| 870 | |
---|
| 871 | // Fail if CM energy below beam masses. |
---|
| 872 | if (eCM < mA + mB) { |
---|
| 873 | info.errorMsg("Error in Pythia::initKinematics: too low energy"); |
---|
| 874 | return false; |
---|
| 875 | } |
---|
| 876 | |
---|
| 877 | // Set up CM-frame kinematics with beams along +-z axis. |
---|
| 878 | pzAcm = 0.5 * sqrtpos( (eCM + mA + mB) * (eCM - mA - mB) |
---|
| 879 | * (eCM - mA + mB) * (eCM + mA - mB) ) / eCM; |
---|
| 880 | pzBcm = -pzAcm; |
---|
| 881 | eA = sqrt(mA*mA + pzAcm*pzAcm); |
---|
| 882 | eB = sqrt(mB*mB + pzBcm*pzBcm); |
---|
| 883 | |
---|
| 884 | // If in CM frame then store beam four-vectors (else already done above). |
---|
| 885 | if (boostType != 2 && boostType != 3) { |
---|
| 886 | pAinit = Vec4( 0., 0., pzAcm, eA); |
---|
| 887 | pBinit = Vec4( 0., 0., pzBcm, eB); |
---|
| 888 | } |
---|
| 889 | |
---|
| 890 | // Store main info for access in process generation. |
---|
| 891 | info.setBeamA( idA, pzAcm, eA, mA); |
---|
| 892 | info.setBeamB( idB, pzBcm, eB, mB); |
---|
| 893 | info.setECM( eCM); |
---|
| 894 | |
---|
| 895 | // Must allow for generic boost+rotation when beam momentum spread. |
---|
| 896 | if (doMomentumSpread) boostType = 3; |
---|
| 897 | |
---|
| 898 | // Done. |
---|
| 899 | return true; |
---|
| 900 | |
---|
| 901 | } |
---|
| 902 | |
---|
| 903 | //-------------------------------------------------------------------------- |
---|
| 904 | |
---|
| 905 | // Set up pointers to PDFs. |
---|
| 906 | |
---|
| 907 | bool Pythia::initPDFs() { |
---|
| 908 | |
---|
| 909 | // Delete any PDF's created in a previous init call. |
---|
| 910 | if (useNewPdfHard) { |
---|
| 911 | if (pdfHardAPtr != pdfAPtr) { |
---|
| 912 | delete pdfHardAPtr; |
---|
| 913 | pdfHardAPtr = 0; |
---|
| 914 | } |
---|
| 915 | if (pdfHardBPtr != pdfBPtr) { |
---|
| 916 | delete pdfHardBPtr; |
---|
| 917 | pdfHardBPtr = 0; |
---|
| 918 | } |
---|
| 919 | useNewPdfHard = false; |
---|
| 920 | } |
---|
| 921 | if (useNewPdfA) { |
---|
| 922 | delete pdfAPtr; |
---|
| 923 | useNewPdfA = false; |
---|
| 924 | pdfAPtr = 0; |
---|
| 925 | } |
---|
| 926 | if (useNewPdfB) { |
---|
| 927 | delete pdfBPtr; |
---|
| 928 | useNewPdfB = false; |
---|
| 929 | pdfBPtr = 0; |
---|
| 930 | } |
---|
| 931 | if (useNewPdfPomA) { |
---|
| 932 | delete pdfPomAPtr; |
---|
| 933 | useNewPdfPomA = false; |
---|
| 934 | pdfPomAPtr = 0; |
---|
| 935 | } |
---|
| 936 | if (useNewPdfPomB) { |
---|
| 937 | delete pdfPomBPtr; |
---|
| 938 | useNewPdfPomB = false; |
---|
| 939 | pdfPomBPtr = 0; |
---|
| 940 | } |
---|
| 941 | |
---|
| 942 | // Set up the PDF's, if not already done. |
---|
| 943 | if (pdfAPtr == 0) { |
---|
| 944 | pdfAPtr = getPDFPtr(idA); |
---|
| 945 | if (pdfAPtr == 0 || !pdfAPtr->isSetup()) { |
---|
| 946 | info.errorMsg("Error in Pythia::init: " |
---|
| 947 | "could not set up PDF for beam A"); |
---|
| 948 | return false; |
---|
| 949 | } |
---|
| 950 | pdfHardAPtr = pdfAPtr; |
---|
| 951 | useNewPdfA = true; |
---|
| 952 | } |
---|
| 953 | if (pdfBPtr == 0) { |
---|
| 954 | pdfBPtr = getPDFPtr(idB); |
---|
| 955 | if (pdfBPtr == 0 || !pdfBPtr->isSetup()) { |
---|
| 956 | info.errorMsg("Error in Pythia::init: " |
---|
| 957 | "could not set up PDF for beam B"); |
---|
| 958 | return false; |
---|
| 959 | } |
---|
| 960 | pdfHardBPtr = pdfBPtr; |
---|
| 961 | useNewPdfB = true; |
---|
| 962 | } |
---|
| 963 | |
---|
| 964 | // Optionally set up separate PDF's for hard process. |
---|
| 965 | if (settings.flag("PDF:useHard") && useNewPdfA && useNewPdfB) { |
---|
| 966 | pdfHardAPtr = getPDFPtr(idA, 2); |
---|
| 967 | if (!pdfHardAPtr->isSetup()) return false; |
---|
| 968 | pdfHardBPtr = getPDFPtr(idB, 2); |
---|
| 969 | if (!pdfHardBPtr->isSetup()) return false; |
---|
| 970 | useNewPdfHard = true; |
---|
| 971 | } |
---|
| 972 | |
---|
| 973 | // Optionally set up Pomeron PDF's for diffractive physics. |
---|
| 974 | if ( doDiffraction) { |
---|
| 975 | if (pdfPomAPtr == 0) { |
---|
| 976 | pdfPomAPtr = getPDFPtr(990); |
---|
| 977 | useNewPdfPomA = true; |
---|
| 978 | } |
---|
| 979 | if (pdfPomBPtr == 0) { |
---|
| 980 | pdfPomBPtr = getPDFPtr(990); |
---|
| 981 | useNewPdfPomB = true; |
---|
| 982 | } |
---|
| 983 | } |
---|
| 984 | |
---|
| 985 | // Done. |
---|
| 986 | return true; |
---|
| 987 | |
---|
| 988 | } |
---|
| 989 | |
---|
| 990 | //-------------------------------------------------------------------------- |
---|
| 991 | |
---|
| 992 | // Main routine to generate the next event, using internal machinery. |
---|
| 993 | |
---|
| 994 | bool Pythia::next() { |
---|
| 995 | |
---|
| 996 | // Regularly print how many events have been generated. |
---|
| 997 | int nPrevious = info.getCounter(3); |
---|
| 998 | if (nCount > 0 && nPrevious > 0 && nPrevious%nCount == 0) |
---|
| 999 | cout << "\n Pythia::next(): " << nPrevious |
---|
| 1000 | << " events have been generated " << endl; |
---|
| 1001 | |
---|
| 1002 | // Set/reset info counters specific to each event. |
---|
| 1003 | info.addCounter(3); |
---|
| 1004 | for (int i = 10; i < 13; ++i) info.setCounter(i); |
---|
| 1005 | |
---|
| 1006 | // Simpler option when no hard process, i.e. mainly hadron level. |
---|
| 1007 | if (!doProcessLevel) { |
---|
| 1008 | |
---|
| 1009 | // Optionally fetch in resonance decays from LHA interface. |
---|
| 1010 | if (doLHA && !processLevel.nextLHAdec( event)) { |
---|
| 1011 | if (info.atEndOfFile()) info.errorMsg("Abort from Pythia::next:" |
---|
| 1012 | " reached end of Les Houches Events File"); |
---|
| 1013 | return false; |
---|
| 1014 | } |
---|
| 1015 | |
---|
| 1016 | // Reset info array (while event record contains data). |
---|
| 1017 | info.clear(); |
---|
| 1018 | |
---|
| 1019 | // Set correct energy for system. |
---|
| 1020 | Vec4 pSum = 0.; |
---|
| 1021 | for (int i = 1; i < event.size(); ++i) |
---|
| 1022 | if (event[i].isFinal()) pSum += event[i].p(); |
---|
| 1023 | event[0].p( pSum ); |
---|
| 1024 | event[0].m( pSum.mCalc() ); |
---|
| 1025 | |
---|
| 1026 | // Generate hadronization and decays. |
---|
| 1027 | bool status = forceHadronLevel(); |
---|
| 1028 | if (status) info.addCounter(4); |
---|
| 1029 | if (status && nPrevious < nShowEvt) event.list(showSaV, showMaD); |
---|
| 1030 | return status; |
---|
| 1031 | } |
---|
| 1032 | |
---|
| 1033 | // Reset arrays. |
---|
| 1034 | info.clear(); |
---|
| 1035 | process.clear(); |
---|
| 1036 | event.clear(); |
---|
| 1037 | partonSystems.clear(); |
---|
| 1038 | beamA.clear(); |
---|
| 1039 | beamB.clear(); |
---|
| 1040 | beamPomA.clear(); |
---|
| 1041 | beamPomB.clear(); |
---|
| 1042 | |
---|
| 1043 | // Can only generate event if initialization worked. |
---|
| 1044 | if (!isInit) { |
---|
| 1045 | info.errorMsg("Abort from Pythia::next: " |
---|
| 1046 | "not properly initialized so cannot generate events"); |
---|
| 1047 | return false; |
---|
| 1048 | } |
---|
| 1049 | |
---|
| 1050 | // Pick beam momentum spread and beam vertex. |
---|
| 1051 | if (doMomentumSpread || doVertexSpread) beamShapePtr->pick(); |
---|
| 1052 | |
---|
| 1053 | // Recalculate kinematics when beam momentum spread. |
---|
| 1054 | if (doMomentumSpread) nextKinematics(); |
---|
| 1055 | |
---|
| 1056 | // Outer loop over hard processes; only relevant for user-set vetoes. |
---|
| 1057 | for ( ; ; ) { |
---|
| 1058 | info.addCounter(10); |
---|
| 1059 | bool hasVetoed = false; |
---|
| 1060 | |
---|
| 1061 | // Provide the hard process that starts it off. Only one try. |
---|
| 1062 | info.clear(); |
---|
| 1063 | process.clear(); |
---|
| 1064 | |
---|
| 1065 | if ( !processLevel.next( process) ) { |
---|
| 1066 | if (doLHA && info.atEndOfFile()) info.errorMsg("Abort from " |
---|
| 1067 | "Pythia::next: reached end of Les Houches Events File"); |
---|
| 1068 | else info.errorMsg("Abort from Pythia::next: " |
---|
| 1069 | "processLevel failed; giving up"); |
---|
| 1070 | return false; |
---|
| 1071 | } |
---|
| 1072 | |
---|
| 1073 | info.addCounter(11); |
---|
| 1074 | |
---|
| 1075 | // Possibility for a user veto of the process-level event. |
---|
| 1076 | if (doVetoProcess) { |
---|
| 1077 | hasVetoed = userHooksPtr->doVetoProcessLevel( process); |
---|
| 1078 | if (hasVetoed) { |
---|
| 1079 | if (abortIfVeto) return false; |
---|
| 1080 | continue; |
---|
| 1081 | } |
---|
| 1082 | } |
---|
| 1083 | |
---|
| 1084 | // Possibility to perform CKKW-L merging on this event. |
---|
| 1085 | if(doMerging) { |
---|
| 1086 | hasVetoed = mergeProcess(); |
---|
| 1087 | if (hasVetoed) { |
---|
| 1088 | if (abortIfVeto) return false; |
---|
| 1089 | continue; |
---|
| 1090 | } |
---|
| 1091 | } |
---|
| 1092 | |
---|
| 1093 | // Possibility to stop the generation at this stage. |
---|
| 1094 | if (!doPartonLevel) { |
---|
| 1095 | boostAndVertex( true, true); |
---|
| 1096 | processLevel.accumulate(); |
---|
| 1097 | info.addCounter(4); |
---|
| 1098 | if (doLHA && nPrevious < nShowLHA) lhaUpPtr->listEvent(); |
---|
| 1099 | if (nPrevious < nShowInfo) info.list(); |
---|
| 1100 | if (nPrevious < nShowProc) process.list(showSaV, showMaD); |
---|
| 1101 | return true; |
---|
| 1102 | } |
---|
| 1103 | |
---|
| 1104 | // Save spare copy of process record in case of problems. |
---|
| 1105 | Event processSave = process; |
---|
| 1106 | int sizeMPI = info.sizeMPIarrays(); |
---|
| 1107 | info.addCounter(12); |
---|
| 1108 | for (int i = 14; i < 19; ++i) info.setCounter(i); |
---|
| 1109 | |
---|
| 1110 | // Allow up to ten tries for parton- and hadron-level processing. |
---|
| 1111 | bool physical = true; |
---|
| 1112 | for (int iTry = 0; iTry < NTRY; ++ iTry) { |
---|
| 1113 | info.addCounter(14); |
---|
| 1114 | physical = true; |
---|
| 1115 | |
---|
| 1116 | // Restore original process record if problems. |
---|
| 1117 | if (iTry > 0) process = processSave; |
---|
| 1118 | if (iTry > 0) info.resizeMPIarrays( sizeMPI); |
---|
| 1119 | |
---|
| 1120 | // Reset event record and (extracted partons from) beam remnants. |
---|
| 1121 | event.clear(); |
---|
| 1122 | beamA.clear(); |
---|
| 1123 | beamB.clear(); |
---|
| 1124 | beamPomA.clear(); |
---|
| 1125 | beamPomB.clear(); |
---|
| 1126 | partonSystems.clear(); |
---|
| 1127 | |
---|
| 1128 | // Parton-level evolution: ISR, FSR, MPI. |
---|
| 1129 | if ( !partonLevel.next( process, event) ) { |
---|
| 1130 | // Skip to next hard process for failure owing to deliberate veto. |
---|
| 1131 | hasVetoed = partonLevel.hasVetoed(); |
---|
| 1132 | if (hasVetoed) { |
---|
| 1133 | if (abortIfVeto) return false; |
---|
| 1134 | break; |
---|
| 1135 | } |
---|
| 1136 | // Else make a new try for other failures. |
---|
| 1137 | info.errorMsg("Error in Pythia::next: " |
---|
| 1138 | "partonLevel failed; try again"); |
---|
| 1139 | physical = false; |
---|
| 1140 | continue; |
---|
| 1141 | } |
---|
| 1142 | info.addCounter(15); |
---|
| 1143 | |
---|
| 1144 | // Possibility for a user veto of the parton-level event. |
---|
| 1145 | if (doVetoPartons) { |
---|
| 1146 | hasVetoed = userHooksPtr->doVetoPartonLevel( event); |
---|
| 1147 | if (hasVetoed) { |
---|
| 1148 | if (abortIfVeto) return false; |
---|
| 1149 | break; |
---|
| 1150 | } |
---|
| 1151 | } |
---|
| 1152 | |
---|
| 1153 | // Boost to lab frame (before decays, for vertices). |
---|
| 1154 | boostAndVertex( true, true); |
---|
| 1155 | |
---|
| 1156 | // Possibility to stop the generation at this stage. |
---|
| 1157 | if (!doHadronLevel) { |
---|
| 1158 | processLevel.accumulate(); |
---|
| 1159 | partonLevel.accumulate(); |
---|
| 1160 | |
---|
| 1161 | // Optionally check final event for problems. |
---|
| 1162 | if (checkEvent && !check()) { |
---|
| 1163 | info.errorMsg("Abort from Pythia::next: " |
---|
| 1164 | "check of event revealed problems"); |
---|
| 1165 | return false; |
---|
| 1166 | } |
---|
| 1167 | info.addCounter(4); |
---|
| 1168 | if (doLHA && nPrevious < nShowLHA) lhaUpPtr->listEvent(); |
---|
| 1169 | if (nPrevious < nShowInfo) info.list(); |
---|
| 1170 | if (nPrevious < nShowProc) process.list(showSaV, showMaD); |
---|
| 1171 | if (nPrevious < nShowEvt) event.list(showSaV, showMaD); |
---|
| 1172 | return true; |
---|
| 1173 | } |
---|
| 1174 | |
---|
| 1175 | // Hadron-level: hadronization, decays. |
---|
| 1176 | info.addCounter(16); |
---|
| 1177 | if ( !hadronLevel.next( event) ) { |
---|
| 1178 | info.errorMsg("Error in Pythia::next: " |
---|
| 1179 | "hadronLevel failed; try again"); |
---|
| 1180 | physical = false; |
---|
| 1181 | continue; |
---|
| 1182 | } |
---|
| 1183 | |
---|
| 1184 | // If R-hadrons have been formed, then (optionally) let them decay. |
---|
| 1185 | if (decayRHadrons && rHadrons.exist() && !doRHadronDecays()) { |
---|
| 1186 | info.errorMsg("Error in Pythia::next: " |
---|
| 1187 | "decayRHadrons failed; try again"); |
---|
| 1188 | physical = false; |
---|
| 1189 | continue; |
---|
| 1190 | } |
---|
| 1191 | info.addCounter(17); |
---|
| 1192 | |
---|
| 1193 | // Optionally check final event for problems. |
---|
| 1194 | if (checkEvent && !check()) { |
---|
| 1195 | info.errorMsg("Error in Pythia::next: " |
---|
| 1196 | "check of event revealed problems"); |
---|
| 1197 | physical = false; |
---|
| 1198 | continue; |
---|
| 1199 | } |
---|
| 1200 | |
---|
| 1201 | // Stop parton- and hadron-level looping if you got this far. |
---|
| 1202 | info.addCounter(18); |
---|
| 1203 | break; |
---|
| 1204 | } |
---|
| 1205 | |
---|
| 1206 | // If event vetoed then to make a new try. |
---|
| 1207 | if (hasVetoed) { |
---|
| 1208 | if (abortIfVeto) return false; |
---|
| 1209 | continue; |
---|
| 1210 | } |
---|
| 1211 | |
---|
| 1212 | // If event failed any other way (after ten tries) then give up. |
---|
| 1213 | if (!physical) { |
---|
| 1214 | info.errorMsg("Abort from Pythia::next: " |
---|
| 1215 | "parton+hadronLevel failed; giving up"); |
---|
| 1216 | return false; |
---|
| 1217 | } |
---|
| 1218 | |
---|
| 1219 | // Process- and parton-level statistics. Event scale. |
---|
| 1220 | processLevel.accumulate(); |
---|
| 1221 | partonLevel.accumulate(); |
---|
| 1222 | event.scale( process.scale() ); |
---|
| 1223 | |
---|
| 1224 | // End of outer loop over hard processes. Done with normal option. |
---|
| 1225 | info.addCounter(13); |
---|
| 1226 | break; |
---|
| 1227 | } |
---|
| 1228 | |
---|
| 1229 | // List events. |
---|
| 1230 | if (doLHA && nPrevious < nShowLHA) lhaUpPtr->listEvent(); |
---|
| 1231 | if (nPrevious < nShowInfo) info.list(); |
---|
| 1232 | if (nPrevious < nShowProc) process.list(showSaV,showMaD); |
---|
| 1233 | if (nPrevious < nShowEvt) event.list(showSaV, showMaD); |
---|
| 1234 | |
---|
| 1235 | // Done. |
---|
| 1236 | info.addCounter(4); |
---|
| 1237 | return true; |
---|
| 1238 | |
---|
| 1239 | } |
---|
| 1240 | |
---|
| 1241 | //-------------------------------------------------------------------------- |
---|
| 1242 | |
---|
| 1243 | // Generate only the hadronization/decay stage, using internal machinery. |
---|
| 1244 | // The "event" instance should already contain a parton-level configuration. |
---|
| 1245 | |
---|
| 1246 | bool Pythia::forceHadronLevel(bool findJunctions) { |
---|
| 1247 | |
---|
| 1248 | // Can only generate event if initialization worked. |
---|
| 1249 | if (!isInit) { |
---|
| 1250 | info.errorMsg("Abort from Pythia::forceHadronLevel: " |
---|
| 1251 | "not properly initialized so cannot generate events"); |
---|
| 1252 | return false; |
---|
| 1253 | } |
---|
| 1254 | |
---|
| 1255 | // Check whether any junctions in system. (Normally done in ProcessLevel.) |
---|
| 1256 | if (findJunctions) { |
---|
| 1257 | event.clearJunctions(); |
---|
| 1258 | processLevel.findJunctions( event); |
---|
| 1259 | } |
---|
| 1260 | |
---|
| 1261 | // Save spare copy of event in case of failure. |
---|
| 1262 | Event spareEvent = event; |
---|
| 1263 | |
---|
| 1264 | // Allow up to ten tries for hadron-level processing. |
---|
| 1265 | bool physical = true; |
---|
| 1266 | for (int iTry = 0; iTry < NTRY; ++ iTry) { |
---|
| 1267 | physical = true; |
---|
| 1268 | |
---|
| 1269 | // Check whether any resonances need to be handled at process level. |
---|
| 1270 | if (doResDec) { |
---|
| 1271 | process = event; |
---|
| 1272 | processLevel.nextDecays( process); |
---|
| 1273 | bool hasDecays = false; |
---|
| 1274 | for (int i = 1; i < process.size(); ++i) |
---|
| 1275 | if (process[i].status() < 0) hasDecays = true; |
---|
| 1276 | |
---|
| 1277 | // Allow for showers if decays happened at process level. |
---|
| 1278 | if (hasDecays) { |
---|
| 1279 | if (doFSRinRes) { |
---|
| 1280 | partonLevel.setupShowerSys( process, event); |
---|
| 1281 | partonLevel.resonanceShowers( process, event, false); |
---|
| 1282 | } else event = process; |
---|
| 1283 | } |
---|
| 1284 | } |
---|
| 1285 | |
---|
| 1286 | // Hadron-level: hadronization, decays. |
---|
| 1287 | if (hadronLevel.next( event)) break; |
---|
| 1288 | |
---|
| 1289 | // If failure then warn, restore original configuration and try again. |
---|
| 1290 | info.errorMsg("Error in Pythia::forceHadronLevel: " |
---|
| 1291 | "hadronLevel failed; try again"); |
---|
| 1292 | physical = false; |
---|
| 1293 | event = spareEvent; |
---|
| 1294 | } |
---|
| 1295 | |
---|
| 1296 | // Done for simpler option. |
---|
| 1297 | if (!physical) { |
---|
| 1298 | info.errorMsg("Abort from Pythia::forceHadronLevel: " |
---|
| 1299 | "hadronLevel failed; giving up"); |
---|
| 1300 | return false; |
---|
| 1301 | } |
---|
| 1302 | |
---|
| 1303 | // Optionally check final event for problems. |
---|
| 1304 | if (checkEvent && !check()) { |
---|
| 1305 | info.errorMsg("Abort from Pythia::forceHadronLevel: " |
---|
| 1306 | "check of event revealed problems"); |
---|
| 1307 | return false; |
---|
| 1308 | } |
---|
| 1309 | |
---|
| 1310 | // Done. |
---|
| 1311 | return true; |
---|
| 1312 | |
---|
| 1313 | } |
---|
| 1314 | |
---|
| 1315 | //-------------------------------------------------------------------------- |
---|
| 1316 | |
---|
| 1317 | // Recalculate kinematics for each event when beam momentum has a spread. |
---|
| 1318 | |
---|
| 1319 | void Pythia::nextKinematics() { |
---|
| 1320 | |
---|
| 1321 | // Read out momentum shift to give current beam momenta. |
---|
| 1322 | pAnow = pAinit + beamShapePtr->deltaPA(); |
---|
| 1323 | pAnow.e( sqrt(pAnow.pAbs2() + mA * mA) ); |
---|
| 1324 | pBnow = pBinit + beamShapePtr->deltaPB(); |
---|
| 1325 | pBnow.e( sqrt(pBnow.pAbs2() + mB * mB) ); |
---|
| 1326 | |
---|
| 1327 | // Construct CM frame kinematics. |
---|
| 1328 | eCM = (pAnow + pBnow).mCalc(); |
---|
| 1329 | pzAcm = 0.5 * sqrtpos( (eCM + mA + mB) * (eCM - mA - mB) |
---|
| 1330 | * (eCM - mA + mB) * (eCM + mA - mB) ) / eCM; |
---|
| 1331 | pzBcm = -pzAcm; |
---|
| 1332 | eA = sqrt(mA*mA + pzAcm*pzAcm); |
---|
| 1333 | eB = sqrt(mB*mB + pzBcm*pzBcm); |
---|
| 1334 | |
---|
| 1335 | // Set relevant info for other classes to use. |
---|
| 1336 | info.setBeamA( idA, pzAcm, eA, mA); |
---|
| 1337 | info.setBeamB( idB, pzBcm, eB, mB); |
---|
| 1338 | info.setECM( eCM); |
---|
| 1339 | beamA.newPzE( pzAcm, eA); |
---|
| 1340 | beamB.newPzE( pzBcm, eB); |
---|
| 1341 | |
---|
| 1342 | // Set boost/rotation matrices from/to CM frame. |
---|
| 1343 | MfromCM.reset(); |
---|
| 1344 | MfromCM.fromCMframe( pAnow, pBnow); |
---|
| 1345 | MtoCM = MfromCM; |
---|
| 1346 | MtoCM.invert(); |
---|
| 1347 | |
---|
| 1348 | } |
---|
| 1349 | |
---|
| 1350 | //-------------------------------------------------------------------------- |
---|
| 1351 | |
---|
| 1352 | // Boost from CM frame to lab frame, or inverse. Set production vertex. |
---|
| 1353 | |
---|
| 1354 | void Pythia::boostAndVertex( bool toLab, bool setVertex) { |
---|
| 1355 | |
---|
| 1356 | // Boost process from CM frame to lab frame. |
---|
| 1357 | if (toLab) { |
---|
| 1358 | if (boostType == 2) process.bst(0., 0., betaZ, gammaZ); |
---|
| 1359 | else if (boostType == 3) process.rotbst(MfromCM); |
---|
| 1360 | |
---|
| 1361 | // Boost nonempty event from CM frame to lab frame. |
---|
| 1362 | if (event.size() > 0) { |
---|
| 1363 | if (boostType == 2) event.bst(0., 0., betaZ, gammaZ); |
---|
| 1364 | else if (boostType == 3) event.rotbst(MfromCM); |
---|
| 1365 | } |
---|
| 1366 | |
---|
| 1367 | // Boost process from lab frame to CM frame. |
---|
| 1368 | } else { |
---|
| 1369 | if (boostType == 2) process.bst(0., 0., -betaZ, gammaZ); |
---|
| 1370 | else if (boostType == 3) process.rotbst(MtoCM); |
---|
| 1371 | |
---|
| 1372 | // Boost nonempty event from lab frame to CM frame. |
---|
| 1373 | if (event.size() > 0) { |
---|
| 1374 | if (boostType == 2) event.bst(0., 0., -betaZ, gammaZ); |
---|
| 1375 | else if (boostType == 3) event.rotbst(MtoCM); |
---|
| 1376 | } |
---|
| 1377 | } |
---|
| 1378 | |
---|
| 1379 | // Set production vertex; assumes particles are in lab frame and at origin. |
---|
| 1380 | if (setVertex && doVertexSpread) { |
---|
| 1381 | Vec4 vertex = beamShapePtr->vertex(); |
---|
| 1382 | for (int i = 0; i < process.size(); ++i) process[i].vProd( vertex); |
---|
| 1383 | for (int i = 0; i < event.size(); ++i) event[i].vProd( vertex); |
---|
| 1384 | } |
---|
| 1385 | |
---|
| 1386 | } |
---|
| 1387 | |
---|
| 1388 | //-------------------------------------------------------------------------- |
---|
| 1389 | |
---|
| 1390 | // Perform R-hadron decays, either as part of normal evolution or forced. |
---|
| 1391 | |
---|
| 1392 | bool Pythia::doRHadronDecays( ) { |
---|
| 1393 | |
---|
| 1394 | // Check if R-hadrons exist to be processed. |
---|
| 1395 | if ( !rHadrons.exist() ) return true; |
---|
| 1396 | |
---|
| 1397 | // Do the R-hadron decay itself. |
---|
| 1398 | if ( !rHadrons.decay( event) ) return false; |
---|
| 1399 | |
---|
| 1400 | // Perform showers in resonance decay chains. |
---|
| 1401 | if ( !partonLevel.resonanceShowers( process, event, false) ) return false; |
---|
| 1402 | |
---|
| 1403 | // Subsequent hadronization and decays. |
---|
| 1404 | if ( !hadronLevel.next( event) ) return false; |
---|
| 1405 | |
---|
| 1406 | // Done. |
---|
| 1407 | return true; |
---|
| 1408 | |
---|
| 1409 | } |
---|
| 1410 | |
---|
| 1411 | //-------------------------------------------------------------------------- |
---|
| 1412 | |
---|
| 1413 | // Print statistics on event generation. |
---|
| 1414 | |
---|
| 1415 | void Pythia::stat() { |
---|
| 1416 | |
---|
| 1417 | // Read out settings for what to include. |
---|
| 1418 | bool showPrL = settings.flag("Stat:showProcessLevel"); |
---|
| 1419 | bool showPaL = settings.flag("Stat:showPartonLevel"); |
---|
| 1420 | bool showErr = settings.flag("Stat:showErrors"); |
---|
| 1421 | bool reset = settings.flag("Stat:reset"); |
---|
| 1422 | |
---|
| 1423 | // Statistics on cross section and number of events. |
---|
| 1424 | if (doProcessLevel) { |
---|
| 1425 | if (showPrL) processLevel.statistics(false); |
---|
| 1426 | if (reset) processLevel.resetStatistics(); |
---|
| 1427 | } |
---|
| 1428 | |
---|
| 1429 | // Statistics from other classes, currently multiparton interactions. |
---|
| 1430 | if (showPaL) partonLevel.statistics(false); |
---|
| 1431 | if (reset) partonLevel.resetStatistics(); |
---|
| 1432 | |
---|
| 1433 | // Warning if significant part of events considerably above the |
---|
| 1434 | // merging scale value. |
---|
| 1435 | if (hasMergingHooks && doMerging && mergingHooksPtr->stats() ) { |
---|
| 1436 | string message="Warning in MergingHooks::stats: Most LH"; |
---|
| 1437 | message+=" Events significantly above Merging:TMS cut. Please check."; |
---|
| 1438 | info.errorMsg(message); |
---|
| 1439 | } |
---|
| 1440 | |
---|
| 1441 | // Summary of which and how many warnings/errors encountered. |
---|
| 1442 | if (showErr) info.errorStatistics(); |
---|
| 1443 | if (reset) info.errorReset(); |
---|
| 1444 | |
---|
| 1445 | } |
---|
| 1446 | |
---|
| 1447 | //-------------------------------------------------------------------------- |
---|
| 1448 | |
---|
| 1449 | // Print statistics on event generation - deprecated version. |
---|
| 1450 | |
---|
| 1451 | void Pythia::statistics(bool all, bool reset) { |
---|
| 1452 | |
---|
| 1453 | // Statistics on cross section and number of events. |
---|
| 1454 | if (doProcessLevel) processLevel.statistics(reset); |
---|
| 1455 | |
---|
| 1456 | // Statistics from other classes, e.g. multiparton interactions. |
---|
| 1457 | if (all) partonLevel.statistics(reset); |
---|
| 1458 | |
---|
| 1459 | // Warning if significant part of events considerably above the |
---|
| 1460 | // merging scale value. |
---|
| 1461 | if (hasMergingHooks && doMerging && mergingHooksPtr->stats() ) { |
---|
| 1462 | string message="Warning in MergingHooks::stats: Most LH"; |
---|
| 1463 | message+=" Events significantly above Merging:TMS cut. Please check."; |
---|
| 1464 | info.errorMsg(message); |
---|
| 1465 | } |
---|
| 1466 | |
---|
| 1467 | // Summary of which and how many warnings/errors encountered. |
---|
| 1468 | info.errorStatistics(); |
---|
| 1469 | if (reset) info.errorReset(); |
---|
| 1470 | |
---|
| 1471 | } |
---|
| 1472 | |
---|
| 1473 | //-------------------------------------------------------------------------- |
---|
| 1474 | |
---|
| 1475 | // Write the Pythia banner, with symbol and version information. |
---|
| 1476 | |
---|
| 1477 | void Pythia::banner(ostream& os) { |
---|
| 1478 | |
---|
| 1479 | // Read in version number and last date of change. |
---|
| 1480 | double versionNumber = settings.parm("Pythia:versionNumber"); |
---|
| 1481 | int versionDate = settings.mode("Pythia:versionDate"); |
---|
| 1482 | string month[12] = {"Jan", "Feb", "Mar", "Apr", "May", "Jun", |
---|
| 1483 | "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"}; |
---|
| 1484 | |
---|
| 1485 | // Get date and time. |
---|
| 1486 | time_t t = time(0); |
---|
| 1487 | char dateNow[12]; |
---|
| 1488 | strftime(dateNow,12,"%d %b %Y",localtime(&t)); |
---|
| 1489 | char timeNow[9]; |
---|
| 1490 | strftime(timeNow,9,"%H:%M:%S",localtime(&t)); |
---|
| 1491 | |
---|
| 1492 | os << "\n" |
---|
| 1493 | << " *-------------------------------------------" |
---|
| 1494 | << "-----------------------------------------* \n" |
---|
| 1495 | << " | " |
---|
| 1496 | << " | \n" |
---|
| 1497 | << " | *----------------------------------------" |
---|
| 1498 | << "--------------------------------------* | \n" |
---|
| 1499 | << " | | " |
---|
| 1500 | << " | | \n" |
---|
| 1501 | << " | | " |
---|
| 1502 | << " | | \n" |
---|
| 1503 | << " | | PPP Y Y TTTTT H H III A " |
---|
| 1504 | << " Welcome to the Lund Monte Carlo! | | \n" |
---|
| 1505 | << " | | P P Y Y T H H I A A " |
---|
| 1506 | << " This is PYTHIA version " << fixed << setprecision(3) |
---|
| 1507 | << setw(5) << versionNumber << " | | \n" |
---|
| 1508 | << " | | PPP Y T HHHHH I AAAAA" |
---|
| 1509 | << " Last date of change: " << setw(2) << versionDate%100 |
---|
| 1510 | << " " << month[ (versionDate/100)%100 - 1 ] |
---|
| 1511 | << " " << setw(4) << versionDate/10000 << " | | \n" |
---|
| 1512 | << " | | P Y T H H I A A" |
---|
| 1513 | << " | | \n" |
---|
| 1514 | << " | | P Y T H H III A A" |
---|
| 1515 | << " Now is " << dateNow << " at " << timeNow << " | | \n" |
---|
| 1516 | << " | | " |
---|
| 1517 | << " | | \n" |
---|
| 1518 | << " | | Torbjorn Sjostrand; Department of As" |
---|
| 1519 | << "tronomy and Theoretical Physics, | | \n" |
---|
| 1520 | << " | | Lund University, Solvegatan 14A, S" |
---|
| 1521 | << "E-223 62 Lund, Sweden; | | \n" |
---|
| 1522 | << " | | phone: + 46 - 46 - 222 48 16; e-ma" |
---|
| 1523 | << "il: torbjorn@thep.lu.se | | \n" |
---|
| 1524 | << " | | Stefan Ask; Department of Physics, U" |
---|
| 1525 | << "niversity of Cambridge, | | \n" |
---|
| 1526 | << " | | Cavendish Laboratory, JJ Thomson A" |
---|
| 1527 | << "ve., Cambridge CB3 0HE, UK; | | \n" |
---|
| 1528 | << " | | phone: + 41 - 22 - 767 6707; e-mai" |
---|
| 1529 | << "l: Stefan.Ask@cern.ch | | \n" |
---|
| 1530 | << " | | Richard Corke; Niels Bohr Institute," |
---|
| 1531 | << " University of Copenhagen, | | \n" |
---|
| 1532 | << " | | Blegdamsvej 17, 2100 Copenhagen, D" |
---|
| 1533 | << "enmark; | | \n" |
---|
| 1534 | << " | | e-mail: richard.corke@thep.lu.se " |
---|
| 1535 | << " | | \n" |
---|
| 1536 | << " | | Stephen Mrenna; Computing Division, " |
---|
| 1537 | << "Simulations Group, | | \n" |
---|
| 1538 | << " | | Fermi National Accelerator Laborat" |
---|
| 1539 | << "ory, MS 234, Batavia, IL 60510, USA; | | \n" |
---|
| 1540 | << " | | phone: + 1 - 630 - 840 - 2556; e-m" |
---|
| 1541 | << "ail: mrenna@fnal.gov | | \n" |
---|
| 1542 | << " | | Stefan Prestel; Department of Astron" |
---|
| 1543 | << "omy and Theoretical Physics, | | \n" |
---|
| 1544 | << " | | Lund University, Solvegatan 14A, S" |
---|
| 1545 | << "E-223 62 Lund, Sweden; | | \n" |
---|
| 1546 | << " | | phone: + 46 - 46 - 222 06 64; e-ma" |
---|
| 1547 | << "il: stefan.prestel@thep.lu.se | | \n" |
---|
| 1548 | << " | | Peter Skands; Theoretical Physics, C" |
---|
| 1549 | << "ERN, CH-1211 Geneva 23, Switzerland; | | \n" |
---|
| 1550 | << " | | phone: + 41 - 22 - 767 2447; e-mai" |
---|
| 1551 | << "l: peter.skands@cern.ch | | \n" |
---|
| 1552 | << " | | " |
---|
| 1553 | << " | | \n" |
---|
| 1554 | << " | | The main program reference is the 'Br" |
---|
| 1555 | << "ief Introduction to PYTHIA 8.1', | | \n" |
---|
| 1556 | << " | | T. Sjostrand, S. Mrenna and P. Skands" |
---|
| 1557 | << ", Comput. Phys. Comm. 178 (2008) 85 | | \n" |
---|
| 1558 | << " | | [arXiv:0710.3820] " |
---|
| 1559 | << " | | \n" |
---|
| 1560 | << " | | " |
---|
| 1561 | << " | | \n" |
---|
| 1562 | << " | | The main physics reference is the 'PY" |
---|
| 1563 | << "THIA 6.4 Physics and Manual', | | \n" |
---|
| 1564 | << " | | T. Sjostrand, S. Mrenna and P. Skands" |
---|
| 1565 | << ", JHEP05 (2006) 026 [hep-ph/0603175]. | | \n" |
---|
| 1566 | << " | | " |
---|
| 1567 | << " | | \n" |
---|
| 1568 | << " | | An archive of program versions and do" |
---|
| 1569 | << "cumentation is found on the web: | | \n" |
---|
| 1570 | << " | | http://www.thep.lu.se/~torbjorn/Pythi" |
---|
| 1571 | << "a.html | | \n" |
---|
| 1572 | << " | | " |
---|
| 1573 | << " | | \n" |
---|
| 1574 | << " | | This program is released under the GN" |
---|
| 1575 | << "U General Public Licence version 2. | | \n" |
---|
| 1576 | << " | | Please respect the MCnet Guidelines f" |
---|
| 1577 | << "or Event Generator Authors and Users. | | \n" |
---|
| 1578 | << " | | " |
---|
| 1579 | << " | | \n" |
---|
| 1580 | << " | | Disclaimer: this program comes withou" |
---|
| 1581 | << "t any guarantees. | | \n" |
---|
| 1582 | << " | | Beware of errors and use common sense" |
---|
| 1583 | << " when interpreting results. | | \n" |
---|
| 1584 | << " | | " |
---|
| 1585 | << " | | \n" |
---|
| 1586 | << " | | Copyright (C) 2012 Torbjorn Sjostrand" |
---|
| 1587 | << " | | \n" |
---|
| 1588 | << " | | " |
---|
| 1589 | << " | | \n" |
---|
| 1590 | << " | | " |
---|
| 1591 | << " | | \n" |
---|
| 1592 | << " | *----------------------------------------" |
---|
| 1593 | << "--------------------------------------* | \n" |
---|
| 1594 | << " | " |
---|
| 1595 | << " | \n" |
---|
| 1596 | << " *-------------------------------------------" |
---|
| 1597 | << "-----------------------------------------* \n" << endl; |
---|
| 1598 | |
---|
| 1599 | } |
---|
| 1600 | |
---|
| 1601 | //-------------------------------------------------------------------------- |
---|
| 1602 | |
---|
| 1603 | // Check for lines in file that mark the beginning of new subrun. |
---|
| 1604 | |
---|
| 1605 | int Pythia::readSubrun(string line, bool warn, ostream& os) { |
---|
| 1606 | |
---|
| 1607 | // If empty line then done. |
---|
| 1608 | int subrunLine = SUBRUNDEFAULT; |
---|
| 1609 | if (line.find_first_not_of(" \n\t\v\b\r\f\a") == string::npos) |
---|
| 1610 | return subrunLine; |
---|
| 1611 | |
---|
| 1612 | // If first character is not a letter, then done. |
---|
| 1613 | string lineNow = line; |
---|
| 1614 | int firstChar = lineNow.find_first_not_of(" \n\t\v\b\r\f\a"); |
---|
| 1615 | if (!isalpha(lineNow[firstChar])) return subrunLine; |
---|
| 1616 | |
---|
| 1617 | // Replace an equal sign by a blank to make parsing simpler. |
---|
| 1618 | while (lineNow.find("=") != string::npos) { |
---|
| 1619 | int firstEqual = lineNow.find_first_of("="); |
---|
| 1620 | lineNow.replace(firstEqual, 1, " "); |
---|
| 1621 | } |
---|
| 1622 | |
---|
| 1623 | // Get first word of a line. |
---|
| 1624 | istringstream splitLine(lineNow); |
---|
| 1625 | string name; |
---|
| 1626 | splitLine >> name; |
---|
| 1627 | |
---|
| 1628 | // Replace two colons by one (:: -> :) to allow for such mistakes. |
---|
| 1629 | while (name.find("::") != string::npos) { |
---|
| 1630 | int firstColonColon = name.find_first_of("::"); |
---|
| 1631 | name.replace(firstColonColon, 2, ":"); |
---|
| 1632 | } |
---|
| 1633 | |
---|
| 1634 | // Convert to lowercase. |
---|
| 1635 | for (int i = 0; i < int(name.length()); ++i) |
---|
| 1636 | name[i] = std::tolower(name[i]); |
---|
| 1637 | |
---|
| 1638 | // If no match then done. |
---|
| 1639 | if (name != "main:subrun") return subrunLine; |
---|
| 1640 | |
---|
| 1641 | // Else find new subrun number and return it. |
---|
| 1642 | splitLine >> subrunLine; |
---|
| 1643 | if (!splitLine) { |
---|
| 1644 | if (warn) os << "\n PYTHIA Warning: Main:subrun number not" |
---|
| 1645 | << " recognized; skip:\n " << line << endl; |
---|
| 1646 | subrunLine = SUBRUNDEFAULT; |
---|
| 1647 | } |
---|
| 1648 | return subrunLine; |
---|
| 1649 | |
---|
| 1650 | } |
---|
| 1651 | |
---|
| 1652 | //-------------------------------------------------------------------------- |
---|
| 1653 | |
---|
| 1654 | // Check that the final event makes sense: no unknown id codes; |
---|
| 1655 | // charge and energy-momentum conserved. |
---|
| 1656 | |
---|
| 1657 | bool Pythia::check(ostream& os) { |
---|
| 1658 | |
---|
| 1659 | // Reset. |
---|
| 1660 | bool physical = true; |
---|
| 1661 | bool listVertices = false; |
---|
| 1662 | bool listHistory = false; |
---|
| 1663 | bool listSystems = false; |
---|
| 1664 | bool listBeams = false; |
---|
| 1665 | iErrId.resize(0); |
---|
| 1666 | iErrCol.resize(0); |
---|
| 1667 | iErrNan.resize(0); |
---|
| 1668 | iErrNanVtx.resize(0); |
---|
| 1669 | Vec4 pSum; |
---|
| 1670 | double chargeSum = 0.; |
---|
| 1671 | |
---|
| 1672 | // Incoming beams counted with negative momentum and charge. |
---|
| 1673 | if (doProcessLevel) { |
---|
| 1674 | pSum = - (event[1].p() + event[2].p()); |
---|
| 1675 | chargeSum = - (event[1].charge() + event[2].charge()); |
---|
| 1676 | |
---|
| 1677 | // If no ProcessLevel then sum final state of process record. |
---|
| 1678 | } else if (process.size() > 0) { |
---|
| 1679 | pSum = - process[0].p(); |
---|
| 1680 | for (int i = 0; i < process.size(); ++i) |
---|
| 1681 | if (process[i].isFinal()) chargeSum -= process[i].charge(); |
---|
| 1682 | |
---|
| 1683 | // If process not filled, then use outgoing primary in event. |
---|
| 1684 | } else { |
---|
| 1685 | pSum = - event[0].p(); |
---|
| 1686 | for (int i = 1; i < event.size(); ++i) |
---|
| 1687 | if (event[i].statusAbs() < 10 || event[i].statusAbs() == 23) |
---|
| 1688 | chargeSum -= event[i].charge(); |
---|
| 1689 | } |
---|
| 1690 | double eLab = abs(pSum.e()); |
---|
| 1691 | |
---|
| 1692 | // Loop over particles in the event. |
---|
| 1693 | for (int i = 0; i < event.size(); ++i) { |
---|
| 1694 | |
---|
| 1695 | // Look for any unrecognized particle codes. |
---|
| 1696 | int id = event[i].id(); |
---|
| 1697 | if (id == 0 || !particleData.isParticle(id)) { |
---|
| 1698 | ostringstream errCode; |
---|
| 1699 | errCode << ", i = " << i << ", id = " << id; |
---|
| 1700 | info.errorMsg("Error in Pythia::check: " |
---|
| 1701 | "unknown particle code", errCode.str()); |
---|
| 1702 | physical = false; |
---|
| 1703 | iErrId.push_back(i); |
---|
| 1704 | |
---|
| 1705 | // Check that colour assignments are the expected ones. |
---|
| 1706 | } else { |
---|
| 1707 | int colType = event[i].colType(); |
---|
| 1708 | int col = event[i].col(); |
---|
| 1709 | int acol = event[i].acol(); |
---|
| 1710 | if ( (colType == 0 && (col > 0 || acol > 0)) |
---|
| 1711 | || (colType == 1 && (col <= 0 || acol > 0)) |
---|
| 1712 | || (colType == -1 && (col > 0 || acol <= 0)) |
---|
| 1713 | || (colType == 2 && (col <= 0 || acol <= 0)) ) { |
---|
| 1714 | ostringstream errCode; |
---|
| 1715 | errCode << ", i = " << i << ", id = " << id << " cols = " << col |
---|
| 1716 | << " " << acol; |
---|
| 1717 | info.errorMsg("Error in Pythia::check: " |
---|
| 1718 | "incorrect colours", errCode.str()); |
---|
| 1719 | physical = false; |
---|
| 1720 | iErrCol.push_back(i); |
---|
| 1721 | } |
---|
| 1722 | } |
---|
| 1723 | |
---|
| 1724 | // Look for particles with not-a-number energy/momentum/mass. |
---|
| 1725 | if (abs(event[i].px()) >= 0. && abs(event[i].py()) >= 0. |
---|
| 1726 | && abs(event[i].pz()) >= 0. && abs(event[i].e()) >= 0. |
---|
| 1727 | && abs(event[i].m()) >= 0.) ; |
---|
| 1728 | else { |
---|
| 1729 | info.errorMsg("Error in Pythia::check: " |
---|
| 1730 | "not-a-number energy/momentum/mass"); |
---|
| 1731 | physical = false; |
---|
| 1732 | iErrNan.push_back(i); |
---|
| 1733 | } |
---|
| 1734 | |
---|
| 1735 | // Look for particles with not-a-number vertex/lifetime. |
---|
| 1736 | if (abs(event[i].xProd()) >= 0. && abs(event[i].yProd()) >= 0. |
---|
| 1737 | && abs(event[i].zProd()) >= 0. && abs(event[i].tProd()) >= 0. |
---|
| 1738 | && abs(event[i].tau()) >= 0.) ; |
---|
| 1739 | else { |
---|
| 1740 | info.errorMsg("Error in Pythia::check: " |
---|
| 1741 | "not-a-number vertex/lifetime"); |
---|
| 1742 | physical = false; |
---|
| 1743 | listVertices = true; |
---|
| 1744 | iErrNanVtx.push_back(i); |
---|
| 1745 | } |
---|
| 1746 | |
---|
| 1747 | // Add final-state four-momentum and charge. |
---|
| 1748 | if (event[i].isFinal()) { |
---|
| 1749 | pSum += event[i].p(); |
---|
| 1750 | chargeSum += event[i].charge(); |
---|
| 1751 | } |
---|
| 1752 | |
---|
| 1753 | // End of particle loop. |
---|
| 1754 | } |
---|
| 1755 | |
---|
| 1756 | // Check energy-momentum/charge conservation. |
---|
| 1757 | double epDev = abs(pSum.e()) + abs(pSum.px()) + abs(pSum.py()) |
---|
| 1758 | + abs(pSum.pz()); |
---|
| 1759 | if (epDev > epTolErr * eLab) { |
---|
| 1760 | info.errorMsg("Error in Pythia::check: energy-momentum not conserved"); |
---|
| 1761 | physical = false; |
---|
| 1762 | } else if (epDev > epTolWarn * eLab) { |
---|
| 1763 | info.errorMsg("Warning in Pythia::check: " |
---|
| 1764 | "energy-momentum not quite conserved"); |
---|
| 1765 | } |
---|
| 1766 | if (abs(chargeSum) > 0.1) { |
---|
| 1767 | info.errorMsg("Error in Pythia::check: charge not conserved"); |
---|
| 1768 | physical = false; |
---|
| 1769 | } |
---|
| 1770 | |
---|
| 1771 | // Check that beams and event records agree on incoming partons. |
---|
| 1772 | // Only meaningful for resolved beams. |
---|
| 1773 | if (info.isResolved()) |
---|
| 1774 | for (int iSys = 0; iSys < beamA.sizeInit(); ++iSys) { |
---|
| 1775 | int eventANw = partonSystems.getInA(iSys); |
---|
| 1776 | int eventBNw = partonSystems.getInB(iSys); |
---|
| 1777 | int beamANw = beamA[iSys].iPos(); |
---|
| 1778 | int beamBNw = beamB[iSys].iPos(); |
---|
| 1779 | if (eventANw != beamANw || eventBNw != beamBNw) { |
---|
| 1780 | info.errorMsg("Error in Pythia::check: " |
---|
| 1781 | "event and beams records disagree"); |
---|
| 1782 | physical = false; |
---|
| 1783 | listSystems = true; |
---|
| 1784 | listBeams = true; |
---|
| 1785 | } |
---|
| 1786 | } |
---|
| 1787 | |
---|
| 1788 | // Check that mother and daughter information match for each particle. |
---|
| 1789 | vector<int> noMot; |
---|
| 1790 | vector<int> noDau; |
---|
| 1791 | vector< pair<int,int> > noMotDau; |
---|
| 1792 | if (checkHistory) { |
---|
| 1793 | |
---|
| 1794 | // Loop through the event and check that there are beam particles. |
---|
| 1795 | bool hasBeams = false; |
---|
| 1796 | for (int i = 0; i < event.size(); ++i) { |
---|
| 1797 | int status = event[i].status(); |
---|
| 1798 | if (abs(status) == 12) hasBeams = true; |
---|
| 1799 | |
---|
| 1800 | // Check that mother and daughter lists not empty where not expected to. |
---|
| 1801 | vector<int> mList = event.motherList(i); |
---|
| 1802 | vector<int> dList = event.daughterList(i); |
---|
| 1803 | if (mList.size() == 0 && abs(status) != 11 && abs(status) != 12) |
---|
| 1804 | noMot.push_back(i); |
---|
| 1805 | if (dList.size() == 0 && status < 0 && status != -11) |
---|
| 1806 | noDau.push_back(i); |
---|
| 1807 | |
---|
| 1808 | // Check that the particle appears in the daughters list of each mother. |
---|
| 1809 | for (int j = 0; j < int(mList.size()); ++j) { |
---|
| 1810 | vector<int> dmList = event.daughterList( mList[j] ); |
---|
| 1811 | bool foundMatch = false; |
---|
| 1812 | for (int k = 0; k < int(dmList.size()); ++k) |
---|
| 1813 | if (dmList[k] == i) { |
---|
| 1814 | foundMatch = true; |
---|
| 1815 | break; |
---|
| 1816 | } |
---|
| 1817 | if (!hasBeams && mList.size() == 1 && mList[0] == 0) foundMatch = true; |
---|
| 1818 | if (!foundMatch) { |
---|
| 1819 | bool oldPair = false; |
---|
| 1820 | for (int k = 0; k < int(noMotDau.size()); ++k) |
---|
| 1821 | if (noMotDau[k].first == mList[j] && noMotDau[k].second == i) { |
---|
| 1822 | oldPair = true; |
---|
| 1823 | break; |
---|
| 1824 | } |
---|
| 1825 | if (!oldPair) noMotDau.push_back( make_pair( mList[j], i) ); |
---|
| 1826 | } |
---|
| 1827 | } |
---|
| 1828 | |
---|
| 1829 | // Check that the particle appears in the mothers list of each daughter. |
---|
| 1830 | for (int j = 0; j < int(dList.size()); ++j) { |
---|
| 1831 | vector<int> mdList = event.motherList( dList[j] ); |
---|
| 1832 | bool foundMatch = false; |
---|
| 1833 | for (int k = 0; k < int(mdList.size()); ++k) |
---|
| 1834 | if (mdList[k] == i) { |
---|
| 1835 | foundMatch = true; |
---|
| 1836 | break; |
---|
| 1837 | } |
---|
| 1838 | if (!foundMatch) { |
---|
| 1839 | bool oldPair = false; |
---|
| 1840 | for (int k = 0; k < int(noMotDau.size()); ++k) |
---|
| 1841 | if (noMotDau[k].first == i && noMotDau[k].second == dList[j]) { |
---|
| 1842 | oldPair = true; |
---|
| 1843 | break; |
---|
| 1844 | } |
---|
| 1845 | if (!oldPair) noMotDau.push_back( make_pair( i, dList[j]) ); |
---|
| 1846 | } |
---|
| 1847 | } |
---|
| 1848 | } |
---|
| 1849 | |
---|
| 1850 | // Warn if any errors were found. |
---|
| 1851 | if (noMot.size() > 0 || noDau.size() > 0 || noMotDau.size() > 0) { |
---|
| 1852 | info.errorMsg("Error in Pythia::check: " |
---|
| 1853 | "mismatch in daughter and mother lists"); |
---|
| 1854 | physical = false; |
---|
| 1855 | listHistory = true; |
---|
| 1856 | } |
---|
| 1857 | } |
---|
| 1858 | |
---|
| 1859 | // Done for sensible events. |
---|
| 1860 | if (physical) return true; |
---|
| 1861 | |
---|
| 1862 | // Print (the first few) flawed events: local info. |
---|
| 1863 | if (nErrEvent < nErrList) { |
---|
| 1864 | os << "\n PYTHIA erroneous event info: \n"; |
---|
| 1865 | if (iErrId.size() > 0) { |
---|
| 1866 | os << " unknown particle codes in lines "; |
---|
| 1867 | for (int i = 0; i < int(iErrId.size()); ++i) |
---|
| 1868 | os << iErrId[i] << " "; |
---|
| 1869 | os << "\n"; |
---|
| 1870 | } |
---|
| 1871 | if (iErrCol.size() > 0) { |
---|
| 1872 | os << " incorrect colour assignments in lines "; |
---|
| 1873 | for (int i = 0; i < int(iErrCol.size()); ++i) |
---|
| 1874 | os << iErrCol[i] << " "; |
---|
| 1875 | os << "\n"; |
---|
| 1876 | } |
---|
| 1877 | if (iErrNan.size() > 0) { |
---|
| 1878 | os << " not-a-number energy/momentum/mass in lines "; |
---|
| 1879 | for (int i = 0; i < int(iErrNan.size()); ++i) |
---|
| 1880 | os << iErrNan[i] << " "; |
---|
| 1881 | os << "\n"; |
---|
| 1882 | } |
---|
| 1883 | if (iErrNanVtx.size() > 0) { |
---|
| 1884 | os << " not-a-number vertex/lifetime in lines "; |
---|
| 1885 | for (int i = 0; i < int(iErrNanVtx.size()); ++i) |
---|
| 1886 | os << iErrNanVtx[i] << " "; |
---|
| 1887 | os << "\n"; |
---|
| 1888 | } |
---|
| 1889 | if (epDev > epTolErr * eLab) os << scientific << setprecision(3) |
---|
| 1890 | << " total energy-momentum non-conservation = " << epDev << "\n"; |
---|
| 1891 | if (abs(chargeSum) > 0.1) os << fixed << setprecision(2) |
---|
| 1892 | << " total charge non-conservation = " << chargeSum << "\n"; |
---|
| 1893 | if (noMot.size() > 0) { |
---|
| 1894 | os << " missing mothers for particles "; |
---|
| 1895 | for (int i = 0; i < int(noMot.size()); ++i) os << noMot[i] << " "; |
---|
| 1896 | os << "\n"; |
---|
| 1897 | } |
---|
| 1898 | if (noDau.size() > 0) { |
---|
| 1899 | os << " missing daughters for particles "; |
---|
| 1900 | for (int i = 0; i < int(noDau.size()); ++i) os << noDau[i] << " "; |
---|
| 1901 | os << "\n"; |
---|
| 1902 | } |
---|
| 1903 | if (noMotDau.size() > 0) { |
---|
| 1904 | os << " inconsistent history for (mother,daughter) pairs "; |
---|
| 1905 | for (int i = 0; i < int(noMotDau.size()); ++i) |
---|
| 1906 | os << "(" << noMotDau[i].first << "," << noMotDau[i].second << ") "; |
---|
| 1907 | os << "\n"; |
---|
| 1908 | } |
---|
| 1909 | |
---|
| 1910 | // Print (the first few) flawed events: standard listings. |
---|
| 1911 | info.list(); |
---|
| 1912 | event.list(listVertices, listHistory); |
---|
| 1913 | if (listSystems) partonSystems.list(); |
---|
| 1914 | if (listBeams) beamA.list(); |
---|
| 1915 | if (listBeams) beamB.list(); |
---|
| 1916 | } |
---|
| 1917 | |
---|
| 1918 | // Update error counter. Done also for flawed event. |
---|
| 1919 | ++nErrEvent; |
---|
| 1920 | return false; |
---|
| 1921 | |
---|
| 1922 | } |
---|
| 1923 | |
---|
| 1924 | //-------------------------------------------------------------------------- |
---|
| 1925 | |
---|
| 1926 | // Routine to set up a PDF pointer. |
---|
| 1927 | |
---|
| 1928 | PDF* Pythia::getPDFPtr(int idIn, int sequence) { |
---|
| 1929 | |
---|
| 1930 | // Temporary pointer to be returned. |
---|
| 1931 | PDF* tempPDFPtr = 0; |
---|
| 1932 | |
---|
| 1933 | // One option is to treat a Pomeron like a pi0. |
---|
| 1934 | if (idIn == 990 && settings.mode("PDF:PomSet") == 2) idIn = 111; |
---|
| 1935 | |
---|
| 1936 | // Proton beam, normal choice. |
---|
| 1937 | if (abs(idIn) == 2212 && sequence == 1) { |
---|
| 1938 | int pSet = settings.mode("PDF:pSet"); |
---|
| 1939 | bool useLHAPDF = settings.flag("PDF:useLHAPDF"); |
---|
| 1940 | |
---|
| 1941 | // Use internal sets. |
---|
| 1942 | if (!useLHAPDF) { |
---|
| 1943 | if (pSet == 1) tempPDFPtr = new GRV94L(idIn); |
---|
| 1944 | else if (pSet == 2) tempPDFPtr = new CTEQ5L(idIn); |
---|
| 1945 | else if (pSet <= 6) |
---|
| 1946 | tempPDFPtr = new MSTWpdf(idIn, pSet - 2, xmlPath, &info); |
---|
| 1947 | else tempPDFPtr = new CTEQ6pdf(idIn, pSet - 6, xmlPath, &info); |
---|
| 1948 | } |
---|
| 1949 | |
---|
| 1950 | // Use sets from LHAPDF. |
---|
| 1951 | else { |
---|
| 1952 | string LHAPDFset = settings.word("PDF:LHAPDFset"); |
---|
| 1953 | int LHAPDFmember = settings.mode("PDF:LHAPDFmember"); |
---|
| 1954 | tempPDFPtr = new LHAPDF(idIn, LHAPDFset, LHAPDFmember, 1, &info); |
---|
| 1955 | |
---|
| 1956 | // Optionally allow extrapolation beyond x and Q2 limits. |
---|
| 1957 | tempPDFPtr->setExtrapolate( settings.flag("PDF:extrapolateLHAPDF") ); |
---|
| 1958 | } |
---|
| 1959 | } |
---|
| 1960 | |
---|
| 1961 | // Proton beam, special choice for the hard process. |
---|
| 1962 | else if (abs(idIn) == 2212) { |
---|
| 1963 | int pSet = settings.mode("PDF:pHardSet"); |
---|
| 1964 | bool useLHAPDF = settings.flag("PDF:useHardLHAPDF"); |
---|
| 1965 | |
---|
| 1966 | // Use internal sets. |
---|
| 1967 | if (!useLHAPDF) { |
---|
| 1968 | if (pSet == 1) tempPDFPtr = new GRV94L(idIn); |
---|
| 1969 | else if (pSet == 2) tempPDFPtr = new CTEQ5L(idIn); |
---|
| 1970 | else if (pSet <= 6) |
---|
| 1971 | tempPDFPtr = new MSTWpdf(idIn, pSet - 2, xmlPath, &info); |
---|
| 1972 | else tempPDFPtr = new CTEQ6pdf(idIn, pSet - 6, xmlPath, &info); |
---|
| 1973 | } |
---|
| 1974 | |
---|
| 1975 | // Use sets from LHAPDF. |
---|
| 1976 | else { |
---|
| 1977 | string LHAPDFset = settings.word("PDF:hardLHAPDFset"); |
---|
| 1978 | int LHAPDFmember = settings.mode("PDF:hardLHAPDFmember"); |
---|
| 1979 | // May need to require LHAPDF to handle two sets simultaneously. |
---|
| 1980 | int nSet = (settings.flag("PDF:useLHAPDF")) ? 2 : 1; |
---|
| 1981 | tempPDFPtr = new LHAPDF(idIn, LHAPDFset, LHAPDFmember, nSet, &info); |
---|
| 1982 | |
---|
| 1983 | // Optionally allow extrapolation beyond x and Q2 limits. |
---|
| 1984 | tempPDFPtr->setExtrapolate( settings.flag("PDF:extrapolateLHAPDF") ); |
---|
| 1985 | } |
---|
| 1986 | } |
---|
| 1987 | |
---|
| 1988 | // Pion beam (or, in one option, Pomeron beam). |
---|
| 1989 | else if (abs(idIn) == 211 || idIn == 111) { |
---|
| 1990 | bool useLHAPDF = settings.flag("PDF:piUseLHAPDF"); |
---|
| 1991 | |
---|
| 1992 | // Use internal sets. |
---|
| 1993 | if (!useLHAPDF) { |
---|
| 1994 | tempPDFPtr = new GRVpiL(idIn); |
---|
| 1995 | } |
---|
| 1996 | |
---|
| 1997 | // Use sets from LHAPDF. |
---|
| 1998 | else { |
---|
| 1999 | string LHAPDFset = settings.word("PDF:piLHAPDFset"); |
---|
| 2000 | int LHAPDFmember = settings.mode("PDF:piLHAPDFmember"); |
---|
| 2001 | tempPDFPtr = new LHAPDF(idIn, LHAPDFset, LHAPDFmember, 1, &info); |
---|
| 2002 | |
---|
| 2003 | // Optionally allow extrapolation beyond x and Q2 limits. |
---|
| 2004 | tempPDFPtr->setExtrapolate( settings.flag("PDF:extrapolateLHAPDF") ); |
---|
| 2005 | } |
---|
| 2006 | } |
---|
| 2007 | |
---|
| 2008 | // Pomeron beam, if not treated like a pi0 beam. |
---|
| 2009 | else if (idIn == 990) { |
---|
| 2010 | int pomSet = settings.mode("PDF:PomSet"); |
---|
| 2011 | double rescale = settings.parm("PDF:PomRescale"); |
---|
| 2012 | |
---|
| 2013 | // A generic Q2-independent parametrization. |
---|
| 2014 | if (pomSet == 1) { |
---|
| 2015 | double gluonA = settings.parm("PDF:PomGluonA"); |
---|
| 2016 | double gluonB = settings.parm("PDF:PomGluonB"); |
---|
| 2017 | double quarkA = settings.parm("PDF:PomQuarkA"); |
---|
| 2018 | double quarkB = settings.parm("PDF:PomQuarkB"); |
---|
| 2019 | double quarkFrac = settings.parm("PDF:PomQuarkFrac"); |
---|
| 2020 | double strangeSupp = settings.parm("PDF:PomStrangeSupp"); |
---|
| 2021 | tempPDFPtr = new PomFix( 990, gluonA, gluonB, quarkA, quarkB, |
---|
| 2022 | quarkFrac, strangeSupp); |
---|
| 2023 | } |
---|
| 2024 | |
---|
| 2025 | // The H1 Q2-dependent parametrizations. Initialization requires files. |
---|
| 2026 | else if (pomSet == 3 || pomSet == 4) |
---|
| 2027 | tempPDFPtr = new PomH1FitAB( 990, pomSet - 2, rescale, xmlPath, &info); |
---|
| 2028 | else if (pomSet == 5) |
---|
| 2029 | tempPDFPtr = new PomH1Jets( 990, rescale, xmlPath, &info); |
---|
| 2030 | else if (pomSet == 6) |
---|
| 2031 | tempPDFPtr = new PomH1FitAB( 990, 3, rescale, xmlPath, &info); |
---|
| 2032 | } |
---|
| 2033 | |
---|
| 2034 | // Lepton beam: neutrino, resolved charged lepton or unresolved ditto. |
---|
| 2035 | else if (abs(idIn) > 10 && abs(idIn) < 17) { |
---|
| 2036 | if (abs(idIn)%2 == 0) tempPDFPtr = new NeutrinoPoint(idIn); |
---|
| 2037 | else if (settings.flag("PDF:lepton")) tempPDFPtr = new Lepton(idIn); |
---|
| 2038 | else tempPDFPtr = new LeptonPoint(idIn); |
---|
| 2039 | } |
---|
| 2040 | |
---|
| 2041 | // Failure for unrecognized particle. |
---|
| 2042 | else tempPDFPtr = 0; |
---|
| 2043 | |
---|
| 2044 | // Done. |
---|
| 2045 | return tempPDFPtr; |
---|
| 2046 | } |
---|
| 2047 | |
---|
| 2048 | //-------------------------------------------------------------------------- |
---|
| 2049 | |
---|
| 2050 | // Initialize SUSY Les Houches Accord data. |
---|
| 2051 | |
---|
| 2052 | bool Pythia::initSLHA() { |
---|
| 2053 | |
---|
| 2054 | // Initial and settings values. |
---|
| 2055 | int ifailLHE = 1; |
---|
| 2056 | int ifailSpc = 1; |
---|
| 2057 | int readFrom = settings.mode("SLHA:readFrom"); |
---|
| 2058 | string lhefFile = settings.word("Beams:LHEF"); |
---|
| 2059 | string lhefHeader = settings.word("Beams:LHEFheader"); |
---|
| 2060 | string slhaFile = settings.word("SLHA:file"); |
---|
| 2061 | int verboseSLHA = settings.mode("SLHA:verbose"); |
---|
| 2062 | bool slhaUseDec = settings.flag("SLHA:useDecayTable"); |
---|
| 2063 | |
---|
| 2064 | // No SUSY by default |
---|
| 2065 | couplingsPtr->isSUSY = false; |
---|
| 2066 | |
---|
| 2067 | // Option with no SLHA read-in at all. |
---|
| 2068 | if (readFrom == 0) return true; |
---|
| 2069 | |
---|
| 2070 | // First check LHEF header (if reading from LHEF) |
---|
| 2071 | if (readFrom == 1) { |
---|
| 2072 | if (lhefHeader != "void") |
---|
| 2073 | ifailLHE = slha.readFile(lhefHeader, verboseSLHA, slhaUseDec ); |
---|
| 2074 | else if (lhefFile != "void") |
---|
| 2075 | ifailLHE = slha.readFile(lhefFile, verboseSLHA, slhaUseDec ); |
---|
| 2076 | } |
---|
| 2077 | |
---|
| 2078 | // If LHEF read successful, everything needed should already be ready |
---|
| 2079 | if (ifailLHE == 0) { |
---|
| 2080 | ifailSpc = 0; |
---|
| 2081 | couplingsPtr->isSUSY = true; |
---|
| 2082 | // If no LHEF file or no SLHA info in header, read from SLHA:file |
---|
| 2083 | } else { |
---|
| 2084 | lhefFile = "void"; |
---|
| 2085 | if ( settings.word("SLHA:file") == "none" |
---|
| 2086 | || settings.word("SLHA:file") == "void" |
---|
| 2087 | || settings.word("SLHA:file") == "" |
---|
| 2088 | || settings.word("SLHA:file") == " ") return true; |
---|
| 2089 | ifailSpc = slha.readFile(slhaFile,verboseSLHA, slhaUseDec); |
---|
| 2090 | } |
---|
| 2091 | |
---|
| 2092 | // In case of problems, print error and fail init. |
---|
| 2093 | if (ifailSpc != 0) { |
---|
| 2094 | info.errorMsg("Error in Pythia::initSLHA: " |
---|
| 2095 | "problem reading SLHA file", slhaFile); |
---|
| 2096 | return false; |
---|
| 2097 | } else { |
---|
| 2098 | couplingsPtr->isSUSY = true; |
---|
| 2099 | } |
---|
| 2100 | |
---|
| 2101 | // Check spectrum for consistency. Switch off SUSY if necessary. |
---|
| 2102 | ifailSpc = slha.checkSpectrum(); |
---|
| 2103 | |
---|
| 2104 | // ifail >= 1 : no MODSEL found -> don't switch on SUSY |
---|
| 2105 | if (ifailSpc == 1) { |
---|
| 2106 | // no SUSY, but MASS ok |
---|
| 2107 | couplingsPtr->isSUSY = false; |
---|
| 2108 | info.errorMsg("Warning in Pythia::initSLHA: " |
---|
| 2109 | "No MODSEL found, keeping internal SUSY switched off"); |
---|
| 2110 | } else if (ifailSpc >= 2) { |
---|
| 2111 | // no SUSY, but problems |
---|
| 2112 | info.errorMsg("Warning in Pythia::initSLHA: " |
---|
| 2113 | "Problem with SLHA MASS or QNUMBERS."); |
---|
| 2114 | couplingsPtr->isSUSY = false; |
---|
| 2115 | } |
---|
| 2116 | // ifail = 0 : MODSEL found, spectrum OK |
---|
| 2117 | else if (ifailSpc == 0) { |
---|
| 2118 | // Print spectrum. Done. |
---|
| 2119 | slha.printSpectrum(0); |
---|
| 2120 | } |
---|
| 2121 | else if (ifailSpc < 0) { |
---|
| 2122 | info.errorMsg("Warning in Pythia::initSLHA: " |
---|
| 2123 | "Problem with SLHA spectrum.", |
---|
| 2124 | "\n Only using masses and switching off SUSY."); |
---|
| 2125 | settings.flag("SUSY:all", false); |
---|
| 2126 | couplingsPtr->isSUSY = false; |
---|
| 2127 | slha.printSpectrum(ifailSpc); |
---|
| 2128 | } |
---|
| 2129 | |
---|
| 2130 | // Import qnumbers |
---|
| 2131 | if ( (ifailSpc == 1 || ifailSpc == 0) && slha.qnumbers.size() > 0) { |
---|
| 2132 | for (int iQnum=0; iQnum < int(slha.qnumbers.size()); iQnum++) { |
---|
| 2133 | // Always use positive id codes |
---|
| 2134 | int id = abs(slha.qnumbers[iQnum](0)); |
---|
| 2135 | ostringstream idCode; |
---|
| 2136 | idCode << id; |
---|
| 2137 | if (particleData.isParticle(id)) { |
---|
| 2138 | info.errorMsg("Warning in Pythia::initSLHA: " |
---|
| 2139 | "ignoring QNUMBERS", "for id = "+idCode.str() |
---|
| 2140 | +" (already exists)", true); |
---|
| 2141 | } else { |
---|
| 2142 | int qEM3 = slha.qnumbers[iQnum](1); |
---|
| 2143 | int nSpins = slha.qnumbers[iQnum](2); |
---|
| 2144 | int colRep = slha.qnumbers[iQnum](3); |
---|
| 2145 | int hasAnti = slha.qnumbers[iQnum](4); |
---|
| 2146 | // Translate colRep to PYTHIA colType |
---|
| 2147 | int colType = 0; |
---|
| 2148 | if (colRep == 3) colType = 1; |
---|
| 2149 | else if (colRep == -3) colType = -1; |
---|
| 2150 | else if (colRep == 8) colType = 2; |
---|
| 2151 | else if (colRep == 6) colType = 3; |
---|
| 2152 | else if (colRep == -6) colType = -3; |
---|
| 2153 | // Default name: PDG code |
---|
| 2154 | string name, antiName; |
---|
| 2155 | ostringstream idStream; |
---|
| 2156 | idStream<<id; |
---|
| 2157 | name = idStream.str(); |
---|
| 2158 | antiName = "-"+name; |
---|
| 2159 | if (iQnum < int(slha.qnumbersName.size())) { |
---|
| 2160 | name = slha.qnumbersName[iQnum]; |
---|
| 2161 | if (iQnum < int(slha.qnumbersAntiName.size())) |
---|
| 2162 | antiName = slha.qnumbersAntiName[iQnum]; |
---|
| 2163 | if (antiName == "") { |
---|
| 2164 | if (name.find("+") != string::npos) { |
---|
| 2165 | antiName = name; |
---|
| 2166 | antiName.replace(antiName.find("+"),1,"-"); |
---|
| 2167 | } else if (name.find("-") != string::npos) { |
---|
| 2168 | antiName = name; |
---|
| 2169 | antiName.replace(antiName.find("-"),1,"+"); |
---|
| 2170 | } else { |
---|
| 2171 | antiName = name+"bar"; |
---|
| 2172 | } |
---|
| 2173 | } |
---|
| 2174 | } |
---|
| 2175 | if ( hasAnti == 0) { |
---|
| 2176 | antiName = ""; |
---|
| 2177 | particleData.addParticle(id, name, nSpins, qEM3, colType); |
---|
| 2178 | } else { |
---|
| 2179 | particleData.addParticle(id, name, antiName, nSpins, qEM3, colType); |
---|
| 2180 | } |
---|
| 2181 | } |
---|
| 2182 | } |
---|
| 2183 | } |
---|
| 2184 | |
---|
| 2185 | // Import mass spectrum. |
---|
| 2186 | bool keepSM = settings.flag("SLHA:keepSM"); |
---|
| 2187 | double minMassSM = settings.parm("SLHA:minMassSM"); |
---|
| 2188 | double massMargin = settings.parm("SLHA:minDecayDeltaM"); |
---|
| 2189 | if (ifailSpc == 1 || ifailSpc == 0) { |
---|
| 2190 | |
---|
| 2191 | // Loop through to update particle data. |
---|
| 2192 | int id = slha.mass.first(); |
---|
| 2193 | for (int i = 1; i <= slha.mass.size() ; i++) { |
---|
| 2194 | double mass = abs(slha.mass(id)); |
---|
| 2195 | |
---|
| 2196 | // Ignore masses for known SM particles or particles with |
---|
| 2197 | // default masses < minMassSM; overwrite masses for rest. |
---|
| 2198 | if (keepSM && (id < 25 || (id > 80 && id < 1000000))) ; |
---|
| 2199 | else if (id < 1000000 && particleData.m0(id) < minMassSM) { |
---|
| 2200 | ostringstream idCode; |
---|
| 2201 | idCode << id; |
---|
| 2202 | info.errorMsg("Warning in Pythia::initSLHA: " |
---|
| 2203 | "ignoring MASS entry", "for id = "+idCode.str() |
---|
| 2204 | +" (m0 < SLHA:minMassSM)", true); |
---|
| 2205 | } |
---|
| 2206 | else particleData.m0(id,mass); |
---|
| 2207 | id = slha.mass.next(); |
---|
| 2208 | }; |
---|
| 2209 | |
---|
| 2210 | } |
---|
| 2211 | |
---|
| 2212 | // Update decay data. |
---|
| 2213 | for (int iTable=0; iTable < int(slha.decays.size()); iTable++) { |
---|
| 2214 | |
---|
| 2215 | // Pointer to this SLHA table |
---|
| 2216 | LHdecayTable* slhaTable=&(slha.decays[iTable]); |
---|
| 2217 | |
---|
| 2218 | // Extract ID and create pointer to corresponding particle data object |
---|
| 2219 | int idRes = slhaTable->getId(); |
---|
| 2220 | ParticleDataEntry* particlePtr |
---|
| 2221 | = particleData.particleDataEntryPtr(idRes); |
---|
| 2222 | |
---|
| 2223 | // Ignore decay channels for known SM particles or particles with |
---|
| 2224 | // default masses < minMassSM; overwrite masses for rest. |
---|
| 2225 | if (keepSM && (idRes < 25 || (idRes > 80 && idRes < 1000000))) continue; |
---|
| 2226 | else if (idRes < 1000000 && particleData.m0(idRes) < minMassSM) { |
---|
| 2227 | ostringstream idCode; |
---|
| 2228 | idCode << idRes; |
---|
| 2229 | info.errorMsg("Warning in Pythia::initSLHA: " |
---|
| 2230 | "ignoring DECAY table", "for id = " + idCode.str() |
---|
| 2231 | + " (m0 < SLHA:minMassSM)", true); |
---|
| 2232 | continue; |
---|
| 2233 | } |
---|
| 2234 | |
---|
| 2235 | // Extract and store total width (absolute value, neg -> switch off) |
---|
| 2236 | double widRes = abs(slhaTable->getWidth()); |
---|
| 2237 | particlePtr->setMWidth(widRes); |
---|
| 2238 | |
---|
| 2239 | // Reset decay table of the particle. Allow decays, treat as resonance. |
---|
| 2240 | if (slhaTable->size() > 0) { |
---|
| 2241 | particlePtr->clearChannels(); |
---|
| 2242 | particleData.mayDecay(idRes,true); |
---|
| 2243 | particleData.isResonance(idRes,true); |
---|
| 2244 | } |
---|
| 2245 | |
---|
| 2246 | // Reset to stable if width <= 0.0 |
---|
| 2247 | if (slhaTable->getWidth() <= 0.0) particleData.mayDecay(idRes,false); |
---|
| 2248 | |
---|
| 2249 | // Set initial minimum mass. |
---|
| 2250 | double brWTsum = 0.; |
---|
| 2251 | double massWTsum = 0.; |
---|
| 2252 | |
---|
| 2253 | // Loop over SLHA channels, import into Pythia, treating channels |
---|
| 2254 | // with negative branching fractions as having the equivalent positive |
---|
| 2255 | // branching fraction, but being switched off for this run |
---|
| 2256 | for (int iChannel=0 ; iChannel<slhaTable->size(); iChannel++) { |
---|
| 2257 | LHdecayChannel slhaChannel = slhaTable->getChannel(iChannel); |
---|
| 2258 | double brat = slhaChannel.getBrat(); |
---|
| 2259 | vector<int> idDa = slhaChannel.getIdDa(); |
---|
| 2260 | if (idDa.size() >= 9) { |
---|
| 2261 | info.errorMsg("Error in Pythia::initSLHA: " |
---|
| 2262 | "max number of decay products is 8."); |
---|
| 2263 | } else if (idDa.size() <= 1) { |
---|
| 2264 | info.errorMsg("Error in Pythia::initSLHA: " |
---|
| 2265 | "min number of decay products is 2."); |
---|
| 2266 | } |
---|
| 2267 | else { |
---|
| 2268 | int onMode = 1; |
---|
| 2269 | if (brat < 0.0) onMode = 0; |
---|
| 2270 | |
---|
| 2271 | // Check phase space, including margin |
---|
| 2272 | double massSum = massMargin; |
---|
| 2273 | for (int jDa=0; jDa<int(idDa.size()); ++jDa) |
---|
| 2274 | massSum += particleData.m0( idDa[jDa] ); |
---|
| 2275 | if (onMode == 1 && brat > 0.0 && massSum > particleData.m0(idRes) ) { |
---|
| 2276 | // String containing decay name |
---|
| 2277 | ostringstream errCode; |
---|
| 2278 | errCode << idRes <<" ->"; |
---|
| 2279 | for (int jDa=0; jDa<int(idDa.size()); ++jDa) errCode<<" "<<idDa[jDa]; |
---|
| 2280 | info.errorMsg("Warning in Pythia::initSLHA: switching off decay", |
---|
| 2281 | errCode.str() + " (mRes - mDa < minDecayDeltaM)" |
---|
| 2282 | "\n (Note: cross sections will be scaled by remaining" |
---|
| 2283 | " open branching fractions!)" , true); |
---|
| 2284 | onMode=0; |
---|
| 2285 | } |
---|
| 2286 | |
---|
| 2287 | // Branching-ratio-weighted average mass in decay. |
---|
| 2288 | brWTsum += abs(brat); |
---|
| 2289 | massWTsum += abs(brat) * massSum; |
---|
| 2290 | |
---|
| 2291 | // Add channel |
---|
| 2292 | int id0 = idDa[0]; |
---|
| 2293 | int id1 = idDa[1]; |
---|
| 2294 | int id2 = (idDa.size() >= 3) ? idDa[2] : 0; |
---|
| 2295 | int id3 = (idDa.size() >= 4) ? idDa[3] : 0; |
---|
| 2296 | int id4 = (idDa.size() >= 5) ? idDa[4] : 0; |
---|
| 2297 | int id5 = (idDa.size() >= 6) ? idDa[5] : 0; |
---|
| 2298 | int id6 = (idDa.size() >= 7) ? idDa[6] : 0; |
---|
| 2299 | int id7 = (idDa.size() >= 8) ? idDa[7] : 0; |
---|
| 2300 | particlePtr->addChannel(onMode,abs(brat),101, |
---|
| 2301 | id0,id1,id2,id3,id4,id5,id6,id7); |
---|
| 2302 | |
---|
| 2303 | } |
---|
| 2304 | } |
---|
| 2305 | |
---|
| 2306 | // Set minimal mass, but always below nominal one. |
---|
| 2307 | if (slhaTable->size() > 0) { |
---|
| 2308 | double massAvg = massWTsum / brWTsum; |
---|
| 2309 | double massMin = min( massAvg, particlePtr->m0()) - massMargin; |
---|
| 2310 | particlePtr->setMMin(massMin); |
---|
| 2311 | } |
---|
| 2312 | } |
---|
| 2313 | |
---|
| 2314 | return true; |
---|
| 2315 | |
---|
| 2316 | } |
---|
| 2317 | |
---|
| 2318 | //-------------------------------------------------------------------------- |
---|
| 2319 | |
---|
| 2320 | // Function to perform CKKW-L merging on the input event. |
---|
| 2321 | |
---|
| 2322 | bool Pythia::mergeProcess() { |
---|
| 2323 | |
---|
| 2324 | // Reset weight of the event. |
---|
| 2325 | mergingHooksPtr->setWeightCKKWL(1.); |
---|
| 2326 | info.setWeightCKKWL(1.); |
---|
| 2327 | double wgt = 1.0; |
---|
| 2328 | |
---|
| 2329 | // Store candidates for the splitting V -> qqbar'. |
---|
| 2330 | mergingHooksPtr->storeHardProcessCandidates( process); |
---|
| 2331 | |
---|
| 2332 | // Calculate number of clustering steps. |
---|
| 2333 | int nSteps = mergingHooksPtr->getNumberOfClusteringSteps( process); |
---|
| 2334 | |
---|
| 2335 | // Check if event passes the merging scale cut. |
---|
| 2336 | double tms = mergingHooksPtr->tms(); |
---|
| 2337 | // Get merging scale in current event. |
---|
| 2338 | double tnow = 0.; |
---|
| 2339 | // Use KT/Durham merging scale definition. |
---|
| 2340 | if ( mergingHooksPtr->doKTMerging() || mergingHooksPtr->doMGMerging() ) |
---|
| 2341 | tnow = mergingHooksPtr->kTms( process ); |
---|
| 2342 | // Use Lund PT merging scale definition. |
---|
| 2343 | else if ( mergingHooksPtr->doPTLundMerging() ) |
---|
| 2344 | tnow = mergingHooksPtr->rhoms( process, false ); |
---|
| 2345 | // Use DeltaR_{ij}, pT_i, Q_{ij} combination merging scale definition. |
---|
| 2346 | else if ( mergingHooksPtr->doCutBasedMerging() ) |
---|
| 2347 | tnow = mergingHooksPtr->cutbasedms( process ); |
---|
| 2348 | // Use user-defined merging scale. |
---|
| 2349 | else |
---|
| 2350 | tnow = mergingHooksPtr->tmsDefinition( process ); |
---|
| 2351 | |
---|
| 2352 | bool enforceCutOnLHE = settings.flag("Merging:enforceCutOnLHE"); |
---|
| 2353 | // Enfore merging scale cut if the event did not pass the merging scale |
---|
| 2354 | // criterion. |
---|
| 2355 | if ( enforceCutOnLHE && nSteps > 0 && tnow < tms ) { |
---|
| 2356 | string message="Warning in Pythia::mergeProcess: Les Houches Event"; |
---|
| 2357 | message+=" fails merging scale cut. Cut by rejecting event."; |
---|
| 2358 | info.errorMsg(message); |
---|
| 2359 | return true; |
---|
| 2360 | } |
---|
| 2361 | |
---|
| 2362 | double TMSMISMATCHPOS = tms / 2. ; |
---|
| 2363 | // Remember if event is significantly above the merging scale. |
---|
| 2364 | if ( enforceCutOnLHE && nSteps > 0 && tms > 0. |
---|
| 2365 | && tnow > tms && (tnow-tms > TMSMISMATCHPOS) ) { |
---|
| 2366 | stringstream tmsmis; |
---|
| 2367 | tmsmis << int(TMSMISMATCHPOS); |
---|
| 2368 | string message="Warning in Pythia::mergeProcess: Les Houches event"; |
---|
| 2369 | message+=" more than "; |
---|
| 2370 | message+=tmsmis.str(); |
---|
| 2371 | message+=" GeV above desired merging scale value."; |
---|
| 2372 | info.errorMsg(message); |
---|
| 2373 | info.addCounter(41); |
---|
| 2374 | } |
---|
| 2375 | |
---|
| 2376 | // For now, prefer construction of ordered histories. |
---|
| 2377 | mergingHooksPtr->orderHistories(true); |
---|
| 2378 | |
---|
| 2379 | // Declare pdfWeight and alpha_s weight. |
---|
| 2380 | double RN = rndm.flat(); |
---|
| 2381 | |
---|
| 2382 | process.scale(0.0); |
---|
| 2383 | // Generate all histories. |
---|
| 2384 | History FullHistory( nSteps, 0.0, process, Clustering(), mergingHooksPtr, |
---|
| 2385 | beamA, beamB, &particleData, &info, true, true, true, true, |
---|
| 2386 | 1.0, 0); |
---|
| 2387 | |
---|
| 2388 | // Project histories onto desired branches, e.g. only ordered paths. |
---|
| 2389 | FullHistory.projectOntoDesiredHistories(); |
---|
| 2390 | |
---|
| 2391 | // Calculate CKKWL weight: |
---|
| 2392 | // Perform reweighting with Sudakov factors, save alpha_s ratios and |
---|
| 2393 | // PDF ratio weights. |
---|
| 2394 | wgt = FullHistory.weightTREE( &trialPartonLevel, |
---|
| 2395 | mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaS_ISR(), RN); |
---|
| 2396 | |
---|
| 2397 | // Event with production scales set for further (trial) showering |
---|
| 2398 | // and starting conditions for the shower. |
---|
| 2399 | FullHistory.getStartingConditions( RN, process ); |
---|
| 2400 | |
---|
| 2401 | // Allow to dampen histories in which the lowest multiplicity reclustered |
---|
| 2402 | // state does not pass the lowest multiplicity cut of the matrix element. |
---|
| 2403 | double dampWeight = mergingHooksPtr->dampenIfFailCuts( |
---|
| 2404 | FullHistory.lowestMultProc(RN) ); |
---|
| 2405 | |
---|
| 2406 | // Save the weight of the event for histogramming. Only change the |
---|
| 2407 | // event weight after trial shower on the matrix element |
---|
| 2408 | // multiplicity event (= in doVetoStep). |
---|
| 2409 | wgt *= dampWeight; |
---|
| 2410 | |
---|
| 2411 | // Set QCD 2->2 starting scale different from arbitrary scale in LHEF! |
---|
| 2412 | // --> Set to minimal mT of partons. |
---|
| 2413 | int nFinal = 0; |
---|
| 2414 | double muf = process[0].e(); |
---|
| 2415 | for ( int i=0; i < process.size(); ++i ) |
---|
| 2416 | if ( process[i].isFinal() |
---|
| 2417 | && (process[i].colType() != 0 || process[i].id() == 22 ) ) { |
---|
| 2418 | nFinal++; |
---|
| 2419 | muf = min( muf, abs(process[i].mT()) ); |
---|
| 2420 | } |
---|
| 2421 | // For pure QCD dijet events (only!), set the process scale to the |
---|
| 2422 | // transverse momentum of the outgoing partons. |
---|
| 2423 | if ( nSteps == 0 && nFinal == 2 |
---|
| 2424 | && ( mergingHooksPtr->getProcessString().compare("pp>jj") == 0 |
---|
| 2425 | || mergingHooksPtr->getProcessString().compare("pp>aj") == 0) ) |
---|
| 2426 | process.scale(muf); |
---|
| 2427 | |
---|
| 2428 | // Save the weight of the event for histogramming. |
---|
| 2429 | mergingHooksPtr->setWeightCKKWL(wgt); |
---|
| 2430 | info.setWeightCKKWL(wgt); |
---|
| 2431 | |
---|
| 2432 | // If the weight of the event is zero, do not continue evolution. |
---|
| 2433 | if(wgt == 0.) return true; |
---|
| 2434 | |
---|
| 2435 | // Done |
---|
| 2436 | return false; |
---|
| 2437 | |
---|
| 2438 | } |
---|
| 2439 | |
---|
| 2440 | //========================================================================== |
---|
| 2441 | |
---|
| 2442 | } // end namespace Pythia8 |
---|