PYTHIA8/pythia8145/include/PhaseSpace.h

   1 // PhaseSpace.h is a part of the PYTHIA event generator.
   2 // Copyright (C) 2010 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 // Header file for phase space generators in kinematics selection.
   7 // PhaseSpace: base class for phase space generators.
   8 // Base class for derived classes> 2 ->1 , 2 -> 2, 2 -> 2 elastic/diffractive,
   9 // 2 -> 2 minbias, 2 -> 3, Les Houches.
  10
  11 #ifndef Pythia8_PhaseSpace_H
  12 #define Pythia8_PhaseSpace_H
  13
  14 #include "Basics.h"
  15 #include "BeamParticle.h"
  16 #include "Info.h"
  17 #include "LesHouches.h"
  18 #include "MultipleInteractions.h"
  19 #include "ParticleData.h"
  20 #include "PartonDistributions.h"
  21 #include "PythiaStdlib.h"
  22 #include "SigmaProcess.h"
  23 #include "SigmaTotal.h"
  24 #include "Settings.h"
  25 #include "StandardModel.h"
  26 #include "UserHooks.h"
  27
  28 namespace Pythia8 {
  29
  30 //==========================================================================
  31
  32 // Forward reference to the UserHooks class.
  33 class UserHooks;
  34
  35 //==========================================================================
  36
  37 // PhaseSpace is a base class for  phase space generators
  38 // used in the selection of hard-process kinematics.
  39
  40 class PhaseSpace {
  41
  42 public:
  43
  44   // Destructor.
  45   virtual ~PhaseSpace() {}
  46
  47   // Perform simple initialization and store pointers.
  48   void init(bool isFirst, SigmaProcess* sigmaProcessPtrIn,
  49     Info* infoPtrIn, Settings* settingsPtrIn, ParticleData* particleDataPtrIn,
  50     Rndm* rndmPtrIn, BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn,
  51     Couplings* couplingsPtrIn, SigmaTotal* sigmaTotPtrIn, UserHooks* userHooksPtrIn);
  52
  53   // Update the CM energy of the event.
  54   void newECM(double eCMin) {eCM = eCMin; s = eCM * eCM;}
  55
  56   // Store or replace Les Houches pointer.
  57   void setLHAPtr(LHAup* lhaUpPtrIn) {lhaUpPtr = lhaUpPtrIn;}
  58
  59   // A pure virtual method, wherein an optimization procedure
  60   // is used to determine how phase space should be sampled.
  61   virtual bool setupSampling() = 0;
  62
  63   // A pure virtual method, wherein a trial event kinematics
  64   // is to be selected in the derived class.
  65   virtual bool trialKin(bool inEvent = true, bool repeatSame = false) = 0;
  66
  67   // A pure virtual method, wherein the accepted event kinematics
  68   // is to be constructed in the derived class.
  69   virtual bool finalKin() = 0;
  70
  71   // Allow for nonisotropic decays when ME's available.
  72   void   decayKinematics( Event& process);
  73
  74   // Give back current or maximum cross section, or set latter.
  75   double sigmaNow() const {return sigmaNw;}
  76   double sigmaMax() const {return sigmaMx;}
  77   bool   newSigmaMax() const {return newSigmaMx;}
  78   void   setSigmaMax(double sigmaMaxIn) {sigmaMx = sigmaMaxIn;}
  79
  80   // For Les Houches with negative event weight needs
  81   virtual double sigmaSumSigned() const {return sigmaMx;}
  82
  83   // Give back constructed four-vectors and known masses.
  84   Vec4   p(int i)   const {return pH[i];}
  85   double m(int i)   const {return mH[i];}
  86
  87   // Give back other event properties.
  88   double ecm()      const {return eCM;}
  89   double x1()       const {return x1H;}
  90   double x2()       const {return x2H;}
  91   double sHat()     const {return sH;}
  92   double tHat()     const {return tH;}
  93   double uHat()     const {return uH;}
  94   double pTHat()    const {return pTH;}
  95   double thetaHat() const {return theta;}
  96   double phiHat()   const {return phi;}
  97   double runBW3()   const {return runBW3H;}
  98   double runBW4()   const {return runBW4H;}
  99   double runBW5()   const {return runBW5H;}
 100
 101   // Inform whether beam particles are resolved in partons or scatter directly.
 102   virtual bool isResolved() const {return true;}
 103
 104 protected:
 105
 106   // Constructor.
 107   PhaseSpace() {}
 108
 109   // Constants: could only be changed in the code itself.
 110   static const int    NMAXTRY, NTRY3BODY;
 111   static const double SAFETYMARGIN, TINY, EVENFRAC, SAMESIGMA, WIDTHMARGIN,
 112                       SAMEMASS, MASSMARGIN, EXTRABWWTMAX, THRESHOLDSIZE,
 113                       THRESHOLDSTEP, YRANGEMARGIN, LEPTONXMIN, LEPTONXMAX,
 114                       LEPTONXLOGMIN, LEPTONXLOGMAX, LEPTONTAUMIN,
 115                       SHATMINZ, PT2RATMINZ, WTCORRECTION[11];
 116
 117   // Pointer to cross section.
 118   SigmaProcess* sigmaProcessPtr;
 119
 120   // Pointer to various information on the generation.
 121   Info*         infoPtr;
 122
 123   // Pointer to the settings database.
 124   Settings*     settingsPtr;
 125
 126   // Pointer to the particle data table.
 127   ParticleData* particleDataPtr;
 128
 129   // Pointer to the random number generator.
 130   Rndm*         rndmPtr;
 131
 132   // Pointers to incoming beams.
 133   BeamParticle* beamAPtr;
 134   BeamParticle* beamBPtr;
 135
 136   // Pointer to Standard Model couplings.
 137   Couplings*         couplingsPtr;
 138
 139   // Pointer to the total/elastic/diffractive cross section object.
 140   SigmaTotal*   sigmaTotPtr;
 141
 142   // Pointer to userHooks object for user interaction with program.
 143   UserHooks*    userHooksPtr;
 144
 145   // Pointer to LHAup for generating external events.
 146   LHAup*        lhaUpPtr;
 147
 148   // Initialization data, normally only set once.
 149   bool   useBreitWigners, doEnergySpread, showSearch, showViolation,
 150          increaseMaximum;
 151   int    gmZmodeGlobal;
 152   double mHatGlobalMin, mHatGlobalMax, pTHatGlobalMin, pTHatGlobalMax,
 153          pTHatMinDiverge, minWidthBreitWigners;
 154
 155   // Information on incoming beams.
 156   int    idA, idB;
 157   double mA, mB, eCM, s;
 158   bool   hasLeptonBeams, hasPointLeptons;
 159
 160  // Cross section information.
 161   bool   newSigmaMx, canModifySigma;
 162   int    gmZmode;
 163   double wtBW, sigmaNw, sigmaMx, sigmaPos, sigmaNeg;
 164
 165   // Process-specific kinematics properties, almost always available.
 166   double mHatMin, mHatMax, sHatMin, sHatMax, pTHatMin, pTHatMax,
 167          pT2HatMin, pT2HatMax;
 168
 169   // Event-specific kinematics properties, almost always available.
 170   double x1H, x2H, m3, m4, m5, s3, s4, s5, mHat, sH, tH, uH, pAbs, p2Abs,
 171          pTH, theta, phi, betaZ;
 172   Vec4   pH[6];
 173   double mH[6];
 174
 175   // Reselect decay products momenta isotropically in phase space.
 176   void decayKinematicsStep( Event& process, int iRes);
 177
 178   // Much common code for normal 2 -> 1, 2 -> 2 and 2 -> 3 cases:
 179
 180   // Determine how phase space should be sampled.
 181   void setup3Body();
 182   bool setupSampling123(bool is2, bool is3, ostream& os = cout);
 183
 184   // Select a trial kinematics phase space point.
 185   bool trialKin123(bool is2, bool is3, bool inEvent = true, ostream& os = cout);
 186
 187   // Presence and properties of any s-channel resonances.
 188   int    idResA, idResB;
 189   double mResA, mResB, GammaResA, GammaResB, tauResA, tauResB, widResA,
 190          widResB;
 191   bool   sameResMass;
 192
 193   // Kinematics properties specific to 2 -> 1/2/3.
 194   bool   useMirrorWeight;
 195   double tau, y, z, tauMin, tauMax, yMax, zMin, zMax, ratio34, unity34,
 196          zNeg, zPos, wtTau, wtY, wtZ, wt3Body, runBW3H, runBW4H, runBW5H,
 197          intTau0, intTau1, intTau2, intTau3, intTau4, intTau5, intTau6,
 198          intY0, intY12, intY34, intY56, mTchan1, sTchan1, mTchan2, sTchan2,
 199          frac3Flat, frac3Pow1, frac3Pow2;
 200   Vec4   p3cm, p4cm, p5cm;
 201
 202   // Coefficients for optimized selection in 2 -> 1/2/3.
 203   int    nTau, nY, nZ;
 204   double tauCoef[8], yCoef[8], zCoef[8], tauCoefSum[8], yCoefSum[8],
 205          zCoefSum[8];
 206
 207   // Calculate kinematical limits for 2 -> 1/2/3.
 208   bool limitTau(bool is2, bool is3);
 209   bool limitY();
 210   bool limitZ();
 211
 212   // Select kinematical variable between defined limits for 2 -> 1/2/3.
 213   void selectTau(int iTau, double tauVal, bool is2);
 214   void selectY(int iY, double yVal);
 215   void selectZ(int iZ, double zVal);
 216   bool select3Body();
 217
 218   // Solve equation system for better phase space coefficients in 2 -> 1/2/3.
 219   void solveSys( int n, int bin[8], double vec[8], double mat[8][8],
 220     double coef[8], ostream& os = cout);
 221
 222   // Properties specific to resonance mass selection in 2 -> 2 and 2 -> 3.
 223   bool   useBW[6];
 224   int    idMass[6];
 225   double mPeak[6], sPeak[6], mWidth[6], mMin[6], mMax[6], mw[6], wmRat[6],
 226          mLower[6], mUpper[6], sLower[6], sUpper[6], fracFlat[6], fracInv[6],
 227          fracInv2[6], atanLower[6], atanUpper[6], intBW[6], intFlat[6],
 228          intInv[6], intInv2[6];
 229
 230   // Setup mass selection for one resonance at a time. Split in two parts.
 231   void   setupMass1(int iM);
 232   void   setupMass2(int iM, double distToThresh);
 233
 234   // Do mass selection and find the associated weight.
 235   void   trialMass(int iM);
 236   double weightMass(int iM);
 237
 238 };
 239
 240 //==========================================================================
 241
 242 // A derived class with 2 -> 1 kinematics set up in tau, y.
 243
 244 class PhaseSpace2to1tauy : public PhaseSpace {
 245
 246 public:
 247
 248   // Constructor.
 249   PhaseSpace2to1tauy() {}
 250
 251   // Optimize subsequent kinematics selection.
 252   virtual bool setupSampling() {if (!setupMass()) return false;
 253     return setupSampling123(false, false);}
 254
 255   // Construct the trial kinematics.
 256   virtual bool trialKin(bool inEvent = true, bool = false) {wtBW = 1.;
 257     return trialKin123(false, false, inEvent);}
 258
 259   // Construct the final event kinematics.
 260   virtual bool finalKin();
 261
 262 private:
 263
 264   // Set up allowed mass range.
 265   bool setupMass();
 266
 267 };
 268
 269 //==========================================================================
 270
 271 // A derived class with 2 -> 2 kinematics set up in tau, y, z = cos(theta).
 272
 273 class PhaseSpace2to2tauyz : public PhaseSpace {
 274
 275 public:
 276
 277   // Constructor.
 278   PhaseSpace2to2tauyz() {}
 279
 280   // Optimize subsequent kinematics selection.
 281   virtual bool setupSampling() {if (!setupMasses()) return false;
 282     return setupSampling123(true, false);}
 283
 284   // Construct the trial kinematics.
 285   virtual bool trialKin(bool inEvent = true, bool = false) {
 286     if (!trialMasses()) return false;
 287     return trialKin123(true, false, inEvent);}
 288
 289   // Construct the final event kinematics.
 290   virtual bool finalKin();
 291
 292 private:
 293
 294   // Set up for fixed or Breit-Wigner mass selection.
 295   bool setupMasses();
 296
 297   // Select fixed or Breit-Wigner-distributed masses.
 298   bool trialMasses();
 299
 300   // Pick off-shell initialization masses when on-shell not allowed.
 301   bool constrainedM3M4();
 302   bool constrainedM3();
 303   bool constrainedM4();
 304
 305 };
 306
 307 //==========================================================================
 308
 309 // A derived class with 2 -> 2 kinematics set up for elastic scattering.
 310
 311 class PhaseSpace2to2elastic : public PhaseSpace {
 312
 313 public:
 314
 315   // Constructor.
 316   PhaseSpace2to2elastic() {}
 317
 318   // Construct the trial or final event kinematics.
 319   virtual bool setupSampling();
 320   virtual bool trialKin(bool inEvent = true, bool = false);
 321   virtual bool finalKin();
 322
 323   // Are beam particles resolved in partons or scatter directly?
 324   virtual bool isResolved() const {return false;}
 325
 326 private:
 327
 328   // Constants: could only be changed in the code itself.
 329   static const double EXPMAX, CONVERTEL;
 330
 331   // Kinematics properties specific to 2 -> 2 elastic.
 332   bool   useCoulomb;
 333   double s1, s2, bSlope, lambda12S, tLow, tUpp, tAux, sigmaTot, rho,
 334          lambda, tAbsMin, phaseCst, alphaEM0, sigmaNuc, sigmaCou, signCou;
 335
 336  // Calculation of alphaElectromagnetic.
 337  AlphaEM alphaEM;
 338
 339 };
 340
 341 //==========================================================================
 342
 343 // A derived class with 2 -> 2 kinematics set up for diffractive scattering.
 344
 345 class PhaseSpace2to2diffractive : public PhaseSpace {
 346
 347 public:
 348
 349   // Constructor.
 350   PhaseSpace2to2diffractive(bool isDiffAin = false, bool isDiffBin = false)
 351     : isDiffA(isDiffAin), isDiffB(isDiffBin) {}
 352
 353   // Construct the trial or final event kinematics.
 354   virtual bool setupSampling();
 355   virtual bool trialKin(bool inEvent = true, bool = false);
 356   virtual bool finalKin();
 357
 358   // Are beam particles resolved in partons or scatter directly?
 359   virtual bool isResolved() const {return false;}
 360
 361 private:
 362
 363   // Constants: could only be changed in the code itself.
 364   static const int    NTRY;
 365   static const double EXPMAX, DIFFMASSMAX;
 366
 367   // Initialization data, in constructor or read from Settings.
 368   bool   isDiffA, isDiffB;
 369   int    PomFlux;
 370   double epsilonPF, alphaPrimePF;
 371
 372   // Initialization: kinematics properties specific to 2 -> 2 diffractive.
 373   double m3ElDiff, m4ElDiff, s1, s2, lambda12, lambda34, tLow, tUpp,
 374          cRes, sResXB, sResAX, sProton, bMin, bSlope, bSlope1, bSlope2,
 375          probSlope1, xIntPF, xtCorPF, mp24DL, coefDL, tAux, tAux1, tAux2;
 376
 377 };
 378
 379 //==========================================================================
 380
 381 // A derived class for minumum bias events. Hardly does anything, since
 382 // the real action is taken care of by the MultipleInteractions class.
 383
 384 class PhaseSpace2to2minbias : public PhaseSpace {
 385
 386 public:
 387
 388   // Constructor.
 389   PhaseSpace2to2minbias() {}
 390
 391   // Construct the trial or final event kinematics.
 392   virtual bool setupSampling() {sigmaNw = sigmaProcessPtr->sigmaHat();
 393     sigmaMx = sigmaNw; return true;}
 394   virtual bool trialKin( bool , bool = false) {return true;}
 395   virtual bool finalKin() {return true;}
 396
 397 private:
 398
 399 };
 400
 401 //==========================================================================
 402
 403 // A derived class with 2 -> 3 kinematics 1 + 2 -> 3 + 4 + 5 set up in
 404 // tau, y, pT2_4, pT2_5, phi_4, phi_5 and y_3 (partial cylindrical symmetry).
 405
 406 class PhaseSpace2to3tauycyl : public PhaseSpace {
 407
 408 public:
 409
 410   // Constructor.
 411   PhaseSpace2to3tauycyl() {}
 412
 413   // Optimize subsequent kinematics selection.
 414   virtual bool setupSampling() {if (!setupMasses()) return false;
 415     setup3Body(); return setupSampling123(false, true);}
 416
 417   // Construct the trial kinematics.
 418   virtual bool trialKin(bool inEvent = true, bool = false) {
 419     if (!trialMasses()) return false;
 420     return trialKin123(false, true, inEvent);}
 421
 422   // Construct the final event kinematics.
 423   virtual bool finalKin();
 424
 425 private:
 426
 427   // Constants: could only be changed in the code itself.
 428   static const int    NITERNR;
 429
 430   // Set up for fixed or Breit-Wigner mass selection.
 431   bool setupMasses();
 432
 433   // Select fixed or Breit-Wigner-distributed masses.
 434   bool trialMasses();
 435
 436 };
 437
 438 //==========================================================================
 439
 440 // A derived class with 2 -> 3 kinematics 1 + 2 -> 3 + 4 + 5 set up in
 441 // y3, y4, y5, pT2_3, pT2_5, phi_3 and phi_5, and with R separation cut.
 442 // Intended specifically for (essentially massless) 2 -> 3 QCD processes.
 443
 444 class PhaseSpace2to3yyycyl : public PhaseSpace {
 445
 446 public:
 447
 448   // Constructor.
 449   PhaseSpace2to3yyycyl() {}
 450
 451   // Optimize subsequent kinematics selection.
 452   virtual bool setupSampling();
 453
 454   // Construct the trial kinematics.
 455   virtual bool trialKin(bool inEvent = true, bool = false);
 456
 457   // Construct the final event kinematics.
 458   virtual bool finalKin();
 459
 460 private:
 461
 462   // Phase space cuts specifically for 2 -> 3 QCD processes.
 463   double pTHat3Min, pTHat3Max, pTHat5Min, pTHat5Max, RsepMin, R2sepMin;
 464   bool   hasBaryonBeams;
 465
 466   // Event kinematics choices.
 467   double pT3Min, pT3Max, pT5Min, pT5Max, y3Max, y4Max, y5Max,
 468          pT3, pT4, pT5, phi3, phi4, phi5, y3, y4, y5, dphi;
 469   Vec4   pInSum;
 470
 471 };
 472
 473 //==========================================================================
 474
 475 // A derived class for Les Houches events.
 476
 477 class PhaseSpaceLHA : public PhaseSpace {
 478
 479 public:
 480
 481   // Constructor.
 482   PhaseSpaceLHA() {idProcSave = 0;}
 483
 484   // Find maximal cross section for comparison with internal processes.
 485   virtual bool setupSampling();
 486
 487   // Construct the next process, by interface to Les Houches class.
 488   virtual bool trialKin( bool , bool repeatSame = false);
 489
 490   // Set scale, alpha_s and alpha_em if not done.
 491   virtual bool finalKin() {sigmaProcessPtr->setScale(); return true;}
 492
 493   // For Les Houches with negative event weight needs
 494   virtual double sigmaSumSigned() const {return sigmaSgn;}
 495
 496 private:
 497
 498   // Constants.
 499   static const double CONVERTPB2MB;
 500
 501   // Local properties.
 502   int    strategy, stratAbs, nProc, idProcSave;
 503   double xMaxAbsSum, xSecSgnSum, sigmaSgn;
 504   vector<int>    idProc;
 505   vector<double> xMaxAbsProc;
 506
 507 };
 508
 509 //==========================================================================
 510
 511 } // end namespace Pythia8
 512
 513 #endif // Pythia8_PhaseSpace_H
 514