Geant4 Cross Reference |
1 //-------------------------------------------- 1 2 // This macrofile was developed by Konstantino 3 // in collaboration with the whole team of mol 4 // Publication: .............................. 5 // For any question please contact through: 6 // chatzipa@cenbg.in2p3.fr (or k.chatzipapas@y 7 //-------------------------------------------- 8 9 // This macro requires the molecular-dna.root 10 // To run this file just insert this command t 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, b 22 system ("hadd -O -f molecular-dna.root molecul 23 24 // Define these parameters of the simulation 25 char ifile[256] = "molecular-dna.root"; // in 26 Double_t r3 = 10575e-9 * 3450e-9 * 10575e-9; 27 Double_t Nbp = 6383.848410; // Mbp // Length o 28 Double_t mass = 997 * 4 * 3.141592 * r3 / 3 ; 29 ////////////////////////////////////////////// 30 //******************************************** 31 32 typedef std::pair <int64_t, int64_t> ipair; 33 bool greaterPair(const ipair &l, const ipair & 34 bool smallerPair(const ipair &l, const ipair & 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"," 45 cfragment->SetLogx(); 46 cfragment->SetLogy(); 47 TH1F *h1fragments = new TH1F("h1fragments","h1 48 BinLogX(h1fragments); 49 50 TCanvas *c1 = new TCanvas("c1", "Molecular DNA 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.5 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 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 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 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 101 Float_t total_EB2, total_ES2, total_OHB2, tota 102 Float_t SD_EB, SD_ES, SD_OHB, SD_OHS, SD_HB, S 103 Float_t SD_SSB, SD_SSBp, SD_SSB2p, SD_sSSB, SD 104 Float_t SD_DSB, SD_DSBp, SD_DSBpp, SD_sDSB, SD 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, tota 115 Float_t total_sSSB2, total_SSBd2, total_SSBi2, 116 Int_t DSBd, DSBi, DSBm, DSBh; 117 Int_t total_sDSB, total_DSBd, total_DSBi, tota 118 Float_t total_sDSB2, total_DSBd2, total_DSBi2, 119 120 Double_t dose = 0; 121 Double_t SD_dose = 0; 122 123 Double_t EB_yield = 0; Double_t ES_yield = 0; 124 Double_t SD_EB_yield = 0; Double_t SD_ES_yield 125 126 Double_t SSB_yield = 0; Double_t SSBp_yield = 127 Double_t SD_SSB_yield = 0; Double_t SD_SSBp_yi 128 Double_t DSB_yield = 0; Double_t DSBp_yield = 129 Double_t SD_DSB_yield = 0; Double_t SD_DSBp_yi 130 131 Double_t sSSB_yield = 0; Double_t SSBi_yield = 132 Double_t SD_sSSB_yield = 0; Double_t SD_SSBi_y 133 Double_t sDSB_yield = 0; Double_t DSBi_yield = 134 Double_t SD_sDSB_yield = 0; Double_t SD_DSBi_y 135 136 total_EB = 0; total_ES = 0; total_OHB = 0; to 137 138 total_SSB = 0; total_SSBp = 0; total_SSB2p = 139 total_SSB2 = 0; total_SSBp2 = 0; total_SSB2p2 140 total_DSB = 0; total_DSBp = 0; total_DSBpp = 141 total_DSB2 = 0; total_DSBp2 = 0; total_DSBpp2 142 143 total_sSSB = 0; total_SSBd = 0; total_SSBi = 0 144 total_sSSB2 = 0; total_SSBd2 = 0; total_SSBi2 145 total_sDSB = 0; total_DSBd = 0; total_DSBi = 0 146 total_sDSB2 = 0; total_DSBd2 = 0; total_DSBi2 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 g 160 TTree* tree = (TTree*) f->Get("tuples/primary_ 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", &P 168 tree->SetBranchAddress("Energy", &E 169 tree->SetBranchAddress("EaqBaseHits", &E 170 tree->SetBranchAddress("EaqStrandHits", &E 171 tree->SetBranchAddress("OHBaseHits", &O 172 tree->SetBranchAddress("OHStrandHits", &O 173 tree->SetBranchAddress("HBaseHits", &H 174 tree->SetBranchAddress("HStrandHits", &H 175 tree->SetBranchAddress("TypeClassification", t 176 tree->SetBranchAddress("BasePair", &B 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+" 197 //cout << "DSB:"<<type<<endl; 198 DSBBPID.push_back(make_pair(i,(int64_t)BPI 199 } 200 201 } 202 203 // Sort DSBs from the one with lower ID value 204 // Then find the number of fragments that have 205 sort(DSBBPID.begin(), DSBBPID.end(), smallerP 206 for(int ie = 0;ie<DSBBPID.size()-1;ie++){ 207 int64_t dsbfragment = DSBBPID[ie+1].second-D 208 209 double val = (double)dsbfragment/1000. 210 double meanw = h1fragments->GetBinCenter 211 double binw = h1fragments->GetBinWidth 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(tot 218 SD_ES = sqrt(((total_ES2 / number) - pow(tot 219 SD_OHB = sqrt(((total_OHB2 / number) - pow(tot 220 SD_OHS = sqrt(((total_OHS2 / number) - pow(tot 221 SD_HB = sqrt(((total_HB2 / number) - pow(tot 222 SD_HS = sqrt(((total_HS2 / number) - pow(tot 223 224 // Read damage classification SSB, SSB+, 2SSB, 225 // As they have been defined in: Nikjoo, H., O 226 // Computational modelling of low-energy elect 227 // and chemical events, International Journal 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 261 SD_SSBp = sqrt(((total_SSBp2 / number) - pow 262 SD_SSB2p = sqrt(((total_SSB2p2 / number) - pow 263 264 SD_DSB = sqrt(((total_DSB2 / number) - pow 265 SD_DSBp = sqrt(((total_DSBp2 / number) - pow 266 SD_DSBpp = sqrt(((total_DSBpp2 / number) - pow 267 268 // Read damage classification SSBd, SSBi, SSBm 269 // As they have been defined in: Nikjoo, H., O 270 // Computational modelling of low-energy elect 271 // and chemical events, International Journal 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(t 311 SD_SSBd = sqrt(((total_SSBd2 / number) - pow(t 312 SD_SSBi = sqrt(((total_SSBi2 / number) - pow(t 313 SD_SSBm = sqrt(((total_SSBm2 / number) - pow(t 314 315 SD_sDSB = sqrt(((total_sDSB2 / number) - pow(t 316 SD_DSBd = sqrt(((total_DSBd2 / number) - pow(t 317 SD_DSBi = sqrt(((total_DSBi2 / number) - pow(t 318 SD_DSBm = sqrt(((total_DSBm2 / number) - pow(t 319 SD_DSBh = sqrt(((total_DSBh2 / number) - pow(t 320 321 322 // Measure the Deposited Energy in the whole v 323 324 tree = (TTree *) f->Get("tuples/chromosome_hit 325 tree->SetBranchAddress("e_chromosome_kev",&Ene 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 *EnergyDepos 331 } 332 tree->SetBranchAddress("e_dna_kev",&EnergyDepo 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 *EnergyDepos 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 t 347 // Default value is 1 to produce the result in 348 // It changes Mbp to Gbp. Some other changes m 349 double norm = 1000; 350 351 // Calculate the yields, together with their s 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 / 368 SSBp_yield = (Double_t) norm * total_SSBp / 369 SSB2p_yield = (Double_t) norm * total_SSB2p / 370 371 DSB_yield = (Double_t) norm * total_DSB / 372 DSBp_yield = (Double_t) norm * total_DSBp / 373 DSBpp_yield = (Double_t) norm * total_DSBpp / 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 / do 385 SSBi_yield = (Double_t) norm * total_SSBi / do 386 SSBd_yield = (Double_t) norm * total_SSBd / do 387 SSBm_yield = (Double_t) norm * total_SSBm / do 388 389 sDSB_yield = (Double_t) norm * total_sDSB / do 390 DSBi_yield = (Double_t) norm * total_DSBi / do 391 DSBd_yield = (Double_t) norm * total_DSBd / do 392 DSBm_yield = (Double_t) norm * total_DSBm / do 393 DSBh_yield = (Double_t) norm * total_DSBh / do 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_ 410 float total_DSB_totalYield = DSB_yield + DSBp_ 411 412 cout<<"\n" <<ifile 413 <<"\nDose Absorbed (Gy): " <<dose 414 <<"Particle : " <<primaryName 415 <<"Energy (MeV) : " <<Energy 416 <<"Number of Primaries : " <<number 417 <<" Output Damage : " 418 <<" Species Hits (Gy-1 Mbp-1) " < 419 <<"EaqBaseHits : " <<EB_yield << 420 <<"EaqStrandHits : " <<ES_yield << 421 <<"OHBaseHits : " <<OHB_yield << 422 <<"OHStrandHits : " <<OHS_yield << 423 <<"HBaseHits : " <<HB_yield << 424 <<"HStrandHits : " <<HS_yield << 425 <<" Damage yield (Gy-1 Gbp-1) " < 426 <<"SSB : " <<SSB_yield << 427 <<"SSB+ : " <<SSBp_yield << 428 <<"2SSB : " <<SSB2p_yield << 429 <<"SSB total : " <<total_SSB_total 430 <<"DSB : " <<DSB_yield << 431 <<"DSB+ : " <<DSBp_yield << 432 <<"DSB++ : " <<DSBpp_yield << 433 <<"DSB total : " <<total_DSB_total 434 <<" Breaks yield (Gy-1 Gbp-1) " < 435 <<"SSB direct : " <<SSBd_yield << 436 <<"SSB indirect : " <<SSBi_yield << 437 <<"SSB mixed : " <<SSBm_yield << 438 <<"SSB total : " <<sSSB_yield << 439 <<"DSB direct : " <<DSBd_yield << 440 <<"DSB indirect : " <<DSBi_yield << 441 <<"DSB mixed : " <<DSBm_yield << 442 <<"DSB hybrid : " <<DSBh_yield << 443 <<"DSB total : " <<sDSB_yield << 444 <<"SSB/DSB : " <<sSSB_yield/sDSB 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 (b 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,O 469 Double_t err_y[n] = {SD_EB_yield,SD_ES_yield,S 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 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_y 488 Double_t err_y2[n] = {SD_SSBp_yield,SD_SSB2p_y 489 TGraph* gr2 = new TGraphErrors(n,x2,y2,0,err_y 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^ 498 gr2->GetYaxis()->SetTitle("Damage yield (Gy^{- 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_y 509 Double_t err_y3[m] = {SD_SSBd_yield,SD_SSBi_yi 510 TGraph* gr3 = new TGraphErrors(m,x3,y3,0,err_y 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^ 517 gr3->GetYaxis()->SetTitle("SSB yield (Gy^{-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_y 528 Double_t err_y4[k] = {SD_DSBd_yield,SD_DSBi_yi 529 TGraph* gr4 = new TGraphErrors(k,x4,y4,0,err_y 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^ 537 gr4->GetYaxis()->SetTitle("DSB yield (Gy^{-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 546 bool greaterPair(const ipair& l, const ipair& 547 bool smallerPair(const ipair& l, const ipair& 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 * 558 } 559 axis->Set(bins, new_bins); 560 delete[] new_bins; 561 562 } 563