00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035 #include "G4VRangeToEnergyConverter.hh"
00036 #include "G4ParticleTable.hh"
00037 #include "G4Material.hh"
00038 #include "G4MaterialTable.hh"
00039 #include "G4PhysicsLogVector.hh"
00040
00041 #include "G4ios.hh"
00042 #include "G4SystemOfUnits.hh"
00043
00044
00045 G4double G4VRangeToEnergyConverter::LowestEnergy = 0.99e-3*MeV;
00046 G4double G4VRangeToEnergyConverter::HighestEnergy = 100.0e6*MeV;
00047
00048
00049 G4double G4VRangeToEnergyConverter::MaxEnergyCut = 10.0*GeV;
00050
00051 G4VRangeToEnergyConverter::G4VRangeToEnergyConverter():
00052 theParticle(0), theLossTable(0), NumberOfElements(0), TotBin(300),
00053 verboseLevel(1)
00054 {
00055 fMaxEnergyCut = 0.;
00056 }
00057
00058 G4VRangeToEnergyConverter::G4VRangeToEnergyConverter(const G4VRangeToEnergyConverter& right) : theParticle(right.theParticle), theLossTable(0), TotBin(right.TotBin)
00059 {
00060 fMaxEnergyCut = 0.;
00061 *this = right;
00062 }
00063
00064 G4VRangeToEnergyConverter & G4VRangeToEnergyConverter::operator=(const G4VRangeToEnergyConverter &right)
00065 {
00066 if (this == &right) return *this;
00067 if (theLossTable) {
00068 theLossTable->clearAndDestroy();
00069 delete theLossTable;
00070 theLossTable=0;
00071 }
00072
00073 NumberOfElements = right.NumberOfElements;
00074 theParticle = right.theParticle;
00075 verboseLevel = right.verboseLevel;
00076
00077
00078 theLossTable = new G4LossTable();
00079 theLossTable->reserve(G4Element::GetNumberOfElements());
00080
00081 for (size_t j=0; j<size_t(NumberOfElements); j++){
00082 G4LossVector* aVector= new
00083 G4LossVector(LowestEnergy, MaxEnergyCut, TotBin);
00084 for (size_t i=0; i<size_t(TotBin); i++) {
00085 G4double Value = (*((*right.theLossTable)[j]))[i];
00086 aVector->PutValue(i,Value);
00087 }
00088 theLossTable->insert(aVector);
00089 }
00090
00091
00092 for (size_t idx=0; idx<fRangeVectorStore.size(); idx++){
00093 delete fRangeVectorStore.at(idx);
00094 }
00095 fRangeVectorStore.clear();
00096
00097
00098 for (size_t j=0; j<((right.fRangeVectorStore).size()); j++){
00099 G4RangeVector* vector = (right.fRangeVectorStore).at(j);
00100 G4RangeVector* rangeVector = 0;
00101 if (vector !=0 ) {
00102 rangeVector = new G4RangeVector(LowestEnergy, MaxEnergyCut, TotBin);
00103 fMaxEnergyCut = MaxEnergyCut;
00104 for (size_t i=0; i<size_t(TotBin); i++) {
00105 G4double Value = (*vector)[i];
00106 rangeVector->PutValue(i,Value);
00107 }
00108 }
00109 fRangeVectorStore.push_back(rangeVector);
00110 }
00111 return *this;
00112 }
00113
00114
00115 G4VRangeToEnergyConverter::~G4VRangeToEnergyConverter()
00116 {
00117 Reset();
00118 }
00119
00120 G4int G4VRangeToEnergyConverter::operator==(const G4VRangeToEnergyConverter &right) const
00121 {
00122 return this == &right;
00123 }
00124
00125 G4int G4VRangeToEnergyConverter::operator!=(const G4VRangeToEnergyConverter &right) const
00126 {
00127 return this != &right;
00128 }
00129
00130
00131
00132
00133
00134 G4double G4VRangeToEnergyConverter::Convert(G4double rangeCut,
00135 const G4Material* material)
00136 {
00137 #ifdef G4VERBOSE
00138 if (GetVerboseLevel()>3) {
00139 G4cout << "G4VRangeToEnergyConverter::Convert() ";
00140 G4cout << "Convert for " << material->GetName()
00141 << " with Range Cut " << rangeCut/mm << "[mm]" << G4endl;
00142 }
00143 #endif
00144
00145 G4double theKineticEnergyCuts = 0.;
00146
00147 if (fMaxEnergyCut != MaxEnergyCut) {
00148 fMaxEnergyCut = MaxEnergyCut;
00149
00150 Reset();
00151 }
00152
00153
00154 BuildLossTable();
00155
00156
00157
00158 G4double tune = 0.025*mm*g/cm3 ,lowen = 30.*keV ;
00159
00160
00161 G4double density = material->GetDensity() ;
00162 if(density <= 0.) {
00163 #ifdef G4VERBOSE
00164 if (GetVerboseLevel()>0) {
00165 G4cout << "G4VRangeToEnergyConverter::Convert() ";
00166 G4cout << material->GetName() << "has zero density "
00167 << "( " << density << ")" << G4endl;
00168 }
00169 #endif
00170 return 0.;
00171 }
00172
00173
00174 const G4MaterialTable* table = G4Material::GetMaterialTable();
00175 G4int ext_size = table->size() - fRangeVectorStore.size();
00176 for (int i=0; i<ext_size; i++) fRangeVectorStore.push_back(0);
00177
00178
00179 G4int idx = material->GetIndex();
00180 G4RangeVector* rangeVector = fRangeVectorStore.at(idx);
00181 if (rangeVector == 0) {
00182 rangeVector = new G4RangeVector(LowestEnergy, MaxEnergyCut, TotBin);
00183 BuildRangeVector(material, rangeVector);
00184 fRangeVectorStore.at(idx) = rangeVector;
00185 }
00186
00187
00188 theKineticEnergyCuts = ConvertCutToKineticEnergy(rangeVector, rangeCut, idx);
00189
00190 if( ((theParticle->GetParticleName()=="e-")||(theParticle->GetParticleName()=="e+"))
00191 && (theKineticEnergyCuts < lowen) ) {
00192
00193 theKineticEnergyCuts /= (1.+(1.-theKineticEnergyCuts/lowen)*
00194 tune/(rangeCut*density));
00195 }
00196
00197 if(theKineticEnergyCuts < LowestEnergy) {
00198 theKineticEnergyCuts = LowestEnergy ;
00199 } else if(theKineticEnergyCuts > MaxEnergyCut) {
00200 theKineticEnergyCuts = MaxEnergyCut;
00201 }
00202
00203 return theKineticEnergyCuts;
00204 }
00205
00206
00207
00208
00209 void G4VRangeToEnergyConverter::SetEnergyRange(G4double lowedge,
00210 G4double highedge)
00211 {
00212
00213 if ( (lowedge<0.0)||(highedge<=lowedge) ){
00214 #ifdef G4VERBOSE
00215 G4cerr << "Error in G4VRangeToEnergyConverter::SetEnergyRange";
00216 G4cerr << " : illegal energy range" << "(" << lowedge/GeV;
00217 G4cerr << "," << highedge/GeV << ") [GeV]" << G4endl;
00218 #endif
00219 G4Exception( "G4VRangeToEnergyConverter::SetEnergyRange()",
00220 "ProcCuts101",
00221 JustWarning, "Illegal energy range ");
00222 } else {
00223 LowestEnergy = lowedge;
00224 HighestEnergy = highedge;
00225 }
00226 }
00227
00228
00229 G4double G4VRangeToEnergyConverter::GetLowEdgeEnergy()
00230 {
00231 return LowestEnergy;
00232 }
00233
00234
00235 G4double G4VRangeToEnergyConverter::GetHighEdgeEnergy()
00236 {
00237 return HighestEnergy;
00238 }
00239
00240
00241
00242
00243 G4double G4VRangeToEnergyConverter::GetMaxEnergyCut()
00244 {
00245 return MaxEnergyCut;
00246 }
00247
00248 void G4VRangeToEnergyConverter::SetMaxEnergyCut(G4double value)
00249 {
00250 MaxEnergyCut = value;
00251 }
00252
00253
00254
00255
00256 void G4VRangeToEnergyConverter::Reset()
00257 {
00258
00259 if (theLossTable) {
00260 theLossTable->clearAndDestroy();
00261 delete theLossTable;
00262 }
00263 theLossTable=0;
00264 NumberOfElements=0;
00265
00266
00267 for (size_t idx=0; idx<fRangeVectorStore.size(); idx++){
00268 delete fRangeVectorStore.at(idx);
00269 }
00270 fRangeVectorStore.clear();
00271 }
00272
00273
00274
00275
00276
00277
00278
00279 void G4VRangeToEnergyConverter::BuildLossTable()
00280 {
00281 if (size_t(NumberOfElements) == G4Element::GetNumberOfElements()) return;
00282
00283
00284 Reset();
00285
00286
00287 NumberOfElements = G4Element::GetNumberOfElements();
00288 theLossTable = new G4LossTable();
00289 theLossTable->reserve(G4Element::GetNumberOfElements());
00290 #ifdef G4VERBOSE
00291 if (GetVerboseLevel()>3) {
00292 G4cout << "G4VRangeToEnergyConverter::BuildLossTable() ";
00293 G4cout << "Create theLossTable[" << theLossTable << "]";
00294 G4cout << " NumberOfElements=" << NumberOfElements <<G4endl;
00295 }
00296 #endif
00297
00298
00299
00300 for (size_t j=0; j<size_t(NumberOfElements); j++){
00301 G4double Value;
00302 G4LossVector* aVector= 0;
00303 aVector= new G4LossVector(LowestEnergy, MaxEnergyCut, TotBin);
00304 for (size_t i=0; i<size_t(TotBin); i++) {
00305 Value = ComputeLoss( (*G4Element::GetElementTable())[j]->GetZ(),
00306 aVector->GetLowEdgeEnergy(i)
00307 );
00308 aVector->PutValue(i,Value);
00309 }
00310 theLossTable->insert(aVector);
00311 }
00312 }
00313
00314
00315
00316
00317 void G4VRangeToEnergyConverter::BuildRangeVector(const G4Material* aMaterial,
00318 G4PhysicsLogVector* rangeVector)
00319 {
00320
00321 const G4ElementVector* elementVector = aMaterial->GetElementVector();
00322 const G4double* atomicNumDensityVector = aMaterial->GetAtomicNumDensityVector();
00323 G4int NumEl = aMaterial->GetNumberOfElements();
00324
00325
00326 size_t i;
00327 std::vector<G4double> lossV;
00328 for ( size_t ib=0; ib<size_t(TotBin); ib++) {
00329 G4double loss=0.;
00330 for (i=0; i<size_t(NumEl); i++) {
00331 G4int IndEl = (*elementVector)[i]->GetIndex();
00332 loss += atomicNumDensityVector[i]*
00333 (*((*theLossTable)[IndEl]))[ib];
00334 }
00335 lossV.push_back(loss);
00336 }
00337
00338
00339 G4double ltt = std::log(MaxEnergyCut/LowestEnergy);
00340 G4double dltau = ltt/TotBin;
00341
00342 G4double s0 = 0.;
00343 G4double Value;
00344 for ( i=0; i<size_t(TotBin); i++) {
00345 G4double t = rangeVector->GetLowEdgeEnergy(i);
00346 G4double q = t/lossV[i];
00347 if (i==0) s0 += 0.5*q;
00348 else s0 += q;
00349
00350 if (i==0) {
00351 Value = (s0 + 0.5*q)*dltau ;
00352 } else {
00353 Value = (s0 - 0.5*q)*dltau ;
00354 }
00355 rangeVector->PutValue(i,Value);
00356 }
00357 }
00358
00359
00360
00361
00362 G4double G4VRangeToEnergyConverter::ConvertCutToKineticEnergy(
00363 G4RangeVector* rangeVector,
00364 G4double theCutInLength,
00365 #ifdef G4VERBOSE
00366 size_t materialIndex
00367 #else
00368 size_t
00369 #endif
00370 ) const
00371 {
00372 const G4double epsilon=0.01;
00373
00374
00375 G4double rmax= -1.e10*mm;
00376
00377 G4double T1 = LowestEnergy;
00378 G4double r1 =(*rangeVector)[0] ;
00379
00380 G4double T2 = MaxEnergyCut;
00381
00382
00383 if ( theCutInLength <= r1 ) { return T1; }
00384
00385
00386
00387 for (size_t ibin=0; ibin<size_t(TotBin); ibin++) {
00388 G4double T=rangeVector->GetLowEdgeEnergy(ibin);
00389 G4double r=(*rangeVector)[ibin];
00390 if ( r>rmax ) rmax=r;
00391 if (r <theCutInLength ) {
00392 T1 = T;
00393 r1 = r;
00394 } else if (r >theCutInLength ) {
00395 T2 = T;
00396 break;
00397 }
00398 }
00399
00400
00401 if ( theCutInLength >= rmax ) {
00402 #ifdef G4VERBOSE
00403 if (GetVerboseLevel()>2) {
00404 G4cout << "G4VRangeToEnergyConverter::ConvertCutToKineticEnergy ";
00405 G4cout << " for " << theParticle->GetParticleName() << G4endl;
00406 G4cout << "The cut in range [" << theCutInLength/mm << " (mm)] ";
00407 G4cout << " is too big " ;
00408 G4cout << " for material idx=" << materialIndex <<G4endl;
00409 }
00410 #endif
00411 return MaxEnergyCut;
00412 }
00413
00414
00415 G4double T3 = std::sqrt(T1*T2);
00416 G4double r3 = rangeVector->Value(T3);
00417 while ( std::fabs(1.-r3/theCutInLength)>epsilon ) {
00418 if ( theCutInLength <= r3 ) {
00419 T2 = T3;
00420 } else {
00421 T1 = T3;
00422 }
00423 T3 = std::sqrt(T1*T2);
00424 r3 = rangeVector->Value(T3);
00425 }
00426
00427 return T3;
00428 }
00429