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