ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RPGAntiXiZeroInelastic.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RPGAntiXiZeroInelastic.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 // NOTE: The FORTRAN version of the cascade, CASAXO, simply called the
29 // routine for the XiZero particle. Hence, the ApplyYourself function
30 // below is just a copy of the ApplyYourself from the XiZero particle.
31 
33 #include "G4Exp.hh"
34 #include "G4PhysicalConstants.hh"
35 #include "G4SystemOfUnits.hh"
36 #include "Randomize.hh"
37 
40  G4Nucleus &targetNucleus )
41 {
42  const G4HadProjectile *originalIncident = &aTrack;
43 
44  // Choose the target particle
45 
46  G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
47 
48  if( verboseLevel > 1 )
49  {
50  const G4Material *targetMaterial = aTrack.GetMaterial();
51  G4cout << "G4RPGAntiXiZeroInelastic::ApplyYourself called" << G4endl;
52  G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
53  G4cout << "target material = " << targetMaterial->GetName() << ", ";
54  G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
55  << G4endl;
56  }
57 
58  // Fermi motion and evaporation
59  // As of Geant3, the Fermi energy calculation had not been Done
60 
61  G4double ek = originalIncident->GetKineticEnergy()/MeV;
62  G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
63  G4ReactionProduct modifiedOriginal;
64  modifiedOriginal = *originalIncident;
65 
66  G4double tkin = targetNucleus.Cinema( ek );
67  ek += tkin;
68  modifiedOriginal.SetKineticEnergy( ek*MeV );
69  G4double et = ek + amas;
70  G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
71  G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
72  if( pp > 0.0 )
73  {
74  G4ThreeVector momentum = modifiedOriginal.GetMomentum();
75  modifiedOriginal.SetMomentum( momentum * (p/pp) );
76  }
77  //
78  // calculate black track energies
79  //
80  tkin = targetNucleus.EvaporationEffects( ek );
81  ek -= tkin;
82  modifiedOriginal.SetKineticEnergy( ek*MeV );
83  et = ek + amas;
84  p = std::sqrt( std::abs((et-amas)*(et+amas)) );
85  pp = modifiedOriginal.GetMomentum().mag()/MeV;
86  if( pp > 0.0 )
87  {
88  G4ThreeVector momentum = modifiedOriginal.GetMomentum();
89  modifiedOriginal.SetMomentum( momentum * (p/pp) );
90  }
91  G4ReactionProduct currentParticle = modifiedOriginal;
92  G4ReactionProduct targetParticle;
93  targetParticle = *originalTarget;
94  currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
95  targetParticle.SetSide( -1 ); // target always goes in backward hemisphere
96  G4bool incidentHasChanged = false;
97  G4bool targetHasChanged = false;
98  G4bool quasiElastic = false;
99  G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec; // vec will contain the secondary particles
100  G4int vecLen = 0;
101  vec.Initialize( 0 );
102 
103  const G4double cutOff = 0.1;
104  const G4double anni = std::min( 1.3*currentParticle.GetTotalMomentum()/GeV, 0.4 );
105  if( (currentParticle.GetKineticEnergy()/MeV > cutOff) || (G4UniformRand() > anni) )
106  Cascade( vec, vecLen,
107  originalIncident, currentParticle, targetParticle,
108  incidentHasChanged, targetHasChanged, quasiElastic );
109 
110  CalculateMomenta( vec, vecLen,
111  originalIncident, originalTarget, modifiedOriginal,
112  targetNucleus, currentParticle, targetParticle,
113  incidentHasChanged, targetHasChanged, quasiElastic );
114 
115  SetUpChange( vec, vecLen,
116  currentParticle, targetParticle,
117  incidentHasChanged );
118 
119  delete originalTarget;
120  return &theParticleChange;
121 }
122 
123 
126  G4int& vecLen,
127  const G4HadProjectile *originalIncident,
128  G4ReactionProduct &currentParticle,
129  G4ReactionProduct &targetParticle,
130  G4bool &incidentHasChanged,
131  G4bool &targetHasChanged,
132  G4bool &quasiElastic )
133 {
134  // Derived from H. Fesefeldt's original FORTRAN code CASAX0
135  // which is just a copy of CASX0 (cascade for Xi0)
136  // AntiXiZero undergoes interaction with nucleon within a nucleus. Check if it is
137  // energetically possible to produce pions/kaons. In not, assume nuclear excitation
138  // occurs and input particle is degraded in energy. No other particles are produced.
139  // If reaction is possible, find the correct number of pions/protons/neutrons
140  // produced using an interpolation to multiplicity data. Replace some pions or
141  // protons/neutrons by kaons or strange baryons according to the average
142  // multiplicity per inelastic reaction.
143 
144  const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass()/MeV;
145  const G4double etOriginal = originalIncident->GetTotalEnergy()/MeV;
146  const G4double targetMass = targetParticle.GetMass()/MeV;
147  G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
148  targetMass*targetMass +
149  2.0*targetMass*etOriginal );
150  G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
151  if (availableEnergy <= G4PionPlus::PionPlus()->GetPDGMass()/MeV) {
152  quasiElastic = true;
153  return;
154  }
155  static G4ThreadLocal G4bool first = true;
156  const G4int numMul = 1200;
157  const G4int numSec = 60;
158  static G4ThreadLocal G4double protmul[numMul], protnorm[numSec]; // proton constants
159  static G4ThreadLocal G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
160 
161  // np = number of pi+, nneg = number of pi-, nz = number of pi0
162  G4int counter, nt=0, np=0, nneg=0, nz=0;
163  G4double test;
164  const G4double c = 1.25;
165  const G4double b[] = { 0.7, 0.7 };
166  if (first) { // Computation of normalization constants will only be done once
167  first = false;
168  G4int i;
169  for( i=0; i<numMul; ++i )protmul[i] = 0.0;
170  for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
171  counter = -1;
172  for (np = 0; np < (numSec/3); ++np) {
173  for( nneg=std::max(0,np-2); nneg<=(np+1); ++nneg )
174  {
175  for( nz=0; nz<numSec/3; ++nz )
176  {
177  if( ++counter < numMul )
178  {
179  nt = np+nneg+nz;
180  if( nt>0 && nt<=numSec )
181  {
182  protmul[counter] = Pmltpc(np,nneg,nz,nt,b[0],c);
183  protnorm[nt-1] += protmul[counter];
184  }
185  }
186  }
187  }
188  }
189  for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
190  for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
191  counter = -1;
192  for( np=0; np<numSec/3; ++np )
193  {
194  for( nneg=std::max(0,np-1); nneg<=(np+2); ++nneg )
195  {
196  for( nz=0; nz<numSec/3; ++nz )
197  {
198  if( ++counter < numMul )
199  {
200  nt = np+nneg+nz;
201  if( nt>0 && nt<=numSec )
202  {
203  neutmul[counter] = Pmltpc(np,nneg,nz,nt,b[1],c);
204  neutnorm[nt-1] += neutmul[counter];
205  }
206  }
207  }
208  }
209  }
210  for( i=0; i<numSec; ++i )
211  {
212  if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
213  if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
214  }
215  } // end of initialization
216 
217  const G4double expxu = 82.; // upper bound for arg. of exp
218  const G4double expxl = -expxu; // lower bound for arg. of exp
224  //
225  // energetically possible to produce pion(s) --> inelastic scattering
226  //
227  G4double n, anpn;
228  GetNormalizationConstant( availableEnergy, n, anpn );
230  G4double dum, excs = 0.0;
231  if( targetParticle.GetDefinition() == aProton )
232  {
233  counter = -1;
234  for( np=0; np<numSec/3 && ran>=excs; ++np )
235  {
236  for( nneg=std::max(0,np-2); nneg<=(np+1) && ran>=excs; ++nneg )
237  {
238  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
239  {
240  if( ++counter < numMul )
241  {
242  nt = np+nneg+nz;
243  if( nt>0 && nt<=numSec )
244  {
245  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
246  dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
247  if( std::fabs(dum) < 1.0 )
248  {
249  if( test >= 1.0e-10 )excs += dum*test;
250  }
251  else
252  excs += dum*test;
253  }
254  }
255  }
256  }
257  }
258  if( ran >= excs ) // 3 previous loops continued to the end
259  {
260  quasiElastic = true;
261  return;
262  }
263  np--; nneg--; nz--;
264  //
265  // number of secondary mesons determined by kno distribution
266  // check for total charge of final state mesons to determine
267  // the kind of baryons to be produced, taking into account
268  // charge and strangeness conservation
269  //
270  if( np < nneg+1 )
271  {
272  if( np != nneg ) // charge mismatch
273  {
274  currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
275  incidentHasChanged = true;
276  //
277  // correct the strangeness by replacing a pi- by a kaon-
278  //
279  vec.Initialize( 1 );
281  p->SetDefinition( aKaonMinus );
282  (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
283  vec.SetElement( vecLen++, p );
284  --nneg;
285  }
286  }
287  else if( np == nneg+1 )
288  {
289  if( G4UniformRand() < 0.5 )
290  {
291  targetParticle.SetDefinitionAndUpdateE( aNeutron );
292  targetHasChanged = true;
293  }
294  else
295  {
296  currentParticle.SetDefinitionAndUpdateE( aXiMinus );
297  incidentHasChanged = true;
298  }
299  }
300  else
301  {
302  currentParticle.SetDefinitionAndUpdateE( aXiMinus );
303  incidentHasChanged = true;
304  targetParticle.SetDefinitionAndUpdateE( aNeutron );
305  targetHasChanged = true;
306  }
307  }
308  else // target must be a neutron
309  {
310  counter = -1;
311  for( np=0; np<numSec/3 && ran>=excs; ++np )
312  {
313  for( nneg=std::max(0,np-1); nneg<=(np+2) && ran>=excs; ++nneg )
314  {
315  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
316  {
317  if( ++counter < numMul )
318  {
319  nt = np+nneg+nz;
320  if( nt>0 && nt<=numSec )
321  {
322  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
323  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
324  if( std::fabs(dum) < 1.0 )
325  {
326  if( test >= 1.0e-10 )excs += dum*test;
327  }
328  else
329  excs += dum*test;
330  }
331  }
332  }
333  }
334  }
335  if( ran >= excs ) // 3 previous loops continued to the end
336  {
337  quasiElastic = true;
338  return;
339  }
340  np--; nneg--; nz--;
341  if( np < nneg )
342  {
343  if( np+1 == nneg )
344  {
345  targetParticle.SetDefinitionAndUpdateE( aProton );
346  targetHasChanged = true;
347  }
348  else // charge mismatch
349  {
350  currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
351  incidentHasChanged = true;
352  targetParticle.SetDefinitionAndUpdateE( aProton );
353  targetHasChanged = true;
354  //
355  // correct the strangeness by replacing a pi- by a kaon-
356  //
357  vec.Initialize( 1 );
359  p->SetDefinition( aKaonMinus );
360  (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
361  vec.SetElement( vecLen++, p );
362  --nneg;
363  }
364  }
365  else if( np == nneg )
366  {
367  if( G4UniformRand() >= 0.5 )
368  {
369  currentParticle.SetDefinitionAndUpdateE( aXiMinus );
370  incidentHasChanged = true;
371  targetParticle.SetDefinitionAndUpdateE( aProton );
372  targetHasChanged = true;
373  }
374  }
375  else
376  {
377  currentParticle.SetDefinitionAndUpdateE( aXiMinus );
378  incidentHasChanged = true;
379  }
380  }
381 
382  SetUpPions(np, nneg, nz, vec, vecLen);
383  return;
384 }
385 
386  /* end of file */
387