[u/mrichter/AliRoot.git] / STEER / AliMagF.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.                  *
 **************************************************************************/


#include <TClass.h>
#include <TFile.h>
#include <TSystem.h>

#include "AliMagF.h"
#include "AliMagWrapCheb.h"
#include "AliLog.h"

ClassImp(AliMagF)

const Double_t AliMagF::fgkSol2DipZ    =  -700.;  

//_______________________________________________________________________
AliMagF::AliMagF():
  TVirtualMagField(),
  fMeasuredMap(0),
  fMapType(k5kG),
  fSolenoid(0),
  fBeamType(kNoBeamField),
  fBeamEnergy(0),
  //
  fInteg(0),
  fPrecInteg(0),
  fFactorSol(1.),
  fFactorDip(1.),
  fMax(15),
  fDipoleOFF(kFALSE),
  //
  fQuadGradient(0),
  fDipoleField(0),
  fCCorrField(0), 
  fACorr1Field(0),
  fACorr2Field(0),
  fParNames("","")
{
  // Default constructor
  //
}

//_______________________________________________________________________
AliMagF::AliMagF(const char *name, const char* title, Int_t integ, 
		 Double_t factorSol, Double_t factorDip, 
		 Double_t fmax, BMap_t maptype, const char* path,
		 BeamType_t bt, Double_t be):
  TVirtualMagField(name),
  fMeasuredMap(0),
  fMapType(maptype),
  fSolenoid(0),
  fBeamType(bt),
  fBeamEnergy(be),
  //
  fInteg(integ),
  fPrecInteg(1),
  fFactorSol(1.),
  fFactorDip(1.),
  fMax(fmax),
  fDipoleOFF(factorDip==0.),
  //
  fQuadGradient(0),
  fDipoleField(0),
  fCCorrField(0), 
  fACorr1Field(0),
  fACorr2Field(0),
  fParNames("","")
{
  // Initialize the field with Geant integration option "integ" and max field "fmax,
  // Impose scaling of parameterized L3 field by factorSol and of dipole by factorDip.
  // The "be" is the energy of the beam in GeV/nucleon
  //
  SetTitle(title);
  if(integ<0 || integ > 2) {
    AliWarning(Form("Invalid magnetic field flag: %5d; Helix tracking chosen instead",integ));
    fInteg = 2;
  }
  if (fInteg == 0) fPrecInteg = 0;
  //
  const char* parname = 0;
  //  
  if      (fMapType == k2kG) parname = fDipoleOFF ? "Sol12_Dip0_Hole":"Sol12_Dip6_Hole";
  else if (fMapType == k5kG) parname = fDipoleOFF ? "Sol30_Dip0_Hole":"Sol30_Dip6_Hole";
  else if (fMapType == k5kGUniform) parname = "Sol30_Dip6_Uniform";
  else AliFatal(Form("Unknown field identifier %d is requested\n",fMapType));
  //
  SetDataFileName(path);
  SetParamName(parname);
  //
  LoadParameterization();
  InitMachineField(fBeamType,fBeamEnergy);
  double xyz[3]={0.,0.,0.};
  fSolenoid = GetBz(xyz);
  SetFactorSol(factorSol);
  SetFactorDip(factorDip);
}

//_______________________________________________________________________
AliMagF::AliMagF(const AliMagF &src):
  TVirtualMagField(src),
  fMeasuredMap(0),
  fMapType(src.fMapType),
  fSolenoid(src.fSolenoid),
  fBeamType(src.fBeamType),
  fBeamEnergy(src.fBeamEnergy),
  fInteg(src.fInteg),
  fPrecInteg(src.fPrecInteg),
  fFactorSol(src.fFactorSol),
  fFactorDip(src.fFactorDip),
  fMax(src.fMax),
  fDipoleOFF(src.fDipoleOFF),
  fQuadGradient(src.fQuadGradient),
  fDipoleField(src.fDipoleField),
  fCCorrField(src.fCCorrField), 
  fACorr1Field(src.fACorr1Field),
  fACorr2Field(src.fACorr2Field),
  fParNames(src.fParNames)
{
  if (src.fMeasuredMap) fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
}

//_______________________________________________________________________
AliMagF::~AliMagF()
{
  delete fMeasuredMap;
}

