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