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