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

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

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