source: trunk/source/processes/hadronic/models/binary_cascade/src/G4Absorber.cc@ 1337

Last change on this file since 1337 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 9.3 KB
Line 
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#include "G4Absorber.hh"
27#include "G4KineticTrack.hh"
28#include "G4PionPlus.hh"
29#include "G4PionMinus.hh"
30#include "G4PionZero.hh"
31#include "G4Proton.hh"
32#include "G4Neutron.hh"
33
34#include "G4LorentzRotation.hh"
35
36G4Absorber::G4Absorber(G4double cutOnP)
37{
38 theCutOnP = cutOnP;
39 theAbsorbers = new G4KineticTrackVector;
40 theProducts = new G4KineticTrackVector;
41}
42
43
44G4Absorber::~G4Absorber()
45{
46 delete theAbsorbers;
47 delete theProducts;
48}
49
50
51bool G4Absorber::WillBeAbsorbed(const G4KineticTrack & kt)
52{
53 // FixMe: actually only for pions
54// if(kt.Get4Momentum().vect().mag() < theCutOnP)
55// Cut on kinetic Energy...
56 if (kt.Get4Momentum().e() - kt.GetActualMass() < theCutOnP)
57 {
58 if(kt.GetDefinition() == G4PionPlus::PionPlus() ||
59 kt.GetDefinition() == G4PionZero::PionZero() ||
60 kt.GetDefinition() == G4PionMinus::PionMinus())
61 {
62 return true;
63 }
64 }
65 return false;
66}
67
68
69
70G4bool G4Absorber::Absorb(G4KineticTrack & kt, G4KineticTrackVector & tgt)
71{
72 if(!FindAbsorbers(kt, tgt))
73 return false;
74 return FindProducts(kt);
75}
76
77
78G4bool G4Absorber::FindAbsorbers(G4KineticTrack & kt,
79 G4KineticTrackVector & tgt)
80{
81 G4KineticTrack * kt1 = NULL;
82 G4KineticTrack * kt2 = NULL;
83 G4double dist1 = DBL_MAX;
84 G4double dist2 = DBL_MAX;
85 G4double charge1 = 0;
86 G4double charge2 = 0;
87 G4double charge0 = kt.GetDefinition()->GetPDGCharge();
88 G4ThreeVector pos = kt.GetPosition();
89
90 std::vector<G4KineticTrack *>::iterator iter;
91 for(iter = tgt.begin(); iter != tgt.end(); ++iter)
92 {
93 G4KineticTrack * curr = *iter;
94 G4double dist = (pos-curr->GetPosition()).mag();
95 if(dist >= dist2)
96 continue;
97 if(dist < dist1)
98 {
99 if(dist1 == DBL_MAX) // accept the candidate
100 {
101 kt1 = curr;
102 charge1 = kt1->GetDefinition()->GetPDGCharge();
103 dist1 = dist;
104 continue;
105 }
106 if(dist2 == DBL_MAX) // accept the candidate put kt1 in kt2
107 {
108 kt2 = kt1;
109 charge2 = charge1;
110 dist2 = dist1;
111 kt1 = curr;
112 charge1 = kt1->GetDefinition()->GetPDGCharge();
113 dist1 = dist;
114 continue;
115 }
116// test the compatibility with charge conservation
117 G4double charge = curr->GetDefinition()->GetPDGCharge();
118 if((charge0+charge1+charge < 0.) ||
119 (charge0+charge1+charge) > 2*eplus)
120 { // incomatible: change kt1 with curr.
121 kt1 = curr;
122 charge1 = charge;
123 dist1 = dist;
124 }
125 else
126 { // compatible: change kt1 with curr and kt2 with kt1
127 kt2 = kt1;
128 charge2 = charge1;
129 dist2 = dist1;
130 kt1 = curr;
131 charge1 = charge;
132 dist1 = dist;
133 }
134 continue;
135 }
136// here if dist1 < dist < dist2
137 if(dist2 == DBL_MAX) // accept the candidate
138 {
139 kt2 = curr;
140 charge2 = kt2->GetDefinition()->GetPDGCharge();
141 dist2 = dist;
142 continue;
143 }
144// test the compatibility with charge conservation
145 G4double charge = curr->GetDefinition()->GetPDGCharge();
146 if((charge0+charge1+charge < 0.) ||
147 (charge0+charge1+charge) > 2*eplus)
148 continue; // incomatible: do nothing
149// compatible: change kt2 with curr
150 kt2 = curr;
151 charge2 = charge;
152 dist2 = dist;
153 }
154
155 theAbsorbers->clear(); // do not delete tracks in theAbsorbers vector!
156 if((kt1 == NULL) || (kt2 == NULL))
157 return false;
158
159 theAbsorbers->push_back(kt1);
160 theAbsorbers->push_back(kt2);
161 return true;
162}
163
164
165
166G4bool G4Absorber::FindProducts(G4KineticTrack & kt)
167{
168// Choose the products type
169 G4ParticleDefinition * prod1;
170 G4ParticleDefinition * prod2;
171 G4KineticTrack * abs1 = (*theAbsorbers)[0];
172 G4KineticTrack * abs2 = (*theAbsorbers)[1];
173
174 G4double charge = kt.GetDefinition()->GetPDGCharge();
175 if(charge == eplus)
176 { // a neutron become proton
177 prod1 = G4Proton::Proton();
178 if(abs1->GetDefinition() == G4Neutron::Neutron())
179 prod2 = abs2->GetDefinition();
180 else
181 prod2 = G4Proton::Proton();
182 }
183 else if(charge == -eplus)
184 { // a proton become neutron
185 prod1 = G4Neutron::Neutron();
186 if(abs1->GetDefinition() == G4Proton::Proton())
187 prod2 = abs2->GetDefinition();
188 else
189 prod2 = G4Neutron::Neutron();
190 }
191 else // charge = 0: leave particle types unchenged
192 {
193 prod1 = abs1->GetDefinition();
194 prod2 = abs2->GetDefinition();
195 }
196
197// Translate to the CMS frame
198 G4LorentzVector momLab = kt.Get4Momentum()+abs1->Get4Momentum()+
199 abs2->Get4Momentum();
200 G4LorentzRotation toCMSFrame((-1)*momLab.boostVector());
201 G4LorentzRotation toLabFrame(momLab.boostVector());
202 G4LorentzVector momCMS = toCMSFrame*momLab;
203
204// Evaluate the final momentum of products
205 G4double m1 = prod1->GetPDGMass();
206 G4double m2 = prod2->GetPDGMass();
207 G4double e0 = momCMS.e();
208 G4double squareP = (e0*e0*e0*e0-2*e0*e0*(m1*m1+m2*m2)+
209 (m2*m2-m1*m1)*(m2*m2-m1*m1))/(4*e0*e0);
210// if(squareP < 0) // should never happen
211// squareP = 0;
212 G4ThreeVector mom1CMS = GetRandomDirection();
213 mom1CMS = std::sqrt(squareP)*mom1CMS;
214 G4LorentzVector final4Mom1CMS(mom1CMS, std::sqrt(squareP+m1*m1));
215 G4LorentzVector final4Mom2CMS((-1)*mom1CMS, std::sqrt(squareP+m2*m2));
216
217// Go back to the lab frame
218 G4LorentzVector mom1 = toLabFrame*final4Mom1CMS;
219 G4LorentzVector mom2 = toLabFrame*final4Mom2CMS;
220
221// ------ debug
222/*
223 G4LorentzVector temp = mom1+mom2;
224
225 cout << (1/MeV)*momLab.x() << " " << (1/MeV)*momLab.y() << " "
226 << (1/MeV)*momLab.z() << " " << (1/MeV)*momLab.t() << " "
227 << (1/MeV)*momLab.vect().mag() << " " << (1/MeV)*momLab.mag() << " "
228 << (1/MeV)*temp.x() << " " << (1/MeV)*temp.y() << " "
229 << (1/MeV)*temp.z() << " " << (1/MeV)*temp.t() << " "
230 << (1/MeV)*temp.vect().mag() << " " << (1/MeV)*temp.mag() << " "
231 << (1/MeV)*std::sqrt(squareP) << endl;
232
233*/
234// ------ end debug
235
236// Build two new kinetic tracks and add to products
237 G4KineticTrack * kt1 = new G4KineticTrack(prod1, 0., abs1->GetPosition(),
238 mom1);
239 G4KineticTrack * kt2 = new G4KineticTrack(prod2, 0., abs2->GetPosition(),
240 mom2);
241// ------ debug
242/*
243 G4LorentzVector initialMom1 = abs1->Get4Momentum();
244 G4LorentzVector initialMom2 = abs2->Get4Momentum();
245 G4LorentzVector pion4MomCMS = toCMSFrame*kt.Get4Momentum();
246 cout << (1/MeV)*initialMom1.x() << " " << (1/MeV)*initialMom1.y() << " "
247 << (1/MeV)*initialMom1.z() << " " << (1/MeV)*initialMom1.e() << " "
248 << (1/MeV)*initialMom1.vect().mag() << " "
249 << (1/MeV)*initialMom2.x() << " " << (1/MeV)*initialMom2.y() << " "
250 << (1/MeV)*initialMom2.z() << " " << (1/MeV)*initialMom2.e() << " "
251 << (1/MeV)*initialMom2.vect().mag() << " "
252 << (1/MeV)*mom1.x() << " " << (1/MeV)*mom1.y() << " "
253 << (1/MeV)*mom1.z() << " " << (1/MeV)*mom1.e() << " "
254 << (1/MeV)*mom1.vect().mag() << " "
255 << (1/MeV)*mom2.x() << " " << (1/MeV)*mom2.y() << " "
256 << (1/MeV)*mom2.z() << " " << (1/MeV)*mom2.e() << " "
257 << (1/MeV)*mom2.vect().mag() << " "
258 << (1/MeV)*pion4MomCMS.x() << " " << (1/MeV)*pion4MomCMS.y() << " "
259 << (1/MeV)*pion4MomCMS.z() << " " << (1/MeV)*pion4MomCMS.e() << " "
260 << (1/MeV)*pion4MomCMS.vect().mag() << " "
261 << (1/MeV)*final4Mom1CMS.x() << " " << (1/MeV)*final4Mom1CMS.y() << " "
262 << (1/MeV)*final4Mom1CMS.z() << " " << (1/MeV)*final4Mom1CMS.e() << " "
263 << (1/MeV)*final4Mom1CMS.vect().mag() << " "
264 << (1/MeV)*final4Mom2CMS.x() << " " << (1/MeV)*final4Mom2CMS.y() << " "
265 << (1/MeV)*final4Mom2CMS.z() << " " << (1/MeV)*final4Mom2CMS.e() << " "
266 << (1/MeV)*final4Mom2CMS.vect().mag() << endl;
267*/
268// ------ end debug
269
270 theProducts->clear();
271 theProducts->push_back(kt1);
272 theProducts->push_back(kt2);
273 return true;
274}
275
276
277
278G4ThreeVector G4Absorber::GetRandomDirection()
279{
280 G4double theta = 2.0*G4UniformRand()-1.0;
281 theta = std::acos(theta);
282 G4double phi = G4UniformRand()*2*pi;
283 G4ThreeVector direction(std::sin(theta)*std::cos(phi), std::sin(theta)*std::sin(phi), std::cos(theta));
284 return direction;
285}
286
287
288
289
290
291
292
Note: See TracBrowser for help on using the repository browser.