[u/mrichter/AliRoot.git] / STEER / AliMagFCheb.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

///////////////////////////////////////////////////////////////////////////////////
//                                                                               //
//  Wrapper for the set of mag.field parameterizations by Chebyshev polinomials  //
//  To obtain the field in cartesian coordinates/components use                  //
//    Field(float* xyz, float* bxyz);                                            //
//  For cylindrical coordinates/components:                                      //
//    FieldCyl(float* rphiz, float* brphiz)                                      //
//                                                                               //
//  For the moment only the solenoid part is parameterized in the volume defined //
//  by R<500, -550<Z<550 cm                                                      //
//                                                                               //
//  The region R<423 cm,  -343.3<Z<481.3 for 30kA and -343.3<Z<481.3 for 12kA    //
//  is parameterized using measured data while outside the Tosca calculation     //
//  is used (matched to data on the boundary of the measurements)                //
//                                                                               //
//  If the querried point is outside the validity region no the return values    //
//  for the field components are set to 0.                                       //
//                                                                               //
///////////////////////////////////////////////////////////////////////////////////

#include "AliLog.h"
#include "AliMagFCheb.h"

ClassImp(AliMagFCheb)


//__________________________________________________________________________________________
AliMagFCheb::AliMagFCheb() :
    TNamed(),
    fNParamsSol(0),
    fNSegZSol(0),
    fNParamsDip(0),
    fSegZSol(0),
    fSegRSol(0),
    fNSegRSol(0),
    fSegZIdSol(0),
    fMinZSol(0.),
    fMaxZSol(0.),
    fMaxRSol(0.),
    fParamsSol(0),
    fParamsDip(0)
{
  Init0();
}

AliMagFCheb::AliMagFCheb(const char* inputFile) :
    TNamed("Field Map", inputFile),
    fNParamsSol(0),
    fNSegZSol(0),
    fNParamsDip(0),
    fSegZSol(0),
    fSegRSol(0),
    fNSegRSol(0),
    fSegZIdSol(0),
    fMinZSol(0.),
    fMaxZSol(0.),
    fMaxRSol(0.),
    fParamsSol(0),
    fParamsDip(0)
{
  Init0();
  LoadData(inputFile);
}

//__________________________________________________________________________________________
AliMagFCheb::AliMagFCheb(const AliMagFCheb& src)
    : TNamed(src), 
      fNParamsSol(src.fNParamsSol),
      fNSegZSol(src.fNSegZSol), 
      fNParamsDip(src.fNParamsDip), 
      fSegZSol(0), 
      fSegRSol(0), 
      fNSegRSol(0), 
      fSegZIdSol(0),
      fMinZSol(src.fMinZSol), 
      fMaxZSol(src.fMaxZSol), 
      fMaxRSol(src.fMaxRSol), 
      fParamsSol(0), 
      fParamsDip(0)
{
// Copy constructor
  if (src.fSegZSol) {
    fSegZSol = new Float_t[fNSegZSol];
    for (int i=fNSegZSol;i--;) fSegZSol[i] = src.fSegZSol[i];
  }
  if (src.fSegRSol) {
    fSegRSol = new Float_t[fNParamsSol];
    for (int i=fNParamsSol;i--;) fSegRSol[i] = src.fSegRSol[i];
  }
  if (src.fNSegRSol) {
    fNSegRSol = new Int_t[fNSegZSol];
    for (int i=fNSegZSol;i--;) fNSegRSol[i] = src.fNSegRSol[i];
  }
  if (src.fSegZIdSol) {
    fSegZIdSol = new Int_t[fNSegZSol];
    for (int i=fNSegZSol;i--;) fSegZIdSol[i] = src.fSegZIdSol[i];
  }
  if (src.fParamsSol) {
    fParamsSol = new TObjArray(fNParamsSol);
    for (int i=0;i<fNParamsSol;i++) {
      AliCheb3D* pr = src.GetParamSol(i);
      if (pr) fParamsSol->AddAtAndExpand(new AliCheb3D(*pr),i);
    }
  }
  if (src.fParamsDip) {
    fParamsDip = new TObjArray(fNParamsDip);
    for (int i=0;i<fNParamsDip;i++) {
      AliCheb3D* pr = src.GetParamDip(i);
      if (pr) fParamsDip->AddAtAndExpand(new AliCheb3D(*pr),i);
    }
  }
  //
}


