PYTHIA8/pythia8130/src/BeamShape.cxx

   1 // BeamShape.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2008 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 BeamShape class.
   7
   8 #include "BeamShape.h"
   9
  10 namespace Pythia8 {
  11
  12 //**************************************************************************
  13
  14 // The BeamShape class.
  15
  16 //*********
  17
  18 // Initialize beam parameters.
  19
  20 void BeamShape::init() {
  21
  22   // Main flags.
  23   allowMomentumSpread = Settings::flag("Beams:allowMomentumSpread");
  24   allowVertexSpread   = Settings::flag("Beams:allowVertexSpread");
  25
  26   // Parameters for beam A momentum spread.
  27   sigmaPxA            = Settings::parm("Beams:sigmaPxA");
  28   sigmaPyA            = Settings::parm("Beams:sigmaPyA");
  29   sigmaPzA            = Settings::parm("Beams:sigmaPzA");
  30   maxDevA             = Settings::parm("Beams:maxDevA");
  31
  32   // Parameters for beam B momentum spread.
  33   sigmaPxB            = Settings::parm("Beams:sigmaPxB");
  34   sigmaPyB            = Settings::parm("Beams:sigmaPyB");
  35   sigmaPzB            = Settings::parm("Beams:sigmaPzB");
  36   maxDevB             = Settings::parm("Beams:maxDevB");
  37
  38   // Parameters for beam vertex spread.
  39   sigmaVertexX        = Settings::parm("Beams:sigmaVertexX");
  40   sigmaVertexY        = Settings::parm("Beams:sigmaVertexY");
  41   sigmaVertexZ        = Settings::parm("Beams:sigmaVertexZ");
  42   maxDevVertex        = Settings::parm("Beams:maxDevVertex");
  43   sigmaTime           = Settings::parm("Beams:sigmaTime");
  44   maxDevTime          = Settings::parm("Beams:maxDevTime");
  45
  46   // Parameters for beam vertex offset.
  47   offsetX             = Settings::parm("Beams:offsetVertexX");
  48   offsetY             = Settings::parm("Beams:offsetVertexY");
  49   offsetZ             = Settings::parm("Beams:offsetVertexZ");
  50   offsetT             = Settings::parm("Beams:offsetTime");
  51
  52 }
  53
  54 //*********
  55
  56 // Set the two beam momentum deviations and the beam vertex.
  57
  58 void BeamShape::pick() {
  59
  60   // Reset all values.
  61   deltaPxA = deltaPyA = deltaPzA = deltaPxB = deltaPyB = deltaPzB
  62     = vertexX = vertexY = vertexZ = vertexT = 0.;
  63
  64   // Set beam A momentum deviation by a three-dimensional Gaussian.
  65   if (allowMomentumSpread) {
  66     double totalDev, gauss;
  67     do {
  68       totalDev = 0.;
  69       if (sigmaPxA > 0.) {
  70         gauss     = Rndm::gauss();
  71         deltaPxA  = sigmaPxA * gauss;
  72         totalDev += gauss * gauss;
  73       }
  74       if (sigmaPyA > 0.) {
  75         gauss     = Rndm::gauss();
  76         deltaPyA  = sigmaPyA * gauss;
  77         totalDev += gauss * gauss;
  78       }
  79       if (sigmaPzA > 0.) {
  80         gauss     = Rndm::gauss();
  81         deltaPzA  = sigmaPzA * gauss;
  82         totalDev += gauss * gauss;
  83       }
  84     } while (totalDev > maxDevA * maxDevA);
  85
  86     // Set beam B momentum deviation by a three-dimensional Gaussian.
  87     do {
  88       totalDev = 0.;
  89       if (sigmaPxB > 0.) {
  90         gauss     = Rndm::gauss();
  91         deltaPxB  = sigmaPxB * gauss;
  92         totalDev += gauss * gauss;
  93       }
  94       if (sigmaPyB > 0.) {
  95         gauss     = Rndm::gauss();
  96         deltaPyB  = sigmaPyB * gauss;
  97         totalDev += gauss * gauss;
  98       }
  99       if (sigmaPzB > 0.) {
 100         gauss     = Rndm::gauss();
 101         deltaPzB  = sigmaPzB * gauss;
 102         totalDev += gauss * gauss;
 103       }
 104     } while (totalDev > maxDevB * maxDevB);
 105   }
 106
 107   // Set beam vertex location by a three-dimensional Gaussian.
 108   if (allowVertexSpread) {
 109     double totalDev, gauss;
 110     do {
 111       totalDev = 0.;
 112       if (sigmaVertexX > 0.) {
 113         gauss     = Rndm::gauss();
 114         vertexX   = sigmaVertexX * gauss;
 115         totalDev += gauss * gauss;
 116       }
 117       if (sigmaVertexY > 0.) {
 118         gauss     = Rndm::gauss();
 119         vertexY   = sigmaVertexY * gauss;
 120         totalDev += gauss * gauss;
 121       }
 122       if (sigmaVertexZ > 0.) {
 123         gauss     = Rndm::gauss();
 124         vertexZ   = sigmaVertexZ * gauss;
 125         totalDev += gauss * gauss;
 126       }
 127     } while (totalDev > maxDevVertex * maxDevVertex);
 128
 129     // Set beam collision time by a Gaussian.
 130     if (sigmaTime > 0.) {
 131       do gauss    = Rndm::gauss();
 132       while (abs(gauss) > maxDevTime);
 133       vertexT     = sigmaTime * gauss;
 134     }
 135
 136     // Add offset to beam vertex.
 137     vertexX      += offsetX;
 138     vertexY      += offsetY;
 139     vertexZ      += offsetZ;
 140     vertexT      += offsetT;
 141   }
 142
 143 }
 144
 145 //**************************************************************************
 146
 147 } // end namespace Pythia8
 148