G4VeLowEnergyLoss.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id$
00028 //
00029 //
00030 // --------------------------------------------------------------
00031 //      GEANT 4 class implementation file
00032 //
00033 //      History: first implementation, based on object model of
00034 //      2nd December 1995, G.Cosmo
00035 // --------------------------------------------------------------
00036 //
00037 // Modifications:
00038 // 20/09/00 update fluctuations V.Ivanchenko
00039 // 22/11/00 minor fix in fluctuations V.Ivanchenko
00040 // 10/05/01  V.Ivanchenko Clean up againist Linux compilation with -Wall
00041 // 22/05/01  V.Ivanchenko Update range calculation
00042 // 23/11/01  V.Ivanchenko Move static member-functions from header to source
00043 // 22/01/03  V.Ivanchenko Cut per region
00044 // 11/02/03  V.Ivanchenko Add limits to fluctuations
00045 // 24/04/03  V.Ivanchenko Fix the problem of table size
00046 //
00047 // --------------------------------------------------------------
00048 
00049 #include "G4VeLowEnergyLoss.hh"
00050 #include "G4PhysicalConstants.hh"
00051 #include "G4SystemOfUnits.hh"
00052 #include "G4ProductionCutsTable.hh"
00053 
00054 G4double     G4VeLowEnergyLoss::ParticleMass ;
00055 G4double     G4VeLowEnergyLoss::taulow       ;
00056 G4double     G4VeLowEnergyLoss::tauhigh      ;
00057 G4double     G4VeLowEnergyLoss::ltaulow       ;
00058 G4double     G4VeLowEnergyLoss::ltauhigh      ;
00059 
00060 
00061 G4bool       G4VeLowEnergyLoss::rndmStepFlag   = false;
00062 G4bool       G4VeLowEnergyLoss::EnlossFlucFlag = true;
00063 G4double     G4VeLowEnergyLoss::dRoverRange    = 20*perCent;
00064 G4double     G4VeLowEnergyLoss::finalRange     = 200*micrometer;
00065 G4double     G4VeLowEnergyLoss::c1lim = dRoverRange ;
00066 G4double     G4VeLowEnergyLoss::c2lim = 2.*(1.-dRoverRange)*finalRange ;
00067 G4double     G4VeLowEnergyLoss::c3lim = -(1.-dRoverRange)*finalRange*finalRange;
00068 
00069 
00070 //
00071 
00072 G4VeLowEnergyLoss::G4VeLowEnergyLoss()
00073                    :G4VContinuousDiscreteProcess("No Name Loss Process"),
00074                     lastMaterial(0),
00075                     imat(-1),
00076                     f1Fluct(0),f2Fluct(0),e1Fluct(0),e2Fluct(0),rateFluct(0),
00077                     ipotFluct(0),e1LogFluct(-1),e2LogFluct(-1),ipotLogFluct(-1),                    
00078                     nmaxCont1(4),
00079                     nmaxCont2(16)
00080 {
00081       G4Exception("G4VeLowEnergyLoss::G4VeLowEnergyLoss()",
00082                     "em1009",FatalException,"default constructor is called");
00083 }
00084 
00085 //
00086 
00087 G4VeLowEnergyLoss::G4VeLowEnergyLoss(const G4String& aName, G4ProcessType aType)
00088                   : G4VContinuousDiscreteProcess(aName, aType),
00089                     lastMaterial(0),
00090                     imat(-1),
00091                     f1Fluct(0),f2Fluct(0),e1Fluct(0),e2Fluct(0),rateFluct(0),
00092                     ipotFluct(0),e1LogFluct(-1),e2LogFluct(-1),ipotLogFluct(-1),            
00093                     nmaxCont1(4),
00094                     nmaxCont2(16)
00095 {;}
00096 
00097 //
00098 
00099 G4VeLowEnergyLoss::~G4VeLowEnergyLoss()
00100 {
00101 }
00102 
00103 //
00104 
00105 G4VeLowEnergyLoss::G4VeLowEnergyLoss(G4VeLowEnergyLoss& right)
00106                   : G4VContinuousDiscreteProcess(right),
00107                     lastMaterial(0),
00108                     imat(-1),
00109                     f1Fluct(0),f2Fluct(0),e1Fluct(0),e2Fluct(0),rateFluct(0),
00110                     ipotFluct(0),e1LogFluct(-1),e2LogFluct(-1),ipotLogFluct(-1),
00111                     nmaxCont1(4),
00112                     nmaxCont2(16)
00113 {;}
00114 
00115 void G4VeLowEnergyLoss::SetRndmStep(G4bool value)
00116 {
00117    rndmStepFlag   = value;
00118 }
00119 
00120 void G4VeLowEnergyLoss::SetEnlossFluc(G4bool value)
00121 {
00122    EnlossFlucFlag = value;
00123 }
00124 
00125 void G4VeLowEnergyLoss::SetStepFunction (G4double c1, G4double c2)
00126 {
00127    dRoverRange = c1;
00128    finalRange = c2;
00129    c1lim=dRoverRange;
00130    c2lim=2.*(1-dRoverRange)*finalRange;
00131    c3lim=-(1.-dRoverRange)*finalRange*finalRange;
00132 }
00133 
00134 G4PhysicsTable* G4VeLowEnergyLoss::BuildRangeTable(G4PhysicsTable* theDEDXTable,
00135                                                    G4PhysicsTable* theRangeTable,
00136                                                    G4double lowestKineticEnergy,
00137                                                    G4double highestKineticEnergy,
00138                                                    G4int TotBin)
00139 // Build range table from the energy loss table
00140 {
00141 
00142    G4int numOfCouples = theDEDXTable->length();
00143 
00144    if(theRangeTable)
00145    { theRangeTable->clearAndDestroy();
00146      delete theRangeTable; }
00147    theRangeTable = new G4PhysicsTable(numOfCouples);
00148 
00149    // loop for materials
00150 
00151    for (G4int J=0;  J<numOfCouples; J++)
00152    {
00153      G4PhysicsLogVector* aVector;
00154      aVector = new G4PhysicsLogVector(lowestKineticEnergy,
00155                               highestKineticEnergy,TotBin);
00156      BuildRangeVector(theDEDXTable,lowestKineticEnergy,highestKineticEnergy,
00157                       TotBin,J,aVector);
00158      theRangeTable->insert(aVector);
00159    }
00160    return theRangeTable ;
00161 }
00162 
00163 //
00164 
00165 void G4VeLowEnergyLoss::BuildRangeVector(G4PhysicsTable* theDEDXTable,
00166                                          G4double lowestKineticEnergy,
00167                                          G4double,
00168                                          G4int TotBin,
00169                                          G4int materialIndex,
00170                                          G4PhysicsLogVector* rangeVector)
00171 
00172 //  create range vector for a material
00173 {
00174   G4bool isOut;
00175   G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex];
00176   G4double energy1 = lowestKineticEnergy;
00177   G4double dedx    = physicsVector->GetValue(energy1,isOut);
00178   G4double range   = 0.5*energy1/dedx;
00179   rangeVector->PutValue(0,range);
00180   G4int n = 100;
00181   G4double del = 1.0/(G4double)n ;
00182 
00183   for (G4int j=1; j<=TotBin; j++) {
00184 
00185     G4double energy2 = rangeVector->GetLowEdgeEnergy(j);
00186     G4double de = (energy2 - energy1) * del ;
00187     G4double dedx1 = dedx ;
00188 
00189     for (G4int i=1; i<n; i++) {
00190       G4double energy = energy1 + i*de ;
00191       G4double dedx2  = physicsVector->GetValue(energy,isOut);
00192       range  += 0.5*de*(1.0/dedx1 + 1.0/dedx2);
00193       dedx1   = dedx2;
00194     }
00195     rangeVector->PutValue(j,range);
00196     dedx = dedx1 ;
00197     energy1 = energy2 ;
00198   }
00199 }
00200 
00201 //
00202 
00203 G4double G4VeLowEnergyLoss::RangeIntLin(G4PhysicsVector* physicsVector,
00204                                         G4int nbin)
00205 //  num. integration, linear binning
00206 {
00207   G4double dtau,Value,taui,ti,lossi,ci;
00208   G4bool isOut;
00209   dtau = (tauhigh-taulow)/nbin;
00210   Value = 0.;
00211 
00212   for (G4int i=0; i<nbin; i++)
00213   {
00214     taui = taulow + dtau*i ;
00215     ti = ParticleMass*taui;
00216     lossi = physicsVector->GetValue(ti,isOut);
00217     if(i==0)
00218       ci=0.5;
00219     else
00220     {
00221       if(i<nbin-1)
00222         ci=1.;
00223       else
00224         ci=0.5;
00225     }
00226     Value += ci/lossi;
00227   }
00228   Value *= ParticleMass*dtau;
00229   return Value;
00230 }
00231 
00232 //
00233 
00234 G4double G4VeLowEnergyLoss::RangeIntLog(G4PhysicsVector* physicsVector,
00235                                         G4int nbin)
00236 //  num. integration, logarithmic binning
00237 {
00238   G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
00239   G4bool isOut;
00240   ltt = ltauhigh-ltaulow;
00241   dltau = ltt/nbin;
00242   Value = 0.;
00243 
00244   for (G4int i=0; i<nbin; i++)
00245   {
00246     ui = ltaulow+dltau*i;
00247     taui = std::exp(ui);
00248     ti = ParticleMass*taui;
00249     lossi = physicsVector->GetValue(ti,isOut);
00250     if(i==0)
00251       ci=0.5;
00252     else
00253     {
00254       if(i<nbin-1)
00255         ci=1.;
00256       else
00257         ci=0.5;
00258     }
00259     Value += ci*taui/lossi;
00260   }
00261   Value *= ParticleMass*dltau;
00262   return Value;
00263 }
00264 
00265 
00266 //
00267 
00268 G4PhysicsTable* G4VeLowEnergyLoss::BuildLabTimeTable(G4PhysicsTable* theDEDXTable,
00269                                                      G4PhysicsTable* theLabTimeTable,
00270                                                      G4double lowestKineticEnergy,
00271                                                      G4double highestKineticEnergy,G4int TotBin)
00272 
00273 {
00274 
00275   G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
00276 
00277   if(theLabTimeTable)
00278   { theLabTimeTable->clearAndDestroy();
00279     delete theLabTimeTable; }
00280   theLabTimeTable = new G4PhysicsTable(numOfCouples);
00281 
00282 
00283   for (G4int J=0;  J<numOfCouples; J++)
00284   {
00285     G4PhysicsLogVector* aVector;
00286 
00287     aVector = new G4PhysicsLogVector(lowestKineticEnergy,
00288                             highestKineticEnergy,TotBin);
00289 
00290     BuildLabTimeVector(theDEDXTable,
00291               lowestKineticEnergy,highestKineticEnergy,TotBin,J,aVector);
00292     theLabTimeTable->insert(aVector);
00293 
00294 
00295   }
00296   return theLabTimeTable ;
00297 }
00298 
00299 //    
00300 
00301 G4PhysicsTable* G4VeLowEnergyLoss::BuildProperTimeTable(G4PhysicsTable* theDEDXTable,
00302                                                         G4PhysicsTable* theProperTimeTable,
00303                                                         G4double lowestKineticEnergy,
00304                                                         G4double highestKineticEnergy,G4int TotBin)
00305                             
00306 {
00307 
00308   G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
00309  
00310   if(theProperTimeTable)
00311   { theProperTimeTable->clearAndDestroy();
00312     delete theProperTimeTable; }
00313   theProperTimeTable = new G4PhysicsTable(numOfCouples);
00314 
00315 
00316   for (G4int J=0;  J<numOfCouples; J++)
00317   {
00318     G4PhysicsLogVector* aVector;
00319 
00320     aVector = new G4PhysicsLogVector(lowestKineticEnergy,
00321                             highestKineticEnergy,TotBin);
00322 
00323     BuildProperTimeVector(theDEDXTable,
00324               lowestKineticEnergy,highestKineticEnergy,TotBin,J,aVector);
00325     theProperTimeTable->insert(aVector);
00326 
00327 
00328   }
00329   return theProperTimeTable ;
00330 }
00331 
00332 //    
00333 
00334 void G4VeLowEnergyLoss::BuildLabTimeVector(G4PhysicsTable* theDEDXTable,
00335                                            G4double, // lowestKineticEnergy
00336                                            G4double /*highestKineticEnergy*/, G4int TotBin,
00337                                            G4int materialIndex, G4PhysicsLogVector* timeVector)
00338 //  create lab time vector for a material
00339 {
00340 
00341   G4int nbin=100;
00342   G4bool isOut;
00343   G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-parlowen ;
00344   G4double losslim,clim,taulim,timelim,
00345            LowEdgeEnergy,tau,Value ;
00346   //G4double ltaulim,ltaumax; 
00347 
00348   G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex];
00349 
00350   // low energy part first...
00351   losslim = physicsVector->GetValue(tlim,isOut);
00352   taulim=tlim/ParticleMass ;
00353   clim=std::sqrt(ParticleMass*tlim/2.)/(c_light*losslim*ppar) ;
00354   //ltaulim = std::log(taulim);
00355   //ltaumax = std::log(highestKineticEnergy/ParticleMass) ;
00356 
00357   G4int i=-1;
00358   G4double oldValue = 0. ;
00359   G4double tauold ;
00360   do
00361   {
00362     i += 1 ;
00363     LowEdgeEnergy = timeVector->GetLowEdgeEnergy(i);
00364     tau = LowEdgeEnergy/ParticleMass ;
00365     if ( tau <= taulim )
00366     {
00367       Value = clim*std::exp(ppar*std::log(tau/taulim)) ;
00368     }
00369     else
00370     {
00371       timelim=clim ;
00372       ltaulow = std::log(taulim);
00373       ltauhigh = std::log(tau);
00374       Value = timelim+LabTimeIntLog(physicsVector,nbin);
00375     }
00376     timeVector->PutValue(i,Value);
00377     oldValue = Value ;
00378     tauold = tau ;
00379   } while (tau<=taulim) ;
00380   i += 1 ;
00381   for (G4int j=i; j<=TotBin; j++)
00382   {
00383     LowEdgeEnergy = timeVector->GetLowEdgeEnergy(j);
00384     tau = LowEdgeEnergy/ParticleMass ;
00385     ltaulow = std::log(tauold);
00386     ltauhigh = std::log(tau);
00387     Value = oldValue+LabTimeIntLog(physicsVector,nbin);
00388     timeVector->PutValue(j,Value);
00389     oldValue = Value ;
00390     tauold = tau ;
00391   }
00392 }
00393 
00394 //    
00395 
00396 void G4VeLowEnergyLoss::BuildProperTimeVector(G4PhysicsTable* theDEDXTable,
00397                                               G4double, // lowestKineticEnergy
00398                                               G4double /*highestKineticEnergy*/, G4int TotBin,
00399                                               G4int materialIndex, G4PhysicsLogVector* timeVector)
00400 //  create proper time vector for a material
00401 {
00402   G4int nbin=100;
00403   G4bool isOut;
00404   G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-parlowen ;
00405   G4double losslim,clim,taulim,timelim,
00406            LowEdgeEnergy,tau,Value ;
00407   //G4double ltaulim,ltaumax;
00408 
00409   G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex];
00410   //const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
00411 
00412   // low energy part first...
00413   losslim = physicsVector->GetValue(tlim,isOut);
00414   taulim=tlim/ParticleMass ;
00415   clim=std::sqrt(ParticleMass*tlim/2.)/(c_light*losslim*ppar) ;
00416   //ltaulim = std::log(taulim);
00417   //ltaumax = std::log(highestKineticEnergy/ParticleMass) ;
00418 
00419   G4int i=-1;
00420   G4double oldValue = 0. ;
00421   G4double tauold ;
00422   do
00423   {
00424     i += 1 ;
00425     LowEdgeEnergy = timeVector->GetLowEdgeEnergy(i);
00426     tau = LowEdgeEnergy/ParticleMass ;
00427     if ( tau <= taulim )
00428     {
00429       Value = clim*std::exp(ppar*std::log(tau/taulim)) ;
00430     }
00431     else
00432     {
00433       timelim=clim ;
00434       ltaulow = std::log(taulim);
00435       ltauhigh = std::log(tau);
00436       Value = timelim+ProperTimeIntLog(physicsVector,nbin);
00437     }
00438     timeVector->PutValue(i,Value);
00439     oldValue = Value ;
00440     tauold = tau ;
00441   } while (tau<=taulim) ;
00442   i += 1 ;
00443   for (G4int j=i; j<=TotBin; j++)
00444   {
00445     LowEdgeEnergy = timeVector->GetLowEdgeEnergy(j);
00446     tau = LowEdgeEnergy/ParticleMass ;
00447     ltaulow = std::log(tauold);
00448     ltauhigh = std::log(tau);
00449     Value = oldValue+ProperTimeIntLog(physicsVector,nbin);
00450     timeVector->PutValue(j,Value);
00451     oldValue = Value ;
00452     tauold = tau ;
00453   }
00454 }
00455 
00456 //    
00457 
00458 G4double G4VeLowEnergyLoss::LabTimeIntLog(G4PhysicsVector* physicsVector,
00459                                           G4int nbin)
00460 //  num. integration, logarithmic binning
00461 {
00462   G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
00463   G4bool isOut;
00464   ltt = ltauhigh-ltaulow;
00465   dltau = ltt/nbin;
00466   Value = 0.;
00467 
00468   for (G4int i=0; i<nbin; i++)
00469   {
00470     ui = ltaulow+dltau*i;
00471     taui = std::exp(ui);
00472     ti = ParticleMass*taui;
00473     lossi = physicsVector->GetValue(ti,isOut);
00474     if(i==0)
00475       ci=0.5;
00476     else
00477     {
00478       if(i<nbin-1)
00479         ci=1.;
00480       else
00481         ci=0.5;
00482     }
00483     Value += ci*taui*(ti+ParticleMass)/(std::sqrt(ti*(ti+2.*ParticleMass))*lossi);
00484   }
00485   Value *= ParticleMass*dltau/c_light;
00486   return Value;
00487 }
00488 
00489 //    
00490 
00491 G4double G4VeLowEnergyLoss::ProperTimeIntLog(G4PhysicsVector* physicsVector,
00492                                              G4int nbin)
00493 //  num. integration, logarithmic binning
00494 {
00495   G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
00496   G4bool isOut;
00497   ltt = ltauhigh-ltaulow;
00498   dltau = ltt/nbin;
00499   Value = 0.;
00500 
00501   for (G4int i=0; i<nbin; i++)
00502   {
00503     ui = ltaulow+dltau*i;
00504     taui = std::exp(ui);
00505     ti = ParticleMass*taui;
00506     lossi = physicsVector->GetValue(ti,isOut);
00507     if(i==0)
00508       ci=0.5;
00509     else
00510     {
00511       if(i<nbin-1)
00512         ci=1.;
00513       else
00514         ci=0.5;
00515     }
00516     Value += ci*taui*ParticleMass/(std::sqrt(ti*(ti+2.*ParticleMass))*lossi);
00517   }
00518   Value *= ParticleMass*dltau/c_light;
00519   return Value;
00520 }
00521 
00522 //    
00523 
00524 G4PhysicsTable* G4VeLowEnergyLoss::BuildInverseRangeTable(G4PhysicsTable* theRangeTable,
00525                                                           G4PhysicsTable*,
00526                                                           G4PhysicsTable*,
00527                                                           G4PhysicsTable*,
00528                                                           G4PhysicsTable* theInverseRangeTable,
00529                                                           G4double, // lowestKineticEnergy,
00530                                                           G4double, // highestKineticEnergy
00531                                                           G4int )   // nbins
00532 // Build inverse table of the range table
00533 {
00534   G4bool b;
00535 
00536   G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
00537 
00538     if(theInverseRangeTable)
00539     { theInverseRangeTable->clearAndDestroy();
00540       delete theInverseRangeTable; }
00541     theInverseRangeTable = new G4PhysicsTable(numOfCouples);
00542 
00543   // loop for materials
00544   for (G4int i=0;  i<numOfCouples; i++)
00545   {
00546 
00547     G4PhysicsVector* pv = (*theRangeTable)[i];
00548     size_t nbins   = pv->GetVectorLength();
00549     G4double elow  = pv->GetLowEdgeEnergy(0);
00550     G4double ehigh = pv->GetLowEdgeEnergy(nbins-1);
00551     G4double rlow  = pv->GetValue(elow, b);
00552     G4double rhigh = pv->GetValue(ehigh, b);
00553 
00554     //rhigh *= std::exp(std::log(rhigh/rlow)/((G4double)(nbins-1)));
00555 
00556     G4PhysicsLogVector* v = new G4PhysicsLogVector(rlow, rhigh, nbins-1);
00557 
00558     v->PutValue(0,elow);
00559     G4double energy1 = elow;
00560     G4double range1  = rlow;
00561     G4double energy2 = elow;
00562     G4double range2  = rlow;
00563     size_t ilow      = 0;
00564     size_t ihigh;
00565 
00566     for (size_t j=1; j<nbins; j++) {
00567 
00568       G4double range = v->GetLowEdgeEnergy(j);
00569 
00570       for (ihigh=ilow+1; ihigh<nbins; ihigh++) {
00571         energy2 = pv->GetLowEdgeEnergy(ihigh);
00572         range2  = pv->GetValue(energy2, b);
00573         if(range2 >= range || ihigh == nbins-1) {
00574           ilow = ihigh - 1;
00575           energy1 = pv->GetLowEdgeEnergy(ilow);
00576           range1  = pv->GetValue(energy1, b); 
00577           break;
00578         }
00579       }
00580 
00581       G4double e = std::log(energy1) + std::log(energy2/energy1)*std::log(range/range1)/std::log(range2/range1);
00582 
00583       v->PutValue(j,std::exp(e));
00584     }
00585     theInverseRangeTable->insert(v);
00586 
00587   }
00588   return theInverseRangeTable ;
00589 }
00590 
00591 //    
00592 
00593 void G4VeLowEnergyLoss::InvertRangeVector(G4PhysicsTable* theRangeTable,
00594                                           G4PhysicsTable* theRangeCoeffATable,
00595                                           G4PhysicsTable* theRangeCoeffBTable,
00596                                           G4PhysicsTable* theRangeCoeffCTable,
00597                                           G4double lowestKineticEnergy,
00598                                           G4double highestKineticEnergy, G4int TotBin,
00599                                           G4int  materialIndex, G4PhysicsLogVector* aVector)
00600 //  invert range vector for a material
00601 {
00602   G4double LowEdgeRange,A,B,C,discr,KineticEnergy ;
00603   G4double RTable = std::exp(std::log(highestKineticEnergy/lowestKineticEnergy)/TotBin) ;
00604   G4double Tbin = lowestKineticEnergy/RTable ;
00605   G4double rangebin = 0.0 ;
00606   G4int binnumber = -1 ;
00607   G4bool isOut ;
00608 
00609   //loop for range values
00610   for( G4int i=0; i<=TotBin; i++)
00611   {
00612     LowEdgeRange = aVector->GetLowEdgeEnergy(i) ;  //i.e. GetLowEdgeValue(i)
00613     if( rangebin < LowEdgeRange )
00614     {
00615       do
00616       {
00617         binnumber += 1 ;
00618         Tbin *= RTable ;
00619         rangebin = (*theRangeTable)(materialIndex)->GetValue(Tbin,isOut) ;
00620       }
00621       while ((rangebin < LowEdgeRange) && (binnumber < TotBin )) ;
00622     }
00623 
00624     if(binnumber == 0)
00625       KineticEnergy = lowestKineticEnergy ;
00626     else if(binnumber == TotBin)
00627       KineticEnergy = highestKineticEnergy ;
00628     else
00629     {
00630       A = (*(*theRangeCoeffATable)(materialIndex))(binnumber-1) ;
00631       B = (*(*theRangeCoeffBTable)(materialIndex))(binnumber-1) ;
00632       C = (*(*theRangeCoeffCTable)(materialIndex))(binnumber-1) ;
00633       if(A==0.)
00634          KineticEnergy = (LowEdgeRange -C )/B ;
00635       else
00636       {
00637          discr = B*B - 4.*A*(C-LowEdgeRange);
00638          discr = discr>0. ? std::sqrt(discr) : 0.;
00639          KineticEnergy = 0.5*(discr-B)/A ;
00640       }
00641     }
00642 
00643     aVector->PutValue(i,KineticEnergy) ;
00644   }
00645 }
00646 
00647 //    
00648 
00649 G4PhysicsTable* G4VeLowEnergyLoss::BuildRangeCoeffATable(G4PhysicsTable* theRangeTable,
00650                                                          G4PhysicsTable* theRangeCoeffATable,
00651                                                          G4double lowestKineticEnergy,
00652                                                          G4double highestKineticEnergy, G4int TotBin)
00653 // Build tables of coefficients for the energy loss calculation
00654 //  create table for coefficients "A"
00655 {
00656 
00657   G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
00658 
00659   if(theRangeCoeffATable)
00660   { theRangeCoeffATable->clearAndDestroy();
00661     delete theRangeCoeffATable; }
00662   theRangeCoeffATable = new G4PhysicsTable(numOfCouples);
00663 
00664   G4double RTable = std::exp(std::log(highestKineticEnergy/lowestKineticEnergy)/TotBin) ;
00665   G4double R2 = RTable*RTable ;
00666   G4double R1 = RTable+1.;
00667   G4double w = R1*(RTable-1.)*(RTable-1.);
00668   G4double w1 = RTable/w , w2 = -RTable*R1/w , w3 = R2/w ;
00669   G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
00670   G4bool isOut;
00671 
00672   //  loop for materials
00673   for (G4int J=0; J<numOfCouples; J++)
00674   {
00675     G4int binmax=TotBin ;
00676     G4PhysicsLinearVector* aVector =
00677                            new G4PhysicsLinearVector(0.,binmax, TotBin);
00678     Ti = lowestKineticEnergy ;
00679     G4PhysicsVector* rangeVector= (*theRangeTable)[J];
00680 
00681     for ( G4int i=0; i<=TotBin; i++)
00682     {
00683       Ri = rangeVector->GetValue(Ti,isOut) ;
00684       if ( i==0 )
00685         Rim = 0. ;
00686       else
00687       {
00688         Tim = Ti/RTable ;
00689         Rim = rangeVector->GetValue(Tim,isOut);
00690       }
00691       if ( i==TotBin)
00692         Rip = Ri ;
00693       else
00694       {
00695         Tip = Ti*RTable ;
00696         Rip = rangeVector->GetValue(Tip,isOut);
00697       }
00698       Value = (w1*Rip + w2*Ri + w3*Rim)/(Ti*Ti) ;
00699 
00700       aVector->PutValue(i,Value);
00701       Ti = RTable*Ti ;
00702     }
00703  
00704     theRangeCoeffATable->insert(aVector);
00705   }
00706   return theRangeCoeffATable ;
00707 }
00708 
00709 //    
00710   
00711 G4PhysicsTable* G4VeLowEnergyLoss::BuildRangeCoeffBTable(G4PhysicsTable* theRangeTable,
00712                                                          G4PhysicsTable* theRangeCoeffBTable,
00713                                                          G4double lowestKineticEnergy,
00714                                                          G4double highestKineticEnergy, G4int TotBin)
00715 // Build tables of coefficients for the energy loss calculation
00716 //  create table for coefficients "B"
00717 {
00718 
00719   G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
00720 
00721   if(theRangeCoeffBTable)
00722   { theRangeCoeffBTable->clearAndDestroy();
00723     delete theRangeCoeffBTable; }
00724   theRangeCoeffBTable = new G4PhysicsTable(numOfCouples);
00725 
00726   G4double RTable = std::exp(std::log(highestKineticEnergy/lowestKineticEnergy)/TotBin) ;
00727   G4double R2 = RTable*RTable ;
00728   G4double R1 = RTable+1.;
00729   G4double w = R1*(RTable-1.)*(RTable-1.);
00730   G4double w1 = -R1/w , w2 = R1*(R2+1.)/w , w3 = -R2*R1/w ;
00731   G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
00732   G4bool isOut;
00733 
00734   //  loop for materials
00735   for (G4int J=0; J<numOfCouples; J++)
00736   {
00737     G4int binmax=TotBin ;
00738     G4PhysicsLinearVector* aVector =
00739                         new G4PhysicsLinearVector(0.,binmax, TotBin);
00740     Ti = lowestKineticEnergy ;
00741     G4PhysicsVector* rangeVector= (*theRangeTable)[J];
00742   
00743     for ( G4int i=0; i<=TotBin; i++)
00744     {
00745       Ri = rangeVector->GetValue(Ti,isOut) ;
00746       if ( i==0 )
00747          Rim = 0. ;
00748       else
00749       {
00750         Tim = Ti/RTable ;
00751         Rim = rangeVector->GetValue(Tim,isOut);
00752       }
00753       if ( i==TotBin)
00754         Rip = Ri ;
00755       else
00756       {
00757         Tip = Ti*RTable ;
00758         Rip = rangeVector->GetValue(Tip,isOut);
00759       }
00760       Value = (w1*Rip + w2*Ri + w3*Rim)/Ti;
00761 
00762       aVector->PutValue(i,Value);
00763       Ti = RTable*Ti ;
00764     }
00765     theRangeCoeffBTable->insert(aVector);
00766   }
00767   return theRangeCoeffBTable ;
00768 }
00769 
00770 //    
00771 
00772 G4PhysicsTable* G4VeLowEnergyLoss::BuildRangeCoeffCTable(G4PhysicsTable* theRangeTable,
00773                                                          G4PhysicsTable* theRangeCoeffCTable,
00774                                                          G4double lowestKineticEnergy,
00775                                                          G4double highestKineticEnergy, G4int TotBin)
00776 // Build tables of coefficients for the energy loss calculation
00777 //  create table for coefficients "C"
00778 {
00779 
00780   G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
00781 
00782   if(theRangeCoeffCTable)
00783   { theRangeCoeffCTable->clearAndDestroy();
00784     delete theRangeCoeffCTable; }
00785   theRangeCoeffCTable = new G4PhysicsTable(numOfCouples);
00786 
00787   G4double RTable = std::exp(std::log(highestKineticEnergy/lowestKineticEnergy)/TotBin) ;
00788   G4double R2 = RTable*RTable ;
00789   G4double R1 = RTable+1.;
00790   G4double w = R1*(RTable-1.)*(RTable-1.);
00791   G4double w1 = 1./w , w2 = -RTable*R1/w , w3 = RTable*R2/w ;
00792   G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
00793   G4bool isOut;
00794 
00795   //  loop for materials
00796   for (G4int J=0; J<numOfCouples; J++)
00797   {
00798     G4int binmax=TotBin ;
00799     G4PhysicsLinearVector* aVector =
00800                       new G4PhysicsLinearVector(0.,binmax, TotBin);
00801     Ti = lowestKineticEnergy ;
00802     G4PhysicsVector* rangeVector= (*theRangeTable)[J];
00803   
00804     for ( G4int i=0; i<=TotBin; i++)
00805     {
00806       Ri = rangeVector->GetValue(Ti,isOut) ;
00807       if ( i==0 )
00808         Rim = 0. ;
00809       else
00810       {
00811         Tim = Ti/RTable ;
00812         Rim = rangeVector->GetValue(Tim,isOut);
00813       }
00814       if ( i==TotBin)
00815         Rip = Ri ;
00816       else
00817       {
00818         Tip = Ti*RTable ;
00819         Rip = rangeVector->GetValue(Tip,isOut);
00820       }
00821       Value = w1*Rip + w2*Ri + w3*Rim ;
00822 
00823       aVector->PutValue(i,Value);
00824       Ti = RTable*Ti ;
00825     }
00826     theRangeCoeffCTable->insert(aVector);
00827   }
00828   return theRangeCoeffCTable ;
00829 }
00830 
00831 //    
00832 
00833 G4double G4VeLowEnergyLoss::GetLossWithFluct(const G4DynamicParticle* aParticle,
00834                                              const G4MaterialCutsCouple* couple,
00835                                              G4double MeanLoss,
00836                                              G4double step)
00837 //  calculate actual loss from the mean loss
00838 //  The model used to get the fluctuation is essentially the same as in Glandz in Geant3.
00839 {
00840    static const G4double minLoss = 1.*eV ;
00841    static const G4double probLim = 0.01 ;
00842    static const G4double sumaLim = -std::log(probLim) ;
00843    static const G4double alim=10.;
00844    static const G4double kappa = 10. ;
00845    static const G4double factor = twopi_mc2_rcl2 ;
00846   const G4Material* aMaterial = couple->GetMaterial();
00847 
00848   // check if the material has changed ( cache mechanism)
00849 
00850   if (aMaterial != lastMaterial)
00851     {
00852       lastMaterial = aMaterial;
00853       imat         = couple->GetIndex();
00854       f1Fluct      = aMaterial->GetIonisation()->GetF1fluct();
00855       f2Fluct      = aMaterial->GetIonisation()->GetF2fluct();
00856       e1Fluct      = aMaterial->GetIonisation()->GetEnergy1fluct();
00857       e2Fluct      = aMaterial->GetIonisation()->GetEnergy2fluct();
00858       e1LogFluct   = aMaterial->GetIonisation()->GetLogEnergy1fluct();
00859       e2LogFluct   = aMaterial->GetIonisation()->GetLogEnergy2fluct();
00860       rateFluct    = aMaterial->GetIonisation()->GetRateionexcfluct();
00861       ipotFluct    = aMaterial->GetIonisation()->GetMeanExcitationEnergy();
00862       ipotLogFluct = aMaterial->GetIonisation()->GetLogMeanExcEnergy();
00863     }
00864   G4double threshold,w1,w2,C,
00865            beta2,suma,e0,loss,lossc,w;
00866   G4double a1,a2,a3;
00867   G4int p1,p2,p3;
00868   G4int nb;
00869   G4double Corrfac, na,alfa,rfac,namean,sa,alfa1,ea,sea;
00870   //  G4double dp1;
00871   G4double dp3;
00872   G4double siga ;
00873 
00874   // shortcut for very very small loss
00875   if(MeanLoss < minLoss) return MeanLoss ;
00876 
00877   // get particle data
00878   G4double Tkin   = aParticle->GetKineticEnergy();
00879 
00880   //  G4cout << "MGP -- Fluc Tkin " << Tkin/keV << " keV "  << " MeanLoss = " << MeanLoss/keV << G4endl;
00881 
00882   threshold = (*((G4ProductionCutsTable::GetProductionCutsTable())
00883                 ->GetEnergyCutsVector(1)))[imat];
00884   G4double rmass = electron_mass_c2/ParticleMass;
00885   G4double tau   = Tkin/ParticleMass, tau1 = tau+1., tau2 = tau*(tau+2.);
00886   G4double Tm    = 2.*electron_mass_c2*tau2/(1.+2.*tau1*rmass+rmass*rmass);
00887 
00888   // G4cout << "MGP Particle mass " << ParticleMass/MeV << " Tm " << Tm << G4endl;
00889 
00890   if(Tm > threshold) Tm = threshold;
00891   beta2 = tau2/(tau1*tau1);
00892 
00893   // Gaussian fluctuation ?
00894   if(MeanLoss >= kappa*Tm || MeanLoss <= kappa*ipotFluct)
00895   {
00896     G4double electronDensity = aMaterial->GetElectronDensity() ;
00897     siga = std::sqrt(Tm*(1.0-0.5*beta2)*step*
00898                 factor*electronDensity/beta2) ;
00899     do {
00900       loss = G4RandGauss::shoot(MeanLoss,siga) ;
00901     } while (loss < 0. || loss > 2.0*MeanLoss);
00902     return loss ;
00903   }
00904 
00905   w1 = Tm/ipotFluct;
00906   w2 = std::log(2.*electron_mass_c2*tau2);
00907 
00908   C = MeanLoss*(1.-rateFluct)/(w2-ipotLogFluct-beta2);
00909 
00910   a1 = C*f1Fluct*(w2-e1LogFluct-beta2)/e1Fluct;
00911   a2 = C*f2Fluct*(w2-e2LogFluct-beta2)/e2Fluct;
00912   a3 = rateFluct*MeanLoss*(Tm-ipotFluct)/(ipotFluct*Tm*std::log(w1));
00913 
00914   suma = a1+a2+a3;
00915 
00916   loss = 0. ;
00917 
00918   if(suma < sumaLim)             // very small Step
00919     {
00920       e0 = aMaterial->GetIonisation()->GetEnergy0fluct();
00921       // G4cout << "MGP e0 = " << e0/keV << G4endl;
00922 
00923       if(Tm == ipotFluct)
00924         {
00925           a3 = MeanLoss/e0;
00926 
00927           if(a3>alim)
00928             {
00929               siga=std::sqrt(a3) ;
00930               p3 = std::max(0,G4int(G4RandGauss::shoot(a3,siga)+0.5));
00931             }
00932           else p3 = G4Poisson(a3);
00933 
00934           loss = p3*e0 ;
00935 
00936           if(p3 > 0) loss += (1.-2.*G4UniformRand())*e0 ;
00937           // G4cout << "MGP very small step " << loss/keV << G4endl;
00938         }
00939       else
00940         {
00941           //      G4cout << "MGP old Tm = " << Tm << " " << ipotFluct << " " << e0 << G4endl;
00942           Tm = Tm-ipotFluct+e0 ;
00943 
00944           // MGP ---- workaround to avoid log argument<0, TO BE CHECKED
00945           if (Tm <= 0.)
00946             {
00947               loss = MeanLoss;
00948               p3 = 0;
00949               // G4cout << "MGP correction loss = MeanLoss " << loss/keV << G4endl;
00950             }
00951           else
00952             {
00953               a3 = MeanLoss*(Tm-e0)/(Tm*e0*std::log(Tm/e0));
00954 
00955               // G4cout << "MGP new Tm = " << Tm << " " << ipotFluct << " " << e0 << " a3= " << a3 << G4endl;
00956               
00957               if(a3>alim)
00958                 {
00959                   siga=std::sqrt(a3) ;
00960                   p3 = std::max(0,G4int(G4RandGauss::shoot(a3,siga)+0.5));
00961                 }
00962               else
00963                 p3 = G4Poisson(a3);
00964               //G4cout << "MGP p3 " << p3 << G4endl;
00965 
00966             }
00967               
00968           if(p3 > 0)
00969             {
00970               w = (Tm-e0)/Tm ;
00971               if(p3 > nmaxCont2)
00972                 {
00973                   // G4cout << "MGP dp3 " << dp3 << " p3 " << p3 << " " << nmaxCont2 << G4endl;
00974                   dp3 = G4double(p3) ;
00975                   Corrfac = dp3/G4double(nmaxCont2) ;
00976                   p3 = nmaxCont2 ;
00977                 }
00978               else
00979                 Corrfac = 1. ;
00980               
00981               for(G4int i=0; i<p3; i++) loss += 1./(1.-w*G4UniformRand()) ;
00982               loss *= e0*Corrfac ; 
00983               // G4cout << "MGP Corrfac = " << Corrfac << " e0 = " << e0/keV << " loss = " << loss/keV << G4endl;
00984             }        
00985         }
00986     }
00987 
00988   else                              // not so small Step
00989     {
00990       // excitation type 1
00991       if(a1>alim)
00992       {
00993         siga=std::sqrt(a1) ;
00994         p1 = std::max(0,int(G4RandGauss::shoot(a1,siga)+0.5));
00995       }
00996       else
00997        p1 = G4Poisson(a1);
00998 
00999       // excitation type 2
01000       if(a2>alim)
01001       {
01002         siga=std::sqrt(a2) ;
01003         p2 = std::max(0,int(G4RandGauss::shoot(a2,siga)+0.5));
01004       }
01005       else
01006         p2 = G4Poisson(a2);
01007 
01008       loss = p1*e1Fluct+p2*e2Fluct;
01009  
01010       // smearing to avoid unphysical peaks
01011       if(p2 > 0)
01012         loss += (1.-2.*G4UniformRand())*e2Fluct;
01013       else if (loss>0.)
01014         loss += (1.-2.*G4UniformRand())*e1Fluct;   
01015       
01016       // ionisation .......................................
01017       if(a3 > 0.)
01018         {
01019           if(a3>alim)
01020             {
01021               siga=std::sqrt(a3) ;
01022               p3 = std::max(0,int(G4RandGauss::shoot(a3,siga)+0.5));
01023             }
01024           else
01025             p3 = G4Poisson(a3);
01026           
01027           lossc = 0.;
01028           if(p3 > 0)
01029             {
01030               na = 0.; 
01031               alfa = 1.;
01032               if (p3 > nmaxCont2)
01033                 {
01034                   dp3        = G4double(p3);
01035                   rfac       = dp3/(G4double(nmaxCont2)+dp3);
01036                   namean     = G4double(p3)*rfac;
01037                   sa         = G4double(nmaxCont1)*rfac;
01038                   na         = G4RandGauss::shoot(namean,sa);
01039                   if (na > 0.)
01040                     {
01041                       alfa   = w1*G4double(nmaxCont2+p3)/(w1*G4double(nmaxCont2)+G4double(p3));
01042                       alfa1  = alfa*std::log(alfa)/(alfa-1.);
01043                       ea     = na*ipotFluct*alfa1;
01044                       sea    = ipotFluct*std::sqrt(na*(alfa-alfa1*alfa1));
01045                       lossc += G4RandGauss::shoot(ea,sea);
01046                     }
01047                 }
01048               
01049               nb = G4int(G4double(p3)-na);
01050               if (nb > 0)
01051                 {
01052                   w2 = alfa*ipotFluct;
01053                   w  = (Tm-w2)/Tm;      
01054                   for (G4int k=0; k<nb; k++) lossc += w2/(1.-w*G4UniformRand());
01055                 }
01056             }        
01057           
01058           loss += lossc;  
01059         }
01060     } 
01061   
01062   return loss ;
01063 }
01064 
01065 //    
01066    

Generated on Mon May 27 17:50:11 2013 for Geant4 by  doxygen 1.4.7