AliMagFCheb& AliMagFCheb::operator=(const AliMagFCheb& rhs)
{
// Assignment operator
    if (this != &rhs) {  
	Clear();
	SetName(rhs.GetName());
	SetTitle(rhs.GetTitle());
	fNParamsSol = rhs.fNParamsSol;
	fNSegZSol   = rhs.fNSegZSol;
	fMinZSol    = rhs.fMinZSol;
	fMaxZSol    = rhs.fMaxZSol;
	fMaxRSol    = rhs.fMaxRSol;
	fNParamsDip = rhs.fNParamsDip;
	fSegZSol    = fSegRSol    = 0; 
	fNSegRSol   = fSegZIdSol  = 0;
	fParamsSol  = fParamsDip  = 0;  
	//
	if (rhs.fSegZSol) {
	    fSegZSol = new Float_t[fNSegZSol];
	    for (int i=fNSegZSol;i--;) fSegZSol[i] = rhs.fSegZSol[i];
	}
	if (rhs.fSegRSol) {
	    fSegRSol = new Float_t[fNParamsSol];
	    for (int i=fNParamsSol;i--;) fSegRSol[i] = rhs.fSegRSol[i];
	}
	if (rhs.fNSegRSol) {
	    fNSegRSol = new Int_t[fNSegZSol];
	    for (int i=fNSegZSol;i--;) fNSegRSol[i] = rhs.fNSegRSol[i];
	}
	if (rhs.fSegZIdSol) {
	    fSegZIdSol = new Int_t[fNSegZSol];
	    for (int i=fNSegZSol;i--;) fSegZIdSol[i] = rhs.fSegZIdSol[i];
	}
	if (rhs.fParamsSol) {
	    fParamsSol = new TObjArray(fNParamsSol);
	    for (int i=0;i<fNParamsSol;i++) {
		AliCheb3D* pr = rhs.GetParamSol(i);
		if (pr) fParamsSol->AddAtAndExpand(new AliCheb3D(*pr),i);
	    }
	}
	if (rhs.fParamsDip) {
	    fParamsDip = new TObjArray(fNParamsDip);
	    for (int i=0;i<fNParamsDip;i++) {
		AliCheb3D* pr = rhs.GetParamDip(i);
		if (pr) fParamsDip->AddAtAndExpand(new AliCheb3D(*pr),i);
	    }
	}
    }
    return *this;  
    //
}


//__________________________________________________________________________________________
AliMagFCheb::~AliMagFCheb()
{
  if (fNParamsSol) {
    delete fParamsSol;
    delete[] fSegZSol;
    delete[] fSegRSol;
    delete[] fNSegRSol;
    delete[] fSegZIdSol;
  }
  //
  // Dipole part ...
  if (fNParamsDip) {
    delete fParamsDip;
  }
}

//__________________________________________________________________________________________
void AliMagFCheb::Init0()
{
  // Solenoid part
  fNParamsSol = 0;
  fNSegZSol   = 0;
  //
  fSegZSol    = 0;
  fSegRSol    = 0;
  //
  fNSegRSol   = 0;
  fSegZIdSol  = 0;
  //
  fMinZSol = fMaxZSol = fMaxRSol = fMaxRSol = 0;
  fParamsSol = 0;
  //
  // Dipole part ...
  fNParamsDip = 0;
  fParamsDip  = 0;
  //
}

//__________________________________________________________________________________________
void AliMagFCheb::AddParamSol(AliCheb3D* param)
{
  // adds new parameterization piece for Sol
  // NOTE: pieces must be added strictly in increasing R then increasing Z order
  //
  if (!fParamsSol) fParamsSol = new TObjArray();
  fParamsSol->Add(param);
  fNParamsSol++;
  //
}

//__________________________________________________________________________________________
void AliMagFCheb::AddParamDip(AliCheb3D* param)
{
  // adds new parameterization piece for Dipole
  //
  if (!fParamsDip) fParamsDip = new TObjArray();
  fParamsDip->Add(param);
  fNParamsDip++;
  //
}

