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