ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LENDElastic.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4LENDElastic.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 
27 #include "G4LENDElastic.hh"
28 #include "G4Pow.hh"
29 #include "G4PhysicalConstants.hh"
30 #include "G4SystemOfUnits.hh"
31 #include "G4Nucleus.hh"
32 #include "G4IonTable.hh"
33 
34 //extern "C" double MyRNG(void*) { return drand48(); }
35 //extern "C" double MyRNG(void*) { return CLHEP::HepRandom::getTheEngine()->flat(); }
36 
38 {
39 
40  G4double temp = aTrack.GetMaterial()->GetTemperature();
41 
42  //G4int iZ = int ( aTarg.GetZ() );
43  //G4int iA = int ( aTarg.GetN() );
44  //migrate to integer A and Z (GetN_asInt returns number of neutrons in the nucleus since this)
45  G4int iZ = aTarg.GetZ_asInt();
46  G4int iA = aTarg.GetA_asInt();
47  G4int iM = 0;
48  if ( aTarg.GetIsotope() != NULL ) {
49  iM = aTarg.GetIsotope()->Getm();
50  }
51 
52  G4double ke = aTrack.GetKineticEnergy();
53 
54  G4HadFinalState* theResult = &theParticleChange;
55  theResult->Clear();
56 
58  if ( aTarget == NULL ) return returnUnchanged( aTrack , theResult );
59 
60  G4double aMu = aTarget->getElasticFinalState( ke*MeV, temp, MyRNG , NULL );
61 
63  G4double theta = std::acos( aMu );
64  //G4double sinth = std::sin( theta );
65 
66  G4ReactionProduct theNeutron( aTrack.GetDefinition() );
67  theNeutron.SetMomentum( aTrack.Get4Momentum().vect() );
68  theNeutron.SetKineticEnergy( ke );
69 
70  //G4ParticleDefinition* pd = G4IonTable::GetIonTable()->GetIon( iZ , iA , iM );
71  //TK 170509 Fix for the case of excited isomer target
72  G4double EE = 0.0;
73  if ( iM != 0 ) {
75  }
76  G4ParticleDefinition* target_pd = G4IonTable::GetIonTable()->GetIon( iZ , iA , EE );
77  G4ReactionProduct theTarget( target_pd );
78 
79  G4double mass = target_pd->GetPDGMass();
80 
81 // add Thermal motion
82  G4double kT = k_Boltzmann*temp;
83  G4ThreeVector v ( G4RandGauss::shoot() * std::sqrt( kT*mass )
84  , G4RandGauss::shoot() * std::sqrt( kT*mass )
85  , G4RandGauss::shoot() * std::sqrt( kT*mass ) );
86  theTarget.SetMomentum( v );
87 
88  G4ThreeVector the3Neutron = theNeutron.GetMomentum();
89  G4double nEnergy = theNeutron.GetTotalEnergy();
90  G4ThreeVector the3Target = theTarget.GetMomentum();
91  G4double tEnergy = theTarget.GetTotalEnergy();
92  G4ReactionProduct theCMS;
93  G4double totE = nEnergy+tEnergy;
94  G4ThreeVector the3CMS = the3Target+the3Neutron;
95  theCMS.SetMomentum(the3CMS);
96  G4double cmsMom = std::sqrt(the3CMS*the3CMS);
97  G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
98  theCMS.SetMass(sqrts);
99  theCMS.SetTotalEnergy(totE);
100 
101  theNeutron.Lorentz(theNeutron, theCMS);
102  theTarget.Lorentz(theTarget, theCMS);
103  G4double en = theNeutron.GetTotalMomentum(); // already in CMS.
104  G4ThreeVector cms3Mom=theNeutron.GetMomentum(); // for neutron direction in CMS
105  G4double cms_theta=cms3Mom.theta();
106  G4double cms_phi=cms3Mom.phi();
107  G4ThreeVector tempVector;
108  tempVector.setX( std::cos(theta)*std::sin(cms_theta)*std::cos(cms_phi)
109  +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::cos(cms_phi)
110  -std::sin(theta)*std::sin(phi)*std::sin(cms_phi) );
111  tempVector.setY( std::cos(theta)*std::sin(cms_theta)*std::sin(cms_phi)
112  +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::sin(cms_phi)
113  +std::sin(theta)*std::sin(phi)*std::cos(cms_phi) );
114  tempVector.setZ( std::cos(theta)*std::cos(cms_theta)
115  -std::sin(theta)*std::cos(phi)*std::sin(cms_theta) );
116  tempVector *= en;
117  theNeutron.SetMomentum(tempVector);
118  theTarget.SetMomentum(-tempVector);
119  G4double tP = theTarget.GetTotalMomentum();
120  G4double tM = theTarget.GetMass();
121  theTarget.SetTotalEnergy(std::sqrt((tP+tM)*(tP+tM)-2.*tP*tM));
122 
123 
124  theNeutron.Lorentz(theNeutron, -1.*theCMS);
125  theTarget.Lorentz(theTarget, -1.*theCMS);
126 
127 //110913 Add Protection for very low energy (1e-6eV) scattering
128  if ( theNeutron.GetKineticEnergy() <= 0 )
129  {
130  theNeutron.SetTotalEnergy ( theNeutron.GetMass() * ( 1 + G4Pow::GetInstance()->powA( 10 , -15.65 ) ) );
131  }
132 
133  if ( theTarget.GetKineticEnergy() < 0 )
134  {
135  theTarget.SetTotalEnergy ( theTarget.GetMass() * ( 1 + G4Pow::GetInstance()->powA( 10 , -15.65 ) ) );
136  }
137 //110913 END
138 
139  theResult->SetEnergyChange(theNeutron.GetKineticEnergy());
140  theResult->SetMomentumChange(theNeutron.GetMomentum().unit());
141  G4DynamicParticle* theRecoil = new G4DynamicParticle;
142 
143  theRecoil->SetDefinition( target_pd );
144  theRecoil->SetMomentum( theTarget.GetMomentum() );
145  theResult->AddSecondary( theRecoil );
146 
147  return theResult;
148 
149 }
150