| Geant4 Cross Reference |
1 //-------------------------------------------------------------------------------//
2 // This macrofile was developed by Konstantinos Chatzipapas at LP2iB (ex. CENBG) //
3 // in collaboration with the whole team of molecularDNA Geant4-DNA example //
4 // Publication: .................................... //
5 // For any question please contact through: //
6 // chatzipa@cenbg.in2p3.fr (or k.chatzipapas@yahoo.com) //
7 //-------------------------------------------------------------------------------//
8
9 // This macro requires the molecular-dna.root file generated from molecularDNA example
10 // To run this file just insert this command to the terminal:
11 // root .X analysis.C
12 // ROOT6.x should be installed
13
14 //***************************************//
15 // Please define the parameters below //
16 // ifile, r3, Nbp (as shown in terminal) //
17 //***************************************//
18
19 {
20 //*******************************************************************************//
21 // If you need to add multiple root outputs, by multithreading, use this command:
22 system ("hadd -O -f molecular-dna.root molecular-dna_t*.root");
23
24 // Define these parameters of the simulation
25 char ifile[256] = "molecular-dna.root"; // input filepath
26 Double_t r3 = 10575e-9 * 3450e-9 * 10575e-9; // a * b * c // Chromosome size, as defined in the mac file, but in meters. If sphere, a=b=c
27 Double_t Nbp = 6383.848410; // Mbp // Length of the DNA chain in Mbp
28 Double_t mass = 997 * 4 * 3.141592 * r3 / 3 ; // waterDensity * 4/3 * pi * r3 in kg
29 ///////////////////////////////////////////////////////////////////////////////////
30 //*******************************************************************************//
31
32 typedef std::pair <int64_t, int64_t> ipair;
33 bool greaterPair(const ipair &l, const ipair &r);
34 bool smallerPair(const ipair &l, const ipair &r);
35
36 void BinLogX(TH1 *h);
37
38 gROOT->Reset();
39 gStyle->SetPalette(1);
40 gROOT->SetStyle("Plain");
41 gStyle->SetOptStat(00000);
42
43 // Initialize output histograms
44 TCanvas *cfragment = new TCanvas("cfragment","DNA Fragments Distribution", 900, 120, 600,400);
45 cfragment->SetLogx();
46 cfragment->SetLogy();
47 TH1F *h1fragments = new TH1F("h1fragments","h1fragments",40,0,5);
48 BinLogX(h1fragments);
49
50 TCanvas *c1 = new TCanvas("c1", "Molecular DNA - Damage Quantification", 60, 120, 800, 800);
51 c1->SetBorderSize(0);
52 c1->SetFillColor(0);
53 c1->SetFillStyle(4000);
54 gPad->SetLeftMargin(0.13);
55
56 TPad* pad1 = new TPad("pad1","Species", 0, 0.51, 0.49, 1);
57 pad1->SetBorderSize(0);
58 pad1->SetFillColor(0);
59 pad1->SetFillStyle(4000);
60 pad1->SetLeftMargin(0.15);
61 pad1->SetRightMargin(0.01);
62 pad1->SetBottomMargin(0.2);
63
64 TPad* pad2 = new TPad("pad2","Damage Yield", 0.51, 0.5, 1, 1);
65 pad2->SetBorderSize(0);
66 pad2->SetFillColor(0);
67 pad2->SetFillStyle(4000);
68 pad2->SetLeftMargin(0.15);
69 pad2->SetRightMargin(0.05);
70 pad2->SetBottomMargin(0.2);
71
72 TPad* pad3 = new TPad("pad3","Breaks Yield SSB", 0, 0, 0.49, 0.49);
73 pad3->SetBorderSize(0);
74 pad3->SetFillColor(0);
75 pad3->SetFillStyle(4000);
76 pad3->SetLeftMargin(0.15);
77 pad3->SetRightMargin(0.01);
78 //pad3->SetTopMargin(0.2);
79 pad3->SetBottomMargin(0.2);
80
81 TPad* pad4 = new TPad("pad4","Breaks Yield DSB", 0.51, 0, 1, 0.49);
82 pad4->SetBorderSize(0);
83 pad4->SetFillColor(0);
84 pad4->SetFillStyle(4000);
85 pad4->SetLeftMargin(0.15);
86 pad4->SetRightMargin(0.05);
87 //pad3->SetTopMargin(0.2);
88 pad4->SetBottomMargin(0.2);
89
90 pad1->Draw();
91 pad2->Draw();
92 pad3->Draw();
93 pad4->Draw();
94
95 // Open root file
96 TFile *f = TFile::Open(ifile);
97
98 // Initialize Variables
99 Int_t EB, ES, OHB, OHS, HB, HS, FL;
100 Int_t total_EB, total_ES, total_OHB, total_OHS, total_HB, total_HS, total_FL;
101 Float_t total_EB2, total_ES2, total_OHB2, total_OHS2, total_HB2, total_HS2, total_FL2;
102 Float_t SD_EB, SD_ES, SD_OHB, SD_OHS, SD_HB, SD_HS;
103 Float_t SD_SSB, SD_SSBp, SD_SSB2p, SD_sSSB, SD_SSBd, SD_SSBi;
104 Float_t SD_DSB, SD_DSBp, SD_DSBpp, SD_sDSB, SD_DSBd, SD_DSBi, SD_DSBm, SD_DSBh;
105
106 Int_t SSB, SSBp, SSB2p;
107 Int_t total_SSB, total_SSBp, total_SSB2p;
108 Float_t total_SSB2, total_SSBp2, total_SSB2p2;
109 Int_t DSB, DSBp, DSBpp;
110 Int_t total_DSB, total_DSBp, total_DSBpp;
111 Float_t total_DSB2, total_DSBp2, total_DSBpp2;
112
113 Int_t SSBd, SSBi, SSBm;
114 Int_t total_sSSB, total_SSBd, total_SSBi, total_SSBm;
115 Float_t total_sSSB2, total_SSBd2, total_SSBi2, total_SSBm2;
116 Int_t DSBd, DSBi, DSBm, DSBh;
117 Int_t total_sDSB, total_DSBd, total_DSBi, total_DSBm, total_DSBh;
118 Float_t total_sDSB2, total_DSBd2, total_DSBi2, total_DSBm2, total_DSBh2;
119
120 Double_t dose = 0;
121 Double_t SD_dose = 0;
122
123 Double_t EB_yield = 0; Double_t ES_yield = 0; Double_t OHB_yield = 0; Double_t OHS_yield = 0; Double_t HB_yield = 0; Double_t HS_yield = 0;
124 Double_t SD_EB_yield = 0; Double_t SD_ES_yield = 0; Double_t SD_OHB_yield = 0; Double_t SD_OHS_yield = 0; Double_t SD_HB_yield = 0; Double_t SD_HS_yield = 0;
125
126 Double_t SSB_yield = 0; Double_t SSBp_yield = 0; Double_t SSB2p_yield = 0;
127 Double_t SD_SSB_yield = 0; Double_t SD_SSBp_yield = 0; Double_t SD_SSB2p_yield = 0;
128 Double_t DSB_yield = 0; Double_t DSBp_yield = 0; Double_t DSBpp_yield = 0;
129 Double_t SD_DSB_yield = 0; Double_t SD_DSBp_yield = 0; Double_t SD_DSBpp_yield = 0;
130
131 Double_t sSSB_yield = 0; Double_t SSBi_yield = 0; Double_t SSBd_yield = 0; Double_t SSBm_yield = 0;
132 Double_t SD_sSSB_yield = 0; Double_t SD_SSBi_yield = 0; Double_t SD_SSBd_yield = 0; Double_t SD_SSBm_yield = 0;
133 Double_t sDSB_yield = 0; Double_t DSBi_yield = 0; Double_t DSBd_yield = 0; Double_t DSBm_yield = 0; Double_t DSBh_yield = 0;
134 Double_t SD_sDSB_yield = 0; Double_t SD_DSBi_yield = 0; Double_t SD_DSBd_yield = 0; Double_t SD_DSBm_yield = 0; Double_t SD_DSBh_yield = 0;
135
136 total_EB = 0; total_ES = 0; total_OHB = 0; total_OHS = 0; total_HB = 0; total_HS = 0;
137
138 total_SSB = 0; total_SSBp = 0; total_SSB2p = 0;
139 total_SSB2 = 0; total_SSBp2 = 0; total_SSB2p2 = 0;
140 total_DSB = 0; total_DSBp = 0; total_DSBpp = 0;
141 total_DSB2 = 0; total_DSBp2 = 0; total_DSBpp2 = 0;
142
143 total_sSSB = 0; total_SSBd = 0; total_SSBi = 0; total_SSBm = 0;
144 total_sSSB2 = 0; total_SSBd2 = 0; total_SSBi2 = 0; total_SSBm2 = 0;
145 total_sDSB = 0; total_DSBd = 0; total_DSBi = 0; total_DSBm = 0; total_DSBh = 0;
146 total_sDSB2 = 0; total_DSBd2 = 0; total_DSBi2 = 0; total_DSBm2 = 0; total_DSBh2 = 0;
147
148 Double_t eVtoJ = 1.60218e-19;
149 Double_t EnergyDeposited_eV = 0;
150 Double_t acc_edep = 0;
151 Double_t acc_edep2 = 0;
152
153 Double_t Energy;
154 Double_t BPID;
155 Char_t Primary;
156 char *primaryName = new char[32];
157 char *type= new char[256];
158
159 // Read trees and leaves from root file, and give values to variables
160 TTree* tree = (TTree*) f->Get("tuples/primary_source");
161 Float_t number = (Float_t) tree->GetEntries();
162
163 vector<pair<int,int64_t>> DSBBPID;
164
165 // For reading species production
166 tree = (TTree*) f->Get("tuples/damage");
167 tree->SetBranchAddress("Primary", &Primary);
168 tree->SetBranchAddress("Energy", &Energy);
169 tree->SetBranchAddress("EaqBaseHits", &EB);
170 tree->SetBranchAddress("EaqStrandHits", &ES);
171 tree->SetBranchAddress("OHBaseHits", &OHB);
172 tree->SetBranchAddress("OHStrandHits", &OHS);
173 tree->SetBranchAddress("HBaseHits", &HB);
174 tree->SetBranchAddress("HStrandHits", &HS);
175 tree->SetBranchAddress("TypeClassification", type);
176 tree->SetBranchAddress("BasePair", &BPID);
177
178
179 Long64_t nentries = tree->GetEntries();
180 for(int i = 0;i<nentries;i++){
181 tree->GetEntry(i);
182
183 total_EB += EB;
184 total_EB2 += pow(EB,2);
185 total_ES += ES;
186 total_ES2 += pow(ES,2);
187 total_OHB += OHB;
188 total_OHB2 += pow(OHB,2);
189 total_OHS += OHS;
190 total_OHS2 += pow(OHS,2);
191 total_HB += HB;
192 total_HB2 += pow(HB,2);
193 total_HS += HS;
194 total_HS2 += pow(HS,2);
195
196 if((string)type=="DSB"||(string)type=="DSB+"||(string)type=="DSB++"){
197 //cout << "DSB:"<<type<<endl;
198 DSBBPID.push_back(make_pair(i,(int64_t)BPID));
199 }
200
201 }
202
203 // Sort DSBs from the one with lower ID value to the one with higher ID value
204 // Then find the number of fragments that have been produced
205 sort(DSBBPID.begin(), DSBBPID.end(), smallerPair);
206 for(int ie = 0;ie<DSBBPID.size()-1;ie++){
207 int64_t dsbfragment = DSBBPID[ie+1].second-DSBBPID[ie].second;
208
209 double val = (double)dsbfragment/1000.;
210 double meanw = h1fragments->GetBinCenter(h1fragments->FindBin(val));
211 double binw = h1fragments->GetBinWidth (h1fragments->FindBin(val));
212 h1fragments->Fill(val,1./binw/1000);//bp-1
213 //cout <<"val:"<<val<<endl;
214 }
215
216 // Calculate the standard deviation of species
217 SD_EB = sqrt(((total_EB2 / number) - pow(total_EB / number,2))/(number -1));
218 SD_ES = sqrt(((total_ES2 / number) - pow(total_ES / number,2))/(number -1));
219 SD_OHB = sqrt(((total_OHB2 / number) - pow(total_OHB / number,2))/(number -1));
220 SD_OHS = sqrt(((total_OHS2 / number) - pow(total_OHS / number,2))/(number -1));
221 SD_HB = sqrt(((total_HB2 / number) - pow(total_HB / number,2))/(number -1));
222 SD_HS = sqrt(((total_HS2 / number) - pow(total_HS / number,2))/(number -1));
223
224 // Read damage classification SSB, SSB+, 2SSB, DSB, DSB+, DSB++
225 // As they have been defined in: Nikjoo, H., O’Neill, O., Goodhead, T., & Terrissol, M. 1997,
226 // Computational modelling of low-energy electron-induced DNA damage by early physical
227 // and chemical events, International Journal of Radiation Biology, 71, 467.
228 tree = (TTree *) f->Get("tuples/classification");
229 tree->SetBranchAddress("Primary",&Primary);
230 tree->SetBranchAddress("Energy", &Energy);
231 tree->SetBranchAddress("SSB", &SSB);
232 tree->SetBranchAddress("SSBp", &SSBp);
233 tree->SetBranchAddress("2SSB", &SSB2p);
234 tree->SetBranchAddress("DSB", &DSB);
235 tree->SetBranchAddress("DSBp", &DSBp);
236 tree->SetBranchAddress("DSBpp", &DSBpp);
237
238
239 Long64_t nentriesC = tree->GetEntries();
240 for(int i = 0;i<nentriesC;i++){
241 tree->GetEntry(i);
242
243 total_SSBp += SSBp;
244 total_SSBp2 += pow(SSBp,2);
245 total_SSB2p += SSB2p;
246 total_SSB2p2 += pow(SSB2p,2);
247 total_SSB += SSB;
248 total_SSB2 += pow(SSB,2);
249
250 total_DSBp += DSBp;
251 total_DSBp2 += pow(DSBp,2);
252 total_DSBpp += DSBpp;
253 total_DSBpp2 += pow(DSBpp,2);
254 total_DSB += DSB;
255 total_DSB2 += pow(DSB,2);
256
257 }
258
259 // Calculate the standard deviation
260 SD_SSB = sqrt(((total_SSB2 / number) - pow(total_SSB / number,2))/(number -1));
261 SD_SSBp = sqrt(((total_SSBp2 / number) - pow(total_SSBp / number,2))/(number -1));
262 SD_SSB2p = sqrt(((total_SSB2p2 / number) - pow(total_SSB2p / number,2))/(number -1));
263
264 SD_DSB = sqrt(((total_DSB2 / number) - pow(total_DSB / number,2))/(number -1));
265 SD_DSBp = sqrt(((total_DSBp2 / number) - pow(total_DSBp / number,2))/(number -1));
266 SD_DSBpp = sqrt(((total_DSBpp2 / number) - pow(total_DSBpp / number,2))/(number -1));
267
268 // Read damage classification SSBd, SSBi, SSBm, DSBd, DSBi, DSBm, DSBh
269 // As they have been defined in: Nikjoo, H., O’Neill, O., Goodhead, T., & Terrissol, M. 1997,
270 // Computational modelling of low-energy electron-induced DNA damage by early physical
271 // and chemical events, International Journal of Radiation Biology, 71, 467.
272 tree = (TTree *) f->Get("tuples/source");
273 tree->SetBranchAddress("Primary",primaryName);
274 tree->SetBranchAddress("Energy", &Energy);
275 tree->SetBranchAddress("SSBd", &SSBd);
276 tree->SetBranchAddress("SSBi", &SSBi);
277 tree->SetBranchAddress("SSBm", &SSBm);
278 tree->SetBranchAddress("DSBd", &DSBd);
279 tree->SetBranchAddress("DSBi", &DSBi);
280 tree->SetBranchAddress("DSBm", &DSBm);
281 tree->SetBranchAddress("DSBh", &DSBh);
282
283 Long64_t nentriesS = tree->GetEntries();
284 for(int i = 0;i<nentriesS;i++){
285 tree->GetEntry(i);
286
287 total_SSBd += SSBd;
288 total_SSBd2 += pow((SSBd),2);
289 total_SSBi += SSBi;
290 total_SSBi2 += pow((SSBi),2);
291 total_SSBm += SSBm;
292 total_SSBm2 += pow((SSBm),2);
293 total_sSSB += SSBd + SSBi + SSBm;
294 total_sSSB2 += pow((SSBd+SSBi+SSBm),2);
295
296 total_DSBd += DSBd;
297 total_DSBd2 += pow(DSBd,2);
298 total_DSBi += DSBi;
299 total_DSBi2 += pow(DSBi,2);
300 total_DSBm += DSBm;
301 total_DSBm2 += pow(DSBm,2);
302 total_DSBh += DSBh;
303 total_DSBh2 += pow(DSBh,2);
304 total_sDSB += DSBd + DSBi + DSBm + DSBh;
305 total_sDSB2 += pow((DSBd+DSBi+DSBm+DSBh),2);
306
307 }
308
309 // Calculate the standard deviation
310 SD_sSSB = sqrt(((total_sSSB2 / number) - pow(total_sSSB / number,2))/(number -1));
311 SD_SSBd = sqrt(((total_SSBd2 / number) - pow(total_SSBd / number,2))/(number -1));
312 SD_SSBi = sqrt(((total_SSBi2 / number) - pow(total_SSBi / number,2))/(number -1));
313 SD_SSBm = sqrt(((total_SSBm2 / number) - pow(total_SSBm / number,2))/(number -1));
314
315 SD_sDSB = sqrt(((total_sDSB2 / number) - pow(total_sDSB / number,2))/(number -1));
316 SD_DSBd = sqrt(((total_DSBd2 / number) - pow(total_DSBd / number,2))/(number -1));
317 SD_DSBi = sqrt(((total_DSBi2 / number) - pow(total_DSBi / number,2))/(number -1));
318 SD_DSBm = sqrt(((total_DSBm2 / number) - pow(total_DSBm / number,2))/(number -1));
319 SD_DSBh = sqrt(((total_DSBh2 / number) - pow(total_DSBh / number,2))/(number -1));
320
321
322 // Measure the Deposited Energy in the whole volume that includes DNA chain (chromosome)
323
324 tree = (TTree *) f->Get("tuples/chromosome_hits");
325 tree->SetBranchAddress("e_chromosome_kev",&EnergyDeposited_eV);
326 nentries = tree->GetEntries();
327 for(int i = 0;i<nentries;i++){
328 tree->GetEntry(i);
329 acc_edep += EnergyDeposited_eV *1e3;
330 acc_edep2 += EnergyDeposited_eV *EnergyDeposited_eV *1e6;
331 }
332 tree->SetBranchAddress("e_dna_kev",&EnergyDeposited_eV);
333 nentries = tree->GetEntries();
334 for(int i = 0;i<nentries;i++){
335 tree->GetEntry(i);
336 acc_edep += EnergyDeposited_eV *1e3;
337 acc_edep2 += EnergyDeposited_eV *EnergyDeposited_eV *1e6;
338 }
339
340 // Close the root file to free space
341 f->Close();
342
343 // Calculate the absorbed dose
344 dose = acc_edep * eVtoJ / mass;
345
346 // This is a normalization factor to produce the output in Gy-1 Gbp-1, or else.
347 // Default value is 1 to produce the result in Gy-1 Mbp-1
348 // It changes Mbp to Gbp. Some other changes may be needed in graphs section (name of axes)
349 double norm = 1000;
350
351 // Calculate the yields, together with their standard deviation
352 EB_yield = (Double_t) total_EB / dose / Nbp;
353 ES_yield = (Double_t) total_ES / dose / Nbp;
354 OHB_yield = (Double_t) total_OHB / dose / Nbp;
355 OHS_yield = (Double_t) total_OHS / dose / Nbp;
356 HB_yield = (Double_t) total_HB / dose / Nbp;
357 HS_yield = (Double_t) total_HS / dose / Nbp;
358
359 SD_EB_yield = SD_EB / dose / Nbp;
360 SD_ES_yield = SD_ES / dose / Nbp;
361 SD_OHB_yield = SD_OHB / dose / Nbp;
362 SD_OHS_yield = SD_OHS / dose / Nbp;
363 SD_HB_yield = SD_HB / dose / Nbp;
364 SD_HS_yield = SD_HS / dose / Nbp;
365
366
367 SSB_yield = (Double_t) norm * total_SSB / dose / Nbp;
368 SSBp_yield = (Double_t) norm * total_SSBp / dose / Nbp;
369 SSB2p_yield = (Double_t) norm * total_SSB2p / dose / Nbp;
370
371 DSB_yield = (Double_t) norm * total_DSB / dose / Nbp;
372 DSBp_yield = (Double_t) norm * total_DSBp / dose / Nbp;
373 DSBpp_yield = (Double_t) norm * total_DSBpp / dose / Nbp;
374
375 SD_SSB_yield = norm * SD_SSB / dose / Nbp;
376 SD_SSBp_yield = norm * SD_SSBp / dose / Nbp;
377 SD_SSB2p_yield = norm * SD_SSB2p / dose / Nbp;
378
379 SD_DSB_yield = norm * SD_DSB / dose / Nbp;
380 SD_DSBp_yield = norm * SD_DSBp / dose / Nbp;
381 SD_DSBpp_yield = norm * SD_DSBpp / dose / Nbp;
382
383
384 sSSB_yield = (Double_t) norm * total_sSSB / dose / Nbp;
385 SSBi_yield = (Double_t) norm * total_SSBi / dose / Nbp;
386 SSBd_yield = (Double_t) norm * total_SSBd / dose / Nbp;
387 SSBm_yield = (Double_t) norm * total_SSBm / dose / Nbp;
388
389 sDSB_yield = (Double_t) norm * total_sDSB / dose / Nbp;
390 DSBi_yield = (Double_t) norm * total_DSBi / dose / Nbp;
391 DSBd_yield = (Double_t) norm * total_DSBd / dose / Nbp;
392 DSBm_yield = (Double_t) norm * total_DSBm / dose / Nbp;
393 DSBh_yield = (Double_t) norm * total_DSBh / dose / Nbp;
394
395 SD_sSSB_yield = norm * SD_sSSB / dose / Nbp;
396 SD_SSBi_yield = norm * SD_SSBi / dose / Nbp;
397 SD_SSBd_yield = norm * SD_SSBd / dose / Nbp;
398 SD_SSBm_yield = norm * SD_SSBm / dose / Nbp;
399
400 SD_sDSB_yield = norm * SD_sDSB / dose / Nbp;
401 SD_DSBi_yield = norm * SD_DSBi / dose / Nbp;
402 SD_DSBd_yield = norm * SD_DSBd / dose / Nbp;
403 SD_DSBm_yield = norm * SD_DSBm / dose / Nbp;
404 SD_DSBh_yield = norm * SD_DSBh / dose / Nbp;
405
406
407 // Print output in terminal
408
409 float total_SSB_totalYield = SSB_yield + SSBp_yield + SSB2p_yield;
410 float total_DSB_totalYield = DSB_yield + DSBp_yield + DSBpp_yield;
411
412 cout<<"\n" <<ifile <<'\n'
413 <<"\nDose Absorbed (Gy): " <<dose <<'\n'
414 <<"Particle : " <<primaryName <<'\t'
415 <<"Energy (MeV) : " <<Energy <<'\t'
416 <<"Number of Primaries : " <<number <<'\n'
417 <<" Output Damage : " <<'\n'<<'\t'
418 <<" Species Hits (Gy-1 Mbp-1) " <<'\n'<<'\t'
419 <<"EaqBaseHits : " <<EB_yield <<" \t" <<" error %: " <<100*SD_EB_yield/EB_yield <<'\n'<<'\t'
420 <<"EaqStrandHits : " <<ES_yield <<" \t" <<" error %: " <<100*SD_ES_yield/ES_yield <<'\n'<<'\t'
421 <<"OHBaseHits : " <<OHB_yield <<" \t" <<" error %: " <<100*SD_OHB_yield/OHB_yield <<'\n'<<'\t'
422 <<"OHStrandHits : " <<OHS_yield <<" \t" <<" error %: " <<100*SD_OHS_yield/OHS_yield <<'\n'<<'\t'
423 <<"HBaseHits : " <<HB_yield <<" \t" <<" error %: " <<100*SD_HB_yield/HB_yield <<'\n'<<'\t'
424 <<"HStrandHits : " <<HS_yield <<" \t" <<" error %: " <<100*SD_HS_yield/HS_yield <<'\n'<<'\n'<<'\t'
425 <<" Damage yield (Gy-1 Gbp-1) " <<'\n'<<'\t'
426 <<"SSB : " <<SSB_yield <<" \t" <<" error %: " <<100*SD_SSB_yield/SSB_yield <<'\n'<<'\t'
427 <<"SSB+ : " <<SSBp_yield <<" \t" <<" error %: " <<100*SD_SSBp_yield/SSBp_yield <<'\n'<<'\t'
428 <<"2SSB : " <<SSB2p_yield <<" \t" <<" error %: " <<100*SD_SSB2p_yield/SSB2p_yield <<'\n'<<'\t'
429 <<"SSB total : " <<total_SSB_totalYield <<'\n'<<'\t'
430 <<"DSB : " <<DSB_yield <<" \t" <<" error %: " <<100*SD_DSB_yield/DSB_yield <<'\n'<<'\t'
431 <<"DSB+ : " <<DSBp_yield <<" \t" <<" error %: " <<100*SD_DSBp_yield/DSBp_yield <<'\n'<<'\t'
432 <<"DSB++ : " <<DSBpp_yield <<" \t" <<" error %: " <<100*SD_DSBpp_yield/DSBpp_yield <<'\n'<<'\t'
433 <<"DSB total : " <<total_DSB_totalYield <<'\n'<<'\n'<<'\t'
434 <<" Breaks yield (Gy-1 Gbp-1) " <<'\n'<<'\t'
435 <<"SSB direct : " <<SSBd_yield <<" \t" <<" error %: " <<100*SD_SSBd_yield/SSBd_yield <<'\n'<<'\t'
436 <<"SSB indirect : " <<SSBi_yield <<" \t" <<" error %: " <<100*SD_SSBi_yield/SSBi_yield <<'\n'<<'\t'
437 <<"SSB mixed : " <<SSBm_yield <<" \t" <<" error %: " <<100*SD_SSBm_yield/SSBi_yield <<'\n'<<'\t'
438 <<"SSB total : " <<sSSB_yield <<" \t" <<" error %: " <<100*SD_sSSB_yield/sSSB_yield <<'\n'<<'\t'
439 <<"DSB direct : " <<DSBd_yield <<" \t" <<" error %: " <<100*SD_DSBd_yield/DSBd_yield <<'\n'<<'\t'
440 <<"DSB indirect : " <<DSBi_yield <<" \t" <<" error %: " <<100*SD_DSBi_yield/DSBi_yield <<'\n'<<'\t'
441 <<"DSB mixed : " <<DSBm_yield <<" \t" <<" error %: " <<100*SD_DSBm_yield/DSBm_yield <<'\n'<<'\t'
442 <<"DSB hybrid : " <<DSBh_yield <<" \t" <<" error %: " <<100*SD_DSBh_yield/DSBh_yield <<'\n'<<'\t'
443 <<"DSB total : " <<sDSB_yield <<" \t" <<" error %: " <<100*SD_sDSB_yield/sDSB_yield <<'\n'<<'\n'<<'\t'
444 <<"SSB/DSB : " <<sSSB_yield/sDSB_yield <<'\n'<<'\n';
445
446
447 // Plot Histograms
448
449 cfragment->GetCanvas()->cd();
450 h1fragments->SetStats(false);
451 h1fragments->SetMarkerSize(0.1);
452 h1fragments->SetMarkerColor(kRed);
453 h1fragments->SetLineColor (kRed);
454 h1fragments->Scale(1./(Nbp*1e6)); //bp^-1
455 h1fragments->SetTitle("");
456 h1fragments->SetYTitle("Number of Fragments (bp^{-2})");
457 h1fragments->SetXTitle("Fragment Length (kbp)");
458 h1fragments->SetAxisRange(10,1e4);
459 h1fragments->SetMaximum(3e-11);
460 h1fragments->SetMinimum(1e-15);
461 h1fragments->Draw();
462
463
464 c1->GetCanvas()->cd();
465 pad1->cd();
466 const Int_t n = 6;
467 Double_t x[n] = {1,2,3,4,5,6};
468 Double_t y[n] = {EB_yield,ES_yield,OHB_yield,OHS_yield,HB_yield,HS_yield};
469 Double_t err_y[n] = {SD_EB_yield,SD_ES_yield,SD_OHB_yield,SD_OHS_yield,SD_HB_yield,SD_HS_yield};
470 TGraph* gr = new TGraphErrors(n,x,y,0,err_y);
471 gr->SetTitle("Species");
472 gr->GetXaxis()->SetBinLabel(9, "EaqBaseHits");
473 gr->GetXaxis()->SetBinLabel(25,"EaqStrandHits");
474 gr->GetXaxis()->SetBinLabel(42,"OHBaseHits");
475 gr->GetXaxis()->SetBinLabel(58,"OHStrandHits");
476 gr->GetXaxis()->SetBinLabel(75,"HBaseHits");
477 gr->GetXaxis()->SetBinLabel(92,"HStrandHits");
478 gr->GetYaxis()->SetTitle("Species Hits (Gy^{-1} Mbp^{-1})");
479 gr->GetYaxis()->SetTitleOffset(2);
480
481 gr->SetFillColor(49);
482 gr->Draw("ba");
483
484
485 pad2->cd();
486 Double_t x2[n] = {1,2,3,4,5,6};
487 Double_t y2[n] = {SSBp_yield,SSB2p_yield,SSB_yield,DSBp_yield,DSBpp_yield,DSB_yield};
488 Double_t err_y2[n] = {SD_SSBp_yield,SD_SSB2p_yield,SD_SSB_yield,SD_DSBp_yield,SD_DSBpp_yield,SD_DSB_yield};
489 TGraph* gr2 = new TGraphErrors(n,x2,y2,0,err_y2);
490 gr2->SetTitle("Damage Yield");
491 gr2->GetXaxis()->SetBinLabel(9, "SSB+");
492 gr2->GetXaxis()->SetBinLabel(25,"2SSB");
493 gr2->GetXaxis()->SetBinLabel(42,"SSB");
494 gr2->GetXaxis()->SetBinLabel(58,"DSB+");
495 gr2->GetXaxis()->SetBinLabel(75,"DSB++");
496 gr2->GetXaxis()->SetBinLabel(92,"DSB");
497 //gr2->GetYaxis()->SetTitle("Damage yield (Gy^{-1} Mbp^{-1})");
498 gr2->GetYaxis()->SetTitle("Damage yield (Gy^{-1} Gbp^{-1})");
499 gr2->GetYaxis()->SetTitleOffset(2);
500
501 gr2->SetFillColor(8);
502 gr2->Draw("ba");
503
504
505 pad3->cd();
506 const Int_t m = 4;
507 Double_t x3[m] = {1,2,3,4};
508 Double_t y3[m] = {SSBd_yield,SSBi_yield,SSBm_yield,sSSB_yield};
509 Double_t err_y3[m] = {SD_SSBd_yield,SD_SSBi_yield,SD_SSBm_yield,SD_sSSB_yield};
510 TGraph* gr3 = new TGraphErrors(m,x3,y3,0,err_y3);
511 gr3->SetTitle("Breaks Yield");
512 gr3->GetXaxis()->SetBinLabel(8, "SSB direct");
513 gr3->GetXaxis()->SetBinLabel(35,"SSB indirect");
514 gr3->GetXaxis()->SetBinLabel(64,"SSB mixed");
515 gr3->GetXaxis()->SetBinLabel(92,"SSB all");
516 //gr3->GetYaxis()->SetTitle("Breaks yield (Gy^{-1} Mbp^{-1})");
517 gr3->GetYaxis()->SetTitle("SSB yield (Gy^{-1} Gbp^{-1})");
518 gr3->GetYaxis()->SetTitleOffset(2);
519
520 gr3->SetFillColor(7);
521 gr3->Draw("ba");
522
523
524 pad4->cd();
525 const Int_t k = 5;
526 Double_t x4[k] = {1,2,3,4,5};
527 Double_t y4[k] = {DSBd_yield,DSBi_yield,DSBm_yield,DSBh_yield,sDSB_yield};
528 Double_t err_y4[k] = {SD_DSBd_yield,SD_DSBi_yield,SD_DSBm_yield,SD_DSBh_yield,SD_sDSB_yield};
529 TGraph* gr4 = new TGraphErrors(k,x4,y4,0,err_y4);
530 gr4->SetTitle("Breaks Yield");
531 gr4->GetXaxis()->SetBinLabel(8,"DSB direct");
532 gr4->GetXaxis()->SetBinLabel(29,"DSB indirect");
533 gr4->GetXaxis()->SetBinLabel(50,"DSB mixed");
534 gr4->GetXaxis()->SetBinLabel(71,"DSB hybrid");
535 gr4->GetXaxis()->SetBinLabel(92,"DSB all");
536 //gr4->GetYaxis()->SetTitle("Breaks yield (Gy^{-1} Mbp^{-1})");
537 gr4->GetYaxis()->SetTitle("DSB yield (Gy^{-1} Gbp^{-1})");
538 gr4->GetYaxis()->SetTitleOffset(2);
539
540 gr4->SetFillColor(4);
541 gr4->Draw("ba");
542
543 }
544
545 // Some important bools that are needed to run the root macro file
546 bool greaterPair(const ipair& l, const ipair& r){return l.second > r.second;}
547 bool smallerPair(const ipair& l, const ipair& r){return l.second < r.second;}
548
549 void BinLogX(TH1 *h) {
550 TAxis *axis = h->GetXaxis();
551 int bins = axis->GetNbins();
552 Axis_t from = axis->GetXmin();
553 Axis_t to = axis->GetXmax();
554 Axis_t width = (to - from) / bins;
555 Axis_t *new_bins = new Axis_t[bins + 1];
556 for (int i = 0; i <= bins; i++) {
557 new_bins[i] = TMath::Power(10, from + i * width);
558 }
559 axis->Set(bins, new_bins);
560 delete[] new_bins;
561
562 }
563