G4SandiaTable.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 
00027 // $Id: G4SandiaTable.cc 67044 2013-01-30 08:50:06Z gcosmo $
00028 
00029 //
00030 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
00031 //
00032 // 10.06.97 created. V. Grichine
00033 // 18.11.98 simplified public interface; new methods for materials.  mma
00034 // 31.01.01 redesign of ComputeMatSandiaMatrix().  mma
00035 // 16.02.01 adapted for STL.  mma
00036 // 22.02.01 GetsandiaCofForMaterial(energy) return 0 below lowest interval  mma  
00037 // 03.04.01 fnulcof returned if energy < emin
00038 // 10.07.01 Migration to STL. M. Verderi.
00039 // 03.02.04 Update distructor V.Ivanchenko
00040 // 05.03.04 New methods for old sorting algorithm for PAI model. V.Grichine
00041 // 26.10.11 new scheme for G4Exception  (mma) 
00042 //
00043 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
00044 
00045 
00046 #include "G4SandiaTable.hh"
00047 #include "G4StaticSandiaData.hh"
00048 #include "G4Material.hh"
00049 #include "G4MaterialTable.hh"
00050 #include "G4PhysicalConstants.hh"
00051 #include "G4SystemOfUnits.hh"
00052 
00053 G4int    G4SandiaTable::fCumulInterval[101]  = {0};
00054 G4double G4SandiaTable::fSandiaCofPerAtom[4] = {0.0};
00055 G4double const G4SandiaTable::funitc[5] = {keV,
00056                                            cm2*keV/g,     
00057                                            cm2*keV*keV/g,     
00058                                            cm2*keV*keV*keV/g,     
00059                                            cm2*keV*keV*keV*keV/g};
00060  
00061 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
00062 
00063 G4SandiaTable::G4SandiaTable(G4Material* material)
00064   : fMaterial(material)
00065 {
00066   fMatSandiaMatrix    = 0; 
00067   fMatSandiaMatrixPAI = 0;
00068   fPhotoAbsorptionCof = 0;
00069 
00070   fMatNbOfIntervals   = 0;
00071 
00072  
00073   fMaxInterval        = 0;
00074   fVerbose            = 0;  
00075 
00076   //build the CumulInterval array
00077 
00078   fCumulInterval[0] = 1;
00079 
00080   for (G4int Z=1; Z<101; ++Z) {
00081     fCumulInterval[Z] = fCumulInterval[Z-1] + fNbOfIntervals[Z];
00082   }
00083   
00084   //initialisation of fnulcof
00085   fnulcof[0] = fnulcof[1] = fnulcof[2] = fnulcof[3] = 0.;
00086 
00087   fMaxInterval = 0;
00088 
00089   //compute macroscopic Sandia coefs for a material   
00090   ComputeMatSandiaMatrix(); // mma
00091 
00092   // ComputeMatTable();  // vmg
00093 }
00094                                                         
00095 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
00096 
00097 // Fake default constructor - sets only member data and allocates memory
00098 //                            for usage restricted to object persistency
00099 
00100 G4SandiaTable::G4SandiaTable(__void__&)
00101   : fMaterial(0),fMatSandiaMatrix(0),fMatSandiaMatrixPAI(0),fPhotoAbsorptionCof(0)
00102 {
00103   fnulcof[0] = fnulcof[1] = fnulcof[2] = fnulcof[3] = 0.;
00104   fMaxInterval = 0;
00105   fMatNbOfIntervals = 0;
00106   fVerbose          = 0;  
00107 }
00108 
00109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
00110 
00111 G4SandiaTable::~G4SandiaTable()
00112 { 
00113   if(fMatSandiaMatrix) {
00114     fMatSandiaMatrix->clearAndDestroy();
00115     delete fMatSandiaMatrix;
00116   }
00117   if(fMatSandiaMatrixPAI) {
00118     fMatSandiaMatrixPAI->clearAndDestroy();
00119     delete fMatSandiaMatrixPAI;
00120   }
00121   if(fPhotoAbsorptionCof)
00122   {
00123     delete [] fPhotoAbsorptionCof;
00124   }
00125 }
00126 
00127 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
00128 
00129 G4double*
00130 G4SandiaTable::GetSandiaCofPerAtom(G4int Z, G4double energy)
00131 {
00132   G4double Emin  = fSandiaTable[fCumulInterval[Z-1]][0]*keV;
00133   G4double Iopot = fIonizationPotentials[Z]*eV;
00134   if (Iopot > Emin) Emin = Iopot;
00135    
00136   G4int interval = fNbOfIntervals[Z] - 1;
00137   G4int row = fCumulInterval[Z-1] + interval;
00138   while ((interval>0) && (energy<fSandiaTable[row][0]*keV)) {
00139     --interval;
00140     row = fCumulInterval[Z-1] + interval;
00141   }
00142   if (energy >= Emin)
00143     {        
00144       G4double AoverAvo = Z*amu/fZtoAratio[Z];
00145          
00146       fSandiaCofPerAtom[0]=AoverAvo*funitc[1]*fSandiaTable[row][1];     
00147       fSandiaCofPerAtom[1]=AoverAvo*funitc[2]*fSandiaTable[row][2];     
00148       fSandiaCofPerAtom[2]=AoverAvo*funitc[3]*fSandiaTable[row][3];     
00149       fSandiaCofPerAtom[3]=AoverAvo*funitc[4]*fSandiaTable[row][4];
00150     }
00151   else 
00152     {
00153       fSandiaCofPerAtom[0] = fSandiaCofPerAtom[1] = fSandiaCofPerAtom[2] =
00154         fSandiaCofPerAtom[3] = 0.;
00155     }                
00156   return fSandiaCofPerAtom;     
00157 }
00158 
00159 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
00160                                                         
00161 G4double G4SandiaTable::GetZtoA(G4int Z)
00162 {
00163   assert (Z>0 && Z<101);
00164   return fZtoAratio[Z];
00165 }
00166 
00167 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
00168 
00169 void G4SandiaTable::ComputeMatSandiaMatrix()
00170 {  
00171   //get list of elements
00172   const G4int NbElm = fMaterial->GetNumberOfElements();
00173   const G4ElementVector* ElementVector = fMaterial->GetElementVector();
00174   
00175   G4int* Z = new G4int[NbElm];               //Atomic number
00176    
00177   //determine the maximum number of energy-intervals for this material
00178   
00179   G4int MaxIntervals = 0;
00180   G4int elm;
00181     
00182   for ( elm = 0; elm < NbElm; elm++ ) 
00183   {
00184     Z[elm] = (G4int)(*ElementVector)[elm]->GetZ();
00185     MaxIntervals += fNbOfIntervals[Z[elm]];
00186   }  
00187      
00188   //copy the Energy bins in a tmp1 array
00189   //(take care of the Ionization Potential of each element)
00190   
00191   G4double* tmp1 = new G4double[MaxIntervals]; 
00192   G4double IonizationPot;
00193   G4int interval1 = 0;
00194 
00195   for ( elm = 0; elm < NbElm; elm++ ) 
00196   {
00197     IonizationPot = GetIonizationPot(Z[elm]);
00198 
00199     for (G4int row = fCumulInterval[ Z[elm]-1 ]; row < fCumulInterval[Z[elm]]; row++ )
00200     {
00201       tmp1[interval1++] = std::max(fSandiaTable[row][0]*keV,IonizationPot);
00202     }
00203   }   
00204         
00205   //sort the energies in strickly increasing values in a tmp2 array
00206   //(eliminate redondances)
00207   
00208   G4double* tmp2 = new G4double[MaxIntervals];
00209   G4double Emin;
00210   G4int interval2 = 0;
00211   
00212   do 
00213   {
00214     Emin = DBL_MAX;
00215 
00216     for ( G4int i1 = 0; i1 < MaxIntervals; i1++ ) 
00217     {
00218       if (tmp1[i1] < Emin) Emin = tmp1[i1];          //find the minimum
00219     }
00220     if (Emin < DBL_MAX) tmp2[interval2++] = Emin;
00221     //copy Emin in tmp2
00222     for ( G4int j1 = 0; j1 < MaxIntervals; j1++ ) 
00223     {
00224       if (tmp1[j1] <= Emin) tmp1[j1] = DBL_MAX;      //eliminate from tmp1          
00225     }
00226   } while (Emin < DBL_MAX);
00227                                         
00228   //create the sandia matrix for this material
00229     
00230   fMatSandiaMatrix = new G4OrderedTable();
00231   G4int interval;
00232 
00233   for (interval = 0; interval < interval2; interval++ ) 
00234   {
00235     fMatSandiaMatrix->push_back( new G4DataVector(5,0.) );
00236   }
00237   
00238   //ready to compute the Sandia coefs for the material
00239   
00240   const G4double* NbOfAtomsPerVolume = fMaterial->GetVecNbOfAtomsPerVolume();
00241   
00242   const G4double prec = 1.e-03*eV;
00243   G4double coef, oldsum(0.), newsum(0.);
00244   fMatNbOfIntervals = 0;
00245          
00246   for ( interval = 0; interval < interval2; interval++ ) 
00247   {
00248     Emin = (*(*fMatSandiaMatrix)[fMatNbOfIntervals])[0] = tmp2[interval];
00249 
00250     for ( G4int k = 1; k < 5; k++ ) {
00251       (*(*fMatSandiaMatrix)[fMatNbOfIntervals])[k] = 0.;      
00252     }
00253     newsum = 0.;
00254       
00255     for ( elm = 0; elm < NbElm; elm++ ) 
00256     {    
00257       GetSandiaCofPerAtom(Z[elm], Emin+prec);
00258 
00259       for ( G4int j = 1; j < 5; j++ ) 
00260       {
00261         coef = NbOfAtomsPerVolume[elm]*fSandiaCofPerAtom[j-1];
00262         (*(*fMatSandiaMatrix)[fMatNbOfIntervals])[j] += coef;
00263         newsum += std::fabs(coef);
00264       }                                                
00265     }                                                            
00266     //check for null or redondant intervals
00267          
00268     if (newsum != oldsum) { oldsum = newsum; fMatNbOfIntervals++;}
00269   }
00270   delete [] Z;
00271   delete [] tmp1;
00272   delete [] tmp2;
00273 
00274   if ( fVerbose > 0 && fMaterial->GetName() == "G4_Ar" )
00275   {
00276     G4cout<<"mma, G4SandiaTable::ComputeMatSandiaMatrix(), mat = "
00277           <<fMaterial->GetName()<<G4endl;
00278 
00279     for( G4int i = 0; i < fMatNbOfIntervals; i++)
00280     {
00281       G4cout<<i<<"\t"<<GetSandiaCofForMaterial(i,0)/keV<<" keV \t"
00282             <<this->GetSandiaCofForMaterial(i,1)
00283             <<"\t"<<this->GetSandiaCofForMaterial(i,2)
00284             <<"\t"<<this->GetSandiaCofForMaterial(i,3)
00285             <<"\t"<<this->GetSandiaCofForMaterial(i,4)<<G4endl;
00286     }   
00287   }
00288 }
00289 
00290 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
00291 
00292 void G4SandiaTable::ComputeMatSandiaMatrixPAI()
00293 {  
00294   G4int MaxIntervals = 0;
00295   G4int elm, c, i, j, jj, k, k1, k2, c1, n1;    
00296 
00297   const G4int noElm = fMaterial->GetNumberOfElements();
00298   const G4ElementVector* ElementVector = fMaterial->GetElementVector();  
00299   G4int* Z = new G4int[noElm];               //Atomic number
00300 
00301   for ( elm = 0; elm < noElm; elm++ )
00302   { 
00303     Z[elm] = (G4int)(*ElementVector)[elm]->GetZ();
00304     MaxIntervals += fNbOfIntervals[Z[elm]];
00305   }  
00306   fMaxInterval = MaxIntervals + 2;
00307 
00308   if ( fVerbose > 0 && fMaterial->GetName() == "G4_Ar" )
00309   {
00310     G4cout<<"fMaxInterval = "<<fMaxInterval<<G4endl;
00311   } 
00312   G4double* fPhotoAbsorptionCof0 = new G4double[fMaxInterval];
00313   G4double* fPhotoAbsorptionCof1 = new G4double[fMaxInterval];
00314   G4double* fPhotoAbsorptionCof2 = new G4double[fMaxInterval];
00315   G4double* fPhotoAbsorptionCof3 = new G4double[fMaxInterval];
00316   G4double* fPhotoAbsorptionCof4 = new G4double[fMaxInterval];
00317 
00318   for( c = 0; c < fMaxInterval; c++ )   // just in case
00319   {
00320     fPhotoAbsorptionCof0[c] = 0.;
00321     fPhotoAbsorptionCof1[c] = 0.;
00322     fPhotoAbsorptionCof2[c] = 0.;
00323     fPhotoAbsorptionCof3[c] = 0.;
00324     fPhotoAbsorptionCof4[c] = 0.;
00325   }
00326   c = 1;
00327 
00328   for(i = 0; i < noElm; i++)
00329   {
00330     G4double I1 = fIonizationPotentials[Z[i]]*keV;  // I1 in keV
00331     n1          = 1;                                    
00332  
00333     for( j = 1; j < Z[i]; j++ )  n1 += fNbOfIntervals[j];
00334     
00335     G4int n2 = n1 + fNbOfIntervals[Z[i]];
00336     
00337     for( k1 = n1; k1 < n2; k1++ )
00338     {
00339       if( I1  > fSandiaTable[k1][0] )
00340       {
00341          continue;    // no ionization for energies smaller than I1 (first
00342       }                // ionisation potential)              
00343       break;
00344     }
00345     G4int flag = 0;
00346     
00347     for( c1 = 1; c1 < c; c1++ )
00348     {
00349       if( fPhotoAbsorptionCof0[c1] == I1 ) // this value already has existed
00350       {
00351         flag = 1;                      
00352         break;                         
00353       }
00354     }
00355     if(flag == 0)
00356     {
00357       fPhotoAbsorptionCof0[c] = I1;
00358       c++;
00359     }
00360     for( k2 = k1; k2 < n2; k2++ )
00361     {
00362       flag = 0;
00363 
00364       for( c1 = 1; c1 < c; c1++ )
00365       {
00366         if( fPhotoAbsorptionCof0[c1] == fSandiaTable[k2][0] )
00367         {
00368           flag = 1;
00369           break;
00370         }
00371       }
00372       if(flag == 0)
00373       {
00374         fPhotoAbsorptionCof0[c] = fSandiaTable[k2][0];
00375         c++;
00376       }
00377     }       
00378   }   // end for(i)
00379   // sort out
00380 
00381   for( i = 1; i < c; i++ ) 
00382   {
00383     for( j = i + 1; j < c;  j++ )
00384     {
00385       if( fPhotoAbsorptionCof0[i] > fPhotoAbsorptionCof0[j] ) 
00386       {
00387         G4double tmp = fPhotoAbsorptionCof0[i];
00388         fPhotoAbsorptionCof0[i] = fPhotoAbsorptionCof0[j];
00389         fPhotoAbsorptionCof0[j] = tmp;
00390       }
00391     }
00392     if ( fVerbose > 0 && fMaterial->GetName() == "G4_Ar" )
00393     {
00394       G4cout<<i<<"\t energy = "<<fPhotoAbsorptionCof0[i]<<G4endl;
00395     }
00396   } 
00397   fMaxInterval = c;
00398  
00399   const G4double* fractionW = fMaterial->GetFractionVector();
00400 
00401   if ( fVerbose > 0 && fMaterial->GetName() == "G4_Ar" )
00402   {
00403     for( i = 0; i < noElm; i++ )  G4cout<<i<<" = elN, fraction = "<<fractionW[i]<<G4endl;
00404   }      
00405    
00406   for( i = 0; i < noElm; i++ )
00407   {
00408     n1 = 1;
00409     G4double I1 = fIonizationPotentials[Z[i]]*keV;
00410 
00411     for( j = 1; j < Z[i]; j++ )  n1 += fNbOfIntervals[j];
00412         
00413     G4int n2 = n1 + fNbOfIntervals[Z[i]] - 1;
00414 
00415     for(k = n1; k < n2; k++)
00416     {
00417       G4double B1 = fSandiaTable[k][0];
00418       G4double B2 = fSandiaTable[k+1][0];
00419 
00420       for(G4int q = 1; q < fMaxInterval-1; q++)
00421       {
00422         G4double E1 = fPhotoAbsorptionCof0[q];
00423         G4double E2 = fPhotoAbsorptionCof0[q+1];
00424 
00425         if ( fVerbose > 0 && fMaterial->GetName() == "G4_Ar" )
00426         {
00427           G4cout<<"k = "<<k<<", q = "<<q<<", B1 = "<<B1<<", B2 = "<<B2<<", E1 = "<<E1<<", E2 = "<<E2<<G4endl;
00428         }      
00429         if( B1 > E1 || B2 < E2 || E1 < I1 )  
00430         {
00431           if ( fVerbose > 0 && fMaterial->GetName() == "G4_Ar" )
00432           {
00433             G4cout<<"continue for: B1 = "<<B1<<", B2 = "<<B2<<", E1 = "<<E1<<", E2 = "<<E2<<G4endl;
00434           }      
00435           continue;
00436         }               
00437         fPhotoAbsorptionCof1[q] += fSandiaTable[k][1]*fractionW[i];
00438         fPhotoAbsorptionCof2[q] += fSandiaTable[k][2]*fractionW[i];
00439         fPhotoAbsorptionCof3[q] += fSandiaTable[k][3]*fractionW[i];
00440         fPhotoAbsorptionCof4[q] += fSandiaTable[k][4]*fractionW[i];
00441       }  
00442     }   
00443     // Last interval
00444 
00445     fPhotoAbsorptionCof1[fMaxInterval-1] += fSandiaTable[k][1]*fractionW[i];
00446     fPhotoAbsorptionCof2[fMaxInterval-1] += fSandiaTable[k][2]*fractionW[i];
00447     fPhotoAbsorptionCof3[fMaxInterval-1] += fSandiaTable[k][3]*fractionW[i];
00448     fPhotoAbsorptionCof4[fMaxInterval-1] += fSandiaTable[k][4]*fractionW[i];
00449   }     // for(i)  
00450   c = 0;     // Deleting of first intervals where all coefficients = 0
00451 
00452   do                        
00453   {
00454     c++;
00455 
00456     if( fPhotoAbsorptionCof1[c] != 0.0 ||
00457         fPhotoAbsorptionCof2[c] != 0.0 ||
00458         fPhotoAbsorptionCof3[c] != 0.0 || 
00459         fPhotoAbsorptionCof4[c] != 0.0     )  continue;
00460 
00461     if ( fVerbose > 0 && fMaterial->GetName() == "G4_Ar" )
00462     {
00463       G4cout<<c<<" = number with zero cofs"<<G4endl;
00464     }      
00465     for( jj = 2; jj < fMaxInterval; jj++ )
00466     {
00467       fPhotoAbsorptionCof0[jj-1] = fPhotoAbsorptionCof0[jj];
00468       fPhotoAbsorptionCof1[jj-1] = fPhotoAbsorptionCof1[jj];
00469       fPhotoAbsorptionCof2[jj-1] = fPhotoAbsorptionCof2[jj];
00470       fPhotoAbsorptionCof3[jj-1] = fPhotoAbsorptionCof3[jj];
00471       fPhotoAbsorptionCof4[jj-1] = fPhotoAbsorptionCof4[jj];
00472     }
00473     fMaxInterval--;
00474     c--;
00475   }
00476   while( c < fMaxInterval - 1 );
00477         
00478   // create the sandia matrix for this material
00479     
00480   fMatSandiaMatrixPAI = new G4OrderedTable();
00481   G4double density = fMaterial->GetDensity();
00482  
00483   for (i = 0; i < fMaxInterval; i++)  fMatSandiaMatrixPAI->push_back(new G4DataVector(5,0.));
00484                         
00485   for (i = 0; i < fMaxInterval; i++)
00486   {
00487     (*(*fMatSandiaMatrixPAI)[i])[0] = fPhotoAbsorptionCof0[i+1];
00488     (*(*fMatSandiaMatrixPAI)[i])[1] = fPhotoAbsorptionCof1[i+1]*density;
00489     (*(*fMatSandiaMatrixPAI)[i])[2] = fPhotoAbsorptionCof2[i+1]*density;
00490     (*(*fMatSandiaMatrixPAI)[i])[3] = fPhotoAbsorptionCof3[i+1]*density;
00491     (*(*fMatSandiaMatrixPAI)[i])[4] = fPhotoAbsorptionCof4[i+1]*density;
00492   }
00493   if ( fVerbose > 0 && fMaterial->GetName() == "G4_Ar" )
00494   {
00495     G4cout<<"mma, G4SandiaTable::ComputeMatSandiaMatrixPAI(), mat = "<<fMaterial->GetName()<<G4endl;
00496 
00497     for( i = 0; i < fMaxInterval; i++)
00498     {
00499       G4cout<<i<<"\t"<<GetSandiaMatTablePAI(i,0)/keV<<" keV \t"<<this->GetSandiaMatTablePAI(i,1)
00500        <<"\t"<<this->GetSandiaMatTablePAI(i,2)<<"\t"<<this->GetSandiaMatTablePAI(i,3)
00501        <<"\t"<<this->GetSandiaMatTablePAI(i,4)<<G4endl;
00502     }   
00503   }
00504                             
00505   delete [] Z;
00506   delete [] fPhotoAbsorptionCof0;
00507   delete [] fPhotoAbsorptionCof1;
00508   delete [] fPhotoAbsorptionCof2;
00509   delete [] fPhotoAbsorptionCof3;
00510   delete [] fPhotoAbsorptionCof4;
00511   return;
00512 }
00513 
00514 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
00515 
00518 //
00519 // Methods for PAI model
00520 
00521 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
00522 
00523 G4SandiaTable::G4SandiaTable(G4int matIndex)
00524 { 
00525   fMaterial           = 0;
00526   fMatNbOfIntervals   = 0;
00527   fMatSandiaMatrix    = 0; 
00528   fMatSandiaMatrixPAI = 0;
00529   fPhotoAbsorptionCof = 0;
00530 
00531  
00532   fMaxInterval        = 0;
00533   fVerbose            = 0;  
00534 
00535   //initialisation of fnulcof
00536   fnulcof[0] = fnulcof[1] = fnulcof[2] = fnulcof[3] = 0.;
00537 
00538   const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
00539   G4int numberOfMat = G4Material::GetNumberOfMaterials();
00540 
00541   if ( matIndex >= 0 && matIndex < numberOfMat)
00542     {
00543       fMaterial = (*theMaterialTable)[matIndex];
00544       ComputeMatTable();
00545     }
00546   else
00547     {
00548       G4Exception("G4SandiaTable::G4SandiaTable(G4int matIndex)", "mat401",
00549        FatalException, "wrong matIndex");
00550     }
00551 }
00552 
00554 //
00555 // Bubble sorting of left energy interval in SandiaTable in ascening order
00556 //
00557 
00558 void
00559 G4SandiaTable::SandiaSort(G4double** da ,
00560                           G4int sz )
00561 {
00562   for(G4int i = 1;i < sz; i++ ) 
00563    {
00564      for(G4int j = i + 1;j < sz; j++ )
00565      {
00566        if(da[i][0] > da[j][0])   SandiaSwap(da,i,j);      
00567      }
00568    }
00569 }
00570 
00572 //
00573 //  SandiaIntervals 
00574 //
00575 
00576 G4int
00577 G4SandiaTable::SandiaIntervals(G4int Z[],
00578                                G4int el )
00579 {
00580   G4int c,  i, flag = 0, n1 = 1;
00581   G4int j, c1, k1, k2;
00582   G4double I1;
00583   fMaxInterval = 0;
00584 
00585   for( i = 0; i < el; i++ )  fMaxInterval += fNbOfIntervals[ Z[i] ]; 
00586 
00587   fMaxInterval += 2;
00588 
00589   if( fVerbose > 0 ) G4cout<<"begin sanInt, fMaxInterval = "<<fMaxInterval<<G4endl;
00590 
00591   fPhotoAbsorptionCof = new G4double* [fMaxInterval];
00592 
00593   for( i = 0; i < fMaxInterval; i++ )   fPhotoAbsorptionCof[i] = new G4double[5];
00594  
00595   //  for(c = 0; c < fIntervalLimit; c++)   // just in case
00596 
00597   for( c = 0; c < fMaxInterval; c++ )     fPhotoAbsorptionCof[c][0] = 0.;
00598   
00599   c = 1;
00600 
00601   for( i = 0; i < el; i++ )
00602   {
00603     I1 = fIonizationPotentials[ Z[i] ]*keV;  // First ionization
00604     n1 = 1;                                  // potential in keV
00605 
00606     for( j = 1; j < Z[i]; j++ )  n1 += fNbOfIntervals[j];
00607     
00608     G4int n2 = n1 + fNbOfIntervals[Z[i]];
00609     
00610     for( k1 = n1; k1 < n2; k1++ )
00611     {
00612       if( I1  > fSandiaTable[k1][0] )
00613       {
00614          continue;    // no ionization for energies smaller than I1 (first
00615       }                // ionisation potential)              
00616       break;
00617     }
00618     flag = 0;
00619     
00620     for( c1 = 1; c1 < c; c1++ )
00621     {
00622       if( fPhotoAbsorptionCof[c1][0] == I1 ) // this value already has existed
00623       {
00624         flag = 1;                      
00625         break;                         
00626       }
00627     }
00628     if( flag == 0 )
00629     {
00630       fPhotoAbsorptionCof[c][0] = I1;
00631       c++;
00632     }
00633     for( k2 = k1; k2 < n2; k2++ )
00634     {
00635       flag = 0;
00636 
00637       for( c1 = 1; c1 < c; c1++ )
00638       {
00639         if( fPhotoAbsorptionCof[c1][0] == fSandiaTable[k2][0] )
00640         {
00641           flag = 1;
00642           break;
00643         }
00644       }
00645       if( flag == 0 )
00646       {
00647         fPhotoAbsorptionCof[c][0] = fSandiaTable[k2][0];
00648         if( fVerbose > 0 ) G4cout<<"sanInt, c = "<<c<<", E_c = "<<fPhotoAbsorptionCof[c][0]<<G4endl;
00649         c++;
00650       }
00651     }       
00652   }   // end for(i)
00653   
00654   SandiaSort(fPhotoAbsorptionCof,c);
00655   fMaxInterval = c;
00656   if( fVerbose > 0 ) G4cout<<"end SanInt, fMaxInterval = "<<fMaxInterval<<G4endl;
00657   return c;
00658 }   
00659 
00661 //
00662 //  SandiaMixing
00663 //
00664 
00665 G4int
00666 G4SandiaTable::SandiaMixing(         G4int Z[],
00667                                const G4double fractionW[],
00668                                      G4int el,
00669                                      G4int mi     )
00670 {
00671   G4int i, j, n1, k, c=1, jj, kk;
00672   G4double I1, B1, B2, E1, E2;
00673    
00674   for( i = 0; i < mi; i++ )
00675   {
00676     for( j = 1; j < 5; j++ ) fPhotoAbsorptionCof[i][j] = 0.;
00677   }
00678   for( i = 0; i < el; i++ )
00679   {
00680     n1 = 1;
00681     I1 = fIonizationPotentials[Z[i]]*keV;
00682 
00683     for( j = 1; j < Z[i]; j++ )   n1 += fNbOfIntervals[j];
00684       
00685     G4int n2 = n1 + fNbOfIntervals[Z[i]] - 1;
00686 
00687     for( k = n1; k < n2; k++ )
00688     {
00689       B1 = fSandiaTable[k][0];
00690       B2 = fSandiaTable[k+1][0];
00691 
00692       for( c = 1; c < mi-1; c++ )
00693       {
00694         E1 = fPhotoAbsorptionCof[c][0];
00695         E2 = fPhotoAbsorptionCof[c+1][0];
00696 
00697         if( B1 > E1 || B2 < E2 || E1 < I1 )   continue;
00698             
00699         for( j = 1; j < 5; j++ ) 
00700         {
00701           fPhotoAbsorptionCof[c][j] += fSandiaTable[k][j]*fractionW[i];
00702           if( fVerbose > 0 )
00703           {
00704             G4cout<<"c="<<c<<"; j="<<j<<"; fST="<<fSandiaTable[k][j]<<"; frW="<<fractionW[i]<<G4endl;
00705           }
00706         }           
00707       }  
00708     }   
00709     for( j = 1; j < 5; j++ )   // Last interval
00710     {
00711       fPhotoAbsorptionCof[mi-1][j] += fSandiaTable[k][j]*fractionW[i];
00712       if( fVerbose > 0 )
00713       {
00714         G4cout<<"mi-1="<<mi-1<<"; j="<<j<<"; fST="<<fSandiaTable[k][j]<<"; frW="<<fractionW[i]<<G4endl;
00715       }
00716     }
00717   }     // for(i)
00718   c = 0;     // Deleting of first intervals where all coefficients = 0
00719 
00720   do                        
00721   {
00722     c++;
00723 
00724     if( fPhotoAbsorptionCof[c][1] != 0.0 ||
00725         fPhotoAbsorptionCof[c][2] != 0.0 ||
00726         fPhotoAbsorptionCof[c][3] != 0.0 || 
00727         fPhotoAbsorptionCof[c][4] != 0.0     )  continue;
00728        
00729     for( jj = 2; jj < mi; jj++ )
00730     {
00731       for( kk = 0; kk < 5; kk++ ) fPhotoAbsorptionCof[jj-1][kk] = fPhotoAbsorptionCof[jj][kk];    
00732     }
00733     mi--;
00734     c--;
00735   }
00736   while( c < mi - 1 );
00737 
00738   if( fVerbose > 0 ) G4cout<<"end SanMix, mi = "<<mi<<G4endl;
00739     
00740   return mi;
00741 }  
00742 
00743 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
00744 
00745 G4int G4SandiaTable::GetMatNbOfIntervals()  
00746 {
00747   return fMatNbOfIntervals;
00748 }
00749 
00750 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
00751 
00752 G4int G4SandiaTable::GetNbOfIntervals(G4int Z)
00753 {
00754   assert (Z>0 && Z<101);  
00755   return fNbOfIntervals[Z];
00756 }
00757 
00758 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
00759 
00760 G4double
00761 G4SandiaTable::GetSandiaCofPerAtom(G4int Z, G4int interval, G4int j)
00762 {
00763   assert (Z>0 && Z<101 && interval>=0 && interval<fNbOfIntervals[Z]
00764           && j>=0 && j<5);
00765 
00766   G4int row = fCumulInterval[Z-1] + interval;
00767   G4double x = fSandiaTable[row][0]*CLHEP::keV;
00768   if (j > 0) {
00769     x = Z*CLHEP::amu/fZtoAratio[Z]*fSandiaTable[row][j]*funitc[j];     
00770   }
00771   return x;
00772 }
00773 
00774 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
00775 
00776 G4double  
00777 G4SandiaTable::GetSandiaCofForMaterial(G4int interval, G4int j)    
00778 {
00779   assert (interval>=0 && interval<fMatNbOfIntervals && j>=0 && j<5);
00780   return ((*(*fMatSandiaMatrix)[interval])[j]); 
00781 }
00782 
00783 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
00784 
00785 G4double* 
00786 G4SandiaTable::GetSandiaCofForMaterial(G4double energy)
00787 {
00788   G4double* x = fnulcof;
00789   if (energy >= (*(*fMatSandiaMatrix)[0])[0]) {
00790    
00791     G4int interval = fMatNbOfIntervals - 1;
00792     while ((interval>0)&&(energy<(*(*fMatSandiaMatrix)[interval])[0])) 
00793       {interval--;} 
00794     x = &((*(*fMatSandiaMatrix)[interval])[1]);
00795   }
00796   return x;
00797 }
00798 
00799 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
00800 
00801 G4double  
00802 G4SandiaTable::GetSandiaMatTable(G4int interval, G4int j) 
00803 {
00804   assert (interval >= 0 && interval < fMaxInterval && j >= 0 && j < 5 );
00805   return ((*(*fMatSandiaMatrix)[interval])[j])*funitc[j]; 
00806 }
00807 
00808 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
00809 
00810 G4double  
00811 G4SandiaTable::GetSandiaCofForMaterialPAI(G4int interval, G4int j)    
00812 {
00813   assert (interval>=0 && interval<fMatNbOfIntervals && j>=0 && j<5);
00814   if(!fMatSandiaMatrixPAI) ComputeMatSandiaMatrixPAI();
00815   return ((*(*fMatSandiaMatrixPAI)[interval])[j]); 
00816 }
00817 
00818 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
00819 
00820 G4double* 
00821 G4SandiaTable::GetSandiaCofForMaterialPAI(G4double energy)
00822 {
00823   if(!fMatSandiaMatrixPAI) ComputeMatSandiaMatrixPAI();
00824   G4double* x = fnulcof;
00825   if (energy >= (*(*fMatSandiaMatrixPAI)[0])[0]) {
00826    
00827     G4int interval = fMatNbOfIntervals - 1;
00828     while ((interval>0)&&(energy<(*(*fMatSandiaMatrixPAI)[interval])[0])) 
00829       {interval--;} 
00830     x = &((*(*fMatSandiaMatrixPAI)[interval])[1]);
00831   }
00832   return x;
00833 }
00834 
00835 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
00836 
00837 G4double  
00838 G4SandiaTable::GetSandiaMatTablePAI(G4int interval, G4int j) 
00839 {
00840   assert (interval >= 0 && interval < fMaxInterval && j >= 0 && j < 5 );
00841   if(!fMatSandiaMatrixPAI) { ComputeMatSandiaMatrixPAI(); }
00842   return ((*(*fMatSandiaMatrixPAI)[interval])[j])*funitc[j]; 
00843 }
00844 
00845 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
00846 
00847 G4double
00848 G4SandiaTable::GetIonizationPot(G4int Z)
00849 {
00850   assert (Z>0 && Z<101);
00851   return fIonizationPotentials[Z]*CLHEP::eV;
00852 }
00853 
00854 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
00855 
00856 G4OrderedTable*  
00857 G4SandiaTable::GetSandiaMatrixPAI()
00858 {
00859   if(!fMatSandiaMatrixPAI) { ComputeMatSandiaMatrixPAI(); }
00860   return fMatSandiaMatrixPAI;
00861 }
00862 
00864 //
00865 //  Sandia interval and mixing calculations for materialCutsCouple constructor 
00866 //
00867 
00868 void G4SandiaTable::ComputeMatTable()
00869 {
00870   G4int MaxIntervals = 0;
00871   G4int elm, c, i, j, jj, k, kk, k1, k2, c1, n1;    
00872 
00873   const G4int noElm = fMaterial->GetNumberOfElements();
00874   const G4ElementVector* ElementVector = fMaterial->GetElementVector();  
00875   G4int* Z = new G4int[noElm];               //Atomic number
00876 
00877   for (elm = 0; elm<noElm; elm++)
00878   { 
00879     Z[elm] = (G4int)(*ElementVector)[elm]->GetZ();
00880     MaxIntervals += fNbOfIntervals[Z[elm]];
00881   }  
00882   fMaxInterval = 0;
00883 
00884   for(i = 0; i < noElm; i++)  fMaxInterval += fNbOfIntervals[Z[i]]; 
00885   
00886   fMaxInterval += 2;
00887 
00888 //  G4cout<<"fMaxInterval = "<<fMaxInterval<<G4endl;
00889 
00890   fPhotoAbsorptionCof = new G4double* [fMaxInterval];
00891 
00892   for(i = 0; i < fMaxInterval; i++)  
00893   {
00894      fPhotoAbsorptionCof[i] = new G4double[5];
00895   }
00896 
00897   //  for(c = 0; c < fIntervalLimit; c++)   // just in case
00898 
00899   for(c = 0; c < fMaxInterval; c++)   // just in case
00900   {
00901      fPhotoAbsorptionCof[c][0] = 0.;
00902   }
00903   c = 1;
00904 
00905   for(i = 0; i < noElm; i++)
00906   {
00907     G4double I1 = fIonizationPotentials[Z[i]]*keV;  // First ionization
00908     n1 = 1;                                     // potential in keV
00909  
00910     for(j = 1; j < Z[i]; j++)
00911     {
00912       n1 += fNbOfIntervals[j];
00913     }
00914     G4int n2 = n1 + fNbOfIntervals[Z[i]];
00915     
00916     for(k1 = n1; k1 < n2; k1++)
00917     {
00918       if(I1  > fSandiaTable[k1][0])
00919       {
00920          continue;    // no ionization for energies smaller than I1 (first
00921       }                // ionisation potential)              
00922       break;
00923     }
00924     G4int flag = 0;
00925     
00926     for(c1 = 1; c1 < c; c1++)
00927     {
00928       if(fPhotoAbsorptionCof[c1][0] == I1) // this value already has existed
00929       {
00930         flag = 1;                      
00931         break;                         
00932       }
00933     }
00934     if(flag == 0)
00935     {
00936       fPhotoAbsorptionCof[c][0] = I1;
00937       c++;
00938     }
00939     for(k2 = k1; k2 < n2; k2++)
00940     {
00941       flag = 0;
00942 
00943       for(c1 = 1; c1 < c; c1++)
00944       {
00945         if(fPhotoAbsorptionCof[c1][0] == fSandiaTable[k2][0])
00946         {
00947           flag = 1;
00948           break;
00949         }
00950       }
00951       if(flag == 0)
00952       {
00953         fPhotoAbsorptionCof[c][0] = fSandiaTable[k2][0];
00954         c++;
00955       }
00956     }       
00957   }   // end for(i)
00958   
00959   SandiaSort(fPhotoAbsorptionCof,c);
00960   fMaxInterval = c;
00961  
00962   const G4double* fractionW = fMaterial->GetFractionVector();
00963    
00964   for(i = 0; i < fMaxInterval; i++)
00965   {
00966       for(j = 1; j < 5; j++) fPhotoAbsorptionCof[i][j] = 0.;
00967   }
00968   for(i = 0; i < noElm; i++)
00969   {
00970       n1 = 1;
00971       G4double I1 = fIonizationPotentials[Z[i]]*keV;
00972 
00973       for(j = 1; j < Z[i]; j++)
00974       {
00975          n1 += fNbOfIntervals[j];
00976       }
00977       G4int n2 = n1 + fNbOfIntervals[Z[i]] - 1;
00978 
00979       for(k = n1; k < n2; k++)
00980       {
00981          G4double B1 = fSandiaTable[k][0];
00982          G4double B2 = fSandiaTable[k+1][0];
00983          for(G4int q = 1; q < fMaxInterval-1; q++)
00984          {
00985             G4double E1 = fPhotoAbsorptionCof[q][0];
00986             G4double E2 = fPhotoAbsorptionCof[q+1][0];
00987             if(B1 > E1 || B2 < E2 || E1 < I1)
00988             {
00989                continue;
00990             }
00991             for(j = 1; j < 5; j++)
00992             {
00993                fPhotoAbsorptionCof[q][j] += fSandiaTable[k][j]*fractionW[i];
00994             }
00995           }  
00996        }   
00997        for(j = 1; j < 5; j++)   // Last interval
00998        {
00999           fPhotoAbsorptionCof[fMaxInterval-1][j] += fSandiaTable[k][j]*fractionW[i];
01000        }
01001   }     // for(i)
01002 
01003   c = 0;     // Deleting of first intervals where all coefficients = 0
01004 
01005   do                        
01006   {
01007     c++;
01008 
01009     if( fPhotoAbsorptionCof[c][1] != 0.0 ||
01010         fPhotoAbsorptionCof[c][2] != 0.0 ||
01011         fPhotoAbsorptionCof[c][3] != 0.0 || 
01012         fPhotoAbsorptionCof[c][4] != 0.0     )  continue;
01013        
01014     for(jj = 2; jj < fMaxInterval; jj++)
01015     {
01016       for(kk = 0; kk < 5; kk++)
01017       {
01018              fPhotoAbsorptionCof[jj-1][kk]= fPhotoAbsorptionCof[jj][kk];
01019       }
01020     }
01021     fMaxInterval--;
01022     c--;
01023   }
01024   while(c < fMaxInterval - 1);
01025         
01026   // create the sandia matrix for this material
01027 
01028   fMaxInterval--;  // vmg 20.11.10
01029     
01030   fMatSandiaMatrix = new G4OrderedTable();
01031  
01032   for (i = 0; i < fMaxInterval; i++)
01033   {
01034      fMatSandiaMatrix->push_back(new G4DataVector(5,0.));
01035   }                     
01036   for ( i = 0; i < fMaxInterval; i++ )
01037   {
01038     for( j = 0; j < 5; j++ )
01039     {
01040       (*(*fMatSandiaMatrix)[i])[j] = fPhotoAbsorptionCof[i+1][j];
01041     }     
01042   }
01043   fMatNbOfIntervals = fMaxInterval; 
01044                             
01045   if ( fVerbose > 0 )
01046   {
01047     G4cout<<"vmg, G4SandiaTable::ComputeMatTable(), mat = "<<fMaterial->GetName()<<G4endl;
01048 
01049     for ( i = 0; i < fMaxInterval; i++ )
01050     {
01051       // G4cout<<i<<"\t"<<(*(*fMatSandiaMatrix)[i])[0]<<" keV \t"<<(*(*fMatSandiaMatrix)[i])[1]
01052       //       <<"\t"<<(*(*fMatSandiaMatrix)[i])[2]<<"\t"<<(*(*fMatSandiaMatrix)[i])[3]
01053       //   <<"\t"<<(*(*fMatSandiaMatrix)[i])[4]<<G4endl;
01054 
01055       G4cout<<i<<"\t"<<GetSandiaCofForMaterial(i,0)/keV<<" keV \t"<<this->GetSandiaCofForMaterial(i,1)
01056        <<"\t"<<this->GetSandiaCofForMaterial(i,2)<<"\t"<<this->GetSandiaCofForMaterial(i,3)
01057        <<"\t"<<this->GetSandiaCofForMaterial(i,4)<<G4endl;
01058     }
01059   }                         
01060   delete [] Z;
01061   return;
01062 }  
01063 
01064 //     G4SandiaTable class -- end of implementation file
01065 //
01067 
01068 

Generated on Mon May 27 17:49:47 2013 for Geant4 by  doxygen 1.4.7