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

CVS update

Location:
trunk/source/global/management/src
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • 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.