[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8170 / src / MiniStringFragmentation.cxx

// MiniStringFragmentation.cc is a part of the PYTHIA event generator.
// Copyright (C) 2012 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 .
// MiniStringFragmentation class

#include "MiniStringFragmentation.h"

namespace Pythia8 {

//==========================================================================

// The MiniStringFragmentation class.

//--------------------------------------------------------------------------

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

// Since diffractive by definition is > 1 particle, try hard. 
const int MiniStringFragmentation::NTRYDIFFRACTIVE = 200;

// After one-body fragmentation failed, try two-body once more. 
const int MiniStringFragmentation::NTRYLASTRESORT  = 100;

// Loop try to combine available endquarks to valid hadron. 
const int MiniStringFragmentation::NTRYFLAV        = 10;

//--------------------------------------------------------------------------

// Initialize and save pointers.

void MiniStringFragmentation::init(Info* infoPtrIn, Settings& settings,
   ParticleData* particleDataPtrIn, Rndm* rndmPtrIn, 
   StringFlav* flavSelPtrIn, StringPT* pTSelPtrIn, StringZ* zSelPtrIn) {

  // Save pointers.
  infoPtr         = infoPtrIn;
  particleDataPtr = particleDataPtrIn;
  rndmPtr         = rndmPtrIn;
  flavSelPtr      = flavSelPtrIn;
  pTSelPtr        = pTSelPtrIn;
  zSelPtr         = zSelPtrIn;

  // Initialize the MiniStringFragmentation class proper.
  nTryMass        = settings.mode("MiniStringFragmentation:nTry");

  // Initialize the b parameter of the z spectrum, used when joining jets.
  bLund           = zSelPtr->bAreaLund();

}

//--------------------------------------------------------------------------

// Do the fragmentation: driver routine.
  
bool MiniStringFragmentation::fragment(int iSub, ColConfig& colConfig, 
  Event& event, bool isDiff) {

  // Read in info on system to be treated.
  iParton  = colConfig[iSub].iParton;
  flav1    = FlavContainer( event[ iParton.front() ].id() );
  flav2    = FlavContainer( event[ iParton.back() ].id() ); 
  pSum     = colConfig[iSub].pSum;
  mSum     = colConfig[iSub].mass;
  m2Sum    = mSum*mSum;
  isClosed = colConfig[iSub].isClosed;

  // Do not want diffractive systems to easily collapse to one particle.
  int nTryFirst = (isDiff) ? NTRYDIFFRACTIVE : nTryMass;

  // First try to produce two particles from the system.
  if (ministring2two( nTryFirst, event)) return true;  

  // If this fails, then form one hadron and shuffle momentum.
  if (ministring2one( iSub, colConfig, event)) return true;  

  // If also this fails, then try harder to produce two particles.
  if (ministring2two( NTRYLASTRESORT, event)) return true;  

  // Else complete failure.
  infoPtr->errorMsg("Error in MiniStringFragmentation::fragment: "
      "no 1- or 2-body state found above mass threshold"); 
  return false;

}

//--------------------------------------------------------------------------

  // Attempt to produce two particles from the ministring.
  
bool MiniStringFragmentation::ministring2two( int nTry, Event& event) {

  // Properties of the produced hadrons.
  int    idHad1  = 0;
  int    idHad2  = 0;
  double mHad1   = 0.;
  double mHad2   = 0.;
  double mHadSum = 0.;

  // Allow a few attempts to find a particle pair with low enough masses.
  for (int iTry = 0; iTry < nTry; ++iTry) {

    // For closed gluon loop need to pick an initial flavour.
    if (isClosed) do {
      int idStart = flavSelPtr->pickLightQ();
      FlavContainer flavStart(idStart, 1);
      flavStart = flavSelPtr->pick( flavStart);
      flav1 = flavSelPtr->pick( flavStart);
      flav2.anti(flav1);
    } while (flav1.id == 0 || flav1.nPop > 0);
   
    // Create a new q qbar flavour to form two hadrons. 
    // Start from a diquark, if any.
    do {
      FlavContainer flav3 =
        (flav1.isDiquark() || (!flav2.isDiquark() && rndmPtr->flat() < 0.5) )  
        ? flavSelPtr->pick( flav1) : flavSelPtr->pick( flav2).anti();
      idHad1 = flavSelPtr->combine( flav1, flav3);
      idHad2 = flavSelPtr->combine( flav2, flav3.anti()); 
    } while (idHad1 == 0 || idHad2 == 0);

    // Check whether the mass sum fits inside the available phase space.  
    mHad1 = particleDataPtr->mass(idHad1);
    mHad2 = particleDataPtr->mass(idHad2);
    mHadSum = mHad1 + mHad2;
    if (mHadSum < mSum) break;
  } 
  if (mHadSum >= mSum) return false;

  // Define an effective two-parton string, by splitting intermediate
  // gluon momenta in proportion to their closeness to either endpoint.
  Vec4 pSum1 = event[ iParton.front() ].p();
  Vec4 pSum2 = event[ iParton.back() ].p();
  if (iParton.size() > 2) {
    Vec4 pEnd1 = pSum1;
    Vec4 pEnd2 = pSum2;
    Vec4 pEndSum = pEnd1 + pEnd2; 
    for (int i = 1; i < int(iParton.size()) - 1 ; ++i) {
      Vec4 pNow = event[ iParton[i] ].p();
      double ratio = (pEnd2 * pNow) / (pEndSum * pNow);
      pSum1 += ratio * pNow;
      pSum2 += (1. - ratio) * pNow;
    }
  }

  // Set up a string region based on the two effective endpoints.
  StringRegion region;
  region.setUp( pSum1, pSum2);

  // Generate an isotropic decay in the ministring rest frame, 
  // suppressed at large pT by a fragmentation pT Gaussian.
  double pAbs2 = 0.25 * ( pow2(m2Sum - mHad1*mHad1 - mHad2*mHad2)
    - pow2(2. * mHad1 * mHad2) ) / m2Sum; 
  double pT2 = 0.;
  do {
    double cosTheta = rndmPtr->flat();
    pT2 = (1. - pow2(cosTheta)) * pAbs2;
  } while (pTSelPtr->suppressPT2(pT2) < rndmPtr->flat() ); 

  // Construct the forward-backward asymmetry of the two particles.
  double mT21 = mHad1*mHad1 + pT2;
  double mT22 = mHad2*mHad2 + pT2;
  double lambda = sqrtpos( pow2(m2Sum  - mT21 - mT22) - 4. * mT21 * mT22 );
  double probReverse = 1. / (1. + exp( min( 50., bLund * lambda) ) ); 

  // Construct kinematics, as viewed in the transverse rest frame. 
  double xpz1 = 0.5 * lambda/ m2Sum;
  if (probReverse > rndmPtr->flat()) xpz1 = -xpz1; 
  double xmDiff = (mT21 - mT22) / m2Sum;
  double xe1 = 0.5 * (1. + xmDiff);
  double xe2 = 0.5 * (1. - xmDiff ); 

  // Distribute pT isotropically in angle.
  double phi = 2. * M_PI * rndmPtr->flat();
  double pT  = sqrt(pT2);
  double px  = pT * cos(phi);
  double py  = pT * sin(phi);

  // Translate this into kinematics in the string frame.
  Vec4 pHad1 = region.pHad( xe1 + xpz1, xe1 - xpz1,  px,  py);
  Vec4 pHad2 = region.pHad( xe2 - xpz1, xe2 + xpz1, -px, -py);

  // Add produced particles to the event record.
  int iFirst = event.append( idHad1, 82, iParton.front(), iParton.back(), 
    0, 0, 0, 0, pHad1, mHad1);
  int iLast = event.append( idHad2, 82, iParton.front(), iParton.back(), 
    0, 0, 0, 0, pHad2, mHad2);

  // Set decay vertex when this is displaced.
  if (event[iParton.front()].hasVertex()) {
    Vec4 vDec = event[iParton.front()].vDec();
    event[iFirst].vProd( vDec );
    event[iLast].vProd( vDec );
  }

  // Set lifetime of hadrons.
  event[iFirst].tau( event[iFirst].tau0() * rndmPtr->exp() );
  event[iLast].tau( event[iLast].tau0() * rndmPtr->exp() );

  // Mark original partons as hadronized and set their daughter range.
  for (int i = 0; i < int(iParton.size()); ++i) {
    event[ iParton[i] ].statusNeg();
    event[ iParton[i] ].daughters(iFirst, iLast);
  }    

  // Successfully done.
  return true;

}

//--------------------------------------------------------------------------

// Attempt to produce one particle from a ministring.
// Current algorithm: find the system with largest invariant mass
// relative to the existing one, and boost that system appropriately.
// Try more sophisticated alternatives later?? (Z0 mass shifted??)
// Also, if problems, attempt several times to obtain closer mass match??
  
bool MiniStringFragmentation::ministring2one( int iSub, 
  ColConfig& colConfig, Event& event) {

  // Cannot handle qq + qbarqbar system. 
  if (abs(flav1.id) > 100 && abs(flav2.id) > 100) return false;

  // For closed gluon loop need to pick an initial flavour.
  if (isClosed) do {
    int idStart = flavSelPtr->pickLightQ(); 
    FlavContainer flavStart(idStart, 1);
    flav1 = flavSelPtr->pick( flavStart);
    flav2 = flav1.anti();
  } while (abs(flav1.id) > 100);

  // Select hadron flavour from available quark flavours.
  int idHad = 0;
  for (int iTryFlav = 0; iTryFlav < NTRYFLAV; ++iTryFlav) {
    idHad = flavSelPtr->combine( flav1, flav2);
    if (idHad != 0) break;
  } 
  if (idHad == 0) return false;

  // Find mass.  
  double mHad = particleDataPtr->mass(idHad);
  
  // Find the untreated parton system which combines to the largest 
  // squared mass above mimimum required. 
  int iMax = -1;
  double deltaM2 = mHad*mHad - mSum*mSum; 
  double delta2Max = 0.;
  for (int iRec = iSub + 1; iRec < colConfig.size(); ++iRec) {
    double delta2Rec = 2. * (pSum * colConfig[iRec].pSum) - deltaM2 
      - 2. * mHad * colConfig[iRec].mass; 
    if (delta2Rec > delta2Max) { iMax = iRec; delta2Max = delta2Rec;}
  }
  if (iMax == -1) return false;  

  // Construct kinematics of the hadron and recoiling system. 
  Vec4& pRec     = colConfig[iMax].pSum;
  double mRec    = colConfig[iMax].mass;
  double vecProd = pSum * pRec; 
  double coefOld = mSum*mSum + vecProd;
  double coefNew = mHad*mHad + vecProd;
  double coefRec = mRec*mRec + vecProd;
  double coefSum = coefOld + coefNew;
  double sHat    = coefOld + coefRec;
  double root    = sqrtpos( (pow2(coefSum) - 4. * sHat * mHad*mHad)
    / (pow2(vecProd) - pow2(mSum * mRec)) );
  double k2      = 0.5 * (coefOld * root - coefSum) / sHat;
  double k1      = (coefRec * k2 + 0.5 * deltaM2) / coefOld;
  Vec4 pHad      = (1. + k1) * pSum - k2 * pRec;
  Vec4 pRecNew   = (1. + k2) * pRec - k1 * pSum;
  
  // Add the produced particle to the event record.
  int iHad = event.append( idHad, 81, iParton.front(), iParton.back(), 
    0, 0, 0, 0, pHad, mHad);

  // Set decay vertex when this is displaced.
  if (event[iParton.front()].hasVertex()) {
    Vec4 vDec = event[iParton.front()].vDec();
    event[iHad].vProd( vDec );
  }

  // Set lifetime of hadron.
  event[iHad].tau( event[iHad].tau0() * rndmPtr->exp() );

  // Mark original partons as hadronized and set their daughter range.
  for (int i = 0; i < int(iParton.size()); ++i) {
    event[ iParton[i] ].statusNeg();
    event[ iParton[i] ].daughters(iHad, iHad);
  }    
   
  // Copy down recoiling system, with boosted momentum. Update current partons.
  RotBstMatrix M;
  M.bst(pRec, pRecNew); 
  for (int i = 0; i < colConfig[iMax].size(); ++i) {
    int iOld = colConfig[iMax].iParton[i];
    // Do not touch negative iOld = beginning of new junction leg.
    if (iOld >= 0) {
      int iNew = event.copy(iOld, 72);
      event[iNew].rotbst(M);
      colConfig[iMax].iParton[i] = iNew;
    }
  }
  colConfig[iMax].pSum = pRecNew;
  colConfig[iMax].isCollected = true;

  // Successfully done.
  return true;

}

//==========================================================================

} // end namespace Pythia8
Commit	Line	Data
63ba5337	1	// MiniStringFragmentation.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 .
	7	// MiniStringFragmentation class
	8
	9	#include "MiniStringFragmentation.h"
	10
	11	namespace Pythia8 {
	12
	13	//==========================================================================
	14
	15	// The MiniStringFragmentation class.
	16
	17	//--------------------------------------------------------------------------
	18
	19	// Constants: could be changed here if desired, but normally should not.
	20	// These are of technical nature, as described for each.
	21
	22	// Since diffractive by definition is > 1 particle, try hard.
	23	const int MiniStringFragmentation::NTRYDIFFRACTIVE = 200;
	24
	25	// After one-body fragmentation failed, try two-body once more.
	26	const int MiniStringFragmentation::NTRYLASTRESORT = 100;
	27
	28	// Loop try to combine available endquarks to valid hadron.
	29	const int MiniStringFragmentation::NTRYFLAV = 10;
	30
	31	//--------------------------------------------------------------------------
	32
	33	// Initialize and save pointers.
	34
	35	void MiniStringFragmentation::init(Info* infoPtrIn, Settings& settings,
	36	ParticleData* particleDataPtrIn, Rndm* rndmPtrIn,
	37	StringFlav* flavSelPtrIn, StringPT* pTSelPtrIn, StringZ* zSelPtrIn) {
	38
	39	// Save pointers.
	40	infoPtr = infoPtrIn;
	41	particleDataPtr = particleDataPtrIn;
	42	rndmPtr = rndmPtrIn;
	43	flavSelPtr = flavSelPtrIn;
	44	pTSelPtr = pTSelPtrIn;
	45	zSelPtr = zSelPtrIn;
	46
	47	// Initialize the MiniStringFragmentation class proper.
	48	nTryMass = settings.mode("MiniStringFragmentation:nTry");
	49
	50	// Initialize the b parameter of the z spectrum, used when joining jets.
	51	bLund = zSelPtr->bAreaLund();
	52
	53	}
	54
	55	//--------------------------------------------------------------------------
	56
	57	// Do the fragmentation: driver routine.
	58
	59	bool MiniStringFragmentation::fragment(int iSub, ColConfig& colConfig,
	60	Event& event, bool isDiff) {
	61
	62	// Read in info on system to be treated.
	63	iParton = colConfig[iSub].iParton;
	64	flav1 = FlavContainer( event[ iParton.front() ].id() );
65	flav2 = FlavContainer( event[ iParton.back() ].id() );
66	pSum = colConfig[iSub].pSum;
67	mSum = colConfig[iSub].mass;
68	m2Sum = mSum*mSum;
69	isClosed = colConfig[iSub].isClosed;
70
71	// Do not want diffractive systems to easily collapse to one particle.
72	int nTryFirst = (isDiff) ? NTRYDIFFRACTIVE : nTryMass;
73
74	// First try to produce two particles from the system.
75	if (ministring2two( nTryFirst, event)) return true;
76
77	// If this fails, then form one hadron and shuffle momentum.
78	if (ministring2one( iSub, colConfig, event)) return true;
79
80	// If also this fails, then try harder to produce two particles.
81	if (ministring2two( NTRYLASTRESORT, event)) return true;
82
83	// Else complete failure.
84	infoPtr->errorMsg("Error in MiniStringFragmentation::fragment: "
85	"no 1- or 2-body state found above mass threshold");
86	return false;
87
88	}
89
90	//--------------------------------------------------------------------------
91
92	// Attempt to produce two particles from the ministring.
93
94	bool MiniStringFragmentation::ministring2two( int nTry, Event& event) {
95
96	// Properties of the produced hadrons.
97	int idHad1 = 0;
98	int idHad2 = 0;
99	double mHad1 = 0.;
100	double mHad2 = 0.;
101	double mHadSum = 0.;
102
103	// Allow a few attempts to find a particle pair with low enough masses.
104	for (int iTry = 0; iTry < nTry; ++iTry) {
105
106	// For closed gluon loop need to pick an initial flavour.
107	if (isClosed) do {
108	int idStart = flavSelPtr->pickLightQ();
109	FlavContainer flavStart(idStart, 1);
110	flavStart = flavSelPtr->pick( flavStart);
111	flav1 = flavSelPtr->pick( flavStart);
112	flav2.anti(flav1);
113	} while (flav1.id == 0 \|\| flav1.nPop > 0);
114
115	// Create a new q qbar flavour to form two hadrons.
116	// Start from a diquark, if any.
117	do {
118	FlavContainer flav3 =
119	(flav1.isDiquark() \|\| (!flav2.isDiquark() && rndmPtr->flat() < 0.5) )
120	? flavSelPtr->pick( flav1) : flavSelPtr->pick( flav2).anti();
121	idHad1 = flavSelPtr->combine( flav1, flav3);
122	idHad2 = flavSelPtr->combine( flav2, flav3.anti());
123	} while (idHad1 == 0 \|\| idHad2 == 0);
124
125	// Check whether the mass sum fits inside the available phase space.
126	mHad1 = particleDataPtr->mass(idHad1);
127	mHad2 = particleDataPtr->mass(idHad2);
128	mHadSum = mHad1 + mHad2;
129	if (mHadSum < mSum) break;
130	}
131	if (mHadSum >= mSum) return false;
132
133	// Define an effective two-parton string, by splitting intermediate
134	// gluon momenta in proportion to their closeness to either endpoint.
135	Vec4 pSum1 = event[ iParton.front() ].p();
136	Vec4 pSum2 = event[ iParton.back() ].p();
137	if (iParton.size() > 2) {
138	Vec4 pEnd1 = pSum1;
139	Vec4 pEnd2 = pSum2;
140	Vec4 pEndSum = pEnd1 + pEnd2;
141	for (int i = 1; i < int(iParton.size()) - 1 ; ++i) {
142	Vec4 pNow = event[ iParton[i] ].p();
143	double ratio = (pEnd2 * pNow) / (pEndSum * pNow);
144	pSum1 += ratio * pNow;
145	pSum2 += (1. - ratio) * pNow;
146	}
147	}
148
149	// Set up a string region based on the two effective endpoints.
150	StringRegion region;
151	region.setUp( pSum1, pSum2);
152
153	// Generate an isotropic decay in the ministring rest frame,
154	// suppressed at large pT by a fragmentation pT Gaussian.
155	double pAbs2 = 0.25 * ( pow2(m2Sum - mHad1mHad1 - mHad2mHad2)
156	- pow2(2. * mHad1 * mHad2) ) / m2Sum;
157	double pT2 = 0.;
158	do {
159	double cosTheta = rndmPtr->flat();
160	pT2 = (1. - pow2(cosTheta)) * pAbs2;
161	} while (pTSelPtr->suppressPT2(pT2) < rndmPtr->flat() );
162
163	// Construct the forward-backward asymmetry of the two particles.
164	double mT21 = mHad1*mHad1 + pT2;
165	double mT22 = mHad2*mHad2 + pT2;
166	double lambda = sqrtpos( pow2(m2Sum - mT21 - mT22) - 4. * mT21 * mT22 );
167	double probReverse = 1. / (1. + exp( min( 50., bLund * lambda) ) );
168
169	// Construct kinematics, as viewed in the transverse rest frame.
170	double xpz1 = 0.5 * lambda/ m2Sum;
171	if (probReverse > rndmPtr->flat()) xpz1 = -xpz1;
172	double xmDiff = (mT21 - mT22) / m2Sum;
173	double xe1 = 0.5 * (1. + xmDiff);
174	double xe2 = 0.5 * (1. - xmDiff );
175
176	// Distribute pT isotropically in angle.
177	double phi = 2. * M_PI * rndmPtr->flat();
178	double pT = sqrt(pT2);
179	double px = pT * cos(phi);
180	double py = pT * sin(phi);
181
182	// Translate this into kinematics in the string frame.
183	Vec4 pHad1 = region.pHad( xe1 + xpz1, xe1 - xpz1, px, py);
184	Vec4 pHad2 = region.pHad( xe2 - xpz1, xe2 + xpz1, -px, -py);
185
186	// Add produced particles to the event record.
187	int iFirst = event.append( idHad1, 82, iParton.front(), iParton.back(),
188	0, 0, 0, 0, pHad1, mHad1);
189	int iLast = event.append( idHad2, 82, iParton.front(), iParton.back(),
190	0, 0, 0, 0, pHad2, mHad2);
191
192	// Set decay vertex when this is displaced.
193	if (event[iParton.front()].hasVertex()) {
194	Vec4 vDec = event[iParton.front()].vDec();
195	event[iFirst].vProd( vDec );
196	event[iLast].vProd( vDec );
197	}
198
199	// Set lifetime of hadrons.
200	event[iFirst].tau( event[iFirst].tau0() * rndmPtr->exp() );
201	event[iLast].tau( event[iLast].tau0() * rndmPtr->exp() );
202
203	// Mark original partons as hadronized and set their daughter range.
204	for (int i = 0; i < int(iParton.size()); ++i) {
205	event[ iParton[i] ].statusNeg();
206	event[ iParton[i] ].daughters(iFirst, iLast);
207	}
208
209	// Successfully done.
210	return true;
211
212	}
213
214	//--------------------------------------------------------------------------
215
216	// Attempt to produce one particle from a ministring.
217	// Current algorithm: find the system with largest invariant mass
218	// relative to the existing one, and boost that system appropriately.
219	// Try more sophisticated alternatives later?? (Z0 mass shifted??)
220	// Also, if problems, attempt several times to obtain closer mass match??
221
222	bool MiniStringFragmentation::ministring2one( int iSub,
223	ColConfig& colConfig, Event& event) {
224
225	// Cannot handle qq + qbarqbar system.
226	if (abs(flav1.id) > 100 && abs(flav2.id) > 100) return false;
227
228	// For closed gluon loop need to pick an initial flavour.
229	if (isClosed) do {
230	int idStart = flavSelPtr->pickLightQ();
231	FlavContainer flavStart(idStart, 1);
232	flav1 = flavSelPtr->pick( flavStart);
233	flav2 = flav1.anti();
234	} while (abs(flav1.id) > 100);
235
236	// Select hadron flavour from available quark flavours.
237	int idHad = 0;
238	for (int iTryFlav = 0; iTryFlav < NTRYFLAV; ++iTryFlav) {
239	idHad = flavSelPtr->combine( flav1, flav2);
240	if (idHad != 0) break;
241	}
242	if (idHad == 0) return false;
243
244	// Find mass.
245	double mHad = particleDataPtr->mass(idHad);
246
247	// Find the untreated parton system which combines to the largest
248	// squared mass above mimimum required.
249	int iMax = -1;
250	double deltaM2 = mHadmHad - mSummSum;
251	double delta2Max = 0.;
252	for (int iRec = iSub + 1; iRec < colConfig.size(); ++iRec) {
253	double delta2Rec = 2. * (pSum * colConfig[iRec].pSum) - deltaM2
254	- 2. * mHad * colConfig[iRec].mass;
255	if (delta2Rec > delta2Max) { iMax = iRec; delta2Max = delta2Rec;}
256	}
257	if (iMax == -1) return false;
258
259	// Construct kinematics of the hadron and recoiling system.
260	Vec4& pRec = colConfig[iMax].pSum;
261	double mRec = colConfig[iMax].mass;
262	double vecProd = pSum * pRec;
263	double coefOld = mSum*mSum + vecProd;
264	double coefNew = mHad*mHad + vecProd;
265	double coefRec = mRec*mRec + vecProd;
266	double coefSum = coefOld + coefNew;
267	double sHat = coefOld + coefRec;
268	double root = sqrtpos( (pow2(coefSum) - 4. * sHat * mHad*mHad)
269	/ (pow2(vecProd) - pow2(mSum * mRec)) );
270	double k2 = 0.5 * (coefOld * root - coefSum) / sHat;
271	double k1 = (coefRec * k2 + 0.5 * deltaM2) / coefOld;
272	Vec4 pHad = (1. + k1) * pSum - k2 * pRec;
273	Vec4 pRecNew = (1. + k2) * pRec - k1 * pSum;
274
275	// Add the produced particle to the event record.
276	int iHad = event.append( idHad, 81, iParton.front(), iParton.back(),
277	0, 0, 0, 0, pHad, mHad);
278
279	// Set decay vertex when this is displaced.
280	if (event[iParton.front()].hasVertex()) {
281	Vec4 vDec = event[iParton.front()].vDec();
282	event[iHad].vProd( vDec );
283	}
284
285	// Set lifetime of hadron.
286	event[iHad].tau( event[iHad].tau0() * rndmPtr->exp() );
287
288	// Mark original partons as hadronized and set their daughter range.
289	for (int i = 0; i < int(iParton.size()); ++i) {
290	event[ iParton[i] ].statusNeg();
291	event[ iParton[i] ].daughters(iHad, iHad);
292	}
293
294	// Copy down recoiling system, with boosted momentum. Update current partons.
295	RotBstMatrix M;
296	M.bst(pRec, pRecNew);
297	for (int i = 0; i < colConfig[iMax].size(); ++i) {
298	int iOld = colConfig[iMax].iParton[i];
299	// Do not touch negative iOld = beginning of new junction leg.
300	if (iOld >= 0) {
301	int iNew = event.copy(iOld, 72);
302	event[iNew].rotbst(M);
303	colConfig[iMax].iParton[i] = iNew;
304	}
305	}
306	colConfig[iMax].pSum = pRecNew;
307	colConfig[iMax].isCollected = true;
308
309	// Successfully done.
310	return true;
311
312	}
313
314	//==========================================================================
315
316	} // end namespace Pythia8