ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLParticle.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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  theNKaon(0),
77  theHelicity(0.0),
78  emissionTime(0.0),
79  outOfWell(false),
80  theMass(0.)
81  {
82  ID = nextID;
83  nextID++;
84  }
85 
88  : theEnergy(energy),
89  thePropagationEnergy(&theEnergy),
90  theFrozenEnergy(theEnergy),
91  theMomentum(momentum),
92  thePropagationMomentum(&theMomentum),
93  theFrozenMomentum(theMomentum),
94  thePosition(position),
95  nCollisions(0), nDecays(0),
96  thePotentialEnergy(0.),
97  rpCorrelated(false),
98  uncorrelatedMomentum(theMomentum.mag()),
99  theParticleBias(1.),
100  theNKaon(0),
101  theHelicity(0.0),
102  emissionTime(0.0), outOfWell(false)
103  {
105  ID = nextID;
106  nextID++;
107  if(theEnergy <= 0.0) {
108  INCL_WARN("Particle with energy " << theEnergy << " created." << '\n');
109  }
110  setType(t);
112  }
113 
115  ThreeVector const &momentum, ThreeVector const &position)
116  : thePropagationEnergy(&theEnergy),
117  theMomentum(momentum),
118  thePropagationMomentum(&theMomentum),
119  theFrozenMomentum(theMomentum),
120  thePosition(position),
121  nCollisions(0), nDecays(0),
122  thePotentialEnergy(0.),
123  rpCorrelated(false),
124  uncorrelatedMomentum(theMomentum.mag()),
125  theParticleBias(1.),
126  theNKaon(0),
127  theHelicity(0.0),
128  emissionTime(0.0), outOfWell(false)
129  {
131  ID = nextID;
132  nextID++;
133  setType(t);
134  if( isResonance() ) {
135  INCL_ERROR("Cannot create resonance without specifying its momentum four-vector." << '\n');
136  }
137  G4double energy = std::sqrt(theMomentum.mag2() + theMass*theMass);
138  theEnergy = energy;
140  }
141 
143  const G4double p2 = theMomentum.mag2();
145  if( newp2<0.0 ) {
146  INCL_ERROR("Particle has E^2 < m^2." << '\n' << print());
147  newp2 = 0.0;
148  theEnergy = theMass;
149  }
150 
151  theMomentum *= std::sqrt(newp2/p2);
152  return theMomentum;
153  }
154 
156  theEnergy = std::sqrt(theMomentum.mag2() + theMass*theMass);
157  return theEnergy;
158  }
159 
161  for(const_iterator i=begin(), e=end(); i!=e; ++i) {
162  (*i)->rotatePositionAndMomentum(angle, axis);
163  }
164  }
165 
166  void ParticleList::rotatePosition(const G4double angle, const ThreeVector &axis) const {
167  for(const_iterator i=begin(), e=end(); i!=e; ++i) {
168  (*i)->rotatePosition(angle, axis);
169  }
170  }
171 
172  void ParticleList::rotateMomentum(const G4double angle, const ThreeVector &axis) const {
173  for(const_iterator i=begin(), e=end(); i!=e; ++i) {
174  (*i)->rotateMomentum(angle, axis);
175  }
176  }
177 
178  void ParticleList::boost(const ThreeVector &b) const {
179  for(const_iterator i=begin(), e=end(); i!=e; ++i) {
180  (*i)->boost(b);
181  }
182  }
183 
185  if(G4int((*this).size())==0) return 1.;
186  std::vector<G4int> MergedVector;
187  for(ParticleIter i = (*this).begin(), e = (*this).end(); i!=e; ++i){
188  MergedVector = Particle::MergeVectorBias(MergedVector,(*i));
189  }
190  return Particle::getBiasFromVector(MergedVector);
191  }
192 
193  std::vector<G4int> ParticleList::getParticleListBiasVector() const {
194  std::vector<G4int> MergedVector;
195  if(G4int((*this).size())==0) return MergedVector;
196  for(ParticleIter i = (*this).begin(), e = (*this).end(); i!=e; ++i){
197  MergedVector = Particle::MergeVectorBias(MergedVector,(*i));
198  }
199  return MergedVector;
200  }
201 
203 // assert(G4int(Particle::INCLBiasVector.size())==nextBiasedCollisionID);
204  //assert(G4int(Particle::INCLBiasVector.Size())==nextBiasedCollisionID);
205 // assert(std::fabs(newBias - 1.) > 1E-6);
206  Particle::INCLBiasVector.push_back(newBias);
207  //Particle::INCLBiasVector.Push_back(newBias);
209  }
210 
211  G4double Particle::getBiasFromVector(std::vector<G4int> VectorBias) {
212  if(VectorBias.empty()) return 1.;
213 
214  G4double ParticleBias = 1.;
215 
216  for(G4int i=0; i<G4int(VectorBias.size()); i++){
217  ParticleBias *= Particle::INCLBiasVector[G4int(VectorBias[i])];
218  }
219 
220  return ParticleBias;
221  }
222 
223  std::vector<G4int> Particle::MergeVectorBias(Particle const * const p1, Particle const * const p2){
224  std::vector<G4int> MergedVectorBias;
225  std::vector<G4int> VectorBias1 = p1->getBiasCollisionVector();
226  std::vector<G4int> VectorBias2 = p2->getBiasCollisionVector();
227  G4int i = 0;
228  G4int j = 0;
229  if(VectorBias1.size()==0 && VectorBias2.size()==0) return MergedVectorBias;
230  else if(VectorBias1.size()==0) return VectorBias2;
231  else if(VectorBias2.size()==0) return VectorBias1;
232 
233  while(i < G4int(VectorBias1.size()) || j < G4int(VectorBias2.size())){
234  if(VectorBias1[i]==VectorBias2[j]){
235  MergedVectorBias.push_back(VectorBias1[i]);
236  i++;
237  j++;
238  if(i == G4int(VectorBias1.size())){
239  for(;j<G4int(VectorBias2.size());j++) MergedVectorBias.push_back(VectorBias2[j]);
240  }
241  else if(j == G4int(VectorBias2.size())){
242  for(;i<G4int(VectorBias1.size());i++) MergedVectorBias.push_back(VectorBias1[i]);
243  }
244  } else if(VectorBias1[i]<VectorBias2[j]){
245  MergedVectorBias.push_back(VectorBias1[i]);
246  i++;
247  if(i == G4int(VectorBias1.size())){
248  for(;j<G4int(VectorBias2.size());j++) MergedVectorBias.push_back(VectorBias2[j]);
249  }
250  }
251  else {
252  MergedVectorBias.push_back(VectorBias2[j]);
253  j++;
254  if(j == G4int(VectorBias2.size())){
255  for(;i<G4int(VectorBias1.size());i++) MergedVectorBias.push_back(VectorBias1[i]);
256  }
257  }
258  }
259  return MergedVectorBias;
260  }
261 
262  std::vector<G4int> Particle::MergeVectorBias(std::vector<G4int> p1, Particle const * const p2){
263  std::vector<G4int> MergedVectorBias;
264  std::vector<G4int> VectorBias = p2->getBiasCollisionVector();
265  G4int i = 0;
266  G4int j = 0;
267  if(p1.size()==0 && VectorBias.size()==0) return MergedVectorBias;
268  else if(p1.size()==0) return VectorBias;
269  else if(VectorBias.size()==0) return p1;
270 
271  while(i < G4int(p1.size()) || j < G4int(VectorBias.size())){
272  if(p1[i]==VectorBias[j]){
273  MergedVectorBias.push_back(p1[i]);
274  i++;
275  j++;
276  if(i == G4int(p1.size())){
277  for(;j<G4int(VectorBias.size());j++) MergedVectorBias.push_back(VectorBias[j]);
278  }
279  else if(j == G4int(VectorBias.size())){
280  for(;i<G4int(p1.size());i++) MergedVectorBias.push_back(p1[i]);
281  }
282  } else if(p1[i]<VectorBias[j]){
283  MergedVectorBias.push_back(p1[i]);
284  i++;
285  if(i == G4int(p1.size())){
286  for(;j<G4int(VectorBias.size());j++) MergedVectorBias.push_back(VectorBias[j]);
287  }
288  }
289  else {
290  MergedVectorBias.push_back(VectorBias[j]);
291  j++;
292  if(j == G4int(VectorBias.size())){
293  for(;i<G4int(p1.size());i++) MergedVectorBias.push_back(p1[i]);
294  }
295  }
296  }
297  return MergedVectorBias;
298  }
299 
301  G4double TotalBias = 1.;
302  for(G4int i=0; i<G4int(INCLBiasVector.size());i++) TotalBias *= Particle::INCLBiasVector[i];
303  return TotalBias;
304  }
305 
306  void Particle::setINCLBiasVector(std::vector<G4double> NewVector) {
307  Particle::INCLBiasVector = NewVector;
308  }
309 }