//__________________________________________________________________________________________
void AliMagFCheb::BuildTableSol()
{
  // build the indexes for each parameterization of Solenoid
  //
  const float kSafety=0.001;
  //
  fSegRSol   = new Float_t[fNParamsSol];
  float *tmpbufF  = new float[fNParamsSol+1];
  int   *tmpbufI  = new int[fNParamsSol+1];
  int   *tmpbufI1 = new int[fNParamsSol+1];
  //
  // count number of Z slices and number of R slices in each Z slice
  for (int ip=0;ip<fNParamsSol;ip++) {
    if (ip==0 || (GetParamSol(ip)->GetBoundMax(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]++;
  }
  //
  fSegZSol   = new Float_t[fNSegZSol];
  fSegZIdSol = new Int_t[fNSegZSol];
  fNSegRSol  = new Int_t[fNSegZSol];
  for (int iz=0;iz<fNSegZSol;iz++) {
    fSegZSol[iz]   = tmpbufF[iz];
    fNSegRSol[iz]  = tmpbufI[iz];
    fSegZIdSol[iz] = tmpbufI1[iz];
  }
  //
  fMinZSol = GetParamSol(0)->GetBoundMin(2);
  fMaxZSol = GetParamSol(fNParamsSol-1)->GetBoundMax(2);
  fMaxRSol = GetParamSol(fNParamsSol-1)->GetBoundMax(0);
  //
  delete[] tmpbufF;
  delete[] tmpbufI;
  delete[] tmpbufI1;
  //
  //
}

//__________________________________________________________________________________________
void AliMagFCheb::Field(Float_t *xyz, Float_t *b) const
{
  // compute field in cartesian coordinates
  float rphiz[3];
  if (xyz[2]>GetMaxZSol() || xyz[2]<GetMinZSol()) {for (int i=3;i--;) b[i]=0; return;}
  //
  // Sol region
  // convert coordinates to cyl system
  rphiz[0] = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);
  rphiz[1] = TMath::ATan2(xyz[1],xyz[0]);
  rphiz[2] = xyz[2];
  if (rphiz[0]>GetMaxRSol()) {for (int i=3;i--;) b[i]=0; return;}
  //
  FieldCylSol(rphiz,b);
  //
  // convert field to cartesian system
  float btr = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]);
  float psiPLUSphi = rphiz[1] + TMath::ATan2(b[1],b[0]);
  b[0] = btr*TMath::Cos(psiPLUSphi);
  b[1] = btr*TMath::Sin(psiPLUSphi);
  //
}

//__________________________________________________________________________________________
void AliMagFCheb::FieldCylSol(Float_t *rphiz, Float_t *b) const
{
  // compute Solenoid field in Cylindircal coordinates
  // note: the check for the point being inside the parameterized region is done outside
  float &r = rphiz[0];
  float &z = rphiz[2];
  int SolZId = 0;
  while (z>fSegZSol[SolZId]) ++SolZId;    // find Z segment
  int SolRId = fSegZIdSol[SolZId];        // first R segment for this Z
  while (r>fSegRSol[SolRId]) ++SolRId;    // find R segment
  GetParamSol( SolRId )->Eval(rphiz,b);
  //
}

//__________________________________________________________________________________________
void AliMagFCheb::Print(Option_t *) const
{
  printf("Alice magnetic field parameterized by Chebyshev polynomials\n");
  printf("Segmentation for Solenoid (%+.2f<Z<%+.2f cm | R<%.2f cm)\n",fMinZSol,fMaxZSol,fMaxRSol);
  //
  for (int iz=0;iz<fNSegZSol;iz++) {
    AliCheb3D* param = GetParamSol( fSegZIdSol[iz] );
    printf("*** Z Segment %2d (%+7.2f<Z<%+7.2f)\t***\n",iz,param->GetBoundMin(2),param->GetBoundMax(2));
    for (int ir=0;ir<fNSegRSol[iz];ir++) {
      param = GetParamSol( fSegZIdSol[iz]+ir );
      printf("    R Segment %2d (%+7.2f<R<%+7.2f, Precision: %.1e) (ID=%2d)\n",ir, param->GetBoundMin(0),param->GetBoundMax(0),
	     param->GetPrecision(),fSegZIdSol[iz]+ir);
    }
  }
}

//_______________________________________________
#ifdef  _INC_CREATION_ALICHEB3D_
void AliMagFCheb::SaveData(const char* outfile) const
{
  // writes coefficients data to output text file
  TString strf = outfile;
  gSystem->ExpandPathName(strf);
  FILE* stream = fopen(strf,"w+");
  //
  // Sol part
  fprintf(stream,"# Set of Chebyshev parameterizations for ALICE magnetic field\nSTART %s\n",GetName());
  fprintf(stream,"START SOLENOID\n#Number of pieces\n%d\n",fNParamsSol);
  for (int ip=0;ip<fNParamsSol;ip++) GetParamSol(ip)->SaveData(stream);
  fprintf(stream,"#\nEND SOLENOID\n");
  //
  // Dip part
  fprintf(stream,"START DIPOLE\n#Number of pieces\n%d\n",fNParamsDip);
  for (int ip=0;ip<fNParamsDip;ip++) GetParamDip(ip)->SaveData(stream);
  fprintf(stream,"#\nEND DIPOLE\n");
  //
  fprintf(stream,"#\nEND %s\n",GetName());
  //
  fclose(stream);
  //
}
#endif

