ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Absorber.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4Absorber.cc
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 "G4PhysicalConstants.hh"
35 #include "G4SystemOfUnits.hh"
36 #include "G4LorentzRotation.hh"
37 
38 G4Absorber::G4Absorber(G4double cutOnP)
39 {
40  theCutOnP = cutOnP;
41  theAbsorbers = new G4KineticTrackVector;
42  theProducts = new G4KineticTrackVector;
43 }
44 
45 
46 G4Absorber::~G4Absorber()
47 {
48  delete theAbsorbers;
49  delete theProducts;
50 }
51 
52 
53 bool G4Absorber::WillBeAbsorbed(const G4KineticTrack & kt)
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 }
69 
70 
71 
72 G4bool G4Absorber::Absorb(G4KineticTrack & kt, G4KineticTrackVector & tgt)
73 {
74  if(!FindAbsorbers(kt, tgt))
75  return false;
76  return FindProducts(kt);
77 }
78 
79 
80 G4bool G4Absorber::FindAbsorbers(G4KineticTrack & kt,
81  G4KineticTrackVector & tgt)
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 }
171 
172 
173 
174 G4bool G4Absorber::FindProducts(G4KineticTrack & kt)
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 }
283 
284 
285 
286 G4ThreeVector G4Absorber::GetRandomDirection()
287 {
288  G4double theta = 2.0*G4UniformRand()-1.0;
289  theta = std::acos(theta);
290  G4double phi = G4UniformRand()*2*pi;
291  G4ThreeVector direction(std::sin(theta)*std::cos(phi), std::sin(theta)*std::sin(phi), std::cos(theta));
292  return direction;
293 }
294 
295 
296 
297 
298 
299 
300