00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030 #include <string.h>
00031 #include <cmath>
00032 #include <stdio.h>
00033
00034 #include "G4KaonMinusAbsorption.hh"
00035 #include "G4DynamicParticle.hh"
00036 #include "G4ParticleTypes.hh"
00037 #include "Randomize.hh"
00038 #include "G4SystemOfUnits.hh"
00039 #include "G4HadronicProcessStore.hh"
00040 #include "G4HadronicDeprecate.hh"
00041
00042 #define MAX_SECONDARIES 100
00043
00044
00045
00046 G4KaonMinusAbsorption::G4KaonMinusAbsorption(const G4String& processName,
00047 G4ProcessType aType ) :
00048 G4VRestProcess (processName, aType),
00049 massKaonMinus(G4KaonMinus::KaonMinus()->GetPDGMass()/GeV),
00050 massGamma(G4Gamma::Gamma()->GetPDGMass()/GeV),
00051 massPionZero(G4PionZero::PionZero()->GetPDGMass()/GeV),
00052 massProton(G4Proton::Proton()->GetPDGMass()/GeV),
00053 massLambda(G4Lambda::Lambda()->GetPDGMass()/GeV),
00054 pdefKaonMinus(G4KaonMinus::KaonMinus()),
00055 pdefGamma(G4Gamma::Gamma()),
00056 pdefPionZero(G4PionZero::PionZero()),
00057 pdefProton(G4Proton::Proton()),
00058 pdefNeutron(G4Neutron::Neutron()),
00059 pdefLambda(G4Lambda::Lambda()),
00060 pdefDeuteron(G4Deuteron::Deuteron()),
00061 pdefTriton(G4Triton::Triton()),
00062 pdefAlpha(G4Alpha::Alpha())
00063 {
00064 G4HadronicDeprecate("G4KaonMinusAbsorption");
00065 if (verboseLevel>0) {
00066 G4cout << GetProcessName() << " is created "<< G4endl;
00067 }
00068 SetProcessSubType(fHadronAtRest);
00069 pv = new G4GHEKinematicsVector [MAX_SECONDARIES+1];
00070 eve = new G4GHEKinematicsVector [MAX_SECONDARIES];
00071 gkin = new G4GHEKinematicsVector [MAX_SECONDARIES];
00072
00073 G4HadronicProcessStore::Instance()->RegisterExtraProcess(this);
00074 }
00075
00076
00077
00078 G4KaonMinusAbsorption::~G4KaonMinusAbsorption()
00079 {
00080 G4HadronicProcessStore::Instance()->DeRegisterExtraProcess(this);
00081 delete [] pv;
00082 delete [] eve;
00083 delete [] gkin;
00084 }
00085
00086 void G4KaonMinusAbsorption::PreparePhysicsTable(const G4ParticleDefinition& p)
00087 {
00088 G4HadronicProcessStore::Instance()->RegisterParticleForExtraProcess(this, &p);
00089 }
00090
00091 void G4KaonMinusAbsorption::BuildPhysicsTable(const G4ParticleDefinition& p)
00092 {
00093 G4HadronicProcessStore::Instance()->PrintInfo(&p);
00094 }
00095
00096
00097
00098 G4bool G4KaonMinusAbsorption::IsApplicable(
00099 const G4ParticleDefinition& particle
00100 )
00101 {
00102 return ( &particle == pdefKaonMinus );
00103
00104 }
00105
00106
00107 G4int G4KaonMinusAbsorption::GetNumberOfSecondaries()
00108 {
00109 return ( ngkine );
00110
00111 }
00112
00113
00114 G4GHEKinematicsVector* G4KaonMinusAbsorption::GetSecondaryKinematics()
00115 {
00116 return ( &gkin[0] );
00117
00118 }
00119
00120 G4double G4KaonMinusAbsorption::AtRestGetPhysicalInteractionLength(
00121 const G4Track& track,
00122 G4ForceCondition* condition
00123 )
00124 {
00125
00126 ResetNumberOfInteractionLengthLeft();
00127
00128
00129 *condition = NotForced;
00130
00131
00132 currentInteractionLength = GetMeanLifeTime(track, condition);
00133
00134 if ((currentInteractionLength <0.0) || (verboseLevel>2)){
00135 G4cout << "G4KaonMinusAbsorptionProcess::AtRestGetPhysicalInteractionLength ";
00136 G4cout << "[ " << GetProcessName() << "]" <<G4endl;
00137 track.GetDynamicParticle()->DumpInfo();
00138 G4cout << " in Material " << track.GetMaterial()->GetName() <<G4endl;
00139 G4cout << "MeanLifeTime = " << currentInteractionLength/ns << "[ns]" <<G4endl;
00140 }
00141
00142 return theNumberOfInteractionLengthLeft * currentInteractionLength;
00143
00144 }
00145
00146 G4VParticleChange* G4KaonMinusAbsorption::AtRestDoIt(
00147 const G4Track& track,
00148 const G4Step&
00149 )
00150
00151
00152
00153
00154
00155 {
00156
00157
00158
00159
00160
00161 aParticleChange.Initialize(track);
00162
00163
00164
00165 globalTime = track.GetGlobalTime()/s;
00166 G4Material * aMaterial = track.GetMaterial();
00167 const G4int numberOfElements = aMaterial->GetNumberOfElements();
00168 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
00169
00170 const G4double* theAtomicNumberDensity = aMaterial->GetAtomicNumDensityVector();
00171 G4double normalization = 0;
00172 for ( G4int i1=0; i1 < numberOfElements; i1++ )
00173 {
00174 normalization += theAtomicNumberDensity[i1] ;
00175
00176 }
00177 G4double runningSum= 0.;
00178 G4double random = G4UniformRand()*normalization;
00179 for ( G4int i2=0; i2 < numberOfElements; i2++ )
00180 {
00181 runningSum += theAtomicNumberDensity[i2];
00182
00183 if (random<=runningSum)
00184 {
00185 targetCharge = G4double( ((*theElementVector)[i2])->GetZ());
00186 targetAtomicMass = (*theElementVector)[i2]->GetN();
00187 }
00188 }
00189 if (random>runningSum)
00190 {
00191 targetCharge = G4double((*theElementVector)[numberOfElements-1]->GetZ());
00192 targetAtomicMass = (*theElementVector)[numberOfElements-1]->GetN();
00193
00194 }
00195
00196 if (verboseLevel>1) {
00197 G4cout << "G4KaonMinusAbsorption::AtRestDoIt is invoked " <<G4endl;
00198 }
00199
00200 G4ParticleMomentum momentum;
00201 G4float localtime;
00202
00203 G4ThreeVector position = track.GetPosition();
00204
00205 GenerateSecondaries();
00206
00207 aParticleChange.SetNumberOfSecondaries( ngkine );
00208
00209 for ( G4int isec = 0; isec < ngkine; isec++ ) {
00210 G4DynamicParticle* aNewParticle = new G4DynamicParticle;
00211 aNewParticle->SetDefinition( gkin[isec].GetParticleDef() );
00212 aNewParticle->SetMomentum( gkin[isec].GetMomentum() * GeV );
00213
00214 localtime = globalTime + gkin[isec].GetTOF();
00215
00216 G4Track* aNewTrack = new G4Track( aNewParticle, localtime*s, position );
00217 aNewTrack->SetTouchableHandle(track.GetTouchableHandle());
00218 aParticleChange.AddSecondary( aNewTrack );
00219
00220 }
00221
00222 aParticleChange.ProposeLocalEnergyDeposit( 0.0*GeV );
00223
00224 aParticleChange.ProposeTrackStatus(fStopAndKill);
00225
00226
00227
00228 ResetNumberOfInteractionLengthLeft();
00229
00230 return &aParticleChange;
00231
00232 }
00233
00234
00235 void G4KaonMinusAbsorption::GenerateSecondaries()
00236 {
00237 static G4int index;
00238 static G4int l;
00239 static G4int nopt;
00240 static G4int i;
00241
00242
00243 for (i = 1; i <= MAX_SECONDARIES; ++i) {
00244 pv[i].SetZero();
00245 }
00246
00247 ngkine = 0;
00248 ntot = 0;
00249 result.SetZero();
00250 result.SetMass( massKaonMinus );
00251 result.SetKineticEnergyAndUpdate( 0. );
00252 result.SetTOF( 0. );
00253 result.SetParticleDef( pdefKaonMinus );
00254
00255 KaonMinusAbsorption(&nopt);
00256
00257
00258 if (ntot != 0 || result.GetParticleDef() != pdefKaonMinus) {
00259
00260
00261
00262
00263
00264
00265
00266 gkin[0] = result;
00267 gkin[0].SetTOF( result.GetTOF() * 5e-11 );
00268 ngkine = 1;
00269
00270
00271
00272
00273
00274 for (l = 1; l <= ntot; ++l) {
00275 index = l - 1;
00276
00277
00278
00279 if (ngkine < MAX_SECONDARIES) {
00280 gkin[ngkine] = eve[index];
00281 gkin[ngkine].SetTOF( eve[index].GetTOF() * 5e-11 );
00282 ++ngkine;
00283 }
00284 }
00285 }
00286 else {
00287
00288
00289 ngkine = 0;
00290 ntot = 0;
00291 globalTime += result.GetTOF() * G4float(5e-11);
00292 }
00293
00294
00295 ngkine = G4int(std::min(ngkine,G4int(MAX_SECONDARIES)));
00296
00297 }
00298
00299
00300 void G4KaonMinusAbsorption::Poisso(G4float xav, G4int *iran)
00301 {
00302 static G4int i;
00303 static G4float r, p1, p2, p3;
00304 static G4int fivex;
00305 static G4float rr, ran, rrr, ran1;
00306
00307
00308
00309
00310
00311
00312 if (xav > G4float(9.9)) {
00313
00314 Normal(&ran1);
00315 ran1 = xav + ran1 * std::sqrt(xav);
00316 *iran = G4int(ran1);
00317 if (*iran < 0) {
00318 *iran = 0;
00319 }
00320 }
00321 else {
00322 fivex = G4int(xav * G4float(5.));
00323 *iran = 0;
00324 if (fivex > 0) {
00325 r = std::exp(-G4double(xav));
00326 ran1 = G4UniformRand();
00327 if (ran1 > r) {
00328 rr = r;
00329 for (i = 1; i <= fivex; ++i) {
00330 ++(*iran);
00331 if (i <= 5) {
00332 rrr = std::pow(xav, G4float(i)) / NFac(i);
00333 }
00334
00335 if (i > 5) {
00336 rrr = std::exp(i * std::log(xav) -
00337 (i + G4float(.5)) * std::log(i * G4float(1.)) +
00338 i - G4float(.9189385));
00339 }
00340 rr += r * rrr;
00341 if (ran1 <= rr) {
00342 break;
00343 }
00344 }
00345 }
00346 }
00347 else {
00348
00349 p1 = xav * std::exp(-G4double(xav));
00350 p2 = xav * p1 / G4float(2.);
00351 p3 = xav * p2 / G4float(3.);
00352 ran = G4UniformRand();
00353 if (ran >= p3) {
00354 if (ran >= p2) {
00355 if (ran >= p1) {
00356 *iran = 0;
00357 }
00358 else {
00359 *iran = 1;
00360 }
00361 }
00362 else {
00363 *iran = 2;
00364 }
00365 }
00366 else {
00367 *iran = 3;
00368 }
00369 }
00370 }
00371
00372 }
00373
00374
00375 G4int G4KaonMinusAbsorption::NFac(G4int n)
00376 {
00377 G4int ret_val;
00378
00379 static G4int i, j;
00380
00381
00382
00383
00384 ret_val = 1;
00385 j = n;
00386 if (j > 1) {
00387 if (j > 10) {
00388 j = 10;
00389 }
00390 for (i = 2; i <= j; ++i) {
00391 ret_val *= i;
00392 }
00393 }
00394 return ret_val;
00395
00396 }
00397
00398
00399 void G4KaonMinusAbsorption::Normal(G4float *ran)
00400 {
00401 static G4int i;
00402
00403
00404
00405
00406 *ran = G4float(-6.);
00407 for (i = 1; i <= 12; ++i) {
00408 *ran += G4UniformRand();
00409 }
00410
00411 }
00412
00413
00414 void G4KaonMinusAbsorption::KaonMinusAbsorption(G4int *nopt)
00415 {
00416 static G4int i;
00417 static G4int nt, nbl;
00418 static G4float ran, pcm;
00419 static G4int isw;
00420 static G4float tex;
00421 static G4float ran2, tof1, ekin, ekin1, ekin2, black;
00422 static G4float pnrat;
00423 static G4ParticleDefinition* ipa1;
00424 static G4ParticleDefinition* inve;
00425
00426
00427
00428
00429
00430
00431
00432
00433 pv[1].SetZero();
00434 pv[1].SetMass( massKaonMinus );
00435 pv[1].SetKineticEnergyAndUpdate( 0. );
00436 pv[1].SetTOF( result.GetTOF() );
00437 pv[1].SetParticleDef( result.GetParticleDef() );
00438 if (targetAtomicMass <= G4float(1.5)) {
00439 ran = G4UniformRand();
00440 tof1 = std::log(ran) * G4float(-12.5);
00441 tof1 *= G4float(20.);
00442 ran = G4UniformRand();
00443 isw = 1;
00444 if (ran < G4float(.33)) {
00445 isw = 2;
00446 }
00447 *nopt = isw;
00448 pv[3].SetZero();
00449 pv[3].SetMass( massLambda );
00450 pv[3].SetKineticEnergyAndUpdate( 0. );
00451 pv[3].SetTOF( result.GetTOF() + tof1 );
00452 pv[3].SetParticleDef( pdefLambda );
00453 pcm = massKaonMinus + massProton - massLambda;
00454 if (isw != 1) {
00455 pv[2].SetZero();
00456 pv[2].SetMass( massGamma );
00457 pv[2].SetKineticEnergyAndUpdate( pcm );
00458 pv[2].SetTOF( result.GetTOF() + tof1 );
00459 pv[2].SetParticleDef( pdefGamma );
00460 }
00461 else {
00462 pcm = pcm * pcm - massPionZero * massPionZero;
00463 if (pcm <= G4float(0.)) {
00464 pcm = G4float(0.);
00465 }
00466 pv[2].SetZero();
00467 pv[2].SetEnergy( std::sqrt(pcm + massPionZero * massPionZero) );
00468 pv[2].SetMassAndUpdate( massPionZero );
00469 pv[2].SetTOF( result.GetTOF() + tof1 );
00470 pv[2].SetParticleDef( pdefPionZero );
00471 }
00472 result = pv[2];
00473 if (ntot < MAX_SECONDARIES-1) {
00474 eve[ntot++] = pv[3];
00475 }
00476 }
00477 else {
00478
00479
00480
00481 evapEnergy1 = G4float(.3);
00482 evapEnergy3 = G4float(.15);
00483 nt = 1;
00484 tex = evapEnergy1;
00485 black = std::log(targetAtomicMass) * G4float(.5);
00486 Poisso(black, &nbl);
00487 if (nt + nbl > (MAX_SECONDARIES - 2)) {
00488 nbl = (MAX_SECONDARIES - 2) - nt;
00489 }
00490 if (nbl <= 0) {
00491 nbl = 1;
00492 }
00493 ekin = tex / nbl;
00494 ekin2 = G4float(0.);
00495 for (i = 1; i <= nbl; ++i) {
00496 if (nt == (MAX_SECONDARIES - 2)) {
00497 continue;
00498 }
00499 ran2 = G4UniformRand();
00500 ekin1 = -G4double(ekin) * std::log(ran2);
00501 ekin2 += ekin1;
00502 ipa1 = pdefNeutron;
00503 pnrat = G4float(1.) - targetCharge / targetAtomicMass;
00504 if (G4UniformRand() > pnrat) {
00505 ipa1 = pdefProton;
00506 }
00507 ++nt;
00508 pv[nt].SetZero();
00509 pv[nt].SetMass( ipa1->GetPDGMass()/GeV );
00510 pv[nt].SetKineticEnergyAndUpdate( ekin1 );
00511 pv[nt].SetTOF( 2. );
00512 pv[nt].SetParticleDef( ipa1 );
00513 if (ekin2 > tex) {
00514 break;
00515 }
00516 }
00517 tex = evapEnergy3;
00518 black = std::log(targetAtomicMass) * G4float(.5);
00519 Poisso(black, &nbl);
00520 if (nt + nbl > (MAX_SECONDARIES - 2)) {
00521 nbl = (MAX_SECONDARIES - 2) - nt;
00522 }
00523 if (nbl <= 0) {
00524 nbl = 1;
00525 }
00526 ekin = tex / nbl;
00527 ekin2 = G4float(0.);
00528 for (i = 1; i <= nbl; ++i) {
00529 if (nt == (MAX_SECONDARIES - 2)) {
00530 continue;
00531 }
00532 ran2 = G4UniformRand();
00533 ekin1 = -G4double(ekin) * std::log(ran2);
00534 ekin2 += ekin1;
00535 ++nt;
00536 ran = G4UniformRand();
00537 inve = pdefDeuteron;
00538 if (ran > G4float(.6)) {
00539 inve = pdefTriton;
00540 }
00541 if (ran > G4float(.9)) {
00542 inve = pdefAlpha;
00543 }
00544
00545 pv[nt].SetZero();
00546 pv[nt].SetMass( inve->GetPDGMass()/GeV );
00547 pv[nt].SetKineticEnergyAndUpdate( ekin1 );
00548 pv[nt].SetTOF( 2. );
00549 pv[nt].SetParticleDef( inve );
00550 if (ekin2 > tex) {
00551 break;
00552 }
00553 }
00554
00555
00556
00557 ran = G4UniformRand();
00558 tof1 = std::log(ran) * G4float(-12.5);
00559 tof1 *= G4float(20.);
00560 for (i = 2; i <= nt; ++i) {
00561 pv[i].SetTOF( result.GetTOF() + tof1 );
00562 }
00563 result = pv[2];
00564 for (i = 3; i <= nt; ++i) {
00565 if (ntot >= MAX_SECONDARIES) {
00566 break;
00567 }
00568 eve[ntot++] = pv[i];
00569 }
00570 }
00571
00572 }