source: trunk/source/processes/hadronic/models/rpg/src/G4RPGInelastic.cc@ 1038

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

update to geant4.9.2
File size: 20.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//
[962]26// $Id: G4RPGInelastic.cc,v 1.6 2008/03/22 00:03:24 dennis Exp $
[1007]27// GEANT4 tag $Name: geant4-09-02 $
[819]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
[962]170void
171G4RPGInelastic::CalculateMomenta(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)
[819]182{
183 cache = 0;
184 what = originalIncident->Get4Momentum().vect();
185
186 G4ReactionProduct leadingStrangeParticle;
187
[962]188 // strangeProduction.ReactionStage(originalIncident, modifiedOriginal,
189 // incidentHasChanged, originalTarget,
190 // targetParticle, targetHasChanged,
191 // targetNucleus, currentParticle,
192 // vec, vecLen,
193 // false, leadingStrangeParticle);
[819]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
[962]268 // Call fragmentation code if
269 // 1) there is annihilation, or
270 // 2) there are more than 5 secondaries, or
271 // 3) incident KE is > 1 GeV AND
272 // ( incident is a kaon AND rand < 0.5 OR twsup )
273 //
274
275 if( annihilation || vecLen > 5 ||
276 ( modifiedOriginal.GetKineticEnergy()/GeV >= 1.0 &&
277
[819]278 (((originalIncident->GetDefinition() == G4KaonPlus::KaonPlus() ||
279 originalIncident->GetDefinition() == G4KaonMinus::KaonMinus() ||
280 originalIncident->GetDefinition() == G4KaonZeroLong::KaonZeroLong() ||
281 originalIncident->GetDefinition() == G4KaonZeroShort::KaonZeroShort()) &&
[962]282 rand1 < 0.5)
283 || rand2 > twsup[vecLen]) ) )
[819]284
285 finishedGenXPt =
286 fragmentation.ReactionStage(originalIncident, modifiedOriginal,
287 incidentHasChanged, originalTarget,
288 targetParticle, targetHasChanged,
289 targetNucleus, currentParticle,
290 vec, vecLen,
291 leadFlag, leadingStrangeParticle);
292
[962]293 if (finishedGenXPt) return;
[819]294
295 G4bool finishedTwoClu = false;
[962]296
297 if (modifiedOriginal.GetTotalMomentum() < 1.0) {
298 for (G4int i = 0; i < vecLen; i++) delete vec[i];
[819]299 vecLen = 0;
[962]300
301 } else {
[819]302 // Occaisionally, GenerateXandPt will fail in the annihilation channel.
303 // Restore current, target and secondaries to pre-GenerateXandPt state
304 // before trying annihilation in TwoCluster
305
306 if (!finishedGenXPt && annihilation) {
307 currentParticle = saveCurrent;
308 targetParticle = saveTarget;
309 for (G4int i = 0; i < vecLen; i++) delete vec[i];
310 vecLen = 0;
311 vec.Initialize( 0 );
312 for (G4int i = 0; i < G4int(savevec.size()); i++) {
313 G4ReactionProduct* p = new G4ReactionProduct;
314 *p = savevec[i];
315 vec.SetElement( vecLen++, p );
316 }
317 }
318
[962]319 // Big violations of energy conservation in this method - don't use
320 //
321 // pionSuppression.ReactionStage(originalIncident, modifiedOriginal,
322 // incidentHasChanged, originalTarget,
323 // targetParticle, targetHasChanged,
324 // targetNucleus, currentParticle,
325 // vec, vecLen,
326 // false, leadingStrangeParticle);
[819]327
328 try
329 {
330 finishedTwoClu =
331 twoCluster.ReactionStage(originalIncident, modifiedOriginal,
332 incidentHasChanged, originalTarget,
333 targetParticle, targetHasChanged,
334 targetNucleus, currentParticle,
335 vec, vecLen,
336 leadFlag, leadingStrangeParticle);
337 }
338 catch(G4HadReentrentException aC)
339 {
340 aC.Report(G4cout);
341 throw G4HadReentrentException(__FILE__, __LINE__, "Failing to calculate momenta");
342 }
343 }
344
[962]345 if (finishedTwoClu) return;
[819]346
347 twoBody.ReactionStage(originalIncident, modifiedOriginal,
348 incidentHasChanged, originalTarget,
349 targetParticle, targetHasChanged,
350 targetNucleus, currentParticle,
351 vec, vecLen,
352 false, leadingStrangeParticle);
353}
354
[962]355/*
[819]356 void G4RPGInelastic::
357 Rotate(G4FastVector<G4ReactionProduct,256> &vec, G4int &vecLen)
358 {
359 G4double rotation = 2.*pi*G4UniformRand();
360 cache = rotation;
361 G4int i;
362 for( i=0; i<vecLen; ++i )
363 {
364 G4ThreeVector momentum = vec[i]->GetMomentum();
365 momentum = momentum.rotate(rotation, what);
366 vec[i]->SetMomentum(momentum);
367 }
368 }
[962]369*/
[819]370
371void
[962]372G4RPGInelastic::SetUpChange(G4FastVector<G4ReactionProduct,256>& vec,
373 G4int& vecLen,
374 G4ReactionProduct& currentParticle,
375 G4ReactionProduct& targetParticle,
376 G4bool& incidentHasChanged )
[819]377{
378 theParticleChange.Clear();
[962]379 G4ParticleDefinition* aKaonZL = G4KaonZeroLong::KaonZeroLong();
380 G4ParticleDefinition* aKaonZS = G4KaonZeroShort::KaonZeroShort();
[819]381 G4int i;
[962]382
383 if (currentParticle.GetDefinition() == particleDef[k0]) {
384 if (G4UniformRand() < 0.5) {
385 currentParticle.SetDefinitionAndUpdateE(aKaonZL);
[819]386 incidentHasChanged = true;
[962]387 } else {
388 currentParticle.SetDefinitionAndUpdateE(aKaonZS);
[819]389 }
[962]390 } else if (currentParticle.GetDefinition() == particleDef[k0b]) {
391 if (G4UniformRand() < 0.5) {
392 currentParticle.SetDefinitionAndUpdateE(aKaonZL);
393 } else {
394 currentParticle.SetDefinitionAndUpdateE(aKaonZS);
[819]395 incidentHasChanged = true;
396 }
397 }
398
[962]399 if (targetParticle.GetDefinition() == particleDef[k0] ||
400 targetParticle.GetDefinition() == particleDef[k0b] ) {
401 if (G4UniformRand() < 0.5) {
402 targetParticle.SetDefinitionAndUpdateE(aKaonZL);
403 } else {
404 targetParticle.SetDefinitionAndUpdateE(aKaonZS);
405 }
[819]406 }
[962]407
408 for (i = 0; i < vecLen; ++i) {
409 if (vec[i]->GetDefinition() == particleDef[k0] ||
410 vec[i]->GetDefinition() == particleDef[k0b] ) {
411 if (G4UniformRand() < 0.5) {
412 vec[i]->SetDefinitionAndUpdateE(aKaonZL);
413 } else {
414 vec[i]->SetDefinitionAndUpdateE(aKaonZS);
415 }
[819]416 }
417 }
418
[962]419 if (incidentHasChanged) {
[819]420 G4DynamicParticle* p0 = new G4DynamicParticle;
[962]421 p0->SetDefinition(currentParticle.GetDefinition() );
422 p0->SetMomentum(currentParticle.GetMomentum() );
[819]423 theParticleChange.AddSecondary( p0 );
424 theParticleChange.SetStatusChange( stopAndKill );
425 theParticleChange.SetEnergyChange( 0.0 );
[962]426
427 } else {
[819]428 G4double p = currentParticle.GetMomentum().mag()/MeV;
429 G4ThreeVector m = currentParticle.GetMomentum();
[962]430 if (p > DBL_MIN)
[819]431 theParticleChange.SetMomentumChange( m.x()/p, m.y()/p, m.z()/p );
432 else
433 theParticleChange.SetMomentumChange( 0.0, 0.0, 1.0 );
434
435 G4double aE = currentParticle.GetKineticEnergy();
436 if (std::fabs(aE)<.1*eV) aE=.1*eV;
437 theParticleChange.SetEnergyChange( aE );
438 }
439
[962]440 if (targetParticle.GetMass() > 0.0) // Tgt particle can be eliminated in TwoBody
[819]441 {
442 G4ThreeVector momentum = targetParticle.GetMomentum();
443 momentum = momentum.rotate(cache, what);
444 G4double targKE = targetParticle.GetKineticEnergy();
445 G4ThreeVector dir(0.0, 0.0, 1.0);
446 if (targKE < DBL_MIN)
447 targKE = DBL_MIN;
448 else
449 dir = momentum/momentum.mag();
450
451 G4DynamicParticle* p1 =
452 new G4DynamicParticle(targetParticle.GetDefinition(), dir, targKE);
453
454 theParticleChange.AddSecondary( p1 );
455 }
456
457 G4DynamicParticle* p;
[962]458 for (i = 0; i < vecLen; ++i) {
[819]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}
[962]472
473
474std::pair<G4int, G4double>
475G4RPGInelastic::interpolateEnergy(G4double e) const
476{
477 G4int index = 29;
478 G4double fraction = 0.0;
479
480 for (G4int i = 1; i < 30; i++) {
481 if (e < energyScale[i]) {
482 index = i-1;
483 fraction = (e - energyScale[index]) / (energyScale[i] - energyScale[index]);
484 break;
485 }
486 }
487 return std::pair<G4int, G4double>(index, fraction);
488}
489
490
491G4int
492G4RPGInelastic::sampleFlat(std::vector<G4double> sigma) const
493{
494 G4int i;
495 G4double sum(0.);
496 for (i = 0; i < G4int(sigma.size()); i++) sum += sigma[i];
497
498 G4double fsum = sum*G4UniformRand();
499 G4double partialSum = 0.0;
500 G4int channel = 0;
501
502 for (i = 0; i < G4int(sigma.size()); i++) {
503 partialSum += sigma[i];
504 if (fsum < partialSum) {
505 channel = i;
506 break;
507 }
508 }
509
510 return channel;
511}
512
513
514void G4RPGInelastic::CheckQnums(G4FastVector<G4ReactionProduct,256> &vec,
515 G4int &vecLen,
516 G4ReactionProduct &currentParticle,
517 G4ReactionProduct &targetParticle,
518 G4double Q, G4double B, G4double S)
519{
520 G4ParticleDefinition* projDef = currentParticle.GetDefinition();
521 G4ParticleDefinition* targDef = targetParticle.GetDefinition();
522 G4double chargeSum = projDef->GetPDGCharge() + targDef->GetPDGCharge();
523 G4double baryonSum = projDef->GetBaryonNumber() + targDef->GetBaryonNumber();
524 G4double strangenessSum = projDef->GetQuarkContent(3) -
525 projDef->GetAntiQuarkContent(3) +
526 targDef->GetQuarkContent(3) -
527 targDef->GetAntiQuarkContent(3);
528
529 G4ParticleDefinition* secDef = 0;
530 for (G4int i = 0; i < vecLen; i++) {
531 secDef = vec[i]->GetDefinition();
532 chargeSum += secDef->GetPDGCharge();
533 baryonSum += secDef->GetBaryonNumber();
534 strangenessSum += secDef->GetQuarkContent(3)
535 - secDef->GetAntiQuarkContent(3);
536 }
537
538 G4bool OK = true;
539 if (chargeSum != Q) {
540 G4cout << " Charge not conserved " << G4endl;
541 OK = false;
542 }
543 if (baryonSum != B) {
544 G4cout << " Baryon number not conserved " << G4endl;
545 OK = false;
546 }
547 if (strangenessSum != S) {
548 G4cout << " Strangeness not conserved " << G4endl;
549 OK = false;
550 }
551
552 if (!OK) {
553 G4cout << " projectile: " << projDef->GetParticleName()
554 << " target: " << targDef->GetParticleName() << G4endl;
555 for (G4int i = 0; i < vecLen; i++) {
556 secDef = vec[i]->GetDefinition();
557 G4cout << secDef->GetParticleName() << " " ;
558 }
559 G4cout << G4endl;
560 }
561
562}
563
564
565const G4double G4RPGInelastic::energyScale[30] = {
566 0.0, 0.01, 0.013, 0.018, 0.024, 0.032, 0.042, 0.056, 0.075, 0.1,
567 0.13, 0.18, 0.24, 0.32, 0.42, 0.56, 0.75, 1.0, 1.3, 1.8,
568 2.4, 3.2, 4.2, 5.6, 7.5, 10.0, 13.0, 18.0, 24.0, 32.0 };
569
570G4ParticleDefinition* p0 = G4PionZero::PionZero();
571G4ParticleDefinition* p1 = G4PionPlus::PionPlus();
572G4ParticleDefinition* p2 = G4PionMinus::PionMinus();
573G4ParticleDefinition* p3 = G4KaonPlus::KaonPlus();
574G4ParticleDefinition* p4 = G4KaonMinus::KaonMinus();
575G4ParticleDefinition* p5 = G4KaonZero::KaonZero();
576G4ParticleDefinition* p6 = G4AntiKaonZero::AntiKaonZero();
577G4ParticleDefinition* p7 = G4Proton::Proton();
578G4ParticleDefinition* p8 = G4Neutron::Neutron();
579G4ParticleDefinition* p9 = G4Lambda::Lambda();
580G4ParticleDefinition* p10 = G4SigmaPlus::SigmaPlus();
581G4ParticleDefinition* p11 = G4SigmaZero::SigmaZero();
582G4ParticleDefinition* p12 = G4SigmaMinus::SigmaMinus();
583G4ParticleDefinition* p13 = G4XiZero::XiZero();
584G4ParticleDefinition* p14 = G4XiMinus::XiMinus();
585G4ParticleDefinition* p15 = G4OmegaMinus::OmegaMinus();
586G4ParticleDefinition* p16 = G4AntiProton::AntiProton();
587G4ParticleDefinition* p17 = G4AntiNeutron::AntiNeutron();
588
589G4ParticleDefinition* G4RPGInelastic::particleDef[18] = {
590 p0, p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12, p13, p14,
591 p15, p16, p17 };
592
[819]593/* end of file */
Note: See TracBrowser for help on using the repository browser.