G4MuonMinusCaptureAtRest.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 //   G4MuonMinusCaptureAtRest physics process
00029 //   Larry Felawka (TRIUMF) and Art Olin (TRIUMF)
00030 //   April 1998
00031 //---------------------------------------------------------------------
00032 //
00033 // Modifications: 
00034 // 18/08/2000  V.Ivanchenko Update description 
00035 // 12/12/2003  H.P.Wellisch Completly rewrite mu-nuclear part
00036 // 17/05/2006  V.Ivanchenko Cleanup
00037 // 15/11/2006  V.Ivanchenko Review and rewrite all kinematics
00038 // 24/01/2007  V.Ivanchenko Force to work with integer Z and A
00039 // 23/01/2009  V.Ivanchenko Add deregistration
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00088 
00089 G4bool G4MuonMinusCaptureAtRest::IsApplicable(const G4ParticleDefinition& p)
00090 {
00091   return ( &p == G4MuonMinus::MuonMinus() );
00092 }
00093 
00094 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00095 
00096 void G4MuonMinusCaptureAtRest::PreparePhysicsTable(const G4ParticleDefinition& p) 
00097 {
00098   G4HadronicProcessStore::Instance()->RegisterParticleForExtraProcess(this, &p);
00099 }
00100 
00101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00102 
00103 void G4MuonMinusCaptureAtRest::BuildPhysicsTable(const G4ParticleDefinition& p) 
00104 {
00105   G4HadronicProcessStore::Instance()->PrintInfo(&p);
00106 }
00107 
00108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00109 
00110 G4VParticleChange* G4MuonMinusCaptureAtRest::AtRestDoIt(const G4Track& track, 
00111                                                         const G4Step&)
00112 {
00113   //
00114   // Handles MuonMinuss at rest; a MuonMinus can either create secondaries or
00115   // do nothing (in which case it should be sent back to decay-handling
00116   // section
00117   //
00118   aParticleChange.Initialize(track);
00119 
00120   // select element and get Z,A.
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   // Do the electromagnetic cascade of the muon in the nuclear field.
00144   nCascade   = 0;
00145   targetMass = G4NucleiProperties::GetNuclearMass(targetA, targetZ);
00146   nCascade   = pEMCascade->DoCascade(targetZ, targetMass, Cascade);
00147 
00148   // Decide on Decay or Capture, and doit.
00149   G4double lambdac  = pSelector->GetMuonCaptureRate(targetZ, targetA);
00150   G4double lambdad  = pSelector->GetMuonDecayRate(targetZ, targetA);
00151   G4double lambda   = lambdac + lambdad;
00152 
00153   // ===  Throw for capture  time.
00154 
00155   G4double tDelay = -std::log(G4UniformRand()) / lambda;
00156   
00157   G4ReactionProductVector * captureResult = 0;
00158   G4int nEmSecondaries = nCascade;
00159   G4int nSecondaries = nCascade;
00160   /*
00161   G4cout << "lambda= " << lambda << " lambdac= " << lambdac 
00162          << " nem= " << nEmSecondaries << G4endl;
00163   */
00164   if( G4UniformRand()*lambda > lambdac) 
00165     pEMCascade->DoBoundMuonMinusDecay(targetZ, &nEmSecondaries, Cascade);
00166   else 
00167     captureResult = DoMuCapture();
00168   
00169   // fill the final state
00170   if(captureResult) nSecondaries += captureResult->size();
00171   else nSecondaries = nEmSecondaries;
00172   //G4cout << " nsec= " << nSecondaries << " nem= " << nEmSecondaries << G4endl;
00173 
00174   aParticleChange.SetNumberOfSecondaries( nSecondaries );
00175 
00176   G4double globalTime = track.GetGlobalTime();
00177   G4ThreeVector position = track.GetPosition();
00178   // Store nuclear cascade
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   // Store electromagnetic cascade
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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00223 
00224 G4ReactionProductVector* G4MuonMinusCaptureAtRest::DoMuCapture()
00225 {
00226   G4double mumass = G4MuonMinus::MuonMinus()->GetPDGMass();
00227   G4double muBindingEnergy = pEMCascade->GetKShellEnergy(targetZ);
00228   /*
00229   G4cout << "G4MuonMinusCaptureAtRest::DoMuCapture called Emu= "
00230          << muBindingEnergy << G4endl;
00231   */
00232   // Energy on K-shell
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   // p, d, t, 3He or alpha as target
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         //      G4cout << "Extra " << pd->GetParticleName() << G4endl;
00271         availableEnergy -= pd->GetPDGMass();
00272       }
00273     }
00274     //
00275     //  Computation in assumption of CM collision of mu and nucleaon
00276     //  
00277     G4double Enu  = 0.5*(availableEnergy - 
00278                          neutron_mass_c2*neutron_mass_c2/availableEnergy);
00279 
00280     // make the nu, and transform to lab;
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   // pick random proton inside nucleus 
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     // Get the nu momentum in the CMS
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       // make the nu, and transform to lab;
00336       G4ThreeVector nu3Mom = Enu*pEMCascade->GetRandomVec();
00337       G4LorentzVector nuMom(nu3Mom, Enu);
00338 
00339       // nu in lab.
00340       nuMom.boost(bst);
00341       aNu->SetTotalEnergy( nuMom.e() );
00342       aNu->SetMomentum( nuMom.vect() );
00343     
00344       // make residual
00345       momResidual = momInitial - nuMom;
00346 
00347       // Call pre-compound on the rest.
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 //  G4cout << "muonCapture : " << eEx << " " << residualMass 
00358 //         << " A,Z= " << targetA << ", "<< targetZ 
00359 //       << "  " << G4int(targetA) << ", " << G4int(targetZ) << G4endl;
00360 
00361   //
00362   // Start Deexcitation
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   // fill neutrino into result
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 

Generated on Mon May 27 17:48:54 2013 for Geant4 by  doxygen 1.4.7