| Geant4 Cross Reference |
1 //*********************************************************************************
2 // Modified by Sara Zein to calculate the damage probability per plasmid
3 //_________________________________________________________________________________
4 //*********************************************************************************
5
6 // This macro requires the molecular-dna.root file generated from molecularDNA example
7 // To run this file just insert this command to the terminal:
8 // root .X plasmid.C
9
10 //***************************************//
11 // Please define the parameters below //
12 // ifile, r3, Nbp (as shown in terminal) //
13 //***************************************//
14
15 {
16 //*******************************************************************************//
17 // If you need to add multiple root outputs, by multithreading, use this command:
18 system("hadd -O -f molecular-dna.root molecular-dna_t*.root");
19
20 // Define these parameters of the simulation
21 char ifile[256] = "molecular-dna.root"; // input filepath
22 const Int_t numberOfPlasmids = 10144;
23 Double_t r3 = 4.42e-6 * 4.42e-6 * 4.84e-6; // r * r * r (world cube side m)
24 Double_t Nbp = 4.367 * 0.001 * numberOfPlasmids; // Mbp // Length of the DNA chain in Mbp
25 // multiplying by 10142 since we are testing 10142 k plasmids
26 Double_t mass = 997 * r3; // density * r3
27
28 //*******************************************************************************//
29
30 typedef std::pair<int64_t, int64_t> ipair;
31 bool greaterPair(const ipair& l, const ipair& r);
32 bool smallerPair(const ipair& l, const ipair& r);
33
34 void BinLogX(TH1 * h);
35
36 gROOT->Reset();
37 gStyle->SetPalette(1);
38 gROOT->SetStyle("Plain");
39 gStyle->SetOptStat(00000);
40
41 gStyle->SetPadTickX(1);
42 gStyle->SetPadTickY(1);
43
44 // Initialize output histograms
45 TCanvas* cdamage = new TCanvas("cdamage", "DNA Damage Distribution", 900, 120, 600, 400);
46 cdamage->SetLogy();
47 TH1F* h1damage = new TH1F("h1damage", "h1damage", 28, 0, 14);
48 TCanvas* c4damage = new TCanvas("c4damage", "Direct damage per plasmid", 900, 120, 600, 400);
49 // c4damage->SetLogy();
50 TH1F* h4damage = new TH1F("h4damage", "h4damage", 28, 0, 14);
51 TCanvas* cSSB = new TCanvas("cSSB", "SSB Distribution", 900, 120, 600, 400);
52 cSSB->SetLogy();
53 TH1F* h1SSB = new TH1F("h1SSB", "h1SSB", 28, 0, 14);
54 TCanvas* cDSB = new TCanvas("cDSB", "DSB Distribution", 900, 120, 600, 400);
55 cDSB->SetLogy();
56 TH1F* h1DSB = new TH1F("h1DSB", "h1DSB", 40, 0, 10);
57 TCanvas* ccount = new TCanvas("cCount", "Damage per plasmid Distribution", 900, 120, 600, 400);
58 TH1F* h1count = new TH1F("h1count", "h1count", 11000, 0, 11000);
59 TGraph2D* gr1 = new TGraph2D();
60
61 // Open root file
62 TFile* f = TFile::Open(ifile);
63
64 // Initialize Variables
65 Int_t EB, ES, OHB, OHS, HB, HS, FL;
66 Int_t total_EB, total_ES, total_OHB, total_OHS, total_HB, total_HS, total_FL;
67 Float_t total_EB2, total_ES2, total_OHB2, total_OHS2, total_HB2, total_HS2, total_FL2;
68 Float_t SD_EB, SD_ES, SD_OHB, SD_OHS, SD_HB, SD_HS;
69 Float_t SD_SSB, SD_SSBp, SD_SSB2p, SD_sSSB, SD_SSBd, SD_SSBi, SD_SSBm;
70 Float_t SD_DSB, SD_DSBp, SD_DSBpp, SD_sDSB, SD_DSBd, SD_DSBi, SD_DSBm, SD_DSBh;
71
72 Int_t SSB, SSBp, SSB2p;
73 Int_t total_SSB, total_SSBp, total_SSB2p;
74 Float_t total_SSB2, total_SSBp2, total_SSB2p2;
75 Int_t DSB, DSBp, DSBpp;
76 Int_t total_DSB, total_DSBp, total_DSBpp;
77 Float_t total_DSB2, total_DSBp2, total_DSBpp2;
78
79 Int_t SSBd, SSBi, SSBm;
80 Int_t total_sSSB, total_SSBd, total_SSBi, total_SSBm;
81 Float_t total_sSSB2, total_SSBd2, total_SSBi2, total_SSBm2;
82 Int_t DSBd, DSBi, DSBm, DSBh;
83 Int_t total_sDSB, total_DSBd, total_DSBi, total_DSBm, total_DSBh;
84 Float_t total_sDSB2, total_DSBd2, total_DSBi2, total_DSBm2, total_DSBh2;
85
86 Double_t dose = 0;
87 Double_t SD_dose = 0;
88
89 Double_t EB_yield = 0;
90 Double_t ES_yield = 0;
91 Double_t OHB_yield = 0;
92 Double_t OHS_yield = 0;
93 Double_t HB_yield = 0;
94 Double_t HS_yield = 0;
95 Double_t SD_EB_yield = 0;
96 Double_t SD_ES_yield = 0;
97 Double_t SD_OHB_yield = 0;
98 Double_t SD_OHS_yield = 0;
99 Double_t SD_HB_yield = 0;
100 Double_t SD_HS_yield = 0;
101
102 Double_t SSB_yield = 0;
103 Double_t SSBp_yield = 0;
104 Double_t SSB2p_yield = 0;
105 Double_t SD_SSB_yield = 0;
106 Double_t SD_SSBp_yield = 0;
107 Double_t SD_SSB2p_yield = 0;
108 Double_t DSB_yield = 0;
109 Double_t DSBp_yield = 0;
110 Double_t DSBpp_yield = 0;
111 Double_t SD_DSB_yield = 0;
112 Double_t SD_DSBp_yield = 0;
113 Double_t SD_DSBpp_yield = 0;
114
115 Double_t sSSB_yield = 0;
116 Double_t SSBi_yield = 0;
117 Double_t SSBd_yield = 0;
118 Double_t SSBm_yield = 0;
119 Double_t SD_sSSB_yield = 0;
120 Double_t SD_SSBi_yield = 0;
121 Double_t SD_SSBd_yield = 0;
122 Double_t SD_SSBm_yield = 0;
123 Double_t sDSB_yield = 0;
124 Double_t DSBi_yield = 0;
125 Double_t DSBd_yield = 0;
126 Double_t DSBm_yield = 0;
127 Double_t DSBh_yield = 0;
128 Double_t SD_sDSB_yield = 0;
129 Double_t SD_DSBi_yield = 0;
130 Double_t SD_DSBd_yield = 0;
131 Double_t SD_DSBm_yield = 0;
132 Double_t SD_DSBh_yield = 0;
133
134 total_EB = 0;
135 total_ES = 0;
136 total_OHB = 0;
137 total_OHS = 0;
138 total_HB = 0;
139 total_HS = 0;
140
141 total_SSB = 0;
142 total_SSBp = 0;
143 total_SSB2p = 0;
144 total_SSB2 = 0;
145 total_SSBp2 = 0;
146 total_SSB2p2 = 0;
147 total_DSB = 0;
148 total_DSBp = 0;
149 total_DSBpp = 0;
150 total_DSB2 = 0;
151 total_DSBp2 = 0;
152 total_DSBpp2 = 0;
153
154 total_sSSB = 0;
155 total_SSBd = 0;
156 total_SSBi = 0;
157 total_SSBm = 0;
158 total_sSSB2 = 0;
159 total_SSBd2 = 0;
160 total_SSBi2 = 0;
161 total_SSBm2 = 0;
162 total_sDSB = 0;
163 total_DSBd = 0;
164 total_DSBi = 0;
165 total_DSBm = 0;
166 total_DSBh = 0;
167 total_sDSB2 = 0;
168 total_DSBd2 = 0;
169 total_DSBi2 = 0;
170 total_DSBm2 = 0;
171 total_DSBh2 = 0;
172
173 Double_t eVtoJ = 1.60218e-19;
174 Double_t EnergyDeposited_eV = 0;
175 Double_t acc_edep = 0;
176 Double_t acc_edep2 = 0;
177
178 Double_t Energy;
179 Double_t BPID;
180 Int_t Strand;
181 Int_t StrandDamage;
182 Int_t Event1, directB;
183 Char_t Primary;
184 char* primaryName = new char[32];
185 char* type = new char[256];
186 char* sourceClassification = new char[256];
187
188 Double_t xx, yy, zz;
189
190 // Read trees and leaves from root file, and give values to variables
191 TTree* tree = (TTree*)f->Get("tuples/primary_source");
192 Float_t number = (Float_t)tree->GetEntries();
193
194 vector<pair<int, int64_t>> DSBBPID;
195
196 // For reading species production
197 tree = (TTree*)f->Get("tuples/damage");
198 tree->SetBranchAddress("Primary", &Primary);
199 tree->SetBranchAddress("Energy", &Energy);
200 tree->SetBranchAddress("EaqBaseHits", &EB);
201 tree->SetBranchAddress("EaqStrandHits", &ES);
202 tree->SetBranchAddress("OHBaseHits", &OHB);
203 tree->SetBranchAddress("OHStrandHits", &OHS);
204 tree->SetBranchAddress("HBaseHits", &HB);
205 tree->SetBranchAddress("HStrandHits", &HS);
206 tree->SetBranchAddress("TypeClassification", type);
207 tree->SetBranchAddress("BasePair", &BPID);
208 tree->SetBranchAddress("Event", &Event1);
209 tree->SetBranchAddress("Strand", &Strand);
210 tree->SetBranchAddress("StrandDamage", &StrandDamage);
211 tree->SetBranchAddress("Position_x_um", &xx);
212 tree->SetBranchAddress("Position_y_um", &yy);
213 tree->SetBranchAddress("Position_z_um", &zz);
214 tree->SetBranchAddress("DirectBreaks", &directB);
215
216 int primea = 0;
217 Long64_t nentries = tree->GetEntries();
218 int damagecount = 0;
219 int hitcount = 0;
220 // Double_t plasmidsD [25992];
221 Double_t plasmidsD[2599200];
222 for (int i = 0; i < nentries; i++) {
223 tree->GetEntry(i);
224 total_EB += EB;
225 total_EB2 += pow(EB, 2);
226 total_ES += ES;
227 total_ES2 += pow(ES, 2);
228 total_OHB += OHB;
229 total_OHB2 += pow(OHB, 2);
230 total_OHS += OHS;
231 total_OHS2 += pow(OHS, 2);
232 total_HB += HB;
233 total_HB2 += pow(HB, 2);
234 total_HS += HS;
235 total_HS2 += pow(HS, 2);
236 if ((string)type == "DSB" || (string)type == "DSB+" || (string)type == "DSB++") {
237 DSBBPID.push_back(make_pair(i, (int64_t)BPID));
238 }
239 if (StrandDamage != 0) {
240 {
241 damagecount++;
242 h1count->Fill(Strand);
243 }
244 primea++;
245 }
246 }
247 int unbrokenP = 0;
248 int brokenP = 0;
249 for (int i = 0; i < 11000; i++) {
250 plasmidsD[i] = h1count->GetBinContent(i);
251 if (plasmidsD[i] == 0 && i < numberOfPlasmids) unbrokenP++;
252 }
253
254 Double_t totalDs = 0.0;
255 for (int i = 0; i < 11000; i++) {
256 if (plasmidsD[i] != 0) {
257 h4damage->Fill(plasmidsD[i]);
258 totalDs += plasmidsD[i];
259 brokenP++;
260 }
261 }
262
263 Double_t X;
264 for (int i = 0; i < 28; i++) {
265 X = h4damage->GetBinContent(i);
266 h4damage->SetBinContent(i, X * 100 / totalDs);
267 }
268
269 // Sort DSBs from the one with lower ID value to the one with higher ID value
270 // Then find the number of fragments that have been produced
271 sort(DSBBPID.begin(), DSBBPID.end(), smallerPair);
272 for (int ie = 0; ie < DSBBPID.size() - 1; ie++) {
273 int64_t dsbfragment = DSBBPID[ie + 1].second - DSBBPID[ie].second;
274 }
275
276 // Calculate the standard deviation of species
277 SD_EB = sqrt(((total_EB2 / number) - pow(total_EB / number, 2)) / (number - 1));
278 SD_ES = sqrt(((total_ES2 / number) - pow(total_ES / number, 2)) / (number - 1));
279 SD_OHB = sqrt(((total_OHB2 / number) - pow(total_OHB / number, 2)) / (number - 1));
280 SD_OHS = sqrt(((total_OHS2 / number) - pow(total_OHS / number, 2)) / (number - 1));
281 SD_HB = sqrt(((total_HB2 / number) - pow(total_HB / number, 2)) / (number - 1));
282 SD_HS = sqrt(((total_HS2 / number) - pow(total_HS / number, 2)) / (number - 1));
283
284 // Read damage classification SSB, SSB+, 2SSB, DSB, DSB+, DSB++
285 // As they have been defined in: Nikjoo, H., O’Neill, O., Goodhead, T., & Terrissol, M. 1997,
286 // Computational modelling of low-energy electron-induced DNA damage by early physical
287 // and chemical events, International Journal of Radiation Biology, 71, 467.
288 tree = (TTree*)f->Get("tuples/classification");
289 tree->SetBranchAddress("Primary", &Primary);
290 tree->SetBranchAddress("Energy", &Energy);
291 tree->SetBranchAddress("SSB", &SSB);
292 tree->SetBranchAddress("SSBp", &SSBp);
293 tree->SetBranchAddress("2SSB", &SSB2p);
294 tree->SetBranchAddress("DSB", &DSB);
295 tree->SetBranchAddress("DSBp", &DSBp);
296 tree->SetBranchAddress("DSBpp", &DSBpp);
297
298 Long64_t nentriesC = tree->GetEntries();
299 for (int i = 0; i < nentriesC; i++) {
300 tree->GetEntry(i);
301
302 total_SSBp += SSBp;
303 total_SSBp2 += pow(SSBp, 2);
304 total_SSB2p += SSB2p;
305 total_SSB2p2 += pow(SSB2p, 2);
306 total_SSB += SSB;
307 total_SSB2 += pow(SSB, 2);
308
309 total_DSBp += DSBp;
310 total_DSBp2 += pow(DSBp, 2);
311 total_DSBpp += DSBpp;
312 total_DSBpp2 += pow(DSBpp, 2);
313 total_DSB += DSB;
314 total_DSB2 += pow(DSB, 2);
315 }
316
317 // Calculate the standard deviation
318 SD_SSB = sqrt(((total_SSB2 / number) - pow(total_SSB / number, 2)) / (number - 1));
319 SD_SSBp = sqrt(((total_SSBp2 / number) - pow(total_SSBp / number, 2)) / (number - 1));
320 SD_SSB2p = sqrt(((total_SSB2p2 / number) - pow(total_SSB2p / number, 2)) / (number - 1));
321
322 SD_DSB = sqrt(((total_DSB2 / number) - pow(total_DSB / number, 2)) / (number - 1));
323 SD_DSBp = sqrt(((total_DSBp2 / number) - pow(total_DSBp / number, 2)) / (number - 1));
324 SD_DSBpp = sqrt(((total_DSBpp2 / number) - pow(total_DSBpp / number, 2)) / (number - 1));
325
326 // Read damage classification SSBd, SSBi, SSBm, DSBd, DSBi, DSBm, DSBh
327 // As they have been defined in: Nikjoo, H., O’Neill, O., Goodhead, T., & Terrissol, M. 1997,
328 // Computational modelling of low-energy electron-induced DNA damage by early physical
329 // and chemical events, International Journal of Radiation Biology, 71, 467.
330 tree = (TTree*)f->Get("tuples/source");
331 tree->SetBranchAddress("Primary", primaryName);
332 tree->SetBranchAddress("Energy", &Energy);
333 tree->SetBranchAddress("SSBd", &SSBd);
334 tree->SetBranchAddress("SSBi", &SSBi);
335 tree->SetBranchAddress("SSBm", &SSBm);
336 tree->SetBranchAddress("DSBd", &DSBd);
337 tree->SetBranchAddress("DSBi", &DSBi);
338 tree->SetBranchAddress("DSBm", &DSBm);
339 tree->SetBranchAddress("DSBh", &DSBh);
340
341 int iprime = 0;
342
343 Long64_t nentriesS = tree->GetEntries();
344 for (int i = 0; i < nentriesS; i++) {
345 tree->GetEntry(i);
346
347 total_SSBd += SSBd;
348 total_SSBd2 += pow((SSBd), 2);
349 total_SSBi += SSBi;
350 total_SSBi2 += pow((SSBi), 2);
351 total_SSBm += SSBm;
352 total_SSBm2 += pow((SSBm), 2);
353 total_sSSB += SSBd + SSBi + SSBm;
354 total_sSSB2 += pow((SSBd + SSBi + SSBm), 2);
355
356 total_DSBd += DSBd;
357 total_DSBd2 += pow(DSBd, 2);
358 total_DSBi += DSBi;
359 total_DSBi2 += pow(DSBi, 2);
360 total_DSBm += DSBm;
361 total_DSBm2 += pow(DSBm, 2);
362 total_DSBh += DSBh;
363 total_DSBh2 += pow(DSBh, 2);
364 total_sDSB += DSBd + DSBi + DSBm + DSBh;
365 total_sDSB2 += pow((DSBd + DSBi + DSBm + DSBh), 2);
366
367 if (SSBd != 0) h1SSB->Fill(SSBd);
368 if (DSBd != 0) h1DSB->Fill(DSBd);
369 if (SSBd != 0) h1damage->Fill(SSBd);
370 if (DSBd != 0) h1damage->Fill(DSBd);
371 if (SSB != 0 || DSBd != 0) {
372 iprime++;
373 }
374 }
375
376 Double_t Y;
377 for (int i = 0; i < 28; i++) {
378 Y = h1damage->GetBinContent(i);
379 h1damage->SetBinContent(i, Y * 100 / totalDs);
380 }
381
382 // Calculate the standard deviation
383 SD_sSSB = sqrt(((total_sSSB2 / number) - pow(total_sSSB / number, 2)) / (number - 1));
384 SD_SSBd = sqrt(((total_SSBd2 / number) - pow(total_SSBd / number, 2)) / (number - 1));
385 SD_SSBi = sqrt(((total_SSBi2 / number) - pow(total_SSBi / number, 2)) / (number - 1));
386 SD_SSBm = sqrt(((total_SSBm2 / number) - pow(total_SSBm / number, 2)) / (number - 1));
387
388 SD_sDSB = sqrt(((total_sDSB2 / number) - pow(total_sDSB / number, 2)) / (number - 1));
389 SD_DSBd = sqrt(((total_DSBd2 / number) - pow(total_DSBd / number, 2)) / (number - 1));
390 SD_DSBi = sqrt(((total_DSBi2 / number) - pow(total_DSBi / number, 2)) / (number - 1));
391 SD_DSBm = sqrt(((total_DSBm2 / number) - pow(total_DSBm / number, 2)) / (number - 1));
392 SD_DSBh = sqrt(((total_DSBh2 / number) - pow(total_DSBh / number, 2)) / (number - 1));
393
394 // Measure the Deposited Energy in the whole volume that includes DNA chain
395
396 tree = (TTree*)f->Get("tuples/chromosome_hits");
397 tree->SetBranchAddress("e_chromosome_kev", &EnergyDeposited_eV);
398 nentries = tree->GetEntries();
399 for (int i = 0; i < nentries; i++) {
400 tree->GetEntry(i);
401 acc_edep += EnergyDeposited_eV * 1e3;
402 acc_edep2 += EnergyDeposited_eV * EnergyDeposited_eV * 1e6;
403 }
404 tree->SetBranchAddress("e_dna_kev", &EnergyDeposited_eV);
405 nentries = tree->GetEntries();
406 for (int i = 0; i < nentries; i++) {
407 tree->GetEntry(i);
408 acc_edep += (EnergyDeposited_eV * 1e3);
409 acc_edep2 += (EnergyDeposited_eV * EnergyDeposited_eV * 1e6);
410 }
411
412 // Close the root file to free space
413 f->Close();
414 // Calculate the absorbed dose
415 dose = acc_edep * eVtoJ / mass;
416
417 double norm = 1;
418 // Calculate the yields, together with their standard deviation
419 EB_yield = (Double_t)total_EB / dose / Nbp;
420 ES_yield = (Double_t)total_ES / dose / Nbp;
421 OHB_yield = (Double_t)total_OHB / dose / Nbp;
422 OHS_yield = (Double_t)total_OHS / dose / Nbp;
423 HB_yield = (Double_t)total_HB / dose / Nbp;
424 HS_yield = (Double_t)total_HS / dose / Nbp;
425
426 SD_EB_yield = SD_EB / dose / Nbp;
427 SD_ES_yield = SD_ES / dose / Nbp;
428 SD_OHB_yield = SD_OHB / dose / Nbp;
429 SD_OHS_yield = SD_OHS / dose / Nbp;
430 SD_HB_yield = SD_HB / dose / Nbp;
431 SD_HS_yield = SD_HS / dose / Nbp;
432
433 SSB_yield = (Double_t)norm * total_SSB / dose / Nbp;
434 SSBp_yield = (Double_t)norm * total_SSBp / dose / Nbp;
435 SSB2p_yield = (Double_t)norm * total_SSB2p / dose / Nbp;
436
437 DSB_yield = (Double_t)norm * total_DSB / dose / Nbp;
438 DSBp_yield = (Double_t)norm * total_DSBp / dose / Nbp;
439 DSBpp_yield = (Double_t)norm * total_DSBpp / dose / Nbp;
440
441 SD_SSB_yield = norm * SD_SSB / dose / Nbp;
442 SD_SSBp_yield = norm * SD_SSBp / dose / Nbp;
443 SD_SSB2p_yield = norm * SD_SSB2p / dose / Nbp;
444
445 SD_DSB_yield = norm * SD_DSB / dose / Nbp;
446 SD_DSBp_yield = norm * SD_DSBp / dose / Nbp;
447 SD_DSBpp_yield = norm * SD_DSBpp / dose / Nbp;
448
449 sSSB_yield = (Double_t)norm * total_sSSB / dose / Nbp;
450 SSBi_yield = (Double_t)norm * total_SSBi / dose / Nbp;
451 SSBd_yield = (Double_t)norm * total_SSBd / dose / Nbp;
452 SSBm_yield = (Double_t)norm * total_SSBm / dose / Nbp;
453
454 sDSB_yield = (Double_t)norm * total_sDSB / dose / Nbp;
455 DSBi_yield = (Double_t)norm * total_DSBi / dose / Nbp;
456 DSBd_yield = (Double_t)norm * total_DSBd / dose / Nbp;
457 DSBm_yield = (Double_t)norm * total_DSBm / dose / Nbp;
458 DSBh_yield = (Double_t)norm * total_DSBh / dose / Nbp;
459
460 SD_sSSB_yield = norm * SD_sSSB / dose / Nbp;
461 SD_SSBi_yield = norm * SD_SSBi / dose / Nbp;
462 SD_SSBd_yield = norm * SD_SSBd / dose / Nbp;
463 SD_SSBm_yield = norm * SD_SSBm / dose / Nbp;
464
465 SD_sDSB_yield = norm * SD_sDSB / dose / Nbp;
466 SD_DSBi_yield = norm * SD_DSBi / dose / Nbp;
467 SD_DSBd_yield = norm * SD_DSBd / dose / Nbp;
468 SD_DSBm_yield = norm * SD_DSBm / dose / Nbp;
469 SD_DSBh_yield = norm * SD_DSBh / dose / Nbp;
470
471 // Print output in terminal
472
473 float total_SSB_totalYield = SSB_yield + SSBp_yield + SSB2p_yield;
474 float total_DSB_totalYield = DSB_yield + DSBp_yield + DSBpp_yield;
475
476 cout << "\n"
477 << " Output file: " << ifile << '\n'
478 << "\nDose Absorbed (Gy): " << dose << '\n'
479 << "Particle : " << primaryName << '\t' << "Energy (MeV) : " << Energy << '\t'
480 << "Number of Primaries : " << number << '\n'
481
482 << " Output Damage : " << '\n'
483
484 << '\t' << "Number of plasmids: " << numberOfPlasmids << " \t" << '\n'
485 << '\t' << "Unbroken plasmids: " << unbrokenP << " \t" << '\n'
486 << '\t' << "Broken plasmids: " << brokenP << " \t" << '\n'
487 << '\t' << "Rate of damages per plasmid " << totalDs / numberOfPlasmids << " \t" << '\n'
488
489 << '\t' << "Rate of damages per plasmid per dose " << '\t'
490 << totalDs / numberOfPlasmids / dose << " plasmid-1Gy-1" << '\n'
491 << '\t' << "Rate of damages per dose per Mbp " << totalDs / dose / Nbp << " Gy-1Mbp-1"
492 << '\n'
493
494 << '\t' << " Species Hits " << '\n'
495 << '\t' << "EaqBaseHits " << EB_yield * dose * Nbp << " \t" << '\n'
496 << '\t' << "EaqStrandHits " << ES_yield * dose * Nbp << " \t" << '\n'
497 << '\t' << "OHBaseHits " << OHB_yield * dose * Nbp << " \t" << '\n'
498 << '\t' << "OHStrandHits " << OHS_yield * dose * Nbp << " \t" << '\n'
499 << '\t' << "HBaseHits " << HB_yield * dose * Nbp << " \t" << '\n'
500 << '\t' << "HStrandHits " << HS_yield * dose * Nbp << " \t" << '\n'
501 << '\n'
502 << '\t' << " Damage number" << '\n'
503 << '\t' << "SSB " << SSB_yield * dose * Nbp << " \t" << '\n'
504 << '\t' << "SSB+ " << SSBp_yield * dose * Nbp << " \t" << '\n'
505 << '\t' << "2SSB " << SSB2p_yield * dose * Nbp << " \t" << '\n'
506 << '\t' << "SSB total " << total_SSB_totalYield * dose * Nbp << '\n'
507 << '\t' << "DSB " << DSB_yield * dose * Nbp << " \t" << '\n'
508 << '\t' << "DSB+ " << DSBp_yield * dose * Nbp << " \t" << '\n'
509 << '\t' << "DSB++ " << DSBpp_yield * dose * Nbp << " \t" << '\n'
510 << '\t' << "DSB total " << total_DSB_totalYield * dose * Nbp << '\n'
511 << '\n'
512 << '\t' << " Breaks number " << '\n'
513 << '\t' << "SSB direct " << SSBd_yield * dose * Nbp << " \t" << '\n'
514 << '\t' << "SSB indirect " << SSBi_yield * dose * Nbp << " \t" << '\n'
515 << '\t' << "SSB mixed " << SSBm_yield * dose * Nbp << " \t" << '\n'
516 << '\t' << "SSB total " << sSSB_yield * dose * Nbp << " \t" << '\n'
517 << '\t' << "DSB direct " << DSBd_yield * dose * Nbp << " \t" << '\n'
518 << '\t' << "DSB indirect " << DSBi_yield * dose * Nbp << " \t" << '\n'
519 << '\t' << "DSB mixed " << DSBm_yield * dose * Nbp << " \t" << '\n'
520 << '\t' << "DSB hybrid " << DSBh_yield * dose * Nbp << " \t" << '\n'
521 << '\t' << "DSB total " << sDSB_yield * dose * Nbp << " \t" << '\n'
522 << '\n'
523 << '\t' << "SSB/DSB " << sSSB_yield / sDSB_yield << '\n'
524 << '\n';
525
526 // Plot Histograms
527
528 cSSB->GetCanvas()->cd();
529 h1SSB->SetStats(false);
530 h1SSB->SetMarkerSize(0.1);
531 h1SSB->SetMarkerColor(kBlue);
532 h1SSB->SetLineColor(kBlue);
533 h1SSB->SetTitle("");
534 h1SSB->SetYTitle(" ");
535 h1SSB->SetXTitle("Number of direct SSB per event");
536 h1SSB->SetFillColor(kBlue);
537 h1SSB->Draw();
538
539 cDSB->GetCanvas()->cd();
540 h1DSB->SetStats(false);
541 h1DSB->SetMarkerSize(0.1);
542 h1DSB->SetMarkerColor(kGreen + 2);
543 h1DSB->SetLineColor(kGreen + 2);
544 h1DSB->SetTitle("");
545 h1DSB->SetYTitle(" ");
546 h1DSB->SetXTitle("Number of direct DSB per event");
547 h1DSB->SetFillColor(kGreen + 2);
548 h1DSB->Draw();
549
550 cdamage->GetCanvas()->cd();
551 h1damage->SetStats(false);
552 h1damage->SetMarkerSize(0.1);
553 h1damage->SetMarkerColor(kMagenta);
554 h1damage->SetLineColor(kMagenta);
555 h1damage->SetTitle("");
556 h1damage->SetYTitle("Percentage % ");
557 h1damage->SetXTitle("Number of damages per event");
558 h1damage->SetFillColor(kMagenta);
559 h1damage->Draw();
560
561 c4damage->GetCanvas()->cd();
562 h4damage->SetStats(false);
563 h4damage->SetMarkerSize(0.1);
564 h4damage->SetMarkerColor(kRed - 4);
565 h4damage->SetLineColor(kRed - 4);
566 h4damage->SetTitle("");
567 h4damage->SetYTitle("Percentage % ");
568 h4damage->SetXTitle("Number of damages per plasmid");
569 h4damage->SetFillColor(kRed - 4);
570 h4damage->Draw();
571
572 ccount->GetCanvas()->cd();
573 h1count->SetStats(false);
574 h1count->SetMarkerSize(0.1);
575 h1count->SetMarkerColor(kCyan);
576 h1count->SetLineColor(kCyan);
577 h1count->SetTitle("");
578 h1count->SetYTitle("Number of damages");
579 h1count->SetXTitle("Plasmid ID");
580 h1count->Draw();
581 }
582
583 // Some important bools that are needed to run the root macro file
584 bool greaterPair(const ipair& l, const ipair& r)
585 {
586 return l.second > r.second;
587 }
588 bool smallerPair(const ipair& l, const ipair& r)
589 {
590 return l.second < r.second;
591 }
592
593 void BinLogX(TH1* h)
594 {
595 TAxis* axis = h->GetXaxis();
596 int bins = axis->GetNbins();
597 Axis_t from = axis->GetXmin();
598 Axis_t to = axis->GetXmax();
599 Axis_t width = (to - from) / bins;
600 Axis_t* new_bins = new Axis_t[bins + 1];
601 for (int i = 0; i <= bins; i++) {
602 new_bins[i] = TMath::Power(10, from + i * width);
603 }
604 axis->Set(bins, new_bins);
605 delete[] new_bins;
606 }
607