//------------------------------------------------------------------------------- // Modified by Sara Zein to calculate the damage probability per plasmid //------------------------------------------------------------------------------- // // This macro requires the molecular-dna.root file generated from molecularDNA example // To run this file just insert this command to the terminal: // root .X plasmid.C // //***************************************// // Please define the parameters below // // ifile, r3, Nbp (as shown in terminal) // //***************************************// { //*******************************************************************************// // Define these parameters of the simulation char ifile[256] = "molecular-dna.root"; // input filepath const Int_t numberOfPlasmids = 10144; Double_t r3 = 4.42e-6 * 4.42e-6 * 4.84e-6; // r * r * r (world cube side m) Double_t Nbp = 4.367 * 0.001 * numberOfPlasmids; // Mbp // Length of the DNA chain in Mbp // multiplying by 10142 since we are testing 10142 k plasmids Double_t mass = 997 * r3; // density * r3 //*******************************************************************************// typedef std::pair ipair; bool greaterPair(const ipair& l, const ipair& r); bool smallerPair(const ipair& l, const ipair& r); void BinLogX(TH1 * h); gROOT->Reset(); gStyle->SetPalette(1); gROOT->SetStyle("Plain"); gStyle->SetOptStat(00000); gStyle->SetPadTickX(1); gStyle->SetPadTickY(1); // Initialize output histograms TCanvas* cdamage = new TCanvas("cdamage", "DNA Damage Distribution", 900, 120, 600, 400); cdamage->SetLogy(); TH1F* h1damage = new TH1F("h1damage", "h1damage", 28, 0, 14); TCanvas* c4damage = new TCanvas("c4damage", "Direct damage per plasmid", 900, 120, 600, 400); // c4damage->SetLogy(); TH1F* h4damage = new TH1F("h4damage", "h4damage", 28, 0, 14); TCanvas* cSSB = new TCanvas("cSSB", "SSB Distribution", 900, 120, 600, 400); cSSB->SetLogy(); TH1F* h1SSB = new TH1F("h1SSB", "h1SSB", 28, 0, 14); TCanvas* cDSB = new TCanvas("cDSB", "DSB Distribution", 900, 120, 600, 400); cDSB->SetLogy(); TH1F* h1DSB = new TH1F("h1DSB", "h1DSB", 40, 0, 10); TCanvas* ccount = new TCanvas("cCount", "Damage per plasmid Distribution", 900, 120, 600, 400); TH1F* h1count = new TH1F("h1count", "h1count", 11000, 0, 11000); TGraph2D* gr1 = new TGraph2D(); // Open root file TFile* f = TFile::Open(ifile); // Initialize Variables Int_t EB, ES, OHB, OHS, HB, HS, FL; Int_t total_EB, total_ES, total_OHB, total_OHS, total_HB, total_HS, total_FL; Float_t total_EB2, total_ES2, total_OHB2, total_OHS2, total_HB2, total_HS2, total_FL2; Float_t SD_EB, SD_ES, SD_OHB, SD_OHS, SD_HB, SD_HS; Float_t SD_SSB, SD_SSBp, SD_SSB2p, SD_sSSB, SD_SSBd, SD_SSBi, SD_SSBm; Float_t SD_DSB, SD_DSBp, SD_DSBpp, SD_sDSB, SD_DSBd, SD_DSBi, SD_DSBm, SD_DSBh; Int_t SSB, SSBp, SSB2p; Int_t total_SSB, total_SSBp, total_SSB2p; Float_t total_SSB2, total_SSBp2, total_SSB2p2; Int_t DSB, DSBp, DSBpp; Int_t total_DSB, total_DSBp, total_DSBpp; Float_t total_DSB2, total_DSBp2, total_DSBpp2; Int_t SSBd, SSBi, SSBm; Int_t total_sSSB, total_SSBd, total_SSBi, total_SSBm; Float_t total_sSSB2, total_SSBd2, total_SSBi2, total_SSBm2; Int_t DSBd, DSBi, DSBm, DSBh; Int_t total_sDSB, total_DSBd, total_DSBi, total_DSBm, total_DSBh; Float_t total_sDSB2, total_DSBd2, total_DSBi2, total_DSBm2, total_DSBh2; Double_t dose = 0; Double_t SD_dose = 0; 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; 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; Double_t SSB_yield = 0; Double_t SSBp_yield = 0; Double_t SSB2p_yield = 0; Double_t SD_SSB_yield = 0; Double_t SD_SSBp_yield = 0; Double_t SD_SSB2p_yield = 0; Double_t DSB_yield = 0; Double_t DSBp_yield = 0; Double_t DSBpp_yield = 0; Double_t SD_DSB_yield = 0; Double_t SD_DSBp_yield = 0; Double_t SD_DSBpp_yield = 0; Double_t sSSB_yield = 0; Double_t SSBi_yield = 0; Double_t SSBd_yield = 0; Double_t SSBm_yield = 0; Double_t SD_sSSB_yield = 0; Double_t SD_SSBi_yield = 0; Double_t SD_SSBd_yield = 0; Double_t SD_SSBm_yield = 0; 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; 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; total_EB = 0; total_ES = 0; total_OHB = 0; total_OHS = 0; total_HB = 0; total_HS = 0; total_SSB = 0; total_SSBp = 0; total_SSB2p = 0; total_SSB2 = 0; total_SSBp2 = 0; total_SSB2p2 = 0; total_DSB = 0; total_DSBp = 0; total_DSBpp = 0; total_DSB2 = 0; total_DSBp2 = 0; total_DSBpp2 = 0; total_sSSB = 0; total_SSBd = 0; total_SSBi = 0; total_SSBm = 0; total_sSSB2 = 0; total_SSBd2 = 0; total_SSBi2 = 0; total_SSBm2 = 0; total_sDSB = 0; total_DSBd = 0; total_DSBi = 0; total_DSBm = 0; total_DSBh = 0; total_sDSB2 = 0; total_DSBd2 = 0; total_DSBi2 = 0; total_DSBm2 = 0; total_DSBh2 = 0; Double_t eVtoJ = 1.60218e-19; Double_t EnergyDeposited_eV = 0; Double_t acc_edep = 0; Double_t acc_edep2 = 0; Double_t Energy; Double_t BPID; Int_t Strand; Int_t StrandDamage; Int_t Event1, directB; Char_t Primary; char* primaryName = new char[32]; char* type = new char[256]; char* sourceClassification = new char[256]; Double_t xx, yy, zz; // Read trees and leaves from root file, and give values to variables TTree* tree = (TTree*)f->Get("tuples/primary_source"); Float_t number = (Float_t)tree->GetEntries(); if (number<2) { std::cout << "Not enough entries in the \"primary_source\" TTree (" << (long)number << " entries)\n"; gApplication->Terminate(0); } vector> DSBBPID; // For reading species production tree = (TTree*)f->Get("tuples/damage"); tree->SetBranchAddress("Primary", &Primary); tree->SetBranchAddress("Energy", &Energy); tree->SetBranchAddress("EaqBaseHits", &EB); tree->SetBranchAddress("EaqStrandHits", &ES); tree->SetBranchAddress("OHBaseHits", &OHB); tree->SetBranchAddress("OHStrandHits", &OHS); tree->SetBranchAddress("HBaseHits", &HB); tree->SetBranchAddress("HStrandHits", &HS); tree->SetBranchAddress("TypeClassification", type); tree->SetBranchAddress("BasePair", &BPID); tree->SetBranchAddress("Event", &Event1); tree->SetBranchAddress("Strand", &Strand); tree->SetBranchAddress("StrandDamage", &StrandDamage); tree->SetBranchAddress("Position_x_um", &xx); tree->SetBranchAddress("Position_y_um", &yy); tree->SetBranchAddress("Position_z_um", &zz); tree->SetBranchAddress("DirectBreaks", &directB); int primea = 0; Long64_t nentries = tree->GetEntries(); int damagecount = 0; int hitcount = 0; // Double_t plasmidsD [25992]; Double_t plasmidsD[2599200]; for (int i = 0; i < nentries; i++) { tree->GetEntry(i); total_EB += EB; total_EB2 += pow(EB, 2); total_ES += ES; total_ES2 += pow(ES, 2); total_OHB += OHB; total_OHB2 += pow(OHB, 2); total_OHS += OHS; total_OHS2 += pow(OHS, 2); total_HB += HB; total_HB2 += pow(HB, 2); total_HS += HS; total_HS2 += pow(HS, 2); if ((string)type == "DSB" || (string)type == "DSB+" || (string)type == "DSB++") { DSBBPID.push_back(make_pair(i, (int64_t)BPID)); } if (StrandDamage != 0) { { damagecount++; h1count->Fill(Strand); } primea++; } } int unbrokenP = 0; int brokenP = 0; for (int i = 0; i < 11000; i++) { plasmidsD[i] = h1count->GetBinContent(i); if (plasmidsD[i] == 0 && i < numberOfPlasmids) unbrokenP++; } Double_t totalDs = 0.0; for (int i = 0; i < 11000; i++) { if (plasmidsD[i] != 0) { h4damage->Fill(plasmidsD[i]); totalDs += plasmidsD[i]; brokenP++; } } Double_t X; for (int i = 0; i < 28; i++) { X = h4damage->GetBinContent(i); h4damage->SetBinContent(i, X * 100 / totalDs); } // Find the number of fragments that have been produced, but first test if there are enough breaks. // If no more than 2 DSBs exist in the DSBBPID vector ( DSBBPID.size() is 0 or 1 ), // the subtraction DSBBPID.size() - 1 makes the loop condition evaluate to ie < -1 (in unsigned terms, // this becomes a large number, which is incorrect) and leads to undefined behavior (crash). if (DSBBPID.size() < 2) { std::cerr << "Not enough damage to calculate fragments." << std::endl; //return; } if (DSBBPID.size() >= 2) { // Sort DSBs from the one with lower ID value to the one with higher ID value // Then find the number of fragments that have been produced sort(DSBBPID.begin(), DSBBPID.end(), smallerPair); for(int ie = 0; ie < DSBBPID.size() - 1; ie++){ int64_t dsbfragment = DSBBPID[ie + 1].second - DSBBPID[ie].second; } } // Calculate the SEM SD_EB = sqrt(abs(((total_EB2 / nentries) - pow(total_EB / nentries, 2))) / (nentries - 1)); SD_ES = sqrt(abs(((total_ES2 / nentries) - pow(total_ES / nentries, 2))) / (nentries - 1)); SD_OHB = sqrt(abs(((total_OHB2 / nentries) - pow(total_OHB / nentries, 2))) / (nentries - 1)); SD_OHS = sqrt(abs(((total_OHS2 / nentries) - pow(total_OHS / nentries, 2))) / (nentries - 1)); SD_HB = sqrt(abs(((total_HB2 / nentries) - pow(total_HB / nentries, 2))) / (nentries - 1)); SD_HS = sqrt(abs(((total_HS2 / nentries) - pow(total_HS / nentries, 2))) / (nentries - 1)); // Read damage classification SSB, SSB+, 2SSB, DSB, DSB+, DSB++ // As they have been defined in: Nikjoo, H., O’Neill, O., Goodhead, T., & Terrissol, M. 1997, // Computational modelling of low-energy electron-induced DNA damage by early physical // and chemical events, International Journal of Radiation Biology, 71, 467. tree = (TTree*)f->Get("tuples/classification"); tree->SetBranchAddress("Primary", &Primary); tree->SetBranchAddress("Energy", &Energy); tree->SetBranchAddress("SSB", &SSB); tree->SetBranchAddress("SSBp", &SSBp); tree->SetBranchAddress("2SSB", &SSB2p); tree->SetBranchAddress("DSB", &DSB); tree->SetBranchAddress("DSBp", &DSBp); tree->SetBranchAddress("DSBpp", &DSBpp); Long64_t nentriesC = tree->GetEntries(); for (int i = 0; i < nentriesC; i++) { tree->GetEntry(i); total_SSBp += SSBp; total_SSBp2 += pow(SSBp, 2); total_SSB2p += SSB2p; total_SSB2p2 += pow(SSB2p, 2); total_SSB += SSB; total_SSB2 += pow(SSB, 2); total_DSBp += DSBp; total_DSBp2 += pow(DSBp, 2); total_DSBpp += DSBpp; total_DSBpp2 += pow(DSBpp, 2); total_DSB += DSB; total_DSB2 += pow(DSB, 2); } // Calculate the SEM SD_SSB = sqrt(abs(((total_SSB2 / nentriesC) - pow(total_SSB / nentriesC, 2))) / (nentriesC - 1)); SD_SSBp = sqrt(abs(((total_SSBp2 / nentriesC) - pow(total_SSBp / nentriesC, 2))) / (nentriesC - 1)); SD_SSB2p = sqrt(abs(((total_SSB2p2 / nentriesC) - pow(total_SSB2p / nentriesC, 2))) / (nentriesC - 1)); SD_DSB = sqrt(abs(((total_DSB2 / nentriesC) - pow(total_DSB / nentriesC, 2))) / (nentriesC - 1)); SD_DSBp = sqrt(abs(((total_DSBp2 / nentriesC) - pow(total_DSBp / nentriesC, 2))) / (nentriesC - 1)); SD_DSBpp = sqrt(abs(((total_DSBpp2 / nentriesC) - pow(total_DSBpp / nentriesC, 2))) / (nentriesC - 1)); // Read damage classification SSBd, SSBi, SSBm, DSBd, DSBi, DSBm, DSBh // As they have been defined in: Nikjoo, H., O’Neill, O., Goodhead, T., & Terrissol, M. 1997, // Computational modelling of low-energy electron-induced DNA damage by early physical // and chemical events, International Journal of Radiation Biology, 71, 467. tree = (TTree*)f->Get("tuples/source"); tree->SetBranchAddress("Primary", primaryName); tree->SetBranchAddress("Energy", &Energy); tree->SetBranchAddress("SSBd", &SSBd); tree->SetBranchAddress("SSBi", &SSBi); tree->SetBranchAddress("SSBm", &SSBm); tree->SetBranchAddress("DSBd", &DSBd); tree->SetBranchAddress("DSBi", &DSBi); tree->SetBranchAddress("DSBm", &DSBm); tree->SetBranchAddress("DSBh", &DSBh); int iprime = 0; Long64_t nentriesS = tree->GetEntries(); for (int i = 0; i < nentriesS; i++) { tree->GetEntry(i); total_SSBd += SSBd; total_SSBd2 += pow((SSBd), 2); total_SSBi += SSBi; total_SSBi2 += pow((SSBi), 2); total_SSBm += SSBm; total_SSBm2 += pow((SSBm), 2); total_sSSB += SSBd + SSBi + SSBm; total_sSSB2 += pow((SSBd + SSBi + SSBm), 2); total_DSBd += DSBd; total_DSBd2 += pow(DSBd, 2); total_DSBi += DSBi; total_DSBi2 += pow(DSBi, 2); total_DSBm += DSBm; total_DSBm2 += pow(DSBm, 2); total_DSBh += DSBh; total_DSBh2 += pow(DSBh, 2); total_sDSB += DSBd + DSBi + DSBm + DSBh; total_sDSB2 += pow((DSBd + DSBi + DSBm + DSBh), 2); if (SSBd != 0) h1SSB->Fill(SSBd); if (DSBd != 0) h1DSB->Fill(DSBd); if (SSBd != 0) h1damage->Fill(SSBd); if (DSBd != 0) h1damage->Fill(DSBd); if (SSB != 0 || DSBd != 0) { iprime++; } } Double_t Y; for (int i = 0; i < 28; i++) { Y = h1damage->GetBinContent(i); h1damage->SetBinContent(i, Y * 100 / totalDs); } // Calculate the SEM SD_sSSB = sqrt(abs(((total_sSSB2 / nentriesS) - pow(total_sSSB / nentriesS, 2))) / (nentriesS - 1)); SD_SSBd = sqrt(abs(((total_SSBd2 / nentriesS) - pow(total_SSBd / nentriesS, 2))) / (nentriesS - 1)); SD_SSBi = sqrt(abs(((total_SSBi2 / nentriesS) - pow(total_SSBi / nentriesS, 2))) / (nentriesS - 1)); SD_SSBm = sqrt(abs(((total_SSBm2 / nentriesS) - pow(total_SSBm / nentriesS, 2))) / (nentriesS - 1)); SD_sDSB = sqrt(abs(((total_sDSB2 / nentriesS) - pow(total_sDSB / nentriesS, 2))) / (nentriesS - 1)); SD_DSBd = sqrt(abs(((total_DSBd2 / nentriesS) - pow(total_DSBd / nentriesS, 2))) / (nentriesS - 1)); SD_DSBi = sqrt(abs(((total_DSBi2 / nentriesS) - pow(total_DSBi / nentriesS, 2))) / (nentriesS - 1)); SD_DSBm = sqrt(abs(((total_DSBm2 / nentriesS) - pow(total_DSBm / nentriesS, 2))) / (nentriesS - 1)); SD_DSBh = sqrt(abs(((total_DSBh2 / nentriesS) - pow(total_DSBh / nentriesS, 2))) / (nentriesS - 1)); // Measure the Deposited Energy in the whole volume that includes DNA chain Long64_t totnentries =0; tree = (TTree*)f->Get("tuples/chromosome_hits"); tree->SetBranchAddress("e_chromosome_kev", &EnergyDeposited_eV); nentries = tree->GetEntries(); totnentries+=nentries; for (int i = 0; i < nentries; i++) { tree->GetEntry(i); acc_edep += EnergyDeposited_eV * 1e3; acc_edep2 += EnergyDeposited_eV * EnergyDeposited_eV * 1e6; } // Close the root file to free space f->Close(); // Calculate the absorbed dose dose = acc_edep * eVtoJ / mass; SD_dose = sqrt(abs(((acc_edep2 / totnentries) - pow(acc_edep / totnentries,2)))/(totnentries -1))* eVtoJ / mass; double norm = 1; // Calculate the yields, together with their error EB_yield = (Double_t)total_EB / dose / Nbp; ES_yield = (Double_t)total_ES / dose / Nbp; OHB_yield = (Double_t)total_OHB / dose / Nbp; OHS_yield = (Double_t)total_OHS / dose / Nbp; HB_yield = (Double_t)total_HB / dose / Nbp; HS_yield = (Double_t)total_HS / dose / Nbp; SD_EB_yield = total_EB / dose / Nbp * sqrt(pow(SD_EB /total_EB ,2)+pow(SD_dose/dose,2)); SD_ES_yield = total_ES / dose / Nbp * sqrt(pow(SD_ES /total_ES ,2)+pow(SD_dose/dose,2)); SD_OHB_yield = total_OHB / dose / Nbp * sqrt(pow(SD_OHB /total_OHB ,2)+pow(SD_dose/dose,2)); SD_OHS_yield = total_OHS / dose / Nbp * sqrt(pow(SD_OHS /total_OHS ,2)+pow(SD_dose/dose,2)); SD_HB_yield = total_HB / dose / Nbp * sqrt(pow(SD_HB /total_HB ,2)+pow(SD_dose/dose,2)); SD_HS_yield = total_HS / dose / Nbp * sqrt(pow(SD_HS /total_HS ,2)+pow(SD_dose/dose,2)); //Previous method //SD_EB_yield = SD_EB / dose / Nbp; //SD_ES_yield = SD_ES / dose / Nbp; //SD_OHB_yield = SD_OHB / dose / Nbp; //SD_OHS_yield = SD_OHS / dose / Nbp; //SD_HB_yield = SD_HB / dose / Nbp; //SD_HS_yield = SD_HS / dose / Nbp; SSB_yield = (Double_t)norm * total_SSB / dose / Nbp; SSBp_yield = (Double_t)norm * total_SSBp / dose / Nbp; SSB2p_yield = (Double_t)norm * total_SSB2p / dose / Nbp; DSB_yield = (Double_t)norm * total_DSB / dose / Nbp; DSBp_yield = (Double_t)norm * total_DSBp / dose / Nbp; DSBpp_yield = (Double_t)norm * total_DSBpp / dose / Nbp; SD_SSB_yield = norm * total_SSB / dose / Nbp * sqrt(pow(SD_SSB /total_SSB ,2)+pow(SD_dose/dose,2)); SD_SSBp_yield = norm * total_SSBp / dose / Nbp * sqrt(pow(SD_SSBp /total_SSBp ,2)+pow(SD_dose/dose,2)); SD_SSB2p_yield = norm * total_SSB2p / dose / Nbp * sqrt(pow(SD_SSB2p /total_SSB2p ,2)+pow(SD_dose/dose,2)); SD_DSB_yield = norm * total_DSB / dose / Nbp * sqrt(pow(SD_DSB /total_DSB ,2)+pow(SD_dose/dose,2)); SD_DSBp_yield = norm * total_DSBp / dose / Nbp * sqrt(pow(SD_DSBp /total_DSBp ,2)+pow(SD_dose/dose,2)); SD_DSBpp_yield = norm * total_DSBpp / dose / Nbp * sqrt(pow(SD_DSBpp /total_DSBpp ,2)+pow(SD_dose/dose,2)); //Previous method //SD_SSB_yield = norm * SD_SSB / dose / Nbp; //SD_SSBp_yield = norm * SD_SSBp / dose / Nbp; //SD_SSB2p_yield = norm * SD_SSB2p / dose / Nbp; //SD_DSB_yield = norm * SD_DSB / dose / Nbp; //SD_DSBp_yield = norm * SD_DSBp / dose / Nbp; //SD_DSBpp_yield = norm * SD_DSBpp / dose / Nbp; sSSB_yield = (Double_t)norm * total_sSSB / dose / Nbp; SSBi_yield = (Double_t)norm * total_SSBi / dose / Nbp; SSBd_yield = (Double_t)norm * total_SSBd / dose / Nbp; SSBm_yield = (Double_t)norm * total_SSBm / dose / Nbp; sDSB_yield = (Double_t)norm * total_sDSB / dose / Nbp; DSBi_yield = (Double_t)norm * total_DSBi / dose / Nbp; DSBd_yield = (Double_t)norm * total_DSBd / dose / Nbp; DSBm_yield = (Double_t)norm * total_DSBm / dose / Nbp; DSBh_yield = (Double_t)norm * total_DSBh / dose / Nbp; SD_sSSB_yield = norm * total_sSSB/ dose / Nbp* sqrt(pow(SD_sSSB /total_sSSB,2)+pow(SD_dose/dose,2)); SD_SSBi_yield = norm * total_SSBi/ dose / Nbp* sqrt(pow(SD_SSBi /total_SSBi,2)+pow(SD_dose/dose,2)); SD_SSBd_yield = norm * total_SSBd/ dose / Nbp* sqrt(pow(SD_SSBd /total_SSBd,2)+pow(SD_dose/dose,2)); SD_SSBm_yield = norm * total_SSBm/ dose / Nbp* sqrt(pow(SD_SSBm /total_SSBm,2)+pow(SD_dose/dose,2)); SD_sDSB_yield = norm * total_sDSB/ dose / Nbp* sqrt(pow(SD_sDSB /total_sDSB,2)+pow(SD_dose/dose,2)); SD_DSBi_yield = norm * total_DSBi/ dose / Nbp* sqrt(pow(SD_DSBi /total_DSBi,2)+pow(SD_dose/dose,2)); SD_DSBd_yield = norm * total_DSBd/ dose / Nbp* sqrt(pow(SD_DSBd /total_DSBd,2)+pow(SD_dose/dose,2)); SD_DSBm_yield = norm * total_DSBm/ dose / Nbp* sqrt(pow(SD_DSBm /total_DSBm,2)+pow(SD_dose/dose,2)); SD_DSBh_yield = norm * total_DSBh/ dose / Nbp* sqrt(pow(SD_DSBh /total_DSBh,2)+pow(SD_dose/dose,2)); //Previous method //SD_sSSB_yield = norm * SD_sSSB / dose / Nbp; //SD_SSBi_yield = norm * SD_SSBi / dose / Nbp; //SD_SSBd_yield = norm * SD_SSBd / dose / Nbp; //SD_SSBm_yield = norm * SD_SSBm / dose / Nbp; //SD_sDSB_yield = norm * SD_sDSB / dose / Nbp; //SD_DSBi_yield = norm * SD_DSBi / dose / Nbp; //SD_DSBd_yield = norm * SD_DSBd / dose / Nbp; //SD_DSBm_yield = norm * SD_DSBm / dose / Nbp; //SD_DSBh_yield = norm * SD_DSBh / dose / Nbp; // Print output in terminal float total_SSB_totalYield = SSB_yield + SSBp_yield + SSB2p_yield; float total_DSB_totalYield = DSB_yield + DSBp_yield + DSBpp_yield; cout << "\n" << " Output file: " << ifile << '\n' << "\nDose Absorbed (Gy): " << dose << '\n' << "Particle : " << primaryName << '\t' << "Energy (MeV) : " << Energy << '\t' << "Number of Primaries : " << number << '\n' << " Output Damage : " << '\n' << '\t' << "Number of plasmids: " << numberOfPlasmids << " \t" << '\n' << '\t' << "Unbroken plasmids: " << unbrokenP << " \t" << '\n' << '\t' << "Broken plasmids: " << brokenP << " \t" << '\n' << '\t' << "Rate of damages per plasmid " << totalDs / numberOfPlasmids << " \t" << '\n' << '\t' << "Rate of damages per plasmid per dose " << '\t' << totalDs / numberOfPlasmids / dose << " plasmid-1Gy-1" << '\n' << '\t' << "Rate of damages per dose per Mbp " << totalDs / dose / Nbp << " Gy-1Mbp-1" << '\n' << '\t' << " Species Hits " << '\n' << '\t' << "EaqBaseHits " << EB_yield * dose * Nbp << " \t" << '\n' << '\t' << "EaqStrandHits " << ES_yield * dose * Nbp << " \t" << '\n' << '\t' << "OHBaseHits " << OHB_yield * dose * Nbp << " \t" << '\n' << '\t' << "OHStrandHits " << OHS_yield * dose * Nbp << " \t" << '\n' << '\t' << "HBaseHits " << HB_yield * dose * Nbp << " \t" << '\n' << '\t' << "HStrandHits " << HS_yield * dose * Nbp << " \t" << '\n' << '\n' << '\t' << " Damage number" << '\n' << '\t' << "SSB " << SSB_yield * dose * Nbp << " \t" << '\n' << '\t' << "SSB+ " << SSBp_yield * dose * Nbp << " \t" << '\n' << '\t' << "2SSB " << SSB2p_yield * dose * Nbp << " \t" << '\n' << '\t' << "SSB total " << total_SSB_totalYield * dose * Nbp << '\n' << '\t' << "DSB " << DSB_yield * dose * Nbp << " \t" << '\n' << '\t' << "DSB+ " << DSBp_yield * dose * Nbp << " \t" << '\n' << '\t' << "DSB++ " << DSBpp_yield * dose * Nbp << " \t" << '\n' << '\t' << "DSB total " << total_DSB_totalYield * dose * Nbp << '\n' << '\n' << '\t' << " Breaks number " << '\n' << '\t' << "SSB direct " << SSBd_yield * dose * Nbp << " \t" << '\n' << '\t' << "SSB indirect " << SSBi_yield * dose * Nbp << " \t" << '\n' << '\t' << "SSB mixed " << SSBm_yield * dose * Nbp << " \t" << '\n' << '\t' << "SSB total " << sSSB_yield * dose * Nbp << " \t" << '\n' << '\t' << "DSB direct " << DSBd_yield * dose * Nbp << " \t" << '\n' << '\t' << "DSB indirect " << DSBi_yield * dose * Nbp << " \t" << '\n' << '\t' << "DSB mixed " << DSBm_yield * dose * Nbp << " \t" << '\n' << '\t' << "DSB hybrid " << DSBh_yield * dose * Nbp << " \t" << '\n' << '\t' << "DSB total " << sDSB_yield * dose * Nbp << " \t" << '\n' << '\n' << '\t' << "SSB/DSB " << sSSB_yield / sDSB_yield << '\n' << '\n'; // Plot Histograms cSSB->GetCanvas()->cd(); h1SSB->SetStats(false); h1SSB->SetMarkerSize(0.1); h1SSB->SetMarkerColor(kBlue); h1SSB->SetLineColor(kBlue); h1SSB->SetTitle(""); h1SSB->SetYTitle(" "); h1SSB->SetXTitle("Number of direct SSB per event"); h1SSB->SetFillColor(kBlue); h1SSB->Draw(); cDSB->GetCanvas()->cd(); h1DSB->SetStats(false); h1DSB->SetMarkerSize(0.1); h1DSB->SetMarkerColor(kGreen + 2); h1DSB->SetLineColor(kGreen + 2); h1DSB->SetTitle(""); h1DSB->SetYTitle(" "); h1DSB->SetXTitle("Number of direct DSB per event"); h1DSB->SetFillColor(kGreen + 2); h1DSB->Draw(); cdamage->GetCanvas()->cd(); h1damage->SetStats(false); h1damage->SetMarkerSize(0.1); h1damage->SetMarkerColor(kMagenta); h1damage->SetLineColor(kMagenta); h1damage->SetTitle(""); h1damage->SetYTitle("Percentage % "); h1damage->SetXTitle("Number of damages per event"); h1damage->SetFillColor(kMagenta); h1damage->Draw(); c4damage->GetCanvas()->cd(); h4damage->SetStats(false); h4damage->SetMarkerSize(0.1); h4damage->SetMarkerColor(kRed - 4); h4damage->SetLineColor(kRed - 4); h4damage->SetTitle(""); h4damage->SetYTitle("Percentage % "); h4damage->SetXTitle("Number of damages per plasmid"); h4damage->SetFillColor(kRed - 4); h4damage->Draw(); ccount->GetCanvas()->cd(); h1count->SetStats(false); h1count->SetMarkerSize(0.1); h1count->SetMarkerColor(kCyan); h1count->SetLineColor(kCyan); h1count->SetTitle(""); h1count->SetYTitle("Number of damages"); h1count->SetXTitle("Plasmid ID"); h1count->Draw(); } // Some important bools that are needed to run the root macro file bool greaterPair(const ipair& l, const ipair& r) { return l.second > r.second; } bool smallerPair(const ipair& l, const ipair& r) { return l.second < r.second; } void BinLogX(TH1* h) { TAxis* axis = h->GetXaxis(); int bins = axis->GetNbins(); Axis_t from = axis->GetXmin(); Axis_t to = axis->GetXmax(); Axis_t width = (to - from) / bins; Axis_t* new_bins = new Axis_t[bins + 1]; for (int i = 0; i <= bins; i++) { new_bins[i] = TMath::Power(10, from + i * width); } axis->Set(bins, new_bins); delete[] new_bins; }