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

Last change on this file since 1355 was 1340, checked in by garnier, 15 years ago

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