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