| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 // <<
27 // ------------------------------------------- <<
28 // <<
29 // GEANT4 Class <<
30 // File name: G4PAIModel.cc 23 // File name: G4PAIModel.cc
31 // 24 //
32 // Author: Vladimir.Grichine@cern.ch on base o << 25 // Author: Vladimir.Grichine@cern.ch on base of Vladimir Ivanchenko code
33 // 26 //
34 // Creation date: 05.10.2003 27 // Creation date: 05.10.2003
35 // 28 //
36 // Modifications: 29 // Modifications:
37 // 30 //
38 // 17.08.04 V.Grichine, bug fixed for Tkin<=0 31 // 17.08.04 V.Grichine, bug fixed for Tkin<=0 in SampleSecondary
39 // 16.08.04 V.Grichine, bug fixed in massRatio << 32 // 16.08.04 V.Grichine, bug fixed in massRatio for DEDX, CrossSection, SampleSecondary
40 // SampleSecondary <<
41 // 08.04.05 Major optimisation of internal int <<
42 // 26.07.09 Fixed logic to work with several m <<
43 // 21.11.10 V. Grichine verbose flag for proto <<
44 // check sandia table <<
45 // 12.06.13 V. Grichine Bug fixed in SampleSec <<
46 // (fMass -> proton_mass_c2) <<
47 // 19.08.13 V.Ivanchenko extract data handling <<
48 // added sharing of internal data bet <<
49 // 33 //
50 34
51 #include "G4PAIModel.hh" <<
52 <<
53 #include "G4SystemOfUnits.hh" <<
54 #include "G4PhysicalConstants.hh" <<
55 #include "G4Region.hh" 35 #include "G4Region.hh"
>> 36 #include "G4PhysicsLogVector.hh"
>> 37 #include "G4PhysicsFreeVector.hh"
>> 38 #include "G4PhysicsTable.hh"
>> 39 #include "G4ProductionCutsTable.hh"
56 #include "G4MaterialCutsCouple.hh" 40 #include "G4MaterialCutsCouple.hh"
57 #include "G4MaterialTable.hh" 41 #include "G4MaterialTable.hh"
58 #include "G4RegionStore.hh" << 42 #include "G4SandiaTable.hh"
>> 43 #include "G4PAIxSection.hh"
59 44
>> 45 #include "G4PAIModel.hh"
60 #include "Randomize.hh" 46 #include "Randomize.hh"
61 #include "G4Electron.hh" 47 #include "G4Electron.hh"
62 #include "G4Positron.hh" <<
63 #include "G4Poisson.hh" 48 #include "G4Poisson.hh"
64 #include "G4Step.hh" 49 #include "G4Step.hh"
65 #include "G4Material.hh" 50 #include "G4Material.hh"
66 #include "G4DynamicParticle.hh" 51 #include "G4DynamicParticle.hh"
67 #include "G4ParticleDefinition.hh" 52 #include "G4ParticleDefinition.hh"
68 #include "G4ParticleChangeForLoss.hh" <<
69 #include "G4PAIModelData.hh" <<
70 #include "G4DeltaAngle.hh" <<
71 53
72 ////////////////////////////////////////////// 54 ////////////////////////////////////////////////////////////////////////
73 55
74 using namespace std; 56 using namespace std;
75 57
76 G4PAIModel::G4PAIModel(const G4ParticleDefinit 58 G4PAIModel::G4PAIModel(const G4ParticleDefinition* p, const G4String& nam)
77 : G4VEmModel(nam),G4VEmFluctuationModel(nam) 59 : G4VEmModel(nam),G4VEmFluctuationModel(nam),
78 fVerbose(0), << 60 fLowestKineticEnergy(10.0*keV),
79 fModelData(nullptr), << 61 fHighestKineticEnergy(100.*TeV),
80 fParticle(nullptr) << 62 fTotBin(200),
81 { << 63 fMeanNumber(20),
82 fElectron = G4Electron::Electron(); << 64 fParticle(0),
83 fPositron = G4Positron::Positron(); << 65 fHighKinEnergy(100.*TeV),
>> 66 fLowKinEnergy(2.0*MeV),
>> 67 fTwoln10(2.0*log(10.0)),
>> 68 fBg2lim(0.0169),
>> 69 fTaulim(8.4146e-3)
>> 70 {
>> 71 if(p) SetParticle(p);
>> 72 fProtonEnergyVector = new G4PhysicsLogVector(fLowestKineticEnergy,
>> 73 fHighestKineticEnergy,
>> 74 fTotBin);
>> 75 fPAItransferTable = 0;
>> 76 fPAIdEdxTable = 0;
>> 77 fSandiaPhotoAbsCof = 0;
>> 78 fdEdxVector = 0;
>> 79 fLambdaVector = 0;
>> 80 fdNdxCutVector = 0;
>> 81 }
>> 82
>> 83 ////////////////////////////////////////////////////////////////////////////
>> 84
>> 85 G4PAIModel::~G4PAIModel()
>> 86 {
>> 87 if(fProtonEnergyVector) delete fProtonEnergyVector;
>> 88 if(fdEdxVector) delete fdEdxVector ;
>> 89 if ( fLambdaVector) delete fLambdaVector;
>> 90 if ( fdNdxCutVector) delete fdNdxCutVector;
>> 91
>> 92 if( fPAItransferTable )
>> 93 {
>> 94 fPAItransferTable->clearAndDestroy();
>> 95 delete fPAItransferTable ;
>> 96 }
>> 97 if(fSandiaPhotoAbsCof)
>> 98 {
>> 99 for(G4int i=0;i<fSandiaIntervalNumber;i++)
>> 100 {
>> 101 delete[] fSandiaPhotoAbsCof[i];
>> 102 }
>> 103 delete[] fSandiaPhotoAbsCof;
>> 104 }
>> 105 }
>> 106
>> 107 ///////////////////////////////////////////////////////////////////////////////
>> 108
>> 109 void G4PAIModel::SetParticle(const G4ParticleDefinition* p)
>> 110 {
>> 111 fParticle = p;
>> 112 fMass = fParticle->GetPDGMass();
>> 113 fSpin = fParticle->GetPDGSpin();
>> 114 G4double q = fParticle->GetPDGCharge()/eplus;
>> 115 fChargeSquare = q*q;
>> 116 fLowKinEnergy *= fMass/proton_mass_c2;
>> 117 fRatio = electron_mass_c2/fMass;
>> 118 fQc = fMass/fRatio;
>> 119 }
>> 120
>> 121 //////////////////////////////////////////////////////////////////////////////
84 122
85 fParticleChange = nullptr; << 123 G4double G4PAIModel::HighEnergyLimit(const G4ParticleDefinition* p)
>> 124 {
>> 125 if(!fParticle) SetParticle(p);
>> 126 return fHighKinEnergy;
>> 127 }
86 128
87 if(p) { SetParticle(p); } << 129 ///////////////////////////////////////////////////////////////////////////
88 else { SetParticle(fElectron); } <<
89 130
90 // default generator << 131 G4double G4PAIModel::LowEnergyLimit( const G4ParticleDefinition* p )
91 SetAngularDistribution(new G4DeltaAngle()); << 132 {
92 fLowestTcut = 12.5*CLHEP::eV; << 133 if(!fParticle) SetParticle(p);
>> 134 return fLowKinEnergy;
93 } 135 }
94 136
95 ////////////////////////////////////////////// 137 ////////////////////////////////////////////////////////////////////////////
96 138
97 G4PAIModel::~G4PAIModel() << 139 G4double G4PAIModel::MinEnergyCut( const G4ParticleDefinition*,
>> 140 const G4MaterialCutsCouple*)
>> 141 {
>> 142 return 0.*eV; // any positive cut
>> 143 }
>> 144
>> 145 ////////////////////////////////////////////////////////////////////////////
>> 146
>> 147 G4bool G4PAIModel::IsInCharge( const G4ParticleDefinition* p )
98 { 148 {
99 if(IsMaster()) { delete fModelData; } << 149 if(!fParticle) SetParticle(p);
>> 150 return (p->GetPDGCharge() != 0.0 );
100 } 151 }
101 152
102 ////////////////////////////////////////////// 153 ////////////////////////////////////////////////////////////////////////////
103 154
104 void G4PAIModel::Initialise(const G4ParticleDe 155 void G4PAIModel::Initialise(const G4ParticleDefinition* p,
105 const G4DataVector& cuts) << 156 const G4DataVector&)
106 { 157 {
107 if(fVerbose > 1) { << 158 if(!fParticle) SetParticle(p);
108 G4cout<<"G4PAIModel::Initialise for "<<p-> <<
109 } <<
110 SetParticle(p); <<
111 fParticleChange = GetParticleChangeForLoss() <<
112 159
113 if(IsMaster()) { << 160 const G4ProductionCutsTable* theCoupleTable =
114 << 161 G4ProductionCutsTable::GetProductionCutsTable();
115 delete fModelData; << 162
116 fMaterialCutsCoupleVector.clear(); << 163 for(size_t iReg = 0; iReg < fPAIRegionVector.size();++iReg) // region loop
117 if(fVerbose > 1) { << 164 {
118 G4cout << "G4PAIModel instantiates data << 165 const G4Region* curReg = fPAIRegionVector[iReg];
119 << G4endl; << 166
120 } << 167 // (*fPAIRegionVector[iRegion])
121 G4double tmin = LowEnergyLimit()*fRatio; << 168
122 G4double tmax = HighEnergyLimit()*fRatio; << 169 vector<G4Material*>::const_iterator matIter = curReg->GetMaterialIterator();
123 fModelData = new G4PAIModelData(tmin, tmax << 170 size_t jMat;
124 << 171 size_t numOfMat = curReg->GetNumberOfMaterials();
125 // Prepare initialization << 172
>> 173 // for(size_t jMat = 0; jMat < curReg->GetNumberOfMaterials();++jMat){}
126 const G4MaterialTable* theMaterialTable = 174 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
127 std::size_t numOfMat = G4Material::GetNu << 175 size_t numberOfMat = G4Material::GetNumberOfMaterials();
128 std::size_t numRegions = fPAIRegionVector. << 176
>> 177 for(jMat = 0 ; jMat < numOfMat; ++jMat) // region material loop
>> 178 {
>> 179 const G4MaterialCutsCouple* matCouple = theCoupleTable->
>> 180 GetMaterialCutsCouple( *matIter, curReg->GetProductionCuts() );
>> 181 fMaterialCutsCoupleVector.push_back(matCouple);
>> 182
>> 183 size_t iMatGlob;
>> 184 for(iMatGlob = 0 ; iMatGlob < numberOfMat ; iMatGlob++ )
>> 185 {
>> 186 if( *matIter == (*theMaterialTable)[iMatGlob]) break ;
>> 187 }
>> 188 fMatIndex = iMatGlob;
>> 189
>> 190 ComputeSandiaPhotoAbsCof();
>> 191 BuildPAIonisationTable();
>> 192
>> 193 fPAIxscBank.push_back(fPAItransferTable);
>> 194 fPAIdEdxBank.push_back(fPAIdEdxTable);
>> 195 fdEdxTable.push_back(fdEdxVector);
>> 196
>> 197 BuildLambdaVector(matCouple);
>> 198 fdNdxCutTable.push_back(fdNdxCutVector);
>> 199 fLambdaTable.push_back(fLambdaVector);
129 200
130 // protect for unit tests << 201
131 if(0 == numRegions) { << 202 matIter++;
132 G4Exception("G4PAIModel::Initialise()"," <<
133 "no G4Regions are registered <<
134 fPAIRegionVector.push_back(G4RegionStore <<
135 ->GetRegion("DefaultRegionForTheWorld <<
136 numRegions = 1; <<
137 } 203 }
>> 204 }
>> 205 }
138 206
139 if(fVerbose > 1) { << 207 //////////////////////////////////////////////////////////////////
140 G4cout << "G4PAIModel is defined for " < << 208
141 << "; number of materials " << numOfMat << 209 void G4PAIModel::ComputeSandiaPhotoAbsCof()
142 } << 210 {
143 for(std::size_t iReg = 0; iReg<numRegions; << 211 G4int i, j, numberOfElements ;
144 const G4Region* curReg = fPAIRegionVecto << 212 static const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
145 G4Region* reg = const_cast<G4Region*>(cu << 213
146 << 214 G4SandiaTable thisMaterialSandiaTable(fMatIndex) ;
147 for(std::size_t jMat = 0; jMat<numOfMat; << 215 numberOfElements = (*theMaterialTable)[fMatIndex]->
148 G4Material* mat = (*theMaterialTable)[jMat]; << 216 GetNumberOfElements();
149 const G4MaterialCutsCouple* cutCouple = reg- << 217 G4int* thisMaterialZ = new G4int[numberOfElements] ;
150 std::size_t n = fMaterialCutsCoupleVector.si << 218
151 if(nullptr != cutCouple) { << 219 for(i=0;i<numberOfElements;i++)
152 if(fVerbose > 1) { << 220 {
153 G4cout << "Region <" << curReg->GetName( << 221 thisMaterialZ[i] =
154 << mat->GetName() << "> CoupleIndex= " << 222 (G4int)(*theMaterialTable)[fMatIndex]->GetElement(i)->GetZ() ;
155 << cutCouple->GetIndex() << 223 }
156 << " " << p->GetParticleName() << 224 fSandiaIntervalNumber = thisMaterialSandiaTable.SandiaIntervals
157 << " cutsize= " << cuts.size() << G4end << 225 (thisMaterialZ,numberOfElements) ;
158 } << 226
159 // check if this couple is not already ini << 227 fSandiaIntervalNumber = thisMaterialSandiaTable.SandiaMixing
160 G4bool isnew = true; << 228 ( thisMaterialZ ,
161 if(0 < n) { << 229 (*theMaterialTable)[fMatIndex]->GetFractionVector() ,
162 for(std::size_t i=0; i<n; ++i) { << 230 numberOfElements,fSandiaIntervalNumber) ;
163 G4cout << i << G4endl; << 231
164 if(cutCouple == fMaterialCutsCoupleVec << 232 fSandiaPhotoAbsCof = new G4double*[fSandiaIntervalNumber] ;
165 isnew = false; << 233
166 break; << 234 for(i=0;i<fSandiaIntervalNumber;i++) fSandiaPhotoAbsCof[i] = new G4double[5] ;
167 } << 235
168 } << 236 for( i = 0 ; i < fSandiaIntervalNumber ; i++ )
169 } << 237 {
170 // initialise data banks << 238 fSandiaPhotoAbsCof[i][0] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i+1,0) ;
171 // G4cout << " isNew: " << isnew << " " << 239
172 if(isnew) { << 240 for( j = 1; j < 5 ; j++ )
173 fMaterialCutsCoupleVector.push_back(cutC << 241 {
174 fModelData->Initialise(cutCouple, this); << 242 fSandiaPhotoAbsCof[i][j] = thisMaterialSandiaTable.
175 } << 243 GetPhotoAbsorpCof(i+1,j)*
176 } << 244 (*theMaterialTable)[fMatIndex]->GetDensity() ;
177 } <<
178 } 245 }
179 InitialiseElementSelectors(p, cuts); <<
180 } 246 }
>> 247 // delete[] thisMaterialZ ;
181 } 248 }
182 249
183 ////////////////////////////////////////////// << 250 ////////////////////////////////////////////////////////////////////////////
>> 251 //
>> 252 // Build tables for the ionization energy loss
>> 253 // the tables are built for MATERIALS
>> 254 // *********
184 255
185 void G4PAIModel::InitialiseLocal(const G4Parti << 256 void
186 G4VEmModel* masterModel) << 257 G4PAIModel::BuildPAIonisationTable()
187 { 258 {
188 fModelData = static_cast<G4PAIModel*>(master << 259 G4double LowEdgeEnergy , ionloss ;
189 fMaterialCutsCoupleVector = << 260 G4double massRatio, tau, Tmax, Tmin, Tkin, deltaLow, gamma, bg2 ;
190 static_cast<G4PAIModel*>(masterModel)->Get << 261 /*
191 SetElementSelectors(masterModel->GetElementS << 262 if( fPAItransferTable )
>> 263 {
>> 264 fPAItransferTable->clearAndDestroy() ;
>> 265 delete fPAItransferTable ;
>> 266 }
>> 267 */
>> 268 fPAItransferTable = new G4PhysicsTable(fTotBin);
>> 269 /*
>> 270 if( fPAIdEdxTable )
>> 271 {
>> 272 fPAIdEdxTable->clearAndDestroy() ;
>> 273 delete fPAIdEdxTable ;
>> 274 }
>> 275 */
>> 276 fPAIdEdxTable = new G4PhysicsTable(fTotBin);
>> 277
>> 278 // if(fdEdxVector) delete fdEdxVector ;
>> 279 fdEdxVector = new G4PhysicsLogVector( fLowestKineticEnergy,
>> 280 fHighestKineticEnergy,
>> 281 fTotBin ) ;
>> 282 Tmin = fSandiaPhotoAbsCof[0][0] ; // low energy Sandia interval
>> 283 deltaLow = 0.5*eV ;
>> 284
>> 285 for (G4int i = 0 ; i < fTotBin ; i++) //The loop for the kinetic energy
>> 286 {
>> 287 LowEdgeEnergy = fProtonEnergyVector->GetLowEdgeEnergy(i) ;
>> 288 tau = LowEdgeEnergy/proton_mass_c2 ;
>> 289 // if(tau < 0.01) tau = 0.01 ;
>> 290 gamma = tau +1. ;
>> 291 // G4cout<<"gamma = "<<gamma<<endl ;
>> 292 bg2 = tau*(tau + 2. ) ;
>> 293 massRatio = electron_mass_c2/proton_mass_c2 ;
>> 294 Tmax = 2.*electron_mass_c2*bg2/(1.+2.*gamma*massRatio+massRatio*massRatio) ;
>> 295 // G4cout<<"proton Tkin = "<<LowEdgeEnergy/MeV<<" MeV"
>> 296 // <<" Tmax = "<<Tmax/MeV<<" MeV"<<G4endl;
>> 297 // Tkin = DeltaCutInKineticEnergyNow ;
>> 298
>> 299 // if ( DeltaCutInKineticEnergyNow > Tmax) // was <
>> 300 {
>> 301 Tkin = Tmax ;
>> 302 }
>> 303 if ( Tkin < Tmin + deltaLow ) // low energy safety
>> 304 {
>> 305 Tkin = Tmin + deltaLow ;
>> 306 }
>> 307 G4PAIxSection protonPAI( fMatIndex,
>> 308 Tkin,
>> 309 bg2,
>> 310 fSandiaPhotoAbsCof,
>> 311 fSandiaIntervalNumber ) ;
>> 312
>> 313
>> 314 // G4cout<<"ionloss = "<<ionloss*cm/keV<<" keV/cm"<<endl ;
>> 315 // G4cout<<"n1 = "<<protonPAI.GetIntegralPAIxSection(1)*cm<<" 1/cm"<<endl ;
>> 316 // G4cout<<"protonPAI.GetSplineSize() = "<<
>> 317 // protonPAI.GetSplineSize()<<G4endl<<G4endl ;
>> 318
>> 319 G4PhysicsFreeVector* transferVector = new
>> 320 G4PhysicsFreeVector(protonPAI.GetSplineSize()) ;
>> 321 G4PhysicsFreeVector* dEdxVector = new
>> 322 G4PhysicsFreeVector(protonPAI.GetSplineSize()) ;
>> 323
>> 324 for( G4int k = 0 ; k < protonPAI.GetSplineSize() ; k++ )
>> 325 {
>> 326 transferVector->PutValue( k ,
>> 327 protonPAI.GetSplineEnergy(k+1),
>> 328 protonPAI.GetIntegralPAIxSection(k+1) ) ;
>> 329 dEdxVector->PutValue( k ,
>> 330 protonPAI.GetSplineEnergy(k+1),
>> 331 protonPAI.GetIntegralPAIdEdx(k+1) ) ;
>> 332 }
>> 333 ionloss = protonPAI.GetMeanEnergyLoss() ; // total <dE/dx>
>> 334 if ( ionloss <= 0.) ionloss = DBL_MIN ;
>> 335 fdEdxVector->PutValue(i,ionloss) ;
>> 336
>> 337 fPAItransferTable->insertAt(i,transferVector) ;
>> 338 fPAIdEdxTable->insertAt(i,dEdxVector) ;
>> 339
>> 340 // delete[] transferVector ;
>> 341 } // end of Tkin loop
>> 342 // theLossTable->insert(fdEdxVector);
>> 343 // end of material loop
>> 344 // G4cout<<"G4PAIonisation::BuildPAIonisationTable() have been called"<<G4endl ;
>> 345 // G4cout<<"G4PAIonisation::BuildLossTable() have been called"<<G4endl ;
192 } 346 }
193 347
194 ////////////////////////////////////////////// << 348 ///////////////////////////////////////////////////////////////////////
>> 349 //
>> 350 // Build mean free path tables for the delta ray production process
>> 351 // tables are built for MATERIALS
>> 352 //
195 353
196 G4double G4PAIModel::MinEnergyCut(const G4Part << 354 void
197 const G4MaterialCutsCouple*) << 355 G4PAIModel::BuildLambdaVector(const G4MaterialCutsCouple* matCutsCouple)
198 { 356 {
199 return fLowestTcut; << 357 G4int i ;
>> 358 G4double dNdxCut, lambda;
>> 359
>> 360 const G4ProductionCutsTable* theCoupleTable=
>> 361 G4ProductionCutsTable::GetProductionCutsTable();
>> 362
>> 363 size_t numOfCouples = theCoupleTable->GetTableSize();
>> 364 size_t jMatCC;
>> 365
>> 366 for (jMatCC = 0 ; jMatCC < numOfCouples ; jMatCC++ )
>> 367 {
>> 368 if( matCutsCouple == theCoupleTable->GetMaterialCutsCouple(jMatCC) ) break;
>> 369 }
>> 370 if( jMatCC == numOfCouples && jMatCC > 0 ) jMatCC--;
>> 371
>> 372 const vector<G4double>* deltaCutInKineticEnergy = theCoupleTable->
>> 373 GetEnergyCutsVector(idxG4ElectronCut);
>> 374
>> 375 if (fLambdaVector) delete fLambdaVector;
>> 376 if (fdNdxCutVector) delete fdNdxCutVector;
>> 377
>> 378 fLambdaVector = new G4PhysicsLogVector( fLowestKineticEnergy,
>> 379 fHighestKineticEnergy,
>> 380 fTotBin ) ;
>> 381 fdNdxCutVector = new G4PhysicsLogVector( fLowestKineticEnergy,
>> 382 fHighestKineticEnergy,
>> 383 fTotBin ) ;
>> 384 G4double deltaCutInKineticEnergyNow = (*deltaCutInKineticEnergy)[jMatCC] ;
>> 385
>> 386 G4cout<<"PAIModel DeltaCutInKineticEnergyNow = "
>> 387 <<deltaCutInKineticEnergyNow/keV<<" keV"<<G4endl;
>> 388
>> 389 for ( i = 0 ; i < fTotBin ; i++ )
>> 390 {
>> 391 dNdxCut = GetdNdxCut(i,deltaCutInKineticEnergyNow) ;
>> 392 lambda = dNdxCut <= DBL_MIN ? DBL_MAX: 1.0/dNdxCut ;
>> 393
>> 394 if (lambda <= 1000*kCarTolerance) lambda = 1000*kCarTolerance ; // Mmm ???
>> 395
>> 396 fLambdaVector->PutValue(i, lambda) ;
>> 397 fdNdxCutVector->PutValue(i, dNdxCut) ;
>> 398 }
>> 399 }
>> 400
>> 401 ///////////////////////////////////////////////////////////////////////
>> 402 //
>> 403 // Returns integral PAI cross section for energy transfers >= transferCut
>> 404
>> 405 G4double
>> 406 G4PAIModel::GetdNdxCut( G4int iPlace, G4double transferCut)
>> 407 {
>> 408 G4int iTransfer;
>> 409 G4double x1, x2, y1, y2, dNdxCut;
>> 410 // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
>> 411 // G4cout<<"size = "<<G4int((*fPAItransferTable)(iPlace)->GetVectorLength())
>> 412 // <<G4endl;
>> 413 for( iTransfer = 0 ;
>> 414 iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) ;
>> 415 iTransfer++)
>> 416 {
>> 417 if(transferCut <= (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
>> 418 {
>> 419 break ;
>> 420 }
>> 421 }
>> 422 if ( iTransfer >= G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) )
>> 423 {
>> 424 iTransfer = (*fPAItransferTable)(iPlace)->GetVectorLength() - 1 ;
>> 425 }
>> 426 y1 = (*(*fPAItransferTable)(iPlace))(iTransfer-1) ;
>> 427 y2 = (*(*fPAItransferTable)(iPlace))(iTransfer) ;
>> 428 // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
>> 429 x1 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ;
>> 430 x2 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
>> 431 // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
>> 432
>> 433 if ( y1 == y2 ) dNdxCut = y2 ;
>> 434 else
>> 435 {
>> 436 // if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
>> 437 if ( fabs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
>> 438 else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;
>> 439 }
>> 440 // G4cout<<""<<dNdxCut<<G4endl;
>> 441 return dNdxCut ;
>> 442 }
>> 443 ///////////////////////////////////////////////////////////////////////
>> 444 //
>> 445 // Returns integral dEdx for energy transfers >= transferCut
>> 446
>> 447 G4double
>> 448 G4PAIModel::GetdEdxCut( G4int iPlace, G4double transferCut)
>> 449 {
>> 450 G4int iTransfer;
>> 451 G4double x1, x2, y1, y2, dEdxCut;
>> 452 // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
>> 453 // G4cout<<"size = "<<G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength())
>> 454 // <<G4endl;
>> 455 for( iTransfer = 0 ;
>> 456 iTransfer < G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) ;
>> 457 iTransfer++)
>> 458 {
>> 459 if(transferCut <= (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer))
>> 460 {
>> 461 break ;
>> 462 }
>> 463 }
>> 464 if ( iTransfer >= G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) )
>> 465 {
>> 466 iTransfer = (*fPAIdEdxTable)(iPlace)->GetVectorLength() - 1 ;
>> 467 }
>> 468 y1 = (*(*fPAIdEdxTable)(iPlace))(iTransfer-1) ;
>> 469 y2 = (*(*fPAIdEdxTable)(iPlace))(iTransfer) ;
>> 470 // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
>> 471 x1 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ;
>> 472 x2 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
>> 473 // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
>> 474
>> 475 if ( y1 == y2 ) dEdxCut = y2 ;
>> 476 else
>> 477 {
>> 478 // if ( x1 == x2 ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ;
>> 479 if ( fabs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ;
>> 480 else dEdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;
>> 481 }
>> 482 // G4cout<<""<<dEdxCut<<G4endl;
>> 483 return dEdxCut ;
200 } 484 }
201 485
202 ////////////////////////////////////////////// 486 //////////////////////////////////////////////////////////////////////////////
203 487
204 G4double G4PAIModel::ComputeDEDXPerVolume(cons << 488 G4double G4PAIModel::ComputeDEDX(const G4MaterialCutsCouple* matCC,
205 const G4ParticleDefinition* p, << 489 const G4ParticleDefinition* p,
206 G4double kineticEnergy, << 490 G4double kineticEnergy,
207 G4double cutEnergy) << 491 G4double cutEnergy)
208 { << 492 {
209 G4int coupleIndex = FindCoupleIndex(CurrentC << 493 G4int iTkin,iPlace;
210 if(0 > coupleIndex) { return 0.0; } << 494 size_t jMat;
>> 495 G4double massRatio = proton_mass_c2/p->GetPDGMass();
>> 496 G4double scaledTkin = kineticEnergy*massRatio;
>> 497 G4double charge = p->GetPDGCharge();
>> 498 G4double charge2 = charge*charge, dEdx;
>> 499
>> 500 for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
>> 501 {
>> 502 if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
>> 503 }
>> 504 if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;
>> 505
>> 506 fPAIdEdxTable = fPAIdEdxBank[jMat];
>> 507 fdEdxVector = fdEdxTable[jMat];
>> 508 for(iTkin = 0 ; iTkin < fTotBin ; iTkin++)
>> 509 {
>> 510 if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
>> 511 }
>> 512 iPlace = iTkin - 1;
>> 513 if(iPlace < 0) iPlace = 0;
>> 514 dEdx = charge2*( (*fdEdxVector)(iPlace) - GetdEdxCut(iPlace,cutEnergy) ) ;
211 515
212 G4double cut = std::min(MaxSecondaryEnergy(p << 516 if( dEdx < 0.) dEdx = 0.;
213 G4double scaledTkin = kineticEnergy*fRatio; << 517 return dEdx;
214 G4double dedx = fChargeSquare*fModelData->DE <<
215 return dedx; <<
216 } 518 }
217 519
218 ////////////////////////////////////////////// 520 /////////////////////////////////////////////////////////////////////////
219 521
220 G4double G4PAIModel::CrossSectionPerVolume( co << 522 G4double G4PAIModel::CrossSection( const G4MaterialCutsCouple* matCC,
221 const G4ParticleDefinition* p, << 523 const G4ParticleDefinition* p,
222 G4double kineticEnergy, << 524 G4double kineticEnergy,
223 G4double cutEnergy, << 525 G4double cutEnergy,
224 G4double maxEnergy ) << 526 G4double maxEnergy )
225 { << 527 {
226 G4int coupleIndex = FindCoupleIndex(CurrentC << 528 G4int iTkin,iPlace;
227 if(0 > coupleIndex) { return 0.0; } << 529 size_t jMat;
228 << 530 G4double tmax = min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy);
229 G4double tmax = std::min(MaxSecondaryEnergy( << 531 G4double massRatio = proton_mass_c2/p->GetPDGMass();
230 if(tmax <= cutEnergy) { return 0.0; } << 532 G4double scaledTkin = kineticEnergy*massRatio;
231 << 533 G4double charge = p->GetPDGCharge();
232 G4double scaledTkin = kineticEnergy*fRatio; << 534 G4double charge2 = charge*charge, cross, cross1, cross2;
233 G4double xs = fChargeSquare*fModelData->Cros << 535
234 scal << 536 for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
235 return xs; << 537 {
>> 538 if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
>> 539 }
>> 540 if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;
>> 541
>> 542 fPAItransferTable = fPAIxscBank[jMat];
>> 543
>> 544 for(iTkin = 0 ; iTkin < fTotBin ; iTkin++)
>> 545 {
>> 546 if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
>> 547 }
>> 548 iPlace = iTkin - 1;
>> 549 if(iPlace < 0) iPlace = 0;
>> 550
>> 551 // G4cout<<"iPlace = "<<iPlace<<"; tmax = "
>> 552 // <<tmax<<"; cutEnergy = "<<cutEnergy<<G4endl;
>> 553 cross1 = GetdNdxCut(iPlace,tmax) ;
>> 554 // G4cout<<"cross1 = "<<cross1<<G4endl;
>> 555 cross2 = GetdNdxCut(iPlace,cutEnergy) ;
>> 556 // G4cout<<"cross2 = "<<cross2<<G4endl;
>> 557 cross = (cross2-cross1)*charge2;
>> 558 // G4cout<<"cross = "<<cross<<G4endl;
>> 559 if( cross < DBL_MIN) cross = DBL_MIN;
>> 560 // if( cross2 < DBL_MIN) cross2 = DBL_MIN;
>> 561
>> 562 // return cross2;
>> 563 return cross;
236 } 564 }
237 565
238 ////////////////////////////////////////////// 566 ///////////////////////////////////////////////////////////////////////////
239 // 567 //
240 // It is analog of PostStepDoIt in terms of se 568 // It is analog of PostStepDoIt in terms of secondary electron.
241 // 569 //
242 << 570
243 void G4PAIModel::SampleSecondaries(std::vector << 571 G4DynamicParticle*
244 const G4MaterialCutsCouple* matCC, << 572 G4PAIModel::SampleSecondary( const G4MaterialCutsCouple* matCC,
245 const G4DynamicParticle* dp, << 573 const G4DynamicParticle* dp,
246 G4double tmin, << 574 G4double tmin,
247 G4double maxEnergy) << 575 G4double maxEnergy)
248 { 576 {
249 G4int coupleIndex = FindCoupleIndex(matCC); << 577 size_t jMat;
>> 578 for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
>> 579 {
>> 580 if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
>> 581 }
>> 582 if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;
250 583
251 //G4cout << "G4PAIModel::SampleSecondaries: << 584 fPAItransferTable = fPAIxscBank[jMat];
252 if(0 > coupleIndex) { return; } << 585 fdNdxCutVector = fdNdxCutTable[jMat];
253 586
254 SetParticle(dp->GetDefinition()); << 587 G4double tmax = min(MaxSecondaryEnergy(dp), maxEnergy);
>> 588 if( tmin >= tmax )
>> 589 {
>> 590 G4cout<<"G4PAIModel::SampleSecondary: tmin >= tmax "<<G4endl;
>> 591 }
>> 592 G4ThreeVector momentum = dp->GetMomentumDirection();
>> 593 G4double particleMass = dp->GetMass();
255 G4double kineticEnergy = dp->GetKineticEnerg 594 G4double kineticEnergy = dp->GetKineticEnergy();
256 595
257 G4double tmax = MaxSecondaryEnergy(fParticle << 596 G4double massRatio = proton_mass_c2/particleMass;
258 if(maxEnergy < tmax) { tmax = maxEnergy; } << 597 G4double scaledTkin = kineticEnergy*massRatio;
259 if(tmin >= tmax) { return; } << 598 G4double totalEnergy = kineticEnergy + particleMass;
>> 599 G4double pSquare = kineticEnergy*(totalEnergy+particleMass);
>> 600
>> 601 G4double deltaTkin = GetPostStepTransfer(scaledTkin);
>> 602 if( deltaTkin <= 0. )
>> 603 {
>> 604 G4cout<<"Tkin of secondary e- <= 0."<<G4endl;
>> 605 G4cout<<"G4PAIModel::SampleSecondary::deltaTkin = "<<deltaTkin<<G4endl;
>> 606 deltaTkin = 10*eV;
>> 607 G4cout<<"Set G4PAIModel::SampleSecondary::deltaTkin = "<<deltaTkin<<G4endl;
>> 608 }
>> 609 if(deltaTkin > kineticEnergy &&
>> 610 particleMass != electron_mass_c2) deltaTkin = kineticEnergy;
>> 611 if (deltaTkin > 0.5*kineticEnergy &&
>> 612 dp->GetDefinition()->GetParticleName() == "e-") deltaTkin = 0.5*kineticEnergy;
>> 613
>> 614 G4double deltaTotalMomentum = sqrt(deltaTkin*(deltaTkin + 2. * electron_mass_c2 ));
>> 615 G4double totalMomentum = sqrt(pSquare);
>> 616 G4double costheta = deltaTkin*(totalEnergy + electron_mass_c2)
>> 617 /(deltaTotalMomentum * totalMomentum);
>> 618 if( costheta >= 0.99999 ) costheta = 0.99999;
>> 619 G4double sintheta, sin2 = 1. - costheta*costheta;
260 620
261 G4ThreeVector direction= dp->GetMomentumDire << 621 if( sin2 <= 0.) sintheta = 0.;
262 G4double scaledTkin = kineticEnergy*fRati << 622 else sintheta = sqrt(sin2);
263 G4double totalEnergy = kineticEnergy + fMa <<
264 G4double totalMomentum = sqrt(kineticEnergy* <<
265 623
266 G4double deltaTkin = << 624 // direction of the delta electron
267 fModelData->SamplePostStepTransfer(coupleI << 625
>> 626 G4double phi = twopi*G4UniformRand();
>> 627 G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta;
268 628
269 //G4cout<<"G4PAIModel::SampleSecondaries; de << 629 G4ThreeVector deltaDirection(dirx,diry,dirz);
270 // <<" keV "<< " Escaled(MeV)= " << scaledT << 630 deltaDirection.rotateUz(momentum);
>> 631 deltaDirection.unit();
271 632
272 if( !(deltaTkin <= 0.) && !(deltaTkin > 0)) << 633 // create G4DynamicParticle object for e- delta ray
273 G4cout<<"G4PAIModel::SampleSecondaries; de << 634
274 <<" keV "<< " Escaled(MeV)= " << scaledTki << 635 G4DynamicParticle* deltaRay = new G4DynamicParticle;
275 return; << 636 deltaRay->SetDefinition(G4Electron::Electron());
276 } << 637 deltaRay->SetKineticEnergy( deltaTkin ); // !!! trick for last steps /2.0 ???
277 if( deltaTkin <= 0.) { return; } << 638 deltaRay->SetMomentumDirection(deltaDirection);
278 639
279 if( deltaTkin > tmax) { deltaTkin = tmax; } << 640 return deltaRay;
>> 641 }
280 642
281 const G4Element* anElement = SelectTargetAto <<
282 <<
283 643
284 G4int Z = anElement->GetZasInt(); << 644 ///////////////////////////////////////////////////////////////////////
285 << 645 //
286 auto deltaRay = new G4DynamicParticle(fElect << 646 // Returns post step PAI energy transfer > cut electron energy according to passed
287 GetAngularDistribution()->SampleDirectio << 647 // scaled kinetic energy of particle
288 Z, matCC->GetMaterial()), << 648
289 deltaTkin); << 649 G4double
290 << 650 G4PAIModel::GetPostStepTransfer( G4double scaledTkin )
291 // primary change << 651 {
292 kineticEnergy -= deltaTkin; << 652 // G4cout<<"G4PAIModel::GetPostStepTransfer"<<G4endl ;
293 G4ThreeVector dir = totalMomentum*direction << 653
294 direction = dir.unit(); << 654 G4int iTkin, iTransfer, iPlace ;
295 fParticleChange->SetProposedKineticEnergy(ki << 655 G4double transfer = 0.0, position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W ;
296 fParticleChange->SetProposedMomentumDirectio << 656
>> 657 for(iTkin=0;iTkin<fTotBin;iTkin++)
>> 658 {
>> 659 if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
>> 660 }
>> 661 iPlace = iTkin - 1 ;
>> 662 // G4cout<<"from search, iPlace = "<<iPlace<<G4endl ;
>> 663 if(iPlace < 0) iPlace = 0;
>> 664 dNdxCut1 = (*fdNdxCutVector)(iPlace) ;
>> 665 // G4cout<<"dNdxCut1 = "<<dNdxCut1<<G4endl ;
297 666
298 vdp->push_back(deltaRay); << 667
>> 668 if(iTkin == fTotBin) // Fermi plato, try from left
>> 669 {
>> 670 position = dNdxCut1*G4UniformRand() ;
>> 671
>> 672 for( iTransfer = 0;
>> 673 iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()); iTransfer++ )
>> 674 {
>> 675 if(position >= (*(*fPAItransferTable)(iPlace))(iTransfer)) break ;
>> 676 }
>> 677 transfer = GetEnergyTransfer(iPlace,position,iTransfer);
>> 678 }
>> 679 else
>> 680 {
>> 681 dNdxCut2 = (*fdNdxCutVector)(iPlace+1) ;
>> 682 // G4cout<<"dNdxCut2 = "<<dNdxCut2<<G4endl ;
>> 683 if(iTkin == 0) // Tkin is too small, trying from right only
>> 684 {
>> 685 position = dNdxCut2*G4UniformRand() ;
>> 686
>> 687 for( iTransfer = 0;
>> 688 iTransfer < G4int((*fPAItransferTable)(iPlace+1)->GetVectorLength()); iTransfer++ )
>> 689 {
>> 690 if(position >= (*(*fPAItransferTable)(iPlace+1))(iTransfer)) break ;
>> 691 }
>> 692 transfer = GetEnergyTransfer(iPlace+1,position,iTransfer);
>> 693 }
>> 694 else // general case: Tkin between two vectors of the material
>> 695 {
>> 696 E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1) ;
>> 697 E2 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin) ;
>> 698 W = 1.0/(E2 - E1) ;
>> 699 W1 = (E2 - scaledTkin)*W ;
>> 700 W2 = (scaledTkin - E1)*W ;
>> 701
>> 702 position = ( dNdxCut1*W1 + dNdxCut2*W2 )*G4UniformRand() ;
>> 703
>> 704 // G4cout<<position<<"\t" ;
>> 705
>> 706 for( iTransfer = 0;
>> 707 iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()); iTransfer++ )
>> 708 {
>> 709 if( position >=
>> 710 ( (*(*fPAItransferTable)(iPlace))(iTransfer)*W1 +
>> 711 (*(*fPAItransferTable)(iPlace+1))(iTransfer)*W2) ) break ;
>> 712 }
>> 713 transfer = GetEnergyTransfer(iPlace,position,iTransfer);
>> 714 }
>> 715 }
>> 716 // G4cout<<"PAImodel PostStepTransfer = "<<transfer/keV<<" keV"<<G4endl ;
>> 717 if(transfer < 0.0 ) transfer = 0.0 ;
>> 718 return transfer ;
299 } 719 }
300 720
301 ////////////////////////////////////////////// 721 ///////////////////////////////////////////////////////////////////////
>> 722 //
>> 723 // Returns random PAI energy transfer according to passed
>> 724 // indexes of particle kinetic
>> 725
>> 726 G4double
>> 727 G4PAIModel::GetEnergyTransfer( G4int iPlace, G4double position, G4int iTransfer )
>> 728 {
>> 729 G4double x1, x2, y1, y2, energyTransfer ;
>> 730
>> 731 if(iTransfer == 0)
>> 732 {
>> 733 energyTransfer = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
>> 734 }
>> 735 else
>> 736 {
>> 737 if ( iTransfer >= G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) )
>> 738 {
>> 739 iTransfer = (*fPAItransferTable)(iPlace)->GetVectorLength() - 1 ;
>> 740 }
>> 741 y1 = (*(*fPAItransferTable)(iPlace))(iTransfer-1) ;
>> 742 y2 = (*(*fPAItransferTable)(iPlace))(iTransfer) ;
>> 743
>> 744 x1 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ;
>> 745 x2 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
>> 746
>> 747 if ( x1 == x2 ) energyTransfer = x2 ;
>> 748 else
>> 749 {
>> 750 if ( y1 == y2 ) energyTransfer = x1 + (x2 - x1)*G4UniformRand() ;
>> 751 else
>> 752 {
>> 753 energyTransfer = x1 + (position - y1)*(x2 - x1)/(y2 - y1) ;
>> 754 }
>> 755 }
>> 756 }
>> 757 return energyTransfer ;
>> 758 }
302 759
303 G4double G4PAIModel::SampleFluctuations(const << 760
304 const << 761 ////////////////////////////////////////////////////////////////////////////
305 const << 762
306 const G4double, << 763 vector<G4DynamicParticle*>*
307 const << 764 G4PAIModel::SampleSecondaries( const G4MaterialCutsCouple*,
308 const G4double eloss) << 765 const G4DynamicParticle*,
>> 766 G4double,
>> 767 G4double)
309 { 768 {
310 G4int coupleIndex = FindCoupleIndex(matCC); << 769 vector<G4DynamicParticle*>* vdp = 0;
311 if(0 > coupleIndex) { return eloss; } << 770 // vector<G4DynamicParticle*>* vdp = new vector<G4DynamicParticle*>;
>> 771 // G4DynamicParticle* delta = SampleSecondary(couple, dp, tmin, maxEnergy);
>> 772 // vdp->push_back(delta);
>> 773 return vdp;
>> 774 }
312 775
313 SetParticle(aParticle->GetDefinition()); << 776 ///////////////////////////////////////////////////////////////////////
314 777
315 /* << 778 G4double G4PAIModel::SampleFluctuations( const G4Material* material,
316 G4cout << "G4PAIModel::SampleFluctuations st << 779 const G4DynamicParticle* aParticle,
317 << " Eloss(keV)= " << eloss/keV << " in " << 780 G4double&,
318 << matCC->Getmaterial()->GetName() << G4end << 781 G4double& step,
319 */ << 782 G4double&)
>> 783 {
>> 784 size_t jMat;
>> 785 for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
>> 786 {
>> 787 if( material == fMaterialCutsCoupleVector[jMat]->GetMaterial() ) break;
>> 788 }
>> 789 if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;
320 790
321 G4double Tkin = aParticle->GetKineticE << 791 fPAItransferTable = fPAIxscBank[jMat];
322 G4double scaledTkin = Tkin*fRatio; << 792 fdNdxCutVector = fdNdxCutTable[jMat];
323 793
324 G4double loss = fModelData->SampleAlongStepT << 794 G4int iTkin, iTransfer, iPlace ;
325 scaledTkin, tcut, << 795 G4long numOfCollisions=0;
326 step*fChargeSquare); << 796
327 << 797 // G4cout<<"G4PAIModel::SampleFluctuations"<<G4endl ;
328 // G4cout<<"PAIModel AlongStepLoss = "<<loss << 798 // G4cout<<"in: "<<fMaterialCutsCoupleVector[jMat]->GetMaterial()->GetName()<<G4endl ;
329 //<<step/mm<<" mm"<<G4endl; << 799
330 return loss; << 800 G4double loss = 0.0, charge2 ;
>> 801 G4double stepSum = 0., stepDelta, lambda, omega;
>> 802 G4double position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber;
>> 803 G4bool numb = true;
>> 804 G4double Tkin = aParticle->GetKineticEnergy() ;
>> 805 G4double MassRatio = proton_mass_c2/aParticle->GetDefinition()->GetPDGMass() ;
>> 806 G4double charge = aParticle->GetDefinition()->GetPDGCharge() ;
>> 807 charge2 = charge*charge ;
>> 808 G4double TkinScaled = Tkin*MassRatio ;
>> 809
>> 810 for(iTkin=0;iTkin<fTotBin;iTkin++)
>> 811 {
>> 812 if(TkinScaled < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
>> 813 }
>> 814 iPlace = iTkin - 1 ;
>> 815 // G4cout<<"from search, iPlace = "<<iPlace<<G4endl ;
>> 816 dNdxCut1 = (*fdNdxCutVector)(iPlace) ;
>> 817 // G4cout<<"dNdxCut1 = "<<dNdxCut1<<G4endl ;
>> 818
>> 819
>> 820 if(iTkin == fTotBin) // Fermi plato, try from left
>> 821 {
>> 822 meanNumber =((*(*fPAItransferTable)(iPlace))(0)-dNdxCut1)*step*charge2;
>> 823 if(meanNumber < 0.) meanNumber = 0. ;
>> 824 // numOfCollisions = RandPoisson::shoot(meanNumber) ;
>> 825 // numOfCollisions = G4Poisson(meanNumber) ;
>> 826 if( meanNumber > 0.) lambda = step/meanNumber;
>> 827 else lambda = DBL_MAX;
>> 828 while(numb)
>> 829 {
>> 830 stepDelta = RandExponential::shoot(lambda);
>> 831 stepSum += stepDelta;
>> 832 if(stepSum >= step) break;
>> 833 numOfCollisions++;
>> 834 }
>> 835 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl ;
>> 836
>> 837 while(numOfCollisions)
>> 838 {
>> 839 position = dNdxCut1+
>> 840 ((*(*fPAItransferTable)(iPlace))(0)-dNdxCut1)*G4UniformRand() ;
>> 841
>> 842 for( iTransfer = 0;
>> 843 iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()); iTransfer++ )
>> 844 {
>> 845 if(position >= (*(*fPAItransferTable)(iPlace))(iTransfer)) break ;
>> 846 }
>> 847 omega = GetEnergyTransfer(iPlace,position,iTransfer);
>> 848 // G4cout<<omega/keV<<"\t";
>> 849 loss += omega;
>> 850 numOfCollisions-- ;
>> 851 }
>> 852 }
>> 853 else
>> 854 {
>> 855 dNdxCut2 = (*fdNdxCutVector)(iPlace+1) ;
>> 856 // G4cout<<"dNdxCut2 = "<<dNdxCut2<<G4endl ;
>> 857
>> 858 if(iTkin == 0) // Tkin is too small, trying from right only
>> 859 {
>> 860 meanNumber =((*(*fPAItransferTable)(iPlace+1))(0)-dNdxCut2)*step*charge2;
>> 861 if( meanNumber < 0. ) meanNumber = 0. ;
>> 862 // numOfCollisions = RandPoisson::shoot(meanNumber) ;
>> 863 // numOfCollisions = G4Poisson(meanNumber) ;
>> 864 if( meanNumber > 0.) lambda = step/meanNumber;
>> 865 else lambda = DBL_MAX;
>> 866 while(numb)
>> 867 {
>> 868 stepDelta = RandExponential::shoot(lambda);
>> 869 stepSum += stepDelta;
>> 870 if(stepSum >= step) break;
>> 871 numOfCollisions++;
>> 872 }
>> 873
>> 874 // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl ;
>> 875
>> 876 while(numOfCollisions)
>> 877 {
>> 878 position = dNdxCut2+
>> 879 ((*(*fPAItransferTable)(iPlace+1))(0)-dNdxCut2)*G4UniformRand();
>> 880
>> 881 for( iTransfer = 0;
>> 882 iTransfer < G4int((*fPAItransferTable)(iPlace+1)->GetVectorLength()); iTransfer++ )
>> 883 {
>> 884 if(position >= (*(*fPAItransferTable)(iPlace+1))(iTransfer)) break ;
>> 885 }
>> 886 omega = GetEnergyTransfer(iPlace,position,iTransfer);
>> 887 // G4cout<<omega/keV<<"\t";
>> 888 loss += omega;
>> 889 numOfCollisions-- ;
>> 890 }
>> 891 }
>> 892 else // general case: Tkin between two vectors of the material
>> 893 {
>> 894 E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1) ;
>> 895 E2 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin) ;
>> 896 W = 1.0/(E2 - E1) ;
>> 897 W1 = (E2 - TkinScaled)*W ;
>> 898 W2 = (TkinScaled - E1)*W ;
>> 899
>> 900 // G4cout<<"(*(*fPAItransferTable)(iPlace))(0) = "<<
>> 901 // (*(*fPAItransferTable)(iPlace))(0)<<G4endl ;
>> 902 // G4cout<<"(*(*fPAItransferTable)(iPlace+1))(0) = "<<
>> 903 // (*(*fPAItransferTable)(iPlace+1))(0)<<G4endl ;
>> 904
>> 905 meanNumber=( ((*(*fPAItransferTable)(iPlace))(0)-dNdxCut1)*W1 +
>> 906 ((*(*fPAItransferTable)(iPlace+1))(0)-dNdxCut2)*W2 )*step*charge2;
>> 907 if(meanNumber<0.0) meanNumber = 0.0;
>> 908 // numOfCollisions = RandPoisson::shoot(meanNumber) ;
>> 909 // numOfCollisions = G4Poisson(meanNumber) ;
>> 910 if( meanNumber > 0.) lambda = step/meanNumber;
>> 911 else lambda = DBL_MAX;
>> 912 while(numb)
>> 913 {
>> 914 stepDelta = RandExponential::shoot(lambda);
>> 915 stepSum += stepDelta;
>> 916 if(stepSum >= step) break;
>> 917 numOfCollisions++;
>> 918 }
>> 919
>> 920 // G4cout<<"numOfCollisions = "<<numOfCollisions<<endl ;
>> 921
>> 922 while(numOfCollisions)
>> 923 {
>> 924 position = dNdxCut1*W1 + dNdxCut2*W2 +
>> 925 ( ( (*(*fPAItransferTable)(iPlace))(0)-dNdxCut1 )*W1 +
>> 926 dNdxCut2+
>> 927 ( (*(*fPAItransferTable)(iPlace+1))(0)-dNdxCut2 )*W2 )*G4UniformRand();
>> 928
>> 929 // G4cout<<position<<"\t" ;
>> 930
>> 931 for( iTransfer = 0;
>> 932 iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()); iTransfer++ )
>> 933 {
>> 934 if( position >=
>> 935 ( (*(*fPAItransferTable)(iPlace))(iTransfer)*W1 +
>> 936 (*(*fPAItransferTable)(iPlace+1))(iTransfer)*W2) )
>> 937 {
>> 938 break ;
>> 939 }
>> 940 }
>> 941 omega = GetEnergyTransfer(iPlace,position,iTransfer);
>> 942 // G4cout<<omega/keV<<"\t";
>> 943 loss += omega;
>> 944 numOfCollisions-- ;
>> 945 }
>> 946 }
>> 947 }
>> 948 // G4cout<<"PAIModel AlongStepLoss = "<<loss/keV<<" keV"<<G4endl ;
>> 949 if(loss > Tkin) loss=Tkin;
>> 950 return loss ;
331 951
332 } 952 }
333 953
334 ////////////////////////////////////////////// 954 //////////////////////////////////////////////////////////////////////
335 // 955 //
336 // Returns the statistical estimation of the e 956 // Returns the statistical estimation of the energy loss distribution variance
337 // 957 //
338 958
339 959
340 G4double G4PAIModel::Dispersion( const G4Mater 960 G4double G4PAIModel::Dispersion( const G4Material* material,
341 const G4Dynam 961 const G4DynamicParticle* aParticle,
342 const G4double tcut, << 962 G4double& tmax,
343 const G4double tmax, << 963 G4double& step )
344 const G4double step ) <<
345 { 964 {
346 G4double particleMass = aParticle->GetMass( << 965 G4double loss, sumLoss=0., sumLoss2=0., sigma2, meanLoss=0.;
347 G4double electronDensity = material->GetElec << 966 for(G4int i = 0 ; i < fMeanNumber; i++)
348 G4double kineticEnergy = aParticle->GetKinet << 967 {
349 G4double q = aParticle->GetCharge()/eplus; << 968 loss = SampleFluctuations(material,aParticle,tmax,step,meanLoss);
350 G4double etot = kineticEnergy + particleMass << 969 sumLoss += loss;
351 G4double beta2 = kineticEnergy*(kineticEnerg << 970 sumLoss2 += loss*loss;
352 G4double siga = (tmax/beta2 - 0.5*tcut) * t << 971 }
353 * electronDensity * q * q; << 972 meanLoss = sumLoss/fMeanNumber;
354 << 973 sigma2 = meanLoss*meanLoss + (sumLoss2-2*sumLoss*meanLoss)/fMeanNumber;
355 return siga; << 974 return sigma2;
356 } 975 }
357 976
358 ////////////////////////////////////////////// <<
359 977
360 G4double G4PAIModel::MaxSecondaryEnergy( const << 978 //
361 G4double kinEnergy) << 979 //
362 { << 980 /////////////////////////////////////////////////
363 SetParticle(p); << 981
364 G4double tmax = kinEnergy; << 982
365 if(p == fElectron) { tmax *= 0.5; } <<
366 else if(p != fPositron) { <<
367 G4double ratio= electron_mass_c2/fMass; <<
368 G4double gamma= kinEnergy/fMass + 1.0; <<
369 tmax = 2.0*electron_mass_c2*(gamma*gamma - <<
370 (1. + 2.0*gamma*ratio + rati <<
371 } <<
372 return tmax; <<
373 } <<
374 983
375 ////////////////////////////////////////////// <<
376 984
377 void G4PAIModel::DefineForRegion(const G4Regio <<
378 { <<
379 fPAIRegionVector.push_back(r); <<
380 } <<
381 985
382 ////////////////////////////////////////////// <<
383 986
384 987