ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RPGFragmentation.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RPGFragmentation.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 
28 #include <iostream>
29 #include <signal.h>
30 
31 #include "G4Log.hh"
32 #include "G4Pow.hh"
33 #include "G4RPGFragmentation.hh"
34 #include "G4PhysicalConstants.hh"
35 #include "G4SystemOfUnits.hh"
36 #include "G4AntiProton.hh"
37 #include "G4AntiNeutron.hh"
38 #include "Randomize.hh"
39 #include "G4Poisson.hh"
41 
43  : G4RPGReaction()
44 {
45  for (G4int i = 0; i < 20; i++) dndl[i] = 0.0;
46 }
47 
48 
51 {
52  pt = std::max( 0.001, pt );
53  G4double dx = 1./(19.*pt);
54  G4double x;
55  G4double term1;
56  G4double term2;
57 
58  for (G4int i = 1; i < 20; i++) {
59  x = (G4double(i) - 0.5)*dx;
60  term1 = 1. + parMass*parMass*x*x;
61  term2 = pt*x*et*pt*x*et + pt*pt + secMass*secMass;
62  dndl[i] = dx / std::sqrt( term1*term1*term1*term2 )
63  + dndl[i-1];
64  }
65 }
66 
67 
69 ReactionStage(const G4HadProjectile* originalIncident,
70  G4ReactionProduct& modifiedOriginal,
71  G4bool& incidentHasChanged,
72  const G4DynamicParticle* originalTarget,
73  G4ReactionProduct& targetParticle,
74  G4bool& targetHasChanged,
75  const G4Nucleus& targetNucleus,
76  G4ReactionProduct& currentParticle,
78  G4int& vecLen,
79  G4bool leadFlag,
80  G4ReactionProduct& leadingStrangeParticle)
81 {
82  //
83  // Based on H. Fesefeldt's original FORTRAN code GENXPT
84  //
85  // Generation of x- and pT- values for incident, target, and all secondary
86  // particles using a simple single variable description E D3S/DP3= F(Q)
87  // with Q^2 = (M*X)^2 + PT^2. Final state kinematics are produced by an
88  // FF-type iterative cascade method
89  //
90  // Internal units are GeV
91  //
92 
93  // Protection in case no secondary has been created. In that case use
94  // two-body scattering
95  //
96  if (vecLen == 0) return false;
97 
103 
104  G4int i, l;
105  G4bool veryForward = false;
106 
107  const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/GeV;
108  const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
109  const G4double mOriginal = modifiedOriginal.GetMass()/GeV;
110  const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/GeV;
111  G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;
112  G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
113  targetMass*targetMass +
114  2.0*targetMass*etOriginal ); // GeV
115  G4double currentMass = currentParticle.GetMass()/GeV;
116  targetMass = targetParticle.GetMass()/GeV;
117 
118  // Randomize the order of the secondary particles.
119  // Note that the current and target particles are not affected.
120 
121  for (i=0; i<vecLen; ++i) {
122  G4int itemp = G4int( G4UniformRand()*vecLen );
123  G4ReactionProduct pTemp = *vec[itemp];
124  *vec[itemp] = *vec[i];
125  *vec[i] = pTemp;
126  }
127 
128  if (currentMass == 0.0 && targetMass == 0.0) {
129  // Target and projectile have annihilated. Replace them with the first
130  // two secondaries in the list. Current particle KE is maintained.
131 
132  G4double ek = currentParticle.GetKineticEnergy();
133  G4ThreeVector mom = currentParticle.GetMomentum();
134  currentParticle = *vec[0];
135  currentParticle.SetSide(1);
136  targetParticle = *vec[1];
137  targetParticle.SetSide(-1);
138  for( i=0; i<(vecLen-2); ++i )*vec[i] = *vec[i+2];
139  G4ReactionProduct *temp = vec[vecLen-1];
140  delete temp;
141  temp = vec[vecLen-2];
142  delete temp;
143  vecLen -= 2;
144  currentMass = currentParticle.GetMass()/GeV;
145  targetMass = targetParticle.GetMass()/GeV;
146  incidentHasChanged = true;
147  targetHasChanged = true;
148  currentParticle.SetKineticEnergy( ek );
149  currentParticle.SetMomentum(mom);
150  veryForward = true;
151  }
152  const G4double atomicWeight = targetNucleus.GetA_asInt();
153  const G4double atomicNumber = targetNucleus.GetZ_asInt();
154  const G4double protonMass = aProton->GetPDGMass();
155 
156  if (originalIncident->GetDefinition()->GetParticleSubType() == "kaon"
157  && G4UniformRand() >= 0.7) {
158  G4ReactionProduct temp = currentParticle;
159  currentParticle = targetParticle;
160  targetParticle = temp;
161  incidentHasChanged = true;
162  targetHasChanged = true;
163  currentMass = currentParticle.GetMass()/GeV;
164  targetMass = targetParticle.GetMass()/GeV;
165  }
166  const G4double afc = std::min( 0.75,
167  0.312+0.200*G4Log(G4Log(centerofmassEnergy*centerofmassEnergy))+
168  std::pow(centerofmassEnergy*centerofmassEnergy,1.5)/6000.0 );
169 
170  G4double freeEnergy = centerofmassEnergy-currentMass-targetMass;
171  G4double forwardEnergy = freeEnergy/2.;
172  G4int forwardCount = 1; // number of particles in forward hemisphere
173 
174  G4double backwardEnergy = freeEnergy/2.;
175  G4int backwardCount = 1; // number of particles in backward hemisphere
176 
177  if(veryForward) {
178  if(currentParticle.GetSide()==-1)
179  {
180  forwardEnergy += currentMass;
181  forwardCount --;
182  backwardEnergy -= currentMass;
183  backwardCount ++;
184  }
185  if(targetParticle.GetSide()!=-1)
186  {
187  backwardEnergy += targetMass;
188  backwardCount --;
189  forwardEnergy -= targetMass;
190  forwardCount ++;
191  }
192  }
193 
194  for (i=0; i<vecLen; ++i) {
195  if( vec[i]->GetSide() == -1 )
196  {
197  ++backwardCount;
198  backwardEnergy -= vec[i]->GetMass()/GeV;
199  } else {
200  ++forwardCount;
201  forwardEnergy -= vec[i]->GetMass()/GeV;
202  }
203  }
204 
205  // Check that sum of forward particle masses does not exceed forwardEnergy,
206  // and similarly for backward particles. If so, move particles from one
207  // hemisphere to another.
208 
209  if (backwardEnergy < 0.0) {
210  for (i = 0; i < vecLen; ++i) {
211  if (vec[i]->GetSide() == -1) {
212  backwardEnergy += vec[i]->GetMass()/GeV;
213  --backwardCount;
214  vec[i]->SetSide(1);
215  forwardEnergy -= vec[i]->GetMass()/GeV;
216  ++forwardCount;
217  if (backwardEnergy > 0.0) break;
218  }
219  }
220  }
221 
222  if (forwardEnergy < 0.0) {
223  for (i = 0; i < vecLen; ++i) {
224  if (vec[i]->GetSide() == 1) {
225  forwardEnergy += vec[i]->GetMass()/GeV;
226  --forwardCount;
227  vec[i]->SetSide(-1);
228  backwardEnergy -= vec[i]->GetMass()/GeV;
229  ++backwardCount;
230  if (forwardEnergy > 0.0) break;
231  }
232  }
233  }
234 
235  // Special cases for reactions near threshold
236 
237  // 1. There is only one secondary
238  if (forwardEnergy > 0.0 && backwardEnergy < 0.0) {
239  forwardEnergy += backwardEnergy;
240  backwardEnergy = 0;
241  }
242 
243  // 2. Nuclear binding energy is large
244  if (forwardEnergy + backwardEnergy < 0.0) return false;
245 
246 
247  // forwardEnergy now represents the total energy in the forward reaction
248  // hemisphere which is available for kinetic energy and particle creation.
249  // Similarly for backwardEnergy.
250 
251  // Add particles from the intranuclear cascade.
252  // nuclearExcitationCount = number of new secondaries produced by nuclear
253  // excitation
254  // extraCount = number of nucleons within these new secondaries
255  //
256  // Note: eventually have to make sure that enough nucleons are available
257  // in the case of small target nuclei
258 
259  G4double xtarg;
260  G4double a13 = G4Pow::GetInstance()->A13(atomicWeight); // A**(1/3)
261  if (centerofmassEnergy < (2.0+G4UniformRand()) )
262  xtarg = afc * (a13-1.0) * (2.0*backwardCount+vecLen+2)/2.0;
263  else
264  xtarg = afc * (a13-1.0) * (2.0*backwardCount);
265 
266  if( xtarg <= 0.0 )xtarg = 0.01;
267  G4int nuclearExcitationCount = G4Poisson( xtarg );
268  // To do: try reducing nuclearExcitationCount with increasing energy
269  // to simulate cut-off of cascade
270  if(atomicWeight<1.0001) nuclearExcitationCount = 0;
271  G4int extraNucleonCount = 0;
272  G4double extraNucleonMass = 0.0;
273 
274  if (nuclearExcitationCount > 0) {
275  const G4double nucsup[] = { 1.00, 0.7, 0.5, 0.4, 0.35, 0.3 };
276  const G4double psup[] = { 3., 6., 20., 50., 100., 1000. };
277  G4int momentumBin = 0;
278 
279  G4int loop = 0;
281  ed << " While count exceeded " << G4endl;
282  while( (momentumBin < 6) &&
283  (modifiedOriginal.GetTotalMomentum()/GeV > psup[momentumBin]) ) { /* Loop checking, 01.09.2015, D.Wright */
284  ++momentumBin;
285  loop++;
286  if (loop > 1000) {
287  G4Exception("G4RPGFragmentation::ReactionStage()", "HAD_RPG_100", JustWarning, ed);
288  break;
289  }
290  }
291 
292  momentumBin = std::min( 5, momentumBin );
293 
294  // NOTE: in GENXPT, these new particles were given negative codes
295  // here I use NewlyAdded = true instead
296  //
297 
298  for (i = 0; i < nuclearExcitationCount; ++i) {
299  G4ReactionProduct * pVec = new G4ReactionProduct();
300  if (G4UniformRand() < nucsup[momentumBin]) {
301 
302  if (G4UniformRand() > 1.0-atomicNumber/atomicWeight)
303  pVec->SetDefinition( aProton );
304  else
305  pVec->SetDefinition( aNeutron );
306 
307  // nucleon comes from nucleus -
308  // do not subtract its mass from backward energy
309  pVec->SetSide( -2 ); // -2 means backside nucleon
310  ++extraNucleonCount;
311  extraNucleonMass += pVec->GetMass()/GeV;
312  // To do: Remove chosen nucleon from target nucleus
313  pVec->SetNewlyAdded( true );
314  vec.SetElement( vecLen++, pVec );
315  ++backwardCount;
316 
317  } else {
318 
320  if( ran < 0.3181 )
321  pVec->SetDefinition( aPiPlus );
322  else if( ran < 0.6819 )
323  pVec->SetDefinition( aPiZero );
324  else
325  pVec->SetDefinition( aPiMinus );
326 
327  if (backwardEnergy > pVec->GetMass()/GeV) {
328  backwardEnergy -= pVec->GetMass()/GeV; // pion mass comes from free energy
329  ++backwardCount;
330  pVec->SetSide( -1 ); // backside particle, but not a nucleon
331  pVec->SetNewlyAdded( true );
332  vec.SetElement( vecLen++, pVec );
333  }
334 
335  // To do: Change proton to neutron (or vice versa) in target nucleus depending
336  // on pion charge
337  }
338  }
339  }
340 
341  // Define initial state vectors for Lorentz transformations
342  // The pseudoParticles have non-standard masses, hence the "pseudo"
343 
344  G4ReactionProduct pseudoParticle[8];
345  for (i = 0; i < 8; ++i) pseudoParticle[i].SetZero();
346 
347  pseudoParticle[0].SetMass( mOriginal*GeV );
348  pseudoParticle[0].SetMomentum( 0.0, 0.0, pOriginal*GeV );
349  pseudoParticle[0].SetTotalEnergy(
350  std::sqrt( pOriginal*pOriginal + mOriginal*mOriginal )*GeV );
351 
352  pseudoParticle[1].SetMass(protonMass); // this could be targetMass
353  pseudoParticle[1].SetTotalEnergy(protonMass);
354 
355  pseudoParticle[3].SetMass(protonMass*(1+extraNucleonCount) );
356  pseudoParticle[3].SetTotalEnergy(protonMass*(1+extraNucleonCount) );
357 
358  pseudoParticle[2] = pseudoParticle[0] + pseudoParticle[1];
359  pseudoParticle[3] = pseudoParticle[3] + pseudoParticle[0];
360 
361  pseudoParticle[0].Lorentz( pseudoParticle[0], pseudoParticle[2] );
362  pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[2] );
363 
364  // Main loop for 4-momentum generation
365  // See Pitha-report (Aachen) for a detailed description of the method
366 
367  G4double aspar, pt, et, x, pp, pp1, wgt;
368  G4int innerCounter, outerCounter;
369  G4bool eliminateThisParticle, resetEnergies, constantCrossSection;
370 
371  G4double forwardKinetic = 0.0;
372  G4double backwardKinetic = 0.0;
373 
374  // Process the secondary particles in reverse order
375  // The incident particle is done after the secondaries
376  // Nucleons, including the target, in the backward hemisphere are also
377  // done later
378 
379  G4int backwardNucleonCount = 0; // number of nucleons in backward hemisphere
380  G4double totalEnergy, kineticEnergy, vecMass;
381  G4double phi;
382 
383  for (i = vecLen-1; i >= 0; --i) {
384 
385  if (vec[i]->GetNewlyAdded()) { // added from intranuclear cascade
386  if (vec[i]->GetSide() == -2) { // its a nucleon
387  if (backwardNucleonCount < 18) {
388  if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
389  for (G4int j = 0; j < vecLen; j++) delete vec[j];
390  vecLen = 0;
391  throw G4HadReentrentException(__FILE__, __LINE__,
392  "G4RPGFragmentation::ReactionStage : a pion has been counted as a backward nucleon");
393  }
394  vec[i]->SetSide(-3);
395  ++backwardNucleonCount;
396  continue; // Don't generate momenta for the first 17 backward
397  // cascade nucleons. This gets done by the cluster
398  // model later on.
399  }
400  }
401  }
402 
403  // Set pt and phi values, they are changed somewhat in the iteration loop
404  // Set mass parameter for lambda fragmentation model
405 
406  vecMass = vec[i]->GetMass()/GeV;
407  G4double ran = -G4Log(1.0-G4UniformRand())/3.5;
408 
409  if (vec[i]->GetSide() == -2) { // backward nucleon
410  aspar = 0.20;
411  pt = std::sqrt( std::pow( ran, 1.2 ) );
412 
413  } else { // not a backward nucleon
414  if (vec[i]->GetDefinition()->GetParticleSubType() == "pi") {
415  aspar = 0.75;
416  pt = std::sqrt( std::pow( ran, 1.7 ) );
417  } else if (vec[i]->GetDefinition()->GetParticleSubType() == "kaon") {
418  aspar = 0.70;
419  pt = std::sqrt( std::pow( ran, 1.7 ) );
420  } else { // vec[i] must be a baryon or ion
421  aspar = 0.65;
422  pt = std::sqrt( std::pow( ran, 1.5 ) );
423  }
424  }
425 
426  pt = std::max( 0.001, pt );
427  phi = G4UniformRand()*twopi;
428  vec[i]->SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
429  if (vec[i]->GetSide() > 0)
430  et = pseudoParticle[0].GetTotalEnergy()/GeV;
431  else
432  et = pseudoParticle[1].GetTotalEnergy()/GeV;
433 
434  //
435  // Start of outer iteration loop
436  //
437  outerCounter = 0;
438  eliminateThisParticle = true;
439  resetEnergies = true;
440  dndl[0] = 0.0;
441 
442  while (++outerCounter < 3) { /* Loop checking, 01.09.2015, D.Wright */
443 
444  FragmentationIntegral(pt, et, aspar, vecMass);
445  innerCounter = 0;
446  vec[i]->SetMomentum( pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV );
447 
448  // Start of inner iteration loop
449 
450  while (++innerCounter < 7) { /* Loop checking, 01.09.2015, D.Wright */
451 
452  ran = G4UniformRand()*dndl[19];
453  l = 1;
454 
455  G4int loop = 0;
457  ed << " While count exceeded " << G4endl;
458  while( ( ran > dndl[l] ) && ( l < 19 ) ) { /* Loop checking, 01.09.2015, D.Wright */
459  l++;
460  loop++;
461  if (loop > 1000) {
462  G4Exception("G4RPGFragmentation::ReactionStage()", "HAD_RPG_100", JustWarning, ed);
463  break;
464  }
465  }
466 
467  x = (G4double(l-1) + G4UniformRand())/19.;
468  if (vec[i]->GetSide() < 0) x *= -1.;
469  vec[i]->SetMomentum( x*et*GeV ); // set the z-momentum
470  totalEnergy = std::sqrt( x*et*x*et + pt*pt + vecMass*vecMass );
471  vec[i]->SetTotalEnergy( totalEnergy*GeV );
472  kineticEnergy = vec[i]->GetKineticEnergy()/GeV;
473 
474  if (vec[i]->GetSide() > 0) { // forward side
475  if( (forwardKinetic+kineticEnergy) < 0.95*forwardEnergy ) {
476  // Leave at least 5% of the forward free energy for the projectile
477 
478  pseudoParticle[4] = pseudoParticle[4] + (*vec[i]);
479  forwardKinetic += kineticEnergy;
480  outerCounter = 2; // leave outer loop
481  eliminateThisParticle = false; // don't eliminate this particle
482  resetEnergies = false;
483  break; // leave inner loop
484  }
485  if( innerCounter > 5 )break; // leave inner loop
486  if( backwardEnergy >= vecMass ) // switch sides
487  {
488  vec[i]->SetSide(-1);
489  forwardEnergy += vecMass;
490  backwardEnergy -= vecMass;
491  ++backwardCount;
492  }
493  } else { // backward side
494  // if (extraNucleonCount > 19) x = 0.999;
495  // G4double xxx = 0.95+0.05*extraNucleonCount/20.0;
496  // DHW: I think above lines were meant to be as follows:
497  G4double xxx = 0.999;
498  if (extraNucleonCount < 20) xxx = 0.95+0.05*extraNucleonCount/20.0;
499 
500  if ((backwardKinetic+kineticEnergy) < xxx*backwardEnergy) {
501  pseudoParticle[5] = pseudoParticle[5] + (*vec[i]);
502  backwardKinetic += kineticEnergy;
503  outerCounter = 2; // leave outer loop
504  eliminateThisParticle = false; // don't eliminate this particle
505  resetEnergies = false;
506  break; // leave inner loop
507  }
508  if (innerCounter > 5) break; // leave inner loop
509  if (forwardEnergy >= vecMass) { // switch sides
510  vec[i]->SetSide(1);
511  forwardEnergy -= vecMass;
512  backwardEnergy += vecMass;
513  backwardCount--;
514  }
515  }
516  G4ThreeVector momentum = vec[i]->GetMomentum();
517  vec[i]->SetMomentum( momentum.x() * 0.9, momentum.y() * 0.9 );
518  pt *= 0.9;
519  dndl[19] *= 0.9;
520  } // closes inner loop
521 
522  // If we get here, the inner loop has been done 6 times.
523  // If necessary, reduce energies of the previously done particles if
524  // they are lighter than protons or are in the forward hemisphere.
525  // Then continue with outer loop.
526 
527  if (resetEnergies)
528  ReduceEnergiesOfSecondaries(i+1, forwardKinetic, backwardKinetic,
529  vec, vecLen,
530  pseudoParticle[4], pseudoParticle[5],
531  pt);
532 
533  } // closes outer loop
534 
535  if (eliminateThisParticle && vec[i]->GetMayBeKilled()) {
536  // not enough energy, eliminate this particle
537 
538  if (vec[i]->GetSide() > 0) {
539  --forwardCount;
540  forwardEnergy += vecMass;
541  } else {
542  --backwardCount;
543  if (vec[i]->GetSide() == -2) {
544  --extraNucleonCount;
545  extraNucleonMass -= vecMass;
546  } else {
547  backwardEnergy += vecMass;
548  }
549  }
550 
551  for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
552  G4ReactionProduct* temp = vec[vecLen-1];
553  delete temp;
554  // To do: modify target nucleus according to particle eliminated
555 
556  if( --vecLen == 0 ){
557  G4cout << " FALSE RETURN DUE TO ENERGY BALANCE " << G4endl;
558  return false;
559  } // all the secondaries have been eliminated
560  }
561  } // closes main loop over secondaries
562 
563  // Re-balance forward and backward energy if possible and if necessary
564 
565  G4double forwardKEDiff = forwardEnergy - forwardKinetic;
566  G4double backwardKEDiff = backwardEnergy - backwardKinetic;
567 
568  if (forwardKEDiff < 0.0 || backwardKEDiff < 0.0) {
569  ReduceEnergiesOfSecondaries(0, forwardKinetic, backwardKinetic,
570  vec, vecLen,
571  pseudoParticle[4], pseudoParticle[5],
572  pt);
573 
574  forwardKEDiff = forwardEnergy - forwardKinetic;
575  backwardKEDiff = backwardEnergy - backwardKinetic;
576  if (backwardKEDiff < 0.0) {
577  if (forwardKEDiff + backwardKEDiff > 0.0) {
578  backwardEnergy = backwardKinetic;
579  forwardEnergy += backwardKEDiff;
580  forwardKEDiff = forwardEnergy - forwardKinetic;
581  backwardKEDiff = 0.0;
582  } else {
583  G4cout << " False return due to insufficient backward energy " << G4endl;
584  return false;
585  }
586  }
587 
588  if (forwardKEDiff < 0.0) {
589  if (forwardKEDiff + backwardKEDiff > 0.0) {
590  forwardEnergy = forwardKinetic;
591  backwardEnergy += forwardKEDiff;
592  backwardKEDiff = backwardEnergy - backwardKinetic;
593  forwardKEDiff = 0.0;
594  } else {
595  G4cout << " False return due to insufficient forward energy " << G4endl;
596  return false;
597  }
598  }
599  }
600 
601  // Generate momentum for incident (current) particle, which was placed
602  // in the forward hemisphere.
603  // Set mass parameter for lambda fragmentation model.
604  // Set pt and phi values, which are changed somewhat in the iteration loop
605 
606  G4double ran = -G4Log(1.0-G4UniformRand());
607  if (currentParticle.GetDefinition()->GetParticleSubType() == "pi") {
608  aspar = 0.60;
609  pt = std::sqrt( std::pow( ran/6.0, 1.7 ) );
610  } else if (currentParticle.GetDefinition()->GetParticleSubType() == "kaon") {
611  aspar = 0.50;
612  pt = std::sqrt( std::pow( ran/5.0, 1.4 ) );
613  } else {
614  aspar = 0.40;
615  pt = std::sqrt( std::pow( ran/4.0, 1.2 ) );
616  }
617 
618  phi = G4UniformRand()*twopi;
619  currentParticle.SetMomentum(pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV);
620  et = pseudoParticle[0].GetTotalEnergy()/GeV;
621  dndl[0] = 0.0;
622  vecMass = currentParticle.GetMass()/GeV;
623 
624  FragmentationIntegral(pt, et, aspar, vecMass);
625 
626  ran = G4UniformRand()*dndl[19];
627  l = 1;
628 
629  G4int loop = 0;
631  ed << " While count exceeded " << G4endl;
632  while( ( ran > dndl[l] ) && ( l < 19 ) ) { /* Loop checking, 01.09.2015, D.Wright */
633  l++;
634  loop++;
635  if (loop > 1000) {
636  G4Exception("G4RPGFragmentation::ReactionStage()", "HAD_RPG_100", JustWarning, ed);
637  break;
638  }
639  }
640 
641  x = (G4double(l-1) + G4UniformRand())/19.;
642  currentParticle.SetMomentum( x*et*GeV ); // set the z-momentum
643 
644  if (forwardEnergy < forwardKinetic) {
645  totalEnergy = vecMass + 0.04*std::fabs(normal());
646  G4cout << " Not enough forward energy: forwardEnergy = "
647  << forwardEnergy << " forwardKinetic = "
648  << forwardKinetic << " total energy left = "
649  << backwardKEDiff + forwardKEDiff << G4endl;
650  } else {
651  totalEnergy = vecMass + forwardEnergy - forwardKinetic;
652  forwardKinetic = forwardEnergy;
653  }
654  currentParticle.SetTotalEnergy( totalEnergy*GeV );
655  pp = std::sqrt(std::abs( totalEnergy*totalEnergy - vecMass*vecMass) )*GeV;
656  pp1 = currentParticle.GetMomentum().mag();
657 
658  if (pp1 < 1.0e-6*GeV) {
659  G4ThreeVector iso = Isotropic(pp);
660  currentParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
661  } else {
662  currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
663  }
664  pseudoParticle[4] = pseudoParticle[4] + currentParticle;
665 
666  // Current particle now finished
667 
668  // Begin target particle
669 
670  if (backwardNucleonCount < 18) {
671  targetParticle.SetSide(-3);
672  ++backwardNucleonCount;
673 
674  } else {
675  // Set pt and phi values, they are changed somewhat in the iteration loop
676  // Set mass parameter for lambda fragmentation model
677 
678  vecMass = targetParticle.GetMass()/GeV;
679  ran = -G4Log(1.0-G4UniformRand());
680  aspar = 0.40;
681  pt = std::max( 0.001, std::sqrt( std::pow( ran/4.0, 1.2 ) ) );
682  phi = G4UniformRand()*twopi;
683  targetParticle.SetMomentum(pt*std::cos(phi)*GeV, pt*std::sin(phi)*GeV);
684  et = pseudoParticle[1].GetTotalEnergy()/GeV;
685  outerCounter = 0;
686  innerCounter = 0;
687  G4bool marginalEnergy = true;
688  dndl[0] = 0.0;
689  G4double xxx = 0.999;
690  if( extraNucleonCount < 20 ) xxx = 0.95+0.05*extraNucleonCount/20.0;
692 
693  while (++outerCounter < 4) { /* Loop checking, 01.09.2015, D.Wright */
694  FragmentationIntegral(pt, et, aspar, vecMass);
695 
696  for (innerCounter = 0; innerCounter < 6; innerCounter++) {
697  ran = G4UniformRand()*dndl[19];
698  l = 1;
699 
700  G4int loopa = 0;
702  eda << " While count exceeded " << G4endl;
703  while( ( ran > dndl[l] ) && ( l < 19 ) ) { /* Loop checking, 01.09.2015, D.Wright */
704  l++;
705  loopa++;
706  if (loopa > 1000) {
707  G4Exception("G4RPGFragmentation::ReactionStage()", "HAD_RPG_100", JustWarning, eda);
708  break;
709  }
710  }
711 
712  x = -(G4double(l-1) + G4UniformRand())/19.;
713  targetParticle.SetMomentum( x*et*GeV ); // set the z-momentum
714  totalEnergy = std::sqrt(x*et*x*et + pt*pt + vecMass*vecMass);
715  targetParticle.SetTotalEnergy( totalEnergy*GeV );
716 
717  if ((backwardKinetic+totalEnergy-vecMass) < xxx*backwardEnergy) {
718  pseudoParticle[5] = pseudoParticle[5] + targetParticle;
719  backwardKinetic += totalEnergy - vecMass;
720  outerCounter = 3; // leave outer loop
721  marginalEnergy = false;
722  break; // leave inner loop
723  }
724  momentum = targetParticle.GetMomentum();
725  targetParticle.SetMomentum(momentum.x() * 0.9, momentum.y() * 0.9);
726  pt *= 0.9;
727  dndl[19] *= 0.9;
728  }
729  }
730 
731  if (marginalEnergy) {
732  G4cout << " Extra backward kinetic energy = "
733  << 0.999*backwardEnergy - backwardKinetic << G4endl;
734  totalEnergy = vecMass + 0.999*backwardEnergy - backwardKinetic;
735  targetParticle.SetTotalEnergy(totalEnergy*GeV);
736  pp = std::sqrt(std::abs(totalEnergy*totalEnergy - vecMass*vecMass) )*GeV;
737  targetParticle.SetMomentum(momentum.x()/0.9, momentum.y()/0.9);
738  pp1 = targetParticle.GetMomentum().mag();
739  targetParticle.SetMomentum(targetParticle.GetMomentum() * pp/pp1 );
740  pseudoParticle[5] = pseudoParticle[5] + targetParticle;
741  backwardKinetic = 0.999*backwardEnergy;
742  }
743 
744  }
745 
746  if (backwardEnergy < backwardKinetic)
747  G4cout << " Backward Edif = " << backwardEnergy - backwardKinetic << G4endl;
748  if (forwardEnergy != forwardKinetic)
749  G4cout << " Forward Edif = " << forwardEnergy - forwardKinetic << G4endl;
750 
751  // Target particle finished.
752  // Now produce backward nucleons with a cluster model
753  // ps[2] = CM frame of projectile + target
754  // ps[3] = sum of projectile + nucleon cluster in lab frame
755  // ps[6] = proj + cluster 4-vector boosted into CM frame of proj + targ
756  // with secondaries, current and target particles subtracted
757  // = total 4-momentum of backward nucleon cluster
758 
759  pseudoParticle[6].Lorentz( pseudoParticle[3], pseudoParticle[2] );
760  pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[4];
761  pseudoParticle[6] = pseudoParticle[6] - pseudoParticle[5];
762 
763  if (backwardNucleonCount == 1) {
764  // Target particle is the only backward nucleon. Give it the remainder
765  // of the backward kinetic energy.
766 
767  G4double ekin =
768  std::min(backwardEnergy-backwardKinetic, centerofmassEnergy/2.0-protonMass/GeV);
769 
770  if( ekin < 0.04 )ekin = 0.04 * std::fabs( normal() );
771  vecMass = targetParticle.GetMass()/GeV;
772  totalEnergy = ekin + vecMass;
773  targetParticle.SetTotalEnergy( totalEnergy*GeV );
774  pp = std::sqrt(std::abs(totalEnergy*totalEnergy - vecMass*vecMass) )*GeV;
775  pp1 = pseudoParticle[6].GetMomentum().mag();
776  if (pp1 < 1.0e-6*GeV) {
777  G4ThreeVector iso = Isotropic(pp);
778  targetParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
779  } else {
780  targetParticle.SetMomentum( pseudoParticle[6].GetMomentum() * (pp/pp1));
781  }
782  pseudoParticle[5] = pseudoParticle[5] + targetParticle;
783 
784  } else if (backwardNucleonCount > 1) {
785  // Share remaining energy with up to 17 backward nucleons
786 
787  G4int tempCount = 5;
788  if (backwardNucleonCount < 5) tempCount = backwardNucleonCount;
789  tempCount -= 2;
790 
791  G4double clusterMass = 0.;
792  if (targetParticle.GetSide() == -3)
793  clusterMass = targetParticle.GetMass()/GeV;
794  for (i = 0; i < vecLen; ++i)
795  if (vec[i]->GetSide() == -3) clusterMass += vec[i]->GetMass()/GeV;
796  clusterMass += backwardEnergy - backwardKinetic;
797 
798  totalEnergy = pseudoParticle[6].GetTotalEnergy()/GeV;
799  pseudoParticle[6].SetMass(clusterMass*GeV);
800 
801  pp = std::sqrt(std::abs(totalEnergy*totalEnergy -
802  clusterMass*clusterMass) )*GeV;
803  pp1 = pseudoParticle[6].GetMomentum().mag();
804  if (pp1 < 1.0e-6*GeV) {
805  G4ThreeVector iso = Isotropic(pp);
806  pseudoParticle[6].SetMomentum(iso.x(), iso.y(), iso.z());
807  } else {
808  pseudoParticle[6].SetMomentum(pseudoParticle[6].GetMomentum() * (-pp/pp1));
809  }
810 
811  std::vector<G4ReactionProduct*> tempList; // ptrs to backward nucleons
812  if (targetParticle.GetSide() == -3) tempList.push_back(&targetParticle);
813  for (i = 0; i < vecLen; ++i)
814  if (vec[i]->GetSide() == -3) tempList.push_back(vec[i]);
815 
816  constantCrossSection = true;
817 
818  if (tempList.size() > 1) {
819  G4int n_entry = 0;
820  wgt = GenerateNBodyEventT(pseudoParticle[6].GetMass(),
821  constantCrossSection, tempList);
822 
823  if (targetParticle.GetSide() == -3) {
824  targetParticle = *tempList[0];
825  targetParticle.Lorentz(targetParticle, pseudoParticle[6]);
826  n_entry++;
827  }
828 
829  for (i = 0; i < vecLen; ++i) {
830  if (vec[i]->GetSide() == -3) {
831  *vec[i] = *tempList[n_entry];
832  vec[i]->Lorentz(*vec[i], pseudoParticle[6]);
833  n_entry++;
834  }
835  }
836  }
837  } else return false;
838 
839  if (vecLen == 0) return false; // all the secondaries have been eliminated
840 
841  // Lorentz transformation to lab frame
842 
843  currentParticle.Lorentz( currentParticle, pseudoParticle[1] );
844  targetParticle.Lorentz( targetParticle, pseudoParticle[1] );
845  for (i = 0; i < vecLen; ++i) vec[i]->Lorentz(*vec[i], pseudoParticle[1]);
846 
847  // Set leading strange particle flag.
848  // leadFlag will be true if original particle and incident particle are
849  // both strange, in which case the incident particle becomes the leading
850  // particle.
851  // leadFlag will also be true if the target particle is strange, but the
852  // incident particle is not, in which case the target particle becomes the
853  // leading particle.
854 
855  G4bool leadingStrangeParticleHasChanged = true;
856  if (leadFlag)
857  {
858  if (currentParticle.GetDefinition() == leadingStrangeParticle.GetDefinition())
859  leadingStrangeParticleHasChanged = false;
860  if (leadingStrangeParticleHasChanged &&
861  (targetParticle.GetDefinition() == leadingStrangeParticle.GetDefinition()) )
862  leadingStrangeParticleHasChanged = false;
863  if( leadingStrangeParticleHasChanged )
864  {
865  for( i=0; i<vecLen; i++ )
866  {
867  if( vec[i]->GetDefinition() == leadingStrangeParticle.GetDefinition() )
868  {
869  leadingStrangeParticleHasChanged = false;
870  break;
871  }
872  }
873  }
874  if( leadingStrangeParticleHasChanged )
875  {
876  G4bool leadTest =
877  (leadingStrangeParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
878  leadingStrangeParticle.GetDefinition()->GetParticleSubType() == "pi");
879  G4bool targetTest =
880  (targetParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
881  targetParticle.GetDefinition()->GetParticleSubType() == "pi");
882 
883  // following modified by JLC 22-Oct-97
884 
885  if( (leadTest&&targetTest) || !(leadTest||targetTest) ) // both true or both false
886  {
887  targetParticle.SetDefinitionAndUpdateE( leadingStrangeParticle.GetDefinition() );
888  targetHasChanged = true;
889  }
890  else
891  {
892  currentParticle.SetDefinitionAndUpdateE( leadingStrangeParticle.GetDefinition() );
893  incidentHasChanged = false;
894  }
895  }
896  } // end of if( leadFlag )
897 
898  // Get number of final state nucleons and nucleons remaining in
899  // target nucleus
900 
901  std::pair<G4int, G4int> finalStateNucleons =
902  GetFinalStateNucleons(originalTarget, vec, vecLen);
903 
904  G4int protonsInFinalState = finalStateNucleons.first;
905  G4int neutronsInFinalState = finalStateNucleons.second;
906 
907  G4int numberofFinalStateNucleons =
908  protonsInFinalState + neutronsInFinalState;
909 
910  if (currentParticle.GetDefinition()->GetBaryonNumber() == 1 &&
911  targetParticle.GetDefinition()->GetBaryonNumber() == 1 &&
912  originalIncident->GetDefinition()->GetPDGMass() <
914  numberofFinalStateNucleons++;
915 
916  numberofFinalStateNucleons = std::max(1, numberofFinalStateNucleons);
917 
918  G4int PinNucleus = std::max(0,
919  G4int(targetNucleus.GetZ_asInt()) - protonsInFinalState);
920  G4int NinNucleus = std::max(0,
921  G4int(targetNucleus.GetA_asInt()-targetNucleus.GetZ_asInt()) - neutronsInFinalState);
922 
923  pseudoParticle[3].SetMomentum( 0.0, 0.0, pOriginal*GeV );
924  pseudoParticle[3].SetMass( mOriginal*GeV );
925  pseudoParticle[3].SetTotalEnergy(
926  std::sqrt( pOriginal*pOriginal + mOriginal*mOriginal )*GeV );
927 
928  const G4ParticleDefinition* aOrgDef = modifiedOriginal.GetDefinition();
929  G4int diff = 0;
930  if(aOrgDef == G4Proton::Proton() || aOrgDef == G4Neutron::Neutron() ) diff = 1;
931  if(numberofFinalStateNucleons == 1) diff = 0;
932  pseudoParticle[4].SetMomentum( 0.0, 0.0, 0.0 );
933  pseudoParticle[4].SetMass( protonMass*(numberofFinalStateNucleons-diff) );
934  pseudoParticle[4].SetTotalEnergy( protonMass*(numberofFinalStateNucleons-diff) );
935 
936  G4double theoreticalKinetic =
937  pseudoParticle[3].GetTotalEnergy() + pseudoParticle[4].GetTotalEnergy() -
938  currentParticle.GetMass() - targetParticle.GetMass();
939  for (i = 0; i < vecLen; ++i) theoreticalKinetic -= vec[i]->GetMass();
940 
941  G4double simulatedKinetic =
942  currentParticle.GetKineticEnergy() + targetParticle.GetKineticEnergy();
943  for (i = 0; i < vecLen; ++i)
944  simulatedKinetic += vec[i]->GetKineticEnergy();
945 
946  pseudoParticle[5] = pseudoParticle[3] + pseudoParticle[4];
947  pseudoParticle[3].Lorentz( pseudoParticle[3], pseudoParticle[5] );
948  pseudoParticle[4].Lorentz( pseudoParticle[4], pseudoParticle[5] );
949 
950  pseudoParticle[7].SetZero();
951  pseudoParticle[7] = pseudoParticle[7] + currentParticle;
952  pseudoParticle[7] = pseudoParticle[7] + targetParticle;
953  for (i = 0; i < vecLen; ++i)
954  pseudoParticle[7] = pseudoParticle[7] + *vec[i];
955 
956  /*
957  // This code does not appear to do anything to the energy or angle spectra
958  if( vecLen <= 16 && vecLen > 0 )
959  {
960  // must create a new set of ReactionProducts here because GenerateNBody will
961  // modify the momenta for the particles, and we don't want to do this
962  //
963  G4ReactionProduct tempR[130];
964  tempR[0] = currentParticle;
965  tempR[1] = targetParticle;
966  for( i=0; i<vecLen; ++i )tempR[i+2] = *vec[i];
967  G4FastVector<G4ReactionProduct,256> tempV1;
968  tempV1.Initialize( vecLen+2 );
969  G4int tempLen1 = 0;
970  for( i=0; i<vecLen+2; ++i )tempV1.SetElement( tempLen1++, &tempR[i] );
971  constantCrossSection = true;
972 
973  wgt = GenerateNBodyEvent(pseudoParticle[3].GetTotalEnergy() +
974  pseudoParticle[4].GetTotalEnergy(),
975  constantCrossSection, tempV1, tempLen1);
976  if (wgt == -1) {
977  G4double Qvalue = 0;
978  for (i = 0; i < tempLen1; i++) Qvalue += tempV1[i]->GetMass();
979  wgt = GenerateNBodyEvent(Qvalue,
980  constantCrossSection, tempV1, tempLen1);
981  }
982  if(wgt>-.5)
983  {
984  theoreticalKinetic = 0.0;
985  for( i=0; i<tempLen1; ++i )
986  {
987  pseudoParticle[6].Lorentz( *tempV1[i], pseudoParticle[4] );
988  theoreticalKinetic += pseudoParticle[6].GetKineticEnergy();
989  }
990  }
991  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
992  }
993  */
994 
995  //
996  // Make sure that the kinetic energies are correct
997  //
998 
999  if (simulatedKinetic != 0.0) {
1000  wgt = theoreticalKinetic/simulatedKinetic;
1001  theoreticalKinetic = currentParticle.GetKineticEnergy() * wgt;
1002  simulatedKinetic = theoreticalKinetic;
1003  currentParticle.SetKineticEnergy(theoreticalKinetic);
1004  pp = currentParticle.GetTotalMomentum();
1005  pp1 = currentParticle.GetMomentum().mag();
1006  if (pp1 < 1.0e-6*GeV) {
1007  G4ThreeVector iso = Isotropic(pp);
1008  currentParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
1009  } else {
1010  currentParticle.SetMomentum(currentParticle.GetMomentum() * (pp/pp1));
1011  }
1012 
1013  theoreticalKinetic = targetParticle.GetKineticEnergy() * wgt;
1014  targetParticle.SetKineticEnergy(theoreticalKinetic);
1015  simulatedKinetic += theoreticalKinetic;
1016  pp = targetParticle.GetTotalMomentum();
1017  pp1 = targetParticle.GetMomentum().mag();
1018 
1019  if (pp1 < 1.0e-6*GeV) {
1020  G4ThreeVector iso = Isotropic(pp);
1021  targetParticle.SetMomentum(iso.x(), iso.y(), iso.z() );
1022  } else {
1023  targetParticle.SetMomentum(targetParticle.GetMomentum() * (pp/pp1) );
1024  }
1025 
1026  for (i = 0; i < vecLen; ++i ) {
1027  theoreticalKinetic = vec[i]->GetKineticEnergy() * wgt;
1028  simulatedKinetic += theoreticalKinetic;
1029  vec[i]->SetKineticEnergy(theoreticalKinetic);
1030  pp = vec[i]->GetTotalMomentum();
1031  pp1 = vec[i]->GetMomentum().mag();
1032  if( pp1 < 1.0e-6*GeV ) {
1033  G4ThreeVector iso = Isotropic(pp);
1034  vec[i]->SetMomentum(iso.x(), iso.y(), iso.z() );
1035  } else {
1036  vec[i]->SetMomentum(vec[i]->GetMomentum() * (pp/pp1) );
1037  }
1038  }
1039  }
1040 
1041  // Rotate(numberofFinalStateNucleons, pseudoParticle[3].GetMomentum(),
1042  // modifiedOriginal, originalIncident, targetNucleus,
1043  // currentParticle, targetParticle, vec, vecLen );
1044 
1045  // Add black track particles
1046  // the total number of particles produced is restricted to 198
1047  // this may have influence on very high energies
1048 
1049  if( atomicWeight >= 1.5 )
1050  {
1051  // npnb is number of proton/neutron black track particles
1052  // ndta is the number of deuterons, tritons, and alphas produced
1053  // epnb is the kinetic energy available for proton/neutron black track
1054  // particles
1055  // edta is the kinetic energy available for deuteron/triton/alpha particles
1056 
1057  G4int npnb = 0;
1058  G4int ndta = 0;
1059 
1060  G4double epnb, edta;
1061  if (veryForward) {
1062  epnb = targetNucleus.GetAnnihilationPNBlackTrackEnergy();
1063  edta = targetNucleus.GetAnnihilationDTABlackTrackEnergy();
1064  } else {
1065  epnb = targetNucleus.GetPNBlackTrackEnergy();
1066  edta = targetNucleus.GetDTABlackTrackEnergy();
1067  }
1068  /*
1069  G4ReactionProduct* fudge = new G4ReactionProduct();
1070  fudge->SetDefinition( aProton );
1071  G4double TT = epnb + edta;
1072  G4double MM = fudge->GetMass()/GeV;
1073  fudge->SetTotalEnergy(MM*GeV + TT*GeV);
1074  G4double pzz = std::sqrt(TT*(TT + 2.*MM));
1075  fudge->SetMomentum(0.0, 0.0, pzz*GeV);
1076  vec.SetElement(vecLen++, fudge);
1077  // G4cout << " Fudge = " << vec[vecLen-1]->GetKineticEnergy()/GeV
1078  << G4endl;
1079  */
1080 
1081  const G4double pnCutOff = 0.001;
1082  const G4double dtaCutOff = 0.001;
1083  // const G4double kineticMinimum = 1.e-6;
1084  // const G4double kineticFactor = -0.010;
1085  // G4double sprob = 0.0; // sprob = probability of self-absorption in
1086  // heavy molecules
1087  // Not currently used (DHW 9 June 2008) const G4double ekIncident = originalIncident->GetKineticEnergy()/GeV;
1088  // if (ekIncident >= 5.0) sprob = std::min(1.0, 0.6*std::log(ekIncident-4.0));
1089  if (epnb > pnCutOff)
1090  {
1091  npnb = G4Poisson((1.5+1.25*numberofFinalStateNucleons)*epnb/(epnb+edta));
1092  if (numberofFinalStateNucleons + npnb > atomicWeight)
1093  npnb = G4int(atomicWeight+0.00001 - numberofFinalStateNucleons);
1094  npnb = std::min( npnb, 127-vecLen );
1095  }
1096  if( edta >= dtaCutOff )
1097  {
1098  ndta = G4Poisson((1.5+1.25*numberofFinalStateNucleons)*edta/(epnb+edta));
1099  ndta = std::min( ndta, 127-vecLen );
1100  }
1101  if (npnb == 0 && ndta == 0) npnb = 1;
1102 
1103  AddBlackTrackParticles(epnb, npnb, edta, ndta, modifiedOriginal,
1104  PinNucleus, NinNucleus, targetNucleus,
1105  vec, vecLen);
1106  }
1107 
1108  // if( centerofmassEnergy <= (4.0+G4UniformRand()) )
1109  // MomentumCheck( modifiedOriginal, currentParticle, targetParticle,
1110  // vec, vecLen );
1111  //
1112  // calculate time delay for nuclear reactions
1113  //
1114 
1115  if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
1116  currentParticle.SetTOF(
1117  1.0-500.0*G4Exp(-ekOriginal/0.04)*G4Log(G4UniformRand()) );
1118  else
1119  currentParticle.SetTOF( 1.0 );
1120  return true;
1121 
1122 }
1123 
1124 
1127  G4double& forwardKinetic,
1128  G4double& backwardKinetic,
1130  G4int& vecLen,
1131  G4ReactionProduct& forwardPseudoParticle,
1132  G4ReactionProduct& backwardPseudoParticle,
1133  G4double& pt)
1134 {
1135  // Reduce energies and pt of secondaries in order to maintain
1136  // energy conservation
1137 
1138  G4double totalEnergy;
1139  G4double pp;
1140  G4double pp1;
1141  G4double px;
1142  G4double py;
1143  G4double mass;
1144  G4ReactionProduct* pVec;
1145  G4int i;
1146 
1147  forwardKinetic = 0.0;
1148  backwardKinetic = 0.0;
1149  forwardPseudoParticle.SetZero();
1150  backwardPseudoParticle.SetZero();
1151 
1152  for (i = startingIndex; i < vecLen; i++) {
1153  pVec = vec[i];
1154  if (pVec->GetSide() != -3) {
1155  mass = pVec->GetMass();
1156  totalEnergy = 0.95*pVec->GetTotalEnergy() + 0.05*mass;
1157  pVec->SetTotalEnergy(totalEnergy);
1158  pp = std::sqrt( std::abs( totalEnergy*totalEnergy - mass*mass ) );
1159  pp1 = pVec->GetMomentum().mag();
1160  if (pp1 < 1.0e-6*GeV) {
1161  G4ThreeVector iso = Isotropic(pp);
1162  pVec->SetMomentum( iso.x(), iso.y(), iso.z() );
1163  } else {
1164  pVec->SetMomentum(pVec->GetMomentum() * (pp/pp1) );
1165  }
1166 
1167  px = pVec->GetMomentum().x();
1168  py = pVec->GetMomentum().y();
1169  pt = std::max(1.0, std::sqrt( px*px + py*py ) )/GeV;
1170  if (pVec->GetSide() > 0) {
1171  forwardKinetic += pVec->GetKineticEnergy()/GeV;
1172  forwardPseudoParticle = forwardPseudoParticle + (*pVec);
1173  } else {
1174  backwardKinetic += pVec->GetKineticEnergy()/GeV;
1175  backwardPseudoParticle = backwardPseudoParticle + (*pVec);
1176  }
1177  }
1178  }
1179 }
1180 
1181  /* end of file */