source: trunk/source/track/src/G4ParticleChange.cc@ 1273

Last change on this file since 1273 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 18.3 KB
RevLine 
[826]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//
[1196]27// $Id: G4ParticleChange.cc,v 1.31 2009/04/02 02:22:30 kurasige Exp $
28// GEANT4 tag $Name: geant4-09-03-cand-01 $
[826]29//
30//
31// --------------------------------------------------------------
32// GEANT 4 class implementation file
33//
34//
35//
36// ------------------------------------------------------------
37// Implemented for the new scheme 23 Mar. 1998 H.Kurahige
38// Change default debug flag to false 10 May. 1998 H.Kurahige
39// Add Track weight 12 Nov. 1998 H.Kurashige
40// Activate CheckIt method for VERBOSE mode 14 Dec. 1998 H.Kurashige
41// Modified CheckIt method for time 9 Feb. 1999 H.Kurashige
42// Rename SetXXX methods to ProposeXXX DynamicCharge Oct. 2005 H.Kurashige
43// Add get/ProposeMagneticMoment Mar 2007 H.Kurashige
44// --------------------------------------------------------------
45
46#include "G4ParticleChange.hh"
47#include "G4Track.hh"
48#include "G4Step.hh"
49#include "G4TrackFastVector.hh"
50#include "G4DynamicParticle.hh"
51#include "G4ExceptionSeverity.hh"
52
53
54G4ParticleChange::G4ParticleChange():G4VParticleChange()
55{
56 G4VParticleChange::SetSecondaryWeightByProcess(false);
57 G4VParticleChange::SetParentWeightByProcess(false);
58
59}
60
61G4ParticleChange::~G4ParticleChange()
62{
63#ifdef G4VERBOSE
64 if (verboseLevel>2) {
65 G4cout << "G4ParticleChange::~G4ParticleChange() " << G4endl;
66 }
67#endif
68}
69
70// copy constructor
71G4ParticleChange::G4ParticleChange(const G4ParticleChange &right): G4VParticleChange(right)
72{
73 if (verboseLevel>1) {
74 G4cout << "G4ParticleChange:: copy constructor is called " << G4endl;
75 }
76 theMomentumDirectionChange = right.theMomentumDirectionChange;
77 thePolarizationChange = right.thePolarizationChange;
78 thePositionChange = right.thePositionChange;
79 theTimeChange = right.theTimeChange;
80 theEnergyChange = right.theEnergyChange;
81 theMassChange = right.theMassChange;
82 theChargeChange = right.theChargeChange;
83 theMagneticMomentChange = right.theMagneticMomentChange;
84 // theWeightChange = right.theWeightChange;
85 theProperTimeChange = right.theProperTimeChange;
86}
87
88// assignemnt operator
89G4ParticleChange & G4ParticleChange::operator=(const G4ParticleChange &right)
90{
91 if (verboseLevel>1) {
92 G4cout << "G4ParticleChange:: assignment operator is called " << G4endl;
93 }
94 if (this != &right)
95 {
96 theListOfSecondaries = right.theListOfSecondaries;
97 theSizeOftheListOfSecondaries = right.theSizeOftheListOfSecondaries;
98 theNumberOfSecondaries = right.theNumberOfSecondaries;
99 theStatusChange = right.theStatusChange;
100
101 theMomentumDirectionChange = right.theMomentumDirectionChange;
102 thePolarizationChange = right.thePolarizationChange;
103 thePositionChange = right.thePositionChange;
104 theTimeChange = right.theTimeChange;
105 theEnergyChange = right.theEnergyChange;
106 theMassChange = right.theMassChange;
107 theChargeChange = right.theChargeChange;
108 theMagneticMomentChange = right.theMagneticMomentChange;
109 // theWeightChange = right.theWeightChange;
110
111 theTrueStepLength = right.theTrueStepLength;
112 theLocalEnergyDeposit = right.theLocalEnergyDeposit;
113 theSteppingControlFlag = right.theSteppingControlFlag;
114 }
115 return *this;
116}
117
118G4bool G4ParticleChange::operator==(const G4ParticleChange &right) const
119{
120 return ((G4VParticleChange *)this == (G4VParticleChange *) &right);
121}
122
123G4bool G4ParticleChange::operator!=(const G4ParticleChange &right) const
124{
125 return ((G4VParticleChange *)this != (G4VParticleChange *) &right);
126}
127
128
129//----------------------------------------------------------------
130// methods for handling secondaries
131//
132
133void G4ParticleChange::AddSecondary(G4DynamicParticle* aParticle,
134 G4bool IsGoodForTracking )
135{
136 // create track
137 G4Track* aTrack = new G4Track(aParticle, theTimeChange, thePositionChange);
138
139 // set IsGoodGorTrackingFlag
140 if (IsGoodForTracking) aTrack->SetGoodForTrackingFlag();
141
142 // Touchable handle is copied to keep the pointer
143 aTrack->SetTouchableHandle(theCurrentTrack->GetTouchableHandle());
144
145 // add a secondary
146 G4VParticleChange::AddSecondary(aTrack);
147}
148
149void G4ParticleChange::AddSecondary(G4DynamicParticle* aParticle,
150 G4ThreeVector newPosition,
151 G4bool IsGoodForTracking )
152{
153 // create track
154 G4Track* aTrack = new G4Track(aParticle, theTimeChange, newPosition);
155
156 // set IsGoodGorTrackingFlag
157 if (IsGoodForTracking) aTrack->SetGoodForTrackingFlag();
158
159 // Touchable is a temporary object, so you cannot keep the pointer
160 aTrack->SetTouchableHandle((G4VTouchable*)0);
161
162 // add a secondary
163 G4VParticleChange::AddSecondary(aTrack);
164}
165
166void G4ParticleChange::AddSecondary(G4DynamicParticle* aParticle,
167 G4double newTime,
168 G4bool IsGoodForTracking )
169{
170 // create track
171 G4Track* aTrack = new G4Track(aParticle, newTime, thePositionChange);
172
173 // set IsGoodGorTrackingFlag
174 if (IsGoodForTracking) aTrack->SetGoodForTrackingFlag();
175
176 // Touchable handle is copied to keep the pointer
177 aTrack->SetTouchableHandle(theCurrentTrack->GetTouchableHandle());
178
179 // add a secondary
180 G4VParticleChange::AddSecondary(aTrack);
181}
182
183void G4ParticleChange::AddSecondary(G4Track* aTrack)
184{
185 // add a secondary
186 G4VParticleChange::AddSecondary(aTrack);
187}
188
189//----------------------------------------------------------------
190// functions for Initialization
191//
192
193void G4ParticleChange::Initialize(const G4Track& track)
194{
195 // use base class's method at first
196 G4VParticleChange::Initialize(track);
197 theCurrentTrack= &track;
198
199 // set Energy/Momentum etc. equal to those of the parent particle
200 const G4DynamicParticle* pParticle = track.GetDynamicParticle();
201 theEnergyChange = pParticle->GetKineticEnergy();
202 theMomentumDirectionChange = pParticle->GetMomentumDirection();
203 thePolarizationChange = pParticle->GetPolarization();
204 theProperTimeChange = pParticle->GetProperTime();
205
206 // Set mass/charge/MagneticMoment of DynamicParticle
207 theMassChange = pParticle->GetMass();
208 theChargeChange = pParticle->GetCharge();
209 theMagneticMomentChange = pParticle->GetMagneticMoment();
210
211 // set Position/Time etc. equal to those of the parent track
212 thePositionChange = track.GetPosition();
213 theTimeChange = track.GetGlobalTime();
214
215 // theWeightChange = track.GetWeight();
216}
217
218//----------------------------------------------------------------
219// methods for updating G4Step
220//
221
222G4Step* G4ParticleChange::UpdateStepForAlongStep(G4Step* pStep)
223{
224 // A physics process always calculates the final state of the
225 // particle relative to the initial state at the beginning
226 // of the Step, i.e., based on information of G4Track (or
227 // equivalently the PreStepPoint).
228 // So, the differences (delta) between these two states have to be
229 // calculated and be accumulated in PostStepPoint.
230
231 // Take note that the return type of GetMomentumDirectionChange is a
232 // pointer to G4ParticleMometum. Also it is a normalized
233 // momentum vector.
234
235 G4StepPoint* pPreStepPoint = pStep->GetPreStepPoint();
236 G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
237 G4double mass = theMassChange;
238
239 // Set Mass/Charge/MagneticMoment
240 pPostStepPoint->SetMass(theMassChange);
241 pPostStepPoint->SetCharge(theChargeChange);
242 pPostStepPoint->SetMagneticMoment(theMagneticMomentChange);
243
244 // calculate new kinetic energy
245 G4double energy = pPostStepPoint->GetKineticEnergy()
246 + (theEnergyChange - pPreStepPoint->GetKineticEnergy());
247
248 // update kinetic energy and momentum direction
249 if (energy > 0.0) {
250 // calculate new momentum
251 G4ThreeVector pMomentum = pPostStepPoint->GetMomentum()
252 + ( CalcMomentum(theEnergyChange, theMomentumDirectionChange, mass)
253 - pPreStepPoint->GetMomentum());
254 G4double tMomentum = pMomentum.mag();
255 G4ThreeVector direction(1.0,0.0,0.0);
256 if( tMomentum > 0. ){
257 G4double inv_Momentum= 1.0 / tMomentum;
258 direction= pMomentum * inv_Momentum;
259 }
260 pPostStepPoint->SetMomentumDirection(direction);
261 pPostStepPoint->SetKineticEnergy( energy );
262 } else {
263 // stop case
264 //pPostStepPoint->SetMomentumDirection(G4ThreeVector(1., 0., 0.));
265 pPostStepPoint->SetKineticEnergy(0.0);
266 }
267
268 // update polarization
269 pPostStepPoint->AddPolarization( thePolarizationChange
270 - pPreStepPoint->GetPolarization());
271
272 // update position and time
273 pPostStepPoint->AddPosition( thePositionChange
274 - pPreStepPoint->GetPosition() );
275 pPostStepPoint->AddGlobalTime( theTimeChange
276 - pPreStepPoint->GetGlobalTime());
277 pPostStepPoint->AddLocalTime( theTimeChange
278 - pPreStepPoint->GetGlobalTime());
279 pPostStepPoint->AddProperTime( theProperTimeChange
280 - pPreStepPoint->GetProperTime());
281
282 // update weight
283 if (!fSetParentWeightByProcess){
284 G4double newWeight= theParentWeight/(pPreStepPoint->GetWeight())
285 * (pPostStepPoint->GetWeight());
286 pPostStepPoint->SetWeight( newWeight );
287 }
288
289#ifdef G4VERBOSE
290 G4Track* aTrack = pStep->GetTrack();
291 if (debugFlag) CheckIt(*aTrack);
292#endif
293
294 // Update the G4Step specific attributes
295 return UpdateStepInfo(pStep);
296}
297
298G4Step* G4ParticleChange::UpdateStepForPostStep(G4Step* pStep)
299{
300 // A physics process always calculates the final state of the particle
301
302 // Take note that the return type of GetMomentumChange is a
303 // pointer to G4ParticleMometum. Also it is a normalized
304 // momentum vector.
305
306 G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
307 G4Track* aTrack = pStep->GetTrack();
308
309 // Set Mass/Charge
310 pPostStepPoint->SetMass(theMassChange);
311 pPostStepPoint->SetCharge(theChargeChange);
312 pPostStepPoint->SetMagneticMoment(theMagneticMomentChange);
313
314 // update kinetic energy and momentum direction
315 pPostStepPoint->SetMomentumDirection(theMomentumDirectionChange);
316 pPostStepPoint->SetKineticEnergy( theEnergyChange );
317
318 // update polarization
319 pPostStepPoint->SetPolarization( thePolarizationChange );
320
321 // update position and time
322 pPostStepPoint->SetPosition( thePositionChange );
323 pPostStepPoint->SetGlobalTime( theTimeChange );
324 pPostStepPoint->AddLocalTime( theTimeChange
325 - aTrack->GetGlobalTime());
326 pPostStepPoint->SetProperTime( theProperTimeChange );
327
328 // update weight
329 if (!fSetParentWeightByProcess){
330 pPostStepPoint->SetWeight( theParentWeight );
331 }
332
333
334#ifdef G4VERBOSE
335 if (debugFlag) CheckIt(*aTrack);
336#endif
337
338 // Update the G4Step specific attributes
339 return UpdateStepInfo(pStep);
340}
341
342
343G4Step* G4ParticleChange::UpdateStepForAtRest(G4Step* pStep)
344{
345 // A physics process always calculates the final state of the particle
346
347 G4StepPoint* pPostStepPoint = pStep->GetPostStepPoint();
348 G4Track* aTrack = pStep->GetTrack();
349
350 // Set Mass/Charge
351 pPostStepPoint->SetMass(theMassChange);
352 pPostStepPoint->SetCharge(theChargeChange);
353 pPostStepPoint->SetMagneticMoment(theMagneticMomentChange);
354
355 // update kinetic energy and momentum direction
356 pPostStepPoint->SetMomentumDirection(theMomentumDirectionChange);
357 pPostStepPoint->SetKineticEnergy( theEnergyChange );
358
359 // update polarization
360 pPostStepPoint->SetPolarization( thePolarizationChange );
361
362 // update position and time
363 pPostStepPoint->SetPosition( thePositionChange );
364 pPostStepPoint->SetGlobalTime( theTimeChange );
365 pPostStepPoint->AddLocalTime( theTimeChange
366 - aTrack->GetGlobalTime());
367 pPostStepPoint->SetProperTime( theProperTimeChange );
368
369 // update weight
370 if (!fSetParentWeightByProcess){
371 pPostStepPoint->SetWeight( theParentWeight );
372 }
373
374#ifdef G4VERBOSE
375 if (debugFlag) CheckIt(*aTrack);
376#endif
377
378 // Update the G4Step specific attributes
379 return UpdateStepInfo(pStep);
380}
381
382//----------------------------------------------------------------
383// methods for printing messages
384//
385
386void G4ParticleChange::DumpInfo() const
387{
388// use base-class DumpInfo
389 G4VParticleChange::DumpInfo();
390
391 G4cout.precision(3);
392
393 G4cout << " Mass (GeV) : "
394 << std::setw(20) << theMassChange/GeV
395 << G4endl;
396 G4cout << " Charge (eplus) : "
397 << std::setw(20) << theChargeChange/eplus
398 << G4endl;
399 G4cout << " MagneticMoment : "
400 << std::setw(20) << theMagneticMomentChange << G4endl;
401 G4cout << " : = " << std::setw(20)
402 << theMagneticMomentChange*2.*theMassChange/c_squared/eplus/hbar_Planck
403 << "*[e hbar]/[2 m]"
404 << G4endl;
405 G4cout << " Position - x (mm) : "
406 << std::setw(20) << thePositionChange.x()/mm
407 << G4endl;
408 G4cout << " Position - y (mm) : "
409 << std::setw(20) << thePositionChange.y()/mm
410 << G4endl;
411 G4cout << " Position - z (mm) : "
412 << std::setw(20) << thePositionChange.z()/mm
413 << G4endl;
414 G4cout << " Time (ns) : "
415 << std::setw(20) << theTimeChange/ns
416 << G4endl;
417 G4cout << " Proper Time (ns) : "
418 << std::setw(20) << theProperTimeChange/ns
419 << G4endl;
420 G4cout << " Momentum Direct - x : "
421 << std::setw(20) << theMomentumDirectionChange.x()
422 << G4endl;
423 G4cout << " Momentum Direct - y : "
424 << std::setw(20) << theMomentumDirectionChange.y()
425 << G4endl;
426 G4cout << " Momentum Direct - z : "
427 << std::setw(20) << theMomentumDirectionChange.z()
428 << G4endl;
429 G4cout << " Kinetic Energy (MeV): "
430 << std::setw(20) << theEnergyChange/MeV
431 << G4endl;
432 G4cout << " Polarization - x : "
433 << std::setw(20) << thePolarizationChange.x()
434 << G4endl;
435 G4cout << " Polarization - y : "
436 << std::setw(20) << thePolarizationChange.y()
437 << G4endl;
438 G4cout << " Polarization - z : "
439 << std::setw(20) << thePolarizationChange.z()
440 << G4endl;
441}
442
443G4bool G4ParticleChange::CheckIt(const G4Track& aTrack)
444{
445 G4bool exitWithError = false;
446 G4double accuracy;
447
448 // No check in case of "fStopAndKill"
449 if (GetTrackStatus() == fStopAndKill ) {
450 return G4VParticleChange::CheckIt(aTrack);
451 }
452
453 // MomentumDirection should be unit vector
454 G4bool itsOKforMomentum = true;
455 if ( theEnergyChange >0.) {
[1196]456 accuracy = std::fabs(theMomentumDirectionChange.mag2()-1.0);
[826]457 if (accuracy > accuracyForWarning) {
458#ifdef G4VERBOSE
459 G4cout << " G4ParticleChange::CheckIt : ";
460 G4cout << "the Momentum Change is not unit vector !!" << G4endl;
461 G4cout << " Difference: " << accuracy << G4endl;
462#endif
463 itsOKforMomentum = false;
464 if (accuracy > accuracyForException) exitWithError = true;
465 }
466 }
467
468 // Both global and proper time should not go back
469 G4bool itsOKforGlobalTime = true;
470 accuracy = (aTrack.GetGlobalTime()- theTimeChange)/ns;
471 if (accuracy > accuracyForWarning) {
472#ifdef G4VERBOSE
473 G4cout << " G4ParticleChange::CheckIt : ";
474 G4cout << "the global time goes back !!" << G4endl;
475 G4cout << " Difference: " << accuracy << "[ns] " <<G4endl;
476#endif
477 itsOKforGlobalTime = false;
478 if (accuracy > accuracyForException) exitWithError = true;
479 }
480
481 G4bool itsOKforProperTime = true;
482 accuracy = (aTrack.GetProperTime() - theProperTimeChange )/ns;
483 if (accuracy > accuracyForWarning) {
484#ifdef G4VERBOSE
485 G4cout << " G4ParticleChange::CheckIt : ";
486 G4cout << "the proper time goes back !!" << G4endl;
487 G4cout << " Difference: " << accuracy << "[ns] " <<G4endl;
488#endif
489 itsOKforProperTime = false;
490 if (accuracy > accuracyForException) exitWithError = true;
491 }
492
493 // Kinetic Energy should not be negative
494 G4bool itsOKforEnergy = true;
495 accuracy = -1.0*theEnergyChange/MeV;
496 if (accuracy > accuracyForWarning) {
497#ifdef G4VERBOSE
498 G4cout << " G4ParticleChange::CheckIt : ";
499 G4cout << "the kinetic energy is negative !!" << G4endl;
500 G4cout << " Difference: " << accuracy << "[MeV] " <<G4endl;
501#endif
502 itsOKforEnergy = false;
503 if (accuracy > accuracyForException) exitWithError = true;
504 }
505
506 G4bool itsOK = itsOKforMomentum && itsOKforEnergy && itsOKforProperTime && itsOKforGlobalTime;
507 // dump out information of this particle change
508#ifdef G4VERBOSE
509 if (!itsOK) {
510 G4cout << " G4ParticleChange::CheckIt " <<G4endl;
511 DumpInfo();
512 }
513#endif
514
515 // Exit with error
516 if (exitWithError) {
517 G4Exception("G4ParticleChange::CheckIt",
518 "200",
519 EventMustBeAborted,
520 "momentum, energy, and/or time was illegal");
521 }
522 //correction
523 if (!itsOKforMomentum) {
524 G4double vmag = theMomentumDirectionChange.mag();
525 theMomentumDirectionChange = (1./vmag)*theMomentumDirectionChange;
526 }
527 if (!itsOKforGlobalTime) {
528 theTimeChange = aTrack.GetGlobalTime();
529 }
530 if (!itsOKforProperTime) {
531 theProperTimeChange = aTrack.GetProperTime();
532 }
533 if (!itsOKforEnergy) {
534 theEnergyChange = 0.0;
535 }
536
537 itsOK = (itsOK) && G4VParticleChange::CheckIt(aTrack);
538 return itsOK;
539}
540
541
542
Note: See TracBrowser for help on using the repository browser.