Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4INCLParticle.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 /*
39  * Particle.cc
40  *
41  * \date Jun 5, 2009
42  * \author Pekka Kaitaniemi
43  */
44 
45 #include "G4INCLParticle.hh"
46 #include "G4INCLParticleTable.hh"
47 
48 namespace G4INCL {
49 
50 #ifdef INCLXX_IN_GEANT4_MODE
51  std::vector<G4double> Particle::INCLBiasVector;
52 #else
53  G4ThreadLocal std::vector<G4double> Particle::INCLBiasVector;
54  //G4VectorCache<G4double> Particle::INCLBiasVector;
55 #endif
58 
60  : theZ(0), theA(0), theS(0),
61  theParticipantType(TargetSpectator),
62  theType(UnknownParticle),
63  theEnergy(0.0),
64  thePropagationEnergy(&theEnergy),
65  theFrozenEnergy(theEnergy),
66  theMomentum(ThreeVector(0.,0.,0.)),
67  thePropagationMomentum(&theMomentum),
68  theFrozenMomentum(theMomentum),
69  thePosition(ThreeVector(0.,0.,0.)),
70  nCollisions(0),
71  nDecays(0),
72  thePotentialEnergy(0.0),
73  rpCorrelated(false),
74  uncorrelatedMomentum(0.),
75  theParticleBias(1.),
76  theHelicity(0.0),
77  emissionTime(0.0),
78  outOfWell(false),
79  theMass(0.)
80  {
81  ID = nextID;
82  nextID++;
83  }
84 
86  ThreeVector const &momentum, ThreeVector const &position)
87  : theEnergy(energy),
88  thePropagationEnergy(&theEnergy),
89  theFrozenEnergy(theEnergy),
90  theMomentum(momentum),
91  thePropagationMomentum(&theMomentum),
92  theFrozenMomentum(theMomentum),
93  thePosition(position),
94  nCollisions(0), nDecays(0),
95  thePotentialEnergy(0.),
96  rpCorrelated(false),
97  uncorrelatedMomentum(theMomentum.mag()),
98  theParticleBias(1.),
99  theHelicity(0.0),
100  emissionTime(0.0), outOfWell(false)
101  {
103  ID = nextID;
104  nextID++;
105  if(theEnergy <= 0.0) {
106  INCL_WARN("Particle with energy " << theEnergy << " created." << '\n');
107  }
108  setType(t);
110  }
111 
113  ThreeVector const &momentum, ThreeVector const &position)
114  : thePropagationEnergy(&theEnergy),
115  theMomentum(momentum),
116  thePropagationMomentum(&theMomentum),
117  theFrozenMomentum(theMomentum),
118  thePosition(position),
119  nCollisions(0), nDecays(0),
120  thePotentialEnergy(0.),
121  rpCorrelated(false),
122  uncorrelatedMomentum(theMomentum.mag()),
123  theParticleBias(1.),
124  theHelicity(0.0),
125  emissionTime(0.0), outOfWell(false)
126  {
128  ID = nextID;
129  nextID++;
130  setType(t);
131  if( isResonance() ) {
132  INCL_ERROR("Cannot create resonance without specifying its momentum four-vector." << '\n');
133  }
134  G4double energy = std::sqrt(theMomentum.mag2() + theMass*theMass);
135  theEnergy = energy;
137  }
138 
140  const G4double p2 = theMomentum.mag2();
142  if( newp2<0.0 ) {
143  INCL_ERROR("Particle has E^2 < m^2." << '\n' << print());
144  newp2 = 0.0;
145  theEnergy = theMass;
146  }
147 
148  theMomentum *= std::sqrt(newp2/p2);
149  return theMomentum;
150  }
151 
153  theEnergy = std::sqrt(theMomentum.mag2() + theMass*theMass);
154  return theEnergy;
155  }
156 
158  for(const_iterator i=begin(), e=end(); i!=e; ++i) {
159  (*i)->rotatePositionAndMomentum(angle, axis);
160  }
161  }
162 
163  void ParticleList::rotatePosition(const G4double angle, const ThreeVector &axis) const {
164  for(const_iterator i=begin(), e=end(); i!=e; ++i) {
165  (*i)->rotatePosition(angle, axis);
166  }
167  }
168 
169  void ParticleList::rotateMomentum(const G4double angle, const ThreeVector &axis) const {
170  for(const_iterator i=begin(), e=end(); i!=e; ++i) {
171  (*i)->rotateMomentum(angle, axis);
172  }
173  }
174 
175  void ParticleList::boost(const ThreeVector &b) const {
176  for(const_iterator i=begin(), e=end(); i!=e; ++i) {
177  (*i)->boost(b);
178  }
179  }
180 
182  if(G4int((*this).size())==0) return 1.;
183  std::vector<G4int> MergedVector;
184  for(ParticleIter i = (*this).begin(), e = (*this).end(); i!=e; ++i){
185  MergedVector = Particle::MergeVectorBias(MergedVector,(*i));
186  }
187  return Particle::getBiasFromVector(MergedVector);
188  }
189 
190  std::vector<G4int> ParticleList::getParticleListBiasVector() const {
191  std::vector<G4int> MergedVector;
192  if(G4int((*this).size())==0) return MergedVector;
193  for(ParticleIter i = (*this).begin(), e = (*this).end(); i!=e; ++i){
194  MergedVector = Particle::MergeVectorBias(MergedVector,(*i));
195  }
196  return MergedVector;
197  }
198 
200 // assert(G4int(Particle::INCLBiasVector.size())==nextBiasedCollisionID);
201  //assert(G4int(Particle::INCLBiasVector.Size())==nextBiasedCollisionID);
202 // assert(std::fabs(newBias - 1.) > 1E-6);
203  Particle::INCLBiasVector.push_back(newBias);
204  //Particle::INCLBiasVector.Push_back(newBias);
206  }
207 
208  G4double Particle::getBiasFromVector(std::vector<G4int> VectorBias) {
209  if(VectorBias.empty()) return 1.;
210 
211  G4double ParticleBias = 1.;
212 
213  for(G4int i=0; i<G4int(VectorBias.size()); i++){
214  ParticleBias *= Particle::INCLBiasVector[G4int(VectorBias[i])];
215  }
216 
217  return ParticleBias;
218  }
219 
220  std::vector<G4int> Particle::MergeVectorBias(Particle const * const p1, Particle const * const p2){
221  std::vector<G4int> MergedVectorBias;
222  std::vector<G4int> VectorBias1 = p1->getBiasCollisionVector();
223  std::vector<G4int> VectorBias2 = p2->getBiasCollisionVector();
224  G4int i = 0;
225  G4int j = 0;
226  if(VectorBias1.size()==0 && VectorBias2.size()==0) return MergedVectorBias;
227  else if(VectorBias1.size()==0) return VectorBias2;
228  else if(VectorBias2.size()==0) return VectorBias1;
229 
230  while(i < G4int(VectorBias1.size()) || j < G4int(VectorBias2.size())){
231  if(VectorBias1[i]==VectorBias2[j]){
232  MergedVectorBias.push_back(VectorBias1[i]);
233  i++;
234  j++;
235  if(i == G4int(VectorBias1.size())){
236  for(;j<G4int(VectorBias2.size());j++) MergedVectorBias.push_back(VectorBias2[j]);
237  }
238  if(j == G4int(VectorBias2.size())){
239  for(;i<G4int(VectorBias1.size());i++) MergedVectorBias.push_back(VectorBias1[i]);
240  }
241  } else if(VectorBias1[i]<VectorBias2[j]){
242  MergedVectorBias.push_back(VectorBias1[i]);
243  i++;
244  if(i == G4int(VectorBias1.size())){
245  for(;j<G4int(VectorBias2.size());j++) MergedVectorBias.push_back(VectorBias2[j]);
246  }
247  }
248  else {
249  MergedVectorBias.push_back(VectorBias2[j]);
250  j++;
251  if(j == G4int(VectorBias2.size())){
252  for(;i<G4int(VectorBias1.size());i++) MergedVectorBias.push_back(VectorBias1[i]);
253  }
254  }
255  }
256  return MergedVectorBias;
257  }
258 
259  std::vector<G4int> Particle::MergeVectorBias(std::vector<G4int> p1, Particle const * const p2){
260  std::vector<G4int> MergedVectorBias;
261  std::vector<G4int> VectorBias = p2->getBiasCollisionVector();
262  G4int i = 0;
263  G4int j = 0;
264  if(p1.size()==0 && VectorBias.size()==0) return MergedVectorBias;
265  else if(p1.size()==0) return VectorBias;
266  else if(VectorBias.size()==0) return p1;
267 
268  while(i < G4int(p1.size()) || j < G4int(VectorBias.size())){
269  if(p1[i]==VectorBias[j]){
270  MergedVectorBias.push_back(p1[i]);
271  i++;
272  j++;
273  if(i == G4int(p1.size())){
274  for(;j<G4int(VectorBias.size());j++) MergedVectorBias.push_back(VectorBias[j]);
275  }
276  if(j == G4int(VectorBias.size())){
277  for(;i<G4int(p1.size());i++) MergedVectorBias.push_back(p1[i]);
278  }
279  } else if(p1[i]<VectorBias[j]){
280  MergedVectorBias.push_back(p1[i]);
281  i++;
282  if(i == G4int(p1.size())){
283  for(;j<G4int(VectorBias.size());j++) MergedVectorBias.push_back(VectorBias[j]);
284  }
285  }
286  else {
287  MergedVectorBias.push_back(VectorBias[j]);
288  j++;
289  if(j == G4int(VectorBias.size())){
290  for(;i<G4int(p1.size());i++) MergedVectorBias.push_back(p1[i]);
291  }
292  }
293  }
294  return MergedVectorBias;
295  }
296 
298  G4double TotalBias = 1.;
299  for(G4int i=0; i<G4int(INCLBiasVector.size());i++) TotalBias *= Particle::INCLBiasVector[i];
300  return TotalBias;
301  }
302 
303  void Particle::setINCLBiasVector(std::vector<G4double> NewVector) {
304  Particle::INCLBiasVector = NewVector;
305  }
306 }
std::vector< G4int > getBiasCollisionVector() const
Get the vector list of biased vertices on the particle path.
ParticleList::const_iterator ParticleIter
void rotatePositionAndMomentum(const G4double angle, const ThreeVector &axis) const
const ThreeVector & adjustMomentumFromEnergy()
Rescale the momentum to match the total energy.
static std::vector< G4double > INCLBiasVector
Time ordered vector of all bias applied.
static void FillINCLBiasVector(G4double newBias)
void setMass(G4double mass)
G4INCL::ThreeVector theMomentum
ParticipantType theParticipantType
G4bool isResonance() const
Is it a resonance?
static G4double getBiasFromVector(std::vector< G4int > VectorBias)
#define G4ThreadLocal
Definition: tls.hh:69
void boost(const ThreeVector &b) const
static G4double getTotalBias()
General bias vector function.
void rotatePosition(const G4double angle, const ThreeVector &axis) const
double G4double
Definition: G4Types.hh:76
G4double getInvariantMass() const
Get the the particle invariant mass.
static std::vector< G4int > MergeVectorBias(Particle const *const p1, Particle const *const p2)
#define INCL_ERROR(x)
double energy
Definition: plottest35.C:25
std::string print() const
G4double theFrozenEnergy
void setType(ParticleType t)
#define INCL_WARN(x)
static void setINCLBiasVector(std::vector< G4double > NewVector)
int G4int
Definition: G4Types.hh:78
G4double mag2() const
G4double adjustEnergyFromMomentum()
Recompute the energy to match the momentum.
static G4ThreadLocal long nextID
static G4ThreadLocal G4int nextBiasedCollisionID
G4double getParticleListBias() const
std::vector< G4int > getParticleListBiasVector() const
void rotateMomentum(const G4double angle, const ThreeVector &axis) const