Adding simple example to load default debug streamer
[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"
cf585711 39#include "AliFieldMap.h"
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;
47const Double_t AliTPCExBFirst::fgkDriftField=40.e3;
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 //
faf93237 76 ConstructCommon(0,bField);
2abfc1e6 77 SetInstance(this);
faf93237 78}
79
80AliTPCExBFirst::AliTPCExBFirst(const AliFieldMap *bFieldMap,
481f877b 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.),
86 fkNMean(0),
87 fkMeanBx(0),fkMeanBy(0),fkMeanBz(0.) {
faf93237 88 //
89 // The constructor. One has to supply a field map and an (initial)
90 // drift velocity.
91 //
2abfc1e6 92 SetInstance(this);
faf93237 93 fkXMin=bFieldMap->Xmin()
94 -TMath::Ceil( (bFieldMap->Xmin()+250.0)/bFieldMap->DelX())
95 *bFieldMap->DelX();
96 fkXMax=bFieldMap->Xmax()
97 -TMath::Floor((bFieldMap->Xmax()-250.0)/bFieldMap->DelX())
98 *bFieldMap->DelX();
99 fkYMin=bFieldMap->Ymin()
100 -TMath::Ceil( (bFieldMap->Ymin()+250.0)/bFieldMap->DelY())
101 *bFieldMap->DelY();
102 fkYMax=bFieldMap->Ymax()
103 -TMath::Floor((bFieldMap->Ymax()-250.0)/bFieldMap->DelY())
104 *bFieldMap->DelY();
105 fkZMin=bFieldMap->Zmin()
106 -TMath::Ceil( (bFieldMap->Zmin()+250.0)/bFieldMap->DelZ())
107 *bFieldMap->DelZ();
108 fkZMax=bFieldMap->Zmax()
109 -TMath::Floor((bFieldMap->Zmax()-250.0)/bFieldMap->DelZ())
110 *bFieldMap->DelZ();
111
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);
115
116 ConstructCommon(bFieldMap,0);
117}
118
119AliTPCExBFirst::~AliTPCExBFirst() {
120 //
121 // destruct the poor object.
122 //
481f877b 123 delete[] fkMeanBx;
124 delete[] fkMeanBy;
faf93237 125}
126
127void AliTPCExBFirst::Correct(const Double_t *position,Double_t *corrected) {
128 //
129 // correct for the distortion
130 //
cf585711 131 Double_t bx,by;
132 GetMeanFields(position[0],position[1],position[2],&bx,&by);
481f877b 133 if (position[2]>0.) {
cf585711 134 Double_t bxe,bye;
135 GetMeanFields(position[0],position[1],250.,&bxe,&bye);
481f877b 136 if (position[2]!=250.) {
cf585711 137 bx=(500.*bxe-(position[2]+250.)*bx)/(250.-position[2]);
138 by=(500.*bye-(position[2]+250.)*by)/(250.-position[2]);
faf93237 139 }
140 else {
cf585711 141 bx=bxe;
142 by=bye;
faf93237 143 }
faf93237 144 }
481f877b 145
146 Double_t mu=fDriftVelocity/fgkDriftField;
147 Double_t wt=mu*fkMeanBz;
148
cf585711 149 corrected[0]=mu*(wt*bx-by)/(1.+wt*wt);
150 corrected[1]=mu*(wt*by+bx)/(1.+wt*wt);
481f877b 151
152 if (position[2]>0.) {
153 corrected[0]*=(250.-position[2]);
154 corrected[1]*=(250.-position[2]);
155 }
156 else {
157 corrected[0]*=(-250.-position[2]);
158 corrected[1]*=(-250.-position[2]);
159 }
160
161 corrected[0]=position[0]-corrected[0];
162 corrected[1]=position[1]-corrected[1];
163 corrected[2]=position[2];
faf93237 164}
165
166void AliTPCExBFirst::TestThisBeautifulObject(const char* fileName) {
167 //
168 // well, as the name sais...
169 //
170 TTreeSRedirector ts(fileName);
171 Double_t x[3];
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.) {
175 Double_t d[3];
176 Correct(x,d);
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]);
179 Double_t dr=r-rd;
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];
189 ts<<"positions"
190 <<"x0="<<x[0]
191 <<"x1="<<x[1]
192 <<"x2="<<x[2]
193 <<"dx="<<dx
194 <<"dy="<<dy
195 <<"dz="<<dz
196 <<"r="<<r
197 <<"phi="<<phi
198 <<"dr="<<dr
199 <<"drphi="<<drphi
200 <<"\n";
201 }
202}
203
204void AliTPCExBFirst::ConstructCommon(const AliFieldMap *bFieldMap,
205 const AliMagF *bField) {
206 //
207 // THIS IS PRIVATE! (a helper for the constructor)
208 //
481f877b 209 fkNMean=fkNX*fkNY*fkNZ;
cf585711 210 fkMeanBx=new Double_t[fkNMean];
211 fkMeanBy=new Double_t[fkNMean];
faf93237 212
213 Double_t x[3];
214 Double_t nBz=0;
481f877b 215 fkMeanBz=0.;
faf93237 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);
cf585711 220 Double_t r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
221 Double_t bx=0.,by=0.;
faf93237 222 for (int k=0;k<fkNZ;++k) {
223 x[2]=fkZMin+k*(fkZMax-fkZMin)/(fkNZ-1);
cf585711 224 Float_t b[3];
faf93237 225 // the x is not const in the Field function...
226 Float_t xt[3];
227 for (int l=0;l<3;++l) xt[l]=x[l];
228 // that happens due to the lack of a sophisticated class design:
229 if (bFieldMap!=0)
cf585711 230 bFieldMap->Field(xt,b);
faf93237 231 else
cf585711 232 bField->Field(xt,b);
233 bx+=b[0]/10.;
234 by+=b[1]/10.;
235 fkMeanBx[(k*fkNY+j)*fkNX+i]=bx;
236 fkMeanBy[(k*fkNY+j)*fkNX+i]=by;
237 if (90.<=r&&r<=250.) {
238 fkMeanBz+=b[2]/10.;
faf93237 239 ++nBz;
240 }
241 }
242 }
243 }
481f877b 244 fkMeanBz/=nBz;
faf93237 245}
246
247void AliTPCExBFirst::GetMeanFields(Double_t rx,Double_t ry,Double_t rz,
248 Double_t *Bx,Double_t *By) const {
249 //
250 // THIS IS PRIVATE! (calculates the mean field utilising a lookup table)
251 //
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);
255 Int_t xi2=xi1+1;
256 Double_t dx=(x-xi1);
257 Double_t dx1=(xi2-x);
258
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);
262 Int_t yi2=yi1+1;
263 Double_t dy=(y-yi1);
264 Double_t dy1=(yi2-y);
265
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);
269 Int_t zi2=zi1+1;
270 Double_t dz=(z-zi1);
271
481f877b 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;
faf93237 280 Int_t zi0=zi1-1;
281 double snmx,snmy;
282 if (zi0>=0) {
481f877b 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;
faf93237 291 }
292 else
293 snmx=snmy=0.;
481f877b 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;
faf93237 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) {
314 *Bx/=z;
315 *By/=z;
316 }
317}