ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4QMDNucleus.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4QMDNucleus.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 // 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
27 //
28 #include <numeric>
29 
30 #include "G4QMDNucleus.hh"
31 #include "G4Pow.hh"
32 #include "G4SystemOfUnits.hh"
33 #include "G4Proton.hh"
34 #include "G4Neutron.hh"
35 #include "G4NucleiProperties.hh"
36 #include "G4HadronicException.hh"
37 
39 {
41  hbc = parameters->Get_hbc();
42 
43  jj = 0; // will be calcualted in CalEnergyAndAngularMomentumInCM;
44  potentialEnergy = 0.0; // will be set through set method
45  excitationEnergy = 0.0;
46 }
47 
48 
49 
50 //G4QMDNucleus::~G4QMDNucleus()
51 //{
52 // ;
53 //}
54 
55 
57 {
58  G4LorentzVector p( 0 );
59  std::vector< G4QMDParticipant* >::iterator it;
60  for ( it = participants.begin() ; it != participants.end() ; it++ )
61  p += (*it)->Get4Momentum();
62 
63  return p;
64 }
65 
66 
67 
69 {
70 
71  G4int A = 0;
72  std::vector< G4QMDParticipant* >::iterator it;
73  for ( it = participants.begin() ; it != participants.end() ; it++ )
74  {
75  if ( (*it)->GetDefinition() == G4Proton::Proton()
76  || (*it)->GetDefinition() == G4Neutron::Neutron() )
77  A++;
78  }
79 
80  if ( A == 0 ) {
81  throw G4HadronicException(__FILE__, __LINE__, "G4QMDNucleus has the mass number of 0!");
82  }
83 
84  return A;
85 }
86 
87 
88 
90 {
91  G4int Z = 0;
92  std::vector< G4QMDParticipant* >::iterator it;
93  for ( it = participants.begin() ; it != participants.end() ; it++ )
94  {
95  if ( (*it)->GetDefinition() == G4Proton::Proton() )
96  Z++;
97  }
98  return Z;
99 }
100 
101 
102 
104 {
105 
107 
108  if ( mass == 0.0 )
109  {
110 
111  G4int Z = GetAtomicNumber();
112  G4int A = GetMassNumber();
113  G4int N = A - Z;
114 
115 // Weizsacker-Bethe
116 
117  G4double Av = 16*MeV;
118  G4double As = 17*MeV;
119  G4double Ac = 0.7*MeV;
120  G4double Asym = 23*MeV;
121 
122  G4double BE = Av * A
123  - As * G4Pow::GetInstance()->A23 ( G4double ( A ) )
124  - Ac * Z*Z/G4Pow::GetInstance()->A13 ( G4double ( A ) )
125  - Asym * ( N - Z )* ( N - Z ) / A;
126 
127  mass = Z * G4Proton::Proton()->GetPDGMass()
128  + N * G4Neutron::Neutron()->GetPDGMass()
129  - BE;
130 
131  }
132 
133  return mass;
134 }
135 
136 
137 
139 {
140 
141  //G4cout << "CalEnergyAndAngularMomentumInCM " << this->GetAtomicNumber() << " " << GetMassNumber() << G4endl;
142 
143  G4double gamma = Get4Momentum().gamma();
144  G4ThreeVector beta = Get4Momentum().v()/ Get4Momentum().e();
145 
146  G4ThreeVector pcm0( 0.0 ) ;
147 
149  pcm.resize( n );
150 
151  for ( G4int i= 0; i < n ; i++ )
152  {
154 
155  G4double trans = gamma / ( gamma + 1.0 ) * p_i * beta;
156  pcm[i] = p_i - trans*beta;
157 
158  pcm0 += pcm[i];
159  }
160 
161  pcm0 = pcm0 / double ( n );
162 
163  //G4cout << "pcm0 " << pcm0 << G4endl;
164 
165  for ( G4int i= 0; i < n ; i++ )
166  {
167  pcm[i] += -pcm0;
168  //G4cout << "pcm " << i << " " << pcm[i] << G4endl;
169  }
170 
171 
172  G4double tmass = 0;
173  G4ThreeVector rcm0( 0.0 ) ;
174  rcm.resize( n );
175  es.resize( n );
176 
177  for ( G4int i= 0; i < n ; i++ )
178  {
180  G4double trans = gamma / ( gamma + 1.0 ) * ri * beta;
181 
182  es[i] = std::sqrt ( G4Pow::GetInstance()->powN ( GetParticipant( i )->GetMass() , 2 ) + pcm[i]*pcm[i] );
183 
184  rcm[i] = ri + trans*beta;
185 
186  rcm0 += rcm[i]*es[i];
187 
188  tmass += es[i];
189  }
190 
191  rcm0 = rcm0/tmass;
192 
193  for ( G4int i= 0; i < n ; i++ )
194  {
195  rcm[i] += -rcm0;
196  //G4cout << "rcm " << i << " " << rcm[i] << G4endl;
197  }
198 
199 // Angluar momentum
200 
201  G4ThreeVector rl ( 0.0 );
202  for ( G4int i= 0; i < n ; i++ )
203  {
204  rl += rcm[i].cross ( pcm[i] );
205  }
206 
207  jj = int ( std::sqrt ( rl*rl / hbc ) + 0.5 );
208 
209 
210 // kinetic energy per nucleon in CM
211 
212  G4double totalMass = 0.0;
213  for ( G4int i= 0; i < n ; i++ )
214  {
215  // following two lines are equivalent
216  //totalMass += GetParticipant( i )->GetDefinition()->GetPDGMass()/GeV;
217  totalMass += GetParticipant( i )->GetMass();
218  }
219 
220  //G4double kineticEnergyPerNucleon = ( std::accumulate ( es.begin() , es.end() , 0.0 ) - totalMass )/n;
221 
222 // Total (not per nucleion ) Binding Energy
223  G4double bindingEnergy = ( std::accumulate ( es.begin() , es.end() , 0.0 ) -totalMass ) + potentialEnergy;
224 
225  //G4cout << "KineticEnergyPerNucleon in GeV " << kineticEnergyPerNucleon << G4endl;
226  //G4cout << "KineticEnergySum in GeV " << std::accumulate ( es.begin() , es.end() , 0.0 ) - totalMass << G4endl;
227  //G4cout << "PotentialEnergy in GeV " << potentialEnergy << G4endl;
228  //G4cout << "BindingEnergy in GeV " << bindingEnergy << G4endl;
229  //G4cout << "G4BindingEnergy in GeV " << G4NucleiProperties::GetBindingEnergy( GetAtomicNumber() , GetMassNumber() )/GeV << G4endl;
230 
232  //G4cout << "excitationEnergy in GeV " << excitationEnergy << G4endl;
233  if ( excitationEnergy < 0 ) excitationEnergy = 0.0;
234 
235 }