Geant4 Cross Reference |
1 2 void BinLogX(TH1* h) 3 { 4 TAxis *axis = h->GetXaxis(); 5 int bins = axis->GetNbins(); 6 7 Axis_t from = axis->GetXmin(); 8 Axis_t to = axis->GetXmax(); 9 Axis_t width = (to - from) / bins; 10 Axis_t *new_bins = new Axis_t[bins + 1]; 11 12 for (int i = 0; i <= bins; i++) { 13 new_bins[i] = TMath::Power(10, from + i * width); 14 } 15 axis->Set(bins, new_bins); 16 delete []new_bins; 17 } 18 19 20 void plot(){ 21 22 double dens = 1000;//kg/m3 23 24 double R =50; 25 26 //TFile *fin = new TFile("AuNP_Livermore.root"); 27 TFile *fin = new TFile("AuNP.root"); 28 TH1F *h1neve = (TH1F*)fin->Get("h1Events"); 29 TH1F *h1Edep = (TH1F*)fin->Get("h1Edep"); 30 TH1F *h1secnp_cha = (TH1F*)fin->Get("h1SecEnergyNP_charged"); 31 TH1F *h1secnp_nut = (TH1F*)fin->Get("h1SecEnergyNP_nutral"); 32 TH1F *h1secnpsurf_cha= (TH1F*)fin->Get("h1SecEnergyNPSurf_charged"); 33 TH1F *h1secnpsurf_nut= (TH1F*)fin->Get("h1SecEnergyNPSurf_nutral"); 34 TH1F *h1sec_cha = (TH1F*)fin->Get("h1Sec_charged"); 35 TH1F *h1sec_nut = (TH1F*)fin->Get("h1Sec_nutral"); 36 TH1F *h1chem_0 = (TH1F*)fin->Get("h1Chem_0"); 37 TH1F *h1chem_1 = (TH1F*)fin->Get("h1Chem_1"); 38 TH1F *h1chem_2 = (TH1F*)fin->Get("h1Chem_2"); 39 TH1F *h1chem_3 = (TH1F*)fin->Get("h1Chem_3"); 40 TH1F *h1chem_4 = (TH1F*)fin->Get("h1Chem_4"); 41 TH1F *h1chem_5 = (TH1F*)fin->Get("h1Chem_5"); 42 TH1F *h1chem_6 = (TH1F*)fin->Get("h1Chem_6"); 43 TH1F *h1chem_7 = (TH1F*)fin->Get("h1Chem_7"); 44 45 TH2F *h2Edep = (TH2F*)fin->Get("h2Edep"); 46 TH2F *h2sec2_cha = (TH2F*)fin->Get("h2SecEnergyAbs_charged"); 47 TH2F *h2sec2_nut = (TH2F*)fin->Get("h2SecEnergyAbs_nutral"); 48 49 double neve = h1neve->GetBinContent(1); 50 h1Edep ->Scale(1./neve); 51 h1secnp_cha ->Scale(1./neve); 52 h1secnp_nut ->Scale(1./neve); 53 h1secnpsurf_cha->Scale(1./neve); 54 h1secnpsurf_nut->Scale(1./neve); 55 h1sec_cha ->Scale(1./neve); 56 h1sec_nut ->Scale(1./neve); 57 h1chem_0 ->Scale(1./neve); 58 h1chem_1 ->Scale(1./neve); 59 h1chem_2 ->Scale(1./neve); 60 h1chem_3 ->Scale(1./neve); 61 h1chem_4 ->Scale(1./neve); 62 h1chem_5 ->Scale(1./neve); 63 h1chem_6 ->Scale(1./neve); 64 h1chem_7 ->Scale(1./neve); 65 h2Edep ->Scale(1./neve); 66 h2sec2_cha ->Scale(1./neve); 67 h2sec2_nut ->Scale(1./neve); 68 69 double val_cha=0; 70 double err_cha=0; 71 double val_nut=0; 72 double err_nut=0; 73 74 int NR = h1sec_cha -> GetNbinsX(); 75 TH1D *h1sec_tot = new TH1D("h1Sec_tot","h1Sec_tot",NR,1,5); 76 BinLogX(h1sec_tot); 77 for(int i=0;i<NR;i++){ 78 val_cha = h1sec_cha->GetBinContent(i+1); 79 err_cha = h1sec_cha->GetBinError (i+1); 80 val_nut = h1sec_nut->GetBinContent(i+1); 81 err_nut = h1sec_nut->GetBinError (i+1); 82 //val_cha = h1sec_cha->GetBinContent(i+1)/h1sec_cha->GetBinWidth(i+1); 83 //err_cha = h1sec_cha->GetBinError (i+1)/h1sec_cha->GetBinWidth(i+1); 84 //val_nut = h1sec_nut->GetBinContent(i+1)/h1sec_nut->GetBinWidth(i+1); 85 //err_nut = h1sec_nut->GetBinError (i+1)/h1sec_nut->GetBinWidth(i+1); 86 h1sec_cha -> SetBinContent(i+1,val_cha); 87 h1sec_cha -> SetBinError (i+1,err_cha); 88 h1sec_nut -> SetBinContent(i+1,val_nut); 89 h1sec_nut -> SetBinError (i+1,err_nut); 90 h1sec_tot -> SetBinContent(i+1,val_cha+val_nut); 91 h1sec_tot -> SetBinError (i+1,err_cha+err_nut); 92 } 93 94 TH1D *h1chem_tot = new TH1D("h1Chem_tot","h1Chem_tot",NR,1,5); 95 BinLogX(h1chem_tot); 96 for(int i=0;i<NR;i++){ 97 double val_0 = h1chem_0->GetBinContent(i+1)/h1chem_0->GetBinWidth(i+1); 98 double val_1 = h1chem_1->GetBinContent(i+1)/h1chem_1->GetBinWidth(i+1); 99 double val_2 = h1chem_2->GetBinContent(i+1)/h1chem_2->GetBinWidth(i+1); 100 double val_3 = h1chem_3->GetBinContent(i+1)/h1chem_3->GetBinWidth(i+1); 101 double val_4 = h1chem_4->GetBinContent(i+1)/h1chem_4->GetBinWidth(i+1); 102 double val_5 = h1chem_5->GetBinContent(i+1)/h1chem_5->GetBinWidth(i+1); 103 double val_6 = h1chem_6->GetBinContent(i+1)/h1chem_6->GetBinWidth(i+1); 104 double val_7 = h1chem_7->GetBinContent(i+1)/h1chem_7->GetBinWidth(i+1); 105 double err_0 = h1chem_0->GetBinError (i+1)/h1chem_0->GetBinWidth(i+1); 106 double err_1 = h1chem_1->GetBinError (i+1)/h1chem_1->GetBinWidth(i+1); 107 double err_2 = h1chem_2->GetBinError (i+1)/h1chem_2->GetBinWidth(i+1); 108 double err_3 = h1chem_3->GetBinError (i+1)/h1chem_3->GetBinWidth(i+1); 109 double err_4 = h1chem_4->GetBinError (i+1)/h1chem_4->GetBinWidth(i+1); 110 double err_5 = h1chem_5->GetBinError (i+1)/h1chem_5->GetBinWidth(i+1); 111 double err_6 = h1chem_6->GetBinError (i+1)/h1chem_6->GetBinWidth(i+1); 112 double err_7 = h1chem_7->GetBinError (i+1)/h1chem_7->GetBinWidth(i+1); 113 h1chem_tot -> SetBinContent(i+1,val_0+val_1+val_2+val_3+val_4+val_5+val_6+val_7); 114 h1chem_tot -> SetBinError (i+1,err_0+err_1+err_2+err_3+err_4+err_5+err_6+err_7); 115 } 116 117 118 NR = h1secnp_cha -> GetNbinsX(); 119 TH1D *h1secnp_tot = new TH1D("h1SecEnergyNP" ,"Energy Spectra in NP" ,NR,0,6); 120 BinLogX(h1secnp_tot); 121 NR = h1secnpsurf_cha -> GetNbinsX(); 122 TH1D *h1secnpsurf_tot = new TH1D("h1SecEnergyNPSurf","Energy Spectra at NP surface",NR,0,6); 123 BinLogX(h1secnpsurf_tot); 124 for(int i=0;i<NR;i++){ 125 val_cha = h1secnp_cha->GetBinContent(i+1); 126 err_cha = h1secnp_cha->GetBinError (i+1); 127 val_nut = h1secnp_nut->GetBinContent(i+1); 128 err_nut = h1secnp_nut->GetBinError (i+1); 129 //val_cha = h1secnp_cha->GetBinContent(i+1)/h1secnp_cha->GetBinWidth(i+1); 130 //err_cha = h1secnp_cha->GetBinError (i+1)/h1secnp_cha->GetBinWidth(i+1); 131 //val_nut = h1secnp_nut->GetBinContent(i+1)/h1secnp_nut->GetBinWidth(i+1); 132 //err_nut = h1secnp_nut->GetBinError (i+1)/h1secnp_nut->GetBinWidth(i+1); 133 h1secnp_cha -> SetBinContent(i+1,val_cha); 134 h1secnp_cha -> SetBinError (i+1,err_cha); 135 h1secnp_nut -> SetBinContent(i+1,val_nut); 136 h1secnp_nut -> SetBinError (i+1,err_nut); 137 h1secnp_tot -> SetBinContent(i+1,val_cha+val_nut); 138 h1secnp_tot -> SetBinError (i+1,err_cha+err_nut); 139 140 141 val_cha = h1secnpsurf_cha->GetBinContent(i+1); 142 err_cha = h1secnpsurf_cha->GetBinError (i+1); 143 val_nut = h1secnpsurf_nut->GetBinContent(i+1); 144 err_nut = h1secnpsurf_nut->GetBinError (i+1); 145 //val_cha = h1secnpsurf_cha->GetBinContent(i+1)/h1secnpsurf_cha->GetBinWidth(i+1); 146 //err_cha = h1secnpsurf_cha->GetBinError (i+1)/h1secnpsurf_cha->GetBinWidth(i+1); 147 //val_nut = h1secnpsurf_nut->GetBinContent(i+1)/h1secnpsurf_nut->GetBinWidth(i+1); 148 //err_nut = h1secnpsurf_nut->GetBinError (i+1)/h1secnpsurf_nut->GetBinWidth(i+1); 149 h1secnpsurf_cha -> SetBinContent(i+1,val_cha); 150 h1secnpsurf_cha -> SetBinError (i+1,err_cha); 151 h1secnpsurf_nut -> SetBinContent(i+1,val_nut); 152 h1secnpsurf_nut -> SetBinError (i+1,err_nut); 153 h1secnpsurf_tot -> SetBinContent(i+1,val_cha+val_nut); 154 h1secnpsurf_tot -> SetBinError (i+1,err_cha+err_nut); 155 156 } 157 158 NR = h1Edep -> GetNbinsX(); 159 for(int i=0;i<NR;i++){ 160 double rmax = h1Edep->GetXaxis()->GetBinUpEdge (i+1)*1.E-9;//m 161 double rmin = h1Edep->GetXaxis()->GetBinLowEdge(i+1)*1.E-9;//m 162 double vol = 4./3.*TMath::Pi()*(pow(rmax,3)-pow(rmin,3)); 163 double mass = dens*vol; 164 h1Edep -> SetBinContent(i+1,h1Edep->GetBinContent(i+1)/mass); 165 h1Edep -> SetBinError (i+1,h1Edep->GetBinError (i+1)/mass); 166 } 167 168 169 int Nazm =h2Edep->GetXaxis()->GetNbins()-1; 170 NR =h2Edep->GetYaxis()->GetNbins()-1; 171 double Rmax =h2Edep->GetYaxis()->GetBinUpEdge(NR); 172 TH2D *h2Edep_pol = new TH2D("h2Edep_pol","h2Edep_pol",Nazm,0,360,NR,0,Rmax); 173 for(int i=0;i<Nazm;i++){ 174 for(int j=0;j<NR;j++){ 175 double height = 10*1.E-9;//m 176 double rmax = h2Edep_pol->GetYaxis()->GetBinUpEdge (j+1)*1.E-9;//m 177 double rmin = h2Edep_pol->GetYaxis()->GetBinLowEdge (j+1)*1.E-9;//m 178 double vol = (TMath::Pi()*pow(rmax,2)-TMath::Pi()*pow(rmin,2))*height; 179 double mass = dens*vol; 180 h2Edep_pol->SetBinContent(i+1,j+1,h2Edep->GetBinContent(i+1,j+1)/mass); 181 } 182 } 183 int NRX = h2sec2_cha->GetXaxis()->GetNbins(); 184 int NRY = h2sec2_cha->GetYaxis()->GetNbins(); 185 TH2D *h2sec2_tot = (TH2D*)h2sec2_cha->Clone("h2sec2_tot"); 186 for(int i=0;i<NRX;i++){ 187 for(int j=0;j<NRY;j++){ 188 val_cha = h2sec2_cha->GetBinContent(i+1,j+1); 189 val_nut = h2sec2_nut->GetBinContent(i+1,j+1); 190 h2sec2_tot->SetBinContent(i+1,j+1,val_cha+val_nut); 191 } 192 } 193 194 195 196 197 198 199 200 TCanvas *cedep = new TCanvas("cedep","cedep",700,700); 201 cedep->cd(1)->SetTheta(90); 202 cedep->cd(1)->SetPhi(0); 203 cedep->cd(1)->SetLogz(); 204 h2Edep_pol->RebinX(6); 205 h2Edep_pol->RebinY(10); 206 //h2Edep_pol->SetAxisRange(0,1000,"Y"); 207 h2Edep_pol->SetAxisRange(0,150,"Y"); 208 h2Edep_pol->SetTitle(""); 209 h2Edep_pol->SetXTitle("Angle [deg]"); 210 h2Edep_pol->Draw("surf1POL"); 211 212 TCanvas *cedepsec = new TCanvas("cedepsec","cedepsec",1800,700); 213 cedepsec->Divide(2,1); 214 cedepsec->cd(1)->SetLogy(); 215 cedepsec->cd(1)->SetLogx(); 216 h1Edep->SetYTitle("Energy Deposite [Gy/event]"); 217 h1Edep->SetXTitle("Distance from NP center[nm]"); 218 h1Edep->SetMaximum(100); 219 h1Edep->SetMinimum(0.0000000001); 220 h1Edep->SetAxisRange(0,1000000); 221 h1Edep->Draw(); 222 cedepsec->cd(2)->SetLogy(); 223 cedepsec->cd(2)->SetLogx(); 224 h1sec_tot->SetLineColor(kBlack); 225 h1sec_cha->SetLineColor(kRed); 226 h1sec_nut->SetLineColor(kBlue); 227 h1sec_tot->SetMarkerColor(kBlack); 228 h1sec_cha->SetMarkerColor(kRed); 229 h1sec_nut->SetMarkerColor(kBlue); 230 h1sec_tot->SetAxisRange(0,1000000); 231 h1sec_tot->SetMaximum(10); 232 h1sec_tot->SetMinimum(0.000001); 233 h1sec_tot->SetTitle ("Number of Generated Secondaries"); 234 h1sec_tot->SetYTitle("N/nm/event"); 235 h1sec_tot->SetXTitle("Distance from NP center[nm]"); 236 //h1sec_cha->Draw("sameL"); 237 //h1sec_nut->Draw("sameL"); 238 h1sec_tot->Draw(); 239 //cedepsec->cd(3)->SetLogy(); 240 //cedepsec->cd(3)->SetLogx(); 241 //h1chem_tot->SetLineColor (kBlack); 242 //h1chem_tot->SetMarkerColor(kBlack); 243 //h1chem_tot->SetAxisRange(0,1000000); 244 //h1chem_tot->SetMaximum(10); 245 //h1chem_tot->SetMinimum(0.000001); 246 //h1chem_tot->SetTitle ("Number of Generated Chemical Species"); 247 //h1chem_tot->SetYTitle("N/nm/event"); 248 //h1chem_tot->SetXTitle("Distance from NP center[nm]"); 249 //h1chem_tot->Draw(); 250 251 252 TCanvas *csec = new TCanvas("csec","csec",2400,700); 253 csec->Divide(3,1); 254 csec->cd(1)->SetLogx(); 255 csec->cd(1)->SetLogy(); 256 h1secnp_tot->SetLineColor(kBlack); 257 h1secnp_cha->SetLineColor(kRed); 258 h1secnp_nut->SetLineColor(kBlue); 259 h1secnp_tot->SetMarkerColor(kBlack); 260 h1secnp_cha->SetMarkerColor(kRed); 261 h1secnp_nut->SetMarkerColor(kBlue); 262 h1secnp_tot->RebinX(10); 263 h1secnp_tot->Scale(1./10.); 264 h1secnp_tot->SetMaximum(20); 265 h1secnp_tot->SetMinimum(0.000000001); 266 h1secnp_tot->SetAxisRange(1,1.E6); 267 h1secnp_tot->SetTitle ("Secondary Energy produced in NP"); 268 h1secnp_tot->SetYTitle("N/event"); 269 h1secnp_tot->SetXTitle("Secandary Energy [eV]"); 270 //h1secnp_cha->DrawCopy(); 271 //h1secnp_nut->DrawCopy("same"); 272 h1secnp_tot->DrawCopy(); 273 csec->cd(2)->SetLogx(); 274 csec->cd(2)->SetLogy(); 275 h1secnpsurf_tot->SetLineColor(kBlack); 276 h1secnpsurf_cha->SetLineColor(kRed); 277 h1secnpsurf_nut->SetLineColor(kBlue); 278 h1secnpsurf_tot->SetMarkerColor(kBlack); 279 h1secnpsurf_cha->SetMarkerColor(kRed); 280 h1secnpsurf_nut->SetMarkerColor(kBlue); 281 h1secnpsurf_tot->RebinX(10); 282 h1secnpsurf_tot->Scale(1./10.); 283 h1secnpsurf_tot->SetMaximum(20); 284 h1secnpsurf_tot->SetMinimum(0.000000001); 285 h1secnpsurf_tot->SetAxisRange(1,1.E6); 286 h1secnpsurf_tot->SetTitle ("Secondary Energy at NP surface"); 287 h1secnpsurf_tot->SetYTitle("N/event"); 288 h1secnpsurf_tot->SetXTitle("Secandary Energy [eV]"); 289 //h1secnpsurf_cha->DrawCopy(); 290 //h1secnpsurf_nut->DrawCopy("same"); 291 h1secnpsurf_tot->DrawCopy(); 292 csec->cd(3)->SetLogy(); 293 h2sec2_tot->SetTitle ("Secondary Energy vs production point"); 294 h2sec2_tot->SetYTitle("Secondary Energy [eV]"); 295 h2sec2_tot->SetXTitle("Distance from NP center [nm]"); 296 h2sec2_tot->GetXaxis()->SetRange(0,1000); 297 h2sec2_tot->Draw("colz"); 298 299 } 300 301