| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 //
28
29 // Author : Gunter Folger March 2007
30 // Modified by Mikhail Kossov. Apr2009, E/M co
31 // Class Description
32 // Final state production model for theoretica
33 // quasi elastic scattering in geant4;
34 // Class Description - End
35 //
36 // Modified:
37 // 20110805 M. Kelsey -- Follow change to G4V
38 // 20110808 M. Kelsey -- Move #includes from
39
40 #include "G4QuasiElasticChannel.hh"
41
42 #include "G4Fancy3DNucleus.hh"
43 #include "G4DynamicParticle.hh"
44 #include "G4HadTmpUtil.hh" /* lrint */
45 #include "G4KineticTrack.hh"
46 #include "G4KineticTrackVector.hh"
47 #include "G4LorentzVector.hh"
48 #include "G4Neutron.hh"
49 #include "G4Nucleon.hh"
50 #include "G4Nucleus.hh"
51 #include "G4ParticleDefinition.hh"
52 #include "G4ParticleTable.hh"
53 #include "G4IonTable.hh"
54 #include "G4QuasiElRatios.hh"
55 #include "globals.hh"
56 #include <vector>
57 #include "G4PhysicsModelCatalog.hh"
58
59 //#define debug_scatter
60
61
62 G4QuasiElasticChannel::G4QuasiElasticChannel()
63 : G4HadronicInteraction("QuasiElastic"),
64 theQuasiElastic(new G4QuasiElRatios),
65 the3DNucleus(new G4Fancy3DNucleus),
66 secID(-1) {
67 secID = G4PhysicsModelCatalog::GetModelID( "
68 }
69
70 G4QuasiElasticChannel::~G4QuasiElasticChannel(
71 {
72 delete the3DNucleus;
73 delete theQuasiElastic;
74 }
75
76 G4double G4QuasiElasticChannel::GetFraction(G4
77 const G4DynamicParticle & thePrimary)
78 {
79 #ifdef debug_scatter
80 G4cout << "G4QuasiElasticChannel:: P=" <
81 << ", pPDG=" << thePrimary.GetDef
82 << ", Z = " << theNucleus.GetZ_a
83 << ", N = " << theNucleus.GetN_a
84 << ", A = " << theNucleus.GetA_a
85 #endif
86
87 std::pair<G4double,G4double> ratios;
88 ratios=theQuasiElastic->GetRatios(thePrimary
89 thePrimary
90 theNucleus
91 theNucleus
92 #ifdef debug_scatter
93 G4cout << "G4QuasiElasticChannel::ratios
94 << " = " << ratios.first*ratios.
95 #endif
96
97 return ratios.first*ratios.second;
98 }
99
100 G4KineticTrackVector * G4QuasiElasticChannel::
101
102 {
103 G4int A=theNucleus.GetA_asInt();
104 G4int Z=theNucleus.GetZ_asInt();
105 // build Nucleus and choose random nucleon
106 the3DNucleus->Init(theNucleus.GetA_asInt(),t
107 const std::vector<G4Nucleon>& nucleons=the3D
108 G4double targetNucleusMass=the3DNucleus->Get
109 G4LorentzVector targetNucleus4Mom(0.,0.,0.,t
110 G4int index;
111 do {
112 index=G4lrint((A-1)*G4UniformRand());
113 } while (index < 0 || index >= static_cast<G
114
115 const G4ParticleDefinition * pDef= nucleons[
116
117 G4int resA=A - 1;
118 G4int resZ=Z - static_cast<int>(pDef->GetPDG
119 const G4ParticleDefinition* resDef;
120 G4double residualNucleusMass;
121 if(resZ)
122 {
123 resDef=G4ParticleTable::GetParticleTable()
124 residualNucleusMass=resDef->GetPDGMass();
125 }
126 else {
127 resDef=G4Neutron::Neutron();
128 residualNucleusMass=resA * G4Neutron::Neut
129 }
130 #ifdef debug_scatter
131 G4cout<<"G4QElChan::Scatter: neutron - pr
132 <<pDef->GetParticleName()<<G4endl;
133 #endif
134
135 G4LorentzVector pNucleon=nucleons[index].Get
136 G4double residualNucleusEnergy=std::sqrt(sqr
137 pNu
138 pNucleon.setE(targetNucleusMass-residualNucl
139 G4LorentzVector residualNucleus4Mom=targetNu
140
141 std::pair<G4LorentzVector,G4LorentzVector> r
142
143 result=theQuasiElastic->Scatter(pDef->GetPDG
144 thePrimary.G
145 thePrimary.G
146 G4LorentzVector scatteredHadron4Mom;
147 if (result.first.e() > 0.)
148 scatteredHadron4Mom=result.second;
149 else { //scatter failed
150 //G4cout << "Warning - G4QuasiElasticChann
151 //return 0; //no scatter
152 scatteredHadron4Mom=thePrimary.Get4Momentu
153 residualNucleus4Mom=G4LorentzVector(0.,0.,
154 resDef=G4ParticleTable::GetParticleTable()
155 }
156
157 #ifdef debug_scatter
158 G4LorentzVector EpConservation=pNucleon+theP
159 - result.firs
160 if ( (EpConservation.vect().mag2() > .01*M
161 || (std::abs(EpConservation.e()) > 0.1 *
162 {
163 G4cout << "Warning - G4QuasiElasticChannel
164 << EpConservation << G4endl;
165 }
166 #endif
167
168 G4KineticTrackVector * ktv = new G4KineticTr
169 G4KineticTrack * sPrim=new G4KineticTrack(th
170 0.
171 sPrim->SetCreatorModelID( secID );
172 ktv->push_back(sPrim);
173 if (result.first.e() > 0.)
174 {
175 G4KineticTrack * sNuc=new G4KineticTrack(p
176 sNuc->SetCreatorModelID( secID );
177 ktv->push_back(sNuc);
178 }
179 if(resZ || resA==1) // For the only neutron
180 {
181 G4KineticTrack * rNuc=new G4KineticTrack(r
182 0.,G4ThreeVector(0)
183 rNuc->SetCreatorModelID( secID );
184 ktv->push_back(rNuc);
185 }
186 else // The residual nucleus consists of onl
187 {
188 residualNucleus4Mom/=resA; // Split 4-
189 for(G4int in=0; in<resA; in++) // Loop ove
190 {
191 G4KineticTrack* rNuc=new G4KineticTrack(
192 0.,G4ThreeVector(0)
193 rNuc->SetCreatorModelID( secID );
194 ktv->push_back(rNuc);
195 }
196 }
197 #ifdef debug_scatter
198 G4cout<<"G4QElC::Scat: Nucleon: "<<result.fi
199 G4cout<<"G4QElC::Scat: Project: "<<result.se
200 #endif
201 return ktv;
202 }
203