TPC/AliTPCExBExact.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 ////
  17 // This implementation AliTPCExB is using an "exact" calculation of the ExB
  18 // effect. That is, it uses the drift ODE to calculate the distortion
  19 // without any further assumption.
  20 // Due to the numerical integration of the ODE, there are some numerical
  21 // uncertencies involed.
  22 ////
  23
  24 #include "TMath.h"
  25 #include "TTreeStream.h"
  26 #include "AliMagF.h"
  27 #include "AliTPCExBExact.h"
  28
  29 ClassImp(AliTPCExBExact)
  30
  31 const Double_t AliTPCExBExact::fgkEM=1.602176487e-19/9.10938215e-31;
  32 const Double_t AliTPCExBExact::fgkDriftField=-40.e3;
  33
  34 AliTPCExBExact::AliTPCExBExact()
  35   : fDriftVelocity(0),
  36     //fkMap(0),
  37     fkField(0),fkN(0),
  38     fkNX(0),fkNY(0),fkNZ(0),
  39     fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
  40     fkZMin(-250.),fkZMax(250.),
  41     fkNLook(0),fkLook(0) {
  42   //
  43   // purely for I/O
  44   //
  45 }
  46
  47 AliTPCExBExact::AliTPCExBExact(const AliMagF *bField,
  48                                Double_t driftVelocity,
  49                                Int_t nx,Int_t ny,Int_t nz,Int_t n)
  50   : fDriftVelocity(driftVelocity),
  51     //fkMap(0),
  52     fkField(bField),fkN(n),
  53     fkNX(nx),fkNY(ny),fkNZ(nz),
  54     fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
  55     fkZMin(-250.),fkZMax(250.),
  56     fkNLook(0),fkLook(0) {
  57   //
  58   // The constructor. One has to supply a magnetic field and an (initial)
  59   // drift velocity. Since some kind of lookuptable is created the
  60   // number of its meshpoints can be supplied.
  61   // n sets the number of integration steps to be used when integrating
  62   // over the full drift length.
  63   //
  64   CreateLookupTable();
  65 }
  66
  67 /*
  68 AliTPCExBExact::AliTPCExBExact(const AliFieldMap *bFieldMap,
  69                                Double_t driftVelocity,Int_t n)
  70   : fDriftVelocity(driftVelocity),
  71     fkMap(bFieldMap),fkField(0),fkN(n),
  72     fkNX(0),fkNY(0),fkNZ(0),
  73     fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
  74     fkZMin(-250.),fkZMax(250.),
  75     fkNLook(0),fkLook(0) {
  76   //
  77   // The constructor. One has to supply a field map and an (initial)
  78   // drift velocity.
  79   // n sets the number of integration steps to be used when integrating
  80   // over the full drift length.
  81   //
  82
  83   fkXMin=bFieldMap->Xmin()
  84     -TMath::Ceil( (bFieldMap->Xmin()+250.0)/bFieldMap->DelX())
  85     *bFieldMap->DelX();
  86   fkXMax=bFieldMap->Xmax()
  87     -TMath::Floor((bFieldMap->Xmax()-250.0)/bFieldMap->DelX())
  88     *bFieldMap->DelX();
  89   fkYMin=bFieldMap->Ymin()
  90     -TMath::Ceil( (bFieldMap->Ymin()+250.0)/bFieldMap->DelY())
  91     *bFieldMap->DelY();
  92   fkYMax=bFieldMap->Ymax()
  93     -TMath::Floor((bFieldMap->Ymax()-250.0)/bFieldMap->DelY())
  94     *bFieldMap->DelY();
  95   fkZMax=bFieldMap->Zmax()
  96     -TMath::Floor((bFieldMap->Zmax()-250.0)/bFieldMap->DelZ())
  97     *bFieldMap->DelZ();
  98   fkZMax=TMath::Max(0.,fkZMax); // I really hope that this is unnecessary!
  99
 100   fkNX=static_cast<Int_t>((fkXMax-fkXMin)/bFieldMap->DelX()+1.1);
 101   fkNY=static_cast<Int_t>((fkYMax-fkYMin)/bFieldMap->DelY()+1.1);
 102   fkNZ=static_cast<Int_t>((fkZMax-fkZMin)/bFieldMap->DelZ()+1.1);
 103
 104   CreateLookupTable();
 105 }
 106 */
 107
 108 AliTPCExBExact::~AliTPCExBExact() {
 109   //
 110   // destruct the poor object.
 111   //
 112   delete[] fkLook;
 113 }
 114
 115 void AliTPCExBExact::Correct(const Double_t *position, Double_t *corrected) {
 116   //
 117   // correct for the distortion
 118   //
 119   Double_t x=(position[0]-fkXMin)/(fkXMax-fkXMin)*(fkNX-1);
 120   Int_t xi1=static_cast<Int_t>(x);
 121   xi1=TMath::Max(TMath::Min(xi1,fkNX-2),0);
 122   Int_t xi2=xi1+1;
 123   Double_t dx=(x-xi1);
 124   Double_t dx1=(xi2-x);
 125
 126   Double_t y=(position[1]-fkYMin)/(fkYMax-fkYMin)*(fkNY-1);
 127   Int_t yi1=static_cast<Int_t>(y);
 128   yi1=TMath::Max(TMath::Min(yi1,fkNY-2),0);
 129   Int_t yi2=yi1+1;
 130   Double_t dy=(y-yi1);
 131   Double_t dy1=(yi2-y);
 132
 133   Double_t z=position[2]/fkZMax*(fkNZ-1);
 134   Int_t side;
 135   if (z>0) {
 136     side=1;
 137   }
 138   else {
 139     z=-z;
 140     side=0;
 141   }
 142   Int_t zi1=static_cast<Int_t>(z);
 143   zi1=TMath::Max(TMath::Min(zi1,fkNZ-2),0);
 144   Int_t zi2=zi1+1;
 145   Double_t dz=(z-zi1);
 146   Double_t dz1=(zi2-z);
 147
 148   for (int i=0;i<3;++i)
 149     corrected[i]
 150       =fkLook[(((xi1*fkNY+yi1)*fkNZ+zi1)*2+side)*3+i]*dx1*dy1*dz1
 151       +fkLook[(((xi1*fkNY+yi1)*fkNZ+zi2)*2+side)*3+i]*dx1*dy1*dz
 152       +fkLook[(((xi1*fkNY+yi2)*fkNZ+zi1)*2+side)*3+i]*dx1*dy *dz1
 153       +fkLook[(((xi1*fkNY+yi2)*fkNZ+zi2)*2+side)*3+i]*dx1*dy *dz
 154       +fkLook[(((xi2*fkNY+yi2)*fkNZ+zi1)*2+side)*3+i]*dx *dy *dz1
 155       +fkLook[(((xi2*fkNY+yi2)*fkNZ+zi2)*2+side)*3+i]*dx *dy *dz
 156       +fkLook[(((xi2*fkNY+yi1)*fkNZ+zi1)*2+side)*3+i]*dx *dy1*dz1
 157       +fkLook[(((xi2*fkNY+yi1)*fkNZ+zi2)*2+side)*3+i]*dx *dy1*dz ;
 158   //    corrected[2]=position[2];
 159 }
 160
 161 /*
 162 void AliTPCExBExact::TestThisBeautifulObject(const AliFieldMap *bFieldMap,
 163                                              const char* fileName) {
 164   //
 165   // Have a look at the common part "TestThisBeautifulObjectGeneric".
 166   //
 167   fkMap=bFieldMap;
 168   fkField=0;
 169   TestThisBeautifulObjectGeneric(fileName);
 170 }
 171 */
 172
 173 void AliTPCExBExact::TestThisBeautifulObject(const AliMagF *bField,
 174                                              const char* fileName) {
 175   //
 176   // Have a look at the common part "TestThisBeautifulObjectGeneric".
 177   //
 178   fkField=bField;
 179   //fkMap=0;
 180   TestThisBeautifulObjectGeneric(fileName);
 181 }
 182
 183 void AliTPCExBExact::TestThisBeautifulObjectGeneric(const char* fileName) {
 184   //
 185   // Well, as the name sais... it tests the object.
 186   //
 187   TTreeSRedirector ts(fileName);
 188   Double_t x[3];
 189   for (x[0]=-250.;x[0]<=250.;x[0]+=10.)
 190     for (x[1]=-250.;x[1]<=250.;x[1]+=10.)
 191       for (x[2]=-250.;x[2]<=250.;x[2]+=10.) {
 192         Double_t d[3];
 193         Double_t dnl[3];
 194         Correct(x,d);
 195         CalculateDistortion(x,dnl);
 196         Double_t r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
 197         Double_t rd=TMath::Sqrt(d[0]*d[0]+d[1]*d[1]);
 198         Double_t dr=r-rd;
 199         Double_t phi=TMath::ATan2(x[0],x[1]);
 200         Double_t phid=TMath::ATan2(d[0],d[1]);
 201         Double_t dphi=phi-phid;
 202         if (dphi<0.) dphi+=TMath::TwoPi();
 203         if (dphi>TMath::Pi()) dphi=TMath::TwoPi()-dphi;
 204         Double_t drphi=r*dphi;
 205         Double_t dx=x[0]-d[0];
 206         Double_t dy=x[1]-d[1];
 207         Double_t dz=x[2]-d[2];
 208         Double_t dnlx=x[0]-dnl[0];
 209         Double_t dnly=x[1]-dnl[1];
 210         Double_t dnlz=x[2]-dnl[2];
 211         Double_t b[3];
 212         GetB(b,x);
 213         ts<<"positions"
 214           <<"bx="<<b[0]
 215           <<"by="<<b[1]
 216           <<"bz="<<b[2]
 217           <<"x0="<<x[0]
 218           <<"x1="<<x[1]
 219           <<"x2="<<x[2]
 220           <<"dx="<<dx
 221           <<"dy="<<dy
 222           <<"dz="<<dz
 223           <<"dnlx="<<dnlx
 224           <<"dnly="<<dnly
 225           <<"dnlz="<<dnlz
 226           <<"r="<<r
 227           <<"phi="<<phi
 228           <<"dr="<<dr
 229           <<"drphi="<<drphi
 230           <<"\n";
 231       }
 232 }
 233
 234 void AliTPCExBExact::CreateLookupTable() {
 235   //
 236   // Helper function to fill the lookup table.
 237   //
 238   fkNLook=fkNX*fkNY*fkNZ*2*3;
 239   fkLook=new Double_t[fkNLook];
 240   Double_t x[3];
 241   for (int i=0;i<fkNX;++i) {
 242     x[0]=fkXMin+(fkXMax-fkXMin)/(fkNX-1)*i;
 243     for (int j=0;j<fkNY;++j) {
 244       x[1]=fkYMin+(fkYMax-fkYMin)/(fkNY-1)*j;
 245       for (int k=0;k<fkNZ;++k) {
 246         x[2]=1.*fkZMax/(fkNZ-1)*k;
 247         x[2]=TMath::Max((Double_t)0.0001,x[2]); //ugly
 248         CalculateDistortion(x,&fkLook[(((i*fkNY+j)*fkNZ+k)*2+1)*3]);
 249         x[2]=-x[2];
 250         CalculateDistortion(x,&fkLook[(((i*fkNY+j)*fkNZ+k)*2+0)*3]);
 251       }
 252     }
 253   }
 254 }
 255
 256 void AliTPCExBExact::GetE(Double_t *e,const Double_t *x) const {
 257   //
 258   // Helper function returning the E field in SI units (V/m).
 259   //
 260   e[0]=0.;
 261   e[1]=0.;
 262   e[2]=(x[2]<0.?-1.:1.)*fgkDriftField; // in V/m
 263 }
 264
 265 void AliTPCExBExact::GetB(Double_t *b,const Double_t *x) const {
 266   //
 267   // Helper function returning the B field in SI units (T).
 268   //
 269   Double_t xm[3];
 270   // the beautiful m to cm (and the ugly "const_cast") and Double_t
 271   // to Float_t read the NRs introduction!:
 272   for (int i=0;i<3;++i) xm[i]=x[i]*100.;
 273   Double_t bf[3];
 274   //if (fkMap!=0)
 275   //  fkMap->Field(xm,bf);
 276   //else
 277   ((AliMagF*)fkField)->Field(xm,bf);
 278   for (int i=0;i<3;++i) b[i]=bf[i]/10.;
 279 }
 280
 281 void AliTPCExBExact::Motion(const Double_t *x,Double_t,
 282                             Double_t *dxdt) const {
 283   //
 284   // The differential equation of motion of the electrons.
 285   //
 286   Double_t tau=fDriftVelocity/fgkDriftField/fgkEM;
 287   Double_t tau2=tau*tau;
 288   Double_t e[3];
 289   Double_t b[3];
 290   GetE(e,x);
 291   GetB(b,x);
 292   Double_t wx=fgkEM*b[0];
 293   Double_t wy=fgkEM*b[1];
 294   Double_t wz=fgkEM*b[2];
 295   Double_t ex=fgkEM*e[0];
 296   Double_t ey=fgkEM*e[1];
 297   Double_t ez=fgkEM*e[2];
 298   Double_t w2=(wx*wx+wy*wy+wz*wz);
 299   dxdt[0]=(1.+wx*wx*tau2)*ex+(wz*tau+wx*wy*tau2)*ey+(-wy*tau+wx*wz*tau2)*ez;
 300   dxdt[1]=(-wz*tau+wx*wy*tau2)*ex+(1.+wy*wy*tau2)*ey+(wx*tau+wy*wz*tau2)*ez;
 301   dxdt[2]=(wy*tau+wx*wz*tau2)*ex+(-wx*tau+wy*wz*tau2)*ey+(1.+wz*wz*tau2)*ez;
 302   Double_t fac=tau/(1.+w2*tau2);
 303   dxdt[0]*=fac;
 304   dxdt[1]*=fac;
 305   dxdt[2]*=fac;
 306 }
 307
 308 void AliTPCExBExact::CalculateDistortion(const Double_t *x0,
 309                                          Double_t *dist) const {
 310   //
 311   // Helper function that calculates one distortion by integration
 312   // (only used to fill the lookup table).
 313   //
 314   Double_t h=0.01*250./fDriftVelocity/fkN;
 315   Double_t t=0.;
 316   Double_t xt[3];
 317   Double_t xo[3];
 318   for (int i=0;i<3;++i)
 319     xo[i]=xt[i]=x0[i]*0.01;
 320   while (TMath::Abs(xt[2])<250.*0.01) {
 321     for (int i=0;i<3;++i)
 322       xo[i]=xt[i];
 323     DGLStep(xt,t,h);
 324     t+=h;
 325   }
 326   if (t!=0.) {
 327     Double_t p=((xt[2]<0.?-1.:1.)*250.*0.01-xo[2])/(xt[2]-xo[2]);
 328     dist[0]=(xo[0]+p*(xt[0]-xo[0]))*100.;
 329     dist[1]=(xo[1]+p*(xt[1]-xo[1]))*100.;
 330     //    dist[2]=(xo[2]+p*(xt[2]-xo[2]))*100.;
 331     dist[2]=(x0[2]>0.?-1:1.)*(t-h+p*h)*fDriftVelocity*100.;
 332     dist[2]+=(x0[2]<0.?-1:1.)*250.;
 333   }
 334   else {
 335     dist[0]=x0[0];
 336     dist[1]=x0[1];
 337     dist[2]=x0[2];
 338   }
 339   // reverse the distortion, i.e. get the correction
 340   dist[0]=x0[0]-(dist[0]-x0[0]);
 341   dist[1]=x0[1]-(dist[1]-x0[1]);
 342 }
 343
 344 void AliTPCExBExact::DGLStep(Double_t *x,Double_t t,Double_t h) const {
 345   //
 346   // An elementary integration step.
 347   // (simple Euler Method)
 348   //
 349   Double_t dxdt[3];
 350   Motion(x,t,dxdt);
 351   for (int i=0;i<3;++i)
 352     x[i]+=h*dxdt[i];
 353
 354   /* suggestions about how to write it this way are welcome!
 355      void DGLStep(void (*f)(const Double_t *x,Double_t t,Double_t *dxdt),
 356                    Double_t *x,Double_t t,Double_t h,Int_t n) const;
 357   */
 358
 359 }