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
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056 #include "G4PhysicsVector.hh"
00057 #include <iomanip>
00058
00059 G4Allocator<G4PhysicsVector> aPVAllocator;
00060
00061
00062
00063 G4PhysicsVector::G4PhysicsVector(G4bool spline)
00064 : type(T_G4PhysicsVector),
00065 edgeMin(0.), edgeMax(0.), numberOfNodes(0),
00066 useSpline(spline),
00067 dBin(0.), baseBin(0.),
00068 verboseLevel(0)
00069 {
00070 cache = new G4PhysicsVectorCache();
00071 }
00072
00073
00074
00075 G4PhysicsVector::~G4PhysicsVector()
00076 {
00077 delete cache; cache =0;
00078 }
00079
00080
00081
00082 G4PhysicsVector::G4PhysicsVector(const G4PhysicsVector& right)
00083 {
00084 cache = new G4PhysicsVectorCache();
00085 dBin = right.dBin;
00086 baseBin = right.baseBin;
00087 verboseLevel = right.verboseLevel;
00088
00089 DeleteData();
00090 CopyData(right);
00091 }
00092
00093
00094
00095 G4PhysicsVector& G4PhysicsVector::operator=(const G4PhysicsVector& right)
00096 {
00097 if (&right==this) { return *this; }
00098 dBin = right.dBin;
00099 baseBin = right.baseBin;
00100 verboseLevel = right.verboseLevel;
00101
00102 DeleteData();
00103 CopyData(right);
00104 return *this;
00105 }
00106
00107
00108
00109 G4int G4PhysicsVector::operator==(const G4PhysicsVector &right) const
00110 {
00111 return (this == &right);
00112 }
00113
00114
00115
00116 G4int G4PhysicsVector::operator!=(const G4PhysicsVector &right) const
00117 {
00118 return (this != &right);
00119 }
00120
00121
00122
00123 void G4PhysicsVector::DeleteData()
00124 {
00125 secDerivative.clear();
00126 }
00127
00128
00129
00130 void G4PhysicsVector::CopyData(const G4PhysicsVector& vec)
00131 {
00132 type = vec.type;
00133 edgeMin = vec.edgeMin;
00134 edgeMax = vec.edgeMax;
00135 numberOfNodes = vec.numberOfNodes;
00136 cache->lastEnergy = vec.GetLastEnergy();
00137 cache->lastValue = vec.GetLastValue();
00138 cache->lastBin = vec.GetLastBin();
00139 useSpline = vec.useSpline;
00140
00141 size_t i;
00142 dataVector.clear();
00143 for(i=0; i<(vec.dataVector).size(); i++){
00144 dataVector.push_back( (vec.dataVector)[i] );
00145 }
00146 binVector.clear();
00147 for(i=0; i<(vec.binVector).size(); i++){
00148 binVector.push_back( (vec.binVector)[i] );
00149 }
00150 secDerivative.clear();
00151 for(i=0; i<(vec.secDerivative).size(); i++){
00152 secDerivative.push_back( (vec.secDerivative)[i] );
00153 }
00154 }
00155
00156
00157
00158 G4double G4PhysicsVector::GetLowEdgeEnergy(size_t binNumber) const
00159 {
00160 return binVector[binNumber];
00161 }
00162
00163
00164
00165 G4bool G4PhysicsVector::Store(std::ofstream& fOut, G4bool ascii)
00166 {
00167
00168 if (ascii)
00169 {
00170 fOut << *this;
00171 return true;
00172 }
00173
00174
00175
00176 fOut.write((char*)(&edgeMin), sizeof edgeMin);
00177 fOut.write((char*)(&edgeMax), sizeof edgeMax);
00178 fOut.write((char*)(&numberOfNodes), sizeof numberOfNodes);
00179
00180
00181 size_t size = dataVector.size();
00182 fOut.write((char*)(&size), sizeof size);
00183
00184 G4double* value = new G4double[2*size];
00185 for(size_t i = 0; i < size; ++i)
00186 {
00187 value[2*i] = binVector[i];
00188 value[2*i+1]= dataVector[i];
00189 }
00190 fOut.write((char*)(value), 2*size*(sizeof (G4double)));
00191 delete [] value;
00192
00193 return true;
00194 }
00195
00196
00197
00198 G4bool G4PhysicsVector::Retrieve(std::ifstream& fIn, G4bool ascii)
00199 {
00200
00201 cache->lastEnergy=-DBL_MAX;
00202 cache->lastValue =0.;
00203 cache->lastBin =0;
00204 dataVector.clear();
00205 binVector.clear();
00206 secDerivative.clear();
00207
00208
00209 if (ascii){
00210
00211 fIn >> edgeMin >> edgeMax >> numberOfNodes;
00212 if (fIn.fail()) { return false; }
00213
00214 G4int siz=0;
00215 fIn >> siz;
00216 if (fIn.fail()) { return false; }
00217 if (siz<=0)
00218 {
00219 #ifdef G4VERBOSE
00220 G4cerr << "G4PhysicsVector::Retrieve():";
00221 G4cerr << " Invalid vector size: " << siz << G4endl;
00222 #endif
00223 return false;
00224 }
00225
00226 binVector.reserve(siz);
00227 dataVector.reserve(siz);
00228 G4double vBin, vData;
00229
00230 for(G4int i = 0; i < siz ; i++)
00231 {
00232 vBin = 0.;
00233 vData= 0.;
00234 fIn >> vBin >> vData;
00235 if (fIn.fail()) { return false; }
00236 binVector.push_back(vBin);
00237 dataVector.push_back(vData);
00238 }
00239
00240
00241 numberOfNodes = siz;
00242 edgeMin = binVector[0];
00243 edgeMax = binVector[numberOfNodes-1];
00244 return true ;
00245 }
00246
00247
00248
00249 fIn.read((char*)(&edgeMin), sizeof edgeMin);
00250 fIn.read((char*)(&edgeMax), sizeof edgeMax);
00251 fIn.read((char*)(&numberOfNodes), sizeof numberOfNodes );
00252
00253
00254 size_t size;
00255 fIn.read((char*)(&size), sizeof size);
00256
00257 G4double* value = new G4double[2*size];
00258 fIn.read((char*)(value), 2*size*(sizeof(G4double)) );
00259 if (G4int(fIn.gcount()) != G4int(2*size*(sizeof(G4double))) )
00260 {
00261 delete [] value;
00262 return false;
00263 }
00264
00265 binVector.reserve(size);
00266 dataVector.reserve(size);
00267 for(size_t i = 0; i < size; ++i)
00268 {
00269 binVector.push_back(value[2*i]);
00270 dataVector.push_back(value[2*i+1]);
00271 }
00272 delete [] value;
00273
00274
00275 numberOfNodes = size;
00276 edgeMin = binVector[0];
00277 edgeMax = binVector[numberOfNodes-1];
00278
00279 return true;
00280 }
00281
00282
00283
00284 void
00285 G4PhysicsVector::ScaleVector(G4double factorE, G4double factorV)
00286 {
00287 size_t n = dataVector.size();
00288 size_t i;
00289 if(n > 0) {
00290 for(i=0; i<n; ++i) {
00291 binVector[i] *= factorE;
00292 dataVector[i] *= factorV;
00293 }
00294 }
00295
00296
00297 secDerivative.clear();
00298
00299 edgeMin *= factorE;
00300 edgeMax *= factorE;
00301 cache->lastEnergy = factorE*(cache->lastEnergy);
00302 cache->lastValue = factorV*(cache->lastValue);
00303 }
00304
00305
00306
00307 void
00308 G4PhysicsVector::ComputeSecondDerivatives(G4double firstPointDerivative,
00309 G4double endPointDerivative)
00310
00311
00312
00313
00314 {
00315 if(4 > numberOfNodes)
00316 {
00317 ComputeSecDerivatives();
00318 return;
00319 }
00320
00321 if(!SplinePossible()) { return; }
00322
00323 G4int n = numberOfNodes-1;
00324
00325 G4double* u = new G4double [n];
00326
00327 G4double p, sig, un;
00328
00329 u[0] = (6.0/(binVector[1]-binVector[0]))
00330 * ((dataVector[1]-dataVector[0])/(binVector[1]-binVector[0])
00331 - firstPointDerivative);
00332
00333 secDerivative[0] = - 0.5;
00334
00335
00336
00337
00338 for(G4int i=1; i<n; ++i)
00339 {
00340 sig = (binVector[i]-binVector[i-1]) / (binVector[i+1]-binVector[i-1]);
00341 p = sig*(secDerivative[i-1]) + 2.0;
00342 secDerivative[i] = (sig - 1.0)/p;
00343 u[i] = (dataVector[i+1]-dataVector[i])/(binVector[i+1]-binVector[i])
00344 - (dataVector[i]-dataVector[i-1])/(binVector[i]-binVector[i-1]);
00345 u[i] = 6.0*u[i]/(binVector[i+1]-binVector[i-1]) - sig*u[i-1]/p;
00346 }
00347
00348 sig = (binVector[n-1]-binVector[n-2]) / (binVector[n]-binVector[n-2]);
00349 p = sig*secDerivative[n-2] + 2.0;
00350 un = (6.0/(binVector[n]-binVector[n-1]))
00351 *(endPointDerivative -
00352 (dataVector[n]-dataVector[n-1])/(binVector[n]-binVector[n-1])) - u[n-1]/p;
00353 secDerivative[n] = un/(secDerivative[n-1] + 2.0);
00354
00355
00356
00357
00358 for(G4int k=n-1; k>0; --k)
00359 {
00360 secDerivative[k] *=
00361 (secDerivative[k+1] -
00362 u[k]*(binVector[k+1]-binVector[k-1])/(binVector[k+1]-binVector[k]));
00363 }
00364 secDerivative[0] = 0.5*(u[0] - secDerivative[1]);
00365
00366 delete [] u;
00367 }
00368
00369
00370
00371 void G4PhysicsVector::FillSecondDerivatives()
00372
00373
00374
00375 {
00376 if(5 > numberOfNodes)
00377 {
00378 ComputeSecDerivatives();
00379 return;
00380 }
00381
00382 if(!SplinePossible()) { return; }
00383
00384 G4int n = numberOfNodes-1;
00385
00386 G4double* u = new G4double [n];
00387
00388 G4double p, sig;
00389
00390 u[1] = ((dataVector[2]-dataVector[1])/(binVector[2]-binVector[1]) -
00391 (dataVector[1]-dataVector[0])/(binVector[1]-binVector[0]));
00392 u[1] = 6.0*u[1]*(binVector[2]-binVector[1])
00393 / ((binVector[2]-binVector[0])*(binVector[2]-binVector[0]));
00394
00395
00396
00397
00398 secDerivative[1] = (2.0*binVector[1]-binVector[0]-binVector[2])
00399 / (2.0*binVector[2]-binVector[0]-binVector[1]);
00400
00401 for(G4int i=2; i<n-1; ++i)
00402 {
00403 sig = (binVector[i]-binVector[i-1]) / (binVector[i+1]-binVector[i-1]);
00404 p = sig*secDerivative[i-1] + 2.0;
00405 secDerivative[i] = (sig - 1.0)/p;
00406 u[i] = (dataVector[i+1]-dataVector[i])/(binVector[i+1]-binVector[i])
00407 - (dataVector[i]-dataVector[i-1])/(binVector[i]-binVector[i-1]);
00408 u[i] = (6.0*u[i]/(binVector[i+1]-binVector[i-1])) - sig*u[i-1]/p;
00409 }
00410
00411 sig = (binVector[n-1]-binVector[n-2]) / (binVector[n]-binVector[n-2]);
00412 p = sig*secDerivative[n-3] + 2.0;
00413 u[n-1] = (dataVector[n]-dataVector[n-1])/(binVector[n]-binVector[n-1])
00414 - (dataVector[n-1]-dataVector[n-2])/(binVector[n-1]-binVector[n-2]);
00415 u[n-1] = 6.0*sig*u[n-1]/(binVector[n]-binVector[n-2])
00416 - (2.0*sig - 1.0)*u[n-2]/p;
00417
00418 p = (1.0+sig) + (2.0*sig-1.0)*secDerivative[n-2];
00419 secDerivative[n-1] = u[n-1]/p;
00420
00421
00422
00423
00424 for(G4int k=n-2; k>1; --k)
00425 {
00426 secDerivative[k] *=
00427 (secDerivative[k+1] -
00428 u[k]*(binVector[k+1]-binVector[k-1])/(binVector[k+1]-binVector[k]));
00429 }
00430 secDerivative[n] = (secDerivative[n-1] - (1.0-sig)*secDerivative[n-2])/sig;
00431 sig = 1.0 - ((binVector[2]-binVector[1])/(binVector[2]-binVector[0]));
00432 secDerivative[1] *= (secDerivative[2] - u[1]/(1.0-sig));
00433 secDerivative[0] = (secDerivative[1] - sig*secDerivative[2])/(1.0-sig);
00434
00435 delete [] u;
00436 }
00437
00438
00439
00440 void
00441 G4PhysicsVector::ComputeSecDerivatives()
00442
00443 {
00444 if(!SplinePossible()) { return; }
00445
00446 if(3 > numberOfNodes)
00447 {
00448 useSpline = false;
00449 return;
00450 }
00451
00452 size_t n = numberOfNodes-1;
00453
00454 for(size_t i=1; i<n; ++i)
00455 {
00456 secDerivative[i] =
00457 3.0*((dataVector[i+1]-dataVector[i])/(binVector[i+1]-binVector[i]) -
00458 (dataVector[i]-dataVector[i-1])/(binVector[i]-binVector[i-1]))
00459 /(binVector[i+1]-binVector[i-1]);
00460 }
00461 secDerivative[n] = secDerivative[n-1];
00462 secDerivative[0] = secDerivative[1];
00463 }
00464
00465
00466
00467 G4bool G4PhysicsVector::SplinePossible()
00468
00469
00470 {
00471 secDerivative.clear();
00472 if(!useSpline) { return useSpline; }
00473 secDerivative.reserve(numberOfNodes);
00474 for(size_t j=0; j<numberOfNodes; ++j)
00475 {
00476 secDerivative.push_back(0.0);
00477 if(j > 0)
00478 {
00479 if(binVector[j]-binVector[j-1] <= 0.) { useSpline = false; }
00480 }
00481 }
00482 return useSpline;
00483 }
00484
00485
00486
00487 std::ostream& operator<<(std::ostream& out, const G4PhysicsVector& pv)
00488 {
00489
00490 out << std::setprecision(12) << pv.edgeMin << " "
00491 << pv.edgeMax << " " << pv.numberOfNodes << G4endl;
00492
00493
00494 out << pv.dataVector.size() << G4endl;
00495 for(size_t i = 0; i < pv.dataVector.size(); i++)
00496 {
00497 out << pv.binVector[i] << " " << pv.dataVector[i] << G4endl;
00498 }
00499 out << std::setprecision(6);
00500
00501 return out;
00502 }
00503
00504
00505
00506 void G4PhysicsVector::ComputeValue(G4double theEnergy)
00507 {
00508
00509
00510
00511
00512 if( theEnergy < cache->lastEnergy
00513 && theEnergy >= binVector[cache->lastBin]) {
00514 cache->lastEnergy = theEnergy;
00515 Interpolation(cache->lastBin);
00516
00517 } else if( theEnergy <= edgeMin ) {
00518 cache->lastBin = 0;
00519 cache->lastEnergy = edgeMin;
00520 cache->lastValue = dataVector[0];
00521
00522 } else if( theEnergy >= edgeMax ) {
00523 cache->lastBin = numberOfNodes-1;
00524 cache->lastEnergy = edgeMax;
00525 cache->lastValue = dataVector[cache->lastBin];
00526
00527 } else {
00528 cache->lastBin = FindBinLocation(theEnergy);
00529 cache->lastEnergy = theEnergy;
00530 Interpolation(cache->lastBin);
00531 }
00532 }
00533