Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/qmd/src/G4QMDNucleus.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 
 38 G4QMDNucleus::G4QMDNucleus()
 39 {
 40    G4QMDParameters* parameters = G4QMDParameters::GetInstance();
 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 
 56 G4LorentzVector G4QMDNucleus::Get4Momentum()
 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 
 68 G4int G4QMDNucleus::GetMassNumber()
 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 
 89 G4int G4QMDNucleus::GetAtomicNumber()
 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 
103 G4double G4QMDNucleus::GetNuclearMass()
104 {
105 
106    G4double mass = G4NucleiProperties::GetNuclearMass( GetMassNumber() , GetAtomicNumber() );
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 
138 void G4QMDNucleus::CalEnergyAndAngularMomentumInCM()
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 
148    G4int n = GetTotalNumberOfParticipant();
149    pcm.resize( n );
150 
151    for ( G4int i= 0; i < n ; i++ ) 
152    {
153       G4ThreeVector p_i = GetParticipant( i )->GetMomentum();
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    {
179       G4ThreeVector ri = GetParticipant( i )->GetPosition();
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 // Angular 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 // DHW: move hbc outside of sqrt to get correct units
208 //  jj = int ( std::sqrt ( rl*rl / hbc ) + 0.5 );
209 
210    jj = int (std::sqrt(rl*rl)/hbc + 0.5);
211 
212 // kinetic energy per nucleon in CM
213 
214    G4double totalMass = 0.0;
215    for ( G4int i= 0; i < n ; i++ ) 
216    {
217       // following two lines are equivalent
218       //totalMass += GetParticipant( i )->GetDefinition()->GetPDGMass()/GeV;
219       totalMass += GetParticipant( i )->GetMass();
220    }
221 
222    //G4double kineticEnergyPerNucleon = ( std::accumulate ( es.begin() , es.end() , 0.0 ) - totalMass )/n;
223 
224 // Total (not per nucleion ) Binding Energy 
225    G4double bindingEnergy =  ( std::accumulate ( es.begin() , es.end() , 0.0 ) -totalMass ) + potentialEnergy;
226 
227    //G4cout << "KineticEnergyPerNucleon in GeV " << kineticEnergyPerNucleon << G4endl;
228    //G4cout << "KineticEnergySum in GeV " << std::accumulate ( es.begin() , es.end() , 0.0 ) - totalMass << G4endl;
229    //G4cout << "PotentialEnergy in GeV " << potentialEnergy << G4endl;
230    //G4cout << "BindingEnergy in GeV " << bindingEnergy << G4endl;
231    //G4cout << "G4BindingEnergy in GeV " << G4NucleiProperties::GetBindingEnergy( GetAtomicNumber() , GetMassNumber() )/GeV << G4endl;
232 
233    excitationEnergy = bindingEnergy + G4NucleiProperties::GetBindingEnergy( GetMassNumber() , GetAtomicNumber() )/GeV;
234    //G4cout << "excitationEnergy in GeV " << excitationEnergy << G4endl;
235    if ( excitationEnergy < 0 ) excitationEnergy = 0.0; 
236 
237 }
238