ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RPGAntiLambdaInelastic.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RPGAntiLambdaInelastic.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 "G4RPGAntiLambdaInelastic.hh"
29 #include "G4Exp.hh"
30 #include "G4PhysicalConstants.hh"
31 #include "G4SystemOfUnits.hh"
32 #include "Randomize.hh"
33 
34 G4HadFinalState*
35 G4RPGAntiLambdaInelastic::ApplyYourself( const G4HadProjectile &aTrack,
36  G4Nucleus &targetNucleus )
37 {
38  const G4HadProjectile *originalIncident = &aTrack;
39 
40  // Choose the target particle
41 
42  G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
43 
44  if( verboseLevel > 1 )
45  {
46  const G4Material *targetMaterial = aTrack.GetMaterial();
47  G4cout << "G4RPGAntiLambdaInelastic::ApplyYourself called" << G4endl;
48  G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
49  G4cout << "target material = " << targetMaterial->GetName() << ", ";
50  G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
51  << G4endl;
52  }
53 
54  // Fermi motion and evaporation
55  // As of Geant3, the Fermi energy calculation had not been Done
56 
57  G4double ek = originalIncident->GetKineticEnergy()/MeV;
58  G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
59  G4ReactionProduct modifiedOriginal;
60  modifiedOriginal = *originalIncident;
61 
62  G4double tkin = targetNucleus.Cinema( ek );
63  ek += tkin;
64  modifiedOriginal.SetKineticEnergy( ek*MeV );
65  G4double et = ek + amas;
66  G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
67  G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
68  if( pp > 0.0 )
69  {
70  G4ThreeVector momentum = modifiedOriginal.GetMomentum();
71  modifiedOriginal.SetMomentum( momentum * (p/pp) );
72  }
73  //
74  // calculate black track energies
75  //
76  tkin = targetNucleus.EvaporationEffects( ek );
77  ek -= tkin;
78  modifiedOriginal.SetKineticEnergy( ek*MeV );
79  et = ek + amas;
80  p = std::sqrt( std::abs((et-amas)*(et+amas)) );
81  pp = modifiedOriginal.GetMomentum().mag()/MeV;
82  if( pp > 0.0 )
83  {
84  G4ThreeVector momentum = modifiedOriginal.GetMomentum();
85  modifiedOriginal.SetMomentum( momentum * (p/pp) );
86  }
87 
88  G4ReactionProduct currentParticle = modifiedOriginal;
89  G4ReactionProduct targetParticle;
90  targetParticle = *originalTarget;
91  currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
92  targetParticle.SetSide( -1 ); // target always goes in backward hemisphere
93  G4bool incidentHasChanged = false;
94  G4bool targetHasChanged = false;
95  G4bool quasiElastic = false;
96  G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec; // vec will contain the secondary particles
97  G4int vecLen = 0;
98  vec.Initialize( 0 );
99 
100  const G4double cutOff = 0.1;
101  const G4double anni = std::min( 1.3*currentParticle.GetTotalMomentum()/GeV, 0.4 );
102  if( (originalIncident->GetKineticEnergy()/MeV > cutOff) || (G4UniformRand() > anni) )
103  Cascade( vec, vecLen,
104  originalIncident, currentParticle, targetParticle,
105  incidentHasChanged, targetHasChanged, quasiElastic );
106 
107  CalculateMomenta( vec, vecLen,
108  originalIncident, originalTarget, modifiedOriginal,
109  targetNucleus, currentParticle, targetParticle,
110  incidentHasChanged, targetHasChanged, quasiElastic );
111 
112  SetUpChange( vec, vecLen,
113  currentParticle, targetParticle,
114  incidentHasChanged );
115 
116  delete originalTarget;
117  return &theParticleChange;
118 }
119 
120 
121 void G4RPGAntiLambdaInelastic::Cascade(
122  G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
123  G4int &vecLen,
124  const G4HadProjectile *originalIncident,
125  G4ReactionProduct &currentParticle,
126  G4ReactionProduct &targetParticle,
127  G4bool &incidentHasChanged,
128  G4bool &targetHasChanged,
129  G4bool &quasiElastic )
130 {
131  // Derived from H. Fesefeldt's original FORTRAN code CASAL0
132  // AntiLambda undergoes interaction with nucleon within a nucleus. Check if it is
133  // energetically possible to produce pions/kaons. In not, assume nuclear excitation
134  // occurs and input particle is degraded in energy. No other particles are produced.
135  // If reaction is possible, find the correct number of pions/protons/neutrons
136  // produced using an interpolation to multiplicity data. Replace some pions or
137  // protons/neutrons by kaons or strange baryons according to the average
138  // multiplicity per Inelastic reaction.
139 
140  const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass()/MeV;
141  const G4double etOriginal = originalIncident->GetTotalEnergy()/MeV;
142  const G4double targetMass = targetParticle.GetMass()/MeV;
143  const G4double pOriginal = originalIncident->GetTotalMomentum()/GeV;
144  G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
145  targetMass*targetMass +
146  2.0*targetMass*etOriginal );
147  G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
148 
149  static G4ThreadLocal G4bool first = true;
150  const G4int numMul = 1200;
151  const G4int numMulA = 400;
152  const G4int numSec = 60;
153  static G4ThreadLocal G4double protmul[numMul], protnorm[numSec]; // proton constants
154  static G4ThreadLocal G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
155  static G4ThreadLocal G4double protmulA[numMulA], protnormA[numSec]; // proton constants
156  static G4ThreadLocal G4double neutmulA[numMulA], neutnormA[numSec]; // neutron constants
157  // np = number of pi+, nneg = number of pi-, nz = number of pi0
158  G4int nt=0, np=0, nneg=0, nz=0;
159  G4double test;
160  const G4double c = 1.25;
161  const G4double b[] = { 0.7, 0.7 };
162  if( first ) // compute normalization constants, this will only be Done once
163  {
164  first = false;
165  G4int i;
166  for( i=0; i<numMul; ++i )protmul[i] = 0.0;
167  for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
168  G4int counter = -1;
169  for( np=0; np<(numSec/3); ++np )
170  {
171  for( nneg=std::max(0,np-2); nneg<=(np+1); ++nneg )
172  {
173  for( nz=0; nz<numSec/3; ++nz )
174  {
175  if( ++counter < numMul )
176  {
177  nt = np+nneg+nz;
178  if( nt>0 && nt<=numSec )
179  {
180  protmul[counter] = Pmltpc(np,nneg,nz,nt,b[0],c);
181  protnorm[nt-1] += protmul[counter];
182  }
183  }
184  }
185  }
186  }
187  for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
188  for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
189  counter = -1;
190  for( np=0; np<numSec/3; ++np )
191  {
192  for( nneg=std::max(0,np-1); nneg<=(np+2); ++nneg )
193  {
194  for( nz=0; nz<numSec/3; ++nz )
195  {
196  if( ++counter < numMul )
197  {
198  nt = np+nneg+nz;
199  if( nt>0 && nt<=numSec )
200  {
201  neutmul[counter] = Pmltpc(np,nneg,nz,nt,b[1],c);
202  neutnorm[nt-1] += neutmul[counter];
203  }
204  }
205  }
206  }
207  }
208  for( i=0; i<numSec; ++i )
209  {
210  if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
211  if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
212  }
213  //
214  // do the same for annihilation channels
215  //
216  for( i=0; i<numMulA; ++i )protmulA[i] = 0.0;
217  for( i=0; i<numSec; ++i )protnormA[i] = 0.0;
218  counter = -1;
219  for( np=1; np<(numSec/3); ++np )
220  {
221  nneg = np-1;
222  for( nz=0; nz<numSec/3; ++nz )
223  {
224  if( ++counter < numMulA )
225  {
226  nt = np+nneg+nz;
227  if( nt>1 && nt<=numSec )
228  {
229  protmulA[counter] = Pmltpc(np,nneg,nz,nt,b[0],c);
230  protnormA[nt-1] += protmulA[counter];
231  }
232  }
233  }
234  }
235  for( i=0; i<numMulA; ++i )neutmulA[i] = 0.0;
236  for( i=0; i<numSec; ++i )neutnormA[i] = 0.0;
237  counter = -1;
238  for( np=0; np<numSec/3; ++np )
239  {
240  nneg = np;
241  for( nz=0; nz<numSec/3; ++nz )
242  {
243  if( ++counter < numMulA )
244  {
245  nt = np+nneg+nz;
246  if( nt>1 && nt<=numSec )
247  {
248  neutmulA[counter] = Pmltpc(np,nneg,nz,nt,b[1],c);
249  neutnormA[nt-1] += neutmulA[counter];
250  }
251  }
252  }
253  }
254  for( i=0; i<numSec; ++i )
255  {
256  if( protnormA[i] > 0.0 )protnormA[i] = 1.0/protnormA[i];
257  if( neutnormA[i] > 0.0 )neutnormA[i] = 1.0/neutnormA[i];
258  }
259  } // end of initialization
260  const G4double expxu = 82.; // upper bound for arg. of exp
261  const G4double expxl = -expxu; // lower bound for arg. of exp
262 
263  G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
264  G4ParticleDefinition *aProton = G4Proton::Proton();
265  G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
266  G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
267  G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
268  G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus();
269  G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
270  G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong();
271  G4ParticleDefinition *anAntiSigmaZero = G4AntiSigmaZero::AntiSigmaZero();
272  G4ParticleDefinition *anAntiSigmaPlus = G4AntiSigmaPlus::AntiSigmaPlus();
273  G4ParticleDefinition *anAntiSigmaMinus = G4AntiSigmaMinus::AntiSigmaMinus();
274  const G4double anhl[] = {1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,0.97,0.88,
275  0.85,0.81,0.75,0.64,0.64,0.55,0.55,0.45,0.47,0.40,
276  0.39,0.36,0.33,0.10,0.01};
277  G4int iplab = G4int( pOriginal*10.0 );
278  if( iplab > 9 )iplab = G4int( (pOriginal- 1.0)*5.0 ) + 10;
279  if( iplab > 14 )iplab = G4int( pOriginal- 2.0 ) + 15;
280  if( iplab > 22 )iplab = G4int( (pOriginal-10.0)/10.0 ) + 23;
281  if( iplab > 24 )iplab = 24;
282  if( G4UniformRand() > anhl[iplab] )
283  {
284  if( availableEnergy <= aPiPlus->GetPDGMass()/MeV )
285  { // not energetically possible to produce pion(s)
286  quasiElastic = true;
287  return;
288  }
289  G4double n, anpn;
290  GetNormalizationConstant( availableEnergy, n, anpn );
291  G4double ran = G4UniformRand();
292  G4double dum, excs = 0.0;
293  if( targetParticle.GetDefinition() == aProton )
294  {
295  G4int counter = -1;
296  for( np=0; np<numSec/3 && ran>=excs; ++np )
297  {
298  for( nneg=std::max(0,np-2); nneg<=(np+1) && ran>=excs; ++nneg )
299  {
300  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
301  {
302  if( ++counter < numMul )
303  {
304  nt = np+nneg+nz;
305  if( nt>0 && nt<=numSec )
306  {
307  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
308  dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
309  if( std::fabs(dum) < 1.0 )
310  {
311  if( test >= 1.0e-10 )excs += dum*test;
312  }
313  else
314  excs += dum*test;
315  }
316  }
317  }
318  }
319  }
320  if( ran >= excs ) // 3 previous loops continued to the end
321  {
322  quasiElastic = true;
323  return;
324  }
325  np--; nneg--; nz--;
326  G4int ncht = std::min( 4, std::max( 1, np-nneg+2 ) );
327  switch( ncht )
328  {
329  case 1:
330  currentParticle.SetDefinitionAndUpdateE( anAntiSigmaMinus );
331  incidentHasChanged = true;
332  break;
333  case 2:
334  if( G4UniformRand() < 0.5 )
335  {
336  if( G4UniformRand() < 0.5 )
337  {
338  currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
339  incidentHasChanged = true;
340  }
341  else
342  {
343  currentParticle.SetDefinitionAndUpdateE( anAntiSigmaMinus );
344  incidentHasChanged = true;
345  targetParticle.SetDefinitionAndUpdateE( aNeutron );
346  targetHasChanged = true;
347  }
348  }
349  else
350  {
351  if( G4UniformRand() >= 0.5 )
352  {
353  currentParticle.SetDefinitionAndUpdateE( anAntiSigmaMinus );
354  incidentHasChanged = true;
355  targetParticle.SetDefinitionAndUpdateE( aNeutron );
356  targetHasChanged = true;
357  }
358  }
359  break;
360  case 3:
361  if( G4UniformRand() < 0.5 )
362  {
363  if( G4UniformRand() < 0.5 )
364  {
365  currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
366  incidentHasChanged = true;
367  targetParticle.SetDefinitionAndUpdateE( aNeutron );
368  targetHasChanged = true;
369  }
370  else
371  {
372  currentParticle.SetDefinitionAndUpdateE( anAntiSigmaPlus );
373  incidentHasChanged = true;
374  }
375  }
376  else
377  {
378  if( G4UniformRand() < 0.5 )
379  {
380  targetParticle.SetDefinitionAndUpdateE( aNeutron );
381  targetHasChanged = true;
382  }
383  else
384  {
385  currentParticle.SetDefinitionAndUpdateE( anAntiSigmaPlus );
386  incidentHasChanged = true;
387  }
388  }
389  break;
390  case 4:
391  currentParticle.SetDefinitionAndUpdateE( anAntiSigmaPlus );
392  incidentHasChanged = true;
393  targetParticle.SetDefinitionAndUpdateE( aNeutron );
394  targetHasChanged = true;
395  break;
396  }
397  }
398  else // target must be a neutron
399  {
400  G4int counter = -1;
401  for( np=0; np<numSec/3 && ran>=excs; ++np )
402  {
403  for( nneg=std::max(0,np-1); nneg<=(np+2) && ran>=excs; ++nneg )
404  {
405  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
406  {
407  if( ++counter < numMul )
408  {
409  nt = np+nneg+nz;
410  if( nt>0 && nt<=numSec )
411  {
412  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
413  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
414  if( std::fabs(dum) < 1.0 )
415  {
416  if( test >= 1.0e-10 )excs += dum*test;
417  }
418  else
419  excs += dum*test;
420  }
421  }
422  }
423  }
424  }
425  if( ran >= excs ) // 3 previous loops continued to the end
426  {
427  quasiElastic = true;
428  return;
429  }
430  np--; nneg--; nz--;
431  G4int ncht = std::min( 4, std::max( 1, np-nneg+3 ) );
432  switch( ncht )
433  {
434  case 1:
435  currentParticle.SetDefinitionAndUpdateE( anAntiSigmaMinus );
436  incidentHasChanged = true;
437  targetParticle.SetDefinitionAndUpdateE( aProton );
438  targetHasChanged = true;
439  break;
440  case 2:
441  if( G4UniformRand() < 0.5 )
442  {
443  if( G4UniformRand() < 0.5 )
444  {
445  currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
446  incidentHasChanged = true;
447  targetParticle.SetDefinitionAndUpdateE( aProton );
448  targetHasChanged = true;
449  }
450  else
451  {
452  currentParticle.SetDefinitionAndUpdateE( anAntiSigmaMinus );
453  incidentHasChanged = true;
454  }
455  }
456  else
457  {
458  if( G4UniformRand() < 0.5 )
459  {
460  targetParticle.SetDefinitionAndUpdateE( aProton );
461  targetHasChanged = true;
462  }
463  else
464  {
465  currentParticle.SetDefinitionAndUpdateE( anAntiSigmaMinus );
466  incidentHasChanged = true;
467  }
468  }
469  break;
470  case 3:
471  if( G4UniformRand() < 0.5 )
472  {
473  if( G4UniformRand() < 0.5 )
474  {
475  currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
476  incidentHasChanged = true;
477  }
478  else
479  {
480  currentParticle.SetDefinitionAndUpdateE( anAntiSigmaPlus );
481  incidentHasChanged = true;
482  targetParticle.SetDefinitionAndUpdateE( aProton );
483  targetHasChanged = true;
484  }
485  }
486  else
487  {
488  if( G4UniformRand() >= 0.5 )
489  {
490  currentParticle.SetDefinitionAndUpdateE( anAntiSigmaPlus );
491  incidentHasChanged = true;
492  targetParticle.SetDefinitionAndUpdateE( aProton );
493  targetHasChanged = true;
494  }
495  }
496  break;
497  default:
498  currentParticle.SetDefinitionAndUpdateE( anAntiSigmaPlus );
499  incidentHasChanged = true;
500  break;
501  }
502  }
503  }
504  else // random number <= anhl[iplab]
505  {
506  if( centerofmassEnergy <= aPiPlus->GetPDGMass()/MeV+aKaonPlus->GetPDGMass()/MeV )
507  {
508  quasiElastic = true;
509  return;
510  }
511  //
512  // annihilation channels
513  //
514  G4double n, anpn;
515  GetNormalizationConstant( -centerofmassEnergy, n, anpn );
516  G4double ran = G4UniformRand();
517  G4double dum, excs = 0.0;
518  if( targetParticle.GetDefinition() == aProton )
519  {
520  G4int counter = -1;
521  for( np=1; np<numSec/3 && ran>=excs; ++np )
522  {
523  nneg = np-1;
524  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
525  {
526  if( ++counter < numMulA )
527  {
528  nt = np+nneg+nz;
529  if( nt>1 && nt<=numSec )
530  {
531  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
532  dum = (pi/anpn)*nt*protmulA[counter]*protnormA[nt-1]/(2.0*n*n);
533  if( std::fabs(dum) < 1.0 )
534  {
535  if( test >= 1.0e-10 )excs += dum*test;
536  }
537  else
538  excs += dum*test;
539  }
540  }
541  }
542  }
543  }
544  else // target must be a neutron
545  {
546  G4int counter = -1;
547  for( np=0; np<numSec/3 && ran>=excs; ++np )
548  {
549  nneg = np;
550  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
551  {
552  if( ++counter < numMulA )
553  {
554  nt = np+nneg+nz;
555  if( nt>1 && nt<=numSec )
556  {
557  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
558  dum = (pi/anpn)*nt*neutmulA[counter]*neutnormA[nt-1]/(2.0*n*n);
559  if( std::fabs(dum) < 1.0 )
560  {
561  if( test >= 1.0e-10 )excs += dum*test;
562  }
563  else
564  excs += dum*test;
565  }
566  }
567  }
568  }
569  }
570  if( ran >= excs ) // 3 previous loops continued to the end
571  {
572  quasiElastic = true;
573  return;
574  }
575  np--; nz--;
576  currentParticle.SetMass( 0.0 );
577  targetParticle.SetMass( 0.0 );
578  }
579  SetUpPions( np, nneg, nz, vec, vecLen );
580  if( currentParticle.GetMass() == 0.0 )
581  {
582  if( nz == 0 )
583  {
584  if( nneg > 0 )
585  {
586  for( G4int i=0; i<vecLen; ++i )
587  {
588  if( vec[i]->GetDefinition() == aPiMinus )
589  {
590  vec[i]->SetDefinitionAndUpdateE( aKaonMinus );
591  break;
592  }
593  }
594  }
595  }
596  else // nz > 0
597  {
598  if( nneg == 0 )
599  {
600  for( G4int i=0; i<vecLen; ++i )
601  {
602  if( vec[i]->GetDefinition() == aPiZero )
603  {
604  vec[i]->SetDefinitionAndUpdateE( aKaonZL );
605  break;
606  }
607  }
608  }
609  else // nneg > 0
610  {
611  if( G4UniformRand() < 0.5 )
612  {
613  if( nneg > 0 )
614  {
615  for( G4int i=0; i<vecLen; ++i )
616  {
617  if( vec[i]->GetDefinition() == aPiMinus )
618  {
619  vec[i]->SetDefinitionAndUpdateE( aKaonMinus );
620  break;
621  }
622  }
623  }
624  }
625  else // random number >= 0.5
626  {
627  for( G4int i=0; i<vecLen; ++i )
628  {
629  if( vec[i]->GetDefinition() == aPiZero )
630  {
631  vec[i]->SetDefinitionAndUpdateE( aKaonZL );
632  break;
633  }
634  }
635  }
636  }
637  }
638  }
639  return;
640 }
641 
642  /* end of file */
643