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