//_______________________________________________
void AliMagFCheb::LoadData(const char* inpfile)
{
  // read coefficients data from the text file
  //
  TString strf = inpfile;
  gSystem->ExpandPathName(strf);
  FILE* stream = fopen(strf,"r");
  if (!stream) {
    printf("Did not find input file %s\n",strf.Data());
    return;
  }
  //
  TString buffs;
  AliCheb3DCalc::ReadLine(buffs,stream);
  if (!buffs.BeginsWith("START")) {Error("LoadData","Expected: \"START <name>\", found \"%s\"\nStop\n",buffs.Data());exit(1);}
  if (buffs.First(' ')>0) SetName(buffs.Data()+buffs.First(' ')+1);
  //
  // Solenoid part
  AliCheb3DCalc::ReadLine(buffs,stream);
  if (!buffs.BeginsWith("START SOLENOID")) {Error("LoadData","Expected: \"START SOLENOID\", found \"%s\"\nStop\n",buffs.Data());exit(1);}
  AliCheb3DCalc::ReadLine(buffs,stream); // nparam
  int nparSol = buffs.Atoi(); 
  //
  for (int ip=0;ip<nparSol;ip++) {
    AliCheb3D* cheb = new AliCheb3D();
    cheb->LoadData(stream);
    AddParamSol(cheb);
  }
  //
  AliCheb3DCalc::ReadLine(buffs,stream);
  if (!buffs.BeginsWith("END SOLENOID")) {Error("LoadData","Expected \"END SOLENOID\", found \"%s\"\nStop\n",buffs.Data());exit(1);}
  //
  // Dipole part
  AliCheb3DCalc::ReadLine(buffs,stream);
  if (!buffs.BeginsWith("START DIPOLE")) {Error("LoadData","Expected: \"START DIPOLE\", found \"%s\"\nStop\n",buffs.Data());exit(1);}
  AliCheb3DCalc::ReadLine(buffs,stream); // nparam
  int nparDip = buffs.Atoi();  
  //
  for (int ip=0;ip<nparDip;ip++) {
    AliCheb3D* cheb = new AliCheb3D();
    cheb->LoadData(stream);
    AddParamDip(cheb);
  }
  //
  AliCheb3DCalc::ReadLine(buffs,stream);
  if (!buffs.BeginsWith("END DIPOLE")) {Error("LoadData","Expected \"END DIPOLE\", found \"%s\"\nStop\n",GetName(),buffs.Data());exit(1);}
  //
  AliCheb3DCalc::ReadLine(buffs,stream);
  if (!buffs.BeginsWith("END") || !buffs.Contains(GetName())) {Error("LoadData","Expected: \"END %s\", found \"%s\"\nStop\n",GetName(),buffs.Data());exit(1);}
  //
  fclose(stream);
  BuildTableSol();
  //  BuildDipTable();
  printf("Loaded magnetic field \"%s\" from %s\n",GetName(),strf.Data());
  //
}
Commit	Line	Data
0eea9d4d	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/* $Id$ */
	17
	18	///////////////////////////////////////////////////////////////////////////////////
	19	// //
	20	// Wrapper for the set of mag.field parameterizations by Chebyshev polinomials //
	21	// To obtain the field in cartesian coordinates/components use //
	22	// Field(float* xyz, float* bxyz); //
	23	// For cylindrical coordinates/components: //
	24	// FieldCyl(float* rphiz, float* brphiz) //
	25	// //
	26	// For the moment only the solenoid part is parameterized in the volume defined //
	27	// by R<500, -550<Z<550 cm //
	28	// //
	29	// The region R<423 cm, -343.3<Z<481.3 for 30kA and -343.3<Z<481.3 for 12kA //
	30	// is parameterized using measured data while outside the Tosca calculation //
	31	// is used (matched to data on the boundary of the measurements) //
	32	// //
	33	// If the querried point is outside the validity region no the return values //
	34	// for the field components are set to 0. //
	35	// //
	36	///////////////////////////////////////////////////////////////////////////////////
	37
0bc7b414	38	#include "AliLog.h"
0eea9d4d	39	#include "AliMagFCheb.h"
	40
	41	ClassImp(AliMagFCheb)
	42
	43
	44
	45	//__________________________________________________________________________________________
