#include "Detector.h"
#define ABS(x) x>0?x:-x
#define PREDZNAK(x) x>0?1:-1
#define L1 y2[n]=1.;
#define R1 y4[n]=1.;
#define U1 y5[n]=1.;
#define D1 y6[n]=1.;
#define L2 y2[n]=2.;
#define R2 y4[n]=2.;
#define U2 y5[n]=2.;
#define D2 y6[n]=2.;
#define U0 y5[n]=0.;
#define D0 y6[n]=0.;
#define R0 y4[n]=0.;
#define L0 y2[n]=0.;
#define C1 y3[n]=-4.;
#define PI 3.1415927
#define EPS 1.0e-14
#define STEP_DET_CUR 25e-9
ClassImp(Detector)
double **a,*b,*y6,*y2,*y3,*y4,*y5;
double *dvector(long, long);
void free_dvector(double*,long,long);
/**********************************************************
snrm: Calculates the norm of vector
unsigned long n - dimension of the vector
double sx[] - components
itol - <=3 real norm , otherwise max component.
************************************************************/
double snrm(unsigned long n, double sx[], int itol)
{
unsigned long i,isamax;
double ans;
if (itol <= 3) {
ans = 0.0;
for (i=1;i<=n;i++) ans += sx[i]*sx[i];
return sqrt(ans);
} else {
isamax=1;
for (i=1;i<=n;i++) {
if (fabs(sx[i]) > fabs(sx[isamax])) isamax=i;
}
return fabs(sx[isamax]);
}
}
/* (C) Copr. 1986-92 Numerical Recipes Software &c&):)+!. */
/**********************************************************
linbcg: Solves linear eqution sparse system Ax=b
unsigned long n - number of equations
int nx - # of points in x direction
int ny - # of points in y direction
double b[] - right side vector
double x[] - solution of the system
int itol - a way to calculte norm
double tol
int itmax - max # of iterations
int *iter -
double *err -
************************************************************/
/*****************************************************
atimes : How to multiply the vector with the matrices
******************************************************/
void atimes(unsigned long n,int nx, int ny,double x[],double r[],int itrnsp)
{
// This is a function used to multuply the vectors with matrices!
// Used for calculation of the electric and ramo field
int i,j;
long k=0;
for(j=1; j<=ny; j++) /*mnozenje po stolpcu*/
for(i=1; i<=nx; i++) /*mnozenje po vrstici*/
{
k++;
if (i==1 && j==1) r[k]=y3[k]*x[k]+y4[k]*x[k+1]+y5[k]*x[k+nx];
if (i>1 && j==1) r[k]=y3[k]*x[k]+y4[k]*x[k+1]+y5[k]*x[k+nx]+y2[k]*x[k-1];
if (j>1 && j<ny) r[k]=y3[k]*x[k]+y4[k]*x[k+1]+y5[k]*x[k+nx]+y2[k]*x[k-1]+y6[k]*x[k-nx];
if (j==ny && i!=nx) r[k]=y3[k]*x[k]+y4[k]*x[k+1]+y2[k]*x[k-1]+y6[k]*x[k-nx];
if (i==nx && j==ny) r[k]=y3[k]*x[k]+y2[k]*x[k-1]+y6[k]*x[k-nx];
}
if(n!=k) printf("n Nekaj je narobe!");
return;
}
/***************************************************
Solve a very simple system of equations
with the diagonal elements to be the only ones
****************************************************/
void asolve(unsigned long n, double b[], double x[], int itrnsp)
{
//Solve a very simple system of n equations
//with the diagonal elements to be the only ones!
unsigned long i;
for(i=1;i<=n;i++) x[i]=(y3[i] != 0.0 ? b[i]/y3[i] : b[i]);
}
void linbcg(unsigned long n, int nx, int ny, double b[], double x[], int itol, double tol,
int itmax, int *iter, double *err)
{
// The main function for electric field calcualtion
void asolve(unsigned long n, double b[], double x[], int itrnsp);
void atimes(unsigned long n, int nx, int ny,double x[], double r[], int itrnsp);
double snrm(unsigned long n, double sx[], int itol);
double *dvector(long, long);
void free_dvector(double *, long, long);
void nrerror(char error_text[]);
unsigned long j;
double ak,akden,bk,bkden,bknum,bnrm,dxnrm,xnrm,zm1nrm,znrm;
double *p,*pp,*r,*rr,*z,*zz;
p=dvector(1,n);
pp=dvector(1,n);
r=dvector(1,n);
rr=dvector(1,n);
z=dvector(1,n);
zz=dvector(1,n);
*iter=0;
atimes(n,nx,ny,x,r,0);
for (j=1;j<=n;j++) {
r[j]=b[j]-r[j];
rr[j]=r[j];
}
atimes(n,nx,ny,r,rr,0); // minimal residual invariant
znrm=1.0;
if (itol == 1) bnrm=snrm(n,b,itol);
else if (itol == 2) {
asolve(n,b,z,0);
bnrm=snrm(n,z,itol);
}
else if (itol == 3 || itol == 4) {
asolve(n,b,z,0);
bnrm=snrm(n,z,itol);
asolve(n,r,z,0);
znrm=snrm(n,z,itol);
} else nrerror("illegal itol in linbcg");
asolve(n,r,z,0);
while (*iter <= itmax) {
++(*iter);
zm1nrm=znrm;
asolve(n,rr,zz,1);
for (bknum=0.0,j=1;j<=n;j++) bknum += z[j]*rr[j];
if (*iter == 1) {
for (j=1;j<=n;j++) {
p[j]=z[j];
pp[j]=zz[j];
}
}
else {
bk=bknum/bkden;
for (j=1;j<=n;j++) {
p[j]=bk*p[j]+z[j];
pp[j]=bk*pp[j]+zz[j];
}
}
bkden=bknum;
atimes(n,nx,ny,p,z,0);
for (akden=0.0,j=1;j<=n;j++) akden += z[j]*pp[j];
ak=bknum/akden;
atimes(n,nx,ny,pp,zz,1);
for (j=1;j<=n;j++) {
x[j] += ak*p[j];
r[j] -= ak*z[j];
rr[j] -= ak*zz[j];
}
asolve(n,r,z,0);
if (itol == 1 || itol == 2) {
znrm=1.0;
*err=snrm(n,r,itol)/bnrm;
} else if (itol == 3 || itol == 4) {
znrm=snrm(n,z,itol);
if (fabs(zm1nrm-znrm) > EPS*znrm) {
dxnrm=fabs(ak)*snrm(n,p,itol);
*err=znrm/fabs(zm1nrm-znrm)*dxnrm;
} else {
*err=znrm/bnrm;
continue;
}
xnrm=snrm(n,x,itol);
if (*err <= 0.5*xnrm) *err /= xnrm;
else {
*err=znrm/bnrm;
continue;
}
}
if(*iter==itmax) {printf("n Number of iterations exceeded the max. value n");
printf("iter=%4d err=%12.6fn",*iter,*err);}
if (*err <= tol) break;
}
free_dvector(p,1,n);
free_dvector(pp,1,n);
free_dvector(r,1,n);
free_dvector(rr,1,n);
free_dvector(z,1,n);
free_dvector(zz,1,n);
}
double Detector::V(int i, int j, int dowhat)
{
//Sets the boundary conditions for the voltage!
double voltage,fak,overdep=10;
int css,cse,k=0;
if (dowhat==0) {
if(j==1 || j==2) voltage=Voltage;
//if(TMath::Abs(voltage)<fdv()) printf("WARNING FDV<Volatge, Can not calculte field");}
else
voltage=0.;
}
if(dowhat==1){
if(j==1 || j==2) voltage=0;
while (StripPosition[k]!=0) k++; css=StripPosition[(k-1)/2]; cse=StripPosition[(k-1)/2+1]; //printf("%d %d %dn",css,cse,k);
if((j==ny-1 || j==ny) && i>=css && i<=cse ) voltage=1; else voltage=0;}
return voltage;
}
Double_t Detector::kappa(int i,int j, int dowhat )
{
//Sets the effective space charge values for x,y in the equation!
double y;
if (dowhat==0) {y=(Neff->Eval(i*Step,j*Step)*1e6*e_0)/(perm*perm0); /*printf("i=%d,j=%d,y=%en",i,j,y);*/}
//if (dowhat==0) if(j>nc) y=(Step*Step*1e-12)/(Si_mue*Ro*perm*perm0); else y=-(Step*Step*1e-12)/(Si_mue*Ro*perm*perm0);
else
y=0.;
return y;
}
void Detector::Declaration(int dowhat)
{
//Declaration of variables used for calculation of electric field
int i,j,l=2,k,starts=1,ends=nx,ramoth;
long n;
//double kappa(int, int, int, int, float, int);
//double V(int, int, int, int, int *,float, float, int);
//for(i=0;i<NoStrips*2;i+=2) printf("Zacetek %d: konec %d : nx=%d: ny=%dn",StripPosition[i],StripPosition[i+1],nx,ny);
for (j=1;j<=ny;j++)
for(i=1;i<=nx;i++)
{ if (j==1) {
if(i==1) {n=(j-1)*nx+i; R2 U2 C1 L0 D0 b[n]=-kappa(i,j,dowhat);}
if(i!=1 && j==1 && i!=ends+1) {n=(j-1)*nx+i; L1 R1 U2 C1 D0 b[n]=-kappa(i,j,dowhat);}
if (i>=1 && i<=ends){
n=(j-1)*nx+i;
U0 D0 L0 R0 y3[n]=1;
b[n]=V(i,j,dowhat);
if(i==starts && i>2) {n=n-1; L1 R0 D0 U2 C1 b[n]=-kappa(i,j,dowhat)-V(i,j,dowhat);}
if(i==ends && i<nx-1) {n=n+1; R1 L0 D0 U2 C1 b[n]=-kappa(i,j,dowhat)-V(i,j,dowhat);}
if(i==2 && starts==2 ) {n=n-1; L0 R0 D0 U2 C1 b[n]=-kappa(i,j,dowhat)-2*V(i,j,dowhat);}
}
if (i==nx) {n=(j-1)*nx+i; U2 L2 D0 R0 C1 b[n]=-kappa(i,j,dowhat);
if(nx==ends) {U0 D0 L0 R0 y3[n]=1; b[n]=V(i,j,dowhat);}
if(ends==nx-1) {L0 R0 U2 D0 C1 b[n]=-kappa(i,j,dowhat)-2*V(i,j,dowhat);}
}
}
if (i==1 && j!=1) { n=nx*(j-1)+1; D1 U1 L0 R2 C1 b[n]=-kappa(i,j,dowhat);}
if (i==1 && j!=1) { n=nx*(j-1)+1; D1 U1 L0 R2 C1 b[n]=-kappa(i,j,dowhat);}
if (i==nx && j!=1) { n=nx*j; D1 L2 R0 C1 U1 b[n]=-kappa(i,j,dowhat);}
if (i!=1 && i<nx && j!=1 && j<ny ) { n=(j-1)*nx+i; L1 R1 U1 D1 C1 b[n]=-kappa(i,j,dowhat);}
if (j==2 && i>=starts && i<=ends) {
n=(j-1)*nx+i; D0 L1 R1 C1 U1 b[n]=-kappa(i,j,dowhat)-V(i,j,dowhat);;
if(i==1) {D0 L0 R2 C1 U1 b[n]=-kappa(i,j,dowhat)-V(i,j,dowhat);}
if(i==nx) {D0 L2 R0 C1 U1 b[n]=-kappa(i,j,dowhat)-V(i,j,dowhat);}
}
if(j==ny-1 || j==ny)
for(k=0;k<NoStrips*2;k+=2) if(i<StripPosition[k+1]+2) {l=k; break;}// printf("stevilki: i=%d, l=%dn",i,l);}
if (j==ny-1 && i>=StripPosition[l] && i<=StripPosition[l+1]) {
//printf("stevilki: i=%d, j=%d kappa=%f V=%fn",i,j,kappa(i,j,nx,ny,StripPitch,dowhat),V(i,j,nx,ny,StripPitch,overdep,dowhat));
n=(j-1)*nx+i; D1 L1 R1 C1 U0 b[n]=-kappa(i,j,dowhat)-V(i,j,dowhat);
if(i==1) {D1 L0 R2 C1 U0 b[n]=-kappa(i,j,dowhat)-V(i,j,dowhat);}
if(i==nx) {D1 L2 R0 C1 U0 b[n]=-kappa(i,j,dowhat)-V(i,j,dowhat);}
}
if (i==1 && j==ny) {n=(j-1)*nx+i; R2 D2 C1 L0 U0 b[n]=-kappa(i,j,dowhat);}
if (i!=1 && j==ny && i!=StripPosition[l+1]+1) {n=(j-1)*nx+i; L1 R1 D2 C1 U0 b[n]=-kappa(i,j,dowhat);}
if (i>=StripPosition[l] && i<=StripPosition[l+1] && j==ny){
n=(j-1)*nx+i;
U0 D0 L0 R0 y3[n]=1;
b[n]=V(i,j,dowhat);
if(i==StripPosition[l] && i>2) {n=n-1; L1 R0 D2 U0 C1 b[n]=-kappa(i,j,dowhat)-V(i,j,dowhat);}
if(i==StripPosition[l+1] && i<nx-1) {n=n+1; R1 L0 D2 U0 C1 b[n]=-kappa(i,j,dowhat)-V(i,j,dowhat);}
if(i==2 && StripPosition[l]==2 ) {n=n-1; L0 R0 D2 U0 C1 b[n]=-kappa(i,j,dowhat)-2*V(i,j,dowhat);}
}
if (i==nx && j==ny) {n=nx*ny; D2 L2 U0 R0 C1 b[n]=-kappa(i,j,dowhat);
if(nx==StripPosition[NoStrips*2-1]) {U0 D0 L0 R0 y3[n]=1; b[n]=V(i,j,dowhat);}
if(StripPosition[NoStrips*2-1]==nx-1) {L0 R0 D2 U0 C1 b[n]=-kappa(i,j,dowhat)-2*V(i,j,dowhat);}
}
}
return;
}
void Detector::CalDimension()
{
//Calculates different geometrical parameters used in the calculation
// for example: the strip positions and the step size.
int x,y,yy,i,j;
x=(Int_t) (StripPitch*NoStrips);
nx=(int)(x/Step);
ny=(int)(Thickness/Step);
y=(int)(StripPitch/Step);
yy=(int)(StripWidth/Step);
for(i=0;i<NoStrips;i++) {j=i*2; StripPosition[j]=y/2+i*y-yy/2+1; StripPosition[j+1]=y/2+i*y+yy/2;}
return;
}
void Detector::CalPhyField()
{
//Calculates electric field
double err;
int iteracije,i;
int num=nx*ny;
Double_t *x;
b=dvector(1,num);
y6=dvector(1,num);
y2=dvector(1,num);
y3=dvector(1,num);
y4=dvector(1,num);
y5=dvector(1,num);
x=dvector(1,num);
Declaration(0);
//for(i=1;i<=num;i++) printf("%d %f %f %f %f %f %fn",i,b[i],y6[i],y2[i],y3[i],y4[i],y5[i]);
for(i=1;i<=num;i++) x[i]=1.;
linbcg(nx*ny,nx,ny,b,x,1,1e-6,2000,&iteracije,&err);
Real=EField(x,nx,ny,StripPosition);
Real.CalField(Step,1);
//Real=TArrayD(num+1,x);
free_dvector(x, 1,num);
free_dvector(b, 1,num);
free_dvector(y2, 1,num);
free_dvector(y3, 1,num);
free_dvector(y4, 1,num);
free_dvector(y5, 1,num);
free_dvector(y6, 1,num);
//for(i=1;i<=num; i++) { printf("%f ",xx[i]); if(i % nx == 0) printf("n");}
}
void Detector::CalPhyField(Double_t *x,Int_t how)
{
Real=EField(x,nx,ny,StripPosition);
Real.CalField(Step,how);
}
void Detector::CalPhyField(TH2F *uhis,Int_t how, Int_t cutx)
{
//Calculates electric field
double err;
int iteracije,i,j;
//nx=uhis->GetNbinsX()-cutx;
//ny=uhis->GetNbinsY()-cuty;
int num=nx*ny;
Double_t *x;
x=dvector(1,num);
for(i=0;i<ny;i++) for(j=0;j<nx;j++) x[j+nx*i+1]=uhis->GetBinContent(j+1+cutx,i+1);
Real=EField(x,nx,ny,StripPosition);
Real.CalField(Step,how);
//Real=TArrayD(num+1,x);
//for(i=1;i<=num; i++) { printf("%f ",x[i]); if(i % nx == 0) printf("n");}
free_dvector(x, 1,num);
}
void Detector::CalRamoField()
{
//Calculates weighting field
double err;
int iteracije,i;
int num=nx*ny;
Double_t *xr;
b=dvector(1,num);
y6=dvector(1,num);
y2=dvector(1,num);
y3=dvector(1,num);
y4=dvector(1,num);
y5=dvector(1,num);
xr=dvector(1,num);
Declaration(1);
//for(i=1;i<=num;i++) printf("%d %f %f %f %f %f %fn",i,b[i],y6[i],y2[i],y3[i],y4[i],y5[i]);
for(i=1;i<=num;i++) xr[i]=1.;
linbcg(nx*ny,nx,ny,b,xr,1,1e-6,2000,&iteracije,&err);
Ramo=EField(xr,nx,ny,StripPosition);
Ramo.CalField(Step,1);
// Ramo=TArrayD(num+1,xr);
free_dvector(xr, 1,num);
free_dvector(b, 1,num);
free_dvector(y2, 1,num);
free_dvector(y3, 1,num);
free_dvector(y4, 1,num);
free_dvector(y5, 1,num);
free_dvector(y6, 1,num);
// for(i=1;i<=num; i++) { printf("%f ",xxr[i]); if(i % nx == 0) printf("n");}
}
void Detector::CalRamoField(Double_t *x,Int_t how)
{
Ramo=EField(x,nx,ny,StripPosition);
Ramo.CalField(Step,how);
}
void Detector::CalRamoField(TH2F *uhis, Int_t how,Int_t cutx)
{
//Calculates electric field
double err;
int iteracije,i;
//nx=uhis->GetNbinsX();
//ny=uhis->GetNbinsY();
int num=nx*ny,j;
Double_t *x;
x=dvector(1,num);
for(i=0;i<ny;i++) for(j=0;j<nx;j++) x[j+nx*i+1]=uhis->GetBinContent(j+1+cutx,i+1);
Ramo=EField(x,nx,ny,StripPosition);
Ramo.CalField(Step,how);
//Real=TArrayD(num+1,x);
//for(i=1;i<=num; i++) { printf("%f ",x[i]); if(i % nx == 0) printf("n");}
free_dvector(x, 1,num);
}
void Detector::Drift(Double_t sx, Double_t sy,Float_t charg, DStruct *seg, Int_t MobMod, Float_t B, Double_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
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],tempE[3],Eprev,En[3],Erprev,Enr[3],Er[3],pathlen=0;
long st=1,i,ishit=0;
double cx,cy,vel,t=0,sumc=0,cutx=1.,cuty=1.,deltacy;
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);
seg->Efield[st]=E[0];
//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);
// Inclusion of the Magnetic field 28.8.2001, Revised 29.10.2001
if(B!=0)
{
for(Int_t cst=0; cst<3;cst++) tempE[cst]=E[cst];
if(charg>0) {E[2]=tempE[2]-muhh*tempE[1]*B; E[1]=tempE[1]+muhh*tempE[2]*B;}
else {E[2]=tempE[2]-muhe*tempE[1]*B; E[1]=tempE[1]+muhe*tempE[2]*B;}
E[0]=TMath::Sqrt(E[1]*E[1]+E[2]*E[2]);
}
// End of magnetic field inclusion 28.8.2001 , Revised 29.10.2001
while(!ishit)
{
st++;
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, Revised 29.10.2001
if(B!=0)
{
for(Int_t cst=0; cst<3;cst++) tempE[cst]=E[cst];
if(charg>0) {E[2]=tempE[2]-muhh*tempE[1]*B; E[1]=tempE[1]+muhh*tempE[2]*B;}
else {E[2]=tempE[2]-muhe*tempE[1]*B; E[1]=tempE[1]+muhe*tempE[2]*B;}
E[0]=TMath::Sqrt(E[1]*E[1]+E[2]*E[2]);
}
// End of magnetic field inclusion 28.8.2001, Revised 29.10.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; // Changedd 29.10.2001
//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];
//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;
}
(*seg).Xlenght=pathlen; (*seg).Ylenght=pathlen;
(*seg).TTime=t; (*seg).TCharge=sumc; (*seg).Steps=st;
return;
#undef x
}
TH2F *Detector::Draw(Char_t *option)
{
//The function draws weighting and electric field and also the event display
//Char_t *option:
// W -weighting
// E -electric
// P - potential
// F - |field|
// X - x component of the field
// X - y component of the field
// EVENT - event display
Int_t i=0,What=0,Which=1;
TH2F *histo=new TH2F();
if(!strcmp(option,"EVENT"))
{
Real.GetHisto(histo,1); histo->DrawCopy("COLZ");
TLine *line=new TLine(enp[0],enp[1],exp[0],exp[1]);
line->SetLineWidth(2);
line->Draw("SAME");
}
else
{
for(i=0;i<strlen(option);i++)
{
if(option[i]=='W') What=2;
if(option[i]=='E') What=1;
if(option[i]=='P') Which=0;
if(option[i]=='F') Which=1;
if(option[i]=='X') Which=2;
if(option[i]=='Y') Which=3;
}
switch(What)
{
case 1: Real.GetHisto(histo,Which); histo->DrawCopy("SURF2"); break;
case 2: Ramo.GetHisto(histo,Which); histo->DrawCopy("SURF2"); break;
case 0: printf("Can't get histogramn"); return NULL; break;
}
}
return histo;
}
void Detector::MipIR(Int_t div, Int_t MobMod, Float_t B)
{
// 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()); //2.23 -> 3.0
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,B);
seg.GetCH(histop,1);
Drift(sp[0],sp[1],-1,&seg,MobMod,B);
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]);
if(data[0]>1)
{
Drift(data[1],data[2],1,&seg,MobMod,B,data[3]);
seg.GetCH(histop,1,data[0]-1);
Drift(data[1],data[2],-1,&seg,MobMod,B,data[3]);
seg.GetCH(histon,1,data[0]-1);
}
}
}
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;
}
void Detector::ShowMipIR(Int_t div, Int_t MobMod, Float_t B, Int_t color,Int_t how)
{
// 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];
TGraph *gr;
Float_t sp[3];
int i,j;
DStruct seg;
TLine *line;
if(color<2)
{
TH2F *pot = new TH2F("Drift path","pot",nx,0,nx,ny,0,ny);
pot->SetLineColor(1);
pot->SetXTitle("x [ #mum ]");
pot->SetYTitle("y [ #mum ]");
pot->GetXaxis()->SetTitleOffset(1.3);
pot->GetYaxis()->SetTitleOffset(1.3);
pot->SetLabelSize(0.05,"X");
pot->SetLabelSize(0.05,"Y");
pot->Draw();
}
for(i=0;i<NoStrips*2;i+=2)
{
line=new TLine(StripPosition[i],ny,StripPosition[i+1],ny);
line->SetLineWidth(5);
line->Draw("SAME");
}
if(how)
for(i=0;i<NoStrips*2;i+=2)
{
line=new TLine(StripPosition[i],1,StripPosition[i+1],1);
line->SetLineWidth(5);
line->Draw("SAME");
}
for(i=0;i<NoStrips-1;i++)
{
line=new TLine((StripPosition[i*2+1]+StripPosition[i*2+2])/2.,ny,(StripPosition[i*2+1]+StripPosition[i*2+2])/2.,1);
line->SetLineStyle(3);
line->Draw("SAME");
}
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]);
Drift(sp[0],sp[1],1,&seg,MobMod,B);
gr=new TGraph(seg.Steps,&seg.Xtrack[1],&seg.Ytrack[1]);
if(color==1 || color==3) gr->SetLineColor(2); //else
gr->SetLineStyle(1); gr->Draw("L");
// seg.GetGraph(gr); gr->Draw("L");
Drift(sp[0],sp[1],-1,&seg,MobMod,B);
gr=new TGraph(seg.Steps,&seg.Xtrack[1],&seg.Ytrack[1]);
if(color==1 || color==3) gr->SetLineColor(4); //else
gr->SetLineStyle(3); gr->Draw("L");
// seg.GetGraph(gr); gr->Draw("L");
}
}
Detector::Detector(Float_t x0,TF2 *neff,Float_t x1,Float_t x2,Int_t x3,Float_t x4,Float_t x5)
{
//Constructor:
// Float_t x0; Voltage
// TF2 *neff; function describing the distrbution of the effective space cahrge
// Float_t x1; Strip Pitch
// Float_t x2; Strip Width
// Int_t x3; Number of strips
// Float_t x4; Detector thickness
// Float_t x5; Step-size of the mesh for drift and field calculation (defautl = 1um)
Voltage=x0;
StripPitch=x1;
StripWidth=x2;
NoStrips=x3;
Thickness=x4;
Step=x5;
Neff=neff; //printf("TEST=%e %en",neff->Eval(2,2),Neff->Eval(2,2));
StripPosition=TArrayI(NoStrips*2+1);
ran=new TRandom(33);
pos = new TH1F("charge+","Positive Charge",200,0,STEP_DET_CUR);
neg = new TH1F("charge-","Negative Charge",200,0,STEP_DET_CUR);
sum = new TH1F("charge","Total Charge",200,0,STEP_DET_CUR);
enp[0]=((Float_t)NoStrips)/2*StripPitch; enp[1]=1; enp[2]=0;
exp[0]=enp[0]; exp[1]=Thickness; exp[2]=0;
trapping=0;
average=1;
diff=0;
Temperature=263;
Multiplication=0;
CalDimension();
}
Detector::Detector()
{
//Constructor: with the set of parameters (100 , Neff , 50, 10, 3, 301 , 1);
// Neff=zero space charge
Neff=new TF2("Profile","[0]",0,1000,0,1000);
Neff->SetParameter(0,1);
Detector( 100 , Neff , 80, 18, 9, 281 , 1);
}
Detector::~Detector()
{
//destructor
Clear();
}
void Detector::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]-1)/(data[0]-1)>0.8 && DIF[i]>1.02)
{
data[1]+=seg->Xtrack[i]; data[2]+=seg->Ytrack[i];
data[3]+=seg->Time[i]; numreg++;
// printf("::::: %e %e %e %dn",data[1],data[2],data[3],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.