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

Last change on this file was 1340, checked in by garnier, 14 years ago

update ti head
File size: 18.5 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.32 2010/07/21 09:30:15 gcosmo Exp $
28// GEANT4 tag $Name: track-V09-03-09 $
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
53G4ParticleChange::G4ParticleChange()
54  : G4VParticleChange(), theEnergyChange(0.), theTimeChange(0.),
55    theProperTimeChange(0.), theMassChange(0.), theChargeChange(0.),
56    theMagneticMomentChange(0.), theCurrentTrack(0)
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
392  G4int oldprc = G4cout.precision(3);
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;
442  G4cout.precision(oldprc);
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.) {
458    accuracy = std::fabs(theMomentumDirectionChange.mag2()-1.0);
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.