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