source: trunk/source/processes/hadronic/models/rpg/src/G4RPGInelastic.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.0 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: G4RPGInelastic.cc,v 1.2 2007/08/15 20:38:25 dennis Exp $
27// GEANT4 tag $Name: geant4-09-01-patch-02 $
28//
29
30#include "G4RPGInelastic.hh"
31#include "Randomize.hh"
32#include "G4HadReentrentException.hh"
33#include "G4RPGStrangeProduction.hh"
34#include "G4RPGTwoBody.hh"
35
36
37G4double G4RPGInelastic::Pmltpc(G4int np, G4int nm, G4int nz,
38 G4int n, G4double b, G4double c)
39{
40 const G4double expxu = 82.; // upper bound for arg. of exp
41 const G4double expxl = -expxu; // lower bound for arg. of exp
42 G4double npf = 0.0;
43 G4double nmf = 0.0;
44 G4double nzf = 0.0;
45 G4int i;
46 for( i=2; i<=np; i++ )npf += std::log((double)i);
47 for( i=2; i<=nm; i++ )nmf += std::log((double)i);
48 for( i=2; i<=nz; i++ )nzf += std::log((double)i);
49 G4double r;
50 r = std::min( expxu, std::max( expxl, -(np-nm+nz+b)*(np-nm+nz+b)/(2*c*c*n*n)-npf-nmf-nzf ) );
51 return std::exp(r);
52}
53
54
55G4int G4RPGInelastic::Factorial( G4int n )
56{
57 G4int m = std::min(n,10);
58 G4int result = 1;
59 if( m <= 1 )return result;
60 for( G4int i=2; i<=m; ++i )result *= i;
61 return result;
62}
63
64
65G4bool G4RPGInelastic::MarkLeadingStrangeParticle(
66 const G4ReactionProduct &currentParticle,
67 const G4ReactionProduct &targetParticle,
68 G4ReactionProduct &leadParticle )
69{
70 // The following was in GenerateXandPt and TwoCluster.
71 // Add a parameter to the GenerateXandPt function telling it about the
72 // 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
95 void G4RPGInelastic::SetUpPions(const G4int np, const G4int nm,
96 const G4int nz,
97 G4FastVector<G4ReactionProduct,256> &vec,
98 G4int &vecLen)
99 {
100 if( np+nm+nz == 0 )return;
101 G4int i;
102 G4ReactionProduct *p;
103 for( i=0; i<np; ++i )
104 {
105 p = new G4ReactionProduct;
106 p->SetDefinition( G4PionPlus::PionPlus() );
107 (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
108 vec.SetElement( vecLen++, p );
109 }
110 for( i=np; i<np+nm; ++i )
111 {
112 p = new G4ReactionProduct;
113 p->SetDefinition( G4PionMinus::PionMinus() );
114 (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
115 vec.SetElement( vecLen++, p );
116 }
117 for( i=np+nm; i<np+nm+nz; ++i )
118 {
119 p = new G4ReactionProduct;
120 p->SetDefinition( G4PionZero::PionZero() );
121 (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
122 vec.SetElement( vecLen++, p );
123 }
124 }
125
126
127 void G4RPGInelastic::GetNormalizationConstant(
128 const G4double energy, // MeV, <0 means annihilation channels
129 G4double &n,
130 G4double &anpn )
131 {
132 const G4double expxu = 82.; // upper bound for arg. of exp
133 const G4double expxl = -expxu; // lower bound for arg. of exp
134 const G4int numSec = 60;
135 //
136 // the only difference between the calculation for annihilation channels
137 // and normal is the starting value, iBegin, for the loop below
138 //
139 G4int iBegin = 1;
140 G4double en = energy;
141 if( energy < 0.0 )
142 {
143 iBegin = 2;
144 en *= -1.0;
145 }
146 //
147 // number of total particles vs. centre of mass Energy - 2*proton mass
148 //
149 G4double aleab = std::log(en/GeV);
150 n = 3.62567 + aleab*(0.665843 + aleab*(0.336514 + aleab*(0.117712 + 0.0136912*aleab)));
151 n -= 2.0;
152 //
153 // normalization constant for kno-distribution
154 //
155 anpn = 0.0;
156 G4double test, temp;
157 for( G4int i=iBegin; i<=numSec; ++i )
158 {
159 temp = pi*i/(2.0*n*n);
160 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(i*i)/(n*n) ) ) );
161 if( temp < 1.0 )
162 {
163 if( test >= 1.0e-10 )anpn += temp*test;
164 }
165 else
166 anpn += temp*test;
167 }
168 }
169
170void G4RPGInelastic::CalculateMomenta(
171 G4FastVector<G4ReactionProduct,256> &vec,
172 G4int &vecLen,
173 const G4HadProjectile *originalIncident,
174 const G4DynamicParticle *originalTarget,
175 G4ReactionProduct &modifiedOriginal,
176 G4Nucleus &targetNucleus,
177 G4ReactionProduct &currentParticle,
178 G4ReactionProduct &targetParticle,
179 G4bool &incidentHasChanged,
180 G4bool &targetHasChanged,
181 G4bool quasiElastic )
182{
183 cache = 0;
184 what = originalIncident->Get4Momentum().vect();
185
186 G4ReactionProduct leadingStrangeParticle;
187
188 strangeProduction.ReactionStage(originalIncident, modifiedOriginal,
189 incidentHasChanged, originalTarget,
190 targetParticle, targetHasChanged,
191 targetNucleus, currentParticle,
192 vec, vecLen,
193 false, leadingStrangeParticle);
194
195 if( quasiElastic )
196 {
197 twoBody.ReactionStage(originalIncident, modifiedOriginal,
198 incidentHasChanged, originalTarget,
199 targetParticle, targetHasChanged,
200 targetNucleus, currentParticle,
201 vec, vecLen,
202 false, leadingStrangeParticle);
203 return;
204 }
205
206 G4bool leadFlag = MarkLeadingStrangeParticle(currentParticle,
207 targetParticle,
208 leadingStrangeParticle );
209 //
210 // Note: the number of secondaries can be reduced in GenerateXandPt
211 // and TwoCluster
212 //
213 G4bool finishedGenXPt = false;
214 G4bool annihilation = false;
215 if( originalIncident->GetDefinition()->GetPDGEncoding() < 0 &&
216 currentParticle.GetMass() == 0.0 && targetParticle.GetMass() == 0.0 )
217 {
218 // original was an anti-particle and annihilation has taken place
219 annihilation = true;
220 G4double ekcor = 1.0;
221 G4double ek = originalIncident->GetKineticEnergy();
222 G4double ekOrg = ek;
223
224 const G4double tarmas = originalTarget->GetDefinition()->GetPDGMass();
225 if( ek > 1.0*GeV )ekcor = 1./(ek/GeV);
226 const G4double atomicWeight = targetNucleus.GetN();
227 ek = 2*tarmas + ek*(1.+ekcor/atomicWeight);
228 G4double tkin = targetNucleus.Cinema(ek);
229 ek += tkin;
230 ekOrg += tkin;
231 // modifiedOriginal.SetKineticEnergy( ekOrg );
232 //
233 // evaporation -- re-calculate black track energies
234 // this was Done already just before the cascade
235 //
236 tkin = targetNucleus.AnnihilationEvaporationEffects(ek, ekOrg);
237 ekOrg -= tkin;
238 ekOrg = std::max( 0.0001*GeV, ekOrg );
239 modifiedOriginal.SetKineticEnergy( ekOrg );
240 G4double amas = originalIncident->GetDefinition()->GetPDGMass();
241 G4double et = ekOrg + amas;
242 G4double p = std::sqrt( std::abs(et*et-amas*amas) );
243 G4double pp = modifiedOriginal.GetMomentum().mag();
244 if( pp > 0.0 )
245 {
246 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
247 modifiedOriginal.SetMomentum( momentum * (p/pp) );
248 }
249 if( ekOrg <= 0.0001 )
250 {
251 modifiedOriginal.SetKineticEnergy( 0.0 );
252 modifiedOriginal.SetMomentum( 0.0, 0.0, 0.0 );
253 }
254 }
255
256 // twsup gives percentage of time two-cluster model is called
257
258 const G4double twsup[] = { 1.0, 0.7, 0.5, 0.3, 0.2, 0.1 };
259 G4double rand1 = G4UniformRand();
260 G4double rand2 = G4UniformRand();
261
262 // Cache current, target, and secondaries
263 G4ReactionProduct saveCurrent = currentParticle;
264 G4ReactionProduct saveTarget = targetParticle;
265 std::vector<G4ReactionProduct> savevec;
266 for (G4int i = 0; i < vecLen; i++) savevec.push_back(*vec[i]);
267
268 if( annihilation || (vecLen >= 6) ||
269 (modifiedOriginal.GetKineticEnergy()/GeV >= 1.0) &&
270 (((originalIncident->GetDefinition() == G4KaonPlus::KaonPlus() ||
271 originalIncident->GetDefinition() == G4KaonMinus::KaonMinus() ||
272 originalIncident->GetDefinition() == G4KaonZeroLong::KaonZeroLong() ||
273 originalIncident->GetDefinition() == G4KaonZeroShort::KaonZeroShort()) &&
274 rand1 < 0.5) || rand2 > twsup[vecLen]) )
275
276 finishedGenXPt =
277 fragmentation.ReactionStage(originalIncident, modifiedOriginal,
278 incidentHasChanged, originalTarget,
279 targetParticle, targetHasChanged,
280 targetNucleus, currentParticle,
281 vec, vecLen,
282 leadFlag, leadingStrangeParticle);
283
284 if( finishedGenXPt )
285 {
286 Rotate(vec, vecLen);
287 return;
288 }
289
290 G4bool finishedTwoClu = false;
291 if( modifiedOriginal.GetTotalMomentum()/MeV < 1.0 )
292 {
293 for(G4int i=0; i<vecLen; i++) delete vec[i];
294 vecLen = 0;
295 }
296 else
297 {
298 // Occaisionally, GenerateXandPt will fail in the annihilation channel.
299 // Restore current, target and secondaries to pre-GenerateXandPt state
300 // before trying annihilation in TwoCluster
301
302 if (!finishedGenXPt && annihilation) {
303 currentParticle = saveCurrent;
304 targetParticle = saveTarget;
305 for (G4int i = 0; i < vecLen; i++) delete vec[i];
306 vecLen = 0;
307 vec.Initialize( 0 );
308 for (G4int i = 0; i < G4int(savevec.size()); i++) {
309 G4ReactionProduct* p = new G4ReactionProduct;
310 *p = savevec[i];
311 vec.SetElement( vecLen++, p );
312 }
313 }
314
315 pionSuppression.ReactionStage(originalIncident, modifiedOriginal,
316 incidentHasChanged, originalTarget,
317 targetParticle, targetHasChanged,
318 targetNucleus, currentParticle,
319 vec, vecLen,
320 false, leadingStrangeParticle);
321
322 try
323 {
324 finishedTwoClu =
325 twoCluster.ReactionStage(originalIncident, modifiedOriginal,
326 incidentHasChanged, originalTarget,
327 targetParticle, targetHasChanged,
328 targetNucleus, currentParticle,
329 vec, vecLen,
330 leadFlag, leadingStrangeParticle);
331 }
332 catch(G4HadReentrentException aC)
333 {
334 aC.Report(G4cout);
335 throw G4HadReentrentException(__FILE__, __LINE__, "Failing to calculate momenta");
336 }
337 }
338
339 if( finishedTwoClu )
340 {
341 Rotate(vec, vecLen);
342 return;
343 }
344
345 twoBody.ReactionStage(originalIncident, modifiedOriginal,
346 incidentHasChanged, originalTarget,
347 targetParticle, targetHasChanged,
348 targetNucleus, currentParticle,
349 vec, vecLen,
350 false, leadingStrangeParticle);
351}
352
353
354 void G4RPGInelastic::
355 Rotate(G4FastVector<G4ReactionProduct,256> &vec, G4int &vecLen)
356 {
357 G4double rotation = 2.*pi*G4UniformRand();
358 cache = rotation;
359 G4int i;
360 for( i=0; i<vecLen; ++i )
361 {
362 G4ThreeVector momentum = vec[i]->GetMomentum();
363 momentum = momentum.rotate(rotation, what);
364 vec[i]->SetMomentum(momentum);
365 }
366 }
367
368void
369G4RPGInelastic::SetUpChange(G4FastVector<G4ReactionProduct,256> &vec,
370 G4int &vecLen,
371 G4ReactionProduct &currentParticle,
372 G4ReactionProduct &targetParticle,
373 G4bool &incidentHasChanged )
374{
375 theParticleChange.Clear();
376 G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong();
377 G4ParticleDefinition *aKaonZS = G4KaonZeroShort::KaonZeroShort();
378 G4int i;
379 if( currentParticle.GetDefinition() == aKaonZL )
380 {
381 if( G4UniformRand() <= 0.5 )
382 {
383 currentParticle.SetDefinition( aKaonZS );
384 incidentHasChanged = true;
385 }
386 }
387 else if( currentParticle.GetDefinition() == aKaonZS )
388 {
389 if( G4UniformRand() > 0.5 )
390 {
391 currentParticle.SetDefinition( aKaonZL );
392 incidentHasChanged = true;
393 }
394 }
395
396 if( targetParticle.GetDefinition() == aKaonZL )
397 {
398 if( G4UniformRand() <= 0.5 )targetParticle.SetDefinition( aKaonZS );
399 }
400 else if( targetParticle.GetDefinition() == aKaonZS )
401 {
402 if( G4UniformRand() > 0.5 )targetParticle.SetDefinition( aKaonZL );
403 }
404 for( i=0; i<vecLen; ++i )
405 {
406 if( vec[i]->GetDefinition() == aKaonZL )
407 {
408 if( G4UniformRand() <= 0.5 )vec[i]->SetDefinition( aKaonZS );
409 }
410 else if( vec[i]->GetDefinition() == aKaonZS )
411 {
412 if( G4UniformRand() > 0.5 )vec[i]->SetDefinition( aKaonZL );
413 }
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( 0.0, 0.0, 1.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 ) // Tgt particle can be eliminated in TwoBody
440 {
441 G4ThreeVector momentum = targetParticle.GetMomentum();
442 momentum = momentum.rotate(cache, what);
443 G4double targKE = targetParticle.GetKineticEnergy();
444 G4ThreeVector dir(0.0, 0.0, 1.0);
445 if (targKE < DBL_MIN)
446 targKE = DBL_MIN;
447 else
448 dir = momentum/momentum.mag();
449
450 G4DynamicParticle* p1 =
451 new G4DynamicParticle(targetParticle.GetDefinition(), dir, targKE);
452
453 theParticleChange.AddSecondary( p1 );
454 }
455
456 G4DynamicParticle* p;
457 for( i=0; i<vecLen; ++i )
458 {
459 G4double secKE = vec[i]->GetKineticEnergy();
460 G4ThreeVector momentum = vec[i]->GetMomentum();
461 G4ThreeVector dir(0.0, 0.0, 1.0);
462 if (secKE < DBL_MIN)
463 secKE = DBL_MIN;
464 else
465 dir = momentum/momentum.mag();
466
467 p = new G4DynamicParticle(vec[i]->GetDefinition(), dir, secKE);
468 theParticleChange.AddSecondary( p );
469 delete vec[i];
470 }
471}
472
473/* end of file */
Note: See TracBrowser for help on using the repository browser.