1 | #include "G4BinaryLightIonReaction.hh" |
---|
2 | #include "G4Nucleus.hh" |
---|
3 | #include "G4Proton.hh" |
---|
4 | #include "G4PionPlus.hh" |
---|
5 | #include "G4ThreeVector.hh" |
---|
6 | #include "G4LorentzVector.hh" |
---|
7 | #include "G4DynamicParticle.hh" |
---|
8 | |
---|
9 | #include "G4LeptonConstructor.hh" |
---|
10 | #include "G4IonConstructor.hh" |
---|
11 | #include "G4Gamma.hh" |
---|
12 | |
---|
13 | G4V3DNucleus * global3DNucleus; |
---|
14 | |
---|
15 | |
---|
16 | void Init(G4int code, vector<G4ParticleDefinition *> & particles); |
---|
17 | void GetParam(G4int & A, G4int & Z, G4double & pMin, G4double & pMax, |
---|
18 | G4int &nEvents, G4int & particleCode); |
---|
19 | G4ThreeVector GetSpherePoint(G4double r); |
---|
20 | G4ThreeVector GetSurfacePoint(G4double r); |
---|
21 | |
---|
22 | int main(int argc, char * argv[]) |
---|
23 | { |
---|
24 | G4String a=G4Proton::ProtonDefinition()->GetParticleName(); |
---|
25 | G4cout << " At Start " << a << G4endl; |
---|
26 | G4cout.setf( std::ios::scientific, std::ios::floatfield ); |
---|
27 | // totalRes.Start(); |
---|
28 | G4LeptonConstructor leptons; |
---|
29 | leptons.ConstructParticle(); |
---|
30 | G4IonConstructor ions; |
---|
31 | ions.ConstructParticle(); |
---|
32 | |
---|
33 | G4int A, Z, nEvents, particleCode; |
---|
34 | G4double pMin, pMax; |
---|
35 | vector<G4ParticleDefinition *> particles; |
---|
36 | if(argc == 7) |
---|
37 | { |
---|
38 | A = atoi(argv[1]); |
---|
39 | Z = atoi(argv[2]); |
---|
40 | pMin = G4double(atoi(argv[3]))*MeV; |
---|
41 | pMax = G4double(atoi(argv[4]))*MeV; |
---|
42 | nEvents = atoi(argv[5]); |
---|
43 | particleCode = atoi(argv[6]); |
---|
44 | } |
---|
45 | else |
---|
46 | GetParam(A, Z, pMin, pMax, nEvents, particleCode); |
---|
47 | |
---|
48 | if(pMin > pMax) |
---|
49 | { |
---|
50 | G4double pTmp = pMin; |
---|
51 | pMin = pMax; |
---|
52 | pMax = pTmp; |
---|
53 | } |
---|
54 | |
---|
55 | Init(particleCode, particles); |
---|
56 | |
---|
57 | G4BinaryCascade hkm; |
---|
58 | hkm.SetDeExcitation(0); |
---|
59 | |
---|
60 | for(G4int event = 0; event < nEvents; ++event) |
---|
61 | { |
---|
62 | G4cerr << " ---- event " << event << "----------------------"<< G4endl; |
---|
63 | G4V3DNucleus * fancyNucleus = new G4Fancy3DNucleus; // for Propagate() |
---|
64 | fancyNucleus->Init(16, 8); |
---|
65 | G4V3DNucleus * projectile = new G4Fancy3DNucleus; |
---|
66 | projectile->Init(12, 6); |
---|
67 | G4double radius = fancyNucleus->GetOuterRadius(); |
---|
68 | G4double impactMax = fancyNucleus->GetOuterRadius()+projectile->GetOuterRadius(); |
---|
69 | // position |
---|
70 | G4double aX=(2.*G4UniformRand()-1.)*impactMax; |
---|
71 | G4double aY=(2.*G4UniformRand()-1.)*impactMax; |
---|
72 | |
---|
73 | G4ThreeVector pos(aX, aY, -20.*fermi); |
---|
74 | // momentum |
---|
75 | G4double mass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(6, 12); |
---|
76 | G4double p = 280*MeV; |
---|
77 | G4double e = sqrt(mass*mass+p*p); |
---|
78 | G4LorentzVector mom(0, p, 0, e); |
---|
79 | |
---|
80 | // build the projectile |
---|
81 | G4KineticTrackVector * ktv = new G4KineticTrackVector; |
---|
82 | G4ThreeVector beta = mom.beta(); |
---|
83 | projectile->DoLorentzBoost(beta); |
---|
84 | projectile->StartLoop(); |
---|
85 | G4Nucleon * aNuc; |
---|
86 | while( (aNuc=projectile->GetNextNucleon()) ) |
---|
87 | { |
---|
88 | G4KineticTrack * it = |
---|
89 | new G4KineticTrack(aNuc->GetDefinition(), 0, |
---|
90 | aNuc->GetPosition()+pos, const_cast<G4LorentzVector &>(aNuc->GetMomentum()) ); |
---|
91 | ktv->push_back(it); |
---|
92 | } |
---|
93 | |
---|
94 | // do the reactions |
---|
95 | G4ReactionProductVector *result=hkm.Propagate(ktv, fancyNucleus); |
---|
96 | |
---|
97 | delete fancyNucleus; |
---|
98 | delete projectile; |
---|
99 | |
---|
100 | G4ReactionProductVector::iterator result_iter; |
---|
101 | G4int isec=0; |
---|
102 | for(result_iter = result->begin();result_iter != result->end(); |
---|
103 | ++result_iter) |
---|
104 | { |
---|
105 | G4ThreeVector pFinal= (*result_iter)->GetMomentum(); |
---|
106 | |
---|
107 | G4cout << " Final particle # "<< ++isec << ": " |
---|
108 | << pFinal.x() << " " |
---|
109 | << pFinal.y() << " " |
---|
110 | << pFinal.z() << " " |
---|
111 | << (*result_iter)->GetTotalEnergy() << " " |
---|
112 | << G4endl; |
---|
113 | } |
---|
114 | } |
---|
115 | cout.flush(); |
---|
116 | return 0; |
---|
117 | } |
---|
118 | |
---|
119 | /* |
---|
120 | |
---|
121 | void PrintResourceUsage(ResourceMgr & resMgr) |
---|
122 | { |
---|
123 | double userTime = resMgr.GetUserTime(); |
---|
124 | double sysTime = resMgr.GetSysTime(); |
---|
125 | double cpuTime = userTime+sysTime; |
---|
126 | cout << "cpu time: " << cpuTime << " (" << userTime << "+" |
---|
127 | << sysTime << ")" << endl; |
---|
128 | } |
---|
129 | |
---|
130 | |
---|
131 | */ |
---|
132 | |
---|
133 | void Init(G4int code, vector<G4ParticleDefinition *> & particles) |
---|
134 | { |
---|
135 | if(code == 0) |
---|
136 | { |
---|
137 | particles.push_back(G4Gamma::Gamma()); |
---|
138 | return; |
---|
139 | } |
---|
140 | |
---|
141 | G4int count = 0; |
---|
142 | while(code != 0) |
---|
143 | { |
---|
144 | if(code & 1) |
---|
145 | { |
---|
146 | switch (count) |
---|
147 | { |
---|
148 | case 0: |
---|
149 | particles.push_back(G4Proton::Proton()); |
---|
150 | break; |
---|
151 | case 1: |
---|
152 | particles.push_back(G4Neutron::Neutron()); |
---|
153 | break; |
---|
154 | case 2: |
---|
155 | particles.push_back(G4AntiProton::AntiProton()); |
---|
156 | break; |
---|
157 | case 3: |
---|
158 | particles.push_back(G4PionPlus::PionPlus()); |
---|
159 | break; |
---|
160 | case 4: |
---|
161 | particles.push_back(G4PionZero::PionZero()); |
---|
162 | break; |
---|
163 | case 5: |
---|
164 | particles.push_back(G4PionMinus::PionMinus()); |
---|
165 | break; |
---|
166 | case 6: |
---|
167 | particles.push_back(G4KaonPlus::KaonPlus()); |
---|
168 | break; |
---|
169 | case 7: |
---|
170 | particles.push_back(G4KaonZero::KaonZero()); |
---|
171 | break; |
---|
172 | case 8: |
---|
173 | particles.push_back(G4KaonMinus::KaonMinus()); |
---|
174 | break; |
---|
175 | case 9: |
---|
176 | particles.push_back(G4SigmaPlus::SigmaPlus()); |
---|
177 | break; |
---|
178 | case 10: |
---|
179 | particles.push_back(G4SigmaZero::SigmaZero()); |
---|
180 | break; |
---|
181 | case 11: |
---|
182 | particles.push_back(G4SigmaMinus::SigmaMinus()); |
---|
183 | break; |
---|
184 | } |
---|
185 | } |
---|
186 | ++count; |
---|
187 | code = code>>1; |
---|
188 | } |
---|
189 | } |
---|
190 | |
---|
191 | |
---|
192 | void GetParam(G4int & A, G4int & Z, G4double & pMin, G4double & pMax, |
---|
193 | G4int &nEvents, G4int & particleCode) |
---|
194 | { |
---|
195 | // A and Z |
---|
196 | cout << "Input Nucleus A and Z: "; |
---|
197 | cin >> A >> Z; |
---|
198 | // number of events |
---|
199 | cout << "Input number of events: "; |
---|
200 | cin >> nEvents; |
---|
201 | } |
---|
202 | |
---|
203 | |
---|
204 | |
---|
205 | |
---|
206 | G4ThreeVector GetSpherePoint(G4double r) |
---|
207 | { |
---|
208 | // random point INSIDE a sphere |
---|
209 | |
---|
210 | G4double x = r*(2*G4UniformRand()-1); |
---|
211 | G4double tmp = sqrt(r*r-x*x); |
---|
212 | G4double y = tmp*(2*G4UniformRand()-1); |
---|
213 | tmp = sqrt(r*r-x*x-y*y); |
---|
214 | G4double z = tmp*(2*G4UniformRand()-1); |
---|
215 | G4ThreeVector point; |
---|
216 | point.setX(x); |
---|
217 | point.setY(y); |
---|
218 | point.setZ(z); |
---|
219 | return point; |
---|
220 | } |
---|
221 | G4ThreeVector GetSurfacePoint(G4double r) |
---|
222 | { |
---|
223 | // random point oin surface of a sphere |
---|
224 | |
---|
225 | G4double x = r*(2*G4UniformRand()-1); |
---|
226 | G4double tmp = sqrt(r*r-x*x); |
---|
227 | G4double y = tmp*(2*G4UniformRand()-1); |
---|
228 | G4double z = -sqrt(r*r-x*x-y*y); |
---|
229 | G4ThreeVector point; |
---|
230 | point.setX(x); |
---|
231 | point.setY(y); |
---|
232 | point.setZ(z); |
---|
233 | return point; |
---|
234 | } |
---|