Geant4  10.02.p02
HadrontherapyDetectorSD.cc
Go to the documentation of this file.
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 // Hadrontherapy advanced example for Geant4
27 // See more at: https://twiki.cern.ch/twiki/bin/view/Geant4/AdvancedExamplesHadrontherapy
28 
31 #include "G4Step.hh"
32 #include "G4VTouchable.hh"
33 #include "G4TouchableHistory.hh"
34 #include "G4SDManager.hh"
35 #include "HadrontherapyMatrix.hh"
36 #include "HadrontherapyLet.hh"
37 #include "G4Track.hh"
39 #include "G4SystemOfUnits.hh"
40 
44 {
45  G4String HCname;
46  collectionName.insert(HCname="HadrontherapyDetectorHitsCollection");
47  HitsCollection = NULL;
49 
50 }
51 
54 {
55 }
56 
59 {
60 
62  collectionName[0]);
63 }
64 
67 {
68 
69 
70  if (aStep -> GetPreStepPoint() -> GetPhysicalVolume() -> GetName() != "RODetectorZDivisionPhys") return false;
71 
72 
73  // Get kinetic energy
74  G4Track * theTrack = aStep -> GetTrack();
75  G4double kineticEnergy = theTrack -> GetKineticEnergy();
76 
77  G4ParticleDefinition *particleDef = theTrack -> GetDefinition();
78  //Get particle name
79  G4String particleName = particleDef -> GetParticleName();
80 
81  G4int pdg = particleDef ->GetPDGEncoding();
82 
83  // Get unique track_id (in an event)
84  G4int trackID = theTrack -> GetTrackID();
85 
86  G4double energyDeposit = aStep -> GetTotalEnergyDeposit();
87 
88  G4double DX = aStep -> GetStepLength();
89  G4int Z = particleDef-> GetAtomicNumber();
90  G4int A = particleDef-> GetAtomicMass();
91 
92 
93  // Read voxel indexes: i is the x index, k is the z index
94  const G4VTouchable* touchable = aStep->GetPreStepPoint()->GetTouchable();
95  G4int k = touchable->GetReplicaNumber(0);
96  G4int i = touchable->GetReplicaNumber(2);
97  G4int j = touchable->GetReplicaNumber(1);
98 
99 
100 
101 #ifdef G4ANALYSIS_USE_ROOT
103 #endif
104 
107 
108 
109  // ******************** let ***************************
110  if (let)
111  {
112  if ( !(Z==0 && A==1) ) // All but not neutrons
113  {
114  if( energyDeposit>0. && DX >0. )
115  {
116  if (pdg !=22) // not gamma
117  {
118  let -> FillEnergySpectrum(trackID, particleDef,energyDeposit, DX, i, j, k);
119  }
120  else if (kineticEnergy > 50.*keV) // gamma cut
121  {
122  let -> FillEnergySpectrum(trackID, particleDef,energyDeposit, DX, i , j, k);
123  }
124 
125  }
126  }
127  }
128 
129 
130 
131  // ******************** let ***************************
132 
133 
134  if (matrix)
135  {
136 
137  // Increment Fluences & accumulate energy spectra
138  // Hit voxels are marked with track_id throught hitTrack matrix
139  G4int* hitTrack = matrix -> GetHitTrack(i,j,k); // hitTrack MUST BE cleared at every eventAction!
140  if ( *hitTrack != trackID )
141  {
142  *hitTrack = trackID;
143  /*
144  * Fill FLUENCE data for every single nuclide
145  * Exclude e-, neutrons, gamma, ...
146  */
147  if ( Z >= 1) matrix -> Fill(trackID, particleDef, i, j, k, 0, true);
148 
149 
150 #ifdef G4ANALYSIS_USE_ROOT
151 /*
152  // Fragments kinetic energy (ntuple)
153  if (trackID !=1 && Z>=1)
154  {
155  // First step kinetic energy for every fragment
156  analysis -> FillKineticFragmentTuple(i, j, k, A, Z, kineticEnergy/MeV);
157  }
158  // Kinetic energy spectra for primary particles
159 
160  if ( trackID == 1 && i == 0)
161  {
162  // First step kinetic energy for primaries only
163  analysis -> FillKineticEnergyPrimaryNTuple(i, j, k, kineticEnergy/MeV);
164  }
165 */
166 #endif
167  }
168 
169  if(energyDeposit != 0)
170  {
171 /*
172  * This method will fill a dose matrix for every single nuclide.
173  * A method of the HadrontherapyMatrix class (StoreDoseFluenceAscii())
174  * is called automatically at the end of main (or via the macro command /analysis/writeDoseFile.
175  * It permits to store all dose/fluence data into a single plane ASCII file.
176 */
177  // if (A==1 && Z==1) // primary and sec. protons
178  if ( Z>=1 ) // exclude e-, neutrons, gamma, ...
179  matrix -> Fill(trackID, particleDef, i, j, k, energyDeposit);
180  /*
181  * Create a hit with the information of position is in the detector
182  */
184  detectorHit -> SetEdepAndPosition(i, j, k, energyDeposit);
185  HitsCollection -> insert(detectorHit);
186  }
187  }
188 
189 #ifdef G4ANALYSIS_USE_ROOT
190  if(energyDeposit != 0)
191  {
192  if(trackID != 1)
193  {
194  if (particleName == "proton")
195  analysis -> SecondaryProtonEnergyDeposit(i, energyDeposit/MeV);
196 
197  else if (particleName == "neutron")
198  analysis -> SecondaryNeutronEnergyDeposit(i, energyDeposit/MeV);
199 
200  else if (particleName == "alpha")
201  analysis -> SecondaryAlphaEnergyDeposit(i, energyDeposit/MeV);
202 
203  else if (particleName == "gamma")
204  analysis -> SecondaryGammaEnergyDeposit(i, energyDeposit/MeV);
205 
206  else if (particleName == "e-")
207  analysis -> SecondaryElectronEnergyDeposit(i, energyDeposit/MeV);
208 
209  else if (particleName == "triton")
210  analysis -> SecondaryTritonEnergyDeposit(i, energyDeposit/MeV);
211 
212  else if (particleName == "deuteron")
213  analysis -> SecondaryDeuteronEnergyDeposit(i, energyDeposit/MeV);
214 
215  else if (particleName == "pi+" || particleName == "pi-" || particleName == "pi0")
216  analysis -> SecondaryPionEnergyDeposit(i, energyDeposit/MeV);
217  }
218  }
219 #endif
220 
221  return true;
222 }
223 
226 {
227 
228  static G4int HCID = -1;
229  if(HCID < 0)
230  {
231  HCID = GetCollectionID(0);
232  }
233 
234  HCE -> AddHitsCollection(HCID,HitsCollection);
235 }
236 
static HadrontherapyAnalysisManager * GetInstance()
Get the pointer to the analysis manager.
static const double MeV
Definition: G4SIunits.hh:211
static HadrontherapyLet * GetInstance()
G4String name
Definition: TRTMaterials.hh:40
const G4VTouchable * GetTouchable() const
int G4int
Definition: G4Types.hh:78
virtual G4int GetCollectionID(G4int i)
G4StepPoint * GetPreStepPoint() const
static HadrontherapyMatrix * GetInstance()
double A(double temperature)
bool G4bool
Definition: G4Types.hh:79
G4THitsCollection< HadrontherapyDetectorHit > HadrontherapyDetectorHitsCollection
Definition: G4Step.hh:76
virtual G4int GetReplicaNumber(G4int depth=0) const
Definition: G4VTouchable.cc:58
void Initialize(G4HCofThisEvent *)
A class for connecting the simulation to an analysis package.
HadrontherapyDetectorHitsCollection * HitsCollection
void EndOfEvent(G4HCofThisEvent *HCE)
G4CollectionNameVector collectionName
static const double keV
Definition: G4SIunits.hh:213
G4bool ProcessHits(G4Step *aStep, G4TouchableHistory *ROhist)
double G4double
Definition: G4Types.hh:76