source: trunk/source/processes/scoring/src/G4ScoreSplittingProcess.cc@ 1353

Last change on this file since 1353 was 1350, checked in by garnier, 15 years ago

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