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