Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-02-23 09:21:47

0001 
0002 void BinLogX(TH1* h)
0003 {
0004                 TAxis *axis = h->GetXaxis();
0005                 int bins = axis->GetNbins();
0006 
0007                 Axis_t from      = axis->GetXmin();
0008                 Axis_t to        = axis->GetXmax();
0009                 Axis_t width     = (to - from) / bins;
0010                 Axis_t *new_bins = new Axis_t[bins + 1];
0011 
0012                 for (int i = 0; i <= bins; i++) {
0013                                 new_bins[i] = TMath::Power(10, from + i * width);
0014                 }
0015                 axis->Set(bins, new_bins);
0016                 delete []new_bins;
0017 }
0018 
0019 
0020 void plot(){
0021 
0022                 double dens = 1000;//kg/m3
0023 
0024                 double R =50;
0025 
0026                 //TFile *fin           = new TFile("AuNP_Livermore.root");
0027                 TFile *fin           = new TFile("AuNP.root");
0028                 TH1F *h1neve         = (TH1F*)fin->Get("h1Events");
0029                 TH1F *h1Edep         = (TH1F*)fin->Get("h1Edep");
0030                 TH1F *h1secnp_cha    = (TH1F*)fin->Get("h1SecEnergyNP_charged");
0031                 TH1F *h1secnp_nut    = (TH1F*)fin->Get("h1SecEnergyNP_nutral");
0032                 TH1F *h1secnpsurf_cha= (TH1F*)fin->Get("h1SecEnergyNPSurf_charged");
0033                 TH1F *h1secnpsurf_nut= (TH1F*)fin->Get("h1SecEnergyNPSurf_nutral");
0034                 TH1F *h1sec_cha      = (TH1F*)fin->Get("h1Sec_charged");
0035                 TH1F *h1sec_nut      = (TH1F*)fin->Get("h1Sec_nutral");
0036                 TH1F *h1chem_0       = (TH1F*)fin->Get("h1Chem_0");
0037                 TH1F *h1chem_1       = (TH1F*)fin->Get("h1Chem_1");
0038                 TH1F *h1chem_2       = (TH1F*)fin->Get("h1Chem_2");
0039                 TH1F *h1chem_3       = (TH1F*)fin->Get("h1Chem_3");
0040                 TH1F *h1chem_4       = (TH1F*)fin->Get("h1Chem_4");
0041                 TH1F *h1chem_5       = (TH1F*)fin->Get("h1Chem_5");
0042                 TH1F *h1chem_6       = (TH1F*)fin->Get("h1Chem_6");
0043                 TH1F *h1chem_7       = (TH1F*)fin->Get("h1Chem_7");
0044 
0045                 TH2F *h2Edep         = (TH2F*)fin->Get("h2Edep");
0046                 TH2F *h2sec2_cha     = (TH2F*)fin->Get("h2SecEnergyAbs_charged");
0047                 TH2F *h2sec2_nut     = (TH2F*)fin->Get("h2SecEnergyAbs_nutral");
0048 
0049                 double neve = h1neve->GetBinContent(1);
0050                 h1Edep         ->Scale(1./neve);
0051                 h1secnp_cha    ->Scale(1./neve);
0052                 h1secnp_nut    ->Scale(1./neve);
0053                 h1secnpsurf_cha->Scale(1./neve);
0054                 h1secnpsurf_nut->Scale(1./neve);
0055                 h1sec_cha      ->Scale(1./neve);
0056                 h1sec_nut      ->Scale(1./neve);
0057                 h1chem_0       ->Scale(1./neve);
0058                 h1chem_1       ->Scale(1./neve);
0059                 h1chem_2       ->Scale(1./neve);
0060                 h1chem_3       ->Scale(1./neve);
0061                 h1chem_4       ->Scale(1./neve);
0062                 h1chem_5       ->Scale(1./neve);
0063                 h1chem_6       ->Scale(1./neve);
0064                 h1chem_7       ->Scale(1./neve);
0065                 h2Edep         ->Scale(1./neve);
0066                 h2sec2_cha     ->Scale(1./neve);
0067                 h2sec2_nut     ->Scale(1./neve);
0068 
0069                 double val_cha=0;
0070                 double err_cha=0;
0071                 double val_nut=0;
0072                 double err_nut=0;
0073 
0074                 int NR = h1sec_cha  -> GetNbinsX();
0075                 TH1D *h1sec_tot = new TH1D("h1Sec_tot","h1Sec_tot",NR,1,5);
0076         BinLogX(h1sec_tot);
0077                 for(int i=0;i<NR;i++){
0078                                 val_cha = h1sec_cha->GetBinContent(i+1);
0079                                 err_cha = h1sec_cha->GetBinError  (i+1);
0080                                 val_nut = h1sec_nut->GetBinContent(i+1);
0081                                 err_nut = h1sec_nut->GetBinError  (i+1);
0082                                 //val_cha = h1sec_cha->GetBinContent(i+1)/h1sec_cha->GetBinWidth(i+1);
0083                                 //err_cha = h1sec_cha->GetBinError  (i+1)/h1sec_cha->GetBinWidth(i+1);
0084                                 //val_nut = h1sec_nut->GetBinContent(i+1)/h1sec_nut->GetBinWidth(i+1);
0085                                 //err_nut = h1sec_nut->GetBinError  (i+1)/h1sec_nut->GetBinWidth(i+1);
0086                                 h1sec_cha  -> SetBinContent(i+1,val_cha);
0087                                 h1sec_cha  -> SetBinError  (i+1,err_cha);
0088                                 h1sec_nut  -> SetBinContent(i+1,val_nut);
0089                                 h1sec_nut  -> SetBinError  (i+1,err_nut);
0090                                 h1sec_tot  -> SetBinContent(i+1,val_cha+val_nut);
0091                                 h1sec_tot  -> SetBinError  (i+1,err_cha+err_nut);
0092                 }
0093 
0094                 TH1D *h1chem_tot = new TH1D("h1Chem_tot","h1Chem_tot",NR,1,5);
0095         BinLogX(h1chem_tot);
0096                 for(int i=0;i<NR;i++){
0097                                 double val_0 = h1chem_0->GetBinContent(i+1)/h1chem_0->GetBinWidth(i+1);
0098                                 double val_1 = h1chem_1->GetBinContent(i+1)/h1chem_1->GetBinWidth(i+1);
0099                                 double val_2 = h1chem_2->GetBinContent(i+1)/h1chem_2->GetBinWidth(i+1);
0100                                 double val_3 = h1chem_3->GetBinContent(i+1)/h1chem_3->GetBinWidth(i+1);
0101                                 double val_4 = h1chem_4->GetBinContent(i+1)/h1chem_4->GetBinWidth(i+1);
0102                                 double val_5 = h1chem_5->GetBinContent(i+1)/h1chem_5->GetBinWidth(i+1);
0103                                 double val_6 = h1chem_6->GetBinContent(i+1)/h1chem_6->GetBinWidth(i+1);
0104                                 double val_7 = h1chem_7->GetBinContent(i+1)/h1chem_7->GetBinWidth(i+1);
0105                                 double err_0 = h1chem_0->GetBinError  (i+1)/h1chem_0->GetBinWidth(i+1);
0106                                 double err_1 = h1chem_1->GetBinError  (i+1)/h1chem_1->GetBinWidth(i+1);
0107                                 double err_2 = h1chem_2->GetBinError  (i+1)/h1chem_2->GetBinWidth(i+1);
0108                                 double err_3 = h1chem_3->GetBinError  (i+1)/h1chem_3->GetBinWidth(i+1);
0109                                 double err_4 = h1chem_4->GetBinError  (i+1)/h1chem_4->GetBinWidth(i+1);
0110                                 double err_5 = h1chem_5->GetBinError  (i+1)/h1chem_5->GetBinWidth(i+1);
0111                                 double err_6 = h1chem_6->GetBinError  (i+1)/h1chem_6->GetBinWidth(i+1);
0112                                 double err_7 = h1chem_7->GetBinError  (i+1)/h1chem_7->GetBinWidth(i+1);
0113                                 h1chem_tot  -> SetBinContent(i+1,val_0+val_1+val_2+val_3+val_4+val_5+val_6+val_7);
0114                                 h1chem_tot  -> SetBinError  (i+1,err_0+err_1+err_2+err_3+err_4+err_5+err_6+err_7);
0115                 }
0116 
0117 
0118                 NR = h1secnp_cha -> GetNbinsX();
0119                 TH1D *h1secnp_tot     = new TH1D("h1SecEnergyNP"    ,"Energy Spectra in NP"        ,NR,0,6);
0120         BinLogX(h1secnp_tot);
0121                 NR = h1secnpsurf_cha -> GetNbinsX();
0122                 TH1D *h1secnpsurf_tot = new TH1D("h1SecEnergyNPSurf","Energy Spectra at NP surface",NR,0,6);
0123         BinLogX(h1secnpsurf_tot);
0124                 for(int i=0;i<NR;i++){
0125                                 val_cha = h1secnp_cha->GetBinContent(i+1);
0126                                 err_cha = h1secnp_cha->GetBinError  (i+1);
0127                                 val_nut = h1secnp_nut->GetBinContent(i+1);
0128                                 err_nut = h1secnp_nut->GetBinError  (i+1);
0129                                 //val_cha = h1secnp_cha->GetBinContent(i+1)/h1secnp_cha->GetBinWidth(i+1);
0130                                 //err_cha = h1secnp_cha->GetBinError  (i+1)/h1secnp_cha->GetBinWidth(i+1);
0131                                 //val_nut = h1secnp_nut->GetBinContent(i+1)/h1secnp_nut->GetBinWidth(i+1);
0132                                 //err_nut = h1secnp_nut->GetBinError  (i+1)/h1secnp_nut->GetBinWidth(i+1);
0133                                 h1secnp_cha -> SetBinContent(i+1,val_cha);
0134                                 h1secnp_cha -> SetBinError  (i+1,err_cha);
0135                                 h1secnp_nut -> SetBinContent(i+1,val_nut);
0136                                 h1secnp_nut -> SetBinError  (i+1,err_nut);
0137                                 h1secnp_tot -> SetBinContent(i+1,val_cha+val_nut);
0138                                 h1secnp_tot -> SetBinError  (i+1,err_cha+err_nut);
0139 
0140 
0141                                 val_cha = h1secnpsurf_cha->GetBinContent(i+1);
0142                                 err_cha = h1secnpsurf_cha->GetBinError  (i+1);
0143                                 val_nut = h1secnpsurf_nut->GetBinContent(i+1);
0144                                 err_nut = h1secnpsurf_nut->GetBinError  (i+1);
0145                                 //val_cha = h1secnpsurf_cha->GetBinContent(i+1)/h1secnpsurf_cha->GetBinWidth(i+1);
0146                                 //err_cha = h1secnpsurf_cha->GetBinError  (i+1)/h1secnpsurf_cha->GetBinWidth(i+1);
0147                                 //val_nut = h1secnpsurf_nut->GetBinContent(i+1)/h1secnpsurf_nut->GetBinWidth(i+1);
0148                                 //err_nut = h1secnpsurf_nut->GetBinError  (i+1)/h1secnpsurf_nut->GetBinWidth(i+1);
0149                                 h1secnpsurf_cha -> SetBinContent(i+1,val_cha);
0150                                 h1secnpsurf_cha -> SetBinError  (i+1,err_cha);
0151                                 h1secnpsurf_nut -> SetBinContent(i+1,val_nut);
0152                                 h1secnpsurf_nut -> SetBinError  (i+1,err_nut);
0153                                 h1secnpsurf_tot -> SetBinContent(i+1,val_cha+val_nut);
0154                                 h1secnpsurf_tot -> SetBinError  (i+1,err_cha+err_nut);
0155 
0156                 }
0157 
0158                 NR = h1Edep  -> GetNbinsX();
0159                 for(int i=0;i<NR;i++){
0160                                 double rmax = h1Edep->GetXaxis()->GetBinUpEdge (i+1)*1.E-9;//m
0161                                 double rmin = h1Edep->GetXaxis()->GetBinLowEdge(i+1)*1.E-9;//m
0162                                 double vol  = 4./3.*TMath::Pi()*(pow(rmax,3)-pow(rmin,3));
0163                                 double mass = dens*vol;
0164                                 h1Edep  -> SetBinContent(i+1,h1Edep->GetBinContent(i+1)/mass);
0165                                 h1Edep  -> SetBinError  (i+1,h1Edep->GetBinError  (i+1)/mass);
0166                 }
0167 
0168 
0169                 int Nazm =h2Edep->GetXaxis()->GetNbins()-1;
0170                 NR   =h2Edep->GetYaxis()->GetNbins()-1;
0171                 double    Rmax =h2Edep->GetYaxis()->GetBinUpEdge(NR);
0172                 TH2D *h2Edep_pol = new TH2D("h2Edep_pol","h2Edep_pol",Nazm,0,360,NR,0,Rmax);
0173                 for(int i=0;i<Nazm;i++){
0174                                 for(int j=0;j<NR;j++){
0175                                                 double height = 10*1.E-9;//m
0176                                                 double rmax = h2Edep_pol->GetYaxis()->GetBinUpEdge  (j+1)*1.E-9;//m
0177                                                 double rmin = h2Edep_pol->GetYaxis()->GetBinLowEdge (j+1)*1.E-9;//m
0178                                                 double vol  = (TMath::Pi()*pow(rmax,2)-TMath::Pi()*pow(rmin,2))*height;
0179                                         double mass = dens*vol;
0180                                                 h2Edep_pol->SetBinContent(i+1,j+1,h2Edep->GetBinContent(i+1,j+1)/mass);
0181                                 }
0182                 }
0183                 int NRX = h2sec2_cha->GetXaxis()->GetNbins();
0184                 int NRY = h2sec2_cha->GetYaxis()->GetNbins();
0185                 TH2D *h2sec2_tot = (TH2D*)h2sec2_cha->Clone("h2sec2_tot");
0186                 for(int i=0;i<NRX;i++){
0187                                 for(int j=0;j<NRY;j++){
0188                                                 val_cha = h2sec2_cha->GetBinContent(i+1,j+1);
0189                                                 val_nut = h2sec2_nut->GetBinContent(i+1,j+1);
0190                                                 h2sec2_tot->SetBinContent(i+1,j+1,val_cha+val_nut);
0191                                 }
0192                 }
0193 
0194 
0195 
0196 
0197 
0198 
0199 
0200                 TCanvas *cedep = new TCanvas("cedep","cedep",700,700);
0201                 cedep->cd(1)->SetTheta(90);
0202                 cedep->cd(1)->SetPhi(0);
0203                 cedep->cd(1)->SetLogz();
0204                 h2Edep_pol->RebinX(6);
0205                 h2Edep_pol->RebinY(10);
0206                 //h2Edep_pol->SetAxisRange(0,1000,"Y");
0207                 h2Edep_pol->SetAxisRange(0,150,"Y");
0208                 h2Edep_pol->SetTitle("");
0209                 h2Edep_pol->SetXTitle("Angle [deg]");
0210                 h2Edep_pol->Draw("surf1POL");
0211 
0212                 TCanvas *cedepsec = new TCanvas("cedepsec","cedepsec",1800,700);
0213                 cedepsec->Divide(2,1);
0214                 cedepsec->cd(1)->SetLogy();
0215                 cedepsec->cd(1)->SetLogx();
0216                 h1Edep->SetYTitle("Energy Deposite [Gy/event]");
0217                 h1Edep->SetXTitle("Distance from NP center[nm]");
0218                 h1Edep->SetMaximum(100);
0219                 h1Edep->SetMinimum(0.0000000001);
0220                 h1Edep->SetAxisRange(0,1000000);
0221                 h1Edep->Draw();
0222                 cedepsec->cd(2)->SetLogy();
0223                 cedepsec->cd(2)->SetLogx();
0224                 h1sec_tot->SetLineColor(kBlack);
0225                 h1sec_cha->SetLineColor(kRed);
0226                 h1sec_nut->SetLineColor(kBlue);
0227                 h1sec_tot->SetMarkerColor(kBlack);
0228                 h1sec_cha->SetMarkerColor(kRed);
0229                 h1sec_nut->SetMarkerColor(kBlue);
0230                 h1sec_tot->SetAxisRange(0,1000000);
0231                 h1sec_tot->SetMaximum(10);
0232                 h1sec_tot->SetMinimum(0.000001);
0233                 h1sec_tot->SetTitle ("Number of Generated Secondaries");
0234                 h1sec_tot->SetYTitle("N/nm/event");
0235                 h1sec_tot->SetXTitle("Distance from NP center[nm]");
0236                 //h1sec_cha->Draw("sameL");
0237                 //h1sec_nut->Draw("sameL");
0238                 h1sec_tot->Draw();
0239                 //cedepsec->cd(3)->SetLogy();
0240                 //cedepsec->cd(3)->SetLogx();
0241                 //h1chem_tot->SetLineColor  (kBlack);
0242                 //h1chem_tot->SetMarkerColor(kBlack);
0243                 //h1chem_tot->SetAxisRange(0,1000000);
0244                 //h1chem_tot->SetMaximum(10);
0245                 //h1chem_tot->SetMinimum(0.000001);
0246                 //h1chem_tot->SetTitle ("Number of Generated Chemical Species");
0247                 //h1chem_tot->SetYTitle("N/nm/event");
0248                 //h1chem_tot->SetXTitle("Distance from NP center[nm]");
0249                 //h1chem_tot->Draw();
0250 
0251 
0252                 TCanvas *csec = new TCanvas("csec","csec",2400,700);
0253                 csec->Divide(3,1);
0254                 csec->cd(1)->SetLogx();
0255                 csec->cd(1)->SetLogy();
0256                 h1secnp_tot->SetLineColor(kBlack);
0257                 h1secnp_cha->SetLineColor(kRed);
0258                 h1secnp_nut->SetLineColor(kBlue);
0259                 h1secnp_tot->SetMarkerColor(kBlack);
0260                 h1secnp_cha->SetMarkerColor(kRed);
0261                 h1secnp_nut->SetMarkerColor(kBlue);
0262                 h1secnp_tot->RebinX(10);
0263                 h1secnp_tot->Scale(1./10.);
0264                 h1secnp_tot->SetMaximum(20);
0265                 h1secnp_tot->SetMinimum(0.000000001);
0266                 h1secnp_tot->SetAxisRange(1,1.E6);
0267                 h1secnp_tot->SetTitle ("Secondary Energy produced in NP");
0268                 h1secnp_tot->SetYTitle("N/event");
0269                 h1secnp_tot->SetXTitle("Secandary Energy [eV]");
0270                 //h1secnp_cha->DrawCopy();
0271                 //h1secnp_nut->DrawCopy("same");
0272                 h1secnp_tot->DrawCopy();
0273                 csec->cd(2)->SetLogx();
0274                 csec->cd(2)->SetLogy();
0275                 h1secnpsurf_tot->SetLineColor(kBlack);
0276                 h1secnpsurf_cha->SetLineColor(kRed);
0277                 h1secnpsurf_nut->SetLineColor(kBlue);
0278                 h1secnpsurf_tot->SetMarkerColor(kBlack);
0279                 h1secnpsurf_cha->SetMarkerColor(kRed);
0280                 h1secnpsurf_nut->SetMarkerColor(kBlue);
0281                 h1secnpsurf_tot->RebinX(10);
0282                 h1secnpsurf_tot->Scale(1./10.);
0283                 h1secnpsurf_tot->SetMaximum(20);
0284                 h1secnpsurf_tot->SetMinimum(0.000000001);
0285                 h1secnpsurf_tot->SetAxisRange(1,1.E6);
0286                 h1secnpsurf_tot->SetTitle ("Secondary Energy at NP surface");
0287                 h1secnpsurf_tot->SetYTitle("N/event");
0288                 h1secnpsurf_tot->SetXTitle("Secandary Energy [eV]");
0289                 //h1secnpsurf_cha->DrawCopy();
0290                 //h1secnpsurf_nut->DrawCopy("same");
0291                 h1secnpsurf_tot->DrawCopy();
0292                 csec->cd(3)->SetLogy();
0293                 h2sec2_tot->SetTitle ("Secondary Energy vs production point");
0294                 h2sec2_tot->SetYTitle("Secondary Energy [eV]");
0295                 h2sec2_tot->SetXTitle("Distance from NP center [nm]");
0296                 h2sec2_tot->GetXaxis()->SetRange(0,1000);
0297                 h2sec2_tot->Draw("colz");
0298 
0299 }
0300