Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
G4Absorber Class Reference

#include <G4Absorber.hh>

Public Member Functions

 G4Absorber (G4double cutOnP)
 
 ~G4Absorber ()
 
G4bool WillBeAbsorbed (const G4KineticTrack &kt)
 
G4bool Absorb (G4KineticTrack &kt, G4KineticTrackVector &tgt)
 
G4KineticTrackVector * GetAbsorbers ()
 
G4KineticTrackVector * GetProducts ()
 
G4bool FindAbsorbers (G4KineticTrack &kt, G4KineticTrackVector &tgt)
 
G4bool FindProducts (G4KineticTrack &kt)
 

Detailed Description

Definition at line 34 of file G4Absorber.hh.

Constructor & Destructor Documentation

◆ G4Absorber()

G4Absorber::G4Absorber ( G4double  cutOnP)

Definition at line 38 of file G4Absorber.cc.

39{
40 theCutOnP = cutOnP;
41 theAbsorbers = new G4KineticTrackVector;
42 theProducts = new G4KineticTrackVector;
43}

◆ ~G4Absorber()

G4Absorber::~G4Absorber ( )

Definition at line 46 of file G4Absorber.cc.

47{
48 delete theAbsorbers;
49 delete theProducts;
50}

Member Function Documentation

◆ Absorb()

G4bool G4Absorber::Absorb ( G4KineticTrack &  kt,
G4KineticTrackVector &  tgt 
)

Definition at line 72 of file G4Absorber.cc.

73{
74 if(!FindAbsorbers(kt, tgt))
75 return false;
76 return FindProducts(kt);
77}
G4Absorber::FindAbsorbers
G4bool FindAbsorbers(G4KineticTrack &kt, G4KineticTrackVector &tgt)
Definition: G4Absorber.cc:80
G4Absorber::FindProducts
G4bool FindProducts(G4KineticTrack &kt)
Definition: G4Absorber.cc:174

◆ FindAbsorbers()

G4bool G4Absorber::FindAbsorbers ( G4KineticTrack &  kt,
G4KineticTrackVector &  tgt 
)

Definition at line 80 of file G4Absorber.cc.

82{
83// Find a closest ( in space) pair of Nucleons capable to absorb pi+/pi-
84// pi+ can be absorbed on np or nn resulting in pp or np
85// pi- can be absorbed on np or pp resulting in nn or np
86
87// @GF: FindAbsorbers is unused, logic is seriously wrong
88
89 G4KineticTrack * kt1 = NULL;
90 G4KineticTrack * kt2 = NULL;
91 G4double dist1 = DBL_MAX; // dist to closest nucleon
92 G4double dist2 = DBL_MAX; // dist to next close
93 G4double charge1 = 0;
94// G4double charge2 = 0; // charge2 is only assigned to, never used
95 G4double charge0 = kt.GetDefinition()->GetPDGCharge();
96 G4ThreeVector pos = kt.GetPosition();
97
98 std::vector<G4KineticTrack *>::iterator iter;
99 for(iter = tgt.begin(); iter != tgt.end(); ++iter)
100 {
101 G4KineticTrack * curr = *iter;
102 G4double dist = (pos-curr->GetPosition()).mag();
103 if(dist >= dist2)
104 continue;
105 if(dist < dist1)
106 {
107 if(dist1 == DBL_MAX) // accept 1st as a candidate,
108 {
109 kt1 = curr;
110 charge1 = kt1->GetDefinition()->GetPDGCharge();
111 dist1 = dist;
112 continue;
113 }
114 if(dist2 == DBL_MAX) // accept the candidate and shift kt1 to kt2
115 { // @GF: should'nt we check if compatible?
116 kt2 = kt1;
117// charge2 = charge1;
118 dist2 = dist1;
119 kt1 = curr;
120 charge1 = kt1->GetDefinition()->GetPDGCharge();
121 dist1 = dist;
122 continue;
123 }
124// test the compatibility with charge conservation for new config
125 G4double charge = curr->GetDefinition()->GetPDGCharge();
126 if((charge0+charge1+charge < 0.) || //test config (curr,kt1)
127 (charge0+charge1+charge) > 2*eplus)
128 { // incompatible: change kt1 with curr.
129 kt1 = curr;
130 charge1 = charge;
131 dist1 = dist;
132 }
133 else
134 { // compatible: change kt1 with curr and kt2 with kt1
135 kt2 = kt1;
136// charge2 = charge1;
137 dist2 = dist1;
138 kt1 = curr;
139 charge1 = charge;
140 dist1 = dist;
141 }
142 continue;
143 }
144// here if dist1 < dist < dist2
145 if(dist2 == DBL_MAX) // accept the candidate
146 {
147 kt2 = curr;
148// charge2 = kt2->GetDefinition()->GetPDGCharge();
149 dist2 = dist;
150 continue;
151 }
152// test the compatibility with charge conservation
153 G4double charge = curr->GetDefinition()->GetPDGCharge();
154 if((charge0+charge1+charge < 0.) ||
155 (charge0+charge1+charge) > 2*eplus)
156 continue; // incomatible: do nothing
157// compatible: change kt2 with curr
158 kt2 = curr;
159// charge2 = charge;
160 dist2 = dist;
161 }
162
163 theAbsorbers->clear(); // do not delete tracks in theAbsorbers vector!
164 if((kt1 == NULL) || (kt2 == NULL))
165 return false;
166
167 theAbsorbers->push_back(kt1);
168 theAbsorbers->push_back(kt2);
169 return true;
170}
G4double
double G4double
Definition: G4Types.hh:83
G4KineticTrack::GetPosition
const G4ThreeVector & GetPosition() const
Definition: G4KineticTrack.hh:229
G4KineticTrack::GetDefinition
const G4ParticleDefinition * GetDefinition() const
Definition: G4KineticTrack.hh:209
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:112
DBL_MAX
#define DBL_MAX
Definition: templates.hh:62

