Changeset 1193 for trunk/source/global


Ignore:
Timestamp:
Nov 19, 2009, 3:40:00 PM (15 years ago)
Author:
garnier
Message:

CVS update

Location:
trunk/source/global
Files:
21 edited

Legend:

Unmodified
Added
Removed
  • trunk/source/global/History

    r1058 r1193  
    1 $Id: History,v 1.217 2009/01/27 08:49:10 gcosmo Exp $
     1$Id: History,v 1.232 2009/11/18 11:42:03 vnivanch Exp $
    22-------------------------------------------------------------------
    33
     
    1717     * Reverse chronological order (last date on top), please *
    1818     ----------------------------------------------------------
     19
     20November 17, 2009  V.Ivanchenko (global-V09-02-10)
     21- G4PhysicsVector - fixed Value method for the case when energy is
     22                    around edgeMax adding extra protection
     23
     24November 4, 2009  V.Ivanchenko (global-V09-02-09)
     25- G4PhysicsVector: cleaned up initialisation of cached values;
     26  simplified Value() method, now not looking into previous bin;
     27  added method ScaleVector() needed for ICRU'73 data.
     28
     29October 29, 2009  G.Cosmo (global-V09-02-08)
     30- Cleared warning for unused argument in G4Allocator internal method.
     31
     32August 7, 2009  G.Cosmo (global-V09-02-07)
     33- Some improvements to G4String and G4SubString implementation of operators
     34  and comparison stub functions, to reduce generation of temporaries.
     35- Removed obsolete static hash(s) method.
     36
     37July 3, 2009  V.Ivanchenko (global-V09-02-06)
     38- New class G4Pow, a singleton for fast computation of log() and pow() for
     39  integer and double arguments in the interval [1-256].
     40- G4PhysicsVector:
     41  o Added ComputeSecDerivative() and SplinePossible() methods as a
     42    computation of Spline coefficients for small number of bins.
     43  o Removed old obsolete hidden bin approach.
     44  o Use std::vector for second derivatives instead of C arrays.
     45
     46May 26, 2009  A.Bagulya (global-V09-02-05)
     47- G4PhysicsVector: do not change flag "useSpline" during initialisation
     48  to address bug in N02 reported by valgrind
     49
     50May 26, 2009  A.Bagulya (global-V09-02-04)
     51- G4PhysicsVector: set limits of number of bins for spline: 4 for the case
     52  of first derivatives defined; 5 for the case when "Not-a-Knot" algorithm
     53  is used.
     54
     55May 15, 2009  P.Arce (global-V09-02-03)
     56- Added G4ErrorStage enum to G4ErrorPropagatorData class.
     57
     58May 13, 2009  A.Bagulya (global-V09-02-02)
     59- G4PhysicsVector: added method ComputeSecondDerivatives() for the case when
     60  user provides first derivative at endpoits; use "Not-a-Knot" algorithm for
     61  the computation of second derivatives in default method
     62  FillSecondDerivatives().
     63- Changed date for release 9.3-beta.
     64
     65Mar 16, 2009  G.Cosmo (global-V09-02-01)
     66- Changed date for release 9.2.p01.
    1967
    2068Jan 27, 2009  G.Cosmo (global-V09-02-00)
  • trunk/source/global/management/include/G4Allocator.hh

    r1058 r1193  
    2525//
    2626//
    27 // $Id: G4Allocator.hh,v 1.18 2006/06/29 19:01:16 gunter Exp $
    28 // GEANT4 tag $Name: geant4-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 $
    2929//
    3030//
     
    9999      // Returns the address of values
    100100
    101     pointer allocate(size_type n, void* hint = 0)
     101    pointer allocate(size_type n, void* = 0)
    102102    {
    103103      // Allocates space for n elements of type Type, but does not initialise
  • trunk/source/global/management/include/G4ErrorPropagatorData.hh

    r1058 r1193  
    2525//
    2626//
    27 // $Id: G4ErrorPropagatorData.hh,v 1.3 2007/06/08 10:33:47 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-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 $
    2929//
    3030//
     
    3838//   and manager verbosity for the error propagation classes.
    3939
    40 // History:
    4140// - Created. P.Arce, 2004.
    4241// --------------------------------------------------------------------
     
    5655                    G4ErrorState_TargetCloserThanBoundary,
    5756                    G4ErrorState_StoppedAtTarget };
     57
     58enum G4ErrorStage  { G4ErrorStage_Inflation = 1,
     59                     G4ErrorStage_Deflation };
     60
    5861class G4ErrorTarget;
    5962
     
    7275  G4ErrorState GetState() const;
    7376  void SetState( G4ErrorState sta );
     77
     78  G4ErrorStage GetStage() const;
     79  void SetStage( G4ErrorStage sta );
    7480
    7581  const G4ErrorTarget* GetTarget( G4bool mustExist = 0) const;
     
    94100  G4ErrorState theState;
    95101
     102  G4ErrorStage theStage;
     103
    96104  G4ErrorTarget* theTarget;
    97105
  • trunk/source/global/management/include/G4ErrorPropagatorData.icc

    r1058 r1193  
    2525//
    2626//
    27 // $Id: G4ErrorPropagatorData.icc,v 1.4 2007/06/08 10:33:47 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-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 $
    2929//
    3030//
     
    5959
    6060inline
     61void G4ErrorPropagatorData::SetStage( G4ErrorStage sta )
     62{
     63  theStage = sta;
     64}
     65
     66inline
     67G4ErrorStage G4ErrorPropagatorData::GetStage() const
     68{
     69  return theStage;
     70}
     71
     72inline
    6173const G4ErrorTarget* G4ErrorPropagatorData::GetTarget( G4bool mustExist ) const
    6274{
  • trunk/source/global/management/include/G4LPhysicsFreeVector.icc

    r1058 r1193  
    2525//
    2626//
    27 // $Id: G4LPhysicsFreeVector.icc,v 1.13 2008/09/22 08:26:33 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-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 $
    2929//
    3030//
     
    5252   if(binNumber == 0)
    5353     { edgeMin = binValue; }
    54    else if( numberOfBin - 1 == binNumber)
     54   else if( numberOfNodes - 1 == binNumber)
    5555     { edgeMax = binValue; }
    5656}
     
    8484{
    8585   G4int n1 = 0;
    86    G4int n2 = numberOfBin/2;
    87    G4int n3 = numberOfBin - 1;
     86   G4int n2 = numberOfNodes/2;
     87   G4int n3 = numberOfNodes - 1;
    8888   while (n1 != n3 - 1)
    8989   {
  • trunk/source/global/management/include/G4PhysicsFreeVector.hh

    r1058 r1193  
    2525//
    2626//
    27 // $Id: G4PhysicsFreeVector.hh,v 1.12 2008/09/22 14:49:57 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-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 $
    2929//
    3030//
     
    108108
    109109  size_t lowerBound = 0;
    110   size_t upperBound = numberOfBin-1;
     110  size_t upperBound = numberOfNodes-2;
    111111
    112112  while (lowerBound <= upperBound)
  • trunk/source/global/management/include/G4PhysicsOrderedFreeVector.hh

    r1058 r1193  
    2525//
    2626//
    27 // $Id: G4PhysicsOrderedFreeVector.hh,v 1.11 2008/09/22 14:49:57 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-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 $
    2929//
    3030////////////////////////////////////////////////////////////////////////
     
    8686        ////////////
    8787
    88         void InsertValues(G4double energy, G4double value);
     88        void InsertValues(G4double energy, G4double value);
    8989
    9090        G4double GetLowEdgeEnergy(size_t binNumber) const;
  • trunk/source/global/management/include/G4PhysicsOrderedFreeVector.icc

    r1058 r1193  
    2525//
    2626//
    27 // $Id: G4PhysicsOrderedFreeVector.icc,v 1.8 2006/06/29 19:02:31 gunter Exp $
    28 // GEANT4 tag $Name: geant4-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 $
    2929//
    3030////////////////////////////////////////////////////////////////////////
     
    7373void G4PhysicsOrderedFreeVector::DumpValues()
    7474{
    75    for (size_t i = 0; i < numberOfBin; i++)
     75   for (size_t i = 0; i < numberOfNodes; i++)
    7676   {
    7777      G4cout << binVector[i] << "\t" << dataVector[i] << G4endl;
     
    8383{
    8484   G4int n1 = 0;
    85    G4int n2 = numberOfBin/2;
    86    G4int n3 = numberOfBin - 1;
     85   G4int n2 = numberOfNodes/2;
     86   G4int n3 = numberOfNodes - 1;
    8787   while (n1 != n3 - 1)
    8888   {
  • trunk/source/global/management/include/G4PhysicsVector.hh

    r1058 r1193  
    2525//
    2626//
    27 // $Id: G4PhysicsVector.hh,v 1.18 2008/09/22 08:26:33 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-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 $
    2929//
    3030//
     
    5151//    09 Mar. 2001, H.Kurashige : Added G4PhysicsVectorType & Store/Retrieve()
    5252//    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
    5355//
    5456//---------------------------------------------------------------
     
    8486         // destructor
    8587
    86     inline G4double GetValue(G4double theEnergy, G4bool& isOutRange);
     88    inline G4double Value(G4double theEnergy);
    8789         // Get the cross-section/energy-loss value corresponding to the
    8890         // given energy. An appropriate interpolation is used to calculate
    8991         // 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.
    9296
    9397    G4int operator==(const G4PhysicsVector &right) const ;
    9498    G4int operator!=(const G4PhysicsVector &right) const ;
     99
    95100    inline G4double operator[](const size_t binNumber) const ;
    96101         // Returns simply the value in the bin specified by 'binNumber'
    97          // of the dataVector. The boundary check will be Done. If you
    98          // don't want this check, use the operator ().
     102         // of the dataVector. The boundary check will not be done.
     103
    99104    inline G4double operator()(const size_t binNumber) const ;
    100105         // 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);
    105109         // 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 vector
     110         // Take note that the 'index' starts from '0'.
     111         // To fill the vector, you have beforehand to construct a vector
    108112         // 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
    112128    virtual G4double GetLowEdgeEnergy(size_t binNumber) const;
     129         // Obsolete method
    113130         // Get the energy value at the low edge of the specified bin.
    114131         // Take note that the 'binNumber' starts from '0'.
    115          // This value is defined when a physics vector is constructed
    116          // by a constructor of a derived class. Use this function
    117          // 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
    118135    inline size_t GetVectorLength() const;
    119136         // 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
    120155    inline G4bool IsFilledVectorExist() const;
    121156         // Is non-empty physics vector already exist?
     
    124159         // Put a comment to the G4PhysicsVector. This may help to check
    125160         // whether your are accessing to the one you want.
     161
    126162    inline const G4String& GetComment() const;
    127163         // Retrieve the comment of the G4PhysicsVector.
     
    154190    G4PhysicsVectorType type;   // The type of PhysicsVector (enumerator)
    155191
    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;
    159196
    160197    G4double lastEnergy;        // Cache the last input value
     
    163200
    164201    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
    166204
    167205  private:
     206
     207    G4bool SplinePossible();
    168208
    169209    inline G4double LinearInterpolation();
     
    174214    inline void Interpolation();
    175215
    176     void FillSecondDerivatives();
    177       // Initialise second derivatives for spline
    178 
    179     G4double*  secDerivative;
    180 
    181216    G4String   comment;
    182217    G4bool     useSpline;
  • trunk/source/global/management/include/G4PhysicsVector.icc

    r1058 r1193  
    2525//
    2626//
    27 // $Id: G4PhysicsVector.icc,v 1.17 2008/09/22 08:26:33 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-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 $
    2929//
    3030//
     
    5959//---------------------------------------------------------------
    6060
     61inline
     62 G4double G4PhysicsVector::Energy(const size_t binNumber) const
     63{
     64  return binVector[binNumber];
     65}
     66
     67//---------------------------------------------------------------
     68
    6169inline
    6270 size_t G4PhysicsVector::GetVectorLength() const
    6371{
    64   return numberOfBin;
     72  return numberOfNodes;
    6573}
    6674
     
    94102  // numberOfBin-1.
    95103
    96   if( !secDerivative ) { FillSecondDerivatives(); }
     104  if(0 == secDerivative.size() ) { FillSecondDerivatives(); }
    97105
    98106  // 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;
    103113   
    104114  // Final evaluation of cubic spline polynomial for return   
     
    107117
    108118  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;
    111121
    112122  return res;
     
    117127inline
    118128 G4double G4PhysicsVector::GetValue(G4double theEnergy, G4bool&)
     129{
     130  return Value(theEnergy);
     131}
     132
     133//---------------------------------------------------------------
     134
     135inline
     136 G4double G4PhysicsVector::Value(G4double theEnergy)
    119137{
    120138  // Use cache for speed up - check if the value 'theEnergy' is same as the
     
    134152     lastValue  = dataVector[0];
    135153
    136   } else if(theEnergy < lastEnergy && theEnergy >= binVector[lastBin-1]) {
    137      lastBin--;
    138      lastEnergy = theEnergy;
    139      Interpolation();
    140 
    141154  } else if( theEnergy >= edgeMax ){
    142      lastBin = numberOfBin-1;
     155     lastBin = numberOfNodes-2;
    143156     lastEnergy = edgeMax;
    144      lastValue  = dataVector[lastBin];
     157     lastValue  = dataVector[numberOfNodes-1];
    145158
    146159  } else {
    147160     lastBin = FindBinLocation(theEnergy);
     161     if(lastBin >= numberOfNodes-1) {lastBin = numberOfNodes-2;}
    148162     lastEnergy = theEnergy;
    149163     Interpolation();
     
    167181{
    168182  dataVector[binNumber] = theValue;
    169 
    170   // Fill the bin which is hidden to user with theValue. This is to
    171   // handle correctly when Energy=theEmax in getValue.
    172 
    173   if(binNumber==numberOfBin-1) { dataVector[binNumber+1] = theValue; }
    174183}
    175184
     
    181190  G4bool status=false;
    182191
    183   if(numberOfBin > 0) { status=true; }
     192  if(numberOfNodes > 0) { status=true; }
    184193  return status;
    185194}
  • trunk/source/global/management/include/G4String.hh

    r1058 r1193  
    2525//
    2626//
    27 // $Id: G4String.hh,v 1.10 2008/12/08 14:16:05 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-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 $
    2929//
    3030//
     
    7878  inline G4int operator!() const;
    7979
    80   inline G4bool operator==(G4String) const;
     80  inline G4bool operator==(const G4String&) const;
    8181  inline G4bool operator==(const char*) const;
    82   inline G4bool operator!=(G4String) const;
     82  inline G4bool operator!=(const G4String&) const;
    8383  inline G4bool operator!=(const char*) const;
    8484
     
    138138  inline G4bool operator!=(const char*) const;
    139139
    140   //inline G4String operator () (unsigned int, unsigned int);
    141140  inline operator const char*() const;
    142141  inline G4SubString operator()(str_size, str_size);
     
    160159  inline G4int last(char) const;
    161160
    162   inline G4bool contains(std::string) const;
     161  inline G4bool contains(const std::string&) const;
    163162  inline G4bool contains(char) const;
    164163
     
    184183  inline unsigned int hash( caseCompare cmp = exact ) const;
    185184  inline unsigned int stlhash() const;
    186 
    187   // useful for supplying hash functions to template hash collection ctors
    188   //
    189   static inline unsigned int hash(const G4String&);
    190 
    191185};
    192186
  • trunk/source/global/management/include/G4String.icc

    r1058 r1193  
    2525//
    2626//
    27 // $Id: G4String.icc,v 1.11 2008/12/08 14:16:05 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-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 $
    2929//
    3030//
     
    9191G4int G4SubString::operator!() const
    9292{
    93   return extent==0 ? 1 : 0;
    94 }
    95 
    96 G4bool G4SubString::operator==(G4String s) const
    97 {
    98   return mystring->substr(mystart,extent) == s ? true:false;
     93  return (extent==0) ? 1 : 0;
     94}
     95
     96G4bool G4SubString::operator==(const G4String& s) const
     97{
     98  return (mystring->find(s,mystart,extent) != std::string::npos);
    9999}
    100100
    101101G4bool G4SubString::operator==(const char* s) const
    102102{
    103   return mystring->substr(mystart,extent) == std::string(s) ? true:false;
    104 }
    105 
    106 G4bool G4SubString::operator!=(G4String s) const
    107 {
    108   return mystring->substr(mystart,extent) != std::string(s) ? true:false;
     103  return (mystring->find(s,mystart,extent) != std::string::npos);
     104}
     105
     106G4bool G4SubString::operator!=(const G4String& s) const
     107{
     108  return (mystring->find(s,mystart,extent) == std::string::npos);
    109109}
    110110
    111111G4bool G4SubString::operator!=(const char* s) const
    112112{
    113   return mystring->substr(mystart,extent) != std::string(s) ? true:false;
     113  return (mystring->find(s,mystart,extent) == std::string::npos);
    114114}
    115115
     
    126126G4bool G4SubString::isNull() const
    127127{
    128   return extent==0 ? true : false;
     128  return (extent==0);
    129129}
    130130
     
    219219G4bool G4String::operator==(const G4String &s) const
    220220{
    221   const std_string *a=this;
    222   const std_string *b=&s;
    223   return *a==*b;
     221  return (std_string::compare(s) == 0);
    224222}
    225223
    226224G4bool G4String::operator==(const char* s) const
    227225{
    228   const std_string *a=this;
    229   const std_string b=s;
    230   return *a==b;
     226  return (std_string::compare(s) == 0);
    231227}
    232228
    233229G4bool G4String::operator!=(const G4String &s) const
    234230{
    235   const std_string *a=this;
    236   const std_string *b=&s;
    237   return *a!=*b;
     231  return !(*this == s);
    238232}
    239233
    240234G4bool G4String::operator!=(const char* s) const
    241235{
    242   const std_string *a=this;
    243   const std_string b=s;
    244   return *a!=b;
     236  return !(*this == s);
    245237}
    246238
     
    268260G4int G4String::compareTo(const char* s, caseCompare mode)
    269261{
    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);
    274264}
    275265
    276266G4int G4String::compareTo(const G4String& s, caseCompare mode)
    277267{
    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());
    282270}
    283271
     
    298286{
    299287  char tmp[1024];
    300   if ( skipWhite ) {
     288  if ( skipWhite )
     289  {
    301290    s >> std::ws;
    302291    s.getline(tmp,1024);
    303292    *this=tmp;
    304293  }
    305   else {
     294  else
     295  {
    306296    s.getline(tmp,1024);   
    307297    *this=tmp;
     
    325315G4String& G4String::remove(str_size n)
    326316{
    327   if(n<size())
    328     { erase(n,size()-n); }
     317  if(n<size()) { erase(n,size()-n); }
    329318  return *this;
    330319}
     
    346335}
    347336
    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; }
     337G4bool G4String::contains(const std::string& s) const
     338{
     339  return (std_string::find(s) != std_string::npos);
    355340}
    356341
    357342G4bool G4String::contains(char c) const
    358343{
    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);
    364345}
    365346
     
    474455  return str_size(h);
    475456}
    476 
    477 unsigned int G4String::hash(const G4String& s)
    478 {
    479   return s.hash();
    480 }
  • trunk/source/global/management/include/G4Version.hh

    r1058 r1193  
    2525//
    2626//
    27 // $Id: G4Version.hh,v 1.17 2008/12/02 09:15:58 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-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 $
    2929//
    3030// Version information
     
    4747
    4848#ifndef G4VERSION_NUMBER
    49 #define G4VERSION_NUMBER  920
     49#define G4VERSION_NUMBER  930
    5050#endif
    5151
    5252#ifndef G4VERSION_TAG
    53 #define G4VERSION_TAG "$Name: geant4-09-02-ref-02 $"
     53#define G4VERSION_TAG "$Name: global-V09-02-10 $"
    5454#endif
    5555
     
    5858#include "G4String.hh"
    5959
    60 static const G4String G4Version = "$Name: geant4-09-02-ref-02 $";
    61 static const G4String G4Date    = "(19-December-2008)";
     60static const G4String G4Version = "$Name: global-V09-02-10 $";
     61static const G4String G4Date    = "(5-June-2009)";
    6262
    6363#endif
  • trunk/source/global/management/src/G4ErrorPropagatorData.cc

    r1058 r1193  
    2525//
    2626//
    27 // $Id: G4ErrorPropagatorData.cc,v 1.3 2007/06/05 13:04:30 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-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 $
    2929//
    3030//
     
    4343//-------------------------------------------------------------------
    4444
    45 G4ErrorPropagatorData::G4ErrorPropagatorData(){}
    46 G4ErrorPropagatorData::~G4ErrorPropagatorData(){}
     45G4ErrorPropagatorData::~G4ErrorPropagatorData()
     46{
     47}
     48
     49G4ErrorPropagatorData::G4ErrorPropagatorData()
     50{
     51  theStage = G4ErrorStage_Inflation;
     52}
    4753
    4854G4ErrorPropagatorData* G4ErrorPropagatorData::GetErrorPropagatorData()
  • trunk/source/global/management/src/G4LPhysicsFreeVector.cc

    r1058 r1193  
    2525//
    2626//
    27 // $Id: G4LPhysicsFreeVector.cc,v 1.22 2008/09/22 08:26:33 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-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 $
    2929//
    3030//
     
    3535//
    3636// F.W. Jones, TRIUMF, 04-JUN-96
     37//
     38// Modified:
     39//    19 Jun. 2009, V.Ivanchenko : removed hidden bin
    3740//
    3841// --------------------------------------------------------------------
    3942
    4043#include "G4LPhysicsFreeVector.hh"
    41 
    4244#include "G4ios.hh"
    4345
     
    6163   edgeMin = binmin;
    6264   edgeMax = binmax;
    63    numberOfBin = nbin;
    64    binVector.reserve(nbin+1);
    65    dataVector.reserve(nbin+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++)
    6769   {
    6870     binVector.push_back(0.0);
     
    105107void G4LPhysicsFreeVector::DumpValues()
    106108{
    107    for (size_t i = 0; i < numberOfBin; i++)
     109   for (size_t i = 0; i < numberOfNodes; i++)
    108110   {
    109111      G4cout << binVector[i] << "   " << dataVector[i]/millibarn << G4endl;
  • trunk/source/global/management/src/G4PhysicsFreeVector.cc

    r1058 r1193  
    2525//
    2626//
    27 // $Id: G4PhysicsFreeVector.cc,v 1.12 2008/09/22 14:49:57 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-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 $
    2929//
    3030//
     
    4141//    26 Sep. 1996, K.Amako : Constructor with only 'bin size' added
    4242//    11 Nov. 2000, H.Kurashige : use STL vector for dataVector and binVector
     43//    19 Jun. 2009, V.Ivanchenko : removed hidden bin
    4344//
    4445//--------------------------------------------------------------------
     
    5758{
    5859  type = T_G4PhysicsFreeVector;
    59   numberOfBin = theNbin;
     60  numberOfNodes = theNbin;
    6061
    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);
    6564
    66   for (size_t i=0; i<=numberOfBin; i++)
     65  for (size_t i=0; i<numberOfNodes; i++)
    6766  {
    6867     binVector.push_back(0.0);
     
    7574{
    7675  type = T_G4PhysicsFreeVector;
    77   numberOfBin = theBinVector.size();
     76  numberOfNodes = theBinVector.size();
    7877
    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);
    8380
    84   for (size_t i=0; i<numberOfBin; i++)
     81  for (size_t i=0; i<numberOfNodes; i++)
    8582  {
    8683     binVector.push_back(theBinVector[i]);
     
    8885  }
    8986
    90   // Put values to extra hidden bin. For 'binVector', the 'edgeMin' of the
    91   // extra hidden bin is assumed to have the following value. For binary
    92   // search, this value is completely arbitrary if it is greater than
    93   // 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 the
    98   // extra hidden bin is assumed to have the same as the one at 'numberBin-1'.
    99   dataVector.push_back( theDataVector[numberOfBin-1] );
    100 
    10187  edgeMin = binVector[0];
    102   edgeMax = binVector[numberOfBin-1];
     88  edgeMax = binVector[numberOfNodes-1];
    10389
    10490
     
    11399  dataVector[theBinNumber] = theDataValue;
    114100
    115   if( theBinNumber == numberOfBin-1 )
     101  if( theBinNumber == numberOfNodes-1 )
    116102  {
    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];
    129104  }
    130105
  • trunk/source/global/management/src/G4PhysicsLinearVector.cc

    r1058 r1193  
    2525//
    2626//
    27 // $Id: G4PhysicsLinearVector.cc,v 1.14 2008/09/22 14:49:57 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-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 $
    2929//
    3030//
     
    3535//
    3636//  15 Feb 1996 - K.Amako : 1st version
     37//  19 Jun 2009 - V.Ivanchenko : removed hidden bin
    3738//
    3839//--------------------------------------------------------------------
     
    5152  type = T_G4PhysicsLinearVector;
    5253
    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);     
    5757
    58   numberOfBin = theNbin;
    59 
    60   for (size_t i=0; i<=numberOfBin; i++)
     58  for (size_t i=0; i<numberOfNodes; i++)
    6159  {
    6260     binVector.push_back(0.0);
     
    7371  type = T_G4PhysicsLinearVector;
    7472
    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);     
    7976
    80   numberOfBin = theNbin;
    81 
    82   for (size_t i=0; i<numberOfBin+1; i++)
     77  for (size_t i=0; i<numberOfNodes; i++)
    8378  {
    8479    binVector.push_back( theEmin + i*dBin );
     
    8782
    8883  edgeMin = binVector[0];
    89   edgeMax = binVector[numberOfBin-1];
     84  edgeMax = binVector[numberOfNodes-1];
    9085
    9186
  • trunk/source/global/management/src/G4PhysicsLnVector.cc

    r1058 r1193  
    2525//
    2626//
    27 // $Id: G4PhysicsLnVector.cc,v 1.17 2008/09/22 14:49:57 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-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 $
    2929//
    3030//
     
    3535//
    3636//  27 Apr 1999 - M.G.Pia: Created from G4PhysicsLogVector
     37//  19 Jun 2009 - V.Ivanchenko : removed hidden bin
    3738//
    3839// --------------------------------------------------------------
     
    5152  type = T_G4PhysicsLnVector;
    5253
    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);     
    5757
    58   numberOfBin = theNbin;
    59 
    60   for (size_t i=0; i<=numberOfBin; i++)
     58  for (size_t i=0; i<numberOfNodes; i++)
    6159  {
    6260     binVector.push_back(0.0);
     
    7371  type = T_G4PhysicsLnVector;
    7472
    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);     
    7976
    80   numberOfBin = theNbin;
    81 
    82   for (size_t i=0; i<numberOfBin+1; i++)
     77  for (size_t i=0; i<numberOfNodes; i++)
    8378  {
    84     binVector.push_back(std::exp(std::log(theEmin)+i*dBin));
     79    binVector.push_back(std::exp((baseBin+i)*dBin));
    8580    dataVector.push_back(0.0);
    8681  }
    8782
    8883  edgeMin = binVector[0];
    89   edgeMax = binVector[numberOfBin-1];
     84  edgeMax = binVector[numberOfNodes-1];
    9085}
    9186
  • trunk/source/global/management/src/G4PhysicsLogVector.cc

    r1058 r1193  
    2525//
    2626//
    27 // $Id: G4PhysicsLogVector.cc,v 1.21 2008/09/22 08:26:33 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-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 $
    2929//
    3030//
     
    4242//    9  Mar. 2001, H.Kurashige : added PhysicsVector type and Retrieve
    4343//    05 Sep. 2008, V.Ivanchenko : added protections for zero-length vector
     44//    19 Jun. 2009, V.Ivanchenko : removed hidden bin
    4445//
    4546// --------------------------------------------------------------
     
    5859  type = T_G4PhysicsLogVector;
    5960
    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);     
    6464
    65   numberOfBin = theNbin;
    66 
    67   for (size_t i=0; i<=numberOfBin; i++)
     65  for (size_t i=0; i<numberOfNodes; i++)
    6866  {
    6967     binVector.push_back(0.0);
     
    7977  type = T_G4PhysicsLogVector;
    8078
    81   // Add extra one bin (hidden to user) to handle correctly when
    82   // 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.);
    8583
    86   numberOfBin = theNbin;
    87 
    88   for (size_t i=0; i<numberOfBin+1; i++)
     84  for (size_t i=0; i<numberOfNodes; i++)
    8985  {
    90     binVector.push_back(std::pow(10., std::log10(theEmin)+i*dBin));
     86    binVector.push_back(std::exp(g4log10*(baseBin+i)*dBin));
    9187    dataVector.push_back(0.0);
    9288  }
     89
    9390  edgeMin = binVector[0];
    94   edgeMax = binVector[numberOfBin-1];
     91  edgeMax = binVector[numberOfNodes-1];
    9592
    9693
  • trunk/source/global/management/src/G4PhysicsOrderedFreeVector.cc

    r1058 r1193  
    2525//
    2626//
    27 // $Id: G4PhysicsOrderedFreeVector.cc,v 1.12 2008/09/22 14:49:57 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-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 $
    2929//
    3030////////////////////////////////////////////////////////////////////////
     
    4343//              2000-11-11 by H.Kurashige
    4444//              > use STL vector for dataVector and binVector
     45//              19 Jun. 2009-06-19 by V.Ivanchenko
     46//              > removed hidden bin
     47//
    4548// mail:        gum@triumf.ca
    46 //
    4749//
    4850////////////////////////////////////////////////////////////////////////
     
    6567        type = T_G4PhysicsOrderedFreeVector;
    6668
    67         dataVector.reserve(VectorLength+1);
    68         binVector.reserve(VectorLength+1);
    69         numberOfBin = VectorLength;
     69        dataVector.reserve(VectorLength);
     70        binVector.reserve(VectorLength);
     71        numberOfNodes = VectorLength;
    7072
    7173        for (size_t i = 0 ; i < VectorLength ; i++)
     
    7476                dataVector.push_back(Values[i]);
    7577        }
    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];
    8080}
    8181
     
    101101        binVector.push_back(energy);
    102102        dataVector.push_back(value);
    103         numberOfBin++;
     103        numberOfNodes++;
    104104        edgeMin = binVector.front();
    105105        edgeMax = binVector.back();
     
    133133{
    134134   G4int n1 = 0;
    135    G4int n2 = numberOfBin/2;
    136    G4int n3 = numberOfBin - 1;
     135   G4int n2 = numberOfNodes/2;
     136   G4int n3 = numberOfNodes - 1;
    137137   while (n1 != n3 - 1) {
    138138      if (aValue > dataVector[n2])
  • trunk/source/global/management/src/G4PhysicsVector.cc

    r1058 r1193  
    2525//
    2626//
    27 // $Id: G4PhysicsVector.cc,v 1.27 2008/10/16 12:14:36 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-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 $
    2929//
    3030//
     
    4343//    09 Mar. 2001, H.Kurashige : added G4PhysicsVector type
    4444//    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
    4551// --------------------------------------------------------------
    4652
     
    5258G4PhysicsVector::G4PhysicsVector(G4bool spline)
    5359 : 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)
    5762{}
    5863
     
    6065
    6166G4PhysicsVector::~G4PhysicsVector()
    62 {
    63   DeleteData();
    64 }
     67{}
    6568
    6669// --------------------------------------------------------------
     
    7881  if (type != right.type)  { return *this; }
    7982
    80   DeleteData();
     83  //DeleteData();
    8184  CopyData(right);
    8285
     
    102105void G4PhysicsVector::DeleteData()
    103106{
    104   delete [] secDerivative;
    105   secDerivative = 0;
     107  secDerivative.clear();
    106108}
    107109
     
    113115  edgeMin = vec.edgeMin;
    114116  edgeMax = vec.edgeMax;
    115   numberOfBin = vec.numberOfBin;
     117  numberOfNodes = vec.numberOfNodes;
    116118  lastEnergy = vec.lastEnergy;
    117119  lastValue = vec.lastValue;
     
    121123  useSpline = vec.useSpline;
    122124  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;
    135126}
    136127
     
    157148  fOut.write((char*)(&edgeMin), sizeof edgeMin);
    158149  fOut.write((char*)(&edgeMax), sizeof edgeMax);
    159   fOut.write((char*)(&numberOfBin), sizeof numberOfBin);
     150  fOut.write((char*)(&numberOfNodes), sizeof numberOfNodes);
    160151
    161152  // contents
     
    164155
    165156  G4double* value = new G4double[2*size];
    166   for(size_t i = 0; i < size; i++)
     157  for(size_t i = 0; i < size; ++i)
    167158  {
    168159    value[2*i]  =  binVector[i];
     
    180171{
    181172  // clear properties;
    182   lastEnergy=0.;
     173  lastEnergy=-1.0;
    183174  lastValue =0.;
    184175  lastBin   =0;
    185176  dataVector.clear();
    186177  binVector.clear();
     178  secDerivative.clear();
    187179  comment = "";
    188180
     
    191183  {
    192184    // binning
    193     fIn >> edgeMin >> edgeMax >> numberOfBin;
     185    fIn >> edgeMin >> edgeMax >> numberOfNodes;
    194186    if (fIn.fail())  { return false; }
    195187    // contents
     
    218210  fIn.read((char*)(&edgeMin), sizeof edgeMin);
    219211  fIn.read((char*)(&edgeMax), sizeof edgeMax);
    220   fIn.read((char*)(&numberOfBin), sizeof numberOfBin );
     212  fIn.read((char*)(&numberOfNodes), sizeof numberOfNodes );
    221213 
    222214  // contents
     
    234226  binVector.reserve(size);
    235227  dataVector.reserve(size);
    236   for(size_t i = 0; i < size; i++)
     228  for(size_t i = 0; i < size; ++i)
    237229  {
    238230    binVector.push_back(value[2*i]);
     
    245237// --------------------------------------------------------------
    246238
     239void
     240G4PhysicsVector::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
     261void
     262G4PhysicsVector::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
    247325void 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
    248329
    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();
    258333    return;
    259334  }
    260335
    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
     397void
     398G4PhysicsVector::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] =
    264414      3.0*((dataVector[i+1]-dataVector[i])/(binVector[i+1]-binVector[i]) -
    265415           (dataVector[i]-dataVector[i-1])/(binVector[i]-binVector[i-1]))
     
    269419  secDerivative[0] = secDerivative[1];
    270420}
     421
     422// --------------------------------------------------------------
     423
     424G4bool 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}
    271441   
    272442// --------------------------------------------------------------
     
    276446  // binning
    277447  out << std::setprecision(12) << pv.edgeMin;
    278   out <<" " << pv.edgeMax <<" "  << pv.numberOfBin << G4endl;
     448  out <<" " << pv.edgeMax <<" "  << pv.numberOfNodes << G4endl;
    279449
    280450  // contents
Note: See TracChangeset for help on using the changeset viewer.