| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 #include "G4Absorber.hh" 26 #include "G4Absorber.hh"
27 #include "G4KineticTrack.hh" 27 #include "G4KineticTrack.hh"
28 #include "G4PionPlus.hh" 28 #include "G4PionPlus.hh"
29 #include "G4PionMinus.hh" 29 #include "G4PionMinus.hh"
30 #include "G4PionZero.hh" 30 #include "G4PionZero.hh"
31 #include "G4Proton.hh" 31 #include "G4Proton.hh"
32 #include "G4Neutron.hh" 32 #include "G4Neutron.hh"
33 33
34 #include "G4PhysicalConstants.hh" 34 #include "G4PhysicalConstants.hh"
35 #include "G4SystemOfUnits.hh" 35 #include "G4SystemOfUnits.hh"
36 #include "G4LorentzRotation.hh" 36 #include "G4LorentzRotation.hh"
37 37
38 G4Absorber::G4Absorber(G4double cutOnP) 38 G4Absorber::G4Absorber(G4double cutOnP)
39 { 39 {
40 theCutOnP = cutOnP; 40 theCutOnP = cutOnP;
41 theAbsorbers = new G4KineticTrackVector; 41 theAbsorbers = new G4KineticTrackVector;
42 theProducts = new G4KineticTrackVector; 42 theProducts = new G4KineticTrackVector;
43 } 43 }
44 44
45 45
46 G4Absorber::~G4Absorber() 46 G4Absorber::~G4Absorber()
47 { 47 {
48 delete theAbsorbers; 48 delete theAbsorbers;
49 delete theProducts; 49 delete theProducts;
50 } 50 }
51 51
52 52
53 bool G4Absorber::WillBeAbsorbed(const G4Kineti 53 bool G4Absorber::WillBeAbsorbed(const G4KineticTrack & kt)
54 { 54 {
55 // FixMe: actually only for pions 55 // FixMe: actually only for pions
56 // if(kt.Get4Momentum().vect().mag() < theCut 56 // if(kt.Get4Momentum().vect().mag() < theCutOnP)
57 // Cut on kinetic Energy... 57 // Cut on kinetic Energy...
58 if (kt.Get4Momentum().e() - kt.GetActualMass 58 if (kt.Get4Momentum().e() - kt.GetActualMass() < theCutOnP)
59 { 59 {
60 if(kt.GetDefinition() == G4PionPlus::Pio 60 if(kt.GetDefinition() == G4PionPlus::PionPlus() ||
61 kt.GetDefinition() == G4PionZero::PionZero( 61 kt.GetDefinition() == G4PionZero::PionZero() ||
62 kt.GetDefinition() == G4PionMinus::PionMinu 62 kt.GetDefinition() == G4PionMinus::PionMinus())
63 { 63 {
64 return true; 64 return true;
65 } 65 }
66 } 66 }
67 return false; 67 return false;
68 } 68 }
69 69
70 70
71 71
72 G4bool G4Absorber::Absorb(G4KineticTrack & kt, 72 G4bool G4Absorber::Absorb(G4KineticTrack & kt, G4KineticTrackVector & tgt)
73 { 73 {
74 if(!FindAbsorbers(kt, tgt)) 74 if(!FindAbsorbers(kt, tgt))
75 return false; 75 return false;
76 return FindProducts(kt); 76 return FindProducts(kt);
77 } 77 }
78 78
79 79
80 G4bool G4Absorber::FindAbsorbers(G4KineticTrac 80 G4bool G4Absorber::FindAbsorbers(G4KineticTrack & kt,
81 G4KineticTrackVector & tgt) 81 G4KineticTrackVector & tgt)
82 { 82 {
83 // Find a closest ( in space) pair of Nucleon 83 // Find a closest ( in space) pair of Nucleons capable to absorb pi+/pi-
84 // pi+ can be absorbed on np or nn resultin 84 // pi+ can be absorbed on np or nn resulting in pp or np
85 // pi- can be absorbed on np or pp resultin 85 // pi- can be absorbed on np or pp resulting in nn or np
86 86
87 // @GF: FindAbsorbers is unused, logic is seri 87 // @GF: FindAbsorbers is unused, logic is seriously wrong
88 88
89 G4KineticTrack * kt1 = NULL; 89 G4KineticTrack * kt1 = NULL;
90 G4KineticTrack * kt2 = NULL; 90 G4KineticTrack * kt2 = NULL;
91 G4double dist1 = DBL_MAX; // dist to close 91 G4double dist1 = DBL_MAX; // dist to closest nucleon
92 G4double dist2 = DBL_MAX; // dist to next 92 G4double dist2 = DBL_MAX; // dist to next close
93 G4double charge1 = 0; 93 G4double charge1 = 0;
94 // G4double charge2 = 0; // charge2 is only 94 // G4double charge2 = 0; // charge2 is only assigned to, never used
95 G4double charge0 = kt.GetDefinition()->GetPD 95 G4double charge0 = kt.GetDefinition()->GetPDGCharge();
96 G4ThreeVector pos = kt.GetPosition(); 96 G4ThreeVector pos = kt.GetPosition();
97 97
98 std::vector<G4KineticTrack *>::iterator iter 98 std::vector<G4KineticTrack *>::iterator iter;
99 for(iter = tgt.begin(); iter != tgt.end(); + 99 for(iter = tgt.begin(); iter != tgt.end(); ++iter)
100 { 100 {
101 G4KineticTrack * curr = *iter; 101 G4KineticTrack * curr = *iter;
102 G4double dist = (pos-curr->GetPosition()). 102 G4double dist = (pos-curr->GetPosition()).mag();
103 if(dist >= dist2) 103 if(dist >= dist2)
104 continue; 104 continue;
105 if(dist < dist1) 105 if(dist < dist1)
106 { 106 {
107 if(dist1 == DBL_MAX) // accept 1st as a 107 if(dist1 == DBL_MAX) // accept 1st as a candidate,
108 { 108 {
109 kt1 = curr; 109 kt1 = curr;
110 charge1 = kt1->GetDefinition()->GetPDGCharge 110 charge1 = kt1->GetDefinition()->GetPDGCharge();
111 dist1 = dist; 111 dist1 = dist;
112 continue; 112 continue;
113 } 113 }
114 if(dist2 == DBL_MAX) // accept the candi 114 if(dist2 == DBL_MAX) // accept the candidate and shift kt1 to kt2
115 { // @GF: should'nt we check if co 115 { // @GF: should'nt we check if compatible?
116 kt2 = kt1; 116 kt2 = kt1;
117 // charge2 = charge1; 117 // charge2 = charge1;
118 dist2 = dist1; 118 dist2 = dist1;
119 kt1 = curr; 119 kt1 = curr;
120 charge1 = kt1->GetDefinition()->GetPDGCharge 120 charge1 = kt1->GetDefinition()->GetPDGCharge();
121 dist1 = dist; 121 dist1 = dist;
122 continue; 122 continue;
123 } 123 }
124 // test the compatibility with charge conserva 124 // test the compatibility with charge conservation for new config
125 G4double charge = curr->GetDefinition()- 125 G4double charge = curr->GetDefinition()->GetPDGCharge();
126 if((charge0+charge1+charge < 0.) || //te 126 if((charge0+charge1+charge < 0.) || //test config (curr,kt1)
127 (charge0+charge1+charge) > 2*eplus) 127 (charge0+charge1+charge) > 2*eplus)
128 { // incompatible: change kt1 with curr 128 { // incompatible: change kt1 with curr.
129 kt1 = curr; 129 kt1 = curr;
130 charge1 = charge; 130 charge1 = charge;
131 dist1 = dist; 131 dist1 = dist;
132 } 132 }
133 else 133 else
134 { // compatible: change kt1 with curr an 134 { // compatible: change kt1 with curr and kt2 with kt1
135 kt2 = kt1; 135 kt2 = kt1;
136 // charge2 = charge1; 136 // charge2 = charge1;
137 dist2 = dist1; 137 dist2 = dist1;
138 kt1 = curr; 138 kt1 = curr;
139 charge1 = charge; 139 charge1 = charge;
140 dist1 = dist; 140 dist1 = dist;
141 } 141 }
142 continue; 142 continue;
143 } 143 }
144 // here if dist1 < dist < dist2 144 // here if dist1 < dist < dist2
145 if(dist2 == DBL_MAX) // accept the candida 145 if(dist2 == DBL_MAX) // accept the candidate
146 { 146 {
147 kt2 = curr; 147 kt2 = curr;
148 // charge2 = kt2->GetDefinition()->GetPDG 148 // charge2 = kt2->GetDefinition()->GetPDGCharge();
149 dist2 = dist; 149 dist2 = dist;
150 continue; 150 continue;
151 } 151 }
152 // test the compatibility with charge conserva 152 // test the compatibility with charge conservation
153 G4double charge = curr->GetDefinition()->G 153 G4double charge = curr->GetDefinition()->GetPDGCharge();
154 if((charge0+charge1+charge < 0.) || 154 if((charge0+charge1+charge < 0.) ||
155 (charge0+charge1+charge) > 2*eplus) 155 (charge0+charge1+charge) > 2*eplus)
156 continue; // incomatible: do nothing 156 continue; // incomatible: do nothing
157 // compatible: change kt2 with curr 157 // compatible: change kt2 with curr
158 kt2 = curr; 158 kt2 = curr;
159 // charge2 = charge; 159 // charge2 = charge;
160 dist2 = dist; 160 dist2 = dist;
161 } 161 }
162 162
163 theAbsorbers->clear(); // do not delete trac 163 theAbsorbers->clear(); // do not delete tracks in theAbsorbers vector!
164 if((kt1 == NULL) || (kt2 == NULL)) 164 if((kt1 == NULL) || (kt2 == NULL))
165 return false; 165 return false;
166 166
167 theAbsorbers->push_back(kt1); 167 theAbsorbers->push_back(kt1);
168 theAbsorbers->push_back(kt2); 168 theAbsorbers->push_back(kt2);
169 return true; 169 return true;
170 } 170 }
171 171
172 172
173 173
174 G4bool G4Absorber::FindProducts(G4KineticTrack 174 G4bool G4Absorber::FindProducts(G4KineticTrack & kt)
175 { 175 {
176 // Choose the products type 176 // Choose the products type
177 const G4ParticleDefinition * prod1; 177 const G4ParticleDefinition * prod1;
178 const G4ParticleDefinition * prod2; 178 const G4ParticleDefinition * prod2;
179 G4KineticTrack * abs1 = (*theAbsorbers)[0]; 179 G4KineticTrack * abs1 = (*theAbsorbers)[0];
180 G4KineticTrack * abs2 = (*theAbsorbers)[1]; 180 G4KineticTrack * abs2 = (*theAbsorbers)[1];
181 181
182 G4double charge = kt.GetDefinition()->GetPDG 182 G4double charge = kt.GetDefinition()->GetPDGCharge();
183 if(charge == eplus) 183 if(charge == eplus)
184 { // a neutron become proton 184 { // a neutron become proton
185 prod1 = G4Proton::Proton(); 185 prod1 = G4Proton::Proton();
186 if(abs1->GetDefinition() == G4Neutron::Neu 186 if(abs1->GetDefinition() == G4Neutron::Neutron())
187 prod2 = abs2->GetDefinition(); 187 prod2 = abs2->GetDefinition();
188 else 188 else
189 prod2 = G4Proton::Proton(); 189 prod2 = G4Proton::Proton();
190 } 190 }
191 else if(charge == -eplus) 191 else if(charge == -eplus)
192 { // a proton become neutron 192 { // a proton become neutron
193 prod1 = G4Neutron::Neutron(); 193 prod1 = G4Neutron::Neutron();
194 if(abs1->GetDefinition() == G4Proton::Prot 194 if(abs1->GetDefinition() == G4Proton::Proton())
195 prod2 = abs2->GetDefinition(); 195 prod2 = abs2->GetDefinition();
196 else 196 else
197 prod2 = G4Neutron::Neutron(); 197 prod2 = G4Neutron::Neutron();
198 } 198 }
199 else // charge = 0: leave particle types un 199 else // charge = 0: leave particle types unchenged
200 { 200 {
201 prod1 = abs1->GetDefinition(); 201 prod1 = abs1->GetDefinition();
202 prod2 = abs2->GetDefinition(); 202 prod2 = abs2->GetDefinition();
203 } 203 }
204 204
205 // Translate to the CMS frame 205 // Translate to the CMS frame
206 G4LorentzVector momLab = kt.Get4Momentum()+a 206 G4LorentzVector momLab = kt.Get4Momentum()+abs1->Get4Momentum()+
207 abs2->Get4Momentum(); 207 abs2->Get4Momentum();
208 G4LorentzRotation toCMSFrame((-1)*momLab.boo 208 G4LorentzRotation toCMSFrame((-1)*momLab.boostVector());
209 G4LorentzRotation toLabFrame(momLab.boostVec 209 G4LorentzRotation toLabFrame(momLab.boostVector());
210 G4LorentzVector momCMS = toCMSFrame*momLab; 210 G4LorentzVector momCMS = toCMSFrame*momLab;
211 211
212 // Evaluate the final momentum of products 212 // Evaluate the final momentum of products
213 G4double ms1 = prod1->GetPDGMass(); 213 G4double ms1 = prod1->GetPDGMass();
214 G4double ms2 = prod2->GetPDGMass(); 214 G4double ms2 = prod2->GetPDGMass();
215 G4double e0 = momCMS.e(); 215 G4double e0 = momCMS.e();
216 G4double squareP = (e0*e0*e0*e0-2*e0*e0*(ms1 216 G4double squareP = (e0*e0*e0*e0-2*e0*e0*(ms1*ms1+ms2*ms2)+
217 (ms2*ms2-ms1*ms1)*(ms2*ms2-ms1*ms1))/(4*e0 217 (ms2*ms2-ms1*ms1)*(ms2*ms2-ms1*ms1))/(4*e0*e0);
218 // if(squareP < 0) // should never happen 218 // if(squareP < 0) // should never happen
219 // squareP = 0; 219 // squareP = 0;
220 G4ThreeVector mom1CMS = GetRandomDirection() 220 G4ThreeVector mom1CMS = GetRandomDirection();
221 mom1CMS = std::sqrt(squareP)*mom1CMS; 221 mom1CMS = std::sqrt(squareP)*mom1CMS;
222 G4LorentzVector final4Mom1CMS(mom1CMS, std:: 222 G4LorentzVector final4Mom1CMS(mom1CMS, std::sqrt(squareP+ms1*ms1));
223 G4LorentzVector final4Mom2CMS((-1)*mom1CMS, 223 G4LorentzVector final4Mom2CMS((-1)*mom1CMS, std::sqrt(squareP+ms2*ms2));
224 224
225 // Go back to the lab frame 225 // Go back to the lab frame
226 G4LorentzVector mom1 = toLabFrame*final4Mom1 226 G4LorentzVector mom1 = toLabFrame*final4Mom1CMS;
227 G4LorentzVector mom2 = toLabFrame*final4Mom2 227 G4LorentzVector mom2 = toLabFrame*final4Mom2CMS;
228 228
229 // ------ debug 229 // ------ debug
230 /* 230 /*
231 G4LorentzVector temp = mom1+mom2; 231 G4LorentzVector temp = mom1+mom2;
232 232
233 cout << (1/MeV)*momLab.x() << " " << (1/MeV) 233 cout << (1/MeV)*momLab.x() << " " << (1/MeV)*momLab.y() << " "
234 << (1/MeV)*momLab.z() << " " << (1/MeV) 234 << (1/MeV)*momLab.z() << " " << (1/MeV)*momLab.t() << " "
235 << (1/MeV)*momLab.vect().mag() << " " < 235 << (1/MeV)*momLab.vect().mag() << " " << (1/MeV)*momLab.mag() << " "
236 << (1/MeV)*temp.x() << " " << (1/MeV)*t 236 << (1/MeV)*temp.x() << " " << (1/MeV)*temp.y() << " "
237 << (1/MeV)*temp.z() << " " << (1/MeV)*t 237 << (1/MeV)*temp.z() << " " << (1/MeV)*temp.t() << " "
238 << (1/MeV)*temp.vect().mag() << " " << 238 << (1/MeV)*temp.vect().mag() << " " << (1/MeV)*temp.mag() << " "
239 << (1/MeV)*std::sqrt(squareP) << endl; 239 << (1/MeV)*std::sqrt(squareP) << endl;
240 240
241 */ 241 */
242 // ------ end debug 242 // ------ end debug
243 243
244 // Build two new kinetic tracks and add to pro 244 // Build two new kinetic tracks and add to products
245 G4KineticTrack * kt1 = new G4KineticTrack(pr 245 G4KineticTrack * kt1 = new G4KineticTrack(prod1, 0., abs1->GetPosition(),
246 mom1); 246 mom1);
247 G4KineticTrack * kt2 = new G4KineticTrack(pr 247 G4KineticTrack * kt2 = new G4KineticTrack(prod2, 0., abs2->GetPosition(),
248 mom2); 248 mom2);
249 // ------ debug 249 // ------ debug
250 /* 250 /*
251 G4LorentzVector initialMom1 = abs1->Get4Mome 251 G4LorentzVector initialMom1 = abs1->Get4Momentum();
252 G4LorentzVector initialMom2 = abs2->Get4Mome 252 G4LorentzVector initialMom2 = abs2->Get4Momentum();
253 G4LorentzVector pion4MomCMS = toCMSFrame*kt. 253 G4LorentzVector pion4MomCMS = toCMSFrame*kt.Get4Momentum();
254 cout << (1/MeV)*initialMom1.x() << " " << (1 254 cout << (1/MeV)*initialMom1.x() << " " << (1/MeV)*initialMom1.y() << " "
255 << (1/MeV)*initialMom1.z() << " " << (1 255 << (1/MeV)*initialMom1.z() << " " << (1/MeV)*initialMom1.e() << " "
256 << (1/MeV)*initialMom1.vect().mag() << 256 << (1/MeV)*initialMom1.vect().mag() << " "
257 << (1/MeV)*initialMom2.x() << " " << (1 257 << (1/MeV)*initialMom2.x() << " " << (1/MeV)*initialMom2.y() << " "
258 << (1/MeV)*initialMom2.z() << " " << (1 258 << (1/MeV)*initialMom2.z() << " " << (1/MeV)*initialMom2.e() << " "
259 << (1/MeV)*initialMom2.vect().mag() << 259 << (1/MeV)*initialMom2.vect().mag() << " "
260 << (1/MeV)*mom1.x() << " " << (1/MeV)*m 260 << (1/MeV)*mom1.x() << " " << (1/MeV)*mom1.y() << " "
261 << (1/MeV)*mom1.z() << " " << (1/MeV)*m 261 << (1/MeV)*mom1.z() << " " << (1/MeV)*mom1.e() << " "
262 << (1/MeV)*mom1.vect().mag() << " " 262 << (1/MeV)*mom1.vect().mag() << " "
263 << (1/MeV)*mom2.x() << " " << (1/MeV)*m 263 << (1/MeV)*mom2.x() << " " << (1/MeV)*mom2.y() << " "
264 << (1/MeV)*mom2.z() << " " << (1/MeV)*m 264 << (1/MeV)*mom2.z() << " " << (1/MeV)*mom2.e() << " "
265 << (1/MeV)*mom2.vect().mag() << " " 265 << (1/MeV)*mom2.vect().mag() << " "
266 << (1/MeV)*pion4MomCMS.x() << " " << (1 266 << (1/MeV)*pion4MomCMS.x() << " " << (1/MeV)*pion4MomCMS.y() << " "
267 << (1/MeV)*pion4MomCMS.z() << " " << (1 267 << (1/MeV)*pion4MomCMS.z() << " " << (1/MeV)*pion4MomCMS.e() << " "
268 << (1/MeV)*pion4MomCMS.vect().mag() << 268 << (1/MeV)*pion4MomCMS.vect().mag() << " "
269 << (1/MeV)*final4Mom1CMS.x() << " " << 269 << (1/MeV)*final4Mom1CMS.x() << " " << (1/MeV)*final4Mom1CMS.y() << " "
270 << (1/MeV)*final4Mom1CMS.z() << " " << 270 << (1/MeV)*final4Mom1CMS.z() << " " << (1/MeV)*final4Mom1CMS.e() << " "
271 << (1/MeV)*final4Mom1CMS.vect().mag() < 271 << (1/MeV)*final4Mom1CMS.vect().mag() << " "
272 << (1/MeV)*final4Mom2CMS.x() << " " << 272 << (1/MeV)*final4Mom2CMS.x() << " " << (1/MeV)*final4Mom2CMS.y() << " "
273 << (1/MeV)*final4Mom2CMS.z() << " " << 273 << (1/MeV)*final4Mom2CMS.z() << " " << (1/MeV)*final4Mom2CMS.e() << " "
274 << (1/MeV)*final4Mom2CMS.vect().mag() < 274 << (1/MeV)*final4Mom2CMS.vect().mag() << endl;
275 */ 275 */
276 // ------ end debug 276 // ------ end debug
277 277
278 theProducts->clear(); 278 theProducts->clear();
279 theProducts->push_back(kt1); 279 theProducts->push_back(kt1);
280 theProducts->push_back(kt2); 280 theProducts->push_back(kt2);
281 return true; 281 return true;
282 } 282 }
283 283
284 284
285 285
286 G4ThreeVector G4Absorber::GetRandomDirection() 286 G4ThreeVector G4Absorber::GetRandomDirection()
287 { 287 {
288 G4double theta = 2.0*G4UniformRand()-1.0; 288 G4double theta = 2.0*G4UniformRand()-1.0;
289 theta = std::acos(theta); 289 theta = std::acos(theta);
290 G4double phi = G4UniformRand()*2*pi; 290 G4double phi = G4UniformRand()*2*pi;
291 G4ThreeVector direction(std::sin(theta)*std: 291 G4ThreeVector direction(std::sin(theta)*std::cos(phi), std::sin(theta)*std::sin(phi), std::cos(theta));
292 return direction; 292 return direction;
293 } 293 }
294 294
295 295
296 296
297 297
298 298
299 299
300 300
301 301