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