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

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

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