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 #include "G4MuonMinusCaptureAtRest.hh"
00044 #include "G4DynamicParticle.hh"
00045 #include "Randomize.hh"
00046 #include "G4PhysicalConstants.hh"
00047 #include "G4SystemOfUnits.hh"
00048 #include "G4He3.hh"
00049 #include "G4NeutrinoMu.hh"
00050 #include "G4Fragment.hh"
00051 #include "G4ReactionProductVector.hh"
00052 #include "G4Proton.hh"
00053 #include "G4PionPlus.hh"
00054 #include "G4GHEKinematicsVector.hh"
00055 #include "G4Fancy3DNucleus.hh"
00056 #include "G4ExcitationHandler.hh"
00057 #include "G4HadronicProcessStore.hh"
00058
00059
00060
00061 G4MuonMinusCaptureAtRest::G4MuonMinusCaptureAtRest(const G4String& processName,
00062 G4ProcessType aType ) :
00063 G4VRestProcess (processName, aType), nCascade(0), targetZ(0), targetA(0),
00064 isInitialised(false)
00065 {
00066 SetProcessSubType(fHadronAtRest);
00067 Cascade = new G4GHEKinematicsVector [17];
00068 pSelector = new G4StopElementSelector();
00069 pEMCascade = new G4MuMinusCaptureCascade();
00070 theN = new G4Fancy3DNucleus();
00071 theHandler = new G4ExcitationHandler();
00072 G4HadronicProcessStore::Instance()->RegisterExtraProcess(this);
00073 }
00074
00075
00076
00077 G4MuonMinusCaptureAtRest::~G4MuonMinusCaptureAtRest()
00078 {
00079 G4HadronicProcessStore::Instance()->DeRegisterExtraProcess(this);
00080 delete [] Cascade;
00081 delete pSelector;
00082 delete pEMCascade;
00083 delete theN;
00084 delete theHandler;
00085 }
00086
00087
00088
00089 G4bool G4MuonMinusCaptureAtRest::IsApplicable(const G4ParticleDefinition& p)
00090 {
00091 return ( &p == G4MuonMinus::MuonMinus() );
00092 }
00093
00094
00095
00096 void G4MuonMinusCaptureAtRest::PreparePhysicsTable(const G4ParticleDefinition& p)
00097 {
00098 G4HadronicProcessStore::Instance()->RegisterParticleForExtraProcess(this, &p);
00099 }
00100
00101
00102
00103 void G4MuonMinusCaptureAtRest::BuildPhysicsTable(const G4ParticleDefinition& p)
00104 {
00105 G4HadronicProcessStore::Instance()->PrintInfo(&p);
00106 }
00107
00108
00109
00110 G4VParticleChange* G4MuonMinusCaptureAtRest::AtRestDoIt(const G4Track& track,
00111 const G4Step&)
00112 {
00113
00114
00115
00116
00117
00118 aParticleChange.Initialize(track);
00119
00120
00121 G4Element* aEle = pSelector->GetElement(track.GetMaterial());
00122 targetZ = aEle->GetZ();
00123 targetA = G4double(G4int(aEle->GetN()+0.5));
00124 G4int ni = 0;
00125
00126 G4IsotopeVector* isv = aEle->GetIsotopeVector();
00127 if(isv) ni = isv->size();
00128
00129 if(ni == 1) {
00130 targetA = G4double(aEle->GetIsotope(0)->GetN());
00131 } else if(ni > 1) {
00132 G4double* ab = aEle->GetRelativeAbundanceVector();
00133 G4double y = G4UniformRand();
00134 G4int j = -1;
00135 ni--;
00136 do {
00137 j++;
00138 y -= ab[j];
00139 } while (y > 0.0 && j < ni);
00140 targetA = G4double(aEle->GetIsotope(j)->GetN());
00141 }
00142
00143
00144 nCascade = 0;
00145 targetMass = G4NucleiProperties::GetNuclearMass(targetA, targetZ);
00146 nCascade = pEMCascade->DoCascade(targetZ, targetMass, Cascade);
00147
00148
00149 G4double lambdac = pSelector->GetMuonCaptureRate(targetZ, targetA);
00150 G4double lambdad = pSelector->GetMuonDecayRate(targetZ, targetA);
00151 G4double lambda = lambdac + lambdad;
00152
00153
00154
00155 G4double tDelay = -std::log(G4UniformRand()) / lambda;
00156
00157 G4ReactionProductVector * captureResult = 0;
00158 G4int nEmSecondaries = nCascade;
00159 G4int nSecondaries = nCascade;
00160
00161
00162
00163
00164 if( G4UniformRand()*lambda > lambdac)
00165 pEMCascade->DoBoundMuonMinusDecay(targetZ, &nEmSecondaries, Cascade);
00166 else
00167 captureResult = DoMuCapture();
00168
00169
00170 if(captureResult) nSecondaries += captureResult->size();
00171 else nSecondaries = nEmSecondaries;
00172
00173
00174 aParticleChange.SetNumberOfSecondaries( nSecondaries );
00175
00176 G4double globalTime = track.GetGlobalTime();
00177 G4ThreeVector position = track.GetPosition();
00178
00179 if(captureResult) {
00180 G4int n = captureResult->size();
00181 for ( G4int isec = 0; isec < n; isec++ ) {
00182 G4ReactionProduct* aParticle = captureResult->operator[](isec);
00183 G4DynamicParticle * aNewParticle = new G4DynamicParticle();
00184 aNewParticle->SetDefinition( aParticle->GetDefinition() );
00185 G4LorentzVector itV(aParticle->GetTotalEnergy(), aParticle->GetMomentum());
00186 aNewParticle->SetMomentum(itV.vect());
00187 G4double localtime = globalTime + tDelay + aParticle->GetTOF();
00188 G4Track* aNewTrack = new G4Track( aNewParticle, localtime, position);
00189 aNewTrack->SetTouchableHandle(track.GetTouchableHandle());
00190 aParticleChange.AddSecondary( aNewTrack );
00191 delete aParticle;
00192 }
00193 delete captureResult;
00194 }
00195
00196
00197
00198 if(nEmSecondaries > 0) {
00199
00200 for ( G4int isec = 0; isec < nEmSecondaries; isec++ ) {
00201 G4ParticleDefinition* pd = Cascade[isec].GetParticleDef();
00202 G4double localtime = globalTime;
00203 if(isec >= nCascade) localtime += tDelay;
00204 if(pd) {
00205 G4DynamicParticle* aNewParticle = new G4DynamicParticle;
00206 aNewParticle->SetDefinition( pd );
00207 aNewParticle->SetMomentum( Cascade[isec].GetMomentum() );
00208
00209 G4Track* aNewTrack = new G4Track( aNewParticle, localtime, position );
00210 aNewTrack->SetTouchableHandle(track.GetTouchableHandle());
00211 aParticleChange.AddSecondary( aNewTrack );
00212 }
00213 }
00214 }
00215
00216 aParticleChange.ProposeLocalEnergyDeposit(0.0);
00217 aParticleChange.ProposeTrackStatus(fStopAndKill);
00218
00219 return &aParticleChange;
00220 }
00221
00222
00223
00224 G4ReactionProductVector* G4MuonMinusCaptureAtRest::DoMuCapture()
00225 {
00226 G4double mumass = G4MuonMinus::MuonMinus()->GetPDGMass();
00227 G4double muBindingEnergy = pEMCascade->GetKShellEnergy(targetZ);
00228
00229
00230
00231
00232
00233 G4double muEnergy = mumass + muBindingEnergy;
00234 G4double muMom = std::sqrt(muBindingEnergy*(muBindingEnergy + 2.0*mumass));
00235 G4double availableEnergy = targetMass + mumass - muBindingEnergy;
00236 G4LorentzVector momInitial(0.0,0.0,0.0,availableEnergy);
00237 G4LorentzVector momResidual;
00238
00239 G4ThreeVector vmu = muMom*pEMCascade->GetRandomVec();
00240 G4LorentzVector aMuMom(vmu, muEnergy);
00241
00242 G4double residualMass =
00243 G4NucleiProperties::GetNuclearMass(targetA, targetZ - 1.0);
00244
00245 G4ReactionProductVector* aPreResult = 0;
00246 G4ReactionProduct* aNu = new G4ReactionProduct();
00247 aNu->SetDefinition( G4NeutrinoMu::NeutrinoMu() );
00248
00249 G4int iz = G4int(targetZ);
00250 G4int ia = G4int(targetA);
00251
00252
00253 if(iz <= 2) {
00254
00255 if(ia > 1) {
00256 if(iz == 1 && ia == 2) {
00257 availableEnergy -= neutron_mass_c2;
00258 } else if(iz == 1 && ia == 3) {
00259 availableEnergy -= 2.0*neutron_mass_c2;
00260 } else if(iz == 2) {
00261 G4ParticleDefinition* pd = 0;
00262 if (ia == 3) {
00263 pd = G4Deuteron::Deuteron();
00264 } else if(ia == 4) {
00265 pd = G4Triton::Triton();
00266 } else {
00267 pd = G4ParticleTable::GetParticleTable()->FindIon(1,ia-1,0,1);
00268 }
00269
00270
00271 availableEnergy -= pd->GetPDGMass();
00272 }
00273 }
00274
00275
00276
00277 G4double Enu = 0.5*(availableEnergy -
00278 neutron_mass_c2*neutron_mass_c2/availableEnergy);
00279
00280
00281 G4ThreeVector nu3Mom = Enu*pEMCascade->GetRandomVec();
00282
00283 G4ReactionProduct* aN = new G4ReactionProduct();
00284 aN->SetDefinition( G4Neutron::Neutron() );
00285 aN->SetTotalEnergy( availableEnergy - Enu );
00286 aN->SetMomentum( -nu3Mom );
00287
00288 aNu->SetTotalEnergy( Enu );
00289 aNu->SetMomentum( nu3Mom );
00290 aPreResult = new G4ReactionProductVector();
00291
00292 aPreResult->push_back(aN );
00293 aPreResult->push_back(aNu);
00294
00295 if(verboseLevel > 1)
00296 G4cout << "DoMuCapture on H or He"
00297 <<" EkinN(MeV)= " << (availableEnergy - Enu - neutron_mass_c2)/MeV
00298 <<" Enu(MeV)= "<<aNu->GetTotalEnergy()/MeV
00299 <<" n= " << aPreResult->size()
00300 <<G4endl;
00301
00302 return aPreResult;
00303 }
00304
00305
00306 G4double eEx;
00307 do {
00308 theN->Init(ia, iz);
00309 G4LorentzVector thePMom;
00310 G4Nucleon * aNucleon = 0;
00311 G4int theProtonCounter = G4int( targetZ * G4UniformRand() );
00312 G4int counter = 0;
00313 theN->StartLoop();
00314
00315 while( (aNucleon=theN->GetNextNucleon()) ) {
00316
00317 if( aNucleon->GetDefinition() == G4Proton::Proton() ) {
00318 counter++;
00319 if(counter == theProtonCounter) {
00320 thePMom = aNucleon->GetMomentum();
00321 break;
00322 }
00323 }
00324 }
00325
00326
00327 G4LorentzVector theCMS = thePMom + aMuMom;
00328 G4ThreeVector bst = theCMS.boostVector();
00329
00330 G4double Ecms = theCMS.mag();
00331 G4double Enu = 0.5*(Ecms - neutron_mass_c2*neutron_mass_c2/Ecms);
00332 eEx = 0.0;
00333
00334 if(Enu > 0.0) {
00335
00336 G4ThreeVector nu3Mom = Enu*pEMCascade->GetRandomVec();
00337 G4LorentzVector nuMom(nu3Mom, Enu);
00338
00339
00340 nuMom.boost(bst);
00341 aNu->SetTotalEnergy( nuMom.e() );
00342 aNu->SetMomentum( nuMom.vect() );
00343
00344
00345 momResidual = momInitial - nuMom;
00346
00347
00348 eEx = momResidual.mag();
00349 if(verboseLevel > 1)
00350 G4cout << "G4MuonMinusCaptureAtRest::DoMuCapture: "
00351 << " Eex(MeV)= " << (eEx-residualMass)/MeV
00352 << " Enu(MeV)= "<<aNu->GetTotalEnergy()/MeV
00353 <<G4endl;
00354 }
00355 } while(eEx <= residualMass);
00356
00357
00358
00359
00360
00361
00362
00363
00364 G4ThreeVector fromBreit = momResidual.boostVector();
00365 G4LorentzVector fscm(0.0,0.0,0.0, eEx);
00366 G4Fragment anInitialState;
00367 anInitialState.SetA(targetA);
00368 anInitialState.SetZ(G4double(iz - 1));
00369 anInitialState.SetNumberOfParticles(2);
00370 anInitialState.SetNumberOfCharged(0);
00371 anInitialState.SetNumberOfHoles(1);
00372 anInitialState.SetMomentum(fscm);
00373 aPreResult = theHandler->BreakItUp(anInitialState);
00374
00375 G4ReactionProductVector::iterator ires;
00376 G4double eBal = availableEnergy - aNu->GetTotalEnergy();
00377 for(ires=aPreResult->begin(); ires!=aPreResult->end(); ires++) {
00378 G4LorentzVector itV((*ires)->GetTotalEnergy(), (*ires)->GetMomentum());
00379 itV.boost(fromBreit);
00380 (*ires)->SetTotalEnergy(itV.t());
00381 (*ires)->SetMomentum(itV.vect());
00382 eBal -= itV.t();
00383 }
00384
00385
00386
00387 aPreResult->push_back(aNu);
00388
00389 if(verboseLevel > 1)
00390 G4cout << "DoMuCapture: Nsec= "
00391 << aPreResult->size() << " Ebalance(MeV)= " << eBal/MeV
00392 <<" E0(MeV)= " <<availableEnergy/MeV
00393 <<" Mres(GeV)= " <<residualMass/GeV
00394 <<G4endl;
00395
00396 return aPreResult;
00397 }
00398