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