Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
HadrontherapyLet.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 // Hadrontherapy advanced example for Geant4
27 // See more at: https://twiki.cern.ch/twiki/bin/view/Geant4/AdvancedExamplesHadrontherapy
28 
30 #include "HadrontherapyLet.hh"
31 
32 #include "HadrontherapyMatrix.hh"
35 #include "HadrontherapyMatrix.hh"
36 #include "G4RunManager.hh"
37 #include "G4SystemOfUnits.hh"
38 #include <cmath>
39 
42 
44 {
45  if (instance) delete instance;
46  instance = new HadrontherapyLet(pDet);
47  return instance;
48 }
49 
51 {
52  return instance;
53 }
54 
56  :filename("Let.out"),matrix(0) // Default output filename
57 {
58 
60 
61  if (!matrix)
62  G4Exception("HadrontherapyLet::HadrontherapyLet",
63  "Hadrontherapy0005", FatalException,
64  "HadrontherapyMatrix not found. Firstly create an instance of it.");
65 
66  nVoxels = matrix -> GetNvoxel();
67 
68  numberOfVoxelAlongX = matrix -> GetNumberOfVoxelAlongX();
69  numberOfVoxelAlongY = matrix -> GetNumberOfVoxelAlongY();
70  numberOfVoxelAlongZ = matrix -> GetNumberOfVoxelAlongZ();
71 
73  pPGA = (HadrontherapyPrimaryGeneratorAction*)runManager -> GetUserPrimaryGeneratorAction();
74  // Pointer to the detector material
75  detectorMat = pDet -> GetDetectorLogicalVolume() -> GetMaterial();
76  density = detectorMat -> GetDensity();
77  // Instances for Total LET
78  totalLetD = new G4double[nVoxels];
80 
81 }
82 
84 {
85  Clear();
86  delete [] totalLetD;
87  delete [] DtotalLetD;
88 }
89 
90 // Fill energy spectrum for every voxel (local energy spectrum)
92 {
93  for(G4int v=0; v < nVoxels; v++) totalLetD[v] = DtotalLetD[v] = 0.;
94  Clear();
95 }
100 {
101  for (size_t i=0; i < ionLetStore.size(); i++)
102  {
103  delete [] ionLetStore[i].letDN;
104  delete [] ionLetStore[i].letDD;
105  }
106  ionLetStore.clear();
107 }
109  G4ParticleDefinition* particleDef,
110  /*G4double kinEnergy,*/
111  G4double DE,
112  G4double DX,
113  G4int i, G4int j, G4int k)
114 {
115  if (DE <= 0. || DX <=0.) return;
116  if (!doCalculation) return;
117  G4int Z = particleDef -> GetAtomicNumber();
118 
119 
120  G4int PDGencoding = particleDef -> GetPDGEncoding();
121  PDGencoding -= PDGencoding%10;
122 
123  G4int voxel = matrix -> Index(i,j,k);
124  // Total LET calculation...
125  totalLetD[voxel] += DE*(DE/DX);
126  DtotalLetD[voxel] += DE;
127  // Single ion LET
128  if (Z>=1)
129  {
130  // Search for already allocated data...
131  size_t l;
132  for (l=0; l < ionLetStore.size(); l++)
133  {
134  if (ionLetStore[l].PDGencoding == PDGencoding)
135  if ( ((trackID ==1) && (ionLetStore[l].isPrimary)) || ((trackID !=1) && (!ionLetStore[l].isPrimary)))
136  break;
137  }
138 
139  if (l == ionLetStore.size()) // Just another type of ion/particle for our store...
140  {
141 
142  G4int A = particleDef -> GetAtomicMass();
143 
144  G4String fullName = particleDef -> GetParticleName();
145  G4String name = fullName.substr (0, fullName.find("[") ); // cut excitation energy [x.y]
146 
147  ionLet ion =
148  {
149  (trackID == 1) ? true:false, // is it the primary particle?
150  PDGencoding,
151  fullName,
152  name,
153  Z,
154  A,
155  new G4double[nVoxels], // Let Dose Numerator
156  new G4double[nVoxels] // Let Dose Denominator
157  };
158 
159  // Initialize let
160  for(G4int v=0; v < nVoxels; v++) ion.letDN[v] = ion.letDD[v] = 0.;
161  ionLetStore.push_back(ion);
162  //G4cout << "Allocated LET data for " << ion.name << G4endl;
163 
164  }
165  ionLetStore[l].letDN[voxel] += DE*(DE/DX);
166  ionLetStore[l].letDD[voxel] += DE;
167  }
168 }
169 
170 
171 
172 
174 {
175  for(G4int v=0; v < nVoxels; v++) if (DtotalLetD[v]>0.) totalLetD[v] = totalLetD[v]/DtotalLetD[v];
176  // Sort ions by A and then by Z ...
177  std::sort(ionLetStore.begin(), ionLetStore.end());
178  // Compute Let Track and Let Dose for any single ion
179 
180  for(G4int v=0; v < nVoxels; v++)
181  {
182  for (size_t ion=0; ion < ionLetStore.size(); ion++)
183  {
184  if (ionLetStore[ion].letDD[v] >0.) ionLetStore[ion].letDN[v] = ionLetStore[ion].letDN[v] / ionLetStore[ion].letDD[v];
185 
186  }// end loop over ions
187  }
188 
189 }// end loop over voxels
190 
192 {
193 #define width 15L
194  if(ionLetStore.size())
195  {
196  ofs.open(filename, std::ios::out);
197  if (ofs.is_open())
198  {
199 
200  // Write the voxels index and the list of particles/ions
201  ofs << std::setprecision(6) << std::left <<
202  "i\tj\tk\t";
203  ofs << std::setw(width) << "LDT";
204  for (size_t l=0; l < ionLetStore.size(); l++)
205  {
206  G4String a = (ionLetStore[l].isPrimary) ? "_1":"";
207  ofs << std::setw(width) << ionLetStore[l].name + a ;
208  }
209  ofs << G4endl;
210 
211  // Write data
212 
213 
214  G4AnalysisManager* LetFragmentTuple = G4AnalysisManager::Instance();
215 
216  LetFragmentTuple->SetVerboseLevel(1);
217  LetFragmentTuple->SetFirstHistoId(1);
218  LetFragmentTuple->SetFirstNtupleId(1);
219  LetFragmentTuple ->OpenFile("Let");
220 
221 
222  LetFragmentTuple ->CreateNtuple("coordinate", "Let");
223 
224 
225  LetFragmentTuple ->CreateNtupleIColumn("i");//1
226  LetFragmentTuple ->CreateNtupleIColumn("j");//2
227  LetFragmentTuple ->CreateNtupleIColumn("k");//3
228  LetFragmentTuple ->CreateNtupleDColumn("TotalLetD");//4
229  LetFragmentTuple ->CreateNtupleIColumn("A");//5
230  LetFragmentTuple ->CreateNtupleIColumn("Z");//6
231  LetFragmentTuple ->CreateNtupleDColumn("IonLETD");//7
232  LetFragmentTuple ->FinishNtuple();
233 
234 
235  for(G4int i = 0; i < numberOfVoxelAlongX; i++)
236  for(G4int j = 0; j < numberOfVoxelAlongY; j++)
237  for(G4int k = 0; k < numberOfVoxelAlongZ; k++)
238  {
239  LetFragmentTuple->FillNtupleIColumn(1,0, i);
240  LetFragmentTuple->FillNtupleIColumn(1,1, j);
241  LetFragmentTuple->FillNtupleIColumn(1,2, k);
242 
243  G4int v = matrix -> Index(i, j, k);
244 
245  for (size_t l=0; l < ionLetStore.size(); l++)
246  {
247  // Write only not identically null data lines
248 
249 
250  if(ionLetStore[l].letDN)
251  {
252  ofs << G4endl;
253  ofs << i << '\t' << j << '\t' << k << '\t';
254  ofs << std::setw(width) << totalLetD[v]/(keV/um);
255 
256  LetFragmentTuple->FillNtupleDColumn(1,3, totalLetD[v]/(keV/um));
257 
258 
259  for (size_t ll=0; ll < ionLetStore.size(); ll++)
260  {
261  ofs << std::setw(width) << ionLetStore[ll].letDN[v]/(keV/um) ;
262 
263  LetFragmentTuple->FillNtupleIColumn(1,4, ionLetStore[ll].A);
264  LetFragmentTuple->FillNtupleIColumn(1,5, ionLetStore[ll].Z);
265 
266 
267  LetFragmentTuple->FillNtupleDColumn(1,6, ionLetStore[ll].letDN[v]/(keV/um));
268  LetFragmentTuple->AddNtupleRow(1);
269 
270 
271  }
272  break;
273  }
274  }
275  }
276  ofs.close();
277 
278  LetFragmentTuple->Write();
279  LetFragmentTuple->CloseFile();
280  }
281 
282  }
283 
284  }
285 
static G4RunManager * GetRunManager()
Definition: G4RunManager.cc:80
const XML_Char * name
Definition: expat.h:151
HadrontherapyPrimaryGeneratorAction * pPGA
G4Material * detectorMat
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
static constexpr double keV
Definition: G4SIunits.hh:216
static HadrontherapyMatrix * GetInstance()
#define G4endl
Definition: G4ios.hh:61
HadrontherapyMatrix * matrix
G4bool SetFirstHistoId(G4int firstId)
G4bool OpenFile(const G4String &fileName="")
G4bool FillNtupleIColumn(G4int id, G4int value)
std::ofstream ofs
static HadrontherapyLet * GetInstance()
static constexpr double um
Definition: G4SIunits.hh:113
G4int CreateNtupleIColumn(const G4String &name)
static G4bool doCalculation
G4bool FillNtupleDColumn(G4int id, G4double value)
Float_t Z
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
std::vector< ionLet > ionLetStore
#define width
G4double * letDN
void FillEnergySpectrum(G4int trackID, G4ParticleDefinition *particleDef, G4double DE, G4double DX, G4int i, G4int j, G4int k)
double A(double temperature)
HadrontherapyLet(HadrontherapyDetectorConstruction *)
G4double * letDD
G4int CreateNtupleDColumn(const G4String &name)
void SetVerboseLevel(G4int verboseLevel)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
int G4int
Definition: G4Types.hh:78
G4int CreateNtuple(const G4String &name, const G4String &title)
static HadrontherapyLet * instance
G4bool SetFirstNtupleId(G4int firstId)