Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ML2SDWithVoxels.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 // The code was written by :
27 // ^Claudio Andenna claudio.andenna@ispesl.it, claudio.andenna@iss.infn.it
28 // *Barbara Caccia barbara.caccia@iss.it
29 // with the support of Pablo Cirrone (LNS, INFN Catania Italy)
30 // with the contribute of Alessandro Occhigrossi*
31 //
32 // ^INAIL DIPIA - ex ISPESL and INFN Roma, gruppo collegato Sanità, Italy
33 // *Istituto Superiore di Sanità and INFN Roma, gruppo collegato Sanità, Italy
34 // Viale Regina Elena 299, 00161 Roma (Italy)
35 // tel (39) 06 49902246
36 // fax (39) 06 49387075
37 //
38 // more information:
39 // http://g4advancedexamples.lngs.infn.it/Examples/medical-linac
40 //
41 //*******************************************************//
42 
43 #include "ML2SDWithVoxels.hh"
44 #include "ML2ExpVoxels.hh"
45 #include "G4SystemOfUnits.hh"
46 
47 // *************************************************************************************
48 
50  G4int seed, G4String ROGOutFile,
51  G4bool bROG, G4ThreeVector ctr, G4ThreeVector hSiz,
52  G4int NumVX, G4int NumVY, G4int NumVZ)
53 : G4VSensitiveDetector(name),voxelsSum(0), voxelsSingle(0)
54 {
55  saving_in_ROG_Voxels_every_events=voxSave;
56  bSaveROG=bROG;
57  bActive=true;
58  nParticle=0;
59  nParticleValatile=0;
60  nTotalEvents=0;
61  nSingleTotalEvents=0;
62  density=1.;
63  voxelVolume=0.;
64  voxelMass=0.;
65  nRecycling=1;
66 
67  G4String seedName;
68  char a[10];
69  sprintf(a,"%d", seed);
70  seedName=(G4String)a;
71 
72  fullOutFileData=ROGOutFile+"_"+seedName+".txt";
73  fullOutFileDataSingle="";
74  if(bSaveROG)
75  {
76  centre=ctr;
77  halfSize=hSiz;
78  NumberOfVoxelsAlongX=NumVX;
79  NumberOfVoxelsAlongY=NumVY;
80  NumberOfVoxelsAlongZ=NumVZ;
81  halfXVoxelDimensionX=halfSize.getX()/NumberOfVoxelsAlongX;
82  halfXVoxelDimensionY=halfSize.getY()/NumberOfVoxelsAlongY;
83  halfXVoxelDimensionZ=halfSize.getZ()/NumberOfVoxelsAlongZ;
84 
85  voxelVolume=halfXVoxelDimensionX*halfXVoxelDimensionY*halfXVoxelDimensionZ*8.;
86 
87 // voxels to store and save the sum of the geometry configurations
88  voxelsSum=new Svoxel**[NumberOfVoxelsAlongX];
89  for (int ix=0; ix< NumberOfVoxelsAlongX; ix++)
90  {
91  voxelsSum[ix]=new Svoxel*[NumberOfVoxelsAlongY];
92  for (int iy=0; iy< NumberOfVoxelsAlongY; iy++)
93  {
94  voxelsSum[ix][iy]=new Svoxel[NumberOfVoxelsAlongZ];
95  for (int iz=0; iz< NumberOfVoxelsAlongZ; iz++)
96  {
97  voxelsSum[ix][iy][iz].volumeId=-1;
98  voxelsSum[ix][iy][iz].depEnergy=0.;
99  voxelsSum[ix][iy][iz].depEnergy2=0.;
100  voxelsSum[ix][iy][iz].depEnergyNorm=0.;
101  voxelsSum[ix][iy][iz].depEnergyNormError=0.;
102  voxelsSum[ix][iy][iz].expDose=0.;
103  voxelsSum[ix][iy][iz].halfSize.set(halfXVoxelDimensionX, halfXVoxelDimensionY, halfXVoxelDimensionZ);
104  voxelsSum[ix][iy][iz].pos.set(2.*(ix)*halfXVoxelDimensionX -halfSize.getX()+halfXVoxelDimensionX + centre.getX(),
105  2.*(iy)*halfXVoxelDimensionY -halfSize.getY()+halfXVoxelDimensionY + centre.getY(),
106  2.*(iz)*halfXVoxelDimensionZ -halfSize.getZ()+halfXVoxelDimensionZ + centre.getZ());
107  voxelsSum[ix][iy][iz].nEvents=0;
108  }
109  }
110  }
111 
112 
113 // voxels to store and save the single geometry configuration
114  voxelsSingle=new Svoxel**[NumberOfVoxelsAlongX];
115  for (int ix=0; ix< NumberOfVoxelsAlongX; ix++)
116  {
117  voxelsSingle[ix]=new Svoxel*[NumberOfVoxelsAlongY];
118  for (int iy=0; iy< NumberOfVoxelsAlongY; iy++)
119  {
120  voxelsSingle[ix][iy]=new Svoxel[NumberOfVoxelsAlongZ];
121  for (int iz=0; iz< NumberOfVoxelsAlongZ; iz++)
122  {
123  voxelsSingle[ix][iy][iz].volumeId=-1;
124  voxelsSingle[ix][iy][iz].depEnergy=0.;
125  voxelsSingle[ix][iy][iz].depEnergy2=0.;
126  voxelsSingle[ix][iy][iz].depEnergyNorm=0.;
127  voxelsSingle[ix][iy][iz].depEnergyNormError=0.;
128  voxelsSingle[ix][iy][iz].expDose=0.;
129  voxelsSingle[ix][iy][iz].halfSize.set(halfXVoxelDimensionX, halfXVoxelDimensionY, halfXVoxelDimensionZ);
130  voxelsSingle[ix][iy][iz].pos.set(2.*(ix)*halfXVoxelDimensionX -halfSize.getX()+halfXVoxelDimensionX + centre.getX(),
131  2.*(iy)*halfXVoxelDimensionY -halfSize.getY()+halfXVoxelDimensionY + centre.getY(),
132  2.*(iz)*halfXVoxelDimensionZ -halfSize.getZ()+halfXVoxelDimensionZ + centre.getZ());
133  voxelsSingle[ix][iy][iz].nEvents=0;
134  }
135  }
136  }
137  }
138 }
139 
141 {
142  if(bSaveROG)
143  {
144  delete [] voxelsSum;
145  delete [] voxelsSingle;
146  }
147 }
149 {
150  for (int ix=0; ix< NumberOfVoxelsAlongX; ix++)
151  {
152  for (int iy=0; iy< NumberOfVoxelsAlongY; iy++)
153  {
154  for (int iz=0; iz< NumberOfVoxelsAlongZ; iz++)
155  {
156  voxelsSingle[ix][iy][iz].volumeId=-1;
157  voxelsSingle[ix][iy][iz].depEnergy=0.;
158  voxelsSingle[ix][iy][iz].depEnergy2=0.;
159  voxelsSingle[ix][iy][iz].depEnergyNorm=0.;
160  voxelsSingle[ix][iy][iz].depEnergyNormError=0.;
161  voxelsSingle[ix][iy][iz].expDose=0.;
162  voxelsSingle[ix][iy][iz].halfSize.set(halfXVoxelDimensionX, halfXVoxelDimensionY, halfXVoxelDimensionZ);
163  voxelsSingle[ix][iy][iz].pos.set(2.*(ix)*halfXVoxelDimensionX -halfSize.getX()+halfXVoxelDimensionX + centre.getX(),
164  2.*(iy)*halfXVoxelDimensionY -halfSize.getY()+halfXVoxelDimensionY + centre.getY(),
165  2.*(iz)*halfXVoxelDimensionZ -halfSize.getZ()+halfXVoxelDimensionZ + centre.getZ());
166  voxelsSingle[ix][iy][iz].nEvents=0;
167  }
168  }
169  }
170  nSingleTotalEvents=0;
171 }
173 {
174 
175  if (bActive)
176  {
177  G4double energyDep = aStep->GetTotalEnergyDeposit();
178  if (bSaveROG && energyDep>0.)
179  {
180  G4int ix, iy, iz;
181  G4String volumeName;
182 
183  ix=ROHist->GetReplicaNumber(2);
184  iy=ROHist->GetReplicaNumber(0);
185  iz=ROHist->GetReplicaNumber(1);
186 
188 
189  voxelMass=voxelVolume*density;
190  energyDep/=voxelMass*nRecycling;
191  voxelsSum[ix][iy][iz].volumeId=getIdFromVolumeName(aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetMaterial()->GetName());
192  voxelsSum[ix][iy][iz].depEnergy+=energyDep;
193  voxelsSum[ix][iy][iz].depEnergy2+=energyDep*energyDep;
194  voxelsSum[ix][iy][iz].nEvents++;
195 
196  voxelsSingle[ix][iy][iz].volumeId=getIdFromVolumeName(aStep->GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetMaterial()->GetName());
197  voxelsSingle[ix][iy][iz].depEnergy+=energyDep;
198  voxelsSingle[ix][iy][iz].depEnergy2+=energyDep*energyDep;
199  voxelsSingle[ix][iy][iz].nEvents++;
200  nTotalEvents++;
201  nSingleTotalEvents++;
202  if (nTotalEvents%saving_in_ROG_Voxels_every_events==0 && nTotalEvents>0)
203  {
204  save();
205  }
206  }
207  }
208  return true;
209 }
210 G4int CML2SDWithVoxels::getIdFromVolumeName(G4String name)
211 {
212  for (int i=0; i<(int)volumeNameIdLink.size(); i++)
213  {
214  if (volumeNameIdLink[i].volumeName==name)
215  {
216  return volumeNameIdLink[i].volumeId;
217  break;
218  }
219  }
220  return -1;
221 }
223 {
224 std::cout<< "n. of events collected in the whole ROG phantom for all geometries: "<< nTotalEvents<< G4endl;
225 std::cout<< "n. of events collected in the whole ROG phantom for the current geometry: "<< nSingleTotalEvents<< G4endl;
226  if (nTotalEvents>0)
227  {saveData(fullOutFileData, voxelsSum);}
228  if (nSingleTotalEvents>0)
229  {saveData(fullOutFileDataSingle, voxelsSingle);}
230 }
231 void CML2SDWithVoxels::saveData(G4String Filename, Svoxel ***voxels)
232 {
233  std::ofstream out;
234  out.open(Filename, std::ios::out);
235  out << "Sensitive Detector-Voxels. Total number of events, [mm]->centreX centreY centreZ HalfSizeX HalfSizeY HalfSizeZ minX maxX, minY maxY, minZ maxZ, Dx, Dy, Dz, nX, nY, nZ: \n";
236  out <<nTotalEvents<<'\t';
237  out <<centre.getX()/mm << '\t' << centre.getY()/mm<< '\t'<< centre.getZ()/mm<<'\t';
238  out <<halfSize.getX()/mm << '\t' << halfSize.getY()/mm <<'\t'<< halfSize.getZ()/mm<<'\t';
239  out <<(centre.getX()-halfSize.getX())/mm<<'\t'<<(centre.getX()+halfSize.getX())/mm<<'\t';
240  out <<(centre.getY()-halfSize.getY())/mm<<'\t'<<(centre.getY()+halfSize.getY())/mm<<'\t';
241  out <<(centre.getZ()-halfSize.getZ())/mm<<'\t'<<(centre.getZ()+halfSize.getZ())/mm<<'\t';
242  out <<halfXVoxelDimensionX/mm<<'\t'<<halfXVoxelDimensionY/mm<<'\t'<<halfXVoxelDimensionZ/mm<<'\t';
243  out <<NumberOfVoxelsAlongX <<'\t'<<NumberOfVoxelsAlongY <<'\t'<<NumberOfVoxelsAlongZ <<'\n';
244  out << "Number of physical volumes: "<< volumeNameIdLink.size() << '\n';
245  for (int i=0; i<(int)volumeNameIdLink.size(); i++)
246  {
247  out << volumeNameIdLink[i].volumeName <<'\t'<< volumeNameIdLink[i].volumeId << G4endl;
248  }
249 
250 
251  out << "Phys Volume x [mm], y [mm], z [mm], ix, iy, iz, Dose [Gy], Dose2 [Gy^2], nEvents" << G4endl;
252 
253  for (int ix=0; ix< NumberOfVoxelsAlongX; ix++)
254  {
255  for (int iy=0; iy< NumberOfVoxelsAlongY; iy++)
256  {
257  for (int iz=0; iz< NumberOfVoxelsAlongZ; iz++)
258  {
259  if (voxels[ix][iy][iz].nEvents>0)
260  {
261  out << voxels[ix][iy][iz].volumeId << '\t';
262  out << voxels[ix][iy][iz].pos.getX()/mm << '\t';
263  out << voxels[ix][iy][iz].pos.getY()/mm << '\t';
264  out << voxels[ix][iy][iz].pos.getZ()/mm << '\t';
265  out << ix << '\t';
266  out << iy << '\t';
267  out << iz << '\t';
268  out << voxels[ix][iy][iz].depEnergy/(joule/kg) << '\t';
269  out << voxels[ix][iy][iz].depEnergy2/((joule/kg)*(joule/kg)) << '\t';
270  out << voxels[ix][iy][iz].nEvents << G4endl;
271  }
272  }
273  }
274  }
275  out.close();
276 
277 }
279 {
280  unsigned int ind = fullOutFileData.find(".txt");
281  G4String onlyName=fullOutFileData.substr( 0, ind);
282  if (val=="")
283  {
284  static unsigned int indGeom=0;
285  char cT[5];
286  sprintf(cT,"%d",indGeom);
287  fullOutFileDataSingle=onlyName+"Single_"+G4String(cT)+".txt";
288  indGeom++;
289  }
290  else
291  {
292  fullOutFileDataSingle=onlyName+val+".txt";
293  }
294 }
295 
void set(double x, double y, double z)
const XML_Char * name
Definition: expat.h:151
G4int nEvents
G4double depEnergyNorm
G4double depEnergy
static constexpr double mm
Definition: G4SIunits.hh:115
G4Material * GetMaterial() const
G4double depEnergy2
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
G4int volumeId
const G4String & GetName() const
Definition: G4Material.hh:178
G4double GetDensity() const
Definition: G4Material.hh:180
double getY() const
int G4int
Definition: G4Types.hh:78
CML2SDWithVoxels(G4String name, G4int saving_in_ROG_Voxels_every_events, G4int seed, G4String ROGOutFile, G4bool bSaveROG, G4ThreeVector centre, G4ThreeVector halfSize, G4int NumberOfVoxelsAlongX, G4int NumberOfVoxelsAlongY, G4int NumberOfVoxelsAlongZ)
G4double expDose
G4StepPoint * GetPreStepPoint() const
double getX() const
G4VPhysicalVolume * GetPhysicalVolume() const
G4ThreeVector pos
long seed
Definition: chem4.cc:68
bool G4bool
Definition: G4Types.hh:79
G4int GetReplicaNumber(G4int depth=0) const
G4double depEnergyNormError
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
Definition: G4Step.hh:76
static constexpr double kg
Definition: G4SIunits.hh:182
void setFullOutFileDataSingle(G4String val)
G4double GetTotalEnergyDeposit() const
G4bool ProcessHits(G4Step *aStep, G4TouchableHistory *ROHist)
G4LogicalVolume * GetLogicalVolume() const
static constexpr double joule
Definition: G4SIunits.hh:204
double getZ() const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4ThreeVector halfSize