TPC/AliTPCExBExact.cxx

   1 #include "TMath.h"
   2 #include "TTreeStream.h"
   3 #include "AliTPCExBExact.h"
   4
   5 ClassImp(AliTPCExBExact)
   6
   7 const Double_t AliTPCExBExact::fgkEM=1.602176487e-19/9.10938215e-31;
   8 const Double_t AliTPCExBExact::fgkDriftField=40.e3;
   9
  10 AliTPCExBExact::AliTPCExBExact(const AliMagF *bField,
  11                                Double_t driftVelocity,
  12                                Int_t nx,Int_t ny,Int_t nz,Int_t n)
  13   : fkMap(0),fkField(bField),fkN(n),
  14     fkNX(nx),fkNY(ny),fkNZ(nz),
  15     fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
  16     fkZMax(250.) {
  17   //
  18   // The constructor. One has to supply a magnetic field and an (initial)
  19   // drift velocity. Since some kind of lookuptable is created the
  20   // number of its meshpoints can be supplied.
  21   // n sets the number of integration steps to be used when integrating
  22   // over the full drift length.
  23   //
  24   fDriftVelocity=driftVelocity;
  25   CreateLookupTable();
  26 }
  27
  28 AliTPCExBExact::AliTPCExBExact(const AliFieldMap *bFieldMap,
  29                                Double_t driftVelocity,Int_t n)
  30   : fkMap(bFieldMap),fkField(0),fkN(n) {
  31   //
  32   // The constructor. One has to supply a field map and an (initial)
  33   // drift velocity.
  34   // n sets the number of integration steps to be used when integrating
  35   // over the full drift length.
  36   //
  37   fDriftVelocity=driftVelocity;
  38
  39   fkXMin=bFieldMap->Xmin()
  40     -TMath::Ceil( (bFieldMap->Xmin()+250.0)/bFieldMap->DelX())
  41     *bFieldMap->DelX();
  42   fkXMax=bFieldMap->Xmax()
  43     -TMath::Floor((bFieldMap->Xmax()-250.0)/bFieldMap->DelX())
  44     *bFieldMap->DelX();
  45   fkYMin=bFieldMap->Ymin()
  46     -TMath::Ceil( (bFieldMap->Ymin()+250.0)/bFieldMap->DelY())
  47     *bFieldMap->DelY();
  48   fkYMax=bFieldMap->Ymax()
  49     -TMath::Floor((bFieldMap->Ymax()-250.0)/bFieldMap->DelY())
  50     *bFieldMap->DelY();
  51   fkZMax=bFieldMap->Zmax()
  52     -TMath::Floor((bFieldMap->Zmax()-250.0)/bFieldMap->DelZ())
  53     *bFieldMap->DelZ();
  54   fkZMax=TMath::Max(0.,fkZMax); // I really hope that this is unnecessary!
  55
  56   fkNX=static_cast<Int_t>((fkXMax-fkXMin)/bFieldMap->DelX()+1.1);
  57   fkNY=static_cast<Int_t>((fkYMax-fkYMin)/bFieldMap->DelY()+1.1);
  58   fkNZ=static_cast<Int_t>((fkZMax-fkZMin)/bFieldMap->DelZ()+1.1);
  59
  60   CreateLookupTable();
  61 }
  62
  63 AliTPCExBExact::~AliTPCExBExact() {
  64   //
  65   // destruct the poor object.
  66   //
  67   delete[] fLook;
  68 }
  69
  70 void AliTPCExBExact::Correct(const Double_t *position,Double_t *corrected) {
  71   Double_t r=TMath::Sqrt(position[0]*position[0]+position[1]*position[1]);
  72   if (TMath::Abs(position[2])>250.||r<90.||250.<r) {
  73     for (Int_t i=0;i<3;++i) corrected[i]=position[i];
  74   }
  75   else {
  76     Double_t x=(position[0]-fkXMin)/(fkXMax-fkXMin)*(fkNX-1);
  77     Int_t xi1=static_cast<Int_t>(x);
  78     xi1=TMath::Max(TMath::Min(xi1,fkNX-2),0);
  79     Int_t xi2=xi1+1;
  80     Double_t dx=(x-xi1);
  81     Double_t dx1=(xi2-x);
  82
  83     Double_t y=(position[1]-fkYMin)/(fkYMax-fkYMin)*(fkNY-1);
  84     Int_t yi1=static_cast<Int_t>(y);
  85     yi1=TMath::Max(TMath::Min(yi1,fkNY-2),0);
  86     Int_t yi2=yi1+1;
  87     Double_t dy=(y-yi1);
  88     Double_t dy1=(yi2-y);
  89
  90     Double_t z=position[2]/fkZMax*(fkNZ-1);
  91     Int_t side;
  92     if (z>0) {
  93       side=1;
  94     }
  95     else {
  96       z=-z;
  97       side=0;
  98     }
  99     Int_t zi1=static_cast<Int_t>(z);
 100     zi1=TMath::Max(TMath::Min(zi1,fkNZ-2),0);
 101     Int_t zi2=zi1+1;
 102     Double_t dz=(z-zi1);
 103     Double_t dz1=(zi2-z);
 104
 105     for (int i=0;i<3;++i)
 106       corrected[i]
 107         =fLook[(((xi1*fkNY+yi1)*fkNZ+zi1)*2+side)*3+i]*dx1*dy1*dz1
 108         +fLook[(((xi1*fkNY+yi1)*fkNZ+zi2)*2+side)*3+i]*dx1*dy1*dz
 109         +fLook[(((xi1*fkNY+yi2)*fkNZ+zi1)*2+side)*3+i]*dx1*dy *dz1
 110         +fLook[(((xi1*fkNY+yi2)*fkNZ+zi2)*2+side)*3+i]*dx1*dy *dz
 111         +fLook[(((xi2*fkNY+yi2)*fkNZ+zi1)*2+side)*3+i]*dx *dy *dz1
 112         +fLook[(((xi2*fkNY+yi2)*fkNZ+zi2)*2+side)*3+i]*dx *dy *dz
 113         +fLook[(((xi2*fkNY+yi1)*fkNZ+zi1)*2+side)*3+i]*dx *dy1*dz1
 114         +fLook[(((xi2*fkNY+yi1)*fkNZ+zi2)*2+side)*3+i]*dx *dy1*dz ;
 115     //    corrected[2]=position[2];
 116   }
 117 }
 118
 119 void AliTPCExBExact::TestThisBeautifulObject(const char* fileName) {
 120   //
 121   // well, as the name sais...
 122   //
 123   TTreeSRedirector ts(fileName);
 124   Double_t x[3];
 125   for (x[0]=-250.;x[0]<=250.;x[0]+=10.)
 126     for (x[1]=-250.;x[1]<=250.;x[1]+=10.)
 127       for (x[2]=-250.;x[2]<=250.;x[2]+=10.) {
 128         Double_t d[3];
 129         Double_t dnl[3];
 130         Correct(x,d);
 131         CalculateDistortion(x,dnl);
 132         Double_t r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
 133         Double_t rd=TMath::Sqrt(d[0]*d[0]+d[1]*d[1]);
 134         Double_t dr=r-rd;
 135         Double_t phi=TMath::ATan2(x[0],x[1]);
 136         Double_t phid=TMath::ATan2(d[0],d[1]);
 137         Double_t dphi=phi-phid;
 138         if (dphi<0.) dphi+=TMath::TwoPi();
 139         if (dphi>TMath::Pi()) dphi=TMath::TwoPi()-dphi;
 140         Double_t drphi=r*dphi;
 141         Double_t dx=x[0]-d[0];
 142         Double_t dy=x[1]-d[1];
 143         Double_t dz=x[2]-d[2];
 144         Double_t dnlx=x[0]-dnl[0];
 145         Double_t dnly=x[1]-dnl[1];
 146         Double_t dnlz=x[2]-dnl[2];
 147         ts<<"positions"
 148           <<"x0="<<x[0]
 149           <<"x1="<<x[1]
 150           <<"x2="<<x[2]
 151           <<"dx="<<dx
 152           <<"dy="<<dy
 153           <<"dz="<<dz
 154           <<"dnlx="<<dnlx
 155           <<"dnly="<<dnly
 156           <<"dnlz="<<dnlz
 157           <<"r="<<r
 158           <<"phi="<<phi
 159           <<"dr="<<dr
 160           <<"drphi="<<drphi
 161           <<"\n";
 162       }
 163 }
 164
 165 void AliTPCExBExact::CreateLookupTable() {
 166   //
 167   // Helper function to fill the lookup table.
 168   //
 169   fLook=new Double_t[fkNX*fkNY*fkNZ*2*3];
 170   Double_t x[3];
 171   for (int i=0;i<fkNX;++i) {
 172     x[0]=fkXMin+(fkXMax-fkXMin)/(fkNX-1)*i;
 173     for (int j=0;j<fkNY;++j) {
 174       x[1]=fkYMin+(fkYMax-fkYMin)/(fkNY-1)*j;
 175       for (int k=0;k<fkNZ;++k) {
 176         x[2]=1.*fkZMax/(fkNZ-1)*k;
 177         x[2]=TMath::Max((Double_t)0.0001,x[2]); //ugly
 178         CalculateDistortion(x,&fLook[(((i*fkNY+j)*fkNZ+k)*2+1)*3]);
 179         x[2]=-x[2];
 180         CalculateDistortion(x,&fLook[(((i*fkNY+j)*fkNZ+k)*2+0)*3]);
 181       }
 182     }
 183   }
 184 }
 185
 186 void AliTPCExBExact::GetE(Double_t *E,const Double_t *x) const {
 187   //
 188   // Helper function returning the E field in SI units (V/m).
 189   //
 190   E[0]=0.;
 191   E[1]=0.;
 192   E[2]=(x[2]<0.?-1.:1.)*fgkDriftField; // in V/m
 193 }
 194
 195 void AliTPCExBExact::GetB(Double_t *B,const Double_t *x) const {
 196   //
 197   // Helper function returning the B field in SI units (T).
 198   //
 199   Float_t xm[3];
 200   // the beautiful m to cm (and the ugly "const_cast") and Double_t
 201   // to Float_t read the NRs introduction!:
 202   for (int i=0;i<3;++i) xm[i]=x[i]*100.;
 203   Float_t Bf[3];
 204   if (fkMap!=0)
 205     fkMap->Field(xm,Bf);
 206   else
 207     fkField->Field(xm,Bf);
 208   for (int i=0;i<3;++i) B[i]=Bf[i]/10.;
 209 }
 210
 211 void AliTPCExBExact::Motion(const Double_t *x,Double_t,
 212                             Double_t *dxdt) const {
 213   //
 214   // The differential equation of motion of the electrons.
 215   //
 216   const Double_t tau=fDriftVelocity/fgkDriftField/fgkEM;
 217   const Double_t tau2=tau*tau;
 218   Double_t E[3];
 219   Double_t B[3];
 220   GetE(E,x);
 221   GetB(B,x);
 222   Double_t wx=fgkEM*B[0];
 223   Double_t wy=fgkEM*B[1];
 224   Double_t wz=fgkEM*B[2];
 225   Double_t ex=fgkEM*E[0];
 226   Double_t ey=fgkEM*E[1];
 227   Double_t ez=fgkEM*E[2];
 228   Double_t w2=(wx*wx+wy*wy+wz*wz);
 229   dxdt[0]=(1.+wx*wx*tau2)*ex+(wz*tau+wx*wy*tau2)*ey+(-wy*tau+wx*wz*tau2)*ez;
 230   dxdt[1]=(-wz*tau+wx*wy*tau2)*ex+(1.+wy*wy*tau2)*ey+(wx*tau+wy*wz*tau2)*ez;
 231   dxdt[2]=(wy*tau+wx*wz*tau2)*ex+(-wx*tau+wy*wz*tau2)*ey+(1.+wz*wz*tau2)*ez;
 232   Double_t fac=tau/(1.+w2*tau2);
 233   dxdt[0]*=fac;
 234   dxdt[1]*=fac;
 235   dxdt[2]*=fac;
 236 }
 237
 238 void AliTPCExBExact::CalculateDistortion(const Double_t *x0,
 239                                          Double_t *dist) const {
 240   //
 241   // Helper function that calculates one distortion by integration
 242   // (only used to fill the lookup table).
 243   //
 244   const Double_t h=0.01*250./fDriftVelocity/fkN;
 245   Double_t t=0.;
 246   Double_t xt[3];
 247   Double_t xo[3];
 248   for (int i=0;i<3;++i)
 249     xo[i]=xt[i]=x0[i]*0.01;
 250   while (TMath::Abs(xt[2])<250.*0.01) {
 251     for (int i=0;i<3;++i)
 252       xo[i]=xt[i];
 253     DGLStep(xt,t,h);
 254     t+=h;
 255   }
 256   if (t!=0.) {
 257     Double_t p=((xt[2]<0.?-1.:1.)*250.*0.01-xo[2])/(xt[2]-xo[2]);
 258     dist[0]=(xo[0]+p*(xt[0]-xo[0]))*100.;
 259     dist[1]=(xo[1]+p*(xt[1]-xo[1]))*100.;
 260     //    dist[2]=(xo[2]+p*(xt[2]-xo[2]))*100.;
 261     dist[2]=(x0[2]>0.?-1:1.)*(t-h+p*h)*fDriftVelocity*100.;
 262     dist[2]+=(x0[2]<0.?-1:1.)*250.;
 263   }
 264   else {
 265     dist[0]=x0[0];
 266     dist[1]=x0[1];
 267     dist[2]=x0[2];
 268   }
 269   // reverse the distortion, i.e. get the correction
 270   dist[0]=x0[0]-(dist[0]-x0[0]);
 271   dist[1]=x0[1]-(dist[1]-x0[1]);
 272 }
 273
 274 void AliTPCExBExact::DGLStep(Double_t *x,Double_t t,Double_t h) const {
 275   //
 276   // An elementary integration step.
 277   // (simple Euler Method)
 278   //
 279   Double_t dxdt[3];
 280   Motion(x,t,dxdt);
 281   for (int i=0;i<3;++i)
 282     x[i]+=h*dxdt[i];
 283
 284   /* suggestions about how to write it this way are welcome!
 285      void DGLStep(void (*f)(const Double_t *x,Double_t t,Double_t *dxdt),
 286                    Double_t *x,Double_t t,Double_t h,Int_t n) const;
 287   */
 288
 289 }