[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8130 / src / BeamShape.cxx

// BeamShape.cc is a part of the PYTHIA event generator.
// Copyright (C) 2008 Torbjorn Sjostrand.
// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
// Please respect the MCnet Guidelines, see GUIDELINES for details.

// Function definitions (not found in the header) for the BeamShape class. 

#include "BeamShape.h"

namespace Pythia8 {
 
//**************************************************************************

// The BeamShape class.

//*********

// Initialize beam parameters.

void BeamShape::init() { 

  // Main flags.
  allowMomentumSpread = Settings::flag("Beams:allowMomentumSpread");
  allowVertexSpread   = Settings::flag("Beams:allowVertexSpread");

  // Parameters for beam A momentum spread.
  sigmaPxA            = Settings::parm("Beams:sigmaPxA");
  sigmaPyA            = Settings::parm("Beams:sigmaPyA");
  sigmaPzA            = Settings::parm("Beams:sigmaPzA");
  maxDevA             = Settings::parm("Beams:maxDevA");

  // Parameters for beam B momentum spread.
  sigmaPxB            = Settings::parm("Beams:sigmaPxB");
  sigmaPyB            = Settings::parm("Beams:sigmaPyB");
  sigmaPzB            = Settings::parm("Beams:sigmaPzB");
  maxDevB             = Settings::parm("Beams:maxDevB");
 
  // Parameters for beam vertex spread.
  sigmaVertexX        = Settings::parm("Beams:sigmaVertexX");
  sigmaVertexY        = Settings::parm("Beams:sigmaVertexY");
  sigmaVertexZ        = Settings::parm("Beams:sigmaVertexZ");
  maxDevVertex        = Settings::parm("Beams:maxDevVertex");
  sigmaTime           = Settings::parm("Beams:sigmaTime");
  maxDevTime          = Settings::parm("Beams:maxDevTime");
 
  // Parameters for beam vertex offset.
  offsetX             = Settings::parm("Beams:offsetVertexX"); 
  offsetY             = Settings::parm("Beams:offsetVertexY");
  offsetZ             = Settings::parm("Beams:offsetVertexZ");
  offsetT             = Settings::parm("Beams:offsetTime");

}

//*********

// Set the two beam momentum deviations and the beam vertex.

void BeamShape::pick() { 

  // Reset all values.
  deltaPxA = deltaPyA = deltaPzA = deltaPxB = deltaPyB = deltaPzB
    = vertexX = vertexY = vertexZ = vertexT = 0.;

  // Set beam A momentum deviation by a three-dimensional Gaussian.
  if (allowMomentumSpread) {
    double totalDev, gauss;
    do {
      totalDev = 0.;
      if (sigmaPxA > 0.) {
        gauss     = Rndm::gauss();
        deltaPxA  = sigmaPxA * gauss;
        totalDev += gauss * gauss; 
      }
      if (sigmaPyA > 0.) {
        gauss     = Rndm::gauss();
        deltaPyA  = sigmaPyA * gauss;
        totalDev += gauss * gauss; 
      }
      if (sigmaPzA > 0.) {
        gauss     = Rndm::gauss();
        deltaPzA  = sigmaPzA * gauss;
        totalDev += gauss * gauss; 
      }
    } while (totalDev > maxDevA * maxDevA); 

    // Set beam B momentum deviation by a three-dimensional Gaussian.
    do {
      totalDev = 0.;
      if (sigmaPxB > 0.) {
        gauss     = Rndm::gauss();
        deltaPxB  = sigmaPxB * gauss;
        totalDev += gauss * gauss; 
      }
      if (sigmaPyB > 0.) {
        gauss     = Rndm::gauss();
        deltaPyB  = sigmaPyB * gauss;
        totalDev += gauss * gauss; 
      }
      if (sigmaPzB > 0.) {
        gauss     = Rndm::gauss();
        deltaPzB  = sigmaPzB * gauss;
        totalDev += gauss * gauss; 
      }
    } while (totalDev > maxDevB * maxDevB); 
  }  

  // Set beam vertex location by a three-dimensional Gaussian.
  if (allowVertexSpread) {
    double totalDev, gauss;
    do {
      totalDev = 0.;
      if (sigmaVertexX > 0.) {
        gauss     = Rndm::gauss();
        vertexX   = sigmaVertexX * gauss;
        totalDev += gauss * gauss; 
      }
      if (sigmaVertexY > 0.) {
        gauss     = Rndm::gauss();
        vertexY   = sigmaVertexY * gauss;
        totalDev += gauss * gauss; 
      }
      if (sigmaVertexZ > 0.) {
        gauss     = Rndm::gauss();
        vertexZ   = sigmaVertexZ * gauss;
        totalDev += gauss * gauss; 
      }
    } while (totalDev > maxDevVertex * maxDevVertex);

    // Set beam collision time by a Gaussian.
    if (sigmaTime > 0.) {
      do gauss    = Rndm::gauss();
      while (abs(gauss) > maxDevTime);
      vertexT     = sigmaTime * gauss;
    }

    // Add offset to beam vertex.
    vertexX      += offsetX;
    vertexY      += offsetY;
    vertexZ      += offsetZ;
    vertexT      += offsetT;
  }  

}
 
//**************************************************************************

} // end namespace Pythia8
Commit	Line	Data
5ad4eb21	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