source: trunk/source/processes/hadronic/models/low_energy/src/G4LEPionMinusInelastic.cc @ 1358

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

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