Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ML2ExpVoxels.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 <fstream>
44 
45 #include "ML2ExpVoxels.hh"
46 #include "G4SystemOfUnits.hh"
47 
49  G4String FileExperimentalData, G4String FileExperimentalDataOut):startCurve(0),
50  stopCurve(0),chi2Factor(0)
51 {
52  char a[10];
53  sprintf(a,"%d", seed);
54  seedName=(G4String)a;
55  saving_in_Selected_Voxels_every_events=saveEvents;
56  nRecycling=1;
57 
58 
59  fullFileOut=FileExperimentalDataOut+seedName+".m";
60  fullFileIn=FileExperimentalData;
61  nParticle=nTotalEvents=0;
62 
63 // define the extremes of global-volume containing all experimental voxels
64  G4double extr=100000000000.;
65  minZone.set(extr, extr, extr);
66  maxZone.set(-extr, -extr, -extr);
67  bHasExperimentalData=bData;
68 }
69 
71 {
72  delete [] startCurve;
73  delete [] stopCurve;
74  delete [] chi2Factor;
75  delete [] nVoxelsgeometry;
76 }
78 {
79  bHasExperimentalData=true;
80  std::ifstream in;
81 
82  Svoxel voxel; voxel.volumeId=0;
83  G4ThreeVector pos, halfSize;
84  G4double expDose;
85 
86  in.open(fullFileIn, std::ios::in);
87  if (in !=0)
88  {
89  G4String appo;
90  char a[1000];
91  in.getline(a,1000,'\n'); headerText1=(G4String)a;
92  in.getline(a,1000,'\n');
93  in >> nCurves;
94  startCurve=new G4int[nCurves];
95  stopCurve=new G4int[nCurves];
96  chi2Factor=new G4double[nCurves];
97  for (int i=0; i< nCurves; i++)
98  {
99  chi2Factor[i]=0.;
100  in >> startCurve[i];
101  in >> stopCurve[i];
102  in >> chi2Factor[i];
103  }
104  in.getline(a,1000,'\n');
105  in.getline(a,1000,'\n'); headerText2=(G4String)a;
106 
107  while (!in.eof())
108  {
109  in >> pos;
110  in >> halfSize;
111  if (bHasExperimentalData)
112  {
113  in >> expDose;
114  voxel.expDose=expDose/100.*(joule/kg); // input data in cGy
115  }
116  else
117  {
118  voxel.expDose=0.;
119  }
120  voxel.pos=pos;
121  voxel.halfSize=halfSize;
122  voxel.depEnergy=0.;
123  voxel.depEnergy2=0.;
124  voxel.nEvents=0;
125  voxel.depEnergyNorm=0.;
126  voxel.depEnergyNormError=0.;
127  voxels.push_back(voxel);
128 
129 // calculate the actual extremes of the global-volume containing all the experimental data
130  if (minZone.getX()>pos.getX()-halfSize.getX())
131  {minZone.setX(pos.getX()-halfSize.getX());}
132  if (maxZone.getX()<pos.getX()+halfSize.getX())
133  {maxZone.setX(pos.getX()+halfSize.getX());}
134 
135  if (minZone.getY()>pos.getY()-halfSize.getY())
136  {minZone.setY(pos.getY()-halfSize.getY());}
137  if (maxZone.getY()<pos.getY()+halfSize.getY())
138  {maxZone.setY(pos.getY()+halfSize.getY());}
139 
140  if (minZone.getZ()>pos.getZ()-halfSize.getZ())
141  {minZone.setZ(pos.getZ()-halfSize.getZ());}
142  if (maxZone.getZ()<pos.getZ()+halfSize.getZ())
143  {maxZone.setZ(pos.getZ()+halfSize.getZ());}
144 
145  }
146  }
147  else
148  {
149  std::cout << "ERROR I can't find the experimental data file" << G4endl;
150  return false;
151  }
152  in.close();
153 
154  nVoxelsgeometry=new G4int[(G4int) voxels.size()];
156 
157  return true;
158 }
160 {
161  for (int i=0; i<(int) voxels.size(); i++ )
162  {nVoxelsgeometry[i]=0;}
163 }
164 
165 void CML2ExpVoxels::add(const G4Step* aStep)
166 {
167  G4ThreeVector pos;
168  G4double depEnergy, density, voxelVolume;
169 
170  pos=aStep->GetPreStepPoint()->GetPosition();
171  depEnergy=aStep->GetTotalEnergyDeposit();
173 
174  G4ThreeVector minPos, maxPos;
175  G4bool newEvent=false;
176  G4double voxelMass, dose;
177 
178 // check if the event is inside the global-volume
179  if (minZone.getX()<= pos.getX() && pos.getX()<maxZone.getX() &&
180  minZone.getY()<= pos.getY() && pos.getY()<maxZone.getY() &&
181  minZone.getZ()<= pos.getZ() && pos.getZ()<maxZone.getZ())
182  {
183 // look for the voxel containing the event
184  for (int i=0; i<(int)voxels.size(); i++)
185  {
186  minPos=voxels[i].pos-voxels[i].halfSize;
187  maxPos=voxels[i].pos+voxels[i].halfSize;
188  if (minPos.getX()<= pos.getX() && pos.getX()<maxPos.getX() &&
189  minPos.getY()<= pos.getY() && pos.getY()<maxPos.getY() &&
190  minPos.getZ()<= pos.getZ() && pos.getZ()<maxPos.getZ())
191  {
192  voxelVolume=voxels[i].halfSize.getX()*voxels[i].halfSize.getY()*voxels[i].halfSize.getZ()*8.;
193  voxelMass=density*voxelVolume;
194 // calculate the dose
195  dose=depEnergy/(voxelMass*nRecycling);
196  voxels[i].nEvents++;
197  nVoxelsgeometry[i]++;
198  voxels[i].depEnergy+=dose;
199  voxels[i].depEnergy2+=dose*dose;
200  newEvent=true;
201 
202  Sparticle *particle=new Sparticle;
203  particle->dir=aStep->GetPreStepPoint()->GetMomentumDirection();
204  particle->pos=aStep->GetPreStepPoint()->GetPosition();
205  particle->kinEnergy=dose; // I use the same kinEnergy name to store the dose
206  particle->nPrimaryPart=-1;
207  particle->partPDGE=aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
208  particle->primaryParticlePDGE=-1;
209  particle->volumeId=i; // voxel index where the dose is accumulating
210  particle->volumeName="-1";
211  }
212  }
213  if (newEvent)
214  {
215 // save data
216  nTotalEvents++;
217  if (nTotalEvents%saving_in_Selected_Voxels_every_events==0 && nTotalEvents>0)
218  {
219  saveResults();
220  }
221  }
222  }
223 }
224 
226 {
227  int n=voxels[0].nEvents;
228  for (int i=0;i<(int)voxels.size();i++)
229  {
230  if (n>voxels[i].nEvents){n = voxels[i].nEvents;}
231  }
232  return n;
233 }
235 {
236  int n=nVoxelsgeometry[0];
237  for (int i=0;i<(int)voxels.size();i++)
238  {
239  if (n<nVoxelsgeometry[i]){n = nVoxelsgeometry[i];}
240  }
241  return n;
242 }
243 void CML2ExpVoxels::saveHeader()
244 {
245  std::ofstream out;
246  out.open(fullFileOut, std::ios::out);
247  out <<"% "<< headerText1 << G4endl;
248  out <<"n"<< seedName<<"="<< nCurves<<";" << G4endl;
249  out <<"fh"<< seedName<<"=["<< G4endl;
250  for (int i=0; i< nCurves; i++)
251  {
252  out << startCurve[i] << '\t';
253  out << stopCurve[i] << '\t';
254  out << chi2Factor[i]<< G4endl;
255  }
256  out << "];"<<G4endl;
257  out <<"% x [mm], y [mm], z [mm], Dx [mm], Dy [mm], Dz [mm], expDose [Gy], Calculated dose [Gy], Calculated dose2 [Gy^2], nEvents, normDose [Gy], normDoseError [Gy]";
258  out << G4endl;
259  out.close();
260 }
261 
263 {
264  if (nTotalEvents>0)
265  {
266  calculateNormalizedEd(voxels);
267  saveHeader();
268  std::ofstream out;
269  out.open(fullFileOut, std::ios::app);
270  out <<"d"<< seedName<<"=["<< G4endl;
271  for (int i=0; i<(int)voxels.size(); i++)
272  {
273  out <<voxels[i].pos.getX()/mm<<'\t'<<voxels[i].pos.getY()/mm<<'\t'<<voxels[i].pos.getZ()/mm<<'\t';
274  out <<voxels[i].halfSize.getX()/mm<<'\t'<<voxels[i].halfSize.getY()/mm<<'\t'<<voxels[i].halfSize.getZ()/mm<<'\t';
275  out <<voxels[i].expDose/(joule/kg)<<'\t'<<voxels[i].depEnergy/(joule/kg)<<'\t'<<voxels[i].depEnergy2/((joule/kg)*(joule/kg))<<'\t'<<voxels[i].nEvents<< '\t';
276  out <<voxels[i].depEnergyNorm/(joule/kg)<<'\t'<<voxels[i].depEnergyNormError/(joule/kg);
277  out << G4endl;
278  }
279  out << "];"<<G4endl;
280  out.close();
281  }
282 }
283 
284 void CML2ExpVoxels::calculateNormalizedEd(std::vector <Svoxel> &vox)
285 {
286  int i,j;
287  G4double cs, cc;
288  int n;
289  G4double d2, dd;
290  G4double v;
291  for (j=0;j<nCurves;j++)
292  {
293  cs=cc=0.;
294  for (i=startCurve[j]-1;i<stopCurve[j];i++)
295  {
296  cs+=vox[i].depEnergy*vox[i].expDose;
297  cc+=vox[i].depEnergy*vox[i].depEnergy;
298  }
299  if (cc>0.)
300  {
301  chi2Factor[j]=cs/cc;
302  }
303  for (i=startCurve[j]-1;i<stopCurve[j];i++)
304  {
305  dd=vox[i].depEnergy*vox[i].depEnergy;
306  d2=vox[i].depEnergy2;
307  n=vox[i].nEvents;
308  vox[i].depEnergyNorm=chi2Factor[j]*vox[i].depEnergy;
309  v=n*d2-dd;
310  if (v<0.){v=0;}
311  if (n>1){vox[i].depEnergyNormError=chi2Factor[j]*std::sqrt(v/(n-1));}
312  if (n==1){vox[i].depEnergyNormError=vox[i].depEnergyNorm;}
313  }
314  }
315 }