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