source: trunk/source/processes/hadronic/models/rpg/src/G4RPGXiZeroInelastic.cc @ 1196

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

import all except CVS
File size: 13.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// $Id: G4RPGXiZeroInelastic.cc,v 1.1 2007/07/18 21:04:21 dennis Exp $
27//
28 
29#include "G4RPGXiZeroInelastic.hh"
30#include "Randomize.hh"
31 
32 G4HadFinalState *
33  G4RPGXiZeroInelastic::ApplyYourself( const G4HadProjectile &aTrack,
34                                      G4Nucleus &targetNucleus )
35  {   
36    const G4HadProjectile *originalIncident = &aTrack;
37    if (originalIncident->GetKineticEnergy()<= 0.1*MeV) 
38    {
39      theParticleChange.SetStatusChange(isAlive);
40      theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
41      theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 
42      return &theParticleChange;     
43    }
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 << "G4RPGXiZeroInelastic::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    CalculateMomenta( vec, vecLen,
111                      originalIncident, originalTarget, modifiedOriginal,
112                      targetNucleus, currentParticle, targetParticle,
113                      incidentHasChanged, targetHasChanged, quasiElastic );
114   
115    SetUpChange( vec, vecLen,
116                 currentParticle, targetParticle,
117                 incidentHasChanged );
118   
119    delete originalTarget;
120    return &theParticleChange;
121  }
122 
123void G4RPGXiZeroInelastic::Cascade(
124   G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
125   G4int& vecLen,
126   const G4HadProjectile *originalIncident,
127   G4ReactionProduct &currentParticle,
128   G4ReactionProduct &targetParticle,
129   G4bool &incidentHasChanged,
130   G4bool &targetHasChanged,
131   G4bool &quasiElastic )
132{
133  // Derived from H. Fesefeldt's original FORTRAN code CASX0
134  //
135  // XiZero undergoes interaction with nucleon within a nucleus.  Check if it is
136  // energetically possible to produce pions/kaons.  In not, assume nuclear excitation
137  // occurs and input particle is degraded in energy. No other particles are produced.
138  // If reaction is possible, find the correct number of pions/protons/neutrons
139  // produced using an interpolation to multiplicity data.  Replace some pions or
140  // protons/neutrons by kaons or strange baryons according to the average
141  // multiplicity per inelastic reaction.
142
143  const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass()/MeV;
144  const G4double etOriginal = originalIncident->GetTotalEnergy()/MeV;
145  const G4double targetMass = targetParticle.GetMass()/MeV;
146  G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
147                                      targetMass*targetMass +
148                                      2.0*targetMass*etOriginal );
149  G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
150  if( availableEnergy <= G4PionPlus::PionPlus()->GetPDGMass()/MeV )
151  {
152    quasiElastic = true;
153    return;
154  }
155    static G4bool first = true;
156    const G4int numMul = 1200;
157    const G4int numSec = 60;
158    static G4double protmul[numMul], protnorm[numSec]; // proton constants
159    static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
160    // np = number of pi+, nm = number of pi-, nz = number of pi0
161    G4int counter, nt=0, np=0, nm=0, nz=0;
162    G4double test;
163    const G4double c = 1.25;   
164    const G4double b[] = { 0.7, 0.7 };
165    if( first )       // compute normalization constants, this will only be Done once
166    {
167      first = false;
168      G4int i;
169      for( i=0; i<numMul; ++i )protmul[i] = 0.0;
170      for( i=0; i<numSec; ++i )protnorm[i] = 0.0;
171      counter = -1;
172      for( np=0; np<(numSec/3); ++np )
173      {
174        for( nm=std::max(0,np-2); nm<=(np+1); ++nm )
175        {
176          for( nz=0; nz<numSec/3; ++nz )
177          {
178            if( ++counter < numMul )
179            {
180              nt = np+nm+nz;
181              if( nt>0 && nt<=numSec )
182              {
183                protmul[counter] = Pmltpc(np,nm,nz,nt,b[0],c);
184                protnorm[nt-1] += protmul[counter];
185              }
186            }
187          }
188        }
189      }
190      for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
191      for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
192      counter = -1;
193      for( np=0; np<numSec/3; ++np )
194      {
195        for( nm=std::max(0,np-1); nm<=(np+2); ++nm )
196        {
197          for( nz=0; nz<numSec/3; ++nz )
198          {
199            if( ++counter < numMul )
200            {
201              nt = np+nm+nz;
202              if( nt>0 && nt<=numSec )
203              {
204                neutmul[counter] = Pmltpc(np,nm,nz,nt,b[1],c);
205                neutnorm[nt-1] += neutmul[counter];
206              }
207            }
208          }
209        }
210      }
211      for( i=0; i<numSec; ++i )
212      {
213        if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
214        if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
215      }
216    }   // end of initialization
217   
218    const G4double expxu = 82.;           // upper bound for arg. of exp
219    const G4double expxl = -expxu;        // lower bound for arg. of exp
220    G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
221    G4ParticleDefinition *aProton = G4Proton::Proton();
222    G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
223    G4ParticleDefinition *aSigmaPlus = G4SigmaPlus::SigmaPlus();
224    G4ParticleDefinition *aXiMinus = G4XiMinus::XiMinus();
225    //
226    // energetically possible to produce pion(s)  -->  inelastic scattering
227    //
228    G4double n, anpn;
229    GetNormalizationConstant( availableEnergy, n, anpn );
230    G4double ran = G4UniformRand();
231    G4double dum, excs = 0.0;
232    if( targetParticle.GetDefinition() == aProton )
233    {
234      counter = -1;
235      for( np=0; np<numSec/3 && ran>=excs; ++np )
236      {
237        for( nm=std::max(0,np-2); nm<=(np+1) && ran>=excs; ++nm )
238        {
239          for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
240          {
241            if( ++counter < numMul )
242            {
243              nt = np+nm+nz;
244              if( nt>0 && nt<=numSec )
245              {
246                test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
247                dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
248                if( std::fabs(dum) < 1.0 )
249                {
250                  if( test >= 1.0e-10 )excs += dum*test;
251                }
252                else
253                  excs += dum*test;
254              }
255            }
256          }
257        }
258      }
259      if( ran >= excs )  // 3 previous loops continued to the end
260      {
261        quasiElastic = true;
262        return;
263      }
264      np--; nm--; nz--;
265      //
266      // number of secondary mesons determined by kno distribution
267      // check for total charge of final state mesons to determine
268      // the kind of baryons to be produced, taking into account
269      // charge and strangeness conservation
270      //
271      if( np < nm+1 )
272      {
273        if( np != nm )   // charge mismatch
274        {
275          currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
276          incidentHasChanged = true;
277          //
278          // correct the strangeness by replacing a pi- by a kaon-
279          //
280          vec.Initialize( 1 );
281          G4ReactionProduct *p = new G4ReactionProduct;
282          p->SetDefinition( aKaonMinus );
283          (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
284          vec.SetElement( vecLen++, p );
285          --nm;
286        }
287      }
288      else if( np == nm+1 )
289      {
290        if( G4UniformRand() < 0.5 )
291        {
292          targetParticle.SetDefinitionAndUpdateE( aNeutron );
293          targetHasChanged = true;
294        }
295        else
296        {
297          currentParticle.SetDefinitionAndUpdateE( aXiMinus );
298          incidentHasChanged = true;
299        }
300      }
301      else
302      {
303        currentParticle.SetDefinitionAndUpdateE( aXiMinus );
304        incidentHasChanged = true;
305        targetParticle.SetDefinitionAndUpdateE( aNeutron );
306        targetHasChanged = true;
307      }
308    }
309    else  // target must be a neutron
310    {
311      counter = -1;
312      for( np=0; np<numSec/3 && ran>=excs; ++np )
313      {
314        for( nm=std::max(0,np-1); nm<=(np+2) && ran>=excs; ++nm )
315        {
316          for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
317          {
318            if( ++counter < numMul )
319            {
320              nt = np+nm+nz;
321              if( nt>0 && nt<=numSec )
322              {
323                test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
324                dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
325                if( std::fabs(dum) < 1.0 )
326                {
327                  if( test >= 1.0e-10 )excs += dum*test;
328                }
329                else
330                  excs += dum*test;
331              }
332            }
333          }
334        }
335      }
336      if( ran >= excs )  // 3 previous loops continued to the end
337      {
338        quasiElastic = true;
339        return;
340      }
341      np--; nm--; nz--;
342      if( np < nm )
343      {
344        if( np+1 == nm )
345        {
346          targetParticle.SetDefinitionAndUpdateE( aProton );
347          targetHasChanged = true;
348        }
349        else   // charge mismatch
350        {
351          currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
352          incidentHasChanged = true;
353          targetParticle.SetDefinitionAndUpdateE( aProton );
354          targetHasChanged = true;
355          //
356          // correct the strangeness by replacing a pi- by a kaon-
357          //
358          vec.Initialize( 1 );
359          G4ReactionProduct *p = new G4ReactionProduct;
360          p->SetDefinition( aKaonMinus );
361          (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
362          vec.SetElement( vecLen++, p );
363          --nm;
364        }
365      }
366      else if( np == nm )
367      {
368        if( G4UniformRand() >= 0.5 )
369        {
370          currentParticle.SetDefinitionAndUpdateE( aXiMinus );
371          incidentHasChanged = true;
372          targetParticle.SetDefinitionAndUpdateE( aProton );
373          targetHasChanged = true;
374        }
375      }
376      else
377      {
378        currentParticle.SetDefinitionAndUpdateE( aXiMinus );
379        incidentHasChanged = true;
380      }       
381    }
382  SetUpPions( np, nm, nz, vec, vecLen );
383  return;
384}
385
386 /* end of file */
387 
Note: See TracBrowser for help on using the repository browser.