X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPCExB.cxx;h=1ec7ec29a27fe3d9a6159d30de78553829dc2b3d;hb=d42184d72fb96685f5ba8294a08d3d82797c15cd;hp=cff72361408ef0a9f814e3cab5f09d42c302978c;hpb=faf9323738d83434592a89e8b4439f5b9130db9c;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCExB.cxx b/TPC/AliTPCExB.cxx index cff72361408..1ec7ec29a27 100644 --- a/TPC/AliTPCExB.cxx +++ b/TPC/AliTPCExB.cxx @@ -1,4 +1,552 @@ #include "AliTPCExB.h" +#include "TMath.h" +#include "TTreeStream.h" +#include "AliMagF.h" +#include "TLinearFitter.h" +#include "AliTPCcalibDB.h" + + +// +// Abstract class for ExB effect parameterization +// +// +// +// The ExB correction map is stored in the calib DB +// The lookup can be dumped to the tree: +/* + + // + char *storage = "local://OCDBres" + Int_t RunNumber=0; + AliCDBManager::Instance()->SetDefaultStorage(storage); + AliCDBManager::Instance()->SetRun(RunNumber) + AliTPCExBFirst * exb = AliTPCcalibDB::Instance()->GetExB(); + // + // See example macro $ALICE_ROOT/TPC/macros/AliTPCExBdraw.C + // + .L $ALICE_ROOT/TPC/macros/AliTPCExBdraw.C + 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) { + // + // 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]+=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 rd=TMath::Sqrt(d[0]*d[0]+d[1]*d[1]); + Double_t dr=r-rd; + 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 drphi=r*dphi; + Double_t dx=x[0]-d[0]; + Double_t dy=x[1]-d[1]; + Double_t dz=x[2]-d[2]; + // + 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="<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); + 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 bz){ + // + // + // + 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]; + 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 bz){ + // + // + // + 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]; + 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 bz){ + // + // + // + 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]; + 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(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"); + // + + + + + +*/ + +