source: trunk/source/processes/hadronic/models/low_energy/src/G4LEAntiKaonZeroInelastic.cc @ 1340

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

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