Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RKPropagation.cc
Go to the documentation of this file.
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  theMomentumTranfer=G4ThreeVector(0,0,0);
483  // G4cout <<"Stepper input"<<kt->GetTrackingMomentum()<<G4endl;
484  // create the integrator stepper
485  // G4Mag_EqRhs * equation = mapIter->second;
486  G4Mag_EqRhs * equation = (*theEquationMap)[kt->GetDefinition()->GetPDGEncoding()];
487  G4MagIntegratorStepper * stepper = new G4ClassicalRK4(equation);
488 
489  // create the integrator driver
490  G4double hMin = 1.0e-25*second; // arbitrary choice. Means 0.03 fm at c
491  G4MagInt_Driver * driver = new G4MagInt_Driver(hMin, stepper);
492 
493  // Temporary: use driver->AccurateAdvance()
494  // create the G4FieldTrack needed by AccurateAdvance
495  G4double curveLength = 0;
496  G4FieldTrack track(kt->GetPosition(),
497  kt->GetTrackingMomentum().vect().unit(), // momentum direction
498  curveLength, // curvelength
499  kt->GetTrackingMomentum().e()-kt->GetActualMass(), // kinetic energy
500  kt->GetActualMass(), // restmass
501  kt->GetTrackingMomentum().beta()*c_light); // velocity
502  // integrate
503  G4double eps = 0.01;
504  // G4cout << "currTimeStep = " << currTimeStep << G4endl;
505  if(!driver->AccurateAdvance(track, timeStep, eps))
506  { // cannot track this particle
507 #ifdef debug_1_RKPropagation
508  std::cerr << "G4RKPropagation::FieldTransport() warning: integration error."
509  << G4endl << "position " << kt->GetPosition() << " 4mom " <<kt->GetTrackingMomentum()
510  <<G4endl << " timestep " <<timeStep
511  << G4endl;
512 #endif
513  delete driver;
514  delete stepper;
515  return false;
516  }
517  /*
518  G4cout <<" E/Field/Sum be4 : " <<mom.e() << " / "
519  << (*theFieldMap)[encoding]->GetField(pos) << " / "
520  << mom.e()+(*theFieldMap)[encoding]->GetField(pos)
521  << G4endl;
522  */
523 
524  // Correct for momentum ( thus energy) transfered to nucleus, boost particle into moving nuclues frame.
525  G4ThreeVector MomentumTranfer = kt->GetTrackingMomentum().vect() - track.GetMomentum();
526  G4ThreeVector boost= MomentumTranfer / std::sqrt (MomentumTranfer.mag2() +sqr(theNucleus->GetMass()));
527 
528  // update the kt
529  kt->SetPosition(track.GetPosition());
530  G4LorentzVector mom(track.GetMomentum(),std::sqrt(track.GetMomentum().mag2() + sqr(kt->GetActualMass())));
531  mom *= G4LorentzRotation( boost );
532  theMomentumTranfer += ( kt->GetTrackingMomentum() - mom ).vect();
533  kt->SetTrackingMomentum(mom);
534 
535  // G4cout <<"Stepper output"<<kt<<" "<<kt->GetTrackingMomentum()<<" "<<kt->GetPosition()<<G4endl;
536  /*
537  * G4ThreeVector MomentumTranfer2=kt->GetTrackingMomentum().vect() - mom.vect();
538  * G4cout << " MomentumTransfer/corrected" << MomentumTranfer << " " << MomentumTranfer.mag()
539  * << " " << MomentumTranfer2 << " " << MomentumTranfer2.mag() << " "
540  * << MomentumTranfer-MomentumTranfer2 << " "<<
541  * MomentumTranfer-MomentumTranfer2.mag() << " " << G4endl;
542  * G4cout <<" E/Field/Sum aft : " <<mom.e() << " / "
543  * << " / " << (*theFieldMap)[encoding]->GetField(pos)<< " / "
544  * << mom.e()+(*theFieldMap)[encoding]->GetField(pos)
545  * << G4endl;
546  */
547 
548  delete driver;
549  delete stepper;
550  return true;
551 }
552 
553 //----------------------------------------------------------------------------
554 G4bool G4RKPropagation::FreeTransport(G4KineticTrack * kt, const G4double timeStep)
555 //----------------------------------------------------------------------------
556 {
557  G4ThreeVector newpos = kt->GetPosition() +
558  timeStep*c_light/kt->GetTrackingMomentum().e() * kt->GetTrackingMomentum().vect();
559  kt->SetPosition(newpos);
560  return true;
561 }
562 
563 /*
564 G4bool G4RKPropagation::WillBeCaptured(const G4KineticTrack * kt)
565 {
566  G4double radius = theOuterRadius;
567 
568 // evaluate the final energy. Il will be captured if newE or newP < 0
569  G4ParticleDefinition * definition = kt->GetDefinition();
570  G4double mass = definition->GetPDGMass();
571  G4ThreeVector pos = kt->GetPosition();
572  G4LorentzVector mom = kt->GetTrackingMomentum();
573  G4VNuclearField * field = (*theFieldMap)[definition->GetPDGEncoding()];
574  G4ThreeVector newPos(0, 0, radius); // to get the field on the surface
575 
576  G4double newE = mom.e()+field->GetField(pos)-field->GetField(newPos);
577 
578  return ((newE < mass) ? false : true);
579 }
580  */
581 
582 
583 
584 //----------------------------------------------------------------------------
585 G4bool G4RKPropagation::GetSphereIntersectionTimes(const G4double radius,
586  //----------------------------------------------------------------------------
587  const G4ThreeVector & currentPos,
588  const G4LorentzVector & momentum,
589  G4double & t1, G4double & t2)
590 {
591  G4ThreeVector speed = momentum.vect()/momentum.e(); // boost vector
592  G4double scalarProd = currentPos.dot(speed);
593  G4double speedMag2 = speed.mag2();
594  G4double sqrtArg = scalarProd*scalarProd -
595  speedMag2*(currentPos.mag2()-radius*radius);
596  if(sqrtArg <= 0.) // particle will not intersect the sphere
597  {
598  // G4cout << " GetSphereIntersectionTimes sqrtArg negative: " << sqrtArg << G4endl;
599  return false;
600  }
601  t1 = (-scalarProd - std::sqrt(sqrtArg))/speedMag2/c_light;
602  t2 = (-scalarProd + std::sqrt(sqrtArg))/speedMag2/c_light;
603  return true;
604 }
605 
606 //----------------------------------------------------------------------------
607 G4bool G4RKPropagation::GetSphereIntersectionTimes(const G4KineticTrack * kt,
608  G4double & t1, G4double & t2)
609 {
610  G4double radius = theOuterRadius + 3*fermi; // "safety" of 3 fermi
611  G4ThreeVector speed = kt->GetTrackingMomentum().vect()/kt->GetTrackingMomentum().e(); // bost vector
612  G4double scalarProd = kt->GetPosition().dot(speed);
613  G4double speedMag2 = speed.mag2();
614  G4double sqrtArg = scalarProd*scalarProd -
615  speedMag2*(kt->GetPosition().mag2()-radius*radius);
616  if(sqrtArg <= 0.) // particle will not intersect the sphere
617  {
618  return false;
619  }
620  t1 = (-scalarProd - std::sqrt(sqrtArg))/speedMag2/c_light;
621  t2 = (-scalarProd + std::sqrt(sqrtArg))/speedMag2/c_light;
622  return true;
623 }
624 
625 // Implementation methods
626 
627 //----------------------------------------------------------------------------
628 void G4RKPropagation::delete_FieldsAndMap(
629  //----------------------------------------------------------------------------
630  std::map <G4int, G4VNuclearField *, std::less<G4int> > * aMap)
631 {
632  if(aMap)
633  {
634  std::map <G4int, G4VNuclearField *, std::less<G4int> >::iterator cur;
635  for(cur = aMap->begin(); cur != aMap->end(); ++cur)
636  delete (*cur).second;
637 
638  aMap->clear();
639  delete aMap;
640  }
641 
642 }
643 
644 //----------------------------------------------------------------------------
645 void G4RKPropagation::delete_EquationsAndMap(
646  //----------------------------------------------------------------------------
647  std::map <G4int, G4Mag_EqRhs *, std::less<G4int> > * aMap)
648 {
649  if(aMap)
650  {
651  std::map <G4int, G4Mag_EqRhs *, std::less<G4int> >::iterator cur;
652  for(cur = aMap->begin(); cur != aMap->end(); ++cur)
653  delete (*cur).second;
654 
655  aMap->clear();
656  delete aMap;
657  }
658 }
Generated on Sat May 25 2013 14:33:43 for Geant4 by   doxygen 1.8.4