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

// BoseEinstein.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 BoseEinsten class.

#include "BoseEinstein.h"

namespace Pythia8 {

//**************************************************************************

// The BoseEinstein class.

//*********

// Constants: could be changed here if desired, but normally should not.
// These are of technical nature, as described for each.

// Enumeration of id codes and table for particle species considered.
const int    BoseEinstein::IDHADRON[9] = { 211, -211, 111, 321, -321, 
                                           130,  310, 221, 331 };
const int    BoseEinstein::ITABLE[9] = { 0, 0, 0, 1, 1, 1, 1, 2, 3 };

// Distance between table entries, normalized to min( 2*mass, QRef).
const double BoseEinstein::STEPSIZE  = 0.05;

// Skip shift for two extremely close particles, to avoid instabilities. 
const double BoseEinstein::Q2MIN     = 1e-8;

// Parameters of energy compensation procedure: maximally allowed 
// relative energy error, iterative stepsize, and number of iterations. 
const double BoseEinstein::COMPRELERR = 1e-10;
const double BoseEinstein::COMPFACMAX = 1000.;
const int    BoseEinstein::NCOMPSTEP  = 10;

//*********

// Find settings. Precalculate table used to find momentum shifts.

bool BoseEinstein::init(Info* infoPtrIn) {

  // Save pointer.
  infoPtr  = infoPtrIn;

  // Main flags.
  doPion   = Settings::flag("BoseEinstein:Pion");
  doKaon   = Settings::flag("BoseEinstein:Kaon");
  doEta    = Settings::flag("BoseEinstein:Eta");

  // Shape of Bose-Einstein enhancement/suppression.
  lambda   = Settings::parm("BoseEinstein:lambda");
  QRef     = Settings::parm("BoseEinstein:QRef");

  // Multiples and inverses (= "radii") of distance parameters in Q-space.
  QRef2    = 2. * QRef;
  QRef3    = 3. * QRef;
  R2Ref    = 1. / (QRef * QRef);
  R2Ref2   = 1. / (QRef2 * QRef2);
  R2Ref3   = 1. / (QRef3 * QRef3);

  // Masses of particles with Bose-Einstein implemented.
  for (int iSpecies = 0; iSpecies < 9; ++iSpecies) 
    mHadron[iSpecies] = ParticleDataTable::m0( IDHADRON[iSpecies] );

  // Pair pi, K, eta and eta' masses for use in tables.
  mPair[0] = 2. * mHadron[0];   
  mPair[1] = 2. * mHadron[3];   
  mPair[2] = 2. * mHadron[7];   
  mPair[3] = 2. * mHadron[8]; 
  
  // Loop over the four required tables. Local variables.
  double Qnow, Q2now, centerCorr;
  for (int iTab = 0; iTab < 4; ++iTab) {
    m2Pair[iTab]      = mPair[iTab] * mPair[iTab]; 
    
    // Step size and number of steps in normal table.
    deltaQ[iTab]      = STEPSIZE * min(mPair[iTab], QRef);     
    nStep[iTab]       = min( 199, 1 + int(3. * QRef / deltaQ[iTab]) );
    maxQ[iTab]        = (nStep[iTab] - 0.1) * deltaQ[iTab];  
    centerCorr        = deltaQ[iTab] * deltaQ[iTab] / 12.;

    // Construct normal table recursively in Q space.
    shift[iTab][0]    = 0.;
    for (int i = 1; i <= nStep[iTab]; ++i) {
      Qnow            = deltaQ[iTab] * (i - 0.5);
      Q2now           = Qnow * Qnow;
      shift[iTab][i]  = shift[iTab][i - 1] + exp(-Q2now * R2Ref) 
        * deltaQ[iTab] * (Q2now + centerCorr) / sqrt(Q2now + m2Pair[iTab]);
    }       

    // Step size and number of steps in compensation table.
    deltaQ3[iTab]     = STEPSIZE * min(mPair[iTab], QRef3);     
    nStep3[iTab]      = min( 199, 1 + int(9. * QRef / deltaQ3[iTab]) );
    maxQ3[iTab]       = (nStep3[iTab] - 0.1) * deltaQ3[iTab];  
    centerCorr        = deltaQ3[iTab] * deltaQ3[iTab] / 12.;

    // Construct compensation table recursively in Q space.
    shift3[iTab][0]   = 0.;
    for (int i = 1; i <= nStep3[iTab]; ++i) {
      Qnow            = deltaQ3[iTab] * (i - 0.5);
      Q2now           = Qnow * Qnow;
      shift3[iTab][i] = shift3[iTab][i - 1] + exp(-Q2now * R2Ref3) 
        * deltaQ3[iTab] * (Q2now + centerCorr) / sqrt(Q2now + m2Pair[iTab]);
    }       

  }

  // Done.
  return true;

}

//*********

// Perform Bose-Einstein corrections on an event.

bool BoseEinstein::shiftEvent( Event& event) {

  // Reset list of identical particles.
  hadronBE.resize(0);

  // Loop over all hadron species with BE effects.
  nStored[0] = 0;
  for (int iSpecies = 0; iSpecies < 9; ++iSpecies) {
    nStored[iSpecies + 1] = nStored[iSpecies]; 
    if (!doPion && iSpecies <= 2) continue;
    if (!doKaon && iSpecies >= 3 && iSpecies <= 6) continue;
    if (!doEta  && iSpecies >= 7) continue;
   
    // Properties of current hadron species. 
    int idNow = IDHADRON[ iSpecies ];
    int iTab  = ITABLE[ iSpecies ];
   
    // Loop through event record to store copies of current species.
    for (int i = 0; i < event.size(); ++i) 
      if ( event[i].id() == idNow && event[i].isFinal() ) 
        hadronBE.push_back( 
          BoseEinsteinHadron( idNow, i, event[i].p(), event[i].m() ) ); 
    nStored[iSpecies + 1] = hadronBE.size();

    // Loop through pairs of identical particles and find shifts. 
    for (int i1 = nStored[iSpecies]; i1 < nStored[iSpecies+1] - 1; ++i1) 
    for (int i2 = i1 + 1; i2 < nStored[iSpecies+1]; ++i2) 
      shiftPair( i1, i2, iTab);  
  }

  // Must have at least two pairs to carry out compensation.
  if (nStored[9] < 2) return true; 

  // Shift momenta and recalculate energies. 
  double eSumOriginal = 0.;
  double eSumShifted  = 0.;
  double eDiffByComp  = 0.;
  for (int i = 0; i < nStored[9]; ++i) {
    eSumOriginal  += hadronBE[i].p.e();
    hadronBE[i].p += hadronBE[i].pShift;
    hadronBE[i].p.e( sqrt( hadronBE[i].p.pAbs2() + hadronBE[i].m2 ) );  
    eSumShifted   += hadronBE[i].p.e();
    eDiffByComp   += dot3( hadronBE[i].pComp, hadronBE[i].p) 
                     / hadronBE[i].p.e(); 
  }

  // Iterate compensation shift until convergence.
  int iStep = 0;
  while ( abs(eSumShifted - eSumOriginal) > COMPRELERR * eSumOriginal
    && abs(eSumShifted - eSumOriginal) < COMPFACMAX * abs(eDiffByComp)  
    && iStep < NCOMPSTEP ) {
    ++iStep;
    double compFac   = (eSumOriginal - eSumShifted) / eDiffByComp;
    eSumShifted      = 0.;
    eDiffByComp      = 0.;
    for (int i = 0; i < nStored[9]; ++i) {
      hadronBE[i].p += compFac * hadronBE[i].pComp;
      hadronBE[i].p.e( sqrt( hadronBE[i].p.pAbs2() + hadronBE[i].m2 ) );  
      eSumShifted   += hadronBE[i].p.e();
      eDiffByComp   += dot3( hadronBE[i].pComp, hadronBE[i].p) 
                       / hadronBE[i].p.e();
    } 
  }

  // Error if no convergence, and then return without doing BE shift.
  // However, not grave enough to kill event, so return true. 
  if ( abs(eSumShifted - eSumOriginal) > COMPRELERR * eSumOriginal ) {
    infoPtr->errorMsg("Warning in BoseEinstein::shiftEvent: "
      "no consistent BE shift topology found, so skip BE"); 
    return true;
  }
    
  // Store new particle copies with shifted momenta.
  for (int i = 0; i < nStored[9]; ++i) {
    int iNew = event.copy( hadronBE[i].iPos, 99);
    event[ iNew ].p( hadronBE[i].p );  
  }

  // Done.
  return true;

}

//*********

// Calculate shift and (unnormalized) compensation for pair.

void BoseEinstein::shiftPair( int i1, int i2, int iTab) {

  // Calculate old relative momentum.
  double Q2old = m2(hadronBE[i1].p, hadronBE[i2].p) - m2Pair[iTab];
  if (Q2old < Q2MIN) return;
  double Qold  = sqrt(Q2old);
  double psFac = sqrt(Q2old + m2Pair[iTab]) / Q2old;
     
  // Calculate new relative momentum for normal shift.
  double Qmove = 0.;
  if (Qold < deltaQ[iTab]) Qmove = Qold / 3.;
  else if (Qold < maxQ[iTab]) {
    double realQbin = Qold / deltaQ[iTab];
    int    intQbin  = int( realQbin );
    double inter    = (pow3(realQbin) - pow3(intQbin)) 
      / (3 * intQbin * (intQbin + 1) + 1);
    Qmove = ( shift[iTab][intQbin] + inter * (shift[iTab][intQbin + 1]  
      - shift[iTab][intQbin]) ) * psFac;
  } 
  else Qmove = shift[iTab][nStep[iTab]] * psFac; 
  double Q2new = Q2old * pow( Qold / (Qold + 3. * lambda * Qmove), 2. / 3.);

  // Calculate corresponding three-momentum shift.
  double Q2Diff    = Q2new - Q2old;
  double p2DiffAbs = (hadronBE[i1].p - hadronBE[i2].p).pAbs2();
  double p2AbsDiff = hadronBE[i1].p.pAbs2() - hadronBE[i2].p.pAbs2();
  double eSum      = hadronBE[i1].p.e() + hadronBE[i2].p.e();
  double eDiff     = hadronBE[i1].p.e() - hadronBE[i2].p.e();
  double sumQ2E    = Q2Diff + eSum * eSum;
  double rootA     = eSum * eDiff * p2AbsDiff - p2DiffAbs * sumQ2E;
  double rootB     = p2DiffAbs * sumQ2E - p2AbsDiff * p2AbsDiff;
  double factor    = 0.5 * ( rootA + sqrtpos(rootA * rootA 
    + Q2Diff * (sumQ2E - eDiff * eDiff) * rootB) ) / rootB;

  // Add shifts to sum. (Energy component dummy.)
  Vec4   pDiff     = factor * (hadronBE[i1].p - hadronBE[i2].p);
  hadronBE[i1].pShift += pDiff;  
  hadronBE[i2].pShift -= pDiff;  

  // Calculate new relative momentum for compensation shift.
  double Qmove3 = 0.;
  if (Qold < deltaQ3[iTab]) Qmove3 = Qold / 3.;
  else if (Qold < maxQ3[iTab]) {
    double realQbin = Qold / deltaQ3[iTab];
    int    intQbin  = int( realQbin );
    double inter    = (pow3(realQbin) - pow3(intQbin)) 
      / (3 * intQbin * (intQbin + 1) + 1);
    Qmove3 = ( shift3[iTab][intQbin] + inter * (shift3[iTab][intQbin + 1]  
      - shift3[iTab][intQbin]) ) * psFac;
  } 
  else Qmove3 = shift3[iTab][nStep3[iTab]] *psFac; 
  double Q2new3 = Q2old * pow( Qold / (Qold + 3. * lambda * Qmove3), 2. / 3.);

  // Calculate corresponding three-momentum shift.
  Q2Diff    = Q2new3 - Q2old;
  sumQ2E    = Q2Diff + eSum * eSum;
  rootA     = eSum * eDiff * p2AbsDiff - p2DiffAbs * sumQ2E;
  rootB     = p2DiffAbs * sumQ2E - p2AbsDiff * p2AbsDiff;
  factor    = 0.5 * ( rootA + sqrtpos(rootA * rootA 
    + Q2Diff * (sumQ2E - eDiff * eDiff) * rootB) ) / rootB;

  // Extra dampening factor to go from BE_3 to BE_32.
  factor   *= 1. - exp(-Q2old * R2Ref2); 

  // Add shifts to sum. (Energy component dummy.)
  pDiff     = factor * (hadronBE[i1].p - hadronBE[i2].p);
  hadronBE[i1].pComp += pDiff;  
  hadronBE[i2].pComp -= pDiff;  

}

//**************************************************************************

} // end namespace Pythia8
Commit	Line	Data
5ad4eb21	1	// BoseEinstein.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 BoseEinsten class.
	7
	8	#include "BoseEinstein.h"
	9
	10	namespace Pythia8 {
	11
	12	//**************************************************************************
	13
	14	// The BoseEinstein class.
	15
	16	//*********
	17
	18	// Constants: could be changed here if desired, but normally should not.
	19	// These are of technical nature, as described for each.
	20
	21	// Enumeration of id codes and table for particle species considered.
	22	const int BoseEinstein::IDHADRON[9] = { 211, -211, 111, 321, -321,
	23	130, 310, 221, 331 };
	24	const int BoseEinstein::ITABLE[9] = { 0, 0, 0, 1, 1, 1, 1, 2, 3 };
	25
	26	// Distance between table entries, normalized to min( 2*mass, QRef).
	27	const double BoseEinstein::STEPSIZE = 0.05;
	28
	29	// Skip shift for two extremely close particles, to avoid instabilities.
	30	const double BoseEinstein::Q2MIN = 1e-8;
	31
	32	// Parameters of energy compensation procedure: maximally allowed
	33	// relative energy error, iterative stepsize, and number of iterations.
	34	const double BoseEinstein::COMPRELERR = 1e-10;
	35	const double BoseEinstein::COMPFACMAX = 1000.;
	36	const int BoseEinstein::NCOMPSTEP = 10;
	37
	38	//*********
	39
	40	// Find settings. Precalculate table used to find momentum shifts.
	41
	42	bool BoseEinstein::init(Info* infoPtrIn) {
	43
	44	// Save pointer.
	45	infoPtr = infoPtrIn;
	46
	47	// Main flags.
	48	doPion = Settings::flag("BoseEinstein:Pion");
	49	doKaon = Settings::flag("BoseEinstein:Kaon");
	50	doEta = Settings::flag("BoseEinstein:Eta");
	51
	52	// Shape of Bose-Einstein enhancement/suppression.
	53	lambda = Settings::parm("BoseEinstein:lambda");
	54	QRef = Settings::parm("BoseEinstein:QRef");
	55
	56	// Multiples and inverses (= "radii") of distance parameters in Q-space.
	57	QRef2 = 2. * QRef;
	58	QRef3 = 3. * QRef;
	59	R2Ref = 1. / (QRef * QRef);
	60	R2Ref2 = 1. / (QRef2 * QRef2);
	61	R2Ref3 = 1. / (QRef3 * QRef3);
	62
	63	// Masses of particles with Bose-Einstein implemented.
	64	for (int iSpecies = 0; iSpecies < 9; ++iSpecies)
65	mHadron[iSpecies] = ParticleDataTable::m0( IDHADRON[iSpecies] );
66
67	// Pair pi, K, eta and eta' masses for use in tables.
68	mPair[0] = 2. * mHadron[0];
69	mPair[1] = 2. * mHadron[3];
70	mPair[2] = 2. * mHadron[7];
71	mPair[3] = 2. * mHadron[8];
72
73	// Loop over the four required tables. Local variables.
74	double Qnow, Q2now, centerCorr;
75	for (int iTab = 0; iTab < 4; ++iTab) {
76	m2Pair[iTab] = mPair[iTab] * mPair[iTab];
77
78	// Step size and number of steps in normal table.
79	deltaQ[iTab] = STEPSIZE * min(mPair[iTab], QRef);
80	nStep[iTab] = min( 199, 1 + int(3. * QRef / deltaQ[iTab]) );
81	maxQ[iTab] = (nStep[iTab] - 0.1) * deltaQ[iTab];
82	centerCorr = deltaQ[iTab] * deltaQ[iTab] / 12.;
83
84	// Construct normal table recursively in Q space.
85	shift[iTab][0] = 0.;
86	for (int i = 1; i <= nStep[iTab]; ++i) {
87	Qnow = deltaQ[iTab] * (i - 0.5);
88	Q2now = Qnow * Qnow;
89	shift[iTab][i] = shift[iTab][i - 1] + exp(-Q2now * R2Ref)
90	* deltaQ[iTab] * (Q2now + centerCorr) / sqrt(Q2now + m2Pair[iTab]);
91	}
92
93	// Step size and number of steps in compensation table.
94	deltaQ3[iTab] = STEPSIZE * min(mPair[iTab], QRef3);
95	nStep3[iTab] = min( 199, 1 + int(9. * QRef / deltaQ3[iTab]) );
96	maxQ3[iTab] = (nStep3[iTab] - 0.1) * deltaQ3[iTab];
97	centerCorr = deltaQ3[iTab] * deltaQ3[iTab] / 12.;
98
99	// Construct compensation table recursively in Q space.
100	shift3[iTab][0] = 0.;
101	for (int i = 1; i <= nStep3[iTab]; ++i) {
102	Qnow = deltaQ3[iTab] * (i - 0.5);
103	Q2now = Qnow * Qnow;
104	shift3[iTab][i] = shift3[iTab][i - 1] + exp(-Q2now * R2Ref3)
105	* deltaQ3[iTab] * (Q2now + centerCorr) / sqrt(Q2now + m2Pair[iTab]);
106	}
107
108	}
109
110	// Done.
111	return true;
112
113	}
114
115	//*********
116
117	// Perform Bose-Einstein corrections on an event.
118
119	bool BoseEinstein::shiftEvent( Event& event) {
120
121	// Reset list of identical particles.
122	hadronBE.resize(0);
123
124	// Loop over all hadron species with BE effects.
125	nStored[0] = 0;
126	for (int iSpecies = 0; iSpecies < 9; ++iSpecies) {
127	nStored[iSpecies + 1] = nStored[iSpecies];
128	if (!doPion && iSpecies <= 2) continue;
129	if (!doKaon && iSpecies >= 3 && iSpecies <= 6) continue;
130	if (!doEta && iSpecies >= 7) continue;
131
132	// Properties of current hadron species.
133	int idNow = IDHADRON[ iSpecies ];
134	int iTab = ITABLE[ iSpecies ];
135
136	// Loop through event record to store copies of current species.
137	for (int i = 0; i < event.size(); ++i)
138	if ( event[i].id() == idNow && event[i].isFinal() )
139	hadronBE.push_back(
140	BoseEinsteinHadron( idNow, i, event[i].p(), event[i].m() ) );
141	nStored[iSpecies + 1] = hadronBE.size();
142
143	// Loop through pairs of identical particles and find shifts.
144	for (int i1 = nStored[iSpecies]; i1 < nStored[iSpecies+1] - 1; ++i1)
145	for (int i2 = i1 + 1; i2 < nStored[iSpecies+1]; ++i2)
146	shiftPair( i1, i2, iTab);
147	}
148
149	// Must have at least two pairs to carry out compensation.
150	if (nStored[9] < 2) return true;
151
152	// Shift momenta and recalculate energies.
153	double eSumOriginal = 0.;
154	double eSumShifted = 0.;
155	double eDiffByComp = 0.;
156	for (int i = 0; i < nStored[9]; ++i) {
157	eSumOriginal += hadronBE[i].p.e();
158	hadronBE[i].p += hadronBE[i].pShift;
159	hadronBE[i].p.e( sqrt( hadronBE[i].p.pAbs2() + hadronBE[i].m2 ) );
160	eSumShifted += hadronBE[i].p.e();
161	eDiffByComp += dot3( hadronBE[i].pComp, hadronBE[i].p)
162	/ hadronBE[i].p.e();
163	}
164
165	// Iterate compensation shift until convergence.
166	int iStep = 0;
167	while ( abs(eSumShifted - eSumOriginal) > COMPRELERR * eSumOriginal
168	&& abs(eSumShifted - eSumOriginal) < COMPFACMAX * abs(eDiffByComp)
169	&& iStep < NCOMPSTEP ) {
170	++iStep;
171	double compFac = (eSumOriginal - eSumShifted) / eDiffByComp;
172	eSumShifted = 0.;
173	eDiffByComp = 0.;
174	for (int i = 0; i < nStored[9]; ++i) {
175	hadronBE[i].p += compFac * hadronBE[i].pComp;
176	hadronBE[i].p.e( sqrt( hadronBE[i].p.pAbs2() + hadronBE[i].m2 ) );
177	eSumShifted += hadronBE[i].p.e();
178	eDiffByComp += dot3( hadronBE[i].pComp, hadronBE[i].p)
179	/ hadronBE[i].p.e();
180	}
181	}
182
183	// Error if no convergence, and then return without doing BE shift.
184	// However, not grave enough to kill event, so return true.
185	if ( abs(eSumShifted - eSumOriginal) > COMPRELERR * eSumOriginal ) {
186	infoPtr->errorMsg("Warning in BoseEinstein::shiftEvent: "
187	"no consistent BE shift topology found, so skip BE");
188	return true;
189	}
190
191	// Store new particle copies with shifted momenta.
192	for (int i = 0; i < nStored[9]; ++i) {
193	int iNew = event.copy( hadronBE[i].iPos, 99);
194	event[ iNew ].p( hadronBE[i].p );
195	}
196
197	// Done.
198	return true;
199
200	}
201
202	//*********
203
204	// Calculate shift and (unnormalized) compensation for pair.
205
206	void BoseEinstein::shiftPair( int i1, int i2, int iTab) {
207
208	// Calculate old relative momentum.
209	double Q2old = m2(hadronBE[i1].p, hadronBE[i2].p) - m2Pair[iTab];
210	if (Q2old < Q2MIN) return;
211	double Qold = sqrt(Q2old);
212	double psFac = sqrt(Q2old + m2Pair[iTab]) / Q2old;
213
214	// Calculate new relative momentum for normal shift.
215	double Qmove = 0.;
216	if (Qold < deltaQ[iTab]) Qmove = Qold / 3.;
217	else if (Qold < maxQ[iTab]) {
218	double realQbin = Qold / deltaQ[iTab];
219	int intQbin = int( realQbin );
220	double inter = (pow3(realQbin) - pow3(intQbin))
221	/ (3 * intQbin * (intQbin + 1) + 1);
222	Qmove = ( shift[iTab][intQbin] + inter * (shift[iTab][intQbin + 1]
223	- shift[iTab][intQbin]) ) * psFac;
224	}
225	else Qmove = shift[iTab][nStep[iTab]] * psFac;
226	double Q2new = Q2old * pow( Qold / (Qold + 3. * lambda * Qmove), 2. / 3.);
227
228	// Calculate corresponding three-momentum shift.
229	double Q2Diff = Q2new - Q2old;
230	double p2DiffAbs = (hadronBE[i1].p - hadronBE[i2].p).pAbs2();
231	double p2AbsDiff = hadronBE[i1].p.pAbs2() - hadronBE[i2].p.pAbs2();
232	double eSum = hadronBE[i1].p.e() + hadronBE[i2].p.e();
233	double eDiff = hadronBE[i1].p.e() - hadronBE[i2].p.e();
234	double sumQ2E = Q2Diff + eSum * eSum;
235	double rootA = eSum * eDiff * p2AbsDiff - p2DiffAbs * sumQ2E;
236	double rootB = p2DiffAbs * sumQ2E - p2AbsDiff * p2AbsDiff;
237	double factor = 0.5 * ( rootA + sqrtpos(rootA * rootA
238	+ Q2Diff * (sumQ2E - eDiff * eDiff) * rootB) ) / rootB;
239
240	// Add shifts to sum. (Energy component dummy.)
241	Vec4 pDiff = factor * (hadronBE[i1].p - hadronBE[i2].p);
242	hadronBE[i1].pShift += pDiff;
243	hadronBE[i2].pShift -= pDiff;
244
245	// Calculate new relative momentum for compensation shift.
246	double Qmove3 = 0.;
247	if (Qold < deltaQ3[iTab]) Qmove3 = Qold / 3.;
248	else if (Qold < maxQ3[iTab]) {
249	double realQbin = Qold / deltaQ3[iTab];
250	int intQbin = int( realQbin );
251	double inter = (pow3(realQbin) - pow3(intQbin))
252	/ (3 * intQbin * (intQbin + 1) + 1);
253	Qmove3 = ( shift3[iTab][intQbin] + inter * (shift3[iTab][intQbin + 1]
254	- shift3[iTab][intQbin]) ) * psFac;
255	}
256	else Qmove3 = shift3[iTab][nStep3[iTab]] *psFac;
257	double Q2new3 = Q2old * pow( Qold / (Qold + 3. * lambda * Qmove3), 2. / 3.);
258
259	// Calculate corresponding three-momentum shift.
260	Q2Diff = Q2new3 - Q2old;
261	sumQ2E = Q2Diff + eSum * eSum;
262	rootA = eSum * eDiff * p2AbsDiff - p2DiffAbs * sumQ2E;
263	rootB = p2DiffAbs * sumQ2E - p2AbsDiff * p2AbsDiff;
264	factor = 0.5 * ( rootA + sqrtpos(rootA * rootA
265	+ Q2Diff * (sumQ2E - eDiff * eDiff) * rootB) ) / rootB;
266
267	// Extra dampening factor to go from BE_3 to BE_32.
268	factor = 1. - exp(-Q2old R2Ref2);
269
270	// Add shifts to sum. (Energy component dummy.)
271	pDiff = factor * (hadronBE[i1].p - hadronBE[i2].p);
272	hadronBE[i1].pComp += pDiff;
273	hadronBE[i2].pComp -= pDiff;
274
275	}
276
277	//**************************************************************************
278
279	} // end namespace Pythia8
280