source: trunk/source/processes/hadronic/models/binary_cascade/src/G4RKPropagation.cc @ 1347

Last change on this file since 1347 was 1196, checked in by garnier, 15 years ago

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