Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

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

Last change on this file since 1036 was 819, checked in by garnier, 17 years ago
import all except CVS
File size: 26.7 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
78	G4RKPropagation::G4RKPropagation() : theNucleus(0),
79	theFieldMap(0), theEquationMap(0),
80	theField(0)
81	{ }
82
83
84	G4RKPropagation::G4RKPropagation(const G4RKPropagation &) :
85	G4VFieldPropagation()
86	{ }
87
88
89	G4RKPropagation::~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
102	const G4RKPropagation & G4RKPropagation::operator=(const G4RKPropagation &)
103	{
104	throw G4HadronicException(__FILE__, __LINE__, "G4RKPropagation::operator= meant not to be accessible");
105	return *this;
106	}
107
108	G4int G4RKPropagation::operator==(const G4RKPropagation &) const
109	{
110	throw G4HadronicException(__FILE__, __LINE__, "G4RKPropagation::operator== meant not to be accessible");
111	return 0;
112	}
113
114	G4int 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	//----------------------------------------------------------------------------
124	void 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	//----------------------------------------------------------------------------
248	void G4RKPropagation::Transport(G4KineticTrackVector & active,
249	//----------------------------------------------------------------------------
250	const G4KineticTrackVector &,
251	G4double timeStep)
252	{
253	// reset momentum transfer to field
254	theMomentumTranfer=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	std::map <G4int, G4VNuclearField*, std::less<G4int> >::iterator fieldIter= theFieldMap->find(encoding);
265
266	G4VNuclearField* currentField=0;
267	if ( fieldIter != theFieldMap->end() ) currentField=fieldIter->second;
268
269	// debug
270	// if ( timeStep > 1e30 ) {
271	// G4cout << " Name :" << kt->GetDefinition()->GetParticleName() << G4endl;
272	// }
273
274	// Get the time of intersections with the nucleus surface.
275	G4double t_enter, t_leave;
276	// if the particle does not intersecate with the nucleus go to next particle
277	if(!GetSphereIntersectionTimes(kt, t_enter, t_leave))
278	{
279	kt->SetState(G4KineticTrack::miss_nucleus);
280	continue;
281	}
282
283
284	#ifdef debug_1_RKPropagation
285	G4cout <<" kt,timeStep, Intersection times tenter, tleave "
286	<<kt<<" "<< currTimeStep << " / " << t_enter << " / " << t_leave <<G4endl;
287	#endif
288
289	// if the particle is already outside nucleus go to next @@GF should never happen? check!
290	// does happen for particles added as late....
291	// if(t_leave < 0 )
292	// {
293	// throw G4HadronicException(__FILE__, __LINE__, "G4RKPropagation:: Attempt to track particle past a nucleus");
294	// continue;
295	// }
296
297	// Apply a straight line propagation for particle types
298	// not included in the model
299	if( ! currentField )
300	{
301	if(currTimeStep == DBL_MAX)currTimeStep = t_leave*1.05;
302	FreeTransport(kt, currTimeStep);
303	if ( currTimeStep >= t_leave )
304	{
305	if ( kt->GetState() == G4KineticTrack::inside )
306	{ kt->SetState(G4KineticTrack::gone_out); }
307	else
308	{ kt->SetState(G4KineticTrack::miss_nucleus);}
309	}
310	continue;
311	}
312
313	if(t_enter > 0) // the particle is out. Transport free to the surface
314	{
315	if(t_enter > currTimeStep) // the particle won't enter the nucleus
316	{
317	FreeTransport(kt, currTimeStep);
318	continue;
319	}
320	else
321	{
322	FreeTransport(kt, t_enter); // go to surface
323	currTimeStep -= t_enter;
324	t_leave -= t_enter; // time left to leave nucleus
325	// on the surface the particle loose the barrier energy
326	// G4double newE = mom.e()-(*theFieldMap)[encoding]->GetBarrier();
327	// GetField = Barrier + FermiPotential
328	G4double newE = kt->GetTrackingMomentum().e()-currentField->GetField(kt->GetPosition());
329	// G4cout << " enter nucleus, E out/in: " << kt->GetTrackingMomentum().e() << " / " << newE <<G4endl;
330	// G4cout << " the Field "<< currentField->GetField(kt->GetPosition()) << " "<< kt->GetPosition()<<G4endl;
331	// G4cout << " the particle "<<kt->GetDefinition()->GetParticleName()<<G4endl;
332	if(newE <= kt->GetActualMass()) // the particle cannot enter the nucleus
333	{
334	// FixMe: should be "pushed back?"
335	// for the moment take it past the nucleus, so we'll not worry next time..
336	FreeTransport(kt, 1.1*t_leave); // take past nucleus
337	kt->SetState(G4KineticTrack::miss_nucleus);
338	// G4cout << "G4RKPropagation: Warning particle cannot enter Nucleus :" << G4endl;
339	// G4cout << " enter nucleus, E out/in: " << kt->GetTrackingMomentum().e() << " / " << newE <<G4endl;
340	// G4cout << " the Field "<< currentField->GetField(kt->GetPosition()) << " "<< kt->GetPosition()<<G4endl;
341	// G4cout << " the particle "<<kt->GetDefinition()->GetParticleName()<<G4endl;
342	continue;
343	}
344	//
345	G4double newP = std::sqrt(newE*newE- sqr(kt->GetActualMass()));
346	G4LorentzVector new4Mom(newP*kt->GetTrackingMomentum().vect().unit(), newE);
347	G4ThreeVector transfer(kt->GetTrackingMomentum().vect()-new4Mom.vect());
348	G4ThreeVector boost= transfer / std::sqrt(transfer.mag2() + sqr(theNucleus->GetMass()));
349	new4Mom*=G4LorentzRotation(boost);
350	kt->SetTrackingMomentum(new4Mom);
351	kt->SetState(G4KineticTrack::inside);
352	// G4cout <<" Enter Nucleus - E/Field/Sum: " <<kt->GetTrackingMomentum().e() << " / "
353	// << (*theFieldMap)[encoding]->GetField(kt->GetPosition()) << " / "
354	// << kt->GetTrackingMomentum().e()-currentField->GetField(kt->GetPosition())
355	// << G4endl
356	// << " Barrier / field just inside nucleus (0.9999*kt->GetPosition())"
357	// << (*theFieldMap)[encoding]->GetBarrier() << " / "
358	// << (theFieldMap)[encoding]->GetField(0.9999kt->GetPosition())
359	// << G4endl;
360	}
361	}
362
363	// FixMe: should I add a control on theCutOnP here?
364	// Transport the particle into the nucleus
365	// G4cerr << "RKPropagation t_leave, curTimeStep " <<t_leave << " " <<currTimeStep<<G4endl;
366	G4bool is_exiting=false;
367	if(currTimeStep > t_leave) // particle will exit from the nucleus
368	{
369	currTimeStep = t_leave;
370	is_exiting=true;
371	}
372
373	#ifdef debug_1_RKPropagation
374	G4cerr << "RKPropagation is_exiting?, t_leave, curTimeStep " <<is_exiting<<" "<<t_leave << " " <<currTimeStep<<G4endl;
375	G4cout << "RKPropagation Ekin, field, projectile potential, p "
376	<< kt->GetTrackingMomentum().e() - kt->GetTrackingMomentum().mag() << " "
377	<< kt->GetPosition()<<" "
378	<< G4endl << currentField->GetField(kt->GetPosition()) << " "
379	<< kt->GetProjectilePotential()<< G4endl
380	<< kt->GetTrackingMomentum()
381	<< G4endl;
382	#endif
383
384	G4LorentzVector momold=kt->GetTrackingMomentum();
385	G4ThreeVector posold=kt->GetPosition();
386
387	// if (currentField->GetField(kt->GetPosition()) > kt->GetProjectilePotential() \|\|
388	if (currTimeStep > 0 &&
389	! FieldTransport(kt, currTimeStep)) {
390	FreeTransport(kt,currTimeStep);
391	}
392
393	#ifdef debug_1_RKPropagation
394	G4cout << "RKPropagation Ekin, field, p "
395	<< kt->GetTrackingMomentum().e() - kt->GetTrackingMomentum().mag() << " "
396	<< G4endl << currentField->GetField(kt->GetPosition())<< G4endl
397	<< kt->GetTrackingMomentum()
398	// << G4endl;
399	<< "delta p " << momold-kt->GetTrackingMomentum() << G4endl
400	<< "del pos " << posold-kt->GetPosition()
401	<< G4endl;
402	#endif
403
404	// complete the transport
405	// FixMe: in some cases there could be a significant
406	// part to do still in the nucleus, or we stepped to far... depending on
407	// slope of potential
408	G4double t_in=-1, t_out=0; // set onto boundary.
409
410	// should go out, or are already out by a too long step..
411	if(is_exiting \|\|
412	(GetSphereIntersectionTimes(kt, t_in, t_out) &&t_in<0 && t_out<=0 )) // particle is exiting
413	{
414	if(t_in < 0 && t_out >= 0) //still inside, transport safely out.
415	{
416	// transport free to a position that is surely out of the nucleus, to avoid
417	// a new transportation and a new adding the barrier next loop.
418	G4ThreeVector savePos = kt->GetPosition();
419	FreeTransport(kt, t_out);
420	// and evaluate the right the energy
421	G4double newE=kt->GetTrackingMomentum().e();
422
423	// G4cout << " V pos/savePos << "
424	// << (*theFieldMap)[encoding]->GetField(kt->GetPosition())<< " / "
425	// << (*theFieldMap)[encoding]->GetField(savePos)
426	// << G4endl;
427
428	if ( std::abs(currentField->GetField(savePos)) > 0. &&
429	std::abs(currentField->GetField(kt->GetPosition())) > 0.)
430	{ // FixMe GF: savePos/pos may be out of nucleus, where GetField(..)=0
431	// This wrongly adds or subtracts the Barrier here while
432	// this is done later.
433	newE += currentField->GetField(savePos)
434	- currentField->GetField(kt->GetPosition());
435	}
436
437	// G4cout << " go border nucleus, E in/border: " << kt->GetTrackingMomentum() << " / " << newE <<G4endl;
438
439	if(newE < kt->GetActualMass())
440	{
441	#ifdef debug_1_RKPropagation
442	G4cout << "RKPropagation-Transport: problem with particle exiting - ignored" << G4endl;
443	G4cout << " cannot leave nucleus, E in/out: " << kt->GetTrackingMomentum() << " / " << newE <<G4endl;
444	#endif
445	if (kt->GetDefinition() == G4Proton::Proton() \|\|
446	kt->GetDefinition() == G4Neutron::Neutron() ) {
447	kt->SetState(G4KineticTrack::captured);
448	} else {
449	kt->SetState(G4KineticTrack::gone_out); //@@GF tofix
450	}
451	continue; // the particle cannot exit the nucleus
452	}
453	G4double newP = std::sqrt(newE*newE- sqr(kt->GetActualMass()));
454	G4LorentzVector new4Mom(newP*kt->GetTrackingMomentum().vect().unit(), newE);
455	G4ThreeVector transfer(kt->GetTrackingMomentum().vect()-new4Mom.vect());
456	G4ThreeVector boost= transfer / std::sqrt(transfer.mag2() + sqr(theNucleus->GetMass()));
457	new4Mom*=G4LorentzRotation(boost);
458	kt->SetTrackingMomentum(new4Mom);
459	}
460	// add the potential barrier
461	// FixMe the Coulomb field is not parallel to mom, this is simple approximation
462	G4double newE = kt->GetTrackingMomentum().e()+currentField->GetField(kt->GetPosition());
463	if(newE < kt->GetActualMass())
464	{ // the particle cannot exit the nucleus @@@ GF check.
465	#ifdef debug_1_RKPropagation
466	G4cout << " cannot leave nucleus, E in/out: " << kt->GetTrackingMomentum() << " / " << newE <<G4endl;
467	#endif
468	if (kt->GetDefinition() == G4Proton::Proton() \|\|
469	kt->GetDefinition() == G4Neutron::Neutron() ) {
470	kt->SetState(G4KineticTrack::captured);
471	} else {
472	kt->SetState(G4KineticTrack::gone_out); //@@GF tofix
473	}
474	continue;
475	}
476	G4double newP = std::sqrt(newE*newE- sqr(kt->GetActualMass()));
477	G4LorentzVector new4Mom(newP*kt->GetTrackingMomentum().vect().unit(), newE);
478	G4ThreeVector transfer(kt->GetTrackingMomentum().vect()-new4Mom.vect());
479	G4ThreeVector boost= transfer / std::sqrt(transfer.mag2() + sqr(theNucleus->GetMass()));
480	new4Mom*=G4LorentzRotation(boost);
481	kt->SetTrackingMomentum(new4Mom);
482	kt->SetState(G4KineticTrack::gone_out);
483	}
484
485
486
487	}
488
489	}
490
491
492	//----------------------------------------------------------------------------
493	G4ThreeVector G4RKPropagation::GetMomentumTransfer() const
494	//----------------------------------------------------------------------------
495	{
496	return theMomentumTranfer;
497	}
498
499
500	//----------------------------------------------------------------------------
501	G4bool G4RKPropagation::FieldTransport(G4KineticTrack * kt, const G4double timeStep)
502	//----------------------------------------------------------------------------
503	{
504	theMomentumTranfer=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	//----------------------------------------------------------------------------
576	G4bool G4RKPropagation::FreeTransport(G4KineticTrack * kt, const G4double timeStep)
577	//----------------------------------------------------------------------------
578	{
579	G4ThreeVector newpos = kt->GetPosition() +
580	timeStepc_light/kt->GetTrackingMomentum().e() kt->GetTrackingMomentum().vect();
581	kt->SetPosition(newpos);
582	return true;
583	}
584
585	/*
586	G4bool 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	//----------------------------------------------------------------------------
607	G4bool 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	speedMagspeedMag(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	//----------------------------------------------------------------------------
629	G4bool 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()-radiusradius);
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	//----------------------------------------------------------------------------
650	void 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	//----------------------------------------------------------------------------
667	void 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.

Download in other formats: