source: trunk/source/materials/src/G4SandiaTable.cc @ 1315

Last change on this file since 1315 was 1196, checked in by garnier, 15 years ago

update CVS release candidate geant4.9.3.01
File size: 23.0 KB
Line 
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//
27// $Id: G4SandiaTable.cc,v 1.34 2007/10/02 10:13:33 vnivanch Exp $
28// GEANT4 tag $Name: materials-V09-02-18 $
29//
30//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
31//
32// 10.06.97 created. V. Grichine
33// 18.11.98 simplified public interface; new methods for materials.  mma
34// 31.01.01 redesign of ComputeMatSandiaMatrix().  mma
35// 16.02.01 adapted for STL.  mma
36// 22.02.01 GetsandiaCofForMaterial(energy) return 0 below lowest interval  mma 
37// 03.04.01 fnulcof returned if energy < emin
38// 10.07.01 Migration to STL. M. Verderi.
39// 03.02.04 Update distructor V.Ivanchenko
40// 05.03.04 New methods for old sorting algorithm for PAI model. V.Grichine
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
50G4int    G4SandiaTable::fCumulInterval[101];
51G4double G4SandiaTable::fSandiaCofPerAtom[4];
52G4double const G4SandiaTable::funitc[4] = {cm2*keV/g,     
53                                           cm2*keV*keV/g,     
54                                           cm2*keV*keV*keV/g,     
55                                           cm2*keV*keV*keV*keV/g};
56 
57//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
58
59G4SandiaTable::G4SandiaTable(G4Material* material)
60  : fMaterial(material)
61{
62  fMatSandiaMatrix    = 0; 
63  fMatSandiaMatrixPAI = 0;
64  fPhotoAbsorptionCof = 0;
65
66  //build the CumulInterval array
67  fCumulInterval[0] = 1;
68  for (G4int Z=1; Z<101; Z++) {
69    fCumulInterval[Z] = fCumulInterval[Z-1] + fNbOfIntervals[Z];
70  }
71  //compute macroscopic Sandia coefs for a material   
72  ComputeMatSandiaMatrix();
73 
74  //initialisation of fnulcof
75  fnulcof[0] = fnulcof[1] = fnulcof[2] = fnulcof[3] = 0.;   
76}
77                                                       
78//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
79
80// Fake default constructor - sets only member data and allocates memory
81//                            for usage restricted to object persistency
82
83G4SandiaTable::G4SandiaTable(__void__&)
84  : fMaterial(0), fMatSandiaMatrix(0), fPhotoAbsorptionCof(0)
85{
86}
87
88//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
89
90G4SandiaTable::~G4SandiaTable()
91{ 
92  if(fMatSandiaMatrix) {
93    fMatSandiaMatrix->clearAndDestroy();
94    delete fMatSandiaMatrix;
95  }
96  if(fMatSandiaMatrixPAI) {
97    fMatSandiaMatrixPAI->clearAndDestroy();
98    delete fMatSandiaMatrixPAI;
99  }
100  if(fPhotoAbsorptionCof)
101  {
102    delete [] fPhotoAbsorptionCof ;
103  }
104}
105
106//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
107
108G4double G4SandiaTable::GetZtoA(G4int Z)
109{
110  return fZtoAratio[Z];
111}
112
113//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
114
115G4double*
116G4SandiaTable::GetSandiaCofPerAtom(G4int Z, G4double energy)
117{
118  assert (Z > 0 && Z < 101);
119   
120  G4double Emin  = fSandiaTable[fCumulInterval[Z-1]][0]*keV;
121  G4double Iopot = fIonizationPotentials[Z]*eV;
122  if (Iopot > Emin) Emin = Iopot;
123   
124  G4int interval = fNbOfIntervals[Z] - 1;
125  G4int row = fCumulInterval[Z-1] + interval;
126  while ((interval>0) && (energy<fSandiaTable[row][0]*keV)) {
127    --interval;
128    row = fCumulInterval[Z-1] + interval;
129  }
130  if (energy >= Emin)
131    {       
132      G4double AoverAvo = Z*amu/fZtoAratio[Z];
133         
134      fSandiaCofPerAtom[0]=AoverAvo*funitc[0]*fSandiaTable[row][1];     
135      fSandiaCofPerAtom[1]=AoverAvo*funitc[1]*fSandiaTable[row][2];     
136      fSandiaCofPerAtom[2]=AoverAvo*funitc[2]*fSandiaTable[row][3];     
137      fSandiaCofPerAtom[3]=AoverAvo*funitc[3]*fSandiaTable[row][4];
138    }
139  else 
140    {
141      fSandiaCofPerAtom[0] = fSandiaCofPerAtom[1] = fSandiaCofPerAtom[2] =
142        fSandiaCofPerAtom[3] = 0.;
143    }               
144  return fSandiaCofPerAtom;     
145}
146                                                       
147//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
148
149void G4SandiaTable::ComputeMatSandiaMatrix()
150{ 
151  //get list of elements
152  const G4int NbElm = fMaterial->GetNumberOfElements();
153  const G4ElementVector* ElementVector = fMaterial->GetElementVector();
154 
155  G4int* Z = new G4int[NbElm];               //Atomic number
156
157  //   
158  //determine the maximum number of energy-intervals for this material
159  //
160  G4int MaxIntervals = 0;
161  G4int elm;   
162  for (elm=0; elm<NbElm; elm++) {
163    Z[elm] = (G4int)(*ElementVector)[elm]->GetZ();
164    MaxIntervals += fNbOfIntervals[Z[elm]];
165  } 
166     
167  //
168  //copy the Energy bins in a tmp1 array
169  //(take care of the Ionization Potential of each element)
170  //
171  G4double* tmp1 = new G4double[MaxIntervals]; 
172  G4double IonizationPot;
173  G4int interval1=0;
174  for (elm=0; elm<NbElm; elm++) {
175
176    IonizationPot = GetIonizationPot(Z[elm]);
177    for (G4int row=fCumulInterval[Z[elm]-1];row<fCumulInterval[Z[elm]];row++) {
178      tmp1[interval1++] = std::max(fSandiaTable[row][0]*keV,IonizationPot);
179    }
180  }   
181  //       
182  //sort the energies in strickly increasing values in a tmp2 array
183  //(eliminate redondances)
184  //
185  G4double* tmp2 = new G4double[MaxIntervals];
186  G4double Emin;
187  G4int interval2 = 0;
188 
189  do {
190    Emin = DBL_MAX;
191    for (G4int i1=0; i1<MaxIntervals; i1++) {
192      if (tmp1[i1] < Emin) Emin = tmp1[i1];          //find the minimum
193    }
194    if (Emin < DBL_MAX) tmp2[interval2++] = Emin;    //copy Emin in tmp2
195    for (G4int j1=0; j1<MaxIntervals; j1++) {
196      if (tmp1[j1] <= Emin) tmp1[j1] = DBL_MAX;      //eliminate from tmp1         
197    }
198  } while (Emin < DBL_MAX);
199                                   
200  //   
201  //create the sandia matrix for this material
202  // 
203  fMatSandiaMatrix = new G4OrderedTable();
204  G4int interval;
205  for (interval=0; interval<interval2; interval++) {
206    fMatSandiaMatrix->push_back(new G4DataVector(5,0.));
207  }
208  //
209  //ready to compute the Sandia coefs for the material
210  //
211  const G4double* NbOfAtomsPerVolume = fMaterial->GetVecNbOfAtomsPerVolume();
212 
213  const G4double prec = 1.e-03*eV;
214  G4double coef, oldsum(0.), newsum(0.);
215  fMatNbOfIntervals = 0;
216         
217  for (interval=0; interval<interval2; interval++) {
218
219    Emin = (*(*fMatSandiaMatrix)[fMatNbOfIntervals])[0] = tmp2[interval];
220    for (G4int k=1; k<5; k++) {
221      (*(*fMatSandiaMatrix)[fMatNbOfIntervals])[k]=0.;     
222    }
223    newsum = 0.;
224     
225    for (elm=0; elm<NbElm; elm++) {
226       
227      GetSandiaCofPerAtom(Z[elm], Emin+prec);
228      for (G4int j=1; j<5; j++) {
229
230        coef = NbOfAtomsPerVolume[elm]*fSandiaCofPerAtom[j-1];
231        (*(*fMatSandiaMatrix)[fMatNbOfIntervals])[j] += coef;
232        newsum += std::abs(coef);
233      }                                               
234    }         
235                                                         
236    //check for null or redondant intervals     
237    if (newsum != oldsum) { oldsum = newsum; fMatNbOfIntervals++;}
238  }
239  delete [] Z;
240  delete [] tmp1;
241  delete [] tmp2;
242}
243
244//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
245
246void G4SandiaTable::ComputeMatSandiaMatrixPAI()
247{ 
248  G4int MaxIntervals = 0;
249  G4int elm, c, i, j, jj, k, k1, k2, c1, n1;   
250
251  const G4int noElm = fMaterial->GetNumberOfElements();
252  const G4ElementVector* ElementVector = fMaterial->GetElementVector(); 
253  G4int* Z = new G4int[noElm];               //Atomic number
254
255  for (elm = 0; elm<noElm; elm++)
256    { 
257      Z[elm] = (G4int)(*ElementVector)[elm]->GetZ();
258      MaxIntervals += fNbOfIntervals[Z[elm]];
259    } 
260  fMaxInterval = MaxIntervals + 2;
261
262  //  G4cout<<"fMaxInterval = "<<fMaxInterval<<G4endl ;
263 
264  G4double* fPhotoAbsorptionCof0 = new G4double[fMaxInterval];
265  G4double* fPhotoAbsorptionCof1 = new G4double[fMaxInterval];
266  G4double* fPhotoAbsorptionCof2 = new G4double[fMaxInterval];
267  G4double* fPhotoAbsorptionCof3 = new G4double[fMaxInterval];
268  G4double* fPhotoAbsorptionCof4 = new G4double[fMaxInterval];
269
270  for(c = 0 ; c<fMaxInterval ; c++)   // just in case
271  {
272    fPhotoAbsorptionCof0[c] = 0. ;
273    fPhotoAbsorptionCof1[c] = 0. ;
274    fPhotoAbsorptionCof2[c] = 0. ;
275    fPhotoAbsorptionCof3[c] = 0. ;
276    fPhotoAbsorptionCof4[c] = 0. ;
277  }
278  c = 1 ;
279
280  for(i = 0 ; i < noElm ; i++)
281  {
282    G4double I1 = fIonizationPotentials[Z[i]]*keV ;  // First ionization
283    n1 = 1 ;                                     // potential in keV
284 
285    for(j = 1 ; j < Z[i] ; j++)
286    {
287      n1 += fNbOfIntervals[j] ;
288    }
289    G4int n2 = n1 + fNbOfIntervals[Z[i]] ;
290   
291    for(k1 = n1 ; k1 < n2 ; k1++)
292    {
293      if(I1  > fSandiaTable[k1][0])
294      {
295         continue ;    // no ionization for energies smaller than I1 (first
296      }                // ionisation potential)             
297      break ;
298    }
299    G4int flag = 0 ;
300   
301    for(c1 = 1 ; c1 < c ; c1++)
302    {
303      if(fPhotoAbsorptionCof0[c1] == I1) // this value already has existed
304      {
305        flag = 1 ;                     
306        break ;                         
307      }
308    }
309    if(flag == 0)
310    {
311      fPhotoAbsorptionCof0[c] = I1 ;
312      c++ ;
313    }
314    for(k2 = k1 ; k2 < n2 ; k2++)
315    {
316      flag = 0 ;
317
318      for(c1 = 1 ; c1 < c ; c1++)
319      {
320        if(fPhotoAbsorptionCof0[c1] == fSandiaTable[k2][0])
321        {
322          flag = 1 ;
323          break ;
324        }
325      }
326      if(flag == 0)
327      {
328        fPhotoAbsorptionCof0[c] = fSandiaTable[k2][0] ;
329        c++ ;
330      }
331    }       
332  }   // end for(i)
333
334  // sort out
335  for(i = 1; i < c ; i++ ) 
336    {
337      for(j = i + 1 ; j < c;  j++ )
338        {
339          if(fPhotoAbsorptionCof0[i] > fPhotoAbsorptionCof0[j]) 
340            {
341              G4double tmp = fPhotoAbsorptionCof0[i];
342              fPhotoAbsorptionCof0[i] = fPhotoAbsorptionCof0[j];
343              fPhotoAbsorptionCof0[j] = tmp ;
344            }
345        }
346    }
347 
348  fMaxInterval = c ;
349 
350  const G4double* fractionW = fMaterial->GetFractionVector();
351   
352  for(i = 0 ; i < noElm; i++)
353    {
354      n1 = 1 ;
355      G4double I1 = fIonizationPotentials[Z[i]]*keV ;
356
357      for(j = 1 ; j < Z[i] ; j++)
358        {
359          n1 += fNbOfIntervals[j] ;
360        }
361      G4int n2 = n1 + fNbOfIntervals[Z[i]] - 1 ;
362
363      for(k = n1 ; k < n2 ; k++)
364        {
365          G4double B1 = fSandiaTable[k][0];
366          G4double B2 = fSandiaTable[k+1][0];
367          for(G4int c = 1 ; c < fMaxInterval-1 ; c++)
368            {
369              G4double E1 = fPhotoAbsorptionCof0[c];
370              G4double E2 = fPhotoAbsorptionCof0[c+1];
371              if(B1 > E1 || B2 < E2 || E1 < I1)
372                {
373                  continue ;
374                }
375              fPhotoAbsorptionCof1[c] += fSandiaTable[k][1]*fractionW[i] ;
376              fPhotoAbsorptionCof2[c] += fSandiaTable[k][2]*fractionW[i] ;
377              fPhotoAbsorptionCof3[c] += fSandiaTable[k][3]*fractionW[i] ;
378              fPhotoAbsorptionCof4[c] += fSandiaTable[k][4]*fractionW[i] ;
379            } 
380        }   
381      // Last interval
382      fPhotoAbsorptionCof1[fMaxInterval-1] += fSandiaTable[k][1]*fractionW[i] ;
383      fPhotoAbsorptionCof2[fMaxInterval-1] += fSandiaTable[k][2]*fractionW[i] ;
384      fPhotoAbsorptionCof3[fMaxInterval-1] += fSandiaTable[k][3]*fractionW[i] ;
385      fPhotoAbsorptionCof4[fMaxInterval-1] += fSandiaTable[k][4]*fractionW[i] ;
386    }     // for(i)
387 
388  c = 0 ;     // Deleting of first intervals where all coefficients = 0
389
390  do                       
391    {
392      c++ ;
393
394      if( fPhotoAbsorptionCof1[c] != 0.0 ||
395          fPhotoAbsorptionCof2[c] != 0.0 ||
396          fPhotoAbsorptionCof3[c] != 0.0 || 
397          fPhotoAbsorptionCof4[c] != 0.0     )  continue ;
398       
399      for(jj = 2 ; jj < fMaxInterval ; jj++)
400        {
401          fPhotoAbsorptionCof0[jj-1] = fPhotoAbsorptionCof0[jj];
402          fPhotoAbsorptionCof1[jj-1] = fPhotoAbsorptionCof1[jj];
403          fPhotoAbsorptionCof2[jj-1] = fPhotoAbsorptionCof2[jj];
404          fPhotoAbsorptionCof3[jj-1] = fPhotoAbsorptionCof3[jj];
405          fPhotoAbsorptionCof4[jj-1] = fPhotoAbsorptionCof4[jj];
406        }
407      fMaxInterval-- ;
408      c-- ;
409    }
410  while(c < fMaxInterval - 1) ;
411       
412  // create the sandia matrix for this material
413   
414  fMatSandiaMatrixPAI = new G4OrderedTable();
415  G4double density = fMaterial->GetDensity();
416 
417  for (i = 0; i < fMaxInterval; i++)
418    {
419      fMatSandiaMatrixPAI->push_back(new G4DataVector(5,0.));
420    }                   
421  for (i = 0; i < fMaxInterval; i++)
422    {
423      (*(*fMatSandiaMatrixPAI)[i])[0] = fPhotoAbsorptionCof0[i+1];
424      (*(*fMatSandiaMatrixPAI)[i])[1] = fPhotoAbsorptionCof1[i+1]*density;
425      (*(*fMatSandiaMatrixPAI)[i])[2] = fPhotoAbsorptionCof2[i+1]*density;
426      (*(*fMatSandiaMatrixPAI)[i])[3] = fPhotoAbsorptionCof3[i+1]*density;
427      (*(*fMatSandiaMatrixPAI)[i])[4] = fPhotoAbsorptionCof4[i+1]*density;
428
429    }                       
430  delete [] Z;
431  delete [] fPhotoAbsorptionCof0;
432  delete [] fPhotoAbsorptionCof1;
433  delete [] fPhotoAbsorptionCof2;
434  delete [] fPhotoAbsorptionCof3;
435  delete [] fPhotoAbsorptionCof4;
436  return;
437}
438
439//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
440
441////////////////////////////////////////////////////////////////////////
442////////////////////////////////////////////////////////////////////////
443//
444// Methods for PAI model
445
446//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
447
448G4SandiaTable::G4SandiaTable(G4int matIndex)
449{ 
450  fMatSandiaMatrix = 0 ; 
451  fMatSandiaMatrixPAI = 0;
452  fPhotoAbsorptionCof = 0 ;
453
454  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
455  size_t numberOfMat = G4Material::GetNumberOfMaterials();
456
457  if ( matIndex >= 0 && matIndex < G4int(numberOfMat) )
458  {
459    fMaterial = (*theMaterialTable)[matIndex];
460  }
461  else
462  {
463    G4Exception("G4SandiaTable::G4SandiaTable(G4int matIndex): wrong matIndex ");
464  }
465  ComputeMatTable();
466}
467
468///////////////////////////////////////////////////////////////////////
469//
470// Bubble sorting of left energy interval in SandiaTable in ascening order
471//
472
473void
474G4SandiaTable::SandiaSort(G4double** da ,
475                          G4int sz )
476{
477   for(G4int i = 1 ;i < sz ; i++ ) 
478   {
479     for(G4int j = i + 1 ;j < sz ; j++ )
480     {
481       if(da[i][0] > da[j][0]) 
482       {
483          SandiaSwap(da,i,j) ;
484       }
485     }
486   }
487}
488
489////////////////////////////////////////////////////////////////////////////
490//
491//  SandiaIntervals
492//
493
494G4int
495G4SandiaTable::SandiaIntervals(G4int Z[],
496                               G4int el )
497{
498  G4int c,i ;
499
500  fMaxInterval = 0 ;
501
502  for(i=0;i<el;i++)
503  {
504    fMaxInterval += fNbOfIntervals[Z[i]] ; 
505  }
506  fMaxInterval += 2 ;
507
508//  G4cout<<"fMaxInterval = "<<fMaxInterval<<G4endl ;
509
510  fPhotoAbsorptionCof = new G4double* [fMaxInterval] ;
511
512  for(i = 0 ; i < fMaxInterval ; i++) 
513  {
514     fPhotoAbsorptionCof[i] = new G4double[5] ;
515  }
516
517
518  //  for(c = 0 ; c < fIntervalLimit ; c++)   // just in case
519
520  for(c = 0 ; c < fMaxInterval ; c++)   // just in case
521  {
522     fPhotoAbsorptionCof[c][0] = 0. ;
523  }
524  c = 1 ;
525  for(i = 0 ; i < el ; i++)
526  {
527    G4double I1 = fIonizationPotentials[Z[i]]*keV ;  // First ionization
528    G4int n1 = 1 ;                                     // potential in keV
529    G4int j, c1, k1, k2 ;
530    for(j = 1 ; j < Z[i] ; j++)
531    {
532      n1 += fNbOfIntervals[j] ;
533    }
534    G4int n2 = n1 + fNbOfIntervals[Z[i]] ;
535   
536    for(k1 = n1 ; k1 < n2 ; k1++)
537    {
538      if(I1  > fSandiaTable[k1][0])
539      {
540         continue ;    // no ionization for energies smaller than I1 (first
541      }                // ionisation potential)             
542      break ;
543    }
544    G4int flag = 0 ;
545   
546    for(c1 = 1 ; c1 < c ; c1++)
547    {
548      if(fPhotoAbsorptionCof[c1][0] == I1) // this value already has existed
549      {
550        flag = 1 ;                     
551        break ;                         
552      }
553    }
554    if(flag == 0)
555    {
556      fPhotoAbsorptionCof[c][0] = I1 ;
557      c++ ;
558    }
559    for(k2 = k1 ; k2 < n2 ; k2++)
560    {
561      flag = 0 ;
562      for(c1 = 1 ; c1 < c ; c1++)
563      {
564        if(fPhotoAbsorptionCof[c1][0] == fSandiaTable[k2][0])
565        {
566          flag = 1 ;
567          break ;
568        }
569      }
570      if(flag == 0)
571      {
572        fPhotoAbsorptionCof[c][0] = fSandiaTable[k2][0] ;
573        c++ ;
574      }
575    }       
576  }   // end for(i)
577 
578  SandiaSort(fPhotoAbsorptionCof,c) ;
579  fMaxInterval = c ;
580  return c ;
581}   
582
583///////////////////////////////////////////////////////////////////////
584//
585//  SandiaMixing
586//
587
588G4int
589G4SandiaTable::SandiaMixing(         G4int Z[],
590                               const G4double fractionW[],
591                                     G4int el,
592                                     G4int mi     )
593{
594   G4int i;
595   
596   for(i = 0 ; i < mi ; i++)
597   {
598      for(G4int j = 1 ; j < 5 ; j++) fPhotoAbsorptionCof[i][j] = 0. ;
599   }
600   for(i = 0 ; i < el ; i++)
601   {
602      G4int n1 = 1 ;
603      G4int j, k ;
604      G4double I1 = fIonizationPotentials[Z[i]]*keV ;
605      for(j = 1 ; j < Z[i] ; j++)
606      {
607         n1 += fNbOfIntervals[j] ;
608      }
609      G4int n2 = n1 + fNbOfIntervals[Z[i]] - 1 ;
610
611      for(k = n1 ; k < n2 ; k++)
612      {
613         G4double B1 = fSandiaTable[k][0] ;
614         G4double B2 = fSandiaTable[k+1][0] ;
615         for(G4int c = 1 ; c < mi-1 ; c++)
616         {
617            G4double E1 = fPhotoAbsorptionCof[c][0] ;
618            G4double E2 = fPhotoAbsorptionCof[c+1][0] ;
619            if(B1 > E1 || B2 < E2 || E1 < I1)
620            {
621               continue ;
622            }
623            for(j = 1 ; j < 5 ; j++)
624            {
625               fPhotoAbsorptionCof[c][j] += fSandiaTable[k][j]*fractionW[i] ;
626            }
627          } 
628       }   
629       for(j = 1 ; j < 5 ; j++)   // Last interval
630       {
631          fPhotoAbsorptionCof[mi-1][j] += fSandiaTable[k][j]*fractionW[i] ;
632       }
633    }     // for(i)
634    G4int c = 0 ;     // Deleting of first intervals where all coefficients = 0
635    do                       
636    {
637       c++ ;
638       if( fPhotoAbsorptionCof[c][1] != 0.0 ||
639           fPhotoAbsorptionCof[c][2] != 0.0 ||
640           fPhotoAbsorptionCof[c][3] != 0.0 || 
641           fPhotoAbsorptionCof[c][4] != 0.0     )
642       {
643          continue ;
644       }
645       for(G4int jj = 2 ; jj < mi ; jj++)
646       {
647          for(G4int kk = 0 ; kk < 5 ; kk++)
648          {
649             fPhotoAbsorptionCof[jj-1][kk]= fPhotoAbsorptionCof[jj][kk] ;
650          }
651       }
652       mi-- ;
653       c-- ;
654    }
655    while(c < mi - 1) ;
656   
657    return mi ;
658} 
659
660////////////////////////////////////////////////////////////////////////////
661//
662//  Sandia interval and mixing calculations for materialCutsCouple constructor
663//
664
665void G4SandiaTable::ComputeMatTable()
666{
667  G4int MaxIntervals = 0;
668  G4int elm, c, i, j, jj, k, kk, k1, k2, c1, n1;   
669
670  const G4int noElm = fMaterial->GetNumberOfElements();
671  const G4ElementVector* ElementVector = fMaterial->GetElementVector(); 
672  G4int* Z = new G4int[noElm];               //Atomic number
673
674  for (elm = 0; elm<noElm; elm++)
675  { 
676    Z[elm] = (G4int)(*ElementVector)[elm]->GetZ();
677    MaxIntervals += fNbOfIntervals[Z[elm]];
678  } 
679  fMaxInterval = 0 ;
680
681  for(i = 0; i < noElm; i++)  fMaxInterval += fNbOfIntervals[Z[i]] ; 
682 
683  fMaxInterval += 2 ;
684
685//  G4cout<<"fMaxInterval = "<<fMaxInterval<<G4endl ;
686
687  fPhotoAbsorptionCof = new G4double* [fMaxInterval] ;
688
689  for(i = 0 ; i < fMaxInterval ; i++) 
690  {
691     fPhotoAbsorptionCof[i] = new G4double[5] ;
692  }
693
694  //  for(c = 0 ; c < fIntervalLimit ; c++)   // just in case
695
696  for(c = 0 ; c < fMaxInterval ; c++)   // just in case
697  {
698     fPhotoAbsorptionCof[c][0] = 0. ;
699  }
700  c = 1 ;
701
702  for(i = 0 ; i < noElm ; i++)
703  {
704    G4double I1 = fIonizationPotentials[Z[i]]*keV ;  // First ionization
705    n1 = 1 ;                                     // potential in keV
706 
707    for(j = 1 ; j < Z[i] ; j++)
708    {
709      n1 += fNbOfIntervals[j] ;
710    }
711    G4int n2 = n1 + fNbOfIntervals[Z[i]] ;
712   
713    for(k1 = n1 ; k1 < n2 ; k1++)
714    {
715      if(I1  > fSandiaTable[k1][0])
716      {
717         continue ;    // no ionization for energies smaller than I1 (first
718      }                // ionisation potential)             
719      break ;
720    }
721    G4int flag = 0 ;
722   
723    for(c1 = 1 ; c1 < c ; c1++)
724    {
725      if(fPhotoAbsorptionCof[c1][0] == I1) // this value already has existed
726      {
727        flag = 1 ;                     
728        break ;                         
729      }
730    }
731    if(flag == 0)
732    {
733      fPhotoAbsorptionCof[c][0] = I1 ;
734      c++ ;
735    }
736    for(k2 = k1 ; k2 < n2 ; k2++)
737    {
738      flag = 0 ;
739
740      for(c1 = 1 ; c1 < c ; c1++)
741      {
742        if(fPhotoAbsorptionCof[c1][0] == fSandiaTable[k2][0])
743        {
744          flag = 1 ;
745          break ;
746        }
747      }
748      if(flag == 0)
749      {
750        fPhotoAbsorptionCof[c][0] = fSandiaTable[k2][0] ;
751        c++ ;
752      }
753    }       
754  }   // end for(i)
755 
756  SandiaSort(fPhotoAbsorptionCof,c) ;
757  fMaxInterval = c ;
758 
759  const G4double* fractionW = fMaterial->GetFractionVector();
760   
761  for(i = 0 ; i < fMaxInterval ; i++)
762  {
763      for(j = 1 ; j < 5 ; j++) fPhotoAbsorptionCof[i][j] = 0.;
764  }
765  for(i = 0 ; i < noElm; i++)
766  {
767      n1 = 1 ;
768      G4double I1 = fIonizationPotentials[Z[i]]*keV ;
769
770      for(j = 1 ; j < Z[i] ; j++)
771      {
772         n1 += fNbOfIntervals[j] ;
773      }
774      G4int n2 = n1 + fNbOfIntervals[Z[i]] - 1 ;
775
776      for(k = n1 ; k < n2 ; k++)
777      {
778         G4double B1 = fSandiaTable[k][0] ;
779         G4double B2 = fSandiaTable[k+1][0] ;
780         for(G4int c = 1 ; c < fMaxInterval-1 ; c++)
781         {
782            G4double E1 = fPhotoAbsorptionCof[c][0] ;
783            G4double E2 = fPhotoAbsorptionCof[c+1][0] ;
784            if(B1 > E1 || B2 < E2 || E1 < I1)
785            {
786               continue ;
787            }
788            for(j = 1 ; j < 5 ; j++)
789            {
790               fPhotoAbsorptionCof[c][j] += fSandiaTable[k][j]*fractionW[i] ;
791            }
792          } 
793       }   
794       for(j = 1 ; j < 5 ; j++)   // Last interval
795       {
796          fPhotoAbsorptionCof[fMaxInterval-1][j] += fSandiaTable[k][j]*fractionW[i] ;
797       }
798  }     // for(i)
799
800  c = 0 ;     // Deleting of first intervals where all coefficients = 0
801
802  do                       
803  {
804    c++ ;
805
806    if( fPhotoAbsorptionCof[c][1] != 0.0 ||
807        fPhotoAbsorptionCof[c][2] != 0.0 ||
808        fPhotoAbsorptionCof[c][3] != 0.0 || 
809        fPhotoAbsorptionCof[c][4] != 0.0     )  continue ;
810       
811    for(jj = 2 ; jj < fMaxInterval ; jj++)
812    {
813      for(kk = 0 ; kk < 5 ; kk++)
814      {
815             fPhotoAbsorptionCof[jj-1][kk]= fPhotoAbsorptionCof[jj][kk] ;
816      }
817    }
818    fMaxInterval-- ;
819    c-- ;
820  }
821  while(c < fMaxInterval - 1) ;
822       
823  // create the sandia matrix for this material
824   
825  fMatSandiaMatrix = new G4OrderedTable();
826 
827  for (i = 0; i < fMaxInterval; i++)
828  {
829     fMatSandiaMatrix->push_back(new G4DataVector(5,0.));
830  }                     
831  for (i = 0; i < fMaxInterval; i++)
832  {
833    for(j = 0 ; j < 5 ; j++)
834    {
835      (*(*fMatSandiaMatrix)[i])[j] = fPhotoAbsorptionCof[i+1][j];
836    }     
837  }                         
838  delete [] Z;
839  return ;
840} 
841
842//     G4SandiaTable class -- end of implementation file
843//
844////////////////////////////////////////////////////////////////////////////
845
846
Note: See TracBrowser for help on using the repository browser.