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

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

update CVS release candidate geant4.9.3.01
File size: 16.9 KB
RevLine 
[819]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//
[1055]26// $Id: G4InelasticInteraction.cc,v 1.12 2009/01/24 11:56:27 vnivanch Exp $
[1196]27// GEANT4 tag $Name: geant4-09-03-cand-01 $
[819]28//
[1055]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
[819]44
45#include "G4InelasticInteraction.hh"
46#include "Randomize.hh"
47#include "G4HadReentrentException.hh"
[1055]48
49G4InelasticInteraction::G4InelasticInteraction(const G4String& modelName)
50 : G4HadronicInteraction(modelName)
51{ cache = 0.0;}
52
53G4InelasticInteraction::~G4InelasticInteraction()
54{}
[819]55
[1055]56G4double
[819]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
[1055]74G4bool
[819]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
[1055]103void
[819]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
[1055]137void
[819]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
[1055]181void
[819]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
[1055]355void G4InelasticInteraction::
[819]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
[1055]369void
[819]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.