Changeset 1315 for trunk/source/processes/electromagnetic/muons/src/G4MuPairProductionModel.cc

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

File:

• trunk/source/processes/electromagnetic/muons/src/G4MuPairProductionModel.cc (modified) (8 diffs)

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);