| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // This example is provided by the Geant4-DNA 26 // This example is provided by the Geant4-DNA collaboration
27 // Any report or published results obtained us 27 // Any report or published results obtained using the Geant4-DNA software
28 // shall cite the following Geant4-DNA collabo 28 // shall cite the following Geant4-DNA collaboration publication:
29 // Med. Phys. 37 (2010) 4692-4708 29 // Med. Phys. 37 (2010) 4692-4708
30 // The Geant4-DNA web site is available at htt 30 // The Geant4-DNA web site is available at http://geant4-dna.org
31 // 31 //
32 // If you use this example, please cite the fo 32 // If you use this example, please cite the following publication:
33 // Rad. Prot. Dos. 133 (2009) 2-11 33 // Rad. Prot. Dos. 133 (2009) 2-11
34 // 34 //
35 #include "SteppingAction.hh" 35 #include "SteppingAction.hh"
36 #include "RunAction.hh" 36 #include "RunAction.hh"
37 #include "DetectorConstruction.hh" 37 #include "DetectorConstruction.hh"
38 #include "G4AnalysisManager.hh" 38 #include "G4AnalysisManager.hh"
39 #include "G4SystemOfUnits.hh" 39 #include "G4SystemOfUnits.hh"
40 #include "G4SteppingManager.hh" 40 #include "G4SteppingManager.hh"
41 #include "G4Alpha.hh" 41 #include "G4Alpha.hh"
42 #include "G4Electron.hh" 42 #include "G4Electron.hh"
43 43
44 SteppingAction::SteppingAction(RunAction* run, 44 SteppingAction::SteppingAction(RunAction* run,const DetectorConstruction* det)
45 :fRun(run),fDetector(det) 45 :fRun(run),fDetector(det)
46 {} 46 {}
47 47
48 SteppingAction::~SteppingAction() 48 SteppingAction::~SteppingAction()
49 {} 49 {}
50 50
51 void SteppingAction::UserSteppingAction(const 51 void SteppingAction::UserSteppingAction(const G4Step* aStep)
52 52
53 { 53 {
54 // Analysis manager 54 // Analysis manager
55 55
56 G4AnalysisManager* man = G4AnalysisManager:: 56 G4AnalysisManager* man = G4AnalysisManager::Instance();
57 57
58 // Read phantom - Singleton 58 // Read phantom - Singleton
59 59
60 fMyCellParameterisation = CellParameterisati 60 fMyCellParameterisation = CellParameterisation::Instance();
61 61
62 // Material : 1 is cytoplasm, 2 is nucleus 62 // Material : 1 is cytoplasm, 2 is nucleus
63 63
64 G4int matVoxelPRE = -1; 64 G4int matVoxelPRE = -1;
65 G4int matVoxelPOST = -1; 65 G4int matVoxelPOST = -1;
66 66
67 const G4StepPoint* preStep = aStep->GetPreSt 67 const G4StepPoint* preStep = aStep->GetPreStepPoint();
68 const G4StepPoint* postStep = aStep->GetPost 68 const G4StepPoint* postStep = aStep->GetPostStepPoint();
69 const G4Track* track = aStep->GetTrack(); 69 const G4Track* track = aStep->GetTrack();
70 70
71 const G4LogicalVolume* preVolume = 71 const G4LogicalVolume* preVolume =
72 preStep->GetPhysicalVolume()->GetLogicalVo 72 preStep->GetPhysicalVolume()->GetLogicalVolume();
73 73
74 const G4LogicalVolume* postVolume = nullptr; 74 const G4LogicalVolume* postVolume = nullptr;
75 if(postStep->GetPhysicalVolume()) 75 if(postStep->GetPhysicalVolume())
76 { 76 {
77 postVolume = postStep->GetPhysicalVolume() 77 postVolume = postStep->GetPhysicalVolume()->GetLogicalVolume();
78 } 78 }
79 const G4ParticleDefinition* particle = 79 const G4ParticleDefinition* particle =
80 track->GetDynamicParticle()->GetDefinition 80 track->GetDynamicParticle()->GetDefinition();
81 81
82 G4int preReplicaNumber = preStep->GetTouchab 82 G4int preReplicaNumber = preStep->GetTouchableHandle()->GetReplicaNumber();
83 G4double edep = aStep->GetTotalEnergyDeposit 83 G4double edep = aStep->GetTotalEnergyDeposit();
84 84
85 if (preReplicaNumber>0) 85 if (preReplicaNumber>0)
86 { 86 {
87 matVoxelPRE = fMyCellParameterisation->Ge 87 matVoxelPRE = fMyCellParameterisation->GetTissueType(preReplicaNumber);
88 } 88 }
89 89
90 if(postVolume) 90 if(postVolume)
91 { 91 {
92 G4int postReplicaNumber = postStep->GetTou 92 G4int postReplicaNumber = postStep->GetTouchableHandle()->GetReplicaNumber();
93 if (postReplicaNumber>0) 93 if (postReplicaNumber>0)
94 { 94 {
95 matVoxelPOST = fMyCellParameterisation-> 95 matVoxelPOST = fMyCellParameterisation->GetTissueType(postReplicaNumber);
96 } 96 }
97 } 97 }
98 98
99 // COUNT GAS DETECTOR HITS 99 // COUNT GAS DETECTOR HITS
100 100
101 if (particle == G4Alpha::AlphaDefinition()) 101 if (particle == G4Alpha::AlphaDefinition())
102 { 102 {
103 if(postVolume == fDetector->GetLogicalIsob 103 if(postVolume == fDetector->GetLogicalIsobutane() &&
104 ((preVolume == fDetector->GetLogicalCol 104 ((preVolume == fDetector->GetLogicalCollDetYoke())
105 || 105 ||
106 (preVolume == fDetector->GetLogicalCollDetGa 106 (preVolume == fDetector->GetLogicalCollDetGap4())
107 || 107 ||
108 (preVolume == fDetector->GetLogicalCollDetGa 108 (preVolume == fDetector->GetLogicalCollDetGap4())))
109 { 109 {
110 fRun->AddNbOfHitsGas(); 110 fRun->AddNbOfHitsGas();
111 } 111 }
112 112
113 // STOPPING POWER AND BEAM SPOT SIZE AT CE 113 // STOPPING POWER AND BEAM SPOT SIZE AT CELL ENTRANCE
114 if(preVolume == fDetector->GetLogicalPolyp 114 if(preVolume == fDetector->GetLogicalPolyprop() &&
115 ( (postVolume == fDetector->GetLogicalK 115 ( (postVolume == fDetector->GetLogicalKgm()) ||
116 (matVoxelPOST == 1)) ) 116 (matVoxelPOST == 1)) )
117 { 117 {
118 G4double deltaE = preStep->GetKineticEne 118 G4double deltaE = preStep->GetKineticEnergy()
119 - postStep->GetKineticEnergy(); 119 - postStep->GetKineticEnergy();
120 if(deltaE > 0.0) 120 if(deltaE > 0.0)
121 { 121 {
122 //Fill ntupleid=1 122 //Fill ntupleid=1
123 man->FillNtupleDColumn(1,0,preStep->GetKinet 123 man->FillNtupleDColumn(1,0,preStep->GetKineticEnergy()/keV);
124 man->FillNtupleDColumn(1,1,deltaE*micrometer 124 man->FillNtupleDColumn(1,1,deltaE*micrometer/(keV*aStep->GetStepLength()));
125 man->AddNtupleRow(1); 125 man->AddNtupleRow(1);
126 } 126 }
127 127
128 // Average dE over step suggested by Mic 128 // Average dE over step suggested by Michel Maire
129 G4ThreeVector coord1 = preStep->GetPosit 129 G4ThreeVector coord1 = preStep->GetPosition();
130 const G4AffineTransform transformation1 130 const G4AffineTransform transformation1 =
131 preStep->GetTouchable()->GetHistory()->GetTo 131 preStep->GetTouchable()->GetHistory()->GetTopTransform();
132 G4ThreeVector localPosition1 = transform 132 G4ThreeVector localPosition1 = transformation1.TransformPoint(coord1);
133 133
134 G4ThreeVector coord2 = postStep->GetPosi 134 G4ThreeVector coord2 = postStep->GetPosition();
135 const G4AffineTransform transformation2 135 const G4AffineTransform transformation2 =
136 postStep->GetTouchable()->GetHistory()->GetT 136 postStep->GetTouchable()->GetHistory()->GetTopTransform();
137 G4ThreeVector localPosition2 = transform 137 G4ThreeVector localPosition2 = transformation2.TransformPoint(coord2);
138 138
139 G4ThreeVector localPosition = 139 G4ThreeVector localPosition =
140 localPosition1 + G4UniformRand()*(localPosit 140 localPosition1 + G4UniformRand()*(localPosition2-localPosition1);
141 141
142 //Fill ntupleid=2 142 //Fill ntupleid=2
143 man->FillNtupleDColumn(2,0,localPosition 143 man->FillNtupleDColumn(2,0,localPosition.x()/micrometer);
144 man->FillNtupleDColumn(2,1,localPosition 144 man->FillNtupleDColumn(2,1,localPosition.y()/micrometer);
145 man->AddNtupleRow(2); 145 man->AddNtupleRow(2);
146 } 146 }
147 147
148 // ALPHA RANGE 148 // ALPHA RANGE
149 if (postStep->GetKineticEnergy() < eV && 149 if (postStep->GetKineticEnergy() < eV &&
150 ( (matVoxelPOST==1) || 150 ( (matVoxelPOST==1) ||
151 (postVolume == fDetector->GetLogicalKgm()) 151 (postVolume == fDetector->GetLogicalKgm()) ||
152 (matVoxelPOST==2) ) ) 152 (matVoxelPOST==2) ) )
153 { 153 {
154 //Fill ntupleid=3 154 //Fill ntupleid=3
155 man->FillNtupleDColumn(3,0,postStep->Get 155 man->FillNtupleDColumn(3,0,postStep->GetPosition().x()/micrometer);
156 man->FillNtupleDColumn(3,1,postStep->Get 156 man->FillNtupleDColumn(3,1,postStep->GetPosition().y()/micrometer);
157 man->FillNtupleDColumn(3,2,postStep->Get 157 man->FillNtupleDColumn(3,2,postStep->GetPosition().z()/micrometer);
158 man->AddNtupleRow(3); 158 man->AddNtupleRow(3);
159 } 159 }
160 160
161 // TOTAL DOSE DEPOSIT AND DOSE DEPOSIT WIT 161 // TOTAL DOSE DEPOSIT AND DOSE DEPOSIT WITHIN A PHANTOM VOXEL
162 // FOR ALL PARTICLES 162 // FOR ALL PARTICLES
163 } 163 }
164 164
165 if (matVoxelPRE == 2) 165 if (matVoxelPRE == 2)
166 { 166 {
167 G4double dose = (edep/joule)/(fRun->GetMas 167 G4double dose = (edep/joule)/(fRun->GetMassNucleus()/kg);
168 fRun->AddDoseN(dose); 168 fRun->AddDoseN(dose);
169 fRun->AddDoseBox(preReplicaNumber, edep/eV 169 fRun->AddDoseBox(preReplicaNumber, edep/eV);
170 } 170 }
171 else if (matVoxelPRE == 1) 171 else if (matVoxelPRE == 1)
172 { 172 {
173 G4double dose = (edep/joule)/(fRun->GetMas 173 G4double dose = (edep/joule)/(fRun->GetMassCytoplasm()/kg);
174 fRun->AddDoseC(dose); 174 fRun->AddDoseC(dose);
175 fRun->AddDoseBox(preReplicaNumber, edep/eV 175 fRun->AddDoseBox(preReplicaNumber, edep/eV);
176 } 176 }
177 177
178 // PROTECTION AGAINST POSSIBLE MSC LOOPS FOR 178 // PROTECTION AGAINST POSSIBLE MSC LOOPS FOR e-
179 179
180 // if ( edep/MeV<1e-25 && particle == G4Ele 180 // if ( edep/MeV<1e-25 && particle == G4Electron::Electron())
181 // { 181 // {
182 //aStep->GetTrack()->SetTrackStatus(fSto 182 //aStep->GetTrack()->SetTrackStatus(fStopAndKill);
183 /* 183 /*
184 G4cout << "*** Warning *** : msc loop for " 184 G4cout << "*** Warning *** : msc loop for "
185 << track->GetDefinition()->GetParticleName() 185 << track->GetDefinition()->GetParticleName()
186 << " in " << 186 << " in " <<
187 postPoint->GetTouchableHandle()->GetVolume() 187 postPoint->GetTouchableHandle()->GetVolume()->GetName() << G4endl;
188 */ 188 */
189 // } 189 // }
190 } 190 }
191 191