source: trunk/source/processes/parameterisation/src/G4FastStep.cc@ 1066

Last change on this file since 1066 was 1007, checked in by garnier, 17 years ago

update to geant4.9.2
File size: 18.7 KB
RevLine 
[819]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// $Id: G4FastStep.cc,v 1.16 2006/06/29 21:09:32 gunter Exp $
[1007]28// GEANT4 tag $Name: geant4-09-02 $
[819]29//
30//---------------------------------------------------------------
31//
32// G4FastStep.cc
33//
34// Description:
35// Encapsulates a G4ParticleChange and insure friendly interface
36// methods to manage the primary/secondaries final state for
37// Fast Simulation Models.
38//
39// History:
40// Oct 97: Verderi && MoraDeFreitas - First Implementation.
41// Apr 98: MoraDeFreitas - G4FastStep becomes the G4ParticleChange
42// for the Fast Simulation Process.
43//
44//---------------------------------------------------------------
45
46#include "G4FastStep.hh"
47#include "G4Track.hh"
48#include "G4Step.hh"
49#include "G4TrackFastVector.hh"
50#include "G4DynamicParticle.hh"
51
52void G4FastStep::Initialize(const G4FastTrack& fastTrack)
53{
54 // keeps the fastTrack reference
55 fFastTrack=&fastTrack;
56
57 // currentTrack will be used to Initialize the other data members
58 const G4Track& currentTrack = *(fFastTrack->GetPrimaryTrack());
59
60 // use base class's method at first
61 G4VParticleChange::Initialize(currentTrack);
62
63 // set Energy/Momentum etc. equal to those of the parent particle
64 const G4DynamicParticle* pParticle = currentTrack.GetDynamicParticle();
65 theEnergyChange = pParticle->GetKineticEnergy();
66 theMomentumChange = pParticle->GetMomentumDirection();
67 thePolarizationChange = pParticle->GetPolarization();
68 theProperTimeChange = pParticle->GetProperTime();
69
70 // set Position/Time etc. equal to those of the parent track
71 thePositionChange = currentTrack.GetPosition();
72 theTimeChange = currentTrack.GetGlobalTime();
73
74 // switch off stepping hit invokation by default:
75 theSteppingControlFlag = AvoidHitInvocation;
76
77 // event biasing weigth:
78 theWeightChange = currentTrack.GetWeight();
79}
80
81//----------------------------------------
82// -- Set the StopAndKilled signal
83// -- and put kinetic energy to 0.0. in the
84// -- G4ParticleChange.
85//----------------------------------------
86void G4FastStep::KillPrimaryTrack()
87{
88 SetPrimaryTrackFinalKineticEnergy(0.) ;
89 ProposeTrackStatus(fStopAndKill) ;
90}
91
92//--------------------
93//
94//--------------------
95void
96G4FastStep::
97ProposePrimaryTrackFinalPosition(const G4ThreeVector &position,
98 G4bool localCoordinates)
99{
100 // Compute the position coordinate in global
101 // reference system if needed ...
102 G4ThreeVector globalPosition = position;
103 if (localCoordinates)
104 globalPosition = fFastTrack->GetInverseAffineTransformation()->
105 TransformPoint(position);
106 // ...and feed the globalPosition:
107 thePositionChange = globalPosition;
108}
109
110void
111G4FastStep::
112SetPrimaryTrackFinalPosition(const G4ThreeVector &position,
113 G4bool localCoordinates)
114{
115 ProposePrimaryTrackFinalPosition(position, localCoordinates);
116}
117
118//--------------------
119//
120//--------------------
121void
122G4FastStep::
123ProposePrimaryTrackFinalMomentumDirection(const G4ThreeVector &momentum,
124 G4bool localCoordinates)
125{
126 // Compute the momentum in global reference
127 // system if needed ...
128 G4ThreeVector globalMomentum = momentum;
129 if (localCoordinates)
130 globalMomentum = fFastTrack->GetInverseAffineTransformation()->
131 TransformAxis(momentum);
132 // ...and feed the globalMomentum (ensuring unitarity)
133 SetMomentumChange(globalMomentum.unit());
134}
135
136void
137G4FastStep::
138SetPrimaryTrackFinalMomentum(const G4ThreeVector &momentum,
139 G4bool localCoordinates)
140{
141 ProposePrimaryTrackFinalMomentumDirection(momentum, localCoordinates);
142}
143
144//--------------------
145//
146//--------------------
147void
148G4FastStep::
149ProposePrimaryTrackFinalKineticEnergyAndDirection(G4double kineticEnergy,
150 const G4ThreeVector &direction,
151 G4bool localCoordinates)
152{
153 // Compute global direction if needed...
154 G4ThreeVector globalDirection = direction;
155 if (localCoordinates)
156 globalDirection =fFastTrack->GetInverseAffineTransformation()->
157 TransformAxis(direction);
158 // ...and feed the globalMomentum (ensuring unitarity)
159 SetMomentumChange(globalDirection.unit());
160 SetPrimaryTrackFinalKineticEnergy(kineticEnergy);
161}
162
163void
164G4FastStep::
165SetPrimaryTrackFinalKineticEnergyAndDirection(G4double kineticEnergy,
166 const G4ThreeVector &direction,
167 G4bool localCoordinates)
168{
169 ProposePrimaryTrackFinalKineticEnergyAndDirection(kineticEnergy, direction, localCoordinates);
170}
171
172//--------------------
173//
174//--------------------
175void
176G4FastStep::
177ProposePrimaryTrackFinalPolarization(const G4ThreeVector &polarization,
178 G4bool localCoordinates)
179{
180 // Compute polarization in global system if needed:
181 G4ThreeVector globalPolarization(polarization);
182 if (localCoordinates)
183 globalPolarization = fFastTrack->GetInverseAffineTransformation()->
184 TransformAxis(globalPolarization);
185 // Feed the particle globalPolarization:
186 thePolarizationChange = globalPolarization;
187}
188
189void
190G4FastStep::
191SetPrimaryTrackFinalPolarization(const G4ThreeVector &polarization,
192 G4bool localCoordinates)
193{
194 ProposePrimaryTrackFinalPolarization(polarization, localCoordinates);
195}
196
197//--------------------
198//
199//--------------------
200G4Track* G4FastStep::
201CreateSecondaryTrack(const G4DynamicParticle& dynamics,
202 G4ThreeVector polarization,
203 G4ThreeVector position,
204 G4double time,
205 G4bool localCoordinates )
206{
207 G4DynamicParticle dummyDynamics(dynamics);
208
209 // ------------------------------------------
210 // Add the polarization to the dummyDynamics:
211 // ------------------------------------------
212 dummyDynamics.SetPolarization(polarization.x(),
213 polarization.y(),
214 polarization.z());
215
216 return CreateSecondaryTrack(dummyDynamics, position, time, localCoordinates);
217}
218
219//--------------------
220//
221//--------------------
222G4Track* G4FastStep::
223CreateSecondaryTrack(const G4DynamicParticle& dynamics,
224 G4ThreeVector position,
225 G4double time,
226 G4bool localCoordinates )
227{
228 // ----------------------------------------
229 // Quantities in global coordinates system.
230 //
231 // The allocated globalDynamics is deleted
232 // by the destructor of the G4Track.
233 // ----------------------------------------
234 G4DynamicParticle* globalDynamics =
235 new G4DynamicParticle(dynamics);
236 G4ThreeVector globalPosition(position);
237
238 // -----------------------------------
239 // Convert to global system if needed:
240 // -----------------------------------
241 if (localCoordinates)
242 {
243 // -- Momentum Direction:
244 globalDynamics->SetMomentumDirection(fFastTrack->
245 GetInverseAffineTransformation()->
246 TransformAxis(globalDynamics->
247 GetMomentumDirection()));
248 // -- Polarization:
249 G4ThreeVector globalPolarization;
250 globalPolarization = fFastTrack->GetInverseAffineTransformation()->
251 TransformAxis(globalDynamics->GetPolarization());
252 globalDynamics->SetPolarization(
253 globalPolarization.x(),
254 globalPolarization.y(),
255 globalPolarization.z()
256 );
257
258 // -- Position:
259 globalPosition = fFastTrack->GetInverseAffineTransformation()->
260 TransformPoint(globalPosition);
261 }
262
263 //-------------------------------------
264 // Create the G4Track of the secondary:
265 //-------------------------------------
266 G4Track* secondary = new G4Track(
267 globalDynamics,
268 time,
269 globalPosition
270 );
271
272 //-------------------------------
273 // and feed the changes:
274 //-------------------------------
275 AddSecondary(secondary);
276
277 //--------------------------------------
278 // returns the pointer on the secondary:
279 //--------------------------------------
280 return secondary;
281}
282
283// G4FastStep should never be Initialized in this way
284// but we must define it to avoid warnings.
285void G4FastStep::Initialize(const G4Track&)
286{
287 G4Exception("G4FastStep::Initialize()", "InvalidCondition", FatalException,
288 "G4FastStep can be initialised only through G4FastTrack!");
289}
290
291G4FastStep::G4FastStep()
292 : G4VParticleChange()
293{
294 if (verboseLevel>2)
295 {
296 G4cerr << "G4FastStep::G4FastStep() " << G4endl;
297 }
298}
299
300G4FastStep::~G4FastStep()
301{
302 if (verboseLevel>2)
303 {
304 G4cerr << "G4FastStep::~G4FastStep() " << G4endl;
305 }
306}
307
308// copy and assignment operators are implemented as "shallow copy"
309G4FastStep::G4FastStep(const G4FastStep &right)
310 : G4VParticleChange()
311{
312 *this = right;
313}
314
315
316G4FastStep & G4FastStep::operator=(const G4FastStep &right)
317{
318 if (this != &right)
319 {
320 G4VParticleChange::operator=(right);
321 theListOfSecondaries = right.theListOfSecondaries;
322 theSizeOftheListOfSecondaries = right.theSizeOftheListOfSecondaries;
323 theNumberOfSecondaries = right.theNumberOfSecondaries;
324 theStatusChange = right.theStatusChange;
325 theMomentumChange = right.theMomentumChange;
326 thePolarizationChange = right.thePolarizationChange;
327 thePositionChange = right.thePositionChange;
328 theTimeChange = right.theTimeChange;
329 theEnergyChange = right.theEnergyChange;
330 theTrueStepLength = right.theTrueStepLength;
331 theLocalEnergyDeposit = right.theLocalEnergyDeposit;
332 theSteppingControlFlag = right.theSteppingControlFlag;
333 theWeightChange = right.theWeightChange;
334 }
335 return *this;
336}
337
338
339
340
341
342G4bool G4FastStep::operator==(const G4FastStep &right) const
343{
344 return ((G4VParticleChange *)this == (G4VParticleChange *) &right);
345}
346
347G4bool G4FastStep::operator!=(const G4FastStep &right) const
348{
349 return ((G4VParticleChange *)this != (G4VParticleChange *) &right);
350}
351
352//----------------------------------------------------------------
353// methods for updating G4Step
354//
355
356G4Step* G4FastStep::UpdateStepForPostStep(G4Step* pStep)
357{
358 // A physics process always calculates the final state of the particle
359
360 // Take note that the return type of GetMomentumChange is a
361 // pointer to G4ParticleMometum. Also it is a normalized
362 // momentum vector.
363
364 // G4StepPoint* pPreStepPoint = pStep->GetPreStepPoint();
365 G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
366 G4Track* aTrack = pStep->GetTrack();
367 // G4double mass = aTrack->GetDynamicParticle()->GetMass();
368
369 // update kinetic energy and momentum direction
370 pPostStepPoint->SetMomentumDirection(theMomentumChange);
371 pPostStepPoint->SetKineticEnergy( theEnergyChange );
372
373 // update polarization
374 pPostStepPoint->SetPolarization( thePolarizationChange );
375
376 // update position and time
377 pPostStepPoint->SetPosition( thePositionChange );
378 pPostStepPoint->SetGlobalTime( theTimeChange );
379 pPostStepPoint->AddLocalTime( theTimeChange
380 - aTrack->GetGlobalTime());
381 pPostStepPoint->SetProperTime( theProperTimeChange );
382
383 // update weight
384 pPostStepPoint->SetWeight( theWeightChange );
385
386 if (debugFlag) CheckIt(*aTrack);
387
388
389 // Update the G4Step specific attributes
390 return UpdateStepInfo(pStep);
391}
392
393G4Step* G4FastStep::UpdateStepForAtRest(G4Step* pStep)
394{
395 // A physics process always calculates the final state of the particle
396
397 // G4StepPoint* pPreStepPoint = pStep->GetPreStepPoint();
398 G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
399 G4Track* aTrack = pStep->GetTrack();
400 // G4double mass = aTrack->GetDynamicParticle()->GetMass();
401
402 // update kinetic energy and momentum direction
403 pPostStepPoint->SetMomentumDirection(theMomentumChange);
404 pPostStepPoint->SetKineticEnergy( theEnergyChange );
405
406 // update polarization
407 pPostStepPoint->SetPolarization( thePolarizationChange );
408
409 // update position and time
410 pPostStepPoint->SetPosition( thePositionChange );
411 pPostStepPoint->SetGlobalTime( theTimeChange );
412 pPostStepPoint->AddLocalTime( theTimeChange
413 - aTrack->GetGlobalTime());
414 pPostStepPoint->SetProperTime( theProperTimeChange );
415
416 // update weight
417 pPostStepPoint->SetWeight( theWeightChange );
418
419 if (debugFlag) CheckIt(*aTrack);
420
421 // Update the G4Step specific attributes
422 return UpdateStepInfo(pStep);
423}
424
425//----------------------------------------------------------------
426// methods for printing messages
427//
428
429void G4FastStep::DumpInfo() const
430{
431// use base-class DumpInfo
432 G4VParticleChange::DumpInfo();
433
434 G4cout.precision(3);
435 G4cout << " Position - x (mm) : "
436 << std::setw(20) << thePositionChange.x()/mm
437 << G4endl;
438 G4cout << " Position - y (mm) : "
439 << std::setw(20) << thePositionChange.y()/mm
440 << G4endl;
441 G4cout << " Position - z (mm) : "
442 << std::setw(20) << thePositionChange.z()/mm
443 << G4endl;
444 G4cout << " Time (ns) : "
445 << std::setw(20) << theTimeChange/ns
446 << G4endl;
447 G4cout << " Proper Time (ns) : "
448 << std::setw(20) << theProperTimeChange/ns
449 << G4endl;
450 G4cout << " Momentum Direct - x : "
451 << std::setw(20) << theMomentumChange.x()
452 << G4endl;
453 G4cout << " Momentum Direct - y : "
454 << std::setw(20) << theMomentumChange.y()
455 << G4endl;
456 G4cout << " Momentum Direct - z : "
457 << std::setw(20) << theMomentumChange.z()
458 << G4endl;
459 G4cout << " Kinetic Energy (MeV): "
460 << std::setw(20) << theEnergyChange/MeV
461 << G4endl;
462 G4cout << " Polarization - x : "
463 << std::setw(20) << thePolarizationChange.x()
464 << G4endl;
465 G4cout << " Polarization - y : "
466 << std::setw(20) << thePolarizationChange.y()
467 << G4endl;
468 G4cout << " Polarization - z : "
469 << std::setw(20) << thePolarizationChange.z()
470 << G4endl;
471}
472
473G4bool G4FastStep::CheckIt(const G4Track& aTrack)
474{
475 //
476 // In the G4FastStep::CheckIt
477 // We only check a bit
478 //
479 // If the user violates the energy,
480 // We don't care, we agree.
481 //
482 // But for theMomentumDirectionChange,
483 // We do pay attention.
484 // And if too large is its range,
485 // We issue an Exception.
486 //
487 //
488 // It means, the G4FastStep::CheckIt issues an exception only for the
489 // theMomentumDirectionChange which should be an unit vector
490 // and it corrects it because it could cause problems for the ulterior
491 // tracking.For the rest, only warning are issued.
492
493 G4bool itsOK = true;
494 G4bool exitWithError = false;
495 G4double accuracy;
496
497 // Energy should not be larger than the initial value
498 accuracy = ( theEnergyChange - aTrack.GetKineticEnergy())/MeV;
499 if (accuracy > GetAccuracyForWarning()) {
500 G4cout << "ERROR - G4FastStep::CheckIt()" << G4endl
501 << " The energy becomes larger than the initial value !!"
502 << G4endl
503 << " Difference: " << accuracy << "[MeV] " << G4endl;
504 itsOK = false;
505 if (accuracy > GetAccuracyForException()) {exitWithError = true;}
506 }
507
508 G4bool itsOKforMomentum = true;
509 if ( theEnergyChange >0.) {
510 accuracy = std::abs(theMomentumChange.mag2()-1.0);
511 if (accuracy > GetAccuracyForWarning()) {
512 G4cout << "ERROR - G4FastStep::CheckIt()" << G4endl
513 << " The Momentum Change is not unit vector !!" << G4endl
514 << " Difference: " << accuracy << G4endl;
515 itsOK = itsOKforMomentum = false;
516 if (accuracy > GetAccuracyForException()) {exitWithError = true;}
517 }
518 }
519
520 accuracy = (aTrack.GetGlobalTime()- theTimeChange)/ns;
521 if (accuracy > GetAccuracyForWarning()) {
522 G4cout << "ERROR - G4FastStep::CheckIt()" << G4endl
523 << " The global time goes back !!" << G4endl
524 << " Difference: " << accuracy << "[ns] " << G4endl;
525 itsOK = false;
526 }
527
528 accuracy = (aTrack.GetProperTime() - theProperTimeChange )/ns;
529 if (accuracy) {
530 G4cout << "ERROR - G4FastStep::CheckIt()" << G4endl
531 << " The proper time goes back !!" << G4endl
532 << " Difference: " << accuracy << "[ns] " << G4endl;
533 itsOK = false;
534 }
535
536 if (!itsOK) {
537 G4cout << "ERROR - G4FastStep::CheckIt() " << G4endl;
538 G4cout << " Pointer : " << this << G4endl ;
539 DumpInfo();
540 }
541
542 // Exit with error
543 if (exitWithError) {
544 G4Exception("G4FastStep::CheckIt()", "FatalError",
545 FatalException, "Error condition in G4ParticleChange.");
546 }
547
548 //correction for Momentum only.
549 if (!itsOKforMomentum) {
550 G4double vmag = theMomentumChange.mag();
551 theMomentumChange = (1./vmag)*theMomentumChange;
552 }
553
554 itsOK = (itsOK) && G4VParticleChange::CheckIt(aTrack);
555 return itsOK;
556}
Note: See TracBrowser for help on using the repository browser.