File indexing completed on 2025-02-23 09:22:02
0001 {
0002 gROOT->Reset();
0003 gStyle->SetPalette(1);
0004 gROOT->SetStyle("Plain");
0005 gStyle->SetOptStat(00000);
0006
0007 system ("hadd -O -f molecular-dna.root molecular-dna_t*.root");
0008
0009 c1 = new TCanvas("c1", "Damages", 120, 60, 1000, 1000);
0010 c1->SetBorderSize(0);
0011 c1->SetFillColor(0);
0012 c1->SetFillStyle(4000);
0013 c1->Divide(2,1);
0014 gPad->SetLeftMargin(0.13);
0015
0016 Int_t SSBd, SSBi, SSBm, DSBd, DSBi, DSBm, DSBh;
0017 Int_t total_SSBd, total_SSBi, total_SSBm, total_DSBd, total_DSBi, total_DSBm, total_DSBh;
0018 Float_t total_SSBd2, total_SSBi2, total_SSBm2, total_DSBd2, total_DSBi2, total_DSBm2, total_DSBh2;
0019 Float_t mean_SSBd, mean_SSBi, mean_SSBm, mean_DSBd, mean_DSBi, mean_DSBm, mean_DSBh;
0020
0021 Int_t SSB, SSBp, twoSSB, DSB, DSBp, DSBpp;
0022 Int_t total_SSB, total_SSBp, total_twoSSB, total_DSB, total_DSBp, total_DSBpp;
0023 Int_t total_SSB2, total_SSBp2, total_twoSSB2, total_DSB2, total_DSBp2,
0024 total_DSBpp2;
0025 Float_t mean_SSB, mean_SSBp, mean_twoSSB, mean_DSB, mean_DSBp, mean_DSBpp;
0026
0027 total_SSBd = 0;
0028 total_SSBi = 0;
0029 total_SSBm = 0;
0030 total_DSBd = 0;
0031 total_DSBi = 0;
0032 total_DSBm = 0;
0033 total_DSBh = 0;
0034
0035 total_SSBd2 = 0;
0036 total_SSBi2 = 0;
0037 total_SSBm2 = 0;
0038 total_DSBd2 = 0;
0039 total_DSBi2 = 0;
0040 total_DSBm2 = 0;
0041 total_DSBh2 = 0;
0042
0043 total_SSB = 0;
0044 total_SSBp = 0;
0045 total_twoSSB = 0;
0046 total_DSB = 0;
0047 total_DSBp = 0;
0048 total_DSBpp = 0;
0049
0050 total_SSB2 = 0;
0051 total_SSBp2 = 0;
0052 total_twoSSB2 = 0;
0053 total_DSB2 = 0;
0054 total_DSBp2 = 0;
0055 total_DSBpp2 = 0;
0056
0057 Double_t EnergyDeposited_eV = 0;
0058 Double_t acc_edep = 0;
0059 Double_t acc_edep2 = 0;
0060
0061 Double_t Energy;
0062 Char_t Primary;
0063
0064 TFile *f = TFile::Open("molecular-dna.root");
0065 TTree *tree = (TTree *) f->Get("tuples/primary_source");
0066 Float_t number = (Float_t) tree->GetEntries();
0067
0068 tree = (TTree *) f->Get("tuples/source");
0069 tree->SetBranchAddress("Primary",&Primary);
0070 tree->SetBranchAddress("Energy",&Energy);
0071 tree->SetBranchAddress("SSBd",&SSBd);
0072 tree->SetBranchAddress("SSBi",&SSBi);
0073 tree->SetBranchAddress("SSBm",&SSBm);
0074 tree->SetBranchAddress("DSBd",&DSBd);
0075 tree->SetBranchAddress("DSBi",&DSBi);
0076 tree->SetBranchAddress("DSBm",&DSBm);
0077 tree->SetBranchAddress("DSBh",&DSBh);
0078
0079 Long64_t nentries = tree->GetEntries();
0080
0081 for(int i = 0;i<nentries;i++){
0082 tree->GetEntry(i);
0083 total_SSBd += SSBd;
0084 total_SSBd2 += SSBd *SSBd;
0085 total_SSBi += SSBi;
0086 total_SSBi2 += SSBi *SSBi;
0087 total_SSBm += SSBm;
0088 total_SSBm2 += SSBm *SSBm;
0089
0090 total_DSBd += DSBd;
0091 total_DSBd2 += DSBd *DSBd;
0092 total_DSBi += DSBi;
0093 total_DSBi2 += DSBi *DSBi;
0094 total_DSBm += DSBm;
0095 total_DSBm2 += DSBm *DSBm;
0096 total_DSBh += DSBh;
0097 total_DSBh2 += DSBh *DSBh;
0098 }
0099
0100 tree = (TTree *) f->Get("tuples/damage");
0101 tree->SetBranchAddress("EnergyDeposited_eV",&EnergyDeposited_eV);
0102 nentries = tree->GetEntries();
0103 for(int i = 0;i<nentries;i++){
0104 tree->GetEntry(i);
0105 acc_edep += EnergyDeposited_eV;
0106 acc_edep2 += EnergyDeposited_eV *EnergyDeposited_eV;
0107 }
0108
0109 tree = (TTree *) f->Get("tuples/classification");
0110 tree->SetBranchAddress("SSB",&SSB);
0111 tree->SetBranchAddress("SSBp",&SSBp);
0112 tree->SetBranchAddress("2SSB",&twoSSB);
0113 tree->SetBranchAddress("DSB",&DSB);
0114 tree->SetBranchAddress("DSBp",&DSBp);
0115 tree->SetBranchAddress("DSBpp",&DSBpp);
0116
0117 nentries = tree->GetEntries();
0118
0119 for(int i = 0;i<nentries;i++){
0120 tree->GetEntry(i);
0121 total_SSB += SSB;
0122 total_SSB2 += SSB *SSB;
0123 total_SSBp += SSBp;
0124 total_SSBp2 += SSBp *SSBp;
0125 total_twoSSB += twoSSB;
0126 total_twoSSB2 += twoSSB *twoSSB;
0127
0128 total_DSB += DSB;
0129 total_DSB2 += DSB *DSB;
0130 total_DSBp += DSBp;
0131 total_DSBp2 += DSBp *DSBp;
0132 total_DSBpp += DSBpp;
0133 total_DSBpp2 += DSBpp *DSBpp;
0134 }
0135
0136 mean_SSBd = (Float_t) total_SSBd / number;
0137 mean_SSBi = (Float_t) total_SSBi / number;
0138 mean_SSBm = (Float_t) total_SSBm / number;
0139
0140 mean_DSBd = (Float_t) total_DSBd / number;
0141 mean_DSBi = (Float_t) total_DSBi / number;
0142 mean_DSBm = (Float_t) total_DSBm / number;
0143 mean_DSBh = (Float_t) total_DSBh / number;
0144
0145 Double_t SD_SSBd = sqrt(((total_SSBd2 / number) - pow(total_SSBd / number,2))/(number -1));
0146 Double_t SD_SSBi = sqrt(((total_SSBi2 / number) - pow(total_SSBi / number,2))/(number -1));
0147 Double_t SD_SSBm = sqrt(((total_SSBm2 / number) - pow(total_SSBm / number,2))/(number -1));
0148
0149 Double_t SD_DSBd = sqrt(((total_DSBd2 / number) - pow(total_DSBd / number,2))/(number -1));
0150 Double_t SD_DSBi = sqrt(((total_DSBi2 / number) - pow(total_DSBi / number,2))/(number -1));
0151 Double_t SD_DSBm = sqrt(((total_DSBm2 / number) - pow(total_DSBm / number,2))/(number -1));
0152 Double_t SD_DSBh = sqrt(((total_DSBh2 / number) - pow(total_DSBh / number,2))/(number -1));
0153
0154
0155 mean_SSB = (Float_t) total_SSB / number;
0156 mean_SSBp = (Float_t) total_SSBp / number;
0157 mean_twoSSB = (Float_t) total_twoSSB / number;
0158
0159 mean_DSB = (Float_t) total_DSB / number;
0160 mean_DSBp = (Float_t) total_DSBp / number;
0161 mean_DSBpp = (Float_t) total_DSBpp / number;
0162
0163 Double_t SD_SSB = sqrt(((total_SSB2 / number) - pow(total_SSB / number,2))/(number -1));
0164 Double_t SD_SSBp = sqrt(((total_SSBp2 / number) - pow(total_SSBp / number,2))
0165 /(number -1));
0166 Double_t SD_twoSSB = sqrt(((total_twoSSB2 / number) - pow(total_twoSSB /
0167 number,2))
0168 /(number -1));
0169
0170 Double_t SD_DSB = sqrt(((total_DSB2 / number) - pow(total_DSB / number,2))/(number -1));
0171 Double_t SD_DSBp = sqrt(((total_DSBp2 / number) - pow(total_DSBp / number,2))
0172 /(number -1));
0173 Double_t SD_DSBpp = sqrt(((total_DSBpp2 / number) - pow(total_DSBpp / number,2))
0174 /(number -1));
0175
0176 cout<<"Paricle : "<<Primary<<'\t'
0177 <<"Energy [/MeV] : "<<Energy<<'\t'
0178 <<"number : "<<number<<'\n'
0179 <<" Output Damages : "<<'\n'<<
0180 '\t'<<"SSB direct : "<<mean_SSBd<<'\t'<<" error : "<<SD_SSBd<<'\n'<<
0181 '\t'<<"SSB indirect : "<<mean_SSBi<<'\t'<<" error : "<<SD_SSBi<<'\n'<<
0182 '\t'<<"SSB mix : "<<mean_SSBm<<'\t'<<" error : "<<SD_SSBm<<'\n'<<'\n'<<
0183
0184 '\t'<<"DSB direct : "<<mean_DSBd<<'\t'<<" error : "<<SD_DSBd<<'\n'<<
0185 '\t'<<"DSB indirect : "<<mean_DSBi<<'\t'<<" error : "<<SD_DSBi<<'\n'<<
0186 '\t'<<"DSB mix : "<<mean_DSBm<<'\t'<<" error : "<<SD_DSBm<<'\n'<<
0187 '\t'<<"DSB hybrid : "<<mean_DSBh<<'\t'<<" error : "<<SD_DSBh<<'\n'<<'\n'<<
0188
0189 '\t'<<"SSB : "<<mean_SSB<<'\t'<<" error : "<<SD_SSB<<'\n'<<
0190 '\t'<<"SSBp : "<<mean_SSBp<<'\t'<<" error : "<<SD_SSBp<<'\n'<<
0191 '\t'<<"2SSB : "<<mean_twoSSB<<'\t'<<" error : "<<SD_twoSSB<<'\n'<<'\n'<<
0192
0193 '\t'<<"DSB : "<<mean_DSB<<'\t'<<" error : "<<SD_DSB<<'\n'<<
0194 '\t'<<"DSBp : "<<mean_DSBp<<'\t'<<" error : "<<SD_DSBp<<'\n'<<
0195 '\t'<<"DSBpp : "<<mean_DSBpp<<'\t'<<" error : "<<SD_DSBpp<<'\n';
0196
0197
0198 f->Close();
0199 const Int_t n1 = 7;
0200 Double_t x1[n1] = {2,4,6,8,10,12,14};
0201 Double_t _y1[n1] = {mean_SSBd,mean_SSBi,mean_SSBm,mean_DSBd,mean_DSBi,
0202 mean_DSBm,mean_DSBh};
0203 Double_t err_y1[n1] = {SD_SSBd,SD_SSBi,SD_SSBm,SD_DSBd,SD_DSBi,SD_DSBm,SD_DSBh};
0204 TGraph* gr1 = new TGraphErrors(n1,x1,_y1,0,err_y1);
0205 gr1->SetTitle("Break Source Frequency");
0206 gr1->GetXaxis()->SetBinLabel(10,"SSB direct");
0207 gr1->GetXaxis()->SetBinLabel(20,"SSB indirect");
0208 gr1->GetXaxis()->SetBinLabel(35,"SSB mix");
0209 gr1->GetXaxis()->SetBinLabel(50,"DSB direct");
0210 gr1->GetXaxis()->SetBinLabel(65,"DSB indirect");
0211 gr1->GetXaxis()->SetBinLabel(80,"DSB mix");
0212 gr1->GetXaxis()->SetBinLabel(90,"DSB hybrid");
0213 gr1->SetFillColor(49);
0214
0215 const Int_t n2 = 6;
0216 Double_t x2[n2] = {2,4,6,8,10,12};
0217 Double_t _y2[n2] = {mean_SSB,mean_SSBp,mean_twoSSB,mean_DSB,mean_DSBp,
0218 mean_DSBpp};
0219 Double_t err_y2[n2] = {SD_SSB,SD_SSBp,SD_twoSSB,SD_DSB,SD_DSBp,SD_DSBpp};
0220 TGraph* gr2 = new TGraphErrors(n2,x2,_y2,0,err_y2);
0221 gr2->SetTitle("Break Complexity Frequency");
0222 gr2->GetXaxis()->SetBinLabel(10,"SSB");
0223 gr2->GetXaxis()->SetBinLabel(25,"SSBp");
0224 gr2->GetXaxis()->SetBinLabel(45,"2SSB");
0225 gr2->GetXaxis()->SetBinLabel(60,"DSB");
0226 gr2->GetXaxis()->SetBinLabel(75,"DSBp");
0227 gr2->GetXaxis()->SetBinLabel(90,"DSBpp");
0228 gr2->SetFillColor(49);
0229
0230
0231 c1->cd(1);
0232 gr1->Draw("ba");
0233
0234 c1->cd(2);
0235 gr2->Draw("ba");
0236 }
