Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/pii/src/G4PixeCrossSectionHandler.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 //
 27 //
 28 // Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
 29 //
 30 // History:
 31 // -----------
 32 // 16 Jun 2008 MGP   Created; Cross section manager for hadron impact ionization
 33 //                   Documented in:
 34 //                   M.G. Pia et al., PIXE Simulation With Geant4,
 35 //                   IEEE Trans. Nucl. Sci., vol. 56, no. 6, pp. 3614-3649, Dec. 2009
 36 //
 37 // -------------------------------------------------------------------
 38 
 39 #include "G4PixeCrossSectionHandler.hh"
 40 #include "G4PhysicalConstants.hh"
 41 #include "G4IInterpolator.hh"
 42 #include "G4LogLogInterpolator.hh"
 43 #include "G4IDataSet.hh"
 44 #include "G4DataSet.hh"
 45 #include "G4CompositeDataSet.hh"
 46 #include "G4PixeShellDataSet.hh"
 47 #include "G4ProductionCutsTable.hh"
 48 #include "G4Material.hh"
 49 #include "G4Element.hh"
 50 #include "Randomize.hh"
 51 #include "G4SystemOfUnits.hh"
 52 #include "G4ParticleDefinition.hh"
 53 
 54 #include <map>
 55 #include <vector>
 56 #include <fstream>
 57 #include <sstream>
 58 #include <memory>
 59 
 60 
 61 G4PixeCrossSectionHandler::G4PixeCrossSectionHandler()
 62 {
 63   crossSections = 0;
 64   interpolation = 0;
 65   // Initialise with default values
 66   Initialise(0,"","","",1.*keV,0.1*GeV,200,MeV,barn,6,92);
 67   ActiveElements();
 68 }
 69 
 70 
 71 G4PixeCrossSectionHandler::G4PixeCrossSectionHandler(G4IInterpolator* algorithm,
 72                  const G4String& modelK,
 73                  const G4String& modelL,
 74                  const G4String& modelM,
 75                  G4double minE,
 76                  G4double maxE,
 77                  G4int bins,
 78                  G4double unitE,
 79                  G4double unitData,
 80                  G4int minZ, 
 81                  G4int maxZ)
 82   : interpolation(algorithm), eMin(minE), eMax(maxE), nBins(bins),
 83     unit1(unitE), unit2(unitData), zMin(minZ), zMax(maxZ)
 84 {
 85   crossSections = 0;
 86 
 87   crossModel.push_back(modelK);
 88   crossModel.push_back(modelL);
 89   crossModel.push_back(modelM);
 90   
 91   //std::cout << "PixeCrossSectionHandler constructor - crossModel[0] = " 
 92   //    << crossModel[0]
 93   //    << std::endl;
 94 
 95   ActiveElements();
 96 }
 97 
 98 G4PixeCrossSectionHandler::~G4PixeCrossSectionHandler()
 99 {
100   delete interpolation;
101   interpolation = 0;
102   std::map<G4int,G4IDataSet*,std::less<G4int> >::iterator pos;
103 
104   for (pos = dataMap.begin(); pos != dataMap.end(); ++pos)
105     {
106       // The following is a workaround for STL ObjectSpace implementation, 
107       // which does not support the standard and does not accept 
108       // the syntax pos->second
109       // G4IDataSet* dataSet = pos->second;
110       G4IDataSet* dataSet = (*pos).second;
111       delete dataSet;
112     }
113 
114   if (crossSections != 0)
115     {
116       std::size_t n = crossSections->size();
117       for (std::size_t i=0; i<n; ++i)
118   {
119     delete (*crossSections)[i];
120   }
121       delete crossSections;
122       crossSections = 0;
123     }
124 }
125 
126 void G4PixeCrossSectionHandler::Initialise(G4IInterpolator* algorithm,
127              const G4String& modelK,
128              const G4String& modelL,
129              const G4String& modelM,
130              G4double minE, G4double maxE, 
131              G4int numberOfBins,
132              G4double unitE, G4double unitData,
133              G4int minZ, G4int maxZ)
134 {
135   if (algorithm != 0) 
136     {
137       delete interpolation;
138       interpolation = algorithm;
139     }
140   else
141     {
142       interpolation = CreateInterpolation();
143     }
144 
145   eMin = minE;
146   eMax = maxE;
147   nBins = numberOfBins;
148   unit1 = unitE;
149   unit2 = unitData;
150   zMin = minZ;
151   zMax = maxZ;
152 
153   crossModel.push_back(modelK);
154   crossModel.push_back(modelL);
155   crossModel.push_back(modelM);
156 
157 }
158 
159 void G4PixeCrossSectionHandler::PrintData() const
160 {
161   std::map<G4int,G4IDataSet*,std::less<G4int> >::const_iterator pos;
162 
163   for (pos = dataMap.begin(); pos != dataMap.end(); pos++)
164     {
165       // The following is a workaround for STL ObjectSpace implementation, 
166       // which does not support the standard and does not accept 
167       // the syntax pos->first or pos->second
168       // G4int z = pos->first;
169       // G4IDataSet* dataSet = pos->second;
170       G4int z = (*pos).first;
171       G4IDataSet* dataSet = (*pos).second;     
172       G4cout << "---- Data set for Z = "
173        << z
174        << G4endl;
175       dataSet->PrintData();
176       G4cout << "--------------------------------------------------" << G4endl;
177     }
178 }
179 
180 void G4PixeCrossSectionHandler::LoadShellData(const G4String& fileName)
181 {
182   std::size_t nZ = activeZ.size();
183   for (std::size_t i=0; i<nZ; ++i)
184     {
185       G4int Z = (G4int) activeZ[i];     
186       G4IInterpolator* algo = interpolation->Clone();
187       G4IDataSet* dataSet = new G4PixeShellDataSet(Z, algo,crossModel[0],crossModel[1],crossModel[2]);
188 
189       // Degug printing
190       //std::cout << "PixeCrossSectionHandler::Load - "
191       //  << Z
192       //  << ", modelK = "
193       //  << crossModel[0]
194       //  << " fileName = "
195       //  << fileName
196       //  << std::endl;
197 
198       dataSet->LoadData(fileName);
199       dataMap[Z] = dataSet;
200     }
201 
202   // Build cross sections for materials if not already built
203   if (! crossSections)
204     {
205       BuildForMaterials();
206     }
207 
208 }
209 
210 void G4PixeCrossSectionHandler::Clear()
211 {
212   // Reset the map of data sets: remove the data sets from the map 
213   std::map<G4int,G4IDataSet*,std::less<G4int> >::iterator pos;
214 
215   if(! dataMap.empty())
216     {
217       for (pos = dataMap.begin(); pos != dataMap.end(); ++pos)
218   {
219     // The following is a workaround for STL ObjectSpace implementation, 
220     // which does not support the standard and does not accept
221     // the syntax pos->first or pos->second
222     // G4IDataSet* dataSet = pos->second;
223     G4IDataSet* dataSet = (*pos).second;
224     delete dataSet;
225     dataSet = 0;
226     G4int i = (*pos).first;
227     dataMap[i] = 0;
228   }
229       dataMap.clear();
230     }
231 
232   activeZ.clear();
233   ActiveElements();
234 }
235 
236 G4double G4PixeCrossSectionHandler::FindValue(G4int Z, G4double energy) const
237 {
238   G4double value = 0.;
239   
240   std::map<G4int,G4IDataSet*,std::less<G4int> >::const_iterator pos;
241   pos = dataMap.find(Z);
242   if (pos!= dataMap.end())
243     {
244       // The following is a workaround for STL ObjectSpace implementation, 
245       // which does not support the standard and does not accept 
246       // the syntax pos->first or pos->second
247       // G4IDataSet* dataSet = pos->second;
248       G4IDataSet* dataSet = (*pos).second;
249       value = dataSet->FindValue(energy);
250     }
251   else
252     {
253       G4cout << "WARNING: G4PixeCrossSectionHandler::FindValue(Z,e) did not find Z = "
254        << Z << G4endl;
255     }
256   return value;
257 }
258 
259 G4double G4PixeCrossSectionHandler::FindValue(G4int Z, G4double energy, 
260                 G4int shellIndex) const
261 {
262   G4double value = 0.;
263 
264   std::map<G4int,G4IDataSet*,std::less<G4int> >::const_iterator pos;
265   pos = dataMap.find(Z);
266   if (pos!= dataMap.end())
267     {
268       // The following is a workaround for STL ObjectSpace implementation, 
269       // which does not support the standard and does not accept 
270       // the syntax pos->first or pos->second
271       // G4IDataSet* dataSet = pos->second;
272       G4IDataSet* dataSet = (*pos).second;
273       if (shellIndex >= 0) 
274   {
275     G4int nComponents = (G4int)dataSet->NumberOfComponents();
276     if(shellIndex < nComponents)    
277       // The value is the cross section for shell component at given energy
278       value = dataSet->GetComponent(shellIndex)->FindValue(energy);
279     else 
280       {
281         G4cout << "WARNING: G4PixeCrossSectionHandler::FindValue(Z,e,shell) did not find"
282          << " shellIndex= " << shellIndex
283          << " for  Z= "
284          << Z << G4endl;
285       }
286   } else {
287   value = dataSet->FindValue(energy);
288       }
289     }
290   else
291     {
292       G4cout << "WARNING: G4PixeCrossSectionHandler::FindValue did not find Z = "
293        << Z << G4endl;
294     }
295   return value;
296 }
297 
298 
299 G4double G4PixeCrossSectionHandler::ValueForMaterial(const G4Material* material,
300                  G4double energy) const
301 {
302   G4double value = 0.;
303 
304   const G4ElementVector* elementVector = material->GetElementVector();
305   const G4double* nAtomsPerVolume = material->GetVecNbOfAtomsPerVolume();
306   std::size_t nElements = material->GetNumberOfElements();
307 
308   for (std::size_t i=0 ; i<nElements ; ++i)
309     {
310       G4int Z = (G4int) (*elementVector)[i]->GetZ();
311       G4double elementValue = FindValue(Z,energy);
312       G4double nAtomsVol = nAtomsPerVolume[i];
313       value += nAtomsVol * elementValue;
314     }
315 
316   return value;
317 }
318 
319 /*
320   G4IDataSet* G4PixeCrossSectionHandler::BuildMeanFreePathForMaterials(const G4DataVector*  energyCuts )
321   {
322   // Builds a CompositeDataSet containing the mean free path for each material
323   // in the material table
324 
325   G4DataVector energyVector;
326   G4double dBin = std::log10(eMax/eMin) / nBins;
327 
328   for (G4int i=0; i<nBins+1; i++)
329   {
330   energyVector.push_back(std::pow(10., std::log10(eMin)+i*dBin));
331   }
332 
333   // Factory method to build cross sections in derived classes,
334   // related to the type of physics process
335 
336   if (crossSections != 0)
337   {  // Reset the list of cross sections
338   std::vector<G4IDataSet*>::iterator mat;
339   if (! crossSections->empty())
340   {
341   for (mat = crossSections->begin(); mat!= crossSections->end(); ++mat)
342   {
343   G4IDataSet* set = *mat;
344   delete set;
345   set = 0;
346   }
347   crossSections->clear();
348   delete crossSections;
349   crossSections = 0;
350   }
351   }
352 
353   crossSections = BuildCrossSectionsForMaterials(energyVector);
354 
355   if (crossSections == 0)
356   G4Exception("G4PixeCrossSectionHandler::BuildMeanFreePathForMaterials", 
357   "pii00000201",
358   FatalException,
359   "crossSections = 0");
360   
361   G4IInterpolator* algo = CreateInterpolation();
362   G4IDataSet* materialSet = new G4CompositeDataSet(algo);
363 
364   G4DataVector* energies;
365   G4DataVector* data;
366   
367   const G4ProductionCutsTable* theCoupleTable=
368   G4ProductionCutsTable::GetProductionCutsTable();
369   std::size_t numOfCouples = theCoupleTable->GetTableSize();
370 
371 
372   for (std::size_t m=0; m<numOfCouples; m++)
373   {
374   energies = new G4DataVector;
375   data = new G4DataVector;
376   for (G4int bin=0; bin<nBins; bin++)
377   {
378   G4double energy = energyVector[bin];
379   energies->push_back(energy);
380   G4IDataSet* matCrossSet = (*crossSections)[m];
381   G4double materialCrossSection = 0.0;
382   G4int nElm = matCrossSet->NumberOfComponents();
383   for(G4int j=0; j<nElm; j++) {
384   materialCrossSection += matCrossSet->GetComponent(j)->FindValue(energy);
385   }
386 
387   if (materialCrossSection > 0.)
388   {
389   data->push_back(1./materialCrossSection);
390   }
391   else
392   {
393   data->push_back(DBL_MAX);
394   }
395   }
396   G4IInterpolator* algo = CreateInterpolation();
397   G4IDataSet* dataSet = new G4DataSet(m,energies,data,algo,1.,1.);
398   materialSet->AddComponent(dataSet);
399   }
400 
401   return materialSet;
402   }
403 
404 */
405 
406 void G4PixeCrossSectionHandler::BuildForMaterials()
407 {
408   // Builds a CompositeDataSet containing the mean free path for each material
409   // in the material table
410 
411   G4DataVector energyVector;
412   G4double dBin = std::log10(eMax/eMin) / nBins;
413 
414   for (G4int i=0; i<nBins+1; i++)
415     {
416       energyVector.push_back(std::pow(10., std::log10(eMin)+i*dBin));
417     }
418 
419   if (crossSections != 0)
420     {  // Reset the list of cross sections
421       std::vector<G4IDataSet*>::iterator mat;
422       if (! crossSections->empty())
423   {
424     for (mat = crossSections->begin(); mat!= crossSections->end(); ++mat)
425       {
426         G4IDataSet* set = *mat;
427         delete set;
428         set = 0;
429       }
430     crossSections->clear();
431     delete crossSections;
432     crossSections = 0;
433   }
434     }
435 
436   crossSections = BuildCrossSectionsForMaterials(energyVector);
437 
438   if (crossSections == 0)
439     G4Exception("G4PixeCrossSectionHandler::BuildForMaterials",
440     "pii00000210",
441     FatalException,
442     ", crossSections = 0");
443 
444   return;
445 }
446 
447 
448 G4int G4PixeCrossSectionHandler::SelectRandomAtom(const G4Material* material,
449               G4double e) const
450 {
451   // Select randomly an element within the material, according to the weight
452   // determined by the cross sections in the data set
453 
454   G4int nElements = (G4int)material->GetNumberOfElements();
455 
456   // Special case: the material consists of one element
457   if (nElements == 1)
458     {
459       G4int Z = (G4int) material->GetZ();
460       return Z;
461     }
462 
463   // Composite material
464 
465   const G4ElementVector* elementVector = material->GetElementVector();
466   std::size_t materialIndex = material->GetIndex();
467 
468   G4IDataSet* materialSet = (*crossSections)[materialIndex];
469   G4double materialCrossSection0 = 0.0;
470   G4DataVector cross;
471   cross.clear();
472   for ( G4int i=0; i < nElements; ++i )
473     {
474       G4double cr = materialSet->GetComponent(i)->FindValue(e);
475       materialCrossSection0 += cr;
476       cross.push_back(materialCrossSection0);
477     }
478 
479   G4double random = G4UniformRand() * materialCrossSection0;
480 
481   for (G4int k=0 ; k < nElements ; ++k )
482     {
483       if (random <= cross[k]) return (G4int) (*elementVector)[k]->GetZ();
484     }
485   // It should never get here
486   return 0;
487 }
488 
489 /*
490   const G4Element* G4PixeCrossSectionHandler::SelectRandomElement(const G4MaterialCutsCouple* couple,
491   G4double e) const
492   {
493   // Select randomly an element within the material, according to the weight determined
494   // by the cross sections in the data set
495 
496   const G4Material* material = couple->GetMaterial();
497   G4Element* nullElement = 0;
498   G4int nElements = material->GetNumberOfElements();
499   const G4ElementVector* elementVector = material->GetElementVector();
500 
501   // Special case: the material consists of one element
502   if (nElements == 1)
503   {
504   G4Element* element = (*elementVector)[0];
505   return element;
506   }
507   else
508   {
509   // Composite material
510 
511   std::size_t materialIndex = couple->GetIndex();
512 
513   G4IDataSet* materialSet = (*crossSections)[materialIndex];
514   G4double materialCrossSection0 = 0.0;
515   G4DataVector cross;
516   cross.clear();
517   for (G4int i=0; i<nElements; i++)
518   {
519   G4double cr = materialSet->GetComponent(i)->FindValue(e);
520   materialCrossSection0 += cr;
521   cross.push_back(materialCrossSection0);
522   }
523 
524   G4double random = G4UniformRand() * materialCrossSection0;
525 
526   for (G4int k=0 ; k < nElements ; k++ )
527   {
528   if (random <= cross[k]) return (*elementVector)[k];
529   }
530   // It should never end up here
531   G4cout << "G4PixeCrossSectionHandler::SelectRandomElement - no element found" << G4endl;
532   return nullElement;
533   }
534   }
535 */
536 
537 
538 G4int G4PixeCrossSectionHandler::SelectRandomShell(G4int Z, G4double e) const
539 {
540   // Select randomly a shell, according to the weight determined by the cross sections
541   // in the data set
542 
543   // Note for later improvement: it would be useful to add a cache mechanism for already
544   // used shells to improve performance
545 
546   G4int shell = 0;
547 
548   G4double totCrossSection = FindValue(Z,e);
549   G4double random = G4UniformRand() * totCrossSection;
550   G4double partialSum = 0.;
551 
552   G4IDataSet* dataSet = nullptr;
553   auto pos = dataMap.find(Z);
554   // The following is a workaround for STL ObjectSpace implementation,
555   // which does not support the standard and does not accept
556   // the syntax pos->first or pos->second
557   // if (pos != dataMap.end()) dataSet = pos->second;
558   if (pos != dataMap.end()) dataSet = (*pos).second;
559 
560   if (dataSet != nullptr) {
561     G4int nShells = (G4int)dataSet->NumberOfComponents();
562     for (G4int i=0; i<nShells; ++i)
563     {
564       const G4IDataSet* shellDataSet = dataSet->GetComponent(i);
565       if (shellDataSet != 0)
566       {
567         G4double value = shellDataSet->FindValue(e);
568         partialSum += value;
569         if (random <= partialSum) return i;
570       }
571     }
572   }
573   // It should never get here
574   return shell;
575 }
576 
577 void G4PixeCrossSectionHandler::ActiveElements()
578 {
579   const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
580   if (materialTable == 0)
581     G4Exception("G4PixeCrossSectionHandler::ActiveElements",
582           "pii00000220",
583           FatalException,
584           "no MaterialTable found");
585 
586   std::size_t nMaterials = G4Material::GetNumberOfMaterials();
587 
588   for (std::size_t mat=0; mat<nMaterials; ++mat)
589     {
590       const G4Material* material= (*materialTable)[mat];
591       const G4ElementVector* elementVector = material->GetElementVector();
592       const std::size_t nElements = material->GetNumberOfElements();
593 
594       for (std::size_t iEl=0; iEl<nElements; ++iEl)
595   {
596     G4double Z = (*elementVector)[iEl]->GetZ();
597     if (!(activeZ.contains(Z)) && Z >= zMin && Z <= zMax)
598       {
599         activeZ.push_back(Z);
600       }
601   }
602     }
603 }
604 
605 G4IInterpolator* G4PixeCrossSectionHandler::CreateInterpolation()
606 {
607   G4IInterpolator* algorithm = new G4LogLogInterpolator;
608   return algorithm;
609 }
610 
611 G4int G4PixeCrossSectionHandler::NumberOfComponents(G4int Z) const
612 {
613   G4int n = 0;
614 
615   std::map<G4int,G4IDataSet*,std::less<G4int> >::const_iterator pos;
616   pos = dataMap.find(Z);
617   if (pos!= dataMap.end())
618     {
619       G4IDataSet* dataSet = (*pos).second;
620       n = (G4int)dataSet->NumberOfComponents();
621     }
622   else
623     {
624       G4cout << "WARNING: G4PixeCrossSectionHandler::NumberOfComponents did not "
625              << "find Z = "
626              << Z << G4endl;
627     }
628   return n;
629 }
630 
631 
632 std::vector<G4IDataSet*>*
633 G4PixeCrossSectionHandler::BuildCrossSectionsForMaterials(const G4DataVector& energyVector)
634 {
635   G4DataVector* energies;
636   G4DataVector* data;
637 
638   std::vector<G4IDataSet*>* matCrossSections = new std::vector<G4IDataSet*>;
639 
640   //const G4ProductionCutsTable* theCoupleTable=G4ProductionCutsTable::GetProductionCutsTable();
641   //std::size_t numOfCouples = theCoupleTable->GetTableSize();
642 
643   std::size_t nOfBins = energyVector.size();
644   const auto interpolationAlgo = std::unique_ptr<G4IInterpolator>(CreateInterpolation());
645 
646   const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
647   if (materialTable == 0)
648     G4Exception("G4PixeCrossSectionHandler::BuildCrossSectionsForMaterials",
649     "pii00000230",
650     FatalException,
651     "no MaterialTable found");
652 
653   std::size_t nMaterials = G4Material::GetNumberOfMaterials();
654 
655   for (std::size_t mat=0; mat<nMaterials; ++mat)
656     {
657       const G4Material* material = (*materialTable)[mat];
658       G4int nElements = (G4int)material->GetNumberOfElements();
659       const G4ElementVector* elementVector = material->GetElementVector();
660       const G4double* nAtomsPerVolume = material->GetAtomicNumDensityVector();
661 
662       G4IInterpolator* algo = interpolationAlgo->Clone();
663 
664       G4IDataSet* setForMat = new G4CompositeDataSet(algo,1.,1.);
665 
666       for (G4int i=0; i<nElements; ++i) {
667  
668         G4int Z = (G4int) (*elementVector)[i]->GetZ();
669         G4double density = nAtomsPerVolume[i];
670 
671         energies = new G4DataVector;
672         data = new G4DataVector;
673 
674 
675         for (std::size_t bin=0; bin<nOfBins; ++bin)
676     {
677       G4double e = energyVector[bin];
678       energies->push_back(e);
679             G4double cross = 0.;
680       if (Z >= zMin && Z <= zMax) cross = density*FindValue(Z,e);
681       data->push_back(cross);
682     }
683 
684         G4IInterpolator* algo1 = interpolationAlgo->Clone();
685         G4IDataSet* elSet = new G4DataSet(i,energies,data,algo1,1.,1.);
686         setForMat->AddComponent(elSet);
687       }
688 
689       matCrossSections->push_back(setForMat);
690     }
691   return matCrossSections;
692 }
693 
694 
695 G4double G4PixeCrossSectionHandler::MicroscopicCrossSection(const G4ParticleDefinition* particleDef,
696                   G4double kineticEnergy,
697                   G4double Z,
698                   G4double deltaCut) const
699 {
700   // Cross section formula is OK for spin=0, 1/2, 1 only !
701   // Calculates the microscopic cross section in Geant4 internal units
702   // Formula documented in Geant4 Phys. Ref. Manual
703   // ( it is called for elements, AtomicNumber = z )
704 
705     G4double cross = 0.;
706 
707   // Particle mass and energy
708     G4double particleMass = particleDef->GetPDGMass();
709     G4double energy = kineticEnergy + particleMass;
710 
711   // Some kinematics
712   G4double gamma = energy / particleMass;
713   G4double beta2 = 1. - 1. / (gamma * gamma);
714   G4double var = electron_mass_c2 / particleMass;
715   G4double tMax = 2. * electron_mass_c2 * (gamma*gamma - 1.) / (1. + 2.*gamma*var + var*var);
716 
717   // Calculate the total cross section
718 
719   if ( tMax > deltaCut ) 
720     {
721       var = deltaCut / tMax;
722       cross = (1. - var * (1. - beta2 * std::log(var))) / deltaCut;
723       
724       G4double spin = particleDef->GetPDGSpin() ;
725       
726       // +term for spin=1/2 particle
727       if (spin == 0.5) 
728   {
729     cross +=  0.5 * (tMax - deltaCut) / (energy*energy);
730   }
731       // +term for spin=1 particle
732       else if (spin > 0.9 )
733   {
734     cross += -std::log(var) / (3.*deltaCut) + (tMax-deltaCut) * 
735       ((5.+1./var)*0.25 /(energy*energy) - beta2 / (tMax*deltaCut))/3.;
736   }
737       cross *= twopi_mc2_rcl2 * Z / beta2 ;
738     }
739 
740   //std::cout << "Microscopic = " << cross/barn 
741   //    << ", e = " << kineticEnergy/MeV <<std:: endl; 
742 
743   return cross;
744 }
745 
746