#include "AliTPCExB.h"
#include "TMath.h"
-#include "TTreeStream.h"
+//#include "TTreeStream.h"
+#include "AliMagF.h"
+#include "TLinearFitter.h"
+#include "AliTPCcalibDB.h"
//
// See example macro $ALICE_ROOT/TPC/macros/AliTPCExBdraw.C
//
.L $ALICE_ROOT/TPC/macros/AliTPCExBdraw.C
- Draw0(0)
+ Draw(0)
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
+ //
+}
+
+AliTPCExB::AliTPCExB(const AliTPCExB& exb):
+ TObject(exb),
+ fMatBrBz(new TVectorD(*(exb.fMatBrBz))), //param matrix Br/Bz
+ fMatBrfiBz(new TVectorD(*(exb.fMatBrfiBz))), //param matrix Br/Bz
+ fMatBrBzI0(new TVectorD(*(exb.fMatBrBzI0))), //param matrix Br/Bz integral z>0
+ fMatBrBzI1(new TVectorD(*(exb.fMatBrBzI1))), //param matrix Br/Bz integral z<0
+ fMatBrfiBzI0(new TVectorD(*(exb.fMatBrfiBzI0))), //param matrix Br/Bz integral z>0
+ fMatBrfiBzI1(new TVectorD(*(exb.fMatBrfiBzI1))) //param matrix Br/Bz integral z<0
+{
+ //
+ // copy constructor
+ //
+}
+
+AliTPCExB& AliTPCExB::operator=(const AliTPCExB &/*exb*/)
+{
+ //
+ // Dummy assignment
+ //
+ return *this;
+}
+
+
+
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;
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";
}
}
-Double_t AliTPCExB::GetDr(Double_t r, Double_t phi, Double_t z){
+Double_t AliTPCExB::GetDr(Double_t r, Double_t phi, Double_t z, Double_t bz){
//
// Static function
// Posibble to us it for visualization
//
//
- if (!fgInstance) return 0;
- Double_t pos0[3] = {r*TMath::Sin(phi), r*TMath::Cos(phi),z};
+ AliTPCExB *exb = Instance();
+ if (!exb) exb = AliTPCcalibDB::GetExB(bz,kFALSE);
+ if (!exb) return 0;
+ Double_t pos0[3] = {r*TMath::Cos(phi), r*TMath::Sin(phi),z};
Double_t pos1[3];
- fgInstance->Correct(pos0,pos1);
+ exb->Correct(pos0,pos1);
Double_t dx=pos1[0]-pos0[0];
Double_t dy=pos1[1]-pos0[1];
// Double_t dz=pos1[2]-pos0[2];
- return TMath::Sqrt(dx*dx+dy*dy);
+ // return TMath::Sqrt(dx*dx+dy*dy);
+ Float_t dr = (dx*pos0[0]+dy*pos0[1])/r;
+ return dr;
}
-Double_t AliTPCExB::GetDrphi(Double_t r, Double_t phi, Double_t z){
+Double_t AliTPCExB::GetDrphi(Double_t r, Double_t phi, Double_t z, Double_t bz){
//
//
//
- if (!fgInstance) return 0;
- Double_t pos0[3] = {r*TMath::Sin(phi), r*TMath::Cos(phi),z};
+ AliTPCExB *exb = Instance();
+ if (!exb) exb = AliTPCcalibDB::GetExB(bz,kFALSE);
+ if (!exb) return 0;
+ Double_t pos0[3] = {r*TMath::Cos(phi), r*TMath::Sin(phi),z};
Double_t pos1[3];
- fgInstance->Correct(pos0,pos1);
+ exb->Correct(pos0,pos1);
Double_t dphi=TMath::ATan2(pos1[1],pos1[0])-TMath::ATan2(pos0[1],pos0[0]);
+ if (dphi>TMath::Pi()) dphi-=TMath::TwoPi();
+ if (dphi<-TMath::Pi()) dphi+=TMath::TwoPi();
return r*dphi;
}
-Double_t AliTPCExB::GetDphi(Double_t r, Double_t phi, Double_t z){
- //
+Double_t AliTPCExB::GetDphi(Double_t r, Double_t phi, Double_t z, Double_t bz){
//
//
- if (!fgInstance) return 0;
- Double_t pos0[3] = {r*TMath::Sin(phi), r*TMath::Cos(phi),z};
+ //
+ AliTPCExB *exb = Instance();
+ if (!exb) exb = AliTPCcalibDB::GetExB(bz,kFALSE);
+ if (!exb) return 0;
+ Double_t pos0[3] = {r*TMath::Cos(phi), r*TMath::Sin(phi),z};
Double_t pos1[3];
- fgInstance->Correct(pos0,pos1);
+ exb->Correct(pos0,pos1);
Double_t dphi=TMath::ATan2(pos1[1],pos1[0])-TMath::ATan2(pos0[1],pos0[0]);
return dphi;
}
-Double_t AliTPCExB::GetDz(Double_t r, Double_t phi, Double_t z){
+Double_t AliTPCExB::GetDz(Double_t r, Double_t phi, Double_t z, Double_t bz){
//
//
//
- if (!fgInstance) return 0;
- Double_t pos0[3] = {r*TMath::Sin(phi), r*TMath::Cos(phi),z};
+ AliTPCExB *exb = Instance();
+ if (!exb) exb = AliTPCcalibDB::GetExB(bz,kFALSE);
+ if (!exb) return 0;
+ Double_t pos0[3] = {r*TMath::Cos(phi), r*TMath::Sin(phi),z};
Double_t pos1[3];
- fgInstance->Correct(pos0,pos1);
+ exb->Correct(pos0,pos1);
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;
+ Double_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
+ // xyz[1]+=30;
+ Double_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;
+ Double_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
+ // xyz[1]+=30;
+ Double_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;
+ Double_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
+ // xyz[1]+=30;
+ Double_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;
+ Double_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
+ //xyz[1]+=30;
+ Double_t bxyz[3];
+ mag->Field(xyz,bxyz);
+ if (r==0) return 0;
+ Double_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;
+ Double_t xyz[3]={r*TMath::Cos(phi),r*TMath::Sin(phi),z};
+ //xyz[1]+=30;
+ Double_t bxyz[3];
+ mag->Field(xyz,bxyz);
+ if (r==0) return 0;
+ Double_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(const 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");
+ //
+
+
+
+
+
+*/
+
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff