| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // Author: Luciano Pandola
28 //
29 // History:
30 // --------
31 // 03 Dec 2009 L Pandola First implementa
32 // 25 May 2011 L.Pandola Renamed (make v2
33 // 19 Sep 2013 L.Pandola Migration to MT
34 //
35
36 #include "G4PenelopeRayleighModel.hh"
37 #include "G4PhysicalConstants.hh"
38 #include "G4SystemOfUnits.hh"
39 #include "G4PenelopeSamplingData.hh"
40 #include "G4ParticleDefinition.hh"
41 #include "G4MaterialCutsCouple.hh"
42 #include "G4ProductionCutsTable.hh"
43 #include "G4DynamicParticle.hh"
44 #include "G4PhysicsTable.hh"
45 #include "G4ElementTable.hh"
46 #include "G4Element.hh"
47 #include "G4PhysicsFreeVector.hh"
48 #include "G4AutoLock.hh"
49 #include "G4Exp.hh"
50
51 //....oooOO0OOooo........oooOO0OOooo........oo
52
53 const G4int G4PenelopeRayleighModel::fMaxZ;
54 G4PhysicsFreeVector* G4PenelopeRayleighModel::
55 G4PhysicsFreeVector* G4PenelopeRayleighModel::
56
57 //....oooOO0OOooo........oooOO0OOooo........oo
58
59 G4PenelopeRayleighModel::G4PenelopeRayleighMod
60 const G4String& nam)
61 :G4VEmModel(nam),fParticleChange(nullptr),fP
62 fLogFormFactorTable(nullptr),fPMaxTable(nul
63 fIsInitialised(false),fLocalTable(false)
64 {
65 fIntrinsicLowEnergyLimit = 100.0*eV;
66 fIntrinsicHighEnergyLimit = 100.0*GeV;
67 // SetLowEnergyLimit(fIntrinsicLowEnergyLim
68 SetHighEnergyLimit(fIntrinsicHighEnergyLimit
69
70 if (part)
71 SetParticle(part);
72
73 //
74 fVerboseLevel= 0;
75 // Verbosity scale:
76 // 0 = nothing
77 // 1 = warning for energy non-conservation
78 // 2 = details of energy budget
79 // 3 = calculation of cross sections, file o
80 // 4 = entering in methods
81
82 //build the energy grid. It is the same for
83 G4double logenergy = G4Log(fIntrinsicLowEner
84 G4double logmaxenergy = G4Log(1.5*fIntrinsic
85 //finer grid below 160 keV
86 G4double logtransitionenergy = G4Log(160*keV
87 G4double logfactor1 = G4Log(10.)/250.;
88 G4double logfactor2 = logfactor1*10;
89 fLogEnergyGridPMax.push_back(logenergy);
90 do{
91 if (logenergy < logtransitionenergy)
92 logenergy += logfactor1;
93 else
94 logenergy += logfactor2;
95 fLogEnergyGridPMax.push_back(logenergy);
96 }while (logenergy < logmaxenergy);
97 }
98
99 //....oooOO0OOooo........oooOO0OOooo........oo
100
101 G4PenelopeRayleighModel::~G4PenelopeRayleighMo
102 {
103 if (IsMaster() || fLocalTable)
104 {
105
106 for(G4int i=0; i<=fMaxZ; ++i)
107 {
108 if(fLogAtomicCrossSection[i])
109 {
110 delete fLogAtomicCrossSection[i];
111 fLogAtomicCrossSection[i] = nullptr;
112 }
113 if(fAtomicFormFactor[i])
114 {
115 delete fAtomicFormFactor[i];
116 fAtomicFormFactor[i] = nullptr;
117 }
118 }
119 ClearTables();
120 }
121 }
122
123 //....oooOO0OOooo........oooOO0OOooo........oo
124 void G4PenelopeRayleighModel::ClearTables()
125 {
126 if (fLogFormFactorTable)
127 {
128 for (auto& item : (*fLogFormFactorTable
129 if (item.second) delete item.second;
130 delete fLogFormFactorTable;
131 fLogFormFactorTable = nullptr; //zero e
132 }
133 if (fPMaxTable)
134 {
135 for (auto& item : (*fPMaxTable))
136 if (item.second) delete item.second;
137 delete fPMaxTable;
138 fPMaxTable = nullptr; //zero explicitly
139 }
140 if (fSamplingTable)
141 {
142 for (auto& item : (*fSamplingTable))
143 if (item.second) delete item.second;
144 delete fSamplingTable;
145 fSamplingTable = nullptr; //zero explic
146 }
147 return;
148 }
149
150 //....oooOO0OOooo........oooOO0OOooo........oo
151
152 void G4PenelopeRayleighModel::Initialise(const
153 const G4DataVector& )
154 {
155 if (fVerboseLevel > 3)
156 G4cout << "Calling G4PenelopeRayleighModel
157
158 SetParticle(part);
159
160 //Only the master model creates/fills/destro
161 if (IsMaster() && part == fParticle)
162 {
163 //clear tables depending on materials, n
164 ClearTables();
165
166 if (fVerboseLevel > 3)
167 G4cout << "Calling G4PenelopeRayleighModel::
168
169 //create new tables
170 if (!fLogFormFactorTable)
171 fLogFormFactorTable = new std::map<const G4M
172 if (!fPMaxTable)
173 fPMaxTable = new std::map<const G4Material*,
174 if (!fSamplingTable)
175 fSamplingTable = new std::map<const G4Materi
176
177 G4ProductionCutsTable* theCoupleTable =
178 G4ProductionCutsTable::GetProductionCutsTabl
179
180 for (G4int i=0;i<(G4int)theCoupleTable->
181 {
182 const G4Material* material =
183 theCoupleTable->GetMaterialCutsCouple(i)
184 const G4ElementVector* theElementVector =
185
186 for (std::size_t j=0;j<material->GetNumber
187 {
188 G4int iZ = theElementVector->at(j)->Ge
189 //read data files only in the master
190 if (!fLogAtomicCrossSection[iZ])
191 ReadDataFile(iZ);
192 }
193
194 //1) If the table has not been built for t
195 if (!fLogFormFactorTable->count(material))
196 BuildFormFactorTable(material);
197
198 //2) retrieve or build the sampling table
199 if (!(fSamplingTable->count(material)))
200 InitializeSamplingAlgorithm(material);
201
202 //3) retrieve or build the pMax data
203 if (!fPMaxTable->count(material))
204 GetPMaxTable(material);
205 }
206
207 if (fVerboseLevel > 1) {
208 G4cout << "Penelope Rayleigh model v2008 is
209 << "Energy range: "
210 << LowEnergyLimit() / keV << " keV -
211 << HighEnergyLimit() / GeV << " GeV"
212 << G4endl;
213 }
214 }
215
216 if(fIsInitialised) return;
217 fParticleChange = GetParticleChangeForGamma(
218 fIsInitialised = true;
219 }
220
221 //....oooOO0OOooo........oooOO0OOooo........oo
222
223 void G4PenelopeRayleighModel::InitialiseLocal(
224 G4VEmModel *masterModel)
225 {
226 if (fVerboseLevel > 3)
227 G4cout << "Calling G4PenelopeRayleighMode
228 //
229 //Check that particle matches: one might hav
230 //for e+ and e-).
231 //
232 if (part == fParticle)
233 {
234 //Get the const table pointers from the
235 const G4PenelopeRayleighModel* theModel
236 static_cast<G4PenelopeRayleighModel*> (maste
237
238 //Copy pointers to the data tables
239 for(G4int i=0; i<=fMaxZ; ++i)
240 {
241 fLogAtomicCrossSection[i] = theModel->fLog
242 fAtomicFormFactor[i] = theModel->fAtomicFo
243 }
244 fLogFormFactorTable = theModel->fLogForm
245 fPMaxTable = theModel->fPMaxTable;
246 fSamplingTable = theModel->fSamplingTabl
247
248 //copy the G4DataVector with the grid
249 fLogQSquareGrid = theModel->fLogQSquareG
250
251 //Same verbosity for all workers, as the
252 fVerboseLevel = theModel->fVerboseLevel;
253 }
254
255 return;
256 }
257
258
259 //....oooOO0OOooo........oooOO0OOooo........oo
260 namespace { G4Mutex PenelopeRayleighModelMute
261 G4double G4PenelopeRayleighModel::ComputeCross
262 G4double energy,
263 G4double Z,
264 G4double,
265 G4double,
266 G4double)
267 {
268 // Cross section of Rayleigh scattering in P
269 // tabulation, Cuellen et al. (1997), with n
270 // et al. J. Phys. Chem. Ref. Data 4 (1975)
271
272 if (fVerboseLevel > 3)
273 G4cout << "Calling CrossSectionPerAtom() o
274
275 G4int iZ = G4int(Z);
276
277 if (!fLogAtomicCrossSection[iZ])
278 {
279 //If we are here, it means that Initial
280 //not filled up. This can happen in a U
281 if (fVerboseLevel > 0)
282 {
283 //Issue a G4Exception (warning) only in ve
284 G4ExceptionDescription ed;
285 ed << "Unable to retrieve the cross sectio
286 ed << "This can happen only in Unit Tests
287 G4Exception("G4PenelopeRayleighModel::Comp
288 "em2040",JustWarning,ed);
289 }
290 //protect file reading via autolock
291 G4AutoLock lock(&PenelopeRayleighModelM
292 ReadDataFile(iZ);
293 lock.unlock();
294 }
295
296 G4double cross = 0;
297 G4PhysicsFreeVector* atom = fLogAtomicCross
298 if (!atom)
299 {
300 G4ExceptionDescription ed;
301 ed << "Unable to find Z=" << iZ << " in
302 G4Exception("G4PenelopeRayleighModel::C
303 "em2041",FatalException,ed);
304 return 0;
305 }
306 G4double logene = G4Log(energy);
307 G4double logXS = atom->Value(logene);
308 cross = G4Exp(logXS);
309
310 if (fVerboseLevel > 2)
311 {
312 G4cout << "Rayleigh cross section at "
313 " = " << cross/barn << " barn" << G4endl;
314 }
315 return cross;
316 }
317
318
319 //....oooOO0OOooo........oooOO0OOooo........oo
320 void G4PenelopeRayleighModel::BuildFormFactorT
321 {
322 /*
323 1) get composition and equivalent molecula
324 */
325 std::size_t nElements = material->GetNumberO
326 const G4ElementVector* elementVector = mater
327 const G4double* fractionVector = material->G
328
329 std::vector<G4double> *StechiometricFactors
330 for (std::size_t i=0;i<nElements;++i)
331 {
332 G4double fraction = fractionVector[i];
333 G4double atomicWeigth = (*elementVector)
334 StechiometricFactors->push_back(fraction
335 }
336 //Find max
337 G4double MaxStechiometricFactor = 0.;
338 for (std::size_t i=0;i<nElements;++i)
339 {
340 if ((*StechiometricFactors)[i] > MaxStec
341 MaxStechiometricFactor = (*Stechiometr
342 }
343 if (MaxStechiometricFactor<1e-16)
344 {
345 G4ExceptionDescription ed;
346 ed << "Inconsistent data of atomic compo
347 material->GetName() << G4endl;
348 G4Exception("G4PenelopeRayleighModel::Bu
349 "em2042",FatalException,ed);
350 }
351 //Normalize
352 for (std::size_t i=0;i<nElements;++i)
353 (*StechiometricFactors)[i] /= MaxStechiom
354
355 /*
356 CREATE THE FORM FACTOR TABLE
357 */
358 G4PhysicsFreeVector* theFFVec = new G4Physic
359
360 for (std::size_t k=0;k<fLogQSquareGrid.size(
361 {
362 G4double ff2 = 0; //squared form factor
363 for (std::size_t i=0;i<nElements;++i)
364 {
365 G4int iZ = (*elementVector)[i]->GetZasInt(
366 G4PhysicsFreeVector* theAtomVec = fAtomicF
367 G4double f = (*theAtomVec)[k]; //the q-gri
368 ff2 += f*f*(*StechiometricFactors)[i];
369 }
370 if (ff2)
371 theFFVec->PutValue(k,fLogQSquareGrid[k],G4Lo
372 }
373 theFFVec->FillSecondDerivatives(); //vector
374 fLogFormFactorTable->insert(std::make_pair(m
375
376 delete StechiometricFactors;
377 return;
378 }
379
380 //....oooOO0OOooo........oooOO0OOooo........oo
381
382 void G4PenelopeRayleighModel::SampleSecondarie
383 const G4MaterialCutsCouple* couple
384 const G4DynamicParticle* aDynamicG
385 G4double,
386 G4double)
387 {
388 // Sampling of the Rayleigh final state (nam
389 // from the Penelope2008 model. The scatteri
390 // cross section dOmega/d(cosTheta) from Bor
391 // anomalous scattering effects. The Form Fa
392 // analytical cross section is retrieved via
393 // are tabulated for F(Q). Form factor for c
394 // the additivity rule. The sampling from th
395 // Transform with Aliasing (RITA) algorithm;
396 // for each material and managed by G4Penelo
397 // The sampling algorithm (rejection method)
398 // increases with energy. For E=100 keV the
399 // hydrogen and uranium, respectively.
400
401 if (fVerboseLevel > 3)
402 G4cout << "Calling SamplingSecondaries() o
403
404 G4double photonEnergy0 = aDynamicGamma->GetK
405
406 if (photonEnergy0 <= fIntrinsicLowEnergyLimi
407 {
408 fParticleChange->ProposeTrackStatus(fSto
409 fParticleChange->SetProposedKineticEnerg
410 fParticleChange->ProposeLocalEnergyDepos
411 return ;
412 }
413
414 G4ParticleMomentum photonDirection0 = aDynam
415
416 const G4Material* theMat = couple->GetMateri
417
418 //1) Verify if tables are ready
419 //Either Initialize() was not called, or we
420 //not invoked
421 if (!fPMaxTable || !fSamplingTable || !fLogF
422 {
423 //create a **thread-local** version of t
424 //Unit Tests
425 fLocalTable = true;
426 if (!fLogFormFactorTable)
427 fLogFormFactorTable = new std::map<const G4M
428 if (!fPMaxTable)
429 fPMaxTable = new std::map<const G4Material*,
430 if (!fSamplingTable)
431 fSamplingTable = new std::map<const G4Materi
432 }
433
434 if (!fSamplingTable->count(theMat))
435 {
436 //If we are here, it means that Initiali
437 //not filled up. This can happen in a Un
438 if (fVerboseLevel > 0)
439 {
440 //Issue a G4Exception (warning) only in ve
441 G4ExceptionDescription ed;
442 ed << "Unable to find the fSamplingTable d
443 theMat->GetName() << G4endl;
444 ed << "This can happen only in Unit Tests"
445 G4Exception("G4PenelopeRayleighModel::Samp
446 "em2019",JustWarning,ed);
447 }
448 const G4ElementVector* theElementVector
449 //protect file reading via autolock
450 G4AutoLock lock(&PenelopeRayleighModelMu
451 for (std::size_t j=0;j<theMat->GetNumber
452 {
453 G4int iZ = theElementVector->at(j)->GetZas
454 if (!fLogAtomicCrossSection[iZ])
455 {
456 lock.lock();
457 ReadDataFile(iZ);
458 lock.unlock();
459 }
460 }
461 lock.lock();
462 //1) If the table has not been built for
463 if (!fLogFormFactorTable->count(theMat))
464 BuildFormFactorTable(theMat);
465
466 //2) retrieve or build the sampling tabl
467 if (!(fSamplingTable->count(theMat)))
468 InitializeSamplingAlgorithm(theMat);
469
470 //3) retrieve or build the pMax data
471 if (!fPMaxTable->count(theMat))
472 GetPMaxTable(theMat);
473 lock.unlock();
474 }
475
476 //Ok, restart the job
477 G4PenelopeSamplingData* theDataTable = fSamp
478 G4PhysicsFreeVector* thePMax = fPMaxTable->f
479
480 G4double cosTheta = 1.0;
481
482 //OK, ready to go!
483 G4double qmax = 2.0*photonEnergy0/electron_m
484
485 if (qmax < 1e-10) //very low momentum transf
486 {
487 G4bool loopAgain=false;
488 do
489 {
490 loopAgain = false;
491 cosTheta = 1.0-2.0*G4UniformRand();
492 G4double G = 0.5*(1+cosTheta*cosTheta);
493 if (G4UniformRand()>G)
494 loopAgain = true;
495 }while(loopAgain);
496 }
497 else //larger momentum transfer
498 {
499 std::size_t nData = theDataTable->GetNum
500 G4double LastQ2inTheTable = theDataTable
501 G4double q2max = std::min(qmax*qmax,Last
502
503 G4bool loopAgain = false;
504 G4double MaxPValue = thePMax->Value(phot
505 G4double xx=0;
506
507 //Sampling by rejection method. The reje
508 //G = 0.5*(1+cos^2(theta))
509 //
510 do{
511 loopAgain = false;
512 G4double RandomMax = G4UniformRand()*MaxPVal
513 xx = theDataTable->SampleValue(RandomMax);
514 //xx is a random value of q^2 in (0,q2max),s
515 //F(Q^2) via the RITA algorithm
516 if (xx > q2max)
517 loopAgain = true;
518 cosTheta = 1.0-2.0*xx/q2max;
519 G4double G = 0.5*(1+cosTheta*cosTheta);
520 if (G4UniformRand()>G)
521 loopAgain = true;
522 }while(loopAgain);
523 }
524
525 G4double sinTheta = std::sqrt(1-cosTheta*cos
526
527 // Scattered photon angles. ( Z - axis along
528 G4double phi = twopi * G4UniformRand() ;
529 G4double dirX = sinTheta*std::cos(phi);
530 G4double dirY = sinTheta*std::sin(phi);
531 G4double dirZ = cosTheta;
532
533 // Update G4VParticleChange for the scattere
534 G4ThreeVector photonDirection1(dirX, dirY, d
535
536 photonDirection1.rotateUz(photonDirection0);
537 fParticleChange->ProposeMomentumDirection(ph
538 fParticleChange->SetProposedKineticEnergy(ph
539
540 return;
541 }
542
543
544 //....oooOO0OOooo........oooOO0OOooo........oo
545
546 void G4PenelopeRayleighModel::ReadDataFile(con
547 {
548 if (fVerboseLevel > 2)
549 {
550 G4cout << "G4PenelopeRayleighModel::Read
551 G4cout << "Going to read Rayleigh data f
552 }
553 const char* path = G4FindDataDir("G4LEDATA
554 if(!path)
555 {
556 G4String excep = "G4LEDATA environment v
557 G4Exception("G4PenelopeRayleighModel::Re
558 "em0006",FatalException,excep);
559 return;
560 }
561
562 /*
563 Read first the cross section file
564 */
565 std::ostringstream ost;
566 if (Z>9)
567 ost << path << "/penelope/rayleigh/pdgra"
568 else
569 ost << path << "/penelope/rayleigh/pdgra0"
570 std::ifstream file(ost.str().c_str());
571 if (!file.is_open())
572 {
573 G4String excep = "Data file " + G4String
574 G4Exception("G4PenelopeRayleighModel::Re
575 "em0003",FatalException,excep);
576 }
577 G4int readZ =0;
578 std::size_t nPoints= 0;
579 file >> readZ >> nPoints;
580 //check the right file is opened.
581 if (readZ != Z || nPoints <= 0 || nPoints >=
582 {
583 G4ExceptionDescription ed;
584 ed << "Corrupted data file for Z=" << Z
585 G4Exception("G4PenelopeRayleighModel::Re
586 "em0005",FatalException,ed);
587 return;
588 }
589
590 fLogAtomicCrossSection[Z] = new G4PhysicsFre
591 G4double ene=0,f1=0,f2=0,xs=0;
592 for (std::size_t i=0;i<nPoints;++i)
593 {
594 file >> ene >> f1 >> f2 >> xs;
595 //dimensional quantities
596 ene *= eV;
597 xs *= cm2;
598 fLogAtomicCrossSection[Z]->PutValue(i,G4
599 if (file.eof() && i != (nPoints-1)) //fi
600 {
601 G4ExceptionDescription ed ;
602 ed << "Corrupted data file for Z=" << Z <<
603 ed << "Found less than " << nPoints << "en
604 G4Exception("G4PenelopeRayleighModel::Read
605 "em0005",FatalException,ed);
606 }
607 }
608 file.close();
609
610 /*
611 Then read the form factor file
612 */
613 std::ostringstream ost2;
614 if (Z>9)
615 ost2 << path << "/penelope/rayleigh/pdaff"
616 else
617 ost2 << path << "/penelope/rayleigh/pdaff0
618 file.open(ost2.str().c_str());
619 if (!file.is_open())
620 {
621 G4String excep = "Data file " + G4String
622 G4Exception("G4PenelopeRayleighModel::Re
623 "em0003",FatalException,excep);
624 }
625 file >> readZ >> nPoints;
626 //check the right file is opened.
627 if (readZ != Z || nPoints <= 0 || nPoints >=
628 {
629 G4ExceptionDescription ed;
630 ed << "Corrupted data file for Z=" << Z
631 G4Exception("G4PenelopeRayleighModel::Re
632 "em0005",FatalException,ed);
633 return;
634 }
635 fAtomicFormFactor[Z] = new G4PhysicsFreeVect
636 G4double q=0,ff=0,incoh=0;
637 G4bool fillQGrid = false;
638 //fill this vector only the first time.
639 if (!fLogQSquareGrid.size())
640 fillQGrid = true;
641 for (std::size_t i=0;i<nPoints;++i)
642 {
643 file >> q >> ff >> incoh;
644 //q and ff are dimensionless (q is in un
645 fAtomicFormFactor[Z]->PutValue(i,q,ff);
646 if (fillQGrid)
647 {
648 fLogQSquareGrid.push_back(2.0*G4Log(q));
649 }
650 if (file.eof() && i != (nPoints-1)) //fi
651 {
652 G4ExceptionDescription ed;
653 ed << "Corrupted data file for Z=" << Z <<
654 ed << "Found less than " << nPoints << "en
655 G4Exception("G4PenelopeRayleighModel::Read
656 "em0005",FatalException,ed);
657 }
658 }
659 file.close();
660 return;
661 }
662
663 //....oooOO0OOooo........oooOO0OOooo........oo
664
665 G4double G4PenelopeRayleighModel::GetFSquared(
666 {
667 G4double f2 = 0;
668 //Input value QSquared could be zero: protec
669 //the FPE exception
670 //If Q<1e-10, set Q to 1e-10
671 G4double logQSquared = (QSquared>1e-10) ? G4
672 //last value of the table
673 G4double maxlogQ2 = fLogQSquareGrid[fLogQSqu
674
675 //now it should be all right
676 G4PhysicsFreeVector* theVec = fLogFormFactor
677
678 if (!theVec)
679 {
680 G4ExceptionDescription ed;
681 ed << "Unable to retrieve F squared tabl
682 G4Exception("G4PenelopeRayleighModel::Ge
683 "em2046",FatalException,ed);
684 return 0;
685 }
686 if (logQSquared < -20) // Q < 1e-9
687 {
688 G4double logf2 = (*theVec)[0]; //first v
689 f2 = G4Exp(logf2);
690 }
691 else if (logQSquared > maxlogQ2)
692 f2 =0;
693 else
694 {
695 //log(Q^2) vs. log(F^2)
696 G4double logf2 = theVec->Value(logQSquar
697 f2 = G4Exp(logf2);
698
699 }
700 if (fVerboseLevel > 3)
701 {
702 G4cout << "G4PenelopeRayleighModel::GetF
703 G4cout << "Q^2 = " << QSquared << " (un
704 }
705 return f2;
706 }
707
708 //....oooOO0OOooo........oooOO0OOooo........oo
709
710 void G4PenelopeRayleighModel::InitializeSampli
711 {
712 G4double q2min = 0;
713 G4double q2max = 0;
714 const std::size_t np = 150; //hard-coded in
715 //G4cout << "Init N= " << fLogQSquareGrid.si
716 for (std::size_t i=1;i<fLogQSquareGrid.size(
717 {
718 G4double Q2 = G4Exp(fLogQSquareGrid[i]);
719 if (GetFSquared(mat,Q2) > 1e-35)
720 {
721 q2max = G4Exp(fLogQSquareGrid[i-1]);
722 }
723 //G4cout << "Q2= " << Q2 << " q2max= " <
724 }
725
726 std::size_t nReducedPoints = np/4;
727
728 //check for errors
729 if (np < 16)
730 {
731 G4Exception("G4PenelopeRayleighModel::In
732 "em2047",FatalException,
733 "Too few points to initialize the sampli
734 }
735 if (q2min > (q2max-1e-10))
736 {
737 G4cout << "q2min= " << q2min << " q2max=
738 G4Exception("G4PenelopeRayleighModel::In
739 "em2048",FatalException,
740 "Too narrow grid to initialize the sampl
741 }
742
743 //This is subroutine RITAI0 of Penelope
744 //Create an object of type G4PenelopeRayleig
745
746 //temporary vectors --> Then everything is s
747 G4DataVector* x = new G4DataVector();
748
749 /*******************************************
750 Start with a grid of NUNIF points uniforml
751 ********************************************
752 std::size_t NUNIF = std::min(std::max(((std:
753 const G4int nip = 51; //hard-coded in Penelo
754
755 G4double dx = (q2max-q2min)/((G4double) NUNI
756 x->push_back(q2min);
757 for (std::size_t i=1;i<NUNIF-1;++i)
758 {
759 G4double app = q2min + i*dx;
760 x->push_back(app); //increase
761 }
762 x->push_back(q2max);
763
764 if (fVerboseLevel> 3)
765 G4cout << "Vector x has " << x->size() <<
766
767 //Allocate temporary storage vectors
768 G4DataVector* area = new G4DataVector();
769 G4DataVector* a = new G4DataVector();
770 G4DataVector* b = new G4DataVector();
771 G4DataVector* c = new G4DataVector();
772 G4DataVector* err = new G4DataVector();
773
774 for (std::size_t i=0;i<NUNIF-1;++i) //build
775 {
776 //Temporary vectors for this loop
777 G4DataVector* pdfi = new G4DataVector();
778 G4DataVector* pdfih = new G4DataVector()
779 G4DataVector* sumi = new G4DataVector();
780
781 G4double dxi = ((*x)[i+1]-(*x)[i])/(G4do
782 G4double pdfmax = 0;
783 for (G4int k=0;k<nip;k++)
784 {
785 G4double xik = (*x)[i]+k*dxi;
786 G4double pdfk = std::max(GetFSquared(mat,x
787 pdfi->push_back(pdfk);
788 pdfmax = std::max(pdfmax,pdfk);
789 if (k < (nip-1))
790 {
791 G4double xih = xik + 0.5*dxi;
792 G4double pdfIK = std::max(GetFSquared(
793 pdfih->push_back(pdfIK);
794 pdfmax = std::max(pdfmax,pdfIK);
795 }
796 }
797
798 //Simpson's integration
799 G4double cons = dxi*0.5*(1./3.);
800 sumi->push_back(0.);
801 for (G4int k=1;k<nip;k++)
802 {
803 G4double previous = (*sumi)[k-1];
804 G4double next = previous + cons*((*pdfi)[k
805 sumi->push_back(next);
806 }
807
808 G4double lastIntegral = (*sumi)[sumi->si
809 area->push_back(lastIntegral);
810 //Normalize cumulative function
811 G4double factor = 1.0/lastIntegral;
812 for (std::size_t k=0;k<sumi->size();++k)
813 (*sumi)[k] *= factor;
814
815 //When the PDF vanishes at one of the in
816 if ((*pdfi)[0] < 1e-35)
817 (*pdfi)[0] = 1e-5*pdfmax;
818 if ((*pdfi)[pdfi->size()-1] < 1e-35)
819 (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax;
820
821 G4double pli = (*pdfi)[0]*factor;
822 G4double pui = (*pdfi)[pdfi->size()-1]*f
823 G4double B_temp = 1.0-1.0/(pli*pui*dx*dx
824 G4double A_temp = (1.0/(pli*dx))-1.0-B_t
825 G4double C_temp = 1.0+A_temp+B_temp;
826 if (C_temp < 1e-35)
827 {
828 a->push_back(0.);
829 b->push_back(0.);
830 c->push_back(1.);
831 }
832 else
833 {
834 a->push_back(A_temp);
835 b->push_back(B_temp);
836 c->push_back(C_temp);
837 }
838
839 //OK, now get ERR(I), the integral of th
840 //and the true pdf, extended over the in
841 G4int icase = 1; //loop code
842 G4bool reLoop = false;
843 err->push_back(0.);
844 do
845 {
846 reLoop = false;
847 (*err)[i] = 0.; //zero variable
848 for (G4int k=0;k<nip;k++)
849 {
850 G4double rr = (*sumi)[k];
851 G4double pap = (*area)[i]*(1.0+((*a)[i
852 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(*
853 if (k == 0 || k == nip-1)
854 (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k])
855 else
856 (*err)[i] += std::fabs(pap-(*pdfi)[k]);
857 }
858 (*err)[i] *= dxi;
859
860 //If err(I) is too large, the pdf is appro
861 if ((*err)[i] > 0.1*(*area)[i] && icase ==
862 {
863 (*b)[i] = 0;
864 (*a)[i] = 0;
865 (*c)[i] = 1.;
866 icase = 2;
867 reLoop = true;
868 }
869 }while(reLoop);
870 delete pdfi;
871 delete pdfih;
872 delete sumi;
873 } //end of first loop over i
874
875 //Now assign last point
876 (*x)[x->size()-1] = q2max;
877 a->push_back(0.);
878 b->push_back(0.);
879 c->push_back(0.);
880 err->push_back(0.);
881 area->push_back(0.);
882
883 if (x->size() != NUNIF || a->size() != NUNIF
884 err->size() != NUNIF || area->size() !=
885 {
886 G4ExceptionDescription ed;
887 ed << "Problem in building the Table for
888 G4Exception("G4PenelopeRayleighModel::In
889 "em2049",FatalException,ed);
890 }
891
892 /*******************************************
893 New grid points are added by halving the su
894 This is done up to np=150 points in the grid
895 ********************************************
896 do
897 {
898 G4double maxError = 0.0;
899 std::size_t iErrMax = 0;
900 for (std::size_t i=0;i<err->size()-2;++i
901 {
902 //maxError is the lagest of the interval e
903 if ((*err)[i] > maxError)
904 {
905 maxError = (*err)[i];
906 iErrMax = i;
907 }
908 }
909
910 //OK, now I have to insert one new point
911 G4double newx = 0.5*((*x)[iErrMax]+(*x)[
912
913 x->insert(x->begin()+iErrMax+1,newx);
914 //Add place-holders in the other vectors
915 area->insert(area->begin()+iErrMax+1,0.)
916 a->insert(a->begin()+iErrMax+1,0.);
917 b->insert(b->begin()+iErrMax+1,0.);
918 c->insert(c->begin()+iErrMax+1,0.);
919 err->insert(err->begin()+iErrMax+1,0.);
920
921 //Now calculate the other parameters
922 for (std::size_t i=iErrMax;i<=iErrMax+1;
923 {
924 //Temporary vectors for this loop
925 G4DataVector* pdfi = new G4DataVector();
926 G4DataVector* pdfih = new G4DataVector();
927 G4DataVector* sumi = new G4DataVector();
928
929 G4double dxLocal = (*x)[i+1]-(*x)[i];
930 G4double dxi = ((*x)[i+1]-(*x)[i])/(G4doub
931 G4double pdfmax = 0;
932 for (G4int k=0;k<nip;k++)
933 {
934 G4double xik = (*x)[i]+k*dxi;
935 G4double pdfk = std::max(GetFSquared(m
936 pdfi->push_back(pdfk);
937 pdfmax = std::max(pdfmax,pdfk);
938 if (k < (nip-1))
939 {
940 G4double xih = xik + 0.5*dxi;
941 G4double pdfIK = std::max(GetFSquared(ma
942 pdfih->push_back(pdfIK);
943 pdfmax = std::max(pdfmax,pdfIK);
944 }
945 }
946
947 //Simpson's integration
948 G4double cons = dxi*0.5*(1./3.);
949 sumi->push_back(0.);
950 for (G4int k=1;k<nip;k++)
951 {
952 G4double previous = (*sumi)[k-1];
953 G4double next = previous + cons*((*pdf
954 sumi->push_back(next);
955 }
956 G4double lastIntegral = (*sumi)[sumi->size
957 (*area)[i] = lastIntegral;
958
959 //Normalize cumulative function
960 G4double factor = 1.0/lastIntegral;
961 for (std::size_t k=0;k<sumi->size();++k)
962 (*sumi)[k] *= factor;
963
964 //When the PDF vanishes at one of the inte
965 if ((*pdfi)[0] < 1e-35)
966 (*pdfi)[0] = 1e-5*pdfmax;
967 if ((*pdfi)[pdfi->size()-1] < 1e-35)
968 (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax;
969
970 G4double pli = (*pdfi)[0]*factor;
971 G4double pui = (*pdfi)[pdfi->size()-1]*fac
972 G4double B_temp = 1.0-1.0/(pli*pui*dxLocal
973 G4double A_temp = (1.0/(pli*dxLocal))-1.0-
974 G4double C_temp = 1.0+A_temp+B_temp;
975 if (C_temp < 1e-35)
976 {
977 (*a)[i]= 0.;
978 (*b)[i] = 0.;
979 (*c)[i] = 1;
980 }
981 else
982 {
983 (*a)[i]= A_temp;
984 (*b)[i] = B_temp;
985 (*c)[i] = C_temp;
986 }
987 //OK, now get ERR(I), the integral of the
988 //and the true pdf, extended over the inte
989 G4int icase = 1; //loop code
990 G4bool reLoop = false;
991 do
992 {
993 reLoop = false;
994 (*err)[i] = 0.; //zero variable
995 for (G4int k=0;k<nip;k++)
996 {
997 G4double rr = (*sumi)[k];
998 G4double pap = (*area)[i]*(1.0+((*a)[i]+
999 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1
1000 if (k == 0 || k == nip-1)
1001 (*err)[i] += 0.5*std::fabs(pap-(*pdfi
1002 else
1003 (*err)[i] += std::fabs(pap-(*pdfi)[k]
1004 }
1005 (*err)[i] *= dxi;
1006
1007 //If err(I) is too large, the pdf is
1008 if ((*err)[i] > 0.1*(*area)[i] && ica
1009 {
1010 (*b)[i] = 0;
1011 (*a)[i] = 0;
1012 (*c)[i] = 1.;
1013 icase = 2;
1014 reLoop = true;
1015 }
1016 }while(reLoop);
1017 delete pdfi;
1018 delete pdfih;
1019 delete sumi;
1020 }
1021 }while(x->size() < np);
1022
1023 if (x->size() != np || a->size() != np ||
1024 err->size() != np || area->size() != np
1025 {
1026 G4Exception("G4PenelopeRayleighModel::I
1027 "em2050",FatalException,
1028 "Problem in building the extended Table
1029 }
1030
1031 /******************************************
1032 Renormalization
1033 *******************************************
1034 G4double ws = 0;
1035 for (std::size_t i=0;i<np-1;++i)
1036 ws += (*area)[i];
1037 ws = 1.0/ws;
1038 G4double errMax = 0;
1039 for (std::size_t i=0;i<np-1;++i)
1040 {
1041 (*area)[i] *= ws;
1042 (*err)[i] *= ws;
1043 errMax = std::max(errMax,(*err)[i]);
1044 }
1045
1046 //Vector with the normalized cumulative dis
1047 G4DataVector* PAC = new G4DataVector();
1048 PAC->push_back(0.);
1049 for (std::size_t i=0;i<np-1;++i)
1050 {
1051 G4double previous = (*PAC)[i];
1052 PAC->push_back(previous+(*area)[i]);
1053 }
1054 (*PAC)[PAC->size()-1] = 1.;
1055
1056 /******************************************
1057 Pre-calculated limits for the initial binar
1058 *******************************************
1059 std::vector<std::size_t> *ITTL = new std::v
1060 std::vector<std::size_t> *ITTU = new std::v
1061
1062 //Just create place-holders
1063 for (std::size_t i=0;i<np;++i)
1064 {
1065 ITTL->push_back(0);
1066 ITTU->push_back(0);
1067 }
1068
1069 G4double bin = 1.0/(np-1);
1070 (*ITTL)[0]=0;
1071 for (std::size_t i=1;i<(np-1);++i)
1072 {
1073 G4double ptst = i*bin;
1074 G4bool found = false;
1075 for (std::size_t j=(*ITTL)[i-1];j<np &&
1076 {
1077 if ((*PAC)[j] > ptst)
1078 {
1079 (*ITTL)[i] = j-1;
1080 (*ITTU)[i-1] = j;
1081 found = true;
1082 }
1083 }
1084 }
1085 (*ITTU)[ITTU->size()-2] = ITTU->size()-1;
1086 (*ITTU)[ITTU->size()-1] = ITTU->size()-1;
1087 (*ITTL)[ITTL->size()-1] = ITTU->size()-2;
1088
1089 if (ITTU->size() != np || ITTU->size() != n
1090 {
1091 G4Exception("G4PenelopeRayleighModel::I
1092 "em2051",FatalException,
1093 "Problem in building the Limit Tables f
1094 }
1095
1096 /******************************************
1097 Copy tables
1098 *******************************************
1099 G4PenelopeSamplingData* theTable = new G4Pe
1100 for (std::size_t i=0;i<np;++i)
1101 {
1102 theTable->AddPoint((*x)[i],(*PAC)[i],(*
1103 }
1104
1105 if (fVerboseLevel > 2)
1106 {
1107 G4cout << "****************************
1108 G4endl;
1109 G4cout << "Sampling table for Penelope
1110 theTable->DumpTable();
1111 }
1112 fSamplingTable->insert(std::make_pair(mat,t
1113
1114 //Clean up temporary vectors
1115 delete x;
1116 delete a;
1117 delete b;
1118 delete c;
1119 delete err;
1120 delete area;
1121 delete PAC;
1122 delete ITTL;
1123 delete ITTU;
1124
1125 //DONE!
1126 return;
1127 }
1128
1129 //....oooOO0OOooo........oooOO0OOooo........o
1130
1131 void G4PenelopeRayleighModel::GetPMaxTable(co
1132 {
1133 if (!fPMaxTable)
1134 {
1135 G4cout << "G4PenelopeRayleighModel::Bui
1136 G4cout << "Going to instanziate the fPM
1137 G4cout << "That should _not_ be here! "
1138 fPMaxTable = new std::map<const G4Mater
1139 }
1140 //check if the table is already there
1141 if (fPMaxTable->count(mat))
1142 return;
1143
1144 //otherwise build it
1145 if (!fSamplingTable)
1146 {
1147 G4Exception("G4PenelopeRayleighModel::G
1148 "em2052",FatalException,
1149 "SamplingTable is not properly instanti
1150 return;
1151 }
1152
1153 //This should not be: the sampling table is
1154 if (!fSamplingTable->count(mat))
1155 {
1156 G4ExceptionDescription ed;
1157 ed << "Sampling table for material " <
1158 G4Exception("G4PenelopeRayleighModel::
1159 "em2052",FatalException,
1160 ed);
1161 return;
1162 }
1163
1164 G4PenelopeSamplingData *theTable = fSamplin
1165 std::size_t tablePoints = theTable->GetNumb
1166
1167 std::size_t nOfEnergyPoints = fLogEnergyGri
1168 G4PhysicsFreeVector* theVec = new G4Physics
1169
1170 const std::size_t nip = 51; //hard-coded in
1171
1172 for (std::size_t ie=0;ie<fLogEnergyGridPMax
1173 {
1174 G4double energy = G4Exp(fLogEnergyGridP
1175 G4double Qm = 2.0*energy/electron_mass_
1176 G4double Qm2 = Qm*Qm;
1177 G4double firstQ2 = theTable->GetX(0);
1178 G4double lastQ2 = theTable->GetX(tableP
1179 G4double thePMax = 0;
1180
1181 if (Qm2 > firstQ2)
1182 {
1183 if (Qm2 < lastQ2)
1184 {
1185 //bisection to look for the index of
1186 std::size_t lowerBound = 0;
1187 std::size_t upperBound = tablePoints-
1188 while (lowerBound <= upperBound)
1189 {
1190 std::size_t midBin = (lowerBound + uppe
1191 if( Qm2 < theTable->GetX(midBin))
1192 { upperBound = midBin-1; }
1193 else
1194 { lowerBound = midBin+1; }
1195 }
1196 //upperBound is the output (but also
1197 G4double Q1 = theTable->GetX(upperBou
1198 G4double Q2 = Qm2;
1199 G4double DQ = (Q2-Q1)/((G4double)(nip
1200 G4double theA = theTable->GetA(upperB
1201 G4double theB = theTable->GetB(upperB
1202 G4double thePAC = theTable->GetPAC(up
1203 G4DataVector* fun = new G4DataVector(
1204 for (std::size_t k=0;k<nip;++k)
1205 {
1206 G4double qi = Q1 + k*DQ;
1207 G4double tau = (qi-Q1)/
1208 (theTable->GetX(upperBound+1)-Q1);
1209 G4double con1 = 2.0*theB*tau;
1210 G4double ci = 1.0+theA+theB;
1211 G4double con2 = ci-theA*tau;
1212 G4double etap = 0;
1213 if (std::fabs(con1) > 1.0e-16*std::fabs
1214 etap = con2*(1.0-std::sqrt(1.0-2.0*ta
1215 else
1216 etap = tau/con2;
1217 G4double theFun = (theTable->GetPAC(upp
1218 (1.0+(theA+theB*etap)*etap)*(1.0+(the
1219 ((1.0-theB*etap*etap)*ci*(theTable->G
1220 fun->push_back(theFun);
1221 }
1222 //Now intergrate numerically the fun
1223 G4DataVector* sum = new G4DataVector;
1224 G4double CONS = DQ*(1./12.);
1225 G4double HCONS = 0.5*CONS;
1226 sum->push_back(0.);
1227 G4double secondPoint = (*sum)[0] +
1228 (5.0*(*fun)[0]+8.0*(*fun)[1]-(*fun)[2])*C
1229 sum->push_back(secondPoint);
1230 for (std::size_t hh=2;hh<nip-1;++hh)
1231 {
1232 G4double previous = (*sum)[hh-1];
1233 G4double next = previous+(13.0*((*fun)[
1234 (*fun)[hh+1]-(*fun)[hh-2])*HCON
1235 sum->push_back(next);
1236 }
1237 G4double last = (*sum)[nip-2]+(5.0*(*
1238 (*fun)[nip-3])*CONS;
1239 sum->push_back(last);
1240 thePMax = thePAC + (*sum)[sum->size()
1241 delete fun;
1242 delete sum;
1243 }
1244 else
1245 {
1246 thePMax = 1.0;
1247 }
1248 }
1249 else
1250 {
1251 thePMax = theTable->GetPAC(0);
1252 }
1253
1254 //Write number in the table
1255 theVec->PutValue(ie,energy,thePMax);
1256 }
1257
1258 fPMaxTable->insert(std::make_pair(mat,theVe
1259 return;
1260 }
1261
1262 //....oooOO0OOooo........oooOO0OOooo........o
1263
1264 void G4PenelopeRayleighModel::DumpFormFactorT
1265 {
1266 G4cout << "********************************
1267 G4cout << "G4PenelopeRayleighModel: Form Fa
1268 //try to use the same format as Penelope-Fo
1269 G4cout << "Q/(m_e*c) F(Q)
1270 G4cout << "********************************
1271 if (!fLogFormFactorTable->count(mat))
1272 BuildFormFactorTable(mat);
1273
1274 G4PhysicsFreeVector* theVec = fLogFormFacto
1275 for (std::size_t i=0;i<theVec->GetVectorLen
1276 {
1277 G4double logQ2 = theVec->GetLowEdgeEner
1278 G4double Q = G4Exp(0.5*logQ2);
1279 G4double logF2 = (*theVec)[i];
1280 G4double F = G4Exp(0.5*logF2);
1281 G4cout << Q << " " << F <<
1282 }
1283 //DONE
1284 return;
1285 }
1286
1287 //....oooOO0OOooo........oooOO0OOooo........o
1288
1289 void G4PenelopeRayleighModel::SetParticle(con
1290 {
1291 if(!fParticle) {
1292 fParticle = p;
1293 }
1294 }
1295