source: trunk/source/processes/hadronic/models/rpg/src/G4RPGKPlusInelastic.cc @ 962

Last change on this file since 962 was 962, checked in by garnier, 15 years ago

update processes
File size: 14.0 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: G4RPGKPlusInelastic.cc,v 1.1 2007/07/18 21:04:20 dennis Exp $
28// GEANT4 tag $Name: geant4-09-02-ref-02 $
29//
30 
31#include "G4RPGKPlusInelastic.hh"
32#include "Randomize.hh"
33 
34G4HadFinalState*
35G4RPGKPlusInelastic::ApplyYourself( const G4HadProjectile &aTrack,
36                                        G4Nucleus &targetNucleus )
37{
38  const G4HadProjectile *originalIncident = &aTrack;
39  if (originalIncident->GetKineticEnergy()<= 0.1*MeV) 
40  {
41    theParticleChange.SetStatusChange(isAlive);
42    theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
43    theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 
44    return &theParticleChange;     
45  }
46   
47  // create the target particle
48   
49  G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
50  G4ReactionProduct targetParticle( originalTarget->GetDefinition() );
51   
52  if( verboseLevel > 1 )
53  {
54    const G4Material *targetMaterial = aTrack.GetMaterial();
55    G4cout << "G4RPGKPlusInelastic::ApplyYourself called" << G4endl;
56    G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy() << "MeV, ";
57    G4cout << "target material = " << targetMaterial->GetName() << ", ";
58    G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
59           << G4endl;
60  }   
61  G4ReactionProduct currentParticle( const_cast<G4ParticleDefinition *>(originalIncident->GetDefinition()));
62  currentParticle.SetMomentum( originalIncident->Get4Momentum().vect() );
63  currentParticle.SetKineticEnergy( originalIncident->GetKineticEnergy() );
64   
65  // Fermi motion and evaporation
66  // As of Geant3, the Fermi energy calculation had not been Done
67   
68  G4double ek = originalIncident->GetKineticEnergy();
69  G4double amas = originalIncident->GetDefinition()->GetPDGMass();
70   
71  G4double tkin = targetNucleus.Cinema( ek );
72  ek += tkin;
73  currentParticle.SetKineticEnergy( ek );
74  G4double et = ek + amas;
75  G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
76  G4double pp = currentParticle.GetMomentum().mag();
77  if( pp > 0.0 )
78  {
79    G4ThreeVector momentum = currentParticle.GetMomentum();
80    currentParticle.SetMomentum( momentum * (p/pp) );
81  }
82   
83  // calculate black track energies
84   
85  tkin = targetNucleus.EvaporationEffects( ek );
86  ek -= tkin;
87  currentParticle.SetKineticEnergy( ek );
88  et = ek + amas;
89  p = std::sqrt( std::abs((et-amas)*(et+amas)) );
90  pp = currentParticle.GetMomentum().mag();
91  if( pp > 0.0 )
92  {
93    G4ThreeVector momentum = currentParticle.GetMomentum();
94    currentParticle.SetMomentum( momentum * (p/pp) );
95  }
96   
97  G4ReactionProduct modifiedOriginal = currentParticle;
98   
99  currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
100  targetParticle.SetSide( -1 );  // target always goes in backward hemisphere
101  G4bool incidentHasChanged = false;
102  G4bool targetHasChanged = false;
103  G4bool quasiElastic = false;
104  G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec;  // vec will contain the secondary particles
105  G4int vecLen = 0;
106  vec.Initialize( 0 );
107       
108  const G4double cutOff = 0.1*MeV;
109  if( currentParticle.GetKineticEnergy() > cutOff )
110    Cascade( vec, vecLen,
111             originalIncident, currentParticle, targetParticle,
112             incidentHasChanged, targetHasChanged, quasiElastic );
113   
114  CalculateMomenta( vec, vecLen,
115                    originalIncident, originalTarget, modifiedOriginal,
116                    targetNucleus, currentParticle, targetParticle,
117                    incidentHasChanged, targetHasChanged, quasiElastic );
118   
119  SetUpChange( vec, vecLen,
120               currentParticle, targetParticle,
121               incidentHasChanged );
122   
123  delete originalTarget;
124   
125  return &theParticleChange;   
126}
127
128
129void G4RPGKPlusInelastic::Cascade(
130   G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
131   G4int &vecLen,
132   const G4HadProjectile *originalIncident,
133   G4ReactionProduct &currentParticle,
134   G4ReactionProduct &targetParticle,
135   G4bool &incidentHasChanged,
136   G4bool &targetHasChanged,
137   G4bool &quasiElastic )
138{
139  // Derived from H. Fesefeldt's original FORTRAN code CASKP
140  //
141  // K+ undergoes interaction with nucleon within a nucleus.  Check if it is
142  // energetically possible to produce pions/kaons.  In not, assume nuclear excitation
143  // occurs and input particle is degraded in energy. No other particles are produced.
144  // If reaction is possible, find the correct number of pions/protons/neutrons
145  // produced using an interpolation to multiplicity data.  Replace some pions or
146  // protons/neutrons by kaons or strange baryons according to the average
147  // multiplicity per Inelastic reaction.
148  //
149  const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass();
150  const G4double etOriginal = originalIncident->GetTotalEnergy();
151  const G4double targetMass = targetParticle.GetMass();
152  G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
153                                        targetMass*targetMass +
154                                        2.0*targetMass*etOriginal );
155  G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
156  if( availableEnergy < G4PionPlus::PionPlus()->GetPDGMass() )
157  {
158    quasiElastic = true;
159    return;
160  }
161  static G4bool first = true;
162  const G4int numMul = 1200;
163  const G4int numSec = 60;
164  static G4double protmul[numMul], protnorm[numSec]; // proton constants
165  static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
166
167  // np = number of pi+, nm = number of pi-, nz = number of pi0
168
169  G4int nt=0, np=0, nm=0, nz=0;
170  const G4double c = 1.25;   
171  const G4double b[] = { 0.70, 0.70 };
172  if( first )       // compute normalization constants, this will only be Done once
173  {
174      first = false;
175      G4int i;
176      for( i=0; i<numMul; ++i )protmul[i] = 0.0;
177      for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
178      G4int counter = -1;
179      for( np=0; np<(numSec/3); ++np )
180      {
181        for( nm=std::max(0,np-2); nm<=np; ++nm )
182        {
183          for( nz=0; nz<numSec/3; ++nz )
184          {
185            if( ++counter < numMul )
186            {
187              nt = np+nm+nz;
188              if( nt > 0 )
189              {
190                protmul[counter] = Pmltpc(np,nm,nz,nt,b[0],c);
191                protnorm[nt-1] += protmul[counter];
192              }
193            }
194          }
195        }
196      }
197      for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
198      for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
199      counter = -1;
200      for( np=0; np<numSec/3; ++np )
201      {
202        for( nm=std::max(0,np-1); nm<=(np+1); ++nm )
203        {
204          for( nz=0; nz<numSec/3; ++nz )
205          {
206            if( ++counter < numMul )
207            {
208              nt = np+nm+nz;
209              if( (nt>0) && (nt<=numSec) )
210              {
211                neutmul[counter] = Pmltpc(np,nm,nz,nt,b[1],c);
212                neutnorm[nt-1] += neutmul[counter];
213              }
214            }
215          }
216        }
217      }
218      for( i=0; i<numSec; ++i )
219      {
220        if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
221        if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
222      }
223  }   // end of initialization
224       
225  const G4double expxu = 82.;           // upper bound for arg. of exp
226  const G4double expxl = -expxu;        // lower bound for arg. of exp
227  G4ParticleDefinition *aKaonZS = G4KaonZeroShort::KaonZeroShort();
228  G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong();
229  G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
230  G4ParticleDefinition *aProton = G4Proton::Proton();
231  G4int ieab = static_cast<G4int>(availableEnergy*5.0/GeV);
232  const G4double supp[] = {0.,0.4,0.55,0.65,0.75,0.82,0.86,0.90,0.94,0.98};
233  G4double test, w0, wp, wt, wm;
234  if( (availableEnergy < 2.0*GeV) && (G4UniformRand() >= supp[ieab]) )
235  {
236    // suppress high multiplicity events at low momentum
237    // only one pion will be produced
238     
239    nm = np = nz = 0;
240    if( targetParticle.GetDefinition() == aProton )
241    {
242      test = std::exp( std::min( expxu, std::max( expxl, -sqr(1.0+b[0])/(2.0*c*c) ) ) );
243      w0 = test;
244      wp = test*2.0;       
245      if( G4UniformRand() < w0/(w0+wp) )
246        nz = 1;
247      else
248        np = 1;
249    }
250    else  // target is a neutron
251    {
252      test = std::exp( std::min( expxu, std::max( expxl, -sqr(1.0+b[1])/(2.0*c*c) ) ) );
253      w0 = test;
254      wp = test;
255      test = std::exp( std::min( expxu, std::max( expxl, -sqr(-1.0+b[1])/(2.0*c*c) ) ) );
256      wm = test;
257      wt = w0+wp+wm;
258      wp += w0;
259      G4double ran = G4UniformRand();
260      if( ran < w0/wt )
261        nz = 1;
262      else if( ran < wp/wt )
263        np = 1;
264      else
265        nm = 1;
266    }
267  }
268  else
269  {
270      G4double n, anpn;
271      GetNormalizationConstant( availableEnergy, n, anpn );
272      G4double ran = G4UniformRand();
273      G4double dum, excs = 0.0;
274      if( targetParticle.GetDefinition() == aProton )
275      {
276        G4int counter = -1;
277        for( np=0; (np<numSec/3) && (ran>=excs); ++np )
278        {
279          for( nm=std::max(0,np-2); (nm<=np) && (ran>=excs); ++nm )
280          {
281            for( nz=0; (nz<numSec/3) && (ran>=excs); ++nz )
282            {
283              if( ++counter < numMul )
284              {
285                nt = np+nm+nz;
286                if( nt > 0 )
287                {
288                  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
289                  dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
290                  if( std::fabs(dum) < 1.0 )
291                  {
292                    if( test >= 1.0e-10 )excs += dum*test;
293                  }
294                  else
295                    excs += dum*test;
296                }
297              }
298            }
299          }
300        }
301        if( ran >= excs )return;  // 3 previous loops continued to the end
302        np--; nm--; nz--;
303      }
304      else  // target must be a neutron
305      {
306        G4int counter = -1;
307        for( np=0; (np<numSec/3) && (ran>=excs); ++np )
308        {
309          for( nm=std::max(0,np-1); (nm<=(np+1)) && (ran>=excs); ++nm )
310          {
311            for( nz=0; (nz<numSec/3) && (ran>=excs); ++nz )
312            {
313              if( ++counter < numMul )
314              {
315                nt = np+nm+nz;
316                if( (nt>=1) && (nt<=numSec) )
317                {
318                  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
319                  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
320                  if( std::fabs(dum) < 1.0 )
321                  {
322                    if( test >= 1.0e-10 )excs += dum*test;
323                  }
324                  else
325                    excs += dum*test;
326                }
327              }
328            }
329          }
330        }
331        if( ran >= excs )return;  // 3 previous loops continued to the end
332        np--; nm--; nz--;
333      }
334  }
335
336  if( targetParticle.GetDefinition() == aProton )
337  {
338    switch( np-nm )
339    {
340     case 1:
341       if( G4UniformRand() < 0.5 )
342       {
343         if( G4UniformRand() < 0.5 )
344           currentParticle.SetDefinitionAndUpdateE( aKaonZS );
345         else
346           currentParticle.SetDefinitionAndUpdateE( aKaonZL );
347         incidentHasChanged = true;
348       }
349       else
350       {
351         targetParticle.SetDefinitionAndUpdateE( aNeutron );
352         targetHasChanged = true;
353       }
354       break;
355     case 2:
356       if( G4UniformRand() < 0.5 )
357         currentParticle.SetDefinitionAndUpdateE( aKaonZS );
358       else
359         currentParticle.SetDefinitionAndUpdateE( aKaonZL );
360       incidentHasChanged = true;
361       targetParticle.SetDefinitionAndUpdateE( aNeutron );
362       incidentHasChanged = true;
363       targetHasChanged = true;
364       break;
365     default:
366       break;
367    }
368  }
369  else   // target is a neutron
370  {
371    switch( np-nm )
372    {
373     case 0:
374       if( G4UniformRand() < 0.25 )
375       {
376         if( G4UniformRand() < 0.5 )
377           currentParticle.SetDefinitionAndUpdateE( aKaonZS );
378         else
379           currentParticle.SetDefinitionAndUpdateE( aKaonZL );
380         targetParticle.SetDefinitionAndUpdateE( aProton );
381         incidentHasChanged = true;
382         targetHasChanged = true;
383       }
384       break;
385     case 1:
386       if( G4UniformRand() < 0.5 )
387         currentParticle.SetDefinitionAndUpdateE( aKaonZS );
388       else
389         currentParticle.SetDefinitionAndUpdateE( aKaonZL );
390       incidentHasChanged = true;
391       break;
392     default: // assumes nm = np+1 so charge is conserved
393       targetParticle.SetDefinitionAndUpdateE( aProton );
394       targetHasChanged = true;
395       break;
396    }
397  }
398  SetUpPions( np, nm, nz, vec, vecLen );
399  return;
400}
401
402 /* end of file */
403 
Note: See TracBrowser for help on using the repository browser.