| 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 // Author: Luciano Pandola and Gianfranco Pate << 26 // $Id: G4PenelopeRayleighModelMI.hh 75573 2013-11-04 11:48:15Z gcosmo $
>> 27 //
>> 28 // Author: Luciano Pandola and Gianfranco Paternò
27 // 29 //
28 // ------------------------------------------- 30 // -------------------------------------------------------------------
29 // History: 31 // History:
30 // 03 Dec 2009 L. Pandola 1st implementati 32 // 03 Dec 2009 L. Pandola 1st implementation
31 // 25 May 2011 L. Pandola Renamed (make v2 33 // 25 May 2011 L. Pandola Renamed (make v2008 as default Penelope)
32 // 27 Sep 2013 L. Pandola Migration to MT 34 // 27 Sep 2013 L. Pandola Migration to MT paradigm
33 // 20 Aug 2017 G. Paterno Molecular Interf << 35 // 20 Aug 2017 G. Paternò Molecular Interference implementation
34 // 24 Mar 2019 G. Paterno Improved Molecul << 36 // 24 Mar 2019 G. Paternò Improved Molecular Interference implementation
35 // 20 Jun 2020 G. Paterno Read qext separa << 37 // 20 Jun 2020 G. Paternò Read qext separately and leave original atomic form factors
36 // 27 Aug 2020 G. Paterno Further improvem << 38 // 27 Aug 2020 G. Paternò Further improvement of MI implementation
37 // 39 //
38 // ------------------------------------------- 40 // -------------------------------------------------------------------
39 // Class description: 41 // Class description:
40 // Low Energy Electromagnetic Physics, Rayleig 42 // Low Energy Electromagnetic Physics, Rayleigh Scattering
41 // with the model from Penelope, version 2008 43 // with the model from Penelope, version 2008
42 // ------------------------------------------- 44 // -------------------------------------------------------------------
43 // 45 //
44 //....oooOO0OOooo........oooOO0OOooo........oo 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45 47
46 #include "G4PenelopeRayleighModelMI.hh" 48 #include "G4PenelopeRayleighModelMI.hh"
47 49
48 #include "G4PenelopeSamplingData.hh" 50 #include "G4PenelopeSamplingData.hh"
49 #include "G4ParticleDefinition.hh" 51 #include "G4ParticleDefinition.hh"
50 #include "G4MaterialCutsCouple.hh" 52 #include "G4MaterialCutsCouple.hh"
51 #include "G4ProductionCutsTable.hh" 53 #include "G4ProductionCutsTable.hh"
52 #include "G4DynamicParticle.hh" 54 #include "G4DynamicParticle.hh"
53 #include "G4ElementTable.hh" 55 #include "G4ElementTable.hh"
54 #include "G4Element.hh" 56 #include "G4Element.hh"
55 #include "G4PhysicsFreeVector.hh" 57 #include "G4PhysicsFreeVector.hh"
56 #include "G4AutoLock.hh" 58 #include "G4AutoLock.hh"
57 #include "G4Exp.hh" 59 #include "G4Exp.hh"
58 #include "G4ExtendedMaterial.hh" 60 #include "G4ExtendedMaterial.hh"
59 #include "G4CrystalExtension.hh" 61 #include "G4CrystalExtension.hh"
60 #include "G4EmParameters.hh" 62 #include "G4EmParameters.hh"
61 63
62 #include "G4PhysicalConstants.hh" 64 #include "G4PhysicalConstants.hh"
63 #include "G4SystemOfUnits.hh" 65 #include "G4SystemOfUnits.hh"
64 //....oooOO0OOooo........oooOO0OOooo........oo 66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
65 67
66 const G4int G4PenelopeRayleighModelMI::fMaxZ; <<
67 G4PhysicsFreeVector* G4PenelopeRayleighModelMI <<
68 G4PhysicsFreeVector* G4PenelopeRayleighModelMI <<
69 <<
70 //....oooOO0OOooo........oooOO0OOooo........oo <<
71 <<
72 G4PenelopeRayleighModelMI::G4PenelopeRayleighM 68 G4PenelopeRayleighModelMI::G4PenelopeRayleighModelMI(const G4ParticleDefinition* part,
73 const G4String& nam) : 69 const G4String& nam) :
74 G4VEmModel(nam), 70 G4VEmModel(nam),
75 fParticleChange(nullptr),fParticle(nullptr), << 71 fParticleChange(nullptr),fParticle(nullptr),isInitialised(false),logAtomicCrossSection(nullptr),
76 fSamplingTable(nullptr),fMolInterferenceData << 72 atomicFormFactor(nullptr),MolInterferenceData(nullptr),logFormFactorTable(nullptr),
77 fIsInitialised(false),fLocalTable(false),fIs << 73 pMaxTable(nullptr),samplingTable(nullptr),fLocalTable(false),fIsMIActive(true),
>> 74 angularFunction(nullptr), fKnownMaterials(nullptr)
78 { 75 {
79 fIntrinsicLowEnergyLimit = 100.0*eV; 76 fIntrinsicLowEnergyLimit = 100.0*eV;
80 fIntrinsicHighEnergyLimit = 100.0*GeV; 77 fIntrinsicHighEnergyLimit = 100.0*GeV;
81 //SetLowEnergyLimit(fIntrinsicLowEnergyLimit 78 //SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
82 SetHighEnergyLimit(fIntrinsicHighEnergyLimit 79 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
83 80
84 if (part) SetParticle(part); 81 if (part) SetParticle(part);
85 82
86 fVerboseLevel = 0; << 83 verboseLevel = 0;
87 // Verbosity scale: 84 // Verbosity scale:
88 // 0 = nothing 85 // 0 = nothing
89 // 1 = warning for energy non-conservation 86 // 1 = warning for energy non-conservation
90 // 2 = details of energy budget 87 // 2 = details of energy budget
91 // 3 = calculation of FF and CS, file openin 88 // 3 = calculation of FF and CS, file openings, sampling of atoms
92 // 4 = entering in methods 89 // 4 = entering in methods
93 90
94 //build the energy grid. It is the same for 91 //build the energy grid. It is the same for all materials
95 G4double logenergy = G4Log(fIntrinsicLowEner 92 G4double logenergy = G4Log(fIntrinsicLowEnergyLimit/2.);
96 G4double logmaxenergy = G4Log(1.5*fIntrinsic 93 G4double logmaxenergy = G4Log(1.5*fIntrinsicHighEnergyLimit);
97 //finer grid below 160 keV 94 //finer grid below 160 keV
98 G4double logtransitionenergy = G4Log(160*keV 95 G4double logtransitionenergy = G4Log(160*keV);
99 G4double logfactor1 = G4Log(10.)/250.; 96 G4double logfactor1 = G4Log(10.)/250.;
100 G4double logfactor2 = logfactor1*10; 97 G4double logfactor2 = logfactor1*10;
101 fLogEnergyGridPMax.push_back(logenergy); << 98 logEnergyGridPMax.push_back(logenergy);
102 do { 99 do {
103 if (logenergy < logtransitionenergy) 100 if (logenergy < logtransitionenergy)
104 logenergy += logfactor1; 101 logenergy += logfactor1;
105 else 102 else
106 logenergy += logfactor2; 103 logenergy += logfactor2;
107 fLogEnergyGridPMax.push_back(logenergy); << 104 logEnergyGridPMax.push_back(logenergy);
108 } while (logenergy < logmaxenergy); 105 } while (logenergy < logmaxenergy);
109 } 106 }
110 107
111 //....oooOO0OOooo........oooOO0OOooo........oo 108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
112 109
113 G4PenelopeRayleighModelMI::~G4PenelopeRayleigh 110 G4PenelopeRayleighModelMI::~G4PenelopeRayleighModelMI()
114 { 111 {
115 if (IsMaster() || fLocalTable) { 112 if (IsMaster() || fLocalTable) {
116 113
117 for(G4int i=0; i<=fMaxZ; ++i) << 114 if (logAtomicCrossSection) {
118 { << 115 for (auto& item : (*logAtomicCrossSection))
119 if(fLogAtomicCrossSection[i]) <<
120 { <<
121 delete fLogAtomicCrossSection[i]; <<
122 fLogAtomicCrossSection[i] = nullptr; <<
123 } <<
124 if(fAtomicFormFactor[i]) <<
125 { <<
126 delete fAtomicFormFactor[i]; <<
127 fAtomicFormFactor[i] = nullptr; <<
128 } <<
129 } <<
130 if (fMolInterferenceData) { <<
131 for (auto& item : (*fMolInterferenceData <<
132 if (item.second) delete item.second; 116 if (item.second) delete item.second;
133 delete fMolInterferenceData; << 117 delete logAtomicCrossSection;
134 fMolInterferenceData = nullptr; << 118 logAtomicCrossSection = nullptr;
135 } << 119 }
>> 120
>> 121 if (atomicFormFactor) {
>> 122 for (auto& item : (*atomicFormFactor))
>> 123 if (item.second) delete item.second;
>> 124 delete atomicFormFactor;
>> 125 atomicFormFactor = nullptr;
>> 126 }
>> 127
>> 128 if (MolInterferenceData) {
>> 129 for (auto& item : (*MolInterferenceData))
>> 130 if (item.second) delete item.second;
>> 131 delete MolInterferenceData;
>> 132 MolInterferenceData = nullptr;
>> 133 }
>> 134
136 if (fKnownMaterials) 135 if (fKnownMaterials)
137 { 136 {
138 delete fKnownMaterials; 137 delete fKnownMaterials;
139 fKnownMaterials = nullptr; 138 fKnownMaterials = nullptr;
140 } 139 }
141 if (fAngularFunction) << 140
>> 141 if (angularFunction)
142 { 142 {
143 delete fAngularFunction; << 143 delete angularFunction;
144 fAngularFunction = nullptr; << 144 angularFunction = nullptr;
145 } 145 }
146 ClearTables(); << 146
>> 147 ClearTables();
>> 148
147 } 149 }
148 } 150 }
149 151
150 //....oooOO0OOooo........oooOO0OOooo........oo 152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
151 153
152 void G4PenelopeRayleighModelMI::ClearTables() 154 void G4PenelopeRayleighModelMI::ClearTables()
153 { 155 {
154 if (fLogFormFactorTable) { << 156 /*
155 for (auto& item : (*fLogFormFactorTable)) << 157 if (!IsMaster())
>> 158 //Should not be here!
>> 159 G4Exception("G4PenelopeRayleighModelMI::ClearTables()",
>> 160 "em0100",FatalException,"Worker thread in this method");
>> 161 */
>> 162
>> 163 if (logFormFactorTable) {
>> 164 for (auto& item : (*logFormFactorTable))
156 if (item.second) delete item.second; 165 if (item.second) delete item.second;
157 delete fLogFormFactorTable; << 166 delete logFormFactorTable;
158 fLogFormFactorTable = nullptr; //zero expl << 167 logFormFactorTable = nullptr; //zero explicitly
159 } 168 }
160 169
161 if (fPMaxTable) { << 170 if (pMaxTable) {
162 for (auto& item : (*fPMaxTable)) << 171 for (auto& item : (*pMaxTable))
163 if (item.second) delete item.second; 172 if (item.second) delete item.second;
164 delete fPMaxTable; << 173 delete pMaxTable;
165 fPMaxTable = nullptr; //zero explicitly << 174 pMaxTable = nullptr; //zero explicitly
166 } 175 }
167 176
168 if (fSamplingTable) { << 177 if (samplingTable) {
169 for (auto& item : (*fSamplingTable)) << 178 for (auto& item : (*samplingTable))
170 if (item.second) delete item.second; 179 if (item.second) delete item.second;
171 delete fSamplingTable; << 180 delete samplingTable;
172 fSamplingTable = nullptr; //zero explicitl << 181 samplingTable = nullptr; //zero explicitly
173 } 182 }
174 183
175 return; 184 return;
176 } 185 }
177 186
178 //....oooOO0OOooo........oooOO0OOooo........oo 187 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
179 188
180 void G4PenelopeRayleighModelMI::Initialise(con 189 void G4PenelopeRayleighModelMI::Initialise(const G4ParticleDefinition* part,
181 const G4DataVector& ) 190 const G4DataVector& )
182 { 191 {
183 if (fVerboseLevel > 3) << 192 if (verboseLevel > 3)
184 G4cout << "Calling G4PenelopeRayleighModel 193 G4cout << "Calling G4PenelopeRayleighModelMI::Initialise()" << G4endl;
185 194
186 SetParticle(part); 195 SetParticle(part);
187 196
188 if (fVerboseLevel) << 197
>> 198 if (verboseLevel)
189 G4cout << "# Molecular Interference is " < 199 G4cout << "# Molecular Interference is " << (fIsMIActive ? "ON" : "OFF") << " #" << G4endl;
190 200
191 //Only the master model creates/fills/destro 201 //Only the master model creates/fills/destroys the tables
192 if (IsMaster() && part == fParticle) { 202 if (IsMaster() && part == fParticle) {
193 //clear tables depending on materials, not 203 //clear tables depending on materials, not the atomic ones
194 ClearTables(); 204 ClearTables();
195 205
196 //Use here the highest verbosity, from G4E 206 //Use here the highest verbosity, from G4EmParameter or local
197 G4int globVerb = G4EmParameters::Instance( 207 G4int globVerb = G4EmParameters::Instance()->Verbose();
198 if (globVerb > fVerboseLevel) << 208 if (globVerb > verboseLevel)
199 { 209 {
200 fVerboseLevel = globVerb; << 210 verboseLevel = globVerb;
201 if (fVerboseLevel) << 211 if (verboseLevel)
202 G4cout << "Verbosity level of G4PenelopeRa << 212 G4cout << "Verbosity level of G4PenelopeRayleighModelMI set to " << verboseLevel <<
203 " from G4EmParameters()" << G4endl; 213 " from G4EmParameters()" << G4endl;
204 } 214 }
205 if (fVerboseLevel > 3) << 215
>> 216 if (verboseLevel > 3)
206 G4cout << "Calling G4PenelopeRayleighMod 217 G4cout << "Calling G4PenelopeRayleighModelMI::Initialise() [master]" << G4endl;
207 218
208 //Load the list of known materials and the 219 //Load the list of known materials and the DCS integration grid
209 if (fIsMIActive) 220 if (fIsMIActive)
210 { 221 {
211 if (!fKnownMaterials) 222 if (!fKnownMaterials)
212 fKnownMaterials = new std::map<G4String,G4 223 fKnownMaterials = new std::map<G4String,G4String>;
213 if (!fKnownMaterials->size()) 224 if (!fKnownMaterials->size())
214 LoadKnownMIFFMaterials(); 225 LoadKnownMIFFMaterials();
215 if (!fAngularFunction) << 226 if (!angularFunction)
216 { 227 {
217 //Create and fill once 228 //Create and fill once
218 fAngularFunction = new G4PhysicsFreeVect << 229 angularFunction = new G4PhysicsFreeVector(Ntheta);
219 CalculateThetaAndAngFun(); //angular fun 230 CalculateThetaAndAngFun(); //angular funtion for DCS integration
220 } 231 }
221 } 232 }
222 if (fIsMIActive && !fMolInterferenceData) << 233
223 fMolInterferenceData = new std::map<G4St << 234 //create new tables
224 if (!fLogFormFactorTable) << 235 //logAtomicCrossSection and atomicFormFactor are created only once,
225 fLogFormFactorTable = new std::map<const << 236 //since they are never cleared
226 if (!fPMaxTable) << 237 if (!logAtomicCrossSection)
227 fPMaxTable = new std::map<const G4Materi << 238 logAtomicCrossSection = new std::map<G4int,G4PhysicsFreeVector*>;
228 if (!fSamplingTable) << 239
229 fSamplingTable = new std::map<const G4Ma << 240 if (!atomicFormFactor)
>> 241 atomicFormFactor = new std::map<G4int,G4PhysicsFreeVector*>;
>> 242
>> 243 if (fIsMIActive && !MolInterferenceData)
>> 244 MolInterferenceData = new std::map<G4String,G4PhysicsFreeVector*>; //G. Paternò
>> 245
>> 246 if (!logFormFactorTable)
>> 247 logFormFactorTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
>> 248
>> 249 if (!pMaxTable)
>> 250 pMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
>> 251
>> 252 if (!samplingTable)
>> 253 samplingTable = new std::map<const G4Material*,G4PenelopeSamplingData*>;
230 254
231 //loop on the used materials 255 //loop on the used materials
232 G4ProductionCutsTable* theCoupleTable = G4 256 G4ProductionCutsTable* theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable();
233 257
234 for (G4int i=0;i<(G4int)theCoupleTable->Ge << 258 for (size_t i=0;i<theCoupleTable->GetTableSize();i++) {
235 const G4Material* material = 259 const G4Material* material =
236 theCoupleTable->GetMaterialCutsCouple(i)->Ge 260 theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
237 const G4ElementVector* theElementVector 261 const G4ElementVector* theElementVector = material->GetElementVector();
238 262
239 for (std::size_t j=0;j<material->GetNumb << 263 for (size_t j=0;j<material->GetNumberOfElements();j++) {
240 G4int iZ = theElementVector->at(j)->GetZasIn 264 G4int iZ = theElementVector->at(j)->GetZasInt();
241 //read data files only in the master 265 //read data files only in the master
242 if (!fLogAtomicCrossSection[iZ]) << 266 if (!logAtomicCrossSection->count(iZ))
243 ReadDataFile(iZ); 267 ReadDataFile(iZ);
244 } 268 }
245 269
246 //1) Read MI form factors 270 //1) Read MI form factors
247 if (fIsMIActive && !fMolInterferenceData << 271 if (fIsMIActive && !MolInterferenceData->count(material->GetName()))
248 ReadMolInterferenceData(material->GetName()) 272 ReadMolInterferenceData(material->GetName());
249 273
250 //2) If the table has not been built for 274 //2) If the table has not been built for the material, do it!
251 if (!fLogFormFactorTable->count(material << 275 if (!logFormFactorTable->count(material))
252 BuildFormFactorTable(material); 276 BuildFormFactorTable(material);
253 277
254 //3) retrieve or build the sampling tabl 278 //3) retrieve or build the sampling table
255 if (!(fSamplingTable->count(material))) << 279 if (!(samplingTable->count(material)))
256 InitializeSamplingAlgorithm(material); 280 InitializeSamplingAlgorithm(material);
257 281
258 //4) retrieve or build the pMax data 282 //4) retrieve or build the pMax data
259 if (!fPMaxTable->count(material)) << 283 if (!pMaxTable->count(material))
260 GetPMaxTable(material); 284 GetPMaxTable(material);
261 } 285 }
262 286
263 if (fVerboseLevel > 1) { << 287 if (verboseLevel > 1) {
264 G4cout << G4endl << "Penelope Rayleigh m 288 G4cout << G4endl << "Penelope Rayleigh model v2008 is initialized" << G4endl
265 << "Energy range: " 289 << "Energy range: "
266 << LowEnergyLimit() / keV << " keV - " 290 << LowEnergyLimit() / keV << " keV - "
267 << HighEnergyLimit() / GeV << " GeV" 291 << HighEnergyLimit() / GeV << " GeV"
268 << G4endl; 292 << G4endl;
269 } << 293 }
>> 294
270 } 295 }
271 296
272 if (fIsInitialised) << 297 if (isInitialised)
273 return; 298 return;
274 fParticleChange = GetParticleChangeForGamma( 299 fParticleChange = GetParticleChangeForGamma();
275 fIsInitialised = true; << 300 isInitialised = true;
276 } 301 }
277 302
278 //....oooOO0OOooo........oooOO0OOooo........oo 303 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
279 304
280 void G4PenelopeRayleighModelMI::InitialiseLoca 305 void G4PenelopeRayleighModelMI::InitialiseLocal(const G4ParticleDefinition* part,
281 G4VEmModel *masterModel) << 306 G4VEmModel *masterModel)
282 { 307 {
283 if (fVerboseLevel > 3) << 308 if (verboseLevel > 3)
284 G4cout << "Calling G4PenelopeRayleighMode 309 G4cout << "Calling G4PenelopeRayleighModelMI::InitialiseLocal()" << G4endl;
285 310
286 //Check that particle matches: one might hav 311 //Check that particle matches: one might have multiple master models
287 //(e.g. for e+ and e-) 312 //(e.g. for e+ and e-)
288 if (part == fParticle) { 313 if (part == fParticle) {
289 314
290 //Get the const table pointers from the ma 315 //Get the const table pointers from the master to the workers
291 const G4PenelopeRayleighModelMI* theModel 316 const G4PenelopeRayleighModelMI* theModel =
292 static_cast<G4PenelopeRayleighModelMI*> 317 static_cast<G4PenelopeRayleighModelMI*> (masterModel);
293 318
294 //Copy pointers to the data tables 319 //Copy pointers to the data tables
295 for(G4int i=0; i<=fMaxZ; ++i) << 320 logAtomicCrossSection = theModel->logAtomicCrossSection;
296 { << 321 atomicFormFactor = theModel->atomicFormFactor;
297 fLogAtomicCrossSection[i] = theModel->fLogAt << 322 MolInterferenceData = theModel->MolInterferenceData;
298 fAtomicFormFactor[i] = theModel->fAtomicForm << 323 logFormFactorTable = theModel->logFormFactorTable;
299 } << 324 pMaxTable = theModel->pMaxTable;
300 fMolInterferenceData = theModel->fMolInter << 325 samplingTable = theModel->samplingTable;
301 fLogFormFactorTable = theModel->fLogFormFa <<
302 fPMaxTable = theModel->fPMaxTable; <<
303 fSamplingTable = theModel->fSamplingTable; <<
304 fKnownMaterials = theModel->fKnownMaterial 326 fKnownMaterials = theModel->fKnownMaterials;
305 fAngularFunction = theModel->fAngularFunct << 327 angularFunction = theModel->angularFunction;
306 328
307 //Copy the G4DataVector with the grid 329 //Copy the G4DataVector with the grid
308 fLogQSquareGrid = theModel->fLogQSquareGri << 330 logQSquareGrid = theModel->logQSquareGrid;
309 331
310 //Same verbosity for all workers, as the m 332 //Same verbosity for all workers, as the master
311 fVerboseLevel = theModel->fVerboseLevel; << 333 verboseLevel = theModel->verboseLevel;
312 } << 334
>> 335 }
>> 336
313 return; 337 return;
314 } 338 }
315 339
316 //....oooOO0OOooo........oooOO0OOooo........oo 340 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
317 341
318 namespace {G4Mutex PenelopeRayleighModelMutex 342 namespace {G4Mutex PenelopeRayleighModelMutex = G4MUTEX_INITIALIZER;}
319 G4double G4PenelopeRayleighModelMI::ComputeCro 343 G4double G4PenelopeRayleighModelMI::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
320 G4double energy, << 344 G4double energy,
321 G4double Z, << 345 G4double Z,
322 G4double, << 346 G4double,
323 G4double, << 347 G4double,
324 G4double) << 348 G4double)
325 { 349 {
326 //Cross section of Rayleigh scattering in Pe << 350 //Cross section of Rayleigh scattering in Penelope v2008 is calculated by the EPDL97
327 //tabulation, Cuellen et al. (1997), with no << 351 //tabulation, Cuellen et al. (1997), with non-relativistic form factors from Hubbel
328 //et al. J. Phys. Chem. Ref. Data 4 (1975) 4 << 352 //et al. J. Phys. Chem. Ref. Data 4 (1975) 471; Erratum ibid. 6 (1977) 615.
329 << 353
330 if (fVerboseLevel > 3) << 354 if (verboseLevel > 3)
331 G4cout << "Calling CrossSectionPerAtom() o << 355 G4cout << "Calling CrossSectionPerAtom() of G4PenelopeRayleighModelMI" << G4endl;
332 << 356
333 G4int iZ = G4int(Z); << 357 G4int iZ = (G4int) Z;
334 if (!fLogAtomicCrossSection[iZ]) { << 358
335 //If we are here, it means that Initialize << 359 //Either Initialize() was not called, or we are in a slave and InitializeLocal() was
336 //not filled up. This can happen in a Unit << 360 //not invoked
337 if (fVerboseLevel > 0) { << 361 if (!logAtomicCrossSection) {
338 //Issue a G4Exception (warning) only in << 362 //create a **thread-local** version of the table. Used only for G4EmCalculator and
339 G4ExceptionDescription ed; << 363 //Unit Tests
340 ed << "Unable to retrieve the cross sect << 364 fLocalTable = true;
341 ed << "This can happen only in Unit Test << 365 logAtomicCrossSection = new std::map<G4int,G4PhysicsFreeVector*>;
342 G4Exception("G4PenelopeRayleighModelMI:: << 366 }
343 "em2040",JustWarning,ed); <<
344 } <<
345 367
346 //protect file reading via autolock << 368 //now it should be ok
347 G4AutoLock lock(&PenelopeRayleighModelMute << 369 if (!logAtomicCrossSection->count(iZ)) {
348 ReadDataFile(iZ); << 370 //If we are here, it means that Initialize() was inkoved, but the MaterialTable was
349 lock.unlock(); << 371 //not filled up. This can happen in a UnitTest or via G4EmCalculator
350 } << 372 if (verboseLevel > 0) {
>> 373 //Issue a G4Exception (warning) only in verbose mode
>> 374 G4ExceptionDescription ed;
>> 375 ed << "Unable to retrieve the cross section table for Z=" << iZ << G4endl;
>> 376 ed << "This can happen only in Unit Tests or via G4EmCalculator" << G4endl;
>> 377 G4Exception("G4PenelopeRayleighModelMI::ComputeCrossSectionPerAtom()",
>> 378 "em2040",JustWarning,ed);
>> 379 }
351 380
352 G4double cross = 0; << 381 //protect file reading via autolock
353 G4PhysicsFreeVector* atom = fLogAtomicCrossS << 382 G4AutoLock lock(&PenelopeRayleighModelMutex);
354 if (!atom) { << 383 ReadDataFile(iZ);
355 G4ExceptionDescription ed; << 384 lock.unlock();
356 ed << "Unable to find Z=" << iZ << " in th << 385
357 G4Exception("G4PenelopeRayleighModelMI::Co << 386 }
358 "em2041",FatalException,ed); <<
359 return 0; <<
360 } <<
361 387
362 G4double logene = G4Log(energy); << 388 G4double cross = 0;
363 G4double logXS = atom->Value(logene); << 389
364 cross = G4Exp(logXS); << 390 G4PhysicsFreeVector* atom = logAtomicCrossSection->find(iZ)->second;
365 << 391 if (!atom) {
366 if (fVerboseLevel > 2) { << 392 G4ExceptionDescription ed;
367 G4cout << "Rayleigh cross section at " << << 393 ed << "Unable to find Z=" << iZ << " in the atomic cross section table" << G4endl;
368 << Z << " = " << cross/barn << " barn" << << 394 G4Exception("G4PenelopeRayleighModelMI::ComputeCrossSectionPerAtom()",
369 } << 395 "em2041",FatalException,ed);
370 return cross; << 396 return 0;
>> 397 }
>> 398
>> 399 G4double logene = G4Log(energy);
>> 400 G4double logXS = atom->Value(logene);
>> 401 cross = G4Exp(logXS);
>> 402
>> 403 if (verboseLevel > 2) {
>> 404 G4cout << "Rayleigh cross section at " << energy/keV << " keV for Z="
>> 405 << Z << " = " << cross/barn << " barn" << G4endl;
>> 406 }
>> 407
>> 408 return cross;
371 } 409 }
372 410
373 //....oooOO0OOooo........oooOO0OOooo........oo 411 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
374 412
375 void G4PenelopeRayleighModelMI::CalculateTheta 413 void G4PenelopeRayleighModelMI::CalculateThetaAndAngFun()
376 { 414 {
377 G4double theta = 0; 415 G4double theta = 0;
378 for(G4int k=0; k<fNtheta; k++) { << 416
>> 417 for(G4int k=0; k<Ntheta; k++) {
379 theta += fDTheta; 418 theta += fDTheta;
380 G4double value = (1+std::cos(theta)*std::c 419 G4double value = (1+std::cos(theta)*std::cos(theta))*std::sin(theta);
381 fAngularFunction->PutValue(k,theta,value); << 420 angularFunction->PutValue(k,theta,value);
382 if (fVerboseLevel > 3) << 421 if (verboseLevel > 3)
383 G4cout << "theta[" << k << "]: " << fAn << 422 G4cout << "theta[" << k << "]: " << angularFunction->Energy(k)
384 << ", angFun[" << k << "]: " << (*fAngu << 423 << ", angFun[" << k << "]: " << (*angularFunction)[k] << G4endl;
385 } 424 }
386 } 425 }
387 426
388 //....oooOO0OOooo........oooOO0OOooo........oo 427 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
389 428
390 G4double G4PenelopeRayleighModelMI::CalculateQ 429 G4double G4PenelopeRayleighModelMI::CalculateQSquared(G4double angle, G4double energy)
391 { 430 {
392 G4double lambda,x,q,q2 = 0; << 431 G4double lambda,x,q,q2 = 0;
393 432
394 lambda = hbarc*twopi/energy; << 433 lambda = hbarc*twopi/energy;
395 x = 1./lambda*std::sin(angle/2.); << 434 x = 1./lambda*std::sin(angle/2.);
396 q = 2.*h_Planck*x/(electron_mass_c2/c_light) << 435 q = 2.*h_Planck*x/(electron_mass_c2/c_light);
397 436
398 if (fVerboseLevel > 3) { << 437 if (verboseLevel > 3) {
399 G4cout << "E: " << energy/keV << " keV, la << 438 G4cout << "E: " << energy/keV << " keV, lambda: " << lambda/nm << " nm"
400 << ", x: " << x*nm << ", q: " << q << G4e << 439 << ", x: " << x*nm << ", q: " << q << G4endl;
401 } << 440 }
402 q2 = std::pow(q,2); << 441
403 return q2; << 442 q2 = std::pow(q,2);
>> 443
>> 444 return q2;
404 } 445 }
405 446
406 //....oooOO0OOooo........oooOO0OOooo........oo 447 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
407 448
408 //Overriding of parent's (G4VEmModel) method 449 //Overriding of parent's (G4VEmModel) method
409 G4double G4PenelopeRayleighModelMI::CrossSecti 450 G4double G4PenelopeRayleighModelMI::CrossSectionPerVolume(const G4Material* material,
410 const G4ParticleDefinition* p, 451 const G4ParticleDefinition* p,
411 G4double energy, 452 G4double energy,
412 G4double, 453 G4double,
413 G4double) 454 G4double)
414 { 455 {
415 //check if we are in a Unit Test (only for t 456 //check if we are in a Unit Test (only for the first time)
416 static G4bool amInAUnitTest = false; 457 static G4bool amInAUnitTest = false;
417 if (G4ProductionCutsTable::GetProductionCuts 458 if (G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize() == 0 && !amInAUnitTest)
418 { 459 {
419 amInAUnitTest = true; 460 amInAUnitTest = true;
420 G4ExceptionDescription ed; 461 G4ExceptionDescription ed;
421 ed << "The ProductionCuts table is empty 462 ed << "The ProductionCuts table is empty " << G4endl;
422 ed << "This should happen only in Unit T 463 ed << "This should happen only in Unit Tests" << G4endl;
423 G4Exception("G4PenelopeRayleighModelMI:: 464 G4Exception("G4PenelopeRayleighModelMI::CrossSectionPerVolume()",
424 "em2019",JustWarning,ed); 465 "em2019",JustWarning,ed);
425 } 466 }
426 //If the material does not have a MIFF, cont 467 //If the material does not have a MIFF, continue with the old-style calculation
427 G4String matname = material->GetName(); 468 G4String matname = material->GetName();
428 if (amInAUnitTest) 469 if (amInAUnitTest)
429 { 470 {
430 //No need to lock, as this is always cal 471 //No need to lock, as this is always called in a master
431 const G4ElementVector* theElementVector 472 const G4ElementVector* theElementVector = material->GetElementVector();
432 //protect file reading via autolock 473 //protect file reading via autolock
433 for (std::size_t j=0;j<material->GetNumb << 474 for (size_t j=0;j<material->GetNumberOfElements();j++) {
434 G4int iZ = theElementVector->at(j)->GetZasIn 475 G4int iZ = theElementVector->at(j)->GetZasInt();
435 if (!fLogAtomicCrossSection[iZ]) { << 476 if (!logAtomicCrossSection->count(iZ)) {
436 ReadDataFile(iZ); 477 ReadDataFile(iZ);
437 } 478 }
438 } 479 }
439 if (fIsMIActive) 480 if (fIsMIActive)
440 ReadMolInterferenceData(matname); 481 ReadMolInterferenceData(matname);
441 if (!fLogFormFactorTable->count(material << 482 if (!logFormFactorTable->count(material))
442 BuildFormFactorTable(material); 483 BuildFormFactorTable(material);
443 if (!(fSamplingTable->count(material))) << 484 if (!(samplingTable->count(material)))
444 InitializeSamplingAlgorithm(material); 485 InitializeSamplingAlgorithm(material);
445 if (!fPMaxTable->count(material)) << 486 if (!pMaxTable->count(material))
446 GetPMaxTable(material); 487 GetPMaxTable(material);
447 } 488 }
448 G4bool useMIFF = fIsMIActive && (fMolInterfe << 489
>> 490 G4bool useMIFF = fIsMIActive && (MolInterferenceData->count(matname) || matname.find("MedMat") != std::string::npos);
449 if (!useMIFF) 491 if (!useMIFF)
450 { 492 {
451 if (fVerboseLevel > 2) << 493 if (verboseLevel > 2)
452 G4cout << "Rayleigh CS of: " << matname << " 494 G4cout << "Rayleigh CS of: " << matname << " calculated through CSperAtom!" << G4endl;
453 return G4VEmModel::CrossSectionPerVolume 495 return G4VEmModel::CrossSectionPerVolume(material,p,energy);
454 } 496 }
455 497
456 // If we are here, it means that we have to 498 // If we are here, it means that we have to integrate the cross section
457 if (fVerboseLevel > 2) << 499 if (verboseLevel > 2)
458 G4cout << "Rayleigh CS of: " << matname 500 G4cout << "Rayleigh CS of: " << matname
459 << " calculated through integration of th 501 << " calculated through integration of the DCS" << G4endl;
460 << 502
>> 503
461 G4double cs = 0; 504 G4double cs = 0;
462 505
463 //force null cross-section if below the low- 506 //force null cross-section if below the low-energy edge of the table
464 if (energy < LowEnergyLimit()) << 507 if (energy < LowEnergyLimit()) return cs;
465 return cs; << 508
466 509
>> 510
>> 511
467 //if the material is a CRYSTAL, forbid this 512 //if the material is a CRYSTAL, forbid this process
468 if (material->IsExtended() && material->GetN 513 if (material->IsExtended() && material->GetName() != "CustomMat") {
469 G4ExtendedMaterial* extendedmaterial = (G4 514 G4ExtendedMaterial* extendedmaterial = (G4ExtendedMaterial*)material;
470 G4CrystalExtension* crystalExtension = (G4 515 G4CrystalExtension* crystalExtension = (G4CrystalExtension*)extendedmaterial->RetrieveExtension("crystal");
471 if (crystalExtension != 0) { 516 if (crystalExtension != 0) {
472 G4cout << "The material has a crystallin 517 G4cout << "The material has a crystalline structure, a dedicated diffraction model is used!" << G4endl;
473 return 0; 518 return 0;
474 } 519 }
475 } 520 }
476 521
477 //GET MATERIAL INFORMATION << 522
>> 523 //GET MATERIAL INFORMATION
>> 524
478 G4double atomDensity = material->GetTotNbOfA 525 G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
479 std::size_t nElements = material->GetNumberO << 526 G4int nElements = material->GetNumberOfElements();
480 const G4ElementVector* elementVector = mater 527 const G4ElementVector* elementVector = material->GetElementVector();
481 const G4double* fractionVector = material->G 528 const G4double* fractionVector = material->GetFractionVector();
>> 529 //const G4double* theAtomNumDensityVector = material-> GetVecNbOfAtomsPerVolume();
482 530
483 //Stoichiometric factors 531 //Stoichiometric factors
484 std::vector<G4double> *StoichiometricFactors 532 std::vector<G4double> *StoichiometricFactors = new std::vector<G4double>;
485 for (std::size_t i=0;i<nElements;++i) { << 533 for (G4int i=0;i<nElements;i++) {
486 G4double fraction = fractionVector[i]; 534 G4double fraction = fractionVector[i];
487 G4double atomicWeigth = (*elementVector)[i 535 G4double atomicWeigth = (*elementVector)[i]->GetA()/(g/mole);
488 StoichiometricFactors->push_back(fraction/ 536 StoichiometricFactors->push_back(fraction/atomicWeigth);
489 } 537 }
490 G4double MaxStoichiometricFactor = 0.; 538 G4double MaxStoichiometricFactor = 0.;
491 for (std::size_t i=0;i<nElements;++i) { << 539 for (G4int i=0;i<nElements;i++) {
492 if ((*StoichiometricFactors)[i] > MaxStoic 540 if ((*StoichiometricFactors)[i] > MaxStoichiometricFactor)
493 MaxStoichiometricFactor = (*Stoichiometr 541 MaxStoichiometricFactor = (*StoichiometricFactors)[i];
494 } 542 }
495 for (std::size_t i=0;i<nElements;++i) { << 543 for (G4int i=0;i<nElements;i++) {
496 (*StoichiometricFactors)[i] /= MaxStoichi 544 (*StoichiometricFactors)[i] /= MaxStoichiometricFactor;
>> 545 //G4cout << "StoichiometricFactors[" << i+1 << "]: " << (*StoichiometricFactors)[i] << G4endl;
497 } 546 }
498 547
499 //Equivalent atoms per molecule 548 //Equivalent atoms per molecule
500 G4double atPerMol = 0.; << 549 G4double atPerMol = 0;
501 for (std::size_t i=0;i<nElements;++i) << 550 for (G4int i=0;i<nElements;i++)
502 atPerMol += (*StoichiometricFactors)[i]; 551 atPerMol += (*StoichiometricFactors)[i];
503 G4double moleculeDensity = 0.; 552 G4double moleculeDensity = 0.;
504 if (atPerMol != 0.) moleculeDensity = atomDe << 553 if (atPerMol) moleculeDensity = atomDensity/atPerMol;
505 554
506 if (fVerboseLevel > 2) << 555 if (verboseLevel > 2)
507 G4cout << "Material " << material->GetName 556 G4cout << "Material " << material->GetName() << " has " << atPerMol << " atoms "
508 << "per molecule and " << moleculeDensity 557 << "per molecule and " << moleculeDensity/(cm*cm*cm) << " molecule/cm3" << G4endl;
509 558
510 //Equivalent molecular weight (dimensionless 559 //Equivalent molecular weight (dimensionless)
511 G4double MolWeight = 0.; 560 G4double MolWeight = 0.;
512 for (std::size_t i=0;i<nElements;++i) << 561 for (G4int i=0;i<nElements;i++)
513 MolWeight += (*StoichiometricFactors)[i]*( 562 MolWeight += (*StoichiometricFactors)[i]*(*elementVector)[i]->GetA()/(g/mole);
514 563
515 if (fVerboseLevel > 2) << 564 if (verboseLevel > 2)
516 G4cout << "Molecular weight of " << matnam 565 G4cout << "Molecular weight of " << matname << ": " << MolWeight << " g/mol" << G4endl;
517 566
518 G4double IntegrandFun[fNtheta]; << 567 G4double IntegrandFun[Ntheta];
519 for (G4int k=0; k<fNtheta; k++) { << 568 for (G4int k=0; k<Ntheta; k++) {
520 G4double theta = fAngularFunction->Energy( << 569 G4double theta = angularFunction->Energy(k); //the x-value is called "Energy", but is an angle
521 G4double F2 = GetFSquared(material,Calcula 570 G4double F2 = GetFSquared(material,CalculateQSquared(theta,energy));
522 IntegrandFun[k] = (*fAngularFunction)[k]*F << 571 IntegrandFun[k] = (*angularFunction)[k]*F2;
523 } 572 }
524 573
525 G4double constant = pi*classic_electr_radius 574 G4double constant = pi*classic_electr_radius*classic_electr_radius;
526 cs = constant*IntegrateFun(IntegrandFun,fNth << 575 cs = constant*IntegrateFun(IntegrandFun,Ntheta,fDTheta);
527 576
528 //Now cs is the cross section per molecule, 577 //Now cs is the cross section per molecule, let's calculate the cross section per volume
529 G4double csvolume = cs*moleculeDensity; << 578 G4double csvolume = cs*moleculeDensity;
>> 579
530 580
531 //print CS and mfp 581 //print CS and mfp
532 if (fVerboseLevel > 2) << 582 if (verboseLevel > 2)
533 G4cout << "Rayleigh CS of " << matname << 583 G4cout << "Rayleigh CS of " << matname << " at " << energy/keV
534 << " keV: " << cs/barn << " barn" 584 << " keV: " << cs/barn << " barn"
535 << ", mean free path: " << 1./csvolume/mm 585 << ", mean free path: " << 1./csvolume/mm << " mm" << G4endl;
>> 586
536 587
537 delete StoichiometricFactors; << 588 delete StoichiometricFactors;
>> 589
538 //return CS 590 //return CS
539 return csvolume; 591 return csvolume;
540 } 592 }
541 593
542 //....oooOO0OOooo........oooOO0OOooo........oo 594 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
543 595
544 void G4PenelopeRayleighModelMI::BuildFormFacto 596 void G4PenelopeRayleighModelMI::BuildFormFactorTable(const G4Material* material)
545 { 597 {
546 if (fVerboseLevel > 3) << 598 if (verboseLevel > 3)
547 G4cout << "Calling BuildFormFactorTable() 599 G4cout << "Calling BuildFormFactorTable() of G4PenelopeRayleighModelMI" << G4endl;
548 600
>> 601
549 //GET MATERIAL INFORMATION 602 //GET MATERIAL INFORMATION
550 std::size_t nElements = material->GetNumberO << 603 G4int nElements = material->GetNumberOfElements();
551 const G4ElementVector* elementVector = mater 604 const G4ElementVector* elementVector = material->GetElementVector();
552 const G4double* fractionVector = material->G 605 const G4double* fractionVector = material->GetFractionVector();
553 606
554 //Stoichiometric factors 607 //Stoichiometric factors
555 std::vector<G4double> *StoichiometricFactors 608 std::vector<G4double> *StoichiometricFactors = new std::vector<G4double>;
556 for (std::size_t i=0;i<nElements;++i) { << 609 for (G4int i=0;i<nElements;i++) {
557 G4double fraction = fractionVector[i]; 610 G4double fraction = fractionVector[i];
558 G4double atomicWeigth = (*elementVector)[i 611 G4double atomicWeigth = (*elementVector)[i]->GetA()/(g/mole);
559 StoichiometricFactors->push_back(fraction/ 612 StoichiometricFactors->push_back(fraction/atomicWeigth);
560 } 613 }
561 G4double MaxStoichiometricFactor = 0.; 614 G4double MaxStoichiometricFactor = 0.;
562 for (std::size_t i=0;i<nElements;++i) { << 615 for (G4int i=0;i<nElements;i++) {
563 if ((*StoichiometricFactors)[i] > MaxStoic 616 if ((*StoichiometricFactors)[i] > MaxStoichiometricFactor)
564 MaxStoichiometricFactor = (*Stoichiometr 617 MaxStoichiometricFactor = (*StoichiometricFactors)[i];
565 } 618 }
566 if (MaxStoichiometricFactor<1e-16) { 619 if (MaxStoichiometricFactor<1e-16) {
567 G4ExceptionDescription ed; 620 G4ExceptionDescription ed;
568 ed << "Inconsistent data of atomic composi 621 ed << "Inconsistent data of atomic composition for " << material->GetName() << G4endl;
569 G4Exception("G4PenelopeRayleighModelMI::Bu 622 G4Exception("G4PenelopeRayleighModelMI::BuildFormFactorTable()",
570 "em2042",FatalException,ed); 623 "em2042",FatalException,ed);
571 } 624 }
572 for (std::size_t i=0;i<nElements;++i) << 625 for (G4int i=0;i<nElements;i++)
573 (*StoichiometricFactors)[i] /= MaxStoichi 626 (*StoichiometricFactors)[i] /= MaxStoichiometricFactor;
574 627
>> 628 //Equivalent atoms per molecule
>> 629 G4double atomsPerMolecule = 0;
>> 630 for (G4int i=0;i<nElements;i++)
>> 631 atomsPerMolecule += (*StoichiometricFactors)[i];
>> 632
575 //Equivalent molecular weight (dimensionless 633 //Equivalent molecular weight (dimensionless)
576 G4double MolWeight = 0.; 634 G4double MolWeight = 0.;
577 for (std::size_t i=0;i<nElements;++i) << 635 for (G4int i=0;i<nElements;i++)
578 MolWeight += (*StoichiometricFactors)[i]*( 636 MolWeight += (*StoichiometricFactors)[i]*(*elementVector)[i]->GetA()/(g/mole);
579 << 637
>> 638
580 //CREATE THE FORM FACTOR TABLE 639 //CREATE THE FORM FACTOR TABLE
581 //First, the form factors are retrieved [F/s 640 //First, the form factors are retrieved [F/sqrt(W)].
582 //Then, they are squared and multiplied for 641 //Then, they are squared and multiplied for MolWeight -> F2 [dimensionless].
583 //This makes difference for CS calculation, 642 //This makes difference for CS calculation, but not for theta sampling.
584 G4PhysicsFreeVector* theFFVec = new G4Physic << 643 G4PhysicsFreeVector* theFFVec = new G4PhysicsFreeVector(logQSquareGrid.size());
585 /*spline=*/true); << 644 theFFVec->SetSpline(true);
586 645
587 G4String matname = material->GetName(); 646 G4String matname = material->GetName();
588 G4String aFileNameFF = ""; 647 G4String aFileNameFF = "";
589 648
590 //retrieve MIdata (fFileNameFF) 649 //retrieve MIdata (fFileNameFF)
591 G4MIData* dataMI = GetMIData(material); 650 G4MIData* dataMI = GetMIData(material);
592 if (dataMI) 651 if (dataMI)
593 aFileNameFF = dataMI->GetFilenameFF(); 652 aFileNameFF = dataMI->GetFilenameFF();
594 653
595 //read the MIFF from a file passed by the us 654 //read the MIFF from a file passed by the user
596 if (fIsMIActive && aFileNameFF != "") { 655 if (fIsMIActive && aFileNameFF != "") {
597 if (fVerboseLevel) << 656 if (verboseLevel)
598 G4cout << "Read MIFF for " << matname << 657 G4cout << "Read MIFF for " << matname << " from custom file: " << aFileNameFF << G4endl;
599 658
600 ReadMolInterferenceData(matname,aFileNameF 659 ReadMolInterferenceData(matname,aFileNameFF);
601 G4PhysicsFreeVector* Fvector = fMolInterf << 660 G4PhysicsFreeVector* Fvector = MolInterferenceData->find(matname)->second;
>> 661
>> 662 for (size_t k=0;k<logQSquareGrid.size();k++) {
>> 663 G4double ff2 = 0;
>> 664
>> 665 G4double q = std::pow(G4Exp(logQSquareGrid[k]),0.5);
>> 666 G4double f = Fvector->Value(q);
>> 667 //G4cout << "QGrid[" << k+1 << "]: " << q << ", F[" << k+1 << "]: " << f << G4endl;
>> 668
>> 669 ff2 = f*f*MolWeight;
>> 670 if (ff2) theFFVec->PutValue(k,logQSquareGrid[k],G4Log(ff2));
>> 671 }
602 672
603 for (std::size_t k=0;k<fLogQSquareGrid.siz <<
604 G4double q = std::pow(G4Exp(fLogQSquareG <<
605 G4double f = Fvector->Value(q); <<
606 G4double ff2 = f*f*MolWeight; <<
607 if (ff2) <<
608 theFFVec->PutValue(k,fLogQSquareGrid[k],G4Lo <<
609 } <<
610 } 673 }
611 //retrieve the MIFF from the database or use 674 //retrieve the MIFF from the database or use the IAM
612 else { 675 else {
613 //medical material: composition of fat, wa 676 //medical material: composition of fat, water, bonematrix, mineral
614 if (fIsMIActive && matname.find("MedMat") 677 if (fIsMIActive && matname.find("MedMat") != std::string::npos) {
615 if (fVerboseLevel) << 678 if (verboseLevel)
616 G4cout << "Build MIFF from components for " 679 G4cout << "Build MIFF from components for " << matname << G4endl;
617 680
618 //get the material composition from its 681 //get the material composition from its name
619 G4int ki, kf=6, ktot=19; 682 G4int ki, kf=6, ktot=19;
620 G4double comp[4]; 683 G4double comp[4];
621 G4String compstring = matname.substr(kf+ 684 G4String compstring = matname.substr(kf+1, ktot);
622 for (std::size_t j=0; j<4; ++j) { << 685 for (size_t j=0; j<4; j++) {
623 ki = kf+1; 686 ki = kf+1;
624 kf = ki+4; 687 kf = ki+4;
625 compstring = matname.substr(ki, 4); 688 compstring = matname.substr(ki, 4);
626 comp[j] = atof(compstring.c_str()); 689 comp[j] = atof(compstring.c_str());
627 if (fVerboseLevel > 2) << 690 if (verboseLevel > 2)
628 G4cout << " -- MedMat comp[" << j+1 << "]: 691 G4cout << " -- MedMat comp[" << j+1 << "]: " << comp[j] << G4endl;
629 } 692 }
630 693
631 const char* path = G4FindDataDir("G4LEDA << 694 char* path = std::getenv("G4LEDATA");
632 if (!path) { 695 if (!path) {
633 G4String excep = "G4LEDATA environment varia 696 G4String excep = "G4LEDATA environment variable not set!";
634 G4Exception("G4PenelopeRayleighModelMI::Buil 697 G4Exception("G4PenelopeRayleighModelMI::BuildFormFactorTable()",
635 "em0006",FatalException,excep); 698 "em0006",FatalException,excep);
636 } 699 }
637 700
638 if (!fMolInterferenceData->count("Fat_MI << 701 if (!MolInterferenceData->count("Fat_MI"))
639 ReadMolInterferenceData("Fat_MI"); 702 ReadMolInterferenceData("Fat_MI");
640 G4PhysicsFreeVector* fatFF = fMolInterfe << 703 G4PhysicsFreeVector* fatFF = MolInterferenceData->find("Fat_MI")->second;
641 704
642 if (!fMolInterferenceData->count("Water_ << 705 if (!MolInterferenceData->count("Water_MI"))
643 ReadMolInterferenceData("Water_MI"); 706 ReadMolInterferenceData("Water_MI");
644 G4PhysicsFreeVector* waterFF = fMolInter << 707 G4PhysicsFreeVector* waterFF = MolInterferenceData->find("Water_MI")->second;
645 708
646 if (!fMolInterferenceData->count("BoneMa << 709 if (!MolInterferenceData->count("BoneMatrix_MI"))
647 ReadMolInterferenceData("BoneMatrix_MI"); 710 ReadMolInterferenceData("BoneMatrix_MI");
648 G4PhysicsFreeVector* bonematrixFF = fMol << 711 G4PhysicsFreeVector* bonematrixFF = MolInterferenceData->find("BoneMatrix_MI")->second;
649 712
650 if (!fMolInterferenceData->count("Minera << 713 if (!MolInterferenceData->count("Mineral_MI"))
651 ReadMolInterferenceData("Mineral_MI"); 714 ReadMolInterferenceData("Mineral_MI");
652 G4PhysicsFreeVector* mineralFF = fMolInt << 715 G4PhysicsFreeVector* mineralFF = MolInterferenceData->find("Mineral_MI")->second;
653 716
654 //get and combine the molecular form fac 717 //get and combine the molecular form factors with interference effect
655 for (std::size_t k=0;k<fLogQSquareGrid.s << 718 for (size_t k=0;k<logQSquareGrid.size();k++) {
656 G4double ff2 = 0; 719 G4double ff2 = 0;
657 G4double q = std::pow(G4Exp(fLogQSquareGrid[ << 720 G4double q = std::pow(G4Exp(logQSquareGrid[k]),0.5);
658 G4double ffat = fatFF->Value(q); 721 G4double ffat = fatFF->Value(q);
659 G4double fwater = waterFF->Value(q); 722 G4double fwater = waterFF->Value(q);
660 G4double fbonematrix = bonematrixFF->Value(q 723 G4double fbonematrix = bonematrixFF->Value(q);
661 G4double fmineral = mineralFF->Value(q); 724 G4double fmineral = mineralFF->Value(q);
662 ff2 = comp[0]*ffat*ffat+comp[1]*fwater*fwate 725 ff2 = comp[0]*ffat*ffat+comp[1]*fwater*fwater+
663 comp[2]*fbonematrix*fbonematrix+comp[3]*fm 726 comp[2]*fbonematrix*fbonematrix+comp[3]*fmineral*fmineral;
664 ff2 *= MolWeight; 727 ff2 *= MolWeight;
665 if (ff2) theFFVec->PutValue(k,fLogQSquareGri << 728 if (ff2) theFFVec->PutValue(k,logQSquareGrid[k],G4Log(ff2));
666 } 729 }
667 } 730 }
668 //other materials with MIFF (from the data 731 //other materials with MIFF (from the database)
669 else if (fIsMIActive && fMolInterferenceDa << 732 else if (fIsMIActive && MolInterferenceData->count(matname)) {
670 if (fVerboseLevel) << 733 if (verboseLevel)
671 G4cout << "Read MIFF from database " << matn 734 G4cout << "Read MIFF from database " << matname << G4endl;
672 G4PhysicsFreeVector* FF = fMolInterferen << 735 G4PhysicsFreeVector* FF = MolInterferenceData->find(matname)->second;
673 for (std::size_t k=0;k<fLogQSquareGrid.s << 736 for (size_t k=0;k<logQSquareGrid.size();k++) {
674 G4double ff2 = 0; 737 G4double ff2 = 0;
675 G4double q = std::pow(G4Exp(fLogQSquareGrid[ << 738 G4double q = std::pow(G4Exp(logQSquareGrid[k]),0.5);
676 G4double f = FF->Value(q); 739 G4double f = FF->Value(q);
677 ff2 = f*f*MolWeight; 740 ff2 = f*f*MolWeight;
678 if (ff2) theFFVec->PutValue(k,fLogQSquareGri << 741 if (ff2) theFFVec->PutValue(k,logQSquareGrid[k],G4Log(ff2));
679 } << 742 }
680 } << 743
>> 744 }
>> 745
681 //IAM 746 //IAM
682 else { 747 else {
683 if (fVerboseLevel) << 748 if (verboseLevel)
684 G4cout << "FF of " << matname << " calculate 749 G4cout << "FF of " << matname << " calculated according to the IAM" << G4endl;
685 for (std::size_t k=0;k<fLogQSquareGrid.s << 750 for (size_t k=0;k<logQSquareGrid.size();k++) {
686 G4double ff2 = 0; 751 G4double ff2 = 0;
687 for (std::size_t i=0;i<nElements;++i) { << 752 for (G4int i=0;i<nElements;i++) {
688 G4int iZ = (*elementVector)[i]->GetZasInt( 753 G4int iZ = (*elementVector)[i]->GetZasInt();
689 G4PhysicsFreeVector* theAtomVec = fAtomicF << 754 G4PhysicsFreeVector* theAtomVec = atomicFormFactor->find(iZ)->second;
690 G4double q = std::pow(G4Exp(fLogQSquareGri << 755 G4double q = std::pow(G4Exp(logQSquareGrid[k]),0.5);
691 G4double f = theAtomVec->Value(q); 756 G4double f = theAtomVec->Value(q);
692 ff2 += f*f*(*StoichiometricFactors)[i]; 757 ff2 += f*f*(*StoichiometricFactors)[i];
693 } 758 }
694 if (ff2) theFFVec->PutValue(k,fLogQSquareGri << 759 if (ff2) theFFVec->PutValue(k,logQSquareGrid[k],G4Log(ff2));
695 } 760 }
696 } 761 }
697 } 762 }
698 theFFVec->FillSecondDerivatives(); << 763 logFormFactorTable->insert(std::make_pair(material,theFFVec));
699 fLogFormFactorTable->insert(std::make_pair(m <<
700 764
701 if (fVerboseLevel > 3) << 765 if (verboseLevel > 3)
702 DumpFormFactorTable(material); 766 DumpFormFactorTable(material);
703 delete StoichiometricFactors; 767 delete StoichiometricFactors;
704 768
705 return; 769 return;
>> 770
706 } 771 }
707 772
708 //....oooOO0OOooo........oooOO0OOooo........oo 773 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
709 774
710 void G4PenelopeRayleighModelMI::SampleSecondar 775 void G4PenelopeRayleighModelMI::SampleSecondaries(std::vector<G4DynamicParticle*>*,
711 const G4MaterialCutsCouple* coup 776 const G4MaterialCutsCouple* couple,
712 const G4DynamicParticle* aDynami 777 const G4DynamicParticle* aDynamicGamma,
713 G4double, 778 G4double,
714 G4double) 779 G4double)
715 { 780 {
716 // Sampling of the Rayleigh final state (nam << 781 // Sampling of the Rayleigh final state (namely, scattering angle of the photon)
717 // from the Penelope2008 model. The scatteri << 782 // from the Penelope2008 model. The scattering angle is sampled from the atomic
718 // cross section dOmega/d(cosTheta) from Bor << 783 // cross section dOmega/d(cosTheta) from Born ("Atomic Phyisics", 1969), disregarding
719 // anomalous scattering effects. The Form Fa << 784 // anomalous scattering effects. The Form Factor F(Q) function which appears in the
720 // analytical cross section is retrieved via << 785 // analytical cross section is retrieved via the method GetFSquared(); atomic data
721 // are tabulated for F(Q). Form factor for c << 786 // are tabulated for F(Q). Form factor for compounds is calculated according to
722 // the additivity rule. The sampling from th << 787 // the additivity rule. The sampling from the F(Q) is made via a Rational Inverse
723 // Transform with Aliasing (RITA) algorithm; << 788 // Transform with Aliasing (RITA) algorithm; RITA parameters are calculated once
724 // for each material and managed by G4Penelo << 789 // for each material and managed by G4PenelopeSamplingData objects.
725 // The sampling algorithm (rejection method) << 790 // The sampling algorithm (rejection method) has efficiency 67% at low energy, and
726 // increases with energy. For E=100 keV the << 791 // increases with energy. For E=100 keV the efficiency is 100% and 86% for
727 // hydrogen and uranium, respectively. << 792 // hydrogen and uranium, respectively.
728 << 793
729 if (fVerboseLevel > 3) << 794 if (verboseLevel > 3)
730 G4cout << "Calling SamplingSecondaries() o << 795 G4cout << "Calling SamplingSecondaries() of G4PenelopeRayleighModelMI" << G4endl;
731 << 796
732 G4double photonEnergy0 = aDynamicGamma->GetK << 797 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
733 << 798
734 if (photonEnergy0 <= fIntrinsicLowEnergyLimi << 799 if (photonEnergy0 <= fIntrinsicLowEnergyLimit) {
735 fParticleChange->ProposeTrackStatus(fStopA << 800 fParticleChange->ProposeTrackStatus(fStopAndKill);
736 fParticleChange->SetProposedKineticEnergy( << 801 fParticleChange->SetProposedKineticEnergy(0.);
737 fParticleChange->ProposeLocalEnergyDeposit << 802 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
738 return; << 803 return;
739 } << 804 }
740 805
741 G4ParticleMomentum photonDirection0 = aDynam << 806 G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
742 const G4Material* theMat = couple->GetMateri <<
743 807
744 //1) Verify if tables are ready << 808 const G4Material* theMat = couple->GetMaterial();
745 //Either Initialize() was not called, or we <<
746 //not invoked <<
747 if (!fPMaxTable || !fSamplingTable || !fLogF <<
748 //create a **thread-local** version of the <<
749 //Unit Tests <<
750 fLocalTable = true; <<
751 if (!fLogFormFactorTable) <<
752 fLogFormFactorTable = new std::map<const <<
753 if (!fPMaxTable) <<
754 fPMaxTable = new std::map<const G4Materi <<
755 if (!fSamplingTable) <<
756 fSamplingTable = new std::map<const G4Ma <<
757 if (fIsMIActive && !fMolInterferenceData) <<
758 fMolInterferenceData = new std::map<G4St <<
759 } <<
760 <<
761 if (!fSamplingTable->count(theMat)) { <<
762 //If we are here, it means that Initialize <<
763 //not filled up. This can happen in a Unit <<
764 if (fVerboseLevel > 0) { <<
765 //Issue a G4Exception (warning) only in <<
766 G4ExceptionDescription ed; <<
767 ed << "Unable to find the fSamplingTable <<
768 theMat->GetName() << G4endl; <<
769 ed << "This can happen only in Unit Test <<
770 G4Exception("G4PenelopeRayleighModelMI:: <<
771 "em2019",JustWarning,ed); <<
772 } <<
773 const G4ElementVector* theElementVector = <<
774 //protect file reading via autolock <<
775 G4AutoLock lock(&PenelopeRayleighModelMute <<
776 for (std::size_t j=0;j<theMat->GetNumberOf <<
777 G4int iZ = theElementVector->at(j)->GetZ <<
778 if (!fLogAtomicCrossSection[iZ]) { <<
779 lock.lock(); <<
780 ReadDataFile(iZ); <<
781 lock.unlock(); <<
782 } <<
783 } <<
784 lock.lock(); <<
785 809
786 //1) If the table has not been built for t << 810 //1) Verify if tables are ready
787 if (!fLogFormFactorTable->count(theMat)) << 811 //Either Initialize() was not called, or we are in a slave and InitializeLocal() was
788 BuildFormFactorTable(theMat); << 812 //not invoked
>> 813 if (!pMaxTable || !samplingTable || !logAtomicCrossSection || !atomicFormFactor ||
>> 814 !logFormFactorTable) {
>> 815 //create a **thread-local** version of the table. Used only for G4EmCalculator and
>> 816 //Unit Tests
>> 817 fLocalTable = true;
>> 818 if (!logAtomicCrossSection)
>> 819 logAtomicCrossSection = new std::map<G4int,G4PhysicsFreeVector*>;
>> 820 if (!atomicFormFactor)
>> 821 atomicFormFactor = new std::map<G4int,G4PhysicsFreeVector*>;
>> 822 if (!logFormFactorTable)
>> 823 logFormFactorTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
>> 824 if (!pMaxTable)
>> 825 pMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
>> 826 if (!samplingTable)
>> 827 samplingTable = new std::map<const G4Material*,G4PenelopeSamplingData*>;
>> 828 if (fIsMIActive && !MolInterferenceData)
>> 829 MolInterferenceData = new std::map<G4String,G4PhysicsFreeVector*>;
>> 830 }
>> 831
>> 832 if (!samplingTable->count(theMat)) {
>> 833 //If we are here, it means that Initialize() was inkoved, but the MaterialTable was
>> 834 //not filled up. This can happen in a UnitTest
>> 835 if (verboseLevel > 0) {
>> 836 //Issue a G4Exception (warning) only in verbose mode
>> 837 G4ExceptionDescription ed;
>> 838 ed << "Unable to find the samplingTable data for " <<
>> 839 theMat->GetName() << G4endl;
>> 840 ed << "This can happen only in Unit Tests" << G4endl;
>> 841 G4Exception("G4PenelopeRayleighModelMI::SampleSecondaries()",
>> 842 "em2019",JustWarning,ed);
>> 843 }
>> 844 const G4ElementVector* theElementVector = theMat->GetElementVector();
>> 845 //protect file reading via autolock
>> 846 G4AutoLock lock(&PenelopeRayleighModelMutex);
>> 847 for (size_t j=0;j<theMat->GetNumberOfElements();j++) {
>> 848 G4int iZ = theElementVector->at(j)->GetZasInt();
>> 849 if (!logAtomicCrossSection->count(iZ)) {
>> 850 lock.lock();
>> 851 ReadDataFile(iZ);
>> 852 lock.unlock();
>> 853 }
>> 854 }
>> 855 lock.lock();
>> 856
>> 857 //1) If the table has not been built for the material, do it!
>> 858 if (!logFormFactorTable->count(theMat))
>> 859 BuildFormFactorTable(theMat);
789 860
790 //2) retrieve or build the sampling table << 861 //2) retrieve or build the sampling table
791 if (!(fSamplingTable->count(theMat))) << 862 if (!(samplingTable->count(theMat)))
792 InitializeSamplingAlgorithm(theMat); << 863 InitializeSamplingAlgorithm(theMat);
793 864
794 //3) retrieve or build the pMax data << 865 //3) retrieve or build the pMax data
795 if (!fPMaxTable->count(theMat)) << 866 if (!pMaxTable->count(theMat))
796 GetPMaxTable(theMat); << 867 GetPMaxTable(theMat);
797 868
798 lock.unlock(); << 869 lock.unlock();
799 } << 870 }
800 871
801 //Ok, restart the job << 872 //Ok, restart the job
802 G4PenelopeSamplingData* theDataTable = fSamp << 873 G4PenelopeSamplingData* theDataTable = samplingTable->find(theMat)->second;
803 G4PhysicsFreeVector* thePMax = fPMaxTable->f << 874 G4PhysicsFreeVector* thePMax = pMaxTable->find(theMat)->second;
804 G4double cosTheta = 1.0; << 875
805 << 876 G4double cosTheta = 1.0;
806 //OK, ready to go! << 877
807 G4double qmax = 2.0*photonEnergy0/electron_m << 878 //OK, ready to go!
808 << 879 G4double qmax = 2.0*photonEnergy0/electron_mass_c2; //this is non-dimensional now
809 if (qmax < 1e-10) { //very low momentum tran << 880
810 G4bool loopAgain=false; << 881 if (qmax < 1e-10) { //very low momentum transfer
811 do { << 882 G4bool loopAgain=false;
812 loopAgain = false; << 883 do {
813 cosTheta = 1.0-2.0*G4UniformRand(); << 884 loopAgain = false;
814 G4double G = 0.5*(1+cosTheta*cosTheta); << 885 cosTheta = 1.0-2.0*G4UniformRand();
815 if (G4UniformRand()>G) << 886 G4double G = 0.5*(1+cosTheta*cosTheta);
816 loopAgain = true; << 887 if (G4UniformRand()>G)
817 } while(loopAgain); << 888 loopAgain = true;
818 } << 889 } while(loopAgain);
819 else { //larger momentum transfer << 890 }
820 std::size_t nData = theDataTable->GetNumbe <<
821 G4double LastQ2inTheTable = theDataTable-> <<
822 G4double q2max = std::min(qmax*qmax,LastQ2 <<
823 <<
824 G4bool loopAgain = false; <<
825 G4double MaxPValue = thePMax->Value(photon <<
826 G4double xx=0; <<
827 <<
828 //Sampling by rejection method. The reject <<
829 //G = 0.5*(1+cos^2(theta)) <<
830 do { <<
831 loopAgain = false; <<
832 G4double RandomMax = G4UniformRand()*Max <<
833 xx = theDataTable->SampleValue(RandomMax <<
834 <<
835 //xx is a random value of q^2 in (0,q2ma <<
836 //F(Q^2) via the RITA algorithm <<
837 if (xx > q2max) <<
838 loopAgain = true; <<
839 cosTheta = 1.0-2.0*xx/q2max; <<
840 G4double G = 0.5*(1+cosTheta*cosTheta); <<
841 if (G4UniformRand()>G) <<
842 loopAgain = true; <<
843 } while(loopAgain); <<
844 } <<
845 <<
846 G4double sinTheta = std::sqrt(1-cosTheta*cos <<
847 <<
848 //Scattered photon angles. ( Z - axis along <<
849 G4double phi = twopi * G4UniformRand() ; <<
850 G4double dirX = sinTheta*std::cos(phi); <<
851 G4double dirY = sinTheta*std::sin(phi); <<
852 G4double dirZ = cosTheta; <<
853 <<
854 //Update G4VParticleChange for the scattered <<
855 G4ThreeVector photonDirection1(dirX, dirY, d <<
856 <<
857 photonDirection1.rotateUz(photonDirection0); <<
858 fParticleChange->ProposeMomentumDirection(ph <<
859 fParticleChange->SetProposedKineticEnergy(ph <<
860 891
861 return; << 892 else { //larger momentum transfer
>> 893
>> 894 size_t nData = theDataTable->GetNumberOfStoredPoints();
>> 895 G4double LastQ2inTheTable = theDataTable->GetX(nData-1);
>> 896 G4double q2max = std::min(qmax*qmax,LastQ2inTheTable);
>> 897
>> 898 G4bool loopAgain = false;
>> 899 G4double MaxPValue = thePMax->Value(photonEnergy0);
>> 900 G4double xx=0;
>> 901
>> 902 //Sampling by rejection method. The rejection function is
>> 903 //G = 0.5*(1+cos^2(theta))
>> 904 do {
>> 905 loopAgain = false;
>> 906 G4double RandomMax = G4UniformRand()*MaxPValue;
>> 907 xx = theDataTable->SampleValue(RandomMax);
>> 908
>> 909 //xx is a random value of q^2 in (0,q2max),sampled according to
>> 910 //F(Q^2) via the RITA algorithm
>> 911 if (xx > q2max)
>> 912 loopAgain = true;
>> 913 cosTheta = 1.0-2.0*xx/q2max;
>> 914 G4double G = 0.5*(1+cosTheta*cosTheta);
>> 915 if (G4UniformRand()>G)
>> 916 loopAgain = true;
>> 917 } while(loopAgain);
>> 918 }
>> 919
>> 920 G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
>> 921
>> 922 //Scattered photon angles. ( Z - axis along the parent photon)
>> 923 G4double phi = twopi * G4UniformRand() ;
>> 924 G4double dirX = sinTheta*std::cos(phi);
>> 925 G4double dirY = sinTheta*std::sin(phi);
>> 926 G4double dirZ = cosTheta;
>> 927
>> 928 //Update G4VParticleChange for the scattered photon
>> 929 G4ThreeVector photonDirection1(dirX, dirY, dirZ);
>> 930
>> 931 photonDirection1.rotateUz(photonDirection0);
>> 932 fParticleChange->ProposeMomentumDirection(photonDirection1) ;
>> 933 fParticleChange->SetProposedKineticEnergy(photonEnergy0) ;
>> 934
>> 935 return;
862 } 936 }
863 937
864 //....oooOO0OOooo........oooOO0OOooo........oo 938 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
865 939
866 void G4PenelopeRayleighModelMI::ReadDataFile(c 940 void G4PenelopeRayleighModelMI::ReadDataFile(const G4int Z)
867 { 941 {
868 if (fVerboseLevel > 2) { << 942
>> 943 if (verboseLevel > 2) {
869 G4cout << "G4PenelopeRayleighModelMI::Read 944 G4cout << "G4PenelopeRayleighModelMI::ReadDataFile()" << G4endl;
870 G4cout << "Going to read Rayleigh data fil 945 G4cout << "Going to read Rayleigh data files for Z=" << Z << G4endl;
871 } 946 }
872 947
873 const char* path = G4FindDataDir("G4LEDATA") << 948 char* path = std::getenv("G4LEDATA");
874 if (!path) { 949 if (!path) {
875 G4String excep = "G4LEDATA environment var 950 G4String excep = "G4LEDATA environment variable not set!";
876 G4Exception("G4PenelopeRayleighModelMI::Re 951 G4Exception("G4PenelopeRayleighModelMI::ReadDataFile()",
877 "em0006",FatalException,excep); 952 "em0006",FatalException,excep);
878 return; 953 return;
879 } 954 }
880 955
881 //Read first the cross section file (all the 956 //Read first the cross section file (all the files have 250 points)
882 std::ostringstream ost; 957 std::ostringstream ost;
883 if (Z>9) 958 if (Z>9)
884 ost << path << "/penelope/rayleigh/pdgra" 959 ost << path << "/penelope/rayleigh/pdgra" << Z << ".p08";
885 else 960 else
886 ost << path << "/penelope/rayleigh/pdgra0" 961 ost << path << "/penelope/rayleigh/pdgra0" << Z << ".p08";
887 std::ifstream file(ost.str().c_str()); 962 std::ifstream file(ost.str().c_str());
888 963
889 if (!file.is_open()) { 964 if (!file.is_open()) {
890 G4String excep = "Data file " + G4String(o 965 G4String excep = "Data file " + G4String(ost.str()) + " not found!";
891 G4Exception("G4PenelopeRayleighModelMI::Re 966 G4Exception("G4PenelopeRayleighModelMI::ReadDataFile()",
892 "em0003",FatalException,excep); 967 "em0003",FatalException,excep);
893 } 968 }
894 969
895 G4int readZ = 0; 970 G4int readZ = 0;
896 std::size_t nPoints = 0; << 971 size_t nPoints = 0;
897 file >> readZ >> nPoints; 972 file >> readZ >> nPoints;
898 973
899 if (readZ != Z || nPoints <= 0 || nPoints >= 974 if (readZ != Z || nPoints <= 0 || nPoints >= 5000) {
900 G4ExceptionDescription ed; 975 G4ExceptionDescription ed;
901 ed << "Corrupted data file for Z=" << Z << 976 ed << "Corrupted data file for Z=" << Z << G4endl;
902 G4Exception("G4PenelopeRayleighModelMI::Re 977 G4Exception("G4PenelopeRayleighModelMI::ReadDataFile()",
903 "em0005",FatalException,ed); 978 "em0005",FatalException,ed);
904 return; 979 return;
905 } 980 }
906 981
907 fLogAtomicCrossSection[Z] = new G4PhysicsFre << 982 G4PhysicsFreeVector* theVec = new G4PhysicsFreeVector((size_t)nPoints);
908 G4double ene=0,f1=0,f2=0,xs=0; 983 G4double ene=0,f1=0,f2=0,xs=0;
909 for (std::size_t i=0;i<nPoints;++i) { << 984 for (size_t i=0;i<nPoints;i++) {
910 file >> ene >> f1 >> f2 >> xs; 985 file >> ene >> f1 >> f2 >> xs;
911 //dimensional quantities 986 //dimensional quantities
912 ene *= eV; 987 ene *= eV;
913 xs *= cm2; 988 xs *= cm2;
914 fLogAtomicCrossSection[Z]->PutValue(i,G4Lo << 989 theVec->PutValue(i,G4Log(ene),G4Log(xs));
915 if (file.eof() && i != (nPoints-1)) { //fi 990 if (file.eof() && i != (nPoints-1)) { //file ended too early
916 G4ExceptionDescription ed ; 991 G4ExceptionDescription ed ;
917 ed << "Corrupted data file for Z=" << Z 992 ed << "Corrupted data file for Z=" << Z << G4endl;
918 ed << "Found less than " << nPoints << " 993 ed << "Found less than " << nPoints << " entries" << G4endl;
919 G4Exception("G4PenelopeRayleighModelMI:: 994 G4Exception("G4PenelopeRayleighModelMI::ReadDataFile()",
920 "em0005",FatalException,ed); 995 "em0005",FatalException,ed);
921 } 996 }
922 } 997 }
>> 998
>> 999 if (!logAtomicCrossSection) {
>> 1000 G4Exception("G4PenelopeRayleighModelMI::ReadDataFile()",
>> 1001 "em2044",FatalException,"Unable to allocate the atomic cross section table");
>> 1002 delete theVec;
>> 1003 return;
>> 1004 }
>> 1005 logAtomicCrossSection->insert(std::make_pair(Z,theVec));
923 file.close(); 1006 file.close();
924 1007
925 //Then, read the extended momentum transfer 1008 //Then, read the extended momentum transfer file
926 std::ostringstream ost2; 1009 std::ostringstream ost2;
927 ost2 << path << "/penelope/rayleigh/MIFF/qex 1010 ost2 << path << "/penelope/rayleigh/MIFF/qext.dat";
928 file.open(ost2.str().c_str()); 1011 file.open(ost2.str().c_str());
929 1012
930 if (!file.is_open()) { 1013 if (!file.is_open()) {
931 G4String excep = "Data file " + G4String(o 1014 G4String excep = "Data file " + G4String(ost2.str()) + " not found!";
932 G4Exception("G4PenelopeRayleighModelMI::Re 1015 G4Exception("G4PenelopeRayleighModelMI::ReadDataFile()",
933 "em0003",FatalException,excep); 1016 "em0003",FatalException,excep);
934 } << 1017 }
>> 1018
935 G4bool fillQGrid = false; 1019 G4bool fillQGrid = false;
936 if (!fLogQSquareGrid.size()) { << 1020 if (!logQSquareGrid.size()) {
937 fillQGrid = true; 1021 fillQGrid = true;
938 nPoints = 1142; 1022 nPoints = 1142;
939 } 1023 }
>> 1024
940 G4double qext = 0; 1025 G4double qext = 0;
941 if (fillQGrid) { //fLogQSquareGrid filled o << 1026 if (fillQGrid) { //logQSquareGrid filled only one time
942 for (std::size_t i=0;i<nPoints;++i) { << 1027 for (size_t i=0;i<nPoints;i++) {
943 file >> qext; 1028 file >> qext;
944 fLogQSquareGrid.push_back(2.0*G4Log(qext << 1029 logQSquareGrid.push_back(2.0*G4Log(qext));
945 } 1030 }
946 } 1031 }
947 file.close(); 1032 file.close();
948 1033
949 //Finally, read the form factor file 1034 //Finally, read the form factor file
950 std::ostringstream ost3; 1035 std::ostringstream ost3;
951 if (Z>9) 1036 if (Z>9)
952 ost3 << path << "/penelope/rayleigh/pdaff" 1037 ost3 << path << "/penelope/rayleigh/pdaff" << Z << ".p08";
953 else << 1038 else
954 ost3 << path << "/penelope/rayleigh/pdaff0 << 1039 ost3 << path << "/penelope/rayleigh/pdaff0" << Z << ".p08";
955 file.open(ost3.str().c_str()); 1040 file.open(ost3.str().c_str());
956 1041
957 if (!file.is_open()) { 1042 if (!file.is_open()) {
958 G4String excep = "Data file " + G4String(o 1043 G4String excep = "Data file " + G4String(ost3.str()) + " not found!";
959 G4Exception("G4PenelopeRayleighModelMI::Re 1044 G4Exception("G4PenelopeRayleighModelMI::ReadDataFile()",
960 "em0003",FatalException,excep); 1045 "em0003",FatalException,excep);
961 } 1046 }
962 1047
963 file >> readZ >> nPoints; 1048 file >> readZ >> nPoints;
964 1049
965 if (readZ != Z || nPoints <= 0 || nPoints >= 1050 if (readZ != Z || nPoints <= 0 || nPoints >= 5000) {
966 G4ExceptionDescription ed; 1051 G4ExceptionDescription ed;
967 ed << "Corrupted data file for Z=" << Z << 1052 ed << "Corrupted data file for Z=" << Z << G4endl;
968 G4Exception("G4PenelopeRayleighModelMI::Re 1053 G4Exception("G4PenelopeRayleighModelMI::ReadDataFile()",
969 "em0005",FatalException,ed); 1054 "em0005",FatalException,ed);
970 return; 1055 return;
971 } 1056 }
972 1057
973 fAtomicFormFactor[Z] = new G4PhysicsFreeVect << 1058 G4PhysicsFreeVector* theFFVec = new G4PhysicsFreeVector((size_t)nPoints);
974 G4double q=0,ff=0,incoh=0; 1059 G4double q=0,ff=0,incoh=0;
975 for (std::size_t i=0;i<nPoints;++i) { << 1060 for (size_t i=0;i<nPoints;i++) {
976 file >> q >> ff >> incoh; //q and ff are d 1061 file >> q >> ff >> incoh; //q and ff are dimensionless (q is in units of (m_e*c))
977 fAtomicFormFactor[Z]->PutValue(i,q,ff); << 1062 theFFVec->PutValue(i,q,ff);
978 if (file.eof() && i != (nPoints-1)) { //fi 1063 if (file.eof() && i != (nPoints-1)) { //file ended too early
979 G4ExceptionDescription ed; 1064 G4ExceptionDescription ed;
980 ed << "Corrupted data file for Z=" << Z 1065 ed << "Corrupted data file for Z=" << Z << G4endl;
981 ed << "Found less than " << nPoints << " 1066 ed << "Found less than " << nPoints << " entries" << G4endl;
982 G4Exception("G4PenelopeRayleighModelMI:: 1067 G4Exception("G4PenelopeRayleighModelMI::ReadDataFile()",
983 "em0005",FatalException,ed); 1068 "em0005",FatalException,ed);
984 } 1069 }
985 } 1070 }
>> 1071
>> 1072 if (!atomicFormFactor) {
>> 1073 G4Exception("G4PenelopeRayleighModelMI::ReadDataFile()",
>> 1074 "em2045",FatalException,
>> 1075 "Unable to allocate the atomicFormFactor data table");
>> 1076 delete theFFVec;
>> 1077 return;
>> 1078 }
>> 1079
>> 1080 atomicFormFactor->insert(std::make_pair(Z,theFFVec));
986 file.close(); 1081 file.close();
987 return; 1082 return;
988 } 1083 }
989 1084
990 //....oooOO0OOooo........oooOO0OOooo........oo 1085 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
991 1086
992 void G4PenelopeRayleighModelMI::ReadMolInterfe 1087 void G4PenelopeRayleighModelMI::ReadMolInterferenceData(const G4String& matname,
993 const G4String& FFfilename) 1088 const G4String& FFfilename)
994 { 1089 {
995 if (fVerboseLevel > 2) { << 1090
>> 1091 if (verboseLevel > 2) {
996 G4cout << "G4PenelopeRayleighModelMI::Read 1092 G4cout << "G4PenelopeRayleighModelMI::ReadMolInterferenceData() for material " <<
997 matname << G4endl; 1093 matname << G4endl;
998 } 1094 }
999 G4bool isLocalFile = (FFfilename != "NULL"); 1095 G4bool isLocalFile = (FFfilename != "NULL");
1000 1096
1001 const char* path = G4FindDataDir("G4LEDATA" << 1097 char* path = std::getenv("G4LEDATA");
1002 if (!path) { 1098 if (!path) {
1003 G4String excep = "G4LEDATA environment va 1099 G4String excep = "G4LEDATA environment variable not set!";
1004 G4Exception("G4PenelopeRayleighModelMI::R 1100 G4Exception("G4PenelopeRayleighModelMI::ReadMolInterferenceData()",
1005 "em0006",FatalException,excep); 1101 "em0006",FatalException,excep);
1006 return; 1102 return;
1007 } 1103 }
1008 1104
1009 if (!(fKnownMaterials->count(matname)) && ! 1105 if (!(fKnownMaterials->count(matname)) && !isLocalFile) //material not found
1010 return; 1106 return;
1011 1107
1012 G4String aFileName = (isLocalFile) ? FFfile 1108 G4String aFileName = (isLocalFile) ? FFfilename : fKnownMaterials->find(matname)->second;
1013 1109
1014 //if the material has a MIFF, read it from 1110 //if the material has a MIFF, read it from the database
1015 if (aFileName != "") { 1111 if (aFileName != "") {
1016 if (fVerboseLevel > 2) << 1112 if (verboseLevel > 2)
1017 G4cout << "ReadMolInterferenceData(). R 1113 G4cout << "ReadMolInterferenceData(). Read material: " << matname << ", filename: " <<
1018 aFileName << " " << 1114 aFileName << " " <<
1019 (isLocalFile ? "(local)" : "(database)") << 1115 (isLocalFile ? "(local)" : "(database)") << G4endl;
1020 std::ifstream file; 1116 std::ifstream file;
1021 std::ostringstream ostIMFF; 1117 std::ostringstream ostIMFF;
1022 if (isLocalFile) 1118 if (isLocalFile)
1023 ostIMFF << aFileName; 1119 ostIMFF << aFileName;
1024 else 1120 else
1025 ostIMFF << path << "/penelope/rayleigh/ 1121 ostIMFF << path << "/penelope/rayleigh/MIFF/" << aFileName;
1026 file.open(ostIMFF.str().c_str()); 1122 file.open(ostIMFF.str().c_str());
1027 1123
1028 if (!file.is_open()) { 1124 if (!file.is_open()) {
1029 G4String excep = "Data file " + G4Strin 1125 G4String excep = "Data file " + G4String(ostIMFF.str()) + " not found!";
1030 G4Exception("G4PenelopeRayleighModelMI: 1126 G4Exception("G4PenelopeRayleighModelMI::ReadMolInterferenceData()",
1031 "em1031",FatalException,excep); 1127 "em1031",FatalException,excep);
1032 return; 1128 return;
1033 } 1129 }
1034 1130
1035 //check the number of points in the file 1131 //check the number of points in the file
1036 std::size_t nPoints = 0; << 1132 size_t nPoints = 0;
1037 G4double x=0,y=0; 1133 G4double x=0,y=0;
1038 while (file.good()) { 1134 while (file.good()) {
1039 file >> x >> y; 1135 file >> x >> y;
1040 nPoints++; 1136 nPoints++;
1041 } 1137 }
1042 file.close(); 1138 file.close();
1043 nPoints--; 1139 nPoints--;
1044 if (fVerboseLevel > 3) << 1140 if (verboseLevel > 3)
1045 G4cout << "Number of nPoints: " << nPoi 1141 G4cout << "Number of nPoints: " << nPoints << G4endl;
1046 1142
1047 //read the file 1143 //read the file
1048 file.open(ostIMFF.str().c_str()); 1144 file.open(ostIMFF.str().c_str());
1049 G4PhysicsFreeVector* theFFVec = new G4Phy << 1145 G4PhysicsFreeVector* theFFVec = new G4PhysicsFreeVector((size_t)nPoints);
1050 G4double q=0,ff=0; 1146 G4double q=0,ff=0;
1051 for (std::size_t i=0; i<nPoints; ++i) { << 1147 for (size_t i=0; i<nPoints; i++) {
1052 file >> q >> ff; 1148 file >> q >> ff;
1053 1149
1054 //q and ff are dimensionless (q is in u 1150 //q and ff are dimensionless (q is in units of (m_e*c))
1055 theFFVec->PutValue(i,q,ff); 1151 theFFVec->PutValue(i,q,ff);
1056 1152
1057 //file ended too early 1153 //file ended too early
1058 if (file.eof() && i != (nPoints-1)) { 1154 if (file.eof() && i != (nPoints-1)) {
1059 G4ExceptionDescription ed; 1155 G4ExceptionDescription ed;
1060 ed << "Corrupted data file" << G4endl; 1156 ed << "Corrupted data file" << G4endl;
1061 ed << "Found less than " << nPoints << " en 1157 ed << "Found less than " << nPoints << " entries" << G4endl;
1062 G4Exception("G4PenelopeRayleighModelMI::Rea 1158 G4Exception("G4PenelopeRayleighModelMI::ReadMolInterferenceData()",
1063 "em1005",FatalException,ed); 1159 "em1005",FatalException,ed);
1064 } 1160 }
1065 } 1161 }
1066 if (!fMolInterferenceData) { << 1162 if (!MolInterferenceData) {
1067 G4Exception("G4PenelopeRayleighModelMI: 1163 G4Exception("G4PenelopeRayleighModelMI::ReadMolInterferenceData()",
1068 "em2145",FatalException, 1164 "em2145",FatalException,
1069 "Unable to allocate the Molecular Inter 1165 "Unable to allocate the Molecular Interference data table");
1070 delete theFFVec; 1166 delete theFFVec;
1071 return; 1167 return;
1072 } 1168 }
1073 file.close(); 1169 file.close();
1074 fMolInterferenceData->insert(std::make_pa << 1170 MolInterferenceData->insert(std::make_pair(matname,theFFVec));
1075 } << 1171 }
>> 1172
1076 return; 1173 return;
1077 } 1174 }
1078 1175
1079 //....oooOO0OOooo........oooOO0OOooo........o 1176 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1080 1177
1081 G4double G4PenelopeRayleighModelMI::GetFSquar 1178 G4double G4PenelopeRayleighModelMI::GetFSquared(const G4Material* mat, const G4double QSquared)
1082 { 1179 {
1083 G4double f2 = 0; 1180 G4double f2 = 0;
1084 //Input value QSquared could be zero: prote 1181 //Input value QSquared could be zero: protect the log() below against
1085 //the FPE exception 1182 //the FPE exception
1086 1183
1087 //If Q<1e-10, set Q to 1e-10 1184 //If Q<1e-10, set Q to 1e-10
1088 G4double logQSquared = (QSquared>1e-10) ? G 1185 G4double logQSquared = (QSquared>1e-10) ? G4Log(QSquared) : -23.;
1089 //last value of the table 1186 //last value of the table
1090 G4double maxlogQ2 = fLogQSquareGrid[fLogQSq << 1187 G4double maxlogQ2 = logQSquareGrid[logQSquareGrid.size()-1];
1091 1188
1092 //now it should be all right 1189 //now it should be all right
1093 G4PhysicsFreeVector* theVec = fLogFormFacto << 1190 G4PhysicsFreeVector* theVec = logFormFactorTable->find(mat)->second;
1094 1191
1095 if (!theVec) { 1192 if (!theVec) {
1096 G4ExceptionDescription ed; 1193 G4ExceptionDescription ed;
1097 ed << "Unable to retrieve F squared table 1194 ed << "Unable to retrieve F squared table for " << mat->GetName() << G4endl;
1098 G4Exception("G4PenelopeRayleighModelMI::G 1195 G4Exception("G4PenelopeRayleighModelMI::GetFSquared()",
1099 "em2046",FatalException,ed); 1196 "em2046",FatalException,ed);
1100 return 0; 1197 return 0;
1101 } 1198 }
1102 1199
1103 if (logQSquared < -20) { //Q < 1e-9 1200 if (logQSquared < -20) { //Q < 1e-9
1104 G4double logf2 = (*theVec)[0]; //first va 1201 G4double logf2 = (*theVec)[0]; //first value of the table
1105 f2 = G4Exp(logf2); 1202 f2 = G4Exp(logf2);
1106 } 1203 }
1107 else if (logQSquared > maxlogQ2) 1204 else if (logQSquared > maxlogQ2)
1108 f2 =0; 1205 f2 =0;
1109 else { 1206 else {
1110 //log(Q^2) vs. log(F^2) 1207 //log(Q^2) vs. log(F^2)
1111 G4double logf2 = theVec->Value(logQSquare 1208 G4double logf2 = theVec->Value(logQSquared);
1112 f2 = G4Exp(logf2); 1209 f2 = G4Exp(logf2);
1113 } 1210 }
1114 1211
1115 if (fVerboseLevel > 3) { << 1212 if (verboseLevel > 3) {
1116 G4cout << "G4PenelopeRayleighModelMI::Get 1213 G4cout << "G4PenelopeRayleighModelMI::GetFSquared() in " << mat->GetName() << G4endl;
1117 G4cout << "Q^2 = " << QSquared << " (uni 1214 G4cout << "Q^2 = " << QSquared << " (units of 1/(m_e*c)); F^2 = " << f2 << G4endl;
1118 } 1215 }
1119 return f2; 1216 return f2;
1120 } 1217 }
1121 1218
1122 //....oooOO0OOooo........oooOO0OOooo........o 1219 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1123 1220
1124 void G4PenelopeRayleighModelMI::InitializeSam 1221 void G4PenelopeRayleighModelMI::InitializeSamplingAlgorithm(const G4Material* mat)
1125 { 1222 {
>> 1223
1126 G4double q2min = 0; 1224 G4double q2min = 0;
1127 G4double q2max = 0; 1225 G4double q2max = 0;
1128 const std::size_t np = 150; //hard-coded in << 1226 const size_t np = 150; //hard-coded in Penelope
1129 for (std::size_t i=1;i<fLogQSquareGrid.size << 1227 // G4cout << "Init N= " << logQSquareGrid.size() << G4endl;
>> 1228 for (size_t i=1;i<logQSquareGrid.size();i++)
1130 { 1229 {
1131 G4double Q2 = G4Exp(fLogQSquareGrid[i]) << 1230 G4double Q2 = G4Exp(logQSquareGrid[i]);
1132 if (GetFSquared(mat,Q2) > 1e-35) 1231 if (GetFSquared(mat,Q2) > 1e-35)
1133 { 1232 {
1134 q2max = G4Exp(fLogQSquareGrid[i-1]); << 1233 q2max = G4Exp(logQSquareGrid[i-1]);
1135 } 1234 }
>> 1235 //G4cout << "Q2= " << Q2 << " q2max= " << q2max << G4endl;
1136 } 1236 }
1137 std::size_t nReducedPoints = np/4; << 1237
>> 1238 size_t nReducedPoints = np/4;
1138 1239
1139 //check for errors 1240 //check for errors
1140 if (np < 16) 1241 if (np < 16)
1141 { 1242 {
1142 G4Exception("G4PenelopeRayleighModelMI: 1243 G4Exception("G4PenelopeRayleighModelMI::InitializeSamplingAlgorithm()",
1143 "em2047",FatalException, 1244 "em2047",FatalException,
1144 "Too few points to initialize the sampl 1245 "Too few points to initialize the sampling algorithm");
1145 } 1246 }
1146 if (q2min > (q2max-1e-10)) 1247 if (q2min > (q2max-1e-10))
1147 { 1248 {
1148 G4cout << "q2min= " << q2min << " q2max 1249 G4cout << "q2min= " << q2min << " q2max= " << q2max << G4endl;
1149 G4Exception("G4PenelopeRayleighModelMI: 1250 G4Exception("G4PenelopeRayleighModelMI::InitializeSamplingAlgorithm()",
1150 "em2048",FatalException, 1251 "em2048",FatalException,
1151 "Too narrow grid to initialize the samp 1252 "Too narrow grid to initialize the sampling algorithm");
1152 } 1253 }
1153 1254
1154 //This is subroutine RITAI0 of Penelope 1255 //This is subroutine RITAI0 of Penelope
1155 //Create an object of type G4PenelopeRaylei 1256 //Create an object of type G4PenelopeRayleighSamplingData --> store in a map::Material*
1156 1257
1157 //temporary vectors --> Then everything is 1258 //temporary vectors --> Then everything is stored in G4PenelopeSamplingData
1158 G4DataVector* x = new G4DataVector(); 1259 G4DataVector* x = new G4DataVector();
1159 1260
1160 /****************************************** 1261 /*******************************************************************************
1161 Start with a grid of NUNIF points uniform 1262 Start with a grid of NUNIF points uniformly spaced in the interval q2min,q2max
1162 ******************************************* 1263 ********************************************************************************/
1163 std::size_t NUNIF = std::min(std::max(((std << 1264 size_t NUNIF = std::min(std::max(((size_t)8),nReducedPoints),np/2);
1164 const G4int nip = 51; //hard-coded in Penel 1265 const G4int nip = 51; //hard-coded in Penelope
1165 1266
1166 G4double dx = (q2max-q2min)/((G4double) NUN 1267 G4double dx = (q2max-q2min)/((G4double) NUNIF-1);
1167 x->push_back(q2min); 1268 x->push_back(q2min);
1168 for (std::size_t i=1;i<NUNIF-1;++i) << 1269 for (size_t i=1;i<NUNIF-1;i++)
1169 { 1270 {
1170 G4double app = q2min + i*dx; 1271 G4double app = q2min + i*dx;
1171 x->push_back(app); //increase 1272 x->push_back(app); //increase
1172 } 1273 }
1173 x->push_back(q2max); 1274 x->push_back(q2max);
1174 1275
1175 if (fVerboseLevel> 3) << 1276 if (verboseLevel> 3)
1176 G4cout << "Vector x has " << x->size() << 1277 G4cout << "Vector x has " << x->size() << " points, while NUNIF = " << NUNIF << G4endl;
1177 1278
1178 //Allocate temporary storage vectors 1279 //Allocate temporary storage vectors
1179 G4DataVector* area = new G4DataVector(); 1280 G4DataVector* area = new G4DataVector();
1180 G4DataVector* a = new G4DataVector(); 1281 G4DataVector* a = new G4DataVector();
1181 G4DataVector* b = new G4DataVector(); 1282 G4DataVector* b = new G4DataVector();
1182 G4DataVector* c = new G4DataVector(); 1283 G4DataVector* c = new G4DataVector();
1183 G4DataVector* err = new G4DataVector(); 1284 G4DataVector* err = new G4DataVector();
1184 1285
1185 for (std::size_t i=0;i<NUNIF-1;++i) //build << 1286 for (size_t i=0;i<NUNIF-1;i++) //build all points but the last
1186 { 1287 {
1187 //Temporary vectors for this loop 1288 //Temporary vectors for this loop
1188 G4DataVector* pdfi = new G4DataVector() 1289 G4DataVector* pdfi = new G4DataVector();
1189 G4DataVector* pdfih = new G4DataVector( 1290 G4DataVector* pdfih = new G4DataVector();
1190 G4DataVector* sumi = new G4DataVector() 1291 G4DataVector* sumi = new G4DataVector();
1191 1292
1192 G4double dxi = ((*x)[i+1]-(*x)[i])/(G4d 1293 G4double dxi = ((*x)[i+1]-(*x)[i])/(G4double (nip-1));
1193 G4double pdfmax = 0; 1294 G4double pdfmax = 0;
1194 for (G4int k=0;k<nip;k++) 1295 for (G4int k=0;k<nip;k++)
1195 { 1296 {
1196 G4double xik = (*x)[i]+k*dxi; 1297 G4double xik = (*x)[i]+k*dxi;
1197 G4double pdfk = std::max(GetFSquared(mat, 1298 G4double pdfk = std::max(GetFSquared(mat,xik),0.);
1198 pdfi->push_back(pdfk); 1299 pdfi->push_back(pdfk);
1199 pdfmax = std::max(pdfmax,pdfk); 1300 pdfmax = std::max(pdfmax,pdfk);
1200 if (k < (nip-1)) 1301 if (k < (nip-1))
1201 { 1302 {
1202 G4double xih = xik + 0.5*dxi; 1303 G4double xih = xik + 0.5*dxi;
1203 G4double pdfIK = std::max(GetFSquared 1304 G4double pdfIK = std::max(GetFSquared(mat,xih),0.);
1204 pdfih->push_back(pdfIK); 1305 pdfih->push_back(pdfIK);
1205 pdfmax = std::max(pdfmax,pdfIK); 1306 pdfmax = std::max(pdfmax,pdfIK);
1206 } 1307 }
1207 } 1308 }
1208 1309
1209 //Simpson's integration 1310 //Simpson's integration
1210 G4double cons = dxi*0.5*(1./3.); 1311 G4double cons = dxi*0.5*(1./3.);
1211 sumi->push_back(0.); 1312 sumi->push_back(0.);
1212 for (G4int k=1;k<nip;k++) 1313 for (G4int k=1;k<nip;k++)
1213 { 1314 {
1214 G4double previous = (*sumi)[k-1]; 1315 G4double previous = (*sumi)[k-1];
1215 G4double next = previous + cons*((*pdfi)[ 1316 G4double next = previous + cons*((*pdfi)[k-1]+4.0*(*pdfih)[k-1]+(*pdfi)[k]);
1216 sumi->push_back(next); 1317 sumi->push_back(next);
1217 } 1318 }
1218 1319
1219 G4double lastIntegral = (*sumi)[sumi->s 1320 G4double lastIntegral = (*sumi)[sumi->size()-1];
1220 area->push_back(lastIntegral); 1321 area->push_back(lastIntegral);
1221 //Normalize cumulative function 1322 //Normalize cumulative function
1222 G4double factor = 1.0/lastIntegral; 1323 G4double factor = 1.0/lastIntegral;
1223 for (std::size_t k=0;k<sumi->size();++k << 1324 for (size_t k=0;k<sumi->size();k++)
1224 (*sumi)[k] *= factor; 1325 (*sumi)[k] *= factor;
1225 1326
1226 //When the PDF vanishes at one of the i 1327 //When the PDF vanishes at one of the interval end points, its value is modified
1227 if ((*pdfi)[0] < 1e-35) 1328 if ((*pdfi)[0] < 1e-35)
1228 (*pdfi)[0] = 1e-5*pdfmax; 1329 (*pdfi)[0] = 1e-5*pdfmax;
1229 if ((*pdfi)[pdfi->size()-1] < 1e-35) 1330 if ((*pdfi)[pdfi->size()-1] < 1e-35)
1230 (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax; 1331 (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax;
1231 1332
1232 G4double pli = (*pdfi)[0]*factor; 1333 G4double pli = (*pdfi)[0]*factor;
1233 G4double pui = (*pdfi)[pdfi->size()-1]* 1334 G4double pui = (*pdfi)[pdfi->size()-1]*factor;
1234 G4double B_temp = 1.0-1.0/(pli*pui*dx*d 1335 G4double B_temp = 1.0-1.0/(pli*pui*dx*dx);
1235 G4double A_temp = (1.0/(pli*dx))-1.0-B_ 1336 G4double A_temp = (1.0/(pli*dx))-1.0-B_temp;
1236 G4double C_temp = 1.0+A_temp+B_temp; 1337 G4double C_temp = 1.0+A_temp+B_temp;
1237 if (C_temp < 1e-35) 1338 if (C_temp < 1e-35)
1238 { 1339 {
1239 a->push_back(0.); 1340 a->push_back(0.);
1240 b->push_back(0.); 1341 b->push_back(0.);
1241 c->push_back(1.); 1342 c->push_back(1.);
1242 } 1343 }
1243 else 1344 else
1244 { 1345 {
1245 a->push_back(A_temp); 1346 a->push_back(A_temp);
1246 b->push_back(B_temp); 1347 b->push_back(B_temp);
1247 c->push_back(C_temp); 1348 c->push_back(C_temp);
1248 } 1349 }
1249 1350
1250 //OK, now get ERR(I), the integral of t 1351 //OK, now get ERR(I), the integral of the absolute difference between the rational interpolation
1251 //and the true pdf, extended over the i 1352 //and the true pdf, extended over the interval (X(I),X(I+1))
1252 G4int icase = 1; //loop code 1353 G4int icase = 1; //loop code
1253 G4bool reLoop = false; 1354 G4bool reLoop = false;
1254 err->push_back(0.); 1355 err->push_back(0.);
1255 do 1356 do
1256 { 1357 {
1257 reLoop = false; 1358 reLoop = false;
1258 (*err)[i] = 0.; //zero variable 1359 (*err)[i] = 0.; //zero variable
1259 for (G4int k=0;k<nip;k++) 1360 for (G4int k=0;k<nip;k++)
1260 { 1361 {
1261 G4double rr = (*sumi)[k]; 1362 G4double rr = (*sumi)[k];
1262 G4double pap = (*area)[i]*(1.0+((*a)[ 1363 G4double pap = (*area)[i]*(1.0+((*a)[i]+(*b)[i]*rr)*rr)*(1.0+((*a)[i]+(*b)[i]*rr)*rr)/
1263 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-( 1364 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(*x)[i]));
1264 if (k == 0 || k == nip-1) 1365 if (k == 0 || k == nip-1)
1265 (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k] 1366 (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k]);
1266 else 1367 else
1267 (*err)[i] += std::fabs(pap-(*pdfi)[k]); 1368 (*err)[i] += std::fabs(pap-(*pdfi)[k]);
1268 } 1369 }
1269 (*err)[i] *= dxi; 1370 (*err)[i] *= dxi;
1270 1371
1271 //If err(I) is too large, the pdf is appr 1372 //If err(I) is too large, the pdf is approximated by a uniform distribution
1272 if ((*err)[i] > 0.1*(*area)[i] && icase = 1373 if ((*err)[i] > 0.1*(*area)[i] && icase == 1)
1273 { 1374 {
1274 (*b)[i] = 0; 1375 (*b)[i] = 0;
1275 (*a)[i] = 0; 1376 (*a)[i] = 0;
1276 (*c)[i] = 1.; 1377 (*c)[i] = 1.;
1277 icase = 2; 1378 icase = 2;
1278 reLoop = true; 1379 reLoop = true;
1279 } 1380 }
1280 }while(reLoop); 1381 }while(reLoop);
1281 1382
1282 delete pdfi; 1383 delete pdfi;
1283 delete pdfih; 1384 delete pdfih;
1284 delete sumi; 1385 delete sumi;
1285 } //end of first loop over i 1386 } //end of first loop over i
1286 1387
1287 //Now assign last point 1388 //Now assign last point
1288 (*x)[x->size()-1] = q2max; 1389 (*x)[x->size()-1] = q2max;
1289 a->push_back(0.); 1390 a->push_back(0.);
1290 b->push_back(0.); 1391 b->push_back(0.);
1291 c->push_back(0.); 1392 c->push_back(0.);
1292 err->push_back(0.); 1393 err->push_back(0.);
1293 area->push_back(0.); 1394 area->push_back(0.);
1294 1395
1295 if (x->size() != NUNIF || a->size() != NUNI 1396 if (x->size() != NUNIF || a->size() != NUNIF ||
1296 err->size() != NUNIF || area->size() != 1397 err->size() != NUNIF || area->size() != NUNIF)
1297 { 1398 {
1298 G4ExceptionDescription ed; 1399 G4ExceptionDescription ed;
1299 ed << "Problem in building the Table fo 1400 ed << "Problem in building the Table for Sampling: array dimensions do not match" << G4endl;
1300 G4Exception("G4PenelopeRayleighModelMI: 1401 G4Exception("G4PenelopeRayleighModelMI::InitializeSamplingAlgorithm()",
1301 "em2049",FatalException,ed); 1402 "em2049",FatalException,ed);
1302 } 1403 }
1303 1404
1304 /****************************************** 1405 /*******************************************************************************
1305 New grid points are added by halving the s 1406 New grid points are added by halving the sub-intervals with the largest absolute error
1306 This is done up to np=150 points in the gri 1407 This is done up to np=150 points in the grid
1307 ******************************************* 1408 ********************************************************************************/
1308 do 1409 do
1309 { 1410 {
1310 G4double maxError = 0.0; 1411 G4double maxError = 0.0;
1311 std::size_t iErrMax = 0; << 1412 size_t iErrMax = 0;
1312 for (std::size_t i=0;i<err->size()-2;++ << 1413 for (size_t i=0;i<err->size()-2;i++)
1313 { 1414 {
1314 //maxError is the lagest of the interval 1415 //maxError is the lagest of the interval errors err[i]
1315 if ((*err)[i] > maxError) 1416 if ((*err)[i] > maxError)
1316 { 1417 {
1317 maxError = (*err)[i]; 1418 maxError = (*err)[i];
1318 iErrMax = i; 1419 iErrMax = i;
1319 } 1420 }
1320 } 1421 }
1321 1422
1322 //OK, now I have to insert one new poin 1423 //OK, now I have to insert one new point in the position iErrMax
1323 G4double newx = 0.5*((*x)[iErrMax]+(*x) 1424 G4double newx = 0.5*((*x)[iErrMax]+(*x)[iErrMax+1]);
1324 1425
1325 x->insert(x->begin()+iErrMax+1,newx); 1426 x->insert(x->begin()+iErrMax+1,newx);
1326 //Add place-holders in the other vector 1427 //Add place-holders in the other vectors
1327 area->insert(area->begin()+iErrMax+1,0. 1428 area->insert(area->begin()+iErrMax+1,0.);
1328 a->insert(a->begin()+iErrMax+1,0.); 1429 a->insert(a->begin()+iErrMax+1,0.);
1329 b->insert(b->begin()+iErrMax+1,0.); 1430 b->insert(b->begin()+iErrMax+1,0.);
1330 c->insert(c->begin()+iErrMax+1,0.); 1431 c->insert(c->begin()+iErrMax+1,0.);
1331 err->insert(err->begin()+iErrMax+1,0.); 1432 err->insert(err->begin()+iErrMax+1,0.);
1332 1433
1333 //Now calculate the other parameters 1434 //Now calculate the other parameters
1334 for (std::size_t i=iErrMax;i<=iErrMax+1 << 1435 for (size_t i=iErrMax;i<=iErrMax+1;i++)
1335 { 1436 {
1336 //Temporary vectors for this loop 1437 //Temporary vectors for this loop
1337 G4DataVector* pdfi = new G4DataVector(); 1438 G4DataVector* pdfi = new G4DataVector();
1338 G4DataVector* pdfih = new G4DataVector(); 1439 G4DataVector* pdfih = new G4DataVector();
1339 G4DataVector* sumi = new G4DataVector(); 1440 G4DataVector* sumi = new G4DataVector();
1340 1441
1341 G4double dxLocal = (*x)[i+1]-(*x)[i]; 1442 G4double dxLocal = (*x)[i+1]-(*x)[i];
1342 G4double dxi = ((*x)[i+1]-(*x)[i])/(G4dou 1443 G4double dxi = ((*x)[i+1]-(*x)[i])/(G4double (nip-1));
1343 G4double pdfmax = 0; 1444 G4double pdfmax = 0;
1344 for (G4int k=0;k<nip;k++) 1445 for (G4int k=0;k<nip;k++)
1345 { 1446 {
1346 G4double xik = (*x)[i]+k*dxi; 1447 G4double xik = (*x)[i]+k*dxi;
1347 G4double pdfk = std::max(GetFSquared( 1448 G4double pdfk = std::max(GetFSquared(mat,xik),0.);
1348 pdfi->push_back(pdfk); 1449 pdfi->push_back(pdfk);
1349 pdfmax = std::max(pdfmax,pdfk); 1450 pdfmax = std::max(pdfmax,pdfk);
1350 if (k < (nip-1)) 1451 if (k < (nip-1))
1351 { 1452 {
1352 G4double xih = xik + 0.5*dxi; 1453 G4double xih = xik + 0.5*dxi;
1353 G4double pdfIK = std::max(GetFSquared(m 1454 G4double pdfIK = std::max(GetFSquared(mat,xih),0.);
1354 pdfih->push_back(pdfIK); 1455 pdfih->push_back(pdfIK);
1355 pdfmax = std::max(pdfmax,pdfIK); 1456 pdfmax = std::max(pdfmax,pdfIK);
1356 } 1457 }
1357 } 1458 }
1358 1459
1359 //Simpson's integration 1460 //Simpson's integration
1360 G4double cons = dxi*0.5*(1./3.); 1461 G4double cons = dxi*0.5*(1./3.);
1361 sumi->push_back(0.); 1462 sumi->push_back(0.);
1362 for (G4int k=1;k<nip;k++) 1463 for (G4int k=1;k<nip;k++)
1363 { 1464 {
1364 G4double previous = (*sumi)[k-1]; 1465 G4double previous = (*sumi)[k-1];
1365 G4double next = previous + cons*((*pd 1466 G4double next = previous + cons*((*pdfi)[k-1]+4.0*(*pdfih)[k-1]+(*pdfi)[k]);
1366 sumi->push_back(next); 1467 sumi->push_back(next);
1367 } 1468 }
1368 G4double lastIntegral = (*sumi)[sumi->siz 1469 G4double lastIntegral = (*sumi)[sumi->size()-1];
1369 (*area)[i] = lastIntegral; 1470 (*area)[i] = lastIntegral;
1370 1471
1371 //Normalize cumulative function 1472 //Normalize cumulative function
1372 G4double factor = 1.0/lastIntegral; 1473 G4double factor = 1.0/lastIntegral;
1373 for (std::size_t k=0;k<sumi->size();++k) << 1474 for (size_t k=0;k<sumi->size();k++)
1374 (*sumi)[k] *= factor; 1475 (*sumi)[k] *= factor;
1375 1476
1376 //When the PDF vanishes at one of the int 1477 //When the PDF vanishes at one of the interval end points, its value is modified
1377 if ((*pdfi)[0] < 1e-35) 1478 if ((*pdfi)[0] < 1e-35)
1378 (*pdfi)[0] = 1e-5*pdfmax; 1479 (*pdfi)[0] = 1e-5*pdfmax;
1379 if ((*pdfi)[pdfi->size()-1] < 1e-35) 1480 if ((*pdfi)[pdfi->size()-1] < 1e-35)
1380 (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax; 1481 (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax;
1381 1482
1382 G4double pli = (*pdfi)[0]*factor; 1483 G4double pli = (*pdfi)[0]*factor;
1383 G4double pui = (*pdfi)[pdfi->size()-1]*fa 1484 G4double pui = (*pdfi)[pdfi->size()-1]*factor;
1384 G4double B_temp = 1.0-1.0/(pli*pui*dxLoca 1485 G4double B_temp = 1.0-1.0/(pli*pui*dxLocal*dxLocal);
1385 G4double A_temp = (1.0/(pli*dxLocal))-1.0 1486 G4double A_temp = (1.0/(pli*dxLocal))-1.0-B_temp;
1386 G4double C_temp = 1.0+A_temp+B_temp; 1487 G4double C_temp = 1.0+A_temp+B_temp;
1387 if (C_temp < 1e-35) 1488 if (C_temp < 1e-35)
1388 { 1489 {
1389 (*a)[i]= 0.; 1490 (*a)[i]= 0.;
1390 (*b)[i] = 0.; 1491 (*b)[i] = 0.;
1391 (*c)[i] = 1; 1492 (*c)[i] = 1;
1392 } 1493 }
1393 else 1494 else
1394 { 1495 {
1395 (*a)[i]= A_temp; 1496 (*a)[i]= A_temp;
1396 (*b)[i] = B_temp; 1497 (*b)[i] = B_temp;
1397 (*c)[i] = C_temp; 1498 (*c)[i] = C_temp;
1398 } 1499 }
1399 //OK, now get ERR(I), the integral of the 1500 //OK, now get ERR(I), the integral of the absolute difference between the rational interpolation
1400 //and the true pdf, extended over the int 1501 //and the true pdf, extended over the interval (X(I),X(I+1))
1401 G4int icase = 1; //loop code 1502 G4int icase = 1; //loop code
1402 G4bool reLoop = false; 1503 G4bool reLoop = false;
1403 do 1504 do
1404 { 1505 {
1405 reLoop = false; 1506 reLoop = false;
1406 (*err)[i] = 0.; //zero variable 1507 (*err)[i] = 0.; //zero variable
1407 for (G4int k=0;k<nip;k++) 1508 for (G4int k=0;k<nip;k++)
1408 { 1509 {
1409 G4double rr = (*sumi)[k]; 1510 G4double rr = (*sumi)[k];
1410 G4double pap = (*area)[i]*(1.0+((*a)[i] 1511 G4double pap = (*area)[i]*(1.0+((*a)[i]+(*b)[i]*rr)*rr)*(1.0+((*a)[i]+(*b)[i]*rr)*rr)/
1411 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+ 1512 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(*x)[i]));
1412 if (k == 0 || k == nip-1) 1513 if (k == 0 || k == nip-1)
1413 (*err)[i] += 0.5*std::fabs(pap-(*pdfi 1514 (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k]);
1414 else 1515 else
1415 (*err)[i] += std::fabs(pap-(*pdfi)[k] 1516 (*err)[i] += std::fabs(pap-(*pdfi)[k]);
1416 } 1517 }
1417 (*err)[i] *= dxi; 1518 (*err)[i] *= dxi;
1418 1519
1419 //If err(I) is too large, the pdf is 1520 //If err(I) is too large, the pdf is approximated by a uniform distribution
1420 if ((*err)[i] > 0.1*(*area)[i] && ica 1521 if ((*err)[i] > 0.1*(*area)[i] && icase == 1)
1421 { 1522 {
1422 (*b)[i] = 0; 1523 (*b)[i] = 0;
1423 (*a)[i] = 0; 1524 (*a)[i] = 0;
1424 (*c)[i] = 1.; 1525 (*c)[i] = 1.;
1425 icase = 2; 1526 icase = 2;
1426 reLoop = true; 1527 reLoop = true;
1427 } 1528 }
1428 }while(reLoop); 1529 }while(reLoop);
1429 delete pdfi; 1530 delete pdfi;
1430 delete pdfih; 1531 delete pdfih;
1431 delete sumi; 1532 delete sumi;
1432 } 1533 }
1433 }while(x->size() < np); 1534 }while(x->size() < np);
1434 1535
1435 if (x->size() != np || a->size() != np || 1536 if (x->size() != np || a->size() != np ||
1436 err->size() != np || area->size() != np 1537 err->size() != np || area->size() != np)
1437 { 1538 {
1438 G4Exception("G4PenelopeRayleighModelMI: 1539 G4Exception("G4PenelopeRayleighModelMI::InitializeSamplingAlgorithm()",
1439 "em2050",FatalException, 1540 "em2050",FatalException,
1440 "Problem in building the extended Table 1541 "Problem in building the extended Table for Sampling: array dimensions do not match ");
1441 } 1542 }
1442 1543
1443 /****************************************** 1544 /*******************************************************************************
1444 Renormalization 1545 Renormalization
1445 ******************************************* 1546 ********************************************************************************/
1446 G4double ws = 0; 1547 G4double ws = 0;
1447 for (std::size_t i=0;i<np-1;++i) << 1548 for (size_t i=0;i<np-1;i++)
1448 ws += (*area)[i]; 1549 ws += (*area)[i];
1449 ws = 1.0/ws; 1550 ws = 1.0/ws;
1450 G4double errMax = 0; 1551 G4double errMax = 0;
1451 for (std::size_t i=0;i<np-1;++i) << 1552 for (size_t i=0;i<np-1;i++)
1452 { 1553 {
1453 (*area)[i] *= ws; 1554 (*area)[i] *= ws;
1454 (*err)[i] *= ws; 1555 (*err)[i] *= ws;
1455 errMax = std::max(errMax,(*err)[i]); 1556 errMax = std::max(errMax,(*err)[i]);
1456 } 1557 }
1457 1558
1458 //Vector with the normalized cumulative dis 1559 //Vector with the normalized cumulative distribution
1459 G4DataVector* PAC = new G4DataVector(); 1560 G4DataVector* PAC = new G4DataVector();
1460 PAC->push_back(0.); 1561 PAC->push_back(0.);
1461 for (std::size_t i=0;i<np-1;++i) << 1562 for (size_t i=0;i<np-1;i++)
1462 { 1563 {
1463 G4double previous = (*PAC)[i]; 1564 G4double previous = (*PAC)[i];
1464 PAC->push_back(previous+(*area)[i]); 1565 PAC->push_back(previous+(*area)[i]);
1465 } 1566 }
1466 (*PAC)[PAC->size()-1] = 1.; 1567 (*PAC)[PAC->size()-1] = 1.;
1467 1568
1468 /****************************************** 1569 /*******************************************************************************
1469 Pre-calculated limits for the initial binar 1570 Pre-calculated limits for the initial binary search for subsequent sampling
1470 ******************************************* 1571 ********************************************************************************/
1471 std::vector<std::size_t> *ITTL = new std::v << 1572
1472 std::vector<std::size_t> *ITTU = new std::v << 1573 //G4DataVector* ITTL = new G4DataVector();
>> 1574 std::vector<size_t> *ITTL = new std::vector<size_t>;
>> 1575 std::vector<size_t> *ITTU = new std::vector<size_t>;
1473 1576
1474 //Just create place-holders 1577 //Just create place-holders
1475 for (std::size_t i=0;i<np;++i) << 1578 for (size_t i=0;i<np;i++)
1476 { 1579 {
1477 ITTL->push_back(0); 1580 ITTL->push_back(0);
1478 ITTU->push_back(0); 1581 ITTU->push_back(0);
1479 } 1582 }
1480 1583
1481 G4double bin = 1.0/(np-1); 1584 G4double bin = 1.0/(np-1);
1482 (*ITTL)[0]=0; 1585 (*ITTL)[0]=0;
1483 for (std::size_t i=1;i<(np-1);++i) << 1586 for (size_t i=1;i<(np-1);i++)
1484 { 1587 {
1485 G4double ptst = i*bin; 1588 G4double ptst = i*bin;
1486 G4bool found = false; 1589 G4bool found = false;
1487 for (std::size_t j=(*ITTL)[i-1];j<np && << 1590 for (size_t j=(*ITTL)[i-1];j<np && !found;j++)
1488 { 1591 {
1489 if ((*PAC)[j] > ptst) 1592 if ((*PAC)[j] > ptst)
1490 { 1593 {
1491 (*ITTL)[i] = j-1; 1594 (*ITTL)[i] = j-1;
1492 (*ITTU)[i-1] = j; 1595 (*ITTU)[i-1] = j;
1493 found = true; 1596 found = true;
1494 } 1597 }
1495 } 1598 }
1496 } 1599 }
1497 (*ITTU)[ITTU->size()-2] = ITTU->size()-1; 1600 (*ITTU)[ITTU->size()-2] = ITTU->size()-1;
1498 (*ITTU)[ITTU->size()-1] = ITTU->size()-1; 1601 (*ITTU)[ITTU->size()-1] = ITTU->size()-1;
1499 (*ITTL)[ITTL->size()-1] = ITTU->size()-2; 1602 (*ITTL)[ITTL->size()-1] = ITTU->size()-2;
1500 1603
1501 if (ITTU->size() != np || ITTU->size() != n 1604 if (ITTU->size() != np || ITTU->size() != np)
1502 { 1605 {
1503 G4Exception("G4PenelopeRayleighModelMI: 1606 G4Exception("G4PenelopeRayleighModelMI::InitializeSamplingAlgorithm()",
1504 "em2051",FatalException, 1607 "em2051",FatalException,
1505 "Problem in building the Limit Tables f 1608 "Problem in building the Limit Tables for Sampling: array dimensions do not match");
1506 } 1609 }
1507 1610
>> 1611
1508 /****************************************** 1612 /********************************************************************************
1509 Copy tables 1613 Copy tables
1510 ******************************************* 1614 ********************************************************************************/
1511 G4PenelopeSamplingData* theTable = new G4Pe 1615 G4PenelopeSamplingData* theTable = new G4PenelopeSamplingData(np);
1512 for (std::size_t i=0;i<np;++i) << 1616 for (size_t i=0;i<np;i++)
1513 { 1617 {
1514 theTable->AddPoint((*x)[i],(*PAC)[i],(* 1618 theTable->AddPoint((*x)[i],(*PAC)[i],(*a)[i],(*b)[i],(*ITTL)[i],(*ITTU)[i]);
1515 } 1619 }
1516 if (fVerboseLevel > 2) << 1620
>> 1621 if (verboseLevel > 2)
1517 { 1622 {
1518 G4cout << "**************************** 1623 G4cout << "*************************************************************************" <<
1519 G4endl; 1624 G4endl;
1520 G4cout << "Sampling table for Penelope 1625 G4cout << "Sampling table for Penelope Rayleigh scattering in " << mat->GetName() << G4endl;
1521 theTable->DumpTable(); 1626 theTable->DumpTable();
1522 } 1627 }
1523 fSamplingTable->insert(std::make_pair(mat,t << 1628 samplingTable->insert(std::make_pair(mat,theTable));
>> 1629
1524 1630
1525 //Clean up temporary vectors 1631 //Clean up temporary vectors
1526 delete x; 1632 delete x;
1527 delete a; 1633 delete a;
1528 delete b; 1634 delete b;
1529 delete c; 1635 delete c;
1530 delete err; 1636 delete err;
1531 delete area; 1637 delete area;
1532 delete PAC; 1638 delete PAC;
1533 delete ITTL; 1639 delete ITTL;
1534 delete ITTU; 1640 delete ITTU;
1535 1641
1536 //DONE! 1642 //DONE!
1537 return; 1643 return;
>> 1644
1538 } 1645 }
1539 1646
1540 //....oooOO0OOooo........oooOO0OOooo........o 1647 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1541 1648
1542 void G4PenelopeRayleighModelMI::GetPMaxTable( 1649 void G4PenelopeRayleighModelMI::GetPMaxTable(const G4Material* mat)
1543 { 1650 {
1544 if (!fPMaxTable) << 1651
>> 1652 if (!pMaxTable)
1545 { 1653 {
1546 G4cout << "G4PenelopeRayleighModelMI::B 1654 G4cout << "G4PenelopeRayleighModelMI::BuildPMaxTable" << G4endl;
1547 G4cout << "Going to instanziate the fPM << 1655 G4cout << "Going to instanziate the pMaxTable !" << G4endl;
1548 G4cout << "That should _not_ be here! " 1656 G4cout << "That should _not_ be here! " << G4endl;
1549 fPMaxTable = new std::map<const G4Mater << 1657 pMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>;
1550 } 1658 }
1551 //check if the table is already there 1659 //check if the table is already there
1552 if (fPMaxTable->count(mat)) << 1660 if (pMaxTable->count(mat))
1553 return; 1661 return;
1554 1662
1555 //otherwise build it 1663 //otherwise build it
1556 if (!fSamplingTable) << 1664 if (!samplingTable)
1557 { 1665 {
1558 G4Exception("G4PenelopeRayleighModelMI: 1666 G4Exception("G4PenelopeRayleighModelMI::GetPMaxTable()",
1559 "em2052",FatalException, 1667 "em2052",FatalException,
1560 "SamplingTable is not properly instanti 1668 "SamplingTable is not properly instantiated");
1561 return; 1669 return;
1562 } 1670 }
1563 1671
1564 //This should not be: the sampling table is 1672 //This should not be: the sampling table is built before the p-table
1565 if (!fSamplingTable->count(mat)) << 1673 if (!samplingTable->count(mat))
1566 { 1674 {
1567 G4ExceptionDescription ed; << 1675 G4ExceptionDescription ed;
1568 ed << "Sampling table for material " << << 1676 ed << "Sampling table for material " << mat->GetName() << " not found";
1569 G4Exception("G4PenelopeRayleighModelMI: << 1677 G4Exception("G4PenelopeRayleighModelMI::GetPMaxTable()",
1570 "em2052",FatalException, 1678 "em2052",FatalException,
1571 ed); 1679 ed);
1572 return; << 1680 return;
1573 } 1681 }
1574 1682
1575 G4PenelopeSamplingData *theTable = fSamplin << 1683 G4PenelopeSamplingData *theTable = samplingTable->find(mat)->second;
1576 std::size_t tablePoints = theTable->GetNumb << 1684 size_t tablePoints = theTable->GetNumberOfStoredPoints();
1577 std::size_t nOfEnergyPoints = fLogEnergyGri << 1685
>> 1686 size_t nOfEnergyPoints = logEnergyGridPMax.size();
1578 G4PhysicsFreeVector* theVec = new G4Physics 1687 G4PhysicsFreeVector* theVec = new G4PhysicsFreeVector(nOfEnergyPoints);
1579 1688
1580 const std::size_t nip = 51; //hard-coded in << 1689 const size_t nip = 51; //hard-coded in Penelope
1581 1690
1582 for (std::size_t ie=0;ie<fLogEnergyGridPMax << 1691 for (size_t ie=0;ie<logEnergyGridPMax.size();ie++)
1583 { 1692 {
1584 G4double energy = G4Exp(fLogEnergyGridP << 1693 G4double energy = G4Exp(logEnergyGridPMax[ie]);
1585 G4double Qm = 2.0*energy/electron_mass_ 1694 G4double Qm = 2.0*energy/electron_mass_c2; //this is non-dimensional now
1586 G4double Qm2 = Qm*Qm; 1695 G4double Qm2 = Qm*Qm;
1587 G4double firstQ2 = theTable->GetX(0); 1696 G4double firstQ2 = theTable->GetX(0);
1588 G4double lastQ2 = theTable->GetX(tableP 1697 G4double lastQ2 = theTable->GetX(tablePoints-1);
1589 G4double thePMax = 0; 1698 G4double thePMax = 0;
1590 1699
1591 if (Qm2 > firstQ2) 1700 if (Qm2 > firstQ2)
1592 { 1701 {
1593 if (Qm2 < lastQ2) 1702 if (Qm2 < lastQ2)
1594 { 1703 {
1595 //bisection to look for the index of 1704 //bisection to look for the index of Qm
1596 std::size_t lowerBound = 0; << 1705 size_t lowerBound = 0;
1597 std::size_t upperBound = tablePoints- << 1706 size_t upperBound = tablePoints-1;
1598 while (lowerBound <= upperBound) 1707 while (lowerBound <= upperBound)
1599 { 1708 {
1600 std::size_t midBin = (lowerBound + uppe << 1709 size_t midBin = (lowerBound + upperBound)/2;
1601 if( Qm2 < theTable->GetX(midBin)) 1710 if( Qm2 < theTable->GetX(midBin))
1602 { upperBound = midBin-1; } 1711 { upperBound = midBin-1; }
1603 else 1712 else
1604 { lowerBound = midBin+1; } 1713 { lowerBound = midBin+1; }
1605 } 1714 }
1606 //upperBound is the output (but also 1715 //upperBound is the output (but also lowerBounf --> should be the same!)
1607 G4double Q1 = theTable->GetX(upperBou 1716 G4double Q1 = theTable->GetX(upperBound);
1608 G4double Q2 = Qm2; 1717 G4double Q2 = Qm2;
1609 G4double DQ = (Q2-Q1)/((G4double)(nip 1718 G4double DQ = (Q2-Q1)/((G4double)(nip-1));
1610 G4double theA = theTable->GetA(upperB 1719 G4double theA = theTable->GetA(upperBound);
1611 G4double theB = theTable->GetB(upperB 1720 G4double theB = theTable->GetB(upperBound);
1612 G4double thePAC = theTable->GetPAC(up 1721 G4double thePAC = theTable->GetPAC(upperBound);
1613 G4DataVector* fun = new G4DataVector( 1722 G4DataVector* fun = new G4DataVector();
1614 for (std::size_t k=0;k<nip;++k) << 1723 for (size_t k=0;k<nip;k++)
1615 { 1724 {
1616 G4double qi = Q1 + k*DQ; 1725 G4double qi = Q1 + k*DQ;
1617 G4double tau = (qi-Q1)/ 1726 G4double tau = (qi-Q1)/
1618 (theTable->GetX(upperBound+1)-Q1); 1727 (theTable->GetX(upperBound+1)-Q1);
1619 G4double con1 = 2.0*theB*tau; 1728 G4double con1 = 2.0*theB*tau;
1620 G4double ci = 1.0+theA+theB; 1729 G4double ci = 1.0+theA+theB;
1621 G4double con2 = ci-theA*tau; 1730 G4double con2 = ci-theA*tau;
1622 G4double etap = 0; 1731 G4double etap = 0;
1623 if (std::fabs(con1) > 1.0e-16*std::fabs 1732 if (std::fabs(con1) > 1.0e-16*std::fabs(con2))
1624 etap = con2*(1.0-std::sqrt(1.0-2.0*ta 1733 etap = con2*(1.0-std::sqrt(1.0-2.0*tau*con1/(con2*con2)))/con1;
1625 else 1734 else
1626 etap = tau/con2; 1735 etap = tau/con2;
1627 G4double theFun = (theTable->GetPAC(upp 1736 G4double theFun = (theTable->GetPAC(upperBound+1)-thePAC)*
1628 (1.0+(theA+theB*etap)*etap)*(1.0+(the 1737 (1.0+(theA+theB*etap)*etap)*(1.0+(theA+theB*etap)*etap)/
1629 ((1.0-theB*etap*etap)*ci*(theTable->G 1738 ((1.0-theB*etap*etap)*ci*(theTable->GetX(upperBound+1)-Q1));
1630 fun->push_back(theFun); 1739 fun->push_back(theFun);
1631 } 1740 }
1632 //Now intergrate numerically the fun 1741 //Now intergrate numerically the fun Cavalieri-Simpson's method
1633 G4DataVector* sum = new G4DataVector; 1742 G4DataVector* sum = new G4DataVector;
1634 G4double CONS = DQ*(1./12.); 1743 G4double CONS = DQ*(1./12.);
1635 G4double HCONS = 0.5*CONS; 1744 G4double HCONS = 0.5*CONS;
1636 sum->push_back(0.); 1745 sum->push_back(0.);
1637 G4double secondPoint = (*sum)[0] + 1746 G4double secondPoint = (*sum)[0] +
1638 (5.0*(*fun)[0]+8.0*(*fun)[1]-(*fun)[2])*C 1747 (5.0*(*fun)[0]+8.0*(*fun)[1]-(*fun)[2])*CONS;
1639 sum->push_back(secondPoint); 1748 sum->push_back(secondPoint);
1640 for (std::size_t hh=2;hh<nip-1;++hh) << 1749 for (size_t hh=2;hh<nip-1;hh++)
1641 { 1750 {
1642 G4double previous = (*sum)[hh-1]; 1751 G4double previous = (*sum)[hh-1];
1643 G4double next = previous+(13.0*((*fun)[ 1752 G4double next = previous+(13.0*((*fun)[hh-1]+(*fun)[hh])-
1644 (*fun)[hh+1]-(*fun)[hh-2])*HCON 1753 (*fun)[hh+1]-(*fun)[hh-2])*HCONS;
1645 sum->push_back(next); 1754 sum->push_back(next);
1646 } 1755 }
1647 G4double last = (*sum)[nip-2]+(5.0*(* 1756 G4double last = (*sum)[nip-2]+(5.0*(*fun)[nip-1]+8.0*(*fun)[nip-2]-
1648 (*fun)[nip-3])*CONS; 1757 (*fun)[nip-3])*CONS;
1649 sum->push_back(last); 1758 sum->push_back(last);
1650 thePMax = thePAC + (*sum)[sum->size() 1759 thePMax = thePAC + (*sum)[sum->size()-1]; //last point
1651 delete fun; 1760 delete fun;
1652 delete sum; 1761 delete sum;
1653 } 1762 }
1654 else 1763 else
1655 { 1764 {
1656 thePMax = 1.0; 1765 thePMax = 1.0;
1657 } 1766 }
1658 } 1767 }
1659 else 1768 else
1660 { 1769 {
1661 thePMax = theTable->GetPAC(0); 1770 thePMax = theTable->GetPAC(0);
1662 } 1771 }
1663 1772
1664 //Write number in the table 1773 //Write number in the table
1665 theVec->PutValue(ie,energy,thePMax); 1774 theVec->PutValue(ie,energy,thePMax);
1666 } << 1775 }
1667 1776
1668 fPMaxTable->insert(std::make_pair(mat,theVe << 1777 pMaxTable->insert(std::make_pair(mat,theVec));
1669 return; 1778 return;
>> 1779
1670 } 1780 }
1671 1781
1672 //....oooOO0OOooo........oooOO0OOooo........o 1782 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1673 1783
1674 void G4PenelopeRayleighModelMI::DumpFormFacto 1784 void G4PenelopeRayleighModelMI::DumpFormFactorTable(const G4Material* mat)
1675 { 1785 {
1676 G4cout << "******************************** 1786 G4cout << "*****************************************************************" << G4endl;
1677 G4cout << "G4PenelopeRayleighModelMI: Form 1787 G4cout << "G4PenelopeRayleighModelMI: Form Factor Table for " << mat->GetName() << G4endl;
1678 //try to use the same format as Penelope-Fo 1788 //try to use the same format as Penelope-Fortran, namely Q (/m_e*c) and F
1679 G4cout << "Q/(m_e*c) F(Q) 1789 G4cout << "Q/(m_e*c) F(Q) " << G4endl;
1680 G4cout << "******************************** 1790 G4cout << "*****************************************************************" << G4endl;
1681 if (!fLogFormFactorTable->count(mat)) << 1791 if (!logFormFactorTable->count(mat))
1682 BuildFormFactorTable(mat); 1792 BuildFormFactorTable(mat);
1683 1793
1684 G4PhysicsFreeVector* theVec = fLogFormFacto << 1794 G4PhysicsFreeVector* theVec = logFormFactorTable->find(mat)->second;
1685 for (std::size_t i=0;i<theVec->GetVectorLen << 1795 for (size_t i=0;i<theVec->GetVectorLength();i++)
1686 { 1796 {
1687 G4double logQ2 = theVec->GetLowEdgeEner 1797 G4double logQ2 = theVec->GetLowEdgeEnergy(i);
1688 G4double Q = G4Exp(0.5*logQ2); 1798 G4double Q = G4Exp(0.5*logQ2);
1689 G4double logF2 = (*theVec)[i]; 1799 G4double logF2 = (*theVec)[i];
1690 G4double F = G4Exp(0.5*logF2); 1800 G4double F = G4Exp(0.5*logF2);
1691 G4cout << Q << " " << F << 1801 G4cout << Q << " " << F << G4endl;
1692 } 1802 }
1693 //DONE 1803 //DONE
1694 return; 1804 return;
1695 } 1805 }
1696 1806
1697 //....oooOO0OOooo........oooOO0OOooo........o 1807 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
1698 1808
1699 void G4PenelopeRayleighModelMI::SetParticle(c 1809 void G4PenelopeRayleighModelMI::SetParticle(const G4ParticleDefinition* p)
1700 { 1810 {
1701 if(!fParticle) { 1811 if(!fParticle) {
1702 fParticle = p; 1812 fParticle = p;
1703 } 1813 }
1704 } 1814 }
1705 1815
1706 //....oooOO0OOooo........oooOO0OOooo........o 1816 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
1707 void G4PenelopeRayleighModelMI::LoadKnownMIFF 1817 void G4PenelopeRayleighModelMI::LoadKnownMIFFMaterials()
1708 { 1818 {
1709 fKnownMaterials->insert(std::pair<G4String, 1819 fKnownMaterials->insert(std::pair<G4String,G4String>("Fat_MI","FF_fat_Tartari2002.dat"));
1710 fKnownMaterials->insert(std::pair<G4String, 1820 fKnownMaterials->insert(std::pair<G4String,G4String>("Water_MI","FF_water_Tartari2002.dat"));
1711 fKnownMaterials->insert(std::pair<G4String, 1821 fKnownMaterials->insert(std::pair<G4String,G4String>("BoneMatrix_MI","FF_bonematrix_Tartari2002.dat"));
1712 fKnownMaterials->insert(std::pair<G4String, 1822 fKnownMaterials->insert(std::pair<G4String,G4String>("Mineral_MI","FF_mineral_Tartari2002.dat"));
1713 fKnownMaterials->insert(std::pair<G4String, 1823 fKnownMaterials->insert(std::pair<G4String,G4String>("adipose_MI","FF_adipose_Poletti2002.dat"));
1714 fKnownMaterials->insert(std::pair<G4String, 1824 fKnownMaterials->insert(std::pair<G4String,G4String>("glandular_MI","FF_glandular_Poletti2002.dat"));
1715 fKnownMaterials->insert(std::pair<G4String, 1825 fKnownMaterials->insert(std::pair<G4String,G4String>("breast5050_MI","FF_human_breast_Peplow1998.dat"));
1716 fKnownMaterials->insert(std::pair<G4String, 1826 fKnownMaterials->insert(std::pair<G4String,G4String>("carcinoma_MI","FF_carcinoma_Kidane1999.dat"));
1717 fKnownMaterials->insert(std::pair<G4String, 1827 fKnownMaterials->insert(std::pair<G4String,G4String>("muscle_MI","FF_pork_muscle_Peplow1998.dat"));
1718 fKnownMaterials->insert(std::pair<G4String, 1828 fKnownMaterials->insert(std::pair<G4String,G4String>("kidney_MI","FF_pork_kidney_Peplow1998.dat"));
1719 fKnownMaterials->insert(std::pair<G4String, 1829 fKnownMaterials->insert(std::pair<G4String,G4String>("liver_MI","FF_pork_liver_Peplow1998.dat"));
1720 fKnownMaterials->insert(std::pair<G4String, 1830 fKnownMaterials->insert(std::pair<G4String,G4String>("heart_MI","FF_pork_heart_Peplow1998.dat"));
1721 fKnownMaterials->insert(std::pair<G4String, 1831 fKnownMaterials->insert(std::pair<G4String,G4String>("blood_MI","FF_beef_blood_Peplow1998.dat"));
1722 fKnownMaterials->insert(std::pair<G4String, 1832 fKnownMaterials->insert(std::pair<G4String,G4String>("grayMatter_MI","FF_gbrain_DeFelici2008.dat"));
1723 fKnownMaterials->insert(std::pair<G4String, 1833 fKnownMaterials->insert(std::pair<G4String,G4String>("whiteMatter_MI","FF_wbrain_DeFelici2008.dat"));
1724 fKnownMaterials->insert(std::pair<G4String, 1834 fKnownMaterials->insert(std::pair<G4String,G4String>("bone_MI","FF_bone_King2011.dat"));
1725 fKnownMaterials->insert(std::pair<G4String, 1835 fKnownMaterials->insert(std::pair<G4String,G4String>("FatLowX_MI","FF_fat_Tartari2002_joint_lowXdata_ESRF2003.dat"));
1726 fKnownMaterials->insert(std::pair<G4String, 1836 fKnownMaterials->insert(std::pair<G4String,G4String>("BoneMatrixLowX_MI","FF_bonematrix_Tartari2002_joint_lowXdata.dat"));
1727 fKnownMaterials->insert(std::pair<G4String, 1837 fKnownMaterials->insert(std::pair<G4String,G4String>("PMMALowX_MI", "FF_PMMA_Tartari2002_joint_lowXdata_ESRF2003.dat"));
1728 fKnownMaterials->insert(std::pair<G4String, 1838 fKnownMaterials->insert(std::pair<G4String,G4String>("dryBoneLowX_MI","FF_drybone_Tartari2002_joint_lowXdata_ESRF2003.dat"));
1729 fKnownMaterials->insert(std::pair<G4String, 1839 fKnownMaterials->insert(std::pair<G4String,G4String>("CIRS30-70_MI","FF_CIRS30-70_Poletti2002.dat"));
1730 fKnownMaterials->insert(std::pair<G4String, 1840 fKnownMaterials->insert(std::pair<G4String,G4String>("CIRS50-50_MI","FF_CIRS50-50_Poletti2002.dat"));
1731 fKnownMaterials->insert(std::pair<G4String, 1841 fKnownMaterials->insert(std::pair<G4String,G4String>("CIRS70-30_MI","FF_CIRS70-30_Poletti2002.dat"));
1732 fKnownMaterials->insert(std::pair<G4String, 1842 fKnownMaterials->insert(std::pair<G4String,G4String>("RMI454_MI", "FF_RMI454_Poletti2002.dat"));
1733 fKnownMaterials->insert(std::pair<G4String, 1843 fKnownMaterials->insert(std::pair<G4String,G4String>("PMMA_MI","FF_PMMA_Tartari2002.dat"));
1734 fKnownMaterials->insert(std::pair<G4String, 1844 fKnownMaterials->insert(std::pair<G4String,G4String>("Lexan_MI","FF_lexan_Peplow1998.dat"));
1735 fKnownMaterials->insert(std::pair<G4String, 1845 fKnownMaterials->insert(std::pair<G4String,G4String>("Kapton_MI","FF_kapton_Peplow1998.dat"));
1736 fKnownMaterials->insert(std::pair<G4String, 1846 fKnownMaterials->insert(std::pair<G4String,G4String>("Nylon_MI","FF_nylon_Kosanetzky1987.dat"));
1737 fKnownMaterials->insert(std::pair<G4String, 1847 fKnownMaterials->insert(std::pair<G4String,G4String>("Polyethylene_MI","FF_polyethylene_Kosanetzky1987.dat"));
1738 fKnownMaterials->insert(std::pair<G4String, 1848 fKnownMaterials->insert(std::pair<G4String,G4String>("Polystyrene_MI","FF_polystyrene_Kosanetzky1987.dat"));
1739 fKnownMaterials->insert(std::pair<G4String, 1849 fKnownMaterials->insert(std::pair<G4String,G4String>("Formaline_MI","FF_formaline_Peplow1998.dat"));
1740 fKnownMaterials->insert(std::pair<G4String, 1850 fKnownMaterials->insert(std::pair<G4String,G4String>("Acetone_MI","FF_acetone_Cozzini2010.dat"));
1741 fKnownMaterials->insert(std::pair<G4String, 1851 fKnownMaterials->insert(std::pair<G4String,G4String>("Hperoxide_MI","FF_Hperoxide_Cozzini2010.dat"));
>> 1852
1742 } 1853 }
1743 1854
1744 //....oooOO0OOooo........oooOO0OOooo........o 1855 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1745 1856
1746 G4double G4PenelopeRayleighModelMI::Integrate 1857 G4double G4PenelopeRayleighModelMI::IntegrateFun(G4double y[], G4int n, G4double dTheta)
1747 { 1858 {
1748 G4double integral = 0.; 1859 G4double integral = 0.;
1749 for (G4int k=0; k<n-1; k++) { 1860 for (G4int k=0; k<n-1; k++) {
1750 integral += (y[k]+y[k+1]); 1861 integral += (y[k]+y[k+1]);
1751 } 1862 }
1752 integral *= dTheta*0.5; 1863 integral *= dTheta*0.5;
1753 return integral; 1864 return integral;
1754 } 1865 }
1755 1866
1756 //....oooOO0OOooo........oooOO0OOooo........o 1867 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1757 G4MIData* G4PenelopeRayleighModelMI::GetMIDat 1868 G4MIData* G4PenelopeRayleighModelMI::GetMIData(const G4Material* material)
1758 { 1869 {
1759 if (material->IsExtended()) { 1870 if (material->IsExtended()) {
1760 G4ExtendedMaterial* aEM = (G4ExtendedMate 1871 G4ExtendedMaterial* aEM = (G4ExtendedMaterial*)material;
1761 G4MIData* dataMI = (G4MIData*)aEM->Retrie 1872 G4MIData* dataMI = (G4MIData*)aEM->RetrieveExtension("MI");
1762 return dataMI; 1873 return dataMI;
1763 } else { 1874 } else {
1764 return nullptr; 1875 return nullptr;
1765 } 1876 }
1766 } 1877 }
1767 1878