Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/materials/src/G4SandiaTable.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /materials/src/G4SandiaTable.cc (Version 11.3.0) and /materials/src/G4SandiaTable.cc (Version 10.4)


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