| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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 (Al
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& aEnco
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(theFiel
91
92 // free theEquationMap memory
93 if(theEquationMap) delete_EquationsAndMap(t
94
95 if (theField) delete theField;
96 }
97
98 //--------------------------------------------
99 void G4RKPropagation::Init(G4V3DNucleus * nucl
100 //--------------------------------------------
101 {
102
103 // free theFieldMap memory
104 if(theFieldMap) delete_FieldsAndMap(theFiel
105
106 // free theEquationMap memory
107 if(theEquationMap) delete_EquationsAndMap(t
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, G4VNucle
116
117 (*theFieldMap)[G4Proton::Proton()->GetPDGEn
118 (*theFieldMap)[G4Neutron::Neutron()->GetPDG
119 (*theFieldMap)[G4AntiProton::AntiProton()->
120 (*theFieldMap)[G4KaonPlus::KaonPlus()->GetP
121 (*theFieldMap)[G4KaonMinus::KaonMinus()->Ge
122 (*theFieldMap)[G4KaonZero::KaonZero()->GetP
123 (*theFieldMap)[G4PionPlus::PionPlus()->GetP
124 (*theFieldMap)[G4PionMinus::PionMinus()->Ge
125 (*theFieldMap)[G4PionZero::PionZero()->GetP
126 (*theFieldMap)[G4SigmaPlus::SigmaPlus()->Ge
127 (*theFieldMap)[G4SigmaMinus::SigmaMinus()->
128 (*theFieldMap)[G4SigmaZero::SigmaZero()->Ge
129
130 theEquationMap = new std::map <G4int, G4Mag
131
132 // theField needed by the design of G4Mag_e
133 theField = new G4KM_DummyField; //Field
134 G4KM_OpticalEqRhs * opticalEq;
135 G4KM_NucleonEqRhs * nucleonEq;
136 G4double mass;
137 G4double opticalCoeff;
138
139 nucleonEq = new G4KM_NucleonEqRhs(theField,
140 mass = G4Proton::Proton()->GetPDGMass();
141 nucleonEq->SetMass(mass);
142 (*theEquationMap)[G4Proton::Proton()->GetPD
143
144 nucleonEq = new G4KM_NucleonEqRhs(theField,
145 mass = G4Neutron::Neutron()->GetPDGMass();
146 nucleonEq->SetMass(mass);
147 (*theEquationMap)[G4Neutron::Neutron()->Get
148
149 opticalEq = new G4KM_OpticalEqRhs(theField,
150 mass = G4AntiProton::AntiProton()->GetPDGMa
151 opticalCoeff =
152 (*theFieldMap)[G4AntiProton::AntiProt
153 opticalEq->SetFactor(mass,opticalCoeff);
154 (*theEquationMap)[G4AntiProton::AntiProton(
155
156 opticalEq = new G4KM_OpticalEqRhs(theField,
157 mass = G4KaonPlus::KaonPlus()->GetPDGMass()
158 opticalCoeff =
159 (*theFieldMap)[G4KaonPlus::KaonPlus()
160 opticalEq->SetFactor(mass,opticalCoeff);
161 (*theEquationMap)[G4KaonPlus::KaonPlus()->G
162
163 opticalEq = new G4KM_OpticalEqRhs(theField,
164 mass = G4KaonMinus::KaonMinus()->GetPDGMass
165 opticalCoeff =
166 (*theFieldMap)[G4KaonMinus::KaonMinus
167 opticalEq->SetFactor(mass,opticalCoeff);
168 (*theEquationMap)[G4KaonMinus::KaonMinus()-
169
170 opticalEq = new G4KM_OpticalEqRhs(theField,
171 mass = G4KaonZero::KaonZero()->GetPDGMass()
172 opticalCoeff =
173 (*theFieldMap)[G4KaonZero::KaonZero()
174 opticalEq->SetFactor(mass,opticalCoeff);
175 (*theEquationMap)[G4KaonZero::KaonZero()->G
176
177 opticalEq = new G4KM_OpticalEqRhs(theField,
178 mass = G4PionPlus::PionPlus()->GetPDGMass()
179 opticalCoeff =
180 (*theFieldMap)[G4PionPlus::PionPlus()
181 opticalEq->SetFactor(mass,opticalCoeff);
182 (*theEquationMap)[G4PionPlus::PionPlus()->G
183
184 opticalEq = new G4KM_OpticalEqRhs(theField,
185 mass = G4PionMinus::PionMinus()->GetPDGMass
186 opticalCoeff =
187 (*theFieldMap)[G4PionMinus::PionMinus
188 opticalEq->SetFactor(mass,opticalCoeff);
189 (*theEquationMap)[G4PionMinus::PionMinus()-
190
191 opticalEq = new G4KM_OpticalEqRhs(theField,
192 mass = G4PionZero::PionZero()->GetPDGMass()
193 opticalCoeff =
194 (*theFieldMap)[G4PionZero::PionZero()
195 opticalEq->SetFactor(mass,opticalCoeff);
196 (*theEquationMap)[G4PionZero::PionZero()->G
197
198 opticalEq = new G4KM_OpticalEqRhs(theField,
199 mass = G4SigmaPlus::SigmaPlus()->GetPDGMass
200 opticalCoeff =
201 (*theFieldMap)[G4SigmaPlus::SigmaPlus
202 opticalEq->SetFactor(mass,opticalCoeff);
203 (*theEquationMap)[G4SigmaPlus::SigmaPlus()-
204
205 opticalEq = new G4KM_OpticalEqRhs(theField,
206 mass = G4SigmaMinus::SigmaMinus()->GetPDGMa
207 opticalCoeff =
208 (*theFieldMap)[G4SigmaMinus::SigmaMin
209 opticalEq->SetFactor(mass,opticalCoeff);
210 (*theEquationMap)[G4SigmaMinus::SigmaMinus(
211
212 opticalEq = new G4KM_OpticalEqRhs(theField,
213 mass = G4SigmaZero::SigmaZero()->GetPDGMass
214 opticalCoeff =
215 (*theFieldMap)[G4SigmaZero::SigmaZero
216 opticalEq->SetFactor(mass,opticalCoeff);
217 (*theEquationMap)[G4SigmaZero::SigmaZero()-
218 }
219
220
221 //#define debug_1_RKPropagation 1
222 //--------------------------------------------
223 void G4RKPropagation::Transport(G4KineticTrack
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();
235 {
236 G4double currTimeStep = timeStep;
237 G4KineticTrack * kt = *i;
238 G4int encoding = kt->GetDefinition()->Ge
239
240 std::map <G4int, G4VNuclearField*, std::
241
242 G4VNuclearField* currentField=0;
243 if ( fieldIter != theFieldMap->end() ) c
244
245 // debug
246 // if ( timeStep > 1e30 ) {
247 // G4cout << " Name :" << kt->GetDefini
248 // }
249
250 // Get the time of intersections with th
251 G4double t_enter, t_leave;
252 // if the particle does not intersecate
253 if(!GetSphereIntersectionTimes(kt, t_ent
254 {
255 kt->SetState(G4KineticTrack::miss_nuc
256 continue;
257 }
258
259
260 #ifdef debug_1_RKPropagation
261 G4cout <<" kt,timeStep, Intersection tim
262 <<kt<< " / state= " <<kt->GetState
263 #endif
264
265 // if the particle is already outside nu
266 // does happen for particles added as l
267 // if(t_leave < 0 )
268 // {
269 // throw G4HadronicException(__FI
270 // continue;
271 // }
272
273 // Apply a straight line propagation for
274 // not included in the model
275 if( ! currentField )
276 {
277 if(currTimeStep == DBL_MAX)currTimeSt
278 FreeTransport(kt, currTimeStep);
279 if ( currTimeStep >= t_leave )
280 {
281 if ( kt->GetState() == G4KineticTr
282 { kt->SetState(G4KineticTrack::gon
283 else
284 { kt->SetState(G4KineticTrack::mis
285 } else if (kt->GetState() == G4Kineti
286 kt->SetState(G4KineticTrack::insid
287 }
288
289 continue;
290 }
291
292 if(t_enter > 0) // the particle is out.
293 {
294 if(t_enter > currTimeStep) // the pa
295 {
296 FreeTransport(kt, currTimeStep);
297 continue;
298 }
299 else
300 {
301 FreeTransport(kt, t_enter); // g
302 currTimeStep -= t_enter;
303 t_leave -= t_enter; // time
304 // on the surface the particle loo
305 // G4double newE = mom.e()-(*theF
306 // GetField = Barrier + FermiP
307 G4double newE = kt->GetTrackingMom
308
309 if(newE <= kt->GetActualMass()) /
310 {
311 // FixMe: should be "pushed bac
312 // for the moment take it
313 FreeTransport(kt, 1.1*t_leave);
314 kt->SetState(G4KineticTrack::mi
315 // G4cout << "G4RKPropagatio
316 // G4cout << " enter nucleus
317 // G4cout << " the Field "<<
318 // G4cout << " the particl
319 continue;
320 }
321 //
322 G4double newP = std::sqrt(newE*new
323 G4LorentzVector new4Mom(newP*kt->G
324 G4ThreeVector transfer(kt->GetTrac
325 G4ThreeVector boost= transfer / st
326 new4Mom*=G4LorentzRotation(boost);
327 kt->SetTrackingMomentum(new4Mom);
328 kt->SetState(G4KineticTrack::insid
329
330 /*
331 G4cout <<" Enter Nucleus - E/Field/Sum: "
332 << (*theFieldMap)[encoding]->GetField
333 << kt->GetTrackingMomentum().e()-currentF
334 << G4endl
335 << " Barrier / field just inside nucleus
336 << (*theFieldMap)[encoding]->GetBarrier()
337 << (*theFieldMap)[encoding]->GetField(0.9
338 << G4endl;
339 */
340 }
341 }
342
343 // FixMe: should I add a control on theC
344 // Transport the particle into the nucle
345 // G4cerr << "RKPropagation t_leav
346 G4bool is_exiting=false;
347 if(currTimeStep > t_leave) // particle
348 {
349 currTimeStep = t_leave;
350 is_exiting=true;
351 }
352
353 #ifdef debug_1_RKPropagation
354 G4cerr << "RKPropagation is_exiting?, t_
355 G4cout << "RKPropagation Ekin, field, pr
356 << kt->GetTrackingMomentum().e() -
357 << kt->GetPosition()<<" "
358 << G4endl << currentField->GetFiel
359 << kt->GetProjectilePotential()<<
360 << kt->GetTrackingMomentum()
361 << G4endl;
362 #endif
363
364 G4LorentzVector momold=kt->GetTrackingMo
365 G4ThreeVector posold=kt->GetPosition();
366
367 // if (currentField->GetField(kt->Get
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() -
376 << G4endl << currentField->GetFiel
377 << kt->GetTrackingMomentum()
378 << G4endl
379 << "delta p " << momold-kt->GetTra
380 << "del pos " << posold-kt->GetPos
381 << G4endl;
382 #endif
383
384 // complete the transport
385 // FixMe: in some cases there could be a
386 // part to do still in the nucleu
387 // slope of potential
388 G4double t_in=-1, t_out=0; // set onto
389
390 // should go out, or are already out by
391 if(is_exiting ||
392 (GetSphereIntersectionTimes(kt, t_
393 {
394 if(t_in < 0 && t_out >= 0) //still
395 {
396 // transport free to a position th
397 // a new transportation and a new
398 G4ThreeVector savePos = kt->GetPos
399 FreeTransport(kt, t_out);
400 // and evaluate the right the ener
401 G4double newE=kt->GetTrackingMomen
402
403 // G4cout << " V pos/savePos << "
404 // << (*theFieldMap)[encoding]-
405 // << (*theFieldMap)[encoding]-
406 // << G4endl;
407
408 if ( std::abs(currentField->GetFie
409 std::abs(currentField->GetFi
410 { // FixMe GF: savePos/pos may be
411 // This wrongly adds
412 // this is done later
413 newE += currentField->GetField(
414 - currentField
415 }
416
417 // G4cout << " go border nuc
418
419 if(newE < kt->GetActualMass())
420 {
421 #ifdef debug_1_RKPropagation
422 G4cout << "RKPropagation-Transp
423 G4cout << " cannot leave nucleu
424 #endif
425 if (kt->GetDefinition() == G4Pr
426 kt->GetDefinition() == G4
427 kt->SetState(G4KineticTrack:
428 } else {
429 kt->SetState(G4KineticTrack:
430 }
431 continue; // the particle canno
432 }
433 G4double newP = std::sqrt(newE*new
434 G4LorentzVector new4Mom(newP*kt->G
435 G4ThreeVector transfer(kt->GetTrac
436 G4ThreeVector boost= transfer / st
437 new4Mom*=G4LorentzRotation(boost);
438 kt->SetTrackingMomentum(new4Mom);
439 }
440 // add the potential barrier
441 // FixMe the Coulomb field is not par
442 G4double newE = kt->GetTrackingMoment
443 if(newE < kt->GetActualMass())
444 { // the particle cannot exit the nu
445 #ifdef debug_1_RKPropagation
446 G4cout << " cannot leave nucleus,
447 #endif
448 if (kt->GetDefinition() == G4Proto
449 kt->GetDefinition() == G4Neu
450 kt->SetState(G4KineticTrack::ca
451 } else {
452 kt->SetState(G4KineticTrack::go
453 }
454 continue;
455 }
456 G4double newP = std::sqrt(newE*newE-
457 G4LorentzVector new4Mom(newP*kt->GetT
458 G4ThreeVector transfer(kt->GetTrackin
459 G4ThreeVector boost= transfer / std::
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::GetMomentumTran
472 //--------------------------------------------
473 {
474 return theMomentumTranfer;
475 }
476
477
478 //--------------------------------------------
479 G4bool G4RKPropagation::FieldTransport(G4Kinet
480 //--------------------------------------------
481 {
482 // G4cout <<"Stepper input"<<kt->GetTrac
483 // create the integrator stepper
484 // G4Mag_EqRhs * equation = mapIter->sec
485 G4Mag_EqRhs * equation = (*theEquationMap)[
486 G4MagIntegratorStepper * stepper = new G4Cl
487
488 // create the integrator driver
489 G4double hMin = 1.0e-25*second; // arbitr
490 G4MagInt_Driver * driver = new G4MagInt_Dri
491
492 // Temporary: use driver->AccurateAdvance()
493 // create the G4FieldTrack needed by Accura
494 G4double curveLength = 0;
495 G4FieldTrack track(kt->GetPosition(),
496 kt->GetTrackingMomentum().vect().unit
497 curveLength, // curvelength
498 kt->GetTrackingMomentum().e()-kt->Get
499 kt->GetActualMass(), // restmass
500 kt->GetTrackingMomentum().beta()*c_li
501 // integrate
502 G4double eps = 0.01;
503 // G4cout << "currTimeStep = " << currTi
504 if(!driver->AccurateAdvance(track, timeStep
505 { // cannot track this particle
506 #ifdef debug_1_RKPropagation
507 std::cerr << "G4RKPropagation::FieldTran
508 << G4endl << "position " << kt->Ge
509 <<G4endl << " timestep " <<timeSte
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]->GetFie
519 << mom.e()+(*theFieldMap)[encoding]->Ge
520 << G4endl;
521 */
522
523 // Correct for momentum ( thus energy) tran
524 G4ThreeVector MomentumTranfer = kt->GetTrac
525 G4ThreeVector boost= MomentumTranfer / std:
526
527 // update the kt
528 kt->SetPosition(track.GetPosition());
529 G4LorentzVector mom(track.GetMomentum(),std
530 mom *= G4LorentzRotation( boost );
531 theMomentumTranfer += ( kt->GetTrackingMome
532 kt->SetTrackingMomentum(mom);
533
534 // G4cout <<"Stepper output"<<kt<<" "<<k
535 /*
536 * G4ThreeVector MomentumTranfer2=kt->Get
537 * G4cout << " MomentumTransfer/corrected"
538 * << " " << MomentumTranfer2 <<
539 * << MomentumTranfer-MomentumTranfer2
540 * MomentumTranfer-MomentumTranfer2.mag
541 * G4cout <<" E/Field/Sum aft : " <<mo
542 * << " / " << (*theFieldMap)[en
543 * << mom.e()+(*theFieldMap)[encoding]
544 * << G4endl;
545 */
546
547 delete driver;
548 delete stepper;
549 return true;
550 }
551
552 //--------------------------------------------
553 G4bool G4RKPropagation::FreeTransport(G4Kineti
554 //--------------------------------------------
555 {
556 G4ThreeVector newpos = kt->GetPosition() +
557 timeStep*c_light/kt->GetTrackingMomen
558 kt->SetPosition(newpos);
559 return true;
560 }
561
562 /*
563 G4bool G4RKPropagation::WillBeCaptured(const G
564 {
565 G4double radius = theOuterRadius;
566
567 // evaluate the final energy. Il will be captu
568 G4ParticleDefinition * definition = kt->GetD
569 G4double mass = definition->GetPDGMass();
570 G4ThreeVector pos = kt->GetPosition();
571 G4LorentzVector mom = kt->GetTrackingMomentu
572 G4VNuclearField * field = (*theFieldMap)[def
573 G4ThreeVector newPos(0, 0, radius); // to ge
574
575 G4double newE = mom.e()+field->GetField(pos)
576
577 return ((newE < mass) ? false : true);
578 }
579 */
580
581
582
583 //--------------------------------------------
584 G4bool G4RKPropagation::GetSphereIntersectionT
585 //--------------------------------------
586 const G4ThreeVector & currentPos,
587 const G4LorentzVector & momentum,
588 G4double & t1, G4double & t2)
589 {
590 G4ThreeVector speed = momentum.vect()/momen
591 G4double scalarProd = currentPos.dot(speed)
592 G4double speedMag2 = speed.mag2();
593 G4double sqrtArg = scalarProd*scalarProd -
594 speedMag2*(currentPos.mag2()-radius*r
595 if(sqrtArg <= 0.) // particle will not inte
596 {
597 // G4cout << " GetSphereIntersection
598 return false;
599 }
600 t1 = (-scalarProd - std::sqrt(sqrtArg))/spe
601 t2 = (-scalarProd + std::sqrt(sqrtArg))/spe
602 return true;
603 }
604
605 //--------------------------------------------
606 G4bool G4RKPropagation::GetSphereIntersectionT
607 G4double & t1, G4double & t2)
608 {
609 G4double radius = theOuterRadius + 3*fermi;
610 G4ThreeVector speed = kt->GetTrackingMoment
611 G4double scalarProd = kt->GetPosition().dot
612 G4double speedMag2 = speed.mag2();
613 G4double sqrtArg = scalarProd*scalarProd -
614 speedMag2*(kt->GetPosition().mag2()-r
615 if(sqrtArg <= 0.) // particle will not inte
616 {
617 return false;
618 }
619 t1 = (-scalarProd - std::sqrt(sqrtArg))/spe
620 t2 = (-scalarProd + std::sqrt(sqrtArg))/spe
621 return true;
622 }
623
624 // Implementation methods
625
626 //--------------------------------------------
627 void G4RKPropagation::delete_FieldsAndMap(
628 //--------------------------------------
629 std::map <G4int, G4VNuclearField *, std:
630 {
631 if(aMap)
632 {
633 std::map <G4int, G4VNuclearField *, std:
634 for(cur = aMap->begin(); cur != aMap->en
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::les
647 {
648 if(aMap)
649 {
650 std::map <G4int, G4Mag_EqRhs *, std::les
651 for(cur = aMap->begin(); cur != aMap->en
652 delete (*cur).second;
653
654 aMap->clear();
655 delete aMap;
656 }
657 }
658