Changeset 1315 for trunk/source/processes/electromagnetic/muons/src
- Timestamp:
- Jun 18, 2010, 11:42:07 AM (14 years ago)
- Location:
- trunk/source/processes/electromagnetic/muons/src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/electromagnetic/muons/src/G4MuIonisation.cc
r1228 r1315 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4MuIonisation.cc,v 1. 59 2009/02/26 11:04:20 vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 3$26 // $Id: G4MuIonisation.cc,v 1.61 2010/06/04 09:30:40 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-04-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 89 89 #include "G4BohrFluctuations.hh" 90 90 #include "G4UnitsTable.hh" 91 #include "G4ICRU73QOModel.hh" 91 92 92 93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... … … 140 141 141 142 mass = theParticle->GetPDGMass(); 143 G4double q = theParticle->GetPDGCharge(); 144 G4double elow = 0.2*MeV; 142 145 SetSecondaryParticle(G4Electron::Electron()); 143 146 144 147 // Bragg peak model 145 if (!EmModel(1)) SetEmModel(new G4BraggModel(),1); 148 if (!EmModel(1)) { 149 if(q > 0.0) { SetEmModel(new G4BraggModel(),1); } 150 else { 151 SetEmModel(new G4ICRU73QOModel(),1); 152 //elow = 1.0*MeV; 153 } 154 } 146 155 EmModel(1)->SetLowEnergyLimit(MinKinEnergy()); 147 EmModel(1)->SetHighEnergyLimit( 0.2*MeV);156 EmModel(1)->SetHighEnergyLimit(elow); 148 157 AddEmModel(1, EmModel(1), new G4IonFluctuations()); 149 158 150 159 // high energy fluctuation model 151 if (!FluctModel()) SetFluctModel(new G4UniversalFluctuation());160 if (!FluctModel()) { SetFluctModel(new G4UniversalFluctuation()); } 152 161 153 162 // moderate energy model 154 if (!EmModel(2)) SetEmModel(new G4BetheBlochModel(),2);155 EmModel(2)->SetLowEnergyLimit( 0.2*MeV);163 if (!EmModel(2)) { SetEmModel(new G4BetheBlochModel(),2); } 164 EmModel(2)->SetLowEnergyLimit(elow); 156 165 EmModel(2)->SetHighEnergyLimit(1.0*GeV); 157 166 AddEmModel(2, EmModel(2), FluctModel()); 158 167 159 168 // high energy model 160 if (!EmModel(3)) SetEmModel(new G4MuBetheBlochModel(),3);169 if (!EmModel(3)) { SetEmModel(new G4MuBetheBlochModel(),3); } 161 170 EmModel(3)->SetLowEnergyLimit(1.0*GeV); 162 171 EmModel(3)->SetHighEnergyLimit(MaxKinEnergy()); -
trunk/source/processes/electromagnetic/muons/src/G4MuPairProductionModel.cc
r1228 r1315 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4MuPairProductionModel.cc,v 1.4 4 2009/08/11 16:50:07vnivanch Exp $27 // GEANT4 tag $Name: geant4-09-0 3$26 // $Id: G4MuPairProductionModel.cc,v 1.45 2010/06/01 15:21:59 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-04-beta-cand-01 $ 28 28 // 29 29 // ------------------------------------------------------------------- … … 466 466 { 467 467 G4double kineticEnergy = aDynamicParticle->GetKineticEnergy(); 468 G4double totalEnergy = kineticEnergy + particleMass ; 469 G4ParticleMomentum ParticleDirection = 470 aDynamicParticle->GetMomentumDirection(); 468 G4double totalEnergy = kineticEnergy + particleMass; 469 G4double totalMomentum = 470 sqrt(kineticEnergy*(kineticEnergy + 2.0*particleMass)); 471 472 G4ThreeVector partDirection = aDynamicParticle->GetMomentumDirection(); 471 473 472 474 G4int it; 473 for(it=1; it<ntdat; it++) {if(kineticEnergy <= tdat[it]) break;}474 if(it == ntdat) it--;475 for(it=1; it<ntdat; ++it) {if(kineticEnergy <= tdat[it]) break;} 476 if(it == ntdat) { --it; } 475 477 G4double dt = log(kineticEnergy/tdat[it-1])/log(tdat[it]/tdat[it-1]); 476 478 477 479 // select randomly one element constituing the material 478 const G4Element* anElement = SelectRandomAtom(kineticEnergy, dt, it, couple, tmin); 480 const G4Element* anElement = 481 SelectRandomAtom(kineticEnergy, dt, it, couple, tmin); 479 482 SetCurrentElement(anElement->GetZ()); 480 483 … … 484 487 G4double minEnergy = std::max(tmin, minPairEnergy); 485 488 486 if(minEnergy >= maxEnergy) return;489 if(minEnergy >= maxEnergy) { return; } 487 490 //G4cout << "emin= " << minEnergy << " emax= " << maxEnergy 488 491 // << " minPair= " << minPairEnergy << " maxpair= " << maxPairEnergy … … 507 510 G4int iz, iy; 508 511 509 for(iz=1; iz<nzdat; iz++) {if(currentZ <= zdat[iz]) break;}510 if(iz == nzdat) iz--;512 for(iz=1; iz<nzdat; ++iz) { if(currentZ <= zdat[iz]) { break; } } 513 if(iz == nzdat) { --iz; } 511 514 512 515 G4double dz = log(currentZ/zdat[iz-1])/log(zdat[iz]/zdat[iz-1]); … … 520 523 G4double p1 = pmin; 521 524 G4double p2 = pmin; 522 for(iy=iymin+1; iy<=iymax; iy++) {525 for(iy=iymin+1; iy<=iymax; ++iy) { 523 526 p1 = p2; 524 527 p2 = InterpolatedIntegralCrossSection(dt, dz, iz, it, iy, currentZ); … … 532 535 *log(maxPairEnergy/minPairEnergy)); 533 536 534 if(PairEnergy < minEnergy) PairEnergy = minEnergy;535 if(PairEnergy > maxEnergy) PairEnergy = maxEnergy;537 if(PairEnergy < minEnergy) { PairEnergy = minEnergy; } 538 if(PairEnergy > maxEnergy) { PairEnergy = maxEnergy; } 536 539 537 540 // sample r=(E+-E-)/PairEnergy ( uniformly .....) … … 545 548 G4double PositronEnergy = PairEnergy - ElectronEnergy; 546 549 547 // angles of the emitted particles ( Z - axis along the parent particle) 548 // (mean theta for the moment) 549 550 // 551 // scattered electron (positron) angles. ( Z - axis along the parent photon) 552 // 553 // universal distribution suggested by L. Urban 554 // (Geant3 manual (1993) Phys211), 555 // derived from Tsai distribution (Rev Mod Phys 49,421(1977)) 556 // G4cout << "Ee= " << ElectronEnergy << " Ep= " << PositronEnergy << G4endl; 557 G4double u; 558 const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ; 559 560 if (9./(9.+d) >G4UniformRand()) u= - log(G4UniformRand()*G4UniformRand())/a1; 561 else u= - log(G4UniformRand()*G4UniformRand())/a2; 562 563 G4double TetEl = u*electron_mass_c2/ElectronEnergy; 564 G4double TetPo = u*electron_mass_c2/PositronEnergy; 565 G4double Phi = twopi * G4UniformRand(); 566 G4double dxEl= sin(TetEl)*cos(Phi),dyEl= sin(TetEl)*sin(Phi),dzEl=cos(TetEl); 567 G4double dxPo=-sin(TetPo)*cos(Phi),dyPo=-sin(TetPo)*sin(Phi),dzPo=cos(TetPo); 568 569 G4ThreeVector ElectDirection (dxEl, dyEl, dzEl); 570 ElectDirection.rotateUz(ParticleDirection); 550 // The angle of the emitted virtual photon is sampled 551 // according to the muon bremsstrahlung model 552 553 G4double gam = totalEnergy/particleMass; 554 G4double gmax = gam*std::min(1.0, totalEnergy/PairEnergy - 1.0); 555 G4double gmax2= gmax*gmax; 556 G4double x = G4UniformRand()*gmax2/(1.0 + gmax2); 557 558 G4double theta = sqrt(x/(1.0 - x))/gam; 559 G4double sint = sin(theta); 560 G4double phi = twopi * G4UniformRand() ; 561 G4double dirx = sint*cos(phi), diry = sint*sin(phi), dirz = cos(theta) ; 562 563 G4ThreeVector gDirection(dirx, diry, dirz); 564 gDirection.rotateUz(partDirection); 565 566 // the angles of e- and e+ assumed to be the same as virtual gamma 571 567 572 568 // create G4DynamicParticle object for the particle1 573 G4DynamicParticle* aParticle1= new G4DynamicParticle(theElectron, 574 ElectDirection, 575 ElectronEnergy - electron_mass_c2); 576 577 G4ThreeVector PositDirection (dxPo, dyPo, dzPo); 578 PositDirection.rotateUz(ParticleDirection); 569 G4DynamicParticle* aParticle1 = 570 new G4DynamicParticle(theElectron, gDirection, 571 ElectronEnergy - electron_mass_c2); 579 572 580 573 // create G4DynamicParticle object for the particle2 581 574 G4DynamicParticle* aParticle2 = 582 new G4DynamicParticle(thePositron, 583 PositDirection, 575 new G4DynamicParticle(thePositron, gDirection, 584 576 PositronEnergy - electron_mass_c2); 585 577 … … 588 580 fParticleChange->SetProposedKineticEnergy(kineticEnergy); 589 581 582 partDirection *= totalMomentum; 583 partDirection -= (aParticle1->GetMomentum() + aParticle2->GetMomentum()); 584 partDirection = partDirection.unit(); 585 fParticleChange->SetProposedMomentumDirection(partDirection); 586 587 // add secondary 590 588 vdp->push_back(aParticle1); 591 589 vdp->push_back(aParticle2);
Note: See TracChangeset
for help on using the changeset viewer.