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

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

import all except CVS
File size: 18.3 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// $Id: G4ParticleChange.cc,v 1.30 2007/03/30 01:03:53 kurasige Exp $
28// GEANT4 tag $Name:  $
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.) {
456    accuracy = std::abs(theMomentumDirectionChange.mag2()-1.0);
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.