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