0bc7b414	46	AliMagFCheb::AliMagFCheb() :
	47	TNamed(),
	48	fNParamsSol(0),
	49	fNSegZSol(0),
	50	fNParamsDip(0),
	51	fSegZSol(0),
	52	fSegRSol(0),
	53	fNSegRSol(0),
	54	fSegZIdSol(0),
	55	fMinZSol(0.),
	56	fMaxZSol(0.),
	57	fMaxRSol(0.),
	58	fParamsSol(0),
	59	fParamsDip(0)
	60	{
	61	Init0();
	62	}
	63
0eea9d4d	64	AliMagFCheb::AliMagFCheb(const char* inputFile) :
	65	TNamed("Field Map", inputFile),
	66	fNParamsSol(0),
	67	fNSegZSol(0),
	68	fNParamsDip(0),
	69	fSegZSol(0),
	70	fSegRSol(0),
	71	fNSegRSol(0),
	72	fSegZIdSol(0),
	73	fMinZSol(0.),
	74	fMaxZSol(0.),
0bc7b414	75	fMaxRSol(0.),
	76	fParamsSol(0),
	77	fParamsDip(0)
0eea9d4d	78	{
	79	Init0();
	80	LoadData(inputFile);
	81	}
	82
c437b1a5	83	//__________________________________________________________________________________________
	84	AliMagFCheb::AliMagFCheb(const AliMagFCheb& src)
	85	: TNamed(src),
	86	fNParamsSol(src.fNParamsSol),
	87	fNSegZSol(src.fNSegZSol),
	88	fNParamsDip(src.fNParamsDip),
	89	fSegZSol(0),
	90	fSegRSol(0),
	91	fNSegRSol(0),
	92	fSegZIdSol(0),
	93	fMinZSol(src.fMinZSol),
	94	fMaxZSol(src.fMaxZSol),
	95	fMaxRSol(src.fMaxRSol),
	96	fParamsSol(0),
	97	fParamsDip(0)
