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