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