//读取dedx文件
void readfile(string filename, Int_t A,TGraph *g)
{
ifstream in(filename);
string ss;
double a, b, e, dedx;
//g = new TGraph;
if(!in.is_open())
{
cout<<filename<<" does noet exist!"<<endl;
return 0;
}
getline(in,ss);
int i=0;
while(!in.eof())
{
in>>a>>b>>e>>dedx>>a>>b>>a>>b>>a>>b>>a>>b;//e:MeV/u,dedx:MeV/m;
if(!in.eof())
{
g->SetPoint(i,e*A,dedx*1000);//number,MeV,MeV/mm
}
i++;
}
in.close();
}
void energyvsdepth(TGraph *g, TGraph *depth2energy,TGraph *energy2depth)
{
Double_t er = 1000.;//MeV:1000 MeV 粒子;
Double_t dx = 0.001;//mm:穿透每个切片的厚度;
Double_t x=0.0;//mm:穿透总厚度
Int_t j =0;
while(er>0)
{
depth2energy->SetPoint(j,x,er);//写入穿透该切片前位置与dedx
energy2depth->SetPoint(j,er,x);
j++;//下一切片;
x=j*dx;//穿透该切片后位置;
er=er-dx*g->Eval(er);//穿透切片后剩余能量;
}
depth2energy->SetPoint(j,x,0);//写入粒子能量全部沉积后的位置与er;
energy2depth->SetPoint(j,0,x);//写入粒子能量全部沉积后的位置与er;
}
Double_t range(Double_t e)
{
TCanvas c_he4;
TGraph *g_he4,*d2e_he4,*e2d_he4;
g_he4 = new TGraph;
readfile("He4elossinSi.txt",4,g_he4);
d2e_he4= new TGraph;
e2d_he4= new TGraph;
energyvsdepth(g_he4,d2e_he4,e2d_he4);
return e2d_he4->Eval(0)-e2d_he4->Eval(e);
}
Double_t e[4]={5.486,10.,100.,200.};//MeV
for(int i=0;i<4;i++)
{
cout<<"the range of "<<e[i]<<"MeV 4He in Silicon is "<<range(e[i])<<" mm"<<endl;
}
以10 MeV,100 MeV,1000 MeV的4He粒子先后穿过厚度为1000 um,500 um,1000 um的能量损失;
Double_t eloss(Double_t e,Double_t thickness,TGraph *depth2energy,TGraph *energy2depth)
{
Double_t x_depth;
Double_t de;
x_depth=energy2depth->Eval(e);
x_depth=x_depth+thickness;
de=e-depth2energy->Eval(x_depth);
if(de<e)
{
de=de;
}
else
{
de=e;
}
return de;
}
TCanvas c_he4;
TGraph *g_he4,*d2e_he4,*e2d_he4;
g_he4 = new TGraph;
readfile("He4elossinSi.txt",4,g_he4);
d2e_he4= new TGraph;
e2d_he4= new TGraph;
energyvsdepth(g_he4,d2e_he4,e2d_he4);
Double_t e1[3]={10.,100.,200.};//MeV
Double_t thickness[3]={1.,0.5,1.};//mm
for(int i=0;i<3;i++)
{
for(int j=0;j<3;j++)
{
cout<<"the eloss of 4He with "<<e1[i]<<" MeV in "<<thickness[j]*1000<<" um Silicon is "<<eloss(e1[i],thickness[j],d2e_he4,e2d_he4)<<" MeV;"<<endl;
e1[i]=e1[i]-eloss(e1[i],thickness[j],d2e_he4,e2d_he4);
}
cout<<" "<<endl;
}
void PIDRes(string filename, Int_t A,TH2D *d1d2Res,TH2D*d2d3Res)
{
TGraph *g_dedx=new TGraph;
readfile(filename,A,g_dedx);
TGraph *d2e= new TGraph;
TGraph *e2d= new TGraph;
energyvsdepth(g_dedx,d2e,e2d);
TRandom3 *e_random =new TRandom3(0);
Double_t thickness[3]={1.,0.5,1.};//mm
Double_t energy;//MeV
Double_t de[3];//eloss in silicon with different thickness;
Double_t deRes[3];
for(Int_t ientries=0;ientries<10000;ientries++)
{
energy=A*e_random->Uniform(0.,35.);
for(int j=0;j<3;j++)
{
de[j]=eloss(energy,thickness[j],d2e,e2d);
deRes[j]=e_random->Gaus(de[j],de[j]*0.01/2.355);
energy=energy-de[j];
}
d1d2Res->Fill(deRes[1],deRes[0]);
d2d3Res->Fill(deRes[2],deRes[1]);
}
}
//%jsroot on
/*
TCanvas c1;
TH2I*de12Res0= new TH2I("de12Res0","de12Res0",350.,0.,350.,500.,0.,500.);
TH2I*de12Res1= new TH2I("de12Res1","de12Res1",350.,0.,350.,500.,0.,500.);
TH2I*de12Res2= new TH2I("de12Res2","de12Res2",350.,0.,350.,500.,0.,500.);
TH2I*de23Res0= new TH2I("de23Res0","de23Res0",350.,0.,350.,500.,0.,500.);
TH2I*de23Res1= new TH2I("de23Res1","de23Res1",350.,0.,350.,500.,0.,500.);
TH2I*de23Res2= new TH2I("de23Res2","de23Res2",350.,0.,350.,500.,0.,500.);
PIDRes("He4elossinSi.txt",4,de12Res0,de23Res0);
PIDRes("He6elossinSi.txt",6,de12Res1,de23Res1);
PIDRes("He8elossinSi.txt",8,de12Res2,de23Res2);
TH2I *de12Res_t= new TH2I("de12Res_t","de12Res_t",350.,0.,350.,500.,0.,500.);
de12Res_t->Add(de12Res0);
de12Res_t->Add(de12Res1);
de12Res_t->Add(de12Res2);
de12Res_t->Draw("colz");
c1.Draw();
*/
%jsroot on
TCanvas c_t1;
TH2D *de12Res=new TH2D("de12Res","de12Res",350.,0.,350.,500.,0.,500.);
TH2D *de23Res=new TH2D("de23Res","de23Res",350.,0.,350.,500.,0.,500.);
PIDRes("H1elossinSi.txt",1,de12Res,de23Res);
PIDRes("H1elossinSi.txt",2,de12Res,de23Res);
PIDRes("H1elossinSi.txt",3,de12Res,de23Res);
PIDRes("He4elossinSi.txt",3,de12Res,de23Res);
PIDRes("He4elossinSi.txt",4,de12Res,de23Res);
PIDRes("He6elossinSi.txt",6,de12Res,de23Res);
PIDRes("Li6elossinSi.txt",6,de12Res,de23Res);
PIDRes("Li6elossinSi.txt",7,de12Res,de23Res);
PIDRes("Li6elossinSi.txt",8,de12Res,de23Res);
PIDRes("Li6elossinSi.txt",9,de12Res,de23Res);
PIDRes("Be9elossinSi.txt",7,de12Res,de23Res);
PIDRes("Be9elossinSi.txt",9,de12Res,de23Res);
PIDRes("Be9elossinSi.txt",10,de12Res,de23Res);
PIDRes("Be9elossinSi.txt",11,de12Res,de23Res);
PIDRes("B10elossinSi.txt",8,de12Res,de23Res);
PIDRes("B10elossinSi.txt",10,de12Res,de23Res);
PIDRes("B10elossinSi.txt",11,de12Res,de23Res);
PIDRes("B10elossinSi.txt",12,de12Res,de23Res);
PIDRes("B10elossinSi.txt",13,de12Res,de23Res);
PIDRes("C12elossinSi.txt",10,de12Res,de23Res);
PIDRes("C12elossinSi.txt",11,de12Res,de23Res);
PIDRes("C12elossinSi.txt",12,de12Res,de23Res);
PIDRes("C12elossinSi.txt",13,de12Res,de23Res);
PIDRes("C12elossinSi.txt",14,de12Res,de23Res);
PIDRes("C12elossinSi.txt",15,de12Res,de23Res);
PIDRes("N14elossinSi.txt",12,de12Res,de23Res);
PIDRes("N14elossinSi.txt",13,de12Res,de23Res);
PIDRes("N14elossinSi.txt",14,de12Res,de23Res);
PIDRes("N14elossinSi.txt",15,de12Res,de23Res);
PIDRes("N14elossinSi.txt",16,de12Res,de23Res);
PIDRes("N14elossinSi.txt",17,de12Res,de23Res);
PIDRes("O16elossinSi.txt",14,de12Res,de23Res);
PIDRes("O16elossinSi.txt",15,de12Res,de23Res);
PIDRes("O16elossinSi.txt",16,de12Res,de23Res);
PIDRes("O16elossinSi.txt",17,de12Res,de23Res);
PIDRes("O16elossinSi.txt",18,de12Res,de23Res);
PIDRes("O16elossinSi.txt",19,de12Res,de23Res);
PIDRes("O16elossinSi.txt",20,de12Res,de23Res);
de12Res->GetXaxis()->SetTitle("de2 (MeV)");
de12Res->GetYaxis()->SetTitle("de1 (MeV)");
gPad->SetLogz();
//de12Res->Draw("colz");
gStyle->SetPalette(1);
de12Res->Draw("colz");
c_t1.Draw();
%jsroot on
TCanvas c_t2;
de23Res->GetXaxis()->SetTitle("de3 (MeV)");
de23Res->GetYaxis()->SetTitle("de2 (MeV)");
gPad->SetLogz();
//de23Res->Draw("colz");
gStyle->SetPalette(1);
de23Res->Draw("colz");
c_t2.Draw();
!jupyter nbconvert 6.PIDRes-function.ipynb --to html