Changeset 1315 for trunk/source/processes/electromagnetic/muons

Timestamp:

Jun 18, 2010, 11:42:07 AM (14 years ago)

Author:

garnier

Message:

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08

Location:

trunk/source/processes/electromagnetic/muons

Files:

• History (modified) (2 diffs)
• src/G4MuIonisation.cc (modified) (3 diffs)
• src/G4MuPairProductionModel.cc (modified) (8 diffs)
• test/muEnergyLossTest.cc (modified) (1 diff)

trunk/source/processes/electromagnetic/muons/History

@@ -1 @@
-$Id: History,v 1.134 2009/11/09 19:18:01 vnivanch Exp $
+$Id: History,v 1.136 2010/06/04 09:30:40 vnivanch Exp $
 -------------------------------------------------------------------
 
@@ -17 +17 @@
      * Reverse chronological order (last date on top), please *
      ----------------------------------------------------------
+
+4 June 2010: V.Ivant (emmuons-V09-03-01)
+- G4MuIonisation - use G4ICRU73QOModel for mu- for E< 0.2 MeV
+
+1 June 2010: V.Ivant (emmuons-V09-03-00)
+- G4MuPairProductionModel - added sampling recoil of a primary particle
+- G4MuIonisation - added G4ICRU73QOModel for mu- for E< 1 MeV
 
 09 November 09: V.Ivant (emmuons-V09-02-08)

trunk/source/processes/electromagnetic/muons/src/G4MuIonisation.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4MuIonisation.cc,v 1.59 2009/02/26 11:04:20 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4MuIonisation.cc,v 1.61 2010/06/04 09:30:40 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -89 +89 @@
 #include "G4BohrFluctuations.hh"
 #include "G4UnitsTable.hh"
+#include "G4ICRU73QOModel.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
@@ -140 +141 @@
 
     mass = theParticle->GetPDGMass();
+    G4double q = theParticle->GetPDGCharge();
+    G4double elow = 0.2*MeV;
     SetSecondaryParticle(G4Electron::Electron());
 
     // Bragg peak model
-    if (!EmModel(1)) SetEmModel(new G4BraggModel(),1);
+    if (!EmModel(1)) {
+      if(q > 0.0) { SetEmModel(new G4BraggModel(),1); }
+      else { 
+        SetEmModel(new G4ICRU73QOModel(),1); 
+        //elow = 1.0*MeV;
+      }
+    }
     EmModel(1)->SetLowEnergyLimit(MinKinEnergy());
-    EmModel(1)->SetHighEnergyLimit(0.2*MeV);
+    EmModel(1)->SetHighEnergyLimit(elow); 
     AddEmModel(1, EmModel(1), new G4IonFluctuations());
 
     // high energy fluctuation model
-    if (!FluctModel()) SetFluctModel(new G4UniversalFluctuation());
+    if (!FluctModel()) { SetFluctModel(new G4UniversalFluctuation()); }
 
     // moderate energy model
-    if (!EmModel(2)) SetEmModel(new G4BetheBlochModel(),2);
-    EmModel(2)->SetLowEnergyLimit(0.2*MeV);
+    if (!EmModel(2)) { SetEmModel(new G4BetheBlochModel(),2); }
+    EmModel(2)->SetLowEnergyLimit(elow);
     EmModel(2)->SetHighEnergyLimit(1.0*GeV);
     AddEmModel(2, EmModel(2), FluctModel());
 
     // high energy model
-    if (!EmModel(3)) SetEmModel(new G4MuBetheBlochModel(),3);
+    if (!EmModel(3)) { SetEmModel(new G4MuBetheBlochModel(),3); }
     EmModel(3)->SetLowEnergyLimit(1.0*GeV);
     EmModel(3)->SetHighEnergyLimit(MaxKinEnergy());

