Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PenelopeSamplingData.cc
Go to the documentation of this file.
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 // $Id: G4PenelopeSamplingData.cc 66241 2012-12-13 18:34:42Z gunter $
27 //
28 // Author: Luciano Pandola
29 //
30 // History:
31 // --------
32 // 09 Dec 2009 L Pandola First implementation
33 //
35 
36 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
38  np(nPoints)
39 {
40  //create vectors
41  x = new G4DataVector();
42  pac = new G4DataVector();
43  a = new G4DataVector();
44  b = new G4DataVector();
45  ITTL = new std::vector<size_t>;
46  ITTU = new std::vector<size_t>;
47 }
48 
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
51 {
52  if (x) delete x;
53  if (pac) delete pac;
54  if (a) delete a;
55  if (b) delete b;
56  if (ITTL) delete ITTL;
57  if (ITTU) delete ITTU;
58 }
59 
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.......oooOO0OOooo...
62 {
63  size_t points = x->size();
64 
65  //check everything is all right
66  if (pac->size() != points || a->size() != points ||
67  b->size() != points || ITTL->size() != points ||
68  ITTU->size() != points)
69  {
71  ed << "Data vectors look to have different dimensions !" << G4endl;
72  G4Exception("G4PenelopeSamplingData::GetNumberOfStoredPoints()","em2040",
73  FatalException,ed);
74  }
75  return points;
76 }
77 
78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
80 {
81  if (x) delete x;
82  if (pac) delete pac;
83  if (a) delete a;
84  if (b) delete b;
85  if (ITTL) delete ITTL;
86  if (ITTU) delete ITTU;
87  //create vectors
88  x = new G4DataVector();
89  pac = new G4DataVector();
90  a = new G4DataVector();
91  b = new G4DataVector();
92  ITTL = new std::vector<size_t>;
93  ITTU = new std::vector<size_t>;
94 }
95 
96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
98  size_t ITTL0,size_t ITTU0)
99 {
100  x->push_back(x0);
101  pac->push_back(pac0);
102  a->push_back(a0);
103  b->push_back(b0);
104  ITTL->push_back(ITTL0);
105  ITTU->push_back(ITTU0);
106 
107  //check how many points we do have now
108  size_t nOfPoints = GetNumberOfStoredPoints();
109 
110  if (nOfPoints > ((size_t)np))
111  {
112  G4cout << "G4PenelopeSamplingData::AddPoint() " << G4endl;
113  G4cout << "WARNING: Up to now there are " << nOfPoints << " points in the table" << G4endl;
114  G4cout << "while the anticipated (declared) number is " << np << G4endl;
115  }
116  return;
117 }
118 
119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
121 {
122 
123  G4cout << "*************************************************************************" << G4endl;
124  G4cout << GetNumberOfStoredPoints() << " points" << G4endl;
125  G4cout << "*************************************************************************" << G4endl;
126  for (size_t i=0;i<GetNumberOfStoredPoints();i++)
127  {
128  G4cout << i << " " << (*x)[i] << " " << (*pac)[i] << " " << (*a)[i] << " " <<
129  (*b)[i] << " " << (*ITTL)[i] << " " << (*ITTU)[i] << G4endl;
130  }
131  G4cout << "*************************************************************************" << G4endl;
132 }
133 
134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
136 {
137  if (index < x->size())
138  return (*x)[index];
139  else
140  return 0;
141 }
142 
143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
145 {
146  if (index < pac->size())
147  return (*pac)[index];
148  else
149  return 0;
150 }
151 
152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
154 {
155  if (index < a->size())
156  return (*a)[index];
157  else
158  return 0;
159 }
160 
161 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
163 {
164  if (index < b->size())
165  return (*b)[index];
166  else
167  return 0;
168 }
169 
170 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
172 {
173  //One passes here a random number in (0,1).
174  //Notice: it possible that is between (0,b) with b<1
175  size_t points = GetNumberOfStoredPoints();
176 
177  size_t itn = (size_t) (maxRand*(points-1));
178  size_t i = (*ITTL)[itn];
179  size_t j = (*ITTU)[itn];
180 
181  while ((j-i) > 1)
182  {
183  size_t k = (i+j)/2;
184  if (maxRand > (*pac)[k])
185  i = k;
186  else
187  j = k;
188  }
189 
190  //Sampling from the rational inverse cumulative distribution
191  G4double result = 0;
192 
193  G4double rr = maxRand - (*pac)[i];
194  if (rr > 1e-16)
195  {
196  G4double d = (*pac)[i+1]-(*pac)[i];
197  result = (*x)[i]+
198  ((1.0+(*a)[i]+(*b)[i])*d*rr/
199  (d*d+((*a)[i]*d+(*b)[i]*rr)*rr))*((*x)[i+1]-(*x)[i]);
200  }
201  else
202  result = (*x)[i];
203 
204  return result;
205 }
G4double G4ParticleHPJENDLHEData::G4double result
const G4double a0
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetB(size_t index)
G4double GetPAC(size_t index)
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double SampleValue(G4double rndm)
G4double GetA(size_t index)
G4double GetX(size_t index)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
void AddPoint(G4double x0, G4double pac0, G4double a0, G4double b0, size_t ITTL0, size_t ITTU0)
G4PenelopeSamplingData(G4int npoints=150)