Geant4_10
G4ScoreSplittingProcess.cc
Go to the documentation of this file.
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: G4ScoreSplittingProcess.cc 68733 2013-04-05 09:45:28Z gcosmo $
28 //
29 
31 
32 #include "G4ios.hh"
33 #include "G4SystemOfUnits.hh"
34 #include "G4Step.hh"
35 #include "G4VTouchable.hh"
36 #include "G4VPhysicalVolume.hh"
37 #include "G4ParticleChange.hh"
40 #include "G4ParticleChange.hh"
41 #include "G4StepPoint.hh"
42 
43 #include "G4SDManager.hh"
44 #include "G4VSensitiveDetector.hh"
45 
46 #include "G4EnergySplitter.hh"
47 #include "G4TouchableHistory.hh"
48 
49 //--------------------------------
50 // Constructor with name and type:
51 //--------------------------------
53 G4ScoreSplittingProcess(const G4String& processName,G4ProcessType theType)
54  :G4VProcess(processName,theType),
55  fOldTouchableH(), fNewTouchableH(), fInitialTouchableH(), fFinalTouchableH()
56 {
57  pParticleChange = &xParticleChange;
58 
59  fSplitStep = new G4Step();
60  fSplitPreStepPoint = fSplitStep->GetPreStepPoint();
61  fSplitPostStepPoint = fSplitStep->GetPostStepPoint();
62 
63  if (verboseLevel>0)
64  {
65  G4cout << GetProcessName() << " is created " << G4endl;
66  }
67  fpEnergySplitter = new G4EnergySplitter();
68 }
69 
70 // -----------
71 // Destructor:
72 // -----------
74 {
75  delete fSplitStep;
76  delete fpEnergySplitter;
77 }
78 
79 //------------------------------------------------------
80 //
81 // StartTracking
82 //
83 //------------------------------------------------------
85 {
86 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
87 // Setup initial touchables for the first step
88 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
89  const G4Step* pStep= trk->GetStep();
90 
91  fOldTouchableH = trk->GetTouchableHandle();
92  *fSplitPreStepPoint = *(pStep->GetPreStepPoint()); // Best to copy, so as to initialise
93  fSplitPreStepPoint->SetTouchableHandle(fOldTouchableH);
94  fNewTouchableH = fOldTouchableH;
95  *fSplitPostStepPoint= *(pStep->GetPostStepPoint()); // Best to copy, so as to initialise
96  fSplitPostStepPoint->SetTouchableHandle(fNewTouchableH);
97 
99  fSplitPreStepPoint ->SetStepStatus(fUndefined);
100  fSplitPostStepPoint->SetStepStatus(fUndefined);
101 }
102 
103 
104 //----------------------------------------------------------
105 //
106 // PostStepGetPhysicalInteractionLength()
107 //
108 //----------------------------------------------------------
109 G4double
111  const G4Track& /*track*/,
112  G4double /*previousStepSize*/,
114 {
115  // This process must be invoked anyway to score the hit
116  // - to do the scoring if the current volume is a regular structure, or
117  // - else to toggle the flag so that the SteppingManager does the scoring.
118  *condition = StronglyForced;
119 
120  // Future optimisation: check whether in regular structure.
121  // If it is in regular structure, be StronglyForced
122  // If not in regular structure,
123  // ask to be called only if SteppingControl is AvoidHitInvocation
124  // in order to reset it to NormalCondition
125 
126  return DBL_MAX;
127 }
128 
129 //------------------------------------
130 //
131 // PostStepDoIt()
132 //
133 //------------------------------------
135  const G4Track& track,
136  const G4Step& step)
137 {
138  G4VPhysicalVolume* pCurrentVolume= track.GetVolume();
139  G4LogicalVolume* pLogicalVolume= pCurrentVolume->GetLogicalVolume();
140  G4VSensitiveDetector* ptrSD = pLogicalVolume->GetSensitiveDetector();
141 
142  pParticleChange->Initialize(track);
143  if( ( ! pCurrentVolume->IsRegularStructure() ) || ( !ptrSD )
145  // Set the flag to make sure that Stepping Manager does the scoring
147  } else {
148  G4ThreeVector preStepPosition, postStepPosition, direction, finalPostStepPosition;
150 
151  G4double totalEnergyDeposit= step.GetTotalEnergyDeposit();
152  G4StepStatus fullStepStatus= step.GetPostStepPoint()->GetStepStatus();
153 
154  CopyStepStart(step);
155  fSplitPreStepPoint->SetSensitiveDetector(ptrSD);
156  fOldTouchableH = fInitialTouchableH;
157  fNewTouchableH= fOldTouchableH;
158  *fSplitPostStepPoint= *(step.GetPreStepPoint());
159 
160  // Split the energy
161  // ----------------
162  G4int numberVoxelsInStep= fpEnergySplitter->SplitEnergyInVolumes( &step );
163 
164  preStepPosition= step.GetPreStepPoint()->GetPosition();
165  finalPostStepPosition= step.GetPostStepPoint()->GetPosition();
166  direction= (finalPostStepPosition - preStepPosition).unit();
167 
168  fFinalTouchableH= track.GetNextTouchableHandle();
169 
170  postStepPosition= preStepPosition;
171  // Loop over the sub-parts of this step
172  G4int iStep;
173  for ( iStep=0; iStep < numberVoxelsInStep; iStep++ ){
174  G4int idVoxel= -1; // Voxel ID
175  G4double stepLength=0.0, energyLoss= 0.0;
176 
177  *fSplitPreStepPoint = *fSplitPostStepPoint;
178  fOldTouchableH = fNewTouchableH;
179 
180  preStepPosition= postStepPosition;
181  fSplitPreStepPoint->SetPosition( preStepPosition );
182  fSplitPreStepPoint->SetTouchableHandle(fOldTouchableH);
183 
184  fpEnergySplitter->GetLengthAndEnergyDeposited( iStep, idVoxel, stepLength, energyLoss);
185 
186  // Correct the material, so that the track->GetMaterial gives correct answer
187  pLogicalVolume->SetMaterial( fpEnergySplitter->GetVoxelMaterial( iStep) ); // idVoxel) );
188 
189  postStepPosition= preStepPosition + stepLength * direction;
190  fSplitPostStepPoint->SetPosition(postStepPosition);
191 
192  // Load the Step with the new values
193  fSplitStep->SetStepLength(stepLength);
194  fSplitStep->SetTotalEnergyDeposit(energyLoss);
195  if( iStep < numberVoxelsInStep -1 ){
196  fSplitStep->GetPostStepPoint()->SetStepStatus( fGeomBoundary );
197  G4int nextVoxelId= -1;
198  fpEnergySplitter->GetVoxelID( iStep+1, nextVoxelId );
199 
200  // Create new "next" touchable for each section ??
201  G4VTouchable* fNewTouchablePtr=
202  CreateTouchableForSubStep( nextVoxelId, postStepPosition );
203  fNewTouchableH= G4TouchableHandle(fNewTouchablePtr);
204  fSplitPostStepPoint->SetTouchableHandle( fNewTouchableH );
205  } else {
206  fSplitStep->GetPostStepPoint()->SetStepStatus( fullStepStatus );
207  fSplitPostStepPoint->SetTouchableHandle( fFinalTouchableH );
208  }
209 
210 
211  // As first approximation, split the NIEL in the same fractions as the energy deposit
212  G4double eLossFraction;
213  eLossFraction= (totalEnergyDeposit>0.0) ? energyLoss / totalEnergyDeposit : 1.0 ;
214  fSplitStep->SetNonIonizingEnergyDeposit(step.GetNonIonizingEnergyDeposit()*eLossFraction);
215 
216  fSplitPostStepPoint->SetSensitiveDetector( ptrSD );
217 
218  // Call the Sensitive Detector
219  ptrSD->Hit(fSplitStep);
220 
221  if (verboseLevel>1) Verbose(step);
222  }
223  }
224 
225  // This must change the Stepping Control
226  return pParticleChange;
227 }
228 
230 G4ScoreSplittingProcess::CreateTouchableForSubStep( G4int newVoxelNum, G4ThreeVector )
231 {
232  // G4cout << " Creating touchable handle for voxel-no " << newVoxelNum << G4endl;
233 
234  G4TouchableHistory* oldTouchableHistory= dynamic_cast<G4TouchableHistory*>(fOldTouchableH());
236  GetNavigatorForTracking()->CreateTouchableHistory(oldTouchableHistory->GetHistory());
237 
238  // Change the history
239  G4NavigationHistory* ptrNavHistory= const_cast<G4NavigationHistory*>(ptrTouchableHistory->GetHistory());
240  G4VPhysicalVolume* curPhysicalVol= ptrNavHistory->GetTopVolume();
241 
242  EVolume curVolumeType= ptrNavHistory->GetTopVolumeType();
243  if( curVolumeType == kParameterised )
244  {
245  ptrNavHistory->BackLevel();
246  // G4VPVParameterised parameterisedPV= pNewMother
247  G4VPVParameterisation* curParamstn= curPhysicalVol->GetParameterisation();
248 
249  // From G4ParameterisedNavigation::IdentifyAndPlaceSolid() inline method
250  G4VSolid* sampleSolid = curParamstn->ComputeSolid(newVoxelNum, curPhysicalVol);
251  sampleSolid->ComputeDimensions(curParamstn, newVoxelNum, curPhysicalVol);
252  curParamstn->ComputeTransformation(newVoxelNum, curPhysicalVol);
253 
254  ptrNavHistory->NewLevel( curPhysicalVol, kParameterised, newVoxelNum );
255  }
256  else
257  {
258  G4cout << " Current volume type is not Parameterised. " << G4endl;
259  G4Exception("G4ScoreSplittingProcess::CreateTouchableForSubStep",
260  "ErrorRegularParamaterisation", JustWarning,
261  "Score Splitting Process is used for Regular Structure - but did not find one here.");
262  }
263  return ptrTouchableHistory;
264 }
265 
266 void G4ScoreSplittingProcess::CopyStepStart(const G4Step & step)
267 {
268  fSplitStep->SetTrack(step.GetTrack());
269  fSplitStep->SetStepLength(step.GetStepLength());
270  fSplitStep->SetTotalEnergyDeposit(step.GetTotalEnergyDeposit());
272  fSplitStep->SetControlFlag(step.GetControlFlag());
273 
274  *fSplitPreStepPoint = *(step.GetPreStepPoint());
275 
276  fInitialTouchableH= (step.GetPreStepPoint()) ->GetTouchableHandle();
277  fFinalTouchableH = (step.GetPostStepPoint())->GetTouchableHandle();
278 }
279 
281 {
282  G4cout << "In mass geometry ------------------------------------------------" << G4endl;
283  G4cout << " StepLength : " << step.GetStepLength()/mm << " TotalEnergyDeposit : "
284  << step.GetTotalEnergyDeposit()/MeV << G4endl;
285  G4cout << " PreStepPoint : "
286  << step.GetPreStepPoint()->GetPhysicalVolume()->GetName() << " - ";
289  else
290  { G4cout << "NoProcessAssigned"; }
291  G4cout << G4endl;
292  G4cout << " " << step.GetPreStepPoint()->GetPosition() << G4endl;
293  G4cout << " PostStepPoint : ";
294  if(step.GetPostStepPoint()->GetPhysicalVolume())
295  { G4cout << step.GetPostStepPoint()->GetPhysicalVolume()->GetName(); }
296  else
297  { G4cout << "OutOfWorld"; }
298  G4cout << " - ";
301  else
302  { G4cout << "NoProcessAssigned"; }
303  G4cout << G4endl;
304  G4cout << " " << step.GetPostStepPoint()->GetPosition() << G4endl;
305 
306  G4cout << "In ghost geometry ------------------------------------------------" << G4endl;
307  G4cout << " StepLength : " << fSplitStep->GetStepLength()/mm
308  << " TotalEnergyDeposit : "
309  << fSplitStep->GetTotalEnergyDeposit()/MeV << G4endl;
310  G4cout << " PreStepPoint : "
311  << fSplitStep->GetPreStepPoint()->GetPhysicalVolume()->GetName() << " ["
312  << fSplitStep->GetPreStepPoint()->GetTouchable()->GetReplicaNumber()
313  << " ]" << " - ";
314  if(fSplitStep->GetPreStepPoint()->GetProcessDefinedStep())
315  { G4cout << fSplitStep->GetPreStepPoint()->GetProcessDefinedStep()->GetProcessName(); }
316  else
317  { G4cout << "NoProcessAssigned"; }
318  G4cout << G4endl;
319  G4cout << " " << fSplitStep->GetPreStepPoint()->GetPosition() << G4endl;
320  G4cout << " PostStepPoint : ";
321  if(fSplitStep->GetPostStepPoint()->GetPhysicalVolume())
322  {
323  G4cout << fSplitStep->GetPostStepPoint()->GetPhysicalVolume()->GetName() << " ["
324  << fSplitStep->GetPostStepPoint()->GetTouchable()->GetReplicaNumber()
325  << " ]";
326  }
327  else
328  { G4cout << "OutOfWorld"; }
329  G4cout << " - ";
330  if(fSplitStep->GetPostStepPoint()->GetProcessDefinedStep())
331  { G4cout << fSplitStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName(); }
332  else
333  { G4cout << "NoProcessAssigned"; }
334  G4cout << G4endl;
335  G4cout << " " << fSplitStep->GetPostStepPoint()->GetPosition() << " == "
336  << fSplitStep->GetTrack()->GetMomentumDirection()
337  << G4endl;
338 
339 }
340 
341 
342 //----------------------------------------------------------
343 //
344 // AtRestGetPhysicalInteractionLength()
345 //
346 //----------------------------------------------------------
347 G4double
349  const G4Track& /*track*/,
351 {
352  *condition = NotForced; // Was Forced
353  return DBL_MAX;
354 }
355 
356 
357 //---------------------------------------
358 // AlongStepGetPhysicalInteractionLength
359 //---------------------------------------
361  const G4Track& , // track,
362  G4double , // previousStepSize,
363  G4double , // currentMinimumStep,
364  G4double& , // proposedSafety,
365  G4GPILSelection* selection)
366 {
367  *selection = NotCandidateForSelection;
368  return DBL_MAX;
369 }
370 
371 //------------------------------------
372 // AlongStepDoIt()
373 //------------------------------------
374 
376  const G4Track& track, const G4Step& )
377 {
378  // Dummy ParticleChange ie: does nothing
379  // Expecting G4Transportation to move the track
380  dummyParticleChange.Initialize(track);
381  return &dummyParticleChange;
382 }
383 
384 //------------------------------------
385 // AtRestDoIt()
386 //------------------------------------
388  const G4Track& track,
389  const G4Step&)
390 {
391  pParticleChange->Initialize(track);
392  return pParticleChange;
393 }
G4double condition(const G4ErrorSymMatrix &m)
void SetStepLength(G4double value)
virtual void Initialize(const G4Track &)
G4VPhysicalVolume * GetTopVolume() const
void SetPosition(const G4ThreeVector &aValue)
void SetTrack(G4Track *value)
G4int verboseLevel
Definition: G4VProcess.hh:368
void GetLengthAndEnergyDeposited(G4int stepNo, G4int &voxelID, G4double &stepLength, G4double &energyLoss)
G4double GetStepLength() const
virtual void ComputeDimensions(G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
Definition: G4VSolid.cc:137
virtual G4VSolid * ComputeSolid(const G4int, G4VPhysicalVolume *)
G4StepStatus GetStepStatus() const
G4double GetNonIonizingEnergyDeposit() const
G4ScoreSplittingProcess(const G4String &processName="ScoreSplittingProc", G4ProcessType theType=fParameterisation)
G4double AlongStepGetPhysicalInteractionLength(const G4Track &, G4double, G4double, G4double &, G4GPILSelection *)
const std::vector< std::pair< G4int, G4double > > & GetStepLengths()
G4SteppingControl GetControlFlag() const
const G4VTouchable * GetTouchable() const
const G4Step * GetStep() const
int G4int
Definition: G4Types.hh:78
void SetStepStatus(const G4StepStatus aValue)
G4StepStatus
Definition: G4StepStatus.hh:49
G4ReferenceCountedHandle< G4VTouchable > G4TouchableHandle
G4VParticleChange * AlongStepDoIt(const G4Track &, const G4Step &)
G4StepPoint * GetPreStepPoint() const
void SetSensitiveDetector(G4VSensitiveDetector *)
G4GLOB_DLL std::ostream G4cout
G4VPhysicalVolume * GetPhysicalVolume() const
const G4String & GetName() const
const G4ThreeVector & GetPosition() const
void SetControlFlag(G4SteppingControl StepControlFlag)
EVolume GetTopVolumeType() const
virtual G4VPVParameterisation * GetParameterisation() const =0
G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
const G4TouchableHandle & GetNextTouchableHandle() const
G4bool Hit(G4Step *aStep)
void Verbose(const G4Step &) const
Definition: G4Step.hh:76
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
const G4TouchableHandle & GetTouchableHandle() const
G4double GetTotalEnergyDeposit() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4TransportationManager * GetTransportationManager()
void NewLevel(G4VPhysicalVolume *pNewMother, EVolume vType=kNormal, G4int nReplica=-1)
const G4VProcess * GetProcessDefinedStep() const
const G4ThreeVector & GetMomentumDirection() const
G4LogicalVolume * GetLogicalVolume() const
void GetVoxelID(G4int stepNo, G4int &voxelID)
G4StepPoint * GetPostStepPoint() const
virtual void ComputeTransformation(const G4int, G4VPhysicalVolume *) const =0
G4VParticleChange * AtRestDoIt(const G4Track &, const G4Step &)
virtual G4int GetReplicaNumber(G4int depth=0) const
Definition: G4VTouchable.cc:58
void SetNonIonizingEnergyDeposit(G4double value)
G4VParticleChange * pParticleChange
Definition: G4VProcess.hh:283
void ProposeSteppingControl(G4SteppingControl StepControlFlag)
G4VPhysicalVolume * GetVolume() const
#define G4endl
Definition: G4ios.hh:61
void SetTotalEnergyDeposit(G4double value)
const G4NavigationHistory * GetHistory() const
G4double AtRestGetPhysicalInteractionLength(const G4Track &, G4ForceCondition *)
EVolume
Definition: geomdefs.hh:68
double G4double
Definition: G4Types.hh:76
G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)
static G4RegularNavigationHelper * Instance()
virtual G4bool IsRegularStructure() const =0
G4ForceCondition
G4Track * GetTrack() const
void SetTouchableHandle(const G4TouchableHandle &apValue)
void SetMaterial(G4Material *pMaterial)
#define DBL_MAX
Definition: templates.hh:83
G4VSensitiveDetector * GetSensitiveDetector() const
G4int SplitEnergyInVolumes(const G4Step *aStep)
G4GPILSelection
G4Material * GetVoxelMaterial(G4int stepNo)
G4ProcessType