| 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 // ------------------------------------------- 27 // -------------------------------------------------------------------
28 // 28 //
29 // GEANT4 Class 29 // GEANT4 Class
30 // File name: G4PAIModelData.cc 30 // File name: G4PAIModelData.cc
31 // 31 //
32 // Author: V. Ivanchenko based on V.Gri 32 // Author: V. Ivanchenko based on V.Grichine code of G4PAIModel
33 // 33 //
34 // Creation date: 16.08.2013 34 // Creation date: 16.08.2013
35 // 35 //
36 // Modifications: 36 // Modifications:
37 // 37 //
38 38
39 #include "G4PAIModelData.hh" 39 #include "G4PAIModelData.hh"
40 #include "G4PAIModel.hh" 40 #include "G4PAIModel.hh"
41 41
42 #include "G4SystemOfUnits.hh" 42 #include "G4SystemOfUnits.hh"
43 #include "G4PhysicalConstants.hh" 43 #include "G4PhysicalConstants.hh"
44 44
45 #include "G4PhysicsLogVector.hh" 45 #include "G4PhysicsLogVector.hh"
46 #include "G4PhysicsFreeVector.hh" 46 #include "G4PhysicsFreeVector.hh"
47 #include "G4PhysicsTable.hh" 47 #include "G4PhysicsTable.hh"
48 #include "G4MaterialCutsCouple.hh" 48 #include "G4MaterialCutsCouple.hh"
49 #include "G4SandiaTable.hh" 49 #include "G4SandiaTable.hh"
50 #include "Randomize.hh" 50 #include "Randomize.hh"
51 #include "G4Poisson.hh" 51 #include "G4Poisson.hh"
52 52
53 ////////////////////////////////////////////// 53 ////////////////////////////////////////////////////////////////////////
54 54
55 using namespace std; 55 using namespace std;
56 56
57 G4PAIModelData::G4PAIModelData(G4double tmin, 57 G4PAIModelData::G4PAIModelData(G4double tmin, G4double tmax, G4int ver)
58 { 58 {
59 const G4int nPerDecade = 10; 59 const G4int nPerDecade = 10;
60 const G4double lowestTkin = 50*keV; 60 const G4double lowestTkin = 50*keV;
61 const G4double highestTkin = 10*TeV; 61 const G4double highestTkin = 10*TeV;
62 62
63 fPAIySection.SetVerbose(ver); 63 fPAIySection.SetVerbose(ver);
64 64
65 fLowestKineticEnergy = std::max(tmin, lowes 65 fLowestKineticEnergy = std::max(tmin, lowestTkin);
66 fHighestKineticEnergy = tmax; 66 fHighestKineticEnergy = tmax;
67 if(tmax < 10*fLowestKineticEnergy) { 67 if(tmax < 10*fLowestKineticEnergy) {
68 fHighestKineticEnergy = 10*fLowestKineticE 68 fHighestKineticEnergy = 10*fLowestKineticEnergy;
69 } else if(tmax > highestTkin) { 69 } else if(tmax > highestTkin) {
70 fHighestKineticEnergy = std::max(highestTk 70 fHighestKineticEnergy = std::max(highestTkin, 10*fLowestKineticEnergy);
71 } 71 }
72 fTotBin = (G4int)(nPerDecade* 72 fTotBin = (G4int)(nPerDecade*
73 std::log10(fHighestKineticEnergy/fLowe 73 std::log10(fHighestKineticEnergy/fLowestKineticEnergy));
74 74
75 fParticleEnergyVector = new G4PhysicsLogVect 75 fParticleEnergyVector = new G4PhysicsLogVector(fLowestKineticEnergy,
76 fHighestKineticEnergy, 76 fHighestKineticEnergy,
77 fTotBin); 77 fTotBin);
78 if(0 < ver) { 78 if(0 < ver) {
79 G4cout << "### G4PAIModelData: Nbins= " << 79 G4cout << "### G4PAIModelData: Nbins= " << fTotBin
80 << " Tlowest(keV)= " << lowestTkin/keV 80 << " Tlowest(keV)= " << lowestTkin/keV
81 << " Tmin(keV)= " << fLowestKineticEnergy 81 << " Tmin(keV)= " << fLowestKineticEnergy/keV
82 << " Tmax(GeV)= " << fHighestKineticEnerg 82 << " Tmax(GeV)= " << fHighestKineticEnergy/GeV
83 << G4endl; 83 << G4endl;
84 } 84 }
85 } 85 }
86 86
87 ////////////////////////////////////////////// 87 ////////////////////////////////////////////////////////////////////////////
88 88
89 G4PAIModelData::~G4PAIModelData() 89 G4PAIModelData::~G4PAIModelData()
90 { 90 {
91 std::size_t n = fPAIxscBank.size(); << 91 size_t n = fPAIxscBank.size();
92 if(0 < n) { 92 if(0 < n) {
93 for(std::size_t i=0; i<n; ++i) { << 93 for(size_t i=0; i<n; ++i) {
94 if(fPAIxscBank[i]) { 94 if(fPAIxscBank[i]) {
95 fPAIxscBank[i]->clearAndDestroy(); 95 fPAIxscBank[i]->clearAndDestroy();
96 delete fPAIxscBank[i]; 96 delete fPAIxscBank[i];
97 } 97 }
98 if(fPAIdEdxBank[i]) { 98 if(fPAIdEdxBank[i]) {
99 fPAIdEdxBank[i]->clearAndDestroy(); 99 fPAIdEdxBank[i]->clearAndDestroy();
100 delete fPAIdEdxBank[i]; 100 delete fPAIdEdxBank[i];
101 } 101 }
102 delete fdEdxTable[i]; 102 delete fdEdxTable[i];
103 } 103 }
104 } 104 }
105 delete fParticleEnergyVector; 105 delete fParticleEnergyVector;
106 } 106 }
107 107
108 ////////////////////////////////////////////// 108 ///////////////////////////////////////////////////////////////////////////////
109 109
110 void G4PAIModelData::Initialise(const G4Materi 110 void G4PAIModelData::Initialise(const G4MaterialCutsCouple* couple,
111 G4PAIModel* mo 111 G4PAIModel* model)
112 { 112 {
113 const G4Material* mat = couple->GetMaterial( 113 const G4Material* mat = couple->GetMaterial();
114 fSandia.Initialize(const_cast<G4Material*>(m 114 fSandia.Initialize(const_cast<G4Material*>(mat));
115 115
116 auto PAItransferTable = new G4PhysicsTable(f << 116 G4PhysicsTable* PAItransferTable = new G4PhysicsTable(fTotBin+1);
117 auto PAIdEdxTable = new G4PhysicsTable(fTotB << 117 G4PhysicsTable* PAIdEdxTable = new G4PhysicsTable(fTotBin+1);
118 auto dEdxMeanVector = << 118 G4PhysicsLogVector* dEdxMeanVector =
119 new G4PhysicsLogVector(fLowestKineticEnerg 119 new G4PhysicsLogVector(fLowestKineticEnergy,
120 fHighestKineticEnergy, 120 fHighestKineticEnergy,
121 fTotBin); 121 fTotBin);
122 // low energy Sandia interval 122 // low energy Sandia interval
123 G4double Tmin = fSandia.GetSandiaMatTablePAI 123 G4double Tmin = fSandia.GetSandiaMatTablePAI(0,0);
124 124
125 // energy safety 125 // energy safety
126 static const G4double deltaLow = 100.*eV; 126 static const G4double deltaLow = 100.*eV;
127 127
128 for (G4int i = 0; i <= fTotBin; ++i) { 128 for (G4int i = 0; i <= fTotBin; ++i) {
129 129
130 G4double kinEnergy = fParticleEnergyVector 130 G4double kinEnergy = fParticleEnergyVector->Energy(i);
131 G4double Tmax = model->ComputeMaxEnergy(ki 131 G4double Tmax = model->ComputeMaxEnergy(kinEnergy);
132 G4double tau = kinEnergy/proton_mass_c2; 132 G4double tau = kinEnergy/proton_mass_c2;
133 G4double bg2 = tau*( tau + 2. ); 133 G4double bg2 = tau*( tau + 2. );
134 134
135 if (Tmax < Tmin + deltaLow ) { Tmax = Tmin 135 if (Tmax < Tmin + deltaLow ) { Tmax = Tmin + deltaLow; }
136 136
137 fPAIySection.Initialize(mat, Tmax, bg2, &f 137 fPAIySection.Initialize(mat, Tmax, bg2, &fSandia);
138 138
139 //G4cout << i << ". TransferMax(keV)= "<< 139 //G4cout << i << ". TransferMax(keV)= "<< Tmax/keV
140 // << " E(MeV)= " << kinEnergy/MeV << 140 // << " E(MeV)= " << kinEnergy/MeV << G4endl;
141 141
142 G4int n = fPAIySection.GetSplineSize(); 142 G4int n = fPAIySection.GetSplineSize();
143 G4int kmin = 0; 143 G4int kmin = 0;
144 for(G4int k = 0; k < n; ++k) { 144 for(G4int k = 0; k < n; ++k) {
145 if(fPAIySection.GetIntegralPAIySection(k 145 if(fPAIySection.GetIntegralPAIySection(k+1) <= 0.0) {
146 kmin = k; 146 kmin = k;
147 } else { 147 } else {
148 break; 148 break;
149 } 149 }
150 } 150 }
151 n -= kmin; 151 n -= kmin;
152 152
153 auto transferVector = new G4PhysicsFreeVec << 153 G4PhysicsFreeVector* transferVector = new G4PhysicsFreeVector(n);
154 auto dEdxVector = new G4PhysicsFreeVector( << 154 G4PhysicsFreeVector* dEdxVector = new G4PhysicsFreeVector(n);
155 155
156 //G4double tr0 = 0.0; 156 //G4double tr0 = 0.0;
157 G4double tr = 0.0; 157 G4double tr = 0.0;
158 for(G4int k = kmin; k < n; ++k) 158 for(G4int k = kmin; k < n; ++k)
159 { 159 {
160 G4double t = fPAIySection.GetSplineEner 160 G4double t = fPAIySection.GetSplineEnergy(k+1);
161 tr = fPAIySection.GetIntegralPAIySection 161 tr = fPAIySection.GetIntegralPAIySection(k+1);
162 //if(tr >= tr0) { tr0 = tr; } 162 //if(tr >= tr0) { tr0 = tr; }
163 //else { G4cout << "G4PAIModelData::Init 163 //else { G4cout << "G4PAIModelData::Initialise Warning: Ekin(MeV)= "
164 // << t/MeV << " IntegralTransfer 164 // << t/MeV << " IntegralTransfer= " << tr
165 // << " < " << tr0 << G4endl; } 165 // << " < " << tr0 << G4endl; }
166 transferVector->PutValue(k, t, t*tr); 166 transferVector->PutValue(k, t, t*tr);
167 dEdxVector->PutValue(k, t, fPAIySection. 167 dEdxVector->PutValue(k, t, fPAIySection.GetIntegralPAIdEdx(k+1));
168 } 168 }
169 //G4cout << "TransferVector:" << G4endl; 169 //G4cout << "TransferVector:" << G4endl;
170 //G4cout << *transferVector << G4endl; 170 //G4cout << *transferVector << G4endl;
171 //G4cout << "DEDXVector:" << G4endl; 171 //G4cout << "DEDXVector:" << G4endl;
172 //G4cout << *dEdxVector << G4endl; 172 //G4cout << *dEdxVector << G4endl;
173 173
174 G4double ionloss = std::max(fPAIySection.G << 174 G4double ionloss = fPAIySection.GetMeanEnergyLoss();// total <dE/dx>
>> 175
>> 176 if(ionloss < 0.0) ionloss = 0.0;
>> 177
175 dEdxMeanVector->PutValue(i,ionloss); 178 dEdxMeanVector->PutValue(i,ionloss);
176 179
177 PAItransferTable->insertAt(i,transferVecto 180 PAItransferTable->insertAt(i,transferVector);
178 PAIdEdxTable->insertAt(i,dEdxVector); 181 PAIdEdxTable->insertAt(i,dEdxVector);
179 182
>> 183 //transferVector->SetSpline(true);
>> 184 //transferVector->FillSecondDerivatives();
>> 185 //dEdxVector->SetSpline(true);
>> 186 //dEdxVector->FillSecondDerivatives();
>> 187
180 } // end of Tkin loop` 188 } // end of Tkin loop`
181 fPAIxscBank.push_back(PAItransferTable); 189 fPAIxscBank.push_back(PAItransferTable);
182 fPAIdEdxBank.push_back(PAIdEdxTable); 190 fPAIdEdxBank.push_back(PAIdEdxTable);
183 //G4cout << "dEdxMeanVector: " << G4endl; 191 //G4cout << "dEdxMeanVector: " << G4endl;
184 //G4cout << *dEdxMeanVector << G4endl; 192 //G4cout << *dEdxMeanVector << G4endl;
>> 193 /*
>> 194 dEdxMeanVector->SetSpline(true);
>> 195 dEdxMeanVector->FillSecondDerivatives();
>> 196 */
185 fdEdxTable.push_back(dEdxMeanVector); 197 fdEdxTable.push_back(dEdxMeanVector);
186 } 198 }
187 199
188 ////////////////////////////////////////////// 200 //////////////////////////////////////////////////////////////////////////////
189 201
190 G4double G4PAIModelData::DEDXPerVolume(G4int c 202 G4double G4PAIModelData::DEDXPerVolume(G4int coupleIndex, G4double scaledTkin,
191 G4double cut) const 203 G4double cut) const
192 { 204 {
193 // VI: iPlace is the low edge index of the b 205 // VI: iPlace is the low edge index of the bin
194 // iPlace is in interval from 0 to (N-1) 206 // iPlace is in interval from 0 to (N-1)
195 std::size_t iPlace(0); << 207 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
196 G4double dEdx = fdEdxTable[coupleIndex]->Val << 208 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
197 std::size_t nPlace = fParticleEnergyVector-> <<
198 /* 209 /*
199 G4cout << "G4PAIModelData::DEDXPerVolume: co 210 G4cout << "G4PAIModelData::DEDXPerVolume: coupleIdx= " << coupleIndex
200 << " Tscaled= " << scaledTkin << " cut= " < 211 << " Tscaled= " << scaledTkin << " cut= " << cut
201 << " iPlace= " << iPlace << " nPlace= " << 212 << " iPlace= " << iPlace << " nPlace= " << nPlace << G4endl;
202 */ 213 */
203 G4bool one = true; 214 G4bool one = true;
204 if(scaledTkin >= fParticleEnergyVector->Ener 215 if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
205 else if(scaledTkin > fParticleEnergyVector-> 216 else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
206 one = false; 217 one = false;
207 } 218 }
208 219
209 // VI: apply interpolation of the vector 220 // VI: apply interpolation of the vector
>> 221 G4double dEdx = fdEdxTable[coupleIndex]->Value(scaledTkin);
210 G4double del = (*(fPAIdEdxBank[coupleIndex] 222 G4double del = (*(fPAIdEdxBank[coupleIndex]))(iPlace)->Value(cut);
211 //G4cout << "dEdx= " << dEdx << " del= " << 223 //G4cout << "dEdx= " << dEdx << " del= " << del << G4endl;
212 if(!one) { 224 if(!one) {
213 G4double del2 = (*(fPAIdEdxBank[coupleInde 225 G4double del2 = (*(fPAIdEdxBank[coupleIndex]))(iPlace+1)->Value(cut);
214 G4double E1 = fParticleEnergyVector->Energ 226 G4double E1 = fParticleEnergyVector->Energy(iPlace);
215 G4double E2 = fParticleEnergyVector->Energ 227 G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
216 G4double W = 1.0/(E2 - E1); 228 G4double W = 1.0/(E2 - E1);
217 G4double W1 = (E2 - scaledTkin)*W; 229 G4double W1 = (E2 - scaledTkin)*W;
218 G4double W2 = (scaledTkin - E1)*W; 230 G4double W2 = (scaledTkin - E1)*W;
219 del *= W1; 231 del *= W1;
220 del += W2*del2; 232 del += W2*del2;
221 } 233 }
222 dEdx -= del; 234 dEdx -= del;
223 //G4cout << "dEdx= " << dEdx << " del= " << 235 //G4cout << "dEdx= " << dEdx << " del= " << del << G4endl;
224 236
225 dEdx = std::max(dEdx, 0.); 237 dEdx = std::max(dEdx, 0.);
226 return dEdx; 238 return dEdx;
227 } 239 }
228 240
229 ////////////////////////////////////////////// 241 /////////////////////////////////////////////////////////////////////////
230 242
231 G4double 243 G4double
232 G4PAIModelData::CrossSectionPerVolume(G4int co 244 G4PAIModelData::CrossSectionPerVolume(G4int coupleIndex,
233 G4double scaledTkin, 245 G4double scaledTkin,
234 G4double tcut, G4double tmax) co 246 G4double tcut, G4double tmax) const
235 { 247 {
236 G4double cross, cross1, cross2; 248 G4double cross, cross1, cross2;
237 249
238 // iPlace is in interval from 0 to (N-1) 250 // iPlace is in interval from 0 to (N-1)
239 std::size_t iPlace = fParticleEnergyVector-> << 251 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
240 std::size_t nPlace = fParticleEnergyVector-> << 252 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
241 253
242 G4bool one = true; 254 G4bool one = true;
243 if(scaledTkin >= fParticleEnergyVector->Ener 255 if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
244 else if(scaledTkin > fParticleEnergyVector-> 256 else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
245 one = false; 257 one = false;
246 } 258 }
247 G4PhysicsTable* table = fPAIxscBank[coupleIn 259 G4PhysicsTable* table = fPAIxscBank[coupleIndex];
248 260
249 //G4cout<<"iPlace = "<<iPlace<<"; tmax = " 261 //G4cout<<"iPlace = "<<iPlace<<"; tmax = "
250 // <<tmax<<"; cutEnergy = "<<cutEnergy<<G4en 262 // <<tmax<<"; cutEnergy = "<<cutEnergy<<G4endl;
251 cross1 = (*table)(iPlace)->Value(tmax)/tmax; 263 cross1 = (*table)(iPlace)->Value(tmax)/tmax;
252 //G4cout<<"cross1 = "<<cross1<<G4endl; 264 //G4cout<<"cross1 = "<<cross1<<G4endl;
253 cross2 = (*table)(iPlace)->Value(tcut)/tcut; 265 cross2 = (*table)(iPlace)->Value(tcut)/tcut;
254 //G4cout<<"cross2 = "<<cross2<<G4endl; 266 //G4cout<<"cross2 = "<<cross2<<G4endl;
255 cross = (cross2-cross1); 267 cross = (cross2-cross1);
256 //G4cout<<"cross = "<<cross<<G4endl; 268 //G4cout<<"cross = "<<cross<<G4endl;
257 if(!one) { 269 if(!one) {
258 cross2 = (*table)(iPlace+1)->Value(tcut)/t 270 cross2 = (*table)(iPlace+1)->Value(tcut)/tcut
259 - (*table)(iPlace+1)->Value(tmax)/tmax; 271 - (*table)(iPlace+1)->Value(tmax)/tmax;
260 272
261 G4double E1 = fParticleEnergyVector->Energ 273 G4double E1 = fParticleEnergyVector->Energy(iPlace);
262 G4double E2 = fParticleEnergyVector->Energ 274 G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
263 G4double W = 1.0/(E2 - E1); 275 G4double W = 1.0/(E2 - E1);
264 G4double W1 = (E2 - scaledTkin)*W; 276 G4double W1 = (E2 - scaledTkin)*W;
265 G4double W2 = (scaledTkin - E1)*W; 277 G4double W2 = (scaledTkin - E1)*W;
266 cross *= W1; 278 cross *= W1;
267 cross += W2*cross2; 279 cross += W2*cross2;
268 } 280 }
269 281
270 cross = std::max(cross, 0.0); 282 cross = std::max(cross, 0.0);
271 return cross; 283 return cross;
272 } 284 }
273 285
274 ////////////////////////////////////////////// 286 ///////////////////////////////////////////////////////////////////////
275 287
276 G4double G4PAIModelData::SampleAlongStepTransf 288 G4double G4PAIModelData::SampleAlongStepTransfer(G4int coupleIndex,
277 289 G4double kinEnergy,
278 G4double scaledTkin, 290 G4double scaledTkin,
279 G4double tmax, 291 G4double tmax,
280 G4double stepFactor) const 292 G4double stepFactor) const
281 { 293 {
282 //G4cout << "=== G4PAIModelData::SampleAlong 294 //G4cout << "=== G4PAIModelData::SampleAlongStepTransfer" << G4endl;
283 G4double loss = 0.0; 295 G4double loss = 0.0;
284 296
285 std::size_t iPlace = fParticleEnergyVector-> << 297 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
286 std::size_t nPlace = fParticleEnergyVector-> << 298 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
287 299
288 G4bool one = true; 300 G4bool one = true;
289 if(scaledTkin >= fParticleEnergyVector->Ener 301 if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
290 else if(scaledTkin > fParticleEnergyVector-> 302 else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
291 one = false; 303 one = false;
292 } 304 }
293 305
294 G4double meanNumber = 0.0; 306 G4double meanNumber = 0.0;
295 G4double meanN11 = 0.0; 307 G4double meanN11 = 0.0;
296 G4double meanN12 = 0.0; 308 G4double meanN12 = 0.0;
297 G4double meanN21 = 0.0; 309 G4double meanN21 = 0.0;
298 G4double meanN22 = 0.0; 310 G4double meanN22 = 0.0;
299 311
300 G4PhysicsVector* v1 = (*(fPAIxscBank[coupleI 312 G4PhysicsVector* v1 = (*(fPAIxscBank[coupleIndex]))(iPlace);
301 G4PhysicsVector* v2 = nullptr; << 313 G4PhysicsVector* v2 = 0;
302 314
303 G4double e1 = v1->Energy(0); 315 G4double e1 = v1->Energy(0);
304 G4double e2 = std::min(tmax, v1->GetMaxEnerg 316 G4double e2 = std::min(tmax, v1->GetMaxEnergy());
305 317
306 if(e2 >= e1) { 318 if(e2 >= e1) {
307 meanN11 = (*v1)[0]/e1; 319 meanN11 = (*v1)[0]/e1;
308 meanN12 = v1->Value(e2)/e2; 320 meanN12 = v1->Value(e2)/e2;
309 meanNumber = (meanN11 - meanN12)*stepFacto 321 meanNumber = (meanN11 - meanN12)*stepFactor;
310 } 322 }
311 //G4cout<<"iPlace = "<<iPlace<< " meanN11= " 323 //G4cout<<"iPlace = "<<iPlace<< " meanN11= " << meanN11
312 // << " meanN12= " << meanN12 << G4endl; 324 // << " meanN12= " << meanN12 << G4endl;
313 325
314 G4double W1 = 1.0; 326 G4double W1 = 1.0;
315 G4double W2 = 0.0; 327 G4double W2 = 0.0;
316 if(!one) { 328 if(!one) {
317 v2 = (*(fPAIxscBank[coupleIndex]))(iPlace+ 329 v2 = (*(fPAIxscBank[coupleIndex]))(iPlace+1);
318 330
319 e1 = v2->Energy(0); 331 e1 = v2->Energy(0);
320 e2 = std::min(tmax, v2->GetMaxEnergy()); 332 e2 = std::min(tmax, v2->GetMaxEnergy());
321 if(e2 >= e1) { 333 if(e2 >= e1) {
322 meanN21 = (*v2)[0]/e1; 334 meanN21 = (*v2)[0]/e1;
323 meanN22 = v2->Value(e2)/e2; 335 meanN22 = v2->Value(e2)/e2;
324 G4double E1 = fParticleEnergyVector->Ene 336 G4double E1 = fParticleEnergyVector->Energy(iPlace);
325 G4double E2 = fParticleEnergyVector->Ene 337 G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
326 G4double W = 1.0/(E2 - E1); 338 G4double W = 1.0/(E2 - E1);
327 W1 = (E2 - scaledTkin)*W; 339 W1 = (E2 - scaledTkin)*W;
328 W2 = (scaledTkin - E1)*W; 340 W2 = (scaledTkin - E1)*W;
329 meanNumber *= W1; 341 meanNumber *= W1;
330 meanNumber += (meanN21 - meanN22)*stepFa 342 meanNumber += (meanN21 - meanN22)*stepFactor*W2;
331 } 343 }
332 } 344 }
333 345
334 if(meanNumber < 0.0) { return loss; } 346 if(meanNumber < 0.0) { return loss; }
335 G4int numOfCollisions = (G4int)G4Poisson(mea << 347 G4int numOfCollisions = G4Poisson(meanNumber);
336 348
337 //G4cout << "meanNumber= " << meanNumber << 349 //G4cout << "meanNumber= " << meanNumber << " N= " << numOfCollisions << G4endl;
338 350
339 if(0 == numOfCollisions) { return loss; } 351 if(0 == numOfCollisions) { return loss; }
340 352
341 G4double position, omega, omega2; 353 G4double position, omega, omega2;
342 for(G4int i=0; i< numOfCollisions; ++i) { 354 for(G4int i=0; i< numOfCollisions; ++i) {
343 G4double rand = G4UniformRand(); 355 G4double rand = G4UniformRand();
344 position = meanN12 + (meanN11 - meanN12)*r 356 position = meanN12 + (meanN11 - meanN12)*rand;
345 omega = GetEnergyTransfer(coupleIndex, iPl 357 omega = GetEnergyTransfer(coupleIndex, iPlace, position);
346 //G4cout << "omega(keV)= " << omega/keV << 358 //G4cout << "omega(keV)= " << omega/keV << G4endl;
347 if(!one) { 359 if(!one) {
348 position = meanN22 + (meanN21 - meanN22) 360 position = meanN22 + (meanN21 - meanN22)*rand;
349 omega2 = GetEnergyTransfer(coupleIndex, 361 omega2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
350 omega *= W1; 362 omega *= W1;
351 omega += omega2*W2; 363 omega += omega2*W2;
352 } 364 }
353 //G4cout << "omega(keV)= " << omega/keV << 365 //G4cout << "omega(keV)= " << omega/keV << G4endl;
354 366
355 loss += omega; 367 loss += omega;
356 if(loss > kinEnergy) { break; } 368 if(loss > kinEnergy) { break; }
357 } 369 }
358 370
359 //G4cout<<"PAIModelData AlongStepLoss = "<<l 371 //G4cout<<"PAIModelData AlongStepLoss = "<<loss/keV<<" keV"<<G4endl;
360 if(loss > kinEnergy) { loss = kinEnergy; } 372 if(loss > kinEnergy) { loss = kinEnergy; }
361 else if(loss < 0.) { loss = 0.; } 373 else if(loss < 0.) { loss = 0.; }
362 return loss; 374 return loss;
363 } 375 }
364 376
365 ////////////////////////////////////////////// 377 ///////////////////////////////////////////////////////////////////////
366 // 378 //
367 // Returns post step PAI energy transfer > cut 379 // Returns post step PAI energy transfer > cut electron energy
368 // according to passed scaled kinetic energy o 380 // according to passed scaled kinetic energy of particle
369 381
370 G4double G4PAIModelData::SamplePostStepTransfe 382 G4double G4PAIModelData::SamplePostStepTransfer(G4int coupleIndex,
371 G4double scaledTkin, 383 G4double scaledTkin,
372 G4double tmin, 384 G4double tmin,
373 G4double tmax) const 385 G4double tmax) const
374 { 386 {
375 //G4cout<<"=== G4PAIModelData::SamplePostSte 387 //G4cout<<"=== G4PAIModelData::SamplePostStepTransfer idx= "<< coupleIndex
376 // << " Tkin= " << scaledTkin << " Tmax= " 388 // << " Tkin= " << scaledTkin << " Tmax= " << tmax << G4endl;
377 G4double transfer = 0.0; 389 G4double transfer = 0.0;
378 G4double rand = G4UniformRand(); 390 G4double rand = G4UniformRand();
379 391
380 std::size_t nPlace = fParticleEnergyVector-> << 392 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
381 std::size_t iPlace = fParticleEnergyVector-> << 393 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
382 394
383 G4bool one = true; 395 G4bool one = true;
384 if(scaledTkin >= fParticleEnergyVector->Ener 396 if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
385 else if(scaledTkin > fParticleEnergyVector-> 397 else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
386 one = false; 398 one = false;
387 } 399 }
388 G4PhysicsTable* table = fPAIxscBank[coupleIn 400 G4PhysicsTable* table = fPAIxscBank[coupleIndex];
389 G4PhysicsVector* v1 = (*table)[iPlace]; 401 G4PhysicsVector* v1 = (*table)[iPlace];
390 402
391 G4double emin = std::max(tmin, v1->Energy(0) 403 G4double emin = std::max(tmin, v1->Energy(0));
392 G4double emax = std::min(tmax, v1->GetMaxEne 404 G4double emax = std::min(tmax, v1->GetMaxEnergy());
393 if(emax < emin) { return transfer; } 405 if(emax < emin) { return transfer; }
394 406
395 G4double dNdx1 = v1->Value(emin)/emin; 407 G4double dNdx1 = v1->Value(emin)/emin;
396 G4double dNdx2 = v1->Value(emax)/emax; 408 G4double dNdx2 = v1->Value(emax)/emax;
397 /* 409 /*
398 G4cout << "iPlace= " << iPlace << " nPlace= 410 G4cout << "iPlace= " << iPlace << " nPlace= " << nPlace
399 << " emin= " << emin << " emax= " << emax 411 << " emin= " << emin << " emax= " << emax
400 << " dNdx1= " << dNdx1 << " dNdx2= " << dNd 412 << " dNdx1= " << dNdx1 << " dNdx2= " << dNdx2
401 << " one: " << one << G4endl; 413 << " one: " << one << G4endl;
402 */ 414 */
403 G4double position = dNdx2 + (dNdx1 - dNdx2)* 415 G4double position = dNdx2 + (dNdx1 - dNdx2)*rand;
404 transfer = GetEnergyTransfer(coupleIndex, iP 416 transfer = GetEnergyTransfer(coupleIndex, iPlace, position);
405 417
406 //G4cout<<"PAImodel PostStepTransfer = "<<tr 418 //G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"
407 // << " position= " << position << G4endl; 419 // << " position= " << position << G4endl;
408 420
409 if(!one) { 421 if(!one) {
410 422
411 G4PhysicsVector* v2 = (*table)[iPlace+1]; 423 G4PhysicsVector* v2 = (*table)[iPlace+1];
412 emin = std::max(tmin, v2->Energy(0)); 424 emin = std::max(tmin, v2->Energy(0));
413 emax = std::min(tmax, v2->GetMaxEnergy()); 425 emax = std::min(tmax, v2->GetMaxEnergy());
414 if(emin <= emax) { 426 if(emin <= emax) {
415 dNdx1 = v2->Value(emin)/emin; 427 dNdx1 = v2->Value(emin)/emin;
416 dNdx2 = v2->Value(emax)/emax; 428 dNdx2 = v2->Value(emax)/emax;
417 429
418 //G4cout << " emax2= " << emax 430 //G4cout << " emax2= " << emax
419 // << " dNdx2= " << dNdx2 << " dNdx1 431 // << " dNdx2= " << dNdx2 << " dNdx1= " << dNdx1 << G4endl;
420 432
421 G4double E1 = fParticleEnergyVector->Ene 433 G4double E1 = fParticleEnergyVector->Energy(iPlace);
422 G4double E2 = fParticleEnergyVector->Ene 434 G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
423 G4double W = 1.0/(E2 - E1); 435 G4double W = 1.0/(E2 - E1);
424 G4double W1 = (E2 - scaledTkin)*W; 436 G4double W1 = (E2 - scaledTkin)*W;
425 G4double W2 = (scaledTkin - E1)*W; 437 G4double W2 = (scaledTkin - E1)*W;
426 438
427 //G4cout<< "E1= " << E1 << " E2= " << E2 439 //G4cout<< "E1= " << E1 << " E2= " << E2 <<" iPlace= " << iPlace
428 // << " W1= " << W1 << " W2= " << W2 440 // << " W1= " << W1 << " W2= " << W2 <<G4endl;
429 441
430 position = dNdx2 + (dNdx1 - dNdx2)*rand; 442 position = dNdx2 + (dNdx1 - dNdx2)*rand;
431 G4double tr2 = GetEnergyTransfer(coupleI 443 G4double tr2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
432 444
433 //G4cout<<"PAImodel PostStepTransfer1 = 445 //G4cout<<"PAImodel PostStepTransfer1 = "<<tr2/keV<<" keV"
434 // << " position= " << position << G4 446 // << " position= " << position << G4endl;
435 transfer *= W1; 447 transfer *= W1;
436 transfer += tr2*W2; 448 transfer += tr2*W2;
437 } 449 }
438 } 450 }
439 //G4cout<<"PAImodel PostStepTransfer = "<<tr 451 //G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"
440 // << " position= " << position << G4endl; 452 // << " position= " << position << G4endl;
441 transfer = std::max(transfer, 0.0); 453 transfer = std::max(transfer, 0.0);
442 return transfer; 454 return transfer;
443 } 455 }
444 456
445 ////////////////////////////////////////////// 457 ///////////////////////////////////////////////////////////////////////
446 // 458 //
447 // Returns PAI energy transfer according to pa 459 // Returns PAI energy transfer according to passed
448 // indexes of particle kinetic enegry and rand 460 // indexes of particle kinetic enegry and random x-section
449 461
450 G4double G4PAIModelData::GetEnergyTransfer(G4i 462 G4double G4PAIModelData::GetEnergyTransfer(G4int coupleIndex,
451 std::size_t iPlace, << 463 size_t iPlace,
452 G4double position) const 464 G4double position) const
453 { 465 {
454 G4PhysicsVector* v = (*(fPAIxscBank[coupleIn 466 G4PhysicsVector* v = (*(fPAIxscBank[coupleIndex]))(iPlace);
455 if(position*v->Energy(0) >= (*v)[0]) { retur 467 if(position*v->Energy(0) >= (*v)[0]) { return v->Energy(0); }
456 468
457 std::size_t iTransferMax = v->GetVectorLengt << 469 size_t iTransferMax = v->GetVectorLength() - 1;
458 470
459 std::size_t iTransfer; << 471 size_t iTransfer;
460 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), 472 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
461 473
462 //G4cout << "iPlace= " << iPlace << " iTrans 474 //G4cout << "iPlace= " << iPlace << " iTransferMax= " << iTransferMax << G4endl;
463 for(iTransfer=1; iTransfer<=iTransferMax; ++ 475 for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer) {
464 x2 = v->Energy(iTransfer); 476 x2 = v->Energy(iTransfer);
465 y2 = (*v)[iTransfer]/x2; 477 y2 = (*v)[iTransfer]/x2;
466 if(position >= y2) { break; } 478 if(position >= y2) { break; }
467 if(iTransfer == iTransferMax) { return v-> 479 if(iTransfer == iTransferMax) { return v->GetMaxEnergy(); }
468 } 480 }
469 481
470 x1 = v->Energy(iTransfer-1); 482 x1 = v->Energy(iTransfer-1);
471 y1 = (*v)[iTransfer-1]/x1; 483 y1 = (*v)[iTransfer-1]/x1;
472 /* 484 /*
473 G4cout << "i= " << iTransfer << " imax= " << 485 G4cout << "i= " << iTransfer << " imax= " << iTransferMax
474 << " x1= " << x1 << " x2= " << x2 486 << " x1= " << x1 << " x2= " << x2
475 << " y1= " << y1 << " y2= " << y2 << G4en 487 << " y1= " << y1 << " y2= " << y2 << G4endl;
476 */ 488 */
477 energyTransfer = x1; 489 energyTransfer = x1;
478 if ( x1 != x2 ) { 490 if ( x1 != x2 ) {
479 if ( y1 == y2 ) { 491 if ( y1 == y2 ) {
480 energyTransfer += (x2 - x1)*G4UniformRan 492 energyTransfer += (x2 - x1)*G4UniformRand();
481 } else { 493 } else {
482 if(x1*1.1 < x2) { 494 if(x1*1.1 < x2) {
483 const G4int nbins = 5; 495 const G4int nbins = 5;
484 G4double del = (x2 - x1)/G4int(nbins); 496 G4double del = (x2 - x1)/G4int(nbins);
485 x2 = x1; 497 x2 = x1;
486 for(G4int i=1; i<=nbins; ++i) { 498 for(G4int i=1; i<=nbins; ++i) {
487 x2 += del; 499 x2 += del;
488 y2 = v->Value(x2)/x2; 500 y2 = v->Value(x2)/x2;
489 if(position >= y2) { 501 if(position >= y2) {
490 break; 502 break;
491 } 503 }
492 x1 = x2; 504 x1 = x2;
493 y1 = y2; 505 y1 = y2;
494 } 506 }
495 } 507 }
496 //G4cout << "x1(keV)= " << x1/keV << " x 508 //G4cout << "x1(keV)= " << x1/keV << " x2(keV)= " << x2/keV
497 // << " y1= " << y1 << " y2= " << y2 < 509 // << " y1= " << y1 << " y2= " << y2 << " pos= " << position << G4endl;
498 energyTransfer = (y2 - y1)*x1*x2/(positi 510 energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2);
499 } 511 }
500 } 512 }
501 //G4cout << "x1(keV)= " << x1/keV << " x2(ke 513 //G4cout << "x1(keV)= " << x1/keV << " x2(keV)= " << x2/keV
502 // << " y1= " << y1 << " y2= " << y2 << " 514 // << " y1= " << y1 << " y2= " << y2 << " pos= " << position
503 // << " E(keV)= " << energyTransfer/keV << 515 // << " E(keV)= " << energyTransfer/keV << G4endl;
504 return energyTransfer; 516 return energyTransfer;
505 } 517 }
506 518
507 ////////////////////////////////////////////// 519 //////////////////////////////////////////////////////////////////////
508 520
509 521