TPC/AliTPCExB.cxx

   1 #include "AliTPCExB.h"
   2 #include "TMath.h"
   3 #include "TTreeStream.h"
   4 #include "AliMagF.h"
   5 #include "TLinearFitter.h"
   6
   7
   8 //
   9 // Abstract class for ExB effect parameterization
  10 //
  11 //
  12 //
  13 // The ExB correction map is stored in the calib DB
  14 // The lookup can be dumped to the tree:
  15 /*
  16
  17    //
  18   char *storage = "local://OCDBres"
  19   Int_t RunNumber=0;
  20   AliCDBManager::Instance()->SetDefaultStorage(storage);
  21   AliCDBManager::Instance()->SetRun(RunNumber)
  22   AliTPCExBFirst * exb = AliTPCcalibDB::Instance()->GetExB();
  23   //
  24   // See example macro $ALICE_ROOT/TPC/macros/AliTPCExBdraw.C
  25   //
  26   .L $ALICE_ROOT/TPC/macros/AliTPCExBdraw.C
  27   Draw(0)
  28
  29
  30
  31
  32 */
  33
  34 AliTPCExB* AliTPCExB::fgInstance = 0;
  35
  36 TObjArray   AliTPCExB::fgArray;
  37
  38
  39
  40 ClassImp(AliTPCExB)
  41
  42 AliTPCExB::AliTPCExB():
  43   TObject(),
  44   fMatBrBz(0),       //param matrix Br/Bz
  45   fMatBrfiBz(0),     //param matrix Br/Bz
  46   fMatBrBzI0(0),     //param matrix Br/Bz integral  z>0
  47   fMatBrBzI1(0),     //param matrix Br/Bz integral  z<0
  48   fMatBrfiBzI0(0),   //param matrix Br/Bz integral  z>0
  49   fMatBrfiBzI1(0)    //param matrix Br/Bz integral  z<0
  50 {
  51   //
  52   // default constructor
  53   //
  54 }
  55
  56
  57 void AliTPCExB::TestExB(const char* fileName) {
  58   //
  59   // Test ExB  -
  60   // Dump the filed and corrections to the tree in file fileName
  61   //
  62   //
  63   TTreeSRedirector ts(fileName);
  64   Double_t x[3];
  65   for (x[0]=-250.;x[0]<=250.;x[0]+=10.)
  66     for (x[1]=-250.;x[1]<=250.;x[1]+=10.)
  67       for (x[2]=-250.;x[2]<=250.;x[2]+=20.) {
  68         Double_t r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
  69         if (r<20) continue;
  70         if (r>260) continue;
  71         Double_t z = x[2];
  72         Double_t d[3];
  73         Correct(x,d);
  74         Double_t rd=TMath::Sqrt(d[0]*d[0]+d[1]*d[1]);
  75         Double_t dr=r-rd;
  76         Double_t phi=TMath::ATan2(x[1],x[0]);
  77         Double_t phid=TMath::ATan2(d[1],d[0]);
  78         Double_t dphi=phi-phid;
  79         if (dphi<0.) dphi+=TMath::TwoPi();
  80         if (dphi>TMath::Pi()) dphi=TMath::TwoPi()-dphi;
  81         Double_t drphi=r*dphi;
  82         Double_t dx=x[0]-d[0];
  83         Double_t dy=x[1]-d[1];
  84         Double_t dz=x[2]-d[2];
  85         //
  86         Double_t bx    = GetBx(r,phi,z,0);
  87         Double_t by    = GetBy(r,phi,z,0);
  88         Double_t bz    = GetBz(r,phi,z,0);
  89         Double_t br    = GetBr(r,phi,z,0);
  90         Double_t brfi  = GetBrfi(r,phi,z,0);
  91         //
  92         Double_t bxi    = GetBxI(r,phi,z,0);
  93         Double_t byi    = GetByI(r,phi,z,0);
  94         Double_t bzi    = GetBzI(r,phi,z,0);
  95         Double_t bri    = GetBrI(r,phi,z,0);
  96         Double_t brfii  = GetBrfiI(r,phi,z,0);
  97
  98         ts<<"positions"<<
  99           "x0="<<x[0]<<
 100           "x1="<<x[1]<<
 101           "x2="<<x[2]<<
 102           "dx="<<dx<<
 103           "dy="<<dy<<
 104           "dz="<<dz<<
 105           "r="<<r<<
 106           "phi="<<phi<<
 107           "dr="<<dr<<
 108           "drphi="<<drphi<<
 109           //
 110           // B-Field
 111           //
 112           "bx="<<bx<<
 113           "by="<<by<<
 114           "bz="<<bz<<
 115           "br="<< br<<
 116           "brfi="<<brfi<<
 117           // B-field integ
 118           "bxi="<<bxi<<
 119           "byi="<<byi<<
 120           "bzi="<<bzi<<
 121           "bri="<< bri<<
 122           "brfii="<<brfii<<
 123           "\n";
 124       }
 125 }
 126
 127
 128
 129 Double_t AliTPCExB::GetDr(Double_t r, Double_t phi, Double_t z){
 130   //
 131   // Static function
 132   // Posibble to us it for visualization
 133   //
 134   //
 135   if (!fgInstance) return 0;
 136   Double_t pos0[3] = {r*TMath::Cos(phi), r*TMath::Sin(phi),z};
 137   Double_t pos1[3];
 138   fgInstance->Correct(pos0,pos1);
 139   Double_t dx=pos1[0]-pos0[0];
 140   Double_t dy=pos1[1]-pos0[1];
 141   //  Double_t dz=pos1[2]-pos0[2];
 142   return TMath::Sqrt(dx*dx+dy*dy);
 143 }
 144
 145
 146 Double_t AliTPCExB::GetDrphi(Double_t r, Double_t phi, Double_t z){
 147   //
 148   //
 149   //
 150   if (!fgInstance) return 0;
 151   Double_t pos0[3] = {r*TMath::Cos(phi), r*TMath::Sin(phi),z};
 152   Double_t pos1[3];
 153   fgInstance->Correct(pos0,pos1);
 154   Double_t dphi=TMath::ATan2(pos1[1],pos1[0])-TMath::ATan2(pos0[1],pos0[0]);
 155   return r*dphi;
 156
 157 }
 158
 159
 160 Double_t AliTPCExB::GetDphi(Double_t r, Double_t phi, Double_t z){
 161   //
 162   //
 163   //
 164   if (!fgInstance) return 0;
 165   Double_t pos0[3] = {r*TMath::Cos(phi), r*TMath::Sin(phi),z};
 166   Double_t pos1[3];
 167   fgInstance->Correct(pos0,pos1);
 168   Double_t dphi=TMath::ATan2(pos1[1],pos1[0])-TMath::ATan2(pos0[1],pos0[0]);
 169   return dphi;
 170
 171 }
 172
 173 Double_t AliTPCExB::GetDz(Double_t r, Double_t phi, Double_t z){
 174   //
 175   //
 176   //
 177   if (!fgInstance) return 0;
 178   Double_t pos0[3] = {r*TMath::Cos(phi), r*TMath::Sin(phi),z};
 179   Double_t pos1[3];
 180   fgInstance->Correct(pos0,pos1);
 181   Double_t dz=pos1[2]-pos0[2];
 182   return dz;
 183 }
 184
 185 //
 186 // Magnetic field
 187 //
 188
 189
 190
 191
 192 void AliTPCExB::RegisterField(Int_t index, AliMagF * magf){
 193   //
 194   // add the filed to the list
 195   //
 196   fgArray.AddAt(magf,index);
 197 }
 198
 199
 200
 201 Double_t AliTPCExB::GetBx(Double_t r, Double_t phi, Double_t z,Int_t index){
 202   //
 203   //
 204   //
 205   AliMagF *mag = (AliMagF*)fgArray.At(index);
 206   if (!mag) return 0;
 207   Float_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
 208   //  xyz[1]+=30;
 209   Float_t bxyz[3];
 210   mag->Field(xyz,bxyz);
 211   return bxyz[0];
 212 }
 213
 214 Double_t AliTPCExB::GetBy(Double_t r, Double_t phi, Double_t z,Int_t index){
 215   //
 216   //
 217   //
 218   AliMagF *mag = (AliMagF*)fgArray.At(index);
 219   if (!mag) return 0;
 220   Float_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
 221   //  xyz[1]+=30;
 222   Float_t bxyz[3];
 223   mag->Field(xyz,bxyz);
 224   return bxyz[1];
 225 }
 226
 227 Double_t AliTPCExB::GetBz(Double_t r, Double_t phi, Double_t z,Int_t index){
 228   //
 229   //
 230   //
 231   AliMagF *mag = (AliMagF*)fgArray.At(index);
 232   if (!mag) return 0;
 233   Float_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
 234   //  xyz[1]+=30;
 235   Float_t bxyz[3];
 236   mag->Field(xyz,bxyz);
 237   return bxyz[2];
 238 }
 239
 240
 241
 242 Double_t AliTPCExB::GetBr(Double_t r, Double_t phi, Double_t z,Int_t index){
 243   //
 244   //
 245   //
 246   AliMagF *mag = (AliMagF*)fgArray.At(index);
 247   if (!mag) return 0;
 248   Float_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
 249   //xyz[1]+=30;
 250   Float_t bxyz[3];
 251   mag->Field(xyz,bxyz);
 252   if (r==0) return 0;
 253   Float_t br = (bxyz[0]*xyz[0]+bxyz[1]*xyz[1])/r;
 254   return br;
 255 }
 256
 257 Double_t AliTPCExB::GetBrfi(Double_t r, Double_t phi, Double_t z,Int_t index){
 258   //
 259   //
 260   //
 261   AliMagF *mag = (AliMagF*)fgArray.At(index);
 262   if (!mag) return 0;
 263   Float_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
 264   //xyz[1]+=30;
 265   Float_t bxyz[3];
 266   mag->Field(xyz,bxyz);
 267   if (r==0) return 0;
 268   Float_t br = (-bxyz[0]*xyz[1]+bxyz[1]*xyz[0])/r;
 269   return br;
 270 }
 271
 272
 273
 274
 275 Double_t AliTPCExB::GetBxI(Double_t r, Double_t phi, Double_t z,Int_t index)
 276 {
 277   Double_t sumf  =0;
 278   if (z>0 &&z<250){
 279     for (Float_t zi=z;zi<250;zi+=5){
 280       sumf+=GetBx(r,phi,zi,index)/GetBz(r,phi,zi,index);
 281     }
 282   }
 283   if (z<0 &&z>-250){
 284     for (Float_t zi=z;zi>-250;zi-=5){
 285       sumf+=GetBx(r,phi,zi,index)/GetBz(r,phi,zi,index);
 286     }
 287   }
 288   return sumf*5;
 289 }
 290
 291 Double_t AliTPCExB::GetByI(Double_t r, Double_t phi, Double_t z,Int_t index)
 292 {
 293   Double_t sumf =0;
 294   if (z>0 &&z<250){
 295     for (Float_t zi=z;zi<250;zi+=5){
 296       sumf+=GetBy(r,phi,zi,index)/GetBz(r,phi,zi,index);
 297     }
 298   }
 299   if (z<0 &&z>-250){
 300     for (Float_t zi=z;zi>-250;zi-=5){
 301       sumf+=GetBy(r,phi,zi,index)/GetBz(r,phi,zi,index);
 302     }
 303   }
 304   return sumf*5;
 305 }
 306
 307 Double_t AliTPCExB::GetBzI(Double_t r, Double_t phi, Double_t z,Int_t index)
 308 {
 309   Double_t sumf =0;
 310   if (z>0 &&z<250){
 311     for (Float_t zi=z;zi<250;zi+=5){
 312       sumf+=GetBz(r,phi,zi,index);
 313     }
 314   }
 315   if (z<0 &&z>-250){
 316     for (Float_t zi=z;zi>-250;zi-=5){
 317       sumf+=GetBz(r,phi,zi,index);
 318     }
 319   }
 320   return sumf*5;
 321 }
 322
 323
 324 Double_t AliTPCExB::GetBrI(Double_t r, Double_t phi, Double_t z,Int_t index)
 325 {
 326   Double_t sumf =0;
 327   if (z>0 &&z<250){
 328     for (Float_t zi=z;zi<250;zi+=5){
 329       sumf+=GetBr(r,phi,zi,index)/GetBz(r,phi,zi,index);
 330     }
 331   }
 332   if (z<0 &&z>-250){
 333     for (Float_t zi=z;zi>-250;zi-=5){
 334       sumf+=GetBr(r,phi,zi,index)/GetBz(r,phi,zi,index);
 335     }
 336   }
 337   return sumf*5;
 338 }
 339
 340 Double_t AliTPCExB::GetBrfiI(Double_t r, Double_t phi, Double_t z,Int_t index)
 341 {
 342   Double_t sumf =0;
 343   if (z>0 &&z<250){
 344     for (Float_t zi=z;zi<250;zi+=5.){
 345       sumf+=GetBrfi(r,phi,zi,index)/GetBz(r,phi,zi,index);
 346     }
 347   }
 348   if (z<0 &&z>-250){
 349     for (Float_t zi=z;zi>-250;zi-=5){
 350       sumf+=GetBrfi(r,phi,zi,index)/GetBz(r,phi,zi,index);
 351     }
 352   }
 353   return sumf*5;
 354 }
 355
 356
 357 Double_t AliTPCExB::Eval(Int_t type, Double_t r, Double_t phi, Double_t z){
 358   //
 359   // Evaluate parameterization
 360   //
 361   // br integral param
 362   if (type==0) {
 363     if (z>0 && fMatBrBzI0) return EvalMat(*fMatBrBzI0,r,phi,z);
 364     if (z<0 && fMatBrBzI1) return EvalMat(*fMatBrBzI1,r,phi,z);
 365   }
 366   // brfi integral param
 367   if (type==1) {
 368     if (z>0 && fMatBrfiBzI0) return EvalMat(*fMatBrfiBzI0,r,phi,z);
 369     if (z<0 && fMatBrfiBzI1) return EvalMat(*fMatBrfiBzI1,r,phi,z);
 370   }
 371   // brbz param
 372   if (type==2 && fMatBrBz) return EvalMat(*fMatBrBz,r,phi,z);
 373   // brfibz param
 374   if (type==3 && fMatBrfiBz) return EvalMat(*fMatBrfiBz,r,phi,z);
 375   return 0;
 376 }
 377
 378
 379 Double_t AliTPCExB::EvalMat(TVectorD &vec, Double_t r, Double_t phi, Double_t z){
 380   //
 381   // Evaluate taylor expansion in r,phi,z
 382   //
 383   // Variables
 384   //tree->SetAlias("sa","sin(phi+0.0)");
 385   //tree->SetAlias("ca","cos(phi+0.0)");
 386   //tree->SetAlias("sa2","sin(phi*2+0.0)");
 387   //tree->SetAlias("ca2","cos(phi*2+0.0)");
 388   //tree->SetAlias("zn","(x2/250.)");
 389   //tree->SetAlias("rn","(r/250.)")
 390   //   TString fstringSym="";
 391   //   //
 392   //   fstringSym+="zn++";
 393   //   fstringSym+="rn++";
 394   //   fstringSym+="zn*rn++";
 395   //   fstringSym+="zn*zn++";
 396   //   fstringSym+="zn*zn*rn++";
 397   //   fstringSym+="zn*rn*rn++";
 398   //   //
 399   //   fstringSym+="sa++";
 400   //   fstringSym+="ca++";
 401   //   fstringSym+="ca2++";
 402   //   fstringSym+="sa2++";
 403   //   fstringSym+="ca*zn++";
 404   //   fstringSym+="sa*zn++";
 405   //   fstringSym+="ca2*zn++";
 406   //   fstringSym+="sa2*zn++";
 407   //   fstringSym+="ca*zn*zn++";
 408   //   fstringSym+="sa*zn*zn++";
 409   //   fstringSym+="ca*zn*rn++";
 410   //   fstringSym+="sa*zn*rn++";
 411
 412
 413   Double_t sa  = TMath::Sin(phi);
 414   Double_t ca  = TMath::Cos(phi);
 415   Double_t sa2 = TMath::Sin(phi*2);
 416   Double_t ca2 = TMath::Cos(phi*2);
 417   Double_t zn  = z/250.;
 418   Double_t rn  = r/250.;
 419   Int_t  ipoint=0;
 420   Double_t res = vec[ipoint++];
 421   res+=vec[ipoint++]*zn;
 422   res+=vec[ipoint++]*rn;
 423   res+=vec[ipoint++]*zn*rn;
 424   res+=vec[ipoint++]*zn*zn;
 425   res+=vec[ipoint++]*zn*zn*rn;
 426   res+=vec[ipoint++]*zn*rn*rn;
 427   //
 428   res+=vec[ipoint++]*sa;
 429   res+=vec[ipoint++]*ca;
 430   res+=vec[ipoint++]*ca2;
 431   res+=vec[ipoint++]*sa2;
 432   res+=vec[ipoint++]*ca*zn;
 433   res+=vec[ipoint++]*sa*zn;
 434   res+=vec[ipoint++]*ca2*zn;
 435   res+=vec[ipoint++]*sa2*zn;
 436   res+=vec[ipoint++]*ca*zn*zn;
 437   res+=vec[ipoint++]*sa*zn*zn;
 438   res+=vec[ipoint++]*ca*zn*rn;
 439   res+=vec[ipoint++]*sa*zn*rn;
 440   return res;
 441 }
 442
 443
 444
 445
 446
 447
 448
 449
 450
 451
 452
 453
 454
 455
 456 /*
 457
 458  AliTPCExB draw;
 459  draw.RegisterField(0,new AliMagWrapCheb("Maps","Maps", 2, 1., 10., AliMagWrapCheb::k5kG));
 460  draw.RegisterField(1,new AliMagFMaps("Maps","Maps", 2, 1., 10., 2));
 461
 462  TF2 fbz_rz_0pi("fbz_rz_0pi","AliTPCExB::GetBz(x,0*pi,y)",0,250,-250,250);
 463  fbz_rz_0pi->Draw("surf2");
 464
 465  TF1 fbz_z_90_00pi("fbz_z_90_00pi","AliTPCExB::GetBz(90,0*pi,x)",-250,250);
 466   TF1 fbz_z_90_05pi("fbz_z_90_05pi","AliTPCExB::GetBz(90,0.5*pi,x)",-250,250);
 467   TF1 fbz_z_90_10pi("fbz_z_90_10pi","AliTPCExB::GetBz(90,1.0*pi,x)",-250,250);
 468   TF1 fbz_z_90_15pi("fbz_z_90_15pi","AliTPCExB::GetBz(90,1.5*pi,x)",-250,250);
 469   fbz_z_90_00pi->SetLineColor(2);
 470   fbz_z_90_05pi->SetLineColor(3);
 471   fbz_z_90_10pi->SetLineColor(4);
 472   fbz_z_90_15pi->SetLineColor(5);
 473   fbz_z_90_00pi->Draw()
 474   fbz_z_90_05pi->Draw("same")
 475   fbz_z_90_15pi->Draw("same")
 476   fbz_z_90_10pi->Draw("same")
 477
 478
 479   TF1 fbr_z_90_00pi("fbz_z_90_00pi","AliTPCExB::GetBr(90,0*pi,x)",-250,250);
 480   TF1 fbr_z_90_05pi("fbz_z_90_05pi","AliTPCExB::GetBr(90,0.5*pi,x)",-250,250);
 481   TF1 fbr_z_90_10pi("fbz_z_90_10pi","AliTPCExB::GetBr(90,1.0*pi,x)",-250,250);
 482   TF1 fbr_z_90_15pi("fbz_z_90_15pi","AliTPCExB::GetBr(90,1.5*pi,x)",-250,250);
 483   fbr_z_90_00pi->SetLineColor(2);
 484   fbr_z_90_05pi->SetLineColor(3);
 485   fbr_z_90_10pi->SetLineColor(4);
 486   fbr_z_90_15pi->SetLineColor(5);
 487   fbr_z_90_00pi->Draw()
 488   fbr_z_90_05pi->Draw("same")
 489   fbr_z_90_15pi->Draw("same")
 490   fbr_z_90_10pi->Draw("same")
 491
 492   //
 493   TF2 fbz_xy_0z("fbz_xy_0z","AliTPCExB::GetBz(sqrt(x^2+y^2),atan2(y,x),0)",-250,250,-250,250);
 494   fbz_xy_0z.SetNpy(100);
 495   fbz_xy_0z.SetNpx(100);
 496   fbz_xy_0z->Draw("colz");
 497   //
 498   TF2 fbz_xy_250z("fbz_xy_250z","AliTPCExB::GetBz(sqrt(x^2+y^2),atan2(y,x),250)",-250,250,-250,250);
 499   fbz_xy_250z.SetNpy(100);
 500   fbz_xy_250z.SetNpx(100)
 501   fbz_xy_250z->Draw("colz");
 502   //
 503    TF2 fbz_xy_m250z("fbz_xy_m250z","AliTPCExB::GetBz(sqrt(x^2+y^2),atan2(y,x),-250)",-250,250,-250,250);
 504   fbz_xy_m250z.SetNpy(100);
 505   fbz_xy_m250z.SetNpx(100)
 506   fbz_xy_m250z->Draw("colz");
 507   //
 508
 509
 510
 511
 512
 513 */
 514
 515