#include "Rtypes.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "Dioda.h"
#define JMAX 40
#define STEPH 25e-9
ClassImp(Dioda)
//______________________________________________________________________________
Dioda::Dioda( TF1 *x0, Float_t x1,Int_t x2,Int_t x3)
{
//Constructor of the class Dioda:
// TF1 *x0 ; distribution of the Neff(x) given in 10**12! Note that p+ contact starts with x=301 and n+ at x=1
// Int_t x1 ; voltage applied in reverse direction
// Int_t x2 ; width of the diode in um. For the simulation it is not neccessary to put real dimensions
// Int_t x3 ; thickness in um
nx=x2;
ny=x3;
Voltage=x1;
Neff=x0;
//Public Part
ran=new TRandom(33);
pos = new TH1F("charge+","Positive Charge",200,0,STEPH);
neg = new TH1F("charge-","Negative Charge",200,0,STEPH);
sum = new TH1F("charge","Total Charge",200,0,STEPH);
sum->SetXTitle("t [s]"); neg->SetXTitle("t [s]"); pos->SetXTitle("t [s]");
sum->SetYTitle("I [arb.]");neg->SetYTitle("I [arb.]");pos->SetYTitle("I [arb.]");
sum->GetYaxis()->SetTitleOffset(1.4); neg->GetYaxis()->SetTitleOffset(1.4); pos->GetYaxis()->SetTitleOffset(1.4);
// sum->SetLabelSize(0.045,"X"); neg->SetLabelSize(0.045,"X"); pos->SetLabelSize(0.045,"X");
// sum->SetLabelSize(0.045,"Y"); neg->SetLabelSize(0.045,"Y"); pos->SetLabelSize(0.045,"Y");
MobMod=0;
B=0.0;
average=1;
diff=0;
GetRamoField(); GetField();
Temperature=263;
trapping=0;
udregion=0;
material=0;
Multiplication=0;
}
//_________________________________________________________________________________
Float_t Dioda::rtbis(float x1, float x2, float xacc)
{
void nrerror(char error_text[]);
int j;
float dx,f,fmid,xmid,rtb;
f=PoEqSolve(x1);
fmid=PoEqSolve(x2);
if (f*fmid >= 0.0) nrerror("Root must be bracketed for bisection in rtbis");
rtb = f < 0.0 ? (dx=x2-x1,x1) : (dx=x1-x2,x2);
for (j=1;j<=JMAX;j++) {
fmid=PoEqSolve(xmid=rtb+(dx *= 0.5));
if (fmid <= 0.0) rtb=xmid;
if (fabs(dx) < xacc || fmid == 0.0) return rtb;
}
nrerror("Too many bisections in rtbis");
return 0.0;
}
//___________________________________________________________________________________
void Dioda::rk4(float y[], float dydx[], int n, float x, float h, float yout[])
{
float *vector(long,long);
void free_vector(float*,long,long);
int i;
float xh,hh,h6,*dym,*dyt,*yt;
dym=vector(1,n);
dyt=vector(1,n);
yt=vector(1,n);
hh=h*0.5;
h6=h/6.0;
xh=x+hh;
for (i=1;i<=n;i++) yt[i]=y[i]+hh*dydx[i];
Derivs(xh,yt,dyt);
for (i=1;i<=n;i++) yt[i]=y[i]+hh*dyt[i];
Derivs(xh,yt,dym);
for (i=1;i<=n;i++) {
yt[i]=y[i]+h*dym[i];
dym[i] += dyt[i];
}
Derivs(x+h,yt,dyt);
for (i=1;i<=n;i++)
yout[i]=y[i]+h6*(dydx[i]+dyt[i]+2.0*dym[i]);
free_vector(yt,1,n);
free_vector(dyt,1,n);
free_vector(dym,1,n);
}
//_______________________________________________________________________________
Dioda::~Dioda()
{
delete pos;
delete neg;
delete sum;
Clear();
}
void Dioda::Derivs(float x,float y[],float dydx[])
{
Double_t permMat=(material==0)?perm:permDi;
Double_t permV = perm0*1e-6;;
dydx[1]=y[2];
dydx[2]=-Neff->Eval(x)*e_0/(permMat*permV);
}
Float_t Dioda::PoEqSolve(Float_t der)
{
//TArrayF PhyField(ny+2);
//TArrayF PhyPot(ny+2);
Float_t Step=1;
Int_t i;
Float_t h,x=0;
Float_t y[3],dydx[3],yout[3];
h=Step;
y[1]=Voltage;
y[2]=der;
PhyField[0]=y[2]; PhyPot[0]=y[1];
Derivs(x,y,dydx);
for (i=1;i<=ny;i++)
{
rk4(y,dydx,2,x,h,yout);
//printf("%f %f %fn",x+h,yout[1],yout[2]);
y[1]=yout[1];
y[2]=yout[2];
PhyField[i]=y[2]; PhyPot[i]=y[1];
x=x+h;
Derivs(x,y,dydx);
}
// printf("y[1]=%fn",xp1);
return y[1];
}
// void Dioda::fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
// {
// Neff->SetParameters(par);
// GetField();
// for(int i=0;i<PhyField.GetSize();i++) if(PhyField[i]>0) {f=1e5; return;};
// sum->Reset();
// LaserV(entryp,exitp,20);
// //el->preamp(d->sum);
// sum->GetXaxis()->Set(200,0,100);
// f=HDif(sum,neg,0,35);
// sum->GetXaxis()->Set(200,0,100e-9);
// printf("DIff2=%fn",f);
// return;
// }
// void Dioda::Minimize(Int_t npar, Double_t *vstart, Double_t *step)
// {// The z values
// TMinuit *gMinuit = new TMinuit(npar); //initialize TMinuit with a maximum of 5 params
// gMinuit->SetFCN((void *)fcn);
// Double_t arglist[10]; Int_t ierflg = 0;
// arglist[0] = 1;
// gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
// // Set starting values and step sizes for parameters
// // static Double_t vstart[3] = {250, -1. , 1};
// // static Double_t step[3] = {4 , 0.2 , 0.2};
// Char_t parname[3]; parname[0]='a'; parname[1]='x';
// for(Int_t i=0;i<npar;i++)
// {
// parname[1]=48+i;
// gMinuit->mnparm(0, parname, vstart[i], step[i], 0,0,ierflg);
// }
// // Now ready for minimization step
// arglist[0] = 500;
// arglist[1] = 1.;
// gMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);// Print results
// Double_t amin,edm,errdef;
// Int_t nvpar,nparx,icstat;
// gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
// gMinuit->mnprin(3,amin);
// }
void Dioda::GetRamoField()
{
Int_t i,j;
Float_t Step=1;
TArrayF RamoPot(ny+2);
TArrayF RamoField(ny+2);
TArrayD RamoPot2D(nx*ny+1);
TArrayI StripPosition=TArrayI(2); StripPosition[0]=1; StripPosition[1]=nx;
for(j=1;j<=ny;j++)
for(i=1;i<=nx;i++)
RamoPot2D[i+(j-1)*nx]=(Double_t)(j-1)/(Double_t)(ny-1);
// for(i=0;i<nx*ny+1;i++) {printf("%f ",RamoPot2D[i]); if(i%nx==0) printf("n");}
Ramo=EField(RamoPot2D,nx,ny,StripPosition);
Ramo.CalField(Step,1);
}
void Dioda::GetRamoField(TH1F *rf)
{
//Change made by GK on 16.10.2002
Int_t i,j;
Float_t Step=1;
TArrayF RamoPot(ny+2);
TArrayF RamoField(ny+2);
TArrayD RamoPot2D(nx*ny+1);
TArrayI StripPosition=TArrayI(2); StripPosition[0]=1; StripPosition[1]=nx;
for(j=1;j<=ny;j++)
{
for(i=1;i<=nx;i++)
RamoPot2D[i+(j-1)*nx]=rf->GetBinContent(j);
// printf("%e, n",RamoPot2D[i+(j-1)*nx]);
}
// for(i=0;i<nx*ny+1;i++) {printf("%f ",RamoPot2D[i]); if(i%nx==0) printf("n");}
Ramo=EField(RamoPot2D,nx,ny,StripPosition);
Ramo.CalField(Step,1);
}
void Dioda::GetField()
{
Int_t i,j;
Float_t Step=1;
Float_t aa;
PhyPot=TArrayF(ny+2);
PhyField=TArrayF(ny+2);
TArrayD PhyPot2D(nx*ny+1);
TArrayI StripPosition=TArrayI(2); StripPosition[0]=1; StripPosition[1]=nx;
aa=rtbis(-100,100,0.000001);
//if(!Invert && GetDepletionVoltage()>Voltage) aa=rtbis(1,300,0.000001); else aa=rtbis(-10,10,0.000001);
for(j=1;j<=ny;j++)
for(i=1;i<=nx;i++)
{
PhyPot2D[i+(j-1)*nx]=(Double_t)PhyPot[j-1];
}
Real=EField(PhyPot2D,nx,ny,StripPosition);
Real.CalField(Step,1);
// for(i=0;i<nx*ny+1;i++) {printf("%f ",PhyPot2D[i]); if(i%nx==0) printf("n");}
}
void Dioda::GetField(TF1 *potential)
{
//Set the field as defined in the potential function.
Int_t i,j,index;
Float_t aa,val;
Double_t abspot=0;
PhyPot=TArrayF(ny+2);
PhyField=TArrayF(ny+2);
TArrayD PhyPot2D(nx*ny+1);
TArrayI StripPosition=TArrayI(2); StripPosition[0]=1; StripPosition[1]=nx;
Real=EField(PhyPot2D,nx,ny,StripPosition);
abspot=potential->Integral(0,300);// printf("napetost=%en",abspot);
for(j=1;j<=ny;j++)
{
for(i=1;i<=nx;i++)
{
index=i+(j-1)*nx;
Real.SetEfx(index,0);
Real.SetEfy(index,-potential->Eval((Double_t)(ny-j)));
Real.SetEf (index,-potential->Eval((Double_t)(ny-j)));
}
PhyField[j]=-potential->Eval((Double_t)(ny-j));
if(j==1) PhyPot[j]=abspot; else
PhyPot[j]=potential->Integral((Double_t)(ny-j),(Double_t)(ny-j-1))+PhyPot[j-1];
}
//Real=EField(PhyPot2D,nx,ny,StripPosition);
//Real.CalField(Step,1);
// for(i=0;i<nx*ny+1;i++) {printf("%f ",PhyPot2D[i]); if(i%nx==0) printf("n");}
}
TGraph *Dioda::Draw(char *option)
{
//Draws potential "p" or electric field "f"
Char_t Opt[10];
TGraph *gr;
TArrayF xx=TArrayF(ny+1);
for(Int_t i=0;i<=ny;i++) xx[i]=(Float_t)i;
if(strchr(option,'p')!=NULL) gr=new TGraph(ny,xx.GetArray(),PhyPot.GetArray());
if(strchr(option,'f')!=NULL) gr=new TGraph(ny,xx.GetArray(),PhyField.GetArray());
if(strchr(option,'s')!=NULL) sprintf(Opt,"CP"); else sprintf(Opt,"ACP");
//if(! strcmp(option,"p")) gr=new TGraph(ny,xx.GetArray(),PhyPot.GetArray());
//if(! strcmp(option,"f")) gr=new TGraph(ny,xx.GetArray(),PhyField.GetArray());
gr->Draw(Opt);
gr->GetHistogram()->Draw();
gr->Draw(Opt);
return gr;
}
void Dioda::Drift(Double_t sx, Double_t sy,Float_t charg,DStruct *seg,Int_t MobMod,Float_t t0)
{
//Drift simulation for a point charge (charg)
//starting from ( sx,sy)
//DStruct *seg is the structure where the the drift paths, drift times and induced cahrges are stored
// Int_t MobMod; mobility model
// Float_t B; magnetic field
// Float_t t0=start time ;
double th=0;
//Double_t mobility(Double_t *, Double_t *);
//TF2 *mob = new TF2("mobility",mobility,0,10000,0,1e13,2);
Double_t pardif[3]={0,1430,300},sigma,dif,gg,difstepx=0,difstepy=0,difx=0,dify=0;
float *xcv,*ycv,*time,*charge,E[3],Eprev,En[3],Erprev,Enr[3],Er[3],pathlen=0;
long st=1,i,ishit=0;
double cx,cy,vel,t,sumc=0,cutx=1.,cuty=1.,deltacy;
double lastsecpath=0,mfp;
int l;
int nostrips=Real.GetNoStrips();
int *strpos=Real.GetStripPosition();
int nx=Real.GetNX();
int ny=Real.GetNY();
//Inclusion of Magnetic field 28.8.2001
float muhe=1650;
float muhh=310;
//End of inclusion of Magnetic field 28.8.2001
t=t0;
//if(charg>0) { xcv=(*seg).xc.pos; ycv=(*seg).yc.pos; time=(*seg).time.pos; charge=(*seg).charge.pos;}
// else {xcv=(*seg).xc.neg; ycv=(*seg).yc.neg; time=(*seg).time.neg; charge=(*seg).charge.neg;}
xcv=(*seg).Xtrack; ycv=(*seg).Ytrack; time=(*seg).Time; charge=(*seg).Charge;
(*seg).Clear();
(*seg).PCharge=(Int_t)charg;
// for(i=0;i<MAXPOINT;i++) {xcv[i]=0.; ycv[i]=0; time[i]=0; charge[i]=0;}
cx=sx; cy=sy;
xcv[st]=cx;
ycv[st]=cy;
time[st]=0.;
charge[st]=0;
Real.CalFieldXY(cx,cy,E);
Ramo.CalFieldXY(cx,cy,Er);
//printf("n#%d, X=%f, Y=%f, t=%f, charge=%f, sumc=%f, Ex=%f ,Ey=%f, th=%f cutx,y=%f %fn",st,cx,cy,t,charge[st],sumc,E[1],E[2],th,cutx,cuty);
// printf("%d %f %f %f %f %f %f %f %fn",st,cx,cy,t,charge[st],sumc,E[1],E[2],th);
if(cy<udregion && udregion!=0) ishit=1;
while(!ishit)
{
st++;
// Inclusion of the Magnetic field 28.8.2001
if(B!=0)
{
if(charg>0) {E[2]=E[2]-muhh*E[1]*B; E[1]=E[1]+muhh*E[2]*B;} else {E[2]=E[2]-muhe*E[1]*B; E[1]=E[1]+muhe*E[2]*B;}
E[0]=TMath::Sqrt(E[1]*E[1]+E[2]*E[2]);
}
// End of magnetic field inclusion 28.8.2001
//
th=TMath::Abs(atan(E[2]/E[1]));
//if(cy>(float)ny) cut=((float)ny-oldcy); else cut=1;
if(E[2]<=0) {deltacy=charg*sin(th);} else {deltacy=-charg*sin(th);}
if(diff) if(st>2) {gg=cuty*1e-4/vel; sigma=TMath::Sqrt(2*Kboltz*Real.Mobility(E[0],Temperature,charg,TMath::Abs(Neff->Eval(cy)),MobMod)*Temperature*gg); dif=ran->Gaus(0,sigma); dify=dif*1e4; dif=ran->Gaus(0,sigma); difx=dif*1e4;}
if((cy+deltacy+dify)>=ny)
{
if(cy==ny) {cuty=0; cutx=1; difx=0; dify=0;}
else {if(cy+deltacy<ny) {cuty=1; cutx=1; dify=(Double_t) ny-(cy+deltacy);} else {cuty=(ny-cy)/sin(th); cutx=cuty; difx=0; dify=0;} }
} else {cutx=1; cuty=1;}
if((cy+deltacy+dify)<1)
{
if(cy+deltacy<1) {cuty=(cy-1)/sin(th); cutx=cuty; difx=0; dify=0;} else
{cuty=1; cutx=1; dify=(Float_t) cy-(1+deltacy);}
}
//printf("1::::%f %f %f %f %f %f %f %f %f n",cx,cy,deltacy,th,charg,cutx,cuty,difx,dify);
//if(st>2) {gg=cuty*1e-4/vel; sigma=TMath::Sqrt(2*Kboltz*1430*300*gg); dif=ran->Gaus(0,sigma); cuty+=dif*1e4;}// printf("dif=%f gg=%en",dif,gg);}
//gg=cutx*1e-4/vel; sigma=TMath::Sqrt(2*Kboltz*1430*300*gg); dif=ran->Gaus(0,sigma); cutx+=dif*1e4;}// printf("dif=%f gg=%en",dif,gg);}
if(E[1]<=0) {cx=cx+cutx*cos(th)*charg+difx; charge[st]=(cos(th)*cutx+difx)*Er[1];} else {cx=cx-cutx*cos(th)*charg+difx; charge[st]=-(cos(th)*cutx+difx)*Er[1];}
if(E[2]<=0) {cy=cy+cuty*sin(th)*charg+dify; charge[st]+=(sin(th)*cuty+dify)*Er[2];} else {cy=cy-cuty*sin(th)*charg+dify; charge[st]+=-(sin(th)*cuty+dify)*Er[2];}
// printf("2::::%f %f %f %f %f %f %f %f %f n",cx,cy,deltacy,th,charg,cutx,cuty,difx,dify);
//printf("2:::::%f %f %f %f %f %f %f %f %f %f n",cx,cy,E[1],E[2],th,charg,cutx,cuty,difx,dify);
/*
if(E[1]<=0) {cx=cx+cutx*cos(th)*charg+difx;} else {cx=cx-cutx*cos(th)*charg+difx;}
if(E[2]<=0) {cy=cy+cuty*sin(th)*charg+dify;} else {cy=cy-cuty*sin(th)*charg+dify;}
Real.CalFieldXY(cx,cy,En); Ramo.CalFieldXY(cx,cy,Enr);
if(E[1]<=0) {charge[st]=(cos(th)*cutx+difx)*((Er[1]+Enr[1])/2);} else {charge[st]=-(cos(th)*cutx+difx)*((Er[1]+Enr[1])/2);}
if(E[2]<=0) {charge[st]+=(sin(th)*cuty+dify)*((Er[2]+Enr[2])/2);} else {charge[st]+=-(sin(th)*cuty+dify)*((Er[2]+Enr[2])/2);}
*/
sumc+=charge[st];
Eprev=E[0];
Real.CalFieldXY(cx,cy,E);
Ramo.CalFieldXY(cx,cy,Er);
// Inclusion of the Magnetic field 28.8.2001
if(B!=0)
{
if(charg>0) {E[2]=E[2]-muhh*E[1]*B; E[1]=E[1]+muhh*E[2]*B;} else {E[2]=E[2]-muhe*E[1]*B; E[1]=E[1]+muhe*E[2]*B;}
E[0]=TMath::Sqrt(E[1]*E[1]+E[2]*E[2]);
}
// End of magnetic field inclusion 28.8.2001
// vel=Real.DriftVelocity((E[0]+En[0])/2,charg,Temperature,Neff->Eval(cy));
vel=Real.DriftVelocity((E[0]+Eprev)/2,charg,Temperature,TMath::Abs(Neff->Eval(cy)),MobMod);
// for(l=0;l<3;l++) E[l]=En[l];
// printf("Curent velocity (Neff=%f) (Temperature=%f) charg=%f %f %f is=%en",Neff->Eval(cy),Temperature,charg,E[0],Eprev,vel);
t=t+cutx*1e-4/vel;
pardif[0]=(Double_t) t;
//printf("%e %f %fn",pardif[0],pardif[1],pardif[2]);
//ok= difint(pardif,3e-4);
pathlen+=cutx;
xcv[st]=cx;
ycv[st]=cy;
time[st]=t;
seg->Efield[st]=E[0]; //included for avalanche effects GK 14.10.2008
//printf("time %fn",t);
//printf("i=%d;cx=%f; cy=%f; strpos[i]=%d;strpos[i+1]=%d",i,cx,cy,strpos[i],strpos[i+1]);
//printf("#%d, vel=%e , X=%f, Y=%f, t=%e, charge=%f, sumc=%f, Ex=%f ,Ey=%f, th=%f cutx,y=%f %fn",st,vel,cx,cy,t,charge[st],sumc,E[1],E[2],th,cutx,cuty);
//printf("%d %f %f %e %f %f %f %f %f %f:::%f,%f,%f,%fn",st,cx,cy,t,charge[st],sumc,E[1],E[2],th,charg,cutx,cuty,difx,dify);
//if(t>40) {printf("difussionn"); (*seg).difftrack=malloc(sizeof(struct diffusion)); }
// if(charg >0)
for(i=0;i<nostrips*2;i+=2) if(strpos[i]<=cx && strpos[i+1]>=cx && cy>=ny)
{ishit=1; (*seg).DStrip=(int)(i/2)+1; break;}
//if(charg <0)
if(cy<=1) ishit=1;
if(cx<=1 || cx>=nx) ishit=1;
if(st>=MAXPOINT-1) {ishit=1;}// printf("too many iterationsn");}
if(cy<udregion && udregion!=0) ishit=1;
// printf("ishit=%dn",ishit);
}
(*seg).Xlenght=pathlen; (*seg).Ylenght=pathlen;
(*seg).TTime=t; (*seg).TCharge=sumc; (*seg).Steps=st;
return;
#undef x
}
void Dioda::MipIR(float *enp, float *exp, int div)
{
// The simulation of the drift for the minimum ionizing particles.
// A track is devided into Int_ div buckets. Each bucket is drifted in the field. The
// induced currents for each carrier is calculated as the sum all buckets.
// Int_t MobMod; mobility model
// Float_t B; magnetic field
Double_t data[4];
TH1F trn[4],trp[4];
Float_t sp[3];
int i,j;
DStruct seg;
TH1F *histop = new TH1F("ch-","charge+",pos->GetNbinsX(),pos->GetXaxis()->GetXmin(),pos->GetXaxis()->GetXmax());
TH1F *histon = new TH1F("ch+","charge-",neg->GetNbinsX(),neg->GetXaxis()->GetXmin(),neg->GetXaxis()->GetXmax());
sum->Reset(); pos->Reset(); neg->Reset();
for(i=0;i<div;i++)
{
for(j=0;j<3;j++) {sp[j]=((exp[j]-enp[j])/div)*i+enp[j]+(exp[j]-enp[j])/(2*div);}
// printf("#i=%d div=%d, pointx=%f, pointy=%f n",i,div,seg.sp[0],seg.sp[1]);
for(j=0;j<average;j++)
{
Drift(sp[0],sp[1],1,&seg,MobMod);
seg.GetCH(histop,1);
//it is assumed that holes don't contribute to ionization
Drift(sp[0],sp[1],-1,&seg,MobMod);
seg.GetCH(histon,1);
if(Multiplication)
{
CalM(&seg,data);
printf("Multiplication M=%f Xa=%f Ya=%f t0=%en",data[0],data[1],data[2],data[3]);
Drift(data[1],data[2],1,&seg,MobMod,data[3]);
seg.GetCH(histop,1,data[0]);
Drift(data[1],data[2],-1,&seg,MobMod,data[3]);
seg.GetCH(histon,1,data[0]);
}
}
histop->Scale(1/((Float_t)average)); histon->Scale(1/((Float_t)average));
// histop->Scale(0.05); histon->Scale(0.05);
if(trapping) {
ht->Trapping(histop,trp);
et->Trapping(histon,trn);
pos->Add(&trp[3]); neg->Add(&trn[3]);
} else {pos->Add(histop);neg->Add(histon);}
histop->Reset();
histon->Reset();
}
sum->Add(neg);
sum->Add(pos);
delete histop;
delete histon;
}
Double_t laser(Double_t *x, Double_t *par)
{
Double_t dvg=0;
dvg=par[0]*TMath::Exp(- TMath::Power((*x-par[1]),2)/par[2])+par[3]+par[4]*(*x)+par[5]*TMath::Power(*x,2);
return dvg;
}
void Dioda::LaserV(float *enp, float *exip, int div,Float_t lambda)
{
// The simulation of the drift for the red laser illumination!
// A penetration profile is devided into Int_ div buckets. Each bucket is drifted in the field. The
// induced currents for each carrier is calculated as the sum all buckets.
TH1F trn[4],trp[4];
Float_t sp[3];;
DStruct seg;
Double_t MaxPudu=2.7,pudu=0.;
Int_t PuduStep=27;
int i,ii=0,j,k=0,z,cutoff=27,shiftn,shiftp;
//void printsegdrift(struct segdrift *);
Double_t *cc= new Double_t[cutoff];
Double_t *tt= new Double_t[cutoff];
Double_t cd;
float xe,xs;//= new Double_t[cutoff];
TH1F *histop = new TH1F("chl+","charge+",200,0,STEPH);
TH1F *histon = new TH1F("chl-","charge-",200,0,STEPH);
TH1F *outhisp,*outhisn;
Double_t par[6]={-0.038553,1.07634,1.14664,1.54445e-2,-6.08484e-3,6.08802e-4};
TF1 *las=new TF1("laser",laser,0,5,6);
las->SetParameters(par);
while(pudu<MaxPudu)
{
cc[k]=las->Integral(pudu,pudu+MaxPudu/PuduStep)/las->Integral(MaxPudu/PuduStep,MaxPudu);
tt[k]=pudu*1e-9;
// printf("%e %en",cc[k],tt[k]);
pudu+=MaxPudu/PuduStep; k++;
}
sum->Reset(); pos->Reset(); neg->Reset();
for(i=0;i<div;i++)
{
for(j=0;j<3;j++) sp[j]=((exip[j]-enp[j])/div)*i+enp[j]+(exip[j]-enp[j])/(2*div);
xs=((exip[1]-enp[1])/div)*i;
xe=xs+(exip[1]-enp[1])/div;
cd=(TMath::Exp(-fabs(xs)/lambda)-TMath::Exp(-fabs(xe)/lambda));
// printf("%f %f %f %f %fn",sp[0],sp[1],fabs(xs),fabs(xe),cd);
for(ii=0;ii<average;ii++)
{
Drift(sp[0],sp[1],1,&seg,MobMod);
seg.GetCH(histop,1);
Drift(sp[0],sp[1],-1,&seg,MobMod);
seg.GetCH(histon,1);
}
histop->Scale(1/(Float_t)average); histon->Scale(1/(Float_t)average);
// pad2->cd(); histon->Draw(); c1->Update(); printf("%d %en",i,histon->Integral());
if(trapping) {
ht->Trapping(histop,trp); outhisp=&trp[3];
et->Trapping(histon,trn); outhisn=&trn[3];
//pad2->cd(); histon->Draw(); c1->Update();
}
else {outhisp=histop; outhisn=histon;}
shiftp=(Int_t) (MaxPudu*1e-9/outhisp->GetBinWidth(1))+1;
shiftn=(Int_t) (MaxPudu*1e-9/outhisn->GetBinWidth(1))+1;
//printf("%d %dn",shiftp,shiftn);
for(k=0;k<PuduStep;k++)
{
for(j=1;j<=outhisn->GetNbinsX()-shiftn;j++) {z=(Int_t)(tt[k]/histon->GetBinWidth(1)); neg->AddBinContent(j+z,outhisn->GetBinContent(j)*cd*fabs(cc[k]));}
for(j=1;j<=outhisp->GetNbinsX()-shiftp;j++) {z=(Int_t)(tt[k]/histop->GetBinWidth(1)); pos->AddBinContent(j+z,outhisp->GetBinContent(j)*cd*fabs(cc[k]));}
}
histon->Reset(); histop->Reset();
}
sum->Add(neg);
sum->Add(pos);
delete histon;
delete histop;
delete cc;
delete tt;
delete las;
//shapper(taush,cht);
}
Double_t Dioda::dEdx(Double_t E)
{
Double_t konst=0.1535;
Double_t A=28.086; //atomic mass of the material
Double_t Z=14; //athomic number of material
Double_t rho=2.33; //density of silicon
Double_t z=2; //alpha particles;
Double_t mass=3727; //alpha particles;
Double_t me=0.511; //alpha particles;
Double_t C0=-4.44,a=0.1492,m=3.25;
Double_t X0=0.2014,X1=2.87;
Double_t C=0;
Double_t gamma=E/mass+1;
Double_t eta=TMath::Sqrt(gamma*gamma-1);
Double_t beta=TMath::Sqrt(1-1/(gamma*gamma));
Double_t Wmax=2*me*eta*eta; //me<<mass
Double_t X=TMath::Log10(eta);
Double_t delta=0;
Double_t dE;
Double_t I = (9.76 * Z + 58.8 * TMath::Power(Z,-0.19)) * 1e-6;
if(X<X0) delta=0;
if(X0<X && X<X1) delta=4.6052*X+C0+a*TMath::Power((X1-X),m);
if(X>X1) delta=4.6052*X+C0;
Double_t logarg=2*me*eta*eta*Wmax/(I*I);
// printf("logarg=%e %e %e :::::::: ",logarg,TMath::Log(logarg) , I);
dE=konst*Z/A*TMath::Power(z/beta,2)*rho*(TMath::Log(logarg)-2*beta*beta-delta-2*C/Z);
return(-dE);
}
Float_t Dioda::dEX(Double_t E,Double_t *x, Double_t *y,Double_t eps)
{
Int_t k=0;
Float_t E0=E,p=0;
Double_t xx=0,dE;
eps=eps*1e-4;
while(E>0.6)
{
dE=dEdx(E)*eps;
if(dE<0)
E+=dE; else {dE=-E; E=0;}
x[k]=xx*1e4;
xx+=eps;
y[k]=TMath::Abs(dE/E0);
// p+=y[k];
// printf("E=%f , dE=%f , x=%fn",E,y[k],x[k]);
k++;
}
y[k]=TMath::Abs(E/E0);
x[k]=xx*1e4;
printf("E=%f , dE=%f , x=%fn",0.0,y[k],x[k]);
// printf("p=%f",p+y[k]);
return (Float_t) x[k];
}
void Dioda::Alpha(Float_t E, Float_t Angle, Float_t *enp)
{
// The simulation of the drift for the minimum ionizing particles.
// A track is devided into Int_ div buckets. Each bucket is drifted in the field. The
// induced currents for each carrier is calculated as the sum all buckets.
// Int_t MobMod; mobility model
// Float_t B; magnetic field
TH1F trn[4],trp[4];
Double_t x[1000],y[1000];
Float_t sp[3];
int i,j,div;
DStruct seg;
TH1F *histop = new TH1F("ch-","charge+",pos->GetNbinsX(),pos->GetXaxis()->GetXmin(),pos->GetXaxis()->GetXmax());
TH1F *histon = new TH1F("ch+","charge-",neg->GetNbinsX(),neg->GetXaxis()->GetXmin(),neg->GetXaxis()->GetXmax());
sum->Reset(); pos->Reset(); neg->Reset();
div=(Int_t) dEX(E,x,y,1);
for(i=0;i<=div;i++)
{
for(j=0;j<3;j++) {if(j==0) sp[j]=enp[j]+TMath::Sin(Angle)*x[i]; if(j==1) sp[j]=enp[j]-TMath::Cos(Angle)*x[i];}
//if(+enp[j]+(exp[j]-enp[j])/(2*div);}
//printf("#i=%d div=%d, pointx=%f, pointy=%f weight=%fn",i,div,sp[0],sp[1],y[i]);
for(j=0;j<average;j++){
Drift(sp[0],sp[1],1,&seg,MobMod);
seg.GetCH(histop,1);
Drift(sp[0],sp[1],-1,&seg,MobMod);
seg.GetCH(histon,1);
}
histop->Scale(y[i]/((Float_t)average)); histon->Scale(y[i]/((Float_t)average));
// histop->Scale(0.05); histon->Scale(0.05);
if(trapping) {
ht->Trapping(histop,trp);
et->Trapping(histon,trn);
pos->Add(&trp[3]); neg->Add(&trn[3]);
} else {pos->Add(histop);neg->Add(histon);}
histop->Reset();
histon->Reset();
}
sum->Add(neg);
sum->Add(pos);
delete histop;
delete histon;
}
void Dioda::CalM(DStruct *seg, Double_t *data)
{
Int_t i,j,numreg=0;
Double_t *M=new Double_t [seg->Steps+1];
Double_t *DIF=new Double_t [seg->Steps+1];
Double_t dx,dy,sum=1.,sum2=0,sum1=0,dif,dif1,dif2;
printf("Number of Steps=%dn",seg->Steps);
for(i=1; i<seg->Steps;i++)
{
// if(i==seg->Steps) DIF[i]=0;
dx=TMath::Sqrt(TMath::Power((seg->Xtrack[i+1]-seg->Xtrack[i]),2)+TMath::Power((seg->Ytrack[i+1]-seg->Ytrack[i]),2));
dif=Real.alpha(0.5*(seg->Efield[i+1]+seg->Efield[i]));
sum*=(1+dif*dx);
DIF[i]=sum;
}
for(i=1;i<seg->Steps;i++) { printf("Step=%d [%f %f], E=%f , a[i]=%f M(i)=%fn",i,seg->Xtrack[i],seg->Ytrack[i],seg->Efield[i],Real.alpha(seg->Efield[i]),DIF[i]); }
data[0]=DIF[seg->Steps-1]; data[1]=0; data[2]=0; data[3]=0;
printf("KKK %fn",data[0]);
for(i=1;i<seg->Steps;i++)
{
if(DIF[i]>data[0]*0.2)
{
data[1]+=seg->Xtrack[i]; data[2]+=seg->Ytrack[i]; numreg++;
data[3]+=seg->Time[i];
printf("::::: %e %e %e %dn",data[1],data[2],data[3],numreg);
}
}
data[1]/=numreg; data[2]/=numreg; data[3]/=numreg;
// DIF[1]=0;
// for(i=1; i<seg->Steps-1;i++)
// {
// // if(i==seg->Steps) DIF[i]=0;
// dx=TMath::Sqrt(TMath::Power((seg->Xtrack[i+1]-seg->Xtrack[i]),2)+TMath::Power((seg->Ytrack[i+1]-seg->Ytrack[i]),2));
// dif=Real.alpha(0.5*(seg->Efield[i+1]+seg->Efield[i]))-Real.beta(0.5*(seg->Efield[i+1]+seg->Efield[i]));
// sum+=dif*dx;
// DIF[i+1]=sum;
// }
// for(j=1;j<seg->Steps-1;j++) printf("Dif[%d]=%fn",j,DIF[j]);
// for(j=1;j<seg->Steps-1;j++)
// {
// dy=TMath::Sqrt(TMath::Power((seg->Xtrack[j+1]-seg->Xtrack[j]),2)+TMath::Power((seg->Ytrack[j+1]-seg->Ytrack[j]),2));
// dif1=Real.beta(0.5*(seg->Efield[j+1]+seg->Efield[j]))*dy*TMath::Exp(-DIF[j+1]);
// sum1+=dif1;
// }
// printf("sum1=%fn",1/(1-sum1));
// for(j=1;j<seg->Steps-1;j++)
// M[j]=TMath::Exp(-(DIF[j]))/(1-sum1);
// for(i=1; i<seg->Steps-1;i++)
// {
// printf("Step=%d [%f %f], E=%f , M[i]=%f A(i)=%fn",i,seg->Xtrack[i],seg->Ytrack[i],seg->Efield[i],M[i],Real.alpha(seg->Efield[i]));
// }
// data[0]=M[seg->Steps]; data[1]=0; data[2]=0; data[3]=0;
// // printf("Mult=%fn",1/(1-sum1));
// for(seg->Steps;i>1;i++)
// {
// if(M[i]>2.)
// {
// data[1]+=seg->Xtrack[i]; data[2]+=seg->Ytrack[i]; numreg++;
// data[3]+=seg->Time[i];
// // printf("sum1=%f sum2=%f t=%e numreg=%d n",sum1,sum2,seg->Time[i],numreg);
// }
// }
// data[1]/=numreg; data[2]/=numreg; data[3]/=numreg;
}
ROOT page - Class index - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.