ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLIFunction1D.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLIFunction1D.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France
33 //
34 #define INCLXX_IN_GEANT4_MODE 1
35 
36 #include "globals.hh"
37 
45 #include <algorithm>
46 #include <cmath>
47 #include <cstdlib>
48 #include "G4INCLIFunction1D.hh"
49 #include "G4INCLLogger.hh"
51 
52 namespace G4INCL {
53 
55  2.*95.0/288.0,
56  317.0/240.0,
57  23.0/30.0,
58  793.0/720.0,
59  157.0/160.0,
60  157.0/160.0,
61  793.0/720.0,
62  23.0/30.0,
63  317.0/240.0,
64  };
65 
67  G4double xi = std::max(x0, xMin);
68  G4double xa = std::min(x1, xMax);
69  G4double sign;
70 
71  if(x1 <= x0) {
72  sign = -1.0;
73  std::swap(xi, xa);
74  } else
75  sign = 1.0;
76 
77  const G4double interval = xa - xi;
78 
79  G4int nIntervals;
80  if(step<0.) {
81  nIntervals = 45;
82  } else {
83  nIntervals = G4int(interval/step);
84 
85  // Round up nIntervals to the closest multiple of 9
86  G4int remainder = nIntervals % 9;
87  if (remainder != 0)
88  nIntervals += 9 - remainder;
89 
90  nIntervals = std::max(nIntervals, 9);
91  }
92 
93  const G4double dx = interval/nIntervals;
94  G4double result = (operator()(xi) + operator()(xa)) * integrationCoefficients[0]/2;
95  for(G4int j = 1; j<nIntervals; ++j) {
96  const G4double x = xi + interval*G4double(j)/G4double(nIntervals);
97  const unsigned index = j%9;
98  result += operator()(x) * integrationCoefficients[index];
99  }
100 
101  return result*dx*sign;
102 
103  }
104 
106  class Primitive : public IFunction1D {
107  public:
108  Primitive(IFunction1D const * const f) :
109  IFunction1D(f->getXMinimum(), f->getXMaximum()),
110  theFunction(f)
111  {}
112 
113  G4double operator()(const G4double x) const {
114  return theFunction->integrate(xMin,x);
115  }
116  private:
117  IFunction1D const * const theFunction;
118  } *thePrimitive = new Primitive(this);
119 
120  return thePrimitive;
121  }
122 
124  class InverseCDF : public IFunction1D {
125  public:
126  InverseCDF(IFunction1D const * const f, ManipulatorFunc fw) :
127  IFunction1D(f->getXMinimum(), f->getXMaximum()),
128  theFunction(f),
129  normalisation(1./theFunction->integrate(xMin,xMax)),
130  fWrap(fw)
131  {}
132 
133  G4double operator()(const G4double x) const {
134  if(fWrap)
135  return fWrap(std::min(1., normalisation * theFunction->integrate(xMin,x)));
136  else
137  return std::min(1., normalisation * theFunction->integrate(xMin,x));
138  }
139  private:
140  IFunction1D const * const theFunction;
141  const G4double normalisation;
142  ManipulatorFunc fWrap;
143  } *theInverseCDF = new InverseCDF(this, fWrap);
144 
145  InterpolationTable *theTable = new InvFInterpolationTable(*theInverseCDF, nNodes);
146  delete theInverseCDF;
147  return theTable;
148  }
149 }
150