source: trunk/source/processes/electromagnetic/pii/src/G4PixeCrossSectionHandler.cc@ 1357

Last change on this file since 1357 was 1350, checked in by garnier, 15 years ago

update to last version 4.9.4
File size: 21.0 KB
RevLine 
[1350]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// $Id: G4PixeCrossSectionHandler.cc,v 1.2 2010/11/19 17:16:21 pia Exp $
28// GEANT4 tag $Name: geant4-09-04-ref-00 $
29//
30// Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
31//
32// History:
33// -----------
34// 16 Jun 2008 MGP Created; Cross section manager for hadron impact ionization
35// Documented in:
36// M.G. Pia et al., PIXE Simulation With Geant4,
37// IEEE Trans. Nucl. Sci., vol. 56, no. 6, pp. 3614-3649, Dec. 2009
38//
39// -------------------------------------------------------------------
40
41#include "G4PixeCrossSectionHandler.hh"
42#include "G4IInterpolator.hh"
43#include "G4LogLogInterpolator.hh"
44#include "G4IDataSet.hh"
45#include "G4DataSet.hh"
46#include "G4CompositeDataSet.hh"
47#include "G4PixeShellDataSet.hh"
48#include "G4ProductionCutsTable.hh"
49#include "G4Material.hh"
50#include "G4Element.hh"
51#include "Randomize.hh"
52#include "G4SystemOfUnits.hh"
53#include "G4ParticleDefinition.hh"
54
55#include <map>
56#include <vector>
57#include <fstream>
58#include <sstream>
59
60
61G4PixeCrossSectionHandler::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
71G4PixeCrossSectionHandler::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
98G4PixeCrossSectionHandler::~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 size_t n = crossSections->size();
117 for (size_t i=0; i<n; i++)
118 {
119 delete (*crossSections)[i];
120 }
121 delete crossSections;
122 crossSections = 0;
123 }
124}
125
126void 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
159void 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
180void G4PixeCrossSectionHandler::LoadShellData(const G4String& fileName)
181{
182 size_t nZ = activeZ.size();
183 for (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
210void 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
236G4double 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
259G4double 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 = 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
299G4double 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 G4int nElements = material->GetNumberOfElements();
307
308 for (G4int 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, crossSections = 0");
357
358 G4IInterpolator* algo = CreateInterpolation();
359 G4IDataSet* materialSet = new G4CompositeDataSet(algo);
360
361 G4DataVector* energies;
362 G4DataVector* data;
363
364 const G4ProductionCutsTable* theCoupleTable=
365 G4ProductionCutsTable::GetProductionCutsTable();
366 size_t numOfCouples = theCoupleTable->GetTableSize();
367
368
369 for (size_t m=0; m<numOfCouples; m++)
370 {
371 energies = new G4DataVector;
372 data = new G4DataVector;
373 for (G4int bin=0; bin<nBins; bin++)
374 {
375 G4double energy = energyVector[bin];
376 energies->push_back(energy);
377 G4IDataSet* matCrossSet = (*crossSections)[m];
378 G4double materialCrossSection = 0.0;
379 G4int nElm = matCrossSet->NumberOfComponents();
380 for(G4int j=0; j<nElm; j++) {
381 materialCrossSection += matCrossSet->GetComponent(j)->FindValue(energy);
382 }
383
384 if (materialCrossSection > 0.)
385 {
386 data->push_back(1./materialCrossSection);
387 }
388 else
389 {
390 data->push_back(DBL_MAX);
391 }
392 }
393 G4IInterpolator* algo = CreateInterpolation();
394 G4IDataSet* dataSet = new G4DataSet(m,energies,data,algo,1.,1.);
395 materialSet->AddComponent(dataSet);
396 }
397
398 return materialSet;
399 }
400
401*/
402
403void G4PixeCrossSectionHandler::BuildForMaterials()
404{
405 // Builds a CompositeDataSet containing the mean free path for each material
406 // in the material table
407
408 G4DataVector energyVector;
409 G4double dBin = std::log10(eMax/eMin) / nBins;
410
411 for (G4int i=0; i<nBins+1; i++)
412 {
413 energyVector.push_back(std::pow(10., std::log10(eMin)+i*dBin));
414 }
415
416 if (crossSections != 0)
417 { // Reset the list of cross sections
418 std::vector<G4IDataSet*>::iterator mat;
419 if (! crossSections->empty())
420 {
421 for (mat = crossSections->begin(); mat!= crossSections->end(); ++mat)
422 {
423 G4IDataSet* set = *mat;
424 delete set;
425 set = 0;
426 }
427 crossSections->clear();
428 delete crossSections;
429 crossSections = 0;
430 }
431 }
432
433 crossSections = BuildCrossSectionsForMaterials(energyVector);
434
435 if (crossSections == 0)
436 G4Exception("G4PixeCrossSectionHandler::BuildForMaterials, crossSections = 0");
437
438 return;
439}
440
441
442G4int G4PixeCrossSectionHandler::SelectRandomAtom(const G4Material* material,
443 G4double e) const
444{
445 // Select randomly an element within the material, according to the weight
446 // determined by the cross sections in the data set
447
448 G4int nElements = material->GetNumberOfElements();
449
450 // Special case: the material consists of one element
451 if (nElements == 1)
452 {
453 G4int Z = (G4int) material->GetZ();
454 return Z;
455 }
456
457 // Composite material
458
459 const G4ElementVector* elementVector = material->GetElementVector();
460 size_t materialIndex = material->GetIndex();
461
462 G4IDataSet* materialSet = (*crossSections)[materialIndex];
463 G4double materialCrossSection0 = 0.0;
464 G4DataVector cross;
465 cross.clear();
466 for ( G4int i=0; i < nElements; i++ )
467 {
468 G4double cr = materialSet->GetComponent(i)->FindValue(e);
469 materialCrossSection0 += cr;
470 cross.push_back(materialCrossSection0);
471 }
472
473 G4double random = G4UniformRand() * materialCrossSection0;
474
475 for (G4int k=0 ; k < nElements ; k++ )
476 {
477 if (random <= cross[k]) return (G4int) (*elementVector)[k]->GetZ();
478 }
479 // It should never get here
480 return 0;
481}
482
483/*
484 const G4Element* G4PixeCrossSectionHandler::SelectRandomElement(const G4MaterialCutsCouple* couple,
485 G4double e) const
486 {
487 // Select randomly an element within the material, according to the weight determined
488 // by the cross sections in the data set
489
490 const G4Material* material = couple->GetMaterial();
491 G4Element* nullElement = 0;
492 G4int nElements = material->GetNumberOfElements();
493 const G4ElementVector* elementVector = material->GetElementVector();
494
495 // Special case: the material consists of one element
496 if (nElements == 1)
497 {
498 G4Element* element = (*elementVector)[0];
499 return element;
500 }
501 else
502 {
503 // Composite material
504
505 size_t materialIndex = couple->GetIndex();
506
507 G4IDataSet* materialSet = (*crossSections)[materialIndex];
508 G4double materialCrossSection0 = 0.0;
509 G4DataVector cross;
510 cross.clear();
511 for (G4int i=0; i<nElements; i++)
512 {
513 G4double cr = materialSet->GetComponent(i)->FindValue(e);
514 materialCrossSection0 += cr;
515 cross.push_back(materialCrossSection0);
516 }
517
518 G4double random = G4UniformRand() * materialCrossSection0;
519
520 for (G4int k=0 ; k < nElements ; k++ )
521 {
522 if (random <= cross[k]) return (*elementVector)[k];
523 }
524 // It should never end up here
525 G4cout << "G4PixeCrossSectionHandler::SelectRandomElement - no element found" << G4endl;
526 return nullElement;
527 }
528 }
529*/
530
531
532G4int G4PixeCrossSectionHandler::SelectRandomShell(G4int Z, G4double e) const
533{
534 // Select randomly a shell, according to the weight determined by the cross sections
535 // in the data set
536
537 // Note for later improvement: it would be useful to add a cache mechanism for already
538 // used shells to improve performance
539
540 G4int shell = 0;
541
542 G4double totCrossSection = FindValue(Z,e);
543 G4double random = G4UniformRand() * totCrossSection;
544 G4double partialSum = 0.;
545
546 G4IDataSet* dataSet = 0;
547 std::map<G4int,G4IDataSet*,std::less<G4int> >::const_iterator pos;
548 pos = dataMap.find(Z);
549 // The following is a workaround for STL ObjectSpace implementation,
550 // which does not support the standard and does not accept
551 // the syntax pos->first or pos->second
552 // if (pos != dataMap.end()) dataSet = pos->second;
553 if (pos != dataMap.end()) dataSet = (*pos).second;
554
555 size_t nShells = dataSet->NumberOfComponents();
556 for (size_t i=0; i<nShells; i++)
557 {
558 const G4IDataSet* shellDataSet = dataSet->GetComponent(i);
559 if (shellDataSet != 0)
560 {
561 G4double value = shellDataSet->FindValue(e);
562 partialSum += value;
563 if (random <= partialSum) return i;
564 }
565 }
566 // It should never get here
567 return shell;
568}
569
570void G4PixeCrossSectionHandler::ActiveElements()
571{
572 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
573 if (materialTable == 0)
574 G4Exception("G4PixeCrossSectionHandler::ActiveElements - no MaterialTable found)");
575
576 G4int nMaterials = G4Material::GetNumberOfMaterials();
577
578 for (G4int m=0; m<nMaterials; m++)
579 {
580 const G4Material* material= (*materialTable)[m];
581 const G4ElementVector* elementVector = material->GetElementVector();
582 const G4int nElements = material->GetNumberOfElements();
583
584 for (G4int iEl=0; iEl<nElements; iEl++)
585 {
586 G4Element* element = (*elementVector)[iEl];
587 G4double Z = element->GetZ();
588 if (!(activeZ.contains(Z)) && Z >= zMin && Z <= zMax)
589 {
590 activeZ.push_back(Z);
591 }
592 }
593 }
594}
595
596G4IInterpolator* G4PixeCrossSectionHandler::CreateInterpolation()
597{
598 G4IInterpolator* algorithm = new G4LogLogInterpolator;
599 return algorithm;
600}
601
602G4int G4PixeCrossSectionHandler::NumberOfComponents(G4int Z) const
603{
604 G4int n = 0;
605
606 std::map<G4int,G4IDataSet*,std::less<G4int> >::const_iterator pos;
607 pos = dataMap.find(Z);
608 if (pos!= dataMap.end())
609 {
610 G4IDataSet* dataSet = (*pos).second;
611 n = dataSet->NumberOfComponents();
612 }
613 else
614 {
615 G4cout << "WARNING: G4PixeCrossSectionHandler::NumberOfComponents did not "
616 << "find Z = "
617 << Z << G4endl;
618 }
619 return n;
620}
621
622
623std::vector<G4IDataSet*>*
624G4PixeCrossSectionHandler::BuildCrossSectionsForMaterials(const G4DataVector& energyVector)
625{
626 G4DataVector* energies;
627 G4DataVector* data;
628
629 std::vector<G4IDataSet*>* matCrossSections = new std::vector<G4IDataSet*>;
630
631 //const G4ProductionCutsTable* theCoupleTable=G4ProductionCutsTable::GetProductionCutsTable();
632 //size_t numOfCouples = theCoupleTable->GetTableSize();
633
634 size_t nOfBins = energyVector.size();
635 const G4IInterpolator* interpolationAlgo = CreateInterpolation();
636
637 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
638 if (materialTable == 0)
639 G4Exception("G4PixeCrossSectionHandler - no MaterialTable found)");
640
641 G4int nMaterials = G4Material::GetNumberOfMaterials();
642
643 for (G4int m=0; m<nMaterials; m++)
644 {
645 const G4Material* material = (*materialTable)[m];
646 G4int nElements = material->GetNumberOfElements();
647 const G4ElementVector* elementVector = material->GetElementVector();
648 const G4double* nAtomsPerVolume = material->GetAtomicNumDensityVector();
649
650 G4IInterpolator* algo = interpolationAlgo->Clone();
651
652 G4IDataSet* setForMat = new G4CompositeDataSet(algo,1.,1.);
653
654 for (G4int i=0; i<nElements; i++) {
655
656 G4int Z = (G4int) (*elementVector)[i]->GetZ();
657 G4double density = nAtomsPerVolume[i];
658
659 energies = new G4DataVector;
660 data = new G4DataVector;
661
662
663 for (size_t bin=0; bin<nOfBins; bin++)
664 {
665 G4double e = energyVector[bin];
666 energies->push_back(e);
667 G4double cross = 0.;
668 if (Z >= zMin && Z <= zMax) cross = density*FindValue(Z,e);
669 data->push_back(cross);
670 }
671
672 G4IInterpolator* algo1 = interpolationAlgo->Clone();
673 G4IDataSet* elSet = new G4DataSet(i,energies,data,algo1,1.,1.);
674 setForMat->AddComponent(elSet);
675 }
676
677 matCrossSections->push_back(setForMat);
678 }
679 return matCrossSections;
680}
681
682
683G4double G4PixeCrossSectionHandler::MicroscopicCrossSection(const G4ParticleDefinition* particleDef,
684 G4double kineticEnergy,
685 G4double Z,
686 G4double deltaCut) const
687{
688 // Cross section formula is OK for spin=0, 1/2, 1 only !
689 // Calculates the microscopic cross section in Geant4 internal units
690 // Formula documented in Geant4 Phys. Ref. Manual
691 // ( it is called for elements, AtomicNumber = z )
692
693 G4double cross = 0.;
694
695 // Particle mass and energy
696 G4double particleMass = particleDef->GetPDGMass();
697 G4double energy = kineticEnergy + particleMass;
698
699 // Some kinematics
700 G4double gamma = energy / particleMass;
701 G4double beta2 = 1. - 1. / (gamma * gamma);
702 G4double var = electron_mass_c2 / particleMass;
703 G4double tMax = 2. * electron_mass_c2 * (gamma*gamma - 1.) / (1. + 2.*gamma*var + var*var);
704
705 // Calculate the total cross section
706
707 if ( tMax > deltaCut )
708 {
709 var = deltaCut / tMax;
710 cross = (1. - var * (1. - beta2 * std::log(var))) / deltaCut;
711
712 G4double spin = particleDef->GetPDGSpin() ;
713
714 // +term for spin=1/2 particle
715 if (spin == 0.5)
716 {
717 cross += 0.5 * (tMax - deltaCut) / (energy*energy);
718 }
719 // +term for spin=1 particle
720 else if (spin > 0.9 )
721 {
722 cross += -std::log(var) / (3.*deltaCut) + (tMax-deltaCut) *
723 ((5.+1./var)*0.25 /(energy*energy) - beta2 / (tMax*deltaCut))/3.;
724 }
725 cross *= twopi_mc2_rcl2 * Z / beta2 ;
726 }
727
728 //std::cout << "Microscopic = " << cross/barn
729 // << ", e = " << kineticEnergy/MeV <<std:: endl;
730
731 return cross;
732}
733
Note: See TracBrowser for help on using the repository browser.