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