Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4Absorber Class Reference

#include <G4Absorber.hh>

Public Member Functions

 G4Absorber (G4double cutOnP)
 
 ~G4Absorber ()
 
G4bool WillBeAbsorbed (const G4KineticTrack &kt)
 
G4bool Absorb (G4KineticTrack &kt, G4KineticTrackVector &tgt)
 
G4KineticTrackVectorGetAbsorbers ()
 
G4KineticTrackVectorGetProducts ()
 
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 ( 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 ( )

Definition at line 46 of file G4Absorber.cc.

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

Member Function Documentation

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 }
G4bool FindAbsorbers(G4KineticTrack &kt, G4KineticTrackVector &tgt)
Definition: G4Absorber.cc:80
G4bool FindProducts(G4KineticTrack &kt)
Definition: G4Absorber.cc:174

Here is the call graph for this function:

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();
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 }
const G4ThreeVector & GetPosition() const
static constexpr double eplus
Definition: G4SIunits.hh:199
double G4double
Definition: G4Types.hh:76
G4double GetPDGCharge() const
#define DBL_MAX
Definition: templates.hh:83
const G4ParticleDefinition * GetDefinition() const
static const G4double pos

Here is the call graph for this function:

Here is the caller graph for this function:

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 }
Hep3Vector boostVector() const
const G4ThreeVector & GetPosition() const
static constexpr double eplus
Definition: G4SIunits.hh:199
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4double GetPDGMass() const
double G4double
Definition: G4Types.hh:76
const G4LorentzVector & Get4Momentum() const
G4double GetPDGCharge() const
const G4ParticleDefinition * GetDefinition() const

Here is the call graph for this function:

Here is the caller graph for this function:

G4KineticTrackVector * G4Absorber::GetAbsorbers ( )
inline

Definition at line 66 of file G4Absorber.hh.

67 {
68  return theAbsorbers;
69 }
G4KineticTrackVector * G4Absorber::GetProducts ( )
inline

Definition at line 71 of file G4Absorber.hh.

72 {
73  return theProducts;
74 }
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() ||
63  {
64  return true;
65  }
66  }
67  return false;
68 }
G4double GetActualMass() const
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
static G4PionZero * PionZero()
Definition: G4PionZero.cc:108
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
const G4LorentzVector & Get4Momentum() const
const G4ParticleDefinition * GetDefinition() const

Here is the call graph for this function:


The documentation for this class was generated from the following files: