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