Adding fit parameterization of the magnetic field inside of the TPC

(Marian)

diff --git a/TPC/AliTPCExB.cxx b/TPC/AliTPCExB.cxx
index d528fe2542a80280c3c798acbdb924e6dbb0bb4f..38865d996ab1727aaafb329363f0b95ce5d5cf8f 100644 (file)
--- a/TPC/AliTPCExB.cxx
+++ b/TPC/AliTPCExB.cxx
@@ -1,6 +1,8 @@
 #include "AliTPCExB.h"
 #include "TMath.h"
 #include "TTreeStream.h"
+#include "AliMagF.h"
+#include "TLinearFitter.h"
 
 
 //
@@ -31,25 +33,48 @@
 
 AliTPCExB* AliTPCExB::fgInstance = 0;
 
+TObjArray   AliTPCExB::fgArray;
+
+
+
 ClassImp(AliTPCExB)
 
+AliTPCExB::AliTPCExB():
+  TObject(),
+  fMatBrBz(0),       //param matrix Br/Bz
+  fMatBrfiBz(0),     //param matrix Br/Bz
+  fMatBrBzI0(0),     //param matrix Br/Bz integral  z>0 
+  fMatBrBzI1(0),     //param matrix Br/Bz integral  z<0 
+  fMatBrfiBzI0(0),   //param matrix Br/Bz integral  z>0 
+  fMatBrfiBzI1(0)    //param matrix Br/Bz integral  z<0
+{
+  //
+  // default constructor
+  //
+}
+
 
 void AliTPCExB::TestExB(const char* fileName) {
   //
-  // well, as the name sais...
+  // Test ExB  - 
+  // Dump the filed and corrections to the tree in file fileName  
   //
+  // 
   TTreeSRedirector ts(fileName);
   Double_t x[3];
   for (x[0]=-250.;x[0]<=250.;x[0]+=10.)
     for (x[1]=-250.;x[1]<=250.;x[1]+=10.)
-      for (x[2]=-250.;x[2]<=250.;x[2]+=10.) {
+      for (x[2]=-250.;x[2]<=250.;x[2]+=20.) {
+       Double_t r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
+       if (r<20) continue;
+       if (r>260) continue;
+       Double_t z = x[2];
        Double_t d[3];
        Correct(x,d);
-       Double_t r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
        Double_t rd=TMath::Sqrt(d[0]*d[0]+d[1]*d[1]);
        Double_t dr=r-rd;
-       Double_t phi=TMath::ATan2(x[0],x[1]);
-       Double_t phid=TMath::ATan2(d[0],d[1]);
+       Double_t phi=TMath::ATan2(x[1],x[0]);
+       Double_t phid=TMath::ATan2(d[1],d[0]);
        Double_t dphi=phi-phid;
        if (dphi<0.) dphi+=TMath::TwoPi();
        if (dphi>TMath::Pi()) dphi=TMath::TwoPi()-dphi;
@@ -57,18 +82,45 @@ void AliTPCExB::TestExB(const char* fileName) {
        Double_t dx=x[0]-d[0];
        Double_t dy=x[1]-d[1];
        Double_t dz=x[2]-d[2];
-       ts<<"positions"
-         <<"x0="<<x[0]
-         <<"x1="<<x[1]
-         <<"x2="<<x[2]
-         <<"dx="<<dx
-         <<"dy="<<dy
-         <<"dz="<<dz
-         <<"r="<<r
-         <<"phi="<<phi
-         <<"dr="<<dr
-         <<"drphi="<<drphi
-         <<"\n";
+       //
+       Double_t bx    = GetBx(r,phi,z,0);
+       Double_t by    = GetBy(r,phi,z,0);
+       Double_t bz    = GetBz(r,phi,z,0);
+       Double_t br    = GetBr(r,phi,z,0);
+       Double_t brfi  = GetBrfi(r,phi,z,0);
+       //
+       Double_t bxi    = GetBxI(r,phi,z,0);
+       Double_t byi    = GetByI(r,phi,z,0);
+       Double_t bzi    = GetBzI(r,phi,z,0);
+       Double_t bri    = GetBrI(r,phi,z,0);
+       Double_t brfii  = GetBrfiI(r,phi,z,0);
+
+       ts<<"positions"<<
+         "x0="<<x[0]<<
+         "x1="<<x[1]<<
+         "x2="<<x[2]<<
+         "dx="<<dx<<
+         "dy="<<dy<<
+         "dz="<<dz<<
+         "r="<<r<<
+         "phi="<<phi<<
+         "dr="<<dr<<
+         "drphi="<<drphi<<
+         //
+         // B-Field
+         //
+         "bx="<<bx<<
+         "by="<<by<<
+         "bz="<<bz<<
+         "br="<< br<<
+         "brfi="<<brfi<<
+         // B-field integ
+         "bxi="<<bxi<<
+         "byi="<<byi<<
+         "bzi="<<bzi<<
+         "bri="<< bri<<
+         "brfii="<<brfii<<
+         "\n";
       }
 }
 
@@ -129,3 +181,335 @@ Double_t AliTPCExB::GetDz(Double_t r, Double_t phi, Double_t z){
   Double_t dz=pos1[2]-pos0[2];
   return dz;  
 }
+
+//
+// Magnetic field
+//
+
+
+
+
+void AliTPCExB::RegisterField(Int_t index, AliMagF * magf){
+  //
+  // add the filed to the list
+  //
+  fgArray.AddAt(magf,index);
+}
+
+
+
+Double_t AliTPCExB::GetBx(Double_t r, Double_t phi, Double_t z,Int_t index){
+  //
+  // 
+  //
+  AliMagF *mag = (AliMagF*)fgArray.At(index);
+  if (!mag) return 0;
+  Float_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
+  //  xyz[1]+=30;
+  Float_t bxyz[3];
+  mag->Field(xyz,bxyz);
+  return bxyz[0];
+}  
+
+Double_t AliTPCExB::GetBy(Double_t r, Double_t phi, Double_t z,Int_t index){
+  //
+  // 
+  //
+  AliMagF *mag = (AliMagF*)fgArray.At(index);
+  if (!mag) return 0;
+  Float_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
+  //  xyz[1]+=30;
+  Float_t bxyz[3];
+  mag->Field(xyz,bxyz);
+  return bxyz[1];
+}  
+
+Double_t AliTPCExB::GetBz(Double_t r, Double_t phi, Double_t z,Int_t index){
+  //
+  // 
+  //
+  AliMagF *mag = (AliMagF*)fgArray.At(index);
+  if (!mag) return 0;
+  Float_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
+  //  xyz[1]+=30;
+  Float_t bxyz[3];
+  mag->Field(xyz,bxyz);
+  return bxyz[2];
+}  
+
+
+
+Double_t AliTPCExB::GetBr(Double_t r, Double_t phi, Double_t z,Int_t index){
+  //
+  // 
+  //
+  AliMagF *mag = (AliMagF*)fgArray.At(index);
+  if (!mag) return 0;
+  Float_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
+  //xyz[1]+=30;
+  Float_t bxyz[3];
+  mag->Field(xyz,bxyz);
+  if (r==0) return 0;
+  Float_t br = (bxyz[0]*xyz[0]+bxyz[1]*xyz[1])/r;
+  return br;
+}  
+
+Double_t AliTPCExB::GetBrfi(Double_t r, Double_t phi, Double_t z,Int_t index){
+  //
+  // 
+  //
+  AliMagF *mag = (AliMagF*)fgArray.At(index);
+  if (!mag) return 0;
+  Float_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
+  //xyz[1]+=30;
+  Float_t bxyz[3];
+  mag->Field(xyz,bxyz);
+  if (r==0) return 0;
+  Float_t br = (-bxyz[0]*xyz[1]+bxyz[1]*xyz[0])/r;
+  return br;
+}  
+
+
+
+
+Double_t AliTPCExB::GetBxI(Double_t r, Double_t phi, Double_t z,Int_t index)
+{
+  Double_t sumf  =0;
+  if (z>0 &&z<250){
+    for (Float_t zi=z;zi<250;zi+=5){
+      sumf+=GetBx(r,phi,zi,index)/GetBz(r,phi,zi,index);
+    }
+  }
+  if (z<0 &&z>-250){
+    for (Float_t zi=z;zi>-250;zi-=5){
+      sumf+=GetBx(r,phi,zi,index)/GetBz(r,phi,zi,index);
+    }
+  }
+  return sumf*5;
+}
+
+Double_t AliTPCExB::GetByI(Double_t r, Double_t phi, Double_t z,Int_t index)
+{
+  Double_t sumf =0;
+  if (z>0 &&z<250){
+    for (Float_t zi=z;zi<250;zi+=5){
+      sumf+=GetBy(r,phi,zi,index)/GetBz(r,phi,zi,index);
+    }
+  }
+  if (z<0 &&z>-250){
+    for (Float_t zi=z;zi>-250;zi-=5){
+      sumf+=GetBy(r,phi,zi,index)/GetBz(r,phi,zi,index);
+    }
+  }
+  return sumf*5;
+}
+
+Double_t AliTPCExB::GetBzI(Double_t r, Double_t phi, Double_t z,Int_t index)
+{
+  Double_t sumf =0;
+  if (z>0 &&z<250){
+    for (Float_t zi=z;zi<250;zi+=5){
+      sumf+=GetBz(r,phi,zi,index);
+    }
+  }
+  if (z<0 &&z>-250){
+    for (Float_t zi=z;zi>-250;zi-=5){
+      sumf+=GetBz(r,phi,zi,index);
+    }
+  }
+  return sumf*5;
+}
+
+
+Double_t AliTPCExB::GetBrI(Double_t r, Double_t phi, Double_t z,Int_t index)
+{
+  Double_t sumf =0;
+  if (z>0 &&z<250){
+    for (Float_t zi=z;zi<250;zi+=5){
+      sumf+=GetBr(r,phi,zi,index)/GetBz(r,phi,zi,index);
+    }
+  }
+  if (z<0 &&z>-250){
+    for (Float_t zi=z;zi>-250;zi-=5){
+      sumf+=GetBr(r,phi,zi,index)/GetBz(r,phi,zi,index);
+    }
+  }
+  return sumf*5;
+}
+
+Double_t AliTPCExB::GetBrfiI(Double_t r, Double_t phi, Double_t z,Int_t index)
+{
+  Double_t sumf =0;
+  if (z>0 &&z<250){
+    for (Float_t zi=z;zi<250;zi+=5.){
+      sumf+=GetBrfi(r,phi,zi,index)/GetBz(r,phi,zi,index);
+    }
+  }
+  if (z<0 &&z>-250){
+    for (Float_t zi=z;zi>-250;zi-=5){
+      sumf+=GetBrfi(r,phi,zi,index)/GetBz(r,phi,zi,index);
+    }
+  }
+  return sumf*5;
+}
+
+
+Double_t AliTPCExB::Eval(Int_t type, Double_t r, Double_t phi, Double_t z){
+  //
+  // Evaluate parameterization
+  //
+  // br integral param 
+  if (type==0) {
+    if (z>0 && fMatBrBzI0) return EvalMat(*fMatBrBzI0,r,phi,z);
+    if (z<0 && fMatBrBzI1) return EvalMat(*fMatBrBzI1,r,phi,z);
+  }
+  // brfi integral param   
+  if (type==1) {
+    if (z>0 && fMatBrfiBzI0) return EvalMat(*fMatBrfiBzI0,r,phi,z);
+    if (z<0 && fMatBrfiBzI1) return EvalMat(*fMatBrfiBzI1,r,phi,z);
+  }
+  // brbz param
+  if (type==2 && fMatBrBz) return EvalMat(*fMatBrBz,r,phi,z);
+  // brfibz param
+  if (type==3 && fMatBrfiBz) return EvalMat(*fMatBrfiBz,r,phi,z);
+  return 0;
+}
+
+
+Double_t AliTPCExB::EvalMat(TVectorD &vec, Double_t r, Double_t phi, Double_t z){
+  //
+  // Evaluate taylor expansion in r,phi,z
+  //
+  // Variables  
+  //tree->SetAlias("sa","sin(phi+0.0)");
+  //tree->SetAlias("ca","cos(phi+0.0)");
+  //tree->SetAlias("sa2","sin(phi*2+0.0)");
+  //tree->SetAlias("ca2","cos(phi*2+0.0)");
+  //tree->SetAlias("zn","(x2/250.)");
+  //tree->SetAlias("rn","(r/250.)")
+  //   TString fstringSym="";
+  //   //  
+  //   fstringSym+="zn++";
+  //   fstringSym+="rn++";
+  //   fstringSym+="zn*rn++";
+  //   fstringSym+="zn*zn++";
+  //   fstringSym+="zn*zn*rn++";
+  //   fstringSym+="zn*rn*rn++";
+  //   //
+  //   fstringSym+="sa++";
+  //   fstringSym+="ca++";  
+  //   fstringSym+="ca2++";
+  //   fstringSym+="sa2++";
+  //   fstringSym+="ca*zn++";
+  //   fstringSym+="sa*zn++";
+  //   fstringSym+="ca2*zn++";
+  //   fstringSym+="sa2*zn++";
+  //   fstringSym+="ca*zn*zn++";
+  //   fstringSym+="sa*zn*zn++";
+  //   fstringSym+="ca*zn*rn++";
+  //   fstringSym+="sa*zn*rn++";
+
+
+  Double_t sa  = TMath::Sin(phi);
+  Double_t ca  = TMath::Cos(phi);
+  Double_t sa2 = TMath::Sin(phi*2);
+  Double_t ca2 = TMath::Cos(phi*2);
+  Double_t zn  = z/250.;
+  Double_t rn  = r/250.;
+  Int_t  ipoint=0;
+  Double_t res = vec[ipoint++];
+  res+=vec[ipoint++]*zn;
+  res+=vec[ipoint++]*rn;
+  res+=vec[ipoint++]*zn*rn;
+  res+=vec[ipoint++]*zn*zn;
+  res+=vec[ipoint++]*zn*zn*rn;
+  res+=vec[ipoint++]*zn*rn*rn;
+  //
+  res+=vec[ipoint++]*sa;
+  res+=vec[ipoint++]*ca;  
+  res+=vec[ipoint++]*ca2;
+  res+=vec[ipoint++]*sa2;
+  res+=vec[ipoint++]*ca*zn;
+  res+=vec[ipoint++]*sa*zn;
+  res+=vec[ipoint++]*ca2*zn;
+  res+=vec[ipoint++]*sa2*zn;
+  res+=vec[ipoint++]*ca*zn*zn;
+  res+=vec[ipoint++]*sa*zn*zn;
+  res+=vec[ipoint++]*ca*zn*rn;
+  res+=vec[ipoint++]*sa*zn*rn;
+  return res;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/*
+  
+ AliTPCExB draw;
+ draw.RegisterField(0,new AliMagWrapCheb("Maps","Maps", 2, 1., 10., AliMagWrapCheb::k5kG));
+ draw.RegisterField(1,new AliMagFMaps("Maps","Maps", 2, 1., 10., 2));
+
+ TF2 fbz_rz_0pi("fbz_rz_0pi","AliTPCExB::GetBz(x,0*pi,y)",0,250,-250,250);
+ fbz_rz_0pi->Draw("surf2");
+ 
+ TF1 fbz_z_90_00pi("fbz_z_90_00pi","AliTPCExB::GetBz(90,0*pi,x)",-250,250);
+  TF1 fbz_z_90_05pi("fbz_z_90_05pi","AliTPCExB::GetBz(90,0.5*pi,x)",-250,250);
+  TF1 fbz_z_90_10pi("fbz_z_90_10pi","AliTPCExB::GetBz(90,1.0*pi,x)",-250,250);
+  TF1 fbz_z_90_15pi("fbz_z_90_15pi","AliTPCExB::GetBz(90,1.5*pi,x)",-250,250);
+  fbz_z_90_00pi->SetLineColor(2);
+  fbz_z_90_05pi->SetLineColor(3);
+  fbz_z_90_10pi->SetLineColor(4);
+  fbz_z_90_15pi->SetLineColor(5);
+  fbz_z_90_00pi->Draw()
+  fbz_z_90_05pi->Draw("same")
+  fbz_z_90_15pi->Draw("same")
+  fbz_z_90_10pi->Draw("same")
+  
+
+  TF1 fbr_z_90_00pi("fbz_z_90_00pi","AliTPCExB::GetBr(90,0*pi,x)",-250,250);
+  TF1 fbr_z_90_05pi("fbz_z_90_05pi","AliTPCExB::GetBr(90,0.5*pi,x)",-250,250);
+  TF1 fbr_z_90_10pi("fbz_z_90_10pi","AliTPCExB::GetBr(90,1.0*pi,x)",-250,250);
+  TF1 fbr_z_90_15pi("fbz_z_90_15pi","AliTPCExB::GetBr(90,1.5*pi,x)",-250,250);
+  fbr_z_90_00pi->SetLineColor(2);
+  fbr_z_90_05pi->SetLineColor(3);
+  fbr_z_90_10pi->SetLineColor(4);
+  fbr_z_90_15pi->SetLineColor(5);
+  fbr_z_90_00pi->Draw()
+  fbr_z_90_05pi->Draw("same")
+  fbr_z_90_15pi->Draw("same")
+  fbr_z_90_10pi->Draw("same")
+
+  //
+  TF2 fbz_xy_0z("fbz_xy_0z","AliTPCExB::GetBz(sqrt(x^2+y^2),atan2(y,x),0)",-250,250,-250,250);
+  fbz_xy_0z.SetNpy(100);
+  fbz_xy_0z.SetNpx(100);
+  fbz_xy_0z->Draw("colz");
+  //
+  TF2 fbz_xy_250z("fbz_xy_250z","AliTPCExB::GetBz(sqrt(x^2+y^2),atan2(y,x),250)",-250,250,-250,250);
+  fbz_xy_250z.SetNpy(100);
+  fbz_xy_250z.SetNpx(100)
+  fbz_xy_250z->Draw("colz");
+  //
+   TF2 fbz_xy_m250z("fbz_xy_m250z","AliTPCExB::GetBz(sqrt(x^2+y^2),atan2(y,x),-250)",-250,250,-250,250);
+  fbz_xy_m250z.SetNpy(100);
+  fbz_xy_m250z.SetNpx(100)
+  fbz_xy_m250z->Draw("colz");
+  //
+
+
+
+
+
+*/
+
+

diff --git a/TPC/AliTPCExB.h b/TPC/AliTPCExB.h
index 7e03b417fa2f4ada050cd98c73e4cd23d3cfaf56..2a13d0b7a5fb43f10e776d09237669c77ae17b80 100644 (file)
--- a/TPC/AliTPCExB.h
+++ b/TPC/AliTPCExB.h
@@ -1,10 +1,13 @@
 #ifndef ALITPC_EXB
 #define ALITPC_EXB
 
+class AliMagF;
 #include "TObject.h"
+#include "TVectorD.h"
 
 class AliTPCExB:public TObject {
 public:
+  AliTPCExB();
   virtual ~AliTPCExB() {};
   virtual void Correct(const Double_t *position,Double_t *corrected)=0;
   virtual void CorrectInverse(const Double_t *position,Double_t *corrected) {
@@ -13,7 +16,7 @@ public:
       corrected[i]=position[i]-(corrected[i]-position[i]);
   }
   //
-  // Test and visulaization
+  // Test and visualization
   //
   void TestExB(const char* fileName);
   static Double_t GetDr(Double_t r, Double_t phi, Double_t z);
@@ -22,9 +25,40 @@ public:
   static Double_t GetDz(Double_t r, Double_t phi, Double_t z);
   static AliTPCExB*  Instance(){return fgInstance;}
   static void SetInstance(AliTPCExB*param){fgInstance = param;}
- protected:
-  static AliTPCExB*   fgInstance; //! Instance of this class (singleton implementation)
-  ClassDef(AliTPCExB,0)
+  //
+  // Mag field scans
+  //
+  static  void RegisterField(Int_t index, AliMagF * magf);
+  static  Double_t GetBx(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  Double_t GetBz(Double_t r, Double_t phi, Double_t z,Int_t index=0);
+  static  Double_t GetBr(Double_t r, Double_t phi, Double_t z,Int_t index=0);
+  static  Double_t GetBrfi(Double_t r, Double_t phi, Double_t z,Int_t index=0);
+  //
+  static  Double_t GetBxI(Double_t r, Double_t phi, Double_t z,Int_t index=0);
+  static  Double_t GetByI(Double_t r, Double_t phi, Double_t z,Int_t index=0);
+  static  Double_t GetBzI(Double_t r, Double_t phi, Double_t z,Int_t index=0);
+  static  Double_t GetBrI(Double_t r, Double_t phi, Double_t z,Int_t index=0);
+  static  Double_t GetBrfiI(Double_t r, Double_t phi, Double_t z,Int_t index=0);
+  //
+  //
+  Double_t Eval(Int_t type, Double_t r, Double_t phi, Double_t z);
+  Double_t SEval(Int_t type, Double_t r, Double_t phi, Double_t z){return Instance()->Eval(type,r,phi,z);}
+  static Double_t EvalMat(TVectorD &vec, Double_t r, Double_t phi, Double_t z);     // evalute parameterization
+ 
+ public:
+  TVectorD *          fMatBrBz;       //param matrix Br/Bz
+  TVectorD *          fMatBrfiBz;     //param matrix Br/Bz
+  TVectorD *          fMatBrBzI0;     //param matrix Br/Bz integral  z>0 
+  TVectorD *          fMatBrBzI1;     //param matrix Br/Bz integral  z<0 
+  TVectorD *          fMatBrfiBzI0;   //param matrix Br/Bz integral  z>0 
+  TVectorD *          fMatBrfiBzI1;   //param matrix Br/Bz integral  z<0
+  
+ private:
+  static AliTPCExB*   fgInstance;  //! Instance of this class (singleton implementation)
+  static TObjArray    fgArray;     //! array of magnetic fields
+  //
+  ClassDef(AliTPCExB,1)
 };
 
 #endif