Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4PenelopeBremsstrahlungAngular.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 /processes/electromagnetic/lowenergy/src/G4PenelopeBremsstrahlungAngular.cc (Version 11.3.0) and /processes/electromagnetic/lowenergy/src/G4PenelopeBremsstrahlungAngular.cc (Version 8.1)


  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 //                                                 25 //
 26 //                                             <<  26 // $Id: G4PenelopeBremsstrahlungAngular.cc,v 1.7 2006/06/29 19:40:37 gunter Exp $
                                                   >>  27 // GEANT4 tag $Name: geant4-08-01 $
                                                   >>  28 // 
 27 // -------------------------------------------     29 // --------------------------------------------------------------
 28 //                                                 30 //
 29 // File name:     G4PenelopeBremsstrahlungAngu     31 // File name:     G4PenelopeBremsstrahlungAngular
 30 //                                                 32 //
 31 // Author:        Luciano Pandola                  33 // Author:        Luciano Pandola
 32 //                                             <<  34 // 
 33 // Creation date: November 2010                <<  35 // Creation date: February 2003
 34 //                                                 36 //
 35 // History:                                        37 // History:
 36 // -----------                                     38 // -----------
 37 // 23 Nov 2010  L. Pandola       1st implement <<  39 // 04 Feb 2003  L. Pandola       1st implementation
 38 // 24 May 2011  L. Pandola       Renamed (make <<  40 // 19 Mar 2003  L. Pandola       Bugs fixed
 39 // 13 Mar 2012  L. Pandola       Made a derive <<  41 // 07 Nov 2003  L. Pandola       Added GetAtomicNumber method for testing 
 40 //                               and update th <<  42 //                               purposes
 41 // 18 Jul 2012  L. Pandola       Migrated to t << 
 42 //                               Now returns a << 
 43 // 03 Oct 2013  L. Pandola       Migrated to M << 
 44 // 17 Oct 2013  L. Pandola       Partially rev << 
 45 //                                thread-local << 
 46 //                                             << 
 47 //--------------------------------------------     43 //----------------------------------------------------------------
 48                                                    44 
 49 #include "G4PenelopeBremsstrahlungAngular.hh"      45 #include "G4PenelopeBremsstrahlungAngular.hh"
 50                                                <<  46 #include "G4PenelopeInterpolator.hh"
 51 #include "globals.hh"                          << 
 52 #include "G4PhysicalConstants.hh"              << 
 53 #include "G4SystemOfUnits.hh"                  << 
 54 #include "G4PhysicsFreeVector.hh"              << 
 55 #include "G4PhysicsTable.hh"                   << 
 56 #include "G4Material.hh"                       << 
 57 #include "Randomize.hh"                            47 #include "Randomize.hh"
 58 #include "G4Exp.hh"                            <<  48 #include "globals.hh"
 59                                                << 
 60 G4PenelopeBremsstrahlungAngular::G4PenelopeBre << 
 61   G4VEmAngularDistribution("Penelope"), fEffec << 
 62   fLorentzTables1(nullptr),fLorentzTables2(nul << 
 63                                                    49 
                                                   >>  50 G4PenelopeBremsstrahlungAngular::G4PenelopeBremsstrahlungAngular (G4int Zed)
                                                   >>  51   : Zmat(Zed)
 64 {                                                  52 {
 65   fDataRead = false;                           <<  53   InterpolationTableForZ();
 66   fVerbosityLevel = 0;                         <<  54   InterpolationForK();
 67 }                                                  55 }
 68                                                    56 
 69 //....oooOO0OOooo........oooOO0OOooo........oo << 
 70                                                    57 
 71 G4PenelopeBremsstrahlungAngular::~G4PenelopeBr     58 G4PenelopeBremsstrahlungAngular::~G4PenelopeBremsstrahlungAngular()
 72 {                                                  59 {
 73   ClearTables();                               << 
 74 }                                                  60 }
 75                                                    61 
 76 //....oooOO0OOooo........oooOO0OOooo........oo <<  62 G4int G4PenelopeBremsstrahlungAngular::GetAtomicNumber()
 77                                                << 
 78 void G4PenelopeBremsstrahlungAngular::Initiali << 
 79 {                                                  63 {
 80   ClearTables();                               <<  64   return Zmat;
 81 }                                                  65 }
 82                                                    66 
 83 //....oooOO0OOooo........oooOO0OOooo........oo <<  67 void G4PenelopeBremsstrahlungAngular::InterpolationTableForZ()
 84                                                << 
 85 void G4PenelopeBremsstrahlungAngular::ClearTab << 
 86 {                                                  68 {
 87   if (fLorentzTables1)                         <<  69   G4double pZ[NumberofZPoints] = {2.0,8.0,13.0,47.0,79.0,92.0};
 88     {                                          <<  70   G4double pX[NumberofZPoints],pY[NumberofZPoints];
 89       for (auto j=fLorentzTables1->begin(); j  <<  71   G4double QQ1[NumberofZPoints][NumberofEPoints][NumberofKPoints];
 90         {                                      <<  72   G4double QQ2[NumberofZPoints][NumberofEPoints][NumberofKPoints];
 91     G4PhysicsTable* tab = j->second;           <<  73  
 92           tab->clearAndDestroy();              <<  74   //Read information from DataBase file
 93           delete tab;                          <<  75   char* path = getenv("G4LEDATA");
 94         }                                      << 
 95       fLorentzTables1->clear();                << 
 96       delete fLorentzTables1;                  << 
 97       fLorentzTables1 = nullptr;               << 
 98     }                                          << 
 99                                                << 
100   if (fLorentzTables2)                         << 
101     {                                          << 
102       for (auto j=fLorentzTables2->begin(); j  << 
103         {                                      << 
104     G4PhysicsTable* tab = j->second;           << 
105     tab->clearAndDestroy();                    << 
106           delete tab;                          << 
107         }                                      << 
108       fLorentzTables2->clear();                << 
109       delete fLorentzTables2;                  << 
110       fLorentzTables2 = nullptr;               << 
111     }                                          << 
112   if (fEffectiveZSq)                           << 
113     {                                          << 
114       delete fEffectiveZSq;                    << 
115       fEffectiveZSq = nullptr;                 << 
116     }                                          << 
117 }                                              << 
118                                                << 
119 //....oooOO0OOooo........oooOO0OOooo........oo << 
120                                                << 
121 void G4PenelopeBremsstrahlungAngular::ReadData << 
122 {                                              << 
123    //Read information from DataBase file       << 
124   const char* path = G4FindDataDir("G4LEDATA") << 
125   if (!path)                                       76   if (!path)
126     {                                              77     {
127       G4String excep =                         <<  78       G4String excep = "G4PenelopeBremsstrahlungAngular - G4LEDATA environment variable not set!";
128   "G4PenelopeBremsstrahlungAngular - G4LEDATA  <<  79       G4Exception(excep);
129       G4Exception("G4PenelopeBremsstrahlungAng << 
130       "em0006",FatalException,excep);          << 
131       return;                                  << 
132     }                                              80     }
133   G4String pathString(path);                       81   G4String pathString(path);
134   G4String pathFile = pathString + "/penelope/ <<  82   G4String pathFile = pathString + "/penelope/br-ang-pen.dat";
135   std::ifstream file(pathFile);                    83   std::ifstream file(pathFile);
136                                                <<  84   std::filebuf* lsdp = file.rdbuf();
137   if (!file.is_open())                         <<  85   
                                                   >>  86   if (!(lsdp->is_open()))
138     {                                              87     {
139       G4String excep = "G4PenelopeBremsstrahlu     88       G4String excep = "G4PenelopeBremsstrahlungAngular - data file " + pathFile + " not found!";
140       G4Exception("G4PenelopeBremsstrahlungAng <<  89       G4Exception(excep);
141       "em0003",FatalException,excep);          << 
142       return;                                  << 
143     }                                              90     }
144   G4int i=0,j=0,k=0; // i=index for Z, j=index <<  91   G4int i=0,j=0,k=0; // i=index for Z, j=index for E, k=index for K 
145                                                <<  92   G4double a1,a2;
146   for (k=0;k<fNumberofKPoints;k++)             <<  93   while(i != -1) {
147     for (i=0;i<fNumberofZPoints;i++)           <<  94     file >> i >> j >> k >> a1 >> a2; 
148       for (j=0;j<fNumberofEPoints;j++)         <<  95     if (i > -1){
149   {                                            <<  96       QQ1[i][j][k]=a1;
150     G4double a1,a2;                            <<  97       QQ2[i][j][k]=a2;
151     G4int ik1,iz1,ie1;                         <<  98     }
152     G4double zr,er,kr;                         <<  99   } 
153     file >> iz1 >> ie1 >> ik1 >> zr >> er >> k << 
154     //check the data are correct               << 
155     if ((iz1-1 == i) && (ik1-1 == k) && (ie1-1 << 
156       {                                        << 
157         fQQ1[i][j][k]=a1;                      << 
158         fQQ2[i][j][k]=a2;                      << 
159       }                                        << 
160     else                                       << 
161       {                                        << 
162         G4ExceptionDescription ed;             << 
163         ed << "Corrupted data file " << pathFi << 
164         G4Exception("G4PenelopeBremsstrahlungA << 
165       "em0005",FatalException,ed);             << 
166       }                                        << 
167   }                                            << 
168   file.close();                                   100   file.close();
169   fDataRead = true;                            << 
170 }                                              << 
171                                                << 
172 //....oooOO0OOooo........oooOO0OOooo........oo << 
173                                                << 
174 void G4PenelopeBremsstrahlungAngular::PrepareT << 
175 {                                              << 
176   //Unused at the moment: the G4PenelopeBremss << 
177   //builds its own version of the tables.      << 
178                                                   101   
179   //Check if data file has already been read   << 
180   if (!fDataRead)                              << 
181     {                                          << 
182       ReadDataFile();                          << 
183       if (!fDataRead)                          << 
184   G4Exception("G4PenelopeBremsstrahlungAngular << 
185         "em2001",FatalException,"Unable to bui << 
186     }                                          << 
187                                                << 
188   if (!fLorentzTables1)                        << 
189       fLorentzTables1 = new std::map<G4double, << 
190   if (!fLorentzTables2)                        << 
191     fLorentzTables2 = new std::map<G4double,G4 << 
192                                                << 
193   G4double Zmat = CalculateEffectiveZ(material << 
194                                                << 
195   const G4int reducedEnergyGrid=21;            << 
196   //Support arrays.                            << 
197   G4double betas[fNumberofEPoints]; //betas fo << 
198   //tables for interpolation                   << 
199   G4double Q1[fNumberofEPoints][fNumberofKPoin << 
200   G4double Q2[fNumberofEPoints][fNumberofKPoin << 
201   //expanded tables for interpolation          << 
202   G4double Q1E[fNumberofEPoints][reducedEnergy << 
203   G4double Q2E[fNumberofEPoints][reducedEnergy << 
204   G4double pZ[fNumberofZPoints] = {2.0,8.0,13. << 
205                                                   102 
206   G4int i=0,j=0,k=0; // i=index for Z, j=index << 
207   //Interpolation in Z                            103   //Interpolation in Z
208   for (i=0;i<fNumberofEPoints;i++)             << 104   for (i=0;i<NumberofEPoints;i++){
209     {                                          << 105     for (j=0;j<NumberofKPoints;j++){
210       for (j=0;j<fNumberofKPoints;j++)         << 106       for (k=0;k<NumberofZPoints;k++){
211   {                                            << 107   pX[k]=std::log(QQ1[k][i][j]);
212     G4PhysicsFreeVector* QQ1vector =           << 108   pY[k]=QQ2[k][i][j];
213       new G4PhysicsFreeVector(fNumberofZPoints << 109       }
214     G4PhysicsFreeVector* QQ2vector =           << 110       G4PenelopeInterpolator* interpolator1 = new G4PenelopeInterpolator(pZ,pX,NumberofZPoints);
215       new G4PhysicsFreeVector(fNumberofZPoints << 111       Q1[i][j]=std::exp(interpolator1->CubicSplineInterpolation((G4double) Zmat));
216                                                << 112       delete interpolator1;
217     //fill vectors                             << 113       G4PenelopeInterpolator* interpolator2 = new G4PenelopeInterpolator(pZ,pY,NumberofZPoints);    
218     for (k=0;k<fNumberofZPoints;k++)           << 114       Q2[i][j]=interpolator2->CubicSplineInterpolation((G4double) Zmat);
219       {                                        << 115       delete interpolator2;
220         QQ1vector->PutValues(k,pZ[k],G4Log(fQQ << 
221         QQ2vector->PutValues(k,pZ[k],fQQ2[k][i << 
222       }                                        << 
223     //Filled: now calculate derivatives        << 
224     QQ1vector->FillSecondDerivatives();        << 
225     QQ2vector->FillSecondDerivatives();        << 
226                                                << 
227     Q1[i][j]= G4Exp(QQ1vector->Value(Zmat));   << 
228     Q2[i][j]=QQ2vector->Value(Zmat);           << 
229     delete QQ1vector;                          << 
230     delete QQ2vector;                          << 
231   }                                            << 
232     }                                             116     }
233   G4double pE[fNumberofEPoints] = {1.0e-03*MeV << 117   }
234           1.0e-01*MeV,5.0e-01*MeV};            << 118  
235   G4double pK[fNumberofKPoints] = {0.0,0.6,0.8 << 119   
                                                   >> 120   //std::ofstream fil("matrice.dat",std::ios::app);
                                                   >> 121   //fil << "Numero atomico: " << Zmat << G4endl;
                                                   >> 122   //for (i=0;i<NumberofEPoints;i++)
                                                   >> 123   //{
                                                   >> 124   //  fil << Q1[i][0] << " " << Q1[i][1] << " " << Q1[i][2] << " " << Q1[i][3] << G4endl;
                                                   >> 125   //}
                                                   >> 126   //fil.close();
                                                   >> 127   
                                                   >> 128 }
                                                   >> 129 
                                                   >> 130 void G4PenelopeBremsstrahlungAngular::InterpolationForK()
                                                   >> 131 {
                                                   >> 132   G4double pE[NumberofEPoints] = {1.0e-03,5.0e-03,1.0e-02,5.0e-02,1.0e-01,5.0e-01};
                                                   >> 133   G4double pK[NumberofKPoints] = {0.0,0.6,0.8,0.95};
236   G4double ppK[reducedEnergyGrid];                134   G4double ppK[reducedEnergyGrid];
                                                   >> 135   G4double pX[NumberofKPoints];
                                                   >> 136   G4int i,j;
237                                                   137 
238   for(i=0;i<reducedEnergyGrid;i++)             << 138   for(i=0;i<reducedEnergyGrid;i++){
239     ppK[i]=((G4double) i) * 0.05;                 139     ppK[i]=((G4double) i) * 0.05;
                                                   >> 140   }
240                                                   141 
241                                                << 142   for(i=0;i<NumberofEPoints;i++){
242   for(i=0;i<fNumberofEPoints;i++)              << 
243     betas[i]=std::sqrt(pE[i]*(pE[i]+2*electron    143     betas[i]=std::sqrt(pE[i]*(pE[i]+2*electron_mass_c2))/(pE[i]+electron_mass_c2);
                                                   >> 144   }
244                                                   145 
245                                                << 146   for (i=0;i<NumberofEPoints;i++){
246   for (i=0;i<fNumberofEPoints;i++)             << 147     for (j=0;j<NumberofKPoints;j++){
247     {                                          << 148       Q1[i][j]=Q1[i][j]/((G4double) Zmat);
248       for (j=0;j<fNumberofKPoints;j++)         << 
249   Q1[i][j]=Q1[i][j]/Zmat;                      << 
250     }                                             149     }
                                                   >> 150   }
251                                                   151 
252   //Expanded table of distribution parameters     152   //Expanded table of distribution parameters
253   for (i=0;i<fNumberofEPoints;i++)             << 153   for (i=0;i<NumberofEPoints;i++){
254     {                                          << 154     for (j=0;j<NumberofKPoints;j++){
255       G4PhysicsFreeVector* Q1vector = new G4Ph << 155       pX[j]=std::log(Q1[i][j]); //logarithmic 
256       G4PhysicsFreeVector* Q2vector = new G4Ph << 
257                                                << 
258       for (j=0;j<fNumberofKPoints;j++)         << 
259   {                                            << 
260     Q1vector->PutValues(j,pK[j],G4Log(Q1[i][j] << 
261     Q2vector->PutValues(j,pK[j],Q2[i][j]);     << 
262   }                                            << 
263                                                << 
264       for (j=0;j<reducedEnergyGrid;j++)        << 
265   {                                            << 
266     Q1E[i][j]=Q1vector->Value(ppK[j]);         << 
267     Q2E[i][j]=Q2vector->Value(ppK[j]);         << 
268   }                                            << 
269       delete Q1vector;                         << 
270       delete Q2vector;                         << 
271     }                                          << 
272   //                                           << 
273   //TABLES to be stored                        << 
274   //                                           << 
275   G4PhysicsTable* theTable1 = new G4PhysicsTab << 
276   G4PhysicsTable* theTable2 = new G4PhysicsTab << 
277   // the table will contain reducedEnergyGrid  << 
278   // values of k,                              << 
279   // Each of the G4PhysicsFreeVectors has a pr << 
280   // y vs. E                                   << 
281   //                                           << 
282   //reserve space of the vectors.              << 
283   for (j=0;j<reducedEnergyGrid;j++)            << 
284     {                                          << 
285       G4PhysicsFreeVector* thevec = new G4Phys << 
286       theTable1->push_back(thevec);            << 
287       G4PhysicsFreeVector* thevec2 = new G4Phy << 
288       theTable2->push_back(thevec2);           << 
289     }                                             156     }
290                                                << 157     G4PenelopeInterpolator* interpolator = new G4PenelopeInterpolator(pK,pX,NumberofKPoints);
291   for (j=0;j<reducedEnergyGrid;j++)            << 158     for (j=0;j<reducedEnergyGrid;j++){
292     {                                          << 159       Q1E[i][j]=interpolator->CubicSplineInterpolation(ppK[j]);
293       G4PhysicsFreeVector* thevec = (G4Physics << 
294       G4PhysicsFreeVector* thevec2 = (G4Physic << 
295       for (i=0;i<fNumberofEPoints;i++)         << 
296   {                                            << 
297     thevec->PutValues(i,betas[i],Q1E[i][j]);   << 
298     thevec2->PutValues(i,betas[i],Q2E[i][j]);  << 
299   }                                            << 
300       //Vectors are filled: calculate derivati << 
301       thevec->FillSecondDerivatives();         << 
302       thevec2->FillSecondDerivatives();        << 
303     }                                             160     }
304                                                << 161     delete interpolator;
305   if (fLorentzTables1 && fLorentzTables2)      << 162     for (j=0;j<NumberofKPoints;j++){
306     {                                          << 163       pX[j]=Q2[i][j];
307       fLorentzTables1->insert(std::make_pair(Z << 
308       fLorentzTables2->insert(std::make_pair(Z << 
309     }                                             164     }
310   else                                         << 165     G4PenelopeInterpolator* interpolator2 = new G4PenelopeInterpolator(pK,pX,NumberofKPoints);
311     {                                          << 166     for (j=0;j<reducedEnergyGrid;j++){
312       G4ExceptionDescription ed;               << 167       Q2E[i][j]=interpolator2->CubicSplineInterpolation(ppK[j]);
313       ed << "Unable to create tables of Lorent << 
314       ed << "<Z>= "  << Zmat << " in G4Penelop << 
315       delete theTable1;                        << 
316       delete theTable2;                        << 
317       G4Exception("G4PenelopeBremsstrahlungAng << 
318       "em2005",FatalException,ed);             << 
319     }                                             168     }
320                                                << 169     delete interpolator2;
321   return;                                      << 170   }
322 }                                                 171 }
323                                                   172 
324 //....oooOO0OOooo........oooOO0OOooo........oo << 173 G4double G4PenelopeBremsstrahlungAngular::ExtractCosTheta(G4double e1,G4double e2)
325                                                << 
326 G4ThreeVector& G4PenelopeBremsstrahlungAngular << 
327                 G4double eGamma,               << 
328                 G4int,                         << 
329                 const G4Material* material)    << 
330 {                                                 174 {
331   if (!material)                               << 175   //e1 = kinetic energy of the electron
332     {                                          << 176   //e2 = energy of the bremsstrahlung photon
333       G4Exception("G4PenelopeBremsstrahlungAng << 
334       "em2040",FatalException,"The pointer to  << 
335       return fLocalDirection;                  << 
336     }                                          << 
337                                                << 
338   //Retrieve the effective Z                   << 
339   G4double Zmat = 0;                           << 
340                                                << 
341   //The model might be initialized incorrectly << 
342   //G4PenelopeBremsstrahungModel: make sure it << 
343   if (!fEffectiveZSq)                          << 
344     {                                          << 
345       G4Exception("G4PenelopeBremsstrahlungAng << 
346       "em2040",JustWarning,"EffectiveZSq table << 
347       PrepareTables(material,false);           << 
348       //return fLocalDirection;                << 
349     }                                          << 
350                                                << 
351   //found in the table: return it              << 
352   if (fEffectiveZSq->count(material))          << 
353     Zmat = fEffectiveZSq->find(material)->seco << 
354   else //this can happen in unit tests or when << 
355     //models other than G4PenelopeBremsstrahun << 
356     {                                          << 
357       G4Exception("G4PenelopeBremsstrahlungAng << 
358       "em2040",JustWarning,"Material not found << 
359       PrepareTables(material,false);           << 
360       Zmat = fEffectiveZSq->find(material)->se << 
361       //      return fLocalDirection;          << 
362     }                                          << 
363                                                << 
364   if (fVerbosityLevel > 0)                     << 
365     {                                          << 
366       G4cout << "Effective <Z> for material :  << 
367   " = " << Zmat << G4endl;                     << 
368     }                                          << 
369                                                << 
370   G4double ePrimary = dp->GetKineticEnergy();  << 
371                                                << 
372   G4double beta = std::sqrt(ePrimary*(ePrimary << 
373     (ePrimary+electron_mass_c2);               << 
374   G4double cdt = 0;                            << 
375   G4double sinTheta = 0;                       << 
376   G4double phi = 0;                            << 
377                                                   177 
378   //Use a pure dipole distribution for energy  << 178   G4double beta = std::sqrt(e1*(e1+2*electron_mass_c2))/(e1+electron_mass_c2);
379   if (ePrimary > 500*keV)                      << 179   
380     {                                          << 
381       cdt = 2.0*G4UniformRand() - 1.0;         << 
382       if (G4UniformRand() > 0.75)              << 
383   {                                            << 
384     if (cdt<0)                                 << 
385       cdt = -1.0*std::pow(-cdt,1./3.);         << 
386     else                                       << 
387       cdt = std::pow(cdt,1./3.);               << 
388   }                                            << 
389       cdt = (cdt+beta)/(1.0+beta*cdt);         << 
390       //Get primary kinematics                 << 
391       sinTheta = std::sqrt(1. - cdt*cdt);      << 
392       phi  = twopi * G4UniformRand();          << 
393       fLocalDirection.set(sinTheta* std::cos(p << 
394           sinTheta* std::sin(phi),             << 
395           cdt);                                << 
396       //rotate                                 << 
397       fLocalDirection.rotateUz(dp->GetMomentum << 
398       //return                                 << 
399       return fLocalDirection;                  << 
400     }                                          << 
401                                                << 
402   if (!(fLorentzTables1->count(Zmat)) || !(fLo << 
403     {                                          << 
404       G4ExceptionDescription ed;               << 
405       ed << "Unable to retrieve Lorentz tables << 
406       G4Exception("G4PenelopeBremsstrahlungAng << 
407       "em2006",FatalException,ed);             << 
408     }                                          << 
409                                                   180 
410   //retrieve actual tables                     << 
411   const G4PhysicsTable* theTable1 = fLorentzTa << 
412   const G4PhysicsTable* theTable2 = fLorentzTa << 
413                                                   181 
414   G4double RK=20.0*eGamma/ePrimary;            << 182   G4double RK=20.0*e2/e1;
415   G4int ik=std::min((G4int) RK,19);               183   G4int ik=std::min((G4int) RK,19);
416                                                << 184   
417   G4double P10=0,P11=0,P1=0;                      185   G4double P10=0,P11=0,P1=0;
418   G4double P20=0,P21=0,P2=0;                      186   G4double P20=0,P21=0,P2=0;
419                                                << 187   G4double pX[NumberofEPoints];
420   //First coefficient                             188   //First coefficient
421   const G4PhysicsFreeVector* v1 = (G4PhysicsFr << 189   G4int i;
422   const G4PhysicsFreeVector* v2 = (G4PhysicsFr << 190   G4int j = ik;
423   P10 = v1->Value(beta);                       << 191   for (i=0;i<NumberofEPoints;i++){
424   P11 = v2->Value(beta);                       << 192     pX[i]=Q1E[i][j];
                                                   >> 193   }
                                                   >> 194   G4PenelopeInterpolator* interpolator = new G4PenelopeInterpolator(betas,pX,NumberofEPoints);
                                                   >> 195   P10=interpolator->CubicSplineInterpolation(beta);
                                                   >> 196   delete interpolator;
                                                   >> 197   j++; //(j=ik+1)
                                                   >> 198   for (i=0;i<NumberofEPoints;i++){
                                                   >> 199     pX[i]=Q1E[i][j];
                                                   >> 200   }
                                                   >> 201   G4PenelopeInterpolator* interpolator2 = new G4PenelopeInterpolator(betas,pX,NumberofEPoints);
                                                   >> 202   P11=interpolator2->CubicSplineInterpolation(beta);
                                                   >> 203   delete interpolator2;
425   P1=P10+(RK-(G4double) ik)*(P11-P10);            204   P1=P10+(RK-(G4double) ik)*(P11-P10);
426                                                << 205   
427   //Second coefficient                            206   //Second coefficient
428   const G4PhysicsFreeVector* v3 = (G4PhysicsFr << 207   j = ik;
429   const G4PhysicsFreeVector* v4 = (G4PhysicsFr << 208   for (i=0;i<NumberofEPoints;i++){
430   P20=v3->Value(beta);                         << 209     pX[i]=Q2E[i][j];
431   P21=v4->Value(beta);                         << 210   }
                                                   >> 211   G4PenelopeInterpolator* interpolator3 = new G4PenelopeInterpolator(betas,pX,NumberofEPoints);
                                                   >> 212   P20=interpolator3->CubicSplineInterpolation(beta);
                                                   >> 213   delete interpolator3;
                                                   >> 214   j++; //(j=ik+1)
                                                   >> 215   for (i=0;i<NumberofEPoints;i++){
                                                   >> 216     pX[i]=Q2E[i][j];
                                                   >> 217   }
                                                   >> 218   G4PenelopeInterpolator* interpolator4 = new G4PenelopeInterpolator(betas,pX,NumberofEPoints);
                                                   >> 219   P21=interpolator4->CubicSplineInterpolation(beta);
                                                   >> 220   delete interpolator4;
432   P2=P20+(RK-(G4double) ik)*(P21-P20);            221   P2=P20+(RK-(G4double) ik)*(P21-P20);
433                                                << 222   
434   //Sampling from the Lorenz-trasformed dipole    223   //Sampling from the Lorenz-trasformed dipole distributions
435   P1=std::min(G4Exp(P1)/beta,1.0);             << 224   P1=std::min(std::exp(P1)/beta,1.0);
436   G4double betap = std::min(std::max(beta*(1.0    225   G4double betap = std::min(std::max(beta*(1.0+P2/beta),0.0),0.9999);
437                                                << 226   
438   G4double testf=0;                            << 227   G4double cdt=0,testf=0;
439                                                << 228   
440   if (G4UniformRand() < P1)                    << 229   if (G4UniformRand() < P1){
441     {                                          << 230     do{
442       do{                                      << 231       cdt = 2.0*G4UniformRand()-1.0;
443   cdt = 2.0*G4UniformRand()-1.0;               << 232       testf=2.0*G4UniformRand()-(1.0+cdt*cdt);
444   testf=2.0*G4UniformRand()-(1.0+cdt*cdt);     << 233     }while(testf>0);
445       }while(testf>0);                         << 234   }
446     }                                          << 235   else{
447   else                                         << 236     do{
448     {                                          << 237       cdt = 2.0*G4UniformRand()-1.0;
449       do{                                      << 238       testf=G4UniformRand()-(1.0-cdt*cdt);
450   cdt = 2.0*G4UniformRand()-1.0;               << 239     }while(testf>0);
451   testf=G4UniformRand()-(1.0-cdt*cdt);         << 240   }
452       }while(testf>0);                         << 
453     }                                          << 
454   cdt = (cdt+betap)/(1.0+betap*cdt);              241   cdt = (cdt+betap)/(1.0+betap*cdt);
455                                                << 242   return cdt;
456   //Get primary kinematics                     << 
457   sinTheta = std::sqrt(1. - cdt*cdt);          << 
458   phi  = twopi * G4UniformRand();              << 
459   fLocalDirection.set(sinTheta* std::cos(phi), << 
460           sinTheta* std::sin(phi),             << 
461           cdt);                                << 
462   //rotate                                     << 
463   fLocalDirection.rotateUz(dp->GetMomentumDire << 
464   //return                                     << 
465   return fLocalDirection;                      << 
466                                                << 
467 }                                              << 
468                                                << 
469 //....oooOO0OOooo........oooOO0OOooo........oo << 
470                                                << 
471 G4double G4PenelopeBremsstrahlungAngular::Calc << 
472 {                                              << 
473   if (!fEffectiveZSq)                          << 
474     fEffectiveZSq = new std::map<const G4Mater << 
475                                                << 
476   //Already exists: return it                  << 
477   if (fEffectiveZSq->count(material))          << 
478     return fEffectiveZSq->find(material)->seco << 
479                                                << 
480   //Helper for the calculation                 << 
481   std::vector<G4double> *StechiometricFactors  << 
482   G4int nElements = (G4int)material->GetNumber << 
483   const G4ElementVector* elementVector = mater << 
484   const G4double* fractionVector = material->G << 
485   for (G4int i=0;i<nElements;++i)              << 
486     {                                          << 
487       G4double fraction = fractionVector[i];   << 
488       G4double atomicWeigth = (*elementVector) << 
489       StechiometricFactors->push_back(fraction << 
490     }                                          << 
491   //Find max                                   << 
492   G4double MaxStechiometricFactor = 0.;        << 
493   for (G4int i=0;i<nElements;++i)              << 
494     {                                          << 
495       if ((*StechiometricFactors)[i] > MaxStec << 
496         MaxStechiometricFactor = (*Stechiometr << 
497     }                                          << 
498   //Normalize                                  << 
499   for (G4int i=0;i<nElements;++i)              << 
500     (*StechiometricFactors)[i] /=  MaxStechiom << 
501                                                << 
502   G4double sumz2 = 0;                          << 
503   G4double sums = 0;                           << 
504   for (G4int i=0;i<nElements;++i)              << 
505     {                                          << 
506       G4double Z = (*elementVector)[i]->GetZ() << 
507       sumz2 += (*StechiometricFactors)[i]*Z*Z; << 
508       sums  += (*StechiometricFactors)[i];     << 
509     }                                          << 
510   delete StechiometricFactors;                 << 
511                                                << 
512   G4double ZBR = std::sqrt(sumz2/sums);        << 
513   fEffectiveZSq->insert(std::make_pair(materia << 
514                                                << 
515   return ZBR;                                  << 
516 }                                                 243 }
517                                                   244