| 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: G4PAIPhotData.cc 30 // File name: G4PAIPhotData.cc
31 // 31 //
32 // Author: V.Grichine based on G4PAIMod 32 // Author: V.Grichine based on G4PAIModelData MT
33 // 33 //
34 // Creation date: 07.10.2013 34 // Creation date: 07.10.2013
35 // 35 //
36 // Modifications: 36 // Modifications:
37 // 37 //
38 38
39 #include "G4PAIPhotData.hh" 39 #include "G4PAIPhotData.hh"
40 #include "G4PAIPhotModel.hh" 40 #include "G4PAIPhotModel.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 "G4ProductionCutsTable.hh" 49 #include "G4ProductionCutsTable.hh"
50 #include "G4SandiaTable.hh" 50 #include "G4SandiaTable.hh"
51 #include "Randomize.hh" 51 #include "Randomize.hh"
52 #include "G4Poisson.hh" 52 #include "G4Poisson.hh"
53 53
54 ////////////////////////////////////////////// 54 ////////////////////////////////////////////////////////////////////////
55 55
56 using namespace std; 56 using namespace std;
57 57
58 G4PAIPhotData::G4PAIPhotData(G4double tmin, G4 58 G4PAIPhotData::G4PAIPhotData(G4double tmin, G4double tmax, G4int ver)
59 { 59 {
60 const G4int nPerDecade = 10; 60 const G4int nPerDecade = 10;
61 const G4double lowestTkin = 50*keV; 61 const G4double lowestTkin = 50*keV;
62 const G4double highestTkin = 10*TeV; 62 const G4double highestTkin = 10*TeV;
63 63
64 // fPAIxSection.SetVerbose(ver); 64 // fPAIxSection.SetVerbose(ver);
65 65
66 fLowestKineticEnergy = std::max(tmin, lowes 66 fLowestKineticEnergy = std::max(tmin, lowestTkin);
67 fHighestKineticEnergy = tmax; 67 fHighestKineticEnergy = tmax;
68 68
69 if(tmax < 10*fLowestKineticEnergy) 69 if(tmax < 10*fLowestKineticEnergy)
70 { 70 {
71 fHighestKineticEnergy = 10*fLowestKineticE 71 fHighestKineticEnergy = 10*fLowestKineticEnergy;
72 } 72 }
73 else if(tmax > highestTkin) 73 else if(tmax > highestTkin)
74 { 74 {
75 fHighestKineticEnergy = std::max(highestTk 75 fHighestKineticEnergy = std::max(highestTkin, 10*fLowestKineticEnergy);
76 } 76 }
77 fTotBin = (G4int)(nPerDecade* 77 fTotBin = (G4int)(nPerDecade*
78 std::log10(fHighestKineticEnergy/fLowe 78 std::log10(fHighestKineticEnergy/fLowestKineticEnergy));
79 79
80 fParticleEnergyVector = new G4PhysicsLogVect 80 fParticleEnergyVector = new G4PhysicsLogVector(fLowestKineticEnergy,
81 fHighestKineticEnergy, 81 fHighestKineticEnergy,
82 fTotBin); 82 fTotBin);
83 if(0 < ver) { 83 if(0 < ver) {
84 G4cout << "### G4PAIPhotData: Nbins= " << 84 G4cout << "### G4PAIPhotData: Nbins= " << fTotBin
85 << " Tmin(MeV)= " << fLowestKineticEnergy 85 << " Tmin(MeV)= " << fLowestKineticEnergy/MeV
86 << " Tmax(GeV)= " << fHighestKineticEnerg 86 << " Tmax(GeV)= " << fHighestKineticEnergy/GeV
87 << " tmin(keV)= " << tmin/keV << G4endl; 87 << " tmin(keV)= " << tmin/keV << G4endl;
88 } 88 }
89 } 89 }
90 90
91 ////////////////////////////////////////////// 91 ////////////////////////////////////////////////////////////////////////////
92 92
93 G4PAIPhotData::~G4PAIPhotData() 93 G4PAIPhotData::~G4PAIPhotData()
94 { 94 {
95 //G4cout << "G4PAIPhotData::~G4PAIPhotData() 95 //G4cout << "G4PAIPhotData::~G4PAIPhotData() " << this << G4endl;
96 std::size_t n = fPAIxscBank.size(); 96 std::size_t n = fPAIxscBank.size();
97 if(0 < n) 97 if(0 < n)
98 { 98 {
99 for(std::size_t i=0; i<n; ++i) 99 for(std::size_t i=0; i<n; ++i)
100 { 100 {
101 if(fPAIxscBank[i]) 101 if(fPAIxscBank[i])
102 { 102 {
103 fPAIxscBank[i]->clearAndDestroy(); 103 fPAIxscBank[i]->clearAndDestroy();
104 delete fPAIxscBank[i]; 104 delete fPAIxscBank[i];
105 fPAIxscBank[i] = nullptr; 105 fPAIxscBank[i] = nullptr;
106 } 106 }
107 if(fPAIdEdxBank[i]) 107 if(fPAIdEdxBank[i])
108 { 108 {
109 fPAIdEdxBank[i]->clearAndDestroy(); 109 fPAIdEdxBank[i]->clearAndDestroy();
110 delete fPAIdEdxBank[i]; 110 delete fPAIdEdxBank[i];
111 fPAIdEdxBank[i] = nullptr; 111 fPAIdEdxBank[i] = nullptr;
112 } 112 }
113 delete fdEdxTable[i]; 113 delete fdEdxTable[i];
114 delete fdNdxCutTable[i]; 114 delete fdNdxCutTable[i];
115 fdEdxTable[i] = nullptr; 115 fdEdxTable[i] = nullptr;
116 fdNdxCutTable[i] = nullptr; 116 fdNdxCutTable[i] = nullptr;
117 } 117 }
118 } 118 }
119 delete fParticleEnergyVector; 119 delete fParticleEnergyVector;
120 fParticleEnergyVector = nullptr; 120 fParticleEnergyVector = nullptr;
121 //G4cout << "G4PAIPhotData::~G4PAIPhotData() 121 //G4cout << "G4PAIPhotData::~G4PAIPhotData() done for " << this << G4endl;
122 } 122 }
123 123
124 ////////////////////////////////////////////// 124 ///////////////////////////////////////////////////////////////////////////////
125 125
126 void G4PAIPhotData::Initialise(const G4Materia 126 void G4PAIPhotData::Initialise(const G4MaterialCutsCouple* couple,
127 G4double cut, 127 G4double cut, G4PAIPhotModel* model)
128 { 128 {
129 G4ProductionCutsTable* theCoupleTable= 129 G4ProductionCutsTable* theCoupleTable=
130 G4ProductionCutsTable::GetProductionCu 130 G4ProductionCutsTable::GetProductionCutsTable();
131 G4int numOfCouples = (G4int)theCoupleTable-> 131 G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
132 G4int jMatCC; 132 G4int jMatCC;
133 133
134 for (jMatCC = 0; jMatCC < numOfCouples; ++jM 134 for (jMatCC = 0; jMatCC < numOfCouples; ++jMatCC )
135 { 135 {
136 if( couple == theCoupleTable->GetMaterialC 136 if( couple == theCoupleTable->GetMaterialCutsCouple(jMatCC) ) break;
137 } 137 }
138 if( jMatCC == numOfCouples && jMatCC > 0 ) - 138 if( jMatCC == numOfCouples && jMatCC > 0 ) --jMatCC;
139 139
140 const vector<G4double>* deltaCutInKineticEn 140 const vector<G4double>* deltaCutInKineticEnergy = theCoupleTable->GetEnergyCutsVector(idxG4ElectronCut);
141 const vector<G4double>* photonCutInKineticE 141 const vector<G4double>* photonCutInKineticEnergy = theCoupleTable->GetEnergyCutsVector(idxG4GammaCut);
142 G4double deltaCutInKineticEnergyNow = (*del 142 G4double deltaCutInKineticEnergyNow = (*deltaCutInKineticEnergy)[jMatCC];
143 G4double photonCutInKineticEnergyNow = (*pho 143 G4double photonCutInKineticEnergyNow = (*photonCutInKineticEnergy)[jMatCC];
144 /* << 144
145 G4cout<<"G4PAIPhotData::Initialise: "<<"cut 145 G4cout<<"G4PAIPhotData::Initialise: "<<"cut = "<<cut/keV<<" keV; cutEl = "
146 <<deltaCutInKineticEnergyNow/keV<<" ke 146 <<deltaCutInKineticEnergyNow/keV<<" keV; cutPh = "
147 <<photonCutInKineticEnergyNow/keV<<" keV"<<G 147 <<photonCutInKineticEnergyNow/keV<<" keV"<<G4endl;
148 */ << 148
149 // if( deltaCutInKineticEnergyNow != cut ) d 149 // if( deltaCutInKineticEnergyNow != cut ) deltaCutInKineticEnergyNow = cut; // exception??
150 150
151 auto dEdxCutVector = 151 auto dEdxCutVector =
152 new G4PhysicsLogVector(fLowestKineticEnerg 152 new G4PhysicsLogVector(fLowestKineticEnergy,
153 fHighestKineticEnergy, 153 fHighestKineticEnergy,
154 fTotBin); 154 fTotBin);
155 155
156 auto dNdxCutVector = 156 auto dNdxCutVector =
157 new G4PhysicsLogVector(fLowestKineticEnerg 157 new G4PhysicsLogVector(fLowestKineticEnergy,
158 fHighestKineticEnergy, 158 fHighestKineticEnergy,
159 fTotBin); 159 fTotBin);
160 auto dNdxCutPhotonVector = 160 auto dNdxCutPhotonVector =
161 new G4PhysicsLogVector(fLowestKineticEnerg 161 new G4PhysicsLogVector(fLowestKineticEnergy,
162 fHighestKineticEnergy, 162 fHighestKineticEnergy,
163 fTotBin); 163 fTotBin);
164 auto dNdxCutPlasmonVector = 164 auto dNdxCutPlasmonVector =
165 new G4PhysicsLogVector(fLowestKineticEnerg 165 new G4PhysicsLogVector(fLowestKineticEnergy,
166 fHighestKineticEnergy, 166 fHighestKineticEnergy,
167 fTotBin); 167 fTotBin);
168 168
169 const G4Material* mat = couple->GetMaterial( 169 const G4Material* mat = couple->GetMaterial();
170 fSandia.Initialize(const_cast<G4Material*>(m 170 fSandia.Initialize(const_cast<G4Material*>(mat));
171 171
172 auto PAItransferTable = new G4PhysicsTable(f 172 auto PAItransferTable = new G4PhysicsTable(fTotBin+1);
173 auto PAIphotonTable = new G4PhysicsTable(fTo 173 auto PAIphotonTable = new G4PhysicsTable(fTotBin+1);
174 auto PAIplasmonTable = new G4PhysicsTable(fT 174 auto PAIplasmonTable = new G4PhysicsTable(fTotBin+1);
175 175
176 auto PAIdEdxTable = new G4PhysicsTable(fTotB 176 auto PAIdEdxTable = new G4PhysicsTable(fTotBin+1);
177 auto dEdxMeanVector = 177 auto dEdxMeanVector =
178 new G4PhysicsLogVector(fLowestKineticEnerg 178 new G4PhysicsLogVector(fLowestKineticEnergy,
179 fHighestKineticEnergy, 179 fHighestKineticEnergy,
180 fTotBin); 180 fTotBin);
181 181
182 // low energy Sandia interval 182 // low energy Sandia interval
183 G4double Tmin = fSandia.GetSandiaMatTablePAI 183 G4double Tmin = fSandia.GetSandiaMatTablePAI(0,0);
184 184
185 // energy safety 185 // energy safety
186 const G4double deltaLow = 100.*eV; 186 const G4double deltaLow = 100.*eV;
187 187
188 for (G4int i = 0; i <= fTotBin; ++i) 188 for (G4int i = 0; i <= fTotBin; ++i)
189 { 189 {
190 G4double kinEnergy = fParticleEnergyVector 190 G4double kinEnergy = fParticleEnergyVector->Energy(i);
191 G4double Tmax = model->ComputeMaxEnergy(ki 191 G4double Tmax = model->ComputeMaxEnergy(kinEnergy);
192 G4double tau = kinEnergy/proton_mass_c2; 192 G4double tau = kinEnergy/proton_mass_c2;
193 G4double bg2 = tau*( tau + 2. ); 193 G4double bg2 = tau*( tau + 2. );
194 194
195 if ( Tmax < Tmin + deltaLow ) Tmax = Tmin 195 if ( Tmax < Tmin + deltaLow ) Tmax = Tmin + deltaLow;
196 196
197 fPAIxSection.Initialize( mat, Tmax, bg2, & 197 fPAIxSection.Initialize( mat, Tmax, bg2, &fSandia);
198 198
199 //G4cout << i << ". TransferMax(keV)= "<< 199 //G4cout << i << ". TransferMax(keV)= "<< Tmax/keV << " cut(keV)= "
200 // << cut/keV << " E(MeV)= " << kinEn 200 // << cut/keV << " E(MeV)= " << kinEnergy/MeV << G4endl;
201 201
202 G4int n = fPAIxSection.GetSplineSize(); 202 G4int n = fPAIxSection.GetSplineSize();
203 203
204 auto transferVector = new G4PhysicsFreeVec 204 auto transferVector = new G4PhysicsFreeVector(n);
205 auto photonVector = new G4PhysicsFreeVec 205 auto photonVector = new G4PhysicsFreeVector(n);
206 auto plasmonVector = new G4PhysicsFreeVec 206 auto plasmonVector = new G4PhysicsFreeVector(n);
207 207
208 auto dEdxVector = new G4PhysicsFreeVec 208 auto dEdxVector = new G4PhysicsFreeVector(n);
209 209
210 for( G4int k = 0; k < n; k++ ) 210 for( G4int k = 0; k < n; k++ )
211 { 211 {
212 G4double t = fPAIxSection.GetSplineEnerg 212 G4double t = fPAIxSection.GetSplineEnergy(k+1);
213 213
214 transferVector->PutValue(k , t, 214 transferVector->PutValue(k , t,
215 t*fPAIxSection. 215 t*fPAIxSection.GetIntegralPAIxSection(k+1));
216 photonVector->PutValue(k , t, 216 photonVector->PutValue(k , t,
217 t*fPAIxSection. 217 t*fPAIxSection.GetIntegralCerenkov(k+1));
218 plasmonVector->PutValue(k , t, 218 plasmonVector->PutValue(k , t,
219 t*fPAIxSection. 219 t*fPAIxSection.GetIntegralPlasmon(k+1));
220 220
221 dEdxVector->PutValue(k, t, fPAIxSection. 221 dEdxVector->PutValue(k, t, fPAIxSection.GetIntegralPAIdEdx(k+1));
222 } 222 }
223 // G4cout << *transferVector << G4endl; 223 // G4cout << *transferVector << G4endl;
224 224
225 G4double ionloss = std::max(fPAIxSection.G << 225 G4double ionloss = fPAIxSection.GetMeanEnergyLoss();// total <dE/dx>
>> 226
>> 227 if(ionloss < 0.0) ionloss = 0.0;
>> 228
226 dEdxMeanVector->PutValue(i,ionloss); 229 dEdxMeanVector->PutValue(i,ionloss);
227 230
228 G4double dNdxCut = transferVector->Value(d 231 G4double dNdxCut = transferVector->Value(deltaCutInKineticEnergyNow)/deltaCutInKineticEnergyNow;
229 G4double dNdxCutPhoton = photonVector->Val 232 G4double dNdxCutPhoton = photonVector->Value(photonCutInKineticEnergyNow)/photonCutInKineticEnergyNow;
230 G4double dNdxCutPlasmon = plasmonVector->V 233 G4double dNdxCutPlasmon = plasmonVector->Value(deltaCutInKineticEnergyNow)/deltaCutInKineticEnergyNow;
231 234
232 G4double dEdxCut = dEdxVector->Value(cut)/ 235 G4double dEdxCut = dEdxVector->Value(cut)/cut;
233 //G4cout << "i= " << i << " x= " << dNdxCu 236 //G4cout << "i= " << i << " x= " << dNdxCut << G4endl;
234 237
235 if(dNdxCut < 0.0) { dNdxCut = 0.0; } 238 if(dNdxCut < 0.0) { dNdxCut = 0.0; }
236 if(dNdxCutPhoton < 0.0) { dNdxCutPhoton = 239 if(dNdxCutPhoton < 0.0) { dNdxCutPhoton = 0.0; }
237 if(dNdxCutPlasmon < 0.0) { dNdxCutPlasmon 240 if(dNdxCutPlasmon < 0.0) { dNdxCutPlasmon = 0.0; }
238 241
239 dNdxCutVector->PutValue(i, dNdxCut); 242 dNdxCutVector->PutValue(i, dNdxCut);
240 dNdxCutPhotonVector->PutValue(i, dNdxCutPh 243 dNdxCutPhotonVector->PutValue(i, dNdxCutPhoton);
241 dNdxCutPlasmonVector->PutValue(i, dNdxCutP 244 dNdxCutPlasmonVector->PutValue(i, dNdxCutPlasmon);
242 245
243 dEdxCutVector->PutValue(i, dEdxCut); 246 dEdxCutVector->PutValue(i, dEdxCut);
244 247
245 PAItransferTable->insertAt(i,transferVecto 248 PAItransferTable->insertAt(i,transferVector);
246 PAIphotonTable->insertAt(i,photonVector); 249 PAIphotonTable->insertAt(i,photonVector);
247 PAIplasmonTable->insertAt(i,plasmonVector) 250 PAIplasmonTable->insertAt(i,plasmonVector);
248 PAIdEdxTable->insertAt(i,dEdxVector); 251 PAIdEdxTable->insertAt(i,dEdxVector);
249 252
250 } // end of Tkin loop 253 } // end of Tkin loop
251 254
252 fPAIxscBank.push_back(PAItransferTable); 255 fPAIxscBank.push_back(PAItransferTable);
253 fPAIphotonBank.push_back(PAIphotonTable); 256 fPAIphotonBank.push_back(PAIphotonTable);
254 fPAIplasmonBank.push_back(PAIplasmonTable); 257 fPAIplasmonBank.push_back(PAIplasmonTable);
255 258
256 fPAIdEdxBank.push_back(PAIdEdxTable); 259 fPAIdEdxBank.push_back(PAIdEdxTable);
257 fdEdxTable.push_back(dEdxMeanVector); 260 fdEdxTable.push_back(dEdxMeanVector);
258 261
259 fdNdxCutTable.push_back(dNdxCutVector); 262 fdNdxCutTable.push_back(dNdxCutVector);
260 fdNdxCutPhotonTable.push_back(dNdxCutPhotonV 263 fdNdxCutPhotonTable.push_back(dNdxCutPhotonVector);
261 fdNdxCutPlasmonTable.push_back(dNdxCutPlasmo 264 fdNdxCutPlasmonTable.push_back(dNdxCutPlasmonVector);
262 265
263 fdEdxCutTable.push_back(dEdxCutVector); 266 fdEdxCutTable.push_back(dEdxCutVector);
264 } 267 }
265 268
266 ////////////////////////////////////////////// 269 //////////////////////////////////////////////////////////////////////////////
267 270
268 G4double G4PAIPhotData::DEDXPerVolume(G4int co 271 G4double G4PAIPhotData::DEDXPerVolume(G4int coupleIndex, G4double scaledTkin,
269 G4double cut) const 272 G4double cut) const
270 { 273 {
271 // VI: iPlace is the low edge index of the b 274 // VI: iPlace is the low edge index of the bin
272 // iPlace is in interval from 0 to (N-1) 275 // iPlace is in interval from 0 to (N-1)
273 std::size_t iPlace = fParticleEnergyVector-> 276 std::size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
274 std::size_t nPlace = fParticleEnergyVector-> 277 std::size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
275 278
276 G4bool one = true; 279 G4bool one = true;
277 if(scaledTkin >= fParticleEnergyVector->Ener 280 if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; }
278 else if(scaledTkin > fParticleEnergyVector-> 281 else if(scaledTkin > fParticleEnergyVector->Energy(0)) {
279 one = false; 282 one = false;
280 } 283 }
281 284
282 // VI: apply interpolation of the vector 285 // VI: apply interpolation of the vector
283 G4double dEdx = fdEdxTable[coupleIndex]->Val 286 G4double dEdx = fdEdxTable[coupleIndex]->Value(scaledTkin);
284 G4double del = (*(fPAIdEdxBank[coupleIndex] 287 G4double del = (*(fPAIdEdxBank[coupleIndex]))(iPlace)->Value(cut);
285 if(!one) { 288 if(!one) {
286 G4double del2 = (*(fPAIdEdxBank[coupleInde 289 G4double del2 = (*(fPAIdEdxBank[coupleIndex]))(iPlace+1)->Value(cut);
287 G4double E1 = fParticleEnergyVector->Energ 290 G4double E1 = fParticleEnergyVector->Energy(iPlace);
288 G4double E2 = fParticleEnergyVector->Energ 291 G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
289 G4double W = 1.0/(E2 - E1); 292 G4double W = 1.0/(E2 - E1);
290 G4double W1 = (E2 - scaledTkin)*W; 293 G4double W1 = (E2 - scaledTkin)*W;
291 G4double W2 = (scaledTkin - E1)*W; 294 G4double W2 = (scaledTkin - E1)*W;
292 del *= W1; 295 del *= W1;
293 del += W2*del2; 296 del += W2*del2;
294 } 297 }
295 dEdx -= del; 298 dEdx -= del;
296 299
297 if( dEdx < 0.) { dEdx = 0.; } 300 if( dEdx < 0.) { dEdx = 0.; }
298 return dEdx; 301 return dEdx;
299 } 302 }
300 303
301 ////////////////////////////////////////////// 304 /////////////////////////////////////////////////////////////////////////
302 305
303 G4double 306 G4double
304 G4PAIPhotData::CrossSectionPerVolume(G4int cou 307 G4PAIPhotData::CrossSectionPerVolume(G4int coupleIndex,
305 G4double scaledTkin, 308 G4double scaledTkin,
306 G4double tcut, G4double tmax) co 309 G4double tcut, G4double tmax) const
307 { 310 {
308 G4double cross, xscEl, xscEl2, xscPh, xscPh2 311 G4double cross, xscEl, xscEl2, xscPh, xscPh2;
309 312
310 cross=tcut+tmax; 313 cross=tcut+tmax;
311 314
312 // iPlace is in interval from 0 to (N-1) 315 // iPlace is in interval from 0 to (N-1)
313 316
314 std::size_t iPlace = fParticleEnergyVector-> 317 std::size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
315 std::size_t nPlace = fParticleEnergyVector-> 318 std::size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
316 319
317 G4bool one = true; 320 G4bool one = true;
318 321
319 if( scaledTkin >= fParticleEnergyVector 322 if( scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
320 else if( scaledTkin > fParticleEnergyVector- 323 else if( scaledTkin > fParticleEnergyVector->Energy(0) ) one = false;
321 324
322 325
323 xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex] 326 xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace);
324 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex]) 327 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace);
325 328
326 xscPh = xscPh2; 329 xscPh = xscPh2;
327 xscEl = xscEl2; 330 xscEl = xscEl2;
328 331
329 cross = xscPh + xscEl; 332 cross = xscPh + xscEl;
330 333
331 if( !one ) 334 if( !one )
332 { 335 {
333 xscEl2 = (*fdNdxCutPlasmonTable[coupleInde 336 xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace+1);
334 337
335 G4double E1 = fParticleEnergyVector->Energ 338 G4double E1 = fParticleEnergyVector->Energy(iPlace);
336 G4double E2 = fParticleEnergyVector->Energ 339 G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
337 340
338 G4double W = 1.0/(E2 - E1); 341 G4double W = 1.0/(E2 - E1);
339 G4double W1 = (E2 - scaledTkin)*W; 342 G4double W1 = (E2 - scaledTkin)*W;
340 G4double W2 = (scaledTkin - E1)*W; 343 G4double W2 = (scaledTkin - E1)*W;
341 344
342 xscEl *= W1; 345 xscEl *= W1;
343 xscEl += W2*xscEl2; 346 xscEl += W2*xscEl2;
344 347
345 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex 348 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace+1);
346 349
347 E1 = fParticleEnergyVector->Energy(iPlace) 350 E1 = fParticleEnergyVector->Energy(iPlace);
348 E2 = fParticleEnergyVector->Energy(iPlace+ 351 E2 = fParticleEnergyVector->Energy(iPlace+1);
349 352
350 W = 1.0/(E2 - E1); 353 W = 1.0/(E2 - E1);
351 W1 = (E2 - scaledTkin)*W; 354 W1 = (E2 - scaledTkin)*W;
352 W2 = (scaledTkin - E1)*W; 355 W2 = (scaledTkin - E1)*W;
353 356
354 xscPh *= W1; 357 xscPh *= W1;
355 xscPh += W2*xscPh2; 358 xscPh += W2*xscPh2;
356 359
357 cross = xscEl + xscPh; 360 cross = xscEl + xscPh;
358 } 361 }
359 if( cross < 0.0) cross = 0.0; 362 if( cross < 0.0) cross = 0.0;
360 363
361 return cross; 364 return cross;
362 } 365 }
363 366
364 ////////////////////////////////////////////// 367 /////////////////////////////////////////////////////////////////////////
365 368
366 G4double 369 G4double
367 G4PAIPhotData::GetPlasmonRatio(G4int coupleInd 370 G4PAIPhotData::GetPlasmonRatio(G4int coupleIndex, G4double scaledTkin) const
368 { 371 {
369 G4double cross, xscEl, xscEl2, xscPh, xscPh2 372 G4double cross, xscEl, xscEl2, xscPh, xscPh2, plRatio;
370 // iPlace is in interval from 0 to (N-1) 373 // iPlace is in interval from 0 to (N-1)
371 374
372 std::size_t iPlace = fParticleEnergyVector-> 375 std::size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
373 std::size_t nPlace = fParticleEnergyVector-> 376 std::size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
374 377
375 G4bool one = true; 378 G4bool one = true;
376 379
377 if( scaledTkin >= fParticleEnergyVector 380 if( scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
378 else if( scaledTkin > fParticleEnergyVector- 381 else if( scaledTkin > fParticleEnergyVector->Energy(0) ) one = false;
379 382
380 383
381 xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex] 384 xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace);
382 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex]) 385 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace);
383 386
384 xscPh = xscPh2; 387 xscPh = xscPh2;
385 xscEl = xscEl2; 388 xscEl = xscEl2;
386 389
387 cross = xscPh + xscEl; 390 cross = xscPh + xscEl;
388 391
389 if( !one ) 392 if( !one )
390 { 393 {
391 xscEl2 = (*fdNdxCutPlasmonTable[coupleInde 394 xscEl2 = (*fdNdxCutPlasmonTable[coupleIndex])(iPlace+1);
392 395
393 G4double E1 = fParticleEnergyVector->Energ 396 G4double E1 = fParticleEnergyVector->Energy(iPlace);
394 G4double E2 = fParticleEnergyVector->Energ 397 G4double E2 = fParticleEnergyVector->Energy(iPlace+1);
395 398
396 G4double W = 1.0/(E2 - E1); 399 G4double W = 1.0/(E2 - E1);
397 G4double W1 = (E2 - scaledTkin)*W; 400 G4double W1 = (E2 - scaledTkin)*W;
398 G4double W2 = (scaledTkin - E1)*W; 401 G4double W2 = (scaledTkin - E1)*W;
399 402
400 xscEl *= W1; 403 xscEl *= W1;
401 xscEl += W2*xscEl2; 404 xscEl += W2*xscEl2;
402 405
403 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex 406 xscPh2 = (*fdNdxCutPhotonTable[coupleIndex])(iPlace+1);
404 407
405 E1 = fParticleEnergyVector->Energy(iPlace) 408 E1 = fParticleEnergyVector->Energy(iPlace);
406 E2 = fParticleEnergyVector->Energy(iPlace+ 409 E2 = fParticleEnergyVector->Energy(iPlace+1);
407 410
408 W = 1.0/(E2 - E1); 411 W = 1.0/(E2 - E1);
409 W1 = (E2 - scaledTkin)*W; 412 W1 = (E2 - scaledTkin)*W;
410 W2 = (scaledTkin - E1)*W; 413 W2 = (scaledTkin - E1)*W;
411 414
412 xscPh *= W1; 415 xscPh *= W1;
413 xscPh += W2*xscPh2; 416 xscPh += W2*xscPh2;
414 417
415 cross = xscEl + xscPh; 418 cross = xscEl + xscPh;
416 } 419 }
417 if( cross <= 0.0) 420 if( cross <= 0.0)
418 { 421 {
419 plRatio = 2.0; 422 plRatio = 2.0;
420 } 423 }
421 else 424 else
422 { 425 {
423 plRatio = xscEl/cross; 426 plRatio = xscEl/cross;
424 427
425 if( plRatio > 1. || plRatio < 0.) plRatio 428 if( plRatio > 1. || plRatio < 0.) plRatio = 2.0;
426 } 429 }
427 return plRatio; 430 return plRatio;
428 } 431 }
429 432
430 ////////////////////////////////////////////// 433 ///////////////////////////////////////////////////////////////////////
431 434
432 G4double G4PAIPhotData::SampleAlongStepTransfe 435 G4double G4PAIPhotData::SampleAlongStepTransfer(G4int coupleIndex,
433 436 G4double kinEnergy,
434 G4double scaledTkin, 437 G4double scaledTkin,
435 G4double stepFactor) const 438 G4double stepFactor) const
436 { 439 {
437 G4double loss = 0.0; 440 G4double loss = 0.0;
438 G4double omega; 441 G4double omega;
439 G4double position, E1, E2, W1, W2, W, dNdxCu 442 G4double position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber;
440 443
441 std::size_t iPlace = fParticleEnergyVector-> 444 std::size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
442 std::size_t nPlace = fParticleEnergyVector-> 445 std::size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
443 446
444 G4bool one = true; 447 G4bool one = true;
445 448
446 if (scaledTkin >= fParticleEnergyVector- 449 if (scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
447 else if(scaledTkin > fParticleEnergyVector-> 450 else if(scaledTkin > fParticleEnergyVector->Energy(0)) one = false;
448 451
449 G4PhysicsLogVector* vcut = fdNdxCutTable[cou 452 G4PhysicsLogVector* vcut = fdNdxCutTable[coupleIndex];
450 G4PhysicsVector* v1 = (*(fPAIxscBank[co 453 G4PhysicsVector* v1 = (*(fPAIxscBank[coupleIndex]))(iPlace);
451 G4PhysicsVector* v2 = nullptr; 454 G4PhysicsVector* v2 = nullptr;
452 455
453 dNdxCut1 = (*vcut)[iPlace]; 456 dNdxCut1 = (*vcut)[iPlace];
454 G4double e1 = v1->Energy(0); 457 G4double e1 = v1->Energy(0);
455 G4double e2 = e1; 458 G4double e2 = e1;
456 459
457 G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*s 460 G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;
458 461
459 meanNumber = meanN1; 462 meanNumber = meanN1;
460 463
461 // G4cout<<"iPlace = "<<iPlace<< " meanN1= " 464 // G4cout<<"iPlace = "<<iPlace<< " meanN1= " << meanN1
462 // <<" (*v1)[0]/e1 = "<<(*v1)[0]/e1<< " 465 // <<" (*v1)[0]/e1 = "<<(*v1)[0]/e1<< " dNdxCut1= " << dNdxCut1 << G4endl;
463 466
464 dNdxCut2 = dNdxCut1; 467 dNdxCut2 = dNdxCut1;
465 W1 = 1.0; 468 W1 = 1.0;
466 W2 = 0.0; 469 W2 = 0.0;
467 if(!one) 470 if(!one)
468 { 471 {
469 v2 = (*(fPAIxscBank[coupleIndex]))(iPlace+ 472 v2 = (*(fPAIxscBank[coupleIndex]))(iPlace+1);
470 dNdxCut2 = (*vcut)[iPlace+1]; 473 dNdxCut2 = (*vcut)[iPlace+1];
471 e2 = v2->Energy(0); 474 e2 = v2->Energy(0);
472 475
473 G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2) 476 G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;
474 477
475 E1 = fParticleEnergyVector->Energy(iPlace) 478 E1 = fParticleEnergyVector->Energy(iPlace);
476 E2 = fParticleEnergyVector->Energy(iPlace+ 479 E2 = fParticleEnergyVector->Energy(iPlace+1);
477 W = 1.0/(E2 - E1); 480 W = 1.0/(E2 - E1);
478 W1 = (E2 - scaledTkin)*W; 481 W1 = (E2 - scaledTkin)*W;
479 W2 = (scaledTkin - E1)*W; 482 W2 = (scaledTkin - E1)*W;
480 meanNumber = W1*meanN1 + W2*meanN2; 483 meanNumber = W1*meanN1 + W2*meanN2;
481 484
482 //G4cout<<"meanN= " << meanNumber << " me 485 //G4cout<<"meanN= " << meanNumber << " meanN2= " << meanN2
483 // << " dNdxCut2= " << dNdxCut2 << G4en 486 // << " dNdxCut2= " << dNdxCut2 << G4endl;
484 } 487 }
485 if( meanNumber <= 0.0) return 0.0; 488 if( meanNumber <= 0.0) return 0.0;
486 489
487 G4int numOfCollisions = (G4int)G4Poisson(mea 490 G4int numOfCollisions = (G4int)G4Poisson(meanNumber);
488 491
489 //G4cout << "N= " << numOfCollisions << G4en 492 //G4cout << "N= " << numOfCollisions << G4endl;
490 493
491 if( 0 == numOfCollisions) return 0.0; 494 if( 0 == numOfCollisions) return 0.0;
492 495
493 for(G4int i=0; i< numOfCollisions; ++i) 496 for(G4int i=0; i< numOfCollisions; ++i)
494 { 497 {
495 G4double rand = G4UniformRand(); 498 G4double rand = G4UniformRand();
496 position = dNdxCut1 + ((*v1)[0]/e1 - dNdxC 499 position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
497 omega = GetEnergyTransfer(coupleIndex, iPl 500 omega = GetEnergyTransfer(coupleIndex, iPlace, position);
498 501
499 //G4cout << "omega(keV)= " << omega/keV << 502 //G4cout << "omega(keV)= " << omega/keV << G4endl;
500 503
501 if(!one) 504 if(!one)
502 { 505 {
503 position = dNdxCut2 + ((*v2)[0]/e2 - dNd 506 position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
504 G4double omega2 = GetEnergyTransfer(coup 507 G4double omega2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
505 omega = omega*W1 + omega2*W2; 508 omega = omega*W1 + omega2*W2;
506 } 509 }
507 //G4cout << "omega(keV)= " << omega/keV << 510 //G4cout << "omega(keV)= " << omega/keV << G4endl;
508 511
509 loss += omega; 512 loss += omega;
510 if( loss > kinEnergy) { break; } 513 if( loss > kinEnergy) { break; }
511 } 514 }
512 515
513 // G4cout<<"PAIPhotData AlongStepLoss = "<<l 516 // G4cout<<"PAIPhotData AlongStepLoss = "<<loss/keV<<" keV, on step = "
514 //<<step/mm<<" mm"<<G4endl; 517 //<<step/mm<<" mm"<<G4endl;
515 518
516 if ( loss > kinEnergy) loss = kinEnergy; 519 if ( loss > kinEnergy) loss = kinEnergy;
517 else if( loss < 0.) loss = 0.; 520 else if( loss < 0.) loss = 0.;
518 521
519 return loss; 522 return loss;
520 } 523 }
521 524
522 ////////////////////////////////////////////// 525 ////////////////////////////////////////////////////////////////////////
523 526
524 G4double G4PAIPhotData::SampleAlongStepPhotonT 527 G4double G4PAIPhotData::SampleAlongStepPhotonTransfer(G4int coupleIndex,
525 528 G4double kinEnergy,
526 G4double scaledTkin, 529 G4double scaledTkin,
527 G4double stepFactor) const 530 G4double stepFactor) const
528 { 531 {
529 G4double loss = 0.0; 532 G4double loss = 0.0;
530 G4double omega; 533 G4double omega;
531 G4double position, E1, E2, W1, W2, W, dNdxCu 534 G4double position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber;
532 535
533 std::size_t iPlace = fParticleEnergyVector-> 536 std::size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
534 std::size_t nPlace = fParticleEnergyVector-> 537 std::size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
535 538
536 G4bool one = true; 539 G4bool one = true;
537 540
538 if (scaledTkin >= fParticleEnergyVector- 541 if (scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
539 else if(scaledTkin > fParticleEnergyVector-> 542 else if(scaledTkin > fParticleEnergyVector->Energy(0)) one = false;
540 543
541 G4PhysicsLogVector* vcut = fdNdxCutPhotonTab 544 G4PhysicsLogVector* vcut = fdNdxCutPhotonTable[coupleIndex];
542 G4PhysicsVector* v1 = (*(fPAIphotonBank 545 G4PhysicsVector* v1 = (*(fPAIphotonBank[coupleIndex]))(iPlace);
543 G4PhysicsVector* v2 = nullptr; 546 G4PhysicsVector* v2 = nullptr;
544 547
545 dNdxCut1 = (*vcut)[iPlace]; 548 dNdxCut1 = (*vcut)[iPlace];
546 G4double e1 = v1->Energy(0); 549 G4double e1 = v1->Energy(0);
547 G4double e2 = e1; 550 G4double e2 = e1;
548 551
549 G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*s 552 G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;
550 553
551 meanNumber = meanN1; 554 meanNumber = meanN1;
552 555
553 // G4cout<<"iPlace = "<<iPlace<< " meanN1= " 556 // G4cout<<"iPlace = "<<iPlace<< " meanN1= " << meanN1
554 // <<" (*v1)[0]/e1 = "<<(*v1)[0]/e1<< " 557 // <<" (*v1)[0]/e1 = "<<(*v1)[0]/e1<< " dNdxCut1= " << dNdxCut1 << G4endl;
555 558
556 dNdxCut2 = dNdxCut1; 559 dNdxCut2 = dNdxCut1;
557 W1 = 1.0; 560 W1 = 1.0;
558 W2 = 0.0; 561 W2 = 0.0;
559 if(!one) 562 if(!one)
560 { 563 {
561 v2 = (*(fPAIphotonBank[coupleIndex]))(iPla 564 v2 = (*(fPAIphotonBank[coupleIndex]))(iPlace+1);
562 dNdxCut2 = (*vcut)[iPlace+1]; 565 dNdxCut2 = (*vcut)[iPlace+1];
563 e2 = v2->Energy(0); 566 e2 = v2->Energy(0);
564 567
565 G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2) 568 G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;
566 569
567 E1 = fParticleEnergyVector->Energy(iPlace) 570 E1 = fParticleEnergyVector->Energy(iPlace);
568 E2 = fParticleEnergyVector->Energy(iPlace+ 571 E2 = fParticleEnergyVector->Energy(iPlace+1);
569 W = 1.0/(E2 - E1); 572 W = 1.0/(E2 - E1);
570 W1 = (E2 - scaledTkin)*W; 573 W1 = (E2 - scaledTkin)*W;
571 W2 = (scaledTkin - E1)*W; 574 W2 = (scaledTkin - E1)*W;
572 meanNumber = W1*meanN1 + W2*meanN2; 575 meanNumber = W1*meanN1 + W2*meanN2;
573 576
574 //G4cout<<"meanN= " << meanNumber << " me 577 //G4cout<<"meanN= " << meanNumber << " meanN2= " << meanN2
575 // << " dNdxCut2= " << dNdxCut2 << G4en 578 // << " dNdxCut2= " << dNdxCut2 << G4endl;
576 } 579 }
577 if( meanNumber <= 0.0) return 0.0; 580 if( meanNumber <= 0.0) return 0.0;
578 581
579 G4int numOfCollisions = (G4int)G4Poisson(mea 582 G4int numOfCollisions = (G4int)G4Poisson(meanNumber);
580 583
581 //G4cout << "N= " << numOfCollisions << G4en 584 //G4cout << "N= " << numOfCollisions << G4endl;
582 585
583 if( 0 == numOfCollisions) return 0.0; 586 if( 0 == numOfCollisions) return 0.0;
584 587
585 for(G4int i=0; i< numOfCollisions; ++i) 588 for(G4int i=0; i< numOfCollisions; ++i)
586 { 589 {
587 G4double rand = G4UniformRand(); 590 G4double rand = G4UniformRand();
588 position = dNdxCut1 + ((*v1)[0]/e1 - dNdxC 591 position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
589 omega = GetEnergyPhotonTransfer(coupleInde 592 omega = GetEnergyPhotonTransfer(coupleIndex, iPlace, position);
590 593
591 //G4cout << "omega(keV)= " << omega/keV << 594 //G4cout << "omega(keV)= " << omega/keV << G4endl;
592 595
593 if(!one) 596 if(!one)
594 { 597 {
595 position = dNdxCut2 + ((*v2)[0]/e2 - dNd 598 position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
596 G4double omega2 = GetEnergyPhotonTransfe 599 G4double omega2 = GetEnergyPhotonTransfer(coupleIndex, iPlace+1, position);
597 omega = omega*W1 + omega2*W2; 600 omega = omega*W1 + omega2*W2;
598 } 601 }
599 //G4cout << "omega(keV)= " << omega/keV << 602 //G4cout << "omega(keV)= " << omega/keV << G4endl;
600 603
601 loss += omega; 604 loss += omega;
602 if( loss > kinEnergy) { break; } 605 if( loss > kinEnergy) { break; }
603 } 606 }
604 607
605 // G4cout<<"PAIPhotData AlongStepLoss = "<<l 608 // G4cout<<"PAIPhotData AlongStepLoss = "<<loss/keV<<" keV, on step = "
606 //<<step/mm<<" mm"<<G4endl; 609 //<<step/mm<<" mm"<<G4endl;
607 610
608 if ( loss > kinEnergy) loss = kinEnergy; 611 if ( loss > kinEnergy) loss = kinEnergy;
609 else if( loss < 0.) loss = 0.; 612 else if( loss < 0.) loss = 0.;
610 613
611 return loss; 614 return loss;
612 } 615 }
613 616
614 ////////////////////////////////////////////// 617 //////////////////////////////////////////////////////////////////
615 618
616 G4double G4PAIPhotData::SampleAlongStepPlasmon 619 G4double G4PAIPhotData::SampleAlongStepPlasmonTransfer(G4int coupleIndex,
617 620 G4double kinEnergy,
618 G4double scaledTkin, 621 G4double scaledTkin,
619 G4double stepFactor) const 622 G4double stepFactor) const
620 { 623 {
621 G4double loss = 0.0; 624 G4double loss = 0.0;
622 G4double omega; 625 G4double omega;
623 G4double position, E1, E2, W1, W2, W, dNdxCu 626 G4double position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber;
624 627
625 std::size_t iPlace = fParticleEnergyVector-> 628 std::size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
626 std::size_t nPlace = fParticleEnergyVector-> 629 std::size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
627 630
628 G4bool one = true; 631 G4bool one = true;
629 632
630 if (scaledTkin >= fParticleEnergyVector- 633 if (scaledTkin >= fParticleEnergyVector->Energy(nPlace)) iPlace = nPlace;
631 else if(scaledTkin > fParticleEnergyVector-> 634 else if(scaledTkin > fParticleEnergyVector->Energy(0)) one = false;
632 635
633 G4PhysicsLogVector* vcut = fdNdxCutPlasmonTa 636 G4PhysicsLogVector* vcut = fdNdxCutPlasmonTable[coupleIndex];
634 G4PhysicsVector* v1 = (*(fPAIplasmonBan 637 G4PhysicsVector* v1 = (*(fPAIplasmonBank[coupleIndex]))(iPlace);
635 G4PhysicsVector* v2 = nullptr; 638 G4PhysicsVector* v2 = nullptr;
636 639
637 dNdxCut1 = (*vcut)[iPlace]; 640 dNdxCut1 = (*vcut)[iPlace];
638 G4double e1 = v1->Energy(0); 641 G4double e1 = v1->Energy(0);
639 G4double e2 = e1; 642 G4double e2 = e1;
640 643
641 G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*s 644 G4double meanN1 = ((*v1)[0]/e1 - dNdxCut1)*stepFactor;
642 645
643 meanNumber = meanN1; 646 meanNumber = meanN1;
644 647
645 // G4cout<<"iPlace = "<<iPlace<< " meanN1= " 648 // G4cout<<"iPlace = "<<iPlace<< " meanN1= " << meanN1
646 // <<" (*v1)[0]/e1 = "<<(*v1)[0]/e1<< " 649 // <<" (*v1)[0]/e1 = "<<(*v1)[0]/e1<< " dNdxCut1= " << dNdxCut1 << G4endl;
647 650
648 dNdxCut2 = dNdxCut1; 651 dNdxCut2 = dNdxCut1;
649 W1 = 1.0; 652 W1 = 1.0;
650 W2 = 0.0; 653 W2 = 0.0;
651 if(!one) 654 if(!one)
652 { 655 {
653 v2 = (*(fPAIplasmonBank[coupleIndex]))(iPl 656 v2 = (*(fPAIplasmonBank[coupleIndex]))(iPlace+1);
654 dNdxCut2 = (*vcut)[iPlace+1]; 657 dNdxCut2 = (*vcut)[iPlace+1];
655 e2 = v2->Energy(0); 658 e2 = v2->Energy(0);
656 659
657 G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2) 660 G4double meanN2 = ((*v2)[0]/e2 - dNdxCut2)*stepFactor;
658 661
659 E1 = fParticleEnergyVector->Energy(iPlace) 662 E1 = fParticleEnergyVector->Energy(iPlace);
660 E2 = fParticleEnergyVector->Energy(iPlace+ 663 E2 = fParticleEnergyVector->Energy(iPlace+1);
661 W = 1.0/(E2 - E1); 664 W = 1.0/(E2 - E1);
662 W1 = (E2 - scaledTkin)*W; 665 W1 = (E2 - scaledTkin)*W;
663 W2 = (scaledTkin - E1)*W; 666 W2 = (scaledTkin - E1)*W;
664 meanNumber = W1*meanN1 + W2*meanN2; 667 meanNumber = W1*meanN1 + W2*meanN2;
665 668
666 //G4cout<<"meanN= " << meanNumber << " me 669 //G4cout<<"meanN= " << meanNumber << " meanN2= " << meanN2
667 // << " dNdxCut2= " << dNdxCut2 << G4en 670 // << " dNdxCut2= " << dNdxCut2 << G4endl;
668 } 671 }
669 if( meanNumber <= 0.0) return 0.0; 672 if( meanNumber <= 0.0) return 0.0;
670 673
671 G4int numOfCollisions = (G4int)G4Poisson(mea 674 G4int numOfCollisions = (G4int)G4Poisson(meanNumber);
672 675
673 //G4cout << "N= " << numOfCollisions << G4en 676 //G4cout << "N= " << numOfCollisions << G4endl;
674 677
675 if( 0 == numOfCollisions) return 0.0; 678 if( 0 == numOfCollisions) return 0.0;
676 679
677 for(G4int i=0; i< numOfCollisions; ++i) 680 for(G4int i=0; i< numOfCollisions; ++i)
678 { 681 {
679 G4double rand = G4UniformRand(); 682 G4double rand = G4UniformRand();
680 position = dNdxCut1 + ((*v1)[0]/e1 - dNdxC 683 position = dNdxCut1 + ((*v1)[0]/e1 - dNdxCut1)*rand;
681 omega = GetEnergyPlasmonTransfer(coupleInd 684 omega = GetEnergyPlasmonTransfer(coupleIndex, iPlace, position);
682 685
683 //G4cout << "omega(keV)= " << omega/keV << 686 //G4cout << "omega(keV)= " << omega/keV << G4endl;
684 687
685 if(!one) 688 if(!one)
686 { 689 {
687 position = dNdxCut2 + ((*v2)[0]/e2 - dNd 690 position = dNdxCut2 + ((*v2)[0]/e2 - dNdxCut2)*rand;
688 G4double omega2 = GetEnergyPlasmonTransf 691 G4double omega2 = GetEnergyPlasmonTransfer(coupleIndex, iPlace+1, position);
689 omega = omega*W1 + omega2*W2; 692 omega = omega*W1 + omega2*W2;
690 } 693 }
691 //G4cout << "omega(keV)= " << omega/keV << 694 //G4cout << "omega(keV)= " << omega/keV << G4endl;
692 695
693 loss += omega; 696 loss += omega;
694 if( loss > kinEnergy) { break; } 697 if( loss > kinEnergy) { break; }
695 } 698 }
696 699
697 // G4cout<<"PAIPhotData AlongStepLoss = "<<l 700 // G4cout<<"PAIPhotData AlongStepLoss = "<<loss/keV<<" keV, on step = "
698 //<<step/mm<<" mm"<<G4endl; 701 //<<step/mm<<" mm"<<G4endl;
699 702
700 if ( loss > kinEnergy) loss = kinEnergy; 703 if ( loss > kinEnergy) loss = kinEnergy;
701 else if( loss < 0.) loss = 0.; 704 else if( loss < 0.) loss = 0.;
702 705
703 return loss; 706 return loss;
704 } 707 }
705 708
706 ////////////////////////////////////////////// 709 ///////////////////////////////////////////////////////////////////////
707 // 710 //
708 // Returns post step PAI energy transfer > cut 711 // Returns post step PAI energy transfer > cut electron energy
709 // according to passed scaled kinetic energy o 712 // according to passed scaled kinetic energy of particle
710 713
711 G4double G4PAIPhotData::SamplePostStepTransfer 714 G4double G4PAIPhotData::SamplePostStepTransfer(G4int coupleIndex,
712 G4double scaledTkin) const 715 G4double scaledTkin) const
713 { 716 {
714 //G4cout<<"G4PAIPhotData::GetPostStepTransfe 717 //G4cout<<"G4PAIPhotData::GetPostStepTransfer"<<G4endl;
715 G4double transfer = 0.0; 718 G4double transfer = 0.0;
716 G4double rand = G4UniformRand(); 719 G4double rand = G4UniformRand();
717 720
718 std::size_t nPlace = fParticleEnergyVector-> 721 std::size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
719 722
720 // std::size_t iTransfer, iTr1, iTr2; 723 // std::size_t iTransfer, iTr1, iTr2;
721 G4double position, dNdxCut1, dNdxCut2, E1, E 724 G4double position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W;
722 725
723 G4PhysicsVector* cutv = fdNdxCutTable[couple 726 G4PhysicsVector* cutv = fdNdxCutTable[coupleIndex];
724 727
725 // Fermi plato, try from left 728 // Fermi plato, try from left
726 if( scaledTkin >= fParticleEnergyVector->Get 729 if( scaledTkin >= fParticleEnergyVector->GetMaxEnergy())
727 { 730 {
728 position = (*cutv)[nPlace]*rand; 731 position = (*cutv)[nPlace]*rand;
729 transfer = GetEnergyTransfer(coupleIndex, 732 transfer = GetEnergyTransfer(coupleIndex, nPlace, position);
730 } 733 }
731 else if( scaledTkin <= fParticleEnergyVector 734 else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
732 { 735 {
733 position = (*cutv)[0]*rand; 736 position = (*cutv)[0]*rand;
734 transfer = GetEnergyTransfer(coupleIndex, 737 transfer = GetEnergyTransfer(coupleIndex, 0, position);
735 } 738 }
736 else 739 else
737 { 740 {
738 std::size_t iPlace = fParticleEnergyVector 741 std::size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
739 742
740 dNdxCut1 = (*cutv)[iPlace]; 743 dNdxCut1 = (*cutv)[iPlace];
741 dNdxCut2 = (*cutv)[iPlace+1]; 744 dNdxCut2 = (*cutv)[iPlace+1];
742 745
743 E1 = fParticleEnergyVector->Energy(iPlace) 746 E1 = fParticleEnergyVector->Energy(iPlace);
744 E2 = fParticleEnergyVector->Energy(iPlace+ 747 E2 = fParticleEnergyVector->Energy(iPlace+1);
745 W = 1.0/(E2 - E1); 748 W = 1.0/(E2 - E1);
746 W1 = (E2 - scaledTkin)*W; 749 W1 = (E2 - scaledTkin)*W;
747 W2 = (scaledTkin - E1)*W; 750 W2 = (scaledTkin - E1)*W;
748 751
749 //G4cout<<"iPlace= " << " dNdxCut1 = "<<d 752 //G4cout<<"iPlace= " << " dNdxCut1 = "<<dNdxCut1
750 // <<" dNdxCut2 = "<<dNdxCut2<< " W1= " 753 // <<" dNdxCut2 = "<<dNdxCut2<< " W1= " << W1 << " W2= " << W2 <<G4endl;
751 754
752 position = dNdxCut1*rand; 755 position = dNdxCut1*rand;
753 G4double tr1 = GetEnergyTransfer(coupleInd 756 G4double tr1 = GetEnergyTransfer(coupleIndex, iPlace, position);
754 757
755 position = dNdxCut2*rand; 758 position = dNdxCut2*rand;
756 G4double tr2 = GetEnergyTransfer(coupleInd 759 G4double tr2 = GetEnergyTransfer(coupleIndex, iPlace+1, position);
757 760
758 transfer = tr1*W1 + tr2*W2; 761 transfer = tr1*W1 + tr2*W2;
759 } 762 }
760 //G4cout<<"PAImodel PostStepTransfer = "<<tr 763 //G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"<<G4endl;
761 if(transfer < 0.0 ) { transfer = 0.0; } 764 if(transfer < 0.0 ) { transfer = 0.0; }
762 return transfer; 765 return transfer;
763 } 766 }
764 767
765 ////////////////////////////////////////////// 768 ////////////////////////////////////////////////////////////////////////
766 769
767 G4double G4PAIPhotData::SamplePostStepPhotonTr 770 G4double G4PAIPhotData::SamplePostStepPhotonTransfer(G4int coupleIndex,
768 G4double scaledTkin) const 771 G4double scaledTkin) const
769 { 772 {
770 //G4cout<<"G4PAIPhotData::GetPostStepTransfe 773 //G4cout<<"G4PAIPhotData::GetPostStepTransfer"<<G4endl;
771 G4double transfer = 0.0; 774 G4double transfer = 0.0;
772 G4double rand = G4UniformRand(); 775 G4double rand = G4UniformRand();
773 776
774 std::size_t nPlace = fParticleEnergyVector-> 777 std::size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
775 778
776 // std::size_t iTransfer, iTr1, iTr2; 779 // std::size_t iTransfer, iTr1, iTr2;
777 G4double position, dNdxCut1, dNdxCut2, E1, E 780 G4double position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W;
778 781
779 G4PhysicsVector* cutv = fdNdxCutPhotonTable[ 782 G4PhysicsVector* cutv = fdNdxCutPhotonTable[coupleIndex];
780 783
781 // Fermi plato, try from left 784 // Fermi plato, try from left
782 785
783 if( scaledTkin >= fParticleEnergyVector->Get 786 if( scaledTkin >= fParticleEnergyVector->GetMaxEnergy())
784 { 787 {
785 position = (*cutv)[nPlace]*rand; 788 position = (*cutv)[nPlace]*rand;
786 transfer = GetEnergyPhotonTransfer(coupleI 789 transfer = GetEnergyPhotonTransfer(coupleIndex, nPlace, position);
787 } 790 }
788 else if( scaledTkin <= fParticleEnergyVector 791 else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
789 { 792 {
790 position = (*cutv)[0]*rand; 793 position = (*cutv)[0]*rand;
791 transfer = GetEnergyPhotonTransfer(coupleI 794 transfer = GetEnergyPhotonTransfer(coupleIndex, 0, position);
792 } 795 }
793 else 796 else
794 { 797 {
795 std::size_t iPlace = fParticleEnergyVector 798 std::size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
796 799
797 dNdxCut1 = (*cutv)[iPlace]; 800 dNdxCut1 = (*cutv)[iPlace];
798 dNdxCut2 = (*cutv)[iPlace+1]; 801 dNdxCut2 = (*cutv)[iPlace+1];
799 802
800 E1 = fParticleEnergyVector->Energy(iPlace) 803 E1 = fParticleEnergyVector->Energy(iPlace);
801 E2 = fParticleEnergyVector->Energy(iPlace+ 804 E2 = fParticleEnergyVector->Energy(iPlace+1);
802 W = 1.0/(E2 - E1); 805 W = 1.0/(E2 - E1);
803 W1 = (E2 - scaledTkin)*W; 806 W1 = (E2 - scaledTkin)*W;
804 W2 = (scaledTkin - E1)*W; 807 W2 = (scaledTkin - E1)*W;
805 808
806 //G4cout<<"iPlace= " << " dNdxCut1 = "<<d 809 //G4cout<<"iPlace= " << " dNdxCut1 = "<<dNdxCut1
807 // <<" dNdxCut2 = "<<dNdxCut2<< " W1= " 810 // <<" dNdxCut2 = "<<dNdxCut2<< " W1= " << W1 << " W2= " << W2 <<G4endl;
808 811
809 position = dNdxCut1*rand; 812 position = dNdxCut1*rand;
810 813
811 G4double tr1 = GetEnergyPhotonTransfer(cou 814 G4double tr1 = GetEnergyPhotonTransfer(coupleIndex, iPlace, position);
812 815
813 position = dNdxCut2*rand; 816 position = dNdxCut2*rand;
814 G4double tr2 = GetEnergyPhotonTransfer(cou 817 G4double tr2 = GetEnergyPhotonTransfer(coupleIndex, iPlace+1, position);
815 818
816 transfer = tr1*W1 + tr2*W2; 819 transfer = tr1*W1 + tr2*W2;
817 } 820 }
818 //G4cout<<"PAImodel PostStepTransfer = "<<tr 821 //G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"<<G4endl;
819 if(transfer < 0.0 ) { transfer = 0.0; } 822 if(transfer < 0.0 ) { transfer = 0.0; }
820 return transfer; 823 return transfer;
821 } 824 }
822 825
823 ////////////////////////////////////////////// 826 //////////////////////////////////////////////////////////////////////////
824 827
825 G4double G4PAIPhotData::SamplePostStepPlasmonT 828 G4double G4PAIPhotData::SamplePostStepPlasmonTransfer(G4int coupleIndex,
826 G4double scaledTkin) const 829 G4double scaledTkin) const
827 { 830 {
828 //G4cout<<"G4PAIPhotData::GetPostStepTransfe 831 //G4cout<<"G4PAIPhotData::GetPostStepTransfer"<<G4endl;
829 G4double transfer = 0.0; 832 G4double transfer = 0.0;
830 G4double rand = G4UniformRand(); 833 G4double rand = G4UniformRand();
831 834
832 std::size_t nPlace = fParticleEnergyVector-> 835 std::size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1;
833 836
834 // std::size_t iTransfer, iTr1, iTr2; 837 // std::size_t iTransfer, iTr1, iTr2;
835 G4double position, dNdxCut1, dNdxCut2, E1, E 838 G4double position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W;
836 839
837 G4PhysicsVector* cutv = fdNdxCutPlasmonTable 840 G4PhysicsVector* cutv = fdNdxCutPlasmonTable[coupleIndex];
838 841
839 // Fermi plato, try from left 842 // Fermi plato, try from left
840 if( scaledTkin >= fParticleEnergyVector->Get 843 if( scaledTkin >= fParticleEnergyVector->GetMaxEnergy())
841 { 844 {
842 position = (*cutv)[nPlace]*rand; 845 position = (*cutv)[nPlace]*rand;
843 transfer = GetEnergyPlasmonTransfer(couple 846 transfer = GetEnergyPlasmonTransfer(coupleIndex, nPlace, position);
844 } 847 }
845 else if( scaledTkin <= fParticleEnergyVector 848 else if( scaledTkin <= fParticleEnergyVector->Energy(0) )
846 { 849 {
847 position = (*cutv)[0]*rand; 850 position = (*cutv)[0]*rand;
848 transfer = GetEnergyPlasmonTransfer(couple 851 transfer = GetEnergyPlasmonTransfer(coupleIndex, 0, position);
849 } 852 }
850 else 853 else
851 { 854 {
852 std::size_t iPlace = fParticleEnergyVector 855 std::size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0);
853 856
854 dNdxCut1 = (*cutv)[iPlace]; 857 dNdxCut1 = (*cutv)[iPlace];
855 dNdxCut2 = (*cutv)[iPlace+1]; 858 dNdxCut2 = (*cutv)[iPlace+1];
856 859
857 E1 = fParticleEnergyVector->Energy(iPlace) 860 E1 = fParticleEnergyVector->Energy(iPlace);
858 E2 = fParticleEnergyVector->Energy(iPlace+ 861 E2 = fParticleEnergyVector->Energy(iPlace+1);
859 W = 1.0/(E2 - E1); 862 W = 1.0/(E2 - E1);
860 W1 = (E2 - scaledTkin)*W; 863 W1 = (E2 - scaledTkin)*W;
861 W2 = (scaledTkin - E1)*W; 864 W2 = (scaledTkin - E1)*W;
862 865
863 //G4cout<<"iPlace= " << " dNdxCut1 = "<<d 866 //G4cout<<"iPlace= " << " dNdxCut1 = "<<dNdxCut1
864 // <<" dNdxCut2 = "<<dNdxCut2<< " W1= " 867 // <<" dNdxCut2 = "<<dNdxCut2<< " W1= " << W1 << " W2= " << W2 <<G4endl;
865 868
866 position = dNdxCut1*rand; 869 position = dNdxCut1*rand;
867 G4double tr1 = GetEnergyPlasmonTransfer(co 870 G4double tr1 = GetEnergyPlasmonTransfer(coupleIndex, iPlace, position);
868 871
869 position = dNdxCut2*rand; 872 position = dNdxCut2*rand;
870 G4double tr2 = GetEnergyPlasmonTransfer(co 873 G4double tr2 = GetEnergyPlasmonTransfer(coupleIndex, iPlace+1, position);
871 874
872 transfer = tr1*W1 + tr2*W2; 875 transfer = tr1*W1 + tr2*W2;
873 } 876 }
874 //G4cout<<"PAImodel PostStepPlasmonTransfer 877 //G4cout<<"PAImodel PostStepPlasmonTransfer = "<<transfer/keV<<" keV"<<G4endl;
875 878
876 if(transfer < 0.0 ) transfer = 0.0; 879 if(transfer < 0.0 ) transfer = 0.0;
877 880
878 return transfer; 881 return transfer;
879 } 882 }
880 883
881 ////////////////////////////////////////////// 884 ///////////////////////////////////////////////////////////////////////
882 // 885 //
883 // Returns PAI energy transfer according to pa 886 // Returns PAI energy transfer according to passed
884 // indexes of particle kinetic enegry and rand 887 // indexes of particle kinetic enegry and random x-section
885 888
886 G4double G4PAIPhotData::GetEnergyTransfer(G4in 889 G4double G4PAIPhotData::GetEnergyTransfer(G4int coupleIndex,
887 std::size_t iPlace, 890 std::size_t iPlace,
888 G4double position) const 891 G4double position) const
889 { 892 {
890 G4PhysicsVector* v = (*(fPAIxscBank[coupleIn 893 G4PhysicsVector* v = (*(fPAIxscBank[coupleIndex]))(iPlace);
891 if(position*v->Energy(0) >= (*v)[0]) { retur 894 if(position*v->Energy(0) >= (*v)[0]) { return v->Energy(0); }
892 895
893 std::size_t iTransferMax = v->GetVectorLengt 896 std::size_t iTransferMax = v->GetVectorLength() - 1;
894 897
895 std::size_t iTransfer; 898 std::size_t iTransfer;
896 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), 899 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
897 900
898 for(iTransfer=1; iTransfer<=iTransferMax; ++ 901 for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer) {
899 x2 = v->Energy(iTransfer); 902 x2 = v->Energy(iTransfer);
900 y2 = (*v)[iTransfer]/x2; 903 y2 = (*v)[iTransfer]/x2;
901 if(position >= y2) { break; } 904 if(position >= y2) { break; }
902 } 905 }
903 906
904 x1 = v->Energy(iTransfer-1); 907 x1 = v->Energy(iTransfer-1);
905 y1 = (*v)[iTransfer-1]/x1; 908 y1 = (*v)[iTransfer-1]/x1;
906 //G4cout << "i= " << iTransfer << " imax= " 909 //G4cout << "i= " << iTransfer << " imax= " << iTransferMax
907 // << " x1= " << x1 << " x2= " << x2 << G4 910 // << " x1= " << x1 << " x2= " << x2 << G4endl;
908 911
909 energyTransfer = x1; 912 energyTransfer = x1;
910 if ( x1 != x2 ) { 913 if ( x1 != x2 ) {
911 if ( y1 == y2 ) { 914 if ( y1 == y2 ) {
912 energyTransfer += (x2 - x1)*G4UniformRan 915 energyTransfer += (x2 - x1)*G4UniformRand();
913 } else { 916 } else {
914 if(x1*1.1 < x2) { 917 if(x1*1.1 < x2) {
915 const G4int nbins = 5; 918 const G4int nbins = 5;
916 G4double del = (x2 - x1)/G4int(nbins); 919 G4double del = (x2 - x1)/G4int(nbins);
917 x2 = x1; 920 x2 = x1;
918 for(G4int i=1; i<=nbins; ++i) { 921 for(G4int i=1; i<=nbins; ++i) {
919 x2 += del; 922 x2 += del;
920 y2 = v->Value(x2)/x2; 923 y2 = v->Value(x2)/x2;
921 if(position >= y2) { break; } 924 if(position >= y2) { break; }
922 x1 = x2; 925 x1 = x2;
923 y1 = y2; 926 y1 = y2;
924 } 927 }
925 } 928 }
926 energyTransfer = (y2 - y1)*x1*x2/(positi 929 energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2);
927 } 930 }
928 } 931 }
929 // G4cout << "x1(keV)= " << x1/keV << " x2( 932 // G4cout << "x1(keV)= " << x1/keV << " x2(keV)= " << x2/keV
930 // << " y1= " << y1 << " y2= " << y2 << " 933 // << " y1= " << y1 << " y2= " << y2 << " pos= " << position
931 // << " E(keV)= " << energyTransfer/keV << 934 // << " E(keV)= " << energyTransfer/keV << G4endl;
932 return energyTransfer; 935 return energyTransfer;
933 } 936 }
934 937
935 ////////////////////////////////////////////// 938 /////////////////////////////////////////////////////////////////
936 939
937 G4double G4PAIPhotData::GetEnergyPhotonTransfe 940 G4double G4PAIPhotData::GetEnergyPhotonTransfer(G4int coupleIndex,
938 std::size_t iPlace, 941 std::size_t iPlace,
939 G4double position) const 942 G4double position) const
940 { 943 {
941 G4PhysicsVector* v = (*(fPAIphotonBank[coupl 944 G4PhysicsVector* v = (*(fPAIphotonBank[coupleIndex]))(iPlace);
942 if(position*v->Energy(0) >= (*v)[0]) return 945 if(position*v->Energy(0) >= (*v)[0]) return v->Energy(0);
943 946
944 std::size_t iTransferMax = v->GetVectorLengt 947 std::size_t iTransferMax = v->GetVectorLength() - 1;
945 948
946 std::size_t iTransfer; 949 std::size_t iTransfer;
947 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), 950 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
948 951
949 for(iTransfer=1; iTransfer<=iTransferMax; ++ 952 for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer)
950 { 953 {
951 x2 = v->Energy(iTransfer); 954 x2 = v->Energy(iTransfer);
952 y2 = (*v)[iTransfer]/x2; 955 y2 = (*v)[iTransfer]/x2;
953 if(position >= y2) break; 956 if(position >= y2) break;
954 } 957 }
955 x1 = v->Energy(iTransfer-1); 958 x1 = v->Energy(iTransfer-1);
956 y1 = (*v)[iTransfer-1]/x1; 959 y1 = (*v)[iTransfer-1]/x1;
957 960
958 //G4cout << "i= " << iTransfer << " imax= " 961 //G4cout << "i= " << iTransfer << " imax= " << iTransferMax
959 // << " x1= " << x1 << " x2= " << x2 << G4 962 // << " x1= " << x1 << " x2= " << x2 << G4endl;
960 963
961 energyTransfer = x1; 964 energyTransfer = x1;
962 965
963 if ( x1 != x2 ) 966 if ( x1 != x2 )
964 { 967 {
965 if ( y1 == y2 ) 968 if ( y1 == y2 )
966 { 969 {
967 energyTransfer += (x2 - x1)*G4UniformRan 970 energyTransfer += (x2 - x1)*G4UniformRand();
968 } 971 }
969 else 972 else
970 { 973 {
971 if( x1*1.1 < x2 ) 974 if( x1*1.1 < x2 )
972 { 975 {
973 const G4int nbins = 5; 976 const G4int nbins = 5;
974 G4double del = (x2 - x1)/G4int(nbins); 977 G4double del = (x2 - x1)/G4int(nbins);
975 x2 = x1; 978 x2 = x1;
976 979
977 for(G4int i=1; i<=nbins; ++i) 980 for(G4int i=1; i<=nbins; ++i)
978 { 981 {
979 x2 += del; 982 x2 += del;
980 y2 = v->Value(x2)/x2; 983 y2 = v->Value(x2)/x2;
981 if(position >= y2) { break; } 984 if(position >= y2) { break; }
982 x1 = x2; 985 x1 = x2;
983 y1 = y2; 986 y1 = y2;
984 } 987 }
985 } 988 }
986 energyTransfer = (y2 - y1)*x1*x2/(positi 989 energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2);
987 } 990 }
988 } 991 }
989 // G4cout << "x1(keV)= " << x1/keV << " x2( 992 // G4cout << "x1(keV)= " << x1/keV << " x2(keV)= " << x2/keV
990 // << " y1= " << y1 << " y2= " << y2 << " 993 // << " y1= " << y1 << " y2= " << y2 << " pos= " << position
991 // << " E(keV)= " << energyTransfer/keV << 994 // << " E(keV)= " << energyTransfer/keV << G4endl;
992 995
993 return energyTransfer; 996 return energyTransfer;
994 } 997 }
995 998
996 ////////////////////////////////////////////// 999 /////////////////////////////////////////////////////////////////////////
997 1000
998 G4double G4PAIPhotData::GetEnergyPlasmonTransf 1001 G4double G4PAIPhotData::GetEnergyPlasmonTransfer(G4int coupleIndex,
999 std::size_t iPlace, 1002 std::size_t iPlace,
1000 G4double position) const 1003 G4double position) const
1001 { 1004 {
1002 G4PhysicsVector* v = (*(fPAIplasmonBank[cou 1005 G4PhysicsVector* v = (*(fPAIplasmonBank[coupleIndex]))(iPlace);
1003 1006
1004 if( position*v->Energy(0) >= (*v)[0] ) ret 1007 if( position*v->Energy(0) >= (*v)[0] ) return v->Energy(0);
1005 1008
1006 std::size_t iTransferMax = v->GetVectorLeng 1009 std::size_t iTransferMax = v->GetVectorLength() - 1;
1007 1010
1008 std::size_t iTransfer; 1011 std::size_t iTransfer;
1009 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0) 1012 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer;
1010 1013
1011 for(iTransfer = 1; iTransfer <= iTransferMa 1014 for(iTransfer = 1; iTransfer <= iTransferMax; ++iTransfer)
1012 { 1015 {
1013 x2 = v->Energy(iTransfer); 1016 x2 = v->Energy(iTransfer);
1014 y2 = (*v)[iTransfer]/x2; 1017 y2 = (*v)[iTransfer]/x2;
1015 if(position >= y2) break; 1018 if(position >= y2) break;
1016 } 1019 }
1017 x1 = v->Energy(iTransfer-1); 1020 x1 = v->Energy(iTransfer-1);
1018 y1 = (*v)[iTransfer-1]/x1; 1021 y1 = (*v)[iTransfer-1]/x1;
1019 1022
1020 //G4cout << "i= " << iTransfer << " imax= " 1023 //G4cout << "i= " << iTransfer << " imax= " << iTransferMax
1021 // << " x1= " << x1 << " x2= " << x2 << G 1024 // << " x1= " << x1 << " x2= " << x2 << G4endl;
1022 1025
1023 energyTransfer = x1; 1026 energyTransfer = x1;
1024 1027
1025 if ( x1 != x2 ) 1028 if ( x1 != x2 )
1026 { 1029 {
1027 if ( y1 == y2 ) 1030 if ( y1 == y2 )
1028 { 1031 {
1029 energyTransfer += (x2 - x1)*G4UniformRa 1032 energyTransfer += (x2 - x1)*G4UniformRand();
1030 } 1033 }
1031 else 1034 else
1032 { 1035 {
1033 if(x1*1.1 < x2) 1036 if(x1*1.1 < x2)
1034 { 1037 {
1035 const G4int nbins = 5; 1038 const G4int nbins = 5;
1036 G4double del = (x2 - x1)/G4int(nbins) 1039 G4double del = (x2 - x1)/G4int(nbins);
1037 x2 = x1; 1040 x2 = x1;
1038 1041
1039 for( G4int i = 1; i <= nbins; ++i ) 1042 for( G4int i = 1; i <= nbins; ++i )
1040 { 1043 {
1041 x2 += del; 1044 x2 += del;
1042 y2 = v->Value(x2)/x2; 1045 y2 = v->Value(x2)/x2;
1043 1046
1044 if(position >= y2) break; 1047 if(position >= y2) break;
1045 1048
1046 x1 = x2; 1049 x1 = x2;
1047 y1 = y2; 1050 y1 = y2;
1048 } 1051 }
1049 } 1052 }
1050 energyTransfer = (y2 - y1)*x1*x2/(posit 1053 energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2);
1051 } 1054 }
1052 } 1055 }
1053 // G4cout << "x1(keV)= " << x1/keV << " x2 1056 // G4cout << "x1(keV)= " << x1/keV << " x2(keV)= " << x2/keV
1054 // << " y1= " << y1 << " y2= " << y2 << " 1057 // << " y1= " << y1 << " y2= " << y2 << " pos= " << position
1055 // << " E(keV)= " << energyTransfer/keV < 1058 // << " E(keV)= " << energyTransfer/keV << G4endl;
1056 1059
1057 return energyTransfer; 1060 return energyTransfer;
1058 } 1061 }
1059 1062
1060 // 1063 //
1061 // 1064 //
1062 ///////////////////////////////////////////// 1065 //////////////////////////////////////////////////////////////////////
1063 1066
1064 1067