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

Last change on this file since 826 was 822, checked in by garnier, 17 years ago

import all except CVS
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: $
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.