G4MuMinusCapturePrecompound Class Reference

#include <G4MuMinusCapturePrecompound.hh>

Inheritance diagram for G4MuMinusCapturePrecompound:

G4HadronicInteraction

Public Member Functions

 G4MuMinusCapturePrecompound (G4VPreCompoundModel *ptr=0)
 ~G4MuMinusCapturePrecompound ()
G4HadFinalStateApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
void ModelDescription (std::ostream &outFile) const

Detailed Description

Definition at line 65 of file G4MuMinusCapturePrecompound.hh.


Constructor & Destructor Documentation

G4MuMinusCapturePrecompound::G4MuMinusCapturePrecompound ( G4VPreCompoundModel ptr = 0  ) 

Definition at line 64 of file G4MuMinusCapturePrecompound.cc.

References G4HadronicInteractionRegistry::FindModel(), G4ParticleDefinition::GetPDGMass(), G4HadronicInteractionRegistry::Instance(), G4MuonMinus::MuonMinus(), G4Neutron::Neutron(), and G4Proton::Proton().

00066   : G4HadronicInteraction("muMinusNuclearCapture")
00067 { 
00068   fMuMass = G4MuonMinus::MuonMinus()->GetPDGMass(); 
00069   fProton = G4Proton::Proton();
00070   fNeutron = G4Neutron::Neutron();
00071   fThreshold = 10*MeV;
00072   fPreCompound = ptr;
00073   if(!ptr) { 
00074     G4HadronicInteraction* p =
00075       G4HadronicInteractionRegistry::Instance()->FindModel("PRECO");
00076     ptr = static_cast<G4VPreCompoundModel*>(p); 
00077     fPreCompound = ptr;
00078     if(!ptr) { fPreCompound = new G4PreCompoundModel(); }
00079   }
00080 }

G4MuMinusCapturePrecompound::~G4MuMinusCapturePrecompound (  ) 

Definition at line 84 of file G4MuMinusCapturePrecompound.cc.

References G4HadFinalState::Clear().

00085 {
00086   result.Clear();
00087 }


Member Function Documentation

G4HadFinalState * G4MuMinusCapturePrecompound::ApplyYourself ( const G4HadProjectile aTrack,
G4Nucleus targetNucleus 
) [virtual]

Implements G4HadronicInteraction.

Definition at line 92 of file G4MuMinusCapturePrecompound.cc.

References G4HadFinalState::Clear(), G4VPreCompoundModel::DeExcite(), FatalException, G4cout, G4endl, G4Exception(), G4RandomDirection(), G4UniformRand, G4DynamicParticle::Get4Momentum(), G4Nucleus::GetA_asInt(), G4HadProjectile::GetBoundEnergy(), G4HadProjectile::GetGlobalTime(), G4HadronicInteraction::GetModelName(), G4NucleiProperties::GetNuclearMass(), G4Fancy3DNucleus::GetNucleons(), G4HadFinalState::GetNumberOfSecondaries(), G4HadSecondary::GetParticle(), G4HadFinalState::GetSecondary(), G4Nucleus::GetZ_asInt(), G4Fancy3DNucleus::Init(), CLHEP::detail::n, G4NeutrinoMu::NeutrinoMu(), G4HadFinalState::SetStatusChange(), stopAndKill, G4Triton::Triton(), and G4HadronicInteraction::verboseLevel.

