| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // Author: Luciano Pandola 27 // Author: Luciano Pandola
28 // 28 //
29 // History: 29 // History:
30 // -------- 30 // --------
31 // 09 Dec 2009 L Pandola First implementa 31 // 09 Dec 2009 L Pandola First implementation
32 // 32 //
33 #include "G4PenelopeSamplingData.hh" 33 #include "G4PenelopeSamplingData.hh"
34 34
35 //....oooOO0OOooo........oooOO0OOooo........oo 35 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
36 G4PenelopeSamplingData::G4PenelopeSamplingData 36 G4PenelopeSamplingData::G4PenelopeSamplingData(G4int nPoints) :
37 fNP(nPoints) << 37 np(nPoints)
38 { 38 {
39 //create vectors 39 //create vectors
40 fX = new G4DataVector(); << 40 x = new G4DataVector();
41 fPAC = new G4DataVector(); << 41 pac = new G4DataVector();
42 fA = new G4DataVector(); << 42 a = new G4DataVector();
43 fB = new G4DataVector(); << 43 b = new G4DataVector();
44 fITTL = new std::vector<size_t>; << 44 ITTL = new std::vector<size_t>;
45 fITTU = new std::vector<size_t>; << 45 ITTU = new std::vector<size_t>;
46 } 46 }
47 47
48 //....oooOO0OOooo........oooOO0OOooo........oo 48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
49 G4PenelopeSamplingData::~G4PenelopeSamplingDat 49 G4PenelopeSamplingData::~G4PenelopeSamplingData()
50 { 50 {
51 if (fX) delete fX; << 51 if (x) delete x;
52 if (fPAC) delete fPAC; << 52 if (pac) delete pac;
53 if (fA) delete fA; << 53 if (a) delete a;
54 if (fB) delete fB; << 54 if (b) delete b;
55 if (fITTL) delete fITTL; << 55 if (ITTL) delete ITTL;
56 if (fITTU) delete fITTU; << 56 if (ITTU) delete ITTU;
57 } 57 }
58 58
59 //....oooOO0OOooo........oooOO0OOooo........oo 59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.......oooOO0OOooo...
60 <<
61 size_t G4PenelopeSamplingData::GetNumberOfStor 60 size_t G4PenelopeSamplingData::GetNumberOfStoredPoints()
62 { 61 {
63 size_t points = fX->size(); << 62 size_t points = x->size();
64 63
65 //check everything is all right 64 //check everything is all right
66 if (fPAC->size() != points || fA->size() != << 65 if (pac->size() != points || a->size() != points ||
67 fB->size() != points || fITTL->size() != << 66 b->size() != points || ITTL->size() != points ||
68 fITTU->size() != points) << 67 ITTU->size() != points)
69 { 68 {
70 G4ExceptionDescription ed; 69 G4ExceptionDescription ed;
71 ed << "Data vectors look to have differe 70 ed << "Data vectors look to have different dimensions !" << G4endl;
72 G4Exception("G4PenelopeSamplingData::Get 71 G4Exception("G4PenelopeSamplingData::GetNumberOfStoredPoints()","em2040",
73 FatalException,ed); 72 FatalException,ed);
74 } 73 }
75 return points; 74 return points;
76 } 75 }
77 76
78 //....oooOO0OOooo........oooOO0OOooo........oo 77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
79 void G4PenelopeSamplingData::Clear() 78 void G4PenelopeSamplingData::Clear()
80 { 79 {
81 if (fX) delete fX; << 80 if (x) delete x;
82 if (fPAC) delete fPAC; << 81 if (pac) delete pac;
83 if (fA) delete fA; << 82 if (a) delete a;
84 if (fB) delete fB; << 83 if (b) delete b;
85 if (fITTL) delete fITTL; << 84 if (ITTL) delete ITTL;
86 if (fITTU) delete fITTU; << 85 if (ITTU) delete ITTU;
87 //create vectors 86 //create vectors
88 fX = new G4DataVector(); << 87 x = new G4DataVector();
89 fPAC = new G4DataVector(); << 88 pac = new G4DataVector();
90 fA = new G4DataVector(); << 89 a = new G4DataVector();
91 fB = new G4DataVector(); << 90 b = new G4DataVector();
92 fITTL = new std::vector<size_t>; << 91 ITTL = new std::vector<size_t>;
93 fITTU = new std::vector<size_t>; << 92 ITTU = new std::vector<size_t>;
94 } 93 }
95 94
96 //....oooOO0OOooo........oooOO0OOooo........oo 95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
97 void G4PenelopeSamplingData::AddPoint(G4double 96 void G4PenelopeSamplingData::AddPoint(G4double x0,G4double pac0,G4double a0,G4double b0,
98 size_t ITTL0,size_t ITTU0) 97 size_t ITTL0,size_t ITTU0)
99 { 98 {
100 fX->push_back(x0); << 99 x->push_back(x0);
101 fPAC->push_back(pac0); << 100 pac->push_back(pac0);
102 fA->push_back(a0); << 101 a->push_back(a0);
103 fB->push_back(b0); << 102 b->push_back(b0);
104 fITTL->push_back(ITTL0); << 103 ITTL->push_back(ITTL0);
105 fITTU->push_back(ITTU0); << 104 ITTU->push_back(ITTU0);
106 105
107 //check how many points we do have now 106 //check how many points we do have now
108 size_t nOfPoints = GetNumberOfStoredPoints() 107 size_t nOfPoints = GetNumberOfStoredPoints();
109 108
110 if (nOfPoints > ((size_t)fNP)) << 109 if (nOfPoints > ((size_t)np))
111 { 110 {
112 G4cout << "G4PenelopeSamplingData::AddPo 111 G4cout << "G4PenelopeSamplingData::AddPoint() " << G4endl;
113 G4cout << "WARNING: Up to now there are 112 G4cout << "WARNING: Up to now there are " << nOfPoints << " points in the table" << G4endl;
114 G4cout << "while the anticipated (declar << 113 G4cout << "while the anticipated (declared) number is " << np << G4endl;
115 } 114 }
116 return; 115 return;
117 } 116 }
118 117
119 //....oooOO0OOooo........oooOO0OOooo........oo 118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
120 void G4PenelopeSamplingData::DumpTable() 119 void G4PenelopeSamplingData::DumpTable()
121 { 120 {
122 121
123 G4cout << "********************************* 122 G4cout << "*************************************************************************" << G4endl;
124 G4cout << GetNumberOfStoredPoints() << " poi 123 G4cout << GetNumberOfStoredPoints() << " points" << G4endl;
125 G4cout << "********************************* 124 G4cout << "*************************************************************************" << G4endl;
126 for (size_t i=0;i<GetNumberOfStoredPoints(); 125 for (size_t i=0;i<GetNumberOfStoredPoints();i++)
127 { 126 {
128 G4cout << i << " " << (*fX)[i] << " " << << 127 G4cout << i << " " << (*x)[i] << " " << (*pac)[i] << " " << (*a)[i] << " " <<
129 (*fB)[i] << " " << (*fITTL)[i] << " " << (*f << 128 (*b)[i] << " " << (*ITTL)[i] << " " << (*ITTU)[i] << G4endl;
130 } 129 }
131 G4cout << "********************************* 130 G4cout << "*************************************************************************" << G4endl;
132 } 131 }
133 132
134 //....oooOO0OOooo........oooOO0OOooo........oo 133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
135 G4double G4PenelopeSamplingData::GetX(size_t i 134 G4double G4PenelopeSamplingData::GetX(size_t index)
136 { 135 {
137 if (index < fX->size()) << 136 if (index < x->size())
138 return (*fX)[index]; << 137 return (*x)[index];
139 else 138 else
140 return 0; 139 return 0;
141 } 140 }
142 141
143 //....oooOO0OOooo........oooOO0OOooo........oo 142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
144 G4double G4PenelopeSamplingData::GetPAC(size_t 143 G4double G4PenelopeSamplingData::GetPAC(size_t index)
145 { 144 {
146 if (index < fPAC->size()) << 145 if (index < pac->size())
147 return (*fPAC)[index]; << 146 return (*pac)[index];
148 else 147 else
149 return 0; 148 return 0;
150 } 149 }
151 150
152 //....oooOO0OOooo........oooOO0OOooo........oo 151 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
153 G4double G4PenelopeSamplingData::GetA(size_t i 152 G4double G4PenelopeSamplingData::GetA(size_t index)
154 { 153 {
155 if (index < fA->size()) << 154 if (index < a->size())
156 return (*fA)[index]; << 155 return (*a)[index];
157 else 156 else
158 return 0; 157 return 0;
159 } 158 }
160 159
161 //....oooOO0OOooo........oooOO0OOooo........oo 160 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
162 G4double G4PenelopeSamplingData::GetB(size_t i 161 G4double G4PenelopeSamplingData::GetB(size_t index)
163 { 162 {
164 if (index < fB->size()) << 163 if (index < b->size())
165 return (*fB)[index]; << 164 return (*b)[index];
166 else 165 else
167 return 0; 166 return 0;
168 } 167 }
169 168
170 //....oooOO0OOooo........oooOO0OOooo........oo 169 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
171 G4double G4PenelopeSamplingData::SampleValue(G 170 G4double G4PenelopeSamplingData::SampleValue(G4double maxRand)
172 { 171 {
173 //One passes here a random number in (0,1). 172 //One passes here a random number in (0,1).
174 //Notice: it possible that is between (0,b) 173 //Notice: it possible that is between (0,b) with b<1
175 size_t points = GetNumberOfStoredPoints(); 174 size_t points = GetNumberOfStoredPoints();
176 175
177 size_t itn = (size_t) (maxRand*(points-1)); 176 size_t itn = (size_t) (maxRand*(points-1));
178 size_t i = (*fITTL)[itn]; << 177 size_t i = (*ITTL)[itn];
179 size_t j = (*fITTU)[itn]; << 178 size_t j = (*ITTU)[itn];
180 179
181 while ((j-i) > 1) 180 while ((j-i) > 1)
182 { 181 {
183 size_t k = (i+j)/2; 182 size_t k = (i+j)/2;
184 if (maxRand > (*fPAC)[k]) << 183 if (maxRand > (*pac)[k])
185 i = k; 184 i = k;
186 else 185 else
187 j = k; 186 j = k;
188 } 187 }
189 188
190 //Sampling from the rational inverse cumulat 189 //Sampling from the rational inverse cumulative distribution
191 G4double result = 0; 190 G4double result = 0;
192 191
193 G4double rr = maxRand - (*fPAC)[i]; << 192 G4double rr = maxRand - (*pac)[i];
194 if (rr > 1e-16) 193 if (rr > 1e-16)
195 { 194 {
196 G4double d = (*fPAC)[i+1]-(*fPAC)[i]; << 195 G4double d = (*pac)[i+1]-(*pac)[i];
197 result = (*fX)[i]+ << 196 result = (*x)[i]+
198 ((1.0+(*fA)[i]+(*fB)[i])*d*rr/ << 197 ((1.0+(*a)[i]+(*b)[i])*d*rr/
199 (d*d+((*fA)[i]*d+(*fB)[i]*rr)*rr))*((*fX)[i << 198 (d*d+((*a)[i]*d+(*b)[i]*rr)*rr))*((*x)[i+1]-(*x)[i]);
200 } 199 }
201 else 200 else
202 result = (*fX)[i]; << 201 result = (*x)[i];
203 202
204 return result; 203 return result;
205 } 204 }
206 205