| Geant4 Cross Reference |
1 {
2 ////////////////////////////////////////////
3 gROOT->Reset();
4
5 bool saveImage = true;
6 string imageFileName = "microdosimetricSpectra.png";
7
8 bool saveLinearEdep = true;
9 bool saveLogEdep = true;
10 bool saveydySpec = true;
11
12 ////////////////////////////////////////////
13 //Open File where data has been stored
14 ////////////////////////////////////////////
15 TFile* f = new TFile("radioprotection_t0.root");
16
17 TDirectory* dir = (TDirectory*)f->Get("radioprotection_ntuple");
18 TTree * nEdep = (TTree*)dir->Get("102");
19 Double_t edep;
20 nEdep->SetBranchAddress("edep", &edep);
21
22 double meanChordLength = 1; //change based on the thickness of the sensitive volume (make sure to apply a conversion factor if converting to a tissue equivalent material)
23 ////////////////////////////////////////////
24 //Put hits in microdosimeter into linear binned space histogram (typical of experimental setup)
25 ////////////////////////////////////////////
26 int numberLinearBins = 4096;
27 double minLinealEnergy = 0;
28 double maxLinealEnergy = 1000;
29 TH1D* histoLin = new TH1D("h1", "f(y)", numberLinearBins, minLinealEnergy, maxLinealEnergy);
30
31 TH1D* histoYFY = new TH1D("h2", "yf(y)", numberLinearBins, minLinealEnergy, maxLinealEnergy);
32 TH1D* histoDY = new TH1D("h3", "d(y)", numberLinearBins, minLinealEnergy, maxLinealEnergy);
33
34 ////////////////////////////////////////////
35 //Create log bin for microdosimetric plotting
36 ////////////////////////////////////////////
37 //Log bin stuff
38 int B = 20; //the number of increments in each decade
39 //eg if B = 10: 1 2 3 ... 9 10 20 30 ... 90 100 200 ...
40 double yMin = 0.01; // y = lineal energy
41 double yMax = 1000;
42
43 double range = yMax / yMin;
44 double fooRange = range;
45 double fooMag = 10; //any old silly no > 1
46 int order = 0;
47 while (fooMag > 1)
48 {
49 fooMag = fooRange / 10 ;
50 fooRange = fooRange / 10 ;
51 order = order + 1 ;
52 }
53
54 int binsPerDecade = B;
55 const int numberOfLogBins = binsPerDecade * order +1; //get the total number of bins
56
57 vector<double> xbins;
58 for (int b = 0; b < numberOfLogBins; b++)
59 {
60 xbins.push_back((yMin)*pow(10, ((b) / double(B))));
61 // cout << b << " " << xbins[b] << endl;
62 }
63
64 TH1D* histoLog = new TH1D("l1", "edepLin", (numberOfLogBins-1), &xbins[0]);
65
66
67 ////////////////////////////////////////////
68 //loop through the hits in the detector stored in the ntuple
69 ////////////////////////////////////////////
70 int numberHits = nEdep->GetEntries();
71 cout << "Number of hits in detector = " << numberHits << endl;
72
73 for (int i = 0; i < numberHits; i++)
74 {
75 nEdep->GetEntry(i);
76 if (edep <= 0)
77 {
78 cout << "Edep = 0" << endl;
79 }
80
81 histoLin->Fill(edep/meanChordLength);
82 histoLog->Fill(edep/meanChordLength);
83 }
84
85 ////////////////////////////////////////////
86 //Normalise the log histogram based on the bin width
87 ////////////////////////////////////////////
88 for (int bb = 1; bb <= numberOfLogBins-1; bb++)
89 {
90 histoLog->SetBinContent(bb, (histoLog->GetBinContent(bb) / histoLog->GetBinWidth(bb)));
91 }
92
93
94 ////////////////////////////////////////////
95 //---------Save values to file--------------
96 ////////////////////////////////////////////
97 ofstream outFile;
98 if (saveLinearEdep)
99 {
100 outFile.open("linEdepSpec.txt");
101 for (int i = 1 ; i <= numberLinearBins; i++)
102 {
103 outFile << histoLin->GetBinCenter(i) * meanChordLength << "\t" << histoLin->GetBinContent(i) << endl;
104 }
105 outFile.close();
106 }
107 if (saveLogEdep)
108 {
109 outFile.open("logEdepSpec.txt");
110 for (int bb = 1; bb <= numberOfLogBins-1; bb++)
111 {
112 outFile << histoLog->GetBinCenter(bb) * meanChordLength << "\t" << histoLog->GetBinContent(bb) << endl;
113 }
114 outFile.close();
115 }
116
117 ////////////////////////////////////////////
118 //----------Normalise Histo to get f(y)----
119 ////////////////////////////////////////////
120 double intOfHist = histoLin->Integral("width") ; //gets counts in each histogram and multiplies by bin width
121 histoLin->Scale(1./intOfHist);
122 //histoLin->Sumw2();
123
124 double intOfLogHist = histoLog->Integral("width") ;
125 histoLog->Scale(1./intOfLogHist);
126 //histoLog->Sumw2();
127
128
129 //////////////////////////////////////////////////////////////
130 //-------Calculate RBE using the modified MK model-----------
131 /////////////////////////////////////////////////////////////
132 //Technically the use of the modified MKM is applied in
133 //therapeutic ion beams and not used for radioprotection applications
134 //though it is included here for convience.
135 //NOTE: the form used here is for "heavy ions" and not applicable
136 //to the more dose depedent protons
137
138 //---------------------------------
139 //Constants
140 //---------------------------------
141 const double satPar = 150.; //saturation paramater obv. in keV/um
142 double satParSqr = satPar * satPar;
143 const double pi = 3.14159265359;
144 const double SiToMuscleCorrection = 0.57;
145
146 //Cell reltated numbers
147 double alpha0 = 0.13;//0.164;//0.13; //Gy^-1 //0.13 for carbon
148 double beta = 0.05; //Gy^-2
149 double rho = 1.; //g/cm^3
150 double rSV = 0.42; //um, the radius of the HSG cell
151 //-----------------------------------------------------------
152 //In order for consistent units, which in their raw form are:
153 //(1/(Gy^2))(keV/g)(cm^3/um^3)
154 // (1/Gy^2)(keV/um)(cm^3/g)(1/um^2)
155 //------------------------------------------------------------
156 //So convert (keV/g) into Gy it is scaled by (1.602E-16/0.001)(J/kg)
157 //And to get cm^3 into um^3 just scale by 10E12
158 //------------------------------------------------------------
159 double scaleGray = 1.602E-16/0.001;
160 double scaleLength = 1.E12;
161
162 //-------Calculate y*--------
163 double topIntegral = 0.;
164 double botIntegral = 0.;
165
166 //top integral
167 for (Int_t i = 1; i <= numberLinearBins; i++)
168 {
169 topIntegral += (1 - exp(-(histoLin->GetBinCenter(i) * histoLin->GetBinCenter(i)) / satParSqr))*histoLin->GetBinContent(i) * histoLin->GetBinWidth(i);
170 // (1 - exp(-y^2/y_{0}^2)) * f(y) * dy
171 }
172 //bottom integral
173 for (Int_t i = 1; i <= numberLinearBins; i++)
174 {
175 botIntegral += (histoLin->GetBinCenter(i) * histoLin->GetBinContent(i) * histoLin->GetBinWidth(i));
176 // y (size of each being time bin number) * f(y) * dy (bin size)
177 }
178
179 //cout << "Top: " << topIntegral << endl;
180 //cout << "Bot: " << botIntegral << endl;
181
182 //saturation corrected dose averaged linel energy value (y*)
183 double satCor = satParSqr * topIntegral / botIntegral;
184 //(y_{0})^2*top/bottom
185 //cout << "--------------------------------" << endl;
186 //cout << "y* = " << satCor << endl;
187 //cout << "--------------------------------" << endl;
188
189
190 //-----------alpha--------
191 double alpha = alpha0 + (beta / (rho*pi*rSV*rSV))*satCor*scaleGray*scaleLength;
192 //cout << "alpha = " << alpha << endl;
193
194 //-----------finally RBE itself--------
195 double ln01 = -2.302585093;
196 double RBE10;//= (2*beta*5.)((sqrt(alpha*alpha - 4*beta*ln01) - alpha));
197
198 double topRBE = (2 * beta*5.);
199 double botRBE = (sqrt(alpha*alpha - 4 * beta*ln01) - alpha);
200 RBE10 = topRBE / botRBE;
201
202 //cout << "--------------------------------" << endl;
203 cout << "RBE10 = " << RBE10 << endl;
204 //cout << "--------------------------------" << endl;
205
206 ////////////////////////////////////////////
207 //-----------Calculate yF-------------------
208 ////////////////////////////////////////////
209 //yF = int(y * f(y) dy)
210 double linyF = 0.;
211 for (int i = 1 ; i <= numberLinearBins; i++)
212 {
213 linyF += histoLin->GetBinCenter(i) * histoLin->GetBinContent(i) * histoLin->GetBinWidth(i) ;
214 }
215
216 cout << "yF = " << linyF << endl;
217
218 double logyF = 0;
219 for (int i = 1 ; i <= numberOfLogBins; i++)
220 {
221 logyF += histoLog->GetBinCenter(i) * histoLog->GetBinContent(i) * histoLog->GetBinWidth(i) ;
222 }
223 cout << "Log yF = " << logyF << endl;
224
225
226 ////////////////////////////////////////////
227 //--------------Calculate d(y)--------------
228 ////////////////////////////////////////////
229 //d(y) = y*f(y) / yF
230 for (int i = 1 ; i <= numberLinearBins; i++)
231 {
232 histoLin->SetBinContent(i, (histoLin->GetBinCenter(i) * histoLin->GetBinContent(i) / linyF));
233 }
234
235 //check that it's normalised to 1
236 double intOfHistDy = histoLin->Integral("width") ;
237 cout << "Int of d(y) = " << intOfHistDy << endl;
238
239 for (int bb = 1; bb <= numberOfLogBins-1; bb++)
240 {
241 histoLog->SetBinContent(bb, (histoLog->GetBinCenter(bb)* histoLog->GetBinContent(bb) / logyF));
242 }
243
244 double intOfHistLogDy = histoLog->Integral("width") ;
245 cout << "Int of log d(y) = " << intOfHistLogDy << endl;
246
247 ////////////////////////////////////////////
248 //-------Calculate Average Q(y)-----------
249 ////////////////////////////////////////////
250 //Q(y) = (a1/y)[1-exp(-(a2)y*y - (a3)y*y*y)]
251
252 //Quality factor values from ICRP 60 report
253 double qualityCoeff1 = 5510.;
254 double qualityCoeff2 = 5.E-5;
255 double qualityCoeff3 = 2.E-7;
256
257 double aveQuality = 0.;
258 //AveQuality = int(Q(y)d(y)dy)
259
260 for (int i = 1 ; i <= numberLinearBins; i++)
261 {
262 aveQuality += (qualityCoeff1/histoLin->GetBinCenter(i))* (1 - exp(-qualityCoeff2*pow(histoLin->GetBinCenter(i), 2.)-qualityCoeff3*pow(histoLin->GetBinCenter(i), 3.)))* histoLin->GetBinContent(i)* histoLin->GetBinWidth(i) ;
263 }
264 cout << "<Q> = " << aveQuality << endl;
265
266
267 ////////////////////////////////////////////
268 //---------Calculate yD--------------------
269 ////////////////////////////////////////////
270 //yD = int (y * d(y) dy)
271 double linyD = 0.;
272
273 for (int i = 1 ; i <= numberLinearBins; i++)
274 {
275 linyD += histoLin->GetBinCenter(i) * histoLin->GetBinContent(i) * histoLin -> GetBinWidth(i) ;
276 }
277
278 cout << "yD = " << linyD << endl;
279
280 double logyD = 0.;
281
282 for (int bb = 1; bb <= numberOfLogBins-1; bb++)
283 {
284 logyD += histoLog->GetBinCenter(bb) * histoLog->GetBinContent(bb) * histoLog -> GetBinWidth(bb) ;
285 }
286 cout << "log yD = " << logyD << endl;
287
288 ////////////////////////////////////////////
289 //---------Create yd(y)-------------------
290 ////////////////////////////////////////////
291 for (int bb = 1; bb <= numberOfLogBins-1; bb++)
292 {
293 histoLog->SetBinContent(bb, (histoLog->GetBinCenter(bb)* histoLog->GetBinContent(bb)) );
294 }
295 //Scale bins by log(10)/B
296 //See Appendix B of the ICRU 36 report for more details
297 histoLog->Scale((log(10) / B));
298
299 if (saveydySpec)
300 {
301 outFile.open("logYDYspec.txt");
302 for (int bb = 1; bb <= numberOfLogBins-1; bb++)
303 {
304 outFile << histoLog->GetBinCenter(bb) << "\t" << histoLog->GetBinContent(bb) << endl;
305 }
306 outFile.close();
307 }
308
309 if (saveImage)
310 {
311
312 TCanvas *c1 = new TCanvas("c1", "", 1000, 800);
313 //gStyle->SetOptStat(0); //un-comment to remove stat box on top right
314
315 histoLog->SetTitle("");
316 histoLog->GetXaxis()->SetTitle("y (keV/#mum)");
317 histoLog->GetYaxis()->SetTitle("yd(y)");
318
319 histoLog->GetYaxis()->SetLabelFont(42);
320 histoLog->GetXaxis()->SetLabelFont(42);
321 histoLog->GetYaxis()->SetTitleFont(42);
322 histoLog->GetXaxis()->SetTitleFont(42);
323 histoLog->GetYaxis()->CenterTitle(1);
324 histoLog->GetXaxis()->CenterTitle(1);
325 histoLog->GetXaxis()->SetTitleSize(0.045);
326 histoLog->GetYaxis()->SetTitleSize(0.045);
327 //histoLog->GetYaxis()->SetRangeUser(0., 1.8);
328 //histoLog -> GetXaxis()->SetRangeUser(0.01., 1000.);
329 histoLog->GetYaxis()->SetLabelSize(0.04);
330 histoLog->GetXaxis()->SetLabelSize(0.04);
331 histoLog->SetLineColor(1);
332 histoLog->SetLineWidth(4);
333 c1->SetLogx();
334 //c1->SetLogy();
335 //Plotting misbehaves in root6.XX
336 histoLog->Draw("");
337
338 c1->SaveAs(imageFileName.c_str());
339 }
340 }
341