Changeset 1193 for trunk/source/global/management
- Timestamp:
- Nov 19, 2009, 3:40:00 PM (15 years ago)
- Location:
- trunk/source/global/management
- Files:
-
- 20 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/global/management/include/G4Allocator.hh
r1058 r1193 25 25 // 26 26 // 27 // $Id: G4Allocator.hh,v 1.1 8 2006/06/29 19:01:16 gunterExp $28 // GEANT4 tag $Name: g eant4-09-02-ref-02$27 // $Id: G4Allocator.hh,v 1.19 2009/10/29 16:01:28 gcosmo Exp $ 28 // GEANT4 tag $Name: global-V09-02-10 $ 29 29 // 30 30 // … … 99 99 // Returns the address of values 100 100 101 pointer allocate(size_type n, void* hint= 0)101 pointer allocate(size_type n, void* = 0) 102 102 { 103 103 // Allocates space for n elements of type Type, but does not initialise -
trunk/source/global/management/include/G4ErrorPropagatorData.hh
r1058 r1193 25 25 // 26 26 // 27 // $Id: G4ErrorPropagatorData.hh,v 1. 3 2007/06/08 10:33:47 gcosmo Exp $28 // GEANT4 tag $Name: g eant4-09-02-ref-02$27 // $Id: G4ErrorPropagatorData.hh,v 1.5 2009/05/19 13:31:47 gcosmo Exp $ 28 // GEANT4 tag $Name: global-V09-02-10 $ 29 29 // 30 30 // … … 38 38 // and manager verbosity for the error propagation classes. 39 39 40 // History:41 40 // - Created. P.Arce, 2004. 42 41 // -------------------------------------------------------------------- … … 56 55 G4ErrorState_TargetCloserThanBoundary, 57 56 G4ErrorState_StoppedAtTarget }; 57 58 enum G4ErrorStage { G4ErrorStage_Inflation = 1, 59 G4ErrorStage_Deflation }; 60 58 61 class G4ErrorTarget; 59 62 … … 72 75 G4ErrorState GetState() const; 73 76 void SetState( G4ErrorState sta ); 77 78 G4ErrorStage GetStage() const; 79 void SetStage( G4ErrorStage sta ); 74 80 75 81 const G4ErrorTarget* GetTarget( G4bool mustExist = 0) const; … … 94 100 G4ErrorState theState; 95 101 102 G4ErrorStage theStage; 103 96 104 G4ErrorTarget* theTarget; 97 105 -
trunk/source/global/management/include/G4ErrorPropagatorData.icc
r1058 r1193 25 25 // 26 26 // 27 // $Id: G4ErrorPropagatorData.icc,v 1. 4 2007/06/08 10:33:47 gcosmoExp $28 // GEANT4 tag $Name: g eant4-09-02-ref-02$27 // $Id: G4ErrorPropagatorData.icc,v 1.5 2009/05/14 13:46:40 arce Exp $ 28 // GEANT4 tag $Name: global-V09-02-10 $ 29 29 // 30 30 // … … 59 59 60 60 inline 61 void G4ErrorPropagatorData::SetStage( G4ErrorStage sta ) 62 { 63 theStage = sta; 64 } 65 66 inline 67 G4ErrorStage G4ErrorPropagatorData::GetStage() const 68 { 69 return theStage; 70 } 71 72 inline 61 73 const G4ErrorTarget* G4ErrorPropagatorData::GetTarget( G4bool mustExist ) const 62 74 { -
trunk/source/global/management/include/G4LPhysicsFreeVector.icc
r1058 r1193 25 25 // 26 26 // 27 // $Id: G4LPhysicsFreeVector.icc,v 1.1 3 2008/09/22 08:26:33 gcosmoExp $28 // GEANT4 tag $Name: g eant4-09-02-ref-02$27 // $Id: G4LPhysicsFreeVector.icc,v 1.14 2009/06/25 10:05:26 vnivanch Exp $ 28 // GEANT4 tag $Name: global-V09-02-10 $ 29 29 // 30 30 // … … 52 52 if(binNumber == 0) 53 53 { edgeMin = binValue; } 54 else if( numberOf Bin- 1 == binNumber)54 else if( numberOfNodes - 1 == binNumber) 55 55 { edgeMax = binValue; } 56 56 } … … 84 84 { 85 85 G4int n1 = 0; 86 G4int n2 = numberOf Bin/2;87 G4int n3 = numberOf Bin- 1;86 G4int n2 = numberOfNodes/2; 87 G4int n3 = numberOfNodes - 1; 88 88 while (n1 != n3 - 1) 89 89 { -
trunk/source/global/management/include/G4PhysicsFreeVector.hh
r1058 r1193 25 25 // 26 26 // 27 // $Id: G4PhysicsFreeVector.hh,v 1.1 2 2008/09/22 14:49:57 gcosmoExp $28 // GEANT4 tag $Name: g eant4-09-02-ref-02$27 // $Id: G4PhysicsFreeVector.hh,v 1.13 2009/06/25 10:05:26 vnivanch Exp $ 28 // GEANT4 tag $Name: global-V09-02-10 $ 29 29 // 30 30 // … … 108 108 109 109 size_t lowerBound = 0; 110 size_t upperBound = numberOf Bin-1;110 size_t upperBound = numberOfNodes-2; 111 111 112 112 while (lowerBound <= upperBound) -
trunk/source/global/management/include/G4PhysicsOrderedFreeVector.hh
r1058 r1193 25 25 // 26 26 // 27 // $Id: G4PhysicsOrderedFreeVector.hh,v 1.1 1 2008/09/22 14:49:57 gcosmoExp $28 // GEANT4 tag $Name: g eant4-09-02-ref-02$27 // $Id: G4PhysicsOrderedFreeVector.hh,v 1.12 2009/06/25 10:05:26 vnivanch Exp $ 28 // GEANT4 tag $Name: global-V09-02-10 $ 29 29 // 30 30 //////////////////////////////////////////////////////////////////////// … … 86 86 //////////// 87 87 88 void InsertValues(G4double energy, G4double value);88 void InsertValues(G4double energy, G4double value); 89 89 90 90 G4double GetLowEdgeEnergy(size_t binNumber) const; -
trunk/source/global/management/include/G4PhysicsOrderedFreeVector.icc
r1058 r1193 25 25 // 26 26 // 27 // $Id: G4PhysicsOrderedFreeVector.icc,v 1. 8 2006/06/29 19:02:31 gunterExp $28 // GEANT4 tag $Name: g eant4-09-02-ref-02$27 // $Id: G4PhysicsOrderedFreeVector.icc,v 1.9 2009/06/25 10:05:26 vnivanch Exp $ 28 // GEANT4 tag $Name: global-V09-02-10 $ 29 29 // 30 30 //////////////////////////////////////////////////////////////////////// … … 73 73 void G4PhysicsOrderedFreeVector::DumpValues() 74 74 { 75 for (size_t i = 0; i < numberOf Bin; i++)75 for (size_t i = 0; i < numberOfNodes; i++) 76 76 { 77 77 G4cout << binVector[i] << "\t" << dataVector[i] << G4endl; … … 83 83 { 84 84 G4int n1 = 0; 85 G4int n2 = numberOf Bin/2;86 G4int n3 = numberOf Bin- 1;85 G4int n2 = numberOfNodes/2; 86 G4int n3 = numberOfNodes - 1; 87 87 while (n1 != n3 - 1) 88 88 { -
trunk/source/global/management/include/G4PhysicsVector.hh
r1058 r1193 25 25 // 26 26 // 27 // $Id: G4PhysicsVector.hh,v 1. 18 2008/09/22 08:26:33 gcosmoExp $28 // GEANT4 tag $Name: g eant4-09-02-ref-02$27 // $Id: G4PhysicsVector.hh,v 1.25 2009/11/04 11:32:43 vnivanch Exp $ 28 // GEANT4 tag $Name: global-V09-02-10 $ 29 29 // 30 30 // … … 51 51 // 09 Mar. 2001, H.Kurashige : Added G4PhysicsVectorType & Store/Retrieve() 52 52 // 02 Apr. 2008, A.Bagulya : Added SplineInterpolation() and SetSpline() 53 // 11 May 2009, V.Ivanchenko : Added ComputeSecondDerivatives 54 // 19 Jun. 2009, V.Ivanchenko : Removed hidden bin 53 55 // 54 56 //--------------------------------------------------------------- … … 84 86 // destructor 85 87 86 inline G4double GetValue(G4double theEnergy, G4bool& isOutRange);88 inline G4double Value(G4double theEnergy); 87 89 // Get the cross-section/energy-loss value corresponding to the 88 90 // given energy. An appropriate interpolation is used to calculate 89 91 // the value. 90 // [Note] isOutRange is not used anymore. This argument is kept 91 // for the compatibility reason. 92 93 inline G4double GetValue(G4double theEnergy, G4bool& isOutRange); 94 // Obolete method to get value, isOutRange is not used anymore. 95 // This method is kept for the compatibility reason. 92 96 93 97 G4int operator==(const G4PhysicsVector &right) const ; 94 98 G4int operator!=(const G4PhysicsVector &right) const ; 99 95 100 inline G4double operator[](const size_t binNumber) const ; 96 101 // Returns simply the value in the bin specified by 'binNumber' 97 // of the dataVector. The boundary check will be Done. If you98 // don't want this check, use the operator (). 102 // of the dataVector. The boundary check will not be done. 103 99 104 inline G4double operator()(const size_t binNumber) const ; 100 105 // Returns simply the value in the bin specified by 'binNumber' 101 // of the dataVector. The boundary check will not be Done. If 102 // you want this check, use the operator []. 103 104 inline void PutValue(size_t binNumber, G4double theValue); 106 // of the dataVector. The boundary check will not be Done. 107 108 inline void PutValue(size_t index, G4double theValue); 105 109 // Put 'theValue' into the bin specified by 'binNumber'. 106 // Take note that the ' binNumber' starts from '0'.107 // To fill the vector, you have beforehand to Construct a vector110 // Take note that the 'index' starts from '0'. 111 // To fill the vector, you have beforehand to construct a vector 108 112 // by the constructor with Emin, Emax, Nbin. 'theValue' should 109 // be the crosssection/energyloss value corresponding to the low 110 // edge energy of the bin specified by 'binNumber'. You can get 111 // the low edge energy value of a bin by GetLowEdgeEnergy(). 113 // be the crosssection/energyloss value corresponding to the 114 // energy of the index. You can get this energy by the next method 115 // or by the old method GetLowEdgeEnergy(). 116 117 void ScaleVector(G4double factorE, G4double factorV); 118 // Scale all values of the vector and second derivatives 119 // by factorV, energies by vectorE. This method may be applied 120 // for example after Retrieve a vector from an external file to 121 // convert values into Geant4 units 122 123 inline G4double Energy(size_t index) const; 124 // Returns simply the value in the energy specified by 'index' 125 // of the energy vector. The boundary check will not be done. 126 // Use this function when you fill physis vector by PutValue(). 127 112 128 virtual G4double GetLowEdgeEnergy(size_t binNumber) const; 129 // Obsolete method 113 130 // Get the energy value at the low edge of the specified bin. 114 131 // Take note that the 'binNumber' starts from '0'. 115 // This value is defined when a physics vector is constructed116 // by a constructor of a derived class. Use this function117 // when you fill physis vector by PutValue(). 132 // This value should be defined before the call. 133 // The boundary check will not be done. 134 118 135 inline size_t GetVectorLength() const; 119 136 // Get the toal length (bin number) of the vector. 137 138 void FillSecondDerivatives(); 139 // Initialise second derivatives for spline keeping 140 // 3d derivative continues - default algorithm 141 142 void ComputeSecDerivatives(); 143 // Initialise second derivatives for spline using algorithm 144 // which garantee only 1st derivative continues 145 // Warning: this method should be called when the vector 146 // is already filled 147 148 void ComputeSecondDerivatives(G4double firstPointDerivative, 149 G4double endPointDerivative); 150 // Initialise second derivatives for spline using 151 // user defined 1st derivatives at edge points 152 // Warning: this method should be called when the vector 153 // is already filled 154 120 155 inline G4bool IsFilledVectorExist() const; 121 156 // Is non-empty physics vector already exist? … … 124 159 // Put a comment to the G4PhysicsVector. This may help to check 125 160 // whether your are accessing to the one you want. 161 126 162 inline const G4String& GetComment() const; 127 163 // Retrieve the comment of the G4PhysicsVector. … … 154 190 G4PhysicsVectorType type; // The type of PhysicsVector (enumerator) 155 191 156 G4double edgeMin; // Lower edge value of the lowest bin 157 G4double edgeMax; // Lower edge value of the highest bin 158 size_t numberOfBin; 192 G4double edgeMin; // Energy of first point 193 G4double edgeMax; // Energy of the last point 194 195 size_t numberOfNodes; 159 196 160 197 G4double lastEnergy; // Cache the last input value … … 163 200 164 201 G4PVDataVector dataVector; // Vector to keep the crossection/energyloss 165 G4PVDataVector binVector; // Vector to keep the low edge value of bin 202 G4PVDataVector binVector; // Vector to keep energy 203 G4PVDataVector secDerivative; // Vector to keep second derivatives 166 204 167 205 private: 206 207 G4bool SplinePossible(); 168 208 169 209 inline G4double LinearInterpolation(); … … 174 214 inline void Interpolation(); 175 215 176 void FillSecondDerivatives();177 // Initialise second derivatives for spline178 179 G4double* secDerivative;180 181 216 G4String comment; 182 217 G4bool useSpline; -
trunk/source/global/management/include/G4PhysicsVector.icc
r1058 r1193 25 25 // 26 26 // 27 // $Id: G4PhysicsVector.icc,v 1. 17 2008/09/22 08:26:33 gcosmoExp $28 // GEANT4 tag $Name: g eant4-09-02-ref-02$27 // $Id: G4PhysicsVector.icc,v 1.23 2009/11/18 11:42:03 vnivanch Exp $ 28 // GEANT4 tag $Name: global-V09-02-10 $ 29 29 // 30 30 // … … 59 59 //--------------------------------------------------------------- 60 60 61 inline 62 G4double G4PhysicsVector::Energy(const size_t binNumber) const 63 { 64 return binVector[binNumber]; 65 } 66 67 //--------------------------------------------------------------- 68 61 69 inline 62 70 size_t G4PhysicsVector::GetVectorLength() const 63 71 { 64 return numberOf Bin;72 return numberOfNodes; 65 73 } 66 74 … … 94 102 // numberOfBin-1. 95 103 96 if( !secDerivative) { FillSecondDerivatives(); }104 if(0 == secDerivative.size() ) { FillSecondDerivatives(); } 97 105 98 106 // check bin value 99 G4double delta = binVector[lastBin+1] - binVector[lastBin]; 100 101 G4double a = (binVector[lastBin+1] - lastEnergy)/delta; 102 G4double b = (lastEnergy - binVector[lastBin])/delta; 107 G4double x1 = binVector[lastBin]; 108 G4double x2 = binVector[lastBin+1]; 109 G4double delta = x2 - x1; 110 111 G4double a = (x2 - lastEnergy)/delta; 112 G4double b = (lastEnergy - x1)/delta; 103 113 104 114 // Final evaluation of cubic spline polynomial for return … … 107 117 108 118 G4double res = a*y1 + b*y2 + 109 ( (a*a*a - a)*secDerivative[lastBin] +110 (b*b*b - b)*secDerivative[lastBin+1])*delta*delta/6.0;119 ( (a*a*a - a)*secDerivative[lastBin] + 120 (b*b*b - b)*secDerivative[lastBin+1] )*delta*delta/6.0; 111 121 112 122 return res; … … 117 127 inline 118 128 G4double G4PhysicsVector::GetValue(G4double theEnergy, G4bool&) 129 { 130 return Value(theEnergy); 131 } 132 133 //--------------------------------------------------------------- 134 135 inline 136 G4double G4PhysicsVector::Value(G4double theEnergy) 119 137 { 120 138 // Use cache for speed up - check if the value 'theEnergy' is same as the … … 134 152 lastValue = dataVector[0]; 135 153 136 } else if(theEnergy < lastEnergy && theEnergy >= binVector[lastBin-1]) {137 lastBin--;138 lastEnergy = theEnergy;139 Interpolation();140 141 154 } else if( theEnergy >= edgeMax ){ 142 lastBin = numberOf Bin-1;155 lastBin = numberOfNodes-2; 143 156 lastEnergy = edgeMax; 144 lastValue = dataVector[ lastBin];157 lastValue = dataVector[numberOfNodes-1]; 145 158 146 159 } else { 147 160 lastBin = FindBinLocation(theEnergy); 161 if(lastBin >= numberOfNodes-1) {lastBin = numberOfNodes-2;} 148 162 lastEnergy = theEnergy; 149 163 Interpolation(); … … 167 181 { 168 182 dataVector[binNumber] = theValue; 169 170 // Fill the bin which is hidden to user with theValue. This is to171 // handle correctly when Energy=theEmax in getValue.172 173 if(binNumber==numberOfBin-1) { dataVector[binNumber+1] = theValue; }174 183 } 175 184 … … 181 190 G4bool status=false; 182 191 183 if(numberOf Bin> 0) { status=true; }192 if(numberOfNodes > 0) { status=true; } 184 193 return status; 185 194 } -
trunk/source/global/management/include/G4String.hh
r1058 r1193 25 25 // 26 26 // 27 // $Id: G4String.hh,v 1.1 0 2008/12/08 14:16:05gcosmo Exp $28 // GEANT4 tag $Name: g eant4-09-02-ref-02$27 // $Id: G4String.hh,v 1.15 2009/08/10 10:18:09 gcosmo Exp $ 28 // GEANT4 tag $Name: global-V09-02-10 $ 29 29 // 30 30 // … … 78 78 inline G4int operator!() const; 79 79 80 inline G4bool operator==( G4String) const;80 inline G4bool operator==(const G4String&) const; 81 81 inline G4bool operator==(const char*) const; 82 inline G4bool operator!=( G4String) const;82 inline G4bool operator!=(const G4String&) const; 83 83 inline G4bool operator!=(const char*) const; 84 84 … … 138 138 inline G4bool operator!=(const char*) const; 139 139 140 //inline G4String operator () (unsigned int, unsigned int);141 140 inline operator const char*() const; 142 141 inline G4SubString operator()(str_size, str_size); … … 160 159 inline G4int last(char) const; 161 160 162 inline G4bool contains( std::string) const;161 inline G4bool contains(const std::string&) const; 163 162 inline G4bool contains(char) const; 164 163 … … 184 183 inline unsigned int hash( caseCompare cmp = exact ) const; 185 184 inline unsigned int stlhash() const; 186 187 // useful for supplying hash functions to template hash collection ctors188 //189 static inline unsigned int hash(const G4String&);190 191 185 }; 192 186 -
trunk/source/global/management/include/G4String.icc
r1058 r1193 25 25 // 26 26 // 27 // $Id: G4String.icc,v 1.1 1 2008/12/08 14:16:05gcosmo Exp $28 // GEANT4 tag $Name: g eant4-09-02-ref-02$27 // $Id: G4String.icc,v 1.19 2009/08/10 10:30:03 gcosmo Exp $ 28 // GEANT4 tag $Name: global-V09-02-10 $ 29 29 // 30 30 // … … 91 91 G4int G4SubString::operator!() const 92 92 { 93 return extent==0? 1 : 0;94 } 95 96 G4bool G4SubString::operator==( G4Strings) const97 { 98 return mystring->substr(mystart,extent) == s ? true:false;93 return (extent==0) ? 1 : 0; 94 } 95 96 G4bool G4SubString::operator==(const G4String& s) const 97 { 98 return (mystring->find(s,mystart,extent) != std::string::npos); 99 99 } 100 100 101 101 G4bool G4SubString::operator==(const char* s) const 102 102 { 103 return mystring->substr(mystart,extent) == std::string(s) ? true:false;104 } 105 106 G4bool G4SubString::operator!=( G4Strings) const107 { 108 return mystring->substr(mystart,extent) != std::string(s) ? true:false;103 return (mystring->find(s,mystart,extent) != std::string::npos); 104 } 105 106 G4bool G4SubString::operator!=(const G4String& s) const 107 { 108 return (mystring->find(s,mystart,extent) == std::string::npos); 109 109 } 110 110 111 111 G4bool G4SubString::operator!=(const char* s) const 112 112 { 113 return mystring->substr(mystart,extent) != std::string(s) ? true:false;113 return (mystring->find(s,mystart,extent) == std::string::npos); 114 114 } 115 115 … … 126 126 G4bool G4SubString::isNull() const 127 127 { 128 return extent==0 ? true : false;128 return (extent==0); 129 129 } 130 130 … … 219 219 G4bool G4String::operator==(const G4String &s) const 220 220 { 221 const std_string *a=this; 222 const std_string *b=&s; 223 return *a==*b; 221 return (std_string::compare(s) == 0); 224 222 } 225 223 226 224 G4bool G4String::operator==(const char* s) const 227 225 { 228 const std_string *a=this; 229 const std_string b=s; 230 return *a==b; 226 return (std_string::compare(s) == 0); 231 227 } 232 228 233 229 G4bool G4String::operator!=(const G4String &s) const 234 230 { 235 const std_string *a=this; 236 const std_string *b=&s; 237 return *a!=*b; 231 return !(*this == s); 238 232 } 239 233 240 234 G4bool G4String::operator!=(const char* s) const 241 235 { 242 const std_string *a=this; 243 const std_string b=s; 244 return *a!=b; 236 return !(*this == s); 245 237 } 246 238 … … 268 260 G4int G4String::compareTo(const char* s, caseCompare mode) 269 261 { 270 if(mode==exact) 271 { return strcmp(c_str(),s); } 272 else 273 { return strcasecompare(c_str(),s); } 262 return (mode==exact) ? strcmp(c_str(),s) 263 : strcasecompare(c_str(),s); 274 264 } 275 265 276 266 G4int G4String::compareTo(const G4String& s, caseCompare mode) 277 267 { 278 if(mode==exact) 279 { return strcmp(c_str(),s.c_str()); } 280 else 281 { return strcasecompare(c_str(),s.c_str()); } 268 return (mode==exact) ? strcmp(c_str(),s.c_str()) 269 : strcasecompare(c_str(),s.c_str()); 282 270 } 283 271 … … 298 286 { 299 287 char tmp[1024]; 300 if ( skipWhite ) { 288 if ( skipWhite ) 289 { 301 290 s >> std::ws; 302 291 s.getline(tmp,1024); 303 292 *this=tmp; 304 293 } 305 else { 294 else 295 { 306 296 s.getline(tmp,1024); 307 297 *this=tmp; … … 325 315 G4String& G4String::remove(str_size n) 326 316 { 327 if(n<size()) 328 { erase(n,size()-n); } 317 if(n<size()) { erase(n,size()-n); } 329 318 return *this; 330 319 } … … 346 335 } 347 336 348 G4bool G4String::contains(std::string s) const 349 { 350 G4int i=std_string::find(s); 351 if(i != G4int(std_string::npos)) 352 { return true; } 353 else 354 { return false; } 337 G4bool G4String::contains(const std::string& s) const 338 { 339 return (std_string::find(s) != std_string::npos); 355 340 } 356 341 357 342 G4bool G4String::contains(char c) const 358 343 { 359 G4int i=std_string::find(c); 360 if(i != G4int(std_string::npos)) 361 { return true; } 362 else 363 { return false; } 344 return (std_string::find(c) != std_string::npos); 364 345 } 365 346 … … 474 455 return str_size(h); 475 456 } 476 477 unsigned int G4String::hash(const G4String& s)478 {479 return s.hash();480 } -
trunk/source/global/management/include/G4Version.hh
r1058 r1193 25 25 // 26 26 // 27 // $Id: G4Version.hh,v 1. 17 2008/12/02 09:15:58 gcosmoExp $28 // GEANT4 tag $Name: g eant4-09-02-ref-02$27 // $Id: G4Version.hh,v 1.24 2009/06/19 14:29:09 vnivanch Exp $ 28 // GEANT4 tag $Name: global-V09-02-10 $ 29 29 // 30 30 // Version information … … 47 47 48 48 #ifndef G4VERSION_NUMBER 49 #define G4VERSION_NUMBER 9 2049 #define G4VERSION_NUMBER 930 50 50 #endif 51 51 52 52 #ifndef G4VERSION_TAG 53 #define G4VERSION_TAG "$Name: g eant4-09-02-ref-02$"53 #define G4VERSION_TAG "$Name: global-V09-02-10 $" 54 54 #endif 55 55 … … 58 58 #include "G4String.hh" 59 59 60 static const G4String G4Version = "$Name: g eant4-09-02-ref-02$";61 static const G4String G4Date = "( 19-December-2008)";60 static const G4String G4Version = "$Name: global-V09-02-10 $"; 61 static const G4String G4Date = "(5-June-2009)"; 62 62 63 63 #endif -
trunk/source/global/management/src/G4ErrorPropagatorData.cc
r1058 r1193 25 25 // 26 26 // 27 // $Id: G4ErrorPropagatorData.cc,v 1. 3 2007/06/05 13:04:30gcosmo Exp $28 // GEANT4 tag $Name: g eant4-09-02-ref-02$27 // $Id: G4ErrorPropagatorData.cc,v 1.5 2009/05/19 13:31:47 gcosmo Exp $ 28 // GEANT4 tag $Name: global-V09-02-10 $ 29 29 // 30 30 // … … 43 43 //------------------------------------------------------------------- 44 44 45 G4ErrorPropagatorData::G4ErrorPropagatorData(){} 46 G4ErrorPropagatorData::~G4ErrorPropagatorData(){} 45 G4ErrorPropagatorData::~G4ErrorPropagatorData() 46 { 47 } 48 49 G4ErrorPropagatorData::G4ErrorPropagatorData() 50 { 51 theStage = G4ErrorStage_Inflation; 52 } 47 53 48 54 G4ErrorPropagatorData* G4ErrorPropagatorData::GetErrorPropagatorData() -
trunk/source/global/management/src/G4LPhysicsFreeVector.cc
r1058 r1193 25 25 // 26 26 // 27 // $Id: G4LPhysicsFreeVector.cc,v 1.2 2 2008/09/22 08:26:33 gcosmoExp $28 // GEANT4 tag $Name: g eant4-09-02-ref-02$27 // $Id: G4LPhysicsFreeVector.cc,v 1.23 2009/06/25 10:05:26 vnivanch Exp $ 28 // GEANT4 tag $Name: global-V09-02-10 $ 29 29 // 30 30 // … … 35 35 // 36 36 // F.W. Jones, TRIUMF, 04-JUN-96 37 // 38 // Modified: 39 // 19 Jun. 2009, V.Ivanchenko : removed hidden bin 37 40 // 38 41 // -------------------------------------------------------------------- 39 42 40 43 #include "G4LPhysicsFreeVector.hh" 41 42 44 #include "G4ios.hh" 43 45 … … 61 63 edgeMin = binmin; 62 64 edgeMax = binmax; 63 numberOf Bin= nbin;64 binVector.reserve(n bin+1);65 dataVector.reserve(n bin+1);66 for (size_t i=0; i< =numberOfBin; i++)65 numberOfNodes = nbin; 66 binVector.reserve(numberOfNodes); 67 dataVector.reserve(numberOfNodes); 68 for (size_t i=0; i<numberOfNodes; i++) 67 69 { 68 70 binVector.push_back(0.0); … … 105 107 void G4LPhysicsFreeVector::DumpValues() 106 108 { 107 for (size_t i = 0; i < numberOf Bin; i++)109 for (size_t i = 0; i < numberOfNodes; i++) 108 110 { 109 111 G4cout << binVector[i] << " " << dataVector[i]/millibarn << G4endl; -
trunk/source/global/management/src/G4PhysicsFreeVector.cc
r1058 r1193 25 25 // 26 26 // 27 // $Id: G4PhysicsFreeVector.cc,v 1.1 2 2008/09/22 14:49:57 gcosmoExp $28 // GEANT4 tag $Name: g eant4-09-02-ref-02$27 // $Id: G4PhysicsFreeVector.cc,v 1.13 2009/06/25 10:05:26 vnivanch Exp $ 28 // GEANT4 tag $Name: global-V09-02-10 $ 29 29 // 30 30 // … … 41 41 // 26 Sep. 1996, K.Amako : Constructor with only 'bin size' added 42 42 // 11 Nov. 2000, H.Kurashige : use STL vector for dataVector and binVector 43 // 19 Jun. 2009, V.Ivanchenko : removed hidden bin 43 44 // 44 45 //-------------------------------------------------------------------- … … 57 58 { 58 59 type = T_G4PhysicsFreeVector; 59 numberOf Bin= theNbin;60 numberOfNodes = theNbin; 60 61 61 // Add extra one bin (hidden to user) to handle correctly when 62 // Energy=theEmax in getValue. 63 dataVector.reserve(numberOfBin+1); 64 binVector.reserve(numberOfBin+1); 62 dataVector.reserve(numberOfNodes); 63 binVector.reserve(numberOfNodes); 65 64 66 for (size_t i=0; i< =numberOfBin; i++)65 for (size_t i=0; i<numberOfNodes; i++) 67 66 { 68 67 binVector.push_back(0.0); … … 75 74 { 76 75 type = T_G4PhysicsFreeVector; 77 numberOf Bin= theBinVector.size();76 numberOfNodes = theBinVector.size(); 78 77 79 // Add extra one bin (hidden to user) to handle correctly when 80 // Energy=theEmax in getValue. 81 dataVector.reserve(numberOfBin+1); 82 binVector.reserve(numberOfBin+1); 78 dataVector.reserve(numberOfNodes); 79 binVector.reserve(numberOfNodes); 83 80 84 for (size_t i=0; i<numberOf Bin; i++)81 for (size_t i=0; i<numberOfNodes; i++) 85 82 { 86 83 binVector.push_back(theBinVector[i]); … … 88 85 } 89 86 90 // Put values to extra hidden bin. For 'binVector', the 'edgeMin' of the91 // extra hidden bin is assumed to have the following value. For binary92 // search, this value is completely arbitrary if it is greater than93 // the 'edgeMin' at 'numberOfBin-1'.94 binVector.push_back ( theBinVector[numberOfBin-1] + 1.0 );95 96 97 // Put values to extra hidden bin. For 'dataVector', the 'value' of the98 // extra hidden bin is assumed to have the same as the one at 'numberBin-1'.99 dataVector.push_back( theDataVector[numberOfBin-1] );100 101 87 edgeMin = binVector[0]; 102 edgeMax = binVector[numberOf Bin-1];88 edgeMax = binVector[numberOfNodes-1]; 103 89 } 104 90 … … 113 99 dataVector[theBinNumber] = theDataValue; 114 100 115 if( theBinNumber == numberOf Bin-1 )101 if( theBinNumber == numberOfNodes-1 ) 116 102 { 117 edgeMax = binVector[numberOfBin-1]; 118 119 // Put values to extra hidden bin. For 'binVector', the 'edgeMin' 120 // of the extra hidden bin is assumed to have the following value. 121 // For binary search, this value is completely arbitrary if it is 122 // greater than the 'edgeMin' at 'numberOfBin-1'. 123 binVector[numberOfBin] = theBinValue + 1.0; 124 125 // Put values to extra hidden bin. For 'dataVector', the 'value' 126 // of the extra hidden bin is assumed to have the same as the one 127 // at 'numberBin-1'. 128 dataVector[numberOfBin] = theDataValue; 103 edgeMax = binVector[numberOfNodes-1]; 129 104 } 130 105 -
trunk/source/global/management/src/G4PhysicsLinearVector.cc
r1058 r1193 25 25 // 26 26 // 27 // $Id: G4PhysicsLinearVector.cc,v 1.1 4 2008/09/22 14:49:57 gcosmoExp $28 // GEANT4 tag $Name: g eant4-09-02-ref-02$27 // $Id: G4PhysicsLinearVector.cc,v 1.15 2009/06/25 10:05:26 vnivanch Exp $ 28 // GEANT4 tag $Name: global-V09-02-10 $ 29 29 // 30 30 // … … 35 35 // 36 36 // 15 Feb 1996 - K.Amako : 1st version 37 // 19 Jun 2009 - V.Ivanchenko : removed hidden bin 37 38 // 38 39 //-------------------------------------------------------------------- … … 51 52 type = T_G4PhysicsLinearVector; 52 53 53 // Add extra one bin (hidden to user) to handle correctly when 54 // Energy=theEmax in getValue. 55 dataVector.reserve(theNbin+1); 56 binVector.reserve(theNbin+1); 54 numberOfNodes = theNbin + 1; 55 dataVector.reserve(numberOfNodes); 56 binVector.reserve(numberOfNodes); 57 57 58 numberOfBin = theNbin; 59 60 for (size_t i=0; i<=numberOfBin; i++) 58 for (size_t i=0; i<numberOfNodes; i++) 61 59 { 62 60 binVector.push_back(0.0); … … 73 71 type = T_G4PhysicsLinearVector; 74 72 75 // Add extra one bin (hidden to user) to handle correctly when 76 // Energy=theEmax in getValue. 77 dataVector.reserve(theNbin+1); 78 binVector.reserve(theNbin+1); 73 numberOfNodes = theNbin + 1; 74 dataVector.reserve(numberOfNodes); 75 binVector.reserve(numberOfNodes); 79 76 80 numberOfBin = theNbin; 81 82 for (size_t i=0; i<numberOfBin+1; i++) 77 for (size_t i=0; i<numberOfNodes; i++) 83 78 { 84 79 binVector.push_back( theEmin + i*dBin ); … … 87 82 88 83 edgeMin = binVector[0]; 89 edgeMax = binVector[numberOf Bin-1];84 edgeMax = binVector[numberOfNodes-1]; 90 85 91 86 } -
trunk/source/global/management/src/G4PhysicsLnVector.cc
r1058 r1193 25 25 // 26 26 // 27 // $Id: G4PhysicsLnVector.cc,v 1.1 7 2008/09/22 14:49:57 gcosmoExp $28 // GEANT4 tag $Name: g eant4-09-02-ref-02$27 // $Id: G4PhysicsLnVector.cc,v 1.18 2009/06/25 10:05:26 vnivanch Exp $ 28 // GEANT4 tag $Name: global-V09-02-10 $ 29 29 // 30 30 // … … 35 35 // 36 36 // 27 Apr 1999 - M.G.Pia: Created from G4PhysicsLogVector 37 // 19 Jun 2009 - V.Ivanchenko : removed hidden bin 37 38 // 38 39 // -------------------------------------------------------------- … … 51 52 type = T_G4PhysicsLnVector; 52 53 53 // Add extra one bin (hidden to user) to handle correctly when 54 // Energy=theEmax in getValue. 55 dataVector.reserve(theNbin+1); 56 binVector.reserve(theNbin+1); 54 numberOfNodes = theNbin + 1; 55 dataVector.reserve(numberOfNodes); 56 binVector.reserve(numberOfNodes); 57 57 58 numberOfBin = theNbin; 59 60 for (size_t i=0; i<=numberOfBin; i++) 58 for (size_t i=0; i<numberOfNodes; i++) 61 59 { 62 60 binVector.push_back(0.0); … … 73 71 type = T_G4PhysicsLnVector; 74 72 75 // Add extra one bin (hidden to user) to handle correctly when 76 // Energy=theEmax in getValue. 77 dataVector.reserve(theNbin+1); 78 binVector.reserve(theNbin+1); 73 numberOfNodes = theNbin + 1; 74 dataVector.reserve(numberOfNodes); 75 binVector.reserve(numberOfNodes); 79 76 80 numberOfBin = theNbin; 81 82 for (size_t i=0; i<numberOfBin+1; i++) 77 for (size_t i=0; i<numberOfNodes; i++) 83 78 { 84 binVector.push_back(std::exp( std::log(theEmin)+i*dBin));79 binVector.push_back(std::exp((baseBin+i)*dBin)); 85 80 dataVector.push_back(0.0); 86 81 } 87 82 88 83 edgeMin = binVector[0]; 89 edgeMax = binVector[numberOf Bin-1];84 edgeMax = binVector[numberOfNodes-1]; 90 85 } 91 86 -
trunk/source/global/management/src/G4PhysicsLogVector.cc
r1058 r1193 25 25 // 26 26 // 27 // $Id: G4PhysicsLogVector.cc,v 1.2 1 2008/09/22 08:26:33 gcosmoExp $28 // GEANT4 tag $Name: g eant4-09-02-ref-02$27 // $Id: G4PhysicsLogVector.cc,v 1.22 2009/06/25 10:05:26 vnivanch Exp $ 28 // GEANT4 tag $Name: global-V09-02-10 $ 29 29 // 30 30 // … … 42 42 // 9 Mar. 2001, H.Kurashige : added PhysicsVector type and Retrieve 43 43 // 05 Sep. 2008, V.Ivanchenko : added protections for zero-length vector 44 // 19 Jun. 2009, V.Ivanchenko : removed hidden bin 44 45 // 45 46 // -------------------------------------------------------------- … … 58 59 type = T_G4PhysicsLogVector; 59 60 60 // Add extra one bin (hidden to user) to handle correctly when 61 // Energy=theEmax in getValue. 62 dataVector.reserve(theNbin+1); 63 binVector.reserve(theNbin+1); 61 numberOfNodes = theNbin + 1; 62 dataVector.reserve(numberOfNodes); 63 binVector.reserve(numberOfNodes); 64 64 65 numberOfBin = theNbin; 66 67 for (size_t i=0; i<=numberOfBin; i++) 65 for (size_t i=0; i<numberOfNodes; i++) 68 66 { 69 67 binVector.push_back(0.0); … … 79 77 type = T_G4PhysicsLogVector; 80 78 81 // Add extra one bin (hidden to user) to handle correctly when82 // Energy=theEmax in getValue.83 dataVector.reserve(theNbin+1);84 binVector.reserve(theNbin+1);79 numberOfNodes = theNbin + 1; 80 dataVector.reserve(numberOfNodes); 81 binVector.reserve(numberOfNodes); 82 static const G4double g4log10 = std::log(10.); 85 83 86 numberOfBin = theNbin; 87 88 for (size_t i=0; i<numberOfBin+1; i++) 84 for (size_t i=0; i<numberOfNodes; i++) 89 85 { 90 binVector.push_back(std:: pow(10., std::log10(theEmin)+i*dBin));86 binVector.push_back(std::exp(g4log10*(baseBin+i)*dBin)); 91 87 dataVector.push_back(0.0); 92 88 } 89 93 90 edgeMin = binVector[0]; 94 edgeMax = binVector[numberOf Bin-1];91 edgeMax = binVector[numberOfNodes-1]; 95 92 } 96 93 -
trunk/source/global/management/src/G4PhysicsOrderedFreeVector.cc
r1058 r1193 25 25 // 26 26 // 27 // $Id: G4PhysicsOrderedFreeVector.cc,v 1.1 2 2008/09/22 14:49:57 gcosmoExp $28 // GEANT4 tag $Name: g eant4-09-02-ref-02$27 // $Id: G4PhysicsOrderedFreeVector.cc,v 1.13 2009/06/25 10:05:26 vnivanch Exp $ 28 // GEANT4 tag $Name: global-V09-02-10 $ 29 29 // 30 30 //////////////////////////////////////////////////////////////////////// … … 43 43 // 2000-11-11 by H.Kurashige 44 44 // > use STL vector for dataVector and binVector 45 // 19 Jun. 2009-06-19 by V.Ivanchenko 46 // > removed hidden bin 47 // 45 48 // mail: gum@triumf.ca 46 //47 49 // 48 50 //////////////////////////////////////////////////////////////////////// … … 65 67 type = T_G4PhysicsOrderedFreeVector; 66 68 67 dataVector.reserve(VectorLength +1);68 binVector.reserve(VectorLength +1);69 numberOf Bin= VectorLength;69 dataVector.reserve(VectorLength); 70 binVector.reserve(VectorLength); 71 numberOfNodes = VectorLength; 70 72 71 73 for (size_t i = 0 ; i < VectorLength ; i++) … … 74 76 dataVector.push_back(Values[i]); 75 77 } 76 edgeMin = binVector.front(); 77 edgeMax = binVector.back(); 78 binVector.push_back ( binVector[numberOfBin-1] + 1.0 ); 79 dataVector.push_back( dataVector[numberOfBin-1] ); 78 edgeMin = binVector[0]; 79 edgeMax = binVector[numberOfNodes-1]; 80 80 } 81 81 … … 101 101 binVector.push_back(energy); 102 102 dataVector.push_back(value); 103 numberOf Bin++;103 numberOfNodes++; 104 104 edgeMin = binVector.front(); 105 105 edgeMax = binVector.back(); … … 133 133 { 134 134 G4int n1 = 0; 135 G4int n2 = numberOf Bin/2;136 G4int n3 = numberOf Bin- 1;135 G4int n2 = numberOfNodes/2; 136 G4int n3 = numberOfNodes - 1; 137 137 while (n1 != n3 - 1) { 138 138 if (aValue > dataVector[n2]) -
trunk/source/global/management/src/G4PhysicsVector.cc
r1058 r1193 25 25 // 26 26 // 27 // $Id: G4PhysicsVector.cc,v 1. 27 2008/10/16 12:14:36 gcosmoExp $28 // GEANT4 tag $Name: g eant4-09-02-ref-02$27 // $Id: G4PhysicsVector.cc,v 1.40 2009/11/04 11:45:54 vnivanch Exp $ 28 // GEANT4 tag $Name: global-V09-02-10 $ 29 29 // 30 30 // … … 43 43 // 09 Mar. 2001, H.Kurashige : added G4PhysicsVector type 44 44 // 05 Sep. 2008, V.Ivanchenko : added protections for zero-length vector 45 // 11 May 2009, A.Bagulya : added new implementation of methods 46 // ComputeSecondDerivatives - first derivatives at edge points 47 // should be provided by a user 48 // FillSecondDerivatives - default computation base on "not-a-knot" 49 // algorithm 50 // 19 Jun. 2009, V.Ivanchenko : removed hidden bin 45 51 // -------------------------------------------------------------- 46 52 … … 52 58 G4PhysicsVector::G4PhysicsVector(G4bool spline) 53 59 : type(T_G4PhysicsVector), 54 edgeMin(0.), edgeMax(0.), numberOfBin(0), 55 lastEnergy(0.), lastValue(0.), lastBin(0), 56 secDerivative(0), useSpline(spline) 60 edgeMin(0.), edgeMax(0.), numberOfNodes(0), 61 lastEnergy(-1.0), lastValue(0.), lastBin(0), useSpline(spline) 57 62 {} 58 63 … … 60 65 61 66 G4PhysicsVector::~G4PhysicsVector() 62 { 63 DeleteData(); 64 } 67 {} 65 68 66 69 // -------------------------------------------------------------- … … 78 81 if (type != right.type) { return *this; } 79 82 80 DeleteData();83 //DeleteData(); 81 84 CopyData(right); 82 85 … … 102 105 void G4PhysicsVector::DeleteData() 103 106 { 104 delete [] secDerivative; 105 secDerivative = 0; 107 secDerivative.clear(); 106 108 } 107 109 … … 113 115 edgeMin = vec.edgeMin; 114 116 edgeMax = vec.edgeMax; 115 numberOf Bin = vec.numberOfBin;117 numberOfNodes = vec.numberOfNodes; 116 118 lastEnergy = vec.lastEnergy; 117 119 lastValue = vec.lastValue; … … 121 123 useSpline = vec.useSpline; 122 124 comment = vec.comment; 123 if (vec.secDerivative) 124 { 125 secDerivative = new G4double [numberOfBin]; 126 for (size_t i=0; i<numberOfBin; i++) 127 { 128 secDerivative[i] = vec.secDerivative[i]; 129 } 130 } 131 else 132 { 133 secDerivative = 0; 134 } 125 secDerivative = vec.secDerivative; 135 126 } 136 127 … … 157 148 fOut.write((char*)(&edgeMin), sizeof edgeMin); 158 149 fOut.write((char*)(&edgeMax), sizeof edgeMax); 159 fOut.write((char*)(&numberOf Bin), sizeof numberOfBin);150 fOut.write((char*)(&numberOfNodes), sizeof numberOfNodes); 160 151 161 152 // contents … … 164 155 165 156 G4double* value = new G4double[2*size]; 166 for(size_t i = 0; i < size; i++)157 for(size_t i = 0; i < size; ++i) 167 158 { 168 159 value[2*i] = binVector[i]; … … 180 171 { 181 172 // clear properties; 182 lastEnergy= 0.;173 lastEnergy=-1.0; 183 174 lastValue =0.; 184 175 lastBin =0; 185 176 dataVector.clear(); 186 177 binVector.clear(); 178 secDerivative.clear(); 187 179 comment = ""; 188 180 … … 191 183 { 192 184 // binning 193 fIn >> edgeMin >> edgeMax >> numberOf Bin;185 fIn >> edgeMin >> edgeMax >> numberOfNodes; 194 186 if (fIn.fail()) { return false; } 195 187 // contents … … 218 210 fIn.read((char*)(&edgeMin), sizeof edgeMin); 219 211 fIn.read((char*)(&edgeMax), sizeof edgeMax); 220 fIn.read((char*)(&numberOf Bin), sizeof numberOfBin);212 fIn.read((char*)(&numberOfNodes), sizeof numberOfNodes ); 221 213 222 214 // contents … … 234 226 binVector.reserve(size); 235 227 dataVector.reserve(size); 236 for(size_t i = 0; i < size; i++)228 for(size_t i = 0; i < size; ++i) 237 229 { 238 230 binVector.push_back(value[2*i]); … … 245 237 // -------------------------------------------------------------- 246 238 239 void 240 G4PhysicsVector::ScaleVector(G4double factorE, G4double factorV) 241 { 242 size_t n = dataVector.size(); 243 size_t i; 244 if(n > 0) { 245 for(i=0; i<n; ++i) { 246 binVector[i] *= factorE; 247 dataVector[i] *= factorV; 248 } 249 } 250 n = secDerivative.size(); 251 if(n > 0) { for(i=0; i<n; ++i) { secDerivative[i] *= factorV; } } 252 253 edgeMin *= factorE; 254 edgeMax *= factorE; 255 lastEnergy *= factorE; 256 lastValue *= factorV; 257 } 258 259 // -------------------------------------------------------------- 260 261 void 262 G4PhysicsVector::ComputeSecondDerivatives(G4double firstPointDerivative, 263 G4double endPointDerivative) 264 // A standard method of computation of second derivatives 265 // First derivatives at the first and the last point should be provided 266 // See for example W.H. Press et al. "Numerical reciptes and C" 267 // Cambridge University Press, 1997. 268 { 269 if(4 > numberOfNodes) // cannot compute derivatives for less than 4 bins 270 { 271 ComputeSecDerivatives(); 272 return; 273 } 274 275 if(!SplinePossible()) { return; } 276 277 G4int n = numberOfNodes-1; 278 279 G4double* u = new G4double [n]; 280 281 G4double p, sig, un; 282 283 u[0] = (6.0/(binVector[1]-binVector[0])) 284 * ((dataVector[1]-dataVector[0])/(binVector[1]-binVector[0]) 285 - firstPointDerivative); 286 287 secDerivative[0] = - 0.5; 288 289 // Decomposition loop for tridiagonal algorithm. secDerivative[i] 290 // and u[i] are used for temporary storage of the decomposed factors. 291 292 for(G4int i=1; i<n; ++i) 293 { 294 sig = (binVector[i]-binVector[i-1]) / (binVector[i+1]-binVector[i-1]); 295 p = sig*secDerivative[i-1] + 2.0; 296 secDerivative[i] = (sig - 1.0)/p; 297 u[i] = (dataVector[i+1]-dataVector[i])/(binVector[i+1]-binVector[i]) 298 - (dataVector[i]-dataVector[i-1])/(binVector[i]-binVector[i-1]); 299 u[i] = 6.0*u[i]/(binVector[i+1]-binVector[i-1]) - sig*u[i-1]/p; 300 } 301 302 sig = (binVector[n-1]-binVector[n-2]) / (binVector[n]-binVector[n-2]); 303 p = sig*secDerivative[n-2] + 2.0; 304 un = (6.0/(binVector[n]-binVector[n-1])) 305 *(endPointDerivative - 306 (dataVector[n]-dataVector[n-1])/(binVector[n]-binVector[n-1])) - u[n-1]/p; 307 secDerivative[n] = un/(secDerivative[n-1] + 2.0); 308 309 // The back-substitution loop for the triagonal algorithm of solving 310 // a linear system of equations. 311 312 for(G4int k=n-1; k>0; --k) 313 { 314 secDerivative[k] *= 315 (secDerivative[k+1] - 316 u[k]*(binVector[k+1]-binVector[k-1])/(binVector[k+1]-binVector[k])); 317 } 318 secDerivative[0] = 0.5*(u[0] - secDerivative[1]); 319 320 delete [] u; 321 } 322 323 // -------------------------------------------------------------- 324 247 325 void G4PhysicsVector::FillSecondDerivatives() 326 // Computation of second derivatives using "Not-a-knot" endpoint conditions 327 // B.I. Kvasov "Methods of shape-preserving spline approximation" 328 // World Scientific, 2000 248 329 { 249 secDerivative = new G4double [numberOfBin]; 250 251 size_t n = numberOfBin-1; 252 253 // cannot compute derivatives for less than 3 points 254 if(3 > numberOfBin) 255 { 256 secDerivative[0] = 0.0; 257 secDerivative[n] = 0.0; 330 if(5 > numberOfNodes) // cannot compute derivatives for less than 4 points 331 { 332 ComputeSecDerivatives(); 258 333 return; 259 334 } 260 335 261 for(size_t i=1; i<n; i++) 262 { 263 secDerivative[i] = 336 if(!SplinePossible()) { return; } 337 338 G4int n = numberOfNodes-1; 339 340 //G4cout << "G4PhysicsVector::FillSecondDerivatives() n= " << n << G4endl; 341 // G4cout << *this << G4endl; 342 343 G4double* u = new G4double [n]; 344 345 G4double p, sig; 346 347 u[1] = ((dataVector[2]-dataVector[1])/(binVector[2]-binVector[1]) - 348 (dataVector[1]-dataVector[0])/(binVector[1]-binVector[0])); 349 u[1] = 6.0*u[1]*(binVector[2]-binVector[1]) 350 / ((binVector[2]-binVector[0])*(binVector[2]-binVector[0])); 351 352 // Decomposition loop for tridiagonal algorithm. secDerivative[i] 353 // and u[i] are used for temporary storage of the decomposed factors. 354 355 secDerivative[1] = (2.0*binVector[1]-binVector[0]-binVector[2]) 356 / (2.0*binVector[2]-binVector[0]-binVector[1]); 357 358 for(G4int i=2; i<n-1; ++i) 359 { 360 sig = (binVector[i]-binVector[i-1]) / (binVector[i+1]-binVector[i-1]); 361 p = sig*secDerivative[i-1] + 2.0; 362 secDerivative[i] = (sig - 1.0)/p; 363 u[i] = (dataVector[i+1]-dataVector[i])/(binVector[i+1]-binVector[i]) 364 - (dataVector[i]-dataVector[i-1])/(binVector[i]-binVector[i-1]); 365 u[i] = (6.0*u[i]/(binVector[i+1]-binVector[i-1])) - sig*u[i-1]/p; 366 } 367 368 sig = (binVector[n-1]-binVector[n-2]) / (binVector[n]-binVector[n-2]); 369 p = sig*secDerivative[n-3] + 2.0; 370 u[n-1] = (dataVector[n]-dataVector[n-1])/(binVector[n]-binVector[n-1]) 371 - (dataVector[n-1]-dataVector[n-2])/(binVector[n-1]-binVector[n-2]); 372 u[n-1] = 6.0*sig*u[n-1]/(binVector[n]-binVector[n-2]) 373 - (2.0*sig - 1.0)*u[n-2]/p; 374 375 p = (1.0+sig) + (2.0*sig-1.0)*secDerivative[n-2]; 376 secDerivative[n-1] = u[n-1]/p; 377 378 // The back-substitution loop for the triagonal algorithm of solving 379 // a linear system of equations. 380 381 for(G4int k=n-2; k>1; --k) 382 { 383 secDerivative[k] *= 384 (secDerivative[k+1] - 385 u[k]*(binVector[k+1]-binVector[k-1])/(binVector[k+1]-binVector[k])); 386 } 387 secDerivative[n] = (secDerivative[n-1] - (1.0-sig)*secDerivative[n-2])/sig; 388 sig = 1.0 - ((binVector[2]-binVector[1])/(binVector[2]-binVector[0])); 389 secDerivative[1] *= (secDerivative[2] - u[1]/(1.0-sig)); 390 secDerivative[0] = (secDerivative[1] - sig*secDerivative[2])/(1.0-sig); 391 392 delete [] u; 393 } 394 395 // -------------------------------------------------------------- 396 397 void 398 G4PhysicsVector::ComputeSecDerivatives() 399 // A simplified method of computation of second derivatives 400 { 401 if(!SplinePossible()) { return; } 402 403 if(3 > numberOfNodes) // cannot compute derivatives for less than 4 bins 404 { 405 useSpline = false; 406 return; 407 } 408 409 size_t n = numberOfNodes-1; 410 411 for(size_t i=1; i<n; ++i) 412 { 413 secDerivative[i] = 264 414 3.0*((dataVector[i+1]-dataVector[i])/(binVector[i+1]-binVector[i]) - 265 415 (dataVector[i]-dataVector[i-1])/(binVector[i]-binVector[i-1])) … … 269 419 secDerivative[0] = secDerivative[1]; 270 420 } 421 422 // -------------------------------------------------------------- 423 424 G4bool G4PhysicsVector::SplinePossible() 425 // Initialise second derivative array. If neighbor energy coincide 426 // or not ordered than spline cannot be applied 427 { 428 if(!useSpline) return useSpline; 429 secDerivative.clear(); 430 secDerivative.reserve(numberOfNodes); 431 for(size_t j=0; j<numberOfNodes; ++j) 432 { 433 secDerivative.push_back(0.0); 434 if(j > 0) 435 { 436 if(binVector[j]-binVector[j-1] <= 0.) { useSpline = false; } 437 } 438 } 439 return useSpline; 440 } 271 441 272 442 // -------------------------------------------------------------- … … 276 446 // binning 277 447 out << std::setprecision(12) << pv.edgeMin; 278 out <<" " << pv.edgeMax <<" " << pv.numberOf Bin<< G4endl;448 out <<" " << pv.edgeMax <<" " << pv.numberOfNodes << G4endl; 279 449 280 450 // contents
Note: See TracChangeset
for help on using the changeset viewer.