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 180 } // end of Tkin loop` 183 } // end of Tkin loop` 181 fPAIxscBank.push_back(PAItransferTable); 184 fPAIxscBank.push_back(PAItransferTable); 182 fPAIdEdxBank.push_back(PAIdEdxTable); 185 fPAIdEdxBank.push_back(PAIdEdxTable); 183 //G4cout << "dEdxMeanVector: " << G4endl; 186 //G4cout << "dEdxMeanVector: " << G4endl; 184 //G4cout << *dEdxMeanVector << G4endl; 187 //G4cout << *dEdxMeanVector << G4endl; 185 fdEdxTable.push_back(dEdxMeanVector); 188 fdEdxTable.push_back(dEdxMeanVector); 186 } 189 } 187 190 188 ////////////////////////////////////////////// 191 ////////////////////////////////////////////////////////////////////////////// 189 192 190 G4double G4PAIModelData::DEDXPerVolume(G4int c 193 G4double G4PAIModelData::DEDXPerVolume(G4int coupleIndex, G4double scaledTkin, 191 G4double cut) const 194 G4double cut) const 192 { 195 { 193 // VI: iPlace is the low edge index of the b 196 // VI: iPlace is the low edge index of the bin 194 // iPlace is in interval from 0 to (N-1) 197 // iPlace is in interval from 0 to (N-1) 195 std::size_t iPlace(0); << 198 size_t iPlace(0); 196 G4double dEdx = fdEdxTable[coupleIndex]->Val 199 G4double dEdx = fdEdxTable[coupleIndex]->Value(scaledTkin, iPlace); 197 std::size_t nPlace = fParticleEnergyVector-> << 200 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1; 198 /* 201 /* 199 G4cout << "G4PAIModelData::DEDXPerVolume: co 202 G4cout << "G4PAIModelData::DEDXPerVolume: coupleIdx= " << coupleIndex 200 << " Tscaled= " << scaledTkin << " cut= " < 203 << " Tscaled= " << scaledTkin << " cut= " << cut 201 << " iPlace= " << iPlace << " nPlace= " << 204 << " iPlace= " << iPlace << " nPlace= " << nPlace << G4endl; 202 */ 205 */ 203 G4bool one = true; 206 G4bool one = true; 204 if(scaledTkin >= fParticleEnergyVector->Ener 207 if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; } 205 else if(scaledTkin > fParticleEnergyVector-> 208 else if(scaledTkin > fParticleEnergyVector->Energy(0)) { 206 one = false; 209 one = false; 207 } 210 } 208 211 209 // VI: apply interpolation of the vector 212 // VI: apply interpolation of the vector 210 G4double del = (*(fPAIdEdxBank[coupleIndex] 213 G4double del = (*(fPAIdEdxBank[coupleIndex]))(iPlace)->Value(cut); 211 //G4cout << "dEdx= " << dEdx << " del= " << 214 //G4cout << "dEdx= " << dEdx << " del= " << del << G4endl; 212 if(!one) { 215 if(!one) { 213 G4double del2 = (*(fPAIdEdxBank[coupleInde 216 G4double del2 = (*(fPAIdEdxBank[coupleIndex]))(iPlace+1)->Value(cut); 214 G4double E1 = fParticleEnergyVector->Energ 217 G4double E1 = fParticleEnergyVector->Energy(iPlace); 215 G4double E2 = fParticleEnergyVector->Energ 218 G4double E2 = fParticleEnergyVector->Energy(iPlace+1); 216 G4double W = 1.0/(E2 - E1); 219 G4double W = 1.0/(E2 - E1); 217 G4double W1 = (E2 - scaledTkin)*W; 220 G4double W1 = (E2 - scaledTkin)*W; 218 G4double W2 = (scaledTkin - E1)*W; 221 G4double W2 = (scaledTkin - E1)*W; 219 del *= W1; 222 del *= W1; 220 del += W2*del2; 223 del += W2*del2; 221 } 224 } 222 dEdx -= del; 225 dEdx -= del; 223 //G4cout << "dEdx= " << dEdx << " del= " << 226 //G4cout << "dEdx= " << dEdx << " del= " << del << G4endl; 224 227 225 dEdx = std::max(dEdx, 0.); 228 dEdx = std::max(dEdx, 0.); 226 return dEdx; 229 return dEdx; 227 } 230 } 228 231 229 ////////////////////////////////////////////// 232 ///////////////////////////////////////////////////////////////////////// 230 233 231 G4double 234 G4double 232 G4PAIModelData::CrossSectionPerVolume(G4int co 235 G4PAIModelData::CrossSectionPerVolume(G4int coupleIndex, 233 G4double scaledTkin, 236 G4double scaledTkin, 234 G4double tcut, G4double tmax) co 237 G4double tcut, G4double tmax) const 235 { 238 { 236 G4double cross, cross1, cross2; 239 G4double cross, cross1, cross2; 237 240 238 // iPlace is in interval from 0 to (N-1) 241 // iPlace is in interval from 0 to (N-1) 239 std::size_t iPlace = fParticleEnergyVector-> << 242 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0); 240 std::size_t nPlace = fParticleEnergyVector-> << 243 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1; 241 244 242 G4bool one = true; 245 G4bool one = true; 243 if(scaledTkin >= fParticleEnergyVector->Ener 246 if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; } 244 else if(scaledTkin > fParticleEnergyVector-> 247 else if(scaledTkin > fParticleEnergyVector->Energy(0)) { 245 one = false; 248 one = false; 246 } 249 } 247 G4PhysicsTable* table = fPAIxscBank[coupleIn 250 G4PhysicsTable* table = fPAIxscBank[coupleIndex]; 248 251 249 //G4cout<<"iPlace = "<<iPlace<<"; tmax = " 252 //G4cout<<"iPlace = "<<iPlace<<"; tmax = " 250 // <<tmax<<"; cutEnergy = "<<cutEnergy<<G4en 253 // <<tmax<<"; cutEnergy = "<<cutEnergy<<G4endl; 251 cross1 = (*table)(iPlace)->Value(tmax)/tmax; 254 cross1 = (*table)(iPlace)->Value(tmax)/tmax; 252 //G4cout<<"cross1 = "<<cross1<<G4endl; 255 //G4cout<<"cross1 = "<<cross1<<G4endl; 253 cross2 = (*table)(iPlace)->Value(tcut)/tcut; 256 cross2 = (*table)(iPlace)->Value(tcut)/tcut; 254 //G4cout<<"cross2 = "<<cross2<<G4endl; 257 //G4cout<<"cross2 = "<<cross2<<G4endl; 255 cross = (cross2-cross1); 258 cross = (cross2-cross1); 256 //G4cout<<"cross = "<<cross<<G4endl; 259 //G4cout<<"cross = "<<cross<<G4endl; 257 if(!one) { 260 if(!one) { 258 cross2 = (*table)(iPlace+1)->Value(tcut)/t 261 cross2 = (*table)(iPlace+1)->Value(tcut)/tcut 259 - (*table)(iPlace+1)->Value(tmax)/tmax; 262 - (*table)(iPlace+1)->Value(tmax)/tmax; 260 263 261 G4double E1 = fParticleEnergyVector->Energ 264 G4double E1 = fParticleEnergyVector->Energy(iPlace); 262 G4double E2 = fParticleEnergyVector->Energ 265 G4double E2 = fParticleEnergyVector->Energy(iPlace+1); 263 G4double W = 1.0/(E2 - E1); 266 G4double W = 1.0/(E2 - E1); 264 G4double W1 = (E2 - scaledTkin)*W; 267 G4double W1 = (E2 - scaledTkin)*W; 265 G4double W2 = (scaledTkin - E1)*W; 268 G4double W2 = (scaledTkin - E1)*W; 266 cross *= W1; 269 cross *= W1; 267 cross += W2*cross2; 270 cross += W2*cross2; 268 } 271 } 269 272 270 cross = std::max(cross, 0.0); 273 cross = std::max(cross, 0.0); 271 return cross; 274 return cross; 272 } 275 } 273 276 274 ////////////////////////////////////////////// 277 /////////////////////////////////////////////////////////////////////// 275 278 276 G4double G4PAIModelData::SampleAlongStepTransf 279 G4double G4PAIModelData::SampleAlongStepTransfer(G4int coupleIndex, 277 280 G4double kinEnergy, 278 G4double scaledTkin, 281 G4double scaledTkin, 279 G4double tmax, 282 G4double tmax, 280 G4double stepFactor) const 283 G4double stepFactor) const 281 { 284 { 282 //G4cout << "=== G4PAIModelData::SampleAlong 285 //G4cout << "=== G4PAIModelData::SampleAlongStepTransfer" << G4endl; 283 G4double loss = 0.0; 286 G4double loss = 0.0; 284 287 285 std::size_t iPlace = fParticleEnergyVector-> << 288 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0); 286 std::size_t nPlace = fParticleEnergyVector-> << 289 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1; 287 290 288 G4bool one = true; 291 G4bool one = true; 289 if(scaledTkin >= fParticleEnergyVector->Ener 292 if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; } 290 else if(scaledTkin > fParticleEnergyVector-> 293 else if(scaledTkin > fParticleEnergyVector->Energy(0)) { 291 one = false; 294 one = false; 292 } 295 } 293 296 294 G4double meanNumber = 0.0; 297 G4double meanNumber = 0.0; 295 G4double meanN11 = 0.0; 298 G4double meanN11 = 0.0; 296 G4double meanN12 = 0.0; 299 G4double meanN12 = 0.0; 297 G4double meanN21 = 0.0; 300 G4double meanN21 = 0.0; 298 G4double meanN22 = 0.0; 301 G4double meanN22 = 0.0; 299 302 300 G4PhysicsVector* v1 = (*(fPAIxscBank[coupleI 303 G4PhysicsVector* v1 = (*(fPAIxscBank[coupleIndex]))(iPlace); 301 G4PhysicsVector* v2 = nullptr; << 304 G4PhysicsVector* v2 = 0; 302 305 303 G4double e1 = v1->Energy(0); 306 G4double e1 = v1->Energy(0); 304 G4double e2 = std::min(tmax, v1->GetMaxEnerg 307 G4double e2 = std::min(tmax, v1->GetMaxEnergy()); 305 308 306 if(e2 >= e1) { 309 if(e2 >= e1) { 307 meanN11 = (*v1)[0]/e1; 310 meanN11 = (*v1)[0]/e1; 308 meanN12 = v1->Value(e2)/e2; 311 meanN12 = v1->Value(e2)/e2; 309 meanNumber = (meanN11 - meanN12)*stepFacto 312 meanNumber = (meanN11 - meanN12)*stepFactor; 310 } 313 } 311 //G4cout<<"iPlace = "<<iPlace<< " meanN11= " 314 //G4cout<<"iPlace = "<<iPlace<< " meanN11= " << meanN11 312 // << " meanN12= " << meanN12 << G4endl; 315 // << " meanN12= " << meanN12 << G4endl; 313 316 314 G4double W1 = 1.0; 317 G4double W1 = 1.0; 315 G4double W2 = 0.0; 318 G4double W2 = 0.0; 316 if(!one) { 319 if(!one) { 317 v2 = (*(fPAIxscBank[coupleIndex]))(iPlace+ 320 v2 = (*(fPAIxscBank[coupleIndex]))(iPlace+1); 318 321 319 e1 = v2->Energy(0); 322 e1 = v2->Energy(0); 320 e2 = std::min(tmax, v2->GetMaxEnergy()); 323 e2 = std::min(tmax, v2->GetMaxEnergy()); 321 if(e2 >= e1) { 324 if(e2 >= e1) { 322 meanN21 = (*v2)[0]/e1; 325 meanN21 = (*v2)[0]/e1; 323 meanN22 = v2->Value(e2)/e2; 326 meanN22 = v2->Value(e2)/e2; 324 G4double E1 = fParticleEnergyVector->Ene 327 G4double E1 = fParticleEnergyVector->Energy(iPlace); 325 G4double E2 = fParticleEnergyVector->Ene 328 G4double E2 = fParticleEnergyVector->Energy(iPlace+1); 326 G4double W = 1.0/(E2 - E1); 329 G4double W = 1.0/(E2 - E1); 327 W1 = (E2 - scaledTkin)*W; 330 W1 = (E2 - scaledTkin)*W; 328 W2 = (scaledTkin - E1)*W; 331 W2 = (scaledTkin - E1)*W; 329 meanNumber *= W1; 332 meanNumber *= W1; 330 meanNumber += (meanN21 - meanN22)*stepFa 333 meanNumber += (meanN21 - meanN22)*stepFactor*W2; 331 } 334 } 332 } 335 } 333 336 334 if(meanNumber < 0.0) { return loss; } 337 if(meanNumber < 0.0) { return loss; } 335 G4int numOfCollisions = (G4int)G4Poisson(mea << 338 G4int numOfCollisions = G4Poisson(meanNumber); 336 339 337 //G4cout << "meanNumber= " << meanNumber << 340 //G4cout << "meanNumber= " << meanNumber << " N= " << numOfCollisions << G4endl; 338 341 339 if(0 == numOfCollisions) { return loss; } 342 if(0 == numOfCollisions) { return loss; } 340 343 341 G4double position, omega, omega2; 344 G4double position, omega, omega2; 342 for(G4int i=0; i< numOfCollisions; ++i) { 345 for(G4int i=0; i< numOfCollisions; ++i) { 343 G4double rand = G4UniformRand(); 346 G4double rand = G4UniformRand(); 344 position = meanN12 + (meanN11 - meanN12)*r 347 position = meanN12 + (meanN11 - meanN12)*rand; 345 omega = GetEnergyTransfer(coupleIndex, iPl 348 omega = GetEnergyTransfer(coupleIndex, iPlace, position); 346 //G4cout << "omega(keV)= " << omega/keV << 349 //G4cout << "omega(keV)= " << omega/keV << G4endl; 347 if(!one) { 350 if(!one) { 348 position = meanN22 + (meanN21 - meanN22) 351 position = meanN22 + (meanN21 - meanN22)*rand; 349 omega2 = GetEnergyTransfer(coupleIndex, 352 omega2 = GetEnergyTransfer(coupleIndex, iPlace+1, position); 350 omega *= W1; 353 omega *= W1; 351 omega += omega2*W2; 354 omega += omega2*W2; 352 } 355 } 353 //G4cout << "omega(keV)= " << omega/keV << 356 //G4cout << "omega(keV)= " << omega/keV << G4endl; 354 357 355 loss += omega; 358 loss += omega; 356 if(loss > kinEnergy) { break; } 359 if(loss > kinEnergy) { break; } 357 } 360 } 358 361 359 //G4cout<<"PAIModelData AlongStepLoss = "<<l 362 //G4cout<<"PAIModelData AlongStepLoss = "<<loss/keV<<" keV"<<G4endl; 360 if(loss > kinEnergy) { loss = kinEnergy; } 363 if(loss > kinEnergy) { loss = kinEnergy; } 361 else if(loss < 0.) { loss = 0.; } 364 else if(loss < 0.) { loss = 0.; } 362 return loss; 365 return loss; 363 } 366 } 364 367 365 ////////////////////////////////////////////// 368 /////////////////////////////////////////////////////////////////////// 366 // 369 // 367 // Returns post step PAI energy transfer > cut 370 // Returns post step PAI energy transfer > cut electron energy 368 // according to passed scaled kinetic energy o 371 // according to passed scaled kinetic energy of particle 369 372 370 G4double G4PAIModelData::SamplePostStepTransfe 373 G4double G4PAIModelData::SamplePostStepTransfer(G4int coupleIndex, 371 G4double scaledTkin, 374 G4double scaledTkin, 372 G4double tmin, 375 G4double tmin, 373 G4double tmax) const 376 G4double tmax) const 374 { 377 { 375 //G4cout<<"=== G4PAIModelData::SamplePostSte 378 //G4cout<<"=== G4PAIModelData::SamplePostStepTransfer idx= "<< coupleIndex 376 // << " Tkin= " << scaledTkin << " Tmax= " 379 // << " Tkin= " << scaledTkin << " Tmax= " << tmax << G4endl; 377 G4double transfer = 0.0; 380 G4double transfer = 0.0; 378 G4double rand = G4UniformRand(); 381 G4double rand = G4UniformRand(); 379 382 380 std::size_t nPlace = fParticleEnergyVector-> << 383 size_t nPlace = fParticleEnergyVector->GetVectorLength() - 1; 381 std::size_t iPlace = fParticleEnergyVector-> << 384 size_t iPlace = fParticleEnergyVector->FindBin(scaledTkin, 0); 382 385 383 G4bool one = true; 386 G4bool one = true; 384 if(scaledTkin >= fParticleEnergyVector->Ener 387 if(scaledTkin >= fParticleEnergyVector->Energy(nPlace)) { iPlace = nPlace; } 385 else if(scaledTkin > fParticleEnergyVector-> 388 else if(scaledTkin > fParticleEnergyVector->Energy(0)) { 386 one = false; 389 one = false; 387 } 390 } 388 G4PhysicsTable* table = fPAIxscBank[coupleIn 391 G4PhysicsTable* table = fPAIxscBank[coupleIndex]; 389 G4PhysicsVector* v1 = (*table)[iPlace]; 392 G4PhysicsVector* v1 = (*table)[iPlace]; 390 393 391 G4double emin = std::max(tmin, v1->Energy(0) 394 G4double emin = std::max(tmin, v1->Energy(0)); 392 G4double emax = std::min(tmax, v1->GetMaxEne 395 G4double emax = std::min(tmax, v1->GetMaxEnergy()); 393 if(emax < emin) { return transfer; } 396 if(emax < emin) { return transfer; } 394 397 395 G4double dNdx1 = v1->Value(emin)/emin; 398 G4double dNdx1 = v1->Value(emin)/emin; 396 G4double dNdx2 = v1->Value(emax)/emax; 399 G4double dNdx2 = v1->Value(emax)/emax; 397 /* 400 /* 398 G4cout << "iPlace= " << iPlace << " nPlace= 401 G4cout << "iPlace= " << iPlace << " nPlace= " << nPlace 399 << " emin= " << emin << " emax= " << emax 402 << " emin= " << emin << " emax= " << emax 400 << " dNdx1= " << dNdx1 << " dNdx2= " << dNd 403 << " dNdx1= " << dNdx1 << " dNdx2= " << dNdx2 401 << " one: " << one << G4endl; 404 << " one: " << one << G4endl; 402 */ 405 */ 403 G4double position = dNdx2 + (dNdx1 - dNdx2)* 406 G4double position = dNdx2 + (dNdx1 - dNdx2)*rand; 404 transfer = GetEnergyTransfer(coupleIndex, iP 407 transfer = GetEnergyTransfer(coupleIndex, iPlace, position); 405 408 406 //G4cout<<"PAImodel PostStepTransfer = "<<tr 409 //G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV" 407 // << " position= " << position << G4endl; 410 // << " position= " << position << G4endl; 408 411 409 if(!one) { 412 if(!one) { 410 413 411 G4PhysicsVector* v2 = (*table)[iPlace+1]; 414 G4PhysicsVector* v2 = (*table)[iPlace+1]; 412 emin = std::max(tmin, v2->Energy(0)); 415 emin = std::max(tmin, v2->Energy(0)); 413 emax = std::min(tmax, v2->GetMaxEnergy()); 416 emax = std::min(tmax, v2->GetMaxEnergy()); 414 if(emin <= emax) { 417 if(emin <= emax) { 415 dNdx1 = v2->Value(emin)/emin; 418 dNdx1 = v2->Value(emin)/emin; 416 dNdx2 = v2->Value(emax)/emax; 419 dNdx2 = v2->Value(emax)/emax; 417 420 418 //G4cout << " emax2= " << emax 421 //G4cout << " emax2= " << emax 419 // << " dNdx2= " << dNdx2 << " dNdx1 422 // << " dNdx2= " << dNdx2 << " dNdx1= " << dNdx1 << G4endl; 420 423 421 G4double E1 = fParticleEnergyVector->Ene 424 G4double E1 = fParticleEnergyVector->Energy(iPlace); 422 G4double E2 = fParticleEnergyVector->Ene 425 G4double E2 = fParticleEnergyVector->Energy(iPlace+1); 423 G4double W = 1.0/(E2 - E1); 426 G4double W = 1.0/(E2 - E1); 424 G4double W1 = (E2 - scaledTkin)*W; 427 G4double W1 = (E2 - scaledTkin)*W; 425 G4double W2 = (scaledTkin - E1)*W; 428 G4double W2 = (scaledTkin - E1)*W; 426 429 427 //G4cout<< "E1= " << E1 << " E2= " << E2 430 //G4cout<< "E1= " << E1 << " E2= " << E2 <<" iPlace= " << iPlace 428 // << " W1= " << W1 << " W2= " << W2 431 // << " W1= " << W1 << " W2= " << W2 <<G4endl; 429 432 430 position = dNdx2 + (dNdx1 - dNdx2)*rand; 433 position = dNdx2 + (dNdx1 - dNdx2)*rand; 431 G4double tr2 = GetEnergyTransfer(coupleI 434 G4double tr2 = GetEnergyTransfer(coupleIndex, iPlace+1, position); 432 435 433 //G4cout<<"PAImodel PostStepTransfer1 = 436 //G4cout<<"PAImodel PostStepTransfer1 = "<<tr2/keV<<" keV" 434 // << " position= " << position << G4 437 // << " position= " << position << G4endl; 435 transfer *= W1; 438 transfer *= W1; 436 transfer += tr2*W2; 439 transfer += tr2*W2; 437 } 440 } 438 } 441 } 439 //G4cout<<"PAImodel PostStepTransfer = "<<tr 442 //G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV" 440 // << " position= " << position << G4endl; 443 // << " position= " << position << G4endl; 441 transfer = std::max(transfer, 0.0); 444 transfer = std::max(transfer, 0.0); 442 return transfer; 445 return transfer; 443 } 446 } 444 447 445 ////////////////////////////////////////////// 448 /////////////////////////////////////////////////////////////////////// 446 // 449 // 447 // Returns PAI energy transfer according to pa 450 // Returns PAI energy transfer according to passed 448 // indexes of particle kinetic enegry and rand 451 // indexes of particle kinetic enegry and random x-section 449 452 450 G4double G4PAIModelData::GetEnergyTransfer(G4i 453 G4double G4PAIModelData::GetEnergyTransfer(G4int coupleIndex, 451 std::size_t iPlace, << 454 size_t iPlace, 452 G4double position) const 455 G4double position) const 453 { 456 { 454 G4PhysicsVector* v = (*(fPAIxscBank[coupleIn 457 G4PhysicsVector* v = (*(fPAIxscBank[coupleIndex]))(iPlace); 455 if(position*v->Energy(0) >= (*v)[0]) { retur 458 if(position*v->Energy(0) >= (*v)[0]) { return v->Energy(0); } 456 459 457 std::size_t iTransferMax = v->GetVectorLengt << 460 size_t iTransferMax = v->GetVectorLength() - 1; 458 461 459 std::size_t iTransfer; << 462 size_t iTransfer; 460 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), 463 G4double x1(0.0), x2(0.0), y1(0.0), y2(0.0), energyTransfer; 461 464 462 //G4cout << "iPlace= " << iPlace << " iTrans 465 //G4cout << "iPlace= " << iPlace << " iTransferMax= " << iTransferMax << G4endl; 463 for(iTransfer=1; iTransfer<=iTransferMax; ++ 466 for(iTransfer=1; iTransfer<=iTransferMax; ++iTransfer) { 464 x2 = v->Energy(iTransfer); 467 x2 = v->Energy(iTransfer); 465 y2 = (*v)[iTransfer]/x2; 468 y2 = (*v)[iTransfer]/x2; 466 if(position >= y2) { break; } 469 if(position >= y2) { break; } 467 if(iTransfer == iTransferMax) { return v-> 470 if(iTransfer == iTransferMax) { return v->GetMaxEnergy(); } 468 } 471 } 469 472 470 x1 = v->Energy(iTransfer-1); 473 x1 = v->Energy(iTransfer-1); 471 y1 = (*v)[iTransfer-1]/x1; 474 y1 = (*v)[iTransfer-1]/x1; 472 /* 475 /* 473 G4cout << "i= " << iTransfer << " imax= " << 476 G4cout << "i= " << iTransfer << " imax= " << iTransferMax 474 << " x1= " << x1 << " x2= " << x2 477 << " x1= " << x1 << " x2= " << x2 475 << " y1= " << y1 << " y2= " << y2 << G4en 478 << " y1= " << y1 << " y2= " << y2 << G4endl; 476 */ 479 */ 477 energyTransfer = x1; 480 energyTransfer = x1; 478 if ( x1 != x2 ) { 481 if ( x1 != x2 ) { 479 if ( y1 == y2 ) { 482 if ( y1 == y2 ) { 480 energyTransfer += (x2 - x1)*G4UniformRan 483 energyTransfer += (x2 - x1)*G4UniformRand(); 481 } else { 484 } else { 482 if(x1*1.1 < x2) { 485 if(x1*1.1 < x2) { 483 const G4int nbins = 5; 486 const G4int nbins = 5; 484 G4double del = (x2 - x1)/G4int(nbins); 487 G4double del = (x2 - x1)/G4int(nbins); 485 x2 = x1; 488 x2 = x1; 486 for(G4int i=1; i<=nbins; ++i) { 489 for(G4int i=1; i<=nbins; ++i) { 487 x2 += del; 490 x2 += del; 488 y2 = v->Value(x2)/x2; 491 y2 = v->Value(x2)/x2; 489 if(position >= y2) { 492 if(position >= y2) { 490 break; 493 break; 491 } 494 } 492 x1 = x2; 495 x1 = x2; 493 y1 = y2; 496 y1 = y2; 494 } 497 } 495 } 498 } 496 //G4cout << "x1(keV)= " << x1/keV << " x 499 //G4cout << "x1(keV)= " << x1/keV << " x2(keV)= " << x2/keV 497 // << " y1= " << y1 << " y2= " << y2 < 500 // << " y1= " << y1 << " y2= " << y2 << " pos= " << position << G4endl; 498 energyTransfer = (y2 - y1)*x1*x2/(positi 501 energyTransfer = (y2 - y1)*x1*x2/(position*(x1 - x2) - y1*x1 + y2*x2); 499 } 502 } 500 } 503 } 501 //G4cout << "x1(keV)= " << x1/keV << " x2(ke 504 //G4cout << "x1(keV)= " << x1/keV << " x2(keV)= " << x2/keV 502 // << " y1= " << y1 << " y2= " << y2 << " 505 // << " y1= " << y1 << " y2= " << y2 << " pos= " << position 503 // << " E(keV)= " << energyTransfer/keV << 506 // << " E(keV)= " << energyTransfer/keV << G4endl; 504 return energyTransfer; 507 return energyTransfer; 505 } 508 } 506 509 507 ////////////////////////////////////////////// 510 ////////////////////////////////////////////////////////////////////// 508 511 509 512