0bc7b414	98	{
c437b1a5	99	// Copy constructor
	100	if (src.fSegZSol) {
	101	fSegZSol = new Float_t[fNSegZSol];
	102	for (int i=fNSegZSol;i--;) fSegZSol[i] = src.fSegZSol[i];
	103	}
	104	if (src.fSegRSol) {
	105	fSegRSol = new Float_t[fNParamsSol];
	106	for (int i=fNParamsSol;i--;) fSegRSol[i] = src.fSegRSol[i];
	107	}
	108	if (src.fNSegRSol) {
	109	fNSegRSol = new Int_t[fNSegZSol];
	110	for (int i=fNSegZSol;i--;) fNSegRSol[i] = src.fNSegRSol[i];
	111	}
	112	if (src.fSegZIdSol) {
	113	fSegZIdSol = new Int_t[fNSegZSol];
	114	for (int i=fNSegZSol;i--;) fSegZIdSol[i] = src.fSegZIdSol[i];
	115	}
	116	if (src.fParamsSol) {
	117	fParamsSol = new TObjArray(fNParamsSol);
	118	for (int i=0;i<fNParamsSol;i++) {
	119	AliCheb3D* pr = src.GetParamSol(i);
	120	if (pr) fParamsSol->AddAtAndExpand(new AliCheb3D(*pr),i);
	121	}
	122	}
	123	if (src.fParamsDip) {
	124	fParamsDip = new TObjArray(fNParamsDip);
	125	for (int i=0;i<fNParamsDip;i++) {
	126	AliCheb3D* pr = src.GetParamDip(i);
	127	if (pr) fParamsDip->AddAtAndExpand(new AliCheb3D(*pr),i);
	128	}
	129	}
0bc7b414	130	//
0bc7b414	131	}
0bc7b414	132
c437b1a5	133
c437b1a5	134	AliMagFCheb& AliMagFCheb::operator=(const AliMagFCheb& rhs)
0bc7b414	135	{
c437b1a5	136	// Assignment operator
	137	if (this != &rhs) {
	138	Clear();
	139	SetName(rhs.GetName());
	140	SetTitle(rhs.GetTitle());
	141	fNParamsSol = rhs.fNParamsSol;
	142	fNSegZSol = rhs.fNSegZSol;
	143	fMinZSol = rhs.fMinZSol;
	144	fMaxZSol = rhs.fMaxZSol;
	145	fMaxRSol = rhs.fMaxRSol;
	146	fNParamsDip = rhs.fNParamsDip;
	147	fSegZSol = fSegRSol = 0;
	148	fNSegRSol = fSegZIdSol = 0;
	149	fParamsSol = fParamsDip = 0;
	150	//
	151	if (rhs.fSegZSol) {
	152	fSegZSol = new Float_t[fNSegZSol];
	153	for (int i=fNSegZSol;i--;) fSegZSol[i] = rhs.fSegZSol[i];
	154	}
	155	if (rhs.fSegRSol) {
	156	fSegRSol = new Float_t[fNParamsSol];
	157	for (int i=fNParamsSol;i--;) fSegRSol[i] = rhs.fSegRSol[i];
	158	}
	159	if (rhs.fNSegRSol) {
	160	fNSegRSol = new Int_t[fNSegZSol];
	161	for (int i=fNSegZSol;i--;) fNSegRSol[i] = rhs.fNSegRSol[i];
	162	}
	163	if (rhs.fSegZIdSol) {
	164	fSegZIdSol = new Int_t[fNSegZSol];
	165	for (int i=fNSegZSol;i--;) fSegZIdSol[i] = rhs.fSegZIdSol[i];
	166	}
	167	if (rhs.fParamsSol) {
	168	fParamsSol = new TObjArray(fNParamsSol);
	169	for (int i=0;i<fNParamsSol;i++) {
	170	AliCheb3D* pr = rhs.GetParamSol(i);
	171	if (pr) fParamsSol->AddAtAndExpand(new AliCheb3D(*pr),i);
	172	}
	173	}
	174	if (rhs.fParamsDip) {
	175	fParamsDip = new TObjArray(fNParamsDip);
	176	for (int i=0;i<fNParamsDip;i++) {
	177	AliCheb3D* pr = rhs.GetParamDip(i);
	178	if (pr) fParamsDip->AddAtAndExpand(new AliCheb3D(*pr),i);
	179	}
	180	}
	181	}
	182	return *this;
	183	//
0bc7b414	184	}
0bc7b414	185
c437b1a5	186
0eea9d4d	187	//__________________________________________________________________________________________
	188	AliMagFCheb::~AliMagFCheb()
	189	{
	190	if (fNParamsSol) {
	191	delete fParamsSol;
	192	delete[] fSegZSol;
	193	delete[] fSegRSol;
	194	delete[] fNSegRSol;
	195	delete[] fSegZIdSol;
	196	}
	197	//
	198	// Dipole part ...
	199	if (fNParamsDip) {
	200	delete fParamsDip;
	201	}
	202	}
	203
	204	//__________________________________________________________________________________________
	205	void AliMagFCheb::Init0()
	206	{
	207	// Solenoid part
	208	fNParamsSol = 0;
	209	fNSegZSol = 0;
	210	//
	211	fSegZSol = 0;
	212	fSegRSol = 0;
	213	//
	214	fNSegRSol = 0;
	215	fSegZIdSol = 0;
	216	//
	217	fMinZSol = fMaxZSol = fMaxRSol = fMaxRSol = 0;
	218	fParamsSol = 0;
	219	//
	220	// Dipole part ...
	221	fNParamsDip = 0;
	222	fParamsDip = 0;
	223	//
	224	}
	225
	226	//__________________________________________________________________________________________
	227	void AliMagFCheb::AddParamSol(AliCheb3D* param)
	228	{
	229	// adds new parameterization piece for Sol
	230	// NOTE: pieces must be added strictly in increasing R then increasing Z order
	231	//
	232	if (!fParamsSol) fParamsSol = new TObjArray();
	233	fParamsSol->Add(param);
	234	fNParamsSol++;
	235	//
	236	}
	237
	238	//__________________________________________________________________________________________
	239	void AliMagFCheb::AddParamDip(AliCheb3D* param)
	240	{
	241	// adds new parameterization piece for Dipole
	242	//
	243	if (!fParamsDip) fParamsDip = new TObjArray();
	244	fParamsDip->Add(param);
	245	fNParamsDip++;
	246	//
	247	}
	248
	249	//__________________________________________________________________________________________
	250	void AliMagFCheb::BuildTableSol()
