ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RKPropagation.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RKPropagation.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 //
27 // -------------------------------------------------------------------
28 // GEANT 4 class implementation file
29 //
30 // CERN, Geneva, Switzerland
31 //
32 // File name: G4RKPropagation.cc
33 //
34 // Author: Alessandro Brunengo (Alessandro.Brunengo@ge.infn.it)
35 //
36 // Creation date: 6 June 2000
37 // -------------------------------------------------------------------
38 #include "G4RKPropagation.hh"
39 #include "G4PhysicalConstants.hh"
40 #include "G4SystemOfUnits.hh"
41 // nuclear fields
42 #include "G4VNuclearField.hh"
43 #include "G4ProtonField.hh"
44 #include "G4NeutronField.hh"
45 #include "G4AntiProtonField.hh"
46 #include "G4KaonPlusField.hh"
47 #include "G4KaonMinusField.hh"
48 #include "G4KaonZeroField.hh"
49 #include "G4PionPlusField.hh"
50 #include "G4PionMinusField.hh"
51 #include "G4PionZeroField.hh"
52 #include "G4SigmaPlusField.hh"
53 #include "G4SigmaMinusField.hh"
54 #include "G4SigmaZeroField.hh"
55 // particles properties
56 #include "G4Proton.hh"
57 #include "G4Neutron.hh"
58 #include "G4AntiProton.hh"
59 #include "G4KaonPlus.hh"
60 #include "G4KaonMinus.hh"
61 #include "G4KaonZero.hh"
62 #include "G4PionPlus.hh"
63 #include "G4PionMinus.hh"
64 #include "G4PionZero.hh"
65 #include "G4SigmaPlus.hh"
66 #include "G4SigmaMinus.hh"
67 #include "G4SigmaZero.hh"
68 
69 #include "globals.hh"
70 
71 #include "G4KM_OpticalEqRhs.hh"
72 #include "G4KM_NucleonEqRhs.hh"
73 #include "G4ClassicalRK4.hh"
74 #include "G4MagIntegratorDriver.hh"
75 
76 #include "G4LorentzRotation.hh"
77 
78 // unsigned EncodingHashFun(const G4int& aEncoding);
79 
80 G4RKPropagation::G4RKPropagation() :
81 theOuterRadius(0), theNucleus(0),
82 theFieldMap(0), theEquationMap(0),
83 theField(0)
84 { }
85 
86 
87 G4RKPropagation::~G4RKPropagation()
88 {
89  // free theFieldMap memory
90  if(theFieldMap) delete_FieldsAndMap(theFieldMap);
91 
92  // free theEquationMap memory
93  if(theEquationMap) delete_EquationsAndMap(theEquationMap);
94 
95  if (theField) delete theField;
96 }
97 
98 //----------------------------------------------------------------------------
99 void G4RKPropagation::Init(G4V3DNucleus * nucleus)
100 //----------------------------------------------------------------------------
101 {
102 
103  // free theFieldMap memory
104  if(theFieldMap) delete_FieldsAndMap(theFieldMap);
105 
106  // free theEquationMap memory
107  if(theEquationMap) delete_EquationsAndMap(theEquationMap);
108 
109  if (theField) delete theField;
110 
111  // Initialize the nuclear field map.
112  theNucleus = nucleus;
113  theOuterRadius = theNucleus->GetOuterRadius();
114 
115  theFieldMap = new std::map <G4int, G4VNuclearField*, std::less<G4int> >;
116 
117  (*theFieldMap)[G4Proton::Proton()->GetPDGEncoding()] = new G4ProtonField(theNucleus);
118  (*theFieldMap)[G4Neutron::Neutron()->GetPDGEncoding()] = new G4NeutronField(theNucleus);
119  (*theFieldMap)[G4AntiProton::AntiProton()->GetPDGEncoding()] = new G4AntiProtonField(theNucleus);
120  (*theFieldMap)[G4KaonPlus::KaonPlus()->GetPDGEncoding()] = new G4KaonPlusField(theNucleus);
121  (*theFieldMap)[G4KaonMinus::KaonMinus()->GetPDGEncoding()] = new G4KaonMinusField(theNucleus);
122  (*theFieldMap)[G4KaonZero::KaonZero()->GetPDGEncoding()] = new G4KaonZeroField(theNucleus);
123  (*theFieldMap)[G4PionPlus::PionPlus()->GetPDGEncoding()] = new G4PionPlusField(theNucleus);
124  (*theFieldMap)[G4PionMinus::PionMinus()->GetPDGEncoding()] = new G4PionMinusField(theNucleus);
125  (*theFieldMap)[G4PionZero::PionZero()->GetPDGEncoding()] = new G4PionZeroField(theNucleus);
126  (*theFieldMap)[G4SigmaPlus::SigmaPlus()->GetPDGEncoding()] = new G4SigmaPlusField(theNucleus);
127  (*theFieldMap)[G4SigmaMinus::SigmaMinus()->GetPDGEncoding()] = new G4SigmaMinusField(theNucleus);
128  (*theFieldMap)[G4SigmaZero::SigmaZero()->GetPDGEncoding()] = new G4SigmaZeroField(theNucleus);
129 
130  theEquationMap = new std::map <G4int, G4Mag_EqRhs*, std::less<G4int> >;
131 
132  // theField needed by the design of G4Mag_eqRhs
133  theField = new G4KM_DummyField; //Field not needed for integration
134  G4KM_OpticalEqRhs * opticalEq;
135  G4KM_NucleonEqRhs * nucleonEq;
136  G4double mass;
137  G4double opticalCoeff;
138 
139  nucleonEq = new G4KM_NucleonEqRhs(theField, theNucleus);
140  mass = G4Proton::Proton()->GetPDGMass();
141  nucleonEq->SetMass(mass);
142  (*theEquationMap)[G4Proton::Proton()->GetPDGEncoding()] = nucleonEq;
143 
144  nucleonEq = new G4KM_NucleonEqRhs(theField, theNucleus);
145  mass = G4Neutron::Neutron()->GetPDGMass();
146  nucleonEq->SetMass(mass);
147  (*theEquationMap)[G4Neutron::Neutron()->GetPDGEncoding()] = nucleonEq;
148 
149  opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
150  mass = G4AntiProton::AntiProton()->GetPDGMass();
151  opticalCoeff =
152  (*theFieldMap)[G4AntiProton::AntiProton()->GetPDGEncoding()]->GetCoeff();
153  opticalEq->SetFactor(mass,opticalCoeff);
154  (*theEquationMap)[G4AntiProton::AntiProton()->GetPDGEncoding()] = opticalEq;
155 
156  opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
157  mass = G4KaonPlus::KaonPlus()->GetPDGMass();
158  opticalCoeff =
159  (*theFieldMap)[G4KaonPlus::KaonPlus()->GetPDGEncoding()]->GetCoeff();
160  opticalEq->SetFactor(mass,opticalCoeff);
161  (*theEquationMap)[G4KaonPlus::KaonPlus()->GetPDGEncoding()] = opticalEq;
162 
163  opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
164  mass = G4KaonMinus::KaonMinus()->GetPDGMass();
165  opticalCoeff =
166  (*theFieldMap)[G4KaonMinus::KaonMinus()->GetPDGEncoding()]->GetCoeff();
167  opticalEq->SetFactor(mass,opticalCoeff);
168  (*theEquationMap)[G4KaonMinus::KaonMinus()->GetPDGEncoding()] = opticalEq;
169 
170  opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
171  mass = G4KaonZero::KaonZero()->GetPDGMass();
172  opticalCoeff =
173  (*theFieldMap)[G4KaonZero::KaonZero()->GetPDGEncoding()]->GetCoeff();
174  opticalEq->SetFactor(mass,opticalCoeff);
175  (*theEquationMap)[G4KaonZero::KaonZero()->GetPDGEncoding()] = opticalEq;
176 
177  opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
178  mass = G4PionPlus::PionPlus()->GetPDGMass();
179  opticalCoeff =
180  (*theFieldMap)[G4PionPlus::PionPlus()->GetPDGEncoding()]->GetCoeff();
181  opticalEq->SetFactor(mass,opticalCoeff);
182  (*theEquationMap)[G4PionPlus::PionPlus()->GetPDGEncoding()] = opticalEq;
183 
184  opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
185  mass = G4PionMinus::PionMinus()->GetPDGMass();
186  opticalCoeff =
187  (*theFieldMap)[G4PionMinus::PionMinus()->GetPDGEncoding()]->GetCoeff();
188  opticalEq->SetFactor(mass,opticalCoeff);
189  (*theEquationMap)[G4PionMinus::PionMinus()->GetPDGEncoding()] = opticalEq;
190 
191  opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
192  mass = G4PionZero::PionZero()->GetPDGMass();
193  opticalCoeff =
194  (*theFieldMap)[G4PionZero::PionZero()->GetPDGEncoding()]->GetCoeff();
195  opticalEq->SetFactor(mass,opticalCoeff);
196  (*theEquationMap)[G4PionZero::PionZero()->GetPDGEncoding()] = opticalEq;
197 
198  opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
199  mass = G4SigmaPlus::SigmaPlus()->GetPDGMass();
200  opticalCoeff =
201  (*theFieldMap)[G4SigmaPlus::SigmaPlus()->GetPDGEncoding()]->GetCoeff();
202  opticalEq->SetFactor(mass,opticalCoeff);
203  (*theEquationMap)[G4SigmaPlus::SigmaPlus()->GetPDGEncoding()] = opticalEq;
204 
205  opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
206  mass = G4SigmaMinus::SigmaMinus()->GetPDGMass();
207  opticalCoeff =
208  (*theFieldMap)[G4SigmaMinus::SigmaMinus()->GetPDGEncoding()]->GetCoeff();
209  opticalEq->SetFactor(mass,opticalCoeff);
210  (*theEquationMap)[G4SigmaMinus::SigmaMinus()->GetPDGEncoding()] = opticalEq;
211 
212  opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
213  mass = G4SigmaZero::SigmaZero()->GetPDGMass();
214  opticalCoeff =
215  (*theFieldMap)[G4SigmaZero::SigmaZero()->GetPDGEncoding()]->GetCoeff();
216  opticalEq->SetFactor(mass,opticalCoeff);
217  (*theEquationMap)[G4SigmaZero::SigmaZero()->GetPDGEncoding()] = opticalEq;
218 }
219 
220 
221 //#define debug_1_RKPropagation 1
222 //----------------------------------------------------------------------------
223 void G4RKPropagation::Transport(G4KineticTrackVector & active,
224  //----------------------------------------------------------------------------
225  const G4KineticTrackVector &,
226  G4double timeStep)
227 {
228  // reset momentum transfer to field
229  theMomentumTranfer=G4ThreeVector(0,0,0);
230 
231  // Loop over tracks
232 
233  std::vector<G4KineticTrack *>::iterator i;
234  for(i = active.begin(); i != active.end(); ++i)
235  {
236  G4double currTimeStep = timeStep;
237  G4KineticTrack * kt = *i;
238  G4int encoding = kt->GetDefinition()->GetPDGEncoding();
239 
240  std::map <G4int, G4VNuclearField*, std::less<G4int> >::iterator fieldIter= theFieldMap->find(encoding);
241 
242  G4VNuclearField* currentField=0;
243  if ( fieldIter != theFieldMap->end() ) currentField=fieldIter->second;
244 
245  // debug
246  // if ( timeStep > 1e30 ) {
247  // G4cout << " Name :" << kt->GetDefinition()->GetParticleName() << G4endl;
248  // }
249 
250  // Get the time of intersections with the nucleus surface.
251  G4double t_enter, t_leave;
252  // if the particle does not intersecate with the nucleus go to next particle
253  if(!GetSphereIntersectionTimes(kt, t_enter, t_leave))
254  {
255  kt->SetState(G4KineticTrack::miss_nucleus);
256  continue;
257  }
258 
259 
260 #ifdef debug_1_RKPropagation
261  G4cout <<" kt,timeStep, Intersection times tenter, tleave "
262  <<kt<< " / state= " <<kt->GetState() <<" / " <<" "<< currTimeStep << " / " << t_enter << " / " << t_leave <<G4endl;
263 #endif
264 
265  // if the particle is already outside nucleus go to next @@GF should never happen? check!
266  // does happen for particles added as late....
267  // if(t_leave < 0 )
268  // {
269  // throw G4HadronicException(__FILE__, __LINE__, "G4RKPropagation:: Attempt to track particle past a nucleus");
270  // continue;
271  // }
272 
273  // Apply a straight line propagation for particle types
274  // not included in the model
275  if( ! currentField )
276  {
277  if(currTimeStep == DBL_MAX)currTimeStep = t_leave*1.05;
278  FreeTransport(kt, currTimeStep);
279  if ( currTimeStep >= t_leave )
280  {
281  if ( kt->GetState() == G4KineticTrack::inside )
282  { kt->SetState(G4KineticTrack::gone_out); }
283  else
284  { kt->SetState(G4KineticTrack::miss_nucleus);}
285  } else if (kt->GetState() == G4KineticTrack::outside && currTimeStep >= t_enter ){
286  kt->SetState(G4KineticTrack::inside);
287  }
288 
289  continue;
290  }
291 
292  if(t_enter > 0) // the particle is out. Transport free to the surface
293  {
294  if(t_enter > currTimeStep) // the particle won't enter the nucleus
295  {
296  FreeTransport(kt, currTimeStep);
297  continue;
298  }
299  else
300  {
301  FreeTransport(kt, t_enter); // go to surface
302  currTimeStep -= t_enter;
303  t_leave -= t_enter; // time left to leave nucleus
304  // on the surface the particle loose the barrier energy
305  // G4double newE = mom.e()-(*theFieldMap)[encoding]->GetBarrier();
306  // GetField = Barrier + FermiPotential
307  G4double newE = kt->GetTrackingMomentum().e()-currentField->GetField(kt->GetPosition());
308 
309  if(newE <= kt->GetActualMass()) // the particle cannot enter the nucleus
310  {
311  // FixMe: should be "pushed back?"
312  // for the moment take it past the nucleus, so we'll not worry next time..
313  FreeTransport(kt, 1.1*t_leave); // take past nucleus
314  kt->SetState(G4KineticTrack::miss_nucleus);
315  // G4cout << "G4RKPropagation: Warning particle cannot enter Nucleus :" << G4endl;
316  // G4cout << " enter nucleus, E out/in: " << kt->GetTrackingMomentum().e() << " / " << newE <<G4endl;
317  // G4cout << " the Field "<< currentField->GetField(kt->GetPosition()) << " "<< kt->GetPosition()<<G4endl;
318  // G4cout << " the particle "<<kt->GetDefinition()->GetParticleName()<<G4endl;
319  continue;
320  }
321  //
322  G4double newP = std::sqrt(newE*newE- sqr(kt->GetActualMass()));
323  G4LorentzVector new4Mom(newP*kt->GetTrackingMomentum().vect().unit(), newE);
324  G4ThreeVector transfer(kt->GetTrackingMomentum().vect()-new4Mom.vect());
325  G4ThreeVector boost= transfer / std::sqrt(transfer.mag2() + sqr(theNucleus->GetMass()));
326  new4Mom*=G4LorentzRotation(boost);
327  kt->SetTrackingMomentum(new4Mom);
328  kt->SetState(G4KineticTrack::inside);
329 
330  /*
331  G4cout <<" Enter Nucleus - E/Field/Sum: " <<kt->GetTrackingMomentum().e() << " / "
332  << (*theFieldMap)[encoding]->GetField(kt->GetPosition()) << " / "
333  << kt->GetTrackingMomentum().e()-currentField->GetField(kt->GetPosition())
334  << G4endl
335  << " Barrier / field just inside nucleus (0.9999*kt->GetPosition())"
336  << (*theFieldMap)[encoding]->GetBarrier() << " / "
337  << (*theFieldMap)[encoding]->GetField(0.9999*kt->GetPosition())
338  << G4endl;
339  */
340  }
341  }
342 
343  // FixMe: should I add a control on theCutOnP here?
344  // Transport the particle into the nucleus
345  // G4cerr << "RKPropagation t_leave, curTimeStep " <<t_leave << " " <<currTimeStep<<G4endl;
346  G4bool is_exiting=false;
347  if(currTimeStep > t_leave) // particle will exit from the nucleus
348  {
349  currTimeStep = t_leave;
350  is_exiting=true;
351  }
352 
353 #ifdef debug_1_RKPropagation
354  G4cerr << "RKPropagation is_exiting?, t_leave, curTimeStep " <<is_exiting<<" "<<t_leave << " " <<currTimeStep<<G4endl;
355  G4cout << "RKPropagation Ekin, field, projectile potential, p "
356  << kt->GetTrackingMomentum().e() - kt->GetTrackingMomentum().mag() << " "
357  << kt->GetPosition()<<" "
358  << G4endl << currentField->GetField(kt->GetPosition()) << " "
359  << kt->GetProjectilePotential()<< G4endl
360  << kt->GetTrackingMomentum()
361  << G4endl;
362 #endif
363 
364  G4LorentzVector momold=kt->GetTrackingMomentum();
365  G4ThreeVector posold=kt->GetPosition();
366 
367  // if (currentField->GetField(kt->GetPosition()) > kt->GetProjectilePotential() ||
368  if (currTimeStep > 0 &&
369  ! FieldTransport(kt, currTimeStep)) {
370  FreeTransport(kt,currTimeStep);
371  }
372 
373 #ifdef debug_1_RKPropagation
374  G4cout << "RKPropagation Ekin, field, p "
375  << kt->GetTrackingMomentum().e() - kt->GetTrackingMomentum().mag() << " "
376  << G4endl << currentField->GetField(kt->GetPosition())<< G4endl
377  << kt->GetTrackingMomentum()
378  << G4endl
379  << "delta p " << momold-kt->GetTrackingMomentum() << G4endl
380  << "del pos " << posold-kt->GetPosition()
381  << G4endl;
382 #endif
383 
384  // complete the transport
385  // FixMe: in some cases there could be a significant
386  // part to do still in the nucleus, or we stepped to far... depending on
387  // slope of potential
388  G4double t_in=-1, t_out=0; // set onto boundary.
389 
390  // should go out, or are already out by a too long step..
391  if(is_exiting ||
392  (GetSphereIntersectionTimes(kt, t_in, t_out) &&t_in<0 && t_out<=0 )) // particle is exiting
393  {
394  if(t_in < 0 && t_out >= 0) //still inside, transport safely out.
395  {
396  // transport free to a position that is surely out of the nucleus, to avoid
397  // a new transportation and a new adding the barrier next loop.
398  G4ThreeVector savePos = kt->GetPosition();
399  FreeTransport(kt, t_out);
400  // and evaluate the right the energy
401  G4double newE=kt->GetTrackingMomentum().e();
402 
403  // G4cout << " V pos/savePos << "
404  // << (*theFieldMap)[encoding]->GetField(kt->GetPosition())<< " / "
405  // << (*theFieldMap)[encoding]->GetField(savePos)
406  // << G4endl;
407 
408  if ( std::abs(currentField->GetField(savePos)) > 0. &&
409  std::abs(currentField->GetField(kt->GetPosition())) > 0.)
410  { // FixMe GF: savePos/pos may be out of nucleus, where GetField(..)=0
411  // This wrongly adds or subtracts the Barrier here while
412  // this is done later.
413  newE += currentField->GetField(savePos)
414  - currentField->GetField(kt->GetPosition());
415  }
416 
417  // G4cout << " go border nucleus, E in/border: " << kt->GetTrackingMomentum() << " / " << newE <<G4endl;
418 
419  if(newE < kt->GetActualMass())
420  {
421 #ifdef debug_1_RKPropagation
422  G4cout << "RKPropagation-Transport: problem with particle exiting - ignored" << G4endl;
423  G4cout << " cannot leave nucleus, E in/out: " << kt->GetTrackingMomentum() << " / " << newE <<G4endl;
424 #endif
425  if (kt->GetDefinition() == G4Proton::Proton() ||
426  kt->GetDefinition() == G4Neutron::Neutron() ) {
427  kt->SetState(G4KineticTrack::captured);
428  } else {
429  kt->SetState(G4KineticTrack::gone_out); //@@GF tofix
430  }
431  continue; // the particle cannot exit the nucleus
432  }
433  G4double newP = std::sqrt(newE*newE- sqr(kt->GetActualMass()));
434  G4LorentzVector new4Mom(newP*kt->GetTrackingMomentum().vect().unit(), newE);
435  G4ThreeVector transfer(kt->GetTrackingMomentum().vect()-new4Mom.vect());
436  G4ThreeVector boost= transfer / std::sqrt(transfer.mag2() + sqr(theNucleus->GetMass()));
437  new4Mom*=G4LorentzRotation(boost);
438  kt->SetTrackingMomentum(new4Mom);
439  }
440  // add the potential barrier
441  // FixMe the Coulomb field is not parallel to mom, this is simple approximation
442  G4double newE = kt->GetTrackingMomentum().e()+currentField->GetField(kt->GetPosition());
443  if(newE < kt->GetActualMass())
444  { // the particle cannot exit the nucleus @@@ GF check.
445 #ifdef debug_1_RKPropagation
446  G4cout << " cannot leave nucleus, E in/out: " << kt->GetTrackingMomentum() << " / " << newE <<G4endl;
447 #endif
448  if (kt->GetDefinition() == G4Proton::Proton() ||
449  kt->GetDefinition() == G4Neutron::Neutron() ) {
450  kt->SetState(G4KineticTrack::captured);
451  } else {
452  kt->SetState(G4KineticTrack::gone_out); //@@GF tofix
453  }
454  continue;
455  }
456  G4double newP = std::sqrt(newE*newE- sqr(kt->GetActualMass()));
457  G4LorentzVector new4Mom(newP*kt->GetTrackingMomentum().vect().unit(), newE);
458  G4ThreeVector transfer(kt->GetTrackingMomentum().vect()-new4Mom.vect());
459  G4ThreeVector boost= transfer / std::sqrt(transfer.mag2() + sqr(theNucleus->GetMass()));
460  new4Mom*=G4LorentzRotation(boost);
461  kt->SetTrackingMomentum(new4Mom);
462  kt->SetState(G4KineticTrack::gone_out);
463  }
464 
465  }
466 
467 }
468 
469 
470 //----------------------------------------------------------------------------
471 G4ThreeVector G4RKPropagation::GetMomentumTransfer() const
472 //----------------------------------------------------------------------------
473 {
474  return theMomentumTranfer;
475 }
476 
477 
478 //----------------------------------------------------------------------------
479 G4bool G4RKPropagation::FieldTransport(G4KineticTrack * kt, const G4double timeStep)
480 //----------------------------------------------------------------------------
481 {
482  // G4cout <<"Stepper input"<<kt->GetTrackingMomentum()<<G4endl;
483  // create the integrator stepper
484  // G4Mag_EqRhs * equation = mapIter->second;
485  G4Mag_EqRhs * equation = (*theEquationMap)[kt->GetDefinition()->GetPDGEncoding()];
486  G4MagIntegratorStepper * stepper = new G4ClassicalRK4(equation);
487 
488  // create the integrator driver
489  G4double hMin = 1.0e-25*second; // arbitrary choice. Means 0.03 fm at c
490  G4MagInt_Driver * driver = new G4MagInt_Driver(hMin, stepper);
491 
492  // Temporary: use driver->AccurateAdvance()
493  // create the G4FieldTrack needed by AccurateAdvance
494  G4double curveLength = 0;
495  G4FieldTrack track(kt->GetPosition(),
496  kt->GetTrackingMomentum().vect().unit(), // momentum direction
497  curveLength, // curvelength
498  kt->GetTrackingMomentum().e()-kt->GetActualMass(), // kinetic energy
499  kt->GetActualMass(), // restmass
500  kt->GetTrackingMomentum().beta()*c_light); // velocity
501  // integrate
502  G4double eps = 0.01;
503  // G4cout << "currTimeStep = " << currTimeStep << G4endl;
504  if(!driver->AccurateAdvance(track, timeStep, eps))
505  { // cannot track this particle
506 #ifdef debug_1_RKPropagation
507  std::cerr << "G4RKPropagation::FieldTransport() warning: integration error."
508  << G4endl << "position " << kt->GetPosition() << " 4mom " <<kt->GetTrackingMomentum()
509  <<G4endl << " timestep " <<timeStep
510  << G4endl;
511 #endif
512  delete driver;
513  delete stepper;
514  return false;
515  }
516  /*
517  G4cout <<" E/Field/Sum be4 : " <<mom.e() << " / "
518  << (*theFieldMap)[encoding]->GetField(pos) << " / "
519  << mom.e()+(*theFieldMap)[encoding]->GetField(pos)
520  << G4endl;
521  */
522 
523  // Correct for momentum ( thus energy) transfered to nucleus, boost particle into moving nucleus frame.
524  G4ThreeVector MomentumTranfer = kt->GetTrackingMomentum().vect() - track.GetMomentum();
525  G4ThreeVector boost= MomentumTranfer / std::sqrt (MomentumTranfer.mag2() +sqr(theNucleus->GetMass()));
526 
527  // update the kt
528  kt->SetPosition(track.GetPosition());
529  G4LorentzVector mom(track.GetMomentum(),std::sqrt(track.GetMomentum().mag2() + sqr(kt->GetActualMass())));
530  mom *= G4LorentzRotation( boost );
531  theMomentumTranfer += ( kt->GetTrackingMomentum() - mom ).vect();
532  kt->SetTrackingMomentum(mom);
533 
534  // G4cout <<"Stepper output"<<kt<<" "<<kt->GetTrackingMomentum()<<" "<<kt->GetPosition()<<G4endl;
535  /*
536  * G4ThreeVector MomentumTranfer2=kt->GetTrackingMomentum().vect() - mom.vect();
537  * G4cout << " MomentumTransfer/corrected" << MomentumTranfer << " " << MomentumTranfer.mag()
538  * << " " << MomentumTranfer2 << " " << MomentumTranfer2.mag() << " "
539  * << MomentumTranfer-MomentumTranfer2 << " "<<
540  * MomentumTranfer-MomentumTranfer2.mag() << " " << G4endl;
541  * G4cout <<" E/Field/Sum aft : " <<mom.e() << " / "
542  * << " / " << (*theFieldMap)[encoding]->GetField(pos)<< " / "
543  * << mom.e()+(*theFieldMap)[encoding]->GetField(pos)
544  * << G4endl;
545  */
546 
547  delete driver;
548  delete stepper;
549  return true;
550 }
551 
552 //----------------------------------------------------------------------------
553 G4bool G4RKPropagation::FreeTransport(G4KineticTrack * kt, const G4double timeStep)
554 //----------------------------------------------------------------------------
555 {
556  G4ThreeVector newpos = kt->GetPosition() +
557  timeStep*c_light/kt->GetTrackingMomentum().e() * kt->GetTrackingMomentum().vect();
558  kt->SetPosition(newpos);
559  return true;
560 }
561 
562 /*
563 G4bool G4RKPropagation::WillBeCaptured(const G4KineticTrack * kt)
564 {
565  G4double radius = theOuterRadius;
566 
567 // evaluate the final energy. Il will be captured if newE or newP < 0
568  G4ParticleDefinition * definition = kt->GetDefinition();
569  G4double mass = definition->GetPDGMass();
570  G4ThreeVector pos = kt->GetPosition();
571  G4LorentzVector mom = kt->GetTrackingMomentum();
572  G4VNuclearField * field = (*theFieldMap)[definition->GetPDGEncoding()];
573  G4ThreeVector newPos(0, 0, radius); // to get the field on the surface
574 
575  G4double newE = mom.e()+field->GetField(pos)-field->GetField(newPos);
576 
577  return ((newE < mass) ? false : true);
578 }
579  */
580 
581 
582 
583 //----------------------------------------------------------------------------
584 G4bool G4RKPropagation::GetSphereIntersectionTimes(const G4double radius,
585  //----------------------------------------------------------------------------
586  const G4ThreeVector & currentPos,
587  const G4LorentzVector & momentum,
588  G4double & t1, G4double & t2)
589 {
590  G4ThreeVector speed = momentum.vect()/momentum.e(); // boost vector
591  G4double scalarProd = currentPos.dot(speed);
592  G4double speedMag2 = speed.mag2();
593  G4double sqrtArg = scalarProd*scalarProd -
594  speedMag2*(currentPos.mag2()-radius*radius);
595  if(sqrtArg <= 0.) // particle will not intersect the sphere
596  {
597  // G4cout << " GetSphereIntersectionTimes sqrtArg negative: " << sqrtArg << G4endl;
598  return false;
599  }
600  t1 = (-scalarProd - std::sqrt(sqrtArg))/speedMag2/c_light;
601  t2 = (-scalarProd + std::sqrt(sqrtArg))/speedMag2/c_light;
602  return true;
603 }
604 
605 //----------------------------------------------------------------------------
606 G4bool G4RKPropagation::GetSphereIntersectionTimes(const G4KineticTrack * kt,
607  G4double & t1, G4double & t2)
608 {
609  G4double radius = theOuterRadius + 3*fermi; // "safety" of 3 fermi
610  G4ThreeVector speed = kt->GetTrackingMomentum().vect()/kt->GetTrackingMomentum().e(); // bost vector
611  G4double scalarProd = kt->GetPosition().dot(speed);
612  G4double speedMag2 = speed.mag2();
613  G4double sqrtArg = scalarProd*scalarProd -
614  speedMag2*(kt->GetPosition().mag2()-radius*radius);
615  if(sqrtArg <= 0.) // particle will not intersect the sphere
616  {
617  return false;
618  }
619  t1 = (-scalarProd - std::sqrt(sqrtArg))/speedMag2/c_light;
620  t2 = (-scalarProd + std::sqrt(sqrtArg))/speedMag2/c_light;
621  return true;
622 }
623 
624 // Implementation methods
625 
626 //----------------------------------------------------------------------------
627 void G4RKPropagation::delete_FieldsAndMap(
628  //----------------------------------------------------------------------------
629  std::map <G4int, G4VNuclearField *, std::less<G4int> > * aMap)
630 {
631  if(aMap)
632  {
633  std::map <G4int, G4VNuclearField *, std::less<G4int> >::iterator cur;
634  for(cur = aMap->begin(); cur != aMap->end(); ++cur)
635  delete (*cur).second;
636 
637  aMap->clear();
638  delete aMap;
639  }
640 
641 }
642 
643 //----------------------------------------------------------------------------
644 void G4RKPropagation::delete_EquationsAndMap(
645  //----------------------------------------------------------------------------
646  std::map <G4int, G4Mag_EqRhs *, std::less<G4int> > * aMap)
647 {
648  if(aMap)
649  {
650  std::map <G4int, G4Mag_EqRhs *, std::less<G4int> >::iterator cur;
651  for(cur = aMap->begin(); cur != aMap->end(); ++cur)
652  delete (*cur).second;
653 
654  aMap->clear();
655  delete aMap;
656  }
657 }