1. 读取dedx文件

In [1]:
//读取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();   
}

2. 入射粒子的射程

In [2]:
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;
}
In [3]:
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);
}
In [4]:
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;
}
the range of 5.486MeV 4He in Silicon is 0.0277423 mm
the range of 10MeV 4He in Silicon is 0.0692679 mm
the range of 100MeV 4He in Silicon is 3.48517 mm
the range of 200MeV 4He in Silicon is 12.0104 mm

3. 入射粒子在探测器中的能损

以10 MeV,100 MeV,1000 MeV的4He粒子先后穿过厚度为1000 um,500 um,1000 um的能量损失;

In [5]:
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;
}
In [6]:
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;
}
the eloss of 4He with 10 MeV in 1000 um Silicon is 10 MeV;
the eloss of 4He with 0 MeV in 500 um Silicon is 0 MeV;
the eloss of 4He with 0 MeV in 1000 um Silicon is 0 MeV;
 
the eloss of 4He with 100 MeV in 1000 um Silicon is 17.3812 MeV;
the eloss of 4He with 82.6188 MeV in 500 um Silicon is 9.8863 MeV;
the eloss of 4He with 72.7325 MeV in 1000 um Silicon is 24.0471 MeV;
 
the eloss of 4He with 200 MeV in 1000 um Silicon is 9.48408 MeV;
the eloss of 4He with 190.516 MeV in 500 um Silicon is 4.88291 MeV;
the eloss of 4He with 185.633 MeV in 1000 um Silicon is 10.0875 MeV;
 

4. He同位素望远镜粒子鉴别图

In [7]:
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]);
    }
}
In [8]:
//%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();
*/

5. 同为素H,He,Li,Be,B,C,N,O的PID(0~35 MeV/u)

  • 添加1%的能量分辨,并用将结果存入TH2画中;
In [9]:
%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();
In [10]:
%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();
  • 注意:TH2图分bin不要太细,会导致内存占满死机。
In [11]:
!jupyter nbconvert 6.PIDRes-function.ipynb --to html
[NbConvertApp] Converting notebook 6.PIDRes-function.ipynb to html

[NbConvertApp] Writing 325730 bytes to 6.PIDRes-function.html

In [ ]: