source: trunk/source/processes/hadronic/stopping/src/G4KaonMinusAbsorption.cc@ 1036

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

update to geant4.9.2
File size: 15.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// G4KaonMinusAbsorption physics process
27// Larry Felawka (TRIUMF), April 1998
28//---------------------------------------------------------------------
29
30#include "G4KaonMinusAbsorption.hh"
31#include "G4DynamicParticle.hh"
32#include "G4ParticleTypes.hh"
33#include "Randomize.hh"
34#include <string.h>
35#include <cmath>
36#include <stdio.h>
37
38#define MAX_SECONDARIES 100
39
40// constructor
41
42G4KaonMinusAbsorption::G4KaonMinusAbsorption(const G4String& processName,
43 G4ProcessType aType ) :
44 G4VRestProcess (processName, aType), // initialization
45 massKaonMinus(G4KaonMinus::KaonMinus()->GetPDGMass()/GeV),
46 massGamma(G4Gamma::Gamma()->GetPDGMass()/GeV),
47 massPionZero(G4PionZero::PionZero()->GetPDGMass()/GeV),
48 massProton(G4Proton::Proton()->GetPDGMass()/GeV),
49 massLambda(G4Lambda::Lambda()->GetPDGMass()/GeV),
50 pdefKaonMinus(G4KaonMinus::KaonMinus()),
51 pdefGamma(G4Gamma::Gamma()),
52 pdefPionZero(G4PionZero::PionZero()),
53 pdefProton(G4Proton::Proton()),
54 pdefNeutron(G4Neutron::Neutron()),
55 pdefLambda(G4Lambda::Lambda()),
56 pdefDeuteron(G4Deuteron::Deuteron()),
57 pdefTriton(G4Triton::Triton()),
58 pdefAlpha(G4Alpha::Alpha())
59{
60 if (verboseLevel>0) {
61 G4cout << GetProcessName() << " is created "<< G4endl;
62 }
63 SetProcessSubType(fHadronAtRest);
64 pv = new G4GHEKinematicsVector [MAX_SECONDARIES+1];
65 eve = new G4GHEKinematicsVector [MAX_SECONDARIES];
66 gkin = new G4GHEKinematicsVector [MAX_SECONDARIES];
67
68}
69
70// destructor
71
72G4KaonMinusAbsorption::~G4KaonMinusAbsorption()
73{
74 delete [] pv;
75 delete [] eve;
76 delete [] gkin;
77}
78
79
80// methods.............................................................................
81
82G4bool G4KaonMinusAbsorption::IsApplicable(
83 const G4ParticleDefinition& particle
84 )
85{
86 return ( &particle == pdefKaonMinus );
87
88}
89
90// Warning - this method may be optimized away if made "inline"
91G4int G4KaonMinusAbsorption::GetNumberOfSecondaries()
92{
93 return ( ngkine );
94
95}
96
97// Warning - this method may be optimized away if made "inline"
98G4GHEKinematicsVector* G4KaonMinusAbsorption::GetSecondaryKinematics()
99{
100 return ( &gkin[0] );
101
102}
103
104G4double G4KaonMinusAbsorption::AtRestGetPhysicalInteractionLength(
105 const G4Track& track,
106 G4ForceCondition* condition
107 )
108{
109 // beggining of tracking
110 ResetNumberOfInteractionLengthLeft();
111
112 // condition is set to "Not Forced"
113 *condition = NotForced;
114
115 // get mean life time
116 currentInteractionLength = GetMeanLifeTime(track, condition);
117
118 if ((currentInteractionLength <0.0) || (verboseLevel>2)){
119 G4cout << "G4KaonMinusAbsorptionProcess::AtRestGetPhysicalInteractionLength ";
120 G4cout << "[ " << GetProcessName() << "]" <<G4endl;
121 track.GetDynamicParticle()->DumpInfo();
122 G4cout << " in Material " << track.GetMaterial()->GetName() <<G4endl;
123 G4cout << "MeanLifeTime = " << currentInteractionLength/ns << "[ns]" <<G4endl;
124 }
125
126 return theNumberOfInteractionLengthLeft * currentInteractionLength;
127
128}
129
130G4VParticleChange* G4KaonMinusAbsorption::AtRestDoIt(
131 const G4Track& track,
132 const G4Step&
133 )
134//
135// Handles KaonMinus at rest; a KaonMinus can either create secondaries or
136// do nothing (in which case it should be sent back to decay-handling
137// section
138//
139{
140
141// Initialize ParticleChange
142// all members of G4VParticleChange are set to equal to
143// corresponding member in G4Track
144
145 aParticleChange.Initialize(track);
146
147// Store some global quantities that depend on current material and particle
148
149 globalTime = track.GetGlobalTime()/s;
150 G4Material * aMaterial = track.GetMaterial();
151 const G4int numberOfElements = aMaterial->GetNumberOfElements();
152 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
153
154 const G4double* theAtomicNumberDensity = aMaterial->GetAtomicNumDensityVector();
155 G4double normalization = 0;
156 for ( G4int i1=0; i1 < numberOfElements; i1++ )
157 {
158 normalization += theAtomicNumberDensity[i1] ; // change when nucleon specific
159 // probabilities are included.
160 }
161 G4double runningSum= 0.;
162 G4double random = G4UniformRand()*normalization;
163 for ( G4int i2=0; i2 < numberOfElements; i2++ )
164 {
165 runningSum += theAtomicNumberDensity[i2]; // change when nucleon specific
166 // probabilities are included.
167 if (random<=runningSum)
168 {
169 targetCharge = G4double( ((*theElementVector)[i2])->GetZ());
170 targetAtomicMass = (*theElementVector)[i2]->GetN();
171 }
172 }
173 if (random>runningSum)
174 {
175 targetCharge = G4double((*theElementVector)[numberOfElements-1]->GetZ());
176 targetAtomicMass = (*theElementVector)[numberOfElements-1]->GetN();
177
178 }
179
180 if (verboseLevel>1) {
181 G4cout << "G4KaonMinusAbsorption::AtRestDoIt is invoked " <<G4endl;
182 }
183
184 G4ParticleMomentum momentum;
185 G4float localtime;
186
187 G4ThreeVector position = track.GetPosition();
188
189 GenerateSecondaries(); // Generate secondaries
190
191 aParticleChange.SetNumberOfSecondaries( ngkine );
192
193 for ( G4int isec = 0; isec < ngkine; isec++ ) {
194 G4DynamicParticle* aNewParticle = new G4DynamicParticle;
195 aNewParticle->SetDefinition( gkin[isec].GetParticleDef() );
196 aNewParticle->SetMomentum( gkin[isec].GetMomentum() * GeV );
197
198 localtime = globalTime + gkin[isec].GetTOF();
199
200 G4Track* aNewTrack = new G4Track( aNewParticle, localtime*s, position );
201 aNewTrack->SetTouchableHandle(track.GetTouchableHandle());
202 aParticleChange.AddSecondary( aNewTrack );
203
204 }
205
206 aParticleChange.ProposeLocalEnergyDeposit( 0.0*GeV );
207
208 aParticleChange.ProposeTrackStatus(fStopAndKill); // Kill the incident KaonMinus
209
210// clear InteractionLengthLeft
211
212 ResetNumberOfInteractionLengthLeft();
213
214 return &aParticleChange;
215
216}
217
218
219void G4KaonMinusAbsorption::GenerateSecondaries()
220{
221 static G4int index;
222 static G4int l;
223 static G4int nopt;
224 static G4int i;
225 static G4ParticleDefinition* jnd;
226
227 for (i = 1; i <= MAX_SECONDARIES; ++i) {
228 pv[i].SetZero();
229 }
230
231 ngkine = 0; // number of generated secondary particles
232 ntot = 0;
233 result.SetZero();
234 result.SetMass( massKaonMinus );
235 result.SetKineticEnergyAndUpdate( 0. );
236 result.SetTOF( 0. );
237 result.SetParticleDef( pdefKaonMinus );
238
239 KaonMinusAbsorption(&nopt);
240
241 // *** CHECK WHETHER THERE ARE NEW PARTICLES GENERATED ***
242 if (ntot != 0 || result.GetParticleDef() != pdefKaonMinus) {
243 // *** CURRENT PARTICLE IS NOT THE SAME AS IN THE BEGINNING OR/AND ***
244 // *** ONE OR MORE SECONDARIES HAVE BEEN GENERATED ***
245
246 // --- INITIAL PARTICLE TYPE HAS BEEN CHANGED ==> PUT NEW TYPE ON ---
247 // --- THE GEANT TEMPORARY STACK ---
248
249 // --- PUT PARTICLE ON THE STACK ---
250 gkin[0] = result;
251 gkin[0].SetTOF( result.GetTOF() * 5e-11 );
252 ngkine = 1;
253
254 // --- ALL QUANTITIES ARE TAKEN FROM THE GHEISHA STACK WHERE THE ---
255 // --- CONVENTION IS THE FOLLOWING ---
256
257 // --- ONE OR MORE SECONDARIES HAVE BEEN GENERATED ---
258 for (l = 1; l <= ntot; ++l) {
259 index = l - 1;
260 jnd = eve[index].GetParticleDef();
261
262 // --- ADD PARTICLE TO THE STACK IF STACK NOT YET FULL ---
263 if (ngkine < MAX_SECONDARIES) {
264 gkin[ngkine] = eve[index];
265 gkin[ngkine].SetTOF( eve[index].GetTOF() * 5e-11 );
266 ++ngkine;
267 }
268 }
269 }
270 else {
271 // --- NO SECONDARIES GENERATED AND PARTICLE IS STILL THE SAME ---
272 // --- ==> COPY EVERYTHING BACK IN THE CURRENT GEANT STACK ---
273 ngkine = 0;
274 ntot = 0;
275 globalTime += result.GetTOF() * G4float(5e-11);
276 }
277
278 // --- LIMIT THE VALUE OF NGKINE IN CASE OF OVERFLOW ---
279 ngkine = G4int(std::min(ngkine,G4int(MAX_SECONDARIES)));
280
281} // GenerateSecondaries
282
283
284void G4KaonMinusAbsorption::Poisso(G4float xav, G4int *iran)
285{
286 static G4int i;
287 static G4float r, p1, p2, p3;
288 static G4int mm;
289 static G4float rr, ran, rrr, ran1;
290
291 // *** GENERATION OF POISSON DISTRIBUTION ***
292 // *** NVE 16-MAR-1988 CERN GENEVA ***
293 // ORIGIN : H.FESEFELDT (27-OCT-1983)
294
295 // --- USE NORMAL DISTRIBUTION FOR <X> > 9.9 ---
296 if (xav > G4float(9.9)) {
297 // ** NORMAL DISTRIBUTION WITH SIGMA**2 = <X>
298 Normal(&ran1);
299 ran1 = xav + ran1 * std::sqrt(xav);
300 *iran = G4int(ran1);
301 if (*iran < 0) {
302 *iran = 0;
303 }
304 }
305 else {
306 mm = G4int(xav * G4float(5.));
307 *iran = 0;
308 if (mm > 0) {
309 r = std::exp(-G4double(xav));
310 ran1 = G4UniformRand();
311 if (ran1 > r) {
312 rr = r;
313 for (i = 1; i <= mm; ++i) {
314 ++(*iran);
315 if (i <= 5) {
316 rrr = std::pow(xav, G4float(i)) / NFac(i);
317 }
318 // ** STIRLING' S FORMULA FOR LARGE NUMBERS
319 if (i > 5) {
320 rrr = std::exp(i * std::log(xav) -
321 (i + G4float(.5)) * std::log(i * G4float(1.)) +
322 i - G4float(.9189385));
323 }
324 rr += r * rrr;
325 if (ran1 <= rr) {
326 break;
327 }
328 }
329 }
330 }
331 else {
332 // ** FOR VERY SMALL XAV TRY IRAN=1,2,3
333 p1 = xav * std::exp(-G4double(xav));
334 p2 = xav * p1 / G4float(2.);
335 p3 = xav * p2 / G4float(3.);
336 ran = G4UniformRand();
337 if (ran >= p3) {
338 if (ran >= p2) {
339 if (ran >= p1) {
340 *iran = 0;
341 }
342 else {
343 *iran = 1;
344 }
345 }
346 else {
347 *iran = 2;
348 }
349 }
350 else {
351 *iran = 3;
352 }
353 }
354 }
355
356} // Poisso
357
358
359G4int G4KaonMinusAbsorption::NFac(G4int n)
360{
361 G4int ret_val;
362
363 static G4int i, m;
364
365 // *** NVE 16-MAR-1988 CERN GENEVA ***
366 // ORIGIN : H.FESEFELDT (27-OCT-1983)
367
368 ret_val = 1;
369 m = n;
370 if (m > 1) {
371 if (m > 10) {
372 m = 10;
373 }
374 for (i = 2; i <= m; ++i) {
375 ret_val *= i;
376 }
377 }
378 return ret_val;
379
380} // NFac
381
382
383void G4KaonMinusAbsorption::Normal(G4float *ran)
384{
385 static G4int i;
386
387 // *** NVE 14-APR-1988 CERN GENEVA ***
388 // ORIGIN : H.FESEFELDT (27-OCT-1983)
389
390 *ran = G4float(-6.);
391 for (i = 1; i <= 12; ++i) {
392 *ran += G4UniformRand();
393 }
394
395} // Normal
396
397
398void G4KaonMinusAbsorption::KaonMinusAbsorption(G4int *nopt)
399{
400 static G4int i;
401 static G4int nt, nbl;
402 static G4float ran, pcm;
403 static G4int isw;
404 static G4float tex;
405 static G4float ran2, tof1, ekin, ekin1, ekin2, black;
406 static G4float pnrat;
407 static G4ParticleDefinition* ipa1;
408 static G4ParticleDefinition* inve;
409
410 // *** CHARGED KAON ABSORPTION BY A NUCLEUS ***
411 // *** NVE 04-MAR-1988 CERN GENEVA ***
412 // ORIGIN : H.FESEFELDT (09-JULY-1987)
413
414 // PRODUCTION OF A HYPERFRAGMENT WITH SUBSEQUENT DECAY
415 // PANOFSKY RATIO (K- P --> LAMBDA PI0/K- P --> LAMBDA GAMMA) = 3/2
416
417 pv[1].SetZero();
418 pv[1].SetMass( massKaonMinus );
419 pv[1].SetKineticEnergyAndUpdate( 0. );
420 pv[1].SetTOF( result.GetTOF() );
421 pv[1].SetParticleDef( result.GetParticleDef() );
422 if (targetAtomicMass <= G4float(1.5)) {
423 ran = G4UniformRand();
424 tof1 = std::log(ran) * G4float(-12.5);
425 tof1 *= G4float(20.);
426 ran = G4UniformRand();
427 isw = 1;
428 if (ran < G4float(.33)) {
429 isw = 2;
430 }
431 *nopt = isw;
432 pv[3].SetZero();
433 pv[3].SetMass( massLambda );
434 pv[3].SetKineticEnergyAndUpdate( 0. );
435 pv[3].SetTOF( result.GetTOF() + tof1 );
436 pv[3].SetParticleDef( pdefLambda );
437 pcm = massKaonMinus + massProton - massLambda;
438 if (isw != 1) {
439 pv[2].SetZero();
440 pv[2].SetMass( massGamma );
441 pv[2].SetKineticEnergyAndUpdate( pcm );
442 pv[2].SetTOF( result.GetTOF() + tof1 );
443 pv[2].SetParticleDef( pdefGamma );
444 }
445 else {
446 pcm = pcm * pcm - massPionZero * massPionZero;
447 if (pcm <= G4float(0.)) {
448 pcm = G4float(0.);
449 }
450 pv[2].SetZero();
451 pv[2].SetEnergy( std::sqrt(pcm + massPionZero * massPionZero) );
452 pv[2].SetMassAndUpdate( massPionZero );
453 pv[2].SetTOF( result.GetTOF() + tof1 );
454 pv[2].SetParticleDef( pdefPionZero );
455 }
456 result = pv[2];
457 if (ntot < MAX_SECONDARIES-1) {
458 eve[ntot++] = pv[3];
459 }
460 }
461 else {
462 // **
463 // ** STAR PRODUCTION FOR PION ABSORPTION IN HEAVY ELEMENTS
464 // **
465 evapEnergy1 = G4float(.3);
466 evapEnergy3 = G4float(.15);
467 nt = 1;
468 tex = evapEnergy1;
469 black = std::log(targetAtomicMass) * G4float(.5);
470 Poisso(black, &nbl);
471 if (nt + nbl > (MAX_SECONDARIES - 2)) {
472 nbl = (MAX_SECONDARIES - 2) - nt;
473 }
474 if (nbl <= 0) {
475 nbl = 1;
476 }
477 ekin = tex / nbl;
478 ekin2 = G4float(0.);
479 for (i = 1; i <= nbl; ++i) {
480 if (nt == (MAX_SECONDARIES - 2)) {
481 continue;
482 }
483 ran2 = G4UniformRand();
484 ekin1 = -G4double(ekin) * std::log(ran2);
485 ekin2 += ekin1;
486 ipa1 = pdefNeutron;
487 pnrat = G4float(1.) - targetCharge / targetAtomicMass;
488 if (G4UniformRand() > pnrat) {
489 ipa1 = pdefProton;
490 }
491 ++nt;
492 pv[nt].SetZero();
493 pv[nt].SetMass( ipa1->GetPDGMass()/GeV );
494 pv[nt].SetKineticEnergyAndUpdate( ekin1 );
495 pv[nt].SetTOF( 2. );
496 pv[nt].SetParticleDef( ipa1 );
497 if (ekin2 > tex) {
498 break;
499 }
500 }
501 tex = evapEnergy3;
502 black = std::log(targetAtomicMass) * G4float(.5);
503 Poisso(black, &nbl);
504 if (nt + nbl > (MAX_SECONDARIES - 2)) {
505 nbl = (MAX_SECONDARIES - 2) - nt;
506 }
507 if (nbl <= 0) {
508 nbl = 1;
509 }
510 ekin = tex / nbl;
511 ekin2 = G4float(0.);
512 for (i = 1; i <= nbl; ++i) {
513 if (nt == (MAX_SECONDARIES - 2)) {
514 continue;
515 }
516 ran2 = G4UniformRand();
517 ekin1 = -G4double(ekin) * std::log(ran2);
518 ekin2 += ekin1;
519 ++nt;
520 ran = G4UniformRand();
521 inve = pdefDeuteron;
522 if (ran > G4float(.6)) {
523 inve = pdefTriton;
524 }
525 if (ran > G4float(.9)) {
526 inve = pdefAlpha;
527 }
528// PV(5,NT)=(ABS(IPA(NT))-28)*RMASS(14) <-- Wrong! (LF)
529 pv[nt].SetZero();
530 pv[nt].SetMass( inve->GetPDGMass()/GeV );
531 pv[nt].SetKineticEnergyAndUpdate( ekin1 );
532 pv[nt].SetTOF( 2. );
533 pv[nt].SetParticleDef( inve );
534 if (ekin2 > tex) {
535 break;
536 }
537 }
538 // **
539 // ** STORE ON EVENT COMMON
540 // **
541 ran = G4UniformRand();
542 tof1 = std::log(ran) * G4float(-12.5);
543 tof1 *= G4float(20.);
544 for (i = 2; i <= nt; ++i) {
545 pv[i].SetTOF( result.GetTOF() + tof1 );
546 }
547 result = pv[2];
548 for (i = 3; i <= nt; ++i) {
549 if (ntot >= MAX_SECONDARIES) {
550 break;
551 }
552 eve[ntot++] = pv[i];
553 }
554 }
555
556} // KaonMinusAbsorption
Note: See TracBrowser for help on using the repository browser.