AliTPCcalibTimeGain.cxx - Adding the Gamma conversion selected electorns
[u/mrichter/AliRoot.git] / TPC / AliTPCExB.cxx
1 #include "AliTPCExB.h"
2 #include "TMath.h"
3 //#include "TTreeStream.h"
4 #include "AliMagF.h"
5 #include "TLinearFitter.h"
6 #include "AliTPCcalibDB.h"
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 
28   Draw(0)
29
30
31
32
33 */
34
35 AliTPCExB* AliTPCExB::fgInstance = 0;
36
37 TObjArray   AliTPCExB::fgArray;
38
39
40
41 ClassImp(AliTPCExB)
42
43 AliTPCExB::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
57 AliTPCExB::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
71 AliTPCExB& AliTPCExB::operator=(const AliTPCExB &/*exb*/)
72 {
73   //
74   // Dummy  assignment
75   //
76   return *this;
77 }
78
79
80
81
82 void AliTPCExB::TestExB(const char* fileName) {
83   //
84   // Test ExB  - 
85   // Dump the filed and corrections to the tree in file fileName  
86   //
87   // 
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.)
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];
97         Double_t d[3];
98         Correct(x,d);
99         Double_t rd=TMath::Sqrt(d[0]*d[0]+d[1]*d[1]);
100         Double_t dr=r-rd;
101         Double_t phi=TMath::ATan2(x[1],x[0]);
102         Double_t phid=TMath::ATan2(d[1],d[0]);
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];
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";
149       }
150 }
151
152
153
154 Double_t AliTPCExB::GetDr(Double_t r, Double_t phi, Double_t z, Double_t bz){
155   //
156   // Static function
157   // Posibble to us it for visualization 
158   // 
159   //
160   AliTPCExB *exb = Instance();
161   if (!exb) exb = AliTPCcalibDB::GetExB(bz,kFALSE);
162   if (!exb) return 0;
163   Double_t pos0[3] = {r*TMath::Cos(phi), r*TMath::Sin(phi),z};
164   Double_t pos1[3];
165   exb->Correct(pos0,pos1);
166   Double_t dx=pos1[0]-pos0[0];
167   Double_t dy=pos1[1]-pos0[1];
168   //  Double_t dz=pos1[2]-pos0[2];
169   // return TMath::Sqrt(dx*dx+dy*dy);  
170   Float_t dr = (dx*pos0[0]+dy*pos0[1])/r;
171   return dr;
172 }
173
174
175 Double_t AliTPCExB::GetDrphi(Double_t r, Double_t phi, Double_t z, Double_t bz){
176   //
177   //
178   //
179   AliTPCExB *exb = Instance();
180   if (!exb) exb = AliTPCcalibDB::GetExB(bz,kFALSE);
181   if (!exb) return 0;
182   Double_t pos0[3] = {r*TMath::Cos(phi), r*TMath::Sin(phi),z};
183   Double_t pos1[3];
184   exb->Correct(pos0,pos1);
185   Double_t dphi=TMath::ATan2(pos1[1],pos1[0])-TMath::ATan2(pos0[1],pos0[0]);
186   if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
187   if (dphi<-TMath::Pi()) dphi+=TMath::TwoPi();
188   return r*dphi;  
189
190 }
191
192
193 Double_t AliTPCExB::GetDphi(Double_t r, Double_t phi, Double_t z, Double_t bz){
194   //
195   //
196   // 
197   AliTPCExB *exb = Instance();
198   if (!exb) exb = AliTPCcalibDB::GetExB(bz,kFALSE);
199   if (!exb) return 0;
200   Double_t pos0[3] = {r*TMath::Cos(phi), r*TMath::Sin(phi),z};
201   Double_t pos1[3];
202   exb->Correct(pos0,pos1);
203   Double_t dphi=TMath::ATan2(pos1[1],pos1[0])-TMath::ATan2(pos0[1],pos0[0]);
204   return dphi;  
205
206 }
207
208 Double_t AliTPCExB::GetDz(Double_t r, Double_t phi, Double_t z, Double_t bz){
209   //
210   //
211   //
212   AliTPCExB *exb = Instance();
213   if (!exb) exb = AliTPCcalibDB::GetExB(bz,kFALSE);
214   if (!exb) return 0;
215   Double_t pos0[3] = {r*TMath::Cos(phi), r*TMath::Sin(phi),z};
216   Double_t pos1[3];
217   exb->Correct(pos0,pos1);
218   Double_t dz=pos1[2]-pos0[2];
219   return dz;  
220 }
221
222 //
223 // Magnetic field
224 //
225
226
227
228
229 void 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
238 Double_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;
244   Double_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
245   //  xyz[1]+=30;
246   Double_t bxyz[3];
247   mag->Field(xyz,bxyz);
248   return bxyz[0];
249 }  
250
251 Double_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;
257   Double_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
258   //  xyz[1]+=30;
259   Double_t bxyz[3];
260   mag->Field(xyz,bxyz);
261   return bxyz[1];
262 }  
263
264 Double_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;
270   Double_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
271   //  xyz[1]+=30;
272   Double_t bxyz[3];
273   mag->Field(xyz,bxyz);
274   return bxyz[2];
275 }  
276
277
278
279 Double_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;
285   Double_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
286   //xyz[1]+=30;
287   Double_t bxyz[3];
288   mag->Field(xyz,bxyz);
289   if (r==0) return 0;
290   Double_t br = (bxyz[0]*xyz[0]+bxyz[1]*xyz[1])/r;
291   return br;
292 }  
293
294 Double_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;
300   Double_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
301   //xyz[1]+=30;
302   Double_t bxyz[3];
303   mag->Field(xyz,bxyz);
304   if (r==0) return 0;
305   Double_t br = (-bxyz[0]*xyz[1]+bxyz[1]*xyz[0])/r;
306   return br;
307 }  
308
309
310
311
312 Double_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
328 Double_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
344 Double_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
361 Double_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
377 Double_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
394 Double_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
416 Double_t AliTPCExB::EvalMat(const TVectorD &vec, Double_t r, Double_t phi, Double_t z){
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