Referenced by Absorb().

◆ FindProducts()

G4bool G4Absorber::FindProducts ( G4KineticTrack &  kt)

Definition at line 174 of file G4Absorber.cc.

175{
176// Choose the products type
177 const G4ParticleDefinition * prod1;
178 const G4ParticleDefinition * prod2;
179 G4KineticTrack * abs1 = (*theAbsorbers)[0];
180 G4KineticTrack * abs2 = (*theAbsorbers)[1];
181
182 G4double charge = kt.GetDefinition()->GetPDGCharge();
183 if(charge == eplus)
184 { // a neutron become proton
185 prod1 = G4Proton::Proton();
186 if(abs1->GetDefinition() == G4Neutron::Neutron())
187 prod2 = abs2->GetDefinition();
188 else
189 prod2 = G4Proton::Proton();
190 }
191 else if(charge == -eplus)
192 { // a proton become neutron
193 prod1 = G4Neutron::Neutron();
194 if(abs1->GetDefinition() == G4Proton::Proton())
195 prod2 = abs2->GetDefinition();
196 else
197 prod2 = G4Neutron::Neutron();
198 }
199 else // charge = 0: leave particle types unchenged
200 {
201 prod1 = abs1->GetDefinition();
202 prod2 = abs2->GetDefinition();
203 }
204
205// Translate to the CMS frame
206 G4LorentzVector momLab = kt.Get4Momentum()+abs1->Get4Momentum()+
207 abs2->Get4Momentum();
208 G4LorentzRotation toCMSFrame((-1)*momLab.boostVector());
209 G4LorentzRotation toLabFrame(momLab.boostVector());
210 G4LorentzVector momCMS = toCMSFrame*momLab;
211
212// Evaluate the final momentum of products
213 G4double ms1 = prod1->GetPDGMass();
214 G4double ms2 = prod2->GetPDGMass();
215 G4double e0 = momCMS.e();
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*e0);
218// if(squareP < 0) // should never happen
219// squareP = 0;
220 G4ThreeVector mom1CMS = GetRandomDirection();
221 mom1CMS = std::sqrt(squareP)*mom1CMS;
222 G4LorentzVector final4Mom1CMS(mom1CMS, std::sqrt(squareP+ms1*ms1));
223 G4LorentzVector final4Mom2CMS((-1)*mom1CMS, std::sqrt(squareP+ms2*ms2));
224
225// Go back to the lab frame
226 G4LorentzVector mom1 = toLabFrame*final4Mom1CMS;
227 G4LorentzVector mom2 = toLabFrame*final4Mom2CMS;
228
229// ------ debug
230/*
231 G4LorentzVector temp = mom1+mom2;
232
233 cout << (1/MeV)*momLab.x() << " " << (1/MeV)*momLab.y() << " "
234 << (1/MeV)*momLab.z() << " " << (1/MeV)*momLab.t() << " "
235 << (1/MeV)*momLab.vect().mag() << " " << (1/MeV)*momLab.mag() << " "
236 << (1/MeV)*temp.x() << " " << (1/MeV)*temp.y() << " "
237 << (1/MeV)*temp.z() << " " << (1/MeV)*temp.t() << " "
238 << (1/MeV)*temp.vect().mag() << " " << (1/MeV)*temp.mag() << " "
239 << (1/MeV)*std::sqrt(squareP) << endl;
240
241*/
242// ------ end debug
243
244// Build two new kinetic tracks and add to products
245 G4KineticTrack * kt1 = new G4KineticTrack(prod1, 0., abs1->GetPosition(),
246 mom1);
247 G4KineticTrack * kt2 = new G4KineticTrack(prod2, 0., abs2->GetPosition(),
248 mom2);
249// ------ debug
250/*
251 G4LorentzVector initialMom1 = abs1->Get4Momentum();
252 G4LorentzVector initialMom2 = abs2->Get4Momentum();
253 G4LorentzVector pion4MomCMS = toCMSFrame*kt.Get4Momentum();
254 cout << (1/MeV)*initialMom1.x() << " " << (1/MeV)*initialMom1.y() << " "
255 << (1/MeV)*initialMom1.z() << " " << (1/MeV)*initialMom1.e() << " "
256 << (1/MeV)*initialMom1.vect().mag() << " "
257 << (1/MeV)*initialMom2.x() << " " << (1/MeV)*initialMom2.y() << " "
258 << (1/MeV)*initialMom2.z() << " " << (1/MeV)*initialMom2.e() << " "
259 << (1/MeV)*initialMom2.vect().mag() << " "
260 << (1/MeV)*mom1.x() << " " << (1/MeV)*mom1.y() << " "
261 << (1/MeV)*mom1.z() << " " << (1/MeV)*mom1.e() << " "
262 << (1/MeV)*mom1.vect().mag() << " "
263 << (1/MeV)*mom2.x() << " " << (1/MeV)*mom2.y() << " "
264 << (1/MeV)*mom2.z() << " " << (1/MeV)*mom2.e() << " "
265 << (1/MeV)*mom2.vect().mag() << " "
266 << (1/MeV)*pion4MomCMS.x() << " " << (1/MeV)*pion4MomCMS.y() << " "
267 << (1/MeV)*pion4MomCMS.z() << " " << (1/MeV)*pion4MomCMS.e() << " "
268 << (1/MeV)*pion4MomCMS.vect().mag() << " "
269 << (1/MeV)*final4Mom1CMS.x() << " " << (1/MeV)*final4Mom1CMS.y() << " "
270 << (1/MeV)*final4Mom1CMS.z() << " " << (1/MeV)*final4Mom1CMS.e() << " "
271 << (1/MeV)*final4Mom1CMS.vect().mag() << " "
272 << (1/MeV)*final4Mom2CMS.x() << " " << (1/MeV)*final4Mom2CMS.y() << " "
273 << (1/MeV)*final4Mom2CMS.z() << " " << (1/MeV)*final4Mom2CMS.e() << " "
274 << (1/MeV)*final4Mom2CMS.vect().mag() << endl;
275*/
276// ------ end debug
277
278 theProducts->clear();
279 theProducts->push_back(kt1);
280 theProducts->push_back(kt2);
281 return true;
282}
CLHEP::HepLorentzVector::boostVector
Hep3Vector boostVector() const
Definition: LorentzVector.cc:173
G4KineticTrack::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition: G4KineticTrack.hh:239
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:103
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:92

