| 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 // $Id$
26 // 27 //
27 // Author: Luciano Pandola 28 // Author: Luciano Pandola
28 // 29 //
29 // History: 30 // History:
30 // -------- 31 // --------
31 // 09 Dec 2009 L Pandola First implementa 32 // 09 Dec 2009 L Pandola First implementation
32 // 33 //
33 #include "G4PenelopeSamplingData.hh" 34 #include "G4PenelopeSamplingData.hh"
34 35
35 //....oooOO0OOooo........oooOO0OOooo........oo 36 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
36 G4PenelopeSamplingData::G4PenelopeSamplingData 37 G4PenelopeSamplingData::G4PenelopeSamplingData(G4int nPoints) :
37 fNP(nPoints) << 38 np(nPoints)
38 { 39 {
39 //create vectors 40 //create vectors
40 fX = new G4DataVector(); << 41 x = new G4DataVector();
41 fPAC = new G4DataVector(); << 42 pac = new G4DataVector();
42 fA = new G4DataVector(); << 43 a = new G4DataVector();
43 fB = new G4DataVector(); << 44 b = new G4DataVector();
44 fITTL = new std::vector<size_t>; << 45 ITTL = new std::vector<size_t>;
45 fITTU = new std::vector<size_t>; << 46 ITTU = new std::vector<size_t>;
46 } 47 }
47 48
48 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
49 G4PenelopeSamplingData::~G4PenelopeSamplingDat 50 G4PenelopeSamplingData::~G4PenelopeSamplingData()
50 { 51 {
51 if (fX) delete fX; << 52 if (x) delete x;
52 if (fPAC) delete fPAC; << 53 if (pac) delete pac;
53 if (fA) delete fA; << 54 if (a) delete a;
54 if (fB) delete fB; << 55 if (b) delete b;
55 if (fITTL) delete fITTL; << 56 if (ITTL) delete ITTL;
56 if (fITTU) delete fITTU; << 57 if (ITTU) delete ITTU;
57 } 58 }
58 59
59 //....oooOO0OOooo........oooOO0OOooo........oo 60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.......oooOO0OOooo...
60 <<
61 size_t G4PenelopeSamplingData::GetNumberOfStor 61 size_t G4PenelopeSamplingData::GetNumberOfStoredPoints()
62 { 62 {
63 size_t points = fX->size(); << 63 size_t points = x->size();
64 64
65 //check everything is all right 65 //check everything is all right
66 if (fPAC->size() != points || fA->size() != << 66 if (pac->size() != points || a->size() != points ||
67 fB->size() != points || fITTL->size() != << 67 b->size() != points || ITTL->size() != points ||
68 fITTU->size() != points) << 68 ITTU->size() != points)
69 { 69 {
70 G4ExceptionDescription ed; 70 G4ExceptionDescription ed;
71 ed << "Data vectors look to have differe 71 ed << "Data vectors look to have different dimensions !" << G4endl;
72 G4Exception("G4PenelopeSamplingData::Get 72 G4Exception("G4PenelopeSamplingData::GetNumberOfStoredPoints()","em2040",
73 FatalException,ed); 73 FatalException,ed);
74 } 74 }
75 return points; 75 return points;
76 } 76 }
77 77
78 //....oooOO0OOooo........oooOO0OOooo........oo 78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
79 void G4PenelopeSamplingData::Clear() 79 void G4PenelopeSamplingData::Clear()
80 { 80 {
81 if (fX) delete fX; << 81 if (x) delete x;
82 if (fPAC) delete fPAC; << 82 if (pac) delete pac;
83 if (fA) delete fA; << 83 if (a) delete a;
84 if (fB) delete fB; << 84 if (b) delete b;
85 if (fITTL) delete fITTL; << 85 if (ITTL) delete ITTL;
86 if (fITTU) delete fITTU; << 86 if (ITTU) delete ITTU;
87 //create vectors 87 //create vectors
88 fX = new G4DataVector(); << 88 x = new G4DataVector();
89 fPAC = new G4DataVector(); << 89 pac = new G4DataVector();
90 fA = new G4DataVector(); << 90 a = new G4DataVector();
91 fB = new G4DataVector(); << 91 b = new G4DataVector();
92 fITTL = new std::vector<size_t>; << 92 ITTL = new std::vector<size_t>;
93 fITTU = new std::vector<size_t>; << 93 ITTU = new std::vector<size_t>;
94 } 94 }
95 95
96 //....oooOO0OOooo........oooOO0OOooo........oo 96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
97 void G4PenelopeSamplingData::AddPoint(G4double 97 void G4PenelopeSamplingData::AddPoint(G4double x0,G4double pac0,G4double a0,G4double b0,
98 size_t ITTL0,size_t ITTU0) 98 size_t ITTL0,size_t ITTU0)
99 { 99 {
100 fX->push_back(x0); << 100 x->push_back(x0);
101 fPAC->push_back(pac0); << 101 pac->push_back(pac0);
102 fA->push_back(a0); << 102 a->push_back(a0);
103 fB->push_back(b0); << 103 b->push_back(b0);
104 fITTL->push_back(ITTL0); << 104 ITTL->push_back(ITTL0);
105 fITTU->push_back(ITTU0); << 105 ITTU->push_back(ITTU0);
106 106
107 //check how many points we do have now 107 //check how many points we do have now
108 size_t nOfPoints = GetNumberOfStoredPoints() 108 size_t nOfPoints = GetNumberOfStoredPoints();
109 109
110 if (nOfPoints > ((size_t)fNP)) << 110 if (nOfPoints > ((size_t)np))
111 { 111 {
112 G4cout << "G4PenelopeSamplingData::AddPo 112 G4cout << "G4PenelopeSamplingData::AddPoint() " << G4endl;
113 G4cout << "WARNING: Up to now there are 113 G4cout << "WARNING: Up to now there are " << nOfPoints << " points in the table" << G4endl;
114 G4cout << "while the anticipated (declar << 114 G4cout << "while the anticipated (declared) number is " << np << G4endl;
115 } 115 }
116 return; 116 return;
117 } 117 }
118 118
119 //....oooOO0OOooo........oooOO0OOooo........oo 119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
120 void G4PenelopeSamplingData::DumpTable() 120 void G4PenelopeSamplingData::DumpTable()
121 { 121 {
122 122
123 G4cout << "********************************* 123 G4cout << "*************************************************************************" << G4endl;
124 G4cout << GetNumberOfStoredPoints() << " poi 124 G4cout << GetNumberOfStoredPoints() << " points" << G4endl;
125 G4cout << "********************************* 125 G4cout << "*************************************************************************" << G4endl;
126 for (size_t i=0;i<GetNumberOfStoredPoints(); 126 for (size_t i=0;i<GetNumberOfStoredPoints();i++)
127 { 127 {
128 G4cout << i << " " << (*fX)[i] << " " << << 128 G4cout << i << " " << (*x)[i] << " " << (*pac)[i] << " " << (*a)[i] << " " <<
129 (*fB)[i] << " " << (*fITTL)[i] << " " << (*f << 129 (*b)[i] << " " << (*ITTL)[i] << " " << (*ITTU)[i] << G4endl;
130 } 130 }
131 G4cout << "********************************* 131 G4cout << "*************************************************************************" << G4endl;
132 } 132 }
133 133
134 //....oooOO0OOooo........oooOO0OOooo........oo 134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
135 G4double G4PenelopeSamplingData::GetX(size_t i 135 G4double G4PenelopeSamplingData::GetX(size_t index)
136 { 136 {
137 if (index < fX->size()) << 137 if (index < x->size())
138 return (*fX)[index]; << 138 return (*x)[index];
139 else 139 else
140 return 0; 140 return 0;
141 } 141 }
142 142
143 //....oooOO0OOooo........oooOO0OOooo........oo 143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
144 G4double G4PenelopeSamplingData::GetPAC(size_t 144 G4double G4PenelopeSamplingData::GetPAC(size_t index)
145 { 145 {
146 if (index < fPAC->size()) << 146 if (index < pac->size())
147 return (*fPAC)[index]; << 147 return (*pac)[index];
148 else 148 else
149 return 0; 149 return 0;
150 } 150 }
151 151
152 //....oooOO0OOooo........oooOO0OOooo........oo 152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
153 G4double G4PenelopeSamplingData::GetA(size_t i 153 G4double G4PenelopeSamplingData::GetA(size_t index)
154 { 154 {
155 if (index < fA->size()) << 155 if (index < a->size())
156 return (*fA)[index]; << 156 return (*a)[index];
157 else 157 else
158 return 0; 158 return 0;
159 } 159 }
160 160
161 //....oooOO0OOooo........oooOO0OOooo........oo 161 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
162 G4double G4PenelopeSamplingData::GetB(size_t i 162 G4double G4PenelopeSamplingData::GetB(size_t index)
163 { 163 {
164 if (index < fB->size()) << 164 if (index < b->size())
165 return (*fB)[index]; << 165 return (*b)[index];
166 else 166 else
167 return 0; 167 return 0;
168 } 168 }
169 169
170 //....oooOO0OOooo........oooOO0OOooo........oo 170 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
171 G4double G4PenelopeSamplingData::SampleValue(G 171 G4double G4PenelopeSamplingData::SampleValue(G4double maxRand)
172 { 172 {
173 //One passes here a random number in (0,1). 173 //One passes here a random number in (0,1).
174 //Notice: it possible that is between (0,b) 174 //Notice: it possible that is between (0,b) with b<1
175 size_t points = GetNumberOfStoredPoints(); 175 size_t points = GetNumberOfStoredPoints();
176 176
177 size_t itn = (size_t) (maxRand*(points-1)); 177 size_t itn = (size_t) (maxRand*(points-1));
178 size_t i = (*fITTL)[itn]; << 178 size_t i = (*ITTL)[itn];
179 size_t j = (*fITTU)[itn]; << 179 size_t j = (*ITTU)[itn];
180 180
181 while ((j-i) > 1) 181 while ((j-i) > 1)
182 { 182 {
183 size_t k = (i+j)/2; 183 size_t k = (i+j)/2;
184 if (maxRand > (*fPAC)[k]) << 184 if (maxRand > (*pac)[k])
185 i = k; 185 i = k;
186 else 186 else
187 j = k; 187 j = k;
188 } 188 }
189 189
190 //Sampling from the rational inverse cumulat 190 //Sampling from the rational inverse cumulative distribution
191 G4double result = 0; 191 G4double result = 0;
192 192
193 G4double rr = maxRand - (*fPAC)[i]; << 193 G4double rr = maxRand - (*pac)[i];
194 if (rr > 1e-16) 194 if (rr > 1e-16)
195 { 195 {
196 G4double d = (*fPAC)[i+1]-(*fPAC)[i]; << 196 G4double d = (*pac)[i+1]-(*pac)[i];
197 result = (*fX)[i]+ << 197 result = (*x)[i]+
198 ((1.0+(*fA)[i]+(*fB)[i])*d*rr/ << 198 ((1.0+(*a)[i]+(*b)[i])*d*rr/
199 (d*d+((*fA)[i]*d+(*fB)[i]*rr)*rr))*((*fX)[i << 199 (d*d+((*a)[i]*d+(*b)[i]*rr)*rr))*((*x)[i+1]-(*x)[i]);
200 } 200 }
201 else 201 else
202 result = (*fX)[i]; << 202 result = (*x)[i];
203 203
204 return result; 204 return result;
205 } 205 }
206 206