dna/moleculardna/human_cell_alphas.C

0001 //-------------------------------------------------------------------------------//
0002 // This macrofile was developed by Konstantinos Chatzipapas at LP2iB (ex. CENBG) //
0003 // in collaboration with the whole team of molecularDNA Geant4-DNA example       //
0004 // Publication: K. Chatzipapas, et al., Phys. Med. 112 (2023) 102613             //
0005 // For any question please contact through:                                      //
0006 // chatzipa@cenbg.in2p3.fr                                                       //
0007 //-------------------------------------------------------------------------------//
0008
0009 // This macro requires the molecular-dna.root file generated from molecularDNA example
0010
0011 {
0012 //*******************************************************************************//
0013 // If you need to add multiple root outputs, by multithreading, use this command:
0014 // system ("hadd -O -f molecular-dna.root molecular-dna_t*.root");
0015
0016 // Define these parameters of the simulation
0017 char ifile[256] = "molecular-dna.root";  // input filepath to be replaced
0018
0019 //for HTB
0020 Double_t r3 = 8550e-9 * 2500e-9 * 6425e-9; // a * b * c
0021 Double_t Nbp = 6454.227202; // Mbp // Length of the DNA chain in Mbp
0022
0023 // for MCF
0024 //Double_t r3 = 7005e-9 * 2500e-9 * 5300e-9; // a * b * c
0025 //Double_t Nbp = 6424.831996; // Mbp // Length of the DNA chain in Mbp // for MCF
0026
0027 Double_t mass = 997 * 4 * 3.141592 * r3 / 3 ;  // density * 4/3 * pi * r3
0028 ///////////////////////////////////////////////////////////////////////////////////
0029 //*******************************************************************************//
0030
0031 typedef std::pair <int64_t, int64_t> ipair;
0032 bool greaterPair(const ipair &l, const ipair &r);
0033 bool smallerPair(const ipair &l, const ipair &r);
0034
0035 void BinLogX(TH1 *h);
0036
0037 gROOT->Reset();
0038 gStyle->SetPalette(1);
0039 gROOT->SetStyle("Plain");
0040 gStyle->SetOptStat(00000);
0041
0042 // Initialize output histograms
0043 TCanvas *cfragment = new TCanvas("cfragment","DNA Fragments Distribution", 900, 120, 600,400);
0044 cfragment->SetLogx();
0045 cfragment->SetLogy();
0046 TH1F *h1fragments = new TH1F("h1fragments","h1fragments",40,0,5);
0047 BinLogX(h1fragments);
0048
0049 TCanvas *c1 = new TCanvas("c1", "Molecular DNA - Damage Quantification", 60, 120, 800, 800);
0050 c1->SetBorderSize(0);
0051 c1->SetFillColor(0);
0052 c1->SetFillStyle(4000);
0053 gPad->SetLeftMargin(0.13);
0054
0055 TPad* pad1 = new TPad("pad1","Species", 0, 0.51, 0.49, 1);
0056 pad1->SetBorderSize(0);
0057 pad1->SetFillColor(0);
0058 pad1->SetFillStyle(4000);
0059 pad1->SetLeftMargin(0.15);
0060 pad1->SetRightMargin(0.01);
0061 pad1->SetBottomMargin(0.2);
0062
0063 TPad* pad2 = new TPad("pad2","Damage Yield", 0.51, 0.5, 1, 1);
0064 pad2->SetBorderSize(0);
0065 pad2->SetFillColor(0);
0066 pad2->SetFillStyle(4000);
0067 pad2->SetLeftMargin(0.15);
0068 pad2->SetRightMargin(0.05);
0069 pad2->SetBottomMargin(0.2);
0070
0071 TPad* pad3 = new TPad("pad3","Breaks Yield SSB", 0, 0, 0.49, 0.49);
0072 pad3->SetBorderSize(0);
0073 pad3->SetFillColor(0);
0074 pad3->SetFillStyle(4000);
0075 pad3->SetLeftMargin(0.15);
0076 pad3->SetRightMargin(0.01);
0077 //pad3->SetTopMargin(0.2);
0078 pad3->SetBottomMargin(0.2);
0079
0080 TPad* pad4 = new TPad("pad4","Breaks Yield DSB", 0.51, 0, 1, 0.49);
0081 pad4->SetBorderSize(0);
0082 pad4->SetFillColor(0);
0083 pad4->SetFillStyle(4000);
0084 pad4->SetLeftMargin(0.15);
0085 pad4->SetRightMargin(0.05);
0086 //pad3->SetTopMargin(0.2);
0087 pad4->SetBottomMargin(0.2);
0088
0089 pad1->Draw();
0090 pad2->Draw();
0091 pad3->Draw();
0092 pad4->Draw();
0093
0094 // Open root file
0095 TFile *f = TFile::Open(ifile);
0096
0097 // Initialize Variables
0098 Int_t EB, ES, OHB, OHS, HB, HS, FL;
0099 Int_t total_EB, total_ES, total_OHB, total_OHS, total_HB, total_HS, total_FL;
0100 Float_t total_EB2, total_ES2, total_OHB2, total_OHS2, total_HB2, total_HS2, total_FL2;
0101 Float_t SD_EB, SD_ES, SD_OHB, SD_OHS, SD_HB, SD_HS;
0102 Float_t SD_SSB, SD_SSBp, SD_SSB2p, SD_sSSB, SD_SSBd, SD_SSBi;
0103 Float_t SD_DSB, SD_DSBp, SD_DSBpp, SD_sDSB, SD_DSBd, SD_DSBi, SD_DSBm, SD_DSBh;
0104
0105 Int_t SSB, SSBp, SSB2p;
0106 Int_t total_SSB, total_SSBp, total_SSB2p;
0107 Float_t total_SSB2, total_SSBp2, total_SSB2p2;
0108 Int_t DSB, DSBp, DSBpp;
0109 Int_t total_DSB, total_DSBp, total_DSBpp;
0110 Float_t total_DSB2, total_DSBp2, total_DSBpp2;
0111
0112 Int_t SSBd, SSBi, SSBm;
0113 Int_t total_sSSB, total_SSBd, total_SSBi, total_SSBm;
0114 Float_t total_sSSB2, total_SSBd2, total_SSBi2, total_SSBm2;
0115 Int_t DSBd, DSBi, DSBm, DSBh;
0116 Int_t total_sDSB, total_DSBd, total_DSBi, total_DSBm, total_DSBh;
0117 Float_t total_sDSB2, total_DSBd2, total_DSBi2, total_DSBm2, total_DSBh2;
0118
0119 Double_t dose = 0;
0120 Double_t SD_dose = 0;
0121
0122 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;
0123 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;
0124
0125 Double_t SSB_yield = 0; Double_t SSBp_yield = 0; Double_t SSB2p_yield = 0;
0126 Double_t SD_SSB_yield = 0; Double_t SD_SSBp_yield = 0; Double_t SD_SSB2p_yield = 0;
0127 Double_t DSB_yield = 0; Double_t DSBp_yield = 0; Double_t DSBpp_yield = 0;
0128 Double_t SD_DSB_yield = 0; Double_t SD_DSBp_yield = 0; Double_t SD_DSBpp_yield = 0;
0129
0130 Double_t sSSB_yield = 0; Double_t SSBi_yield = 0; Double_t SSBd_yield = 0; Double_t SSBm_yield = 0;
0131 Double_t SD_sSSB_yield = 0; Double_t SD_SSBi_yield = 0; Double_t SD_SSBd_yield = 0; Double_t SD_SSBm_yield = 0;
0132 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;
0133 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;
0134
0135 total_EB  = 0; total_ES = 0; total_OHB = 0; total_OHS = 0; total_HB = 0; total_HS = 0;
0136
0137 total_SSB  = 0; total_SSBp = 0; total_SSB2p = 0;
0138 total_SSB2  = 0; total_SSBp2 = 0; total_SSB2p2 = 0;
0139 total_DSB  = 0; total_DSBp = 0; total_DSBpp = 0;
0140 total_DSB2  = 0; total_DSBp2 = 0; total_DSBpp2 = 0;
0141
0142 total_sSSB = 0; total_SSBd = 0; total_SSBi = 0; total_SSBm = 0;
0143 total_sSSB2 = 0; total_SSBd2 = 0; total_SSBi2 = 0; total_SSBm2 = 0;
0144 total_sDSB = 0; total_DSBd = 0; total_DSBi = 0; total_DSBm = 0; total_DSBh = 0;
0145 total_sDSB2 = 0; total_DSBd2 = 0; total_DSBi2 = 0; total_DSBm2 = 0;  total_DSBh2 = 0;
0146
0147 Double_t eVtoJ = 1.60218e-19;
0148 Double_t EnergyDeposited_eV = 0;
0149 Double_t acc_edep = 0;
0150 Double_t acc_edep2 = 0;
0151
0152 Double_t Energy;
0153 Double_t BPID;
0154 Char_t Primary;
0155 char *primaryName = new char[32];
0156 char *type= new char[256];
0157
0158 // Read trees and leaves from root file, and give values to variables
0159 TTree* tree = (TTree*) f->Get("tuples/primary_source");
0160 Float_t number = (Float_t) tree->GetEntries();
0161
0162 vector<pair<int,int64_t>> DSBBPID;
0163
0164 // For reading species production
0165 tree = (TTree*) f->Get("tuples/damage");
0166 tree->SetBranchAddress("Primary",           &Primary);
0167 tree->SetBranchAddress("Energy",            &Energy);
0168 tree->SetBranchAddress("EaqBaseHits",       &EB);
0169 tree->SetBranchAddress("EaqStrandHits",     &ES);
0170 tree->SetBranchAddress("OHBaseHits",        &OHB);
0171 tree->SetBranchAddress("OHStrandHits",      &OHS);
0172 tree->SetBranchAddress("HBaseHits",         &HB);
0173 tree->SetBranchAddress("HStrandHits",       &HS);
0174 tree->SetBranchAddress("TypeClassification", type);
0175 tree->SetBranchAddress("BasePair",          &BPID);
0176
0177
0178 Long64_t nentries = tree->GetEntries();
0179 for(int i = 0;i<nentries;i++){
0180   tree->GetEntry(i);
0181
0182   total_EB   += EB;
0183   total_EB2  += pow(EB,2);
0184   total_ES   += ES;
0185   total_ES2  += pow(ES,2);
0186   total_OHB  += OHB;
0187   total_OHB2 += pow(OHB,2);
0188   total_OHS  += OHS;
0189   total_OHS2 += pow(OHS,2);
0190   total_HB   += HB;
0191   total_HB2  += pow(HB,2);
0192   total_HS   += HS;
0193   total_HS2  += pow(HS,2);
0194
0195   if((string)type=="DSB"||(string)type=="DSB+"||(string)type=="DSB++"){
0196     //cout << "DSB:"<<type<<endl;
0197     DSBBPID.push_back(make_pair(i,(int64_t)BPID));
0198     }
0199
0200   }
0201
0202 // Sort DSBs from the one with lower ID value to the one with higher ID value
0203 // Then find the number of fragments that have been produced
0204 sort(DSBBPID.begin(),  DSBBPID.end(), smallerPair);
0205 for(int ie = 0;ie<DSBBPID.size()-1;ie++){
0206   int64_t dsbfragment = DSBBPID[ie+1].second-DSBBPID[ie].second;
0207
0208   double val       = (double)dsbfragment/1000.;
0209   double meanw     = h1fragments->GetBinCenter(h1fragments->FindBin(val));
0210   double binw      = h1fragments->GetBinWidth (h1fragments->FindBin(val));
0211   h1fragments->Fill(val,1./binw/1000);//bp-1
0212   //cout <<"val:"<<val<<endl;
0213   }
0214
0215 // Calculate the standard deviation of species
0216 SD_EB  = sqrt(((total_EB2  / number) - pow(total_EB  / number,2))/(number -1));
0217 SD_ES  = sqrt(((total_ES2  / number) - pow(total_ES  / number,2))/(number -1));
0218 SD_OHB = sqrt(((total_OHB2 / number) - pow(total_OHB / number,2))/(number -1));
0219 SD_OHS = sqrt(((total_OHS2 / number) - pow(total_OHS / number,2))/(number -1));
0220 SD_HB  = sqrt(((total_HB2  / number) - pow(total_HB  / number,2))/(number -1));
0221 SD_HS  = sqrt(((total_HS2  / number) - pow(total_HS  / number,2))/(number -1));
0222
0223 // Read damage classification SSB, SSB+, 2SSB, DSB, DSB+, DSB++
0224 // As they have been defined in: Nikjoo, H., O’Neill, O., Goodhead, T., & Terrissol, M. 1997,
0225 // Computational modelling of low-energy electron-induced DNA damage by early physical
0226 // and chemical events, International Journal of Radiation Biology, 71, 467.
0227 tree = (TTree *) f->Get("tuples/classification");
0228 tree->SetBranchAddress("Primary",&Primary);
0229 tree->SetBranchAddress("Energy", &Energy);
0230 tree->SetBranchAddress("SSB",    &SSB);
0231 tree->SetBranchAddress("SSBp",   &SSBp);
0232 tree->SetBranchAddress("SSB2",   &SSB2p);
0233 tree->SetBranchAddress("DSB",    &DSB);
0234 tree->SetBranchAddress("DSBp",   &DSBp);
0235 tree->SetBranchAddress("DSBpp",  &DSBpp);
0236
0237
0238 Long64_t nentriesC = tree->GetEntries();
0239 for(int i = 0;i<nentriesC;i++){
0240   tree->GetEntry(i);
0241
0242   total_SSBp   += SSBp;
0243   total_SSBp2  += pow(SSBp,2);
0244   total_SSB2p  += SSB2p;
0245   total_SSB2p2 += pow(SSB2p,2);
0246   total_SSB    += SSB;
0247   total_SSB2   += pow(SSB,2);
0248
0249   total_DSBp   += DSBp;
0250   total_DSBp2  += pow(DSBp,2);
0251   total_DSBpp  += DSBpp;
0252   total_DSBpp2 += pow(DSBpp,2);
0253   total_DSB    += DSB;
0254   total_DSB2   += pow(DSB,2);
0255
0256   }
0257
0258 // Calculate the standard deviation
0259 SD_SSB   = sqrt(((total_SSB2   / number) - pow(total_SSB   / number,2))/(number -1));
0260 SD_SSBp  = sqrt(((total_SSBp2  / number) - pow(total_SSBp  / number,2))/(number -1));
0261 SD_SSB2p = sqrt(((total_SSB2p2 / number) - pow(total_SSB2p / number,2))/(number -1));
0262
0263 SD_DSB   = sqrt(((total_DSB2   / number) - pow(total_DSB   / number,2))/(number -1));
0264 SD_DSBp  = sqrt(((total_DSBp2  / number) - pow(total_DSBp  / number,2))/(number -1));
0265 SD_DSBpp = sqrt(((total_DSBpp2 / number) - pow(total_DSBpp / number,2))/(number -1));
0266
0267 // Read damage classification SSBd, SSBi, SSBm, DSBd, DSBi, DSBm, DSBh
0268 // As they have been defined in: Nikjoo, H., O’Neill, O., Goodhead, T., & Terrissol, M. 1997,
0269 // Computational modelling of low-energy electron-induced DNA damage by early physical
0270 // and chemical events, International Journal of Radiation Biology, 71, 467.
0271 tree = (TTree *) f->Get("tuples/source");
0272 tree->SetBranchAddress("Primary",primaryName);
0273 tree->SetBranchAddress("Energy", &Energy);
0274 tree->SetBranchAddress("SSBd",   &SSBd);
0275 tree->SetBranchAddress("SSBi",   &SSBi);
0276 tree->SetBranchAddress("SSBm",   &SSBm);
0277 tree->SetBranchAddress("DSBd",   &DSBd);
0278 tree->SetBranchAddress("DSBi",   &DSBi);
0279 tree->SetBranchAddress("DSBm",   &DSBm);
0280 tree->SetBranchAddress("DSBh",   &DSBh);
0281
0282 Long64_t nentriesS = tree->GetEntries();
0283 for(int i = 0;i<nentriesS;i++){
0284   tree->GetEntry(i);
0285
0286   total_SSBd += SSBd;
0287   total_SSBd2 += pow((SSBd),2);
0288   total_SSBi += SSBi;
0289   total_SSBi2 += pow((SSBi),2);
0290   total_SSBm += SSBm;
0291   total_SSBm2 += pow((SSBm),2);
0292   total_sSSB += SSBd + SSBi + SSBm;
0293   total_sSSB2 += pow((SSBd+SSBi+SSBm),2);
0294
0295   total_DSBd  += DSBd;
0296   total_DSBd2 += pow(DSBd,2);
0297   total_DSBi  += DSBi;
0298   total_DSBi2 += pow(DSBi,2);
0299   total_DSBm  += DSBm;
0300   total_DSBm2 += pow(DSBm,2);
0301   total_DSBh  += DSBh;
0302   total_DSBh2 += pow(DSBh,2);
0303   total_sDSB  += DSBd + DSBi + DSBm + DSBh;
0304   total_sDSB2 += pow((DSBd+DSBi+DSBm+DSBh),2);
0305
0306   }
0307
0308 // Calculate the standard deviation
0309 SD_sSSB = sqrt(((total_sSSB2 / number) - pow(total_sSSB / number,2))/(number -1));
0310 SD_SSBd = sqrt(((total_SSBd2 / number) - pow(total_SSBd / number,2))/(number -1));
0311 SD_SSBi = sqrt(((total_SSBi2 / number) - pow(total_SSBi / number,2))/(number -1));
0312 SD_SSBm = sqrt(((total_SSBm2 / number) - pow(total_SSBm / number,2))/(number -1));
0313
0314 SD_sDSB = sqrt(((total_sDSB2 / number) - pow(total_sDSB / number,2))/(number -1));
0315 SD_DSBd = sqrt(((total_DSBd2 / number) - pow(total_DSBd / number,2))/(number -1));
0316 SD_DSBi = sqrt(((total_DSBi2 / number) - pow(total_DSBi / number,2))/(number -1));
0317 SD_DSBm = sqrt(((total_DSBm2 / number) - pow(total_DSBm / number,2))/(number -1));
0318 SD_DSBh = sqrt(((total_DSBh2 / number) - pow(total_DSBh / number,2))/(number -1));
0319
0320
0321 // Measure the Deposited Energy in the whole volume that includes DNA chain
0322
0323 tree = (TTree *) f->Get("tuples/chromosome_hits");
0324 tree->SetBranchAddress("e_chromosome_kev",&EnergyDeposited_eV);
0325 nentries = tree->GetEntries();
0326 for(int i = 0;i<nentries;i++){
0327   tree->GetEntry(i);
0328   acc_edep += EnergyDeposited_eV *1e3;
0329   acc_edep2 += EnergyDeposited_eV *EnergyDeposited_eV *1e6;
0330 }
0331 tree->SetBranchAddress("e_dna_kev",&EnergyDeposited_eV);
0332 nentries = tree->GetEntries();
0333 for(int i = 0;i<nentries;i++){
0334   tree->GetEntry(i);
0335   acc_edep += EnergyDeposited_eV *1e3;
0336   acc_edep2 += EnergyDeposited_eV *EnergyDeposited_eV *1e6;
0337 }
0338
0339 // Close the root file to free space
0340 f->Close();
0341
0342 // Calculate the absorbed dose
0343 dose = acc_edep * eVtoJ / mass;
0344 cout << acc_edep << "\n";
0345
0346 // This is a normalization factor to produce the output in Gy-1 Gbp-1, or else.
0347 // Default value is 1 to produce the result in Gy-1 Mbp-1
0348 // It changes Mbp to Gbp. Some other changes may be needed in graphs section (name of axes)
0349 double norm = 1000;
0350
0351 // Calculate the yields, together with their standard deviation
0352 EB_yield  = (Double_t) total_EB  / dose / Nbp;
0353 ES_yield  = (Double_t) total_ES  / dose / Nbp;
0354 OHB_yield = (Double_t) total_OHB / dose / Nbp;
0355 OHS_yield = (Double_t) total_OHS / dose / Nbp;
0356 HB_yield  = (Double_t) total_HB  / dose / Nbp;
0357 HS_yield  = (Double_t) total_HS  / dose / Nbp;
0358
0359 SD_EB_yield  = SD_EB  / dose / Nbp;
0360 SD_ES_yield  = SD_ES  / dose / Nbp;
0361 SD_OHB_yield = SD_OHB / dose / Nbp;
0362 SD_OHS_yield = SD_OHS / dose / Nbp;
0363 SD_HB_yield  = SD_HB  / dose / Nbp;
0364 SD_HS_yield  = SD_HS  / dose / Nbp;
0365
0366
0367 SSB_yield   = (Double_t) norm * total_SSB   / dose / Nbp;
0368 SSBp_yield  = (Double_t) norm * total_SSBp  / dose / Nbp;
0369 SSB2p_yield = (Double_t) norm * total_SSB2p / dose / Nbp;
0370
0371 DSB_yield   = (Double_t) norm * total_DSB   / dose / Nbp;
0372 DSBp_yield  = (Double_t) norm * total_DSBp  / dose / Nbp;
0373 DSBpp_yield = (Double_t) norm * total_DSBpp / dose / Nbp;
0374
0375 SD_SSB_yield   = norm * SD_SSB   / dose / Nbp;
0376 SD_SSBp_yield  = norm * SD_SSBp  / dose / Nbp;
0377 SD_SSB2p_yield = norm * SD_SSB2p / dose / Nbp;
0378
0379 SD_DSB_yield   = norm * SD_DSB   / dose / Nbp;
0380 SD_DSBp_yield  = norm * SD_DSBp  / dose / Nbp;
0381 SD_DSBpp_yield = norm * SD_DSBpp / dose / Nbp;
0382
0383
0384 sSSB_yield = (Double_t) norm * total_sSSB / dose / Nbp;
0385 SSBi_yield = (Double_t) norm * total_SSBi / dose / Nbp;
0386 SSBd_yield = (Double_t) norm * total_SSBd / dose / Nbp;
0387 SSBm_yield = (Double_t) norm * total_SSBm / dose / Nbp;
0388
0389 sDSB_yield = (Double_t) norm * total_sDSB / dose / Nbp;
0390 DSBi_yield = (Double_t) norm * total_DSBi / dose / Nbp;
0391 DSBd_yield = (Double_t) norm * total_DSBd / dose / Nbp;
0392 DSBm_yield = (Double_t) norm * total_DSBm / dose / Nbp;
0393 DSBh_yield = (Double_t) norm * total_DSBh / dose / Nbp;
0394
0395 SD_sSSB_yield = norm * SD_sSSB / dose / Nbp;
0396 SD_SSBi_yield = norm * SD_SSBi / dose / Nbp;
0397 SD_SSBd_yield = norm * SD_SSBd / dose / Nbp;
0398 SD_SSBm_yield = norm * SD_SSBm / dose / Nbp;
0399
0400 SD_sDSB_yield = norm * SD_sDSB / dose / Nbp;
0401 SD_DSBi_yield = norm * SD_DSBi / dose / Nbp;
0402 SD_DSBd_yield = norm * SD_DSBd / dose / Nbp;
0403 SD_DSBm_yield = norm * SD_DSBm / dose / Nbp;
0404 SD_DSBh_yield = norm * SD_DSBh / dose / Nbp;
0405
0406
0407 // Print output in terminal
0408
0409 float total_SSB_totalYield = SSB_yield + SSBp_yield + SSB2p_yield;
0410 float total_DSB_totalYield = DSB_yield + DSBp_yield + DSBpp_yield;
0411
0412 cout<<"\n"                      <<ifile         <<'\n'
0413     <<"\nDose Absorbed (Gy): "  <<dose          <<'\n'
0414     <<"Particle : "             <<primaryName   <<'\t'
0415     <<"Energy (MeV) : "         <<Energy        <<'\t'
0416     <<"Number of Primaries : "  <<number        <<'\n'
0417     <<"  Output Damage : "                      <<'\n'<<'\t'
0418     <<"     Species Hits (Gy-1 Mbp-1) "      <<'\n'<<'\t'
0419     <<"EaqBaseHits   : "     <<EB_yield     <<"   \t"      <<" error %: "   <<100*SD_EB_yield/EB_yield        <<'\n'<<'\t'
0420     <<"EaqStrandHits : "     <<ES_yield     <<"   \t"      <<" error %: "   <<100*SD_ES_yield/ES_yield        <<'\n'<<'\t'
0421     <<"OHBaseHits    : "     <<OHB_yield    <<"   \t"      <<" error %: "   <<100*SD_OHB_yield/OHB_yield      <<'\n'<<'\t'
0422     <<"OHStrandHits  : "     <<OHS_yield    <<"   \t"      <<" error %: "   <<100*SD_OHS_yield/OHS_yield      <<'\n'<<'\t'
0423     <<"HBaseHits     : "     <<HB_yield     <<"   \t"      <<" error %: "   <<100*SD_HB_yield/HB_yield        <<'\n'<<'\t'
0424     <<"HStrandHits   : "     <<HS_yield     <<"   \t"      <<" error %: "   <<100*SD_HS_yield/HS_yield        <<'\n'<<'\n'<<'\t'
0425     <<"     Damage yield (Gy-1 Gbp-1) "      <<'\n'<<'\t'
0426     <<"SSB           : "     <<SSB_yield    <<"   \t"      <<" error %: "   <<100*SD_SSB_yield/SSB_yield      <<'\n'<<'\t'
0427     <<"SSB+          : "     <<SSBp_yield   <<"   \t"      <<" error %: "   <<100*SD_SSBp_yield/SSBp_yield    <<'\n'<<'\t'
0428     <<"2SSB          : "     <<SSB2p_yield  <<"   \t"      <<" error %: "   <<100*SD_SSB2p_yield/SSB2p_yield  <<'\n'<<'\t'
0429     <<"SSB total     : "     <<total_SSB_totalYield        <<'\n'<<'\t'
0430     <<"DSB           : "     <<DSB_yield    <<"   \t"      <<" error %: "   <<100*SD_DSB_yield/DSB_yield     <<'\n'<<'\t'
0431     <<"DSB+          : "     <<DSBp_yield   <<"   \t"      <<" error %: "   <<100*SD_DSBp_yield/DSBp_yield   <<'\n'<<'\t'
0432     <<"DSB++         : "     <<DSBpp_yield  <<"   \t"      <<" error %: "   <<100*SD_DSBpp_yield/DSBpp_yield <<'\n'<<'\t'
0433     <<"DSB total     : "     <<total_DSB_totalYield        <<'\n'<<'\n'<<'\t'
0434     <<"     Breaks yield (Gy-1 Gbp-1) "      <<'\n'<<'\t'
0435     <<"SSB direct    : "     <<SSBd_yield   <<"   \t"      <<" error %: "   <<100*SD_SSBd_yield/SSBd_yield   <<'\n'<<'\t'
0436     <<"SSB indirect  : "     <<SSBi_yield   <<"   \t"      <<" error %: "   <<100*SD_SSBi_yield/SSBi_yield   <<'\n'<<'\t'
0437     <<"SSB mixed     : "     <<SSBm_yield   <<"   \t"      <<" error %: "   <<100*SD_SSBm_yield/SSBi_yield   <<'\n'<<'\t'
0438     <<"SSB total     : "     <<sSSB_yield   <<"   \t"      <<" error %: "   <<100*SD_sSSB_yield/sSSB_yield   <<'\n'<<'\t'
0439     <<"DSB direct    : "     <<DSBd_yield   <<"   \t"      <<" error %: "   <<100*SD_DSBd_yield/DSBd_yield   <<'\n'<<'\t'
0440     <<"DSB indirect  : "     <<DSBi_yield   <<"   \t"      <<" error %: "   <<100*SD_DSBi_yield/DSBi_yield   <<'\n'<<'\t'
0441     <<"DSB mixed     : "     <<DSBm_yield   <<"   \t"      <<" error %: "   <<100*SD_DSBm_yield/DSBm_yield   <<'\n'<<'\t'
0442     <<"DSB hybrid    : "     <<DSBh_yield   <<"   \t"      <<" error %: "   <<100*SD_DSBh_yield/DSBh_yield   <<'\n'<<'\t'
0443     <<"DSB total     : "     <<sDSB_yield   <<"   \t"      <<" error %: "   <<100*SD_sDSB_yield/sDSB_yield   <<'\n'<<'\n'<<'\t'
0444     <<"SSB/DSB       : "     <<sSSB_yield/sDSB_yield       <<'\n'<<'\n';
0445
0446
0447 // Plot Histograms
0448
0449 cfragment->GetCanvas()->cd();
0450 h1fragments->SetStats(false);
0451 h1fragments->SetMarkerSize(0.1);
0452 h1fragments->SetMarkerColor(kRed);
0453 h1fragments->SetLineColor  (kRed);
0454 h1fragments->Scale(1./(Nbp*1e6)); //bp^-1
0455 h1fragments->SetTitle("");
0456 h1fragments->SetYTitle("Number of Fragments (bp^{-2})");
0457 h1fragments->SetXTitle("Fragment Length (kbp)");
0458 h1fragments->SetAxisRange(10,1e4);
0459 h1fragments->SetMaximum(3e-11);
0460 h1fragments->SetMinimum(1e-15);
0461 h1fragments->Draw();
0462
0463
0464 c1->GetCanvas()->cd();
0465 pad1->cd();
0466 const Int_t n = 6;
0467 Double_t x[n] = {1,2,3,4,5,6};
0468 Double_t y[n] = {EB_yield,ES_yield,OHB_yield,OHS_yield,HB_yield,HS_yield};
0469 Double_t err_y[n] = {SD_EB_yield,SD_ES_yield,SD_OHB_yield,SD_OHS_yield,SD_HB_yield,SD_HS_yield};
0470 TGraph* gr = new TGraphErrors(n,x,y,0,err_y);
0471 gr->SetTitle("Species");
0472 gr->GetXaxis()->SetBinLabel(9, "EaqBaseHits");
0473 gr->GetXaxis()->SetBinLabel(25,"EaqStrandHits");
0474 gr->GetXaxis()->SetBinLabel(42,"OHBaseHits");
0475 gr->GetXaxis()->SetBinLabel(58,"OHStrandHits");
0476 gr->GetXaxis()->SetBinLabel(75,"HBaseHits");
0477 gr->GetXaxis()->SetBinLabel(92,"HStrandHits");
0478 gr->GetYaxis()->SetTitle("Species Hits (Gy^{-1} Mbp^{-1})");
0479 gr->GetYaxis()->SetTitleOffset(2);
0480
0481 gr->SetFillColor(49);
0482 gr->Draw("ba");
0483
0484
0485 pad2->cd();
0486 Double_t x2[n] = {1,2,3,4,5,6};
0487 Double_t y2[n] = {SSBp_yield,SSB2p_yield,SSB_yield,DSBp_yield,DSBpp_yield,DSB_yield};
0488 Double_t err_y2[n] = {SD_SSBp_yield,SD_SSB2p_yield,SD_SSB_yield,SD_DSBp_yield,SD_DSBpp_yield,SD_DSB_yield};
0489 TGraph* gr2 = new TGraphErrors(n,x2,y2,0,err_y2);
0490 gr2->SetTitle("Damage Yield");
0491 gr2->GetXaxis()->SetBinLabel(9, "SSB+");
0492 gr2->GetXaxis()->SetBinLabel(25,"2SSB");
0493 gr2->GetXaxis()->SetBinLabel(42,"SSB");
0494 gr2->GetXaxis()->SetBinLabel(58,"DSB+");
0495 gr2->GetXaxis()->SetBinLabel(75,"DSB++");
0496 gr2->GetXaxis()->SetBinLabel(92,"DSB");
0497 //gr2->GetYaxis()->SetTitle("Damage yield (Gy^{-1} Mbp^{-1})");
0498 gr2->GetYaxis()->SetTitle("Damage yield (Gy^{-1} Gbp^{-1})");
0499 gr2->GetYaxis()->SetTitleOffset(2);
0500
0501 gr2->SetFillColor(8);
0502 gr2->Draw("ba");
0503
0504
0505 pad3->cd();
0506 const Int_t m = 4;
0507 Double_t x3[m] = {1,2,3,4};
0508 Double_t y3[m] = {SSBd_yield,SSBi_yield,SSBm_yield,sSSB_yield};
0509 Double_t err_y3[m] = {SD_SSBd_yield,SD_SSBi_yield,SD_SSBm_yield,SD_sSSB_yield};
0510 TGraph* gr3 = new TGraphErrors(m,x3,y3,0,err_y3);
0511 gr3->SetTitle("Breaks Yield");
0512 gr3->GetXaxis()->SetBinLabel(8, "SSB direct");
0513 gr3->GetXaxis()->SetBinLabel(35,"SSB indirect");
0514 gr3->GetXaxis()->SetBinLabel(64,"SSB mixed");
0515 gr3->GetXaxis()->SetBinLabel(92,"SSB all");
0516 //gr3->GetYaxis()->SetTitle("Breaks yield (Gy^{-1} Mbp^{-1})");
0517 gr3->GetYaxis()->SetTitle("SSB yield (Gy^{-1} Gbp^{-1})");
0518 gr3->GetYaxis()->SetTitleOffset(2);
0519
0520 gr3->SetFillColor(7);
0521 gr3->Draw("ba");
0522
0523
0524 pad4->cd();
0525 const Int_t k = 5;
0526 Double_t x4[k] = {1,2,3,4,5};
0527 Double_t y4[k] = {DSBd_yield,DSBi_yield,DSBm_yield,DSBh_yield,sDSB_yield};
0528 Double_t err_y4[k] = {SD_DSBd_yield,SD_DSBi_yield,SD_DSBm_yield,SD_DSBh_yield,SD_sDSB_yield};
0529 TGraph* gr4 = new TGraphErrors(k,x4,y4,0,err_y4);
0530 gr4->SetTitle("Breaks Yield");
0531 gr4->GetXaxis()->SetBinLabel(8,"DSB direct");
0532 gr4->GetXaxis()->SetBinLabel(29,"DSB indirect");
0533 gr4->GetXaxis()->SetBinLabel(50,"DSB mixed");
0534 gr4->GetXaxis()->SetBinLabel(71,"DSB hybrid");
0535 gr4->GetXaxis()->SetBinLabel(92,"DSB all");
0536 //gr4->GetYaxis()->SetTitle("Breaks yield (Gy^{-1} Mbp^{-1})");
0537 gr4->GetYaxis()->SetTitle("DSB yield (Gy^{-1} Gbp^{-1})");
0538 gr4->GetYaxis()->SetTitleOffset(2);
0539
0540 gr4->SetFillColor(4);
0541 gr4->Draw("ba");
0542
0543 }
0544
0545 // Some important bools that are needed to run the root macro file
0546 bool greaterPair(const ipair& l, const ipair& r){return l.second > r.second;}
0547 bool smallerPair(const ipair& l, const ipair& r){return l.second < r.second;}
0548
0549 void BinLogX(TH1 *h) {
0550   TAxis *axis = h->GetXaxis();
0551   int bins = axis->GetNbins();
0552   Axis_t from = axis->GetXmin();
0553   Axis_t to = axis->GetXmax();
0554   Axis_t width = (to - from) / bins;
0555   Axis_t *new_bins = new Axis_t[bins + 1];
0556   for (int i = 0; i <= bins; i++) {
0557     new_bins[i] = TMath::Power(10, from + i * width);
0558   }
0559   axis->Set(bins, new_bins);
0560   delete[] new_bins;
0561
0562 }