X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliMagWrapCheb.cxx;h=b194ce2ed7ac44fdbad299a661b7901decaa33a7;hb=1e122ba8ea49d0934052b7fba98235f2f440c60d;hp=dbc939d3de7ef24b71e356adf0acadf8556ce3c6;hpb=82d1e5564c07302098d0119ca821b93178aab3e0;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliMagWrapCheb.cxx b/STEER/AliMagWrapCheb.cxx index dbc939d3de7..b194ce2ed7a 100644 --- a/STEER/AliMagWrapCheb.cxx +++ b/STEER/AliMagWrapCheb.cxx @@ -22,52 +22,17 @@ ClassImp(AliMagWrapCheb) //__________________________________________________________________________________________ AliMagWrapCheb::AliMagWrapCheb() : - fNParamsSol(0), - fNSegZSol(0), - fNParamsTPCInt(0), - fNSegZTPCInt(0), - fNParamsDip(0), +fNParamsSol(0),fNZSegSol(0),fNPSegSol(0),fNRSegSol(0), + fSegZSol(0),fSegPSol(0),fSegRSol(0), + fBegSegPSol(0),fNSegPSol(0),fBegSegRSol(0),fNSegRSol(0),fSegIDSol(0),fMinZSol(1.e6),fMaxZSol(-1.e6),fParamsSol(0),fMaxRSol(0), // - fNZSegDip(0), - fNYSegDip(0), - fNXSegDip(0), + fNParamsTPC(0),fNZSegTPC(0),fNPSegTPC(0),fNRSegTPC(0), + fSegZTPC(0),fSegPTPC(0),fSegRTPC(0), + fBegSegPTPC(0),fNSegPTPC(0),fBegSegRTPC(0),fNSegRTPC(0),fSegIDTPC(0),fMinZTPC(1.e6),fMaxZTPC(-1.e6),fParamsTPC(0),fMaxRTPC(0), // - fSegZSol(0), - fSegRSol(0), -// - fSegZTPCInt(0), - fSegRTPCInt(0), -// - fSegZDip(0), - fSegYDip(0), - fSegXDip(0), -// - fNSegRSol(0), - fSegZIdSol(0), -// - fNSegRTPCInt(0), - fSegZIdTPCInt(0), -// - fBegSegYDip(0), - fNSegYDip(0), - fBegSegXDip(0), - fNSegXDip(0), - fSegIDDip(0), -// - fMinZSol(1e6), - fMaxZSol(-1e6), - fMaxRSol(-1e6), -// - fMinZDip(1e6), - fMaxZDip(-1e6), -// - fMinZTPCInt(1e6), - fMaxZTPCInt(-1e6), - fMaxRTPCInt(-1e6), -// - fParamsSol(0), - fParamsDip(0), - fParamsTPCInt(0) + fNParamsDip(0),fNZSegDip(0),fNYSegDip(0),fNXSegDip(0), + fSegZDip(0),fSegYDip(0),fSegXDip(0), + fBegSegYDip(0),fNSegYDip(0),fBegSegXDip(0),fNSegXDip(0),fSegIDDip(0),fMinZDip(1.e6),fMaxZDip(-1.e6),fParamsDip(0) // { // default constructor @@ -76,52 +41,17 @@ AliMagWrapCheb::AliMagWrapCheb() : //__________________________________________________________________________________________ AliMagWrapCheb::AliMagWrapCheb(const AliMagWrapCheb& src) : TNamed(src), - fNParamsSol(0), - fNSegZSol(0), - fNParamsTPCInt(0), - fNSegZTPCInt(0), - fNParamsDip(0), -// - fNZSegDip(0), - fNYSegDip(0), - fNXSegDip(0), -// - fSegZSol(0), - fSegRSol(0), -// - fSegZTPCInt(0), - fSegRTPCInt(0), + fNParamsSol(0),fNZSegSol(0),fNPSegSol(0),fNRSegSol(0), + fSegZSol(0),fSegPSol(0),fSegRSol(0), + fBegSegPSol(0),fNSegPSol(0),fBegSegRSol(0),fNSegRSol(0),fSegIDSol(0),fMinZSol(1.e6),fMaxZSol(-1.e6),fParamsSol(0),fMaxRSol(0), // - fSegZDip(0), - fSegYDip(0), - fSegXDip(0), + fNParamsTPC(0),fNZSegTPC(0),fNPSegTPC(0),fNRSegTPC(0), + fSegZTPC(0),fSegPTPC(0),fSegRTPC(0), + fBegSegPTPC(0),fNSegPTPC(0),fBegSegRTPC(0),fNSegRTPC(0),fSegIDTPC(0),fMinZTPC(1.e6),fMaxZTPC(-1.e6),fParamsTPC(0),fMaxRTPC(0), // - fNSegRSol(0), - fSegZIdSol(0), -// - fNSegRTPCInt(0), - fSegZIdTPCInt(0), -// - fBegSegYDip(0), - fNSegYDip(0), - fBegSegXDip(0), - fNSegXDip(0), - fSegIDDip(0), -// - fMinZSol(1e6), - fMaxZSol(-1e6), - fMaxRSol(-1e6), -// - fMinZDip(1e6), - fMaxZDip(-1e6), -// - fMinZTPCInt(1e6), - fMaxZTPCInt(-1e6), - fMaxRTPCInt(-1e6), -// - fParamsSol(0), - fParamsDip(0), - fParamsTPCInt(0) + fNParamsDip(0),fNZSegDip(0),fNYSegDip(0),fNXSegDip(0), + fSegZDip(0),fSegYDip(0),fSegXDip(0), + fBegSegYDip(0),fNSegYDip(0),fBegSegXDip(0),fNSegXDip(0),fSegIDDip(0),fMinZDip(1.e6),fMaxZDip(-1.e6),fParamsDip(0) { // copy constructor CopyFrom(src); @@ -134,36 +64,53 @@ void AliMagWrapCheb::CopyFrom(const AliMagWrapCheb& src) Clear(); SetName(src.GetName()); SetTitle(src.GetTitle()); + // fNParamsSol = src.fNParamsSol; - fNSegZSol = src.fNSegZSol; - fNParamsTPCInt = src.fNParamsTPCInt; - fNSegZTPCInt = src.fNSegZTPCInt; - fNParamsDip = src.fNParamsDip; + fNZSegSol = src.fNZSegSol; + fNPSegSol = src.fNPSegSol; + fNRSegSol = src.fNRSegSol; + fMinZSol = src.fMinZSol; + fMaxZSol = src.fMaxZSol; + fMaxRSol = src.fMaxRSol; + if (src.fNParamsSol) { + memcpy(fSegZSol = new Float_t[fNZSegSol], src.fSegZSol, sizeof(Float_t)*fNZSegSol); + memcpy(fSegPSol = new Float_t[fNPSegSol], src.fSegPSol, sizeof(Float_t)*fNPSegSol); + memcpy(fSegRSol = new Float_t[fNRSegSol], src.fSegRSol, sizeof(Float_t)*fNRSegSol); + memcpy(fBegSegPSol= new Int_t[fNZSegSol], src.fBegSegPSol, sizeof(Int_t)*fNZSegSol); + memcpy(fNSegPSol = new Int_t[fNZSegSol], src.fNSegPSol, sizeof(Int_t)*fNZSegSol); + memcpy(fBegSegRSol= new Int_t[fNPSegSol], src.fBegSegRSol, sizeof(Int_t)*fNPSegSol); + memcpy(fNSegRSol = new Int_t[fNPSegSol], src.fNSegRSol, sizeof(Int_t)*fNPSegSol); + memcpy(fSegIDSol = new Int_t[fNRSegSol], src.fSegIDSol, sizeof(Int_t)*fNRSegSol); + fParamsSol = new TObjArray(fNParamsSol); + for (int i=0;iAddAtAndExpand(new AliCheb3D(*src.GetParamSol(i)),i); + } // + fNParamsTPC = src.fNParamsTPC; + fNZSegTPC = src.fNZSegTPC; + fNPSegTPC = src.fNPSegTPC; + fNRSegTPC = src.fNRSegTPC; + fMinZTPC = src.fMinZTPC; + fMaxZTPC = src.fMaxZTPC; + fMaxRTPC = src.fMaxRTPC; + if (src.fNParamsTPC) { + memcpy(fSegZTPC = new Float_t[fNZSegTPC], src.fSegZTPC, sizeof(Float_t)*fNZSegTPC); + memcpy(fSegPTPC = new Float_t[fNPSegTPC], src.fSegPTPC, sizeof(Float_t)*fNPSegTPC); + memcpy(fSegRTPC = new Float_t[fNRSegTPC], src.fSegRTPC, sizeof(Float_t)*fNRSegTPC); + memcpy(fBegSegPTPC= new Int_t[fNZSegTPC], src.fBegSegPTPC, sizeof(Int_t)*fNZSegTPC); + memcpy(fNSegPTPC = new Int_t[fNZSegTPC], src.fNSegPTPC, sizeof(Int_t)*fNZSegTPC); + memcpy(fBegSegRTPC= new Int_t[fNPSegTPC], src.fBegSegRTPC, sizeof(Int_t)*fNPSegTPC); + memcpy(fNSegRTPC = new Int_t[fNPSegTPC], src.fNSegRTPC, sizeof(Int_t)*fNPSegTPC); + memcpy(fSegIDTPC = new Int_t[fNRSegTPC], src.fSegIDTPC, sizeof(Int_t)*fNRSegTPC); + fParamsTPC = new TObjArray(fNParamsTPC); + for (int i=0;iAddAtAndExpand(new AliCheb3D(*src.GetParamTPCInt(i)),i); + } + // + fNParamsDip = src.fNParamsDip; fNZSegDip = src.fNZSegDip; fNYSegDip = src.fNYSegDip; fNXSegDip = src.fNXSegDip; - // - fMinZSol = src.fMinZSol; - fMaxZSol = src.fMaxZSol; - fMaxRSol = src.fMaxRSol; - // fMinZDip = src.fMinZDip; fMaxZDip = src.fMaxZDip; - // - fMinZTPCInt = src.fMinZTPCInt; - fMaxZTPCInt = src.fMaxZTPCInt; - fMaxRTPCInt = src.fMaxRTPCInt; - // - if (src.fNParamsSol) { - memcpy(fSegZSol = new Float_t[fNSegZSol], src.fSegZSol, sizeof(Float_t)*fNSegZSol); - memcpy(fSegRSol = new Float_t[fNParamsSol], src.fSegRSol, sizeof(Float_t)*fNParamsSol); - memcpy(fNSegRSol = new Int_t[fNSegZSol], src.fNSegRSol, sizeof(Int_t)*fNSegZSol); - memcpy(fSegZIdSol= new Int_t[fNSegZSol], src.fSegZIdSol, sizeof(Int_t)*fNSegZSol); - fParamsSol = new TObjArray(fNParamsSol); - for (int i=0;iAddAtAndExpand(new AliCheb3D(*src.GetParamSol(i)),i); - } - // if (src.fNParamsDip) { memcpy(fSegZDip = new Float_t[fNZSegDip], src.fSegZDip, sizeof(Float_t)*fNZSegDip); memcpy(fSegYDip = new Float_t[fNYSegDip], src.fSegYDip, sizeof(Float_t)*fNYSegDip); @@ -177,15 +124,6 @@ void AliMagWrapCheb::CopyFrom(const AliMagWrapCheb& src) for (int i=0;iAddAtAndExpand(new AliCheb3D(*src.GetParamDip(i)),i); } // - if (src.fNParamsTPCInt) { - memcpy(fSegZTPCInt = new Float_t[fNSegZTPCInt], src.fSegZTPCInt, sizeof(Float_t)*fNSegZTPCInt); - memcpy(fSegRTPCInt = new Float_t[fNParamsTPCInt], src.fSegRTPCInt, sizeof(Float_t)*fNParamsTPCInt); - memcpy(fNSegRTPCInt = new Int_t[fNSegZTPCInt], src.fNSegRTPCInt, sizeof(Int_t)*fNSegZTPCInt); - memcpy(fSegZIdTPCInt= new Int_t[fNSegZTPCInt], src.fSegZIdTPCInt, sizeof(Int_t)*fNSegZTPCInt); - fParamsTPCInt = new TObjArray(fNParamsTPCInt); - for (int i=0;iAddAtAndExpand(new AliCheb3D(*src.GetParamTPCInt(i)),i); - } - // } //__________________________________________________________________________________________ @@ -205,36 +143,51 @@ void AliMagWrapCheb::Clear(const Option_t *) { // clear all dynamic parts if (fNParamsSol) { - delete fParamsSol; - delete[] fSegZSol; - delete[] fSegRSol; - delete[] fNSegRSol; - delete[] fSegZIdSol; - } + delete fParamsSol; fParamsSol = 0; + delete[] fSegZSol; fSegZSol = 0; + delete[] fSegPSol; fSegPSol = 0; + delete[] fSegRSol; fSegRSol = 0; + delete[] fBegSegPSol; fBegSegPSol = 0; + delete[] fNSegPSol; fNSegPSol = 0; + delete[] fBegSegRSol; fBegSegRSol = 0; + delete[] fNSegRSol; fNSegRSol = 0; + delete[] fSegIDSol; fSegIDSol = 0; + } + fNParamsSol = fNZSegSol = fNPSegSol = fNRSegSol = 0; + fMinZSol = 1e6; + fMaxZSol = -1e6; + fMaxRSol = 0; + // + if (fNParamsTPC) { + delete fParamsTPC; fParamsTPC = 0; + delete[] fSegZTPC; fSegZTPC = 0; + delete[] fSegPTPC; fSegPTPC = 0; + delete[] fSegRTPC; fSegRTPC = 0; + delete[] fBegSegPTPC; fBegSegPTPC = 0; + delete[] fNSegPTPC; fNSegPTPC = 0; + delete[] fBegSegRTPC; fBegSegRTPC = 0; + delete[] fNSegRTPC; fNSegRTPC = 0; + delete[] fSegIDTPC; fSegIDTPC = 0; + } + fNParamsTPC = fNZSegTPC = fNPSegTPC = fNRSegTPC = 0; + fMinZTPC = 1e6; + fMaxZTPC = -1e6; + fMaxRTPC = 0; // - if (fNParamsTPCInt) { - delete fParamsTPCInt; - delete[] fSegZTPCInt; - delete[] fSegRTPCInt; - delete[] fNSegRTPCInt; - delete[] fSegZIdTPCInt; - } - // if (fNParamsDip) { - delete fParamsDip; - delete[] fSegZDip; - delete[] fSegYDip; - delete[] fSegXDip; - delete[] fBegSegYDip; - delete[] fNSegYDip; - delete[] fBegSegXDip; - delete[] fNSegXDip; - delete[] fSegIDDip; - } - fNParamsSol = fNParamsTPCInt = fNParamsDip = fNZSegDip = fNYSegDip = fNXSegDip = 0; - fNSegZSol = fNSegZTPCInt = 0; - fMinZSol = fMinZDip = fMinZTPCInt = 1e6; - fMaxZSol = fMaxZDip = fMaxZTPCInt = fMaxRSol = fMaxRTPCInt = -1e6; + delete fParamsDip; fParamsDip = 0; + delete[] fSegZDip; fSegZDip = 0; + delete[] fSegYDip; fSegYDip = 0; + delete[] fSegXDip; fSegXDip = 0; + delete[] fBegSegYDip; fBegSegYDip = 0; + delete[] fNSegYDip; fNSegYDip = 0; + delete[] fBegSegXDip; fBegSegXDip = 0; + delete[] fNSegXDip; fNSegXDip = 0; + delete[] fSegIDDip; fSegIDDip = 0; + } + fNParamsDip = fNZSegDip = fNYSegDip = fNXSegDip = 0; + fMinZDip = 1e6; + fMaxZDip = -1e6; // } @@ -243,33 +196,27 @@ void AliMagWrapCheb::Field(const Double_t *xyz, Double_t *b) const { // compute field in cartesian coordinates. If point is outside of the parameterized region // get it at closest valid point - static Double_t rphiz[3]; + Double_t rphiz[3]; // #ifndef _BRING_TO_BOUNDARY_ // exact matching to fitted volume is requested - if ( !(xyz[2]>=GetMinZSol()&&xyz[2]<=GetMaxZSol()) && - !(xyz[2]>=GetMinZDip()&&xyz[2]<=GetMaxZDip()) ) {for (int i=3;i--;) b[i]=0; return;} + b[0] = b[1] = b[2] = 0; #endif // - if (xyz[2]IsInside(xyz)) {par->Eval(xyz,b); return;} - for (int i=3;i--;) b[i]=0; return; -#else - GetParamDip(FindDipSegment(xyz))->Eval(xyz,b); return; -#endif + if (xyz[2]>fMinZSol) { + CartToCyl(xyz,rphiz); + FieldCylSol(rphiz,b); + // convert field to cartesian system + CylToCartCylB(rphiz, b,b); + return; } // - // Sol region: convert coordinates to cyl system - CartToCyl(xyz,rphiz); + int iddip = FindDipSegment(xyz); + if (iddip<0) return; + AliCheb3D* par = GetParamDip(iddip); #ifndef _BRING_TO_BOUNDARY_ - if (rphiz[0]>GetMaxRSol()) {for (int i=3;i--;) b[i]=0; return;} + if (!par->IsInside(xyz)) return; #endif - // - FieldCylSol(rphiz,b); - // - // convert field to cartesian system - CylToCartCylB(rphiz, b,b); + par->Eval(xyz,b); // } @@ -278,29 +225,20 @@ Double_t AliMagWrapCheb::GetBz(const Double_t *xyz) const { // compute Bz for the point in cartesian coordinates. If point is outside of the parameterized region // get it at closest valid point - static Double_t rphiz[3]; + Double_t rphiz[3]; // -#ifndef _BRING_TO_BOUNDARY_ // exact matching to fitted volume is requested - if ( !(xyz[2]>=GetMinZSol()&&xyz[2]<=GetMaxZSol()) && - !(xyz[2]>=GetMinZDip()&&xyz[2]<=GetMaxZDip()) ) return 0.; -#endif - // - if (xyz[2]IsInside(xyz)) return par->Eval(xyz,2); - else return 0.; -#else - return GetParamDip(FindDipSegment(xyz))->Eval(xyz,2); -#endif + if (xyz[2]>fMinZSol) { + CartToCyl(xyz,rphiz); + return FieldCylSolBz(rphiz); } - // Sol region: convert coordinates to cyl system - CartToCyl(xyz,rphiz); + // + int iddip = FindDipSegment(xyz); + if (iddip<0) return 0.; + AliCheb3D* par = GetParamDip(iddip); #ifndef _BRING_TO_BOUNDARY_ - if (rphiz[0]>GetMaxRSol()) return 0.; + if (!par->IsInside(xyz)) return 0.; #endif - // - return FieldCylSolBz(rphiz); + return par->Eval(xyz,2); } @@ -317,39 +255,16 @@ void AliMagWrapCheb::Print(Option_t *) const GetParamSol(i)->Print(); } } - /* - for (int iz=0;izGetBoundMin(2),param->GetBoundMax(2)); - for (int ir=0;irGetBoundMin(0), - param->GetBoundMax(0),param->GetPrecision(),fSegZIdSol[iz]+ir); - } - } - */ // - printf("Segmentation for TPC field integral (%+.2fPrint(); } } // - /* - for (int iz=0;izGetBoundMin(2),param->GetBoundMax(2)); - for (int ir=0;irGetBoundMin(0), - param->GetBoundMax(0),param->GetPrecision(),fSegZIdTPCInt[iz]+ir); - } - } - */ - // printf("Segmentation for Dipole (%+.2ffSegZSol[SolZId] && SolZIdfSegRSol[SolRId] && SolRIdEval(rphiz,b); + int id = FindSolSegment(rphiz); + if (id<0) return; + AliCheb3D* par = GetParamSol(id); +#ifndef _BRING_TO_BOUNDARY_ // exact matching to fitted volume is requested + if (!par->IsInside(rphiz)) return; +#endif + par->Eval(rphiz,b); + return; // } @@ -421,12 +379,13 @@ Double_t AliMagWrapCheb::FieldCylSolBz(const Double_t *rphiz) const { // compute Solenoid field in Cylindircal coordinates // note: if the point is outside the volume get the field in closest parameterized point - int SolZId = 0; - while (rphiz[2]>fSegZSol[SolZId] && SolZIdfSegRSol[SolRId] && SolRIdEval(rphiz,2); + int id = FindSolSegment(rphiz); + if (id<0) return 0.; + AliCheb3D* par = GetParamSol(id); +#ifndef _BRING_TO_BOUNDARY_ + return par->IsInside(rphiz) ? par->Eval(rphiz,2) : 0; +#endif + return par->Eval(rphiz,2); // } @@ -435,12 +394,15 @@ void AliMagWrapCheb::GetTPCIntCyl(const Double_t *rphiz, Double_t *b) const { // compute field integral in TPC region in Cylindircal coordinates // note: the check for the point being inside the parameterized region is done outside - int tpcIntZId = 0; - while (rphiz[2]>fSegZTPCInt[tpcIntZId] && tpcIntZIdfSegRTPCInt[tpcIntRId] && tpcIntRIdEval(rphiz,b); + int id = FindTPCSegment(rphiz); + if (id<0) { + b[0] = b[1] = b[2] = 0; + return; + } + AliCheb3D* par = GetParamTPCInt(id); + if (par->IsInside(rphiz)) {par->Eval(rphiz,b); return;} + b[0] = b[1] = b[2] = 0; + return; // } @@ -539,11 +501,49 @@ void AliMagWrapCheb::LoadData(const char* inpfile) // --------------------------------------------------------------------------- fclose(stream); BuildTableSol(); - BuildTableTPCInt(); BuildTableDip(); + BuildTableTPCInt(); + // printf("Loaded magnetic field \"%s\" from %s\n",GetName(),strf.Data()); // } + +//__________________________________________________________________________________________ +void AliMagWrapCheb::BuildTableSol() +{ + BuildTable(fNParamsSol,fParamsSol, + fNZSegSol,fNPSegSol,fNRSegSol, + fMinZSol,fMaxZSol, + &fSegZSol,&fSegPSol,&fSegRSol, + &fBegSegPSol,&fNSegPSol, + &fBegSegRSol,&fNSegRSol, + &fSegIDSol); +} + +//__________________________________________________________________________________________ +void AliMagWrapCheb::BuildTableDip() +{ + BuildTable(fNParamsDip,fParamsDip, + fNZSegDip,fNYSegDip,fNXSegDip, + fMinZDip,fMaxZDip, + &fSegZDip,&fSegYDip,&fSegXDip, + &fBegSegYDip,&fNSegYDip, + &fBegSegXDip,&fNSegXDip, + &fSegIDDip); +} + +//__________________________________________________________________________________________ +void AliMagWrapCheb::BuildTableTPCInt() +{ + BuildTable(fNParamsTPC,fParamsTPC, + fNZSegTPC,fNPSegTPC,fNRSegTPC, + fMinZTPC,fMaxZTPC, + &fSegZTPC,&fSegPTPC,&fSegRTPC, + &fBegSegPTPC,&fNSegPTPC, + &fBegSegRTPC,&fNSegRTPC, + &fSegIDTPC); +} + #endif //_______________________________________________ @@ -551,53 +551,17 @@ void AliMagWrapCheb::LoadData(const char* inpfile) //__________________________________________________________________________________________ AliMagWrapCheb::AliMagWrapCheb(const char* inputFile) : - fNParamsSol(0), - fNSegZSol(0), - fNParamsTPCInt(0), - fNSegZTPCInt(0), - fNParamsDip(0), -// - fNZSegDip(0), - fNYSegDip(0), - fNXSegDip(0), -// - fSegZSol(0), - fSegRSol(0), -// - fSegZTPCInt(0), - fSegRTPCInt(0), -// - fSegZDip(0), - fSegYDip(0), - fSegXDip(0), -// - fNSegRSol(0), - fSegZIdSol(0), + fNParamsSol(0),fNZSegSol(0),fNPSegSol(0),fNRSegSol(0), + fSegZSol(0),fSegPSol(0),fSegRSol(0), + fBegSegPSol(0),fNSegPSol(0),fBegSegRSol(0),fNSegRSol(0),fSegIDSol(0),fMinZSol(1.e6),fMaxZSol(-1.e6),fParamsSol(0),fMaxRSol(0), // - fNSegRTPCInt(0), - fSegZIdTPCInt(0), -// - fBegSegYDip(0), - fNSegYDip(0), - fBegSegXDip(0), - fNSegXDip(0), - fSegIDDip(0), -// - fMinZSol(0), - fMaxZSol(0), - fMaxRSol(0), -// - fMinZDip(0), - fMaxZDip(0), -// - fMinZTPCInt(0), - fMaxZTPCInt(0), - fMaxRTPCInt(0), -// - fParamsSol(0), - fParamsDip(0), - fParamsTPCInt(0) + fNParamsTPC(0),fNZSegTPC(0),fNPSegTPC(0),fNRSegTPC(0), + fSegZTPC(0),fSegPTPC(0),fSegRTPC(0), + fBegSegPTPC(0),fNSegPTPC(0),fBegSegRTPC(0),fNSegRTPC(0),fSegIDTPC(0),fMinZTPC(1.e6),fMaxZTPC(-1.e6),fParamsTPC(0),fMaxRTPC(0), // + fNParamsDip(0),fNZSegDip(0),fNYSegDip(0),fNXSegDip(0), + fSegZDip(0),fSegYDip(0),fSegXDip(0), + fBegSegYDip(0),fNSegYDip(0),fBegSegXDip(0),fNSegXDip(0),fSegIDDip(0),fMinZDip(1.e6),fMaxZDip(-1.e6),fParamsDip(0) // { // construct from coeffs from the text file @@ -613,6 +577,7 @@ void AliMagWrapCheb::AddParamSol(const AliCheb3D* param) if (!fParamsSol) fParamsSol = new TObjArray(); fParamsSol->Add( (AliCheb3D*)param ); fNParamsSol++; + if (fMaxRSolGetBoundMax(0)) fMaxRSol = param->GetBoundMax(0); // } @@ -622,9 +587,10 @@ void AliMagWrapCheb::AddParamTPCInt(const AliCheb3D* param) // adds new parameterization piece for TPCInt // NOTE: pieces must be added strictly in increasing R then increasing Z order // - if (!fParamsTPCInt) fParamsTPCInt = new TObjArray(); - fParamsTPCInt->Add( (AliCheb3D*)param); - fNParamsTPCInt++; + if (!fParamsTPC) fParamsTPC = new TObjArray(); + fParamsTPC->Add( (AliCheb3D*)param); + fNParamsTPC++; + if (fMaxRTPCGetBoundMax(0)) fMaxRTPC = param->GetBoundMax(0); // } @@ -643,26 +609,123 @@ void AliMagWrapCheb::AddParamDip(const AliCheb3D* param) void AliMagWrapCheb::ResetTPCInt() { // clean TPC field integral (used for update) - if (!fNParamsTPCInt) return; - delete fParamsTPCInt; - delete[] fSegZTPCInt; - delete[] fSegRTPCInt; - delete[] fNSegRTPCInt; - delete[] fSegZIdTPCInt; - // - fNParamsTPCInt = 0; - fNSegZTPCInt = 0; - fSegZTPCInt = 0; - fSegRTPCInt = 0; - fNSegRTPCInt = 0; - fSegZIdTPCInt = 0; - fMinZTPCInt = 0; - fMaxZTPCInt = 0; - fMaxRTPCInt = 0; - fParamsTPCInt = 0; + if (fNParamsTPC) { + delete fParamsTPC; fParamsTPC = 0; + delete[] fSegZTPC; fSegZTPC = 0; + delete[] fSegPTPC; fSegPTPC = 0; + delete[] fSegRTPC; fSegRTPC = 0; + delete[] fBegSegPTPC; fBegSegPTPC = 0; + delete[] fNSegPTPC; fNSegPTPC = 0; + delete[] fBegSegRTPC; fBegSegRTPC = 0; + delete[] fNSegRTPC; fNSegRTPC = 0; + delete[] fSegIDTPC; fSegIDTPC = 0; + } + fNParamsTPC = fNZSegTPC = fNPSegTPC = fNRSegTPC = 0; + fMinZTPC = 1e6; + fMaxZTPC = -1e6; + fMaxRTPC = 0; + // +} + + +//__________________________________________________ +void AliMagWrapCheb::BuildTable(Int_t npar,TObjArray *parArr, Int_t &nZSeg, Int_t &nYSeg, Int_t &nXSeg, + Float_t &minZ,Float_t &maxZ, + Float_t **segZ,Float_t **segY,Float_t **segX, + Int_t **begSegY,Int_t **nSegY, + Int_t **begSegX,Int_t **nSegX, + Int_t **segID) +{ + // build lookup table for dipole + // + if (npar<1) return; + TArrayF segYArr,segXArr; + TArrayI begSegYDipArr,begSegXDipArr; + TArrayI nSegYDipArr,nSegXDipArr; + TArrayI segIDArr; + float *tmpSegZ,*tmpSegY,*tmpSegX; + // + // create segmentation in Z + nZSeg = SegmentDimension(&tmpSegZ, parArr, npar, 2, 1,-1, 1,-1, 1,-1) - 1; + nYSeg = 0; + nXSeg = 0; + // + // for each Z slice create segmentation in Y + begSegYDipArr.Set(nZSeg); + nSegYDipArr.Set(nZSeg); + float xyz[3]; + for (int iz=0;izAt(ipar); + if (!cheb->IsInside(xyz)) continue; + segIDArr[nXSeg+ix] = ipar; + break; + } + } + nXSeg += nx; + // + delete[] tmpSegX; + } + delete[] tmpSegY; + nYSeg += ny; + } + // + minZ = tmpSegZ[0]; + maxZ = tmpSegZ[nZSeg]; + (*segZ) = new Float_t[nZSeg]; + for (int i=nZSeg;i--;) (*segZ)[i] = tmpSegZ[i]; + delete[] tmpSegZ; + // + (*segY) = new Float_t[nYSeg]; + (*segX) = new Float_t[nXSeg]; + (*begSegY) = new Int_t[nZSeg]; + (*nSegY) = new Int_t[nZSeg]; + (*begSegX) = new Int_t[nYSeg]; + (*nSegX) = new Int_t[nYSeg]; + (*segID) = new Int_t[nXSeg]; + // + for (int i=nYSeg;i--;) (*segY)[i] = segYArr[i]; + for (int i=nXSeg;i--;) (*segX)[i] = segXArr[i]; + for (int i=nZSeg;i--;) {(*begSegY)[i] = begSegYDipArr[i]; (*nSegY)[i] = nSegYDipArr[i];} + for (int i=nYSeg;i--;) {(*begSegX)[i] = begSegXDipArr[i]; (*nSegX)[i] = nSegXDipArr[i];} + for (int i=nXSeg;i--;) {(*segID)[i] = segIDArr[i];} // } +/* //__________________________________________________ void AliMagWrapCheb::BuildTableDip() { @@ -676,7 +739,7 @@ void AliMagWrapCheb::BuildTableDip() float *tmpSegZ,*tmpSegY,*tmpSegX; // // create segmentation in Z - fNZSegDip = SegmentDipDimension(&tmpSegZ, fParamsDip, fNParamsDip, 2, 1,-1, 1,-1, 1,-1) - 1; + fNZSegDip = SegmentDimension(&tmpSegZ, fParamsDip, fNParamsDip, 2, 1,-1, 1,-1, 1,-1) - 1; fNYSegDip = 0; fNXSegDip = 0; // @@ -686,7 +749,7 @@ void AliMagWrapCheb::BuildTableDip() float xyz[3]; for (int iz=0;izGetBoundMax(2)-GetParamSol(ip-1)->GetBoundMax(2))>kSafety) { // new Z slice - tmpbufF[fNSegZSol] = GetParamSol(ip)->GetBoundMax(2); // - tmpbufI[fNSegZSol] = 0; - tmpbufI1[fNSegZSol++] = ip; - } - fSegRSol[ip] = GetParamSol(ip)->GetBoundMax(0); // upper R - tmpbufI[fNSegZSol-1]++; - if (fMaxRSolGetBoundMin(2); - fMaxZSol = GetParamSol(fNParamsSol-1)->GetBoundMax(2); - // - delete[] tmpbufF; - delete[] tmpbufI; - delete[] tmpbufI1; - // - // -} - -//__________________________________________________________________________________________ -void AliMagWrapCheb::BuildTableTPCInt() -{ - // build the indexes for each parameterization of TPC field integral - // - const float kSafety=0.001; - // - if (fNParamsTPCInt<1) return; - fNSegZTPCInt = 0; - fSegRTPCInt = new Float_t[fNParamsTPCInt]; - float *tmpbufF = new float[fNParamsTPCInt+1]; - int *tmpbufI = new int[fNParamsTPCInt+1]; - int *tmpbufI1 = new int[fNParamsTPCInt+1]; - // - // count number of Z slices and number of R slices in each Z slice - for (int ip=0;ipGetBoundMax(2)-GetParamTPCInt(ip-1)->GetBoundMax(2))>kSafety) { // new Z slice - tmpbufF[fNSegZTPCInt] = GetParamTPCInt(ip)->GetBoundMax(2); // - tmpbufI[fNSegZTPCInt] = 0; - tmpbufI1[fNSegZTPCInt++] = ip; - } - fSegRTPCInt[ip] = GetParamTPCInt(ip)->GetBoundMax(0); // upper R - tmpbufI[fNSegZTPCInt-1]++; - } - // - fSegZTPCInt = new Float_t[fNSegZTPCInt]; - fSegZIdTPCInt = new Int_t[fNSegZTPCInt]; - fNSegRTPCInt = new Int_t[fNSegZTPCInt]; - for (int iz=0;izGetBoundMin(2); - fMaxZTPCInt = GetParamTPCInt(fNParamsTPCInt-1)->GetBoundMax(2); - fMaxRTPCInt = GetParamTPCInt(fNParamsTPCInt-1)->GetBoundMax(0); - // - delete[] tmpbufF; - delete[] tmpbufI; - delete[] tmpbufI1; - // - // -} - +//________________________________________________________________ void AliMagWrapCheb::SaveData(const char* outfile) const { // writes coefficients data to output text file @@ -861,8 +835,8 @@ void AliMagWrapCheb::SaveData(const char* outfile) const // // TPCInt part --------------------------------------------------------- fprintf(stream,"# Set of Chebyshev parameterizations for ALICE magnetic field\nSTART %s\n",GetName()); - fprintf(stream,"START TPCINT\n#Number of pieces\n%d\n",fNParamsTPCInt); - for (int ip=0;ipSaveData(stream); + fprintf(stream,"START TPCINT\n#Number of pieces\n%d\n",fNParamsTPC); + for (int ip=0;ipSaveData(stream); fprintf(stream,"#\nEND TPCINT\n"); // // Dip part --------------------------------------------------------- @@ -876,7 +850,7 @@ void AliMagWrapCheb::SaveData(const char* outfile) const // } -Int_t AliMagWrapCheb::SegmentDipDimension(float** seg,const TObjArray* par,int npar, int dim, +Int_t AliMagWrapCheb::SegmentDimension(float** seg,const TObjArray* par,int npar, int dim, float xmn,float xmx,float ymn,float ymx,float zmn,float zmx) { // find all boundaries in deimension dim for boxes in given region.