Referenced by Absorb().

◆ GetAbsorbers()

G4KineticTrackVector * G4Absorber::GetAbsorbers ( )
inline

Definition at line 66 of file G4Absorber.hh.

67{
68 return theAbsorbers;
69}

◆ GetProducts()

G4KineticTrackVector * G4Absorber::GetProducts ( )
inline

Definition at line 71 of file G4Absorber.hh.

72{
73 return theProducts;
74}

◆ WillBeAbsorbed()

bool G4Absorber::WillBeAbsorbed ( const G4KineticTrack &  kt)

Definition at line 53 of file G4Absorber.cc.

54{
55 // FixMe: actually only for pions
56// if(kt.Get4Momentum().vect().mag() < theCutOnP)
57// Cut on kinetic Energy...
58 if (kt.Get4Momentum().e() - kt.GetActualMass() < theCutOnP)
59 {
60 if(kt.GetDefinition() == G4PionPlus::PionPlus() ||
61 kt.GetDefinition() == G4PionZero::PionZero() ||
62 kt.GetDefinition() == G4PionMinus::PionMinus())
63 {
64 return true;
65 }
66 }
67 return false;
68}
G4KineticTrack::GetActualMass
G4double GetActualMass() const
Definition: G4KineticTrack.hh:318
G4PionMinus::PionMinus
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:97
G4PionPlus::PionPlus
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:97
G4PionZero::PionZero
static G4PionZero * PionZero()
Definition: G4PionZero.cc:107
The documentation for this class was generated from the following files: