Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
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 
29 
30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 G4VLEPTSModel::G4VLEPTSModel(const G4String& modelName) : G4VEmModel(modelName),isInitialised(false)
32 {
34 
35  theNumbBinTable=100;
36 
37  verboseLevel = 0;
38 }
39 
40 
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
43 {
44 
47  delete theMeanFreePathTable;
48  }
49 }
50 
51 
52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
54 {
57  //t theHighestEnergyLimit = 15.0*CLHEP::MeV;
60  theNumbBinTable = 100;
61 
62 }
63 
64 
65 
66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
68  const G4ParticleDefinition* ,
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......
88 {
89  //CHECK IF PATH VARIABLE IS DEFINED
90  char* path = getenv("G4LEDATA");
91  if( !path ) {
92  G4Exception("G4VLEPTSModel",
93  "",
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) {
104  delete theMeanFreePathTable;
105  }
106 
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;
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  "",
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  "",
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();
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......
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......
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);
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......
294 {
295  G4double cosTeta = std::cos(x); //*twopi/360.0);
296  G4double sinTeta = std::sqrt(1.0-cosTeta*cosTeta);
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......
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......
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 
Float_t x
Definition: compare.C:6
static size_t GetNumberOfMaterials()
Definition: G4Material.cc:600
std::map< const G4Material *, G4LEPTSDistribution * > theRMTDistr
G4double GetLowEdgeEnergy(size_t binNumber) const
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:707
std::map< const G4Material *, G4LEPTSElossDistr * > theElostDistr
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:593
G4double SampleAngle(const G4Material *aMaterial, G4double e, G4double el)
static constexpr double g
size_t GetIndex() const
Definition: G4Material.hh:262
#define G4endl
Definition: G4ios.hh:61
static constexpr double mole
const G4String & GetParticleName() const
G4PhysicsTable * theMeanFreePathTable
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
void insertAt(size_t, G4PhysicsVector *)
const G4int * GetAtomsVector() const
Definition: G4Material.hh:199
G4double GetMeanFreePath(const G4Material *mate, const G4ParticleDefinition *aParticle, G4double kineticEnergy)
void BuildPhysicsTable(const G4ParticleDefinition &aParticleType)
G4VLEPTSModel(const G4String &processName)
std::map< const G4Material *, G4double > theIonisPotInt
fin
Definition: AddMC.C:2
const G4String & GetName() const
Definition: G4Material.hh:179
static ulg bb
Definition: csz_inflate.cc:344
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
static constexpr double Avogadro
virtual std::map< G4int, std::vector< G4double > > ReadIXS(G4String fileName, const G4Material *aMaterial)
std::map< const G4Material *, G4LEPTSDiffXS * > theDiffXS
std::map< const G4Material *, G4int > theNXSdat
static constexpr double MeV
#define G4UniformRand()
Definition: Randomize.hh:53
static constexpr double eV
G4GLOB_DLL std::ostream G4cerr
std::vector< G4Material * > G4MaterialTable
static constexpr double cm2
static constexpr double twopi
Definition: SystemOfUnits.h:55
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
int G4int
Definition: G4Types.hh:78
const G4String & GetName() const
Definition: G4VEmModel.hh:777
std::map< const G4Material *, G4double > theIonisPot
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:714
G4double GetA() const
Definition: G4Element.hh:139
void BuildMeanFreePathTable(const G4Material *aMaterial, std::map< G4int, std::vector< G4double > > &integralXS)
#define DBL_MIN
Definition: templates.hh:75
G4double theHighestEnergyLimit
G4GLOB_DLL std::ostream G4cout
std::map< const G4Material *, G4int > theNXSsub
G4ThreeVector SampleNewDirection(const G4Material *aMaterial, G4ThreeVector Dir, G4double e, G4double el)
G4bool ReadParam(G4String fileName, const G4Material *aMaterial)
void clearAndDestroy()
G4double theLowestEnergyLimit
G4int first(char) const
G4double SampleEnergyLoss(const G4Material *aMaterial, G4double eMin, G4double eMax)
#define DBL_MAX
Definition: templates.hh:83
const G4Element * GetElement(G4int iel) const
Definition: G4Material.hh:203
G4int theNumbBinTable
std::map< const G4Material *, G4double > theMolecularMass
void PutValue(size_t index, G4double theValue)
size_t GetNumberOfElements() const
Definition: G4Material.hh:187
G4double GetDensity() const
Definition: G4Material.hh:181