//_______________________________________________________________________
Bool_t AliMagF::LoadParameterization()
{
  if (fMeasuredMap) {
    AliError(Form("Field data %s are already loaded from %s\n",GetParamName(),GetDataFileName()));
    return kTRUE;
  }
  //
  char* fname = gSystem->ExpandPathName(GetDataFileName());
  TFile* file = TFile::Open(fname);
  if (!file) {
    AliError(Form("Failed to open magnetic field data file %s\n",fname)); 
    return kFALSE;
  }
  //
  fMeasuredMap = dynamic_cast<AliMagWrapCheb*>(file->Get(GetParamName()));
  if (!fMeasuredMap) {
    AliError(Form("Did not find field %s in %s\n",GetParamName(),fname)); 
    return kFALSE;
  }
  file->Close();
  delete file;
  return kTRUE;
}


//_______________________________________________________________________
void AliMagF::Field(const Double_t *xyz, Double_t *b)
{
  // Method to calculate the field at point  xyz
  //
  //  b[0]=b[1]=b[2]=0.0;
  if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
    fMeasuredMap->Field(xyz,b);
    if (xyz[2]>fgkSol2DipZ || fDipoleOFF) for (int i=3;i--;) b[i] *= fFactorSol;
    else                                  for (int i=3;i--;) b[i] *= fFactorDip;    
  }
  else MachineField(xyz, b);
  //
}

//_______________________________________________________________________
Double_t AliMagF::GetBz(const Double_t *xyz) const
{
  // Method to calculate the field at point  xyz
  //
  if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
    double bz = fMeasuredMap->GetBz(xyz);
    return (xyz[2]>fgkSol2DipZ || fDipoleOFF) ? bz*fFactorSol : bz*fFactorDip;    
  }
  else return 0.;
}

//_______________________________________________________________________
AliMagF& AliMagF::operator=(const AliMagF& src)
{
  if (this != &src && src.fMeasuredMap) { 
    if (fMeasuredMap) delete fMeasuredMap;
    fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
    SetName(src.GetName());
    fSolenoid    = src.fSolenoid;
    fBeamType    = src.fBeamType;
    fBeamEnergy  = src.fBeamEnergy;
    fInteg       = src.fInteg;
    fPrecInteg   = src.fPrecInteg;
    fFactorSol   = src.fFactorSol;
    fFactorDip   = src.fFactorDip;
    fMax         = src.fMax;
    fDipoleOFF   = src.fDipoleOFF;
    fParNames    = src.fParNames;
  }
  return *this;
}

//_______________________________________________________________________
void AliMagF::InitMachineField(BeamType_t btype, Double_t benergy)
{
  if (btype==kNoBeamField || benergy<1.) {
    fQuadGradient = fDipoleField = fCCorrField = fACorr1Field = fACorr2Field = 0.;
    return;
  }
  //
  double rigScale = benergy/7000.;   // scale according to ratio of E/Enominal
  // for ions assume PbPb (with energy provided per nucleon) and account for A/Z
  if (btype == kBeamTypeAA) rigScale *= 208./82.;
  //
  fQuadGradient = 22.0002*rigScale;
  fDipoleField  = 37.8781*rigScale;
  //
  // SIDE C
  fCCorrField   = -9.6980;
  // SIDE A
  fACorr1Field  = -13.2247;
  fACorr2Field  =  11.7905;
  //
}

