| Geant4 Cross Reference |
1 {
2 gROOT->Reset();
3 gStyle->SetPalette(1);
4 gROOT->SetStyle("Plain");
5 gStyle->SetOptStat(00000);
6
7 system ("hadd -O -f molecular-dna.root molecular-dna_t*.root");
8
9 c1 = new TCanvas("c1", "Damages", 120, 60, 1000, 1000);
10 c1->SetBorderSize(0);
11 c1->SetFillColor(0);
12 c1->SetFillStyle(4000);
13 c1->Divide(2,1);
14 gPad->SetLeftMargin(0.13);
15
16 Int_t SSBd, SSBi, SSBm, DSBd, DSBi, DSBm, DSBh;
17 Int_t total_SSBd, total_SSBi, total_SSBm, total_DSBd, total_DSBi, total_DSBm, total_DSBh;
18 Float_t total_SSBd2, total_SSBi2, total_SSBm2, total_DSBd2, total_DSBi2, total_DSBm2, total_DSBh2;
19 Float_t mean_SSBd, mean_SSBi, mean_SSBm, mean_DSBd, mean_DSBi, mean_DSBm, mean_DSBh;
20
21 Int_t SSB, SSBp, twoSSB, DSB, DSBp, DSBpp;
22 Int_t total_SSB, total_SSBp, total_twoSSB, total_DSB, total_DSBp, total_DSBpp;
23 Int_t total_SSB2, total_SSBp2, total_twoSSB2, total_DSB2, total_DSBp2,
24 total_DSBpp2;
25 Float_t mean_SSB, mean_SSBp, mean_twoSSB, mean_DSB, mean_DSBp, mean_DSBpp;
26
27 total_SSBd = 0;
28 total_SSBi = 0;
29 total_SSBm = 0;
30 total_DSBd = 0;
31 total_DSBi = 0;
32 total_DSBm = 0;
33 total_DSBh = 0;
34
35 total_SSBd2 = 0;
36 total_SSBi2 = 0;
37 total_SSBm2 = 0;
38 total_DSBd2 = 0;
39 total_DSBi2 = 0;
40 total_DSBm2 = 0;
41 total_DSBh2 = 0;
42
43 total_SSB = 0;
44 total_SSBp = 0;
45 total_twoSSB = 0;
46 total_DSB = 0;
47 total_DSBp = 0;
48 total_DSBpp = 0;
49
50 total_SSB2 = 0;
51 total_SSBp2 = 0;
52 total_twoSSB2 = 0;
53 total_DSB2 = 0;
54 total_DSBp2 = 0;
55 total_DSBpp2 = 0;
56
57 Double_t EnergyDeposited_eV = 0;
58 Double_t acc_edep = 0;
59 Double_t acc_edep2 = 0;
60
61 Double_t Energy;
62 Char_t Primary;
63
64 TFile *f = TFile::Open("molecular-dna.root");
65 TTree *tree = (TTree *) f->Get("tuples/primary_source");
66 Float_t number = (Float_t) tree->GetEntries();
67
68 tree = (TTree *) f->Get("tuples/source");
69 tree->SetBranchAddress("Primary",&Primary);
70 tree->SetBranchAddress("Energy",&Energy);
71 tree->SetBranchAddress("SSBd",&SSBd);
72 tree->SetBranchAddress("SSBi",&SSBi);
73 tree->SetBranchAddress("SSBm",&SSBm);
74 tree->SetBranchAddress("DSBd",&DSBd);
75 tree->SetBranchAddress("DSBi",&DSBi);
76 tree->SetBranchAddress("DSBm",&DSBm);
77 tree->SetBranchAddress("DSBh",&DSBh);
78
79 Long64_t nentries = tree->GetEntries();
80
81 for(int i = 0;i<nentries;i++){
82 tree->GetEntry(i);
83 total_SSBd += SSBd;
84 total_SSBd2 += SSBd *SSBd;
85 total_SSBi += SSBi;
86 total_SSBi2 += SSBi *SSBi;
87 total_SSBm += SSBm;
88 total_SSBm2 += SSBm *SSBm;
89
90 total_DSBd += DSBd;
91 total_DSBd2 += DSBd *DSBd;
92 total_DSBi += DSBi;
93 total_DSBi2 += DSBi *DSBi;
94 total_DSBm += DSBm;
95 total_DSBm2 += DSBm *DSBm;
96 total_DSBh += DSBh;
97 total_DSBh2 += DSBh *DSBh;
98 }
99
100 tree = (TTree *) f->Get("tuples/damage");
101 tree->SetBranchAddress("EnergyDeposited_eV",&EnergyDeposited_eV);
102 nentries = tree->GetEntries();
103 for(int i = 0;i<nentries;i++){
104 tree->GetEntry(i);
105 acc_edep += EnergyDeposited_eV;
106 acc_edep2 += EnergyDeposited_eV *EnergyDeposited_eV;
107 }
108
109 tree = (TTree *) f->Get("tuples/classification");
110 tree->SetBranchAddress("SSB",&SSB);
111 tree->SetBranchAddress("SSBp",&SSBp);
112 tree->SetBranchAddress("2SSB",&twoSSB);
113 tree->SetBranchAddress("DSB",&DSB);
114 tree->SetBranchAddress("DSBp",&DSBp);
115 tree->SetBranchAddress("DSBpp",&DSBpp);
116
117 nentries = tree->GetEntries();
118
119 for(int i = 0;i<nentries;i++){
120 tree->GetEntry(i);
121 total_SSB += SSB;
122 total_SSB2 += SSB *SSB;
123 total_SSBp += SSBp;
124 total_SSBp2 += SSBp *SSBp;
125 total_twoSSB += twoSSB;
126 total_twoSSB2 += twoSSB *twoSSB;
127
128 total_DSB += DSB;
129 total_DSB2 += DSB *DSB;
130 total_DSBp += DSBp;
131 total_DSBp2 += DSBp *DSBp;
132 total_DSBpp += DSBpp;
133 total_DSBpp2 += DSBpp *DSBpp;
134 }
135
136 mean_SSBd = (Float_t) total_SSBd / number;
137 mean_SSBi = (Float_t) total_SSBi / number;
138 mean_SSBm = (Float_t) total_SSBm / number;
139
140 mean_DSBd = (Float_t) total_DSBd / number;
141 mean_DSBi = (Float_t) total_DSBi / number;
142 mean_DSBm = (Float_t) total_DSBm / number;
143 mean_DSBh = (Float_t) total_DSBh / number;
144
145 Double_t SD_SSBd = sqrt(((total_SSBd2 / number) - pow(total_SSBd / number,2))/(number -1));
146 Double_t SD_SSBi = sqrt(((total_SSBi2 / number) - pow(total_SSBi / number,2))/(number -1));
147 Double_t SD_SSBm = sqrt(((total_SSBm2 / number) - pow(total_SSBm / number,2))/(number -1));
148
149 Double_t SD_DSBd = sqrt(((total_DSBd2 / number) - pow(total_DSBd / number,2))/(number -1));
150 Double_t SD_DSBi = sqrt(((total_DSBi2 / number) - pow(total_DSBi / number,2))/(number -1));
151 Double_t SD_DSBm = sqrt(((total_DSBm2 / number) - pow(total_DSBm / number,2))/(number -1));
152 Double_t SD_DSBh = sqrt(((total_DSBh2 / number) - pow(total_DSBh / number,2))/(number -1));
153
154
155 mean_SSB = (Float_t) total_SSB / number;
156 mean_SSBp = (Float_t) total_SSBp / number;
157 mean_twoSSB = (Float_t) total_twoSSB / number;
158
159 mean_DSB = (Float_t) total_DSB / number;
160 mean_DSBp = (Float_t) total_DSBp / number;
161 mean_DSBpp = (Float_t) total_DSBpp / number;
162
163 Double_t SD_SSB = sqrt(((total_SSB2 / number) - pow(total_SSB / number,2))/(number -1));
164 Double_t SD_SSBp = sqrt(((total_SSBp2 / number) - pow(total_SSBp / number,2))
165 /(number -1));
166 Double_t SD_twoSSB = sqrt(((total_twoSSB2 / number) - pow(total_twoSSB /
167 number,2))
168 /(number -1));
169
170 Double_t SD_DSB = sqrt(((total_DSB2 / number) - pow(total_DSB / number,2))/(number -1));
171 Double_t SD_DSBp = sqrt(((total_DSBp2 / number) - pow(total_DSBp / number,2))
172 /(number -1));
173 Double_t SD_DSBpp = sqrt(((total_DSBpp2 / number) - pow(total_DSBpp / number,2))
174 /(number -1));
175
176 cout<<"Paricle : "<<Primary<<'\t'
177 <<"Energy [/MeV] : "<<Energy<<'\t'
178 <<"number : "<<number<<'\n'
179 <<" Output Damages : "<<'\n'<<
180 '\t'<<"SSB direct : "<<mean_SSBd<<'\t'<<" error : "<<SD_SSBd<<'\n'<<
181 '\t'<<"SSB indirect : "<<mean_SSBi<<'\t'<<" error : "<<SD_SSBi<<'\n'<<
182 '\t'<<"SSB mix : "<<mean_SSBm<<'\t'<<" error : "<<SD_SSBm<<'\n'<<'\n'<<
183
184 '\t'<<"DSB direct : "<<mean_DSBd<<'\t'<<" error : "<<SD_DSBd<<'\n'<<
185 '\t'<<"DSB indirect : "<<mean_DSBi<<'\t'<<" error : "<<SD_DSBi<<'\n'<<
186 '\t'<<"DSB mix : "<<mean_DSBm<<'\t'<<" error : "<<SD_DSBm<<'\n'<<
187 '\t'<<"DSB hybrid : "<<mean_DSBh<<'\t'<<" error : "<<SD_DSBh<<'\n'<<'\n'<<
188
189 '\t'<<"SSB : "<<mean_SSB<<'\t'<<" error : "<<SD_SSB<<'\n'<<
190 '\t'<<"SSBp : "<<mean_SSBp<<'\t'<<" error : "<<SD_SSBp<<'\n'<<
191 '\t'<<"2SSB : "<<mean_twoSSB<<'\t'<<" error : "<<SD_twoSSB<<'\n'<<'\n'<<
192
193 '\t'<<"DSB : "<<mean_DSB<<'\t'<<" error : "<<SD_DSB<<'\n'<<
194 '\t'<<"DSBp : "<<mean_DSBp<<'\t'<<" error : "<<SD_DSBp<<'\n'<<
195 '\t'<<"DSBpp : "<<mean_DSBpp<<'\t'<<" error : "<<SD_DSBpp<<'\n';
196
197
198 f->Close();
199 const Int_t n1 = 7;
200 Double_t x1[n1] = {2,4,6,8,10,12,14};
201 Double_t _y1[n1] = {mean_SSBd,mean_SSBi,mean_SSBm,mean_DSBd,mean_DSBi,
202 mean_DSBm,mean_DSBh};
203 Double_t err_y1[n1] = {SD_SSBd,SD_SSBi,SD_SSBm,SD_DSBd,SD_DSBi,SD_DSBm,SD_DSBh};
204 TGraph* gr1 = new TGraphErrors(n1,x1,_y1,0,err_y1);
205 gr1->SetTitle("Break Source Frequency");
206 gr1->GetXaxis()->SetBinLabel(10,"SSB direct");
207 gr1->GetXaxis()->SetBinLabel(20,"SSB indirect");
208 gr1->GetXaxis()->SetBinLabel(35,"SSB mix");
209 gr1->GetXaxis()->SetBinLabel(50,"DSB direct");
210 gr1->GetXaxis()->SetBinLabel(65,"DSB indirect");
211 gr1->GetXaxis()->SetBinLabel(80,"DSB mix");
212 gr1->GetXaxis()->SetBinLabel(90,"DSB hybrid");
213 gr1->SetFillColor(49);
214
215 const Int_t n2 = 6;
216 Double_t x2[n2] = {2,4,6,8,10,12};
217 Double_t _y2[n2] = {mean_SSB,mean_SSBp,mean_twoSSB,mean_DSB,mean_DSBp,
218 mean_DSBpp};
219 Double_t err_y2[n2] = {SD_SSB,SD_SSBp,SD_twoSSB,SD_DSB,SD_DSBp,SD_DSBpp};
220 TGraph* gr2 = new TGraphErrors(n2,x2,_y2,0,err_y2);
221 gr2->SetTitle("Break Complexity Frequency");
222 gr2->GetXaxis()->SetBinLabel(10,"SSB");
223 gr2->GetXaxis()->SetBinLabel(25,"SSBp");
224 gr2->GetXaxis()->SetBinLabel(45,"2SSB");
225 gr2->GetXaxis()->SetBinLabel(60,"DSB");
226 gr2->GetXaxis()->SetBinLabel(75,"DSBp");
227 gr2->GetXaxis()->SetBinLabel(90,"DSBpp");
228 gr2->SetFillColor(49);
229
230 //Draw
231 c1->cd(1);
232 gr1->Draw("ba");
233
234 c1->cd(2);
235 gr2->Draw("ba");
236 }
237