The Frisch-grid ionization chamber (GIC) 信号模拟¶
致谢¶
- 感谢北京大学重离子物理研究所张国辉教授研究组的刘杰、崔增琪等同学提供照片及实验数据，感谢刘杰同学的答疑解惑。
- 实验装置参见：J.Liu etc., Nucl. Instrum. Methods Phys. Res. A 1014 (2021) 165751
探测器设置¶
GIC Setup：¶
The electrodes of the GIC are squares with both sides of 156.0 mm in length. The distances of cathode-grid and anode-grid are 60.35 and 15.24 mm, respectively. The grid consists of gold-coated tungsten wires with the radius of 0.05 mm in parallel and the distance between two adjacent wires is 2.0 mm.
气体：0.922 bar Ar+3.45%CO2
- 阴极电压：-1500V
阳极电压：+750V
Preamplifier：MPR-1
- RC=27$\mu s$, In electron-sensitive operation
探测器照片 & 波形¶
alpha source:¶
- A back-to-back compound alpha source consisting of 234U, 239Pu, 238Pu, 244Cm, with alpha energies of 4.775 MeV, 5.155 MeV, 5.499 MeV, and 5.805 MeV, respectively, is placed at the center of the common cathode to generate electron–ion pairs in the working gas. The activity of the alpha source is about 16 Bq (30%, 14%, 42% and 14% for the four isotopes, respectively).
模拟步骤. 1¶
- 读入有SRIM软件计算的alpha粒子的能损文件，生成dE/dx~ E曲线
- 由dE/dx ~ E 计算给定能量入射粒子的bragg curve
读入srim软件计算的能损文件，生成dE/dx ~ E曲线¶
In [1]:
TGraph* GetdEdxSrim(string elossfile)
{
ifstream inf(elossfile);//读入SRIM软件生成的12C入射到H2O的数据
if(!inf.is_open()){ //文件不存时。。
cout<<"Data File: "<<elossfile<< " does not exist!"<<endl;
return 0;
}
auto *g=new TGraph;//定义新的二维散点图(TGraph) g
string ss;
string sunit;//能量单位
double a, e,ededx, ndedx;//ededx-电子阻止本领，ndedx-核阻止本领
for(int j=0;j<27;j++) getline(inf,ss);//读入文件中开始的24行描述性文本
int i=0;
while(inf >> e >> sunit >> ededx >> ndedx >> a >> ss >> a >> ss >> a >> ss)
{
if(abs(e)<1.0e-4) break;//对于非数据，e值为零，此时停止读数据
if(sunit=="keV") e /=1e3;
if(sunit=="GeV") e *= 1e3;
g->SetPoint(i++, e, (ededx+ndedx)* 1.6499E-01);
}
inf.close();
return g;
}
In [2]:
string elossfile="alpha_0.922barAr+3.45%CO2.txt";
TGraph *gdEdx_E=GetdEdxSrim(elossfile);
float Emax=gdEdx_E->GetPointX(gdEdx_E->GetN()-1);
cout<<"Maximum Ein ="<<Emax<<" MeV"<<endl;
gdEdx_E->SetTitle("dE/dx(MeV/mm) vs. MeV for Alpha in 0.922bar Ar+3.45%CO2");//设定图的标题
gdEdx_E->SetLineColor(kGreen);
In [3]:
//%jsroot on
TCanvas *c1=new TCanvas;
gdEdx_E->Draw();
c1->Draw();
计算bragg curve (dE/dx ~ x for a given $E_{in}$)¶
In [4]:
TGraph* GetdEdx_x(float E, TGraph *gdEdx_E)
{
TGraph *g1=new TGraph;
g1->SetTitle(Form(" dE/dx(MeV/mm) vs. x for alpha with Ein= %.2f MeV",E));
float dE=0.0001;//MeV
float dx;
float x=0;
int i=0;
while(E>0) {
g1->SetPoint(i,x,gdEdx_E->Eval(E));
E-=dE;
dx=dE/gdEdx_E->Eval(E);
x+=dx;
i++;
}
return g1;
}
In [5]:
float Ein=5.805;//MeV
TGraph *gdEdx_x=GetdEdx_x(Ein,gdEdx_E);
In [6]:
c1->Clear();
gdEdx_x->Draw();
c1->Draw();
模拟步骤 2¶
- 计算探测器电流脉冲，电荷脉冲
- 计算电荷灵敏前置放大器的电压脉冲
计算探测器的电流脉冲信号¶
- 粒子从坐标原点沿$\theta$方向出射，电离发生在下图的黄色区域。在x方向上的原始电离产生的电子数目沿x的分布可由bragg curve得到。
- 即Nx(x)=gdEdx_x->Eval(x)
- 假设原始电离在t=0时刻产生，在t=t1时漂移电子的最前沿到达栅极，t=t2时漂移电子后延到达栅极。
- 利用Schottky-Ramo Theorem计算电流脉冲
In [7]:
float theta=0.;
In [8]:
void GIC_CurrentPulse( float theta,TGraph *gdEdx_x, TGraph *gdcurrent, TGraph *gacurrent)
{
float dcg=60.35;//mm
float dag=15.24;//mm
float vcg=40;//mm/us, drift velocity of electron in a given condition.
float vag=40;//mm/us
gacurrent->SetTitle(Form("Anode current pulse for #theta=%.1f degree",theta));
gdcurrent->SetTitle(Form("Cathod current pulse for #theta=%.1f degree",theta));
gacurrent->SetMarkerColor(kRed);
gdcurrent->SetMarkerColor(kGreen);
gacurrent->SetLineColor(kRed);
gdcurrent->SetLineColor(kGreen);
theta *=3.14159/180.;
float acurrent[2000]={0.};
float dcurrent[2000]={0.};
float xmax=gdEdx_x->GetPointX(gdEdx_x->GetN()-1);
float dx=xmax/1000;
for(int i=0;i<=xmax/dx;i++) {
float x=i*dx;
float Nx=gdEdx_x->Eval(x)*dx;// number of primary ionization, 30eV per ion pair
float d=dcg-x*cos(theta);//drift path
float t=d/vcg*1000;//drift time(ns)
for(int j=0;j<t;j++) {
dcurrent[j]+=Nx/(dcg/vcg)/1000.;
}
for(int j=t;j<t+dag/vag*1000;j++) {
acurrent[j]+=Nx/(dag/vag)/1000.;
}
}
for(int i=0;i<2000;i++) {
gdcurrent->SetPoint(i,i,dcurrent[i]);
gacurrent->SetPoint(i,i,acurrent[i]);
}
}
In [9]:
void GIC_ChargePulse(TGraph *gdcurrent, TGraph *gacurrent, TGraph *gdcharge, TGraph *gacharge)
{
double tmpx,tmpy;
float dcharge=0;
float acharge=0;
gacharge->SetTitle(Form("Anode charge pulse for #theta=%.1f degree",theta));
gdcharge->SetTitle(Form("Cathod charge pulse for #theta=%.1f degree",theta));
gacharge->SetMarkerColor(kRed);
gdcharge->SetMarkerColor(kGreen);
gacharge->SetLineColor(kRed);
gdcharge->SetLineColor(kGreen);
for(int i=0;i<gdcurrent->GetN();i++)
{
gdcurrent->GetPoint(i,tmpx,tmpy);
dcharge+=tmpy;
gdcharge->SetPoint(i,tmpx,dcharge);
}
for(int i=0;i<gacurrent->GetN();i++)
{
gacurrent->GetPoint(i,tmpx,tmpy);
acharge+=tmpy;
gacharge->SetPoint(i,tmpx,acharge);
}
}
In [10]:
TGraph *gdcurrent=new TGraph;
TGraph *gacurrent=new TGraph;
TGraph *gdcharge=new TGraph;
TGraph *gacharge=new TGraph;
In [11]:
GIC_CurrentPulse(theta,gdEdx_x,gdcurrent,gacurrent);
GIC_ChargePulse(gdcurrent,gacurrent,gdcharge,gacharge);
double aqmax = TMath::MaxElement(gacharge->GetN(),gacharge->GetY());//get x minimum for TGraph
double dqmax = TMath::MaxElement(gdcharge->GetN(),gdcharge->GetY());
cout<<"cathod Q:"<<dqmax<<", anode Q:"<<aqmax<<endl;
current pulse , Green-cathod, Red-Anode¶
In [12]:
c1->Clear();
gacurrent->Draw();
gdcurrent->Draw("same");
c1->Draw();
charge pulse, Green-cathod, Red-Anode¶
In [13]:
c1->Clear();
gacharge->Draw();
gdcharge->Draw("same");
c1->Draw();
the Response of a Charge-Sensitive Preamplifier¶
- through convolution
$\tau=RC$
$$ I(t)=0, U(t)=0 \qquad for \quad t<0\\ U(t)=-\frac{1}{C_f}\int_0^t{I(t')\cdot e^{-(t-t')/\tau}dt'} $$In [14]:
float getUval(float t,TGraph *gcurrent,float tau)
{
float usum=0;
float dt=10;
float tmax=gcurrent->GetPointX(gcurrent->GetN()-1);
for(int i=0;i<=t/dt;i++) {//1ns per bin;
float t1=i*dt;
float Ival= t1>tmax ? 0 : gcurrent->Eval(t1);
usum += Ival*exp(-(t-t1)/tau)*dt;
}
return usum;
}
In [15]:
void PreampPulse(float tau, TGraph* gdcurrent, TGraph *gacurrent,TGraph *gdu, TGraph *gau)
{
gau->SetMarkerColor(kRed);
gdu->SetMarkerColor(kGreen);
gau->SetLineColor(kRed);
gdu->SetLineColor(kGreen);
int ig=0;
float dt=10;
for(int i=-5000./dt;i<30000/dt;i++){//up to 30us
float t=i*dt;
if(t<0){
gdu->SetPoint(ig,t,0);
gau->SetPoint(ig,t,0);
}
else {
float ud=getUval(t,gdcurrent,tau);
gdu->SetPoint(ig,t,ud);
float ua=getUval(t,gacurrent,tau);
gau->SetPoint(ig,t,ua);
}
ig++;
}
}
In [16]:
float GetRiseTime(TGraph *g1)//debug
{
double ymax = TMath::MaxElement(g1->GetN(),g1->GetY());
double ya=0.1*ymax;
double yb=0.9*ymax;
TGraph *gt=new TGraph;
int igt=0;
for(int i=0;i<g1->GetN();i++){
double x,y;
g1->GetPoint(i,x,y);
if(y>0.99*ymax) break;
if(y>0) {
gt->SetPoint(igt++,y,x);
}
}
return gt->Eval(yb)-gt->Eval(ya);
}
voltage pulse of preamplifier, Green-cathod, Red-Anode¶
In [17]:
float tau=27*1000.;//us; RC for preamplifier
TGraph *gdu=new TGraph;
TGraph *gau=new TGraph;
TGraph *gdut=new TGraph;
TGraph *gaut=new TGraph;
PreampPulse(tau,gdcurrent,gacurrent,gdu,gau);//return gdu,gau;
double aumax = TMath::MaxElement(gau->GetN(),gau->GetY());//get x minimum for TGraph
double dumax = TMath::MaxElement(gdu->GetN(),gdu->GetY());
cout<<"cathod Vmax:"<<dumax<<", anode Vmax:"<<aumax<<endl;
double durt=GetRiseTime(gdu);
double aurt=GetRiseTime(gau);
cout<<"rise time of cathod signal:"<<durt<<", anode signal:"<<aurt<<endl;
In [18]:
c1->Clear();
gau->SetMinimum(-0.5);
gau->Draw();
gdu->Draw("same");
c1->Draw();
TTree output¶
In [19]:
double ein[9]={1,2,3,4,4.775, 5.155, 5.499,5.805,6};
double Ein,theta;
double dcurrent[2000],acurrent[2000];
double dcharge[2000],acharge[2000];
double dqmax,aqmax;
double du[3500],au[3500];
double dumax,aumax;
double durt,aurt;
TFile *f=new TFile("gic.root","recreate");
TTree *tree=new TTree("tree","gic tree");
tree->Branch("e",&Ein,"e/D");
tree->Branch("theta",&theta,"theta/D");
tree->Branch("ai",&acurrent,"ai[2000]/D");
tree->Branch("ci",&dcurrent,"ci[2000]/D");
tree->Branch("aq",&acharge,"aq[2000]/D");
tree->Branch("cq",&dcharge,"cq[2000]/D");
tree->Branch("au",&au,"au[3500]/D");
tree->Branch("cu",&du,"cu[3500]/D");
tree->Branch("aqmax",&aqmax,"aqmax/D");
tree->Branch("cqmax",&dqmax,"cqmax/D");
tree->Branch("aumax",&aumax,"aumax/D");
tree->Branch("cumax",&dumax,"cumax/D");
tree->Branch("aurt",&aurt,"aurt/D");
tree->Branch("curt",&durt,"curt/D");
for(int ie=0;ie<9;ie++) {
for(int ia=0;ia<=16;ia++) {
Ein=ein[ie];
theta=ia*5;
TGraph *gdEdx_x=GetdEdx_x(Ein,gdEdx_E);
TGraph *gdcurrent=new TGraph;//2000
TGraph *gacurrent=new TGraph;
TGraph *gdcharge=new TGraph;
TGraph *gacharge=new TGraph;
GIC_CurrentPulse(theta,gdEdx_x,gdcurrent,gacurrent);
GIC_ChargePulse(gdcurrent,gacurrent,gdcharge,gacharge);
for(int i=0;i<2000;i++) {
dcurrent[i]=gdcurrent->GetY()[i];
acurrent[i]=gacurrent->GetY()[i];
dcharge[i]=gdcharge->GetY()[i];
acharge[i]=gacharge->GetY()[i];
}
aqmax = TMath::MaxElement(gacharge->GetN(),gacharge->GetY());//get x minimum for TGraph
dqmax = TMath::MaxElement(gdcharge->GetN(),gdcharge->GetY());
float tau=27.*1000.;//us; RC for preamplifier
TGraph *gdu=new TGraph;//3500
TGraph *gau=new TGraph;
PreampPulse(tau,gdcurrent,gacurrent,gdu,gau);//return gdu,gau;
for(int i=0;i<3500;i++) {
du[i]=gdu->GetY()[i];
au[i]=gau->GetY()[i];
}
aumax = TMath::MaxElement(gacharge->GetN(),gacharge->GetY());//get x minimum for TGraph
dumax = TMath::MaxElement(gdcharge->GetN(),gdcharge->GetY());
durt=GetRiseTime(gdu);
aurt=GetRiseTime(gau);
cout<<ie<<","<<ia<<";";
tree->Fill();
}
cout<<endl;
}
tree->Write();
f->Close();