source: trunk/source/processes/hadronic/models/low_energy/src/G4LEAntiNeutronInelastic.cc @ 819

Last change on this file since 819 was 819, checked in by garnier, 16 years ago

import all except CVS
File size: 18.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: G4LEAntiNeutronInelastic.cc,v 1.14 2006/06/29 20:44:43 gunter Exp $
28// GEANT4 tag $Name:  $
29//
30 // Hadronic Process: AntiNeutron Inelastic Process
31 // J.L. Chuma, TRIUMF, 18-Feb-1997
32 // Last modified: 27-Mar-1997
33 // J.P.Wellisch: 23-Apr-97: Added theNucleus.SetParameters call
34 // J.P. Wellisch: 23-Apr-97: nm = np+1; in line 392
35 // Modified by J.L.Chuma 30-Apr-97: added originalTarget for CalculateMomenta
36 
37#include "G4LEAntiNeutronInelastic.hh"
38#include "Randomize.hh"
39
40 G4HadFinalState *
41  G4LEAntiNeutronInelastic::ApplyYourself( const G4HadProjectile &aTrack,
42                                           G4Nucleus &targetNucleus )
43  { 
44    const G4HadProjectile *originalIncident = &aTrack;
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 << "G4LEAntiNeutronInelastic::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   
94    G4ReactionProduct currentParticle = modifiedOriginal;
95    G4ReactionProduct targetParticle;
96    targetParticle = *originalTarget;
97    currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
98    targetParticle.SetSide( -1 );  // target always goes in backward hemisphere
99    G4bool incidentHasChanged = false;
100    G4bool targetHasChanged = false;
101    G4bool quasiElastic = false;
102    G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec;  // vec will contain the secondary particles
103    G4int vecLen = 0;
104    vec.Initialize( 0 );
105   
106    const G4double cutOff = 0.1*MeV;
107    const G4double anni = std::min( 1.3*currentParticle.GetTotalMomentum()/GeV, 0.4 );
108   
109    if( (currentParticle.GetKineticEnergy()/MeV > cutOff) ||
110        (G4UniformRand() > anni) )
111      Cascade( vec, vecLen,
112               originalIncident, currentParticle, targetParticle,
113               incidentHasChanged, targetHasChanged, quasiElastic );
114    else
115      quasiElastic = true;
116   
117    CalculateMomenta( vec, vecLen,
118                      originalIncident, originalTarget, modifiedOriginal,
119                      targetNucleus, currentParticle, targetParticle,
120                      incidentHasChanged, targetHasChanged, quasiElastic );
121   
122    SetUpChange( vec, vecLen,
123                 currentParticle, targetParticle,
124                 incidentHasChanged );
125   
126    delete originalTarget;
127    return &theParticleChange;
128  }
129 
130 void
131  G4LEAntiNeutronInelastic::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 CASNB by H. Fesefeldt (13-Sep-1987)
142    //
143    // AntiNeutron 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->GetTotalEnergy()/MeV;
153    const G4double pOriginal = originalIncident->GetTotalMomentum()/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 numMulA = 400;
163    const G4int numSec = 60;
164    static G4double protmul[numMul], protnorm[numSec]; // proton constants
165    static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
166    static G4double protmulA[numMulA], protnormA[numSec]; // proton constants
167    static G4double neutmulA[numMulA], neutnormA[numSec]; // neutron constants
168    // np = number of pi+, nm = number of pi-, nz = number of pi0
169    G4int counter, nt=0, np=0, nm=0, nz=0;
170    G4double test;
171    const G4double c = 1.25;   
172    const G4double b[] = { 0.70, 0.70 };
173    if( first )       // compute normalization constants, this will only be Done once
174    {
175      first = false;
176      G4int i;
177      for( i=0; i<numMul; ++i )protmul[i] = 0.0;
178      for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
179      counter = -1;
180      for( np=0; np<(numSec/3); ++np )
181      {
182        for( nm=std::max(0,np-2); nm<=np; ++nm )
183        {
184          for( nz=0; nz<numSec/3; ++nz )
185          {
186            if( ++counter < numMul )
187            {
188              nt = np+nm+nz;
189              if( nt>0 && nt<=numSec )
190              {
191                protmul[counter] = Pmltpc(np,nm,nz,nt,b[0],c);
192                protnorm[nt-1] += protmul[counter];
193              }
194            }
195          }
196        }
197      }
198      for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
199      for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
200      counter = -1;
201      for( np=0; np<numSec/3; ++np )
202      {
203        for( nm=std::max(0,np-1); nm<=(np+1); ++nm )
204        {
205          for( nz=0; nz<numSec/3; ++nz )
206          {
207            if( ++counter < numMul )
208            {
209              nt = np+nm+nz;
210              if( (nt>0) && (nt<=numSec) )
211              {
212                neutmul[counter] = Pmltpc(np,nm,nz,nt,b[1],c);
213                neutnorm[nt-1] += neutmul[counter];
214              }
215            }
216          }
217        }
218      }
219      for( i=0; i<numSec; ++i )
220      {
221        if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
222        if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
223      }
224      //
225      // do the same for annihilation channels
226      //
227      for( i=0; i<numMulA; ++i )protmulA[i] = 0.0;
228      for( i=0; i<numSec; ++i )protnormA[i] = 0.0;
229      counter = -1;
230      for( np=1; np<(numSec/3); ++np )
231      {
232        nm = np-1;
233        for( nz=0; nz<numSec/3; ++nz )
234        {
235          if( ++counter < numMulA )
236          {
237            nt = np+nm+nz;
238            if( nt>1 && nt<=numSec )
239            {
240              protmulA[counter] = Pmltpc(np,nm,nz,nt,b[0],c);
241              protnormA[nt-1] += protmulA[counter];
242            }
243          }
244        }
245      }
246      for( i=0; i<numMulA; ++i )neutmulA[i] = 0.0;
247      for( i=0; i<numSec; ++i )neutnormA[i] = 0.0;
248      counter = -1;
249      for( np=0; np<numSec/3; ++np )
250      {
251        nm = np;
252        for( nz=0; nz<numSec/3; ++nz )
253        {
254          if( ++counter < numMulA )
255          {
256            nt = np+nm+nz;
257            if( nt>1 && nt<=numSec )
258            {
259              neutmulA[counter] = Pmltpc(np,nm,nz,nt,b[1],c);
260              neutnormA[nt-1] += neutmulA[counter];
261            }
262          }
263        }
264      }
265      for( i=0; i<numSec; ++i )
266      {
267        if( protnormA[i] > 0.0 )protnormA[i] = 1.0/protnormA[i];
268        if( neutnormA[i] > 0.0 )neutnormA[i] = 1.0/neutnormA[i];
269      }
270    }   // end of initialization
271    const G4double expxu = 82.;           // upper bound for arg. of exp
272    const G4double expxl = -expxu;        // lower bound for arg. of exp
273    G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
274    G4ParticleDefinition *aProton = G4Proton::Proton();
275    G4ParticleDefinition *anAntiProton = G4AntiProton::AntiProton();
276    G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
277   
278    // energetically possible to produce pion(s)  -->  inelastic scattering
279    //                                   otherwise quasi-elastic scattering
280   
281    const G4double anhl[] = {1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,0.97,0.88,
282                             0.85,0.81,0.75,0.64,0.64,0.55,0.55,0.45,0.47,0.40,
283                             0.39,0.36,0.33,0.10,0.01};
284    G4int iplab = G4int( pOriginal/GeV*10.0 );
285    if( iplab >  9 )iplab = G4int( (pOriginal/GeV- 1.0)*5.0 ) + 10;
286    if( iplab > 14 )iplab = G4int(  pOriginal/GeV- 2.0 ) + 15;
287    if( iplab > 22 )iplab = G4int( (pOriginal/GeV-10.0)/10.0 ) + 23;
288    if( iplab > 24 )iplab = 24;
289    if( G4UniformRand() > anhl[iplab] )
290    {
291      if( availableEnergy <= aPiPlus->GetPDGMass()/MeV )
292      {
293        quasiElastic = true;
294        return;
295      }
296      G4int ieab = static_cast<G4int>(availableEnergy*5.0/GeV);
297      const G4double supp[] = {0.,0.4,0.55,0.65,0.75,0.82,0.86,0.90,0.94,0.98};
298      G4double w0, wp, wt, wm;
299      if( (availableEnergy < 2.0*GeV) && (G4UniformRand() >= supp[ieab]) )
300      {
301        // suppress high multiplicity events at low momentum
302        // only one pion will be produced
303        //
304        np = nm = nz = 0;
305        if( targetParticle.GetDefinition() == aProton )
306        {
307          test = std::exp( std::min( expxu, std::max( expxl, -(1.0+b[0])*(1.0+b[0])/(2.0*c*c) ) ) );
308          w0 = test;
309          wp = test;
310          if( G4UniformRand() < w0/(w0+wp) )
311            nz = 1;
312          else
313            np = 1;
314        }
315        else  // target is a neutron
316        {
317          test = std::exp( std::min( expxu, std::max( expxl, -(1.0+b[1])*(1.0+b[1])/(2.0*c*c) ) ) );
318          w0 = test;
319          wp = test;
320          test = std::exp( std::min( expxu, std::max( expxl, -(-1.0+b[1])*(-1.0+b[1])/(2.0*c*c) ) ) );
321          wm = test;
322          wt = w0+wp+wm;
323          wp += w0;
324          G4double ran = G4UniformRand();
325          if( ran < w0/wt )
326            nz = 1;
327          else if( ran < wp/wt )
328            np = 1;
329          else
330            nm = 1;
331        }
332      }
333      else   // (availableEnergy >= 2.0*GeV) || (random number < supp[ieab])
334      {
335        G4double n, anpn;
336        GetNormalizationConstant( availableEnergy, n, anpn );
337        G4double ran = G4UniformRand();
338        G4double dum, excs = 0.0;
339        if( targetParticle.GetDefinition() == aProton )
340        {
341          counter = -1;
342          for( np=0; np<numSec/3 && ran>=excs; ++np )
343          {
344            for( nm=std::max(0,np-2); nm<=np && ran>=excs; ++nm )
345            {
346              for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
347              {
348                if( ++counter < numMul )
349                {
350                  nt = np+nm+nz;
351                  if( nt > 0 )
352                  {
353                    test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
354                    dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
355                    if( std::fabs(dum) < 1.0 )
356                    {
357                      if( test >= 1.0e-10 )excs += dum*test;
358                    }
359                    else
360                      excs += dum*test;
361                  }
362                }
363              }
364            }
365          }
366          if( ran >= excs )  // 3 previous loops continued to the end
367          {
368            quasiElastic = true;
369            return;
370          }
371          np--; nm--; nz--;
372        }
373        else  // target must be a neutron
374        {
375          counter = -1;
376          for( np=0; np<numSec/3 && ran>=excs; ++np )
377          {
378            for( nm=std::max(0,np-1); nm<=(np+1) && ran>=excs; ++nm )
379            {
380              for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
381              {
382                if( ++counter < numMul )
383                {
384                  nt = np+nm+nz;
385                  if( (nt>=1) && (nt<=numSec) )
386                  {
387                    test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
388                    dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
389                    if( std::fabs(dum) < 1.0 )
390                    {
391                      if( test >= 1.0e-10 )excs += dum*test;
392                    }
393                    else
394                      excs += dum*test;
395                  }
396                }
397              }
398            }
399          }
400          if( ran >= excs )  // 3 previous loops continued to the end
401          {
402            quasiElastic = true;
403            return;
404          }
405          np--; nm--; nz--;
406        }
407      }
408      if( targetParticle.GetDefinition() == aProton )
409      {
410        switch( np-nm )
411        {
412         case 1:
413           if( G4UniformRand() < 0.5 )
414           {
415             currentParticle.SetDefinitionAndUpdateE( anAntiProton );
416             incidentHasChanged = true;
417           }
418           else
419           {
420             targetParticle.SetDefinitionAndUpdateE( aNeutron );
421             targetHasChanged = true;
422           }
423           break;
424         case 2:
425           currentParticle.SetDefinitionAndUpdateE( anAntiProton );
426           targetParticle.SetDefinitionAndUpdateE( aNeutron );
427           incidentHasChanged = true;
428           targetHasChanged = true;
429           break;
430         default:
431           break;
432        }
433      }
434      else  // target must be a neutron
435      {
436        switch( np-nm )
437        {
438         case 0:
439           if( G4UniformRand() < 0.33 )
440           {
441             currentParticle.SetDefinitionAndUpdateE( anAntiProton );
442             targetParticle.SetDefinitionAndUpdateE( aProton );
443             incidentHasChanged = true;
444             targetHasChanged = true;
445           }
446           break;
447         case 1:
448           currentParticle.SetDefinitionAndUpdateE( anAntiProton );
449           incidentHasChanged = true;
450           break;
451         default:
452           targetParticle.SetDefinitionAndUpdateE( aProton );
453           targetHasChanged = true;
454           break;
455        }
456      }
457    }
458    else   // random number <= anhl[iplab]
459    {
460      if( centerofmassEnergy <= 2*aPiPlus->GetPDGMass()/MeV )
461      {
462        quasiElastic = true;
463        return;
464      }
465      //
466      // annihilation channels
467      //
468      G4double n, anpn;
469      GetNormalizationConstant( -centerofmassEnergy, n, anpn );
470      G4double ran = G4UniformRand();
471      G4double dum, excs = 0.0;
472      if( targetParticle.GetDefinition() == aProton )
473      {
474        counter = -1;
475        for( np=1; (np<numSec/3) && (ran>=excs); ++np )
476        {
477          nm = np-1;
478          for( nz=0; (nz<numSec/3) && (ran>=excs); ++nz )
479          {
480            if( ++counter < numMulA )
481            {
482              nt = np+nm+nz;
483              if( nt>1 && nt<=numSec )
484              {
485                test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
486                dum = (pi/anpn)*nt*protmulA[counter]*protnormA[nt-1]/(2.0*n*n);
487                if( std::fabs(dum) < 1.0 )
488                {
489                  if( test >= 1.0e-10 )excs += dum*test;
490                }
491                else
492                  excs += dum*test;
493              }
494            }
495          }
496        }
497      }
498      else  // target must be a neutron
499      {
500        counter = -1;
501        for( np=0; (np<numSec/3) && (ran>=excs); ++np )
502        {
503          nm = np;
504          for( nz=0; (nz<numSec/3) && (ran>=excs); ++nz )
505          {
506            if( ++counter < numMulA )
507            {
508              nt = np+nm+nz;
509              if( (nt>1) && (nt<=numSec) )
510              {
511                test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
512                dum = (pi/anpn)*nt*neutmulA[counter]*neutnormA[nt-1]/(2.0*n*n);
513                if( std::fabs(dum) < 1.0 )
514                {
515                  if( test >= 1.0e-10 )excs += dum*test;
516                }
517                else
518                  excs += dum*test;
519              }
520            }
521          }
522        }
523      }
524      if( ran >= excs )  // 3 previous loops continued to the end
525      {
526        quasiElastic = true;
527        return;
528      }
529      np--; nz--;
530      currentParticle.SetMass( 0.0 );
531      targetParticle.SetMass( 0.0 );
532    }
533    while(np+nm+nz<3) nz++;
534    SetUpPions( np, nm, nz, vec, vecLen );
535    return;
536  }
537
538 /* end of file */
539 
Note: See TracBrowser for help on using the repository browser.