ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4VLEPTSModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4VLEPTSModel.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 #include "G4VLEPTSModel.hh"
27 
28 #include "CLHEP/Units/SystemOfUnits.h"
29 
30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 G4VLEPTSModel::G4VLEPTSModel(const G4String& modelName) : G4VEmModel(modelName),isInitialised(false)
32 {
33  theMeanFreePathTable=NULL;
34 
35  theNumbBinTable=100;
36 
37  verboseLevel = 0;
38 }
39 
40 
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
42 G4VLEPTSModel::~G4VLEPTSModel()
43 {
44 
45  if(theMeanFreePathTable) {
46  theMeanFreePathTable->clearAndDestroy();
47  delete theMeanFreePathTable;
48  }
49 }
50 
51 
52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
53 void G4VLEPTSModel::Init()
54 {
55  theLowestEnergyLimit = 0.5*CLHEP::eV;
56  theHighestEnergyLimit = 1.0*CLHEP::MeV;
57  //t theHighestEnergyLimit = 15.0*CLHEP::MeV;
58  SetLowEnergyLimit(theLowestEnergyLimit);
59  SetHighEnergyLimit(theHighestEnergyLimit);
60  theNumbBinTable = 100;
61 
62 }
63 
64 
65 
66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
67 G4double G4VLEPTSModel::GetMeanFreePath(const G4Material* aMaterial,
68  const G4ParticleDefinition* ,
69  G4double kineticEnergy )
70 {
71 
72  G4double MeanFreePath;
73  G4bool isOutRange ;
74 
75  if( verboseLevel >= 3 ) G4cout << aMaterial->GetIndex() << " G4VLEPTSModel::GetMeanFreePath " << kineticEnergy << " > " << theHighestEnergyLimit << " < " << theLowestEnergyLimit << G4endl;
76  if (kineticEnergy > theHighestEnergyLimit || kineticEnergy < theLowestEnergyLimit)
77  MeanFreePath = DBL_MAX;
78  else
79  MeanFreePath = (*theMeanFreePathTable)(aMaterial->GetIndex())->
80  GetValue(kineticEnergy, isOutRange);
81 
82  return MeanFreePath;
83 }
84 
85 
86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
87 void G4VLEPTSModel::BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
88 {
89  //CHECK IF PATH VARIABLE IS DEFINED
90  char* path = std::getenv("G4LEDATA");
91  if( !path ) {
92  G4Exception("G4VLEPTSModel",
93  "",
94  FatalException,
95  "variable G4LEDATA not defined");
96  }
97 
98  // Build microscopic cross section table and mean free path table
99 
100  G4String aParticleName = aParticleType.GetParticleName();
101 
102  if (theMeanFreePathTable) {
103  theMeanFreePathTable->clearAndDestroy();
104  delete theMeanFreePathTable;
105  }
106 
107  theMeanFreePathTable = new G4PhysicsTable(G4Material::GetNumberOfMaterials());
108 
109  //LOOP TO MATERIALS IN GEOMETRY
110  const G4MaterialTable * materialTable = G4Material::GetMaterialTable() ;
111  std::vector<G4Material*>::const_iterator matite;
112  for( matite = materialTable->begin(); matite != materialTable->end(); matite++ ) {
113  const G4Material * aMaterial = (*matite);
114  G4String mateName = aMaterial->GetName();
115 
116  //READ PARAMETERS FOR THIS MATERIAL
117  std::string dirName = std::string(path) + "/lepts/";
118  std::string fnParam = dirName + mateName + "." + aParticleName + ".param.dat";
119  std::string baseName = std::string(path) + "/lepts/" + mateName + "." + aParticleName;
120  G4bool bData = ReadParam( fnParam, aMaterial );
121  if( !bData ) continue; // MATERIAL NOT EXISTING, DO NOT READ OTHER FILES
122 
123  //READ INTEGRAL CROSS SECTION FOR THIS MATERIAL
124  std::string fnIXS = baseName + ".IXS.dat";
125 
126  std::map< G4int, std::vector<G4double> > integralXS = ReadIXS(fnIXS, aMaterial);
127  if( verboseLevel >= 2 ) G4cout << GetName() << " : " << theXSType << " " << mateName << " INTEGRALXS " << integralXS.size() << G4endl;
128 
129  if( integralXS.size() == 0 ) {
130  G4cerr << " Integral cross sections will be set to 0. for material " << mateName << G4endl;
131  G4PhysicsLogVector* ptrVector = new G4PhysicsLogVector(theLowestEnergyLimit, theHighestEnergyLimit, 2);
132  ptrVector->PutValue(0, DBL_MAX);
133  ptrVector->PutValue(1, DBL_MAX);
134 
135  unsigned int matIdx = aMaterial->GetIndex();
136  theMeanFreePathTable->insertAt( matIdx , ptrVector ) ;
137 
138  } else {
139 
140  if( verboseLevel >= 2 ) {
141  std::map< G4int, std::vector<G4double> >::const_iterator itei;
142  for( itei = integralXS.begin(); itei != integralXS.end(); itei++ ){
143  G4cout << GetName() << " : " << (*itei).first << " INTEGRALXS NDATA " << (*itei).second.size() << G4endl;
144  }
145  }
146 
147  BuildMeanFreePathTable( aMaterial, integralXS );
148 
149  std::string fnDXS = baseName + ".DXS.dat";
150  std::string fnRMT = baseName + ".RMT.dat";
151  std::string fnEloss = baseName + ".Eloss.dat";
152  std::string fnEloss2 = baseName + ".Eloss2.dat";
153 
154  theDiffXS[aMaterial] = new G4LEPTSDiffXS(fnDXS);
155  if( !theDiffXS[aMaterial]->IsFileFound() ) {
156  G4Exception("G4VLEPTSModel::BuildPhysicsTable",
157  "",
158  FatalException,
159  (G4String("File not found :" + fnDXS).c_str()));
160  }
161 
162  theRMTDistr[aMaterial] = new G4LEPTSDistribution();
163  theRMTDistr[aMaterial]->ReadFile(fnRMT);
164 
165  theElostDistr[aMaterial] = new G4LEPTSElossDistr(fnEloss);
166  if( !theElostDistr[aMaterial]->IsFileFound() ) {
167  G4Exception("G4VLEPTSModel::BuildPhysicsTable",
168  "",
169  FatalException,
170  (G4String("File not found :" + fnEloss).c_str()));
171  }
172  }
173 
174  }
175 
176 }
177 
178 void G4VLEPTSModel::BuildMeanFreePathTable( const G4Material* aMaterial, std::map< G4int, std::vector<G4double> >& integralXS )
179 {
180  G4double LowEdgeEnergy, fValue;
181 
182  //BUILD MEAN FREE PATH TABLE FROM INTEGRAL CROSS SECTION
183  unsigned int matIdx = aMaterial->GetIndex();
184  G4PhysicsLogVector* ptrVector = new G4PhysicsLogVector(theLowestEnergyLimit, theHighestEnergyLimit, theNumbBinTable);
185 
186  for (G4int ii=0; ii < theNumbBinTable; ii++) {
187  LowEdgeEnergy = ptrVector->GetLowEdgeEnergy(ii);
188  if( verboseLevel >= 2 ) G4cout << GetName() << " " << ii << " LowEdgeEnergy " << LowEdgeEnergy << " > " << theLowestEnergyLimit << " < " << theHighestEnergyLimit << G4endl;
189  //- fValue = ComputeMFP(LowEdgeEnergy, material, aParticleName);
190  fValue = 0.;
191  if( LowEdgeEnergy >= theLowestEnergyLimit &&
192  LowEdgeEnergy <= theHighestEnergyLimit) {
193  G4double NbOfMoleculesPerVolume = aMaterial->GetDensity()/theMolecularMass[aMaterial]*CLHEP::Avogadro;
194 
195  G4double SIGMA = 0. ;
196  //- for ( size_t elm=0 ; elm < aMaterial->GetNumberOfElements() ; elm++ ) {
197  G4double crossSection = 0.;
198 
199  G4double eVEnergy = LowEdgeEnergy/CLHEP::eV;
200 
201  //- if( verboseLevel >= 2 ) G4cout << " eVEnergy " << eVEnergy << " LowEdgeE " << LowEdgeEnergy << " " << integralXS[theXSType][1] << G4endl;
202 
203  if(eVEnergy < integralXS[0][1] ) {
204  crossSection = 0.;
205  } else {
206  G4int Bin = 0; // locate bin
207  G4double aa, bb;
208  for( G4int jj=1; jj<theNXSdat[aMaterial]; jj++) { // Extrapolate for E > Emax !!!
209  if( verboseLevel >= 3 ) G4cout << " GET BIN " << jj << " "<< eVEnergy << " > " << integralXS[0][jj] << G4endl;
210  if( eVEnergy > integralXS[0][jj]) {
211  Bin = jj;
212  } else {
213  break;
214  }
215  }
216  aa = integralXS[0][Bin];
217  bb = integralXS[0][Bin+1];
218  crossSection = (integralXS[theXSType][Bin] + (integralXS[theXSType][Bin+1]-integralXS[theXSType][Bin])/(bb-aa)*(eVEnergy-aa) ) * 1.e-16*CLHEP::cm2;
219 
220  if( verboseLevel >= 3 ) G4cout << " crossSection " << crossSection << " " <<integralXS[theXSType][Bin] << " + " << (integralXS[theXSType][Bin+1]-integralXS[theXSType][Bin]) << " / " << (bb-aa) << " *" << (eVEnergy-aa) << " * " << 1.e-16*CLHEP::cm2 << G4endl;;
221 
222  // SIGMA += NbOfAtomsPerVolume[elm] * crossSection;
223  SIGMA = NbOfMoleculesPerVolume * crossSection;
224  if( verboseLevel >= 2 ) G4cout << GetName() << " ADDING SIGMA " << SIGMA << " NAtoms " << NbOfMoleculesPerVolume
225  << " Bin " << Bin << " TOTAL " << aa << " " << bb
226  << " XS " << integralXS[theXSType][Bin] << " " << integralXS[theXSType][Bin+1] << G4endl;
227  }
228 
229  //-}
230 
231  fValue = SIGMA > DBL_MIN ? 1./SIGMA : DBL_MAX;
232  }
233 
234  ptrVector->PutValue(ii, fValue);
235  if( verboseLevel >= 2 ) G4cout << GetName() << " BUILDXS " << ii << " : " << LowEdgeEnergy << " = " << fValue << G4endl;
236  }
237 
238  theMeanFreePathTable->insertAt( matIdx , ptrVector ) ;
239 }
240 
241 
242 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
243 G4double G4VLEPTSModel::SampleAngle(const G4Material* aMaterial, G4double e, G4double el)
244 {
245  G4double x;
246 
247  if( e < 10001) {
248  x = theDiffXS[aMaterial]->SampleAngleMT(e, el);
249  }
250  else {
251  G4double Ei = e; //incidente
252  G4double Ed = e -el; //dispersado
253 
254  G4double Pi = std::sqrt( std::pow( (Ei/27.2/137),2) +2*Ei/27.2); //incidente
255  G4double Pd = std::sqrt( std::pow( (Ed/27.2/137),2) +2*Ed/27.2); //dispersado
256 
257  G4double Kmin = Pi - Pd;
258  G4double Kmax = Pi + Pd;
259 
260  G4double KR = theRMTDistr[aMaterial]->Sample(Kmin, Kmax); //sorteo mom. transf.
261 
262  G4double co = (Pi*Pi + Pd*Pd - KR*KR) / (2*Pi*Pd); //cos ang. disp.
263  if( co > 1. ) co = 1.;
264  x = std::acos(co); //*360/twopi; //ang. dispers.
265  }
266  return(x);
267 }
268 
269 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
270 G4ThreeVector G4VLEPTSModel::SampleNewDirection(const G4Material* aMaterial, G4ThreeVector P0Dir, G4double e, G4double el) {
271 
272  G4double x = SampleAngle(aMaterial, e, el);
273 
274  G4double cosTeta = std::cos(x); //*twopi/360.0);
275  G4double sinTeta = std::sqrt(1.0-cosTeta*cosTeta);
276  G4double Phi = CLHEP::twopi * G4UniformRand() ;
277  G4double dirx = sinTeta*std::cos(Phi) , diry = sinTeta*std::sin(Phi) , dirz = cosTeta ;
278 
279  G4ThreeVector P1Dir(dirx, diry, dirz);
280 #ifdef DEBUG_LEPTS
281  if( verboseLevel >= 2 ) G4cout << " G4VLEPTSModel::SampleNewDirection " <<P1Dir << G4endl;
282 #endif
283  P1Dir.rotateUz(P0Dir);
284 #ifdef DEBUG_LEPTS
285  if( verboseLevel >= 2 ) G4cout << " G4VLEPTSModel::SampleNewDirection rotated " <<P1Dir << " " << P0Dir << G4endl;
286 #endif
287 
288  return(P1Dir);
289 }
290 
291 
292 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
293 G4ThreeVector G4VLEPTSModel::SampleNewDirection(G4ThreeVector P0Dir, G4double x)
294 {
295  G4double cosTeta = std::cos(x); //*twopi/360.0);
296  G4double sinTeta = std::sqrt(1.0-cosTeta*cosTeta);
297  G4double Phi = CLHEP::twopi * G4UniformRand() ;
298  G4double dirx = sinTeta*std::cos(Phi) , diry = sinTeta*std::sin(Phi) , dirz = cosTeta ;
299 
300  G4ThreeVector P1Dir( dirx,diry,dirz );
301  P1Dir.rotateUz(P0Dir);
302 
303  return(P1Dir);
304 }
305 
306 
307 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
308 G4double G4VLEPTSModel::SampleEnergyLoss(const G4Material* aMaterial, G4double eMin, G4double eMax)
309 {
310  G4double el;
311  el = theElostDistr[aMaterial]->Sample(eMin/CLHEP::eV, eMax/CLHEP::eV)*CLHEP::eV;
312 
313 #ifdef DEBUG_LEPTS
314  if( verboseLevel >= 2 ) G4cout << aMaterial->GetName() <<"SampleEnergyLoss/eV " << el/CLHEP::eV << " eMax/eV " << eMax/CLHEP::eV << " "
315  << " " << GetName() << G4endl;
316 #endif
317  return el;
318 }
319 
320 
321 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
322 G4bool G4VLEPTSModel::ReadParam(G4String fnParam, const G4Material* aMaterial )
323 {
324  std::ifstream fin(fnParam);
325  if (!fin.is_open()) {
326  G4Exception("G4VLEPTSModel::ReadParam",
327  "",
328  JustWarning,
329  (G4String("File not found: ")+ fnParam).c_str());
330  return false;
331  }
332 
333  G4double IonisPot, IonisPotInt;
334 
335  fin >> IonisPot >> IonisPotInt;
336  if( verboseLevel >= 1 ) G4cout << "Read param (" << fnParam << ")\t IonisPot: " << IonisPot
337  << " IonisPotInt: " << IonisPotInt << G4endl;
338 
339  theIonisPot[aMaterial] = IonisPot * CLHEP::eV;
340  theIonisPotInt[aMaterial] = IonisPotInt * CLHEP::eV;
341 
342  G4double MolecularMass = 0;
343  size_t nelem = aMaterial->GetNumberOfElements();
344  const G4int* atomsV = aMaterial->GetAtomsVector();
345  for( size_t ii = 0; ii < nelem; ii++ ) {
346  MolecularMass += aMaterial->GetElement(ii)->GetA()*atomsV[ii]/CLHEP::g;
347  // G4cout << " MMASS1 " << mmass/CLHEP::g << " " << aMaterial->GetElement(ii)->GetName() << " " << aMaterial->GetElement(ii)->GetA()/CLHEP::g << G4endl;
348  }
349  // G4cout << " MMASS " << MolecularMass << " " << MolecularMass*CLHEP::g << " ME " << mmass << " " << mmass/CLHEP::g << G4endl;
350  theMolecularMass[aMaterial] = MolecularMass* CLHEP::g/CLHEP::mole;
351  // theMolecularMass[aMaterial] = aMaterial->GetMassOfMolecule()*CLHEP::Avogadro; // Material mixtures do not calculate molecular mass
352 
353  if( verboseLevel >= 1) G4cout << " IonisPot: " << IonisPot/CLHEP::eV << " eV "
354  << " IonisPotInt: " << IonisPotInt/CLHEP::eV << " eV"
355  << " MolecularMass " << MolecularMass/(CLHEP::g/CLHEP::mole) << " g/mole" << G4endl;
356 
357  return true;
358 }
359 
360 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
361 std::map< G4int, std::vector<G4double> > G4VLEPTSModel::ReadIXS(G4String fnIXS, const G4Material* aMaterial )
362 {
363  std::map< G4int, std::vector<G4double> > integralXS; // process type - energy
364  //G4cout << "fnIXS (" << fnIXS << ")" << G4endl;
365 
366  std::ifstream fin(fnIXS);
367  if (!fin.is_open()) {
368  G4Exception("G4VLEPTSModel::ReadIXS",
369  "",
370  JustWarning,
371  (G4String("File not found: ")+ fnIXS).c_str());
372  return integralXS;
373  }
374 
375  G4int nXSdat, nXSsub;
376  fin >> nXSdat >> nXSsub;
377  if( verboseLevel >= 1 ) G4cout << "Read IXS (" << fnIXS << ")\t nXSdat: " << nXSdat
378  << " nXSsub: " << nXSsub << G4endl;
379  theNXSdat[aMaterial] = nXSdat;
380  theNXSsub[aMaterial] = nXSsub;
381 
382  G4double xsdat;
383  for (G4int ip=0; ip<=nXSsub; ip++) {
384  integralXS[ip].push_back(0.);
385  }
386  for (G4int ie=1; ie<=nXSdat; ie++) {
387  for (G4int ip=0; ip<=nXSsub; ip++) {
388  fin >> xsdat;
389  integralXS[ip].push_back(xsdat);
390  if( verboseLevel >= 3 ) G4cout << GetName() << " FILL IXS " << ip << " " << ie << " = " << integralXS[ip][ie] << " " << xsdat << G4endl;
391  // xsdat 1e-16*cm2
392  }
393  }
394  fin.close();
395 
396  return integralXS;
397 }
398