00094 {
00095   result.Clear();
00096   result.SetStatusChange(stopAndKill);
00097   fTime = projectile.GetGlobalTime();
00098   G4double time0 = fTime;
00099 
00100   G4double muBindingEnergy = projectile.GetBoundEnergy();
00101 
00102   G4int Z = targetNucleus.GetZ_asInt(); 
00103   G4int A = targetNucleus.GetA_asInt();
00104   G4double massA = G4NucleiProperties::GetNuclearMass(A, Z);
00105 
00106   /*
00107   G4cout << "G4MuMinusCapturePrecompound::ApplyYourself: Emu= "
00108          << muBindingEnergy << G4endl;
00109   */
00110   // Energy on K-shell
00111   G4double muEnergy = fMuMass + muBindingEnergy;
00112   G4double muMom = std::sqrt(muBindingEnergy*(muBindingEnergy + 2.0*fMuMass));
00113   G4double availableEnergy = massA + fMuMass - muBindingEnergy;
00114   G4double residualMass = G4NucleiProperties::GetNuclearMass(A, Z - 1);
00115 
00116   G4ThreeVector vmu = muMom*G4RandomDirection();
00117   G4LorentzVector aMuMom(vmu, muEnergy);
00118 
00119   // p or 3He as a target 
00120   // two body reaction mu- + A(Z,A) -> nuMu + A(Z-1,A)
00121   if((1 == Z && 1 == A) || (2 == Z && 3 == A)) {
00122 
00123     G4ParticleDefinition* pd = 0;
00124     if(1 == Z) { pd = fNeutron; }
00125     else { pd = G4Triton::Triton(); }
00126 
00127     //
00128     //  Computation in assumption of CM reaction
00129     //  
00130     G4double e = 0.5*(availableEnergy - 
00131                       residualMass*residualMass/availableEnergy);
00132 
00133     G4ThreeVector nudir = G4RandomDirection();
00134     AddNewParticle(G4NeutrinoMu::NeutrinoMu(), nudir, e);
00135     nudir *= -1.0;
00136     AddNewParticle(pd, nudir, availableEnergy - e - residualMass);
00137 
00138 
00139   } else {
00140     // sample mu- + p -> nuMu + n reaction in CM of muonic atom
00141 
00142     // muon
00143 //    
00144 // NOTE by K.Genser and J.Yarba:
00145 // The code below isn't working because emu always turns smaller than fMuMass
00146 // For this reason the sqrt is producing a NaN    
00147 //
00148 //    G4double emu = (availableEnergy*availableEnergy - massA*massA
00149 //                  + fMuMass*fMuMass)/(2*availableEnergy);
00150 //    G4ThreeVector mudir = G4RandomDirection();
00151 //    G4LorentzVector momMuon(std::sqrt(emu*emu - fMuMass*fMuMass)*mudir, emu);
00152 
00153     // nucleus
00154     G4LorentzVector momInitial(0.0,0.0,0.0,availableEnergy);
00155     G4LorentzVector momResidual, momNu;
00156 
00157     // pick random proton inside nucleus 
00158     G4double eEx;
00159     fNucleus.Init(A, Z);
00160     const std::vector<G4Nucleon>& nucleons= fNucleus.GetNucleons();
00161     G4ParticleDefinition* pDef;
00162 
00163     G4int nneutrons = 1;
00164     G4int reentryCount = 0;
00165   
00166     do {
00167       ++reentryCount;
00168       G4int index = 0;
00169       do {
00170         index=G4int(A*G4UniformRand());
00171         pDef = nucleons[index].GetDefinition();
00172       } while(pDef != fProton);
00173       G4LorentzVector momP = nucleons[index].Get4Momentum();
00174 
00175       // Get CMS kinematics
00176       G4LorentzVector theCMS = momP + aMuMom;
00177       G4ThreeVector bst = theCMS.boostVector();
00178 
00179       G4double Ecms = theCMS.mag();
00180       G4double Enu  = 0.5*(Ecms - neutron_mass_c2*neutron_mass_c2/Ecms);
00181       eEx = 0.0;
00182 
00183       if(Enu > 0.0) {
00184         // make the nu, and transform to lab;
00185         momNu.set(Enu*G4RandomDirection(), Enu);
00186 
00187         // nu in lab.
00188         momNu.boost(bst);
00189         momResidual = momInitial - momNu;
00190         eEx = momResidual.mag() - residualMass;
00191 
00192         // release neutron
00193         
00194         if(eEx > 0.0) {
00195           G4double eth = residualMass - massA + fThreshold + 2*neutron_mass_c2;
00196           if(Ecms - Enu > eth) {
00197             theCMS -= momNu;
00198             G4double ekin =  theCMS.e() - eth;
00199             G4ThreeVector dir = theCMS.vect().unit();
00200             AddNewParticle(fNeutron, dir, ekin);
00201             momResidual -= 
00202               result.GetSecondary(0)->GetParticle()->Get4Momentum();
00203             --Z;
00204             --A; 
00205             residualMass = G4NucleiProperties::GetNuclearMass(A, Z);
00206             nneutrons = 0;  
00207           }
00208         }
00209       }
00210       if(Enu <= 0.0 && eEx <= 0.0 && reentryCount > 100) {
00211         G4ExceptionDescription ed;
00212         ed << "Call for " << GetModelName() << G4endl;
00213         ed << "Target  Z= " << Z  
00214            << "  A= " << A << G4endl;
00215         ed << " ApplyYourself does not completed after 100 attempts" << G4endl;
00216         G4Exception("G4MuMinusCapturePrecompound::AtRestDoIt", "had006", 
00217                     FatalException, ed);        
00218       }
00219     } while(eEx <= 0.0);
00220 
00221     G4ThreeVector dir = momNu.vect().unit();
00222     AddNewParticle(G4NeutrinoMu::NeutrinoMu(), dir, momNu.e());
00223 
00224     G4Fragment initialState(A, Z, momResidual);
00225     initialState.SetNumberOfExcitedParticle(nneutrons,0);
00226     initialState.SetNumberOfHoles(1,1);
00227 
00228     // decay time for pre-compound/de-excitation starts from zero
00229     G4ReactionProductVector* rpv = fPreCompound->DeExcite(initialState);
00230     size_t n = rpv->size();
00231     for(size_t i=0; i<n; ++i) {
00232       G4ReactionProduct* rp = (*rpv)[i];
00233 
00234       // reaction time
00235       fTime = time0 + rp->GetTOF();
00236       G4ThreeVector direction = rp->GetMomentum().unit();
00237       AddNewParticle(rp->GetDefinition(), direction, rp->GetKineticEnergy());
00238       delete rp;
00239     }
00240     delete rpv;
00241   } 
00242   if(verboseLevel > 1)
00243     G4cout << "G4MuMinusCapturePrecompound::ApplyYourself:  Nsec= " 
00244            << result.GetNumberOfSecondaries() 
00245            <<" E0(MeV)= " <<availableEnergy/MeV
00246            <<" Mres(GeV)= " <<residualMass/GeV
00247            <<G4endl;
00248 
00249   return &result;
00250 }

void G4MuMinusCapturePrecompound::ModelDescription ( std::ostream &  outFile  )  const [virtual]

Reimplemented from G4HadronicInteraction.

Definition at line 254 of file G4MuMinusCapturePrecompound.cc.

00255 {
00256   outFile << "Sampling of mu- capture by atomic nucleus from K-shell"
00257           << " mesoatom orbit.\n"
00258           << "Primary reaction mu- + p -> n + neutrino, neutron providing\n"
00259           << "  initial excitation of the target nucleus and PreCompound"
00260           << " model samples final state\n";
00261 }


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:52:32 2013 for Geant4 by  doxygen 1.4.7