Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RPGTwoCluster.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$
27 //
28 
29 #include <iostream>
30 #include <signal.h>
31 
32 #include "G4RPGTwoCluster.hh"
33 #include "G4PhysicalConstants.hh"
34 #include "G4SystemOfUnits.hh"
35 #include "Randomize.hh"
36 #include "G4Poisson.hh"
38 
40  : G4RPGReaction() {}
41 
42 
44 ReactionStage(const G4HadProjectile* originalIncident,
45  G4ReactionProduct& modifiedOriginal,
46  G4bool& incidentHasChanged,
47  const G4DynamicParticle* originalTarget,
48  G4ReactionProduct& targetParticle,
49  G4bool& targetHasChanged,
50  const G4Nucleus& targetNucleus,
51  G4ReactionProduct& currentParticle,
53  G4int& vecLen,
54  G4bool leadFlag,
55  G4ReactionProduct& leadingStrangeParticle)
56 {
57  // Derived from H. Fesefeldt's FORTRAN code TWOCLU
58  //
59  // A simple two cluster model is used to generate x- and pt- values for
60  // incident, target, and all secondary particles.
61  // This should be sufficient for low energy interactions.
62 
63  G4int i;
69  G4bool veryForward = false;
70 
71  const G4double protonMass = aProton->GetPDGMass()/MeV;
72  const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/GeV;
73  const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
74  const G4double mOriginal = modifiedOriginal.GetMass()/GeV;
75  const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/GeV;
76  G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV;
77  G4double centerofmassEnergy = std::sqrt(mOriginal*mOriginal +
78  targetMass*targetMass +
79  2.0*targetMass*etOriginal); // GeV
80  G4double currentMass = currentParticle.GetMass()/GeV;
81  targetMass = targetParticle.GetMass()/GeV;
82 
83  if (currentMass == 0.0 && targetMass == 0.0) {
84  G4double ek = currentParticle.GetKineticEnergy();
85  G4ThreeVector mom = currentParticle.GetMomentum();
86  currentParticle = *vec[0];
87  targetParticle = *vec[1];
88  for (i = 0; i < (vecLen-2); ++i) *vec[i] = *vec[i+2];
89  if (vecLen < 2) {
90  for (G4int j = 0; j < vecLen; j++) delete vec[j];
91  vecLen = 0;
92  throw G4HadReentrentException(__FILE__, __LINE__,
93  "G4RPGTwoCluster::ReactionStage : Negative number of particles");
94  }
95  delete vec[vecLen-1];
96  delete vec[vecLen-2];
97  vecLen -= 2;
98  currentMass = currentParticle.GetMass()/GeV;
99  targetMass = targetParticle.GetMass()/GeV;
100  incidentHasChanged = true;
101  targetHasChanged = true;
102  currentParticle.SetKineticEnergy(ek);
103  currentParticle.SetMomentum(mom);
104  veryForward = true;
105  }
106 
107  const G4double atomicWeight = targetNucleus.GetA_asInt();
108  const G4double atomicNumber = targetNucleus.GetZ_asInt();
109 
110  // particles have been distributed in forward and backward hemispheres
111  // in center of mass system of the hadron nucleon interaction
112 
113  // Incident particle always in forward hemisphere
114 
115  G4int forwardCount = 1; // number of particles in forward hemisphere
116  currentParticle.SetSide( 1 );
117  G4double forwardMass = currentParticle.GetMass()/GeV;
118  G4double cMass = forwardMass;
119 
120  // Target particle always in backward hemisphere
121  G4int backwardCount = 1; // number of particles in backward hemisphere
122  targetParticle.SetSide( -1 );
123  G4double backwardMass = targetParticle.GetMass()/GeV;
124  G4double bMass = backwardMass;
125 
126  // G4int backwardNucleonCount = 1; // number of nucleons in backward hemisphere
127  for (i = 0; i < vecLen; ++i) {
128  if (vec[i]->GetSide() < 0) vec[i]->SetSide(-1); // to take care of
129  // case where vec has been preprocessed by GenerateXandPt
130  // and some of them have been set to -2 or -3
131  if (vec[i]->GetSide() == -1) {
132  ++backwardCount;
133  backwardMass += vec[i]->GetMass()/GeV;
134  } else {
135  ++forwardCount;
136  forwardMass += vec[i]->GetMass()/GeV;
137  }
138  }
139 
140  // Add nucleons and some pions from intra-nuclear cascade
141  G4double term1 = std::log(centerofmassEnergy*centerofmassEnergy);
142  if(term1 < 0) term1 = 0.0001; // making sure xtarg<0;
143  const G4double afc = 0.312 + 0.2 * std::log(term1);
144  G4double xtarg;
145  if( centerofmassEnergy < 2.0+G4UniformRand() ) // added +2 below, JLC 4Jul97
146  xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2*backwardCount+vecLen+2)/2.0;
147  else
148  xtarg = afc * (std::pow(atomicWeight,0.33)-1.0) * (2*backwardCount);
149  if( xtarg <= 0.0 )xtarg = 0.01;
150  G4int nuclearExcitationCount = G4Poisson( xtarg );
151 
152  if(atomicWeight<1.0001) nuclearExcitationCount = 0;
153  // G4int extraNucleonCount = 0;
154  // G4double extraMass = 0.0;
155  // G4double extraNucleonMass = 0.0;
156  if( nuclearExcitationCount > 0 )
157  {
158  G4int momentumBin = std::min( 4, G4int(pOriginal/3.0) );
159  const G4double nucsup[] = { 1.0, 0.8, 0.6, 0.5, 0.4 };
160  //
161  // NOTE: in TWOCLU, these new particles were given negative codes
162  // here we use NewlyAdded = true instead
163  //
164  for( i=0; i<nuclearExcitationCount; ++i )
165  {
166  G4ReactionProduct* pVec = new G4ReactionProduct();
167  if( G4UniformRand() < nucsup[momentumBin] ) // add proton or neutron
168  {
169  if( G4UniformRand() > 1.0-atomicNumber/atomicWeight )
170  pVec->SetDefinition( aProton );
171  else
172  pVec->SetDefinition( aNeutron );
173  // Not used ++backwardNucleonCount;
174  // Not used ++extraNucleonCount;
175  // Not used extraNucleonMass += pVec->GetMass()/GeV;
176  }
177  else
178  { // add a pion
179  G4double ran = G4UniformRand();
180  if( ran < 0.3181 )
181  pVec->SetDefinition( aPiPlus );
182  else if( ran < 0.6819 )
183  pVec->SetDefinition( aPiZero );
184  else
185  pVec->SetDefinition( aPiMinus );
186 
187  // DHW: add following two lines to correct energy balance
188  // ++backwardCount;
189  // backwardMass += pVec->GetMass()/GeV;
190  }
191  pVec->SetSide( -2 ); // backside particle
192  // Not used extraMass += pVec->GetMass()/GeV;
193  pVec->SetNewlyAdded( true );
194  vec.SetElement( vecLen++, pVec );
195  }
196  }
197 
198  // Masses of particles added from cascade not included in energy balance.
199  // That's correct for nucleons from the intra-nuclear cascade but not for
200  // pions from the cascade.
201 
202  G4double forwardEnergy = (centerofmassEnergy-cMass-bMass)/2.0 +cMass - forwardMass;
203  G4double backwardEnergy = (centerofmassEnergy-cMass-bMass)/2.0 +bMass - backwardMass;
204  G4double eAvailable = centerofmassEnergy - (forwardMass+backwardMass);
205  G4bool secondaryDeleted;
206  G4double pMass;
207 
208  while( eAvailable <= 0.0 ) // must eliminate a particle
209  {
210  secondaryDeleted = false;
211  for( i=(vecLen-1); i>=0; --i )
212  {
213  if( vec[i]->GetSide() == 1 && vec[i]->GetMayBeKilled())
214  {
215  pMass = vec[i]->GetMass()/GeV;
216  for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
217  --forwardCount;
218  forwardEnergy += pMass;
219  forwardMass -= pMass;
220  secondaryDeleted = true;
221  break;
222  }
223  else if( vec[i]->GetSide() == -1 && vec[i]->GetMayBeKilled())
224  {
225  pMass = vec[i]->GetMass()/GeV;
226  for( G4int j=i; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
227  --backwardCount;
228  backwardEnergy += pMass;
229  backwardMass -= pMass;
230  secondaryDeleted = true;
231  break;
232  }
233  } // breaks go down to here
234 
235  if( secondaryDeleted )
236  {
237  delete vec[vecLen-1];
238  --vecLen;
239  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
240  }
241  else
242  {
243  if( vecLen == 0 ) return false; // all secondaries have been eliminated
244  if( targetParticle.GetSide() == -1 )
245  {
246  pMass = targetParticle.GetMass()/GeV;
247  targetParticle = *vec[0];
248  for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
249  --backwardCount;
250  backwardEnergy += pMass;
251  backwardMass -= pMass;
252  secondaryDeleted = true;
253  }
254  else if( targetParticle.GetSide() == 1 )
255  {
256  pMass = targetParticle.GetMass()/GeV;
257  targetParticle = *vec[0];
258  for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
259  --forwardCount;
260  forwardEnergy += pMass;
261  forwardMass -= pMass;
262  secondaryDeleted = true;
263  }
264 
265  if( secondaryDeleted )
266  {
267  delete vec[vecLen-1];
268  --vecLen;
269  }
270  else
271  {
272  if( currentParticle.GetSide() == -1 )
273  {
274  pMass = currentParticle.GetMass()/GeV;
275  currentParticle = *vec[0];
276  for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
277  --backwardCount;
278  backwardEnergy += pMass;
279  backwardMass -= pMass;
280  secondaryDeleted = true;
281  }
282  else if( currentParticle.GetSide() == 1 )
283  {
284  pMass = currentParticle.GetMass()/GeV;
285  currentParticle = *vec[0];
286  for( G4int j=0; j<(vecLen-1); ++j )*vec[j] = *vec[j+1]; // shift up
287  --forwardCount;
288  forwardEnergy += pMass;
289  forwardMass -= pMass;
290  secondaryDeleted = true;
291  }
292  if( secondaryDeleted )
293  {
294  delete vec[vecLen-1];
295  --vecLen;
296  }
297  else break;
298 
299  } // secondary not deleted
300  } // secondary not deleted
301 
302  eAvailable = centerofmassEnergy - (forwardMass+backwardMass);
303  } // while
304 
305  //
306  // This is the start of the TwoCluster function
307  // Choose multi-particle resonance masses by sampling
308  // P(M) = gc[g(M-M0)]**(c-1) *exp[-(g(M-M0))**c]
309  // for M > M0
310  //
311  // Use for the forward and backward clusters, but not
312  // the cascade cluster
313 
314  const G4double cpar[] = { 1.60, 1.35, 1.15, 1.10 };
315  const G4double gpar[] = { 2.60, 1.80, 1.30, 1.20 };
316  G4int ntc = 0;
317 
318  if (forwardCount < 1 || backwardCount < 1) return false; // array bounds protection
319 
320  G4double rmc = forwardMass;
321  if (forwardCount > 1) {
322  ntc = std::min(3,forwardCount-2);
323  rmc += std::pow(-std::log(1.0-G4UniformRand()),1./cpar[ntc])/gpar[ntc];
324  }
325  G4double rmd = backwardMass;
326  if( backwardCount > 1 ) {
327  ntc = std::min(3,backwardCount-2);
328  rmd += std::pow(-std::log(1.0-G4UniformRand()),1./cpar[ntc])/gpar[ntc];
329  }
330 
331  while( rmc+rmd > centerofmassEnergy )
332  {
333  if( (rmc <= forwardMass) && (rmd <= backwardMass) )
334  {
335  G4double temp = 0.999*centerofmassEnergy/(rmc+rmd);
336  rmc *= temp;
337  rmd *= temp;
338  }
339  else
340  {
341  rmc = 0.1*forwardMass + 0.9*rmc;
342  rmd = 0.1*backwardMass + 0.9*rmd;
343  }
344  }
345 
346  G4ReactionProduct pseudoParticle[8];
347  for( i=0; i<8; ++i )pseudoParticle[i].SetZero();
348 
349  pseudoParticle[1].SetMass( mOriginal*GeV );
350  pseudoParticle[1].SetTotalEnergy( etOriginal*GeV );
351  pseudoParticle[1].SetMomentum( 0.0, 0.0, pOriginal*GeV );
352 
353  pseudoParticle[2].SetMass( protonMass*MeV );
354  pseudoParticle[2].SetTotalEnergy( protonMass*MeV );
355  pseudoParticle[2].SetMomentum( 0.0, 0.0, 0.0 );
356  //
357  // transform into center of mass system
358  //
359  pseudoParticle[0] = pseudoParticle[1] + pseudoParticle[2];
360  pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[0] );
361  pseudoParticle[2].Lorentz( pseudoParticle[2], pseudoParticle[0] );
362 
363  // Calculate cm momentum for forward and backward masses
364  // W = sqrt(pf*pf + rmc*rmc) + sqrt(pf*pf + rmd*rmd)
365  // Solve for pf
366 
367  const G4double pfMin = 0.0001;
368  G4double pf = (centerofmassEnergy*centerofmassEnergy+rmd*rmd-rmc*rmc);
369  pf *= pf;
370  pf -= 4*centerofmassEnergy*centerofmassEnergy*rmd*rmd;
371  pf = std::sqrt( std::max(pf,pfMin) )/(2.0*centerofmassEnergy);
372  //
373  // set final state masses and energies in centre of mass system
374  //
375  pseudoParticle[3].SetMass( rmc*GeV );
376  pseudoParticle[3].SetTotalEnergy( std::sqrt(pf*pf+rmc*rmc)*GeV );
377 
378  pseudoParticle[4].SetMass( rmd*GeV );
379  pseudoParticle[4].SetTotalEnergy( std::sqrt(pf*pf+rmd*rmd)*GeV );
380 
381  //
382  // Get cm scattering angle by sampling t from tmin to tmax
383  //
384  const G4double bMin = 0.01;
385  const G4double b1 = 4.0;
386  const G4double b2 = 1.6;
387  G4double pin = pseudoParticle[1].GetMomentum().mag()/GeV;
388  G4double dtb = 4.0*pin*pf*std::max( bMin, b1+b2*std::log(pOriginal) );
389  G4double factor = 1.0 - std::exp(-dtb);
390  G4double costheta = 1.0 + 2.0*std::log(1.0 - G4UniformRand()*factor) / dtb;
391 
392  costheta = std::max(-1.0, std::min(1.0, costheta));
393  G4double sintheta = std::sqrt((1.0-costheta)*(1.0+costheta));
394  G4double phi = G4UniformRand() * twopi;
395  //
396  // calculate final state momenta in centre of mass system
397  //
398  pseudoParticle[3].SetMomentum( pf*sintheta*std::cos(phi)*GeV,
399  pf*sintheta*std::sin(phi)*GeV,
400  pf*costheta*GeV );
401  pseudoParticle[4].SetMomentum( -pseudoParticle[3].GetMomentum());
402 
403  // Backward cluster of nucleons and pions from intra-nuclear cascade
404  // Set up in lab system and transform to cms
405 
406  G4double pp, pp1;
407  if( nuclearExcitationCount > 0 )
408  {
409  const G4double ga = 1.2;
410  G4double ekit1 = 0.04;
411  G4double ekit2 = 0.6; // Max KE of cascade particle
412  if( ekOriginal <= 5.0 )
413  {
414  ekit1 *= ekOriginal*ekOriginal/25.0;
415  ekit2 *= ekOriginal*ekOriginal/25.0;
416  }
417  G4double scale = std::pow(ekit2/ekit1, 1.0-ga) - 1.0;
418  for( i=0; i<vecLen; ++i )
419  {
420  if( vec[i]->GetSide() == -2 )
421  {
422  G4double kineticE = ekit1*std::pow((1.0 + G4UniformRand()*scale), 1.0/(1.0-ga) );
423  vec[i]->SetKineticEnergy( kineticE*GeV );
424  G4double vMass = vec[i]->GetMass()/MeV;
425  G4double totalE = kineticE*GeV + vMass;
426  pp = std::sqrt( std::abs(totalE*totalE-vMass*vMass) );
427  G4double cost = std::min( 1.0, std::max( -1.0, std::log(2.23*G4UniformRand()+0.383)/0.96 ) );
428  G4double sint = std::sqrt(1.0-cost*cost);
429  phi = twopi*G4UniformRand();
430  vec[i]->SetMomentum( pp*sint*std::cos(phi)*MeV,
431  pp*sint*std::sin(phi)*MeV,
432  pp*cost*MeV );
433  vec[i]->Lorentz( *vec[i], pseudoParticle[0] );
434  }
435  }
436  }
437 
438  //
439  // Fragmentation of forward and backward clusters
440  //
441 
442  currentParticle.SetMomentum( pseudoParticle[3].GetMomentum() );
443  currentParticle.SetTotalEnergy( pseudoParticle[3].GetTotalEnergy() );
444 
445  targetParticle.SetMomentum( pseudoParticle[4].GetMomentum() );
446  targetParticle.SetTotalEnergy( pseudoParticle[4].GetTotalEnergy() );
447 
448  pseudoParticle[5].SetMomentum( pseudoParticle[3].GetMomentum() * (-1.0) );
449  pseudoParticle[5].SetMass( pseudoParticle[3].GetMass() );
450  pseudoParticle[5].SetTotalEnergy( pseudoParticle[3].GetTotalEnergy() );
451 
452  pseudoParticle[6].SetMomentum( pseudoParticle[4].GetMomentum() * (-1.0) );
453  pseudoParticle[6].SetMass( pseudoParticle[4].GetMass() );
454  pseudoParticle[6].SetTotalEnergy( pseudoParticle[4].GetTotalEnergy() );
455 
456  G4double wgt;
457  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
458  if( forwardCount > 1 ) // tempV will contain the forward particles
459  {
461  tempV.Initialize( forwardCount );
462  G4bool constantCrossSection = true;
463  G4int tempLen = 0;
464  if( currentParticle.GetSide() == 1 )
465  tempV.SetElement( tempLen++, &currentParticle );
466  if( targetParticle.GetSide() == 1 )
467  tempV.SetElement( tempLen++, &targetParticle );
468  for( i=0; i<vecLen; ++i )
469  {
470  if( vec[i]->GetSide() == 1 )
471  {
472  if( tempLen < 18 )
473  tempV.SetElement( tempLen++, vec[i] );
474  else
475  {
476  vec[i]->SetSide( -1 );
477  continue;
478  }
479  }
480  }
481  if( tempLen >= 2 )
482  {
483  wgt = GenerateNBodyEvent( pseudoParticle[3].GetMass()/MeV,
484  constantCrossSection, tempV, tempLen );
485  if( currentParticle.GetSide() == 1 )
486  currentParticle.Lorentz( currentParticle, pseudoParticle[5] );
487  if( targetParticle.GetSide() == 1 )
488  targetParticle.Lorentz( targetParticle, pseudoParticle[5] );
489  for( i=0; i<vecLen; ++i )
490  {
491  if( vec[i]->GetSide() == 1 )vec[i]->Lorentz( *vec[i], pseudoParticle[5] );
492  }
493  }
494  }
495  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
496  if( backwardCount > 1 ) // tempV will contain the backward particles,
497  { // but not those created from the intranuclear cascade
499  tempV.Initialize( backwardCount );
500  G4bool constantCrossSection = true;
501  G4int tempLen = 0;
502  if( currentParticle.GetSide() == -1 )
503  tempV.SetElement( tempLen++, &currentParticle );
504  if( targetParticle.GetSide() == -1 )
505  tempV.SetElement( tempLen++, &targetParticle );
506  for( i=0; i<vecLen; ++i )
507  {
508  if( vec[i]->GetSide() == -1 )
509  {
510  if( tempLen < 18 )
511  tempV.SetElement( tempLen++, vec[i] );
512  else
513  {
514  vec[i]->SetSide( -2 );
515  vec[i]->SetKineticEnergy( 0.0 );
516  vec[i]->SetMomentum( 0.0, 0.0, 0.0 );
517  continue;
518  }
519  }
520  }
521  if( tempLen >= 2 )
522  {
523  wgt = GenerateNBodyEvent( pseudoParticle[4].GetMass()/MeV,
524  constantCrossSection, tempV, tempLen );
525  if( currentParticle.GetSide() == -1 )
526  currentParticle.Lorentz( currentParticle, pseudoParticle[6] );
527  if( targetParticle.GetSide() == -1 )
528  targetParticle.Lorentz( targetParticle, pseudoParticle[6] );
529  for( i=0; i<vecLen; ++i )
530  {
531  if( vec[i]->GetSide() == -1 )vec[i]->Lorentz( *vec[i], pseudoParticle[6] );
532  }
533  }
534  }
535 
536  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
537  //
538  // Lorentz transformation in lab system
539  //
540  currentParticle.Lorentz( currentParticle, pseudoParticle[2] );
541  targetParticle.Lorentz( targetParticle, pseudoParticle[2] );
542  for( i=0; i<vecLen; ++i ) vec[i]->Lorentz( *vec[i], pseudoParticle[2] );
543 
544  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
545  //
546  // sometimes the leading strange particle is lost, set it back
547  //
548  G4bool dum = true;
549  if( leadFlag )
550  {
551  // leadFlag will be true
552  // iff original particle is strange AND if incident particle is strange
553  // leadFlag is set to the incident particle
554  // or
555  // target particle is strange leadFlag is set to the target particle
556 
557  if( currentParticle.GetDefinition() == leadingStrangeParticle.GetDefinition() )
558  dum = false;
559  else if( targetParticle.GetDefinition() == leadingStrangeParticle.GetDefinition() )
560  dum = false;
561  else
562  {
563  for( i=0; i<vecLen; ++i )
564  {
565  if( vec[i]->GetDefinition() == leadingStrangeParticle.GetDefinition() )
566  {
567  dum = false;
568  break;
569  }
570  }
571  }
572  if( dum )
573  {
574  G4double leadMass = leadingStrangeParticle.GetMass()/MeV;
575  G4double ekin;
576  if( ((leadMass < protonMass) && (targetParticle.GetMass()/MeV < protonMass)) ||
577  ((leadMass >= protonMass) && (targetParticle.GetMass()/MeV >= protonMass)) )
578  {
579  ekin = targetParticle.GetKineticEnergy()/GeV;
580  pp1 = targetParticle.GetMomentum().mag()/MeV; // old momentum
581  targetParticle.SetDefinition( leadingStrangeParticle.GetDefinition() );
582  targetParticle.SetKineticEnergy( ekin*GeV );
583  pp = targetParticle.GetTotalMomentum()/MeV; // new momentum
584  if( pp1 < 1.0e-3 ) {
585  G4ThreeVector iso = Isotropic(pp);
586  targetParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
587  } else {
588  targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
589  }
590  targetHasChanged = true;
591  }
592  else
593  {
594  ekin = currentParticle.GetKineticEnergy()/GeV;
595  pp1 = currentParticle.GetMomentum().mag()/MeV;
596  currentParticle.SetDefinition( leadingStrangeParticle.GetDefinition() );
597  currentParticle.SetKineticEnergy( ekin*GeV );
598  pp = currentParticle.GetTotalMomentum()/MeV;
599  if( pp1 < 1.0e-3 ) {
600  G4ThreeVector iso = Isotropic(pp);
601  currentParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
602  } else {
603  currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
604  }
605  incidentHasChanged = true;
606  }
607  }
608  } // end of if( leadFlag )
609 
610  // Get number of final state nucleons and nucleons remaining in
611  // target nucleus
612 
613  std::pair<G4int, G4int> finalStateNucleons =
614  GetFinalStateNucleons(originalTarget, vec, vecLen);
615 
616  G4int protonsInFinalState = finalStateNucleons.first;
617  G4int neutronsInFinalState = finalStateNucleons.second;
618 
619  G4int numberofFinalStateNucleons =
620  protonsInFinalState + neutronsInFinalState;
621 
622  if (currentParticle.GetDefinition()->GetBaryonNumber() == 1 &&
623  targetParticle.GetDefinition()->GetBaryonNumber() == 1 &&
624  originalIncident->GetDefinition()->GetPDGMass() <
626  numberofFinalStateNucleons++;
627 
628  numberofFinalStateNucleons = std::max(1, numberofFinalStateNucleons);
629 
630  G4int PinNucleus = std::max(0,
631  G4int(targetNucleus.GetZ_asInt()) - protonsInFinalState);
632  G4int NinNucleus = std::max(0,
633  G4int(targetNucleus.GetA_asInt()-targetNucleus.GetZ_asInt()) - neutronsInFinalState);
634  //
635  // for various reasons, the energy balance is not sufficient,
636  // check that, energy balance, angle of final system, etc.
637  //
638  pseudoParticle[4].SetMass( mOriginal*GeV );
639  pseudoParticle[4].SetTotalEnergy( etOriginal*GeV );
640  pseudoParticle[4].SetMomentum( 0.0, 0.0, pOriginal*GeV );
641 
642  G4ParticleDefinition * aOrgDef = modifiedOriginal.GetDefinition();
643  G4int diff = 0;
644  if(aOrgDef == G4Proton::Proton() || aOrgDef == G4Neutron::Neutron() ) diff = 1;
645  if(numberofFinalStateNucleons == 1) diff = 0;
646  pseudoParticle[5].SetMomentum( 0.0, 0.0, 0.0 );
647  pseudoParticle[5].SetMass( protonMass*(numberofFinalStateNucleons-diff)*MeV);
648  pseudoParticle[5].SetTotalEnergy( protonMass*(numberofFinalStateNucleons-diff)*MeV);
649 
650  G4double theoreticalKinetic =
651  pseudoParticle[4].GetTotalEnergy()/GeV + pseudoParticle[5].GetTotalEnergy()/GeV;
652 
653  pseudoParticle[6] = pseudoParticle[4] + pseudoParticle[5];
654  pseudoParticle[4].Lorentz( pseudoParticle[4], pseudoParticle[6] );
655  pseudoParticle[5].Lorentz( pseudoParticle[5], pseudoParticle[6] );
656 
657  if( vecLen < 16 )
658  {
659  G4ReactionProduct tempR[130];
660  tempR[0] = currentParticle;
661  tempR[1] = targetParticle;
662  for( i=0; i<vecLen; ++i )tempR[i+2] = *vec[i];
663 
665  tempV.Initialize( vecLen+2 );
666  G4bool constantCrossSection = true;
667  G4int tempLen = 0;
668  for( i=0; i<vecLen+2; ++i )tempV.SetElement( tempLen++, &tempR[i] );
669 
670  if( tempLen >= 2 )
671  {
672  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
673  wgt = GenerateNBodyEvent( pseudoParticle[4].GetTotalEnergy()/MeV +
674  pseudoParticle[5].GetTotalEnergy()/MeV,
675  constantCrossSection, tempV, tempLen );
676  if (wgt == -1) {
677  G4double Qvalue = 0;
678  for (i = 0; i < tempLen; i++) Qvalue += tempV[i]->GetMass();
679  wgt = GenerateNBodyEvent( Qvalue/MeV,
680  constantCrossSection, tempV, tempLen );
681  }
682  theoreticalKinetic = 0.0;
683  for( i=0; i<vecLen+2; ++i )
684  {
685  pseudoParticle[7].SetMomentum( tempV[i]->GetMomentum() );
686  pseudoParticle[7].SetMass( tempV[i]->GetMass() );
687  pseudoParticle[7].SetTotalEnergy( tempV[i]->GetTotalEnergy() );
688  pseudoParticle[7].Lorentz( pseudoParticle[7], pseudoParticle[5] );
689  theoreticalKinetic += pseudoParticle[7].GetKineticEnergy()/GeV;
690  }
691  }
692  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
693  }
694  else
695  {
696  theoreticalKinetic -=
697  ( currentParticle.GetMass()/GeV + targetParticle.GetMass()/GeV );
698  for( i=0; i<vecLen; ++i )theoreticalKinetic -= vec[i]->GetMass()/GeV;
699  }
700  G4double simulatedKinetic =
701  currentParticle.GetKineticEnergy()/GeV + targetParticle.GetKineticEnergy()/GeV;
702  for( i=0; i<vecLen; ++i )simulatedKinetic += vec[i]->GetKineticEnergy()/GeV;
703 
704  // make sure that kinetic energies are correct
705  // the backward nucleon cluster is not produced within proper kinematics!!!
706 
707  if( simulatedKinetic != 0.0 )
708  {
709  wgt = (theoreticalKinetic)/simulatedKinetic;
710  currentParticle.SetKineticEnergy( wgt*currentParticle.GetKineticEnergy() );
711  pp = currentParticle.GetTotalMomentum()/MeV;
712  pp1 = currentParticle.GetMomentum().mag()/MeV;
713  if( pp1 < 0.001*MeV ) {
714  G4ThreeVector iso = Isotropic(pp);
715  currentParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
716  } else {
717  currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) );
718  }
719 
720  targetParticle.SetKineticEnergy( wgt*targetParticle.GetKineticEnergy() );
721  pp = targetParticle.GetTotalMomentum()/MeV;
722  pp1 = targetParticle.GetMomentum().mag()/MeV;
723  if( pp1 < 0.001*MeV ) {
724  G4ThreeVector iso = Isotropic(pp);
725  targetParticle.SetMomentum( iso.x(), iso.y(), iso.z() );
726  } else {
727  targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) );
728  }
729 
730  for( i=0; i<vecLen; ++i )
731  {
732  vec[i]->SetKineticEnergy( wgt*vec[i]->GetKineticEnergy() );
733  pp = vec[i]->GetTotalMomentum()/MeV;
734  pp1 = vec[i]->GetMomentum().mag()/MeV;
735  if( pp1 < 0.001 ) {
736  G4ThreeVector iso = Isotropic(pp);
737  vec[i]->SetMomentum( iso.x(), iso.y(), iso.z() );
738  } else {
739  vec[i]->SetMomentum( vec[i]->GetMomentum() * (pp/pp1) );
740  }
741  }
742  }
743  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
744 
745  Rotate( numberofFinalStateNucleons, pseudoParticle[4].GetMomentum(),
746  modifiedOriginal, originalIncident, targetNucleus,
747  currentParticle, targetParticle, vec, vecLen );
748 
749  // Add black track particles
750  // the total number of particles produced is restricted to 198
751  // this may have influence on very high energies
752 
753  if( atomicWeight >= 1.5 )
754  {
755  // npnb is number of proton/neutron black track particles
756  // ndta is the number of deuterons, tritons, and alphas produced
757  // epnb is the kinetic energy available for proton/neutron black track
758  // particles
759  // edta is the kinetic energy available for deuteron/triton/alpha
760  // particles
761 
762  G4int npnb = 0;
763  G4int ndta = 0;
764 
765  G4double epnb, edta;
766  if (veryForward) {
767  epnb = targetNucleus.GetAnnihilationPNBlackTrackEnergy();
768  edta = targetNucleus.GetAnnihilationDTABlackTrackEnergy();
769  } else {
770  epnb = targetNucleus.GetPNBlackTrackEnergy();
771  edta = targetNucleus.GetDTABlackTrackEnergy();
772  }
773 
774  const G4double pnCutOff = 0.001; // GeV
775  const G4double dtaCutOff = 0.001; // GeV
776  // const G4double kineticMinimum = 1.e-6;
777  // const G4double kineticFactor = -0.005;
778 
779  // G4double sprob = 0.0; // sprob = probability of self-absorption in
780  // heavy molecules
781  // Not currently used (DHW 9 June 2008) const G4double ekIncident = originalIncident->GetKineticEnergy()/GeV;
782  // if( ekIncident >= 5.0 )sprob = std::min( 1.0, 0.6*std::log(ekIncident-4.0) );
783 
784  if( epnb >= pnCutOff )
785  {
786  npnb = G4Poisson((1.5+1.25*numberofFinalStateNucleons)*epnb/(epnb+edta));
787  if( numberofFinalStateNucleons + npnb > atomicWeight )
788  npnb = G4int(atomicWeight - numberofFinalStateNucleons);
789  npnb = std::min( npnb, 127-vecLen );
790  }
791  if( edta >= dtaCutOff )
792  {
793  ndta = G4Poisson( (1.5+1.25*numberofFinalStateNucleons)*edta/(epnb+edta) );
794  ndta = std::min( ndta, 127-vecLen );
795  }
796  if (npnb == 0 && ndta == 0) npnb = 1;
797 
798  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
799 
800  AddBlackTrackParticles(epnb, npnb, edta, ndta, modifiedOriginal,
801  PinNucleus, NinNucleus, targetNucleus,
802  vec, vecLen );
803  // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
804  }
805 
806  //if( centerofmassEnergy <= (4.0+G4UniformRand()) )
807  // MomentumCheck( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
808  //
809  // calculate time delay for nuclear reactions
810  //
811  if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) )
812  currentParticle.SetTOF( 1.0-500.0*std::exp(-ekOriginal/0.04)*std::log(G4UniformRand()) );
813  else
814  currentParticle.SetTOF( 1.0 );
815 
816  return true;
817 }
818 
819  /* end of file */