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