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

Last change on this file since 1330 was 1055, checked in by garnier, 17 years ago

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