source: trunk/source/processes/hadronic/models/management/src/G4InelasticInteraction.cc @ 1340

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

update ti head
File size: 16.9 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: G4InelasticInteraction.cc,v 1.12 2009/01/24 11:56:27 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-03-ref-09 $
28//
29// Hadronic Process: Inelastic Interaction
30// original by H.P. Wellisch
31// modified by J.L. Chuma, TRIUMF, 22-Nov-1996
32// Last modified: 27-Mar-1997
33// J.P. Wellisch: 23-Apr-97: throw G4HadronicException(__FILE__, __LINE__,  removed
34// J.P. Wellisch: 24-Apr-97: correction for SetUpPions
35// Modified by J.L. Chuma, 30-Apr-97: added originalTarget to CalculateMomenta
36//                                    since TwoBody needed to reset the target particle
37// J.L. Chuma, 20-Jun-97: Modified CalculateMomenta to correct the decision process
38//                        for whether to use GenerateXandPt or TwoCluster
39// J.L. Chuma, 06-Aug-97: added original incident particle, before Fermi motion and
40//                        evaporation effects are included, needed for calculating
41//                        self absorption and corrections for single particle spectra
42// HPW removed misunderstanding of LocalEnergyDeposit, 11.04.98.
43// 23-Jan-2009 V.Ivanchenko move constructor and destructor to the body
44 
45#include "G4InelasticInteraction.hh"
46#include "Randomize.hh"
47#include "G4HadReentrentException.hh"
48
49G4InelasticInteraction::G4InelasticInteraction(const G4String& modelName) 
50  : G4HadronicInteraction(modelName)
51{ cache = 0.0;}
52   
53G4InelasticInteraction::~G4InelasticInteraction()
54{}
55 
56G4double
57  G4InelasticInteraction::Pmltpc(      // used in Cascade functions
58   G4int np, G4int nm, G4int nz, G4int n, G4double b, G4double c )
59  {
60    const G4double expxu =  82.;           // upper bound for arg. of exp
61    const G4double expxl = -expxu;         // lower bound for arg. of exp
62    G4double npf = 0.0;
63    G4double nmf = 0.0;
64    G4double nzf = 0.0;
65    G4int i;
66    for( i=2; i<=np; i++ )npf += std::log((double)i);
67    for( i=2; i<=nm; i++ )nmf += std::log((double)i);
68    for( i=2; i<=nz; i++ )nzf += std::log((double)i);
69    G4double r;
70    r = std::min( expxu, std::max( expxl, -(np-nm+nz+b)*(np-nm+nz+b)/(2*c*c*n*n)-npf-nmf-nzf ) );
71    return std::exp(r);
72  }
73
74G4bool
75  G4InelasticInteraction::MarkLeadingStrangeParticle(
76   const G4ReactionProduct &currentParticle,
77   const G4ReactionProduct &targetParticle,
78   G4ReactionProduct &leadParticle )
79  {
80    // the following was in GenerateXandPt and TwoCluster
81    // add a parameter to the GenerateXandPt function telling it about the strange particle
82    //
83    // assumes that the original particle was a strange particle
84    //
85    G4bool lead = false;
86    if( (currentParticle.GetMass() >= G4KaonPlus::KaonPlus()->GetPDGMass()) &&
87        (currentParticle.GetDefinition() != G4Proton::Proton()) &&
88        (currentParticle.GetDefinition() != G4Neutron::Neutron()) )
89    {
90      lead = true;
91      leadParticle = currentParticle;              //  set lead to the incident particle
92    }
93    else if( (targetParticle.GetMass() >= G4KaonPlus::KaonPlus()->GetPDGMass()) &&
94             (targetParticle.GetDefinition() != G4Proton::Proton()) &&
95             (targetParticle.GetDefinition() != G4Neutron::Neutron()) )
96    {
97      lead = true;
98      leadParticle = targetParticle;              //   set lead to the target particle
99    }
100    return lead;
101  }
102 
103void
104  G4InelasticInteraction::SetUpPions(
105   const G4int np,
106   const G4int nm,
107   const G4int nz,
108   G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
109   G4int &vecLen )
110  {
111    if( np+nm+nz == 0 )return;
112    G4int i;
113    G4ReactionProduct *p;
114    for( i=0; i<np; ++i )
115    {
116      p = new G4ReactionProduct;
117      p->SetDefinition( G4PionPlus::PionPlus() );
118      (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
119      vec.SetElement( vecLen++, p );
120    }
121    for( i=np; i<np+nm; ++i )
122    {
123      p = new G4ReactionProduct;
124      p->SetDefinition( G4PionMinus::PionMinus() );
125      (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
126      vec.SetElement( vecLen++, p );
127    }
128    for( i=np+nm; i<np+nm+nz; ++i )
129    {
130      p = new G4ReactionProduct;
131      p->SetDefinition( G4PionZero::PionZero() );
132      (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
133      vec.SetElement( vecLen++, p );
134    }
135  }
136 
137void
138  G4InelasticInteraction::GetNormalizationConstant(
139   const G4double energy,  // MeV,  <0 means annihilation channels
140   G4double &n,
141   G4double &anpn )
142  {
143    const G4double expxu =  82.;          // upper bound for arg. of exp
144    const G4double expxl = -expxu;        // lower bound for arg. of exp
145    const G4int numSec = 60;
146    //
147    // the only difference between the calculation for annihilation channels
148    // and normal is the starting value, iBegin, for the loop below
149    //
150    G4int iBegin = 1;
151    G4double en = energy;
152    if( energy < 0.0 )
153    {
154      iBegin = 2;
155      en *= -1.0;
156    }
157    //
158    // number of total particles vs. centre of mass Energy - 2*proton mass
159    //
160    G4double aleab = std::log(en/GeV);
161    n = 3.62567 + aleab*(0.665843 + aleab*(0.336514 + aleab*(0.117712 + 0.0136912*aleab)));
162    n -= 2.0;
163    //
164    // normalization constant for kno-distribution
165    //
166    anpn = 0.0;
167    G4double test, temp;
168    for( G4int i=iBegin; i<=numSec; ++i )
169    {
170      temp = pi*i/(2.0*n*n);
171      test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(i*i)/(n*n) ) ) );
172      if( temp < 1.0 )
173      {
174        if( test >= 1.0e-10 )anpn += temp*test;
175      }
176      else
177        anpn += temp*test;
178    }
179  }
180 
181void
182  G4InelasticInteraction::CalculateMomenta(
183   G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
184   G4int &vecLen,
185   const G4HadProjectile *originalIncident,   // the original incident particle
186   const G4DynamicParticle *originalTarget,
187   G4ReactionProduct &modifiedOriginal,   // Fermi motion and evap. effects included
188   G4Nucleus &targetNucleus,
189   G4ReactionProduct &currentParticle,
190   G4ReactionProduct &targetParticle,
191   G4bool &incidentHasChanged,
192   G4bool &targetHasChanged,
193   G4bool quasiElastic )
194  {
195    cache = 0;
196    what = originalIncident->Get4Momentum().vect();
197
198
199    theReactionDynamics.ProduceStrangeParticlePairs( vec, vecLen,
200                                                     modifiedOriginal, originalTarget,
201                                                     currentParticle, targetParticle,
202                                                     incidentHasChanged, targetHasChanged );
203
204    if( quasiElastic )
205    {
206      theReactionDynamics.TwoBody( vec, vecLen,
207                                   modifiedOriginal, originalTarget,
208                                   currentParticle, targetParticle,
209                                   targetNucleus, targetHasChanged );
210      return;
211    }
212    G4ReactionProduct leadingStrangeParticle;
213    G4bool leadFlag = MarkLeadingStrangeParticle( currentParticle,
214                                                  targetParticle,
215                                                  leadingStrangeParticle );
216    //
217    // Note: the number of secondaries can be reduced in GenerateXandPt and TwoCluster
218    //
219    G4bool finishedGenXPt = false;
220    G4bool annihilation = false;
221    if( originalIncident->GetDefinition()->GetPDGEncoding() < 0 &&
222        currentParticle.GetMass() == 0.0 && targetParticle.GetMass() == 0.0 )
223    {
224      // original was an anti-particle and annihilation has taken place
225      annihilation = true;
226      G4double ekcor = 1.0;
227      G4double ek = originalIncident->GetKineticEnergy();
228      G4double ekOrg = ek;
229     
230      const G4double tarmas = originalTarget->GetDefinition()->GetPDGMass();
231      if( ek > 1.0*GeV )ekcor = 1./(ek/GeV);
232      const G4double atomicWeight = targetNucleus.GetN();
233      ek = 2*tarmas + ek*(1.+ekcor/atomicWeight);
234      G4double tkin = targetNucleus.Cinema(ek);
235      ek += tkin;
236      ekOrg += tkin;
237      //      modifiedOriginal.SetKineticEnergy( ekOrg );
238      //
239      // evaporation --  re-calculate black track energies
240      //                 this was Done already just before the cascade
241      //
242      tkin = targetNucleus.AnnihilationEvaporationEffects(ek, ekOrg);
243      ekOrg -= tkin;
244      ekOrg = std::max( 0.0001*GeV, ekOrg );
245      modifiedOriginal.SetKineticEnergy( ekOrg );
246      G4double amas = originalIncident->GetDefinition()->GetPDGMass();
247      G4double et = ekOrg + amas;
248      G4double p = std::sqrt( std::abs(et*et-amas*amas) );
249      G4double pp = modifiedOriginal.GetMomentum().mag();
250      if( pp > 0.0 )
251      {
252        G4ThreeVector momentum = modifiedOriginal.GetMomentum();
253        modifiedOriginal.SetMomentum( momentum * (p/pp) );
254      }
255      if( ekOrg <= 0.0001 )
256      {
257        modifiedOriginal.SetKineticEnergy( 0.0 );
258        modifiedOriginal.SetMomentum( 0.0, 0.0, 0.0 );
259      }
260    }
261    const G4double twsup[] = { 1.0, 0.7, 0.5, 0.3, 0.2, 0.1 };
262    G4double rand1 = G4UniformRand();
263    G4double rand2 = G4UniformRand();
264
265    // Cache current, target, and secondaries
266    G4ReactionProduct saveCurrent = currentParticle;
267    G4ReactionProduct saveTarget = targetParticle;
268    std::vector<G4ReactionProduct> savevec;
269    for (G4int i = 0; i < vecLen; i++) savevec.push_back(*vec[i]);
270
271    if (annihilation || 
272        vecLen >= 6 ||
273        ( modifiedOriginal.GetKineticEnergy()/GeV >= 1.0 &&
274          ( ( (originalIncident->GetDefinition() == G4KaonPlus::KaonPlus() ||
275               originalIncident->GetDefinition() == G4KaonMinus::KaonMinus() ||
276               originalIncident->GetDefinition() == G4KaonZeroLong::KaonZeroLong() ||
277               originalIncident->GetDefinition() == G4KaonZeroShort::KaonZeroShort() ) 
278               &&
279               rand1 < 0.5 ) 
280            || rand2 > twsup[vecLen] )  )  )
281
282      finishedGenXPt =
283        theReactionDynamics.GenerateXandPt( vec, vecLen,
284                                            modifiedOriginal, originalIncident,
285                                            currentParticle, targetParticle,
286                                            originalTarget,
287                                            targetNucleus, incidentHasChanged,
288                                            targetHasChanged, leadFlag,
289                                            leadingStrangeParticle );
290    if( finishedGenXPt )
291    {
292      Rotate(vec, vecLen);
293      return;
294    }
295
296    G4bool finishedTwoClu = false;
297    if( modifiedOriginal.GetTotalMomentum()/MeV < 1.0 )
298    {
299      for(G4int i=0; i<vecLen; i++) delete vec[i];
300      vecLen = 0;
301    }
302    else
303    {
304      // Occaisionally, GenerateXandPt will fail in the annihilation channel.
305      // Restore current, target and secondaries to pre-GenerateXandPt state
306      // before trying annihilation in TwoCluster
307
308      if (!finishedGenXPt && annihilation) {
309        currentParticle = saveCurrent;
310        targetParticle = saveTarget;
311        for (G4int i = 0; i < vecLen; i++) delete vec[i];
312        vecLen = 0;
313        vec.Initialize( 0 );
314        for (G4int i = 0; i < G4int(savevec.size()); i++) {
315          G4ReactionProduct* p = new G4ReactionProduct;
316          *p = savevec[i];
317          vec.SetElement( vecLen++, p );
318        }
319      }
320
321      theReactionDynamics.SuppressChargedPions( vec, vecLen,
322                                      modifiedOriginal, currentParticle,
323                                      targetParticle, targetNucleus,
324                                      incidentHasChanged, targetHasChanged );
325      try
326      {
327      finishedTwoClu = theReactionDynamics.TwoCluster( vec, vecLen,
328                                      modifiedOriginal, originalIncident,
329                                      currentParticle, targetParticle,
330                                      originalTarget,
331                                      targetNucleus, incidentHasChanged,
332                                      targetHasChanged, leadFlag,
333                                      leadingStrangeParticle );
334       }
335       catch(G4HadReentrentException aC)
336       {
337         aC.Report(G4cout);
338         throw G4HadReentrentException(__FILE__, __LINE__, "Failing to calculate momenta");
339       }
340    }
341
342    if( finishedTwoClu )
343    {
344      Rotate(vec, vecLen);
345      return;
346    }
347
348    theReactionDynamics.TwoBody( vec, vecLen,
349                                 modifiedOriginal, originalTarget,
350                                 currentParticle, targetParticle,
351                                 targetNucleus, targetHasChanged );
352  }
353
354 
355void G4InelasticInteraction::
356 Rotate(G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec, G4int &vecLen)
357 {
358   G4double rotation = 2.*pi*G4UniformRand();
359   cache = rotation;
360   G4int i;
361   for( i=0; i<vecLen; ++i )
362   {
363     G4ThreeVector momentum = vec[i]->GetMomentum();
364     momentum = momentum.rotate(rotation, what);
365     vec[i]->SetMomentum(momentum);
366   }
367 }     
368
369void
370  G4InelasticInteraction::SetUpChange(
371   G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
372   G4int &vecLen,
373   G4ReactionProduct &currentParticle,
374   G4ReactionProduct &targetParticle,
375   G4bool &incidentHasChanged )
376  {
377    theParticleChange.Clear();
378    G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong();
379    G4ParticleDefinition *aKaonZS = G4KaonZeroShort::KaonZeroShort();
380    G4int i;
381    if( currentParticle.GetDefinition() == aKaonZL )
382    {
383      if( G4UniformRand() <= 0.5 )
384      {
385        currentParticle.SetDefinition( aKaonZS );
386        incidentHasChanged = true;
387      }
388    }
389    else if( currentParticle.GetDefinition() == aKaonZS )
390    {
391      if( G4UniformRand() > 0.5 )
392      {
393        currentParticle.SetDefinition( aKaonZL );
394        incidentHasChanged = true;
395      }
396    }
397    if( targetParticle.GetDefinition() == aKaonZL )
398    {
399      if( G4UniformRand() <= 0.5 )targetParticle.SetDefinition( aKaonZS );
400    }
401    else if( targetParticle.GetDefinition() == aKaonZS )
402    {
403      if( G4UniformRand() > 0.5 )targetParticle.SetDefinition( aKaonZL );
404    }
405    for( i=0; i<vecLen; ++i )
406    {
407      if( vec[i]->GetDefinition() == aKaonZL )
408      {
409        if( G4UniformRand() <= 0.5 )vec[i]->SetDefinition( aKaonZS );
410      }
411      else if( vec[i]->GetDefinition() == aKaonZS )
412      {
413        if( G4UniformRand() > 0.5 )vec[i]->SetDefinition( aKaonZL );
414      }
415    }     
416    if( incidentHasChanged )
417    {
418      G4DynamicParticle* p0 = new G4DynamicParticle;
419      p0->SetDefinition( currentParticle.GetDefinition() );
420      p0->SetMomentum( currentParticle.GetMomentum() );
421      theParticleChange.AddSecondary( p0 );
422      theParticleChange.SetStatusChange( stopAndKill );
423      theParticleChange.SetEnergyChange( 0.0 );
424    }
425    else
426    {
427      G4double p = currentParticle.GetMomentum().mag()/MeV;
428      G4ThreeVector m = currentParticle.GetMomentum();
429      if( p > DBL_MIN )
430        theParticleChange.SetMomentumChange( m.x()/p, m.y()/p, m.z()/p );
431      else
432        theParticleChange.SetMomentumChange( 1.0, 0.0, 0.0 );
433
434      G4double aE = currentParticle.GetKineticEnergy();
435      if (std::fabs(aE)<.1*eV) aE=.1*eV;
436      theParticleChange.SetEnergyChange( aE );
437    }
438
439    if( targetParticle.GetMass() > 0.0 )  // targetParticle can be eliminated in TwoBody
440    {
441      G4DynamicParticle *p1 = new G4DynamicParticle;
442      p1->SetDefinition( targetParticle.GetDefinition() );
443      G4ThreeVector momentum = targetParticle.GetMomentum();
444      momentum = momentum.rotate(cache, what);
445      p1->SetMomentum( momentum );
446      theParticleChange.AddSecondary( p1 );
447    }
448
449    G4DynamicParticle *p;
450    for( i=0; i<vecLen; ++i )
451    {
452      p = new G4DynamicParticle();
453      p->SetDefinition( vec[i]->GetDefinition() );
454      p->SetMomentum( vec[i]->GetMomentum() );
455      theParticleChange.AddSecondary( p );
456      delete vec[i];
457    }
458  }
459 
460  /* end of file */
461 
Note: See TracBrowser for help on using the repository browser.