source: trunk/source/processes/hadronic/models/low_energy/src/G4LELambdaInelastic.cc @ 1337

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

tag geant4.9.4 beta 1 + modifs locales
File size: 13.8 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//
27// $Id: G4LELambdaInelastic.cc,v 1.11 2006/06/29 20:45:05 gunter Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-01 $
29//
30 // Hadronic Process: Lambda Inelastic Process
31 // J.L. Chuma, TRIUMF, 18-Feb-1997
32 // Last modified: 27-Mar-1997
33 // Modified by J.L.Chuma 30-Apr-97: added originalTarget for CalculateMomenta
34 
35#include "G4LELambdaInelastic.hh"
36#include "Randomize.hh"
37
38 G4HadFinalState *
39  G4LELambdaInelastic::ApplyYourself( const G4HadProjectile &aTrack,
40                                      G4Nucleus &targetNucleus )
41  {
42    const G4HadProjectile *originalIncident = &aTrack;
43    //
44    // create the target particle
45    //
46    G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
47   
48    if( verboseLevel > 1 )
49    {
50      const G4Material *targetMaterial = aTrack.GetMaterial();
51      G4cout << "G4LELambdaInelastic::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   
92    G4ReactionProduct currentParticle = modifiedOriginal;
93    G4ReactionProduct targetParticle;
94    targetParticle = *originalTarget;
95    currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
96    targetParticle.SetSide( -1 );  // target always goes in backward hemisphere
97    G4bool incidentHasChanged = false;
98    G4bool targetHasChanged = false;
99    G4bool quasiElastic = false;
100    G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec;  // vec will contain the secondary particles
101    G4int vecLen = 0;
102    vec.Initialize( 0 );
103   
104    const G4double cutOff = 0.1;
105    if( currentParticle.GetKineticEnergy()/MeV > cutOff )
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 void
124  G4LELambdaInelastic::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 original FORTRAN code CASL0 by H. Fesefeldt (13-Sep-1987)
135    //
136    // Lambda 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.70, 0.35 };
166    if( first ) {      // compute normalization constants, this 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( nm=std::max(0,np-2); nm<=(np+1); ++nm ) {
174          for( nz=0; nz<numSec/3; ++nz ) {
175            if( ++counter < numMul ) {
176              nt = np+nm+nz;
177              if( nt>0 && nt<=numSec ) {
178                protmul[counter] = Pmltpc(np,nm,nz,nt,b[0],c);
179                protnorm[nt-1] += protmul[counter];
180              }
181            }
182          }
183        }
184      }
185      for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
186      for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
187      counter = -1;
188      for( np=0; np<numSec/3; ++np ) {
189        for( nm=std::max(0,np-1); nm<=(np+2); ++nm ) {
190          for( nz=0; nz<numSec/3; ++nz ) {
191            if( ++counter < numMul ) {
192              nt = np+nm+nz;
193              if( nt>0 && nt<=numSec ) {
194                neutmul[counter] = Pmltpc(np,nm,nz,nt,b[1],c);
195                neutnorm[nt-1] += neutmul[counter];
196              }
197            }
198          }
199        }
200      }
201      for( i=0; i<numSec; ++i ) {
202        if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
203        if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
204      }
205    }   // end of initialization
206   
207    const G4double expxu = 82.;           // upper bound for arg. of exp
208    const G4double expxl = -expxu;        // lower bound for arg. of exp
209    G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
210    G4ParticleDefinition *aProton = G4Proton::Proton();
211    G4ParticleDefinition *aSigmaPlus = G4SigmaPlus::SigmaPlus();
212    G4ParticleDefinition *aSigmaMinus = G4SigmaMinus::SigmaMinus();
213    G4ParticleDefinition *aSigmaZero = G4SigmaZero::SigmaZero();
214   
215    // energetically possible to produce pion(s)  -->  inelastic scattering
216    //                                   otherwise quasi-elastic scattering
217   
218    G4double n, anpn;
219    GetNormalizationConstant( availableEnergy, n, anpn );
220    G4double ran = G4UniformRand();
221    G4double dum, excs = 0.0;
222    if( targetParticle.GetDefinition() == aProton ) {
223      counter = -1;
224      for( np=0; np<numSec/3 && ran>=excs; ++np ) {
225        for( nm=std::max(0,np-2); nm<=(np+1) && ran>=excs; ++nm ) {
226          for( nz=0; nz<numSec/3 && ran>=excs; ++nz ) {
227            if( ++counter < numMul ) {
228              nt = np+nm+nz;
229              if( nt>0 && nt<=numSec ) {
230                test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
231                dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
232                if( std::fabs(dum) < 1.0 ) {
233                  if( test >= 1.0e-10 )excs += dum*test;
234                } else {
235                  excs += dum*test;
236                }
237              }
238            }
239          }
240        }
241      }
242      if( ran >= excs )  // 3 previous loops continued to the end
243      {
244        quasiElastic = true;
245        return;
246      }
247      np--; nm--; nz--;
248      G4int ncht = std::max( 1, np-nm );
249      switch( ncht ) {
250       case 1:
251         currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
252         incidentHasChanged = true;
253         break;
254       case 2:
255         if( G4UniformRand() < 0.5 ) {
256           if( G4UniformRand() < 0.5 ) {
257             currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
258             incidentHasChanged = true;
259           }
260         } else {
261           currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
262           incidentHasChanged = true;
263           targetParticle.SetDefinitionAndUpdateE( aNeutron );
264           targetHasChanged = true;
265         }
266         break;
267       case 3:
268         if( G4UniformRand() < 0.5 ) {
269           if( G4UniformRand() < 0.5 ) {
270             currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
271             incidentHasChanged = true;
272           }
273           targetParticle.SetDefinitionAndUpdateE( aNeutron );
274           targetHasChanged = true;
275         } else {
276           currentParticle.SetDefinitionAndUpdateE( aSigmaMinus );
277           incidentHasChanged = true;
278         }
279         break;
280       default:
281         currentParticle.SetDefinitionAndUpdateE( aSigmaMinus );
282         incidentHasChanged = true;
283         targetParticle.SetDefinitionAndUpdateE( aNeutron );
284         targetHasChanged = true;
285         break;
286      }
287    }
288    else  // target must be a neutron
289    {
290      counter = -1;
291      for( np=0; np<numSec/3 && ran>=excs; ++np ) {
292        for( nm=std::max(0,np-1); nm<=(np+2) && ran>=excs; ++nm ) {
293          for( nz=0; nz<numSec/3 && ran>=excs; ++nz ) {
294            if( ++counter < numMul ) {
295              nt = np+nm+nz;
296              if( nt>0 && nt<=numSec ) {
297                test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
298                dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
299                if( std::fabs(dum) < 1.0 ) {
300                  if( test >= 1.0e-10 )excs += dum*test;
301                } else {
302                  excs += dum*test;
303                }
304              }
305            }
306          }
307        }
308      }
309      if( ran >= excs )  // 3 previous loops continued to the end
310      {
311        quasiElastic = true;
312        return;
313      }
314      np--; nm--; nz--;
315      G4int ncht = std::max( 1, np-nm+3 );
316      switch( ncht ) {
317       case 1:
318         currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
319         incidentHasChanged = true;
320         targetParticle.SetDefinitionAndUpdateE( aProton );
321         targetHasChanged = true;
322         break;
323       case 2:
324         if( G4UniformRand() < 0.5 ) {
325           if( G4UniformRand() < 0.5 ) {
326             currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
327             incidentHasChanged = true;
328           }
329           targetParticle.SetDefinitionAndUpdateE( aProton );
330           targetHasChanged = true;
331         } else {
332           currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
333           incidentHasChanged = true;
334         }
335         break;
336       case 3:
337         if( G4UniformRand() < 0.5 ) {
338           if( G4UniformRand() < 0.5 ) {
339             currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
340             incidentHasChanged = true;
341           }
342         } else {
343           currentParticle.SetDefinitionAndUpdateE( aSigmaMinus );
344           incidentHasChanged = true;
345           targetParticle.SetDefinitionAndUpdateE( aProton );
346           targetHasChanged = true;
347         }
348         break;
349       default:
350         currentParticle.SetDefinitionAndUpdateE( aSigmaMinus );
351         incidentHasChanged = true;
352         break;
353      }
354    }
355    SetUpPions( np, nm, nz, vec, vecLen );
356    return;
357  }
358
359 /* end of file */
360 
Note: See TracBrowser for help on using the repository browser.