251	{
252	// build the indexes for each parameterization of Solenoid
253	//
254	const float kSafety=0.001;
255	//
256	fSegRSol = new Float_t[fNParamsSol];
257	float *tmpbufF = new float[fNParamsSol+1];
258	int *tmpbufI = new int[fNParamsSol+1];
259	int *tmpbufI1 = new int[fNParamsSol+1];
260	//
261	// count number of Z slices and number of R slices in each Z slice
262	for (int ip=0;ip<fNParamsSol;ip++) {
263	if (ip==0 \|\| (GetParamSol(ip)->GetBoundMax(2)-GetParamSol(ip-1)->GetBoundMax(2))>kSafety) { // new Z slice
264	tmpbufF[fNSegZSol] = GetParamSol(ip)->GetBoundMax(2); //
265	tmpbufI[fNSegZSol] = 0;
266	tmpbufI1[fNSegZSol++] = ip;
267	}
268	fSegRSol[ip] = GetParamSol(ip)->GetBoundMax(0); // upper R
269	tmpbufI[fNSegZSol-1]++;
270	}
271	//
272	fSegZSol = new Float_t[fNSegZSol];
273	fSegZIdSol = new Int_t[fNSegZSol];
274	fNSegRSol = new Int_t[fNSegZSol];
275	for (int iz=0;iz<fNSegZSol;iz++) {
276	fSegZSol[iz] = tmpbufF[iz];
277	fNSegRSol[iz] = tmpbufI[iz];
278	fSegZIdSol[iz] = tmpbufI1[iz];
279	}
280	//
281	fMinZSol = GetParamSol(0)->GetBoundMin(2);
282	fMaxZSol = GetParamSol(fNParamsSol-1)->GetBoundMax(2);
283	fMaxRSol = GetParamSol(fNParamsSol-1)->GetBoundMax(0);
284	//
285	delete[] tmpbufF;
286	delete[] tmpbufI;
287	delete[] tmpbufI1;
288	//
289	//
290	}
291
292	//__________________________________________________________________________________________
293	void AliMagFCheb::Field(Float_t xyz, Float_t b) const
294	{
295	// compute field in cartesian coordinates
296	float rphiz[3];
297	if (xyz[2]>GetMaxZSol() \|\| xyz[2]<GetMinZSol()) {for (int i=3;i--;) b[i]=0; return;}
298	//
299	// Sol region
300	// convert coordinates to cyl system
301	rphiz[0] = TMath::Sqrt(xyz[0]xyz[0]+xyz[1]xyz[1]);
302	rphiz[1] = TMath::ATan2(xyz[1],xyz[0]);
303	rphiz[2] = xyz[2];
304	if (rphiz[0]>GetMaxRSol()) {for (int i=3;i--;) b[i]=0; return;}
305	//
306	FieldCylSol(rphiz,b);
307	//
308	// convert field to cartesian system
309	float btr = TMath::Sqrt(b[0]b[0]+b[1]b[1]);
310	float psiPLUSphi = rphiz[1] + TMath::ATan2(b[1],b[0]);
311	b[0] = btr*TMath::Cos(psiPLUSphi);
312	b[1] = btr*TMath::Sin(psiPLUSphi);
313	//
314	}
315
316	//__________________________________________________________________________________________
317	void AliMagFCheb::FieldCylSol(Float_t rphiz, Float_t b) const
318	{
319	// compute Solenoid field in Cylindircal coordinates
320	// note: the check for the point being inside the parameterized region is done outside
321	float &r = rphiz[0];
322	float &z = rphiz[2];
323	int SolZId = 0;
324	while (z>fSegZSol[SolZId]) ++SolZId; // find Z segment
325	int SolRId = fSegZIdSol[SolZId]; // first R segment for this Z
326	while (r>fSegRSol[SolRId]) ++SolRId; // find R segment
327	GetParamSol( SolRId )->Eval(rphiz,b);
328	//
329	}
330
331	//__________________________________________________________________________________________
332	void AliMagFCheb::Print(Option_t *) const
333	{
334	printf("Alice magnetic field parameterized by Chebyshev polynomials\n");
335	printf("Segmentation for Solenoid (%+.2f<Z<%+.2f cm \| R<%.2f cm)\n",fMinZSol,fMaxZSol,fMaxRSol);
336	//
337	for (int iz=0;iz<fNSegZSol;iz++) {
338	AliCheb3D* param = GetParamSol( fSegZIdSol[iz] );
339	printf("* Z Segment %2d (%+7.2f<Z<%+7.2f)\t*\n",iz,param->GetBoundMin(2),param->GetBoundMax(2));
340	for (int ir=0;ir<fNSegRSol[iz];ir++) {
341	param = GetParamSol( fSegZIdSol[iz]+ir );
342	printf(" R Segment %2d (%+7.2f<R<%+7.2f, Precision: %.1e) (ID=%2d)\n",ir, param->GetBoundMin(0),param->GetBoundMax(0),
343	param->GetPrecision(),fSegZIdSol[iz]+ir);
344	}
345	}
346	}
347
348	//_______________________________________________
349	#ifdef _INC_CREATION_ALICHEB3D_
350	void AliMagFCheb::SaveData(const char* outfile) const
351	{
352	// writes coefficients data to output text file
353	TString strf = outfile;
354	gSystem->ExpandPathName(strf);
355	FILE* stream = fopen(strf,"w+");
356	//
357	// Sol part
358	fprintf(stream,"# Set of Chebyshev parameterizations for ALICE magnetic field\nSTART %s\n",GetName());
359	fprintf(stream,"START SOLENOID\n#Number of pieces\n%d\n",fNParamsSol);
360	for (int ip=0;ip<fNParamsSol;ip++) GetParamSol(ip)->SaveData(stream);
361	fprintf(stream,"#\nEND SOLENOID\n");
362	//
363	// Dip part
364	fprintf(stream,"START DIPOLE\n#Number of pieces\n%d\n",fNParamsDip);
365	for (int ip=0;ip<fNParamsDip;ip++) GetParamDip(ip)->SaveData(stream);
366	fprintf(stream,"#\nEND DIPOLE\n");
367	//
368	fprintf(stream,"#\nEND %s\n",GetName());
369	//
370	fclose(stream);
371	//
372	}
373	#endif
374
375	//_______________________________________________
376	void AliMagFCheb::LoadData(const char* inpfile)
377	{
378	// read coefficients data from the text file
379	//
380	TString strf = inpfile;
381	gSystem->ExpandPathName(strf);
382	FILE* stream = fopen(strf,"r");
383	if (!stream) {
384	printf("Did not find input file %s\n",strf.Data());
385	return;
386	}
387	//
388	TString buffs;
389	AliCheb3DCalc::ReadLine(buffs,stream);
390	if (!buffs.BeginsWith("START")) {Error("LoadData","Expected: \"START <name>\", found \"%s\"\nStop\n",buffs.Data());exit(1);}
391	if (buffs.First(' ')>0) SetName(buffs.Data()+buffs.First(' ')+1);
392	//
393	// Solenoid part
394	AliCheb3DCalc::ReadLine(buffs,stream);
395	if (!buffs.BeginsWith("START SOLENOID")) {Error("LoadData","Expected: \"START SOLENOID\", found \"%s\"\nStop\n",buffs.Data());exit(1);}
396	AliCheb3DCalc::ReadLine(buffs,stream); // nparam
397	int nparSol = buffs.Atoi();
398	//
399	for (int ip=0;ip<nparSol;ip++) {
400	AliCheb3D* cheb = new AliCheb3D();
401	cheb->LoadData(stream);
402	AddParamSol(cheb);
403	}
404	//
405	AliCheb3DCalc::ReadLine(buffs,stream);
406	if (!buffs.BeginsWith("END SOLENOID")) {Error("LoadData","Expected \"END SOLENOID\", found \"%s\"\nStop\n",buffs.Data());exit(1);}
407	//
408	// Dipole part
409	AliCheb3DCalc::ReadLine(buffs,stream);
410	if (!buffs.BeginsWith("START DIPOLE")) {Error("LoadData","Expected: \"START DIPOLE\", found \"%s\"\nStop\n",buffs.Data());exit(1);}
411	AliCheb3DCalc::ReadLine(buffs,stream); // nparam
412	int nparDip = buffs.Atoi();
413	//
414	for (int ip=0;ip<nparDip;ip++) {
415	AliCheb3D* cheb = new AliCheb3D();
416	cheb->LoadData(stream);
417	AddParamDip(cheb);
418	}
419	//
420	AliCheb3DCalc::ReadLine(buffs,stream);
421	if (!buffs.BeginsWith("END DIPOLE")) {Error("LoadData","Expected \"END DIPOLE\", found \"%s\"\nStop\n",GetName(),buffs.Data());exit(1);}
422	//
423	AliCheb3DCalc::ReadLine(buffs,stream);
424	if (!buffs.BeginsWith("END") \|\| !buffs.Contains(GetName())) {Error("LoadData","Expected: \"END %s\", found \"%s\"\nStop\n",GetName(),buffs.Data());exit(1);}
425	//
426	fclose(stream);
427	BuildTableSol();
428	// BuildDipTable();
429	printf("Loaded magnetic field \"%s\" from %s\n",GetName(),strf.Data());
430	//
431	}