]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCExB.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
CORR does not need a runloader
[u/mrichter/AliRoot.git] / TPC / AliTPCExB.cxx
CommitLineData
faf93237 1#include "AliTPCExB.h"
2abfc1e6 2#include "TMath.h"
3#include "TTreeStream.h"
9170bb1b 4#include "AliMagF.h"
5#include "TLinearFitter.h"
2abfc1e6 6
7
8//
9// Abstract class for ExB effect parameterization
10//
11//
12//
13// The ExB correction map is stored in the calib DB
14// The lookup can be dumped to the tree:
15/*
16
17 //
18 char *storage = "local://OCDBres"
19 Int_t RunNumber=0;
20 AliCDBManager::Instance()->SetDefaultStorage(storage);
21 AliCDBManager::Instance()->SetRun(RunNumber)
22 AliTPCExBFirst * exb = AliTPCcalibDB::Instance()->GetExB();
23 //
24 // See example macro $ALICE_ROOT/TPC/macros/AliTPCExBdraw.C
25 //
26 .L $ALICE_ROOT/TPC/macros/AliTPCExBdraw.C
74e16bc1 27 Draw(0)
2abfc1e6 28
29
30
31
32*/
33
34AliTPCExB* AliTPCExB::fgInstance = 0;
faf93237 35
9170bb1b 36TObjArray AliTPCExB::fgArray;
37
38
39
faf93237 40ClassImp(AliTPCExB)
41
9170bb1b 42AliTPCExB::AliTPCExB():
43 TObject(),
44 fMatBrBz(0), //param matrix Br/Bz
45 fMatBrfiBz(0), //param matrix Br/Bz
46 fMatBrBzI0(0), //param matrix Br/Bz integral z>0
47 fMatBrBzI1(0), //param matrix Br/Bz integral z<0
48 fMatBrfiBzI0(0), //param matrix Br/Bz integral z>0
49 fMatBrfiBzI1(0) //param matrix Br/Bz integral z<0
50{
51 //
52 // default constructor
53 //
54}
55
2abfc1e6 56
57void AliTPCExB::TestExB(const char* fileName) {
58 //
9170bb1b 59 // Test ExB -
60 // Dump the filed and corrections to the tree in file fileName
2abfc1e6 61 //
9170bb1b 62 //
2abfc1e6 63 TTreeSRedirector ts(fileName);
64 Double_t x[3];
65 for (x[0]=-250.;x[0]<=250.;x[0]+=10.)
66 for (x[1]=-250.;x[1]<=250.;x[1]+=10.)
9170bb1b 67 for (x[2]=-250.;x[2]<=250.;x[2]+=20.) {
68 Double_t r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
69 if (r<20) continue;
70 if (r>260) continue;
71 Double_t z = x[2];
2abfc1e6 72 Double_t d[3];
73 Correct(x,d);
2abfc1e6 74 Double_t rd=TMath::Sqrt(d[0]*d[0]+d[1]*d[1]);
75 Double_t dr=r-rd;
9170bb1b 76 Double_t phi=TMath::ATan2(x[1],x[0]);
77 Double_t phid=TMath::ATan2(d[1],d[0]);
2abfc1e6 78 Double_t dphi=phi-phid;
79 if (dphi<0.) dphi+=TMath::TwoPi();
80 if (dphi>TMath::Pi()) dphi=TMath::TwoPi()-dphi;
81 Double_t drphi=r*dphi;
82 Double_t dx=x[0]-d[0];
83 Double_t dy=x[1]-d[1];
84 Double_t dz=x[2]-d[2];
9170bb1b 85 //
86 Double_t bx = GetBx(r,phi,z,0);
87 Double_t by = GetBy(r,phi,z,0);
88 Double_t bz = GetBz(r,phi,z,0);
89 Double_t br = GetBr(r,phi,z,0);
90 Double_t brfi = GetBrfi(r,phi,z,0);
91 //
92 Double_t bxi = GetBxI(r,phi,z,0);
93 Double_t byi = GetByI(r,phi,z,0);
94 Double_t bzi = GetBzI(r,phi,z,0);
95 Double_t bri = GetBrI(r,phi,z,0);
96 Double_t brfii = GetBrfiI(r,phi,z,0);
97
98 ts<<"positions"<<
99 "x0="<<x[0]<<
100 "x1="<<x[1]<<
101 "x2="<<x[2]<<
102 "dx="<<dx<<
103 "dy="<<dy<<
104 "dz="<<dz<<
105 "r="<<r<<
106 "phi="<<phi<<
107 "dr="<<dr<<
108 "drphi="<<drphi<<
109 //
110 // B-Field
111 //
112 "bx="<<bx<<
113 "by="<<by<<
114 "bz="<<bz<<
115 "br="<< br<<
116 "brfi="<<brfi<<
117 // B-field integ
118 "bxi="<<bxi<<
119 "byi="<<byi<<
120 "bzi="<<bzi<<
121 "bri="<< bri<<
122 "brfii="<<brfii<<
123 "\n";
2abfc1e6 124 }
125}
126
127
128
129Double_t AliTPCExB::GetDr(Double_t r, Double_t phi, Double_t z){
130 //
131 // Static function
132 // Posibble to us it for visualization
133 //
134 //
135 if (!fgInstance) return 0;
136 Double_t pos0[3] = {r*TMath::Sin(phi), r*TMath::Cos(phi),z};
137 Double_t pos1[3];
138 fgInstance->Correct(pos0,pos1);
139 Double_t dx=pos1[0]-pos0[0];
140 Double_t dy=pos1[1]-pos0[1];
141 // Double_t dz=pos1[2]-pos0[2];
142 return TMath::Sqrt(dx*dx+dy*dy);
143}
144
145
146Double_t AliTPCExB::GetDrphi(Double_t r, Double_t phi, Double_t z){
147 //
148 //
149 //
150 if (!fgInstance) return 0;
151 Double_t pos0[3] = {r*TMath::Sin(phi), r*TMath::Cos(phi),z};
152 Double_t pos1[3];
153 fgInstance->Correct(pos0,pos1);
154 Double_t dphi=TMath::ATan2(pos1[1],pos1[0])-TMath::ATan2(pos0[1],pos0[0]);
155 return r*dphi;
156
157}
158
159
160Double_t AliTPCExB::GetDphi(Double_t r, Double_t phi, Double_t z){
161 //
162 //
163 //
164 if (!fgInstance) return 0;
165 Double_t pos0[3] = {r*TMath::Sin(phi), r*TMath::Cos(phi),z};
166 Double_t pos1[3];
167 fgInstance->Correct(pos0,pos1);
168 Double_t dphi=TMath::ATan2(pos1[1],pos1[0])-TMath::ATan2(pos0[1],pos0[0]);
169 return dphi;
170
171}
172
173Double_t AliTPCExB::GetDz(Double_t r, Double_t phi, Double_t z){
174 //
175 //
176 //
177 if (!fgInstance) return 0;
178 Double_t pos0[3] = {r*TMath::Sin(phi), r*TMath::Cos(phi),z};
179 Double_t pos1[3];
180 fgInstance->Correct(pos0,pos1);
181 Double_t dz=pos1[2]-pos0[2];
182 return dz;
183}
9170bb1b 184
185//
186// Magnetic field
187//
188
189
190
191
192void AliTPCExB::RegisterField(Int_t index, AliMagF * magf){
193 //
194 // add the filed to the list
195 //
196 fgArray.AddAt(magf,index);
197}
198
199
200
201Double_t AliTPCExB::GetBx(Double_t r, Double_t phi, Double_t z,Int_t index){
202 //
203 //
204 //
205 AliMagF *mag = (AliMagF*)fgArray.At(index);
206 if (!mag) return 0;
207 Float_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
208 // xyz[1]+=30;
209 Float_t bxyz[3];
210 mag->Field(xyz,bxyz);
211 return bxyz[0];
212}
213
214Double_t AliTPCExB::GetBy(Double_t r, Double_t phi, Double_t z,Int_t index){
215 //
216 //
217 //
218 AliMagF *mag = (AliMagF*)fgArray.At(index);
219 if (!mag) return 0;
220 Float_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
221 // xyz[1]+=30;
222 Float_t bxyz[3];
223 mag->Field(xyz,bxyz);
224 return bxyz[1];
225}
226
227Double_t AliTPCExB::GetBz(Double_t r, Double_t phi, Double_t z,Int_t index){
228 //
229 //
230 //
231 AliMagF *mag = (AliMagF*)fgArray.At(index);
232 if (!mag) return 0;
233 Float_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
234 // xyz[1]+=30;
235 Float_t bxyz[3];
236 mag->Field(xyz,bxyz);
237 return bxyz[2];
238}
239
240
241
242Double_t AliTPCExB::GetBr(Double_t r, Double_t phi, Double_t z,Int_t index){
243 //
244 //
245 //
246 AliMagF *mag = (AliMagF*)fgArray.At(index);
247 if (!mag) return 0;
248 Float_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
249 //xyz[1]+=30;
250 Float_t bxyz[3];
251 mag->Field(xyz,bxyz);
252 if (r==0) return 0;
253 Float_t br = (bxyz[0]*xyz[0]+bxyz[1]*xyz[1])/r;
254 return br;
255}
256
257Double_t AliTPCExB::GetBrfi(Double_t r, Double_t phi, Double_t z,Int_t index){
258 //
259 //
260 //
261 AliMagF *mag = (AliMagF*)fgArray.At(index);
262 if (!mag) return 0;
263 Float_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
264 //xyz[1]+=30;
265 Float_t bxyz[3];
266 mag->Field(xyz,bxyz);
267 if (r==0) return 0;
268 Float_t br = (-bxyz[0]*xyz[1]+bxyz[1]*xyz[0])/r;
269 return br;
270}
271
272
273
274
275Double_t AliTPCExB::GetBxI(Double_t r, Double_t phi, Double_t z,Int_t index)
276{
277 Double_t sumf =0;
278 if (z>0 &&z<250){
279 for (Float_t zi=z;zi<250;zi+=5){
280 sumf+=GetBx(r,phi,zi,index)/GetBz(r,phi,zi,index);
281 }
282 }
283 if (z<0 &&z>-250){
284 for (Float_t zi=z;zi>-250;zi-=5){
285 sumf+=GetBx(r,phi,zi,index)/GetBz(r,phi,zi,index);
286 }
287 }
288 return sumf*5;
289}
290
291Double_t AliTPCExB::GetByI(Double_t r, Double_t phi, Double_t z,Int_t index)
292{
293 Double_t sumf =0;
294 if (z>0 &&z<250){
295 for (Float_t zi=z;zi<250;zi+=5){
296 sumf+=GetBy(r,phi,zi,index)/GetBz(r,phi,zi,index);
297 }
298 }
299 if (z<0 &&z>-250){
300 for (Float_t zi=z;zi>-250;zi-=5){
301 sumf+=GetBy(r,phi,zi,index)/GetBz(r,phi,zi,index);
302 }
303 }
304 return sumf*5;
305}
306
307Double_t AliTPCExB::GetBzI(Double_t r, Double_t phi, Double_t z,Int_t index)
308{
309 Double_t sumf =0;
310 if (z>0 &&z<250){
311 for (Float_t zi=z;zi<250;zi+=5){
312 sumf+=GetBz(r,phi,zi,index);
313 }
314 }
315 if (z<0 &&z>-250){
316 for (Float_t zi=z;zi>-250;zi-=5){
317 sumf+=GetBz(r,phi,zi,index);
318 }
319 }
320 return sumf*5;
321}
322
323
324Double_t AliTPCExB::GetBrI(Double_t r, Double_t phi, Double_t z,Int_t index)
325{
326 Double_t sumf =0;
327 if (z>0 &&z<250){
328 for (Float_t zi=z;zi<250;zi+=5){
329 sumf+=GetBr(r,phi,zi,index)/GetBz(r,phi,zi,index);
330 }
331 }
332 if (z<0 &&z>-250){
333 for (Float_t zi=z;zi>-250;zi-=5){
334 sumf+=GetBr(r,phi,zi,index)/GetBz(r,phi,zi,index);
335 }
336 }
337 return sumf*5;
338}
339
340Double_t AliTPCExB::GetBrfiI(Double_t r, Double_t phi, Double_t z,Int_t index)
341{
342 Double_t sumf =0;
343 if (z>0 &&z<250){
344 for (Float_t zi=z;zi<250;zi+=5.){
345 sumf+=GetBrfi(r,phi,zi,index)/GetBz(r,phi,zi,index);
346 }
347 }
348 if (z<0 &&z>-250){
349 for (Float_t zi=z;zi>-250;zi-=5){
350 sumf+=GetBrfi(r,phi,zi,index)/GetBz(r,phi,zi,index);
351 }
352 }
353 return sumf*5;
354}
355
356
357Double_t AliTPCExB::Eval(Int_t type, Double_t r, Double_t phi, Double_t z){
358 //
359 // Evaluate parameterization
360 //
361 // br integral param
362 if (type==0) {
363 if (z>0 && fMatBrBzI0) return EvalMat(*fMatBrBzI0,r,phi,z);
364 if (z<0 && fMatBrBzI1) return EvalMat(*fMatBrBzI1,r,phi,z);
365 }
366 // brfi integral param
367 if (type==1) {
368 if (z>0 && fMatBrfiBzI0) return EvalMat(*fMatBrfiBzI0,r,phi,z);
369 if (z<0 && fMatBrfiBzI1) return EvalMat(*fMatBrfiBzI1,r,phi,z);
370 }
371 // brbz param
372 if (type==2 && fMatBrBz) return EvalMat(*fMatBrBz,r,phi,z);
373 // brfibz param
374 if (type==3 && fMatBrfiBz) return EvalMat(*fMatBrfiBz,r,phi,z);
375 return 0;
376}
377
378
379Double_t AliTPCExB::EvalMat(TVectorD &vec, Double_t r, Double_t phi, Double_t z){
380 //
381 // Evaluate taylor expansion in r,phi,z
382 //
383 // Variables
384 //tree->SetAlias("sa","sin(phi+0.0)");
385 //tree->SetAlias("ca","cos(phi+0.0)");
386 //tree->SetAlias("sa2","sin(phi*2+0.0)");
387 //tree->SetAlias("ca2","cos(phi*2+0.0)");
388 //tree->SetAlias("zn","(x2/250.)");
389 //tree->SetAlias("rn","(r/250.)")
390 // TString fstringSym="";
391 // //
392 // fstringSym+="zn++";
393 // fstringSym+="rn++";
394 // fstringSym+="zn*rn++";
395 // fstringSym+="zn*zn++";
396 // fstringSym+="zn*zn*rn++";
397 // fstringSym+="zn*rn*rn++";
398 // //
399 // fstringSym+="sa++";
400 // fstringSym+="ca++";
401 // fstringSym+="ca2++";
402 // fstringSym+="sa2++";
403 // fstringSym+="ca*zn++";
404 // fstringSym+="sa*zn++";
405 // fstringSym+="ca2*zn++";
406 // fstringSym+="sa2*zn++";
407 // fstringSym+="ca*zn*zn++";
408 // fstringSym+="sa*zn*zn++";
409 // fstringSym+="ca*zn*rn++";
410 // fstringSym+="sa*zn*rn++";
411
412
413 Double_t sa = TMath::Sin(phi);
414 Double_t ca = TMath::Cos(phi);
415 Double_t sa2 = TMath::Sin(phi*2);
416 Double_t ca2 = TMath::Cos(phi*2);
417 Double_t zn = z/250.;
418 Double_t rn = r/250.;
419 Int_t ipoint=0;
420 Double_t res = vec[ipoint++];
421 res+=vec[ipoint++]*zn;
422 res+=vec[ipoint++]*rn;
423 res+=vec[ipoint++]*zn*rn;
424 res+=vec[ipoint++]*zn*zn;
425 res+=vec[ipoint++]*zn*zn*rn;
426 res+=vec[ipoint++]*zn*rn*rn;
427 //
428 res+=vec[ipoint++]*sa;
429 res+=vec[ipoint++]*ca;
430 res+=vec[ipoint++]*ca2;
431 res+=vec[ipoint++]*sa2;
432 res+=vec[ipoint++]*ca*zn;
433 res+=vec[ipoint++]*sa*zn;
434 res+=vec[ipoint++]*ca2*zn;
435 res+=vec[ipoint++]*sa2*zn;
436 res+=vec[ipoint++]*ca*zn*zn;
437 res+=vec[ipoint++]*sa*zn*zn;
438 res+=vec[ipoint++]*ca*zn*rn;
439 res+=vec[ipoint++]*sa*zn*rn;
440 return res;
441}
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456/*
457
458 AliTPCExB draw;
459 draw.RegisterField(0,new AliMagWrapCheb("Maps","Maps", 2, 1., 10., AliMagWrapCheb::k5kG));
460 draw.RegisterField(1,new AliMagFMaps("Maps","Maps", 2, 1., 10., 2));
461
462 TF2 fbz_rz_0pi("fbz_rz_0pi","AliTPCExB::GetBz(x,0*pi,y)",0,250,-250,250);
463 fbz_rz_0pi->Draw("surf2");
464
465 TF1 fbz_z_90_00pi("fbz_z_90_00pi","AliTPCExB::GetBz(90,0*pi,x)",-250,250);
466 TF1 fbz_z_90_05pi("fbz_z_90_05pi","AliTPCExB::GetBz(90,0.5*pi,x)",-250,250);
467 TF1 fbz_z_90_10pi("fbz_z_90_10pi","AliTPCExB::GetBz(90,1.0*pi,x)",-250,250);
468 TF1 fbz_z_90_15pi("fbz_z_90_15pi","AliTPCExB::GetBz(90,1.5*pi,x)",-250,250);
469 fbz_z_90_00pi->SetLineColor(2);
470 fbz_z_90_05pi->SetLineColor(3);
471 fbz_z_90_10pi->SetLineColor(4);
472 fbz_z_90_15pi->SetLineColor(5);
473 fbz_z_90_00pi->Draw()
474 fbz_z_90_05pi->Draw("same")
475 fbz_z_90_15pi->Draw("same")
476 fbz_z_90_10pi->Draw("same")
477
478
479 TF1 fbr_z_90_00pi("fbz_z_90_00pi","AliTPCExB::GetBr(90,0*pi,x)",-250,250);
480 TF1 fbr_z_90_05pi("fbz_z_90_05pi","AliTPCExB::GetBr(90,0.5*pi,x)",-250,250);
481 TF1 fbr_z_90_10pi("fbz_z_90_10pi","AliTPCExB::GetBr(90,1.0*pi,x)",-250,250);
482 TF1 fbr_z_90_15pi("fbz_z_90_15pi","AliTPCExB::GetBr(90,1.5*pi,x)",-250,250);
483 fbr_z_90_00pi->SetLineColor(2);
484 fbr_z_90_05pi->SetLineColor(3);
485 fbr_z_90_10pi->SetLineColor(4);
486 fbr_z_90_15pi->SetLineColor(5);
487 fbr_z_90_00pi->Draw()
488 fbr_z_90_05pi->Draw("same")
489 fbr_z_90_15pi->Draw("same")
490 fbr_z_90_10pi->Draw("same")
491
492 //
493 TF2 fbz_xy_0z("fbz_xy_0z","AliTPCExB::GetBz(sqrt(x^2+y^2),atan2(y,x),0)",-250,250,-250,250);
494 fbz_xy_0z.SetNpy(100);
495 fbz_xy_0z.SetNpx(100);
496 fbz_xy_0z->Draw("colz");
497 //
498 TF2 fbz_xy_250z("fbz_xy_250z","AliTPCExB::GetBz(sqrt(x^2+y^2),atan2(y,x),250)",-250,250,-250,250);
499 fbz_xy_250z.SetNpy(100);
500 fbz_xy_250z.SetNpx(100)
501 fbz_xy_250z->Draw("colz");
502 //
503 TF2 fbz_xy_m250z("fbz_xy_m250z","AliTPCExB::GetBz(sqrt(x^2+y^2),atan2(y,x),-250)",-250,250,-250,250);
504 fbz_xy_m250z.SetNpy(100);
505 fbz_xy_m250z.SetNpx(100)
506 fbz_xy_m250z->Draw("colz");
507 //
508
509
510
511
512
513*/
514
515