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