source: trunk/source/processes/hadronic/models/low_energy/src/G4LEKaonMinusInelastic.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: 19.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: G4LEKaonMinusInelastic.cc,v 1.15 2006/06/29 20:44:59 gunter Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-01 $
29//
30 // Hadronic Process: Low Energy KaonMinus Inelastic Process
31 // J.L. Chuma, TRIUMF, 12-Feb-1997
32 // Last modified: 27-Mar-1997
33 // J.P.Wellisch 23-Apr-97: bug-hunting (missing initialization of np,nm,nz fixed)
34 // Modified by J.L.Chuma 30-Apr-97: added originalTarget for CalculateMomenta
35 
36#include "G4LEKaonMinusInelastic.hh"
37 #include "Randomize.hh"
38
39 G4HadFinalState *
40  G4LEKaonMinusInelastic::ApplyYourself( const G4HadProjectile &aTrack,
41                                         G4Nucleus &targetNucleus )
42  {
43    const G4HadProjectile *originalIncident = &aTrack;
44    if (originalIncident->GetKineticEnergy()<= 0.1*MeV) 
45    {
46      theParticleChange.SetStatusChange(isAlive);
47      theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
48      theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 
49      return &theParticleChange;     
50    }
51
52    // create the target particle
53   
54    G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
55    G4ReactionProduct targetParticle( originalTarget->GetDefinition() );
56   
57    if( verboseLevel > 1 )
58    {
59      const G4Material *targetMaterial = aTrack.GetMaterial();
60      G4cout << "G4LEKaonMinusInelastic::ApplyYourself called" << G4endl;
61      G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy() << "MeV, ";
62      G4cout << "target material = " << targetMaterial->GetName() << ", ";
63      G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
64           << G4endl;
65    }
66    G4ReactionProduct currentParticle( const_cast<G4ParticleDefinition *>(originalIncident->GetDefinition()) );
67    currentParticle.SetMomentum( originalIncident->Get4Momentum().vect() );
68    currentParticle.SetKineticEnergy( originalIncident->GetKineticEnergy() );
69   
70    // Fermi motion and evaporation
71    // As of Geant3, the Fermi energy calculation had not been Done
72   
73    G4double ek = originalIncident->GetKineticEnergy();
74    G4double amas = originalIncident->GetDefinition()->GetPDGMass();
75   
76    G4double tkin = targetNucleus.Cinema( ek );
77    ek += tkin;
78    currentParticle.SetKineticEnergy( ek );
79    G4double et = ek + amas;
80    G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
81    G4double pp = currentParticle.GetMomentum().mag();
82    if( pp > 0.0 )
83    {
84      G4ThreeVector momentum = currentParticle.GetMomentum();
85      currentParticle.SetMomentum( momentum * (p/pp) );
86    }
87   
88    // calculate black track energies
89   
90    tkin = targetNucleus.EvaporationEffects( ek );
91    ek -= tkin;
92    currentParticle.SetKineticEnergy( ek );
93    et = ek + amas;
94    p = std::sqrt( std::abs((et-amas)*(et+amas)) );
95    pp = currentParticle.GetMomentum().mag();
96    if( pp > 0.0 )
97    {
98      G4ThreeVector momentum = currentParticle.GetMomentum();
99      currentParticle.SetMomentum( momentum * (p/pp) );
100    }
101
102    G4ReactionProduct modifiedOriginal = currentParticle;
103
104    currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
105    targetParticle.SetSide( -1 );  // target always goes in backward hemisphere
106    G4bool incidentHasChanged = false;
107    G4bool targetHasChanged = false;
108    G4bool quasiElastic = false;
109    G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec;  // vec will contain the secondary particles
110    G4int vecLen = 0;
111    vec.Initialize( 0 );
112   
113    const G4double cutOff = 0.1*MeV;
114    if( currentParticle.GetKineticEnergy() > cutOff )
115      Cascade( vec, vecLen,
116               originalIncident, currentParticle, targetParticle,
117               incidentHasChanged, targetHasChanged, quasiElastic );
118   
119    CalculateMomenta( vec, vecLen,
120                      originalIncident, originalTarget, modifiedOriginal,
121                      targetNucleus, currentParticle, targetParticle,
122                      incidentHasChanged, targetHasChanged, quasiElastic );
123   
124    SetUpChange( vec, vecLen,
125                 currentParticle, targetParticle,
126                 incidentHasChanged );
127   
128    delete originalTarget;
129    return &theParticleChange;   
130  }
131 
132 void
133  G4LEKaonMinusInelastic::Cascade(
134   G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
135   G4int& vecLen,
136   const G4HadProjectile *originalIncident,
137   G4ReactionProduct &currentParticle,
138   G4ReactionProduct &targetParticle,
139   G4bool &incidentHasChanged,
140   G4bool &targetHasChanged,
141   G4bool &quasiElastic )
142  {
143    // derived from original FORTRAN code CASKM by H. Fesefeldt (13-Sep-1987)
144    //
145    // K- undergoes interaction with nucleon within a nucleus.  Check if it is
146    // energetically possible to produce pions/kaons.  In not, assume nuclear excitation
147    // occurs and input particle is degraded in energy. No other particles are produced.
148    // If reaction is possible, find the correct number of pions/protons/neutrons
149    // produced using an interpolation to multiplicity data.  Replace some pions or
150    // protons/neutrons by kaons or strange baryons according to the average
151    // multiplicity per Inelastic reaction.
152    //
153    const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass();
154    const G4double etOriginal = originalIncident->GetTotalEnergy();
155    const G4double pOriginal = originalIncident->GetTotalMomentum();
156    const G4double targetMass = targetParticle.GetMass();
157    G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
158                                        targetMass*targetMass +
159                                        2.0*targetMass*etOriginal );
160    G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
161   
162    static G4bool first = true;
163    const G4int numMul = 1200;
164    const G4int numSec = 60;
165    static G4double protmul[numMul], protnorm[numSec]; // proton constants
166    static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
167    // np = number of pi+, nm = number of pi-, nz = number of pi0
168    G4int nt(0), np(0), nm(0), nz(0);
169    const G4double c = 1.25;   
170    const G4double b[] = { 0.70, 0.70 };
171    if( first )       // compute normalization constants, this will only be Done once
172    {
173      first = false;
174      G4int i;
175      for( i=0; i<numMul; ++i )protmul[i] = 0.0;
176      for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
177      G4int counter = -1;
178      for( np=0; np<(numSec/3); ++np )
179      {
180        for( nm=std::max(0,np-1); nm<=(np+1); ++nm )
181        {
182          for( nz=0; nz<numSec/3; ++nz )
183          {
184            if( ++counter < numMul )
185            {
186              nt = np+nm+nz;
187              if( (nt>0) && (nt<=numSec) )
188              {
189                protmul[counter] = Pmltpc(np,nm,nz,nt,b[0],c);
190                protnorm[nt-1] += protmul[counter];
191              }
192            }
193          }
194        }
195      }
196      for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
197      for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
198      counter = -1;
199      for( np=0; np<numSec/3; ++np )
200      {
201        for( nm=np; nm<=(np+2); ++nm )
202        {
203          for( nz=0; nz<numSec/3; ++nz )
204          {
205            if( ++counter < numMul )
206            {
207              nt = np+nm+nz;
208              if( (nt>0) && (nt<=numSec) )
209              {
210                neutmul[counter] = Pmltpc(np,nm,nz,nt,b[1],c);
211                neutnorm[nt-1] += neutmul[counter];
212              }
213            }
214          }
215        }
216      }
217      for( i=0; i<numSec; ++i )
218      {
219        if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
220        if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
221      }
222    }   // end of initialization
223   
224    const G4double expxu = 82.;           // upper bound for arg. of exp
225    const G4double expxl = -expxu;        // lower bound for arg. of exp
226    G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
227    G4ParticleDefinition *aKaonZS = G4KaonZeroShort::KaonZeroShort();
228    G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong();
229    G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
230    G4ParticleDefinition *aProton = G4Proton::Proton();
231    G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
232    G4ParticleDefinition *aPiMinus = G4PionMinus::PionMinus();
233    G4ParticleDefinition *aPiZero = G4PionZero::PionZero();
234    G4ParticleDefinition *aLambda = G4Lambda::Lambda();
235    G4ParticleDefinition *aSigmaPlus = G4SigmaPlus::SigmaPlus();
236    G4ParticleDefinition *aSigmaMinus = G4SigmaMinus::SigmaMinus();
237    G4ParticleDefinition *aSigmaZero = G4SigmaZero::SigmaZero();
238    const G4double cech[] = {1.,1.,1.,0.70,0.60,0.55,0.35,0.25,0.18,0.15};
239    G4int iplab = G4int(std::min( 9.0, pOriginal/GeV*5.0 ));
240    if( (pOriginal <= 2.0*GeV) && (G4UniformRand() < cech[iplab]) )
241    {
242      np = nm = nz = nt = 0;
243      iplab = G4int(std::min( 19.0, pOriginal/GeV*10.0 ));
244      const G4double cnk0[] = {0.17,0.18,0.17,0.24,0.26,0.20,0.22,0.21,0.34,0.45,
245                               0.58,0.55,0.36,0.29,0.29,0.32,0.32,0.33,0.33,0.33};
246      if( G4UniformRand() <= cnk0[iplab] )
247      {
248        quasiElastic = true;
249        if( targetParticle.GetDefinition() == aProton )
250        {
251          currentParticle.SetDefinitionAndUpdateE( aKaonZL );
252          incidentHasChanged = true;
253          targetParticle.SetDefinitionAndUpdateE( aNeutron );
254          targetHasChanged = true;
255        }
256      }
257      else  // random number > cnk0[iplab]
258      {
259        G4double ran = G4UniformRand();
260        if( ran < 0.25 )         // k- p --> pi- s+
261        {
262          if( targetParticle.GetDefinition() == aProton )
263          {
264            currentParticle.SetDefinitionAndUpdateE( aPiMinus );
265            targetParticle.SetDefinitionAndUpdateE( aSigmaPlus );
266            incidentHasChanged = true;
267            targetHasChanged = true;
268          }
269        }
270        else if( ran < 0.50 )  // k- p --> pi0 s0  or  k- n --> pi- s0
271        {
272          if( targetParticle.GetDefinition() == aNeutron )
273            currentParticle.SetDefinitionAndUpdateE( aPiMinus );
274          else
275            currentParticle.SetDefinitionAndUpdateE( aPiZero );
276          targetParticle.SetDefinitionAndUpdateE( aSigmaZero );
277          incidentHasChanged = true;
278          targetHasChanged = true;
279        }
280        else if( ran < 0.75 )  // k- p --> pi+ s-  or  k- n --> pi0 s-
281        {
282          if( targetParticle.GetDefinition() == aNeutron )
283            currentParticle.SetDefinitionAndUpdateE( aPiZero );
284          else
285            currentParticle.SetDefinitionAndUpdateE( aPiPlus );
286          targetParticle.SetDefinitionAndUpdateE( aSigmaMinus );
287          incidentHasChanged = true;
288          targetHasChanged = true;
289        }
290        else                   // k- p --> pi0 L  or  k- n --> pi- L
291        {
292          if( targetParticle.GetDefinition() == aNeutron )
293            currentParticle.SetDefinitionAndUpdateE( aPiMinus );
294          else
295            currentParticle.SetDefinitionAndUpdateE( aPiZero );
296          targetParticle.SetDefinitionAndUpdateE( aLambda );
297          incidentHasChanged = true;
298          targetHasChanged = true;
299        }
300      }
301    }
302    else  // (pOriginal > 2.0*GeV) || (random number >= cech[iplab])
303    {
304      if( availableEnergy < aPiPlus->GetPDGMass() )
305      {               // not energetically possible to produce pion(s)
306        quasiElastic = true;
307        return;
308      }
309      G4double n, anpn;
310      GetNormalizationConstant( availableEnergy, n, anpn );
311      G4double ran = G4UniformRand();
312      G4double dum, test, excs = 0.0;
313      if( targetParticle.GetDefinition() == aProton )
314      {
315        G4int counter = -1;
316        for( np=0; np<numSec/3 && ran>=excs; ++np )
317        {
318          for( nm=std::max(0,np-1); nm<=(np+1) && ran>=excs; ++nm )
319          {
320            for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
321            {
322              if( ++counter < numMul )
323              {
324                nt = np+nm+nz;
325                if( nt > 0 )
326                {
327                  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
328                  dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
329                  if( std::fabs(dum) < 1.0 )
330                  {
331                    if( test >= 1.0e-10 )excs += dum*test;
332                  }
333                  else
334                    excs += dum*test;
335                }
336              }
337            }
338          }
339        }
340        if( ran >= excs )  // 3 previous loops continued to the end
341        {
342          quasiElastic = true;
343          return;
344        }
345        np--; nm--; nz--;
346        if( np == nm )
347        {
348          if( G4UniformRand() >= 0.75 )
349          {
350            currentParticle.SetDefinitionAndUpdateE( aKaonZL );
351            targetParticle.SetDefinitionAndUpdateE( aNeutron );
352            incidentHasChanged = true;
353            targetHasChanged = true;
354          }
355        }
356        else if( np == nm+1 )
357        {
358          targetParticle.SetDefinitionAndUpdateE( aNeutron );
359          targetHasChanged = true;
360        }
361        else
362        {
363          currentParticle.SetDefinitionAndUpdateE( aKaonZL );
364          incidentHasChanged = true;
365        }
366      }
367      else   // target must be a neutron
368      {
369        G4int counter = -1;
370        for( np=0; np<numSec/3 && ran>=excs; ++np )
371        {
372          for( nm=np; nm<=(np+2) && ran>=excs; ++nm )
373          {
374            for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
375            {
376              if( ++counter < numMul )
377              {
378                nt = np+nm+nz;
379                if( (nt>=1) && (nt<=numSec) )
380                {
381                  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
382                  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
383                  if( std::fabs(dum) < 1.0 )
384                  {
385                    if( test >= 1.0e-10 )excs += dum*test;
386                  }
387                  else
388                    excs += dum*test;
389                }
390              }
391            }
392          }
393        }
394        if( ran >= excs )  // 3 previous loops continued to the end
395        {
396          quasiElastic = true;
397          return;
398        }
399        np--; nm--; nz--;
400        if( np == nm-1 )
401        {
402          if( G4UniformRand() < 0.5 )
403          {
404            targetParticle.SetDefinitionAndUpdateE( aProton );
405            targetHasChanged = true;
406          }
407          else
408          {
409            currentParticle.SetDefinitionAndUpdateE( aKaonZL );
410            incidentHasChanged = true;
411          }
412        }
413        else if( np != nm )
414        {
415          currentParticle.SetDefinitionAndUpdateE( aKaonZL );
416          incidentHasChanged = true;
417        }
418      }
419      if( G4UniformRand() >= 0.5 )
420      {
421        if( (currentParticle.GetDefinition() == aKaonMinus &&
422             targetParticle.GetDefinition() == aNeutron )     ||
423            (currentParticle.GetDefinition() == aKaonZL &&
424             targetParticle.GetDefinition() == aProton ) )
425        {
426          ran = G4UniformRand();
427          if( ran < 0.68 )
428          {
429            if( targetParticle.GetDefinition() == aProton )
430            {
431              currentParticle.SetDefinitionAndUpdateE( aPiPlus );
432              targetParticle.SetDefinitionAndUpdateE( aLambda );
433              incidentHasChanged = true;
434              targetHasChanged = true;
435            }
436            else
437            {
438              currentParticle.SetDefinitionAndUpdateE( aPiMinus );
439              targetParticle.SetDefinitionAndUpdateE( aLambda );
440              incidentHasChanged = true;
441              targetHasChanged = true;
442            }
443          }
444          else if( ran < 0.84 )
445          {
446            if( targetParticle.GetDefinition() == aProton )
447            {
448              currentParticle.SetDefinitionAndUpdateE( aPiZero );
449              targetParticle.SetDefinitionAndUpdateE( aSigmaPlus );
450              incidentHasChanged = true;
451              targetHasChanged = true;
452            }
453            else
454            {
455              currentParticle.SetDefinitionAndUpdateE( aPiMinus );
456              targetParticle.SetDefinitionAndUpdateE( aSigmaZero );
457              incidentHasChanged = true;
458              targetHasChanged = true;
459            }
460          }
461          else
462          {
463            if( targetParticle.GetDefinition() == aProton )
464            {
465              currentParticle.SetDefinitionAndUpdateE( aPiPlus );
466              targetParticle.SetDefinitionAndUpdateE( aSigmaZero );
467              incidentHasChanged = true;
468              targetHasChanged = true;
469            }
470            else
471            {
472              currentParticle.SetDefinitionAndUpdateE( aPiZero );
473              targetParticle.SetDefinitionAndUpdateE( aSigmaMinus );
474              incidentHasChanged = true;
475              targetHasChanged = true;
476            }
477          }
478        }
479        else  // ( current != aKaonMinus || target != aNeutron ) &&
480              // ( current != aKaonZL    || target != aProton  )
481        {
482          ran = G4UniformRand();
483          if( ran < 0.67 )
484          {
485            currentParticle.SetDefinitionAndUpdateE( aPiZero );
486            targetParticle.SetDefinitionAndUpdateE( aLambda );
487            incidentHasChanged = true;
488            targetHasChanged = true;
489          }
490          else if( ran < 0.78 )
491          {
492            currentParticle.SetDefinitionAndUpdateE( aPiMinus );
493            targetParticle.SetDefinitionAndUpdateE( aSigmaPlus );
494            incidentHasChanged = true;
495            targetHasChanged = true;
496          }
497          else if( ran < 0.89 )
498          {
499            currentParticle.SetDefinitionAndUpdateE( aPiZero );
500            targetParticle.SetDefinitionAndUpdateE( aSigmaZero );
501            incidentHasChanged = true;
502            targetHasChanged = true;
503          }
504          else
505          {
506            currentParticle.SetDefinitionAndUpdateE( aPiPlus );
507            targetParticle.SetDefinitionAndUpdateE( aSigmaMinus );
508            incidentHasChanged = true;
509            targetHasChanged = true;
510          }
511        }
512      }
513    }
514    if( currentParticle.GetDefinition() == aKaonZL )
515    {
516      if( G4UniformRand() >= 0.5 )
517      {
518        currentParticle.SetDefinitionAndUpdateE( aKaonZS );
519        incidentHasChanged = true;
520      }
521    }
522    if( targetParticle.GetDefinition() == aKaonZL )
523    {
524      if( G4UniformRand() >= 0.5 )
525      {
526        targetParticle.SetDefinitionAndUpdateE( aKaonZS );
527        targetHasChanged = true;
528      }
529    }
530    SetUpPions( np, nm, nz, vec, vecLen );
531    return;
532  }
533
534 /* end of file */
535 
536
537
Note: See TracBrowser for help on using the repository browser.