ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLNuclearDensity.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLNuclearDensity.cc
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 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France
33 //
34 #define INCLXX_IN_GEANT4_MODE 1
35 
36 #include "globals.hh"
37 
38 #include "G4INCLNuclearDensity.hh"
39 #include "G4INCLParticleTable.hh"
40 #include "G4INCLGlobals.hh"
41 #include <algorithm>
42 
43 namespace G4INCL {
44 
45  NuclearDensity::NuclearDensity(const G4int A, const G4int Z, const G4int S, InterpolationTable const * const rpCorrelationTableProton, InterpolationTable const * const rpCorrelationTableNeutron, InterpolationTable const * const rpCorrelationTableLambda) :
46  theA(A),
47  theZ(Z),
48  theS(S),
49  theMaximumRadius(std::min((*rpCorrelationTableProton)(1.), (*rpCorrelationTableNeutron)(1.))),
50  theProtonNuclearRadius(ParticleTable::getNuclearRadius(Proton,theA,theZ))
51  {
52  std::fill(rFromP, rFromP + UnknownParticle, static_cast<InterpolationTable*>(NULL));
53  rFromP[Proton] = rpCorrelationTableProton;
54  rFromP[Neutron] = rpCorrelationTableNeutron;
55  rFromP[Lambda] = rpCorrelationTableLambda;
56  rFromP[DeltaPlusPlus] = rpCorrelationTableProton;
57  rFromP[DeltaPlus] = rpCorrelationTableProton;
58  rFromP[DeltaZero] = rpCorrelationTableNeutron;
59  rFromP[DeltaMinus] = rpCorrelationTableNeutron;
60  // The interpolation table for local-energy look-ups is simply obtained by
61  // inverting the r-p correlation table.
62  std::fill(pFromR, pFromR + UnknownParticle, static_cast<InterpolationTable*>(NULL));
63  pFromR[Proton] = new InterpolationTable(rFromP[Proton]->getNodeValues(), rFromP[Proton]->getNodeAbscissae());
64  pFromR[Neutron] = new InterpolationTable(rFromP[Neutron]->getNodeValues(), rFromP[Neutron]->getNodeAbscissae());
65  pFromR[Lambda] = new InterpolationTable(rFromP[Lambda]->getNodeValues(), rFromP[Lambda]->getNodeAbscissae());
66  pFromR[DeltaPlusPlus] = new InterpolationTable(rFromP[DeltaPlusPlus]->getNodeValues(), rFromP[DeltaPlusPlus]->getNodeAbscissae());
67  pFromR[DeltaPlus] = new InterpolationTable(rFromP[DeltaPlus]->getNodeValues(), rFromP[DeltaPlus]->getNodeAbscissae());
68  pFromR[DeltaZero] = new InterpolationTable(rFromP[DeltaZero]->getNodeValues(), rFromP[DeltaZero]->getNodeAbscissae());
69  pFromR[DeltaMinus] = new InterpolationTable(rFromP[DeltaMinus]->getNodeValues(), rFromP[DeltaMinus]->getNodeAbscissae());
70  INCL_DEBUG("Interpolation table for proton local energy (A=" << theA << ", Z=" << theZ << ") initialised:"
71  << '\n'
72  << pFromR[Proton]->print()
73  << '\n'
74  << "Interpolation table for neutron local energy (A=" << theA << ", Z=" << theZ << ") initialised:"
75  << '\n'
76  << pFromR[Neutron]->print()
77  << '\n'
78  << "Interpolation table for lambda local energy (A=" << theA << ", Z=" << theZ << ", S=" << theS << ") initialised:"
79  << '\n'
80  << pFromR[Lambda]->print()
81  << '\n'
82  << "Interpolation table for delta++ local energy (A=" << theA << ", Z=" << theZ << ") initialised:"
83  << '\n'
85  << '\n'
86  << "Interpolation table for delta+ local energy (A=" << theA << ", Z=" << theZ << ") initialised:"
87  << '\n'
88  << pFromR[DeltaPlus]->print()
89  << '\n'
90  << "Interpolation table for delta0 local energy (A=" << theA << ", Z=" << theZ << ") initialised:"
91  << '\n'
92  << pFromR[DeltaZero]->print()
93  << '\n'
94  << "Interpolation table for delta- local energy (A=" << theA << ", Z=" << theZ << ") initialised:"
95  << '\n'
96  << pFromR[DeltaMinus]->print()
97  << '\n');
99  }
100 
102  // We don't delete the rFromP tables, which are cached in the
103  // NuclearDensityFactory
104  delete pFromR[Proton];
105  delete pFromR[Neutron];
106  delete pFromR[Lambda];
107  delete pFromR[DeltaPlusPlus];
108  delete pFromR[DeltaPlus];
109  delete pFromR[DeltaZero];
110  delete pFromR[DeltaMinus];
111  }
112 
114  theA(rhs.theA),
115  theZ(rhs.theZ),
116  theS(rhs.theS),
117  theMaximumRadius(rhs.theMaximumRadius),
118  theProtonNuclearRadius(rhs.theProtonNuclearRadius)
119  {
120  // rFromP is owned by NuclearDensityFactory, so shallow copy is sufficient
121  std::fill(rFromP, rFromP + UnknownParticle, static_cast<InterpolationTable*>(NULL));
122  rFromP[Proton] = rhs.rFromP[Proton];
123  rFromP[Neutron] = rhs.rFromP[Neutron];
124  rFromP[Lambda] = rhs.rFromP[Lambda];
129  // deep copy for pFromR
130  std::fill(pFromR, pFromR + UnknownParticle, static_cast<InterpolationTable*>(NULL));
131  pFromR[Proton] = new InterpolationTable(*(rhs.pFromR[Proton]));
133  pFromR[Lambda] = new InterpolationTable(*(rhs.pFromR[Lambda]));
139  }
140 
142  NuclearDensity temporaryDensity(rhs);
143  swap(temporaryDensity);
144  return *this;
145  }
146 
148  std::swap(theA, rhs.theA);
149  std::swap(theZ, rhs.theZ);
150  std::swap(theS, rhs.theS);
154  std::swap(rFromP[Proton], rhs.rFromP[Proton]);
155  std::swap(rFromP[Neutron], rhs.rFromP[Neutron]);
156  std::swap(rFromP[Lambda], rhs.rFromP[Lambda]);
157  std::swap(rFromP[DeltaPlusPlus], rhs.rFromP[DeltaPlusPlus]);
158  std::swap(rFromP[DeltaPlus], rhs.rFromP[DeltaPlus]);
159  std::swap(rFromP[DeltaZero], rhs.rFromP[DeltaZero]);
160  std::swap(rFromP[DeltaMinus], rhs.rFromP[DeltaMinus]);
161  std::swap(pFromR[Proton], rhs.pFromR[Proton]);
162  std::swap(pFromR[Neutron], rhs.pFromR[Neutron]);
163  std::swap(pFromR[DeltaPlusPlus], rhs.pFromR[DeltaPlusPlus]);
164  std::swap(pFromR[DeltaPlus], rhs.pFromR[DeltaPlus]);
165  std::swap(pFromR[DeltaZero], rhs.pFromR[DeltaZero]);
166  std::swap(pFromR[DeltaMinus], rhs.pFromR[DeltaMinus]);
167  }
168 
170  const G4double theProtonRadius = 0.88; // fm
171  const G4double theProtonTransmissionRadius = theProtonNuclearRadius + theProtonRadius;
172 
173  transmissionRadius[Proton] = theProtonTransmissionRadius;
176  transmissionRadius[DeltaPlusPlus] = theProtonTransmissionRadius;
177  transmissionRadius[DeltaPlus] = theProtonTransmissionRadius;
178  transmissionRadius[DeltaMinus] = theProtonTransmissionRadius;
180  transmissionRadius[SigmaPlus] = theProtonTransmissionRadius;
181  transmissionRadius[SigmaMinus] = theProtonTransmissionRadius;
184 
185  // transmission radii for neutral particles intentionally left uninitialised
186  }
187 
189 // assert(t==Proton || t==Neutron || t==Lambda || t==DeltaPlusPlus || t==DeltaPlus || t==DeltaZero || t==DeltaMinus);
190  return (*(rFromP[t]))(p);
191  }
192 
194 // assert(t==Proton || t==Neutron || t==Lambda || t==DeltaPlusPlus || t==DeltaPlus || t==DeltaZero || t==DeltaMinus);
195  return (*(pFromR[t]))(r);
196  }
197 
198 }