Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-10-13 08:25:11 

0001 //-------------------------------------------------------------------------------
0002 // Modified by Sara Zein to calculate the damage probability per plasmid
0003 //-------------------------------------------------------------------------------
0004 //
0005 // This macro requires the molecular-dna.root file generated from molecularDNA example
0006 // To run this file just insert this command to the terminal:
0007 // root .X plasmid.C
0008 //
0009 //***************************************//
0010 // Please define the parameters below    //
0011 // ifile, r3, Nbp (as shown in terminal) //
0012 //***************************************//
0013 
0014 {
0015   //*******************************************************************************//
0016   // If you need to add multiple root outputs, by multithreading, use this command:
0017   system("hadd -O -f molecular-dna.root molecular-dna_t*.root");
0018 
0019   // Define these parameters of the simulation
0020   char ifile[256] = "molecular-dna.root";  // input filepath
0021   const Int_t numberOfPlasmids = 10144;
0022   Double_t r3 = 4.42e-6 * 4.42e-6 * 4.84e-6;  // r * r * r  (world cube side m)
0023   Double_t Nbp = 4.367 * 0.001 * numberOfPlasmids;  // Mbp // Length of the DNA chain in Mbp
0024   // multiplying by 10142 since we are testing 10142 k plasmids
0025   Double_t mass = 997 * r3;  // density * r3
0026 
0027   //*******************************************************************************//
0028 
0029   typedef std::pair<int64_t, int64_t> ipair;
0030   bool greaterPair(const ipair& l, const ipair& r);
0031   bool smallerPair(const ipair& l, const ipair& r);
0032 
0033   void BinLogX(TH1 * h);
0034 
0035   gROOT->Reset();
0036   gStyle->SetPalette(1);
0037   gROOT->SetStyle("Plain");
0038   gStyle->SetOptStat(00000);
0039 
0040   gStyle->SetPadTickX(1);
0041   gStyle->SetPadTickY(1);
0042 
0043   // Initialize output histograms
0044   TCanvas* cdamage = new TCanvas("cdamage", "DNA Damage Distribution", 900, 120, 600, 400);
0045   cdamage->SetLogy();
0046   TH1F* h1damage = new TH1F("h1damage", "h1damage", 28, 0, 14);
0047   TCanvas* c4damage = new TCanvas("c4damage", "Direct damage per plasmid", 900, 120, 600, 400);
0048   // c4damage->SetLogy();
0049   TH1F* h4damage = new TH1F("h4damage", "h4damage", 28, 0, 14);
0050   TCanvas* cSSB = new TCanvas("cSSB", "SSB Distribution", 900, 120, 600, 400);
0051   cSSB->SetLogy();
0052   TH1F* h1SSB = new TH1F("h1SSB", "h1SSB", 28, 0, 14);
0053   TCanvas* cDSB = new TCanvas("cDSB", "DSB Distribution", 900, 120, 600, 400);
0054   cDSB->SetLogy();
0055   TH1F* h1DSB = new TH1F("h1DSB", "h1DSB", 40, 0, 10);
0056   TCanvas* ccount = new TCanvas("cCount", "Damage per plasmid Distribution", 900, 120, 600, 400);
0057   TH1F* h1count = new TH1F("h1count", "h1count", 11000, 0, 11000);
0058   TGraph2D* gr1 = new TGraph2D();
0059 
0060   // Open root file
0061   TFile* f = TFile::Open(ifile);
0062 
0063   // Initialize Variables
0064   Int_t EB, ES, OHB, OHS, HB, HS, FL;
0065   Int_t total_EB, total_ES, total_OHB, total_OHS, total_HB, total_HS, total_FL;
0066   Float_t total_EB2, total_ES2, total_OHB2, total_OHS2, total_HB2, total_HS2, total_FL2;
0067   Float_t SD_EB, SD_ES, SD_OHB, SD_OHS, SD_HB, SD_HS;
0068   Float_t SD_SSB, SD_SSBp, SD_SSB2p, SD_sSSB, SD_SSBd, SD_SSBi, SD_SSBm;
0069   Float_t SD_DSB, SD_DSBp, SD_DSBpp, SD_sDSB, SD_DSBd, SD_DSBi, SD_DSBm, SD_DSBh;
0070 
0071   Int_t SSB, SSBp, SSB2p;
0072   Int_t total_SSB, total_SSBp, total_SSB2p;
0073   Float_t total_SSB2, total_SSBp2, total_SSB2p2;
0074   Int_t DSB, DSBp, DSBpp;
0075   Int_t total_DSB, total_DSBp, total_DSBpp;
0076   Float_t total_DSB2, total_DSBp2, total_DSBpp2;
0077 
0078   Int_t SSBd, SSBi, SSBm;
0079   Int_t total_sSSB, total_SSBd, total_SSBi, total_SSBm;
0080   Float_t total_sSSB2, total_SSBd2, total_SSBi2, total_SSBm2;
0081   Int_t DSBd, DSBi, DSBm, DSBh;
0082   Int_t total_sDSB, total_DSBd, total_DSBi, total_DSBm, total_DSBh;
0083   Float_t total_sDSB2, total_DSBd2, total_DSBi2, total_DSBm2, total_DSBh2;
0084 
0085   Double_t dose = 0;
0086   Double_t SD_dose = 0;
0087 
0088   Double_t EB_yield = 0;
0089   Double_t ES_yield = 0;
0090   Double_t OHB_yield = 0;
0091   Double_t OHS_yield = 0;
0092   Double_t HB_yield = 0;
0093   Double_t HS_yield = 0;
0094   Double_t SD_EB_yield = 0;
0095   Double_t SD_ES_yield = 0;
0096   Double_t SD_OHB_yield = 0;
0097   Double_t SD_OHS_yield = 0;
0098   Double_t SD_HB_yield = 0;
0099   Double_t SD_HS_yield = 0;
0100 
0101   Double_t SSB_yield = 0;
0102   Double_t SSBp_yield = 0;
0103   Double_t SSB2p_yield = 0;
0104   Double_t SD_SSB_yield = 0;
0105   Double_t SD_SSBp_yield = 0;
0106   Double_t SD_SSB2p_yield = 0;
0107   Double_t DSB_yield = 0;
0108   Double_t DSBp_yield = 0;
0109   Double_t DSBpp_yield = 0;
0110   Double_t SD_DSB_yield = 0;
0111   Double_t SD_DSBp_yield = 0;
0112   Double_t SD_DSBpp_yield = 0;
0113 
0114   Double_t sSSB_yield = 0;
0115   Double_t SSBi_yield = 0;
0116   Double_t SSBd_yield = 0;
0117   Double_t SSBm_yield = 0;
0118   Double_t SD_sSSB_yield = 0;
0119   Double_t SD_SSBi_yield = 0;
0120   Double_t SD_SSBd_yield = 0;
0121   Double_t SD_SSBm_yield = 0;
0122   Double_t sDSB_yield = 0;
0123   Double_t DSBi_yield = 0;
0124   Double_t DSBd_yield = 0;
0125   Double_t DSBm_yield = 0;
0126   Double_t DSBh_yield = 0;
0127   Double_t SD_sDSB_yield = 0;
0128   Double_t SD_DSBi_yield = 0;
0129   Double_t SD_DSBd_yield = 0;
0130   Double_t SD_DSBm_yield = 0;
0131   Double_t SD_DSBh_yield = 0;
0132 
0133   total_EB = 0;
0134   total_ES = 0;
0135   total_OHB = 0;
0136   total_OHS = 0;
0137   total_HB = 0;
0138   total_HS = 0;
0139 
0140   total_SSB = 0;
0141   total_SSBp = 0;
0142   total_SSB2p = 0;
0143   total_SSB2 = 0;
0144   total_SSBp2 = 0;
0145   total_SSB2p2 = 0;
0146   total_DSB = 0;
0147   total_DSBp = 0;
0148   total_DSBpp = 0;
0149   total_DSB2 = 0;
0150   total_DSBp2 = 0;
0151   total_DSBpp2 = 0;
0152 
0153   total_sSSB = 0;
0154   total_SSBd = 0;
0155   total_SSBi = 0;
0156   total_SSBm = 0;
0157   total_sSSB2 = 0;
0158   total_SSBd2 = 0;
0159   total_SSBi2 = 0;
0160   total_SSBm2 = 0;
0161   total_sDSB = 0;
0162   total_DSBd = 0;
0163   total_DSBi = 0;
0164   total_DSBm = 0;
0165   total_DSBh = 0;
0166   total_sDSB2 = 0;
0167   total_DSBd2 = 0;
0168   total_DSBi2 = 0;
0169   total_DSBm2 = 0;
0170   total_DSBh2 = 0;
0171 
0172   Double_t eVtoJ = 1.60218e-19;
0173   Double_t EnergyDeposited_eV = 0;
0174   Double_t acc_edep = 0;
0175   Double_t acc_edep2 = 0;
0176 
0177   Double_t Energy;
0178   Double_t BPID;
0179   Int_t Strand;
0180   Int_t StrandDamage;
0181   Int_t Event1, directB;
0182   Char_t Primary;
0183   char* primaryName = new char[32];
0184   char* type = new char[256];
0185   char* sourceClassification = new char[256];
0186 
0187   Double_t xx, yy, zz;
0188 
0189   // Read trees and leaves from root file, and give values to variables
0190   TTree* tree = (TTree*)f->Get("tuples/primary_source");
0191   Float_t number = (Float_t)tree->GetEntries();
0192 
0193   if (number<2) {
0194     std::cout << "Not enough entries in the \"primary_source\" TTree (" << (long)number << " entries)\n";
0195     gApplication->Terminate(0);
0196   }
0197   
0198   vector<pair<int, int64_t>> DSBBPID;
0199 
0200   // For reading species production
0201   tree = (TTree*)f->Get("tuples/damage");
0202   tree->SetBranchAddress("Primary", &Primary);
0203   tree->SetBranchAddress("Energy", &Energy);
0204   tree->SetBranchAddress("EaqBaseHits", &EB);
0205   tree->SetBranchAddress("EaqStrandHits", &ES);
0206   tree->SetBranchAddress("OHBaseHits", &OHB);
0207   tree->SetBranchAddress("OHStrandHits", &OHS);
0208   tree->SetBranchAddress("HBaseHits", &HB);
0209   tree->SetBranchAddress("HStrandHits", &HS);
0210   tree->SetBranchAddress("TypeClassification", type);
0211   tree->SetBranchAddress("BasePair", &BPID);
0212   tree->SetBranchAddress("Event", &Event1);
0213   tree->SetBranchAddress("Strand", &Strand);
0214   tree->SetBranchAddress("StrandDamage", &StrandDamage);
0215   tree->SetBranchAddress("Position_x_um", &xx);
0216   tree->SetBranchAddress("Position_y_um", &yy);
0217   tree->SetBranchAddress("Position_z_um", &zz);
0218   tree->SetBranchAddress("DirectBreaks", &directB);
0219 
0220   int primea = 0;
0221   Long64_t nentries = tree->GetEntries();
0222   int damagecount = 0;
0223   int hitcount = 0;
0224   // Double_t plasmidsD [25992];
0225   Double_t plasmidsD[2599200];
0226   for (int i = 0; i < nentries; i++) {
0227     tree->GetEntry(i);
0228     total_EB += EB;
0229     total_EB2 += pow(EB, 2);
0230     total_ES += ES;
0231     total_ES2 += pow(ES, 2);
0232     total_OHB += OHB;
0233     total_OHB2 += pow(OHB, 2);
0234     total_OHS += OHS;
0235     total_OHS2 += pow(OHS, 2);
0236     total_HB += HB;
0237     total_HB2 += pow(HB, 2);
0238     total_HS += HS;
0239     total_HS2 += pow(HS, 2);
0240     if ((string)type == "DSB" || (string)type == "DSB+" || (string)type == "DSB++") {
0241       DSBBPID.push_back(make_pair(i, (int64_t)BPID));
0242     }
0243     if (StrandDamage != 0) {
0244       {
0245         damagecount++;
0246         h1count->Fill(Strand);
0247       }
0248       primea++;
0249     }
0250   }
0251   int unbrokenP = 0;
0252   int brokenP = 0;
0253   for (int i = 0; i < 11000; i++) {
0254     plasmidsD[i] = h1count->GetBinContent(i);
0255     if (plasmidsD[i] == 0 && i < numberOfPlasmids) unbrokenP++;
0256   }
0257 
0258   Double_t totalDs = 0.0;
0259   for (int i = 0; i < 11000; i++) {
0260     if (plasmidsD[i] != 0) {
0261       h4damage->Fill(plasmidsD[i]);
0262       totalDs += plasmidsD[i];
0263       brokenP++;
0264     }
0265   }
0266 
0267   Double_t X;
0268   for (int i = 0; i < 28; i++) {
0269     X = h4damage->GetBinContent(i);
0270     h4damage->SetBinContent(i, X * 100 / totalDs);
0271   }
0272 
0273   // Find the number of fragments that have been produced, but first test if there are enough breaks.
0274   // If no more than 2 DSBs exist in the DSBBPID vector ( DSBBPID.size() is 0 or 1 ),
0275   // the subtraction DSBBPID.size() - 1 makes the loop condition evaluate to ie < -1 (in unsigned terms,
0276   // this becomes a large number, which is incorrect) and leads to undefined behavior (crash).
0277   if (DSBBPID.size() < 2) {
0278       std::cerr << "Not enough damage to calculate fragments." << std::endl;
0279       //return;
0280   }
0281   if (DSBBPID.size() >= 2) {
0282     // Sort DSBs from the one with lower ID value to the one with higher ID value
0283     // Then find the number of fragments that have been produced
0284     sort(DSBBPID.begin(), DSBBPID.end(), smallerPair);
0285     for(int ie = 0; ie < DSBBPID.size() - 1; ie++){
0286       int64_t dsbfragment = DSBBPID[ie + 1].second - DSBBPID[ie].second;
0287     }
0288   }
0289 
0290   // Calculate the SEM
0291   SD_EB = sqrt(abs(((total_EB2 / number) - pow(total_EB / number, 2))) / (number - 1));
0292   SD_ES = sqrt(abs(((total_ES2 / number) - pow(total_ES / number, 2))) / (number - 1));
0293   SD_OHB = sqrt(abs(((total_OHB2 / number) - pow(total_OHB / number, 2))) / (number - 1));
0294   SD_OHS = sqrt(abs(((total_OHS2 / number) - pow(total_OHS / number, 2))) / (number - 1));
0295   SD_HB = sqrt(abs(((total_HB2 / number) - pow(total_HB / number, 2))) / (number - 1));
0296   SD_HS = sqrt(abs(((total_HS2 / number) - pow(total_HS / number, 2))) / (number - 1));
0297 
0298   // Read damage classification SSB, SSB+, 2SSB, DSB, DSB+, DSB++
0299   // As they have been defined in: Nikjoo, H., O’Neill, O., Goodhead, T., & Terrissol, M. 1997,
0300   // Computational modelling of low-energy electron-induced DNA damage by early physical
0301   // and chemical events, International Journal of Radiation Biology, 71, 467.
0302   tree = (TTree*)f->Get("tuples/classification");
0303   tree->SetBranchAddress("Primary", &Primary);
0304   tree->SetBranchAddress("Energy", &Energy);
0305   tree->SetBranchAddress("SSB", &SSB);
0306   tree->SetBranchAddress("SSBp", &SSBp);
0307   tree->SetBranchAddress("2SSB", &SSB2p);
0308   tree->SetBranchAddress("DSB", &DSB);
0309   tree->SetBranchAddress("DSBp", &DSBp);
0310   tree->SetBranchAddress("DSBpp", &DSBpp);
0311 
0312   Long64_t nentriesC = tree->GetEntries();
0313   for (int i = 0; i < nentriesC; i++) {
0314     tree->GetEntry(i);
0315 
0316     total_SSBp += SSBp;
0317     total_SSBp2 += pow(SSBp, 2);
0318     total_SSB2p += SSB2p;
0319     total_SSB2p2 += pow(SSB2p, 2);
0320     total_SSB += SSB;
0321     total_SSB2 += pow(SSB, 2);
0322 
0323     total_DSBp += DSBp;
0324     total_DSBp2 += pow(DSBp, 2);
0325     total_DSBpp += DSBpp;
0326     total_DSBpp2 += pow(DSBpp, 2);
0327     total_DSB += DSB;
0328     total_DSB2 += pow(DSB, 2);
0329   }
0330 
0331   // Calculate the SEM
0332   SD_SSB = sqrt(abs(((total_SSB2 / number) - pow(total_SSB / number, 2))) / (number - 1));
0333   SD_SSBp = sqrt(abs(((total_SSBp2 / number) - pow(total_SSBp / number, 2))) / (number - 1));
0334   SD_SSB2p = sqrt(abs(((total_SSB2p2 / number) - pow(total_SSB2p / number, 2))) / (number - 1));
0335 
0336   SD_DSB = sqrt(abs(((total_DSB2 / number) - pow(total_DSB / number, 2))) / (number - 1));
0337   SD_DSBp = sqrt(abs(((total_DSBp2 / number) - pow(total_DSBp / number, 2))) / (number - 1));
0338   SD_DSBpp = sqrt(abs(((total_DSBpp2 / number) - pow(total_DSBpp / number, 2))) / (number - 1));
0339 
0340   // Read damage classification SSBd, SSBi, SSBm, DSBd, DSBi, DSBm, DSBh
0341   // As they have been defined in: Nikjoo, H., O’Neill, O., Goodhead, T., & Terrissol, M. 1997,
0342   // Computational modelling of low-energy electron-induced DNA damage by early physical
0343   // and chemical events, International Journal of Radiation Biology, 71, 467.
0344   tree = (TTree*)f->Get("tuples/source");
0345   tree->SetBranchAddress("Primary", primaryName);
0346   tree->SetBranchAddress("Energy", &Energy);
0347   tree->SetBranchAddress("SSBd", &SSBd);
0348   tree->SetBranchAddress("SSBi", &SSBi);
0349   tree->SetBranchAddress("SSBm", &SSBm);
0350   tree->SetBranchAddress("DSBd", &DSBd);
0351   tree->SetBranchAddress("DSBi", &DSBi);
0352   tree->SetBranchAddress("DSBm", &DSBm);
0353   tree->SetBranchAddress("DSBh", &DSBh);
0354 
0355   int iprime = 0;
0356 
0357   Long64_t nentriesS = tree->GetEntries();
0358   for (int i = 0; i < nentriesS; i++) {
0359     tree->GetEntry(i);
0360 
0361     total_SSBd += SSBd;
0362     total_SSBd2 += pow((SSBd), 2);
0363     total_SSBi += SSBi;
0364     total_SSBi2 += pow((SSBi), 2);
0365     total_SSBm += SSBm;
0366     total_SSBm2 += pow((SSBm), 2);
0367     total_sSSB += SSBd + SSBi + SSBm;
0368     total_sSSB2 += pow((SSBd + SSBi + SSBm), 2);
0369 
0370     total_DSBd += DSBd;
0371     total_DSBd2 += pow(DSBd, 2);
0372     total_DSBi += DSBi;
0373     total_DSBi2 += pow(DSBi, 2);
0374     total_DSBm += DSBm;
0375     total_DSBm2 += pow(DSBm, 2);
0376     total_DSBh += DSBh;
0377     total_DSBh2 += pow(DSBh, 2);
0378     total_sDSB += DSBd + DSBi + DSBm + DSBh;
0379     total_sDSB2 += pow((DSBd + DSBi + DSBm + DSBh), 2);
0380 
0381     if (SSBd != 0) h1SSB->Fill(SSBd);
0382     if (DSBd != 0) h1DSB->Fill(DSBd);
0383     if (SSBd != 0) h1damage->Fill(SSBd);
0384     if (DSBd != 0) h1damage->Fill(DSBd);
0385     if (SSB != 0 || DSBd != 0) {
0386       iprime++;
0387     }
0388   }
0389 
0390   Double_t Y;
0391   for (int i = 0; i < 28; i++) {
0392     Y = h1damage->GetBinContent(i);
0393     h1damage->SetBinContent(i, Y * 100 / totalDs);
0394   }
0395 
0396   // Calculate the SEM
0397   SD_sSSB = sqrt(abs(((total_sSSB2 / number) - pow(total_sSSB / number, 2))) / (number - 1));
0398   SD_SSBd = sqrt(abs(((total_SSBd2 / number) - pow(total_SSBd / number, 2))) / (number - 1));
0399   SD_SSBi = sqrt(abs(((total_SSBi2 / number) - pow(total_SSBi / number, 2))) / (number - 1));
0400   SD_SSBm = sqrt(abs(((total_SSBm2 / number) - pow(total_SSBm / number, 2))) / (number - 1));
0401 
0402   SD_sDSB = sqrt(abs(((total_sDSB2 / number) - pow(total_sDSB / number, 2))) / (number - 1));
0403   SD_DSBd = sqrt(abs(((total_DSBd2 / number) - pow(total_DSBd / number, 2))) / (number - 1));
0404   SD_DSBi = sqrt(abs(((total_DSBi2 / number) - pow(total_DSBi / number, 2))) / (number - 1));
0405   SD_DSBm = sqrt(abs(((total_DSBm2 / number) - pow(total_DSBm / number, 2))) / (number - 1));
0406   SD_DSBh = sqrt(abs(((total_DSBh2 / number) - pow(total_DSBh / number, 2))) / (number - 1));
0407 
0408   // Measure the Deposited Energy in the whole volume that includes DNA chain
0409 
0410   tree = (TTree*)f->Get("tuples/chromosome_hits");
0411   tree->SetBranchAddress("e_chromosome_kev", &EnergyDeposited_eV);
0412   nentries = tree->GetEntries();
0413   for (int i = 0; i < nentries; i++) {
0414     tree->GetEntry(i);
0415     acc_edep += EnergyDeposited_eV * 1e3;
0416     acc_edep2 += EnergyDeposited_eV * EnergyDeposited_eV * 1e6;
0417   }
0418   tree->SetBranchAddress("e_dna_kev", &EnergyDeposited_eV);
0419   nentries = tree->GetEntries();
0420   for (int i = 0; i < nentries; i++) {
0421     tree->GetEntry(i);
0422     acc_edep += (EnergyDeposited_eV * 1e3);
0423     acc_edep2 += (EnergyDeposited_eV * EnergyDeposited_eV * 1e6);
0424   }
0425 
0426   // Close the root file to free space
0427   f->Close();
0428   // Calculate the absorbed dose
0429   dose = acc_edep * eVtoJ / mass;
0430 
0431   double norm = 1;
0432   // Calculate the yields, together with their error
0433   EB_yield = (Double_t)total_EB / dose / Nbp;
0434   ES_yield = (Double_t)total_ES / dose / Nbp;
0435   OHB_yield = (Double_t)total_OHB / dose / Nbp;
0436   OHS_yield = (Double_t)total_OHS / dose / Nbp;
0437   HB_yield = (Double_t)total_HB / dose / Nbp;
0438   HS_yield = (Double_t)total_HS / dose / Nbp;
0439 
0440   SD_EB_yield = SD_EB / dose / Nbp;
0441   SD_ES_yield = SD_ES / dose / Nbp;
0442   SD_OHB_yield = SD_OHB / dose / Nbp;
0443   SD_OHS_yield = SD_OHS / dose / Nbp;
0444   SD_HB_yield = SD_HB / dose / Nbp;
0445   SD_HS_yield = SD_HS / dose / Nbp;
0446 
0447   SSB_yield = (Double_t)norm * total_SSB / dose / Nbp;
0448   SSBp_yield = (Double_t)norm * total_SSBp / dose / Nbp;
0449   SSB2p_yield = (Double_t)norm * total_SSB2p / dose / Nbp;
0450 
0451   DSB_yield = (Double_t)norm * total_DSB / dose / Nbp;
0452   DSBp_yield = (Double_t)norm * total_DSBp / dose / Nbp;
0453   DSBpp_yield = (Double_t)norm * total_DSBpp / dose / Nbp;
0454 
0455   SD_SSB_yield = norm * SD_SSB / dose / Nbp;
0456   SD_SSBp_yield = norm * SD_SSBp / dose / Nbp;
0457   SD_SSB2p_yield = norm * SD_SSB2p / dose / Nbp;
0458 
0459   SD_DSB_yield = norm * SD_DSB / dose / Nbp;
0460   SD_DSBp_yield = norm * SD_DSBp / dose / Nbp;
0461   SD_DSBpp_yield = norm * SD_DSBpp / dose / Nbp;
0462 
0463   sSSB_yield = (Double_t)norm * total_sSSB / dose / Nbp;
0464   SSBi_yield = (Double_t)norm * total_SSBi / dose / Nbp;
0465   SSBd_yield = (Double_t)norm * total_SSBd / dose / Nbp;
0466   SSBm_yield = (Double_t)norm * total_SSBm / dose / Nbp;
0467 
0468   sDSB_yield = (Double_t)norm * total_sDSB / dose / Nbp;
0469   DSBi_yield = (Double_t)norm * total_DSBi / dose / Nbp;
0470   DSBd_yield = (Double_t)norm * total_DSBd / dose / Nbp;
0471   DSBm_yield = (Double_t)norm * total_DSBm / dose / Nbp;
0472   DSBh_yield = (Double_t)norm * total_DSBh / dose / Nbp;
0473 
0474   SD_sSSB_yield = norm * SD_sSSB / dose / Nbp;
0475   SD_SSBi_yield = norm * SD_SSBi / dose / Nbp;
0476   SD_SSBd_yield = norm * SD_SSBd / dose / Nbp;
0477   SD_SSBm_yield = norm * SD_SSBm / dose / Nbp;
0478 
0479   SD_sDSB_yield = norm * SD_sDSB / dose / Nbp;
0480   SD_DSBi_yield = norm * SD_DSBi / dose / Nbp;
0481   SD_DSBd_yield = norm * SD_DSBd / dose / Nbp;
0482   SD_DSBm_yield = norm * SD_DSBm / dose / Nbp;
0483   SD_DSBh_yield = norm * SD_DSBh / dose / Nbp;
0484 
0485   // Print output in terminal
0486 
0487   float total_SSB_totalYield = SSB_yield + SSBp_yield + SSB2p_yield;
0488   float total_DSB_totalYield = DSB_yield + DSBp_yield + DSBpp_yield;
0489 
0490   cout << "\n"
0491        << " Output file: " << ifile << '\n'
0492        << "\nDose Absorbed (Gy): " << dose << '\n'
0493        << "Particle : " << primaryName << '\t' << "Energy (MeV) : " << Energy << '\t'
0494        << "Number of Primaries : " << number << '\n'
0495 
0496        << "  Output Damage : " << '\n'
0497 
0498        << '\t' << "Number of plasmids:         " << numberOfPlasmids << "   \t" << '\n'
0499        << '\t' << "Unbroken plasmids:          " << unbrokenP << "   \t" << '\n'
0500        << '\t' << "Broken plasmids:            " << brokenP << "   \t" << '\n'
0501        << '\t' << "Rate of damages per plasmid " << totalDs / numberOfPlasmids << "   \t" << '\n'
0502 
0503        << '\t' << "Rate of damages per plasmid per dose " << '\t'
0504        << totalDs / numberOfPlasmids / dose << " plasmid-1Gy-1" << '\n'
0505        << '\t' << "Rate of damages per dose per Mbp " << totalDs / dose / Nbp << " Gy-1Mbp-1"
0506        << '\n'
0507 
0508        << '\t' << "     Species Hits " << '\n'
0509        << '\t' << "EaqBaseHits    " << EB_yield * dose * Nbp << "   \t" << '\n'
0510        << '\t' << "EaqStrandHits  " << ES_yield * dose * Nbp << "   \t" << '\n'
0511        << '\t' << "OHBaseHits     " << OHB_yield * dose * Nbp << "   \t" << '\n'
0512        << '\t' << "OHStrandHits   " << OHS_yield * dose * Nbp << "   \t" << '\n'
0513        << '\t' << "HBaseHits      " << HB_yield * dose * Nbp << "   \t" << '\n'
0514        << '\t' << "HStrandHits    " << HS_yield * dose * Nbp << "   \t" << '\n'
0515        << '\n'
0516        << '\t' << "     Damage number" << '\n'
0517        << '\t' << "SSB            " << SSB_yield * dose * Nbp << "   \t" << '\n'
0518        << '\t' << "SSB+           " << SSBp_yield * dose * Nbp << "   \t" << '\n'
0519        << '\t' << "2SSB           " << SSB2p_yield * dose * Nbp << "   \t" << '\n'
0520        << '\t' << "SSB total      " << total_SSB_totalYield * dose * Nbp << '\n'
0521        << '\t' << "DSB            " << DSB_yield * dose * Nbp << "   \t" << '\n'
0522        << '\t' << "DSB+           " << DSBp_yield * dose * Nbp << "   \t" << '\n'
0523        << '\t' << "DSB++          " << DSBpp_yield * dose * Nbp << "   \t" << '\n'
0524        << '\t' << "DSB total      " << total_DSB_totalYield * dose * Nbp << '\n'
0525        << '\n'
0526        << '\t' << "     Breaks number " << '\n'
0527        << '\t' << "SSB direct     " << SSBd_yield * dose * Nbp << "   \t" << '\n'
0528        << '\t' << "SSB indirect   " << SSBi_yield * dose * Nbp << "   \t" << '\n'
0529        << '\t' << "SSB mixed      " << SSBm_yield * dose * Nbp << "   \t" << '\n'
0530        << '\t' << "SSB total      " << sSSB_yield * dose * Nbp << "   \t" << '\n'
0531        << '\t' << "DSB direct     " << DSBd_yield * dose * Nbp << "   \t" << '\n'
0532        << '\t' << "DSB indirect   " << DSBi_yield * dose * Nbp << "   \t" << '\n'
0533        << '\t' << "DSB mixed      " << DSBm_yield * dose * Nbp << "   \t" << '\n'
0534        << '\t' << "DSB hybrid     " << DSBh_yield * dose * Nbp << "   \t" << '\n'
0535        << '\t' << "DSB total      " << sDSB_yield * dose * Nbp << "   \t" << '\n'
0536        << '\n'
0537        << '\t' << "SSB/DSB        " << sSSB_yield / sDSB_yield << '\n'
0538        << '\n';
0539 
0540   // Plot Histograms
0541 
0542   cSSB->GetCanvas()->cd();
0543   h1SSB->SetStats(false);
0544   h1SSB->SetMarkerSize(0.1);
0545   h1SSB->SetMarkerColor(kBlue);
0546   h1SSB->SetLineColor(kBlue);
0547   h1SSB->SetTitle("");
0548   h1SSB->SetYTitle(" ");
0549   h1SSB->SetXTitle("Number of direct SSB per event");
0550   h1SSB->SetFillColor(kBlue);
0551   h1SSB->Draw();
0552 
0553   cDSB->GetCanvas()->cd();
0554   h1DSB->SetStats(false);
0555   h1DSB->SetMarkerSize(0.1);
0556   h1DSB->SetMarkerColor(kGreen + 2);
0557   h1DSB->SetLineColor(kGreen + 2);
0558   h1DSB->SetTitle("");
0559   h1DSB->SetYTitle(" ");
0560   h1DSB->SetXTitle("Number of direct DSB per event");
0561   h1DSB->SetFillColor(kGreen + 2);
0562   h1DSB->Draw();
0563 
0564   cdamage->GetCanvas()->cd();
0565   h1damage->SetStats(false);
0566   h1damage->SetMarkerSize(0.1);
0567   h1damage->SetMarkerColor(kMagenta);
0568   h1damage->SetLineColor(kMagenta);
0569   h1damage->SetTitle("");
0570   h1damage->SetYTitle("Percentage % ");
0571   h1damage->SetXTitle("Number of damages per event");
0572   h1damage->SetFillColor(kMagenta);
0573   h1damage->Draw();
0574 
0575   c4damage->GetCanvas()->cd();
0576   h4damage->SetStats(false);
0577   h4damage->SetMarkerSize(0.1);
0578   h4damage->SetMarkerColor(kRed - 4);
0579   h4damage->SetLineColor(kRed - 4);
0580   h4damage->SetTitle("");
0581   h4damage->SetYTitle("Percentage % ");
0582   h4damage->SetXTitle("Number of damages per plasmid");
0583   h4damage->SetFillColor(kRed - 4);
0584   h4damage->Draw();
0585 
0586   ccount->GetCanvas()->cd();
0587   h1count->SetStats(false);
0588   h1count->SetMarkerSize(0.1);
0589   h1count->SetMarkerColor(kCyan);
0590   h1count->SetLineColor(kCyan);
0591   h1count->SetTitle("");
0592   h1count->SetYTitle("Number of damages");
0593   h1count->SetXTitle("Plasmid ID");
0594   h1count->Draw();
0595 }
0596 
0597 // Some important bools that are needed to run the root macro file
0598 bool greaterPair(const ipair& l, const ipair& r)
0599 {
0600   return l.second > r.second;
0601 }
0602 bool smallerPair(const ipair& l, const ipair& r)
0603 {
0604   return l.second < r.second;
0605 }
0606 
0607 void BinLogX(TH1* h)
0608 {
0609   TAxis* axis = h->GetXaxis();
0610   int bins = axis->GetNbins();
0611   Axis_t from = axis->GetXmin();
0612   Axis_t to = axis->GetXmax();
0613   Axis_t width = (to - from) / bins;
0614   Axis_t* new_bins = new Axis_t[bins + 1];
0615   for (int i = 0; i <= bins; i++) {
0616     new_bins[i] = TMath::Power(10, from + i * width);
0617   }
0618   axis->Set(bins, new_bins);
0619   delete[] new_bins;
0620 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!