ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4hCoulombScatteringModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4hCoulombScatteringModel.cc
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 using namespace std;
65 
66 G4hCoulombScatteringModel::G4hCoulombScatteringModel(G4bool combined)
67  : G4VEmModel("hCoulombScattering"),
68  cosThetaMin(1.0),
69  cosThetaMax(-1.0),
70  isCombined(combined)
71 {
72  fParticleChange = nullptr;
73  fNistManager = G4NistManager::Instance();
74  theIonTable = G4ParticleTable::GetParticleTable()->GetIonTable();
75  theProton = G4Proton::Proton();
76  currentMaterial = nullptr;
77  fixedCut = -1.0;
78 
79  pCuts = nullptr;
80 
81  recoilThreshold = 0.0; // by default does not work
82 
83  particle = nullptr;
84  currentCouple = nullptr;
85  wokvi = new G4WentzelVIRelXSection();
86 
87  currentMaterialIndex = 0;
88  mass = CLHEP::proton_mass_c2;
89  elecRatio = 0.0;
90 }
91 
92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
93 
94 G4hCoulombScatteringModel::~G4hCoulombScatteringModel()
95 {
96  delete wokvi;
97 }
98 
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100 
101 void G4hCoulombScatteringModel::Initialise(const G4ParticleDefinition* part,
102  const G4DataVector& cuts)
103 {
104  SetupParticle(part);
105  currentCouple = nullptr;
106 
107  // defined theta limit between single and multiple scattering
108  isCombined = true;
109  G4double tet = PolarAngleLimit();
110 
111  if(tet <= 0.0) {
112  cosThetaMin = 1.0;
113  isCombined = false;
114  } else if(tet >= CLHEP::pi) {
115  cosThetaMin = -1.0;
116  } else {
117  cosThetaMin = cos(tet);
118  }
119 
120  wokvi->Initialise(part, cosThetaMin);
121  /*
122  G4cout << "G4hCoulombScatteringModel: " << particle->GetParticleName()
123  << " 1-cos(ThetaLimit)= " << 1 - cosThetaMin
124  << " cos(thetaMax)= " << cosThetaMax
125  << G4endl;
126  */
127  pCuts = &cuts;
128  //G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(3);
129  /*
130  G4cout << "!!! G4hCoulombScatteringModel::Initialise for "
131  << part->GetParticleName() << " cos(TetMin)= " << cosThetaMin
132  << " cos(TetMax)= " << cosThetaMax <<G4endl;
133  G4cout << "cut= " << (*pCuts)[0] << " cut1= " << (*pCuts)[1] << G4endl;
134  */
135  if(!fParticleChange) {
136  fParticleChange = GetParticleChangeForGamma();
137  }
138  if(IsMaster() && mass < GeV && part->GetParticleName() != "GenericIon") {
139  InitialiseElementSelectors(part, cuts);
140  }
141 }
142 
143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
144 
145 void G4hCoulombScatteringModel::InitialiseLocal(const G4ParticleDefinition*,
146  G4VEmModel* masterModel)
147 {
148  SetElementSelectors(masterModel->GetElementSelectors());
149 }
150 
151 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
152 
153 G4double
154 G4hCoulombScatteringModel::MinPrimaryEnergy(const G4Material* material,
155  const G4ParticleDefinition* part,
156  G4double)
157 {
158  SetupParticle(part);
159 
160  // define cut using cuts for proton
161  G4double cut =
162  std::max(recoilThreshold, (*pCuts)[CurrentCouple()->GetIndex()]);
163 
164  // find out lightest element
165  const G4ElementVector* theElementVector = material->GetElementVector();
166  G4int nelm = material->GetNumberOfElements();
167 
168  // select lightest element
169  G4int Z = 300;
170  for (G4int j=0; j<nelm; ++j) {
171  Z = std::min(Z,(*theElementVector)[j]->GetZasInt());
172  }
173  G4int A = G4lrint(fNistManager->GetAtomicMassAmu(Z));
174  G4double targetMass = G4NucleiProperties::GetNuclearMass(A, Z);
175  G4double t = std::max(cut, 0.5*(cut + sqrt(2*cut*targetMass)));
176 
177  return t;
178 }
179 
180 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
181 
182 G4double G4hCoulombScatteringModel::ComputeCrossSectionPerAtom(
183  const G4ParticleDefinition* p,
184  G4double kinEnergy,
185  G4double Z, G4double,
186  G4double cutEnergy, G4double)
187 {
188  //G4cout << "### G4hCoulombScatteringModel::ComputeCrossSectionPerAtom for "
189  //<< p->GetParticleName()<<" Z= "<<Z<<" e(MeV)= "<< kinEnergy/MeV << G4endl;
190  G4double cross = 0.0;
191  elecRatio = 0.0;
192  if(p != particle) { SetupParticle(p); }
193 
194  // cross section is set to zero to avoid problems in sample secondary
195  if(kinEnergy <= 0.0) { return cross; }
196  DefineMaterial(CurrentCouple());
197 
198  G4int iz = G4lrint(Z);
199  G4double tmass = (1 == iz) ? proton_mass_c2 :
200  fNistManager->GetAtomicMassAmu(iz)*amu_c2;
201  wokvi->SetTargetMass(tmass);
202 
203  G4double costmin =
204  wokvi->SetupKinematic(kinEnergy, currentMaterial);
205 
206  if(cosThetaMax < costmin) {
207  G4double cut = (0.0 < fixedCut) ? fixedCut : cutEnergy;
208  costmin = wokvi->SetupTarget(iz, cut);
209  G4double costmax =
210  (1 == iz && particle == theProton && cosThetaMax < 0.0)
211  ? 0.0 : cosThetaMax;
212  if(costmin > costmax) {
213  cross = wokvi->ComputeNuclearCrossSection(costmin, costmax)
214  + wokvi->ComputeElectronCrossSection(costmin, costmax);
215  }
216  /*
217  if(p->GetParticleName() == "mu+")
218  G4cout << "e(MeV)= " << kinEnergy/MeV << " cross(b)= " << cross/barn
219  << " 1-costmin= " << 1-costmin
220  << " 1-costmax= " << 1-costmax
221  << " 1-cosThetaMax= " << 1-cosThetaMax
222  << G4endl;
223  */
224  }
225  return cross;
226 }
227 
228 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
229 
230 void G4hCoulombScatteringModel::SampleSecondaries(
231  std::vector<G4DynamicParticle*>* fvect,
232  const G4MaterialCutsCouple* couple,
233  const G4DynamicParticle* dp,
234  G4double cutEnergy,
235  G4double)
236 {
237  G4double kinEnergy = dp->GetKineticEnergy();
238  SetupParticle(dp->GetDefinition());
239  DefineMaterial(couple);
240 
241  // Choose nucleus
242  G4double cut = (0.0 < fixedCut) ? fixedCut : cutEnergy;
243 
244  const G4Element* elm = SelectRandomAtom(couple,particle,
245  kinEnergy,cut,kinEnergy);
246 
247  G4int iz = elm->GetZasInt();
248  G4int ia = SelectIsotopeNumber(elm);
249  G4double mass2 = G4NucleiProperties::GetNuclearMass(ia, iz);
250 
251  wokvi->SetTargetMass(mass2);
252  wokvi->SetupKinematic(kinEnergy, currentMaterial);
253  G4double costmin = wokvi->SetupTarget(iz, cut);
254  G4double costmax = (1 == iz && particle == theProton && cosThetaMax < 0.0)
255  ? 0.0 : cosThetaMax;
256  if(costmin <= costmax) { return; }
257 
258  G4double cross = wokvi->ComputeNuclearCrossSection(costmin, costmax);
259  G4double ecross = wokvi->ComputeElectronCrossSection(costmin, costmax);
260  G4double ratio = ecross/(cross + ecross);
261 
262  G4ThreeVector newDirection =
263  wokvi->SampleSingleScattering(costmin, costmax, ratio);
264 
265  // kinematics in the Lab system
266  G4double ptot = sqrt(kinEnergy*(kinEnergy + 2.0*mass));
267  G4double e1 = mass + kinEnergy;
268 
269  // Lab. system kinematics along projectile direction
270  G4LorentzVector v0 = G4LorentzVector(0, 0, ptot, e1+mass2);
271  G4LorentzVector v1 = G4LorentzVector(0, 0, ptot, e1);
272  G4ThreeVector bst = v0.boostVector();
273  v1.boost(-bst);
274  // CM projectile
275  G4double momCM = v1.pz();
276 
277  // Momentum after scattering of incident particle
278  v1.setX(momCM*newDirection.x());
279  v1.setY(momCM*newDirection.y());
280  v1.setZ(momCM*newDirection.z());
281 
282  // CM--->Lab
283  v1.boost(bst);
284 
285  G4ThreeVector dir = dp->GetMomentumDirection();
286  newDirection = v1.vect().unit();
287  newDirection.rotateUz(dir);
288 
289  fParticleChange->ProposeMomentumDirection(newDirection);
290 
291  // recoil
292  v0 -= v1;
293  G4double trec = std::max(v0.e() - mass2, 0.0);
294  G4double edep = 0.0;
295 
296  G4double tcut = recoilThreshold;
297  if(pCuts) { tcut= std::max(tcut,(*pCuts)[currentMaterialIndex]); }
298 
299  if(trec > tcut) {
300  G4ParticleDefinition* ion = theIonTable->GetIon(iz, ia, 0);
301  newDirection = v0.vect().unit();
302  newDirection.rotateUz(dir);
303  G4DynamicParticle* newdp = new G4DynamicParticle(ion, newDirection, trec);
304  fvect->push_back(newdp);
305  } else if(trec > 0.0) {
306  edep = trec;
307  fParticleChange->ProposeNonIonizingEnergyDeposit(edep);
308  }
309 
310  // finelize primary energy and energy balance
311  G4double finalT = v1.e() - mass;
312  if(finalT < 0.0) {
313  edep += finalT;
314  finalT = 0.0;
315  }
316  edep = std::max(edep, 0.0);
317  fParticleChange->SetProposedKineticEnergy(finalT);
318  fParticleChange->ProposeLocalEnergyDeposit(edep);
319 }
320 
321 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......