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

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

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