G4INCLIFunction1D.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // INCL++ intra-nuclear cascade model
00027 // Pekka Kaitaniemi, CEA and Helsinki Institute of Physics
00028 // Davide Mancusi, CEA
00029 // Alain Boudard, CEA
00030 // Sylvie Leray, CEA
00031 // Joseph Cugnon, University of Liege
00032 //
00033 #define INCLXX_IN_GEANT4_MODE 1
00034 
00035 #include "globals.hh"
00036 
00044 #include <algorithm>
00045 #include <cmath>
00046 #include <cstdlib>
00047 #include "G4INCLIFunction1D.hh"
00048 #include "G4INCLLogger.hh"
00049 #include "G4INCLInverseInterpolationTable.hh"
00050 
00051 namespace G4INCL {
00052 
00053   const G4double IFunction1D::integrationCoefficients[] = {
00054       2.*95.0/288.0,
00055       317.0/240.0,
00056       23.0/30.0,
00057       793.0/720.0,
00058       157.0/160.0,
00059       157.0/160.0,
00060       793.0/720.0,
00061       23.0/30.0,
00062       317.0/240.0,
00063   };
00064 
00065   G4double IFunction1D::integrate(const G4double x0, const G4double x1, const G4double step) const {
00066     G4double xi = std::max(x0, xMin);
00067     G4double xa = std::min(x1, xMax);
00068     G4double sign;
00069 
00070     if(x1 <= x0) {
00071       sign = -1.0;
00072       std::swap(xi, xa);
00073     } else
00074       sign = 1.0;
00075 
00076     const G4double interval = xa - xi;
00077 
00078     G4int nIntervals;
00079     if(step<0.) {
00080       nIntervals = 45;
00081     } else {
00082       nIntervals = G4int(interval/step);
00083 
00084       // Round up nIntervals to the closest multiple of 9
00085       G4int remainder = nIntervals % 9;
00086       if (remainder != 0)
00087         nIntervals += 9 - remainder;
00088 
00089       nIntervals = std::max(nIntervals, 9);
00090     }
00091 
00092     const G4double dx = interval/nIntervals;
00093     G4double result = (operator()(xi) + operator()(xa)) * integrationCoefficients[0]/2;
00094     for(G4int j = 1; j<nIntervals; ++j) {
00095       const G4double x = xi + interval*G4double(j)/G4double(nIntervals);
00096       const unsigned index = j%9;
00097       result += operator()(x) * integrationCoefficients[index];
00098     }
00099 
00100     return result*dx*sign;
00101 
00102   }
00103 
00104   IFunction1D *IFunction1D::primitive() const {
00105     class Primitive : public IFunction1D {
00106       public:
00107         Primitive(IFunction1D const * const f) :
00108           IFunction1D(f->getXMinimum(), f->getXMaximum()),
00109           theFunction(f)
00110       {}
00111 
00112         G4double operator()(const G4double x) const {
00113           return theFunction->integrate(xMin,x);
00114         }
00115       private:
00116         IFunction1D const * const theFunction;
00117     } *thePrimitive = new Primitive(this);
00118 
00119     return thePrimitive;
00120   }
00121 
00122   InverseInterpolationTable *IFunction1D::inverseCDFTable(const G4int nNodes) const {
00123     class InverseCDF : public IFunction1D {
00124       public:
00125         InverseCDF(IFunction1D const * const f) :
00126           IFunction1D(f->getXMinimum(), f->getXMaximum()),
00127           theFunction(f),
00128           normalisation(1./theFunction->integrate(xMin,xMax))
00129       {}
00130 
00131         G4double operator()(const G4double x) const {
00132           return std::min(1., normalisation * theFunction->integrate(xMin,x));
00133         }
00134       private:
00135         IFunction1D const * const theFunction;
00136         const G4double normalisation;
00137     } *theInverseCDF = new InverseCDF(this);
00138 
00139     InverseInterpolationTable *theTable = new InverseInterpolationTable(*theInverseCDF, nNodes);
00140     delete theInverseCDF;
00141     return theTable;
00142   }
00143 
00144 }
00145 

Generated on Mon May 27 17:48:35 2013 for Geant4 by  doxygen 1.4.7