Fixes for building of DA (Anshul)
[u/mrichter/AliRoot.git] / TPC / macros / AliMagFDraw.cxx
1 /*
2
3   .L $ALICE_ROOT/TPC/macros/AliMagFDraw.cxx+
4   AliMagFDraw draw;
5   draw.RegisterField(0,new AliMagWrapCheb("Maps","Maps", 2, 1., 10., AliMagWrapCheb::k5kG));
6   draw.RegisterField(1,new AliMagFMaps("Maps","Maps", 2, 1., 10., 2));
7
8   TF2 fbz_rz_0pi("fbz_rz_0pi","AliMagFDraw::GetBz(x,0*pi,y)",0,250,-250,250);
9   fbz_rz_0pi->Draw("surf2");
10   
11   TF1 fbz_z_90_00pi("fbz_z_90_00pi","AliMagFDraw::GetBz(90,0*pi,x)",-250,250);
12   TF1 fbz_z_90_05pi("fbz_z_90_05pi","AliMagFDraw::GetBz(90,0.5*pi,x)",-250,250);
13   TF1 fbz_z_90_10pi("fbz_z_90_10pi","AliMagFDraw::GetBz(90,1.0*pi,x)",-250,250);
14   TF1 fbz_z_90_15pi("fbz_z_90_15pi","AliMagFDraw::GetBz(90,1.5*pi,x)",-250,250);
15   fbz_z_90_00pi->SetLineColor(2);
16   fbz_z_90_05pi->SetLineColor(3);
17   fbz_z_90_10pi->SetLineColor(4);
18   fbz_z_90_15pi->SetLineColor(5);
19   fbz_z_90_00pi->Draw()
20   fbz_z_90_05pi->Draw("same")
21   fbz_z_90_15pi->Draw("same")
22   fbz_z_90_10pi->Draw("same")
23   
24
25   TF1 fbr_z_90_00pi("fbz_z_90_00pi","AliMagFDraw::GetBr(90,0*pi,x)",-250,250);
26   TF1 fbr_z_90_05pi("fbz_z_90_05pi","AliMagFDraw::GetBr(90,0.5*pi,x)",-250,250);
27   TF1 fbr_z_90_10pi("fbz_z_90_10pi","AliMagFDraw::GetBr(90,1.0*pi,x)",-250,250);
28   TF1 fbr_z_90_15pi("fbz_z_90_15pi","AliMagFDraw::GetBr(90,1.5*pi,x)",-250,250);
29   fbr_z_90_00pi->SetLineColor(2);
30   fbr_z_90_05pi->SetLineColor(3);
31   fbr_z_90_10pi->SetLineColor(4);
32   fbr_z_90_15pi->SetLineColor(5);
33   fbr_z_90_00pi->Draw()
34   fbr_z_90_05pi->Draw("same")
35   fbr_z_90_15pi->Draw("same")
36   fbr_z_90_10pi->Draw("same")
37
38   //
39   TF2 fbz_xy_0z("fbz_xy_0z","AliMagFDraw::GetBz(sqrt(x^2+y^2),atan2(y,x),0)",-250,250,-250,250);
40   fbz_xy_0z.SetNpy(100);
41   fbz_xy_0z.SetNpx(100);
42   fbz_xy_0z->Draw("colz");
43   //
44   TF2 fbz_xy_250z("fbz_xy_250z","AliMagFDraw::GetBz(sqrt(x^2+y^2),atan2(y,x),250)",-250,250,-250,250);
45   fbz_xy_250z.SetNpy(100);
46   fbz_xy_250z.SetNpx(100)
47   fbz_xy_250z->Draw("colz");
48   //
49    TF2 fbz_xy_m250z("fbz_xy_m250z","AliMagFDraw::GetBz(sqrt(x^2+y^2),atan2(y,x),-250)",-250,250,-250,250);
50   fbz_xy_m250z.SetNpy(100);
51   fbz_xy_m250z.SetNpx(100)
52   fbz_xy_m250z->Draw("colz");
53   //
54
55
56 */
57 #include "TObjArray.h"
58 #include "TMath.h"
59 #include "AliMagF.h"
60 #include "TLinearFitter.h"
61 #include "TString.h"
62
63 class AliMagFDraw  : public AliMagF
64 {
65 public:
66   AliMagFDraw():AliMagF(){}
67   static  void RegisterField(Int_t index, AliMagF * magf);
68   static  Double_t GetBx(Double_t r, Double_t phi, Double_t z,Int_t index=0);
69   static  Double_t GetBy(Double_t r, Double_t phi, Double_t z,Int_t index=0);
70   static  Double_t GetBz(Double_t r, Double_t phi, Double_t z,Int_t index=0);
71   static  Double_t GetBr(Double_t r, Double_t phi, Double_t z,Int_t index=0);
72   static  Double_t GetBrfi(Double_t r, Double_t phi, Double_t z,Int_t index=0);
73   //static  Double_t GetBr2(Double_t r, Double_t phi, Double_t z,Int_t index=0);
74   //static  Double_t GetBrfi2(Double_t r, Double_t phi, Double_t z,Int_t index=0);
75   static  TObjArray *Fit(const char *formula, Int_t index=0);
76 public:
77   static TObjArray   fgArray;
78   ClassDef(AliMagFDraw,2) 
79 };
80
81
82 ClassImp(AliMagFDraw)
83
84
85 TObjArray   AliMagFDraw::fgArray;
86
87 void AliMagFDraw::RegisterField(Int_t index, AliMagF * magf){
88   //
89   // add the filed to the list
90   //
91   fgArray.AddAt(magf,index);
92 }
93
94 Double_t AliMagFDraw::GetBz(Double_t r, Double_t phi, Double_t z,Int_t index){
95   //
96   // 
97   //
98   AliMagF *mag = (AliMagF*)fgArray.At(index);
99   if (!mag) return 0;
100   Float_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
101   //  xyz[1]+=30;
102   Float_t bxyz[3];
103   mag->Field(xyz,bxyz);
104   return bxyz[2];
105 }  
106
107 Double_t AliMagFDraw::GetBy(Double_t r, Double_t phi, Double_t z,Int_t index){
108   //
109   // 
110   //
111   AliMagF *mag = (AliMagF*)fgArray.At(index);
112   if (!mag) return 0;
113   Float_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
114   //  xyz[1]+=30;
115   Float_t bxyz[3];
116   mag->Field(xyz,bxyz);
117   return bxyz[1];
118 }  
119
120
121 Double_t AliMagFDraw::GetBx(Double_t r, Double_t phi, Double_t z,Int_t index){
122   //
123   // 
124   //
125   AliMagF *mag = (AliMagF*)fgArray.At(index);
126   if (!mag) return 0;
127   Float_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
128   //  xyz[1]+=30;
129   Float_t bxyz[3];
130   mag->Field(xyz,bxyz);
131   return bxyz[0];
132 }  
133
134
135
136
137 Double_t AliMagFDraw::GetBr(Double_t r, Double_t phi, Double_t z,Int_t index){
138   //
139   // 
140   //
141   AliMagF *mag = (AliMagF*)fgArray.At(index);
142   if (!mag) return 0;
143   Float_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
144   //xyz[1]+=30;
145   Float_t bxyz[3];
146   mag->Field(xyz,bxyz);
147   if (r==0) return 0;
148   Float_t br = bxyz[0]*xyz[0]/r+bxyz[1]*xyz[1]/r;
149   return br;
150 }  
151
152 Double_t AliMagFDraw::GetBrfi(Double_t r, Double_t phi, Double_t z,Int_t index){
153   //
154   // 
155   //
156   AliMagF *mag = (AliMagF*)fgArray.At(index);
157   if (!mag) return 0;
158   Float_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
159   //xyz[1]+=30;
160   Float_t bxyz[3];
161   mag->Field(xyz,bxyz);
162   if (r==0) return 0;
163   Float_t br = -bxyz[0]*xyz[1]/r+bxyz[1]*xyz[0]/r;
164   return br;
165 }  
166
167
168 TObjArray * AliMagFDraw::Fit(const char *formula, Int_t index){
169   //
170   /*
171     formula="1++x+x^2++cos(y)++cos(y)^2++z++z^2"
172     index=0
173   */
174   //
175   TObjArray *fstrings = TString(formula).Tokenize("++");
176   Int_t ndim = fstrings->GetEntries();
177   TObjArray *formulas = new TObjArray(ndim);
178   for (Int_t i=0;i<ndim;i++){
179     formulas->AddAt(new TFormula(Form("fff_%d",i),fstrings->At(i)->GetName()),i);                   
180   }  
181   TLinearFitter * fitR   = new TLinearFitter(ndim+1,Form("hyp%d",ndim));
182   TLinearFitter * fitRFI = new TLinearFitter(ndim+1,Form("hyp%d",ndim));
183   TLinearFitter * fitZ   = new TLinearFitter(ndim+1,Form("hyp%d",ndim));
184   Double_t x[ndim];    
185   for (Float_t r=20; r<250;r+=20){
186     for (Float_t fi=0; fi<TMath::Pi()*2;fi+=0.2){
187       for (Float_t z=-250; z<250;z+=20){
188         for (Int_t ifor=0;ifor<ndim;ifor++){
189           x[ifor]= ((TFormula*)formulas->At(ifor))->Eval(r/250.,fi,z/250.);
190         }
191         fitR->AddPoint(x,AliMagFDraw::GetBr(r,fi,z,index));
192         fitRFI->AddPoint(x,AliMagFDraw::GetBrfi(r,fi,z,index));
193         fitZ->AddPoint(x,AliMagFDraw::GetBz(r,fi,z,index));     
194       }
195     }
196   }
197   fitR->Eval();  
198   fitRFI->Eval();  
199   fitZ->Eval();  
200   TObjArray *res = new TObjArray; 
201   res->AddAt(fitR,0);
202   res->AddAt(fitRFI,1);
203   res->AddAt(fitZ,2);
204   printf("\tchi2\tn\tRMS\n");
205   printf("\t%f\t%d\t%f\n",fitR->GetChisquare(),fitR->GetNpoints(),TMath::Sqrt(fitR->GetChisquare()/fitR->GetNpoints()));
206   printf("\t%f\t%d\t%f\n",fitRFI->GetChisquare(),fitRFI->GetNpoints(),TMath::Sqrt(fitRFI->GetChisquare()/fitRFI->GetNpoints()));
207   printf("\t%f\t%d\t%f\n",fitZ->GetChisquare(),fitZ->GetNpoints(),TMath::Sqrt(fitZ->GetChisquare()/fitZ->GetNpoints()));
208
209   TFormula * funBZ = new TFormula("funBZ",formula);
210   TFormula * funBR = new TFormula("funBR",formula);
211   TFormula * funBRFI = new TFormula("funBRFI",formula);
212   TVectorD vec;
213   fitR->GetParameters(vec);
214   funBR->SetParameters(vec.GetMatrixArray());
215   fitRFI->GetParameters(vec);
216   funBRFI->SetParameters(vec.GetMatrixArray());
217   fitZ->GetParameters(vec);
218   funBZ->SetParameters(vec.GetMatrixArray());
219
220   return res;
221 }
222
223
224 TString MakeString(){
225   TString str="";
226   {
227     Int_t counter=0;
228     for (Int_t ix=0;ix<3;ix++)
229       for (Int_t iz=0;iz<3;iz++){
230         if (ix+iz>0) {
231           str+=Form("x^%d*z^%d++",ix,iz);
232           printf("x^%d*z^%d++\n",ix,iz);
233         }
234         for (Int_t iy=1;iy<4;iy++){
235           if (ix+iz+iy==0) continue;
236           str+=Form("x^%d*z^%d*sin(y*%d)++",ix,iz,iy);
237           str+=Form("x^%d*z^%d*cos(y*%d)++",ix,iz,iy);
238           printf(   "x^%d*z^%d*sin(y*%d)++\n",ix,iz,iy);
239           printf(   "x^%d*z^%d*cos(y*%d)++\n",ix,iz,iy);
240           counter++;
241         }
242       }
243     str[str.Length()-2]=0;
244   }
245   return str;
246 }
247
248
249
250 /*
251
252 TObjArray * array = AliMagFDraw::Fit("x",0)
253 //
254 chi2    n       RMS
255 15.534116       9600    0.040226
256 4.422160        9600    0.021463
257 8.378622        9600    0.029543
258
259 TObjArray * array = AliMagFDraw::Fit("x++z++z^2++x*sin(y)++x*cos(y)++z*x*sin(y)++z*x*cos(y)++z^2*x*sin(y)++z^2*x*cos(y)++z^2*x*sin(y)^2++z^2*x*cos(y)^2++z^3*x*sin(y)^2++z^3*x*cos(y)^2",0);
260 chi2    n       RMS
261 1.842443        9600    0.013854
262 0.957056        9600    0.009985
263 0.097799        9600    0.003192
264
265 TObjArray * array = AliMagFDraw::Fit("x++z++z^2++x*sin(y)++x*cos(y)++z*x*sin(y)++z*x*cos(y)++z^2*x*sin(y)++z^2*x*cos(y)++z^2*x*sin(y)^2++z^2*x*cos(y)^2++z^3*x*sin(y)^2++z^3*x*cos(y)++z^3*sin(y)++z^3*cos(y)++z^3*x*cos(y)++z^2*sin(y)++z^2*cos(y)++z*sin(y)++z*cos(y)++sin(y)++cos(y)",0);
266
267 chi2    n       RMS
268 1.564687        9600    0.012767
269 0.002291        9600    0.000489
270 0.097063        9600    0.003180
271
272
273 TObjArray * array = AliMagFDraw::Fit(MakeString(),0);
274
275 chi2    n       RMS
276 0.000303        9600    0.000178
277 0.000066        9600    0.000083
278 0.003110        9600    0.000569
279
280
281 }
282 */
283