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
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00060
00061 #include "G4EMDissociation.hh"
00062 #include "G4PhysicalConstants.hh"
00063 #include "G4SystemOfUnits.hh"
00064 #include "G4Evaporation.hh"
00065 #include "G4FermiBreakUp.hh"
00066 #include "G4StatMF.hh"
00067 #include "G4ParticleDefinition.hh"
00068 #include "G4LorentzVector.hh"
00069 #include "G4PhysicsFreeVector.hh"
00070 #include "G4EMDissociationCrossSection.hh"
00071 #include "G4Proton.hh"
00072 #include "G4Neutron.hh"
00073 #include "G4ParticleTable.hh"
00074 #include "G4IonTable.hh"
00075 #include "G4GeneralPhaseSpaceDecay.hh"
00076 #include "G4DecayProducts.hh"
00077 #include "G4DynamicParticle.hh"
00078 #include "G4Fragment.hh"
00079 #include "G4ReactionProductVector.hh"
00080 #include "Randomize.hh"
00081 #include "globals.hh"
00082
00083 G4EMDissociation::G4EMDissociation():G4HadronicInteraction("EMDissociation") {
00084
00085
00086
00087 PrintWelcomeMessage();
00088
00089
00090 theExcitationHandler = new G4ExcitationHandler;
00091 G4Evaporation* theEvaporation = new G4Evaporation;
00092 G4FermiBreakUp* theFermiBreakUp = new G4FermiBreakUp;
00093 G4StatMF* theMF = new G4StatMF;
00094 theExcitationHandler->SetEvaporation(theEvaporation);
00095 theExcitationHandler->SetFermiModel(theFermiBreakUp);
00096 theExcitationHandler->SetMultiFragmentation(theMF);
00097 theExcitationHandler->SetMaxAandZForFermiBreakUp(12, 6);
00098 theExcitationHandler->SetMinEForMultiFrag(5.0*MeV);
00099 handlerDefinedInternally = true;
00100
00101
00102
00103 dissociationCrossSection = new G4EMDissociationCrossSection;
00104 thePhotonSpectrum = new G4EMDissociationSpectrum;
00105
00106
00107
00108 SetMinEnergy(100.0*MeV);
00109 SetMaxEnergy(500.0*GeV);
00110
00111
00112 verboseLevel = 0;
00113 }
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137 G4EMDissociation::G4EMDissociation (G4ExcitationHandler *aExcitationHandler)
00138 {
00139
00140
00141
00142
00143
00144 PrintWelcomeMessage();
00145
00146 theExcitationHandler = aExcitationHandler;
00147 handlerDefinedInternally = false;
00148
00149
00150
00151
00152
00153 dissociationCrossSection = new G4EMDissociationCrossSection;
00154 thePhotonSpectrum = new G4EMDissociationSpectrum;
00155
00156
00157
00158
00159
00160 SetMinEnergy(100.0*MeV);
00161 SetMaxEnergy(500.0*GeV);
00162
00163
00164
00165
00166 verboseLevel = 0;
00167 }
00168
00169
00170 G4EMDissociation::~G4EMDissociation() {
00171 if (handlerDefinedInternally) delete theExcitationHandler;
00172
00173
00174
00175 delete thePhotonSpectrum;
00176 }
00177
00178
00179 G4HadFinalState *G4EMDissociation::ApplyYourself
00180 (const G4HadProjectile &theTrack, G4Nucleus &theTarget)
00181 {
00182
00183
00184
00185
00186
00187 theParticleChange.Clear();
00188 theParticleChange.SetStatusChange(stopAndKill);
00189
00190
00191
00192
00193
00194
00195 const G4ParticleDefinition *definitionP = theTrack.GetDefinition();
00196 const G4double AP = definitionP->GetBaryonNumber();
00197 const G4double ZP = definitionP->GetPDGCharge();
00198 G4LorentzVector pP = theTrack.Get4Momentum();
00199 G4double E = theTrack.GetKineticEnergy()/AP;
00200 G4double MP = theTrack.GetTotalEnergy() - E*AP;
00201 G4double b = pP.beta();
00202 G4double AT = theTarget.GetA_asInt();
00203 G4double ZT = theTarget.GetZ_asInt();
00204 G4double MT = G4NucleiProperties::GetNuclearMass(AT,ZT);
00205
00206
00207
00208
00209
00210 if (verboseLevel >= 2)
00211 {
00212 G4cout.precision(6);
00213 G4cout <<"########################################"
00214 <<"########################################"
00215 <<G4endl;
00216 G4cout <<"IN G4EMDissociation" <<G4endl;
00217 G4cout <<"Initial projectile A=" <<AP
00218 <<", Z=" <<ZP
00219 <<G4endl;
00220 G4cout <<"Initial target A=" <<AT
00221 <<", Z=" <<ZT
00222 <<G4endl;
00223 G4cout <<"Projectile momentum and Energy/nuc = " <<pP <<" ," <<E <<G4endl;
00224 }
00225
00226
00227
00228
00229
00230 G4ParticleDefinition *typeNucleon = NULL;
00231 G4ParticleDefinition *typeDaughter = NULL;
00232 G4double Eg = 0.0;
00233 G4double mass = 0.0;
00234 G4ThreeVector boost = G4ThreeVector(0.0, 0.0, 0.0);
00235
00236
00237
00238
00239
00240
00241
00242 G4double bmin = thePhotonSpectrum->GetClosestApproach(AP, ZP, AT, ZT, b);
00243 G4PhysicsFreeVector *crossSectionP = dissociationCrossSection->
00244 GetCrossSectionForProjectile(AP, ZP, AT, ZT, b, bmin);
00245 G4PhysicsFreeVector *crossSectionT = dissociationCrossSection->
00246 GetCrossSectionForTarget(AP, ZP, AT, ZT, b, bmin);
00247
00248 G4double totCrossSectionP = (*crossSectionP)[0]+(*crossSectionP)[1];
00249 G4double totCrossSectionT = (*crossSectionT)[0]+(*crossSectionT)[1];
00250
00251
00252
00253
00254
00255 if (G4UniformRand() <
00256 totCrossSectionP / (totCrossSectionP + totCrossSectionT))
00257 {
00258
00259
00260
00261
00262
00263
00264
00265
00266 mass = MP;
00267 if (G4UniformRand() < dissociationCrossSection->
00268 GetWilsonProbabilityForProtonDissociation (AP, ZP))
00269 {
00270 if (verboseLevel >= 2)
00271 G4cout <<"Projectile underwent EM dissociation producing a proton"
00272 <<G4endl;
00273 typeNucleon = G4Proton::ProtonDefinition();
00274 typeDaughter = G4ParticleTable::GetParticleTable()->
00275 GetIon((G4int) ZP-1, (G4int) AP-1, 0.0);
00276 }
00277 else
00278 {
00279 if (verboseLevel >= 2)
00280 G4cout <<"Projectile underwent EM dissociation producing a neutron"
00281 <<G4endl;
00282 typeNucleon = G4Neutron::NeutronDefinition();
00283 typeDaughter = G4ParticleTable::GetParticleTable()->
00284 GetIon((G4int) ZP, (G4int) AP-1, 0.0);
00285 }
00286 if (G4UniformRand() < (*crossSectionP)[0]/totCrossSectionP)
00287 {
00288 Eg = crossSectionP->GetLowEdgeEnergy(0);
00289 if (verboseLevel >= 2)
00290 G4cout <<"Transition type was E1" <<G4endl;
00291 }
00292 else
00293 {
00294 Eg = crossSectionP->GetLowEdgeEnergy(1);
00295 if (verboseLevel >= 2)
00296 G4cout <<"Transition type was E2" <<G4endl;
00297 }
00298
00299
00300
00301
00302
00303
00304 pP.setE(pP.e()+Eg);
00305 boost = pP.findBoostToCM();
00306 }
00307 else
00308 {
00309
00310
00311
00312
00313
00314
00315
00316 mass = MT;
00317 if (G4UniformRand() < dissociationCrossSection->
00318 GetWilsonProbabilityForProtonDissociation (AT, ZT))
00319 {
00320 if (verboseLevel >= 2)
00321 G4cout <<"Target underwent EM dissociation producing a proton"
00322 <<G4endl;
00323 typeNucleon = G4Proton::ProtonDefinition();
00324 typeDaughter = G4ParticleTable::GetParticleTable()->
00325 GetIon((G4int) ZT-1, (G4int) AT-1, 0.0);
00326 }
00327 else
00328 {
00329 if (verboseLevel >= 2)
00330 G4cout <<"Target underwent EM dissociation producing a neutron"
00331 <<G4endl;
00332 typeNucleon = G4Neutron::NeutronDefinition();
00333 typeDaughter = G4ParticleTable::GetParticleTable()->
00334 GetIon((G4int) ZT, (G4int) AT-1, 0.0);
00335 }
00336 if (G4UniformRand() < (*crossSectionT)[0]/totCrossSectionT)
00337 {
00338 Eg = crossSectionT->GetLowEdgeEnergy(0);
00339 if (verboseLevel >= 2)
00340 G4cout <<"Transition type was E1" <<G4endl;
00341 }
00342 else
00343 {
00344 Eg = crossSectionT->GetLowEdgeEnergy(1);
00345 if (verboseLevel >= 2)
00346 G4cout <<"Transition type was E2" <<G4endl;
00347 }
00348
00349
00350
00351
00352
00353
00354 G4ThreeVector v = pP.vect();
00355 v.setMag(1.0);
00356 G4DynamicParticle *changedP = new G4DynamicParticle
00357 (const_cast<G4ParticleDefinition*>(definitionP), v, E*AP-Eg);
00358 theParticleChange.AddSecondary (changedP);
00359 if (verboseLevel >= 2)
00360 {
00361 G4cout <<"Projectile change:" <<G4endl;
00362 changedP->DumpInfo();
00363 }
00364 }
00365
00366
00367
00368
00369
00370
00371
00372
00373 G4double e = mass + Eg;
00374 G4double mass1 = typeNucleon->GetPDGMass();
00375 G4double mass2 = typeDaughter->GetPDGMass();
00376 G4double pp = (e+mass1+mass2)*(e+mass1-mass2)*
00377 (e-mass1+mass2)*(e-mass1-mass2)/(4.0*e*e);
00378 if (pp < 0.0)
00379 {
00380 pp = 1.0*eV;
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390 }
00391 else
00392 pp = std::sqrt(pp);
00393 G4double costheta = 2.*G4UniformRand()-1.0;
00394 G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
00395 G4double phi = 2.0*pi*G4UniformRand()*rad;
00396 G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
00397 G4DynamicParticle *dynamicNucleon =
00398 new G4DynamicParticle(typeNucleon, direction*pp);
00399 dynamicNucleon->Set4Momentum(dynamicNucleon->Get4Momentum().boost(-boost));
00400 G4DynamicParticle *dynamicDaughter =
00401 new G4DynamicParticle(typeDaughter, -direction*pp);
00402 dynamicDaughter->Set4Momentum(dynamicDaughter->Get4Momentum().boost(-boost));
00403
00404
00405
00406
00407
00408 theParticleChange.AddSecondary (dynamicNucleon);
00409 if (verboseLevel >= 2)
00410 {
00411 G4cout <<"Nucleon from the EMD process:" <<G4endl;
00412 dynamicNucleon->DumpInfo();
00413 }
00414
00415 G4Fragment *theFragment = new
00416 G4Fragment((G4int) typeDaughter->GetBaryonNumber(),
00417 (G4int) typeDaughter->GetPDGCharge(), dynamicDaughter->Get4Momentum());
00418 if (verboseLevel >= 2)
00419 {
00420 G4cout <<"Dynamic properties of the prefragment:" <<G4endl;
00421 G4cout.precision(6);
00422 dynamicDaughter->DumpInfo();
00423 G4cout <<"Nuclear properties of the prefragment:" <<G4endl;
00424 G4cout <<theFragment <<G4endl;
00425 }
00426 G4ReactionProductVector *products =
00427 theExcitationHandler->BreakItUp(*theFragment);
00428 delete theFragment;
00429 theFragment = NULL;
00430
00431 G4ReactionProductVector::iterator iter;
00432 for (iter = products->begin(); iter != products->end(); ++iter)
00433 {
00434 G4DynamicParticle *secondary =
00435 new G4DynamicParticle((*iter)->GetDefinition(),
00436 (*iter)->GetTotalEnergy(), (*iter)->GetMomentum());
00437 theParticleChange.AddSecondary (secondary);
00438 }
00439
00440 if (verboseLevel >= 2)
00441 G4cout <<"########################################"
00442 <<"########################################"
00443 <<G4endl;
00444
00445 return &theParticleChange;
00446 }
00448
00449 void G4EMDissociation::PrintWelcomeMessage ()
00450 {
00451 G4cout <<G4endl;
00452 G4cout <<" ****************************************************************"
00453 <<G4endl;
00454 G4cout <<" EM dissociation model for nuclear-nuclear interactions activated"
00455 <<G4endl;
00456 G4cout <<" (Written by QinetiQ Ltd for the European Space Agency)"
00457 <<G4endl;
00458 G4cout <<" ****************************************************************"
00459 <<G4endl;
00460 G4cout << G4endl;
00461
00462 return;
00463 }
00465