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
  58 void AliTPCExB::TestExB(const char* fileName) {
  59   //
  60   // Test ExB  -
  61   // Dump the filed and corrections to the tree in file fileName
  62   //
  63   //
  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.)
  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];
  73         Double_t d[3];
  74         Correct(x,d);
  75         Double_t rd=TMath::Sqrt(d[0]*d[0]+d[1]*d[1]);
  76         Double_t dr=r-rd;
  77         Double_t phi=TMath::ATan2(x[1],x[0]);
  78         Double_t phid=TMath::ATan2(d[1],d[0]);
  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];
  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";
 125       }
 126 }
 127
 128
 129
 130 Double_t AliTPCExB::GetDr(Double_t r, Double_t phi, Double_t z, Double_t bz){
 131   //
 132   // Static function
 133   // Posibble to us it for visualization
 134   //
 135   //
 136   AliTPCExB *exb = AliTPCcalibDB::GetExB(bz,kFALSE);
 137   if (!exb) return 0;
 138   Double_t pos0[3] = {r*TMath::Cos(phi), r*TMath::Sin(phi),z};
 139   Double_t pos1[3];
 140   exb->Correct(pos0,pos1);
 141   Double_t dx=pos1[0]-pos0[0];
 142   Double_t dy=pos1[1]-pos0[1];
 143   //  Double_t dz=pos1[2]-pos0[2];
 144   // return TMath::Sqrt(dx*dx+dy*dy);
 145   Float_t dr = (dx*pos0[0]+dy*pos0[1])/r;
 146   return dr;
 147 }
 148
 149
 150 Double_t AliTPCExB::GetDrphi(Double_t r, Double_t phi, Double_t z, Double_t bz){
 151   //
 152   //
 153   //
 154   AliTPCExB *exb = AliTPCcalibDB::GetExB(bz,kFALSE);
 155   if (!exb) return 0;
 156   Double_t pos0[3] = {r*TMath::Cos(phi), r*TMath::Sin(phi),z};
 157   Double_t pos1[3];
 158   exb->Correct(pos0,pos1);
 159   Double_t dphi=TMath::ATan2(pos1[1],pos1[0])-TMath::ATan2(pos0[1],pos0[0]);
 160   return r*dphi;
 161
 162 }
 163
 164
 165 Double_t AliTPCExB::GetDphi(Double_t r, Double_t phi, Double_t z, Double_t bz){
 166   //
 167   //
 168   //
 169   AliTPCExB *exb = AliTPCcalibDB::GetExB(bz,kFALSE);
 170   if (!exb) return 0;
 171   Double_t pos0[3] = {r*TMath::Cos(phi), r*TMath::Sin(phi),z};
 172   Double_t pos1[3];
 173   exb->Correct(pos0,pos1);
 174   Double_t dphi=TMath::ATan2(pos1[1],pos1[0])-TMath::ATan2(pos0[1],pos0[0]);
 175   return dphi;
 176
 177 }
 178
 179 Double_t AliTPCExB::GetDz(Double_t r, Double_t phi, Double_t z, Double_t bz){
 180   //
 181   //
 182   //
 183   AliTPCExB *exb = AliTPCcalibDB::GetExB(bz,kFALSE);
 184   if (!exb) return 0;
 185   Double_t pos0[3] = {r*TMath::Cos(phi), r*TMath::Sin(phi),z};
 186   Double_t pos1[3];
 187   exb->Correct(pos0,pos1);
 188   Double_t dz=pos1[2]-pos0[2];
 189   return dz;
 190 }
 191
 192 //
 193 // Magnetic field
 194 //
 195
 196
 197
 198
 199 void 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
 208 Double_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
 221 Double_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
 234 Double_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
 249 Double_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
 264 Double_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
 282 Double_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
 298 Double_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
 314 Double_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
 331 Double_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
 347 Double_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
 364 Double_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
 386 Double_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