trunk/source/processes/electromagnetic/muons/src/G4MuPairProductionModel.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4MuPairProductionModel.cc,v 1.44 2009/08/11 16:50:07 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4MuPairProductionModel.cc,v 1.45 2010/06/01 15:21:59 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -466 +466 @@
 {
   G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
-  G4double totalEnergy   = kineticEnergy + particleMass ;
-  G4ParticleMomentum ParticleDirection = 
-    aDynamicParticle->GetMomentumDirection();
+  G4double totalEnergy   = kineticEnergy + particleMass;
+  G4double totalMomentum = 
+    sqrt(kineticEnergy*(kineticEnergy + 2.0*particleMass));
+
+  G4ThreeVector partDirection = aDynamicParticle->GetMomentumDirection();
 
   G4int it;
-  for(it=1; it<ntdat; it++) {if(kineticEnergy <= tdat[it]) break;}
-  if(it == ntdat) it--;
+  for(it=1; it<ntdat; ++it) {if(kineticEnergy <= tdat[it]) break;}
+  if(it == ntdat) { --it; }
   G4double dt = log(kineticEnergy/tdat[it-1])/log(tdat[it]/tdat[it-1]);
 
   // select randomly one element constituing the material
-  const G4Element* anElement = SelectRandomAtom(kineticEnergy, dt, it, couple, tmin);
+  const G4Element* anElement = 
+    SelectRandomAtom(kineticEnergy, dt, it, couple, tmin);
   SetCurrentElement(anElement->GetZ());
 
@@ -484 +487 @@
   G4double minEnergy     = std::max(tmin, minPairEnergy);
 
-  if(minEnergy >= maxEnergy) return;
+  if(minEnergy >= maxEnergy) { return; }
   //G4cout << "emin= " << minEnergy << " emax= " << maxEnergy 
   //     << " minPair= " << minPairEnergy << " maxpair= " << maxPairEnergy 
@@ -507 +510 @@
   G4int iz, iy;
 
-  for(iz=1; iz<nzdat; iz++) {if(currentZ <= zdat[iz]) break;}
-  if(iz == nzdat) iz--;
+  for(iz=1; iz<nzdat; ++iz) { if(currentZ <= zdat[iz]) { break; } }
+  if(iz == nzdat) { --iz; }
 
   G4double dz = log(currentZ/zdat[iz-1])/log(zdat[iz]/zdat[iz-1]);
@@ -520 +523 @@
   G4double p1 = pmin;
   G4double p2 = pmin;
-  for(iy=iymin+1; iy<=iymax; iy++) {
+  for(iy=iymin+1; iy<=iymax; ++iy) {
     p1 = p2;
     p2 = InterpolatedIntegralCrossSection(dt, dz, iz, it, iy, currentZ);
@@ -532 +535 @@
                        *log(maxPairEnergy/minPairEnergy));
                        
-  if(PairEnergy < minEnergy) PairEnergy = minEnergy;
-  if(PairEnergy > maxEnergy) PairEnergy = maxEnergy;
+  if(PairEnergy < minEnergy) { PairEnergy = minEnergy; }
+  if(PairEnergy > maxEnergy) { PairEnergy = maxEnergy; }
 
   // sample r=(E+-E-)/PairEnergy  ( uniformly .....)
@@ -545 +548 @@
   G4double PositronEnergy = PairEnergy - ElectronEnergy;
 
-  //  angles of the emitted particles ( Z - axis along the parent particle)
-  //      (mean theta for the moment)
-
-  //
-  // scattered electron (positron) angles. ( Z - axis along the parent photon)
-  //
-  //  universal distribution suggested by L. Urban
-  // (Geant3 manual (1993) Phys211),
-  //  derived from Tsai distribution (Rev Mod Phys 49,421(1977))
-  //  G4cout << "Ee= " << ElectronEnergy << " Ep= " << PositronEnergy << G4endl;
-  G4double u;
-  const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
-
-  if (9./(9.+d) >G4UniformRand()) u= - log(G4UniformRand()*G4UniformRand())/a1;
-  else                            u= - log(G4UniformRand()*G4UniformRand())/a2;
-
-  G4double TetEl = u*electron_mass_c2/ElectronEnergy;
-  G4double TetPo = u*electron_mass_c2/PositronEnergy;
-  G4double Phi  = twopi * G4UniformRand();
-  G4double dxEl= sin(TetEl)*cos(Phi),dyEl= sin(TetEl)*sin(Phi),dzEl=cos(TetEl);
-  G4double dxPo=-sin(TetPo)*cos(Phi),dyPo=-sin(TetPo)*sin(Phi),dzPo=cos(TetPo);
-
-  G4ThreeVector ElectDirection (dxEl, dyEl, dzEl);
-  ElectDirection.rotateUz(ParticleDirection);
+  // The angle of the emitted virtual photon is sampled
+  // according to the muon bremsstrahlung model
+ 
+  G4double gam  = totalEnergy/particleMass;
+  G4double gmax = gam*std::min(1.0, totalEnergy/PairEnergy - 1.0);
+  G4double gmax2= gmax*gmax;
+  G4double x = G4UniformRand()*gmax2/(1.0 + gmax2);
+
+  G4double theta = sqrt(x/(1.0 - x))/gam;
+  G4double sint  = sin(theta);
+  G4double phi   = twopi * G4UniformRand() ;
+  G4double dirx  = sint*cos(phi), diry = sint*sin(phi), dirz = cos(theta) ;
+
+  G4ThreeVector gDirection(dirx, diry, dirz);
+  gDirection.rotateUz(partDirection);
+
+  // the angles of e- and e+ assumed to be the same as virtual gamma
 
   // create G4DynamicParticle object for the particle1
-  G4DynamicParticle* aParticle1= new G4DynamicParticle(theElectron,
-                                                       ElectDirection,
-                                             ElectronEnergy - electron_mass_c2);
-
-  G4ThreeVector PositDirection (dxPo, dyPo, dzPo);
-  PositDirection.rotateUz(ParticleDirection);
+  G4DynamicParticle* aParticle1 = 
+    new G4DynamicParticle(theElectron, gDirection, 
+                          ElectronEnergy - electron_mass_c2);
 
   // create G4DynamicParticle object for the particle2
   G4DynamicParticle* aParticle2 = 
-    new G4DynamicParticle(thePositron,
-                          PositDirection,
+    new G4DynamicParticle(thePositron, gDirection,
                           PositronEnergy - electron_mass_c2);
 
@@ -588 +580 @@
   fParticleChange->SetProposedKineticEnergy(kineticEnergy);
 
+  partDirection *= totalMomentum;
+  partDirection -= (aParticle1->GetMomentum() + aParticle2->GetMomentum());
+  partDirection = partDirection.unit();
+  fParticleChange->SetProposedMomentumDirection(partDirection);
+
+  // add secondary
   vdp->push_back(aParticle1);
   vdp->push_back(aParticle2);

trunk/source/processes/electromagnetic/muons/test/muEnergyLossTest.cc

@@ -26 +26 @@
 //
 // $Id: muEnergyLossTest.cc,v 1.7 2006/06/29 19:49:56 gunter Exp $
-// GEANT4 tag $Name: geant4-09-03-cand-01 $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 //-----------------------------------------------------------------