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

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

update CVS release candidate geant4.9.3.01
File size: 20.4 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.7 2009/05/25 19:07:15 dennis Exp $
27// GEANT4 tag $Name: geant4-09-03-cand-01 $
28//
29
30#include "G4RPGInelastic.hh"
31#include "Randomize.hh"
32#include "G4HadReentrentException.hh"
33#include "G4RPGStrangeProduction.hh"
34#include "G4RPGTwoBody.hh"
35
36
37G4RPGInelastic::G4RPGInelastic(const G4String& modelName)
38 : G4HadronicInteraction(modelName)
39{
40 cache = 0.0;
41 particleDef[0] = G4PionZero::PionZero();
42 particleDef[1] = G4PionPlus::PionPlus();
43 particleDef[2] = G4PionMinus::PionMinus();
44 particleDef[3] = G4KaonPlus::KaonPlus();
45 particleDef[4] = G4KaonMinus::KaonMinus();
46 particleDef[5] = G4KaonZero::KaonZero();
47 particleDef[6] = G4AntiKaonZero::AntiKaonZero();
48 particleDef[7] = G4Proton::Proton();
49 particleDef[8] = G4Neutron::Neutron();
50 particleDef[9] = G4Lambda::Lambda();
51 particleDef[10] = G4SigmaPlus::SigmaPlus();
52 particleDef[11] = G4SigmaZero::SigmaZero();
53 particleDef[12] = G4SigmaMinus::SigmaMinus();
54 particleDef[13] = G4XiZero::XiZero();
55 particleDef[14] = G4XiMinus::XiMinus();
56 particleDef[15] = G4OmegaMinus::OmegaMinus();
57 particleDef[16] = G4AntiProton::AntiProton();
58 particleDef[17] = G4AntiNeutron::AntiNeutron();
59}
60
61
62G4double G4RPGInelastic::Pmltpc(G4int np, G4int nm, G4int nz,
63 G4int n, G4double b, G4double c)
64{
65 const G4double expxu = 82.; // upper bound for arg. of exp
66 const G4double expxl = -expxu; // lower bound for arg. of exp
67 G4double npf = 0.0;
68 G4double nmf = 0.0;
69 G4double nzf = 0.0;
70 G4int i;
71 for( i=2; i<=np; i++ )npf += std::log((double)i);
72 for( i=2; i<=nm; i++ )nmf += std::log((double)i);
73 for( i=2; i<=nz; i++ )nzf += std::log((double)i);
74 G4double r;
75 r = std::min( expxu, std::max( expxl, -(np-nm+nz+b)*(np-nm+nz+b)/(2*c*c*n*n)-npf-nmf-nzf ) );
76 return std::exp(r);
77}
78
79
80G4int G4RPGInelastic::Factorial( G4int n )
81{
82 G4int m = std::min(n,10);
83 G4int result = 1;
84 if( m <= 1 )return result;
85 for( G4int i=2; i<=m; ++i )result *= i;
86 return result;
87}
88
89
90G4bool G4RPGInelastic::MarkLeadingStrangeParticle(
91 const G4ReactionProduct &currentParticle,
92 const G4ReactionProduct &targetParticle,
93 G4ReactionProduct &leadParticle )
94{
95 // The following was in GenerateXandPt and TwoCluster.
96 // Add a parameter to the GenerateXandPt function telling it about the
97 // strange particle.
98 //
99 // Assumes that the original particle was a strange particle
100 //
101 G4bool lead = false;
102 if( (currentParticle.GetMass() >= G4KaonPlus::KaonPlus()->GetPDGMass()) &&
103 (currentParticle.GetDefinition() != G4Proton::Proton()) &&
104 (currentParticle.GetDefinition() != G4Neutron::Neutron()) )
105 {
106 lead = true;
107 leadParticle = currentParticle; // set lead to the incident particle
108 }
109 else if( (targetParticle.GetMass() >= G4KaonPlus::KaonPlus()->GetPDGMass()) &&
110 (targetParticle.GetDefinition() != G4Proton::Proton()) &&
111 (targetParticle.GetDefinition() != G4Neutron::Neutron()) )
112 {
113 lead = true;
114 leadParticle = targetParticle; // set lead to the target particle
115 }
116 return lead;
117}
118
119
120 void G4RPGInelastic::SetUpPions(const G4int np, const G4int nm,
121 const G4int nz,
122 G4FastVector<G4ReactionProduct,256> &vec,
123 G4int &vecLen)
124 {
125 if( np+nm+nz == 0 )return;
126 G4int i;
127 G4ReactionProduct *p;
128 for( i=0; i<np; ++i )
129 {
130 p = new G4ReactionProduct;
131 p->SetDefinition( G4PionPlus::PionPlus() );
132 (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
133 vec.SetElement( vecLen++, p );
134 }
135 for( i=np; i<np+nm; ++i )
136 {
137 p = new G4ReactionProduct;
138 p->SetDefinition( G4PionMinus::PionMinus() );
139 (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
140 vec.SetElement( vecLen++, p );
141 }
142 for( i=np+nm; i<np+nm+nz; ++i )
143 {
144 p = new G4ReactionProduct;
145 p->SetDefinition( G4PionZero::PionZero() );
146 (G4UniformRand() < 0.5) ? p->SetSide( -1 ) : p->SetSide( 1 );
147 vec.SetElement( vecLen++, p );
148 }
149 }
150
151
152 void G4RPGInelastic::GetNormalizationConstant(
153 const G4double energy, // MeV, <0 means annihilation channels
154 G4double &n,
155 G4double &anpn )
156 {
157 const G4double expxu = 82.; // upper bound for arg. of exp
158 const G4double expxl = -expxu; // lower bound for arg. of exp
159 const G4int numSec = 60;
160 //
161 // the only difference between the calculation for annihilation channels
162 // and normal is the starting value, iBegin, for the loop below
163 //
164 G4int iBegin = 1;
165 G4double en = energy;
166 if( energy < 0.0 )
167 {
168 iBegin = 2;
169 en *= -1.0;
170 }
171 //
172 // number of total particles vs. centre of mass Energy - 2*proton mass
173 //
174 G4double aleab = std::log(en/GeV);
175 n = 3.62567 + aleab*(0.665843 + aleab*(0.336514 + aleab*(0.117712 + 0.0136912*aleab)));
176 n -= 2.0;
177 //
178 // normalization constant for kno-distribution
179 //
180 anpn = 0.0;
181 G4double test, temp;
182 for( G4int i=iBegin; i<=numSec; ++i )
183 {
184 temp = pi*i/(2.0*n*n);
185 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(i*i)/(n*n) ) ) );
186 if( temp < 1.0 )
187 {
188 if( test >= 1.0e-10 )anpn += temp*test;
189 }
190 else
191 anpn += temp*test;
192 }
193 }
194
195void
196G4RPGInelastic::CalculateMomenta(G4FastVector<G4ReactionProduct,256>& vec,
197 G4int& vecLen,
198 const G4HadProjectile* originalIncident,
199 const G4DynamicParticle* originalTarget,
200 G4ReactionProduct& modifiedOriginal,
201 G4Nucleus& targetNucleus,
202 G4ReactionProduct& currentParticle,
203 G4ReactionProduct& targetParticle,
204 G4bool& incidentHasChanged,
205 G4bool& targetHasChanged,
206 G4bool quasiElastic)
207{
208 cache = 0;
209 what = originalIncident->Get4Momentum().vect();
210
211 G4ReactionProduct leadingStrangeParticle;
212
213 // strangeProduction.ReactionStage(originalIncident, modifiedOriginal,
214 // incidentHasChanged, originalTarget,
215 // targetParticle, targetHasChanged,
216 // targetNucleus, currentParticle,
217 // vec, vecLen,
218 // false, leadingStrangeParticle);
219
220 if( quasiElastic )
221 {
222 twoBody.ReactionStage(originalIncident, modifiedOriginal,
223 incidentHasChanged, originalTarget,
224 targetParticle, targetHasChanged,
225 targetNucleus, currentParticle,
226 vec, vecLen,
227 false, leadingStrangeParticle);
228 return;
229 }
230
231 G4bool leadFlag = MarkLeadingStrangeParticle(currentParticle,
232 targetParticle,
233 leadingStrangeParticle );
234 //
235 // Note: the number of secondaries can be reduced in GenerateXandPt
236 // and TwoCluster
237 //
238 G4bool finishedGenXPt = false;
239 G4bool annihilation = false;
240 if( originalIncident->GetDefinition()->GetPDGEncoding() < 0 &&
241 currentParticle.GetMass() == 0.0 && targetParticle.GetMass() == 0.0 )
242 {
243 // original was an anti-particle and annihilation has taken place
244 annihilation = true;
245 G4double ekcor = 1.0;
246 G4double ek = originalIncident->GetKineticEnergy();
247 G4double ekOrg = ek;
248
249 const G4double tarmas = originalTarget->GetDefinition()->GetPDGMass();
250 if( ek > 1.0*GeV )ekcor = 1./(ek/GeV);
251 const G4double atomicWeight = targetNucleus.GetN();
252 ek = 2*tarmas + ek*(1.+ekcor/atomicWeight);
253 G4double tkin = targetNucleus.Cinema(ek);
254 ek += tkin;
255 ekOrg += tkin;
256 // modifiedOriginal.SetKineticEnergy( ekOrg );
257 //
258 // evaporation -- re-calculate black track energies
259 // this was Done already just before the cascade
260 //
261 tkin = targetNucleus.AnnihilationEvaporationEffects(ek, ekOrg);
262 ekOrg -= tkin;
263 ekOrg = std::max( 0.0001*GeV, ekOrg );
264 modifiedOriginal.SetKineticEnergy( ekOrg );
265 G4double amas = originalIncident->GetDefinition()->GetPDGMass();
266 G4double et = ekOrg + amas;
267 G4double p = std::sqrt( std::abs(et*et-amas*amas) );
268 G4double pp = modifiedOriginal.GetMomentum().mag();
269 if( pp > 0.0 )
270 {
271 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
272 modifiedOriginal.SetMomentum( momentum * (p/pp) );
273 }
274 if( ekOrg <= 0.0001 )
275 {
276 modifiedOriginal.SetKineticEnergy( 0.0 );
277 modifiedOriginal.SetMomentum( 0.0, 0.0, 0.0 );
278 }
279 }
280
281 // twsup gives percentage of time two-cluster model is called
282
283 const G4double twsup[] = { 1.0, 0.7, 0.5, 0.3, 0.2, 0.1 };
284 G4double rand1 = G4UniformRand();
285 G4double rand2 = G4UniformRand();
286
287 // Cache current, target, and secondaries
288 G4ReactionProduct saveCurrent = currentParticle;
289 G4ReactionProduct saveTarget = targetParticle;
290 std::vector<G4ReactionProduct> savevec;
291 for (G4int i = 0; i < vecLen; i++) savevec.push_back(*vec[i]);
292
293 // Call fragmentation code if
294 // 1) there is annihilation, or
295 // 2) there are more than 5 secondaries, or
296 // 3) incident KE is > 1 GeV AND
297 // ( incident is a kaon AND rand < 0.5 OR twsup )
298 //
299
300 if( annihilation || vecLen > 5 ||
301 ( modifiedOriginal.GetKineticEnergy()/GeV >= 1.0 &&
302
303 (((originalIncident->GetDefinition() == G4KaonPlus::KaonPlus() ||
304 originalIncident->GetDefinition() == G4KaonMinus::KaonMinus() ||
305 originalIncident->GetDefinition() == G4KaonZeroLong::KaonZeroLong() ||
306 originalIncident->GetDefinition() == G4KaonZeroShort::KaonZeroShort()) &&
307 rand1 < 0.5)
308 || rand2 > twsup[vecLen]) ) )
309
310 finishedGenXPt =
311 fragmentation.ReactionStage(originalIncident, modifiedOriginal,
312 incidentHasChanged, originalTarget,
313 targetParticle, targetHasChanged,
314 targetNucleus, currentParticle,
315 vec, vecLen,
316 leadFlag, leadingStrangeParticle);
317
318 if (finishedGenXPt) return;
319
320 G4bool finishedTwoClu = false;
321
322 if (modifiedOriginal.GetTotalMomentum() < 1.0) {
323 for (G4int i = 0; i < vecLen; i++) delete vec[i];
324 vecLen = 0;
325
326 } else {
327 // Occaisionally, GenerateXandPt will fail in the annihilation channel.
328 // Restore current, target and secondaries to pre-GenerateXandPt state
329 // before trying annihilation in TwoCluster
330
331 if (!finishedGenXPt && annihilation) {
332 currentParticle = saveCurrent;
333 targetParticle = saveTarget;
334 for (G4int i = 0; i < vecLen; i++) delete vec[i];
335 vecLen = 0;
336 vec.Initialize( 0 );
337 for (G4int i = 0; i < G4int(savevec.size()); i++) {
338 G4ReactionProduct* p = new G4ReactionProduct;
339 *p = savevec[i];
340 vec.SetElement( vecLen++, p );
341 }
342 }
343
344 // Big violations of energy conservation in this method - don't use
345 //
346 // pionSuppression.ReactionStage(originalIncident, modifiedOriginal,
347 // incidentHasChanged, originalTarget,
348 // targetParticle, targetHasChanged,
349 // targetNucleus, currentParticle,
350 // vec, vecLen,
351 // false, leadingStrangeParticle);
352
353 try
354 {
355 finishedTwoClu =
356 twoCluster.ReactionStage(originalIncident, modifiedOriginal,
357 incidentHasChanged, originalTarget,
358 targetParticle, targetHasChanged,
359 targetNucleus, currentParticle,
360 vec, vecLen,
361 leadFlag, leadingStrangeParticle);
362 }
363 catch(G4HadReentrentException aC)
364 {
365 aC.Report(G4cout);
366 throw G4HadReentrentException(__FILE__, __LINE__, "Failing to calculate momenta");
367 }
368 }
369
370 if (finishedTwoClu) return;
371
372 twoBody.ReactionStage(originalIncident, modifiedOriginal,
373 incidentHasChanged, originalTarget,
374 targetParticle, targetHasChanged,
375 targetNucleus, currentParticle,
376 vec, vecLen,
377 false, leadingStrangeParticle);
378}
379
380/*
381 void G4RPGInelastic::
382 Rotate(G4FastVector<G4ReactionProduct,256> &vec, G4int &vecLen)
383 {
384 G4double rotation = 2.*pi*G4UniformRand();
385 cache = rotation;
386 G4int i;
387 for( i=0; i<vecLen; ++i )
388 {
389 G4ThreeVector momentum = vec[i]->GetMomentum();
390 momentum = momentum.rotate(rotation, what);
391 vec[i]->SetMomentum(momentum);
392 }
393 }
394*/
395
396void
397G4RPGInelastic::SetUpChange(G4FastVector<G4ReactionProduct,256>& vec,
398 G4int& vecLen,
399 G4ReactionProduct& currentParticle,
400 G4ReactionProduct& targetParticle,
401 G4bool& incidentHasChanged )
402{
403 theParticleChange.Clear();
404 G4ParticleDefinition* aKaonZL = G4KaonZeroLong::KaonZeroLong();
405 G4ParticleDefinition* aKaonZS = G4KaonZeroShort::KaonZeroShort();
406 G4int i;
407
408 if (currentParticle.GetDefinition() == particleDef[k0]) {
409 if (G4UniformRand() < 0.5) {
410 currentParticle.SetDefinitionAndUpdateE(aKaonZL);
411 incidentHasChanged = true;
412 } else {
413 currentParticle.SetDefinitionAndUpdateE(aKaonZS);
414 }
415 } else if (currentParticle.GetDefinition() == particleDef[k0b]) {
416 if (G4UniformRand() < 0.5) {
417 currentParticle.SetDefinitionAndUpdateE(aKaonZL);
418 } else {
419 currentParticle.SetDefinitionAndUpdateE(aKaonZS);
420 incidentHasChanged = true;
421 }
422 }
423
424 if (targetParticle.GetDefinition() == particleDef[k0] ||
425 targetParticle.GetDefinition() == particleDef[k0b] ) {
426 if (G4UniformRand() < 0.5) {
427 targetParticle.SetDefinitionAndUpdateE(aKaonZL);
428 } else {
429 targetParticle.SetDefinitionAndUpdateE(aKaonZS);
430 }
431 }
432
433 for (i = 0; i < vecLen; ++i) {
434 if (vec[i]->GetDefinition() == particleDef[k0] ||
435 vec[i]->GetDefinition() == particleDef[k0b] ) {
436 if (G4UniformRand() < 0.5) {
437 vec[i]->SetDefinitionAndUpdateE(aKaonZL);
438 } else {
439 vec[i]->SetDefinitionAndUpdateE(aKaonZS);
440 }
441 }
442 }
443
444 if (incidentHasChanged) {
445 G4DynamicParticle* p0 = new G4DynamicParticle;
446 p0->SetDefinition(currentParticle.GetDefinition() );
447 p0->SetMomentum(currentParticle.GetMomentum() );
448 theParticleChange.AddSecondary( p0 );
449 theParticleChange.SetStatusChange( stopAndKill );
450 theParticleChange.SetEnergyChange( 0.0 );
451
452 } else {
453 G4double p = currentParticle.GetMomentum().mag()/MeV;
454 G4ThreeVector m = currentParticle.GetMomentum();
455 if (p > DBL_MIN)
456 theParticleChange.SetMomentumChange( m.x()/p, m.y()/p, m.z()/p );
457 else
458 theParticleChange.SetMomentumChange( 0.0, 0.0, 1.0 );
459
460 G4double aE = currentParticle.GetKineticEnergy();
461 if (std::fabs(aE)<.1*eV) aE=.1*eV;
462 theParticleChange.SetEnergyChange( aE );
463 }
464
465 if (targetParticle.GetMass() > 0.0) // Tgt particle can be eliminated in TwoBody
466 {
467 G4ThreeVector momentum = targetParticle.GetMomentum();
468 momentum = momentum.rotate(cache, what);
469 G4double targKE = targetParticle.GetKineticEnergy();
470 G4ThreeVector dir(0.0, 0.0, 1.0);
471 if (targKE < DBL_MIN)
472 targKE = DBL_MIN;
473 else
474 dir = momentum/momentum.mag();
475
476 G4DynamicParticle* p1 =
477 new G4DynamicParticle(targetParticle.GetDefinition(), dir, targKE);
478
479 theParticleChange.AddSecondary( p1 );
480 }
481
482 G4DynamicParticle* p;
483 for (i = 0; i < vecLen; ++i) {
484 G4double secKE = vec[i]->GetKineticEnergy();
485 G4ThreeVector momentum = vec[i]->GetMomentum();
486 G4ThreeVector dir(0.0, 0.0, 1.0);
487 if (secKE < DBL_MIN)
488 secKE = DBL_MIN;
489 else
490 dir = momentum/momentum.mag();
491
492 p = new G4DynamicParticle(vec[i]->GetDefinition(), dir, secKE);
493 theParticleChange.AddSecondary( p );
494 delete vec[i];
495 }
496}
497
498
499std::pair<G4int, G4double>
500G4RPGInelastic::interpolateEnergy(G4double e) const
501{
502 G4int index = 29;
503 G4double fraction = 0.0;
504
505 for (G4int i = 1; i < 30; i++) {
506 if (e < energyScale[i]) {
507 index = i-1;
508 fraction = (e - energyScale[index]) / (energyScale[i] - energyScale[index]);
509 break;
510 }
511 }
512 return std::pair<G4int, G4double>(index, fraction);
513}
514
515
516G4int
517G4RPGInelastic::sampleFlat(std::vector<G4double> sigma) const
518{
519 G4int i;
520 G4double sum(0.);
521 for (i = 0; i < G4int(sigma.size()); i++) sum += sigma[i];
522
523 G4double fsum = sum*G4UniformRand();
524 G4double partialSum = 0.0;
525 G4int channel = 0;
526
527 for (i = 0; i < G4int(sigma.size()); i++) {
528 partialSum += sigma[i];
529 if (fsum < partialSum) {
530 channel = i;
531 break;
532 }
533 }
534
535 return channel;
536}
537
538
539void G4RPGInelastic::CheckQnums(G4FastVector<G4ReactionProduct,256> &vec,
540 G4int &vecLen,
541 G4ReactionProduct &currentParticle,
542 G4ReactionProduct &targetParticle,
543 G4double Q, G4double B, G4double S)
544{
545 G4ParticleDefinition* projDef = currentParticle.GetDefinition();
546 G4ParticleDefinition* targDef = targetParticle.GetDefinition();
547 G4double chargeSum = projDef->GetPDGCharge() + targDef->GetPDGCharge();
548 G4double baryonSum = projDef->GetBaryonNumber() + targDef->GetBaryonNumber();
549 G4double strangenessSum = projDef->GetQuarkContent(3) -
550 projDef->GetAntiQuarkContent(3) +
551 targDef->GetQuarkContent(3) -
552 targDef->GetAntiQuarkContent(3);
553
554 G4ParticleDefinition* secDef = 0;
555 for (G4int i = 0; i < vecLen; i++) {
556 secDef = vec[i]->GetDefinition();
557 chargeSum += secDef->GetPDGCharge();
558 baryonSum += secDef->GetBaryonNumber();
559 strangenessSum += secDef->GetQuarkContent(3)
560 - secDef->GetAntiQuarkContent(3);
561 }
562
563 G4bool OK = true;
564 if (chargeSum != Q) {
565 G4cout << " Charge not conserved " << G4endl;
566 OK = false;
567 }
568 if (baryonSum != B) {
569 G4cout << " Baryon number not conserved " << G4endl;
570 OK = false;
571 }
572 if (strangenessSum != S) {
573 G4cout << " Strangeness not conserved " << G4endl;
574 OK = false;
575 }
576
577 if (!OK) {
578 G4cout << " projectile: " << projDef->GetParticleName()
579 << " target: " << targDef->GetParticleName() << G4endl;
580 for (G4int i = 0; i < vecLen; i++) {
581 secDef = vec[i]->GetDefinition();
582 G4cout << secDef->GetParticleName() << " " ;
583 }
584 G4cout << G4endl;
585 }
586
587}
588
589
590const G4double G4RPGInelastic::energyScale[30] = {
591 0.0, 0.01, 0.013, 0.018, 0.024, 0.032, 0.042, 0.056, 0.075, 0.1,
592 0.13, 0.18, 0.24, 0.32, 0.42, 0.56, 0.75, 1.0, 1.3, 1.8,
593 2.4, 3.2, 4.2, 5.6, 7.5, 10.0, 13.0, 18.0, 24.0, 32.0 };
594
595
596/* end of file */
Note: See TracBrowser for help on using the repository browser.