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