Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/standard/src/G4hCoulombScatteringModel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 // -------------------------------------------------------------------
 27 //
 28 // GEANT4 Class file
 29 //
 30 //
 31 // File name:     G4hCoulombScatteringModel
 32 //
 33 // Author:        Vladimir Ivanchenko 
 34 //
 35 // Creation date: 08.06.2012 from G4eCoulombScatteringModel
 36 //
 37 // Modifications:
 38 //
 39 //
 40 // Class Description:
 41 //
 42 // -------------------------------------------------------------------
 43 //
 44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 46 
 47 #include "G4hCoulombScatteringModel.hh"
 48 #include "G4PhysicalConstants.hh"
 49 #include "G4SystemOfUnits.hh"
 50 #include "Randomize.hh"
 51 #include "G4DataVector.hh"
 52 #include "G4ElementTable.hh"
 53 #include "G4ParticleChangeForGamma.hh"
 54 #include "G4Proton.hh"
 55 #include "G4ParticleTable.hh"
 56 #include "G4IonTable.hh"
 57 #include "G4ProductionCutsTable.hh"
 58 #include "G4NucleiProperties.hh"
 59 #include "G4Pow.hh"
 60 #include "G4NistManager.hh"
 61 
 62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 63 
 64 G4hCoulombScatteringModel::G4hCoulombScatteringModel(G4bool combined)
 65   : G4VEmModel("hCoulombScattering"),
 66     cosThetaMin(1.0),
 67     cosThetaMax(-1.0),
 68     isCombined(combined)
 69 {
 70   fParticleChange = nullptr;
 71   fNistManager = G4NistManager::Instance();
 72   theIonTable  = G4ParticleTable::GetParticleTable()->GetIonTable();
 73   theProton    = G4Proton::Proton();
 74   currentMaterial = nullptr; 
 75   fixedCut = -1.0;
 76 
 77   pCuts = nullptr;
 78 
 79   recoilThreshold = 0.0; // by default does not work
 80 
 81   particle = nullptr;
 82   currentCouple = nullptr;
 83   wokvi = new G4WentzelVIRelXSection();
 84 
 85   currentMaterialIndex = 0;
 86   mass = CLHEP::proton_mass_c2;
 87   elecRatio = 0.0;
 88 }
 89 
 90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 91 
 92 G4hCoulombScatteringModel::~G4hCoulombScatteringModel()
 93 {
 94   delete wokvi;
 95 }
 96 
 97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 98 
 99 void G4hCoulombScatteringModel::Initialise(const G4ParticleDefinition* part,
100              const G4DataVector& cuts)
101 {
102   SetupParticle(part);
103   currentCouple = nullptr;
104 
105   // defined theta limit between single and multiple scattering 
106   isCombined = true;
107   G4double tet = PolarAngleLimit();
108 
109   if(tet <= 0.0) { 
110     cosThetaMin = 1.0; 
111     isCombined = false; 
112   } else if(tet >= CLHEP::pi) { 
113     cosThetaMin = -1.0; 
114   } else { 
115     cosThetaMin = std::cos(tet); 
116   }
117 
118   wokvi->Initialise(part, cosThetaMin);
119   /*
120   G4cout << "G4hCoulombScatteringModel: " << particle->GetParticleName()
121          << "  1-cos(ThetaLimit)= " << 1 - cosThetaMin
122    << "  cos(thetaMax)= " <<  cosThetaMax
123    << G4endl;
124   */
125   pCuts = &cuts;
126   //G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(3);
127   /*
128   G4cout << "!!! G4hCoulombScatteringModel::Initialise for " 
129      << part->GetParticleName() << "  cos(TetMin)= " << cosThetaMin 
130      << "  cos(TetMax)= " << cosThetaMax <<G4endl;
131   G4cout << "cut= " << (*pCuts)[0] << "  cut1= " << (*pCuts)[1] << G4endl;
132   */
133   if(!fParticleChange) {
134     fParticleChange = GetParticleChangeForGamma();
135   }
136   if(IsMaster() && mass < CLHEP::GeV && part->GetParticleName() != "GenericIon") {
137     InitialiseElementSelectors(part, cuts);
138   } 
139 }
140 
141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
142 
143 void G4hCoulombScatteringModel::InitialiseLocal(const G4ParticleDefinition*, 
144             G4VEmModel* masterModel)
145 {
146   SetElementSelectors(masterModel->GetElementSelectors());
147 }
148 
149 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
150 
151 G4double 
152 G4hCoulombScatteringModel::MinPrimaryEnergy(const G4Material* material,
153               const G4ParticleDefinition* part,
154               G4double)
155 {
156   SetupParticle(part);
157 
158   // define cut using cuts for proton
159   G4double cut = 
160     std::max(recoilThreshold, (*pCuts)[CurrentCouple()->GetIndex()]);
161 
162   // find out lightest element
163   const G4ElementVector* theElementVector = material->GetElementVector();
164   std::size_t nelm = material->GetNumberOfElements();
165 
166   // select lightest element
167   G4int Z = 300;
168   for (std::size_t j=0; j<nelm; ++j) {
169     Z = std::min(Z,(*theElementVector)[j]->GetZasInt());
170   }
171   G4int A = G4lrint(fNistManager->GetAtomicMassAmu(Z));
172   G4double targetMass = G4NucleiProperties::GetNuclearMass(A, Z);
173   G4double t = std::max(cut, 0.5*(cut + std::sqrt(2*cut*targetMass)));
174 
175   return t;
176 }
177 
178 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
179 
180 G4double G4hCoulombScatteringModel::ComputeCrossSectionPerAtom(
181                 const G4ParticleDefinition* p,
182     G4double kinEnergy,
183     G4double Z, G4double,
184     G4double cutEnergy, G4double)
185 {
186   //G4cout << "### G4hCoulombScatteringModel::ComputeCrossSectionPerAtom  for " 
187   //<< p->GetParticleName()<<" Z= "<<Z<<" e(MeV)= "<< kinEnergy/MeV << G4endl; 
188   G4double cross = 0.0;
189   elecRatio = 0.0;
190   if(p != particle) { SetupParticle(p); }
191 
192   // cross section is set to zero to avoid problems in sample secondary
193   if(kinEnergy <= 0.0) { return cross; }
194   DefineMaterial(CurrentCouple());
195 
196   G4int iz = G4lrint(Z);
197   G4double tmass = (1 == iz) ? proton_mass_c2 :
198     fNistManager->GetAtomicMassAmu(iz)*amu_c2; 
199   wokvi->SetTargetMass(tmass);
200 
201   G4double costmin = 
202     wokvi->SetupKinematic(kinEnergy, currentMaterial);
203 
204   if(cosThetaMax < costmin) {
205     G4double cut = (0.0 < fixedCut) ? fixedCut : cutEnergy;
206     costmin = wokvi->SetupTarget(iz, cut);
207     G4double costmax = 
208       (1 == iz && particle == theProton && cosThetaMax < 0.0) 
209       ? 0.0 : cosThetaMax; 
210     if(costmin > costmax) {
211       cross = wokvi->ComputeNuclearCrossSection(costmin, costmax)
212   + wokvi->ComputeElectronCrossSection(costmin, costmax);
213     }
214     /*  
215   if(p->GetParticleName() == "mu+") 
216   G4cout << "e(MeV)= " << kinEnergy/MeV << " cross(b)= " << cross/barn  
217    << " 1-costmin= " << 1-costmin
218    << " 1-costmax= " << 1-costmax
219    << " 1-cosThetaMax= " << 1-cosThetaMax
220    << G4endl;
221     */
222   }
223   return cross;  
224 }
225 
226 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
227 
228 void G4hCoulombScatteringModel::SampleSecondaries(
229                 std::vector<G4DynamicParticle*>* fvect,
230     const G4MaterialCutsCouple* couple,
231     const G4DynamicParticle* dp,
232     G4double cutEnergy,
233     G4double)
234 {
235   G4double kinEnergy = dp->GetKineticEnergy();
236   SetupParticle(dp->GetDefinition());
237   DefineMaterial(couple);
238 
239   // Choose nucleus
240   G4double cut = (0.0 < fixedCut) ? fixedCut : cutEnergy;
241 
242   const G4Element* elm = SelectRandomAtom(couple,particle,
243             kinEnergy,cut,kinEnergy);
244 
245   G4int iz = elm->GetZasInt();
246   G4int ia = SelectIsotopeNumber(elm);
247   G4double mass2 = G4NucleiProperties::GetNuclearMass(ia, iz);
248 
249   wokvi->SetTargetMass(mass2);
250   wokvi->SetupKinematic(kinEnergy, currentMaterial);
251   G4double costmin = wokvi->SetupTarget(iz, cut);
252   G4double costmax = (1 == iz && particle == theProton && cosThetaMax < 0.0) 
253     ? 0.0 :  cosThetaMax; 
254   if(costmin <= costmax) { return; }
255 
256   G4double cross = wokvi->ComputeNuclearCrossSection(costmin, costmax);
257   G4double ecross = wokvi->ComputeElectronCrossSection(costmin, costmax);
258   G4double ratio = ecross/(cross + ecross);
259 
260   G4ThreeVector newDirection = 
261     wokvi->SampleSingleScattering(costmin, costmax, ratio);
262 
263   // kinematics in the Lab system
264   G4double ptot = std::sqrt(kinEnergy*(kinEnergy + 2.0*mass));
265   G4double e1   = mass + kinEnergy;
266   
267   // Lab. system kinematics along projectile direction
268   G4LorentzVector v0 = G4LorentzVector(0, 0, ptot, e1+mass2);
269   G4LorentzVector v1 = G4LorentzVector(0, 0, ptot, e1);
270   G4ThreeVector bst = v0.boostVector();
271   v1.boost(-bst);
272   // CM projectile
273   G4double momCM = v1.pz(); 
274   
275   // Momentum after scattering of incident particle
276   v1.setX(momCM*newDirection.x());
277   v1.setY(momCM*newDirection.y());
278   v1.setZ(momCM*newDirection.z());
279 
280   // CM--->Lab
281   v1.boost(bst);
282 
283   G4ThreeVector dir = dp->GetMomentumDirection(); 
284   newDirection = v1.vect().unit();
285   newDirection.rotateUz(dir);   
286   
287   fParticleChange->ProposeMomentumDirection(newDirection);   
288   
289   // recoil
290   v0 -= v1; 
291   G4double trec = std::max(v0.e() - mass2, 0.0);
292   G4double edep = 0.0;
293 
294   G4double tcut = recoilThreshold;
295   if(pCuts) { tcut= std::max(tcut,(*pCuts)[currentMaterialIndex]); }
296  
297   if(trec > tcut) {
298     G4ParticleDefinition* ion = theIonTable->GetIon(iz, ia, 0);
299     newDirection = v0.vect().unit();
300     newDirection.rotateUz(dir);   
301     auto newdp = new G4DynamicParticle(ion, newDirection, trec);
302     fvect->push_back(newdp);
303   } else if(trec > 0.0) {
304     edep = trec;
305     fParticleChange->ProposeNonIonizingEnergyDeposit(edep);
306   }
307 
308   // finelize primary energy and energy balance
309   G4double finalT = v1.e() - mass;
310   if(finalT < 0.0) { 
311     edep += finalT;
312     finalT = 0.0;
313   } 
314   edep = std::max(edep, 0.0);
315   fParticleChange->SetProposedKineticEnergy(finalT);
316   fParticleChange->ProposeLocalEnergyDeposit(edep);
317 }
318 
319 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
320