Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/eRosita/physics/src/G4RDVeLowEnergyLoss.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 /examples/advanced/eRosita/physics/src/G4RDVeLowEnergyLoss.cc (Version 11.3.0) and /examples/advanced/eRosita/physics/src/G4RDVeLowEnergyLoss.cc (Version 10.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 //
                                                   >>  27 // $Id$
                                                   >>  28 // GEANT4 tag $Name: geant4-09-01-ref-00 $
 27 //                                                 29 //
 28 //                                                 30 //
 29 // -------------------------------------------     31 // --------------------------------------------------------------
 30 //  GEANT 4 class implementation file              32 //  GEANT 4 class implementation file
 31 //                                                 33 //
 32 //  History: first implementation, based on ob     34 //  History: first implementation, based on object model of
 33 //  2nd December 1995, G.Cosmo                     35 //  2nd December 1995, G.Cosmo
 34 // -------------------------------------------     36 // --------------------------------------------------------------
 35 //                                                 37 //
 36 // Modifications:                                  38 // Modifications:
 37 // 20/09/00 update fluctuations V.Ivanchenko       39 // 20/09/00 update fluctuations V.Ivanchenko
 38 // 22/11/00 minor fix in fluctuations V.Ivanch     40 // 22/11/00 minor fix in fluctuations V.Ivanchenko
 39 // 10/05/01  V.Ivanchenko Clean up againist Li     41 // 10/05/01  V.Ivanchenko Clean up againist Linux compilation with -Wall
 40 // 22/05/01  V.Ivanchenko Update range calcula     42 // 22/05/01  V.Ivanchenko Update range calculation
 41 // 23/11/01  V.Ivanchenko Move static member-f     43 // 23/11/01  V.Ivanchenko Move static member-functions from header to source
 42 // 22/01/03  V.Ivanchenko Cut per region           44 // 22/01/03  V.Ivanchenko Cut per region
 43 // 11/02/03  V.Ivanchenko Add limits to fluctu     45 // 11/02/03  V.Ivanchenko Add limits to fluctuations
 44 // 24/04/03  V.Ivanchenko Fix the problem of t     46 // 24/04/03  V.Ivanchenko Fix the problem of table size
 45 //                                                 47 //
 46 // -------------------------------------------     48 // --------------------------------------------------------------
 47                                                    49 
 48 #include "G4RDVeLowEnergyLoss.hh"                  50 #include "G4RDVeLowEnergyLoss.hh"
 49 #include "G4PhysicalConstants.hh"                  51 #include "G4PhysicalConstants.hh"
 50 #include "G4SystemOfUnits.hh"                      52 #include "G4SystemOfUnits.hh"
 51 #include "G4ProductionCutsTable.hh"                53 #include "G4ProductionCutsTable.hh"
 52                                                    54 
 53 G4double     G4RDVeLowEnergyLoss::ParticleMass     55 G4double     G4RDVeLowEnergyLoss::ParticleMass ;
 54 G4double     G4RDVeLowEnergyLoss::taulow           56 G4double     G4RDVeLowEnergyLoss::taulow       ;
 55 G4double     G4RDVeLowEnergyLoss::tauhigh          57 G4double     G4RDVeLowEnergyLoss::tauhigh      ;
 56 G4double     G4RDVeLowEnergyLoss::ltaulow          58 G4double     G4RDVeLowEnergyLoss::ltaulow       ;
 57 G4double     G4RDVeLowEnergyLoss::ltauhigh         59 G4double     G4RDVeLowEnergyLoss::ltauhigh      ;
 58                                                    60 
 59                                                    61 
 60 G4bool       G4RDVeLowEnergyLoss::rndmStepFlag     62 G4bool       G4RDVeLowEnergyLoss::rndmStepFlag   = false;
 61 G4bool       G4RDVeLowEnergyLoss::EnlossFlucFl     63 G4bool       G4RDVeLowEnergyLoss::EnlossFlucFlag = true;
 62 G4double     G4RDVeLowEnergyLoss::dRoverRange      64 G4double     G4RDVeLowEnergyLoss::dRoverRange    = 20*perCent;
 63 G4double     G4RDVeLowEnergyLoss::finalRange       65 G4double     G4RDVeLowEnergyLoss::finalRange     = 200*micrometer;
 64 G4double     G4RDVeLowEnergyLoss::c1lim = dRov     66 G4double     G4RDVeLowEnergyLoss::c1lim = dRoverRange ;
 65 G4double     G4RDVeLowEnergyLoss::c2lim = 2.*(     67 G4double     G4RDVeLowEnergyLoss::c2lim = 2.*(1.-dRoverRange)*finalRange ;
 66 G4double     G4RDVeLowEnergyLoss::c3lim = -(1.     68 G4double     G4RDVeLowEnergyLoss::c3lim = -(1.-dRoverRange)*finalRange*finalRange;
 67                                                    69 
 68                                                    70 
 69 //                                                 71 //
 70                                                    72 
 71 G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss()         73 G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss()
 72                    :G4VContinuousDiscreteProce     74                    :G4VContinuousDiscreteProcess("No Name Loss Process"),
 73      lastMaterial(0),                              75      lastMaterial(0),
 74      nmaxCont1(4),                                 76      nmaxCont1(4),
 75      nmaxCont2(16)                                 77      nmaxCont2(16)
 76 {                                                  78 {
 77   G4Exception("G4RDVeLowEnergyLoss::G4RDVeLowE     79   G4Exception("G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss()", "InvalidCall",
 78               FatalException, "Default constru     80               FatalException, "Default constructor is called!");
 79 }                                                  81 }
 80                                                    82 
 81 //                                                 83 //
 82                                                    84 
 83 G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss(const     85 G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss(const G4String& aName, G4ProcessType aType)
 84                   : G4VContinuousDiscreteProce     86                   : G4VContinuousDiscreteProcess(aName, aType),
 85      lastMaterial(0),                              87      lastMaterial(0),
 86      nmaxCont1(4),                                 88      nmaxCont1(4),
 87      nmaxCont2(16)                                 89      nmaxCont2(16)
 88 {                                                  90 {
 89 }                                                  91 }
 90                                                    92 
 91 //                                                 93 //
 92                                                    94 
 93 G4RDVeLowEnergyLoss::~G4RDVeLowEnergyLoss()        95 G4RDVeLowEnergyLoss::~G4RDVeLowEnergyLoss()
 94 {                                                  96 {
 95 }                                                  97 }
 96                                                    98 
 97 //                                                 99 //
 98                                                   100 
 99 G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss(G4RDV    101 G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss(G4RDVeLowEnergyLoss& right)
100                   : G4VContinuousDiscreteProce    102                   : G4VContinuousDiscreteProcess(right),
101      lastMaterial(0),                             103      lastMaterial(0),
102      nmaxCont1(4),                                104      nmaxCont1(4),
103      nmaxCont2(16)                                105      nmaxCont2(16)
104 {                                                 106 {
105 }                                                 107 }
106                                                   108 
107 void G4RDVeLowEnergyLoss::SetRndmStep(G4bool v    109 void G4RDVeLowEnergyLoss::SetRndmStep(G4bool value)
108 {                                                 110 {
109    rndmStepFlag   = value;                        111    rndmStepFlag   = value;
110 }                                                 112 }
111                                                   113 
112 void G4RDVeLowEnergyLoss::SetEnlossFluc(G4bool    114 void G4RDVeLowEnergyLoss::SetEnlossFluc(G4bool value)
113 {                                                 115 {
114    EnlossFlucFlag = value;                        116    EnlossFlucFlag = value;
115 }                                                 117 }
116                                                   118 
117 void G4RDVeLowEnergyLoss::SetStepFunction (G4d    119 void G4RDVeLowEnergyLoss::SetStepFunction (G4double c1, G4double c2)
118 {                                                 120 {
119    dRoverRange = c1;                              121    dRoverRange = c1;
120    finalRange = c2;                               122    finalRange = c2;
121    c1lim=dRoverRange;                             123    c1lim=dRoverRange;
122    c2lim=2.*(1-dRoverRange)*finalRange;           124    c2lim=2.*(1-dRoverRange)*finalRange;
123    c3lim=-(1.-dRoverRange)*finalRange*finalRan    125    c3lim=-(1.-dRoverRange)*finalRange*finalRange;
124 }                                                 126 }
125                                                   127 
126 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRang    128 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRangeTable(G4PhysicsTable* theDEDXTable,
127                G4PhysicsTable* theRangeTable,     129                G4PhysicsTable* theRangeTable,
128                G4double lowestKineticEnergy,      130                G4double lowestKineticEnergy,
129                G4double highestKineticEnergy,     131                G4double highestKineticEnergy,
130                G4int TotBin)                      132                G4int TotBin)
131 // Build range table from the energy loss tabl    133 // Build range table from the energy loss table
132 {                                                 134 {
133                                                   135 
134    G4int numOfCouples = theDEDXTable->length()    136    G4int numOfCouples = theDEDXTable->length();
135                                                   137 
136    if(theRangeTable)                              138    if(theRangeTable)
137    { theRangeTable->clearAndDestroy();            139    { theRangeTable->clearAndDestroy();
138      delete theRangeTable; }                      140      delete theRangeTable; }
139    theRangeTable = new G4PhysicsTable(numOfCou    141    theRangeTable = new G4PhysicsTable(numOfCouples);
140                                                   142 
141    // loop for materials                          143    // loop for materials
142                                                   144 
143    for (G4int J=0;  J<numOfCouples; J++)          145    for (G4int J=0;  J<numOfCouples; J++)
144    {                                              146    {
145      G4PhysicsLogVector* aVector;                 147      G4PhysicsLogVector* aVector;
146      aVector = new G4PhysicsLogVector(lowestKi    148      aVector = new G4PhysicsLogVector(lowestKineticEnergy,
147                               highestKineticEn    149                               highestKineticEnergy,TotBin);
148      BuildRangeVector(theDEDXTable,lowestKinet    150      BuildRangeVector(theDEDXTable,lowestKineticEnergy,highestKineticEnergy,
149                       TotBin,J,aVector);          151                       TotBin,J,aVector);
150      theRangeTable->insert(aVector);              152      theRangeTable->insert(aVector);
151    }                                              153    }
152    return theRangeTable ;                         154    return theRangeTable ;
153 }                                                 155 }
154                                                   156 
155 //                                                157 //
156                                                   158 
157 void G4RDVeLowEnergyLoss::BuildRangeVector(G4P    159 void G4RDVeLowEnergyLoss::BuildRangeVector(G4PhysicsTable* theDEDXTable,
158                                          G4dou    160                                          G4double lowestKineticEnergy,
159                                          G4dou    161                                          G4double,
160                                          G4int    162                                          G4int TotBin,
161                                          G4int    163                                          G4int materialIndex,
162                                          G4Phy    164                                          G4PhysicsLogVector* rangeVector)
163                                                   165 
164 //  create range vector for a material            166 //  create range vector for a material
165 {                                                 167 {
166   G4bool isOut;                                   168   G4bool isOut;
167   G4PhysicsVector* physicsVector= (*theDEDXTab    169   G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex];
168   G4double energy1 = lowestKineticEnergy;         170   G4double energy1 = lowestKineticEnergy;
169   G4double dedx    = physicsVector->GetValue(e    171   G4double dedx    = physicsVector->GetValue(energy1,isOut);
170   G4double range   = 0.5*energy1/dedx;            172   G4double range   = 0.5*energy1/dedx;
171   rangeVector->PutValue(0,range);                 173   rangeVector->PutValue(0,range);
172   G4int n = 100;                                  174   G4int n = 100;
173   G4double del = 1.0/(G4double)n ;                175   G4double del = 1.0/(G4double)n ;
174                                                   176 
175   for (G4int j=1; j<TotBin; j++) {                177   for (G4int j=1; j<TotBin; j++) {
176                                                   178 
177     G4double energy2 = rangeVector->GetLowEdge    179     G4double energy2 = rangeVector->GetLowEdgeEnergy(j);
178     G4double de = (energy2 - energy1) * del ;     180     G4double de = (energy2 - energy1) * del ;
179     G4double dedx1 = dedx ;                       181     G4double dedx1 = dedx ;
180                                                   182 
181     for (G4int i=1; i<n; i++) {                   183     for (G4int i=1; i<n; i++) {
182       G4double energy = energy1 + i*de ;          184       G4double energy = energy1 + i*de ;
183       G4double dedx2  = physicsVector->GetValu    185       G4double dedx2  = physicsVector->GetValue(energy,isOut);
184       range  += 0.5*de*(1.0/dedx1 + 1.0/dedx2)    186       range  += 0.5*de*(1.0/dedx1 + 1.0/dedx2);
185       dedx1   = dedx2;                            187       dedx1   = dedx2;
186     }                                             188     }
187     rangeVector->PutValue(j,range);               189     rangeVector->PutValue(j,range);
188     dedx = dedx1 ;                                190     dedx = dedx1 ;
189     energy1 = energy2 ;                           191     energy1 = energy2 ;
190   }                                               192   }
191 }                                                 193 }
192                                                   194 
193 //                                                195 //
194                                                   196 
195 G4double G4RDVeLowEnergyLoss::RangeIntLin(G4Ph    197 G4double G4RDVeLowEnergyLoss::RangeIntLin(G4PhysicsVector* physicsVector,
196           G4int nbin)                             198           G4int nbin)
197 //  num. integration, linear binning              199 //  num. integration, linear binning
198 {                                                 200 {
199   G4double dtau,Value,taui,ti,lossi,ci;           201   G4double dtau,Value,taui,ti,lossi,ci;
200   G4bool isOut;                                   202   G4bool isOut;
201   dtau = (tauhigh-taulow)/nbin;                   203   dtau = (tauhigh-taulow)/nbin;
202   Value = 0.;                                     204   Value = 0.;
203                                                   205 
204   for (G4int i=0; i<=nbin; i++)                   206   for (G4int i=0; i<=nbin; i++)
205   {                                               207   {
206     taui = taulow + dtau*i ;                      208     taui = taulow + dtau*i ;
207     ti = ParticleMass*taui;                       209     ti = ParticleMass*taui;
208     lossi = physicsVector->GetValue(ti,isOut);    210     lossi = physicsVector->GetValue(ti,isOut);
209     if(i==0)                                      211     if(i==0)
210       ci=0.5;                                     212       ci=0.5;
211     else                                          213     else
212     {                                             214     {
213       if(i<nbin)                                  215       if(i<nbin)
214         ci=1.;                                    216         ci=1.;
215       else                                        217       else
216         ci=0.5;                                   218         ci=0.5;
217     }                                             219     }
218     Value += ci/lossi;                            220     Value += ci/lossi;
219   }                                               221   }
220   Value *= ParticleMass*dtau;                     222   Value *= ParticleMass*dtau;
221   return Value;                                   223   return Value;
222 }                                                 224 }
223                                                   225 
224 //                                                226 //
225                                                   227 
226 G4double G4RDVeLowEnergyLoss::RangeIntLog(G4Ph    228 G4double G4RDVeLowEnergyLoss::RangeIntLog(G4PhysicsVector* physicsVector,
227           G4int nbin)                             229           G4int nbin)
228 //  num. integration, logarithmic binning         230 //  num. integration, logarithmic binning
229 {                                                 231 {
230   G4double ltt,dltau,Value,ui,taui,ti,lossi,ci    232   G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
231   G4bool isOut;                                   233   G4bool isOut;
232   ltt = ltauhigh-ltaulow;                         234   ltt = ltauhigh-ltaulow;
233   dltau = ltt/nbin;                               235   dltau = ltt/nbin;
234   Value = 0.;                                     236   Value = 0.;
235                                                   237 
236   for (G4int i=0; i<=nbin; i++)                   238   for (G4int i=0; i<=nbin; i++)
237   {                                               239   {
238     ui = ltaulow+dltau*i;                         240     ui = ltaulow+dltau*i;
239     taui = std::exp(ui);                          241     taui = std::exp(ui);
240     ti = ParticleMass*taui;                       242     ti = ParticleMass*taui;
241     lossi = physicsVector->GetValue(ti,isOut);    243     lossi = physicsVector->GetValue(ti,isOut);
242     if(i==0)                                      244     if(i==0)
243       ci=0.5;                                     245       ci=0.5;
244     else                                          246     else
245     {                                             247     {
246       if(i<nbin)                                  248       if(i<nbin)
247         ci=1.;                                    249         ci=1.;
248       else                                        250       else
249         ci=0.5;                                   251         ci=0.5;
250     }                                             252     }
251     Value += ci*taui/lossi;                       253     Value += ci*taui/lossi;
252   }                                               254   }
253   Value *= ParticleMass*dltau;                    255   Value *= ParticleMass*dltau;
254   return Value;                                   256   return Value;
255 }                                                 257 }
256                                                   258 
257                                                   259 
258 //                                                260 //
259                                                   261 
260 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildLabT    262 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildLabTimeTable(G4PhysicsTable* theDEDXTable,
261                  G4PhysicsTable* theLabTimeTab    263                  G4PhysicsTable* theLabTimeTable,
262                  G4double lowestKineticEnergy,    264                  G4double lowestKineticEnergy,
263                  G4double highestKineticEnergy    265                  G4double highestKineticEnergy,G4int TotBin)
264                                                   266 
265 {                                                 267 {
266                                                   268 
267   G4int numOfCouples = G4ProductionCutsTable::    269   G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
268                                                   270 
269   if(theLabTimeTable)                             271   if(theLabTimeTable)
270   { theLabTimeTable->clearAndDestroy();           272   { theLabTimeTable->clearAndDestroy();
271     delete theLabTimeTable; }                     273     delete theLabTimeTable; }
272   theLabTimeTable = new G4PhysicsTable(numOfCo    274   theLabTimeTable = new G4PhysicsTable(numOfCouples);
273                                                   275 
274                                                   276 
275   for (G4int J=0;  J<numOfCouples; J++)           277   for (G4int J=0;  J<numOfCouples; J++)
276   {                                               278   {
277     G4PhysicsLogVector* aVector;                  279     G4PhysicsLogVector* aVector;
278                                                   280 
279     aVector = new G4PhysicsLogVector(lowestKin    281     aVector = new G4PhysicsLogVector(lowestKineticEnergy,
280                             highestKineticEner    282                             highestKineticEnergy,TotBin);
281                                                   283 
282     BuildLabTimeVector(theDEDXTable,              284     BuildLabTimeVector(theDEDXTable,
283               lowestKineticEnergy,highestKinet    285               lowestKineticEnergy,highestKineticEnergy,TotBin,J,aVector);
284     theLabTimeTable->insert(aVector);             286     theLabTimeTable->insert(aVector);
285                                                   287 
286                                                   288 
287   }                                               289   }
288   return theLabTimeTable ;                        290   return theLabTimeTable ;
289 }                                                 291 }
290                                                   292 
291 //                                                293 //    
292                                                   294 
293 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildProp    295 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildProperTimeTable(G4PhysicsTable* theDEDXTable,
294               G4PhysicsTable* theProperTimeTab    296               G4PhysicsTable* theProperTimeTable,
295               G4double lowestKineticEnergy,       297               G4double lowestKineticEnergy,
296               G4double highestKineticEnergy,G4    298               G4double highestKineticEnergy,G4int TotBin)
297                                                   299                             
298 {                                                 300 {
299                                                   301 
300   G4int numOfCouples = G4ProductionCutsTable::    302   G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
301                                                   303  
302   if(theProperTimeTable)                          304   if(theProperTimeTable)
303   { theProperTimeTable->clearAndDestroy();        305   { theProperTimeTable->clearAndDestroy();
304     delete theProperTimeTable; }                  306     delete theProperTimeTable; }
305   theProperTimeTable = new G4PhysicsTable(numO    307   theProperTimeTable = new G4PhysicsTable(numOfCouples);
306                                                   308 
307                                                   309 
308   for (G4int J=0;  J<numOfCouples; J++)           310   for (G4int J=0;  J<numOfCouples; J++)
309   {                                               311   {
310     G4PhysicsLogVector* aVector;                  312     G4PhysicsLogVector* aVector;
311                                                   313 
312     aVector = new G4PhysicsLogVector(lowestKin    314     aVector = new G4PhysicsLogVector(lowestKineticEnergy,
313                             highestKineticEner    315                             highestKineticEnergy,TotBin);
314                                                   316 
315     BuildProperTimeVector(theDEDXTable,           317     BuildProperTimeVector(theDEDXTable,
316               lowestKineticEnergy,highestKinet    318               lowestKineticEnergy,highestKineticEnergy,TotBin,J,aVector);
317     theProperTimeTable->insert(aVector);          319     theProperTimeTable->insert(aVector);
318                                                   320 
319                                                   321 
320   }                                               322   }
321   return theProperTimeTable ;                     323   return theProperTimeTable ;
322 }                                                 324 }
323                                                   325 
324 //                                                326 //    
325                                                   327 
326 void G4RDVeLowEnergyLoss::BuildLabTimeVector(G    328 void G4RDVeLowEnergyLoss::BuildLabTimeVector(G4PhysicsTable* theDEDXTable,
327              G4double, // lowestKineticEnergy,    329              G4double, // lowestKineticEnergy,
328              G4double, // highestKineticEnergy    330              G4double, // highestKineticEnergy,
329                                            G4i    331                                            G4int TotBin,
330              G4int materialIndex, G4PhysicsLog    332              G4int materialIndex, G4PhysicsLogVector* timeVector)
331 //  create lab time vector for a material         333 //  create lab time vector for a material
332 {                                                 334 {
333                                                   335 
334   G4int nbin=100;                                 336   G4int nbin=100;
335   G4bool isOut;                                   337   G4bool isOut;
336   G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-p    338   G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-parlowen ;
337   G4double losslim,clim,taulim,timelim,           339   G4double losslim,clim,taulim,timelim,
338            LowEdgeEnergy,tau,Value ;              340            LowEdgeEnergy,tau,Value ;
339                                                   341 
340   G4PhysicsVector* physicsVector= (*theDEDXTab    342   G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex];
341                                                   343 
342   // low energy part first...                     344   // low energy part first...
343   losslim = physicsVector->GetValue(tlim,isOut    345   losslim = physicsVector->GetValue(tlim,isOut);
344   taulim=tlim/ParticleMass ;                      346   taulim=tlim/ParticleMass ;
345   clim=std::sqrt(ParticleMass*tlim/2.)/(c_ligh    347   clim=std::sqrt(ParticleMass*tlim/2.)/(c_light*losslim*ppar) ;
346                                                   348 
347   G4int i=-1;                                     349   G4int i=-1;
348   G4double oldValue = 0. ;                        350   G4double oldValue = 0. ;
349   G4double tauold ;                               351   G4double tauold ;
350   do                                              352   do
351   {                                               353   {
352     i += 1 ;                                      354     i += 1 ;
353     LowEdgeEnergy = timeVector->GetLowEdgeEner    355     LowEdgeEnergy = timeVector->GetLowEdgeEnergy(i);
354     tau = LowEdgeEnergy/ParticleMass ;            356     tau = LowEdgeEnergy/ParticleMass ;
355     if ( tau <= taulim )                          357     if ( tau <= taulim )
356     {                                             358     {
357       Value = clim*std::exp(ppar*std::log(tau/    359       Value = clim*std::exp(ppar*std::log(tau/taulim)) ;
358     }                                             360     }
359     else                                          361     else
360     {                                             362     {
361       timelim=clim ;                              363       timelim=clim ;
362       ltaulow = std::log(taulim);                 364       ltaulow = std::log(taulim);
363       ltauhigh = std::log(tau);                   365       ltauhigh = std::log(tau);
364       Value = timelim+LabTimeIntLog(physicsVec    366       Value = timelim+LabTimeIntLog(physicsVector,nbin);
365     }                                             367     }
366     timeVector->PutValue(i,Value);                368     timeVector->PutValue(i,Value);
367     oldValue = Value ;                            369     oldValue = Value ;
368     tauold = tau ;                                370     tauold = tau ;
369   } while (tau<=taulim) ;                         371   } while (tau<=taulim) ;
370   i += 1 ;                                        372   i += 1 ;
371   for (G4int j=i; j<TotBin; j++)                  373   for (G4int j=i; j<TotBin; j++)
372   {                                               374   {
373     LowEdgeEnergy = timeVector->GetLowEdgeEner    375     LowEdgeEnergy = timeVector->GetLowEdgeEnergy(j);
374     tau = LowEdgeEnergy/ParticleMass ;            376     tau = LowEdgeEnergy/ParticleMass ;
375     ltaulow = std::log(tauold);                   377     ltaulow = std::log(tauold);
376     ltauhigh = std::log(tau);                     378     ltauhigh = std::log(tau);
377     Value = oldValue+LabTimeIntLog(physicsVect    379     Value = oldValue+LabTimeIntLog(physicsVector,nbin);
378     timeVector->PutValue(j,Value);                380     timeVector->PutValue(j,Value);
379     oldValue = Value ;                            381     oldValue = Value ;
380     tauold = tau ;                                382     tauold = tau ;
381   }                                               383   }
382 }                                                 384 }
383                                                   385 
384 //                                                386 //    
385                                                   387 
386 void G4RDVeLowEnergyLoss::BuildProperTimeVecto    388 void G4RDVeLowEnergyLoss::BuildProperTimeVector(G4PhysicsTable* theDEDXTable,
387                 G4double, // lowestKineticEner    389                 G4double, // lowestKineticEnergy,
388                 G4double, // highestKineticEne    390                 G4double, // highestKineticEnergy,
389                                                   391                                               G4int TotBin,
390                 G4int materialIndex, G4Physics    392                 G4int materialIndex, G4PhysicsLogVector* timeVector)
391 //  create proper time vector for a material      393 //  create proper time vector for a material
392 {                                                 394 {
393   G4int nbin=100;                                 395   G4int nbin=100;
394   G4bool isOut;                                   396   G4bool isOut;
395   G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-p    397   G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-parlowen ;
396   G4double losslim,clim,taulim,timelim,           398   G4double losslim,clim,taulim,timelim,
397            LowEdgeEnergy,tau,Value ;              399            LowEdgeEnergy,tau,Value ;
398                                                   400 
399   G4PhysicsVector* physicsVector= (*theDEDXTab    401   G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex];
400   //const G4MaterialTable* theMaterialTable =     402   //const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
401                                                   403 
402   // low energy part first...                     404   // low energy part first...
403   losslim = physicsVector->GetValue(tlim,isOut    405   losslim = physicsVector->GetValue(tlim,isOut);
404   taulim=tlim/ParticleMass ;                      406   taulim=tlim/ParticleMass ;
405   clim=std::sqrt(ParticleMass*tlim/2.)/(c_ligh    407   clim=std::sqrt(ParticleMass*tlim/2.)/(c_light*losslim*ppar) ;
406                                                   408 
407   G4int i=-1;                                     409   G4int i=-1;
408   G4double oldValue = 0. ;                        410   G4double oldValue = 0. ;
409   G4double tauold ;                               411   G4double tauold ;
410   do                                              412   do
411   {                                               413   {
412     i += 1 ;                                      414     i += 1 ;
413     LowEdgeEnergy = timeVector->GetLowEdgeEner    415     LowEdgeEnergy = timeVector->GetLowEdgeEnergy(i);
414     tau = LowEdgeEnergy/ParticleMass ;            416     tau = LowEdgeEnergy/ParticleMass ;
415     if ( tau <= taulim )                          417     if ( tau <= taulim )
416     {                                             418     {
417       Value = clim*std::exp(ppar*std::log(tau/    419       Value = clim*std::exp(ppar*std::log(tau/taulim)) ;
418     }                                             420     }
419     else                                          421     else
420     {                                             422     {
421       timelim=clim ;                              423       timelim=clim ;
422       ltaulow = std::log(taulim);                 424       ltaulow = std::log(taulim);
423       ltauhigh = std::log(tau);                   425       ltauhigh = std::log(tau);
424       Value = timelim+ProperTimeIntLog(physics    426       Value = timelim+ProperTimeIntLog(physicsVector,nbin);
425     }                                             427     }
426     timeVector->PutValue(i,Value);                428     timeVector->PutValue(i,Value);
427     oldValue = Value ;                            429     oldValue = Value ;
428     tauold = tau ;                                430     tauold = tau ;
429   } while (tau<=taulim) ;                         431   } while (tau<=taulim) ;
430   i += 1 ;                                        432   i += 1 ;
431   for (G4int j=i; j<TotBin; j++)                  433   for (G4int j=i; j<TotBin; j++)
432   {                                               434   {
433     LowEdgeEnergy = timeVector->GetLowEdgeEner    435     LowEdgeEnergy = timeVector->GetLowEdgeEnergy(j);
434     tau = LowEdgeEnergy/ParticleMass ;            436     tau = LowEdgeEnergy/ParticleMass ;
435     ltaulow = std::log(tauold);                   437     ltaulow = std::log(tauold);
436     ltauhigh = std::log(tau);                     438     ltauhigh = std::log(tau);
437     Value = oldValue+ProperTimeIntLog(physicsV    439     Value = oldValue+ProperTimeIntLog(physicsVector,nbin);
438     timeVector->PutValue(j,Value);                440     timeVector->PutValue(j,Value);
439     oldValue = Value ;                            441     oldValue = Value ;
440     tauold = tau ;                                442     tauold = tau ;
441   }                                               443   }
442 }                                                 444 }
443                                                   445 
444 //                                                446 //    
445                                                   447 
446 G4double G4RDVeLowEnergyLoss::LabTimeIntLog(G4    448 G4double G4RDVeLowEnergyLoss::LabTimeIntLog(G4PhysicsVector* physicsVector,
447             G4int nbin)                           449             G4int nbin)
448 //  num. integration, logarithmic binning         450 //  num. integration, logarithmic binning
449 {                                                 451 {
450   G4double ltt,dltau,Value,ui,taui,ti,lossi,ci    452   G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
451   G4bool isOut;                                   453   G4bool isOut;
452   ltt = ltauhigh-ltaulow;                         454   ltt = ltauhigh-ltaulow;
453   dltau = ltt/nbin;                               455   dltau = ltt/nbin;
454   Value = 0.;                                     456   Value = 0.;
455                                                   457 
456   for (G4int i=0; i<=nbin; i++)                   458   for (G4int i=0; i<=nbin; i++)
457   {                                               459   {
458     ui = ltaulow+dltau*i;                         460     ui = ltaulow+dltau*i;
459     taui = std::exp(ui);                          461     taui = std::exp(ui);
460     ti = ParticleMass*taui;                       462     ti = ParticleMass*taui;
461     lossi = physicsVector->GetValue(ti,isOut);    463     lossi = physicsVector->GetValue(ti,isOut);
462     if(i==0)                                      464     if(i==0)
463       ci=0.5;                                     465       ci=0.5;
464     else                                          466     else
465     {                                             467     {
466       if(i<nbin)                                  468       if(i<nbin)
467         ci=1.;                                    469         ci=1.;
468       else                                        470       else
469         ci=0.5;                                   471         ci=0.5;
470     }                                             472     }
471     Value += ci*taui*(ti+ParticleMass)/(std::s    473     Value += ci*taui*(ti+ParticleMass)/(std::sqrt(ti*(ti+2.*ParticleMass))*lossi);
472   }                                               474   }
473   Value *= ParticleMass*dltau/c_light;            475   Value *= ParticleMass*dltau/c_light;
474   return Value;                                   476   return Value;
475 }                                                 477 }
476                                                   478 
477 //                                                479 //    
478                                                   480 
479 G4double G4RDVeLowEnergyLoss::ProperTimeIntLog    481 G4double G4RDVeLowEnergyLoss::ProperTimeIntLog(G4PhysicsVector* physicsVector,
480                G4int nbin)                        482                G4int nbin)
481 //  num. integration, logarithmic binning         483 //  num. integration, logarithmic binning
482 {                                                 484 {
483   G4double ltt,dltau,Value,ui,taui,ti,lossi,ci    485   G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
484   G4bool isOut;                                   486   G4bool isOut;
485   ltt = ltauhigh-ltaulow;                         487   ltt = ltauhigh-ltaulow;
486   dltau = ltt/nbin;                               488   dltau = ltt/nbin;
487   Value = 0.;                                     489   Value = 0.;
488                                                   490 
489   for (G4int i=0; i<=nbin; i++)                   491   for (G4int i=0; i<=nbin; i++)
490   {                                               492   {
491     ui = ltaulow+dltau*i;                         493     ui = ltaulow+dltau*i;
492     taui = std::exp(ui);                          494     taui = std::exp(ui);
493     ti = ParticleMass*taui;                       495     ti = ParticleMass*taui;
494     lossi = physicsVector->GetValue(ti,isOut);    496     lossi = physicsVector->GetValue(ti,isOut);
495     if(i==0)                                      497     if(i==0)
496       ci=0.5;                                     498       ci=0.5;
497     else                                          499     else
498     {                                             500     {
499       if(i<nbin)                                  501       if(i<nbin)
500         ci=1.;                                    502         ci=1.;
501       else                                        503       else
502         ci=0.5;                                   504         ci=0.5;
503     }                                             505     }
504     Value += ci*taui*ParticleMass/(std::sqrt(t    506     Value += ci*taui*ParticleMass/(std::sqrt(ti*(ti+2.*ParticleMass))*lossi);
505   }                                               507   }
506   Value *= ParticleMass*dltau/c_light;            508   Value *= ParticleMass*dltau/c_light;
507   return Value;                                   509   return Value;
508 }                                                 510 }
509                                                   511 
510 //                                                512 //    
511                                                   513 
512 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildInve    514 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildInverseRangeTable(G4PhysicsTable* theRangeTable,
513                 G4PhysicsTable*,                  515                 G4PhysicsTable*,
514                 G4PhysicsTable*,                  516                 G4PhysicsTable*,
515                 G4PhysicsTable*,                  517                 G4PhysicsTable*,
516                 G4PhysicsTable* theInverseRang    518                 G4PhysicsTable* theInverseRangeTable,
517                 G4double, // lowestKineticEner    519                 G4double, // lowestKineticEnergy,
518                 G4double, // highestKineticEne    520                 G4double, // highestKineticEnergy
519                 G4int )   // nbins                521                 G4int )   // nbins
520 // Build inverse table of the range table         522 // Build inverse table of the range table
521 {                                                 523 {
522   G4bool b;                                       524   G4bool b;
523                                                   525 
524   G4int numOfCouples = G4ProductionCutsTable::    526   G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
525                                                   527 
526     if(theInverseRangeTable)                      528     if(theInverseRangeTable)
527     { theInverseRangeTable->clearAndDestroy();    529     { theInverseRangeTable->clearAndDestroy();
528       delete theInverseRangeTable; }              530       delete theInverseRangeTable; }
529     theInverseRangeTable = new G4PhysicsTable(    531     theInverseRangeTable = new G4PhysicsTable(numOfCouples);
530                                                   532 
531   // loop for materials                           533   // loop for materials
532   for (G4int i=0;  i<numOfCouples; i++)           534   for (G4int i=0;  i<numOfCouples; i++)
533   {                                               535   {
534                                                   536 
535     G4PhysicsVector* pv = (*theRangeTable)[i];    537     G4PhysicsVector* pv = (*theRangeTable)[i];
536     size_t nbins   = pv->GetVectorLength();       538     size_t nbins   = pv->GetVectorLength();
537     G4double elow  = pv->GetLowEdgeEnergy(0);     539     G4double elow  = pv->GetLowEdgeEnergy(0);
538     G4double ehigh = pv->GetLowEdgeEnergy(nbin    540     G4double ehigh = pv->GetLowEdgeEnergy(nbins-1);
539     G4double rlow  = pv->GetValue(elow, b);       541     G4double rlow  = pv->GetValue(elow, b);
540     G4double rhigh = pv->GetValue(ehigh, b);      542     G4double rhigh = pv->GetValue(ehigh, b);
541                                                   543 
542     rhigh *= std::exp(std::log(rhigh/rlow)/((G    544     rhigh *= std::exp(std::log(rhigh/rlow)/((G4double)(nbins-1)));
543                                                   545 
544     G4PhysicsLogVector* v = new G4PhysicsLogVe    546     G4PhysicsLogVector* v = new G4PhysicsLogVector(rlow, rhigh, nbins);
545                                                   547 
546     v->PutValue(0,elow);                          548     v->PutValue(0,elow);
547     G4double energy1 = elow;                      549     G4double energy1 = elow;
548     G4double range1  = rlow;                      550     G4double range1  = rlow;
549     G4double energy2 = elow;                      551     G4double energy2 = elow;
550     G4double range2  = rlow;                      552     G4double range2  = rlow;
551     size_t ilow      = 0;                         553     size_t ilow      = 0;
552     size_t ihigh;                                 554     size_t ihigh;
553                                                   555 
554     for (size_t j=1; j<nbins; j++) {              556     for (size_t j=1; j<nbins; j++) {
555                                                   557 
556       G4double range = v->GetLowEdgeEnergy(j);    558       G4double range = v->GetLowEdgeEnergy(j);
557                                                   559 
558       for (ihigh=ilow+1; ihigh<nbins; ihigh++)    560       for (ihigh=ilow+1; ihigh<nbins; ihigh++) {
559         energy2 = pv->GetLowEdgeEnergy(ihigh);    561         energy2 = pv->GetLowEdgeEnergy(ihigh);
560         range2  = pv->GetValue(energy2, b);       562         range2  = pv->GetValue(energy2, b);
561         if(range2 >= range || ihigh == nbins-1    563         if(range2 >= range || ihigh == nbins-1) {
562           ilow = ihigh - 1;                       564           ilow = ihigh - 1;
563           energy1 = pv->GetLowEdgeEnergy(ilow)    565           energy1 = pv->GetLowEdgeEnergy(ilow);
564           range1  = pv->GetValue(energy1, b);     566           range1  = pv->GetValue(energy1, b); 
565           break;                                  567           break;
566   }                                               568   }
567       }                                           569       }
568                                                   570 
569       G4double e = std::log(energy1) + std::lo    571       G4double e = std::log(energy1) + std::log(energy2/energy1)*std::log(range/range1)/std::log(range2/range1);
570                                                   572 
571       v->PutValue(j,std::exp(e));                 573       v->PutValue(j,std::exp(e));
572     }                                             574     }
573     theInverseRangeTable->insert(v);              575     theInverseRangeTable->insert(v);
574                                                   576 
575   }                                               577   }
576   return theInverseRangeTable ;                   578   return theInverseRangeTable ;
577 }                                                 579 }
578                                                   580 
579 //                                                581 //    
580                                                   582 
581 void G4RDVeLowEnergyLoss::InvertRangeVector(G4    583 void G4RDVeLowEnergyLoss::InvertRangeVector(G4PhysicsTable* theRangeTable,
582             G4PhysicsTable* theRangeCoeffATabl    584             G4PhysicsTable* theRangeCoeffATable,
583             G4PhysicsTable* theRangeCoeffBTabl    585             G4PhysicsTable* theRangeCoeffBTable,
584             G4PhysicsTable* theRangeCoeffCTabl    586             G4PhysicsTable* theRangeCoeffCTable,
585             G4double lowestKineticEnergy,         587             G4double lowestKineticEnergy,
586             G4double highestKineticEnergy, G4i    588             G4double highestKineticEnergy, G4int TotBin,
587             G4int  materialIndex, G4PhysicsLog    589             G4int  materialIndex, G4PhysicsLogVector* aVector)
588 //  invert range vector for a material            590 //  invert range vector for a material
589 {                                                 591 {
590   G4double LowEdgeRange,A,B,C,discr,KineticEne    592   G4double LowEdgeRange,A,B,C,discr,KineticEnergy ;
591   G4double RTable = std::exp(std::log(highestK    593   G4double RTable = std::exp(std::log(highestKineticEnergy/lowestKineticEnergy)/TotBin) ;
592   G4double Tbin = lowestKineticEnergy/RTable ;    594   G4double Tbin = lowestKineticEnergy/RTable ;
593   G4double rangebin = 0.0 ;                       595   G4double rangebin = 0.0 ;
594   G4int binnumber = -1 ;                          596   G4int binnumber = -1 ;
595   G4bool isOut ;                                  597   G4bool isOut ;
596                                                   598 
597   //loop for range values                         599   //loop for range values
598   for( G4int i=0; i<TotBin; i++)                  600   for( G4int i=0; i<TotBin; i++)
599   {                                               601   {
600     LowEdgeRange = aVector->GetLowEdgeEnergy(i    602     LowEdgeRange = aVector->GetLowEdgeEnergy(i) ;  //i.e. GetLowEdgeValue(i)
601     if( rangebin < LowEdgeRange )                 603     if( rangebin < LowEdgeRange )
602     {                                             604     {
603       do                                          605       do
604       {                                           606       {
605         binnumber += 1 ;                          607         binnumber += 1 ;
606         Tbin *= RTable ;                          608         Tbin *= RTable ;
607         rangebin = (*theRangeTable)(materialIn    609         rangebin = (*theRangeTable)(materialIndex)->GetValue(Tbin,isOut) ;
608       }                                           610       }
609       while ((rangebin < LowEdgeRange) && (bin    611       while ((rangebin < LowEdgeRange) && (binnumber < TotBin )) ;
610     }                                             612     }
611                                                   613 
612     if(binnumber == 0)                            614     if(binnumber == 0)
613       KineticEnergy = lowestKineticEnergy ;       615       KineticEnergy = lowestKineticEnergy ;
614     else if(binnumber == TotBin-1)                616     else if(binnumber == TotBin-1)
615       KineticEnergy = highestKineticEnergy ;      617       KineticEnergy = highestKineticEnergy ;
616     else                                          618     else
617     {                                             619     {
618       A = (*(*theRangeCoeffATable)(materialInd    620       A = (*(*theRangeCoeffATable)(materialIndex))(binnumber-1) ;
619       B = (*(*theRangeCoeffBTable)(materialInd    621       B = (*(*theRangeCoeffBTable)(materialIndex))(binnumber-1) ;
620       C = (*(*theRangeCoeffCTable)(materialInd    622       C = (*(*theRangeCoeffCTable)(materialIndex))(binnumber-1) ;
621       if(A==0.)                                   623       if(A==0.)
622          KineticEnergy = (LowEdgeRange -C )/B     624          KineticEnergy = (LowEdgeRange -C )/B ;
623       else                                        625       else
624       {                                           626       {
625          discr = B*B - 4.*A*(C-LowEdgeRange);     627          discr = B*B - 4.*A*(C-LowEdgeRange);
626          discr = discr>0. ? std::sqrt(discr) :    628          discr = discr>0. ? std::sqrt(discr) : 0.;
627          KineticEnergy = 0.5*(discr-B)/A ;        629          KineticEnergy = 0.5*(discr-B)/A ;
628       }                                           630       }
629     }                                             631     }
630                                                   632 
631     aVector->PutValue(i,KineticEnergy) ;          633     aVector->PutValue(i,KineticEnergy) ;
632   }                                               634   }
633 }                                                 635 }
634                                                   636 
635 //                                                637 //    
636                                                   638 
637 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRang    639 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRangeCoeffATable(G4PhysicsTable* theRangeTable,
638                G4PhysicsTable* theRangeCoeffAT    640                G4PhysicsTable* theRangeCoeffATable,
639                G4double lowestKineticEnergy,      641                G4double lowestKineticEnergy,
640                G4double highestKineticEnergy,     642                G4double highestKineticEnergy, G4int TotBin)
641 // Build tables of coefficients for the energy    643 // Build tables of coefficients for the energy loss calculation
642 //  create table for coefficients "A"             644 //  create table for coefficients "A"
643 {                                                 645 {
644                                                   646 
645   G4int numOfCouples = G4ProductionCutsTable::    647   G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
646                                                   648 
647   if(theRangeCoeffATable)                         649   if(theRangeCoeffATable)
648   { theRangeCoeffATable->clearAndDestroy();       650   { theRangeCoeffATable->clearAndDestroy();
649     delete theRangeCoeffATable; }                 651     delete theRangeCoeffATable; }
650   theRangeCoeffATable = new G4PhysicsTable(num    652   theRangeCoeffATable = new G4PhysicsTable(numOfCouples);
651                                                   653 
652   G4double RTable = std::exp(std::log(highestK    654   G4double RTable = std::exp(std::log(highestKineticEnergy/lowestKineticEnergy)/TotBin) ;
653   G4double R2 = RTable*RTable ;                   655   G4double R2 = RTable*RTable ;
654   G4double R1 = RTable+1.;                        656   G4double R1 = RTable+1.;
655   G4double w = R1*(RTable-1.)*(RTable-1.);        657   G4double w = R1*(RTable-1.)*(RTable-1.);
656   G4double w1 = RTable/w , w2 = -RTable*R1/w ,    658   G4double w1 = RTable/w , w2 = -RTable*R1/w , w3 = R2/w ;
657   G4double Ti , Tim , Tip , Ri , Rim , Rip , V    659   G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
658   G4bool isOut;                                   660   G4bool isOut;
659                                                   661 
660   //  loop for materials                          662   //  loop for materials
661   for (G4int J=0; J<numOfCouples; J++)            663   for (G4int J=0; J<numOfCouples; J++)
662   {                                               664   {
663     G4int binmax=TotBin ;                         665     G4int binmax=TotBin ;
664     G4PhysicsLinearVector* aVector =              666     G4PhysicsLinearVector* aVector =
665                            new G4PhysicsLinear    667                            new G4PhysicsLinearVector(0.,binmax, TotBin);
666     Ti = lowestKineticEnergy ;                    668     Ti = lowestKineticEnergy ;
667     G4PhysicsVector* rangeVector= (*theRangeTa    669     G4PhysicsVector* rangeVector= (*theRangeTable)[J];
668                                                   670 
669     for ( G4int i=0; i<TotBin; i++)               671     for ( G4int i=0; i<TotBin; i++)
670     {                                             672     {
671       Ri = rangeVector->GetValue(Ti,isOut) ;      673       Ri = rangeVector->GetValue(Ti,isOut) ;
672       if ( i==0 )                                 674       if ( i==0 )
673         Rim = 0. ;                                675         Rim = 0. ;
674       else                                        676       else
675       {                                           677       {
676         Tim = Ti/RTable ;                         678         Tim = Ti/RTable ;
677         Rim = rangeVector->GetValue(Tim,isOut)    679         Rim = rangeVector->GetValue(Tim,isOut);
678       }                                           680       }
679       if ( i==(TotBin-1))                         681       if ( i==(TotBin-1))
680         Rip = Ri ;                                682         Rip = Ri ;
681       else                                        683       else
682       {                                           684       {
683         Tip = Ti*RTable ;                         685         Tip = Ti*RTable ;
684         Rip = rangeVector->GetValue(Tip,isOut)    686         Rip = rangeVector->GetValue(Tip,isOut);
685       }                                           687       }
686       Value = (w1*Rip + w2*Ri + w3*Rim)/(Ti*Ti    688       Value = (w1*Rip + w2*Ri + w3*Rim)/(Ti*Ti) ;
687                                                   689 
688       aVector->PutValue(i,Value);                 690       aVector->PutValue(i,Value);
689       Ti = RTable*Ti ;                            691       Ti = RTable*Ti ;
690     }                                             692     }
691                                                   693  
692     theRangeCoeffATable->insert(aVector);         694     theRangeCoeffATable->insert(aVector);
693   }                                               695   }
694   return theRangeCoeffATable ;                    696   return theRangeCoeffATable ;
695 }                                                 697 }
696                                                   698 
697 //                                                699 //    
698                                                   700   
699 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRang    701 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRangeCoeffBTable(G4PhysicsTable* theRangeTable,
700                G4PhysicsTable* theRangeCoeffBT    702                G4PhysicsTable* theRangeCoeffBTable,
701                G4double lowestKineticEnergy,      703                G4double lowestKineticEnergy,
702                G4double highestKineticEnergy,     704                G4double highestKineticEnergy, G4int TotBin)
703 // Build tables of coefficients for the energy    705 // Build tables of coefficients for the energy loss calculation
704 //  create table for coefficients "B"             706 //  create table for coefficients "B"
705 {                                                 707 {
706                                                   708 
707   G4int numOfCouples = G4ProductionCutsTable::    709   G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
708                                                   710 
709   if(theRangeCoeffBTable)                         711   if(theRangeCoeffBTable)
710   { theRangeCoeffBTable->clearAndDestroy();       712   { theRangeCoeffBTable->clearAndDestroy();
711     delete theRangeCoeffBTable; }                 713     delete theRangeCoeffBTable; }
712   theRangeCoeffBTable = new G4PhysicsTable(num    714   theRangeCoeffBTable = new G4PhysicsTable(numOfCouples);
713                                                   715 
714   G4double RTable = std::exp(std::log(highestK    716   G4double RTable = std::exp(std::log(highestKineticEnergy/lowestKineticEnergy)/TotBin) ;
715   G4double R2 = RTable*RTable ;                   717   G4double R2 = RTable*RTable ;
716   G4double R1 = RTable+1.;                        718   G4double R1 = RTable+1.;
717   G4double w = R1*(RTable-1.)*(RTable-1.);        719   G4double w = R1*(RTable-1.)*(RTable-1.);
718   G4double w1 = -R1/w , w2 = R1*(R2+1.)/w , w3    720   G4double w1 = -R1/w , w2 = R1*(R2+1.)/w , w3 = -R2*R1/w ;
719   G4double Ti , Tim , Tip , Ri , Rim , Rip , V    721   G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
720   G4bool isOut;                                   722   G4bool isOut;
721                                                   723 
722   //  loop for materials                          724   //  loop for materials
723   for (G4int J=0; J<numOfCouples; J++)            725   for (G4int J=0; J<numOfCouples; J++)
724   {                                               726   {
725     G4int binmax=TotBin ;                         727     G4int binmax=TotBin ;
726     G4PhysicsLinearVector* aVector =              728     G4PhysicsLinearVector* aVector =
727                         new G4PhysicsLinearVec    729                         new G4PhysicsLinearVector(0.,binmax, TotBin);
728     Ti = lowestKineticEnergy ;                    730     Ti = lowestKineticEnergy ;
729     G4PhysicsVector* rangeVector= (*theRangeTa    731     G4PhysicsVector* rangeVector= (*theRangeTable)[J];
730                                                   732   
731     for ( G4int i=0; i<TotBin; i++)               733     for ( G4int i=0; i<TotBin; i++)
732     {                                             734     {
733       Ri = rangeVector->GetValue(Ti,isOut) ;      735       Ri = rangeVector->GetValue(Ti,isOut) ;
734       if ( i==0 )                                 736       if ( i==0 )
735          Rim = 0. ;                               737          Rim = 0. ;
736       else                                        738       else
737       {                                           739       {
738         Tim = Ti/RTable ;                         740         Tim = Ti/RTable ;
739         Rim = rangeVector->GetValue(Tim,isOut)    741         Rim = rangeVector->GetValue(Tim,isOut);
740       }                                           742       }
741       if ( i==(TotBin-1))                         743       if ( i==(TotBin-1))
742         Rip = Ri ;                                744         Rip = Ri ;
743       else                                        745       else
744       {                                           746       {
745         Tip = Ti*RTable ;                         747         Tip = Ti*RTable ;
746         Rip = rangeVector->GetValue(Tip,isOut)    748         Rip = rangeVector->GetValue(Tip,isOut);
747       }                                           749       }
748       Value = (w1*Rip + w2*Ri + w3*Rim)/Ti;       750       Value = (w1*Rip + w2*Ri + w3*Rim)/Ti;
749                                                   751 
750       aVector->PutValue(i,Value);                 752       aVector->PutValue(i,Value);
751       Ti = RTable*Ti ;                            753       Ti = RTable*Ti ;
752     }                                             754     }
753     theRangeCoeffBTable->insert(aVector);         755     theRangeCoeffBTable->insert(aVector);
754   }                                               756   }
755   return theRangeCoeffBTable ;                    757   return theRangeCoeffBTable ;
756 }                                                 758 }
757                                                   759 
758 //                                                760 //    
759                                                   761 
760 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRang    762 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRangeCoeffCTable(G4PhysicsTable* theRangeTable,
761                G4PhysicsTable* theRangeCoeffCT    763                G4PhysicsTable* theRangeCoeffCTable,
762                G4double lowestKineticEnergy,      764                G4double lowestKineticEnergy,
763                G4double highestKineticEnergy,     765                G4double highestKineticEnergy, G4int TotBin)
764 // Build tables of coefficients for the energy    766 // Build tables of coefficients for the energy loss calculation
765 //  create table for coefficients "C"             767 //  create table for coefficients "C"
766 {                                                 768 {
767                                                   769 
768   G4int numOfCouples = G4ProductionCutsTable::    770   G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
769                                                   771 
770   if(theRangeCoeffCTable)                         772   if(theRangeCoeffCTable)
771   { theRangeCoeffCTable->clearAndDestroy();       773   { theRangeCoeffCTable->clearAndDestroy();
772     delete theRangeCoeffCTable; }                 774     delete theRangeCoeffCTable; }
773   theRangeCoeffCTable = new G4PhysicsTable(num    775   theRangeCoeffCTable = new G4PhysicsTable(numOfCouples);
774                                                   776 
775   G4double RTable = std::exp(std::log(highestK    777   G4double RTable = std::exp(std::log(highestKineticEnergy/lowestKineticEnergy)/TotBin) ;
776   G4double R2 = RTable*RTable ;                   778   G4double R2 = RTable*RTable ;
777   G4double R1 = RTable+1.;                        779   G4double R1 = RTable+1.;
778   G4double w = R1*(RTable-1.)*(RTable-1.);        780   G4double w = R1*(RTable-1.)*(RTable-1.);
779   G4double w1 = 1./w , w2 = -RTable*R1/w , w3     781   G4double w1 = 1./w , w2 = -RTable*R1/w , w3 = RTable*R2/w ;
780   G4double Ti , Tim , Tip , Ri , Rim , Rip , V    782   G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
781   G4bool isOut;                                   783   G4bool isOut;
782                                                   784 
783   //  loop for materials                          785   //  loop for materials
784   for (G4int J=0; J<numOfCouples; J++)            786   for (G4int J=0; J<numOfCouples; J++)
785   {                                               787   {
786     G4int binmax=TotBin ;                         788     G4int binmax=TotBin ;
787     G4PhysicsLinearVector* aVector =              789     G4PhysicsLinearVector* aVector =
788                       new G4PhysicsLinearVecto    790                       new G4PhysicsLinearVector(0.,binmax, TotBin);
789     Ti = lowestKineticEnergy ;                    791     Ti = lowestKineticEnergy ;
790     G4PhysicsVector* rangeVector= (*theRangeTa    792     G4PhysicsVector* rangeVector= (*theRangeTable)[J];
791                                                   793   
792     for ( G4int i=0; i<TotBin; i++)               794     for ( G4int i=0; i<TotBin; i++)
793     {                                             795     {
794       Ri = rangeVector->GetValue(Ti,isOut) ;      796       Ri = rangeVector->GetValue(Ti,isOut) ;
795       if ( i==0 )                                 797       if ( i==0 )
796         Rim = 0. ;                                798         Rim = 0. ;
797       else                                        799       else
798       {                                           800       {
799         Tim = Ti/RTable ;                         801         Tim = Ti/RTable ;
800         Rim = rangeVector->GetValue(Tim,isOut)    802         Rim = rangeVector->GetValue(Tim,isOut);
801       }                                           803       }
802       if ( i==(TotBin-1))                         804       if ( i==(TotBin-1))
803         Rip = Ri ;                                805         Rip = Ri ;
804       else                                        806       else
805       {                                           807       {
806         Tip = Ti*RTable ;                         808         Tip = Ti*RTable ;
807         Rip = rangeVector->GetValue(Tip,isOut)    809         Rip = rangeVector->GetValue(Tip,isOut);
808       }                                           810       }
809       Value = w1*Rip + w2*Ri + w3*Rim ;           811       Value = w1*Rip + w2*Ri + w3*Rim ;
810                                                   812 
811       aVector->PutValue(i,Value);                 813       aVector->PutValue(i,Value);
812       Ti = RTable*Ti ;                            814       Ti = RTable*Ti ;
813     }                                             815     }
814     theRangeCoeffCTable->insert(aVector);         816     theRangeCoeffCTable->insert(aVector);
815   }                                               817   }
816   return theRangeCoeffCTable ;                    818   return theRangeCoeffCTable ;
817 }                                                 819 }
818                                                   820 
819 //                                                821 //    
820                                                   822 
821 G4double G4RDVeLowEnergyLoss::GetLossWithFluct    823 G4double G4RDVeLowEnergyLoss::GetLossWithFluct(const G4DynamicParticle* aParticle,
822                                              c    824                                              const G4MaterialCutsCouple* couple,
823                G4double MeanLoss,                 825                G4double MeanLoss,
824                G4double step)                     826                G4double step)
825 //  calculate actual loss from the mean loss      827 //  calculate actual loss from the mean loss
826 //  The model used to get the fluctuation is e    828 //  The model used to get the fluctuation is essentially the same as in Glandz in Geant3.
827 {                                                 829 {
828    static const G4double minLoss = 1.*eV ;        830    static const G4double minLoss = 1.*eV ;
829    static const G4double probLim = 0.01 ;         831    static const G4double probLim = 0.01 ;
830    static const G4double sumaLim = -std::log(p    832    static const G4double sumaLim = -std::log(probLim) ;
831    static const G4double alim=10.;                833    static const G4double alim=10.;
832    static const G4double kappa = 10. ;            834    static const G4double kappa = 10. ;
833    static const G4double factor = twopi_mc2_rc    835    static const G4double factor = twopi_mc2_rcl2 ;
834   const G4Material* aMaterial = couple->GetMat    836   const G4Material* aMaterial = couple->GetMaterial();
835                                                   837 
836   // check if the material has changed ( cache    838   // check if the material has changed ( cache mechanism)
837                                                   839 
838   if (aMaterial != lastMaterial)                  840   if (aMaterial != lastMaterial)
839     {                                             841     {
840       lastMaterial = aMaterial;                   842       lastMaterial = aMaterial;
841       imat         = couple->GetIndex();          843       imat         = couple->GetIndex();
842       f1Fluct      = aMaterial->GetIonisation(    844       f1Fluct      = aMaterial->GetIonisation()->GetF1fluct();
843       f2Fluct      = aMaterial->GetIonisation(    845       f2Fluct      = aMaterial->GetIonisation()->GetF2fluct();
844       e1Fluct      = aMaterial->GetIonisation(    846       e1Fluct      = aMaterial->GetIonisation()->GetEnergy1fluct();
845       e2Fluct      = aMaterial->GetIonisation(    847       e2Fluct      = aMaterial->GetIonisation()->GetEnergy2fluct();
846       e1LogFluct   = aMaterial->GetIonisation(    848       e1LogFluct   = aMaterial->GetIonisation()->GetLogEnergy1fluct();
847       e2LogFluct   = aMaterial->GetIonisation(    849       e2LogFluct   = aMaterial->GetIonisation()->GetLogEnergy2fluct();
848       rateFluct    = aMaterial->GetIonisation(    850       rateFluct    = aMaterial->GetIonisation()->GetRateionexcfluct();
849       ipotFluct    = aMaterial->GetIonisation(    851       ipotFluct    = aMaterial->GetIonisation()->GetMeanExcitationEnergy();
850       ipotLogFluct = aMaterial->GetIonisation(    852       ipotLogFluct = aMaterial->GetIonisation()->GetLogMeanExcEnergy();
851     }                                             853     }
852   G4double threshold,w1,w2,C,                     854   G4double threshold,w1,w2,C,
853            beta2,suma,e0,loss,lossc,w;            855            beta2,suma,e0,loss,lossc,w;
854   G4double a1,a2,a3;                              856   G4double a1,a2,a3;
855   G4int p1,p2,p3;                                 857   G4int p1,p2,p3;
856   G4int nb;                                       858   G4int nb;
857   G4double Corrfac, na,alfa,rfac,namean,sa,alf    859   G4double Corrfac, na,alfa,rfac,namean,sa,alfa1,ea,sea;
858   //  G4double dp1;                               860   //  G4double dp1;
859   G4double dp3;                                   861   G4double dp3;
860   G4double siga ;                                 862   G4double siga ;
861                                                   863 
862   // shortcut for very very small loss            864   // shortcut for very very small loss
863   if(MeanLoss < minLoss) return MeanLoss ;        865   if(MeanLoss < minLoss) return MeanLoss ;
864                                                   866 
865   // get particle data                            867   // get particle data
866   G4double Tkin   = aParticle->GetKineticEnerg    868   G4double Tkin   = aParticle->GetKineticEnergy();
867                                                   869 
868   //  G4cout << "MGP -- Fluc Tkin " << Tkin/ke    870   //  G4cout << "MGP -- Fluc Tkin " << Tkin/keV << " keV "  << " MeanLoss = " << MeanLoss/keV << G4endl;
869                                                   871 
870   threshold = (*((G4ProductionCutsTable::GetPr    872   threshold = (*((G4ProductionCutsTable::GetProductionCutsTable())
871                 ->GetEnergyCutsVector(1)))[ima    873                 ->GetEnergyCutsVector(1)))[imat];
872   G4double rmass = electron_mass_c2/ParticleMa    874   G4double rmass = electron_mass_c2/ParticleMass;
873   G4double tau   = Tkin/ParticleMass, tau1 = t    875   G4double tau   = Tkin/ParticleMass, tau1 = tau+1., tau2 = tau*(tau+2.);
874   G4double Tm    = 2.*electron_mass_c2*tau2/(1    876   G4double Tm    = 2.*electron_mass_c2*tau2/(1.+2.*tau1*rmass+rmass*rmass);
875                                                   877 
876   // G4cout << "MGP Particle mass " << Particl    878   // G4cout << "MGP Particle mass " << ParticleMass/MeV << " Tm " << Tm << G4endl;
877                                                   879 
878   if(Tm > threshold) Tm = threshold;              880   if(Tm > threshold) Tm = threshold;
879   beta2 = tau2/(tau1*tau1);                       881   beta2 = tau2/(tau1*tau1);
880                                                   882 
881   // Gaussian fluctuation ?                       883   // Gaussian fluctuation ?
882   if(MeanLoss >= kappa*Tm || MeanLoss <= kappa    884   if(MeanLoss >= kappa*Tm || MeanLoss <= kappa*ipotFluct)
883   {                                               885   {
884     G4double electronDensity = aMaterial->GetE    886     G4double electronDensity = aMaterial->GetElectronDensity() ;
885     siga = std::sqrt(Tm*(1.0-0.5*beta2)*step*     887     siga = std::sqrt(Tm*(1.0-0.5*beta2)*step*
886                 factor*electronDensity/beta2)     888                 factor*electronDensity/beta2) ;
887     do {                                          889     do {
888       loss = G4RandGauss::shoot(MeanLoss,siga)    890       loss = G4RandGauss::shoot(MeanLoss,siga) ;
889     } while (loss < 0. || loss > 2.0*MeanLoss)    891     } while (loss < 0. || loss > 2.0*MeanLoss);
890     return loss ;                                 892     return loss ;
891   }                                               893   }
892                                                   894 
893   w1 = Tm/ipotFluct;                              895   w1 = Tm/ipotFluct;
894   w2 = std::log(2.*electron_mass_c2*tau2);        896   w2 = std::log(2.*electron_mass_c2*tau2);
895                                                   897 
896   C = MeanLoss*(1.-rateFluct)/(w2-ipotLogFluct    898   C = MeanLoss*(1.-rateFluct)/(w2-ipotLogFluct-beta2);
897                                                   899 
898   a1 = C*f1Fluct*(w2-e1LogFluct-beta2)/e1Fluct    900   a1 = C*f1Fluct*(w2-e1LogFluct-beta2)/e1Fluct;
899   a2 = C*f2Fluct*(w2-e2LogFluct-beta2)/e2Fluct    901   a2 = C*f2Fluct*(w2-e2LogFluct-beta2)/e2Fluct;
900   a3 = rateFluct*MeanLoss*(Tm-ipotFluct)/(ipot    902   a3 = rateFluct*MeanLoss*(Tm-ipotFluct)/(ipotFluct*Tm*std::log(w1));
901                                                   903 
902   suma = a1+a2+a3;                                904   suma = a1+a2+a3;
903                                                   905 
904   loss = 0. ;                                     906   loss = 0. ;
905                                                   907 
906   if(suma < sumaLim)             // very small    908   if(suma < sumaLim)             // very small Step
907     {                                             909     {
908       e0 = aMaterial->GetIonisation()->GetEner    910       e0 = aMaterial->GetIonisation()->GetEnergy0fluct();
909       // G4cout << "MGP e0 = " << e0/keV << G4    911       // G4cout << "MGP e0 = " << e0/keV << G4endl;
910                                                   912 
911       if(Tm == ipotFluct)                         913       if(Tm == ipotFluct)
912   {                                               914   {
913     a3 = MeanLoss/e0;                             915     a3 = MeanLoss/e0;
914                                                   916 
915     if(a3>alim)                                   917     if(a3>alim)
916       {                                           918       {
917         siga=std::sqrt(a3) ;                      919         siga=std::sqrt(a3) ;
918         p3 = std::max(0,G4int(G4RandGauss::sho    920         p3 = std::max(0,G4int(G4RandGauss::shoot(a3,siga)+0.5));
919       }                                           921       }
920     else p3 = G4Poisson(a3);                      922     else p3 = G4Poisson(a3);
921                                                   923 
922     loss = p3*e0 ;                                924     loss = p3*e0 ;
923                                                   925 
924     if(p3 > 0) loss += (1.-2.*G4UniformRand())    926     if(p3 > 0) loss += (1.-2.*G4UniformRand())*e0 ;
925     // G4cout << "MGP very small step " << los    927     // G4cout << "MGP very small step " << loss/keV << G4endl;
926   }                                               928   }
927       else                                        929       else
928   {                                               930   {
929     //    G4cout << "MGP old Tm = " << Tm << "    931     //    G4cout << "MGP old Tm = " << Tm << " " << ipotFluct << " " << e0 << G4endl;
930     Tm = Tm-ipotFluct+e0 ;                        932     Tm = Tm-ipotFluct+e0 ;
931                                                   933 
932     // MGP ---- workaround to avoid log argume    934     // MGP ---- workaround to avoid log argument<0, TO BE CHECKED
933     if (Tm <= 0.)                                 935     if (Tm <= 0.)
934       {                                           936       {
935         loss = MeanLoss;                          937         loss = MeanLoss;
936         p3 = 0;                                   938         p3 = 0;
937         // G4cout << "MGP correction loss = Me    939         // G4cout << "MGP correction loss = MeanLoss " << loss/keV << G4endl;
938       }                                           940       }
939     else                                          941     else
940       {                                           942       {
941         a3 = MeanLoss*(Tm-e0)/(Tm*e0*std::log(    943         a3 = MeanLoss*(Tm-e0)/(Tm*e0*std::log(Tm/e0));
942                                                   944 
943         // G4cout << "MGP new Tm = " << Tm <<     945         // G4cout << "MGP new Tm = " << Tm << " " << ipotFluct << " " << e0 << " a3= " << a3 << G4endl;
944                                                   946         
945         if(a3>alim)                               947         if(a3>alim)
946     {                                             948     {
947       siga=std::sqrt(a3) ;                        949       siga=std::sqrt(a3) ;
948       p3 = std::max(0,G4int(G4RandGauss::shoot    950       p3 = std::max(0,G4int(G4RandGauss::shoot(a3,siga)+0.5));
949     }                                             951     }
950         else                                      952         else
951     p3 = G4Poisson(a3);                           953     p3 = G4Poisson(a3);
952         //G4cout << "MGP p3 " << p3 << G4endl;    954         //G4cout << "MGP p3 " << p3 << G4endl;
953                                                   955 
954       }                                           956       }
955                                                   957         
956     if(p3 > 0)                                    958     if(p3 > 0)
957       {                                           959       {
958         w = (Tm-e0)/Tm ;                          960         w = (Tm-e0)/Tm ;
959         if(p3 > nmaxCont2)                        961         if(p3 > nmaxCont2)
960     {                                             962     {
961       // G4cout << "MGP dp3 " << dp3 << " p3 "    963       // G4cout << "MGP dp3 " << dp3 << " p3 " << p3 << " " << nmaxCont2 << G4endl;
962       dp3 = G4double(p3) ;                        964       dp3 = G4double(p3) ;
963       Corrfac = dp3/G4double(nmaxCont2) ;         965       Corrfac = dp3/G4double(nmaxCont2) ;
964       p3 = nmaxCont2 ;                            966       p3 = nmaxCont2 ;
965     }                                             967     }
966         else                                      968         else
967     Corrfac = 1. ;                                969     Corrfac = 1. ;
968                                                   970         
969         for(G4int i=0; i<p3; i++) loss += 1./(    971         for(G4int i=0; i<p3; i++) loss += 1./(1.-w*G4UniformRand()) ;
970         loss *= e0*Corrfac ;                      972         loss *= e0*Corrfac ; 
971         // G4cout << "MGP Corrfac = " << Corrf    973         // G4cout << "MGP Corrfac = " << Corrfac << " e0 = " << e0/keV << " loss = " << loss/keV << G4endl;
972       }                                           974       }        
973   }                                               975   }
974     }                                             976     }
975                                                   977 
976   else                              // not so     978   else                              // not so small Step
977     {                                             979     {
978       // excitation type 1                        980       // excitation type 1
979       if(a1>alim)                                 981       if(a1>alim)
980       {                                           982       {
981         siga=std::sqrt(a1) ;                      983         siga=std::sqrt(a1) ;
982         p1 = std::max(0,int(G4RandGauss::shoot    984         p1 = std::max(0,int(G4RandGauss::shoot(a1,siga)+0.5));
983       }                                           985       }
984       else                                        986       else
985        p1 = G4Poisson(a1);                        987        p1 = G4Poisson(a1);
986                                                   988 
987       // excitation type 2                        989       // excitation type 2
988       if(a2>alim)                                 990       if(a2>alim)
989       {                                           991       {
990         siga=std::sqrt(a2) ;                      992         siga=std::sqrt(a2) ;
991         p2 = std::max(0,int(G4RandGauss::shoot    993         p2 = std::max(0,int(G4RandGauss::shoot(a2,siga)+0.5));
992       }                                           994       }
993       else                                        995       else
994         p2 = G4Poisson(a2);                       996         p2 = G4Poisson(a2);
995                                                   997 
996       loss = p1*e1Fluct+p2*e2Fluct;               998       loss = p1*e1Fluct+p2*e2Fluct;
997                                                   999  
998       // smearing to avoid unphysical peaks       1000       // smearing to avoid unphysical peaks
999       if(p2 > 0)                                  1001       if(p2 > 0)
1000         loss += (1.-2.*G4UniformRand())*e2Flu    1002         loss += (1.-2.*G4UniformRand())*e2Fluct;
1001       else if (loss>0.)                          1003       else if (loss>0.)
1002         loss += (1.-2.*G4UniformRand())*e1Flu    1004         loss += (1.-2.*G4UniformRand())*e1Fluct;   
1003                                                  1005       
1004       // ionisation .........................    1006       // ionisation .......................................
1005       if(a3 > 0.)                                1007       if(a3 > 0.)
1006   {                                              1008   {
1007     if(a3>alim)                                  1009     if(a3>alim)
1008       {                                          1010       {
1009         siga=std::sqrt(a3) ;                     1011         siga=std::sqrt(a3) ;
1010         p3 = std::max(0,int(G4RandGauss::shoo    1012         p3 = std::max(0,int(G4RandGauss::shoot(a3,siga)+0.5));
1011       }                                          1013       }
1012     else                                         1014     else
1013       p3 = G4Poisson(a3);                        1015       p3 = G4Poisson(a3);
1014                                                  1016     
1015     lossc = 0.;                                  1017     lossc = 0.;
1016     if(p3 > 0)                                   1018     if(p3 > 0)
1017       {                                          1019       {
1018         na = 0.;                                 1020         na = 0.; 
1019         alfa = 1.;                               1021         alfa = 1.;
1020         if (p3 > nmaxCont2)                      1022         if (p3 > nmaxCont2)
1021     {                                            1023     {
1022       dp3        = G4double(p3);                 1024       dp3        = G4double(p3);
1023       rfac       = dp3/(G4double(nmaxCont2)+d    1025       rfac       = dp3/(G4double(nmaxCont2)+dp3);
1024       namean     = G4double(p3)*rfac;            1026       namean     = G4double(p3)*rfac;
1025       sa         = G4double(nmaxCont1)*rfac;     1027       sa         = G4double(nmaxCont1)*rfac;
1026       na         = G4RandGauss::shoot(namean,    1028       na         = G4RandGauss::shoot(namean,sa);
1027       if (na > 0.)                               1029       if (na > 0.)
1028         {                                        1030         {
1029           alfa   = w1*G4double(nmaxCont2+p3)/    1031           alfa   = w1*G4double(nmaxCont2+p3)/(w1*G4double(nmaxCont2)+G4double(p3));
1030           alfa1  = alfa*std::log(alfa)/(alfa-    1032           alfa1  = alfa*std::log(alfa)/(alfa-1.);
1031           ea     = na*ipotFluct*alfa1;           1033           ea     = na*ipotFluct*alfa1;
1032           sea    = ipotFluct*std::sqrt(na*(al    1034           sea    = ipotFluct*std::sqrt(na*(alfa-alfa1*alfa1));
1033           lossc += G4RandGauss::shoot(ea,sea)    1035           lossc += G4RandGauss::shoot(ea,sea);
1034         }                                        1036         }
1035     }                                            1037     }
1036                                                  1038         
1037         nb = G4int(G4double(p3)-na);             1039         nb = G4int(G4double(p3)-na);
1038         if (nb > 0)                              1040         if (nb > 0)
1039     {                                            1041     {
1040       w2 = alfa*ipotFluct;                       1042       w2 = alfa*ipotFluct;
1041       w  = (Tm-w2)/Tm;                           1043       w  = (Tm-w2)/Tm;      
1042       for (G4int k=0; k<nb; k++) lossc += w2/    1044       for (G4int k=0; k<nb; k++) lossc += w2/(1.-w*G4UniformRand());
1043     }                                            1045     }
1044       }                                          1046       }        
1045                                                  1047     
1046     loss += lossc;                               1048     loss += lossc;  
1047   }                                              1049   }
1048     }                                            1050     } 
1049                                                  1051   
1050   return loss ;                                  1052   return loss ;
1051 }                                                1053 }
1052                                                  1054 
1053 //                                               1055 //    
1054                                                  1056    
1055                                                  1057