| 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 // CPA100 elastic model class for electrons 26 // CPA100 elastic model class for electrons
27 // 27 //
28 // Based on the work of M. Terrissol and M. C. 28 // Based on the work of M. Terrissol and M. C. Bordage
29 // 29 //
30 // Users are requested to cite the following p 30 // Users are requested to cite the following papers:
31 // - M. Terrissol, A. Baudre, Radiat. Prot. Do 31 // - M. Terrissol, A. Baudre, Radiat. Prot. Dosim. 31 (1990) 175-177
32 // - M.C. Bordage, J. Bordes, S. Edel, M. Terr << 32 // - M.C. Bordage, J. Bordes, S. Edel, M. Terrissol, X. Franceries,
33 // M. Bardies, N. Lampe, S. Incerti, Phys. M 33 // M. Bardies, N. Lampe, S. Incerti, Phys. Med. 32 (2016) 1833-1840
34 // 34 //
35 // Authors of this class: << 35 // Authors of this class:
36 // M.C. Bordage, M. Terrissol, S. Edel, J. Bor 36 // M.C. Bordage, M. Terrissol, S. Edel, J. Bordes, S. Incerti
37 // 37 //
38 // 15.01.2014: creation 38 // 15.01.2014: creation
39 // 39 //
40 // Based on the study by S. Zein et. al. Nucl. <<
41 // 1/2/2023 : Hoang added modification <<
42 40
43 #include "G4DNACPA100ElasticModel.hh" 41 #include "G4DNACPA100ElasticModel.hh"
44 << 42 #include "G4PhysicalConstants.hh"
45 #include "G4DNAMaterialManager.hh" <<
46 #include "G4DNAMolecularMaterial.hh" <<
47 #include "G4EnvironmentUtils.hh" <<
48 #include "G4SystemOfUnits.hh" 43 #include "G4SystemOfUnits.hh"
>> 44 #include "G4DNAMolecularMaterial.hh"
>> 45
>> 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 47
49 using namespace std; 48 using namespace std;
50 G4DNACPA100ElasticModel::G4DNACPA100ElasticMod <<
51 : G4VDNAModel(nam, "all") <<
52 { <<
53 fpGuanine = G4Material::GetMaterial("G4_GUAN <<
54 fpG4_WATER = G4Material::GetMaterial("G4_WAT <<
55 fpDeoxyribose = G4Material::GetMaterial("G4_ <<
56 fpCytosine = G4Material::GetMaterial("G4_CYT <<
57 fpThymine = G4Material::GetMaterial("G4_THYM <<
58 fpAdenine = G4Material::GetMaterial("G4_ADEN <<
59 fpPhosphate = G4Material::GetMaterial("G4_PH <<
60 fpParticle = G4Electron::ElectronDefinition( <<
61 } <<
62 49
63 void G4DNACPA100ElasticModel::Initialise(const << 50 // #define CPA100_VERBOSE
64 const << 51
65 { << 52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
66 if (isInitialised) { << 53
67 return; << 54 G4DNACPA100ElasticModel::G4DNACPA100ElasticModel(const G4ParticleDefinition*,
68 } << 55 const G4String& nam)
69 if (verboseLevel > 3) { << 56 :G4VEmModel(nam),isInitialised(false)
70 G4cout << "Calling G4DNACPA100ExcitationMo << 57 {
>> 58
>> 59 SetLowEnergyLimit(11*eV);
>> 60 SetHighEnergyLimit(255955*eV);
>> 61
>> 62 verboseLevel= 0;
>> 63 // Verbosity scale:
>> 64 // 0 = nothing
>> 65 // 1 = warning for energy non-conservation
>> 66 // 2 = details of energy budget
>> 67 // 3 = calculation of cross sections, file openings, sampling of atoms
>> 68 // 4 = entering in methods
>> 69
>> 70 #ifdef UEHARA_VERBOSE
>> 71 if( verboseLevel>0 )
>> 72 {
>> 73 G4cout << "CPA100 Elastic model is constructed " << G4endl
>> 74 << "Energy range: "
>> 75 << LowEnergyLimit()/eV << " eV - "
>> 76 << HighEnergyLimit()/ keV << " keV"
>> 77 << G4endl;
71 } 78 }
>> 79 #endif
>> 80
>> 81 fParticleChangeForGamma = 0;
>> 82 fpMolWaterDensity = 0;
72 83
73 if (!G4DNAMaterialManager::Instance()->IsLoc << 84 // Selection of stationary mode
74 if (p != fpParticle) { <<
75 std::ostringstream oss; <<
76 oss << " Model is not applied for this p <<
77 G4Exception("G4DNACPA100ElasticModel::G4 <<
78 oss.str().c_str()); <<
79 } <<
80 85
81 const char* path = G4FindDataDir("G4LEDATA << 86 statCode = false;
>> 87 }
82 88
83 if (path == nullptr) { << 89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
84 G4Exception("G4DNACPA100ElasticModel::In <<
85 "G4LEDATA environment variab <<
86 return; <<
87 } <<
88 90
89 std::size_t index; << 91 G4DNACPA100ElasticModel::~G4DNACPA100ElasticModel()
90 if (fpG4_WATER != nullptr) { << 92 {
91 index = fpG4_WATER->GetIndex(); << 93 // For total cross section
92 fLevels[index] = 1.214e-4; << 94
93 AddCrossSectionData(index, p, "dna/sigma << 95 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
94 "dna/sigmadiff_cumul << 96 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
95 SetLowELimit(index, p, 11. * eV); << 97 {
96 SetHighELimit(index, p, 255955. * eV); << 98 G4DNACrossSectionDataSet* table = pos->second;
97 } << 99 delete table;
98 if (fpGuanine != nullptr) { << 100 }
99 index = fpGuanine->GetIndex(); <<
100 fLevels[index] = 1.4504480e-05; <<
101 AddCrossSectionData(index, p, "dna/sigma <<
102 "dna/sigmadiff_cumul <<
103 SetLowELimit(index, p, 11 * eV); <<
104 SetHighELimit(index, p, 1 * MeV); <<
105 } <<
106 if (fpDeoxyribose != nullptr) { <<
107 index = fpDeoxyribose->GetIndex(); <<
108 fLevels[index] = 1.6343100e-05; <<
109 AddCrossSectionData(index, p, "dna/sigma <<
110 "dna/sigmadiff_cumul <<
111 SetLowELimit(index, p, 11 * eV); <<
112 SetHighELimit(index, p, 1 * MeV); <<
113 } <<
114 if (fpCytosine != nullptr) { <<
115 index = fpCytosine->GetIndex(); <<
116 fLevels[index] = 1.9729660e-05; <<
117 AddCrossSectionData(index, p, "dna/sigma <<
118 "dna/sigmadiff_cumul <<
119 SetLowELimit(index, p, 11 * eV); <<
120 SetHighELimit(index, p, 1 * MeV); <<
121 } <<
122 if (fpThymine != nullptr) { <<
123 index = fpThymine->GetIndex(); <<
124 fLevels[index] = 1.7381300e-05; <<
125 AddCrossSectionData(index, p, "dna/sigma <<
126 "dna/sigmadiff_cumul <<
127 SetLowELimit(index, p, 11 * eV); <<
128 SetHighELimit(index, p, 1 * MeV); <<
129 } <<
130 if (fpAdenine != nullptr) { <<
131 index = fpAdenine->GetIndex(); <<
132 fLevels[index] = 1.6221800e-05; <<
133 AddCrossSectionData(index, p, "dna/sigma <<
134 "dna/sigmadiff_cumul <<
135 SetLowELimit(index, p, 11 * eV); <<
136 SetHighELimit(index, p, 1 * MeV); <<
137 } <<
138 if (fpPhosphate != nullptr) { <<
139 index = fpPhosphate->GetIndex(); <<
140 fLevels[index] = 2.2369600e-05; <<
141 AddCrossSectionData(index, p, "dna/sigma <<
142 "dna/sigmadiff_cumul <<
143 SetLowELimit(index, p, 11 * eV); <<
144 SetHighELimit(index, p, 1 * MeV); <<
145 } <<
146 101
147 // Load data << 102 // For final state
148 LoadCrossSectionData(p); << 103 eVecm.clear();
149 G4DNAMaterialManager::Instance()->SetMaste << 104
150 fpModelData = this; << 105 }
>> 106
>> 107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 108
>> 109 void G4DNACPA100ElasticModel::Initialise(const G4ParticleDefinition*
>> 110 particle,
>> 111 const G4DataVector& /*cuts*/)
>> 112 {
>> 113
>> 114 #ifdef UEHARA_VERBOSE
>> 115 if (verboseLevel > 3)
>> 116 G4cout << "Calling G4DNACPA100ElasticModel::Initialise()" << G4endl;
>> 117 #endif
>> 118
>> 119 if(particle->GetParticleName() != "e-")
>> 120 {
>> 121 G4Exception("*** WARNING: the G4DNACPA100ElasticModel is "
>> 122 "not intented to be used with another particle than the electron",
>> 123 "",FatalException,"") ;
>> 124 }
>> 125
>> 126 // Energy limits
>> 127
>> 128 if (LowEnergyLimit() < 11.*eV)
>> 129 {
>> 130 G4cout << "G4DNACPA100ElasticModel: low energy limit increased from " <<
>> 131 LowEnergyLimit()/eV << " eV to " << 11 << " eV" << G4endl;
>> 132 SetLowEnergyLimit(11.*eV);
>> 133 }
>> 134
>> 135 if (HighEnergyLimit() > 255955.*eV)
>> 136 {
>> 137 G4cout << "G4DNACPA100ElasticModel: high energy limit decreased from " <<
>> 138 HighEnergyLimit()/keV << " keV to " << 255.955 << " keV"
>> 139 << G4endl;
>> 140 SetHighEnergyLimit(255955.*eV);
>> 141 }
>> 142
>> 143 // Reading of data files
>> 144
>> 145 G4double scaleFactor = 1e-20*m*m;
>> 146
>> 147 G4String fileElectron("dna/sigma_elastic_e_cpa100");
>> 148
>> 149 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
>> 150 G4String electron;
>> 151
>> 152 // *** ELECTRON
>> 153
>> 154 // For total cross section
>> 155
>> 156 electron = electronDef->GetParticleName();
>> 157
>> 158 tableFile[electron] = fileElectron;
>> 159
>> 160 G4DNACrossSectionDataSet* tableE =
>> 161 new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
>> 162 eV,scaleFactor );
>> 163
>> 164 /*
>> 165 G4DNACrossSectionDataSet* tableE =
>> 166 new G4DNACrossSectionDataSet(new G4DNACPA100LogLogInterpolation,
>> 167 eV,scaleFactor );
>> 168 */
>> 169
>> 170 tableE->LoadData(fileElectron);
>> 171
>> 172 tableData[electron] = tableE;
>> 173
>> 174 // For final state
>> 175
>> 176 char *path = getenv("G4LEDATA");
>> 177
>> 178 if (!path)
>> 179 {
>> 180 G4Exception("G4DNACPA100ElasticModel::Initialise","em0006",
>> 181 FatalException,"G4LEDATA environment variable not set.");
>> 182 return;
151 } 183 }
152 else { << 184
153 auto dataModel = dynamic_cast<G4DNACPA100E << 185 std::ostringstream eFullFileName;
154 G4DNAMaterialManager::Instance()->GetMod << 186
155 if (dataModel == nullptr) { << 187 eFullFileName << path
156 G4cout << "G4DNACPA100ElasticModel::Cros << 188 << "/dna/sigmadiff_cumulated_elastic_e_cpa100.dat";
157 G4Exception("G4DNACPA100ElasticModel::Cr << 189
158 "no modelData is registered" << 190 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
159 } << 191
160 else { << 192 if (!eDiffCrossSection)
161 fpModelData = dataModel; << 193 G4Exception("G4DNACPA100ElasticModel::Initialise","em0003",
162 } << 194 FatalException,
>> 195 "Missing data file:/dna/sigmadiff_cumulated_elastic_e_cpa100.dat");
>> 196
>> 197 // March 25th, 2014 - Vaclav Stepan, Sebastien Incerti
>> 198 // Added clear for MT
>> 199
>> 200 eTdummyVec.clear();
>> 201 eVecm.clear();
>> 202 eDiffCrossSectionData.clear();
>> 203
>> 204 //
>> 205
>> 206 eTdummyVec.push_back(0.);
>> 207
>> 208 while(!eDiffCrossSection.eof())
>> 209 {
>> 210 G4double tDummy;
>> 211 G4double eDummy;
>> 212 eDiffCrossSection>>tDummy>>eDummy;
>> 213
>> 214 // SI : mandatory eVecm initialization
>> 215
>> 216 if (tDummy != eTdummyVec.back())
>> 217 {
>> 218 eTdummyVec.push_back(tDummy);
>> 219 eVecm[tDummy].push_back(0.);
>> 220 }
>> 221
>> 222 eDiffCrossSection>>eDiffCrossSectionData[tDummy][eDummy];
>> 223
>> 224 if (eDummy != eVecm[tDummy].back()) eVecm[tDummy].push_back(eDummy);
>> 225 }
>> 226
>> 227 // End final state
>> 228
>> 229 #ifdef UEHARA_VERBOSE
>> 230 if (verboseLevel > 2)
>> 231 G4cout << "Loaded cross section files for CPA100 Elastic model" << G4endl;
>> 232 #endif
>> 233
>> 234 #ifdef UEHARA_VERBOSE
>> 235 if( verboseLevel>0 )
>> 236 {
>> 237 G4cout << "CPA100 Elastic model is initialized " << G4endl
>> 238 << "Energy range: "
>> 239 << LowEnergyLimit() / eV << " eV - "
>> 240 << HighEnergyLimit() / keV << " keV"
>> 241 << G4endl;
163 } 242 }
>> 243 #endif
>> 244
>> 245 // Initialize water density pointer
>> 246 fpMolWaterDensity = G4DNAMolecularMaterial::Instance()
>> 247 ->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
164 248
>> 249 if (isInitialised) { return; }
165 fParticleChangeForGamma = GetParticleChangeF 250 fParticleChangeForGamma = GetParticleChangeForGamma();
166 isInitialised = true; 251 isInitialised = true;
>> 252
167 } 253 }
168 254
169 //....oooOO0OOooo........oooOO0OOooo........oo 255 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
170 256
171 G4double G4DNACPA100ElasticModel::CrossSection << 257 G4double G4DNACPA100ElasticModel::CrossSectionPerVolume
172 << 258 (const G4Material* material,
173 << 259 const G4ParticleDefinition* p,
>> 260 G4double ekin,
>> 261 G4double,
>> 262 G4double)
174 { 263 {
175 // Get the name of the current particle <<
176 const G4String& particleName = p->GetParticl <<
177 auto materialID = pMaterial->GetIndex(); <<
178 264
179 // set killBelowEnergy value for current mat << 265 #ifdef UEHARA_VERBOSE
180 fKillBelowEnergy = fpModelData->GetLowELimit << 266 if (verboseLevel > 3)
>> 267 G4cout <<
>> 268 "Calling CrossSectionPerVolume() of G4DNACPA100ElasticModel" << G4endl;
>> 269 #endif
181 270
182 G4double sigma = 0.; << 271 // Calculate total cross section for model
183 272
184 if (ekin < fpModelData->GetHighELimit(materi << 273 G4double sigma=0;
185 if (ekin < fKillBelowEnergy) { <<
186 return DBL_MAX; <<
187 } <<
188 274
189 auto tableData = fpModelData->GetData(); << 275 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
190 276
191 if ((*tableData)[materialID][p] == nullptr << 277 const G4String& particleName = p->GetParticleName();
192 G4Exception("G4DNACPA100ElasticModel::Cr <<
193 "No model is registered"); <<
194 } <<
195 sigma = (*tableData)[materialID][p]->FindV <<
196 } <<
197 278
198 if (verboseLevel > 2) { << 279 if (ekin <= HighEnergyLimit() && ekin >= LowEnergyLimit())
199 auto MolDensity = << 280 {
200 (*G4DNAMolecularMaterial::Instance()->Ge << 281 //SI : XS must not be zero otherwise sampling of secondaries
>> 282 // method ignored
>> 283
>> 284 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
>> 285 pos = tableData.find(particleName);
>> 286
>> 287 if (pos != tableData.end())
>> 288 {
>> 289 G4DNACrossSectionDataSet* table = pos->second;
>> 290 if (table != 0)
>> 291 {
>> 292 sigma = table->FindValue(ekin);
>> 293 //
>> 294 //Dump in non-MT mode
>> 295 //
>> 296 /*
>> 297 G4double minEnergy = 10.481 * eV;
>> 298 G4double maxEnergy = 255955. * eV;
>> 299 G4int nEnergySteps = 1000;
>> 300 G4double energy(minEnergy);
>> 301 G4double
>> 302 stpEnergy(std::pow(maxEnergy/energy,
>> 303 1./static_cast<G4double>(nEnergySteps-1)));
>> 304 G4int step(nEnergySteps);
>> 305 system ("rm -rf elastic-cpa100.out");
>> 306 FILE* myFile=fopen("elastic-cpa100.out","a");
>> 307 while (step>0)
>> 308 {
>> 309 step--;
>> 310 fprintf (myFile,"%16.9le %16.9le\n",
>> 311 energy/eV,
>> 312 table->FindValue(energy)/(1e-20*m*m));
>> 313 energy*=stpEnergy;
>> 314 }
>> 315 fclose (myFile);
>> 316 abort();
>> 317 */
>> 318 //
>> 319 // end of dump
>> 320 //
>> 321 }
>> 322 }
>> 323 else
>> 324 {
>> 325 G4Exception("G4DNACPA100ElasticModel::ComputeCrossSectionPerVolume",
>> 326 "em0002",
>> 327 FatalException,"Model not applicable to particle type.");
>> 328 }
>> 329 }
>> 330
>> 331 #ifdef UEHARA_VERBOSE
>> 332 if (verboseLevel > 2)
>> 333 {
201 G4cout << "_______________________________ 334 G4cout << "__________________________________" << G4endl;
202 G4cout << "°°° G4DNACPA100ElasticModel << 335 G4cout << "G4DNACPA100ElasticModel - XS INFO START" << G4endl;
203 G4cout << "°°° Kinetic energy(eV)=" << << 336 G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
204 G4cout << "°°° lowLim (eV) = " << GetLo << 337 G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
205 << " highLim (eV) : " << GetHighELi << 338 G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
206 G4cout << "°°° Materials = " << (*G4Mat << 339 // G4cout << " - Cross section per water molecule (cm^-1)="
207 << G4endl; << 340 // << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
208 G4cout << "°°° Cross section per molecu << 341 G4cout << "G4DNACPA100ElasticModel - XS INFO END" << G4endl;
209 G4cout << "°°° Cross section per Phosph << 342 }
210 << G4endl; << 343 #endif
211 G4cout << "°°° G4DNACPA100ElasticModel << 344
212 } << 345 return sigma*waterDensity;
213 <<
214 auto MolDensity = <<
215 (*G4DNAMolecularMaterial::Instance()->GetN <<
216 return sigma * MolDensity; <<
217 } 346 }
218 347
219 //....oooOO0OOooo........oooOO0OOooo........oo 348 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
220 349
221 void G4DNACPA100ElasticModel::SampleSecondarie 350 void G4DNACPA100ElasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
222 << 351 const G4MaterialCutsCouple* /*couple*/,
223 << 352 const G4DynamicParticle* aDynamicElectron,
224 << 353 G4double,
225 { << 354 G4double)
226 G4double electronEnergy0 = aDynamicElectron- << 355 {
227 auto materialID = couple->GetMaterial()->Get << 356 #ifdef UEHARA_VERBOSE
228 auto p = aDynamicElectron->GetParticleDefini << 357 if (verboseLevel > 3)
>> 358 G4cout << "Calling SampleSecondaries() of G4DNACPA100ElasticModel" << G4endl;
>> 359 #endif
229 360
230 if (p != fpParticle) { << 361 G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
231 G4Exception("G4DNACPA100ElasticModel::Samp << 362
232 "This particle is not applied << 363 G4double cosTheta = RandomizeCosTheta(electronEnergy0);
233 } << 364 G4double phi = 2. * pi * G4UniformRand();
234 if (electronEnergy0 < fKillBelowEnergy) { << 365
235 return; << 366 G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
236 } << 367
237 G4double cosTheta = fpModelData->RandomizeCo << 368 //G4ThreeVector xVers = zVers.orthogonal();
238 G4double phi = 2. * CLHEP::pi * G4UniformRan << 369 //G4ThreeVector yVers = zVers.cross(xVers);
>> 370 //G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
>> 371 //G4double yDir = xDir;
>> 372 //xDir *= std::cos(phi);
>> 373 //yDir *= std::sin(phi);
239 374
240 const G4ThreeVector& zVers = aDynamicElectro << 375 // Computation of scattering angles (from Subroutine DIRAN in CPA100)
241 376
242 G4double CT1, ST1, CF1, SF1, CT2, ST2, CF2, 377 G4double CT1, ST1, CF1, SF1, CT2, ST2, CF2, SF2;
243 G4double sinTheta = std::sqrt(1 - cosTheta * << 378 G4double sinTheta = std::sqrt (1-cosTheta*cosTheta);
>> 379
>> 380 CT1=0;
>> 381 ST1=0;
>> 382 CF1=0;
>> 383 SF1=0;
>> 384 CT2=0;
>> 385 ST2=0;
>> 386 CF2=0;
>> 387 SF2=0;
244 388
245 CT1 = zVers.z(); 389 CT1 = zVers.z();
246 ST1 = std::sqrt(1. - CT1 * CT1); << 390 ST1=std::sqrt(1.-CT1*CT1);
247 391
248 if (ST1 != 0) << 392 if (ST1!=0) CF1 = zVers.x()/ST1; else CF1 = std::cos(2. * pi * G4UniformRand());
249 CF1 = zVers.x() / ST1; << 393 if (ST1!=0) SF1 = zVers.y()/ST1; else SF1 = std::sqrt(1.-CF1*CF1);
250 else <<
251 CF1 = std::cos(2. * CLHEP::pi * G4UniformR <<
252 if (ST1 != 0) <<
253 SF1 = zVers.y() / ST1; <<
254 else <<
255 SF1 = std::sqrt(1. - CF1 * CF1); <<
256 394
257 G4double A3, A4, A5, A2, A1; 395 G4double A3, A4, A5, A2, A1;
>> 396 A3=0;
>> 397 A4=0;
>> 398 A5=0;
>> 399 A2=0;
>> 400 A1=0;
258 401
259 A3 = sinTheta * std::cos(phi); << 402 A3 = sinTheta*std::cos(phi);
260 A4 = A3 * CT1 + ST1 * cosTheta; << 403 A4 = A3*CT1 + ST1*cosTheta;
261 A5 = sinTheta * std::sin(phi); 404 A5 = sinTheta * std::sin(phi);
262 A2 = A4 * SF1 + A5 * CF1; 405 A2 = A4 * SF1 + A5 * CF1;
263 A1 = A4 * CF1 - A5 * SF1; 406 A1 = A4 * CF1 - A5 * SF1;
264 407
265 CT2 = CT1 * cosTheta - ST1 * A3; << 408 CT2 = CT1*cosTheta - ST1*A3;
266 ST2 = std::sqrt(1. - CT2 * CT2); << 409 ST2 = std::sqrt(1.-CT2*CT2);
267 410
268 if (ST2 == 0) ST2 = 1E-6; << 411 if (ST2==0) ST2=1E-6;
269 CF2 = A1 / ST2; << 412 CF2 = A1/ST2;
270 SF2 = A2 / ST2; << 413 SF2 = A2/ST2;
271 G4ThreeVector zPrimeVers(ST2 * CF2, ST2 * SF <<
272 <<
273 fParticleChangeForGamma->ProposeMomentumDire <<
274 <<
275 auto EnergyDeposit = fpModelData->GetElastic <<
276 fParticleChangeForGamma->ProposeLocalEnergyD <<
277 if (statCode) { <<
278 fParticleChangeForGamma->SetProposedKineti <<
279 } <<
280 else { <<
281 auto newEnergy = electronEnergy0 - EnergyD <<
282 fParticleChangeForGamma->SetProposedKineti <<
283 } <<
284 } <<
285 414
286 G4double G4DNACPA100ElasticModel::Theta(const << 415 /*
287 G4doub << 416 G4cout << "CT1=" << CT1 << G4endl;
288 { << 417 G4cout << "ST1=" << ST1 << G4endl;
289 G4double theta, valueT1, valueT2, valueE21, << 418 G4cout << "CF1=" << CF1 << G4endl;
290 G4double xs11 = 0; << 419 G4cout << "SF1=" << SF1 << G4endl;
291 G4double xs12 = 0; << 420 G4cout << "cosTheta=" << cosTheta << G4endl;
292 G4double xs21 = 0; << 421 G4cout << "sinTheta=" << sinTheta << G4endl;
293 G4double xs22 = 0; << 422 G4cout << "cosPhi=" << std::cos(phi) << G4endl;
294 if (p == G4Electron::ElectronDefinition()) { << 423 G4cout << "sinPhi=" << std::sin(phi) << G4endl;
295 if (k == tValuesVec[materialID][p].back()) << 424 G4cout << "CT2=" << CT2 << G4endl;
296 k = k * (1. - 1e-12); << 425 G4cout << "ST2=" << ST2 << G4endl;
297 } << 426 G4cout << "CF2=" << CF2 << G4endl;
298 auto t2 = << 427 G4cout << "SF2=" << SF2 << G4endl;
299 std::upper_bound(tValuesVec[materialID][ << 428 */
300 auto t1 = t2 - 1; <<
301 <<
302 auto e12 = std::upper_bound(eValuesVect[ma <<
303 eValuesVect[ma <<
304 auto e11 = e12 - 1; <<
305 <<
306 auto e22 = std::upper_bound(eValuesVect[ma <<
307 eValuesVect[ma <<
308 auto e21 = e22 - 1; <<
309 <<
310 valueT1 = *t1; <<
311 valueT2 = *t2; <<
312 valueE21 = *e21; <<
313 valueE22 = *e22; <<
314 valueE12 = *e12; <<
315 valueE11 = *e11; <<
316 <<
317 xs11 = diffCrossSectionData[materialID][p] <<
318 xs12 = diffCrossSectionData[materialID][p] <<
319 xs21 = diffCrossSectionData[materialID][p] <<
320 xs22 = diffCrossSectionData[materialID][p] <<
321 } <<
322 429
323 if (xs11 == 0 && xs12 == 0 && xs21 == 0 && x << 430 G4ThreeVector zPrimeVers(ST2*CF2,ST2*SF2,CT2);
324 return (0.); << 431
325 } << 432 //
>> 433
>> 434 fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit()) ;
326 435
327 theta = QuadInterpolator(valueE11, valueE12, << 436 if (!statCode)
328 valueT2, k, integrD << 437
>> 438 fParticleChangeForGamma->SetProposedKineticEnergy
>> 439 (electronEnergy0-1.214E-4*(1.-cosTheta)*electronEnergy0);
>> 440
>> 441 else fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
329 442
>> 443 //
>> 444
>> 445 fParticleChangeForGamma->ProposeLocalEnergyDeposit(1.214E-4*(1.-cosTheta)*electronEnergy0);
>> 446
>> 447 }
>> 448
>> 449 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 450
>> 451 G4double G4DNACPA100ElasticModel::Theta
>> 452 (G4ParticleDefinition *, G4double k, G4double integrDiff)
>> 453 {
>> 454
>> 455 G4double theta = 0.;
>> 456 G4double valueT1 = 0;
>> 457 G4double valueT2 = 0;
>> 458 G4double valueE21 = 0;
>> 459 G4double valueE22 = 0;
>> 460 G4double valueE12 = 0;
>> 461 G4double valueE11 = 0;
>> 462 G4double xs11 = 0;
>> 463 G4double xs12 = 0;
>> 464 G4double xs21 = 0;
>> 465 G4double xs22 = 0;
>> 466
>> 467 // Protection against out of boundary access
>> 468 if (k==eTdummyVec.back()) k=k*(1.-1e-12);
>> 469 //
>> 470
>> 471 std::vector<G4double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
>> 472 std::vector<G4double>::iterator t1 = t2-1;
>> 473
>> 474 std::vector<G4double>::iterator e12 = std::upper_bound(eVecm[(*t1)].begin(),eVecm[(*t1)].end(),
>> 475 integrDiff);
>> 476 std::vector<G4double>::iterator e11 = e12-1;
>> 477
>> 478 std::vector<G4double>::iterator e22 = std::upper_bound(eVecm[(*t2)].begin(),eVecm[(*t2)].end(),
>> 479 integrDiff);
>> 480 std::vector<G4double>::iterator e21 = e22-1;
>> 481
>> 482 valueT1 =*t1;
>> 483 valueT2 =*t2;
>> 484 valueE21 =*e21;
>> 485 valueE22 =*e22;
>> 486 valueE12 =*e12;
>> 487 valueE11 =*e11;
>> 488
>> 489 xs11 = eDiffCrossSectionData[valueT1][valueE11];
>> 490 xs12 = eDiffCrossSectionData[valueT1][valueE12];
>> 491 xs21 = eDiffCrossSectionData[valueT2][valueE21];
>> 492 xs22 = eDiffCrossSectionData[valueT2][valueE22];
>> 493
>> 494 //TEST CPA100
>> 495 //if(k==valueT1) xs22 = eDiffCrossSectionData[valueT1][valueE12];
>> 496
>> 497 if (xs11==0 && xs12==0 && xs21==0 && xs22==0) return (0.);
>> 498
>> 499 theta = QuadInterpolator(
>> 500 valueE11, valueE12,
>> 501 valueE21, valueE22,
>> 502 xs11, xs12,
>> 503 xs21, xs22,
>> 504 valueT1, valueT2,
>> 505 k, integrDiff);
>> 506
330 return theta; 507 return theta;
>> 508
>> 509 //TEST CPA100
>> 510 //return xs22;
331 } 511 }
332 512
333 //....oooOO0OOooo........oooOO0OOooo........oo 513 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
334 514
335 G4double G4DNACPA100ElasticModel::LinLogInterp << 515 G4double G4DNACPA100ElasticModel::LinLogInterpolate(G4double e1,
336 << 516 G4double e2,
>> 517 G4double e,
>> 518 G4double xs1,
>> 519 G4double xs2)
337 { 520 {
338 G4double d1 = std::log(xs1); 521 G4double d1 = std::log(xs1);
339 G4double d2 = std::log(xs2); 522 G4double d2 = std::log(xs2);
340 G4double value = std::exp(d1 + (d2 - d1) * ( << 523 G4double value = std::exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
341 return value; 524 return value;
342 } 525 }
343 526
344 //....oooOO0OOooo........oooOO0OOooo........oo 527 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
345 528
346 G4double G4DNACPA100ElasticModel::LinLinInterp << 529 G4double G4DNACPA100ElasticModel::LinLinInterpolate(G4double e1,
347 << 530 G4double e2,
>> 531 G4double e,
>> 532 G4double xs1,
>> 533 G4double xs2)
348 { 534 {
349 G4double d1 = xs1; 535 G4double d1 = xs1;
350 G4double d2 = xs2; 536 G4double d2 = xs2;
351 G4double value = (d1 + (d2 - d1) * (e - e1) << 537 G4double value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
352 return value; 538 return value;
353 } 539 }
354 540
355 //....oooOO0OOooo........oooOO0OOooo........oo 541 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
356 542
357 G4double G4DNACPA100ElasticModel::LogLogInterp << 543 G4double G4DNACPA100ElasticModel::LogLogInterpolate(G4double e1,
358 << 544 G4double e2,
359 { << 545 G4double e,
360 G4double a = (std::log10(xs2) - std::log10(x << 546 G4double xs1,
361 G4double b = std::log10(xs2) - a * std::log1 << 547 G4double xs2)
362 G4double sigma = a * std::log10(e) + b; << 548 {
363 G4double value = (std::pow(10., sigma)); << 549 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
>> 550 G4double b = std::log10(xs2) - a*std::log10(e2);
>> 551 G4double sigma = a*std::log10(e) + b;
>> 552 G4double value = (std::pow(10.,sigma));
364 return value; 553 return value;
365 } 554 }
366 555
367 //....oooOO0OOooo........oooOO0OOooo........oo 556 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
368 557
369 G4double G4DNACPA100ElasticModel::QuadInterpol << 558
370 << 559 G4double G4DNACPA100ElasticModel::QuadInterpolator(G4double e11, G4double e12,
371 << 560 G4double e21, G4double e22,
372 << 561 G4double xs11, G4double xs12,
>> 562 G4double xs21, G4double xs22,
>> 563 G4double t1, G4double t2,
>> 564 G4double t, G4double e)
373 { 565 {
374 // Log-Log 566 // Log-Log
375 /* << 567 /*
376 G4double interpolatedvalue1 = LogLogInterp << 568 G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
377 G4double interpolatedvalue2 = LogLogInterp << 569 G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
378 G4double value = LogLogInterpolate(t1, t2, << 570 G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
379 571
380 572
381 // Lin-Log << 573 // Lin-Log
382 G4double interpolatedvalue1 = LinLogInterp << 574 G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12);
383 G4double interpolatedvalue2 = LinLogInterp << 575 G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22);
384 G4double value = LinLogInterpolate(t1, t2, << 576 G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
385 */ << 577 */
386 578
387 // Lin-Lin 579 // Lin-Lin
388 G4double interpolatedvalue1 = LinLinInterpol 580 G4double interpolatedvalue1 = LinLinInterpolate(e11, e12, e, xs11, xs12);
389 G4double interpolatedvalue2 = LinLinInterpol 581 G4double interpolatedvalue2 = LinLinInterpolate(e21, e22, e, xs21, xs22);
390 G4double value = LinLinInterpolate(t1, t2, t 582 G4double value = LinLinInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
391 583
392 return value; 584 return value;
393 } 585 }
394 586
395 //....oooOO0OOooo........oooOO0OOooo........oo 587 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
396 588
397 G4double G4DNACPA100ElasticModel::RandomizeCos << 589 G4double G4DNACPA100ElasticModel::RandomizeCosTheta(G4double k)
398 { 590 {
399 G4double integrdiff = 0; // PROBABILITY bet << 591
400 G4double uniformRand = G4UniformRand(); << 592 G4double integrdiff=0; // PROBABILITY between 0 and 1.
>> 593 G4double uniformRand=G4UniformRand();
401 integrdiff = uniformRand; 594 integrdiff = uniformRand;
402 G4double cosTheta = 0.; << 595
403 cosTheta = 1 - Theta(G4Electron::ElectronDef << 596 G4double cosTheta=0.;
404 // cosTheta = std::cos(theta * CLHEP::pi / 1 <<
405 return cosTheta; <<
406 } <<
407 //....oooOO0OOooo........oooOO0OOooo........oo <<
408 <<
409 void G4DNACPA100ElasticModel::ReadDiffCSFile(c <<
410 c <<
411 c <<
412 { <<
413 const char* path = G4FindDataDir("G4LEDATA") <<
414 if (path == nullptr) { <<
415 G4Exception("G4DNACPA100ElasticModel::Read <<
416 "G4LEDATA environment variable <<
417 return; <<
418 } <<
419 <<
420 std::ostringstream fullFileName; <<
421 fullFileName << path << "/" << file << ".dat <<
422 <<
423 std::ifstream diffCrossSection(fullFileName. <<
424 // error if file is not there <<
425 std::stringstream endPath; <<
426 if (!diffCrossSection) { <<
427 endPath << "Missing data file: " << file; <<
428 G4Exception("G4DNACPA100ElasticModel::Init <<
429 endPath.str().c_str()); <<
430 } <<
431 <<
432 tValuesVec[materialName][particleName].push_ <<
433 597
434 G4String line; << 598 // 1 - COS THETA is read from the data file
435 while (std::getline(diffCrossSection, line)) << 599 cosTheta = 1 - Theta(G4Electron::ElectronDefinition(),k/eV,integrdiff);
436 // <<
437 std::istringstream testIss(line); <<
438 G4String test; <<
439 testIss >> test; <<
440 if (test == "#") { <<
441 continue; <<
442 } <<
443 // check if line is empty <<
444 if (line.empty()) { <<
445 continue; <<
446 } <<
447 std::istringstream iss(line); <<
448 <<
449 G4double tDummy; <<
450 G4double eDummy; <<
451 600
452 iss >> tDummy >> eDummy; << 601 //
453 << 602 //
454 if (tDummy != tValuesVec[materialName][par << 603 //Dump
455 // Add the current T value << 604 //
456 tValuesVec[materialName][particleName].p << 605 //G4cout << "theta=" << theta << G4endl;
457 // Make it correspond to a default zero << 606 //G4cout << "cos theta=" << std::cos(theta*pi/180) << G4endl;
458 eValuesVect[materialName][particleName][ << 607 //G4cout << "sin theta=" << std::sin(theta*pi/180) << G4endl;
459 } << 608 //G4cout << "acos(cos theta)=" << std::acos(cosTheta) << G4endl;
460 iss >> diffCrossSectionData[materialName][ << 609 //G4cout << "cos theta="<< cosTheta << G4endl;
>> 610 //G4cout << "1 - cos theta="<< 1. - cosTheta << G4endl;
>> 611 //G4cout << "sin theta=" << std::sqrt(1-cosTheta*cosTheta) << G4endl;
>> 612 //
>> 613 /*
>> 614 G4double minProb = 0; // we scan probability between 0 and one
>> 615 G4double maxProb = 1;
>> 616 G4int nProbSteps = 100;
>> 617 G4double prob(minProb);
>> 618 G4double stepProb((maxProb-minProb)/static_cast<G4double>(nProbSteps));
>> 619 G4int step(nProbSteps);
>> 620 system ("rm -rf elastic-cumul-cpa100-100keV.out");
>> 621 FILE* myFile=fopen("elastic-cumul-cpa100-100keV.out","a");
>> 622 while (step>=0)
>> 623 {
>> 624 step--;
>> 625 fprintf (myFile,"%16.9le %16.9le\n",
>> 626 prob,
>> 627 Theta(G4Electron::ElectronDefinition(),100000,prob)); // SELECT NRJ IN eV !!!
>> 628 prob=prob+stepProb;
>> 629 }
>> 630 fclose (myFile);
>> 631 abort();
>> 632 */
>> 633 //
>> 634 // end of dump
>> 635 //
461 636
462 if (eDummy != eValuesVect[materialName][pa << 637 return cosTheta;
463 eValuesVect[materialName][particleName][ <<
464 } <<
465 } <<
466 } 638 }
467 639