| 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 // This example is provided by the Geant4-DNA collaboration
27 // Any report or published results obtained using the Geant4-DNA software
28 // shall cite the following Geant4-DNA collaboration publication:
29 // Med. Phys. 37 (2010) 4692-4708
30 // The Geant4-DNA web site is available at http://geant4-dna.org
31 //
32 // If you use this example, please cite the following publication:
33 // Rad. Prot. Dos. 133 (2009) 2-11
34 //
35 #include "SteppingAction.hh"
36 #include "RunAction.hh"
37 #include "DetectorConstruction.hh"
38 #include "G4AnalysisManager.hh"
39 #include "G4SystemOfUnits.hh"
40 #include "G4SteppingManager.hh"
41 #include "G4Alpha.hh"
42 #include "G4Electron.hh"
43
44 SteppingAction::SteppingAction(RunAction* run,const DetectorConstruction* det)
45 :fRun(run),fDetector(det)
46 {}
47
48 SteppingAction::~SteppingAction()
49 {}
50
51 void SteppingAction::UserSteppingAction(const G4Step* aStep)
52
53 {
54 // Analysis manager
55
56 G4AnalysisManager* man = G4AnalysisManager::Instance();
57
58 // Read phantom - Singleton
59
60 fMyCellParameterisation = CellParameterisation::Instance();
61
62 // Material : 1 is cytoplasm, 2 is nucleus
63
64 G4int matVoxelPRE = -1;
65 G4int matVoxelPOST = -1;
66
67 const G4StepPoint* preStep = aStep->GetPreStepPoint();
68 const G4StepPoint* postStep = aStep->GetPostStepPoint();
69 const G4Track* track = aStep->GetTrack();
70
71 const G4LogicalVolume* preVolume =
72 preStep->GetPhysicalVolume()->GetLogicalVolume();
73
74 const G4LogicalVolume* postVolume = nullptr;
75 if(postStep->GetPhysicalVolume())
76 {
77 postVolume = postStep->GetPhysicalVolume()->GetLogicalVolume();
78 }
79 const G4ParticleDefinition* particle =
80 track->GetDynamicParticle()->GetDefinition();
81
82 G4int preReplicaNumber = preStep->GetTouchableHandle()->GetReplicaNumber();
83 G4double edep = aStep->GetTotalEnergyDeposit();
84
85 if (preReplicaNumber>0)
86 {
87 matVoxelPRE = fMyCellParameterisation->GetTissueType(preReplicaNumber);
88 }
89
90 if(postVolume)
91 {
92 G4int postReplicaNumber = postStep->GetTouchableHandle()->GetReplicaNumber();
93 if (postReplicaNumber>0)
94 {
95 matVoxelPOST = fMyCellParameterisation->GetTissueType(postReplicaNumber);
96 }
97 }
98
99 // COUNT GAS DETECTOR HITS
100
101 if (particle == G4Alpha::AlphaDefinition())
102 {
103 if(postVolume == fDetector->GetLogicalIsobutane() &&
104 ((preVolume == fDetector->GetLogicalCollDetYoke())
105 ||
106 (preVolume == fDetector->GetLogicalCollDetGap4())
107 ||
108 (preVolume == fDetector->GetLogicalCollDetGap4())))
109 {
110 fRun->AddNbOfHitsGas();
111 }
112
113 // STOPPING POWER AND BEAM SPOT SIZE AT CELL ENTRANCE
114 if(preVolume == fDetector->GetLogicalPolyprop() &&
115 ( (postVolume == fDetector->GetLogicalKgm()) ||
116 (matVoxelPOST == 1)) )
117 {
118 G4double deltaE = preStep->GetKineticEnergy()
119 - postStep->GetKineticEnergy();
120 if(deltaE > 0.0)
121 {
122 //Fill ntupleid=1
123 man->FillNtupleDColumn(1,0,preStep->GetKineticEnergy()/keV);
124 man->FillNtupleDColumn(1,1,deltaE*micrometer/(keV*aStep->GetStepLength()));
125 man->AddNtupleRow(1);
126 }
127
128 // Average dE over step suggested by Michel Maire
129 G4ThreeVector coord1 = preStep->GetPosition();
130 const G4AffineTransform transformation1 =
131 preStep->GetTouchable()->GetHistory()->GetTopTransform();
132 G4ThreeVector localPosition1 = transformation1.TransformPoint(coord1);
133
134 G4ThreeVector coord2 = postStep->GetPosition();
135 const G4AffineTransform transformation2 =
136 postStep->GetTouchable()->GetHistory()->GetTopTransform();
137 G4ThreeVector localPosition2 = transformation2.TransformPoint(coord2);
138
139 G4ThreeVector localPosition =
140 localPosition1 + G4UniformRand()*(localPosition2-localPosition1);
141
142 //Fill ntupleid=2
143 man->FillNtupleDColumn(2,0,localPosition.x()/micrometer);
144 man->FillNtupleDColumn(2,1,localPosition.y()/micrometer);
145 man->AddNtupleRow(2);
146 }
147
148 // ALPHA RANGE
149 if (postStep->GetKineticEnergy() < eV &&
150 ( (matVoxelPOST==1) ||
151 (postVolume == fDetector->GetLogicalKgm()) ||
152 (matVoxelPOST==2) ) )
153 {
154 //Fill ntupleid=3
155 man->FillNtupleDColumn(3,0,postStep->GetPosition().x()/micrometer);
156 man->FillNtupleDColumn(3,1,postStep->GetPosition().y()/micrometer);
157 man->FillNtupleDColumn(3,2,postStep->GetPosition().z()/micrometer);
158 man->AddNtupleRow(3);
159 }
160
161 // TOTAL DOSE DEPOSIT AND DOSE DEPOSIT WITHIN A PHANTOM VOXEL
162 // FOR ALL PARTICLES
163 }
164
165 if (matVoxelPRE == 2)
166 {
167 G4double dose = (edep/joule)/(fRun->GetMassNucleus()/kg);
168 fRun->AddDoseN(dose);
169 fRun->AddDoseBox(preReplicaNumber, edep/eV);
170 }
171 else if (matVoxelPRE == 1)
172 {
173 G4double dose = (edep/joule)/(fRun->GetMassCytoplasm()/kg);
174 fRun->AddDoseC(dose);
175 fRun->AddDoseBox(preReplicaNumber, edep/eV);
176 }
177
178 // PROTECTION AGAINST POSSIBLE MSC LOOPS FOR e-
179
180 // if ( edep/MeV<1e-25 && particle == G4Electron::Electron())
181 // {
182 //aStep->GetTrack()->SetTrackStatus(fStopAndKill);
183 /*
184 G4cout << "*** Warning *** : msc loop for "
185 << track->GetDefinition()->GetParticleName()
186 << " in " <<
187 postPoint->GetTouchableHandle()->GetVolume()->GetName() << G4endl;
188 */
189 // }
190 }
191