Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4SandiaTable.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 // $Id: G4SandiaTable.cc 96794 2016-05-09 10:09:30Z gcosmo $
27 
28 //
29 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
30 //
31 // 10.06.97 created. V. Grichine
32 // 18.11.98 simplified public interface; new methods for materials. mma
33 // 31.01.01 redesign of ComputeMatSandiaMatrix(). mma
34 // 16.02.01 adapted for STL. mma
35 // 22.02.01 GetsandiaCofForMaterial(energy) return 0 below lowest interval mma
36 // 03.04.01 fnulcof returned if energy < emin
37 // 10.07.01 Migration to STL. M. Verderi.
38 // 03.02.04 Update distructor V.Ivanchenko
39 // 05.03.04 New methods for old sorting algorithm for PAI model. V.Grichine
40 // 26.10.11 new scheme for G4Exception (mma)
41 // 22.05.13 preparation of material table without dynamical arrays. V. Grichine
42 // 09.07.14 modify low limit in GetSandiaCofPerAtom() and material. mma
43 // 10.07.14 modify low limit for water. VI
44 //
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
46 
47 
48 #include "G4SandiaTable.hh"
49 #include "G4StaticSandiaData.hh"
50 #include "G4Material.hh"
51 #include "G4MaterialTable.hh"
52 #include "G4PhysicalConstants.hh"
53 #include "G4SystemOfUnits.hh"
54 
56 {
57  CLHEP::keV,
62 };
63 
65 
66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
67 
69  : fMaterial(material)
70 {
71  fMatSandiaMatrix = nullptr;
72  fMatSandiaMatrixPAI = nullptr;
73  fPhotoAbsorptionCof = nullptr;
74 
76 
77  fMaxInterval = 0;
78  fVerbose = 0;
79 
80  //build the CumulInterval array
81  if(0 == fCumulInterval[0]) {
82  fCumulInterval[0] = 1;
83 
84  for (G4int Z=1; Z<101; ++Z) {
86  }
87  }
88 
89  fMaxInterval = 0;
90  fSandiaCofPerAtom.resize(4,0.0);
91  fLowerI1 = false;
92  //compute macroscopic Sandia coefs for a material
93  ComputeMatSandiaMatrix(); // mma
94 }
95 
96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
97 
98 // Fake default constructor - sets only member data and allocates memory
99 // for usage restricted to object persistency
100 
102  : fMaterial(nullptr),fMatSandiaMatrix(nullptr),
103  fMatSandiaMatrixPAI(nullptr),fPhotoAbsorptionCof(nullptr)
104 {
105  fMaxInterval = 0;
106  fMatNbOfIntervals = 0;
107  fLowerI1 = false;
108  fVerbose = 0;
109  fSandiaCofPerAtom.resize(4,0.0);
110 }
111 
112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
113 
115 {
116  if(fMatSandiaMatrix)
117  {
119  delete fMatSandiaMatrix;
120  }
122  {
124  delete fMatSandiaMatrixPAI;
125  }
127  {
128  delete [] fPhotoAbsorptionCof;
129  }
130 }
131 
132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
133 
134 void
136  std::vector<G4double>& coeff) const
137 {
138 #ifdef G4VERBOSE
139  if(Z < 1 || Z > 100) {
140  Z = PrintErrorZ(Z, "GetSandiaCofPerAtom");
141  }
142  if(4 > coeff.size()) {
143  PrintErrorV("GetSandiaCofPerAtom(): input vector is resized");
144  coeff.resize(4);
145  }
146 #endif
148  //G4double Iopot = fIonizationPotentials[Z]*eV;
149  //if (Emin < Iopot) Emin = Iopot;
150 
151  G4int row = 0;
152  if (energy <= Emin) {
153  energy = Emin;
154 
155  } else {
156  G4int interval = fNbOfIntervals[Z] - 1;
157  row = fCumulInterval[Z-1] + interval;
158  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
159  while ((interval>0) && (energy<fSandiaTable[row][0]*CLHEP::keV)) {
160  --interval;
161  row = fCumulInterval[Z-1] + interval;
162  }
163  }
164 
165  G4double AoverAvo = Z*amu/fZtoAratio[Z];
166 
167  coeff[0]=AoverAvo*funitc[1]*fSandiaTable[row][1];
168  coeff[1]=AoverAvo*funitc[2]*fSandiaTable[row][2];
169  coeff[2]=AoverAvo*funitc[3]*fSandiaTable[row][3];
170  coeff[3]=AoverAvo*funitc[4]*fSandiaTable[row][4];
171 }
172 
173 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
174 
175 void
177  std::vector<G4double>& coeff) const
178 {
179 #ifdef G4VERBOSE
180  if(4 > coeff.size()) {
181  PrintErrorV("GetSandiaCofWater: input vector is resized");
182  coeff.resize(4);
183  }
184 #endif
185  G4int i = 0;
186  if(energy > fH2OlowerI1[0][0]*CLHEP::keV) {
187  i = fH2OlowerInt - 1;
188  for(; i>0; --i) {
189  if(energy >= fH2OlowerI1[i][0]*CLHEP::keV) { break; }
190  }
191  }
192  coeff[0]=funitc[1]*fH2OlowerI1[i][1];
193  coeff[1]=funitc[2]*fH2OlowerI1[i][2];
194  coeff[2]=funitc[3]*fH2OlowerI1[i][3];
195  coeff[3]=funitc[4]*fH2OlowerI1[i][4];
196 }
197 
198 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
199 
201 {
202  return fH2OlowerI1[fH2OlowerInt - 1][0]*CLHEP::keV;
203 }
204 
205 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
206 
208 {
209  return fH2OlowerI1[i][j]*funitc[j];
210 }
211 
212 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
213 
215 {
216 #ifdef G4VERBOSE
217  if(Z < 1 || Z > 100) {
218  Z = PrintErrorZ(Z, "GetSandiaCofPerAtom");
219  }
220 #endif
221  return fZtoAratio[Z];
222 }
223 
224 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
225 
226 #ifdef G4VERBOSE
227 
228 G4int G4SandiaTable::PrintErrorZ(G4int Z, const G4String& ss)
229 {
230  G4String sss = "G4SandiaTable::"+ss+"()";
232  ed << "Atomic number out of range Z= " << Z << "; closest value is used";
233  G4Exception(sss,"mat060",JustWarning,ed,"");
234  return (Z > 100) ? 100 : 1;
235 }
236 
237 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
238 
239 void G4SandiaTable::PrintErrorV(const G4String& ss)
240 {
241  G4String sss = "G4SandiaTable::"+ss;
243  G4Exception(sss,"mat061",JustWarning,"Wrong input parameters");
244 }
245 #endif
246 
247 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
248 
250 {
251  //get list of elements
252  const G4int NbElm = fMaterial->GetNumberOfElements();
253  const G4ElementVector* ElementVector = fMaterial->GetElementVector();
254 
255  G4int* Z = new G4int[NbElm]; //Atomic number
256 
257  //determine the maximum number of energy-intervals for this material
258  G4int MaxIntervals = 0;
259  G4int elm, z;
260 
261  // here we compute only for a mixture, so no waring or exception
262  // if z is out of validity interval
263  for (elm = 0; elm < NbElm; ++elm)
264  {
265  z = G4lrint((*ElementVector)[elm]->GetZ());
266  if(z < 1) { z = 1; }
267  else if(z > 100) { z = 100; }
268  Z[elm] = z;
269  MaxIntervals += fNbOfIntervals[z];
270  }
271 
272  //copy the Energy bins in a tmp1 array
273  //(take care of the Ionization Potential of each element)
274  G4double* tmp1 = new G4double[MaxIntervals];
275  G4double IonizationPot;
276  G4int interval1 = 0;
277 
278  for (elm = 0; elm < NbElm; ++elm)
279  {
280  z = Z[elm];
281  IonizationPot = fIonizationPotentials[z]*CLHEP::eV;
282  for(G4int row = fCumulInterval[z-1]; row<fCumulInterval[z]; ++row)
283  {
284  tmp1[interval1] = std::max(fSandiaTable[row][0]*CLHEP::keV,
285  IonizationPot);
286  ++interval1;
287  }
288  }
289  //sort the energies in strickly increasing values in a tmp2 array
290  //(eliminate redondances)
291 
292  G4double* tmp2 = new G4double[MaxIntervals];
293  G4double Emin;
294  G4int interval2 = 0;
295 
296  do
297  {
298  Emin = DBL_MAX;
299 
300  for ( G4int i1 = 0; i1 < MaxIntervals; ++i1)
301  {
302  Emin = std::min(Emin, tmp1[i1]); //find the minimum
303  }
304  if (Emin < DBL_MAX) {
305  tmp2[interval2] = Emin;
306  ++interval2;
307  }
308  //copy Emin in tmp2
309  for ( G4int j1 = 0; j1 < MaxIntervals; ++j1)
310  {
311  if (tmp1[j1] <= Emin) { tmp1[j1] = DBL_MAX; } //eliminate from tmp1
312  }
313  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
314  } while (Emin < DBL_MAX);
315 
316  //create the sandia matrix for this material
317 
319  G4int interval;
320 
321  for (interval = 0; interval < interval2; ++interval)
322  {
323  fMatSandiaMatrix->push_back( new G4DataVector(5,0.) );
324  }
325 
326  //ready to compute the Sandia coefs for the material
327 
328  const G4double* NbOfAtomsPerVolume = fMaterial->GetVecNbOfAtomsPerVolume();
329 
330  static const G4double prec = 1.e-03*CLHEP::eV;
331  G4double coef, oldsum(0.), newsum(0.);
332  fMatNbOfIntervals = 0;
333 
334  for ( interval = 0; interval < interval2; ++interval)
335  {
336  Emin = (*(*fMatSandiaMatrix)[fMatNbOfIntervals])[0] = tmp2[interval];
337 
338  for ( G4int k = 1; k < 5; ++k ) {
339  (*(*fMatSandiaMatrix)[fMatNbOfIntervals])[k] = 0.;
340  }
341  newsum = 0.;
342 
343  for ( elm = 0; elm < NbElm; elm++ )
344  {
345  GetSandiaCofPerAtom(Z[elm], Emin+prec, fSandiaCofPerAtom);
346 
347  for ( G4int j = 1; j < 5; ++j )
348  {
349  coef = NbOfAtomsPerVolume[elm]*fSandiaCofPerAtom[j-1];
350  (*(*fMatSandiaMatrix)[fMatNbOfIntervals])[j] += coef;
351  newsum += std::abs(coef);
352  }
353  }
354  //check for null or redondant intervals
355 
356  if (newsum != oldsum) { oldsum = newsum; ++fMatNbOfIntervals;}
357  }
358  delete [] Z;
359  delete [] tmp1;
360  delete [] tmp2;
361 
362  if ( fVerbose > 0 )
363  {
364  G4cout<<"G4SandiaTable::ComputeMatSandiaMatrix(), mat = "
365  <<fMaterial->GetName()<<G4endl;
366 
367  for( G4int i = 0; i < fMatNbOfIntervals; ++i)
368  {
369  G4cout<<i<<"\t"<<GetSandiaCofForMaterial(i,0)/keV<<" keV \t"
371  <<"\t"<< GetSandiaCofForMaterial(i,2)
372  <<"\t"<< GetSandiaCofForMaterial(i,3)
373  <<"\t"<< GetSandiaCofForMaterial(i,4)<<G4endl;
374  }
375  }
376 }
377 
379 //
380 // Sandia matrix for PAI models based on vectors ...
381 
383 {
384  G4int MaxIntervals = 0;
385  G4int elm, c, i, j, jj, k, k1, k2, c1, n1, z;
386 
387  const G4int noElm = fMaterial->GetNumberOfElements();
388  const G4ElementVector* ElementVector = fMaterial->GetElementVector();
389 
390  std::vector<G4int> Z(noElm); //Atomic number
391 
392  for ( elm = 0; elm < noElm; elm++ )
393  {
394  z = G4lrint((*ElementVector)[elm]->GetZ());
395  if(z < 1) { z = 1; }
396  else if(z > 100) { z = 100; }
397  Z[elm] = z;
398  MaxIntervals += fNbOfIntervals[Z[elm]];
399  }
400  fMaxInterval = MaxIntervals + 2;
401 
402  if ( fVerbose > 0 )
403  {
404  G4cout<<"G4SandiaTable::ComputeMatSandiaMatrixPAI: fMaxInterval = "
405  <<fMaxInterval<<G4endl;
406  }
407 
408  G4DataVector fPhotoAbsorptionCof0(fMaxInterval);
409  G4DataVector fPhotoAbsorptionCof1(fMaxInterval);
410  G4DataVector fPhotoAbsorptionCof2(fMaxInterval);
411  G4DataVector fPhotoAbsorptionCof3(fMaxInterval);
412  G4DataVector fPhotoAbsorptionCof4(fMaxInterval);
413 
414  for( c = 0; c < fMaxInterval; ++c ) // just in case
415  {
416  fPhotoAbsorptionCof0[c] = 0.;
417  fPhotoAbsorptionCof1[c] = 0.;
418  fPhotoAbsorptionCof2[c] = 0.;
419  fPhotoAbsorptionCof3[c] = 0.;
420  fPhotoAbsorptionCof4[c] = 0.;
421  }
422  c = 1;
423 
424  for(i = 0; i < noElm; ++i)
425  {
426  G4double I1 = fIonizationPotentials[Z[i]]*CLHEP::keV; // I1 in keV
427  n1 = 1;
428 
429  for( j = 1; j < Z[i]; ++j ) n1 += fNbOfIntervals[j];
430 
431  G4int n2 = n1 + fNbOfIntervals[Z[i]];
432 
433  for( k1 = n1; k1 < n2; k1++ )
434  {
435  if( I1 > fSandiaTable[k1][0] )
436  {
437  continue; // no ionization for energies smaller than I1 (first
438  } // ionisation potential)
439  break;
440  }
441  G4int flag = 0;
442 
443  for( c1 = 1; c1 < c; c1++ )
444  {
445  if( fPhotoAbsorptionCof0[c1] == I1 ) // this value already has existed
446  {
447  flag = 1;
448  break;
449  }
450  }
451  if(flag == 0)
452  {
453  fPhotoAbsorptionCof0[c] = I1;
454  ++c;
455  }
456  for( k2 = k1; k2 < n2; k2++ )
457  {
458  flag = 0;
459 
460  for( c1 = 1; c1 < c; c1++ )
461  {
462  if( fPhotoAbsorptionCof0[c1] == fSandiaTable[k2][0] )
463  {
464  flag = 1;
465  break;
466  }
467  }
468  if(flag == 0)
469  {
470  fPhotoAbsorptionCof0[c] = fSandiaTable[k2][0];
471  ++c;
472  }
473  }
474  } // end for(i)
475  // sort out
476 
477  for( i = 1; i < c; ++i )
478  {
479  for( j = i + 1; j < c; ++j )
480  {
481  if( fPhotoAbsorptionCof0[i] > fPhotoAbsorptionCof0[j] )
482  {
483  G4double tmp = fPhotoAbsorptionCof0[i];
484  fPhotoAbsorptionCof0[i] = fPhotoAbsorptionCof0[j];
485  fPhotoAbsorptionCof0[j] = tmp;
486  }
487  }
488  if ( fVerbose > 0)
489  {
490  G4cout<<i<<"\t energy = "<<fPhotoAbsorptionCof0[i]<<G4endl;
491  }
492  }
493  fMaxInterval = c;
494 
495  const G4double* fractionW = fMaterial->GetFractionVector();
496 
497  if ( fVerbose > 0)
498  {
499  for( i = 0; i < noElm; ++i )
500  G4cout<<i<<" = elN, fraction = "<<fractionW[i]<<G4endl;
501  }
502 
503  for( i = 0; i < noElm; ++i )
504  {
505  n1 = 1;
507 
508  for( j = 1; j < Z[i]; ++j ) n1 += fNbOfIntervals[j];
509 
510  G4int n2 = n1 + fNbOfIntervals[Z[i]] - 1;
511 
512  for(k = n1; k < n2; ++k)
513  {
514  G4double B1 = fSandiaTable[k][0];
515  G4double B2 = fSandiaTable[k+1][0];
516 
517  for(G4int q = 1; q < fMaxInterval-1; q++)
518  {
519  G4double E1 = fPhotoAbsorptionCof0[q];
520  G4double E2 = fPhotoAbsorptionCof0[q+1];
521 
522  if ( fVerbose > 0 )
523  {
524  G4cout<<"k = "<<k<<", q = "<<q<<", B1 = "<<B1<<", B2 = "<<B2
525  <<", E1 = "<<E1<<", E2 = "<<E2<<G4endl;
526  }
527  if( B1 > E1 || B2 < E2 || E1 < I1 )
528  {
529  if ( fVerbose > 0 )
530  {
531  G4cout<<"continue for: B1 = "<<B1<<", B2 = "<<B2<<", E1 = "
532  <<E1<<", E2 = "<<E2<<G4endl;
533  }
534  continue;
535  }
536  fPhotoAbsorptionCof1[q] += fSandiaTable[k][1]*fractionW[i];
537  fPhotoAbsorptionCof2[q] += fSandiaTable[k][2]*fractionW[i];
538  fPhotoAbsorptionCof3[q] += fSandiaTable[k][3]*fractionW[i];
539  fPhotoAbsorptionCof4[q] += fSandiaTable[k][4]*fractionW[i];
540  }
541  }
542  // Last interval
543 
544  fPhotoAbsorptionCof1[fMaxInterval-1] += fSandiaTable[k][1]*fractionW[i];
545  fPhotoAbsorptionCof2[fMaxInterval-1] += fSandiaTable[k][2]*fractionW[i];
546  fPhotoAbsorptionCof3[fMaxInterval-1] += fSandiaTable[k][3]*fractionW[i];
547  fPhotoAbsorptionCof4[fMaxInterval-1] += fSandiaTable[k][4]*fractionW[i];
548  } // for(i)
549  c = 0; // Deleting of first intervals where all coefficients = 0
550 
551  do
552  {
553  ++c;
554 
555  if( fPhotoAbsorptionCof1[c] != 0.0 ||
556  fPhotoAbsorptionCof2[c] != 0.0 ||
557  fPhotoAbsorptionCof3[c] != 0.0 ||
558  fPhotoAbsorptionCof4[c] != 0.0 ) continue;
559 
560  if ( fVerbose > 0 )
561  {
562  G4cout<<c<<" = number with zero cofs"<<G4endl;
563  }
564  for( jj = 2; jj < fMaxInterval; ++jj )
565  {
566  fPhotoAbsorptionCof0[jj-1] = fPhotoAbsorptionCof0[jj];
567  fPhotoAbsorptionCof1[jj-1] = fPhotoAbsorptionCof1[jj];
568  fPhotoAbsorptionCof2[jj-1] = fPhotoAbsorptionCof2[jj];
569  fPhotoAbsorptionCof3[jj-1] = fPhotoAbsorptionCof3[jj];
570  fPhotoAbsorptionCof4[jj-1] = fPhotoAbsorptionCof4[jj];
571  }
572  --fMaxInterval;
573  --c;
574  }
575  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
576  while( c < fMaxInterval - 1 );
577 
578  if( fPhotoAbsorptionCof0[fMaxInterval-1] == 0.0 ) fMaxInterval--;
579 
580  // create the sandia matrix for this material
581 
583 
584  G4double density = fMaterial->GetDensity();
585 
586  for (i = 0; i < fMaxInterval; ++i) // -> G4units
587  {
588  fPhotoAbsorptionCof0[i+1] *= funitc[0];
589  fPhotoAbsorptionCof1[i+1] *= funitc[1]*density;
590  fPhotoAbsorptionCof2[i+1] *= funitc[2]*density;
591  fPhotoAbsorptionCof3[i+1] *= funitc[3]*density;
592  fPhotoAbsorptionCof4[i+1] *= funitc[4]*density;
593  }
594  if(fLowerI1)
595  {
596  if( fMaterial->GetName() == "G4_WATER")
597  {
598  fMaxInterval += fH2OlowerInt;
599 
600  for (i = 0; i < fMaxInterval; ++i) // init vector table
601  {
602  fMatSandiaMatrixPAI->push_back( new G4DataVector(5,0.) );
603  }
604  for (i = 0; i < fH2OlowerInt; ++i)
605  {
606  (*(*fMatSandiaMatrixPAI)[i])[0] = fH2OlowerI1[i][0];
607  (*(*fMatSandiaMatrixPAI)[i])[1] = fH2OlowerI1[i][1];
608  (*(*fMatSandiaMatrixPAI)[i])[2] = fH2OlowerI1[i][2];
609  (*(*fMatSandiaMatrixPAI)[i])[3] = fH2OlowerI1[i][3];
610  (*(*fMatSandiaMatrixPAI)[i])[4] = fH2OlowerI1[i][4];
611  }
612  for (i = fH2OlowerInt; i < fMaxInterval; ++i)
613  {
614  (*(*fMatSandiaMatrixPAI)[i])[0] = fPhotoAbsorptionCof0[i+1-fH2OlowerInt];
615  (*(*fMatSandiaMatrixPAI)[i])[1] = fPhotoAbsorptionCof1[i+1-fH2OlowerInt];
616  (*(*fMatSandiaMatrixPAI)[i])[2] = fPhotoAbsorptionCof2[i+1-fH2OlowerInt];
617  (*(*fMatSandiaMatrixPAI)[i])[3] = fPhotoAbsorptionCof3[i+1-fH2OlowerInt];
618  (*(*fMatSandiaMatrixPAI)[i])[4] = fPhotoAbsorptionCof4[i+1-fH2OlowerInt];
619  }
620  }
621  }
622  else
623  {
624  for (i = 0; i < fMaxInterval; ++i) // init vector table
625  {
626  fMatSandiaMatrixPAI->push_back( new G4DataVector(5,0.) );
627  }
628  for (i = 0; i < fMaxInterval; ++i)
629  {
630  (*(*fMatSandiaMatrixPAI)[i])[0] = fPhotoAbsorptionCof0[i+1];
631  (*(*fMatSandiaMatrixPAI)[i])[1] = fPhotoAbsorptionCof1[i+1]; // *density;
632  (*(*fMatSandiaMatrixPAI)[i])[2] = fPhotoAbsorptionCof2[i+1]; // *density;
633  (*(*fMatSandiaMatrixPAI)[i])[3] = fPhotoAbsorptionCof3[i+1]; // *density;
634  (*(*fMatSandiaMatrixPAI)[i])[4] = fPhotoAbsorptionCof4[i+1]; // *density;
635  }
636  }
637  // --fMaxInterval;
638  // to avoid duplicate at 500 keV or extra zeros in last interval
639 
640  if ( fVerbose > 0 )
641  {
642  G4cout<<"G4SandiaTable::ComputeMatSandiaMatrixPAI(), mat = "
643  <<fMaterial->GetName()<<G4endl;
644 
645  for( i = 0; i < fMaxInterval; ++i)
646  {
647  G4cout<<i<<"\t"<<GetSandiaMatTablePAI(i,0)/keV<<" keV \t"
648  << GetSandiaMatTablePAI(i,1)
649  <<"\t"<<GetSandiaMatTablePAI(i,2)
650  <<"\t"<<GetSandiaMatTablePAI(i,3)
651  <<"\t"<<GetSandiaMatTablePAI(i,4)<<G4endl;
652  }
653  }
654  return;
655 }
656 
658 // Methods for PAI model only
659 //
660 
662 {
663  fMaterial = 0;
664  fMatNbOfIntervals = 0;
665  fMatSandiaMatrix = 0;
668 
669  fMaxInterval = 0;
670  fVerbose = 0;
671  fLowerI1 = false;
672 
673  fSandiaCofPerAtom.resize(4,0.0);
674 
675  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
676  G4int numberOfMat = G4Material::GetNumberOfMaterials();
677 
678  if ( matIndex >= 0 && matIndex < numberOfMat)
679  {
680  fMaterial = (*theMaterialTable)[matIndex];
681  }
682  else
683  {
684  G4Exception("G4SandiaTable::G4SandiaTable(G4int matIndex)", "mat401",
685  FatalException, "wrong matIndex");
686  }
687 }
688 
690 
692 {
693  fMaterial = 0;
694  fMatNbOfIntervals = 0;
695  fMatSandiaMatrix = 0;
698 
699  fMaxInterval = 0;
700  fVerbose = 0;
701  fLowerI1 = false;
702 
703  fSandiaCofPerAtom.resize(4,0.0);
704 }
705 
707 
709 {
710  fMaterial = mat;
712 }
713 
715 
717 {
718  return fMaxInterval;
719 }
720 
722 
724 {
726  return fPhotoAbsorptionCof;
727 }
728 
730 
732  G4int i,
733  G4int j )
734 {
735  G4double tmp = da[i][0] ;
736  da[i][0] = da[j][0] ;
737  da[j][0] = tmp ;
738 }
739 
741 
743 {
744  return fPhotoAbsorptionCof[i][j]*funitc[j];
745 }
746 
748 //
749 // Bubble sorting of left energy interval in SandiaTable in ascening order
750 //
751 
752 void
754 {
755  for(G4int i = 1;i < sz; ++i )
756  {
757  for(G4int j = i + 1;j < sz; ++j )
758  {
759  if(da[i][0] > da[j][0]) SandiaSwap(da,i,j);
760  }
761  }
762 }
763 
765 //
766 // SandiaIntervals
767 //
768 
770 {
771  G4int c, i, flag = 0, n1 = 1;
772  G4int j, c1, k1, k2;
773  G4double I1;
774  fMaxInterval = 0;
775 
776  for( i = 0; i < el; ++i ) fMaxInterval += fNbOfIntervals[ Z[i] ];
777 
778  fMaxInterval += 2;
779 
780  if( fVerbose > 0 ) {
781  G4cout<<"begin sanInt, fMaxInterval = "<<fMaxInterval<<G4endl;
782  }
783 
785 
786  for( i = 0; i < fMaxInterval; ++i ) {
787  fPhotoAbsorptionCof[i] = new G4double[5];
788  }
789  // for(c = 0; c < fIntervalLimit; ++c) // just in case
790 
791  for( c = 0; c < fMaxInterval; ++c ) { fPhotoAbsorptionCof[c][0] = 0.; }
792 
793  c = 1;
794 
795  for( i = 0; i < el; ++i )
796  {
797  I1 = fIonizationPotentials[ Z[i] ]*keV; // First ionization
798  n1 = 1; // potential in keV
799 
800  for( j = 1; j < Z[i]; ++j ) n1 += fNbOfIntervals[j];
801 
802  G4int n2 = n1 + fNbOfIntervals[Z[i]];
803 
804  for( k1 = n1; k1 < n2; k1++ )
805  {
806  if( I1 > fSandiaTable[k1][0] )
807  {
808  continue; // no ionization for energies smaller than I1 (first
809  } // ionisation potential)
810  break;
811  }
812  flag = 0;
813 
814  for( c1 = 1; c1 < c; c1++ )
815  {
816  if( fPhotoAbsorptionCof[c1][0] == I1 ) // this value already has existed
817  {
818  flag = 1;
819  break;
820  }
821  }
822  if( flag == 0 )
823  {
824  fPhotoAbsorptionCof[c][0] = I1;
825  ++c;
826  }
827  for( k2 = k1; k2 < n2; k2++ )
828  {
829  flag = 0;
830 
831  for( c1 = 1; c1 < c; c1++ )
832  {
833  if( fPhotoAbsorptionCof[c1][0] == fSandiaTable[k2][0] )
834  {
835  flag = 1;
836  break;
837  }
838  }
839  if( flag == 0 )
840  {
841  fPhotoAbsorptionCof[c][0] = fSandiaTable[k2][0];
842  if( fVerbose > 0 ) {
843  G4cout<<"sanInt, c = "<<c<<", E_c = "<<fPhotoAbsorptionCof[c][0]
844  <<G4endl;
845  }
846  ++c;
847  }
848  }
849  } // end for(i)
850 
852  fMaxInterval = c;
853  if( fVerbose > 0 ) {
854  G4cout<<"end SanInt, fMaxInterval = "<<fMaxInterval<<G4endl;
855  }
856  return c;
857 }
858 
860 //
861 // SandiaMixing
862 //
863 
864 G4int
866  const G4double fractionW[],
867  G4int el,
868  G4int mi )
869 {
870  G4int i, j, n1, k, c=1, jj, kk;
871  G4double I1, B1, B2, E1, E2;
872 
873  for( i = 0; i < mi; ++i )
874  {
875  for( j = 1; j < 5; ++j ) fPhotoAbsorptionCof[i][j] = 0.;
876  }
877  for( i = 0; i < el; ++i )
878  {
879  n1 = 1;
880  I1 = fIonizationPotentials[Z[i]]*keV;
881 
882  for( j = 1; j < Z[i]; ++j ) n1 += fNbOfIntervals[j];
883 
884  G4int n2 = n1 + fNbOfIntervals[Z[i]] - 1;
885 
886  for( k = n1; k < n2; ++k )
887  {
888  B1 = fSandiaTable[k][0];
889  B2 = fSandiaTable[k+1][0];
890 
891  for( c = 1; c < mi-1; ++c )
892  {
893  E1 = fPhotoAbsorptionCof[c][0];
894  E2 = fPhotoAbsorptionCof[c+1][0];
895 
896  if( B1 > E1 || B2 < E2 || E1 < I1 ) continue;
897 
898  for( j = 1; j < 5; ++j )
899  {
900  fPhotoAbsorptionCof[c][j] += fSandiaTable[k][j]*fractionW[i];
901  if( fVerbose > 0 )
902  {
903  G4cout<<"c="<<c<<"; j="<<j<<"; fST="<<fSandiaTable[k][j]
904  <<"; frW="<<fractionW[i]<<G4endl;
905  }
906  }
907  }
908  }
909  for( j = 1; j < 5; ++j ) // Last interval
910  {
911  fPhotoAbsorptionCof[mi-1][j] += fSandiaTable[k][j]*fractionW[i];
912  if( fVerbose > 0 )
913  {
914  G4cout<<"mi-1="<<mi-1<<"; j="<<j<<"; fST="<<fSandiaTable[k][j]
915  <<"; frW="<<fractionW[i]<<G4endl;
916  }
917  }
918  } // for(i)
919  c = 0; // Deleting of first intervals where all coefficients = 0
920 
921  do
922  {
923  ++c;
924 
925  if( fPhotoAbsorptionCof[c][1] != 0.0 ||
926  fPhotoAbsorptionCof[c][2] != 0.0 ||
927  fPhotoAbsorptionCof[c][3] != 0.0 ||
928  fPhotoAbsorptionCof[c][4] != 0.0 ) continue;
929 
930  for( jj = 2; jj < mi; ++jj )
931  {
932  for( kk = 0; kk < 5; ++kk ) {
933  fPhotoAbsorptionCof[jj-1][kk] = fPhotoAbsorptionCof[jj][kk];
934  }
935  }
936  mi--;
937  c--;
938  }
939  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
940  while( c < mi - 1 );
941 
942  if( fVerbose > 0 ) G4cout<<"end SanMix, mi = "<<mi<<G4endl;
943 
944  return mi;
945 }
946 
948 
950 {
951  return fMatNbOfIntervals;
952 }
953 
955 
956 G4double
958 {
959 #ifdef G4VERBOSE
960  if(Z < 1 || Z > 100) {
961  Z = PrintErrorZ(Z, "GetSandiaPerAtom");
962  }
963  if(interval<0 || interval>=fNbOfIntervals[Z]) {
964  PrintErrorV("GetSandiaPerAtom");
965  interval = (interval<0) ? 0 : fNbOfIntervals[Z]-1;
966  }
967  if(j<0 || j>4) {
968  PrintErrorV("GetSandiaPerAtom");
969  j = (j<0) ? 0 : 4;
970  }
971 #endif
972  G4int row = fCumulInterval[Z-1] + interval;
973  G4double x = fSandiaTable[row][0]*CLHEP::keV;
974  if (j > 0) {
975  x = Z*CLHEP::amu/fZtoAratio[Z]*fSandiaTable[row][j]*funitc[j];
976  }
977  return x;
978 }
979 
981 
982 G4double
984 {
985 #ifdef G4VERBOSE
986  if(interval<0 || interval>=fMatNbOfIntervals) {
987  PrintErrorV("GetSandiaCofForMaterial");
988  interval = (interval<0) ? 0 : fMatNbOfIntervals-1;
989  }
990  if(j<0 || j>4) {
991  PrintErrorV("GetSandiaCofForMaterial");
992  j = (j<0) ? 0 : 4;
993  }
994 #endif
995  return ((*(*fMatSandiaMatrix)[interval])[j]);
996 }
997 
999 
1000 const G4double*
1002 {
1003  G4int interval = 0;
1004  if (energy > (*(*fMatSandiaMatrix)[0])[0]) {
1005  interval = fMatNbOfIntervals - 1;
1006  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
1007  while ((interval>0)&&(energy<(*(*fMatSandiaMatrix)[interval])[0]))
1008  { --interval; }
1009  }
1010  return &((*(*fMatSandiaMatrix)[interval])[1]);
1011 }
1012 
1014 
1015 G4double
1017 {
1018 #ifdef G4VERBOSE
1019  if(interval<0 || interval>=fMatNbOfIntervals) {
1020  PrintErrorV("GetSandiaCofForMaterial");
1021  interval = (interval<0) ? 0 : fMatNbOfIntervals-1;
1022  }
1023  if(j<0 || j>4) {
1024  PrintErrorV("GetSandiaCofForMaterial");
1025  j = (j<0) ? 0 : 4;
1026  }
1027 #endif
1028  return ((*(*fMatSandiaMatrix)[interval])[j])*funitc[j];
1029 }
1030 
1032 
1033 G4double
1035 {
1036 #ifdef G4VERBOSE
1037  if(interval<0 || interval>=fMaxInterval) {
1038  PrintErrorV("GetSandiaCofForMaterialPAI");
1039  interval = (interval<0) ? 0 : fMaxInterval-1;
1040  }
1041  if(j<0 || j>4) {
1042  PrintErrorV("GetSandiaCofForMaterialPAI");
1043  j = (j<0) ? 0 : 4;
1044  }
1045 #endif
1046  return ((*(*fMatSandiaMatrixPAI)[interval])[j]);
1047 }
1048 
1050 //
1051 // Sandia interval and mixing calculations for materialCutsCouple constructor
1052 //
1053 
1055 {
1056  G4int MaxIntervals = 0;
1057  G4int elm, c, i, j, jj, k, kk, k1, k2, c1, n1;
1058 
1059  const G4int noElm = fMaterial->GetNumberOfElements();
1060  const G4ElementVector* ElementVector = fMaterial->GetElementVector();
1061  G4int* Z = new G4int[noElm]; //Atomic number
1062 
1063  for (elm = 0; elm<noElm; ++elm)
1064  {
1065  Z[elm] = (G4int)(*ElementVector)[elm]->GetZ();
1066  MaxIntervals += fNbOfIntervals[Z[elm]];
1067  }
1068  fMaxInterval = 0;
1069 
1070  for(i = 0; i < noElm; ++i) fMaxInterval += fNbOfIntervals[Z[i]];
1071 
1072  fMaxInterval += 2;
1073 
1074  // G4cout<<"fMaxInterval = "<<fMaxInterval<<G4endl;
1075 
1077 
1078  for(i = 0; i < fMaxInterval; ++i)
1079  {
1080  fPhotoAbsorptionCof[i] = new G4double[5];
1081  }
1082 
1083  // for(c = 0; c < fIntervalLimit; ++c) // just in case
1084 
1085  for(c = 0; c < fMaxInterval; ++c) // just in case
1086  {
1087  fPhotoAbsorptionCof[c][0] = 0.;
1088  }
1089  c = 1;
1090 
1091  for(i = 0; i < noElm; ++i)
1092  {
1093  G4double I1 = fIonizationPotentials[Z[i]]*keV; // First ionization
1094  n1 = 1; // potential in keV
1095 
1096  for(j = 1; j < Z[i]; ++j)
1097  {
1098  n1 += fNbOfIntervals[j];
1099  }
1100  G4int n2 = n1 + fNbOfIntervals[Z[i]];
1101 
1102  for(k1 = n1; k1 < n2; ++k1)
1103  {
1104  if(I1 > fSandiaTable[k1][0])
1105  {
1106  continue; // no ionization for energies smaller than I1 (first
1107  } // ionisation potential)
1108  break;
1109  }
1110  G4int flag = 0;
1111 
1112  for(c1 = 1; c1 < c; ++c1)
1113  {
1114  if(fPhotoAbsorptionCof[c1][0] == I1) // this value already has existed
1115  {
1116  flag = 1;
1117  break;
1118  }
1119  }
1120  if(flag == 0)
1121  {
1122  fPhotoAbsorptionCof[c][0] = I1;
1123  ++c;
1124  }
1125  for(k2 = k1; k2 < n2; ++k2)
1126  {
1127  flag = 0;
1128 
1129  for(c1 = 1; c1 < c; ++c1)
1130  {
1131  if(fPhotoAbsorptionCof[c1][0] == fSandiaTable[k2][0])
1132  {
1133  flag = 1;
1134  break;
1135  }
1136  }
1137  if(flag == 0)
1138  {
1139  fPhotoAbsorptionCof[c][0] = fSandiaTable[k2][0];
1140  ++c;
1141  }
1142  }
1143  } // end for(i)
1144 
1146  fMaxInterval = c;
1147 
1148  const G4double* fractionW = fMaterial->GetFractionVector();
1149 
1150  for(i = 0; i < fMaxInterval; ++i)
1151  {
1152  for(j = 1; j < 5; ++j) fPhotoAbsorptionCof[i][j] = 0.;
1153  }
1154  for(i = 0; i < noElm; ++i)
1155  {
1156  n1 = 1;
1157  G4double I1 = fIonizationPotentials[Z[i]]*keV;
1158 
1159  for(j = 1; j < Z[i]; ++j)
1160  {
1161  n1 += fNbOfIntervals[j];
1162  }
1163  G4int n2 = n1 + fNbOfIntervals[Z[i]] - 1;
1164 
1165  for(k = n1; k < n2; ++k)
1166  {
1167  G4double B1 = fSandiaTable[k][0];
1168  G4double B2 = fSandiaTable[k+1][0];
1169  for(G4int q = 1; q < fMaxInterval-1; q++)
1170  {
1171  G4double E1 = fPhotoAbsorptionCof[q][0];
1172  G4double E2 = fPhotoAbsorptionCof[q+1][0];
1173  if(B1 > E1 || B2 < E2 || E1 < I1)
1174  {
1175  continue;
1176  }
1177  for(j = 1; j < 5; ++j)
1178  {
1179  fPhotoAbsorptionCof[q][j] += fSandiaTable[k][j]*fractionW[i];
1180  }
1181  }
1182  }
1183  for(j = 1; j < 5; ++j) // Last interval
1184  {
1185  fPhotoAbsorptionCof[fMaxInterval-1][j] +=
1186  fSandiaTable[k][j]*fractionW[i];
1187  }
1188  } // for(i)
1189 
1190  c = 0; // Deleting of first intervals where all coefficients = 0
1191 
1192  do
1193  {
1194  ++c;
1195 
1196  if( fPhotoAbsorptionCof[c][1] != 0.0 ||
1197  fPhotoAbsorptionCof[c][2] != 0.0 ||
1198  fPhotoAbsorptionCof[c][3] != 0.0 ||
1199  fPhotoAbsorptionCof[c][4] != 0.0 ) continue;
1200 
1201  for(jj = 2; jj < fMaxInterval; ++jj)
1202  {
1203  for(kk = 0; kk < 5; ++kk)
1204  {
1205  fPhotoAbsorptionCof[jj-1][kk]= fPhotoAbsorptionCof[jj][kk];
1206  }
1207  }
1208  --fMaxInterval;
1209  --c;
1210  }
1211  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
1212  while( c < fMaxInterval - 1 );
1213 
1214  // create the sandia matrix for this material
1215 
1216  --fMaxInterval; // vmg 20.11.10
1217 
1219 
1220  for (i = 0; i < fMaxInterval; ++i)
1221  {
1222  fMatSandiaMatrix->push_back(new G4DataVector(5,0.));
1223  }
1224  for ( i = 0; i < fMaxInterval; ++i )
1225  {
1226  for( j = 0; j < 5; ++j )
1227  {
1228  (*(*fMatSandiaMatrix)[i])[j] = fPhotoAbsorptionCof[i+1][j];
1229  }
1230  }
1232 
1233  if ( fVerbose > 0 )
1234  {
1235  G4cout<<"vmg, G4SandiaTable::ComputeMatTable(), mat = "
1236  <<fMaterial->GetName()<<G4endl;
1237 
1238  for ( i = 0; i < fMaxInterval; ++i )
1239  {
1240  // G4cout<<i<<"\t"<<(*(*fMatSandiaMatrix)[i])[0]<<" keV \t"
1241  // <<(*(*fMatSandiaMatrix)[i])[1]
1242  // <<"\t"<<(*(*fMatSandiaMatrix)[i])[2]<<"\t"
1243  // <<(*(*fMatSandiaMatrix)[i])[3]
1244  // <<"\t"<<(*(*fMatSandiaMatrix)[i])[4]<<G4endl;
1245 
1246  G4cout<<i<<"\t"<<GetSandiaCofForMaterial(i,0)/keV
1247  <<" keV \t"<<this->GetSandiaCofForMaterial(i,1)
1248  <<"\t"<<this->GetSandiaCofForMaterial(i,2)
1249  <<"\t"<<this->GetSandiaCofForMaterial(i,3)
1250  <<"\t"<<this->GetSandiaCofForMaterial(i,4)<<G4endl;
1251  }
1252  }
1253  delete [] Z;
1254  return;
1255 }
1256 
1257 //
1258 //
G4int SandiaMixing(G4int Z[], const G4double *fractionW, G4int el, G4int mi)
Float_t x
Definition: compare.C:6
void clearAndDestroy()
T max(const T t1, const T t2)
brief Return the largest of the two arguments
static const double prec
Definition: RanecuEngine.cc:58
static const G4double Emin
const G4double * GetFractionVector() const
Definition: G4Material.hh:195
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:45
G4double ** GetPointerToCof()
static size_t GetNumberOfMaterials()
Definition: G4Material.cc:600
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:593
static constexpr double keV
Definition: G4SIunits.hh:216
static constexpr double g
#define G4endl
Definition: G4ios.hh:61
Double_t z
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:191
static const G4double fIonizationPotentials[101]
std::vector< G4double > fSandiaCofPerAtom
void SandiaSwap(G4double **da, G4int i, G4int j)
void ComputeMatSandiaMatrix()
const G4double * GetVecNbOfAtomsPerVolume() const
Definition: G4Material.hh:207
Float_t tmp
G4double GetSandiaPerAtom(G4int Z, G4int, G4int) const
static const G4int fNbOfIntervals[101]
G4int GetMaxInterval() const
void ComputeMatTable()
void SandiaSort(G4double **da, G4int sz)
G4int GetMatNbOfIntervals() const
const G4String & GetName() const
Definition: G4Material.hh:179
Float_t Z
static const char sss[MAX_N_PAR+2]
Definition: Evaluator.cc:64
void GetSandiaCofWater(G4double energy, std::vector< G4double > &coeff) const
static G4double GetZtoA(G4int Z)
double G4double
Definition: G4Types.hh:76
void Initialize(G4Material *)
double energy
Definition: plottest35.C:25
G4double ** fPhotoAbsorptionCof
void GetSandiaCofPerAtom(G4int Z, G4double energy, std::vector< G4double > &coeff) const
G4int SandiaIntervals(G4int Z[], G4int el)
G4Material * fMaterial
static constexpr double eV
Float_t mat
static const G4int fH2OlowerInt
G4double GetSandiaMatTablePAI(G4int, G4int) const
G4double GetSandiaCofForMaterial(G4int, G4int) const
G4int fMatNbOfIntervals
static const G4double fZtoAratio[101]
std::vector< G4Material * > G4MaterialTable
static constexpr double cm2
G4double GetPhotoAbsorpCof(G4int i, G4int j) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
G4double GetWaterEnergyLimit() const
void ComputeMatSandiaMatrixPAI()
int G4lrint(double ad)
Definition: templates.hh:151
std::vector< G4Element * > G4ElementVector
int G4int
Definition: G4Types.hh:78
G4OrderedTable * fMatSandiaMatrixPAI
static const G4double fSandiaTable[981][5]
G4GLOB_DLL std::ostream G4cout
static G4int fCumulInterval[101]
static const G4double fH2OlowerI1[23][5]
G4OrderedTable * fMatSandiaMatrix
static constexpr double amu
G4double GetSandiaMatTable(G4int, G4int) const
static const G4double funitc[5]
static constexpr double keV
#define DBL_MAX
Definition: templates.hh:83
G4double GetWaterCofForMaterial(G4int, G4int) const
size_t GetNumberOfElements() const
Definition: G4Material.hh:187
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double GetDensity() const
Definition: G4Material.hh:181