AliRoot Core  3dc7879 (3dc7879)
AliMagFDraw.cxx
Go to the documentation of this file.
1 
54 #include "TObjArray.h"
55 #include "TMath.h"
56 #include "AliMagF.h"
57 #include "TLinearFitter.h"
58 #include "TString.h"
59 
60 class AliMagFDraw : public AliMagF
61 {
62 public:
64  static void RegisterField(Int_t index, AliMagF * magf);
65  static Double_t GetBx(Double_t r, Double_t phi, Double_t z,Int_t index=0);
66  static Double_t GetBy(Double_t r, Double_t phi, Double_t z,Int_t index=0);
67  static Double_t GetBz(Double_t r, Double_t phi, Double_t z,Int_t index=0);
68  static Double_t GetBr(Double_t r, Double_t phi, Double_t z,Int_t index=0);
69  static Double_t GetBrfi(Double_t r, Double_t phi, Double_t z,Int_t index=0);
70  //static Double_t GetBr2(Double_t r, Double_t phi, Double_t z,Int_t index=0);
71  //static Double_t GetBrfi2(Double_t r, Double_t phi, Double_t z,Int_t index=0);
72  static TObjArray *Fit(const char *formula, Int_t index=0);
73 public:
75  ClassDef(AliMagFDraw,2)
76 };
77 
78 
80 ClassImp(AliMagFDraw)
82 
83 
85 
86 void AliMagFDraw::RegisterField(Int_t index, AliMagF * magf){
88 
89  fgArray.AddAt(magf,index);
90 }
91 
92 Double_t AliMagFDraw::GetBz(Double_t r, Double_t phi, Double_t z,Int_t index){
94 
95  AliMagF *mag = (AliMagF*)fgArray.At(index);
96  if (!mag) return 0;
97  Float_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
98  // xyz[1]+=30;
99  Float_t bxyz[3];
100  mag->Field(xyz,bxyz);
101  return bxyz[2];
102 }
103 
104 Double_t AliMagFDraw::GetBy(Double_t r, Double_t phi, Double_t z,Int_t index){
106 
107  AliMagF *mag = (AliMagF*)fgArray.At(index);
108  if (!mag) return 0;
109  Float_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
110  // xyz[1]+=30;
111  Float_t bxyz[3];
112  mag->Field(xyz,bxyz);
113  return bxyz[1];
114 }
115 
116 
117 Double_t AliMagFDraw::GetBx(Double_t r, Double_t phi, Double_t z,Int_t index){
119 
120  AliMagF *mag = (AliMagF*)fgArray.At(index);
121  if (!mag) return 0;
122  Float_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
123  // xyz[1]+=30;
124  Float_t bxyz[3];
125  mag->Field(xyz,bxyz);
126  return bxyz[0];
127 }
128 
129 
130 
131 
132 Double_t AliMagFDraw::GetBr(Double_t r, Double_t phi, Double_t z,Int_t index){
134 
135  AliMagF *mag = (AliMagF*)fgArray.At(index);
136  if (!mag) return 0;
137  Float_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
138  //xyz[1]+=30;
139  Float_t bxyz[3];
140  mag->Field(xyz,bxyz);
141  if (r==0) return 0;
142  Float_t br = bxyz[0]*xyz[0]/r+bxyz[1]*xyz[1]/r;
143  return br;
144 }
145 
146 Double_t AliMagFDraw::GetBrfi(Double_t r, Double_t phi, Double_t z,Int_t index){
148 
149  AliMagF *mag = (AliMagF*)fgArray.At(index);
150  if (!mag) return 0;
151  Float_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
152  //xyz[1]+=30;
153  Float_t bxyz[3];
154  mag->Field(xyz,bxyz);
155  if (r==0) return 0;
156  Float_t br = -bxyz[0]*xyz[1]/r+bxyz[1]*xyz[0]/r;
157  return br;
158 }
159 
160 
161 TObjArray * AliMagFDraw::Fit(const char *formula, Int_t index){
164  TObjArray *fstrings = TString(formula).Tokenize("++");
165  Int_t ndim = fstrings->GetEntries();
166  TObjArray *formulas = new TObjArray(ndim);
167  for (Int_t i=0;i<ndim;i++){
168  formulas->AddAt(new TFormula(Form("fff_%d",i),fstrings->At(i)->GetName()),i);
169  }
170  TLinearFitter * fitR = new TLinearFitter(ndim+1,Form("hyp%d",ndim));
171  TLinearFitter * fitRFI = new TLinearFitter(ndim+1,Form("hyp%d",ndim));
172  TLinearFitter * fitZ = new TLinearFitter(ndim+1,Form("hyp%d",ndim));
173  Double_t x[ndim];
174  for (Float_t r=20; r<250;r+=20){
175  for (Float_t fi=0; fi<TMath::Pi()*2;fi+=0.2){
176  for (Float_t z=-250; z<250;z+=20){
177  for (Int_t ifor=0;ifor<ndim;ifor++){
178  x[ifor]= ((TFormula*)formulas->At(ifor))->Eval(r/250.,fi,z/250.);
179  }
180  fitR->AddPoint(x,AliMagFDraw::GetBr(r,fi,z,index));
181  fitRFI->AddPoint(x,AliMagFDraw::GetBrfi(r,fi,z,index));
182  fitZ->AddPoint(x,AliMagFDraw::GetBz(r,fi,z,index));
183  }
184  }
185  }
186  fitR->Eval();
187  fitRFI->Eval();
188  fitZ->Eval();
189  TObjArray *res = new TObjArray;
190  res->AddAt(fitR,0);
191  res->AddAt(fitRFI,1);
192  res->AddAt(fitZ,2);
193  printf("\tchi2\tn\tRMS\n");
194  printf("\t%f\t%d\t%f\n",fitR->GetChisquare(),fitR->GetNpoints(),TMath::Sqrt(fitR->GetChisquare()/fitR->GetNpoints()));
195  printf("\t%f\t%d\t%f\n",fitRFI->GetChisquare(),fitRFI->GetNpoints(),TMath::Sqrt(fitRFI->GetChisquare()/fitRFI->GetNpoints()));
196  printf("\t%f\t%d\t%f\n",fitZ->GetChisquare(),fitZ->GetNpoints(),TMath::Sqrt(fitZ->GetChisquare()/fitZ->GetNpoints()));
197 
198  TFormula * funBZ = new TFormula("funBZ",formula);
199  TFormula * funBR = new TFormula("funBR",formula);
200  TFormula * funBRFI = new TFormula("funBRFI",formula);
201  TVectorD vec;
202  fitR->GetParameters(vec);
203  funBR->SetParameters(vec.GetMatrixArray());
204  fitRFI->GetParameters(vec);
205  funBRFI->SetParameters(vec.GetMatrixArray());
206  fitZ->GetParameters(vec);
207  funBZ->SetParameters(vec.GetMatrixArray());
208 
209  return res;
210 }
211 
212 
213 TString MakeString(){
214  TString str="";
215  {
216  Int_t counter=0;
217  for (Int_t ix=0;ix<3;ix++)
218  for (Int_t iz=0;iz<3;iz++){
219  if (ix+iz>0) {
220  str+=Form("x^%d*z^%d++",ix,iz);
221  printf("x^%d*z^%d++\n",ix,iz);
222  }
223  for (Int_t iy=1;iy<4;iy++){
224  if (ix+iz+iy==0) continue;
225  str+=Form("x^%d*z^%d*sin(y*%d)++",ix,iz,iy);
226  str+=Form("x^%d*z^%d*cos(y*%d)++",ix,iz,iy);
227  printf( "x^%d*z^%d*sin(y*%d)++\n",ix,iz,iy);
228  printf( "x^%d*z^%d*cos(y*%d)++\n",ix,iz,iy);
229  counter++;
230  }
231  }
232  str[str.Length()-2]=0;
233  }
234  return str;
235 }
236 
237 
238 
239 /*
240 
241 TObjArray * array = AliMagFDraw::Fit("x",0)
242 //
243 chi2 n RMS
244 15.534116 9600 0.040226
245 4.422160 9600 0.021463
246 8.378622 9600 0.029543
247 
248 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);
249 chi2 n RMS
250 1.842443 9600 0.013854
251 0.957056 9600 0.009985
252 0.097799 9600 0.003192
253 
254 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);
255 
256 chi2 n RMS
257 1.564687 9600 0.012767
258 0.002291 9600 0.000489
259 0.097063 9600 0.003180
260 
261 
262 TObjArray * array = AliMagFDraw::Fit(MakeString(),0);
263 
264 chi2 n RMS
265 0.000303 9600 0.000178
266 0.000066 9600 0.000083
267 0.003110 9600 0.000569
268 
269 
270 }
271 */
272 
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
#define TObjArray
TVectorD vec
Definition: AnalyzeLaser.C:8
static Double_t GetBx(Double_t r, Double_t phi, Double_t z, Int_t index=0)
virtual void Field(const Double_t *x, Double_t *b)
Definition: AliMagF.cxx:280
static TObjArray * Fit(const char *formula, Int_t index=0)
static void RegisterField(Int_t index, AliMagF *magf)
Definition: AliMagFDraw.cxx:86
static Double_t GetBrfi(Double_t r, Double_t phi, Double_t z, Int_t index=0)
static Double_t GetBy(Double_t r, Double_t phi, Double_t z, Int_t index=0)
static AliMagF::BMap_t mag
static Double_t GetBz(Double_t r, Double_t phi, Double_t z, Int_t index=0)
Definition: AliMagFDraw.cxx:92
TString MakeString()
static Double_t GetBr(Double_t r, Double_t phi, Double_t z, Int_t index=0)
void res(Char_t i)
Definition: Resolution.C:2
static TObjArray fgArray
Definition: AliMagFDraw.cxx:74