Update TPCCEda to write output file in parts (to avoid too big files produced in...
[u/mrichter/AliRoot.git] / TPC / AliTPCExBFirst.cxx
CommitLineData
cf585711 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
2abfc1e6 16//
cf585711 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.
2abfc1e6 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
cf585711 37
faf93237 38#include "TMath.h"
f7a1cc68 39//#include "AliFieldMap.h"
cf585711 40#include "AliMagF.h"
faf93237 41#include "TTreeStream.h"
42#include "AliTPCExBFirst.h"
43
44ClassImp(AliTPCExBFirst)
45
46const Double_t AliTPCExBFirst::fgkEM=1.602176487e-19/9.10938215e-31;
a187c771 47const Double_t AliTPCExBFirst::fgkDriftField=-40.e3;
faf93237 48
481f877b 49AliTPCExBFirst::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 //
2abfc1e6 59 SetInstance(this);
481f877b 60}
61
faf93237 62AliTPCExBFirst::AliTPCExBFirst(const AliMagF *bField,
63 Double_t driftVelocity,
64 Int_t nx,Int_t ny,Int_t nz)
481f877b 65 : fDriftVelocity(driftVelocity),
66 fkNX(nx),fkNY(ny),fkNZ(nz),
faf93237 67 fkXMin(-250.),fkXMax(250.),fkYMin(-250.),fkYMax(250.),
481f877b 68 fkZMin(-250.),fkZMax(250.),
69 fkNMean(0),
70 fkMeanBx(0),fkMeanBy(0),fkMeanBz(0.) {
faf93237 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 //
f7a1cc68 76 // ConstructCommon(0,bField);
77 ConstructCommon(bField);
2abfc1e6 78 SetInstance(this);
faf93237 79}
80
f7a1cc68 81/*
faf93237 82AliTPCExBFirst::AliTPCExBFirst(const AliFieldMap *bFieldMap,
481f877b 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.) {
faf93237 90 //
91 // The constructor. One has to supply a field map and an (initial)
92 // drift velocity.
93 //
2abfc1e6 94 SetInstance(this);
faf93237 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}
f7a1cc68 120*/
faf93237 121
122AliTPCExBFirst::~AliTPCExBFirst() {
123 //
124 // destruct the poor object.
125 //
481f877b 126 delete[] fkMeanBx;
127 delete[] fkMeanBy;
faf93237 128}
129
130void AliTPCExBFirst::Correct(const Double_t *position,Double_t *corrected) {
131 //
132 // correct for the distortion
133 //
cf585711 134 Double_t bx,by;
135 GetMeanFields(position[0],position[1],position[2],&bx,&by);
481f877b 136 if (position[2]>0.) {
cf585711 137 Double_t bxe,bye;
138 GetMeanFields(position[0],position[1],250.,&bxe,&bye);
481f877b 139 if (position[2]!=250.) {
cf585711 140 bx=(500.*bxe-(position[2]+250.)*bx)/(250.-position[2]);
141 by=(500.*bye-(position[2]+250.)*by)/(250.-position[2]);
faf93237 142 }
143 else {
cf585711 144 bx=bxe;
145 by=bye;
faf93237 146 }
faf93237 147 }
481f877b 148
149 Double_t mu=fDriftVelocity/fgkDriftField;
150 Double_t wt=mu*fkMeanBz;
151
cf585711 152 corrected[0]=mu*(wt*bx-by)/(1.+wt*wt);
153 corrected[1]=mu*(wt*by+bx)/(1.+wt*wt);
481f877b 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];
faf93237 167}
168
169void 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
f7a1cc68 207
208void AliTPCExBFirst::ConstructCommon(//const AliFieldMap *bFieldMap,
faf93237 209 const AliMagF *bField) {
210 //
211 // THIS IS PRIVATE! (a helper for the constructor)
212 //
481f877b 213 fkNMean=fkNX*fkNY*fkNZ;
cf585711 214 fkMeanBx=new Double_t[fkNMean];
215 fkMeanBy=new Double_t[fkNMean];
faf93237 216
217 Double_t x[3];
218 Double_t nBz=0;
481f877b 219 fkMeanBz=0.;
faf93237 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);
cf585711 224 Double_t r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
225 Double_t bx=0.,by=0.;
faf93237 226 for (int k=0;k<fkNZ;++k) {
227 x[2]=fkZMin+k*(fkZMax-fkZMin)/(fkNZ-1);
f7a1cc68 228 Double_t b[3];
faf93237 229 // the x is not const in the Field function...
f7a1cc68 230 Double_t xt[3];
faf93237 231 for (int l=0;l<3;++l) xt[l]=x[l];
232 // that happens due to the lack of a sophisticated class design:
f7a1cc68 233 // if (bFieldMap!=0)
234 // bFieldMap->Field(xt,b);
235 // else
236 ((AliMagF*)bField)->Field(xt,b);
cf585711 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.;
faf93237 243 ++nBz;
244 }
245 }
246 }
247 }
481f877b 248 fkMeanBz/=nBz;
faf93237 249}
250
f7a1cc68 251
faf93237 252void 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
481f877b 277 double s0x=fkMeanBx[yi1*fkNX+xi1]*dx1*dy1
278 +fkMeanBx[yi2*fkNX+xi1]*dx1*dy
82e0477c 279 +fkMeanBx[yi1*fkNX+xi2]*dx *dy1
280 +fkMeanBx[yi2*fkNX+xi2]*dx *dy;
481f877b 281 double s0y=fkMeanBy[yi1*fkNX+xi1]*dx1*dy1
282 +fkMeanBy[yi2*fkNX+xi1]*dx1*dy
82e0477c 283 +fkMeanBy[yi1*fkNX+xi2]*dx *dy1
284 +fkMeanBy[yi2*fkNX+xi2]*dx *dy;
faf93237 285 Int_t zi0=zi1-1;
286 double snmx,snmy;
287 if (zi0>=0) {
481f877b 288 snmx=fkMeanBx[(zi0*fkNY+yi1)*fkNX+xi1]*dx1*dy1
289 +fkMeanBx[(zi0*fkNY+yi2)*fkNX+xi1]*dx1*dy
82e0477c 290 +fkMeanBx[(zi0*fkNY+yi1)*fkNX+xi2]*dx *dy1
291 +fkMeanBx[(zi0*fkNY+yi2)*fkNX+xi2]*dx *dy;
481f877b 292 snmy=fkMeanBy[(zi0*fkNY+yi1)*fkNX+xi1]*dx1*dy1
293 +fkMeanBy[(zi0*fkNY+yi2)*fkNX+xi1]*dx1*dy
82e0477c 294 +fkMeanBy[(zi0*fkNY+yi1)*fkNX+xi2]*dx *dy1
295 +fkMeanBy[(zi0*fkNY+yi2)*fkNX+xi2]*dx *dy;
faf93237 296 }
297 else
298 snmx=snmy=0.;
481f877b 299 double snx=fkMeanBx[(zi1*fkNY+yi1)*fkNX+xi1]*dx1*dy1
300 +fkMeanBx[(zi1*fkNY+yi2)*fkNX+xi1]*dx1*dy
82e0477c 301 +fkMeanBx[(zi1*fkNY+yi1)*fkNX+xi2]*dx *dy1
302 +fkMeanBx[(zi1*fkNY+yi2)*fkNX+xi2]*dx *dy;
481f877b 303 double sny=fkMeanBy[(zi1*fkNY+yi1)*fkNX+xi1]*dx1*dy1
304 +fkMeanBy[(zi1*fkNY+yi2)*fkNX+xi1]*dx1*dy
82e0477c 305 +fkMeanBy[(zi1*fkNY+yi1)*fkNX+xi2]*dx *dy1
306 +fkMeanBy[(zi1*fkNY+yi2)*fkNX+xi2]*dx *dy;
481f877b 307 double snpx=fkMeanBx[(zi2*fkNY+yi1)*fkNX+xi1]*dx1*dy1
308 +fkMeanBx[(zi2*fkNY+yi2)*fkNX+xi1]*dx1*dy
82e0477c 309 +fkMeanBx[(zi2*fkNY+yi1)*fkNX+xi2]*dx *dy1
310 +fkMeanBx[(zi2*fkNY+yi2)*fkNX+xi2]*dx *dy;
481f877b 311 double snpy=fkMeanBy[(zi2*fkNY+yi1)*fkNX+xi1]*dx1*dy1
312 +fkMeanBy[(zi2*fkNY+yi2)*fkNX+xi1]*dx1*dy
82e0477c 313 +fkMeanBy[(zi2*fkNY+yi1)*fkNX+xi2]*dx *dy1
314 +fkMeanBy[(zi2*fkNY+yi2)*fkNX+xi2]*dx *dy;
315
316
317
faf93237 318 *Bx=0.5*(((snpx-2.*snx+snmx)*dz+2.*(snx-snmx))*dz+snx-s0x+snmx);
319 *By=0.5*(((snpy-2.*sny+snmy)*dz+2.*(sny-snmy))*dz+sny-s0y+snmy);
320 //TODO: make this nice
321 if (TMath::Abs(z)>0.001) {
322 *Bx/=z;
323 *By/=z;
324 }
325}