Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-31 09:22:09 

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 //
0028 //
0029 // --------------------------------------------------------------
0030 //  GEANT 4 class implementation file
0031 //
0032 //  History: first implementation, based on object model of
0033 //  2nd December 1995, G.Cosmo
0034 // --------------------------------------------------------------
0035 //
0036 // Modifications:
0037 // 20/09/00 update fluctuations V.Ivanchenko
0038 // 22/11/00 minor fix in fluctuations V.Ivanchenko
0039 // 10/05/01  V.Ivanchenko Clean up againist Linux compilation with -Wall
0040 // 22/05/01  V.Ivanchenko Update range calculation
0041 // 23/11/01  V.Ivanchenko Move static member-functions from header to source
0042 // 22/01/03  V.Ivanchenko Cut per region
0043 // 11/02/03  V.Ivanchenko Add limits to fluctuations
0044 // 24/04/03  V.Ivanchenko Fix the problem of table size
0045 //
0046 // --------------------------------------------------------------
0047 
0048 #include "G4RDVeLowEnergyLoss.hh"
0049 #include "G4PhysicalConstants.hh"
0050 #include "G4SystemOfUnits.hh"
0051 #include "G4ProductionCutsTable.hh"
0052 
0053 G4double     G4RDVeLowEnergyLoss::ParticleMass ;
0054 G4double     G4RDVeLowEnergyLoss::taulow       ;
0055 G4double     G4RDVeLowEnergyLoss::tauhigh      ;
0056 G4double     G4RDVeLowEnergyLoss::ltaulow       ;
0057 G4double     G4RDVeLowEnergyLoss::ltauhigh      ;
0058 
0059 
0060 G4bool       G4RDVeLowEnergyLoss::rndmStepFlag   = false;
0061 G4bool       G4RDVeLowEnergyLoss::EnlossFlucFlag = true;
0062 G4double     G4RDVeLowEnergyLoss::dRoverRange    = 20*perCent;
0063 G4double     G4RDVeLowEnergyLoss::finalRange     = 200*micrometer;
0064 G4double     G4RDVeLowEnergyLoss::c1lim = dRoverRange ;
0065 G4double     G4RDVeLowEnergyLoss::c2lim = 2.*(1.-dRoverRange)*finalRange ;
0066 G4double     G4RDVeLowEnergyLoss::c3lim = -(1.-dRoverRange)*finalRange*finalRange;
0067 
0068 
0069 //
0070 
0071 G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss()
0072                    :G4VContinuousDiscreteProcess("No Name Loss Process"),
0073      lastMaterial(0),
0074      nmaxCont1(4),
0075      nmaxCont2(16)
0076 {
0077   G4Exception("G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss()", "InvalidCall",
0078               FatalException, "Default constructor is called!");
0079 }
0080 
0081 //
0082 
0083 G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss(const G4String& aName, G4ProcessType aType)
0084                   : G4VContinuousDiscreteProcess(aName, aType),
0085      lastMaterial(0),
0086      nmaxCont1(4),
0087      nmaxCont2(16)
0088 {
0089 }
0090 
0091 //
0092 
0093 G4RDVeLowEnergyLoss::~G4RDVeLowEnergyLoss()
0094 {
0095 }
0096 
0097 //
0098 
0099 G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss(G4RDVeLowEnergyLoss& right)
0100                   : G4VContinuousDiscreteProcess(right),
0101      lastMaterial(0),
0102      nmaxCont1(4),
0103      nmaxCont2(16)
0104 {
0105 }
0106 
0107 void G4RDVeLowEnergyLoss::SetRndmStep(G4bool value)
0108 {
0109    rndmStepFlag   = value;
0110 }
0111 
0112 void G4RDVeLowEnergyLoss::SetEnlossFluc(G4bool value)
0113 {
0114    EnlossFlucFlag = value;
0115 }
0116 
0117 void G4RDVeLowEnergyLoss::SetStepFunction (G4double c1, G4double c2)
0118 {
0119    dRoverRange = c1;
0120    finalRange = c2;
0121    c1lim=dRoverRange;
0122    c2lim=2.*(1-dRoverRange)*finalRange;
0123    c3lim=-(1.-dRoverRange)*finalRange*finalRange;
0124 }
0125 
0126 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRangeTable(G4PhysicsTable* theDEDXTable,
0127                            G4PhysicsTable* theRangeTable,
0128                            G4double lowestKineticEnergy,
0129                            G4double highestKineticEnergy,
0130                            G4int TotBin)
0131 // Build range table from the energy loss table
0132 {
0133 
0134    G4int numOfCouples = theDEDXTable->length();
0135 
0136    if(theRangeTable)
0137    { theRangeTable->clearAndDestroy();
0138      delete theRangeTable; }
0139    theRangeTable = new G4PhysicsTable(numOfCouples);
0140 
0141    // loop for materials
0142 
0143    for (G4int J=0;  J<numOfCouples; J++)
0144    {
0145      G4PhysicsLogVector* aVector;
0146      aVector = new G4PhysicsLogVector(lowestKineticEnergy,
0147                               highestKineticEnergy,TotBin);
0148      BuildRangeVector(theDEDXTable,lowestKineticEnergy,highestKineticEnergy,
0149                       TotBin,J,aVector);
0150      theRangeTable->insert(aVector);
0151    }
0152    return theRangeTable ;
0153 }
0154 
0155 //
0156 
0157 void G4RDVeLowEnergyLoss::BuildRangeVector(G4PhysicsTable* theDEDXTable,
0158                                          G4double lowestKineticEnergy,
0159                                          G4double,
0160                                          G4int TotBin,
0161                                          G4int materialIndex,
0162                                          G4PhysicsLogVector* rangeVector)
0163 
0164 //  create range vector for a material
0165 {
0166   G4bool isOut;
0167   G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex];
0168   G4double energy1 = lowestKineticEnergy;
0169   G4double dedx    = physicsVector->GetValue(energy1,isOut);
0170   G4double range   = 0.5*energy1/dedx;
0171   rangeVector->PutValue(0,range);
0172   G4int n = 100;
0173   G4double del = 1.0/(G4double)n ;
0174 
0175   for (G4int j=1; j<TotBin; j++) {
0176 
0177     G4double energy2 = rangeVector->GetLowEdgeEnergy(j);
0178     G4double de = (energy2 - energy1) * del ;
0179     G4double dedx1 = dedx ;
0180 
0181     for (G4int i=1; i<n; i++) {
0182       G4double energy = energy1 + i*de ;
0183       G4double dedx2  = physicsVector->GetValue(energy,isOut);
0184       range  += 0.5*de*(1.0/dedx1 + 1.0/dedx2);
0185       dedx1   = dedx2;
0186     }
0187     rangeVector->PutValue(j,range);
0188     dedx = dedx1 ;
0189     energy1 = energy2 ;
0190   }
0191 }
0192 
0193 //
0194 
0195 G4double G4RDVeLowEnergyLoss::RangeIntLin(G4PhysicsVector* physicsVector,
0196                     G4int nbin)
0197 //  num. integration, linear binning
0198 {
0199   G4double dtau,Value,taui,ti,lossi,ci;
0200   G4bool isOut;
0201   dtau = (tauhigh-taulow)/nbin;
0202   Value = 0.;
0203 
0204   for (G4int i=0; i<=nbin; i++)
0205   {
0206     taui = taulow + dtau*i ;
0207     ti = ParticleMass*taui;
0208     lossi = physicsVector->GetValue(ti,isOut);
0209     if(i==0)
0210       ci=0.5;
0211     else
0212     {
0213       if(i<nbin)
0214         ci=1.;
0215       else
0216         ci=0.5;
0217     }
0218     Value += ci/lossi;
0219   }
0220   Value *= ParticleMass*dtau;
0221   return Value;
0222 }
0223 
0224 //
0225 
0226 G4double G4RDVeLowEnergyLoss::RangeIntLog(G4PhysicsVector* physicsVector,
0227                     G4int nbin)
0228 //  num. integration, logarithmic binning
0229 {
0230   G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
0231   G4bool isOut;
0232   ltt = ltauhigh-ltaulow;
0233   dltau = ltt/nbin;
0234   Value = 0.;
0235 
0236   for (G4int i=0; i<=nbin; i++)
0237   {
0238     ui = ltaulow+dltau*i;
0239     taui = std::exp(ui);
0240     ti = ParticleMass*taui;
0241     lossi = physicsVector->GetValue(ti,isOut);
0242     if(i==0)
0243       ci=0.5;
0244     else
0245     {
0246       if(i<nbin)
0247         ci=1.;
0248       else
0249         ci=0.5;
0250     }
0251     Value += ci*taui/lossi;
0252   }
0253   Value *= ParticleMass*dltau;
0254   return Value;
0255 }
0256 
0257 
0258 //
0259 
0260 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildLabTimeTable(G4PhysicsTable* theDEDXTable,
0261                              G4PhysicsTable* theLabTimeTable,
0262                              G4double lowestKineticEnergy,
0263                              G4double highestKineticEnergy,G4int TotBin)
0264 
0265 {
0266 
0267   G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
0268 
0269   if(theLabTimeTable)
0270   { theLabTimeTable->clearAndDestroy();
0271     delete theLabTimeTable; }
0272   theLabTimeTable = new G4PhysicsTable(numOfCouples);
0273 
0274 
0275   for (G4int J=0;  J<numOfCouples; J++)
0276   {
0277     G4PhysicsLogVector* aVector;
0278 
0279     aVector = new G4PhysicsLogVector(lowestKineticEnergy,
0280                             highestKineticEnergy,TotBin);
0281 
0282     BuildLabTimeVector(theDEDXTable,
0283               lowestKineticEnergy,highestKineticEnergy,TotBin,J,aVector);
0284     theLabTimeTable->insert(aVector);
0285 
0286 
0287   }
0288   return theLabTimeTable ;
0289 }
0290 
0291 //    
0292 
0293 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildProperTimeTable(G4PhysicsTable* theDEDXTable,
0294                             G4PhysicsTable* theProperTimeTable,
0295                             G4double lowestKineticEnergy,
0296                             G4double highestKineticEnergy,G4int TotBin)
0297                             
0298 {
0299 
0300   G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
0301  
0302   if(theProperTimeTable)
0303   { theProperTimeTable->clearAndDestroy();
0304     delete theProperTimeTable; }
0305   theProperTimeTable = new G4PhysicsTable(numOfCouples);
0306 
0307 
0308   for (G4int J=0;  J<numOfCouples; J++)
0309   {
0310     G4PhysicsLogVector* aVector;
0311 
0312     aVector = new G4PhysicsLogVector(lowestKineticEnergy,
0313                             highestKineticEnergy,TotBin);
0314 
0315     BuildProperTimeVector(theDEDXTable,
0316               lowestKineticEnergy,highestKineticEnergy,TotBin,J,aVector);
0317     theProperTimeTable->insert(aVector);
0318 
0319 
0320   }
0321   return theProperTimeTable ;
0322 }
0323 
0324 //    
0325 
0326 void G4RDVeLowEnergyLoss::BuildLabTimeVector(G4PhysicsTable* theDEDXTable,
0327                        G4double, // lowestKineticEnergy,
0328                        G4double, // highestKineticEnergy,
0329                                            G4int TotBin,
0330                        G4int materialIndex, G4PhysicsLogVector* timeVector)
0331 //  create lab time vector for a material
0332 {
0333 
0334   G4int nbin=100;
0335   G4bool isOut;
0336   G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-parlowen ;
0337   G4double losslim,clim,taulim,timelim,
0338            LowEdgeEnergy,tau,Value ;
0339 
0340   G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex];
0341 
0342   // low energy part first...
0343   losslim = physicsVector->GetValue(tlim,isOut);
0344   taulim=tlim/ParticleMass ;
0345   clim=std::sqrt(ParticleMass*tlim/2.)/(c_light*losslim*ppar) ;
0346 
0347   G4int i=-1;
0348   G4double oldValue = 0. ;
0349   G4double tauold ;
0350   do
0351   {
0352     i += 1 ;
0353     LowEdgeEnergy = timeVector->GetLowEdgeEnergy(i);
0354     tau = LowEdgeEnergy/ParticleMass ;
0355     if ( tau <= taulim )
0356     {
0357       Value = clim*std::exp(ppar*std::log(tau/taulim)) ;
0358     }
0359     else
0360     {
0361       timelim=clim ;
0362       ltaulow = std::log(taulim);
0363       ltauhigh = std::log(tau);
0364       Value = timelim+LabTimeIntLog(physicsVector,nbin);
0365     }
0366     timeVector->PutValue(i,Value);
0367     oldValue = Value ;
0368     tauold = tau ;
0369   } while (tau<=taulim) ;
0370   i += 1 ;
0371   for (G4int j=i; j<TotBin; j++)
0372   {
0373     LowEdgeEnergy = timeVector->GetLowEdgeEnergy(j);
0374     tau = LowEdgeEnergy/ParticleMass ;
0375     ltaulow = std::log(tauold);
0376     ltauhigh = std::log(tau);
0377     Value = oldValue+LabTimeIntLog(physicsVector,nbin);
0378     timeVector->PutValue(j,Value);
0379     oldValue = Value ;
0380     tauold = tau ;
0381   }
0382 }
0383 
0384 //    
0385 
0386 void G4RDVeLowEnergyLoss::BuildProperTimeVector(G4PhysicsTable* theDEDXTable,
0387                           G4double, // lowestKineticEnergy,
0388                           G4double, // highestKineticEnergy,
0389                                               G4int TotBin,
0390                           G4int materialIndex, G4PhysicsLogVector* timeVector)
0391 //  create proper time vector for a material
0392 {
0393   G4int nbin=100;
0394   G4bool isOut;
0395   G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-parlowen ;
0396   G4double losslim,clim,taulim,timelim,
0397            LowEdgeEnergy,tau,Value ;
0398 
0399   G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex];
0400   //const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
0401 
0402   // low energy part first...
0403   losslim = physicsVector->GetValue(tlim,isOut);
0404   taulim=tlim/ParticleMass ;
0405   clim=std::sqrt(ParticleMass*tlim/2.)/(c_light*losslim*ppar) ;
0406 
0407   G4int i=-1;
0408   G4double oldValue = 0. ;
0409   G4double tauold ;
0410   do
0411   {
0412     i += 1 ;
0413     LowEdgeEnergy = timeVector->GetLowEdgeEnergy(i);
0414     tau = LowEdgeEnergy/ParticleMass ;
0415     if ( tau <= taulim )
0416     {
0417       Value = clim*std::exp(ppar*std::log(tau/taulim)) ;
0418     }
0419     else
0420     {
0421       timelim=clim ;
0422       ltaulow = std::log(taulim);
0423       ltauhigh = std::log(tau);
0424       Value = timelim+ProperTimeIntLog(physicsVector,nbin);
0425     }
0426     timeVector->PutValue(i,Value);
0427     oldValue = Value ;
0428     tauold = tau ;
0429   } while (tau<=taulim) ;
0430   i += 1 ;
0431   for (G4int j=i; j<TotBin; j++)
0432   {
0433     LowEdgeEnergy = timeVector->GetLowEdgeEnergy(j);
0434     tau = LowEdgeEnergy/ParticleMass ;
0435     ltaulow = std::log(tauold);
0436     ltauhigh = std::log(tau);
0437     Value = oldValue+ProperTimeIntLog(physicsVector,nbin);
0438     timeVector->PutValue(j,Value);
0439     oldValue = Value ;
0440     tauold = tau ;
0441   }
0442 }
0443 
0444 //    
0445 
0446 G4double G4RDVeLowEnergyLoss::LabTimeIntLog(G4PhysicsVector* physicsVector,
0447                       G4int nbin)
0448 //  num. integration, logarithmic binning
0449 {
0450   G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
0451   G4bool isOut;
0452   ltt = ltauhigh-ltaulow;
0453   dltau = ltt/nbin;
0454   Value = 0.;
0455 
0456   for (G4int i=0; i<=nbin; i++)
0457   {
0458     ui = ltaulow+dltau*i;
0459     taui = std::exp(ui);
0460     ti = ParticleMass*taui;
0461     lossi = physicsVector->GetValue(ti,isOut);
0462     if(i==0)
0463       ci=0.5;
0464     else
0465     {
0466       if(i<nbin)
0467         ci=1.;
0468       else
0469         ci=0.5;
0470     }
0471     Value += ci*taui*(ti+ParticleMass)/(std::sqrt(ti*(ti+2.*ParticleMass))*lossi);
0472   }
0473   Value *= ParticleMass*dltau/c_light;
0474   return Value;
0475 }
0476 
0477 //    
0478 
0479 G4double G4RDVeLowEnergyLoss::ProperTimeIntLog(G4PhysicsVector* physicsVector,
0480                          G4int nbin)
0481 //  num. integration, logarithmic binning
0482 {
0483   G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
0484   G4bool isOut;
0485   ltt = ltauhigh-ltaulow;
0486   dltau = ltt/nbin;
0487   Value = 0.;
0488 
0489   for (G4int i=0; i<=nbin; i++)
0490   {
0491     ui = ltaulow+dltau*i;
0492     taui = std::exp(ui);
0493     ti = ParticleMass*taui;
0494     lossi = physicsVector->GetValue(ti,isOut);
0495     if(i==0)
0496       ci=0.5;
0497     else
0498     {
0499       if(i<nbin)
0500         ci=1.;
0501       else
0502         ci=0.5;
0503     }
0504     Value += ci*taui*ParticleMass/(std::sqrt(ti*(ti+2.*ParticleMass))*lossi);
0505   }
0506   Value *= ParticleMass*dltau/c_light;
0507   return Value;
0508 }
0509 
0510 //    
0511 
0512 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildInverseRangeTable(G4PhysicsTable* theRangeTable,
0513                               G4PhysicsTable*,
0514                               G4PhysicsTable*,
0515                               G4PhysicsTable*,
0516                               G4PhysicsTable* theInverseRangeTable,
0517                               G4double, // lowestKineticEnergy,
0518                               G4double, // highestKineticEnergy
0519                               G4int )   // nbins
0520 // Build inverse table of the range table
0521 {
0522   G4bool b;
0523 
0524   G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
0525 
0526     if(theInverseRangeTable)
0527     { theInverseRangeTable->clearAndDestroy();
0528       delete theInverseRangeTable; }
0529     theInverseRangeTable = new G4PhysicsTable(numOfCouples);
0530 
0531   // loop for materials
0532   for (G4int i=0;  i<numOfCouples; i++)
0533   {
0534 
0535     G4PhysicsVector* pv = (*theRangeTable)[i];
0536     size_t nbins   = pv->GetVectorLength();
0537     G4double elow  = pv->GetLowEdgeEnergy(0);
0538     G4double ehigh = pv->GetLowEdgeEnergy(nbins-1);
0539     G4double rlow  = pv->GetValue(elow, b);
0540     G4double rhigh = pv->GetValue(ehigh, b);
0541 
0542     rhigh *= std::exp(std::log(rhigh/rlow)/((G4double)(nbins-1)));
0543 
0544     G4PhysicsLogVector* v = new G4PhysicsLogVector(rlow, rhigh, nbins);
0545 
0546     v->PutValue(0,elow);
0547     G4double energy1 = elow;
0548     G4double range1  = rlow;
0549     G4double energy2 = elow;
0550     G4double range2  = rlow;
0551     size_t ilow      = 0;
0552     size_t ihigh;
0553 
0554     for (size_t j=1; j<nbins; j++) {
0555 
0556       G4double range = v->GetLowEdgeEnergy(j);
0557 
0558       for (ihigh=ilow+1; ihigh<nbins; ihigh++) {
0559         energy2 = pv->GetLowEdgeEnergy(ihigh);
0560         range2  = pv->GetValue(energy2, b);
0561         if(range2 >= range || ihigh == nbins-1) {
0562           ilow = ihigh - 1;
0563           energy1 = pv->GetLowEdgeEnergy(ilow);
0564           range1  = pv->GetValue(energy1, b); 
0565           break;
0566     }
0567       }
0568 
0569       G4double e = std::log(energy1) + std::log(energy2/energy1)*std::log(range/range1)/std::log(range2/range1);
0570 
0571       v->PutValue(j,std::exp(e));
0572     }
0573     theInverseRangeTable->insert(v);
0574 
0575   }
0576   return theInverseRangeTable ;
0577 }
0578 
0579 //    
0580 
0581 void G4RDVeLowEnergyLoss::InvertRangeVector(G4PhysicsTable* theRangeTable,
0582                       G4PhysicsTable* theRangeCoeffATable,
0583                       G4PhysicsTable* theRangeCoeffBTable,
0584                       G4PhysicsTable* theRangeCoeffCTable,
0585                       G4double lowestKineticEnergy,
0586                       G4double highestKineticEnergy, G4int TotBin,
0587                       G4int  materialIndex, G4PhysicsLogVector* aVector)
0588 //  invert range vector for a material
0589 {
0590   G4double LowEdgeRange,A,B,C,discr,KineticEnergy ;
0591   G4double RTable = std::exp(std::log(highestKineticEnergy/lowestKineticEnergy)/TotBin) ;
0592   G4double Tbin = lowestKineticEnergy/RTable ;
0593   G4double rangebin = 0.0 ;
0594   G4int binnumber = -1 ;
0595   G4bool isOut ;
0596 
0597   //loop for range values
0598   for( G4int i=0; i<TotBin; i++)
0599   {
0600     LowEdgeRange = aVector->GetLowEdgeEnergy(i) ;  //i.e. GetLowEdgeValue(i)
0601     if( rangebin < LowEdgeRange )
0602     {
0603       do
0604       {
0605         binnumber += 1 ;
0606         Tbin *= RTable ;
0607         rangebin = (*theRangeTable)(materialIndex)->GetValue(Tbin,isOut) ;
0608       }
0609       while ((rangebin < LowEdgeRange) && (binnumber < TotBin )) ;
0610     }
0611 
0612     if(binnumber == 0)
0613       KineticEnergy = lowestKineticEnergy ;
0614     else if(binnumber == TotBin-1)
0615       KineticEnergy = highestKineticEnergy ;
0616     else
0617     {
0618       A = (*(*theRangeCoeffATable)(materialIndex))(binnumber-1) ;
0619       B = (*(*theRangeCoeffBTable)(materialIndex))(binnumber-1) ;
0620       C = (*(*theRangeCoeffCTable)(materialIndex))(binnumber-1) ;
0621       if(A==0.)
0622          KineticEnergy = (LowEdgeRange -C )/B ;
0623       else
0624       {
0625          discr = B*B - 4.*A*(C-LowEdgeRange);
0626          discr = discr>0. ? std::sqrt(discr) : 0.;
0627          KineticEnergy = 0.5*(discr-B)/A ;
0628       }
0629     }
0630 
0631     aVector->PutValue(i,KineticEnergy) ;
0632   }
0633 }
0634 
0635 //    
0636 
0637 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRangeCoeffATable(G4PhysicsTable* theRangeTable,
0638                              G4PhysicsTable* theRangeCoeffATable,
0639                              G4double lowestKineticEnergy,
0640                              G4double highestKineticEnergy, G4int TotBin)
0641 // Build tables of coefficients for the energy loss calculation
0642 //  create table for coefficients "A"
0643 {
0644 
0645   G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
0646 
0647   if(theRangeCoeffATable)
0648   { theRangeCoeffATable->clearAndDestroy();
0649     delete theRangeCoeffATable; }
0650   theRangeCoeffATable = new G4PhysicsTable(numOfCouples);
0651 
0652   G4double RTable = std::exp(std::log(highestKineticEnergy/lowestKineticEnergy)/TotBin) ;
0653   G4double R2 = RTable*RTable ;
0654   G4double R1 = RTable+1.;
0655   G4double w = R1*(RTable-1.)*(RTable-1.);
0656   G4double w1 = RTable/w , w2 = -RTable*R1/w , w3 = R2/w ;
0657   G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
0658   G4bool isOut;
0659 
0660   //  loop for materials
0661   for (G4int J=0; J<numOfCouples; J++)
0662   {
0663     G4int binmax=TotBin ;
0664     G4PhysicsLinearVector* aVector =
0665                            new G4PhysicsLinearVector(0.,binmax, TotBin);
0666     Ti = lowestKineticEnergy ;
0667     G4PhysicsVector* rangeVector= (*theRangeTable)[J];
0668 
0669     for ( G4int i=0; i<TotBin; i++)
0670     {
0671       Ri = rangeVector->GetValue(Ti,isOut) ;
0672       if ( i==0 )
0673         Rim = 0. ;
0674       else
0675       {
0676         Tim = Ti/RTable ;
0677         Rim = rangeVector->GetValue(Tim,isOut);
0678       }
0679       if ( i==(TotBin-1))
0680         Rip = Ri ;
0681       else
0682       {
0683         Tip = Ti*RTable ;
0684         Rip = rangeVector->GetValue(Tip,isOut);
0685       }
0686       Value = (w1*Rip + w2*Ri + w3*Rim)/(Ti*Ti) ;
0687 
0688       aVector->PutValue(i,Value);
0689       Ti = RTable*Ti ;
0690     }
0691  
0692     theRangeCoeffATable->insert(aVector);
0693   }
0694   return theRangeCoeffATable ;
0695 }
0696 
0697 //    
0698   
0699 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRangeCoeffBTable(G4PhysicsTable* theRangeTable,
0700                              G4PhysicsTable* theRangeCoeffBTable,
0701                              G4double lowestKineticEnergy,
0702                              G4double highestKineticEnergy, G4int TotBin)
0703 // Build tables of coefficients for the energy loss calculation
0704 //  create table for coefficients "B"
0705 {
0706 
0707   G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
0708 
0709   if(theRangeCoeffBTable)
0710   { theRangeCoeffBTable->clearAndDestroy();
0711     delete theRangeCoeffBTable; }
0712   theRangeCoeffBTable = new G4PhysicsTable(numOfCouples);
0713 
0714   G4double RTable = std::exp(std::log(highestKineticEnergy/lowestKineticEnergy)/TotBin) ;
0715   G4double R2 = RTable*RTable ;
0716   G4double R1 = RTable+1.;
0717   G4double w = R1*(RTable-1.)*(RTable-1.);
0718   G4double w1 = -R1/w , w2 = R1*(R2+1.)/w , w3 = -R2*R1/w ;
0719   G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
0720   G4bool isOut;
0721 
0722   //  loop for materials
0723   for (G4int J=0; J<numOfCouples; J++)
0724   {
0725     G4int binmax=TotBin ;
0726     G4PhysicsLinearVector* aVector =
0727                         new G4PhysicsLinearVector(0.,binmax, TotBin);
0728     Ti = lowestKineticEnergy ;
0729     G4PhysicsVector* rangeVector= (*theRangeTable)[J];
0730   
0731     for ( G4int i=0; i<TotBin; i++)
0732     {
0733       Ri = rangeVector->GetValue(Ti,isOut) ;
0734       if ( i==0 )
0735          Rim = 0. ;
0736       else
0737       {
0738         Tim = Ti/RTable ;
0739         Rim = rangeVector->GetValue(Tim,isOut);
0740       }
0741       if ( i==(TotBin-1))
0742         Rip = Ri ;
0743       else
0744       {
0745         Tip = Ti*RTable ;
0746         Rip = rangeVector->GetValue(Tip,isOut);
0747       }
0748       Value = (w1*Rip + w2*Ri + w3*Rim)/Ti;
0749 
0750       aVector->PutValue(i,Value);
0751       Ti = RTable*Ti ;
0752     }
0753     theRangeCoeffBTable->insert(aVector);
0754   }
0755   return theRangeCoeffBTable ;
0756 }
0757 
0758 //    
0759 
0760 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRangeCoeffCTable(G4PhysicsTable* theRangeTable,
0761                              G4PhysicsTable* theRangeCoeffCTable,
0762                              G4double lowestKineticEnergy,
0763                              G4double highestKineticEnergy, G4int TotBin)
0764 // Build tables of coefficients for the energy loss calculation
0765 //  create table for coefficients "C"
0766 {
0767 
0768   G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
0769 
0770   if(theRangeCoeffCTable)
0771   { theRangeCoeffCTable->clearAndDestroy();
0772     delete theRangeCoeffCTable; }
0773   theRangeCoeffCTable = new G4PhysicsTable(numOfCouples);
0774 
0775   G4double RTable = std::exp(std::log(highestKineticEnergy/lowestKineticEnergy)/TotBin) ;
0776   G4double R2 = RTable*RTable ;
0777   G4double R1 = RTable+1.;
0778   G4double w = R1*(RTable-1.)*(RTable-1.);
0779   G4double w1 = 1./w , w2 = -RTable*R1/w , w3 = RTable*R2/w ;
0780   G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
0781   G4bool isOut;
0782 
0783   //  loop for materials
0784   for (G4int J=0; J<numOfCouples; J++)
0785   {
0786     G4int binmax=TotBin ;
0787     G4PhysicsLinearVector* aVector =
0788                       new G4PhysicsLinearVector(0.,binmax, TotBin);
0789     Ti = lowestKineticEnergy ;
0790     G4PhysicsVector* rangeVector= (*theRangeTable)[J];
0791   
0792     for ( G4int i=0; i<TotBin; i++)
0793     {
0794       Ri = rangeVector->GetValue(Ti,isOut) ;
0795       if ( i==0 )
0796         Rim = 0. ;
0797       else
0798       {
0799         Tim = Ti/RTable ;
0800         Rim = rangeVector->GetValue(Tim,isOut);
0801       }
0802       if ( i==(TotBin-1))
0803         Rip = Ri ;
0804       else
0805       {
0806         Tip = Ti*RTable ;
0807         Rip = rangeVector->GetValue(Tip,isOut);
0808       }
0809       Value = w1*Rip + w2*Ri + w3*Rim ;
0810 
0811       aVector->PutValue(i,Value);
0812       Ti = RTable*Ti ;
0813     }
0814     theRangeCoeffCTable->insert(aVector);
0815   }
0816   return theRangeCoeffCTable ;
0817 }
0818 
0819 //    
0820 
0821 G4double G4RDVeLowEnergyLoss::GetLossWithFluct(const G4DynamicParticle* aParticle,
0822                                              const G4MaterialCutsCouple* couple,
0823                          G4double MeanLoss,
0824                          G4double step)
0825 //  calculate actual loss from the mean loss
0826 //  The model used to get the fluctuation is essentially the same as in Glandz in Geant3.
0827 {
0828    static const G4double minLoss = 1.*eV ;
0829    static const G4double probLim = 0.01 ;
0830    static const G4double sumaLim = -std::log(probLim) ;
0831    static const G4double alim=10.;
0832    static const G4double kappa = 10. ;
0833    static const G4double factor = twopi_mc2_rcl2 ;
0834   const G4Material* aMaterial = couple->GetMaterial();
0835 
0836   // check if the material has changed ( cache mechanism)
0837 
0838   if (aMaterial != lastMaterial)
0839     {
0840       lastMaterial = aMaterial;
0841       imat         = couple->GetIndex();
0842       f1Fluct      = aMaterial->GetIonisation()->GetF1fluct();
0843       f2Fluct      = aMaterial->GetIonisation()->GetF2fluct();
0844       e1Fluct      = aMaterial->GetIonisation()->GetEnergy1fluct();
0845       e2Fluct      = aMaterial->GetIonisation()->GetEnergy2fluct();
0846       e1LogFluct   = aMaterial->GetIonisation()->GetLogEnergy1fluct();
0847       e2LogFluct   = aMaterial->GetIonisation()->GetLogEnergy2fluct();
0848       rateFluct    = aMaterial->GetIonisation()->GetRateionexcfluct();
0849       ipotFluct    = aMaterial->GetIonisation()->GetMeanExcitationEnergy();
0850       ipotLogFluct = aMaterial->GetIonisation()->GetLogMeanExcEnergy();
0851     }
0852   G4double threshold,w1,w2,C,
0853            beta2,suma,e0,loss,lossc,w;
0854   G4double a1,a2,a3;
0855   G4int p1,p2,p3;
0856   G4int nb;
0857   G4double Corrfac, na,alfa,rfac,namean,sa,alfa1,ea,sea;
0858   //  G4double dp1;
0859   G4double dp3;
0860   G4double siga ;
0861 
0862   // shortcut for very very small loss
0863   if(MeanLoss < minLoss) return MeanLoss ;
0864 
0865   // get particle data
0866   G4double Tkin   = aParticle->GetKineticEnergy();
0867 
0868   //  G4cout << "MGP -- Fluc Tkin " << Tkin/keV << " keV "  << " MeanLoss = " << MeanLoss/keV << G4endl;
0869 
0870   threshold = (*((G4ProductionCutsTable::GetProductionCutsTable())
0871                 ->GetEnergyCutsVector(1)))[imat];
0872   G4double rmass = electron_mass_c2/ParticleMass;
0873   G4double tau   = Tkin/ParticleMass, tau1 = tau+1., tau2 = tau*(tau+2.);
0874   G4double Tm    = 2.*electron_mass_c2*tau2/(1.+2.*tau1*rmass+rmass*rmass);
0875 
0876   // G4cout << "MGP Particle mass " << ParticleMass/MeV << " Tm " << Tm << G4endl;
0877 
0878   if(Tm > threshold) Tm = threshold;
0879   beta2 = tau2/(tau1*tau1);
0880 
0881   // Gaussian fluctuation ?
0882   if(MeanLoss >= kappa*Tm || MeanLoss <= kappa*ipotFluct)
0883   {
0884     G4double electronDensity = aMaterial->GetElectronDensity() ;
0885     siga = std::sqrt(Tm*(1.0-0.5*beta2)*step*
0886                 factor*electronDensity/beta2) ;
0887     do {
0888       loss = G4RandGauss::shoot(MeanLoss,siga) ;
0889     } while (loss < 0. || loss > 2.0*MeanLoss);
0890     return loss ;
0891   }
0892 
0893   w1 = Tm/ipotFluct;
0894   w2 = std::log(2.*electron_mass_c2*tau2);
0895 
0896   C = MeanLoss*(1.-rateFluct)/(w2-ipotLogFluct-beta2);
0897 
0898   a1 = C*f1Fluct*(w2-e1LogFluct-beta2)/e1Fluct;
0899   a2 = C*f2Fluct*(w2-e2LogFluct-beta2)/e2Fluct;
0900   a3 = rateFluct*MeanLoss*(Tm-ipotFluct)/(ipotFluct*Tm*std::log(w1));
0901 
0902   suma = a1+a2+a3;
0903 
0904   loss = 0. ;
0905 
0906   if(suma < sumaLim)             // very small Step
0907     {
0908       e0 = aMaterial->GetIonisation()->GetEnergy0fluct();
0909       // G4cout << "MGP e0 = " << e0/keV << G4endl;
0910 
0911       if(Tm == ipotFluct)
0912     {
0913       a3 = MeanLoss/e0;
0914 
0915       if(a3>alim)
0916         {
0917           siga=std::sqrt(a3) ;
0918           p3 = std::max(0,G4int(G4RandGauss::shoot(a3,siga)+0.5));
0919         }
0920       else p3 = G4Poisson(a3);
0921 
0922       loss = p3*e0 ;
0923 
0924       if(p3 > 0) loss += (1.-2.*G4UniformRand())*e0 ;
0925       // G4cout << "MGP very small step " << loss/keV << G4endl;
0926     }
0927       else
0928     {
0929       //      G4cout << "MGP old Tm = " << Tm << " " << ipotFluct << " " << e0 << G4endl;
0930       Tm = Tm-ipotFluct+e0 ;
0931 
0932       // MGP ---- workaround to avoid log argument<0, TO BE CHECKED
0933       if (Tm <= 0.)
0934         {
0935           loss = MeanLoss;
0936           p3 = 0;
0937           // G4cout << "MGP correction loss = MeanLoss " << loss/keV << G4endl;
0938         }
0939       else
0940         {
0941           a3 = MeanLoss*(Tm-e0)/(Tm*e0*std::log(Tm/e0));
0942 
0943           // G4cout << "MGP new Tm = " << Tm << " " << ipotFluct << " " << e0 << " a3= " << a3 << G4endl;
0944           
0945           if(a3>alim)
0946         {
0947           siga=std::sqrt(a3) ;
0948           p3 = std::max(0,G4int(G4RandGauss::shoot(a3,siga)+0.5));
0949         }
0950           else
0951         p3 = G4Poisson(a3);
0952           //G4cout << "MGP p3 " << p3 << G4endl;
0953 
0954         }
0955           
0956       if(p3 > 0)
0957         {
0958           w = (Tm-e0)/Tm ;
0959           if(p3 > nmaxCont2)
0960         {
0961           // G4cout << "MGP dp3 " << dp3 << " p3 " << p3 << " " << nmaxCont2 << G4endl;
0962           dp3 = G4double(p3) ;
0963           Corrfac = dp3/G4double(nmaxCont2) ;
0964           p3 = nmaxCont2 ;
0965         }
0966           else
0967         Corrfac = 1. ;
0968           
0969           for(G4int i=0; i<p3; i++) loss += 1./(1.-w*G4UniformRand()) ;
0970           loss *= e0*Corrfac ; 
0971           // G4cout << "MGP Corrfac = " << Corrfac << " e0 = " << e0/keV << " loss = " << loss/keV << G4endl;
0972         }        
0973     }
0974     }
0975 
0976   else                              // not so small Step
0977     {
0978       // excitation type 1
0979       if(a1>alim)
0980       {
0981         siga=std::sqrt(a1) ;
0982         p1 = std::max(0,int(G4RandGauss::shoot(a1,siga)+0.5));
0983       }
0984       else
0985        p1 = G4Poisson(a1);
0986 
0987       // excitation type 2
0988       if(a2>alim)
0989       {
0990         siga=std::sqrt(a2) ;
0991         p2 = std::max(0,int(G4RandGauss::shoot(a2,siga)+0.5));
0992       }
0993       else
0994         p2 = G4Poisson(a2);
0995 
0996       loss = p1*e1Fluct+p2*e2Fluct;
0997  
0998       // smearing to avoid unphysical peaks
0999       if(p2 > 0)
1000         loss += (1.-2.*G4UniformRand())*e2Fluct;
1001       else if (loss>0.)
1002         loss += (1.-2.*G4UniformRand())*e1Fluct;   
1003       
1004       // ionisation .......................................
1005       if(a3 > 0.)
1006     {
1007       if(a3>alim)
1008         {
1009           siga=std::sqrt(a3) ;
1010           p3 = std::max(0,int(G4RandGauss::shoot(a3,siga)+0.5));
1011         }
1012       else
1013         p3 = G4Poisson(a3);
1014       
1015       lossc = 0.;
1016       if(p3 > 0)
1017         {
1018           na = 0.; 
1019           alfa = 1.;
1020           if (p3 > nmaxCont2)
1021         {
1022           dp3        = G4double(p3);
1023           rfac       = dp3/(G4double(nmaxCont2)+dp3);
1024           namean     = G4double(p3)*rfac;
1025           sa         = G4double(nmaxCont1)*rfac;
1026           na         = G4RandGauss::shoot(namean,sa);
1027           if (na > 0.)
1028             {
1029               alfa   = w1*G4double(nmaxCont2+p3)/(w1*G4double(nmaxCont2)+G4double(p3));
1030               alfa1  = alfa*std::log(alfa)/(alfa-1.);
1031               ea     = na*ipotFluct*alfa1;
1032               sea    = ipotFluct*std::sqrt(na*(alfa-alfa1*alfa1));
1033               lossc += G4RandGauss::shoot(ea,sea);
1034             }
1035         }
1036           
1037           nb = G4int(G4double(p3)-na);
1038           if (nb > 0)
1039         {
1040           w2 = alfa*ipotFluct;
1041           w  = (Tm-w2)/Tm;      
1042           for (G4int k=0; k<nb; k++) lossc += w2/(1.-w*G4UniformRand());
1043         }
1044         }        
1045       
1046       loss += lossc;  
1047     }
1048     } 
1049   
1050   return loss ;
1051 }
1052 
1053 //    
1054
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!