source: trunk/source/processes/hadronic/models/rpg/src/G4RPGSigmaPlusInelastic.cc @ 1228

Last change on this file since 1228 was 1228, checked in by garnier, 14 years ago

update geant4.9.3 tag
File size: 13.3 KB
Line 
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: G4RPGSigmaPlusInelastic.cc,v 1.1 2007/07/18 21:04:21 dennis Exp $
27// GEANT4 tag $Name: geant4-09-03 $
28//
29 
30#include "G4RPGSigmaPlusInelastic.hh"
31#include "Randomize.hh"
32 
33G4HadFinalState*
34G4RPGSigmaPlusInelastic::ApplyYourself( const G4HadProjectile &aTrack,
35                                         G4Nucleus &targetNucleus )
36{
37  const G4HadProjectile *originalIncident = &aTrack;
38  if (originalIncident->GetKineticEnergy()<= 0.1*MeV) 
39  {
40    theParticleChange.SetStatusChange(isAlive);
41    theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
42    theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 
43    return &theParticleChange;     
44  }
45
46  // create the target particle
47
48  G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
49   
50  if( verboseLevel > 1 )
51  {
52    const G4Material *targetMaterial = aTrack.GetMaterial();
53    G4cout << "G4RPGSigmaPlusInelastic::ApplyYourself called" << G4endl;
54    G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
55    G4cout << "target material = " << targetMaterial->GetName() << ", ";
56    G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
57           << G4endl;
58  }
59
60  // Fermi motion and evaporation
61  // As of Geant3, the Fermi energy calculation had not been Done
62
63    G4double ek = originalIncident->GetKineticEnergy()/MeV;
64    G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
65    G4ReactionProduct modifiedOriginal;
66    modifiedOriginal = *originalIncident;
67   
68    G4double tkin = targetNucleus.Cinema( ek );
69    ek += tkin;
70    modifiedOriginal.SetKineticEnergy( ek*MeV );
71    G4double et = ek + amas;
72    G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
73    G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
74    if( pp > 0.0 )
75    {
76      G4ThreeVector momentum = modifiedOriginal.GetMomentum();
77      modifiedOriginal.SetMomentum( momentum * (p/pp) );
78    }
79    //
80    // calculate black track energies
81    //
82    tkin = targetNucleus.EvaporationEffects( ek );
83    ek -= tkin;
84    modifiedOriginal.SetKineticEnergy( ek*MeV );
85    et = ek + amas;
86    p = std::sqrt( std::abs((et-amas)*(et+amas)) );
87    pp = modifiedOriginal.GetMomentum().mag()/MeV;
88    if( pp > 0.0 )
89    {
90      G4ThreeVector momentum = modifiedOriginal.GetMomentum();
91      modifiedOriginal.SetMomentum( momentum * (p/pp) );
92    }
93    G4ReactionProduct currentParticle = modifiedOriginal;
94    G4ReactionProduct targetParticle;
95    targetParticle = *originalTarget;
96    currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
97    targetParticle.SetSide( -1 );  // target always goes in backward hemisphere
98    G4bool incidentHasChanged = false;
99    G4bool targetHasChanged = false;
100    G4bool quasiElastic = false;
101    G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec;  // vec will contain the secondary particles
102    G4int vecLen = 0;
103    vec.Initialize( 0 );
104   
105    const G4double cutOff = 0.1;
106    if( currentParticle.GetKineticEnergy()/MeV > cutOff )
107      Cascade( vec, vecLen,
108               originalIncident, currentParticle, targetParticle,
109               incidentHasChanged, targetHasChanged, quasiElastic );
110   
111    CalculateMomenta( vec, vecLen,
112                      originalIncident, originalTarget, modifiedOriginal,
113                      targetNucleus, currentParticle, targetParticle,
114                      incidentHasChanged, targetHasChanged, quasiElastic );
115   
116    SetUpChange( vec, vecLen,
117                 currentParticle, targetParticle,
118                 incidentHasChanged );
119   
120  delete originalTarget;
121  return &theParticleChange;
122}
123 
124void G4RPGSigmaPlusInelastic::Cascade(
125   G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
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 CASSP
135  //
136  // SigmaPlus 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  {
153    quasiElastic = true;
154    return;
155  }
156    static G4bool first = true;
157    const G4int numMul = 1200;
158    const G4int numSec = 60;
159    static G4double protmul[numMul], protnorm[numSec]; // proton constants
160    static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
161    // np = number of pi+, nm = number of pi-, nz = number of pi0
162    G4int counter, nt=0, np=0, nm=0, nz=0;
163    G4double test;
164    const G4double c = 1.25;   
165    const G4double b[] = { 0.7, 0.7 };
166    if( first )       // compute normalization constants, this will only be Done once
167    {
168      first = false;
169      G4int i;
170      for( i=0; i<numMul; ++i )protmul[i] = 0.0;
171      for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
172      counter = -1;
173      for( np=0; np<(numSec/3); ++np )
174      {
175        for( nm=np; nm<=(np+2); ++nm )
176        {
177          for( nz=0; nz<numSec/3; ++nz )
178          {
179            if( ++counter < numMul )
180            {
181              nt = np+nm+nz;
182              if( nt>0 && nt<=numSec )
183              {
184                protmul[counter] = Pmltpc(np,nm,nz,nt,b[0],c);
185                protnorm[nt-1] += protmul[counter];
186              }
187            }
188          }
189        }
190      }
191      for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
192      for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
193      counter = -1;
194      for( np=0; np<numSec/3; ++np )
195      {
196        for( nm=std::max(0,np-1); nm<=(np+1); ++nm )
197        {
198          for( nz=0; nz<numSec/3; ++nz )
199          {
200            if( ++counter < numMul )
201            {
202              nt = np+nm+nz;
203              if( nt>0 && nt<=numSec )
204              {
205                neutmul[counter] = Pmltpc(np,nm,nz,nt,b[1],c);
206                neutnorm[nt-1] += neutmul[counter];
207              }
208            }
209          }
210        }
211      }
212      for( i=0; i<numSec; ++i )
213      {
214        if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
215        if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
216      }
217    }   // end of initialization
218   
219    const G4double expxu = 82.;           // upper bound for arg. of exp
220    const G4double expxl = -expxu;        // lower bound for arg. of exp
221    G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
222    G4ParticleDefinition *aProton = G4Proton::Proton();
223    G4ParticleDefinition *aLambda = G4Lambda::Lambda();
224    G4ParticleDefinition *aSigmaZero = G4SigmaZero::SigmaZero();
225    //
226    // energetically possible to produce pion(s)  -->  inelastic scattering
227    //
228    G4double n, anpn;
229    GetNormalizationConstant( availableEnergy, n, anpn );
230    G4double ran = G4UniformRand();
231    G4double dum, excs = 0.0;
232    if( targetParticle.GetDefinition() == aProton )
233    {
234      counter = -1;
235      for( np=0; np<numSec/3 && ran>=excs; ++np )
236      {
237        for( nm=np; nm<=(np+2) && ran>=excs; ++nm )
238        {
239          for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
240          {
241            if( ++counter < numMul )
242            {
243              nt = np+nm+nz;
244              if( nt>0 && nt<=numSec )
245              {
246                test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
247                dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
248                if( std::fabs(dum) < 1.0 )
249                {
250                  if( test >= 1.0e-10 )excs += dum*test;
251                }
252                else
253                  excs += dum*test;
254              }
255            }
256          }
257        }
258      }
259      if( ran >= excs )  // 3 previous loops continued to the end
260      {
261        quasiElastic = true;
262        return;
263      }
264      np--; nm--; nz--;
265      switch( std::min( 3, std::max( 1, np-nm+3 ) ) )
266      {
267       case 1:
268         if( G4UniformRand() < 0.5 )
269           currentParticle.SetDefinitionAndUpdateE( aLambda );
270         else
271           currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
272         incidentHasChanged = true;
273         targetParticle.SetDefinitionAndUpdateE( aNeutron );
274         targetHasChanged = true;
275         break;
276       case 2:
277         if( G4UniformRand() < 0.5 )
278         {
279           targetParticle.SetDefinitionAndUpdateE( aNeutron );
280           targetHasChanged = true;
281         }
282         else
283         {
284           if( G4UniformRand() < 0.5 )
285             currentParticle.SetDefinitionAndUpdateE( aLambda );
286           else
287             currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
288           incidentHasChanged = true;
289         }             
290         break;
291       case 3:
292         break;
293      }
294    }
295    else  // target must be a neutron
296    {
297      counter = -1;
298      for( np=0; np<numSec/3 && ran>=excs; ++np )
299      {
300        for( nm=std::max(0,np-1); nm<=(np+1) && ran>=excs; ++nm )
301        {
302          for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
303          {
304            if( ++counter < numMul )
305            {
306              nt = np+nm+nz;
307              if( nt>0 && nt<=numSec )
308              {
309                test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
310                dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
311                if( std::fabs(dum) < 1.0 )
312                {
313                  if( test >= 1.0e-10 )excs += dum*test;
314                }
315                else
316                  excs += dum*test;
317              }
318            }
319          }
320        }
321      }
322      if( ran >= excs )  // 3 previous loops continued to the end
323      {
324        quasiElastic = true;
325        return;
326      }
327      np--; nm--; nz--;
328      switch( std::min( 3, std::max( 1, np-nm+2 ) ) )
329      {
330       case 1:
331         targetParticle.SetDefinitionAndUpdateE( aProton );
332         targetHasChanged = true;
333         break;
334       case 2:
335         if( G4UniformRand() < 0.5 )
336         {
337           if( G4UniformRand() < 0.5 )
338           {
339             currentParticle.SetDefinitionAndUpdateE( aLambda );
340             incidentHasChanged = true;
341             targetParticle.SetDefinitionAndUpdateE( aProton );
342             targetHasChanged = true;
343           }
344           else
345           {
346             targetParticle.SetDefinitionAndUpdateE( aNeutron );
347             targetHasChanged = true;
348           }
349         }
350         else
351         {
352           if( G4UniformRand() < 0.5 )
353           {
354             currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
355             incidentHasChanged = true;
356             targetParticle.SetDefinitionAndUpdateE( aProton );
357             targetHasChanged = true;
358           }
359         }
360         break;
361       case 3:
362         if( G4UniformRand() < 0.5 )
363           currentParticle.SetDefinitionAndUpdateE( aLambda );
364         else
365           currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
366         incidentHasChanged = true;
367         break;
368      }
369    }
370  SetUpPions( np, nm, nz, vec, vecLen );
371  return;
372}
373
374 /* end of file */
375 
Note: See TracBrowser for help on using the repository browser.