//_______________________________________________________________________
void AliMagF::MachineField(const Double_t *x, Double_t *b) const
{
  // ---- This is the ZDC part
  // Compansators for Alice Muon Arm Dipole
  const Double_t kBComp1CZ = 1075., kBComp1hDZ = 260./2., kBComp1SqR = 4.0*4.0; 
  const Double_t kBComp2CZ = 2049., kBComp2hDZ = 153./2., kBComp2SqR = 4.5*4.5; 
  //  
  const Double_t kTripQ1CZ = 2615., kTripQ1hDZ = 637./2., kTripQ1SqR = 3.5*3.5;
  const Double_t kTripQ2CZ = 3480., kTripQ2hDZ = 550./2., kTripQ2SqR = 3.5*3.5;
  const Double_t kTripQ3CZ = 4130., kTripQ3hDZ = 550./2., kTripQ3SqR = 3.5*3.5;
  const Double_t kTripQ4CZ = 5015., kTripQ4hDZ = 637./2., kTripQ4SqR = 3.5*3.5;
  //
  const Double_t kDip1CZ = 6310.8,  kDip1hDZ = 945./2., kDip1SqRC = 4.5*4.5, kDip1SqRA = 3.375*3.375;
  const Double_t kDip2CZ = 12640.3, kDip2hDZ = 945./2., kDip2SqRC = 4.5*4.5, kDip2SqRA = 3.75*3.75;
  const Double_t kDip2DXC = 9.7, kDip2DXA = 9.4;
  //
  double rad2 = x[0] * x[0] + x[1] * x[1];
  //
  b[0] = b[1] = b[2] = 0;
  //
  // SIDE C **************************************************
  if(x[2]<0.){  
    if(TMath::Abs(x[2]+kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
      b[0] = fCCorrField*fFactorDip;
    } 
    else if(TMath::Abs(x[2]+kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
      b[0] = fQuadGradient*x[1];
      b[1] = fQuadGradient*x[0];
    }
    else if(TMath::Abs(x[2]+kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
      b[0] = -fQuadGradient*x[1];
      b[1] = -fQuadGradient*x[0];
    }
    else if(TMath::Abs(x[2]+kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
      b[0] = -fQuadGradient*x[1];
      b[1] = -fQuadGradient*x[0];
    }
    else if(TMath::Abs(x[2]+kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
      b[0] = fQuadGradient*x[1];
      b[1] = fQuadGradient*x[0];
    }
    else if(TMath::Abs(x[2]+kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRC){
      b[1] = fDipoleField;
    }
    else if(TMath::Abs(x[2]+kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRC) {
      double dxabs = TMath::Abs(x[0])-kDip2DXC;
      if ( (dxabs*dxabs + x[1]*x[1])<kDip2SqRC) {
	b[1] = -fDipoleField;
      }
    }
  }
  //
  // SIDE A **************************************************
  else{        
    if(TMath::Abs(x[2]-kBComp1CZ)<kBComp1hDZ && rad2 < kBComp1SqR) {
      // Compensator magnet at z = 1075 m 
      b[0] = fACorr1Field*fFactorDip;
    }
    //
    if(TMath::Abs(x[2]-kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
      b[0] = fACorr2Field*fFactorDip;
    }
    else if(TMath::Abs(x[2]-kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
      b[0] = -fQuadGradient*x[1];
      b[1] = -fQuadGradient*x[0];
    }
    else if(TMath::Abs(x[2]-kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
      b[0] =  fQuadGradient*x[1];
      b[1] =  fQuadGradient*x[0];
    }
    else if(TMath::Abs(x[2]-kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
      b[0] =  fQuadGradient*x[1];
      b[1] =  fQuadGradient*x[0];
    }
    else if(TMath::Abs(x[2]-kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
      b[0] = -fQuadGradient*x[1];
      b[1] = -fQuadGradient*x[0];
    }
    else if(TMath::Abs(x[2]-kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRA){
      b[1] = -fDipoleField;
    }
    else if(TMath::Abs(x[2]-kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRA) {
      double dxabs = TMath::Abs(x[0])-kDip2DXA;
      if ( (dxabs*dxabs + x[1]*x[1])<kDip2SqRA) {
	b[1] = fDipoleField;
      }
    }
  }
  //
}

//_______________________________________________________________________
void AliMagF::GetTPCInt(const Double_t *xyz, Double_t *b) const
{
  // Method to calculate the integral of magnetic integral from xyz to nearest cathode plane
  b[0]=b[1]=b[2]=0.0;
  if (fMeasuredMap) {
    fMeasuredMap->GetTPCInt(xyz,b);
    for (int i=3;i--;) b[i] *= fFactorSol;
  }
}

//_______________________________________________________________________
void AliMagF::GetTPCIntCyl(const Double_t *rphiz, Double_t *b) const
{
  // Method to calculate the integral of magnetic integral from point to nearest cathode plane
  // in cylindrical coordiates ( -pi<phi<pi convention )
  b[0]=b[1]=b[2]=0.0;
  if (fMeasuredMap) {
    fMeasuredMap->GetTPCIntCyl(rphiz,b);
    for (int i=3;i--;) b[i] *= fFactorSol;
  }
}
Commit	Line	Data
4c039060	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
4c039060	16
db83d72f	17	#include <TClass.h>
	18	#include <TFile.h>
	19	#include <TSystem.h>
fe4da5cc	20
fe4da5cc	21	#include "AliMagF.h"
db83d72f	22	#include "AliMagWrapCheb.h"
db83d72f	23	#include "AliLog.h"
972ca52f	24
fe4da5cc	25	ClassImp(AliMagF)
fe4da5cc	26
9251fceb	27	const Double_t AliMagF::fgkSol2DipZ = -700.;
db83d72f	28
e2afb3b6	29	//_______________________________________________________________________
e2afb3b6	30	AliMagF::AliMagF():
db83d72f	31	TVirtualMagField(),
	32	fMeasuredMap(0),
	33	fMapType(k5kG),
	34	fSolenoid(0),
	35	fBeamType(kNoBeamField),
	36	fBeamEnergy(0),
db83d72f	37	//
e2afb3b6	38	fInteg(0),
db83d72f	39	fPrecInteg(0),
	40	fFactorSol(1.),
	41	fFactorDip(1.),
	42	fMax(15),
	43	fDipoleOFF(kFALSE),
e2afb3b6	44	//
db83d72f	45	fQuadGradient(0),
	46	fDipoleField(0),
	47	fCCorrField(0),
	48	fACorr1Field(0),
	49	fACorr2Field(0),
	50	fParNames("","")
	51	{
e2afb3b6	52	// Default constructor
	53	//
	54	}
	55
	56	//_______________________________________________________________________
db83d72f	57	AliMagF::AliMagF(const char name, const char title, Int_t integ,
	58	Double_t factorSol, Double_t factorDip,
	59	Double_t fmax, BMap_t maptype, const char* path,
9251fceb	60	BeamType_t bt, Double_t be):
db83d72f	61	TVirtualMagField(name),
	62	fMeasuredMap(0),
	63	fMapType(maptype),
	64	fSolenoid(0),
	65	fBeamType(bt),
	66	fBeamEnergy(be),
db83d72f	67	//
db83d72f	68	fInteg(integ),
604e0531	69	fPrecInteg(1),
db83d72f	70	fFactorSol(1.),
db83d72f	71	fFactorDip(1.),
972ca52f	72	fMax(fmax),
db83d72f	73	fDipoleOFF(factorDip==0.),
	74	//
	75	fQuadGradient(0),
	76	fDipoleField(0),
	77	fCCorrField(0),
	78	fACorr1Field(0),
	79	fACorr2Field(0),
	80	fParNames("","")
fe4da5cc	81	{
9251fceb	82	// Initialize the field with Geant integration option "integ" and max field "fmax,
	83	// Impose scaling of parameterized L3 field by factorSol and of dipole by factorDip.
	84	// The "be" is the energy of the beam in GeV/nucleon
aee8290b	85	//
db83d72f	86	SetTitle(title);
	87	if(integ<0 \|\| integ > 2) {
	88	AliWarning(Form("Invalid magnetic field flag: %5d; Helix tracking chosen instead",integ));
	89	fInteg = 2;
	90	}
	91	if (fInteg == 0) fPrecInteg = 0;
aee8290b	92	//
db83d72f	93	const char* parname = 0;
db83d72f	94	//
f04e7f5f	95	if (fMapType == k2kG) parname = fDipoleOFF ? "Sol12_Dip0_Hole":"Sol12_Dip6_Hole";
	96	else if (fMapType == k5kG) parname = fDipoleOFF ? "Sol30_Dip0_Hole":"Sol30_Dip6_Hole";
	97	else if (fMapType == k5kGUniform) parname = "Sol30_Dip6_Uniform";
	98	else AliFatal(Form("Unknown field identifier %d is requested\n",fMapType));
db83d72f	99	//
	100	SetDataFileName(path);
	101	SetParamName(parname);
	102	//
db83d72f	103	LoadParameterization();
db83d72f	104	InitMachineField(fBeamType,fBeamEnergy);
f04e7f5f	105	double xyz[3]={0.,0.,0.};
	106	fSolenoid = GetBz(xyz);
	107	SetFactorSol(factorSol);
	108	SetFactorDip(factorDip);
fe4da5cc	109	}
fe4da5cc	110
eeda4611	111	//_______________________________________________________________________
eeda4611	112	AliMagF::AliMagF(const AliMagF &src):
db83d72f	113	TVirtualMagField(src),
	114	fMeasuredMap(0),
	115	fMapType(src.fMapType),
	116	fSolenoid(src.fSolenoid),
	117	fBeamType(src.fBeamType),
	118	fBeamEnergy(src.fBeamEnergy),
eeda4611	119	fInteg(src.fInteg),
eeda4611	120	fPrecInteg(src.fPrecInteg),
db83d72f	121	fFactorSol(src.fFactorSol),
db83d72f	122	fFactorDip(src.fFactorDip),
eeda4611	123	fMax(src.fMax),
db83d72f	124	fDipoleOFF(src.fDipoleOFF),
	125	fQuadGradient(src.fQuadGradient),
	126	fDipoleField(src.fDipoleField),
	127	fCCorrField(src.fCCorrField),
	128	fACorr1Field(src.fACorr1Field),
	129	fACorr2Field(src.fACorr2Field),
	130	fParNames(src.fParNames)
eeda4611	131	{
db83d72f	132	if (src.fMeasuredMap) fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
eeda4611	133	}
eeda4611	134
e2afb3b6	135	//_______________________________________________________________________
db83d72f	136	AliMagF::~AliMagF()
ff66b122	137	{
db83d72f	138	delete fMeasuredMap;
	139	}
	140
	141	//_______________________________________________________________________
	142	Bool_t AliMagF::LoadParameterization()
	143	{
	144	if (fMeasuredMap) {
	145	AliError(Form("Field data %s are already loaded from %s\n",GetParamName(),GetDataFileName()));
	146	return kTRUE;
	147	}
ff66b122	148	//
db83d72f	149	char* fname = gSystem->ExpandPathName(GetDataFileName());
	150	TFile* file = TFile::Open(fname);
	151	if (!file) {
	152	AliError(Form("Failed to open magnetic field data file %s\n",fname));
	153	return kFALSE;
	154	}
ff66b122	155	//
db83d72f	156	fMeasuredMap = dynamic_cast<AliMagWrapCheb*>(file->Get(GetParamName()));
	157	if (!fMeasuredMap) {
	158	AliError(Form("Did not find field %s in %s\n",GetParamName(),fname));
	159	return kFALSE;
	160	}
	161	file->Close();
	162	delete file;
	163	return kTRUE;
ff66b122	164	}
ff66b122	165
db83d72f	166
ff66b122	167	//_______________________________________________________________________
db83d72f	168	void AliMagF::Field(const Double_t xyz, Double_t b)
fe4da5cc	169	{
db83d72f	170	// Method to calculate the field at point xyz
aee8290b	171	//
9251fceb	172	// b[0]=b[1]=b[2]=0.0;
9251fceb	173	if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
db83d72f	174	fMeasuredMap->Field(xyz,b);
db83d72f	175	if (xyz[2]>fgkSol2DipZ \|\| fDipoleOFF) for (int i=3;i--;) b[i] *= fFactorSol;
9251fceb	176	else for (int i=3;i--;) b[i] *= fFactorDip;
db83d72f	177	}
9251fceb	178	else MachineField(xyz, b);
aee8290b	179	//
fe4da5cc	180	}
eeda4611	181
eeda4611	182	//_______________________________________________________________________
db83d72f	183	Double_t AliMagF::GetBz(const Double_t *xyz) const
eeda4611	184	{
db83d72f	185	// Method to calculate the field at point xyz
db83d72f	186	//
9251fceb	187	if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
	188	double bz = fMeasuredMap->GetBz(xyz);
	189	return (xyz[2]>fgkSol2DipZ \|\| fDipoleOFF) ? bzfFactorSol : bzfFactorDip;
db83d72f	190	}
9251fceb	191	else return 0.;
eeda4611	192	}
	193
	194	//_______________________________________________________________________
db83d72f	195	AliMagF& AliMagF::operator=(const AliMagF& src)
eeda4611	196	{
db83d72f	197	if (this != &src && src.fMeasuredMap) {
	198	if (fMeasuredMap) delete fMeasuredMap;
	199	fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
	200	SetName(src.GetName());
	201	fSolenoid = src.fSolenoid;
	202	fBeamType = src.fBeamType;
	203	fBeamEnergy = src.fBeamEnergy;
db83d72f	204	fInteg = src.fInteg;
	205	fPrecInteg = src.fPrecInteg;
	206	fFactorSol = src.fFactorSol;
	207	fFactorDip = src.fFactorDip;
	208	fMax = src.fMax;
	209	fDipoleOFF = src.fDipoleOFF;
	210	fParNames = src.fParNames;
	211	}
	212	return *this;
eeda4611	213	}
	214
	215	//_______________________________________________________________________
db83d72f	216	void AliMagF::InitMachineField(BeamType_t btype, Double_t benergy)
eeda4611	217	{
9251fceb	218	if (btype==kNoBeamField \|\| benergy<1.) {
db83d72f	219	fQuadGradient = fDipoleField = fCCorrField = fACorr1Field = fACorr2Field = 0.;
9251fceb	220	return;
db83d72f	221	}
db83d72f	222	//
9251fceb	223	double rigScale = benergy/7000.; // scale according to ratio of E/Enominal
	224	// for ions assume PbPb (with energy provided per nucleon) and account for A/Z
	225	if (btype == kBeamTypeAA) rigScale *= 208./82.;
	226	//
	227	fQuadGradient = 22.0002*rigScale;
	228	fDipoleField = 37.8781*rigScale;
	229	//
	230	// SIDE C
	231	fCCorrField = -9.6980;
	232	// SIDE A
	233	fACorr1Field = -13.2247;
	234	fACorr2Field = 11.7905;
db83d72f	235	//
eeda4611	236	}
eed8a1a2	237
db83d72f	238	//_______________________________________________________________________
db83d72f	239	void AliMagF::MachineField(const Double_t x, Double_t b) const
eed8a1a2	240	{
db83d72f	241	// ---- This is the ZDC part
9251fceb	242	// Compansators for Alice Muon Arm Dipole
	243	const Double_t kBComp1CZ = 1075., kBComp1hDZ = 260./2., kBComp1SqR = 4.0*4.0;
	244	const Double_t kBComp2CZ = 2049., kBComp2hDZ = 153./2., kBComp2SqR = 4.5*4.5;
	245	//
	246	const Double_t kTripQ1CZ = 2615., kTripQ1hDZ = 637./2., kTripQ1SqR = 3.5*3.5;
90ae20c9	247	const Double_t kTripQ2CZ = 3480., kTripQ2hDZ = 550./2., kTripQ2SqR = 3.5*3.5;
9251fceb	248	const Double_t kTripQ3CZ = 4130., kTripQ3hDZ = 550./2., kTripQ3SqR = 3.5*3.5;
9251fceb	249	const Double_t kTripQ4CZ = 5015., kTripQ4hDZ = 637./2., kTripQ4SqR = 3.5*3.5;
db83d72f	250	//
9251fceb	251	const Double_t kDip1CZ = 6310.8, kDip1hDZ = 945./2., kDip1SqRC = 4.54.5, kDip1SqRA = 3.3753.375;
	252	const Double_t kDip2CZ = 12640.3, kDip2hDZ = 945./2., kDip2SqRC = 4.54.5, kDip2SqRA = 3.753.75;
	253	const Double_t kDip2DXC = 9.7, kDip2DXA = 9.4;
db83d72f	254	//
	255	double rad2 = x[0] * x[0] + x[1] * x[1];
	256	//
9251fceb	257	b[0] = b[1] = b[2] = 0;
9251fceb	258	//
db83d72f	259	// SIDE C **************************************************
db83d72f	260	if(x[2]<0.){
9251fceb	261	if(TMath::Abs(x[2]+kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
9251fceb	262	b[0] = fCCorrField*fFactorDip;
db83d72f	263	}
9251fceb	264	else if(TMath::Abs(x[2]+kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
db83d72f	265	b[0] = fQuadGradient*x[1];
db83d72f	266	b[1] = fQuadGradient*x[0];
db83d72f	267	}
9251fceb	268	else if(TMath::Abs(x[2]+kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
db83d72f	269	b[0] = -fQuadGradient*x[1];
db83d72f	270	b[1] = -fQuadGradient*x[0];
db83d72f	271	}
9251fceb	272	else if(TMath::Abs(x[2]+kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
db83d72f	273	b[0] = -fQuadGradient*x[1];
db83d72f	274	b[1] = -fQuadGradient*x[0];
db83d72f	275	}
9251fceb	276	else if(TMath::Abs(x[2]+kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
db83d72f	277	b[0] = fQuadGradient*x[1];
db83d72f	278	b[1] = fQuadGradient*x[0];
db83d72f	279	}
9251fceb	280	else if(TMath::Abs(x[2]+kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRC){
db83d72f	281	b[1] = fDipoleField;
db83d72f	282	}
9251fceb	283	else if(TMath::Abs(x[2]+kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRC) {
	284	double dxabs = TMath::Abs(x[0])-kDip2DXC;
	285	if ( (dxabsdxabs + x[1]x[1])<kDip2SqRC) {
db83d72f	286	b[1] = -fDipoleField;
db83d72f	287	}
	288	}
	289	}
	290	//
	291	// SIDE A **************************************************
	292	else{
9251fceb	293	if(TMath::Abs(x[2]-kBComp1CZ)<kBComp1hDZ && rad2 < kBComp1SqR) {
db83d72f	294	// Compensator magnet at z = 1075 m
9251fceb	295	b[0] = fACorr1Field*fFactorDip;
db83d72f	296	}
db83d72f	297	//
9251fceb	298	if(TMath::Abs(x[2]-kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
	299	b[0] = fACorr2Field*fFactorDip;
	300	}
	301	else if(TMath::Abs(x[2]-kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
db83d72f	302	b[0] = -fQuadGradient*x[1];
db83d72f	303	b[1] = -fQuadGradient*x[0];
eed8a1a2	304	}
9251fceb	305	else if(TMath::Abs(x[2]-kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
	306	b[0] = fQuadGradient*x[1];
	307	b[1] = fQuadGradient*x[0];
eed8a1a2	308	}
9251fceb	309	else if(TMath::Abs(x[2]-kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
	310	b[0] = fQuadGradient*x[1];
	311	b[1] = fQuadGradient*x[0];
db83d72f	312	}
9251fceb	313	else if(TMath::Abs(x[2]-kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
db83d72f	314	b[0] = -fQuadGradient*x[1];
db83d72f	315	b[1] = -fQuadGradient*x[0];
db83d72f	316	}
9251fceb	317	else if(TMath::Abs(x[2]-kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRA){
db83d72f	318	b[1] = -fDipoleField;
db83d72f	319	}
9251fceb	320	else if(TMath::Abs(x[2]-kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRA) {
	321	double dxabs = TMath::Abs(x[0])-kDip2DXA;
	322	if ( (dxabsdxabs + x[1]x[1])<kDip2SqRA) {
db83d72f	323	b[1] = fDipoleField;
	324	}
	325	}
	326	}
9251fceb	327	//
db83d72f	328	}
	329
	330	//_______________________________________________________________________
	331	void AliMagF::GetTPCInt(const Double_t xyz, Double_t b) const
	332	{
	333	// Method to calculate the integral of magnetic integral from xyz to nearest cathode plane
	334	b[0]=b[1]=b[2]=0.0;
	335	if (fMeasuredMap) {
	336	fMeasuredMap->GetTPCInt(xyz,b);
	337	for (int i=3;i--;) b[i] *= fFactorSol;
	338	}
	339	}
	340
	341	//_______________________________________________________________________
	342	void AliMagF::GetTPCIntCyl(const Double_t rphiz, Double_t b) const
	343	{
	344	// Method to calculate the integral of magnetic integral from point to nearest cathode plane
	345	// in cylindrical coordiates ( -pi<phi<pi convention )
	346	b[0]=b[1]=b[2]=0.0;
	347	if (fMeasuredMap) {
	348	fMeasuredMap->GetTPCIntCyl(rphiz,b);
	349	for (int i=3;i--;) b[i] *= fFactorSol;
	350	}
eed8a1a2	351	}