Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/dna/moleculardna/human_cell_alphas.C

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

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