ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4eCoulombScatteringModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4eCoulombScatteringModel.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 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4eCoulombScatteringModel
33 //
34 // Author: Vladimir Ivanchenko
35 //
36 // Creation date: 22.08.2005
37 //
38 // Modifications: V.Ivanchenko
39 //
40 //
41 //
42 // Class Description:
43 //
44 // -------------------------------------------------------------------
45 //
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48 
49 #include "G4eCoulombScatteringModel.hh"
50 #include "G4PhysicalConstants.hh"
51 #include "G4SystemOfUnits.hh"
52 #include "Randomize.hh"
53 #include "G4DataVector.hh"
54 #include "G4ElementTable.hh"
55 #include "G4ParticleChangeForGamma.hh"
56 #include "G4Proton.hh"
57 #include "G4ParticleTable.hh"
58 #include "G4IonTable.hh"
59 #include "G4ProductionCutsTable.hh"
60 #include "G4NucleiProperties.hh"
61 #include "G4Pow.hh"
62 #include "G4NistManager.hh"
63 
64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
65 
66 using namespace std;
67 
68 G4eCoulombScatteringModel::G4eCoulombScatteringModel(G4bool combined)
69  : G4VEmModel("eCoulombScattering"),
70  cosThetaMin(1.0),
71  cosThetaMax(-1.0),
72  isCombined(combined)
73 {
74  fParticleChange = nullptr;
75  fNistManager = G4NistManager::Instance();
76  theIonTable = G4ParticleTable::GetParticleTable()->GetIonTable();
77  theProton = G4Proton::Proton();
78  currentMaterial = nullptr;
79  fixedCut = -1.0;
80 
81  pCuts = nullptr;
82 
83  recoilThreshold = 0.0; // by default does not work
84 
85  particle = nullptr;
86  currentCouple = nullptr;
87 
88  wokvi = new G4WentzelOKandVIxSection(isCombined);
89 
90  currentMaterialIndex = 0;
91  mass = CLHEP::proton_mass_c2;
92  elecRatio = 0.0;
93 }
94 
95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
96 
97 G4eCoulombScatteringModel::~G4eCoulombScatteringModel()
98 {
99  delete wokvi;
100 }
101 
102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
103 
104 void G4eCoulombScatteringModel::Initialise(const G4ParticleDefinition* part,
105  const G4DataVector& cuts)
106 {
107  SetupParticle(part);
108  currentCouple = nullptr;
109 
110  // defined theta limit between single and multiple scattering
111  if(isCombined) {
112  cosThetaMin = 1.0;
113  G4double tet = PolarAngleLimit();
114  if(tet >= pi) { cosThetaMin = -1.0; }
115  else if(tet > 0.0) { cosThetaMin = cos(tet); }
116  }
117 
118  wokvi->Initialise(part, cosThetaMin);
119  pCuts = &cuts;
120  /*
121  G4cout << "G4eCoulombScatteringModel::Initialise for "
122  << part->GetParticleName() << " 1-cos(TetMin)= " << 1.0 - cosThetaMin
123  << " 1-cos(TetMax)= " << 1. - cosThetaMax << G4endl;
124  G4cout << "cut[0]= " << (*pCuts)[0] << G4endl;
125  */
126  if(!fParticleChange) {
127  fParticleChange = GetParticleChangeForGamma();
128  }
129  if(IsMaster() && mass < GeV && part->GetParticleName() != "GenericIon") {
130  InitialiseElementSelectors(part, cuts);
131  }
132 }
133 
134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
135 
136 void G4eCoulombScatteringModel::InitialiseLocal(const G4ParticleDefinition*,
137  G4VEmModel* masterModel)
138 {
139  SetElementSelectors(masterModel->GetElementSelectors());
140 }
141 
142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
143 
144 G4double
145 G4eCoulombScatteringModel::MinPrimaryEnergy(const G4Material* material,
146  const G4ParticleDefinition* part,
147  G4double)
148 {
149  SetupParticle(part);
150 
151  // define cut using cuts for proton
152  G4double cut =
153  std::max(recoilThreshold, (*pCuts)[CurrentCouple()->GetIndex()]);
154 
155  // find out lightest element
156  const G4ElementVector* theElementVector = material->GetElementVector();
157  G4int nelm = material->GetNumberOfElements();
158 
159  // select lightest element
160  G4int Z = 300;
161  for (G4int j=0; j<nelm; ++j) {
162  Z = std::min(Z,(*theElementVector)[j]->GetZasInt());
163  }
164  G4int A = G4lrint(fNistManager->GetAtomicMassAmu(Z));
165  G4double targetMass = G4NucleiProperties::GetNuclearMass(A, Z);
166  G4double t = std::max(cut, 0.5*(cut + sqrt(2*cut*targetMass)));
167 
168  return t;
169 }
170 
171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
172 
173 G4double G4eCoulombScatteringModel::ComputeCrossSectionPerAtom(
174  const G4ParticleDefinition* p,
175  G4double kinEnergy,
176  G4double Z, G4double,
177  G4double cutEnergy, G4double)
178 {
179  /*
180  G4cout << "### G4eCoulombScatteringModel::ComputeCrossSectionPerAtom for "
181  << p->GetParticleName()<<" Z= "<<Z<<" e(MeV)= "<< kinEnergy/MeV
182  << G4endl;
183  */
184  G4double cross = 0.0;
185  elecRatio = 0.0;
186  if(p != particle) { SetupParticle(p); }
187 
188  // cross section is set to zero to avoid problems in sample secondary
189  if(kinEnergy <= 0.0) { return cross; }
190  DefineMaterial(CurrentCouple());
191  G4double costmin = wokvi->SetupKinematic(kinEnergy, currentMaterial);
192 
193  //G4cout << "cosThetaMax= "<<cosThetaMax<<" costmin= "<<costmin<< G4endl;
194 
195  if(cosThetaMax < costmin) {
196  G4int iz = G4lrint(Z);
197  G4double cut = (0.0 < fixedCut) ? fixedCut : cutEnergy;
198  costmin = wokvi->SetupTarget(iz, cut);
199  //G4cout << "SetupTarget: Z= " << iz << " cut= " << cut << " "
200  // << costmin << G4endl;
201  G4double costmax = (1 == iz && particle == theProton && cosThetaMax < 0.0)
202  ? 0.0 : cosThetaMax;
203  if(costmin > costmax) {
204  cross = wokvi->ComputeNuclearCrossSection(costmin, costmax)
205  + wokvi->ComputeElectronCrossSection(costmin, costmax);
206  }
207  /*
208  if(p->GetParticleName() == "e-")
209  G4cout << "Z= " << Z << " e(MeV)= " << kinEnergy/MeV
210  << " cross(b)= " << cross/barn << " 1-costmin= " << 1-costmin
211  << " 1-costmax= " << 1-costmax
212  << " 1-cosThetaMax= " << 1-cosThetaMax
213  << " " << currentMaterial->GetName()
214  << G4endl;
215  */
216  }
217  //G4cout << "====== cross= " << cross << G4endl;
218  return cross;
219 }
220 
221 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
222 
223 void G4eCoulombScatteringModel::SampleSecondaries(
224  std::vector<G4DynamicParticle*>* fvect,
225  const G4MaterialCutsCouple* couple,
226  const G4DynamicParticle* dp,
227  G4double cutEnergy,
228  G4double)
229 {
230  G4double kinEnergy = dp->GetKineticEnergy();
231  SetupParticle(dp->GetDefinition());
232  DefineMaterial(couple);
233  /*
234  G4cout << "G4eCoulombScatteringModel::SampleSecondaries e(MeV)= "
235  << kinEnergy << " " << particle->GetParticleName()
236  << " cut= " << cutEnergy<< G4endl;
237  */
238  // Choose nucleus
239  G4double cut = (0.0 < fixedCut) ? fixedCut : cutEnergy;
240 
241  wokvi->SetupKinematic(kinEnergy, currentMaterial);
242 
243  const G4Element* currentElement = SelectTargetAtom(couple,particle,kinEnergy,
244  dp->GetLogKineticEnergy(),cut,kinEnergy);
245  G4int iz = currentElement->GetZasInt();
246 
247  G4double costmin = wokvi->SetupTarget(iz, cut);
248  G4double costmax = (1 == iz && particle == theProton && cosThetaMax < 0.0)
249  ? 0.0 : cosThetaMax;
250  if(costmin <= costmax) { return; }
251 
252  G4double cross = wokvi->ComputeNuclearCrossSection(costmin, costmax);
253  G4double ecross = wokvi->ComputeElectronCrossSection(costmin, costmax);
254  G4double ratio = ecross/(cross + ecross);
255 
256  G4int ia = SelectIsotopeNumber(currentElement);
257  G4double targetMass = G4NucleiProperties::GetNuclearMass(ia, iz);
258  wokvi->SetTargetMass(targetMass);
259 
260  G4ThreeVector newDirection =
261  wokvi->SampleSingleScattering(costmin, costmax, ratio);
262  G4double cost = newDirection.z();
263  /*
264  G4cout << "SampleSec: e(MeV)= " << kinEnergy/MeV
265  << " 1-costmin= " << 1-costmin
266  << " 1-costmax= " << 1-costmax
267  << " 1-cost= " << 1-cost
268  << " ratio= " << ratio
269  << G4endl;
270  */
271  G4ThreeVector direction = dp->GetMomentumDirection();
272  newDirection.rotateUz(direction);
273 
274  fParticleChange->ProposeMomentumDirection(newDirection);
275 
276  // recoil sampling assuming a small recoil
277  // and first order correction to primary 4-momentum
278  G4double mom2 = wokvi->GetMomentumSquare();
279  G4double trec = mom2*(1.0 - cost)
280  /(targetMass + (mass + kinEnergy)*(1.0 - cost));
281 
282  // the check likely not needed
283  trec = std::min(trec, kinEnergy);
284  G4double finalT = kinEnergy - trec;
285  G4double edep = 0.0;
286  /*
287  G4cout<<"G4eCoulombScatteringModel: finalT= "<<finalT<<" Trec= "
288  <<trec << " Z= " << iz << " A= " << ia
289  << " tcut(keV)= " << (*pCuts)[currentMaterialIndex]/keV << G4endl;
290  */
291  G4double tcut = recoilThreshold;
292  if(pCuts) { tcut= std::max(tcut,(*pCuts)[currentMaterialIndex]); }
293 
294  if(trec > tcut) {
295  G4ParticleDefinition* ion = theIonTable->GetIon(iz, ia, 0);
296  G4ThreeVector dir = (direction*sqrt(mom2) -
297  newDirection*sqrt(finalT*(2*mass + finalT))).unit();
298  G4DynamicParticle* newdp = new G4DynamicParticle(ion, dir, trec);
299  fvect->push_back(newdp);
300  } else {
301  edep = trec;
302  fParticleChange->ProposeNonIonizingEnergyDeposit(edep);
303  }
304 
305  // finelize primary energy and energy balance
306  // this threshold may be applied only because for low-enegry
307  // e+e- msc model is applied
308  if(finalT < 0.0) {
309  edep += finalT;
310  finalT = 0.0;
311  }
312  edep = std::max(edep, 0.0);
313  fParticleChange->SetProposedKineticEnergy(finalT);
314  fParticleChange->ProposeLocalEnergyDeposit(edep);
315 }
316 
317 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......