1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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.
22 // The ExB correction map is stored in the calib DB
23 // To test current version:
25 char *storage = "local://OCDBres"
27 AliCDBManager::Instance()->SetDefaultStorage(storage);
28 AliCDBManager::Instance()->SetRun(RunNumber)
29 AliTPCExBFirst * exb = AliTPCcalibDB::Instance()->GetExB();
31 // exb->TestExB("exb.root");
32 // TFile f("exb.root");
33 //positions->Draw("drphi");
39 #include "AliFieldMap.h"
41 #include "TTreeStream.h"
42 #include "AliTPCExBFirst.h"
44 ClassImp(AliTPCExBFirst)
46 const Double_t AliTPCExBFirst::fgkEM=1.602176487e-19/9.10938215e-31;
47 const Double_t AliTPCExBFirst::fgkDriftField=40.e3;
49 AliTPCExBFirst::AliTPCExBFirst()
51 fkNX(0),fkNY(0),fkNZ(0),
52 fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
53 fkZMin(-250.),fkZMax(250.),
55 fkMeanBx(0),fkMeanBy(0),fkMeanBz(0.) {
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.),
70 fkMeanBx(0),fkMeanBy(0),fkMeanBz(0.) {
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.
76 ConstructCommon(0,bField);
80 AliTPCExBFirst::AliTPCExBFirst(const AliFieldMap *bFieldMap,
81 Double_t driftVelocity)
82 : fDriftVelocity(driftVelocity),
83 fkNX(0),fkNY(0),fkNZ(0),
84 fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
85 fkZMin(-250.),fkZMax(250.),
87 fkMeanBx(0),fkMeanBy(0),fkMeanBz(0.) {
89 // The constructor. One has to supply a field map and an (initial)
93 fkXMin=bFieldMap->Xmin()
94 -TMath::Ceil( (bFieldMap->Xmin()+250.0)/bFieldMap->DelX())
96 fkXMax=bFieldMap->Xmax()
97 -TMath::Floor((bFieldMap->Xmax()-250.0)/bFieldMap->DelX())
99 fkYMin=bFieldMap->Ymin()
100 -TMath::Ceil( (bFieldMap->Ymin()+250.0)/bFieldMap->DelY())
102 fkYMax=bFieldMap->Ymax()
103 -TMath::Floor((bFieldMap->Ymax()-250.0)/bFieldMap->DelY())
105 fkZMin=bFieldMap->Zmin()
106 -TMath::Ceil( (bFieldMap->Zmin()+250.0)/bFieldMap->DelZ())
108 fkZMax=bFieldMap->Zmax()
109 -TMath::Floor((bFieldMap->Zmax()-250.0)/bFieldMap->DelZ())
112 fkNX=static_cast<Int_t>((fkXMax-fkXMin)/bFieldMap->DelX()+1.1);
113 fkNY=static_cast<Int_t>((fkYMax-fkYMin)/bFieldMap->DelY()+1.1);
114 fkNZ=static_cast<Int_t>((fkZMax-fkZMin)/bFieldMap->DelZ()+1.1);
116 ConstructCommon(bFieldMap,0);
119 AliTPCExBFirst::~AliTPCExBFirst() {
121 // destruct the poor object.
127 void AliTPCExBFirst::Correct(const Double_t *position,Double_t *corrected) {
129 // correct for the distortion
132 GetMeanFields(position[0],position[1],position[2],&bx,&by);
133 if (position[2]>0.) {
135 GetMeanFields(position[0],position[1],250.,&bxe,&bye);
136 if (position[2]!=250.) {
137 bx=(500.*bxe-(position[2]+250.)*bx)/(250.-position[2]);
138 by=(500.*bye-(position[2]+250.)*by)/(250.-position[2]);
146 Double_t mu=fDriftVelocity/fgkDriftField;
147 Double_t wt=mu*fkMeanBz;
149 corrected[0]=mu*(wt*bx-by)/(1.+wt*wt);
150 corrected[1]=mu*(wt*by+bx)/(1.+wt*wt);
152 if (position[2]>0.) {
153 corrected[0]*=(250.-position[2]);
154 corrected[1]*=(250.-position[2]);
157 corrected[0]*=(-250.-position[2]);
158 corrected[1]*=(-250.-position[2]);
161 corrected[0]=position[0]-corrected[0];
162 corrected[1]=position[1]-corrected[1];
163 corrected[2]=position[2];
166 void AliTPCExBFirst::TestThisBeautifulObject(const char* fileName) {
168 // well, as the name sais...
170 TTreeSRedirector ts(fileName);
172 for (x[0]=-250.;x[0]<=250.;x[0]+=10.)
173 for (x[1]=-250.;x[1]<=250.;x[1]+=10.)
174 for (x[2]=-250.;x[2]<=250.;x[2]+=10.) {
177 Double_t r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
178 Double_t rd=TMath::Sqrt(d[0]*d[0]+d[1]*d[1]);
180 Double_t phi=TMath::ATan2(x[0],x[1]);
181 Double_t phid=TMath::ATan2(d[0],d[1]);
182 Double_t dphi=phi-phid;
183 if (dphi<0.) dphi+=TMath::TwoPi();
184 if (dphi>TMath::Pi()) dphi=TMath::TwoPi()-dphi;
185 Double_t drphi=r*dphi;
186 Double_t dx=x[0]-d[0];
187 Double_t dy=x[1]-d[1];
188 Double_t dz=x[2]-d[2];
204 void AliTPCExBFirst::ConstructCommon(const AliFieldMap *bFieldMap,
205 const AliMagF *bField) {
207 // THIS IS PRIVATE! (a helper for the constructor)
209 fkNMean=fkNX*fkNY*fkNZ;
210 fkMeanBx=new Double_t[fkNMean];
211 fkMeanBy=new Double_t[fkNMean];
216 for (int i=0;i<fkNX;++i) {
217 x[0]=fkXMin+i*(fkXMax-fkXMin)/(fkNX-1);
218 for (int j=0;j<fkNY;++j) {
219 x[1]=fkYMin+j*(fkYMax-fkYMin)/(fkNY-1);
220 Double_t r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
221 Double_t bx=0.,by=0.;
222 for (int k=0;k<fkNZ;++k) {
223 x[2]=fkZMin+k*(fkZMax-fkZMin)/(fkNZ-1);
225 // the x is not const in the Field function...
227 for (int l=0;l<3;++l) xt[l]=x[l];
228 // that happens due to the lack of a sophisticated class design:
230 bFieldMap->Field(xt,b);
235 fkMeanBx[(k*fkNY+j)*fkNX+i]=bx;
236 fkMeanBy[(k*fkNY+j)*fkNX+i]=by;
237 if (90.<=r&&r<=250.) {
247 void AliTPCExBFirst::GetMeanFields(Double_t rx,Double_t ry,Double_t rz,
248 Double_t *Bx,Double_t *By) const {
250 // THIS IS PRIVATE! (calculates the mean field utilising a lookup table)
252 Double_t x=(fkNX-1)*(rx-fkXMin)/(fkXMax-fkXMin);
253 Int_t xi1=static_cast<Int_t>(x);
254 xi1=TMath::Max(TMath::Min(xi1,fkNX-2),0);
257 Double_t dx1=(xi2-x);
259 Double_t y=(fkNY-1)*(ry-fkYMin)/(fkYMax-fkYMin);
260 Int_t yi1=static_cast<Int_t>(y);
261 yi1=TMath::Max(TMath::Min(yi1,fkNY-2),0);
264 Double_t dy1=(yi2-y);
266 Double_t z=(fkNZ-1)*(rz-fkZMin)/(fkZMax-fkZMin);
267 Int_t zi1=static_cast<Int_t>(z);
268 zi1=TMath::Max(TMath::Min(zi1,fkNZ-2),0);
272 double s0x=fkMeanBx[yi1*fkNX+xi1]*dx1*dy1
273 +fkMeanBx[yi2*fkNX+xi1]*dx1*dy
274 +fkMeanBx[yi1*fkNX+xi2]*dx *dy
275 +fkMeanBx[yi2*fkNX+xi2]*dx *dy1;
276 double s0y=fkMeanBy[yi1*fkNX+xi1]*dx1*dy1
277 +fkMeanBy[yi2*fkNX+xi1]*dx1*dy
278 +fkMeanBy[yi1*fkNX+xi2]*dx *dy
279 +fkMeanBy[yi2*fkNX+xi2]*dx *dy1;
283 snmx=fkMeanBx[(zi0*fkNY+yi1)*fkNX+xi1]*dx1*dy1
284 +fkMeanBx[(zi0*fkNY+yi2)*fkNX+xi1]*dx1*dy
285 +fkMeanBx[(zi0*fkNY+yi1)*fkNX+xi2]*dx *dy
286 +fkMeanBx[(zi0*fkNY+yi2)*fkNX+xi2]*dx *dy1;
287 snmy=fkMeanBy[(zi0*fkNY+yi1)*fkNX+xi1]*dx1*dy1
288 +fkMeanBy[(zi0*fkNY+yi2)*fkNX+xi1]*dx1*dy
289 +fkMeanBy[(zi0*fkNY+yi1)*fkNX+xi2]*dx *dy
290 +fkMeanBy[(zi0*fkNY+yi2)*fkNX+xi2]*dx *dy1;
294 double snx=fkMeanBx[(zi1*fkNY+yi1)*fkNX+xi1]*dx1*dy1
295 +fkMeanBx[(zi1*fkNY+yi2)*fkNX+xi1]*dx1*dy
296 +fkMeanBx[(zi1*fkNY+yi1)*fkNX+xi2]*dx *dy
297 +fkMeanBx[(zi1*fkNY+yi2)*fkNX+xi2]*dx *dy1;
298 double sny=fkMeanBy[(zi1*fkNY+yi1)*fkNX+xi1]*dx1*dy1
299 +fkMeanBy[(zi1*fkNY+yi2)*fkNX+xi1]*dx1*dy
300 +fkMeanBy[(zi1*fkNY+yi1)*fkNX+xi2]*dx *dy
301 +fkMeanBy[(zi1*fkNY+yi2)*fkNX+xi2]*dx *dy1;
302 double snpx=fkMeanBx[(zi2*fkNY+yi1)*fkNX+xi1]*dx1*dy1
303 +fkMeanBx[(zi2*fkNY+yi2)*fkNX+xi1]*dx1*dy
304 +fkMeanBx[(zi2*fkNY+yi1)*fkNX+xi2]*dx *dy
305 +fkMeanBx[(zi2*fkNY+yi2)*fkNX+xi2]*dx *dy1;
306 double snpy=fkMeanBy[(zi2*fkNY+yi1)*fkNX+xi1]*dx1*dy1
307 +fkMeanBy[(zi2*fkNY+yi2)*fkNX+xi1]*dx1*dy
308 +fkMeanBy[(zi2*fkNY+yi1)*fkNX+xi2]*dx *dy
309 +fkMeanBy[(zi2*fkNY+yi2)*fkNX+xi2]*dx *dy1;
310 *Bx=0.5*(((snpx-2.*snx+snmx)*dz+2.*(snx-snmx))*dz+snx-s0x+snmx);
311 *By=0.5*(((snpy-2.*sny+snmy)*dz+2.*(sny-snmy))*dz+sny-s0y+snmy);
312 //TODO: make this nice
313 if (TMath::Abs(z)>0.001) {