File indexing completed on 2025-10-13 08:25:09
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 auto 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 if (number<2) {
0069 std::cout << "Not enough entries in the \"primary_source\" TTree (" << (long)number << " entries)\n";
0070 gApplication->Terminate(0);
0071 }
0072
0073 tree = (TTree *) f->Get("tuples/source");
0074 tree->SetBranchAddress("Primary",&Primary);
0075 tree->SetBranchAddress("Energy",&Energy);
0076 tree->SetBranchAddress("SSBd",&SSBd);
0077 tree->SetBranchAddress("SSBi",&SSBi);
0078 tree->SetBranchAddress("SSBm",&SSBm);
0079 tree->SetBranchAddress("DSBd",&DSBd);
0080 tree->SetBranchAddress("DSBi",&DSBi);
0081 tree->SetBranchAddress("DSBm",&DSBm);
0082 tree->SetBranchAddress("DSBh",&DSBh);
0083
0084 Long64_t nentries = tree->GetEntries();
0085
0086 for(int i = 0;i<nentries;i++){
0087 tree->GetEntry(i);
0088 total_SSBd += SSBd;
0089 total_SSBd2 += SSBd *SSBd;
0090 total_SSBi += SSBi;
0091 total_SSBi2 += SSBi *SSBi;
0092 total_SSBm += SSBm;
0093 total_SSBm2 += SSBm *SSBm;
0094
0095 total_DSBd += DSBd;
0096 total_DSBd2 += DSBd *DSBd;
0097 total_DSBi += DSBi;
0098 total_DSBi2 += DSBi *DSBi;
0099 total_DSBm += DSBm;
0100 total_DSBm2 += DSBm *DSBm;
0101 total_DSBh += DSBh;
0102 total_DSBh2 += DSBh *DSBh;
0103 }
0104
0105 tree = (TTree *) f->Get("tuples/damage");
0106 tree->SetBranchAddress("EnergyDeposited_eV",&EnergyDeposited_eV);
0107 nentries = tree->GetEntries();
0108 for(int i = 0;i<nentries;i++){
0109 tree->GetEntry(i);
0110 acc_edep += EnergyDeposited_eV;
0111 acc_edep2 += EnergyDeposited_eV *EnergyDeposited_eV;
0112 }
0113
0114 tree = (TTree *) f->Get("tuples/classification");
0115 tree->SetBranchAddress("SSB",&SSB);
0116 tree->SetBranchAddress("SSBp",&SSBp);
0117 tree->SetBranchAddress("2SSB",&twoSSB);
0118 tree->SetBranchAddress("DSB",&DSB);
0119 tree->SetBranchAddress("DSBp",&DSBp);
0120 tree->SetBranchAddress("DSBpp",&DSBpp);
0121
0122 nentries = tree->GetEntries();
0123
0124 for(int i = 0;i<nentries;i++){
0125 tree->GetEntry(i);
0126 total_SSB += SSB;
0127 total_SSB2 += SSB *SSB;
0128 total_SSBp += SSBp;
0129 total_SSBp2 += SSBp *SSBp;
0130 total_twoSSB += twoSSB;
0131 total_twoSSB2 += twoSSB *twoSSB;
0132
0133 total_DSB += DSB;
0134 total_DSB2 += DSB *DSB;
0135 total_DSBp += DSBp;
0136 total_DSBp2 += DSBp *DSBp;
0137 total_DSBpp += DSBpp;
0138 total_DSBpp2 += DSBpp *DSBpp;
0139 }
0140
0141 mean_SSBd = (Float_t) total_SSBd / number;
0142 mean_SSBi = (Float_t) total_SSBi / number;
0143 mean_SSBm = (Float_t) total_SSBm / number;
0144
0145 mean_DSBd = (Float_t) total_DSBd / number;
0146 mean_DSBi = (Float_t) total_DSBi / number;
0147 mean_DSBm = (Float_t) total_DSBm / number;
0148 mean_DSBh = (Float_t) total_DSBh / number;
0149
0150 Double_t SD_SSBd = sqrt(abs(((total_SSBd2 / number) - pow(total_SSBd / number,2)))/(number -1));
0151 Double_t SD_SSBi = sqrt(abs(((total_SSBi2 / number) - pow(total_SSBi / number,2)))/(number -1));
0152 Double_t SD_SSBm = sqrt(abs(((total_SSBm2 / number) - pow(total_SSBm / number,2)))/(number -1));
0153
0154 Double_t SD_DSBd = sqrt(abs(((total_DSBd2 / number) - pow(total_DSBd / number,2)))/(number -1));
0155 Double_t SD_DSBi = sqrt(abs(((total_DSBi2 / number) - pow(total_DSBi / number,2)))/(number -1));
0156 Double_t SD_DSBm = sqrt(abs(((total_DSBm2 / number) - pow(total_DSBm / number,2)))/(number -1));
0157 Double_t SD_DSBh = sqrt(abs(((total_DSBh2 / number) - pow(total_DSBh / number,2)))/(number -1));
0158
0159 mean_SSB = (Float_t) total_SSB / number;
0160 mean_SSBp = (Float_t) total_SSBp / number;
0161 mean_twoSSB = (Float_t) total_twoSSB / number;
0162
0163 mean_DSB = (Float_t) total_DSB / number;
0164 mean_DSBp = (Float_t) total_DSBp / number;
0165 mean_DSBpp = (Float_t) total_DSBpp / number;
0166
0167 Double_t SD_SSB = sqrt(abs(((total_SSB2 / number) - pow(total_SSB / number,2))/(number -1)));
0168 Double_t SD_SSBp = sqrt(abs(((total_SSBp2 / number) - pow(total_SSBp / number,2))
0169 /(number -1)));
0170 Double_t SD_twoSSB = sqrt(abs(((total_twoSSB2 / number) - pow(total_twoSSB /
0171 number,2))
0172 /(number -1)));
0173
0174 Double_t SD_DSB = sqrt(abs(((total_DSB2 / number) - pow(total_DSB / number,2))/(number -1)));
0175 Double_t SD_DSBp = sqrt(abs(((total_DSBp2 / number) - pow(total_DSBp / number,2))
0176 /(number -1)));
0177 Double_t SD_DSBpp = sqrt(abs(((total_DSBpp2 / number) - pow(total_DSBpp / number,2))
0178 /(number -1)));
0179
0180 cout<<"Paricle : "<<Primary<<'\t'
0181 <<"Energy [/MeV] : "<<Energy<<'\t'
0182 <<"number : "<<number<<'\n'
0183 <<" Output Damages : "<<'\n'<<
0184 '\t'<<"SSB direct : "<<mean_SSBd<<'\t'<<" error : "<<SD_SSBd<<'\n'<<
0185 '\t'<<"SSB indirect : "<<mean_SSBi<<'\t'<<" error : "<<SD_SSBi<<'\n'<<
0186 '\t'<<"SSB mix : "<<mean_SSBm<<'\t'<<" error : "<<SD_SSBm<<'\n'<<'\n'<<
0187
0188 '\t'<<"DSB direct : "<<mean_DSBd<<'\t'<<" error : "<<SD_DSBd<<'\n'<<
0189 '\t'<<"DSB indirect : "<<mean_DSBi<<'\t'<<" error : "<<SD_DSBi<<'\n'<<
0190 '\t'<<"DSB mix : "<<mean_DSBm<<'\t'<<" error : "<<SD_DSBm<<'\n'<<
0191 '\t'<<"DSB hybrid : "<<mean_DSBh<<'\t'<<" error : "<<SD_DSBh<<'\n'<<'\n'<<
0192
0193 '\t'<<"SSB : "<<mean_SSB<<'\t'<<" error : "<<SD_SSB<<'\n'<<
0194 '\t'<<"SSBp : "<<mean_SSBp<<'\t'<<" error : "<<SD_SSBp<<'\n'<<
0195 '\t'<<"2SSB : "<<mean_twoSSB<<'\t'<<" error : "<<SD_twoSSB<<'\n'<<'\n'<<
0196
0197 '\t'<<"DSB : "<<mean_DSB<<'\t'<<" error : "<<SD_DSB<<'\n'<<
0198 '\t'<<"DSBp : "<<mean_DSBp<<'\t'<<" error : "<<SD_DSBp<<'\n'<<
0199 '\t'<<"DSBpp : "<<mean_DSBpp<<'\t'<<" error : "<<SD_DSBpp<<'\n';
0200
0201
0202 f->Close();
0203 const Int_t n1 = 7;
0204 Double_t x1[n1] = {2,4,6,8,10,12,14};
0205 Double_t _y1[n1] = {mean_SSBd,mean_SSBi,mean_SSBm,mean_DSBd,mean_DSBi,
0206 mean_DSBm,mean_DSBh};
0207 Double_t err_y1[n1] = {SD_SSBd,SD_SSBi,SD_SSBm,SD_DSBd,SD_DSBi,SD_DSBm,SD_DSBh};
0208 TGraph* gr1 = new TGraphErrors(n1,x1,_y1,0,err_y1);
0209 gr1->SetTitle("Break Source Frequency");
0210 gr1->GetXaxis()->SetBinLabel(10,"SSB direct");
0211 gr1->GetXaxis()->SetBinLabel(20,"SSB indirect");
0212 gr1->GetXaxis()->SetBinLabel(35,"SSB mix");
0213 gr1->GetXaxis()->SetBinLabel(50,"DSB direct");
0214 gr1->GetXaxis()->SetBinLabel(65,"DSB indirect");
0215 gr1->GetXaxis()->SetBinLabel(80,"DSB mix");
0216 gr1->GetXaxis()->SetBinLabel(90,"DSB hybrid");
0217 gr1->SetFillColor(49);
0218
0219 const Int_t n2 = 6;
0220 Double_t x2[n2] = {2,4,6,8,10,12};
0221 Double_t _y2[n2] = {mean_SSB,mean_SSBp,mean_twoSSB,mean_DSB,mean_DSBp,
0222 mean_DSBpp};
0223 Double_t err_y2[n2] = {SD_SSB,SD_SSBp,SD_twoSSB,SD_DSB,SD_DSBp,SD_DSBpp};
0224 TGraph* gr2 = new TGraphErrors(n2,x2,_y2,0,err_y2);
0225 gr2->SetTitle("Break Complexity Frequency");
0226 gr2->GetXaxis()->SetBinLabel(10,"SSB");
0227 gr2->GetXaxis()->SetBinLabel(25,"SSBp");
0228 gr2->GetXaxis()->SetBinLabel(45,"2SSB");
0229 gr2->GetXaxis()->SetBinLabel(60,"DSB");
0230 gr2->GetXaxis()->SetBinLabel(75,"DSBp");
0231 gr2->GetXaxis()->SetBinLabel(90,"DSBpp");
0232 gr2->SetFillColor(49);
0233
0234 c1->cd(1);
0235 gr1->Draw("ba");
0236
0237 c1->cd(2);
0238 gr2->Draw("ba");
0239 }
