TPC/AliTPCExBFirst.cxx

   1 #include "TMath.h"
   2 #include "TTreeStream.h"
   3 #include "AliTPCExBFirst.h"
   4
   5 ClassImp(AliTPCExBFirst)
   6
   7 const Double_t AliTPCExBFirst::fgkEM=1.602176487e-19/9.10938215e-31;
   8 const Double_t AliTPCExBFirst::fgkDriftField=40.e3;
   9
  10 AliTPCExBFirst::AliTPCExBFirst()
  11   : fDriftVelocity(0),
  12     fkNX(0),fkNY(0),fkNZ(0),
  13     fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
  14     fkZMin(-250.),fkZMax(250.),
  15     fkNMean(0),
  16     fkMeanBx(0),fkMeanBy(0),fkMeanBz(0.) {
  17   //
  18   // purely for I/O
  19   //
  20 }
  21
  22 AliTPCExBFirst::AliTPCExBFirst(const AliMagF *bField,
  23                                Double_t driftVelocity,
  24                                Int_t nx,Int_t ny,Int_t nz)
  25   : fDriftVelocity(driftVelocity),
  26     fkNX(nx),fkNY(ny),fkNZ(nz),
  27     fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
  28     fkZMin(-250.),fkZMax(250.),
  29     fkNMean(0),
  30     fkMeanBx(0),fkMeanBy(0),fkMeanBz(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   //
  36   ConstructCommon(0,bField);
  37 }
  38
  39 AliTPCExBFirst::AliTPCExBFirst(const AliFieldMap *bFieldMap,
  40                                Double_t driftVelocity)
  41   : fDriftVelocity(driftVelocity),
  42     fkNX(0),fkNY(0),fkNZ(0),
  43     fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
  44     fkZMin(-250.),fkZMax(250.),
  45     fkNMean(0),
  46     fkMeanBx(0),fkMeanBy(0),fkMeanBz(0.) {
  47   //
  48   // The constructor. One has to supply a field map and an (initial)
  49   // drift velocity.
  50   //
  51
  52   fkXMin=bFieldMap->Xmin()
  53     -TMath::Ceil( (bFieldMap->Xmin()+250.0)/bFieldMap->DelX())
  54     *bFieldMap->DelX();
  55   fkXMax=bFieldMap->Xmax()
  56     -TMath::Floor((bFieldMap->Xmax()-250.0)/bFieldMap->DelX())
  57     *bFieldMap->DelX();
  58   fkYMin=bFieldMap->Ymin()
  59     -TMath::Ceil( (bFieldMap->Ymin()+250.0)/bFieldMap->DelY())
  60     *bFieldMap->DelY();
  61   fkYMax=bFieldMap->Ymax()
  62     -TMath::Floor((bFieldMap->Ymax()-250.0)/bFieldMap->DelY())
  63     *bFieldMap->DelY();
  64   fkZMin=bFieldMap->Zmin()
  65     -TMath::Ceil( (bFieldMap->Zmin()+250.0)/bFieldMap->DelZ())
  66     *bFieldMap->DelZ();
  67   fkZMax=bFieldMap->Zmax()
  68     -TMath::Floor((bFieldMap->Zmax()-250.0)/bFieldMap->DelZ())
  69     *bFieldMap->DelZ();
  70
  71   fkNX=static_cast<Int_t>((fkXMax-fkXMin)/bFieldMap->DelX()+1.1);
  72   fkNY=static_cast<Int_t>((fkYMax-fkYMin)/bFieldMap->DelY()+1.1);
  73   fkNZ=static_cast<Int_t>((fkZMax-fkZMin)/bFieldMap->DelZ()+1.1);
  74
  75   ConstructCommon(bFieldMap,0);
  76 }
  77
  78 AliTPCExBFirst::~AliTPCExBFirst() {
  79   //
  80   // destruct the poor object.
  81   //
  82   delete[] fkMeanBx;
  83   delete[] fkMeanBy;
  84 }
  85
  86 void AliTPCExBFirst::Correct(const Double_t *position,Double_t *corrected) {
  87   //
  88   // correct for the distortion
  89   //
  90   Double_t Bx,By;
  91   GetMeanFields(position[0],position[1],position[2],&Bx,&By);
  92   if (position[2]>0.) {
  93     Double_t Bxe,Bye;
  94     GetMeanFields(position[0],position[1],250.,&Bxe,&Bye);
  95     if (position[2]!=250.) {
  96       Bx=(500.*Bxe-(position[2]+250.)*Bx)/(250.-position[2]);
  97       By=(500.*Bye-(position[2]+250.)*By)/(250.-position[2]);
  98     }
  99     else {
 100       Bx=Bxe;
 101       By=Bye;
 102     }
 103   }
 104
 105   Double_t mu=fDriftVelocity/fgkDriftField;
 106   Double_t wt=mu*fkMeanBz;
 107
 108   corrected[0]=mu*(wt*Bx-By)/(1.+wt*wt);
 109   corrected[1]=mu*(wt*By+Bx)/(1.+wt*wt);
 110
 111   if (position[2]>0.) {
 112     corrected[0]*=(250.-position[2]);
 113     corrected[1]*=(250.-position[2]);
 114   }
 115   else {
 116     corrected[0]*=(-250.-position[2]);
 117     corrected[1]*=(-250.-position[2]);
 118   }
 119
 120   corrected[0]=position[0]-corrected[0];
 121   corrected[1]=position[1]-corrected[1];
 122   corrected[2]=position[2];
 123 }
 124
 125 void AliTPCExBFirst::TestThisBeautifulObject(const char* fileName) {
 126   //
 127   // well, as the name sais...
 128   //
 129   TTreeSRedirector ts(fileName);
 130   Double_t x[3];
 131   for (x[0]=-250.;x[0]<=250.;x[0]+=10.)
 132     for (x[1]=-250.;x[1]<=250.;x[1]+=10.)
 133       for (x[2]=-250.;x[2]<=250.;x[2]+=10.) {
 134         Double_t d[3];
 135         Correct(x,d);
 136         Double_t r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
 137         Double_t rd=TMath::Sqrt(d[0]*d[0]+d[1]*d[1]);
 138         Double_t dr=r-rd;
 139         Double_t phi=TMath::ATan2(x[0],x[1]);
 140         Double_t phid=TMath::ATan2(d[0],d[1]);
 141         Double_t dphi=phi-phid;
 142         if (dphi<0.) dphi+=TMath::TwoPi();
 143         if (dphi>TMath::Pi()) dphi=TMath::TwoPi()-dphi;
 144         Double_t drphi=r*dphi;
 145         Double_t dx=x[0]-d[0];
 146         Double_t dy=x[1]-d[1];
 147         Double_t dz=x[2]-d[2];
 148         ts<<"positions"
 149           <<"x0="<<x[0]
 150           <<"x1="<<x[1]
 151           <<"x2="<<x[2]
 152           <<"dx="<<dx
 153           <<"dy="<<dy
 154           <<"dz="<<dz
 155           <<"r="<<r
 156           <<"phi="<<phi
 157           <<"dr="<<dr
 158           <<"drphi="<<drphi
 159           <<"\n";
 160       }
 161 }
 162
 163 void AliTPCExBFirst::ConstructCommon(const AliFieldMap *bFieldMap,
 164                                      const AliMagF *bField) {
 165   //
 166   // THIS IS PRIVATE! (a helper for the constructor)
 167   //
 168   fkNMean=fkNX*fkNY*fkNZ;
 169   fkMeanBx=new Double_t[fkNMean];
 170   fkMeanBy=new Double_t[fkNMean];
 171
 172   Double_t x[3];
 173   Double_t nBz=0;
 174   fkMeanBz=0.;
 175   for (int i=0;i<fkNX;++i) {
 176     x[0]=fkXMin+i*(fkXMax-fkXMin)/(fkNX-1);
 177     for (int j=0;j<fkNY;++j) {
 178       x[1]=fkYMin+j*(fkYMax-fkYMin)/(fkNY-1);
 179       Double_t R=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
 180       Double_t Bx=0.,By=0.;
 181       for (int k=0;k<fkNZ;++k) {
 182         x[2]=fkZMin+k*(fkZMax-fkZMin)/(fkNZ-1);
 183         Float_t B[3];
 184         // the x is not const in the Field function...
 185         Float_t xt[3];
 186         for (int l=0;l<3;++l) xt[l]=x[l];
 187         // that happens due to the lack of a sophisticated class design:
 188         if (bFieldMap!=0)
 189           bFieldMap->Field(xt,B);
 190         else
 191           bField->Field(xt,B);
 192         Bx+=B[0]/10.;
 193         By+=B[1]/10.;
 194         /*old
 195         fkMeanBx[(k*fkNY+j)*fkNX+i]=Bx/(k+1);
 196         fkMeanBy[(k*fkNY+j)*fkNX+i]=By/(k+1);
 197         */
 198         fkMeanBx[(k*fkNY+j)*fkNX+i]=Bx;
 199         fkMeanBy[(k*fkNY+j)*fkNX+i]=By;
 200         if (90.<=R&&R<=250.) {
 201           fkMeanBz+=B[2]/10.;
 202           ++nBz;
 203         }
 204       }
 205     }
 206   }
 207   fkMeanBz/=nBz;
 208 }
 209
 210 void AliTPCExBFirst::GetMeanFields(Double_t rx,Double_t ry,Double_t rz,
 211                                    Double_t *Bx,Double_t *By) const {
 212   //
 213   // THIS IS PRIVATE! (calculates the mean field utilising a lookup table)
 214   //
 215   Double_t x=(fkNX-1)*(rx-fkXMin)/(fkXMax-fkXMin);
 216   Int_t xi1=static_cast<Int_t>(x);
 217   xi1=TMath::Max(TMath::Min(xi1,fkNX-2),0);
 218   Int_t xi2=xi1+1;
 219   Double_t dx=(x-xi1);
 220   Double_t dx1=(xi2-x);
 221
 222   Double_t y=(fkNY-1)*(ry-fkYMin)/(fkYMax-fkYMin);
 223   Int_t yi1=static_cast<Int_t>(y);
 224   yi1=TMath::Max(TMath::Min(yi1,fkNY-2),0);
 225   Int_t yi2=yi1+1;
 226   Double_t dy=(y-yi1);
 227   Double_t dy1=(yi2-y);
 228
 229   Double_t z=(fkNZ-1)*(rz-fkZMin)/(fkZMax-fkZMin);
 230   Int_t zi1=static_cast<Int_t>(z);
 231   zi1=TMath::Max(TMath::Min(zi1,fkNZ-2),0);
 232   Int_t zi2=zi1+1;
 233   Double_t dz=(z-zi1);
 234
 235   double s0x=fkMeanBx[yi1*fkNX+xi1]*dx1*dy1
 236             +fkMeanBx[yi2*fkNX+xi1]*dx1*dy
 237             +fkMeanBx[yi1*fkNX+xi2]*dx *dy
 238             +fkMeanBx[yi2*fkNX+xi2]*dx *dy1;
 239   double s0y=fkMeanBy[yi1*fkNX+xi1]*dx1*dy1
 240             +fkMeanBy[yi2*fkNX+xi1]*dx1*dy
 241             +fkMeanBy[yi1*fkNX+xi2]*dx *dy
 242             +fkMeanBy[yi2*fkNX+xi2]*dx *dy1;
 243   Int_t zi0=zi1-1;
 244   double snmx,snmy;
 245   if (zi0>=0) {
 246     snmx=fkMeanBx[(zi0*fkNY+yi1)*fkNX+xi1]*dx1*dy1
 247         +fkMeanBx[(zi0*fkNY+yi2)*fkNX+xi1]*dx1*dy
 248         +fkMeanBx[(zi0*fkNY+yi1)*fkNX+xi2]*dx *dy
 249         +fkMeanBx[(zi0*fkNY+yi2)*fkNX+xi2]*dx *dy1;
 250     snmy=fkMeanBy[(zi0*fkNY+yi1)*fkNX+xi1]*dx1*dy1
 251         +fkMeanBy[(zi0*fkNY+yi2)*fkNX+xi1]*dx1*dy
 252         +fkMeanBy[(zi0*fkNY+yi1)*fkNX+xi2]*dx *dy
 253         +fkMeanBy[(zi0*fkNY+yi2)*fkNX+xi2]*dx *dy1;
 254   }
 255   else
 256     snmx=snmy=0.;
 257   double snx=fkMeanBx[(zi1*fkNY+yi1)*fkNX+xi1]*dx1*dy1
 258             +fkMeanBx[(zi1*fkNY+yi2)*fkNX+xi1]*dx1*dy
 259             +fkMeanBx[(zi1*fkNY+yi1)*fkNX+xi2]*dx *dy
 260             +fkMeanBx[(zi1*fkNY+yi2)*fkNX+xi2]*dx *dy1;
 261   double sny=fkMeanBy[(zi1*fkNY+yi1)*fkNX+xi1]*dx1*dy1
 262             +fkMeanBy[(zi1*fkNY+yi2)*fkNX+xi1]*dx1*dy
 263             +fkMeanBy[(zi1*fkNY+yi1)*fkNX+xi2]*dx *dy
 264             +fkMeanBy[(zi1*fkNY+yi2)*fkNX+xi2]*dx *dy1;
 265   double snpx=fkMeanBx[(zi2*fkNY+yi1)*fkNX+xi1]*dx1*dy1
 266              +fkMeanBx[(zi2*fkNY+yi2)*fkNX+xi1]*dx1*dy
 267              +fkMeanBx[(zi2*fkNY+yi1)*fkNX+xi2]*dx *dy
 268              +fkMeanBx[(zi2*fkNY+yi2)*fkNX+xi2]*dx *dy1;
 269   double snpy=fkMeanBy[(zi2*fkNY+yi1)*fkNX+xi1]*dx1*dy1
 270              +fkMeanBy[(zi2*fkNY+yi2)*fkNX+xi1]*dx1*dy
 271              +fkMeanBy[(zi2*fkNY+yi1)*fkNX+xi2]*dx *dy
 272              +fkMeanBy[(zi2*fkNY+yi2)*fkNX+xi2]*dx *dy1;
 273   *Bx=0.5*(((snpx-2.*snx+snmx)*dz+2.*(snx-snmx))*dz+snx-s0x+snmx);
 274   *By=0.5*(((snpy-2.*sny+snmy)*dz+2.*(sny-snmy))*dz+sny-s0y+snmy);
 275   //TODO: make this nice
 276   if (TMath::Abs(z)>0.001) {
 277     *Bx/=z;
 278     *By/=z;
 279   }
 280 }