[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) {
    fSolenoid = 2.;
    parname = fDipoleOFF ? "Sol12_Dip0_Hole":"Sol12_Dip6_Hole";
  } else if (fMapType == k5kG) {
    fSolenoid = 5.;
    parname = fDipoleOFF ? "Sol30_Dip0_Hole":"Sol30_Dip6_Hole";
  } else if (fMapType == k5kGUniform) {
    fSolenoid = 5.;
    parname = "Sol30_Dip6_Uniform";
  } else {
    AliFatal(Form("Unknown field identifier %d is requested\n",fMapType)); 
  }
  //
  SetDataFileName(path);
  SetParamName(parname);
  //
  SetFactorSol(factorSol);
  SetFactorDip(factorDip);
  LoadParameterization();
  InitMachineField(fBeamType,fBeamEnergy);
}

//_______________________________________________________________________
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;
	94	//
	95	if (fMapType == k2kG) {
	96	fSolenoid = 2.;
	97	parname = fDipoleOFF ? "Sol12_Dip0_Hole":"Sol12_Dip6_Hole";
	98	} else if (fMapType == k5kG) {
	99	fSolenoid = 5.;
	100	parname = fDipoleOFF ? "Sol30_Dip0_Hole":"Sol30_Dip6_Hole";
	101	} else if (fMapType == k5kGUniform) {
	102	fSolenoid = 5.;
	103	parname = "Sol30_Dip6_Uniform";
	104	} else {
	105	AliFatal(Form("Unknown field identifier %d is requested\n",fMapType));
	106	}
	107	//
	108	SetDataFileName(path);
	109	SetParamName(parname);
	110	//
	111	SetFactorSol(factorSol);
	112	SetFactorDip(factorDip);
	113	LoadParameterization();
	114	InitMachineField(fBeamType,fBeamEnergy);
fe4da5cc	115	}
fe4da5cc	116
eeda4611	117	//_______________________________________________________________________
eeda4611	118	AliMagF::AliMagF(const AliMagF &src):
db83d72f	119	TVirtualMagField(src),
	120	fMeasuredMap(0),
	121	fMapType(src.fMapType),
	122	fSolenoid(src.fSolenoid),
	123	fBeamType(src.fBeamType),
	124	fBeamEnergy(src.fBeamEnergy),
eeda4611	125	fInteg(src.fInteg),
eeda4611	126	fPrecInteg(src.fPrecInteg),
db83d72f	127	fFactorSol(src.fFactorSol),
db83d72f	128	fFactorDip(src.fFactorDip),
eeda4611	129	fMax(src.fMax),
db83d72f	130	fDipoleOFF(src.fDipoleOFF),
	131	fQuadGradient(src.fQuadGradient),
	132	fDipoleField(src.fDipoleField),
	133	fCCorrField(src.fCCorrField),
	134	fACorr1Field(src.fACorr1Field),
	135	fACorr2Field(src.fACorr2Field),
	136	fParNames(src.fParNames)
eeda4611	137	{
db83d72f	138	if (src.fMeasuredMap) fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
eeda4611	139	}
eeda4611	140
e2afb3b6	141	//_______________________________________________________________________
db83d72f	142	AliMagF::~AliMagF()
ff66b122	143	{
db83d72f	144	delete fMeasuredMap;
	145	}
	146
	147	//_______________________________________________________________________
	148	Bool_t AliMagF::LoadParameterization()
	149	{
	150	if (fMeasuredMap) {
	151	AliError(Form("Field data %s are already loaded from %s\n",GetParamName(),GetDataFileName()));
	152	return kTRUE;
	153	}
ff66b122	154	//
db83d72f	155	char* fname = gSystem->ExpandPathName(GetDataFileName());
	156	TFile* file = TFile::Open(fname);
	157	if (!file) {
	158	AliError(Form("Failed to open magnetic field data file %s\n",fname));
	159	return kFALSE;
	160	}
ff66b122	161	//
db83d72f	162	fMeasuredMap = dynamic_cast<AliMagWrapCheb*>(file->Get(GetParamName()));
	163	if (!fMeasuredMap) {
	164	AliError(Form("Did not find field %s in %s\n",GetParamName(),fname));
	165	return kFALSE;
	166	}
	167	file->Close();
	168	delete file;
	169	return kTRUE;
ff66b122	170	}
ff66b122	171
db83d72f	172
ff66b122	173	//_______________________________________________________________________
db83d72f	174	void AliMagF::Field(const Double_t xyz, Double_t b)
fe4da5cc	175	{
db83d72f	176	// Method to calculate the field at point xyz
aee8290b	177	//
9251fceb	178	// b[0]=b[1]=b[2]=0.0;
9251fceb	179	if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
db83d72f	180	fMeasuredMap->Field(xyz,b);
db83d72f	181	if (xyz[2]>fgkSol2DipZ \|\| fDipoleOFF) for (int i=3;i--;) b[i] *= fFactorSol;
9251fceb	182	else for (int i=3;i--;) b[i] *= fFactorDip;
db83d72f	183	}
9251fceb	184	else MachineField(xyz, b);
aee8290b	185	//
fe4da5cc	186	}
eeda4611	187
eeda4611	188	//_______________________________________________________________________
db83d72f	189	Double_t AliMagF::GetBz(const Double_t *xyz) const
eeda4611	190	{
db83d72f	191	// Method to calculate the field at point xyz
db83d72f	192	//
9251fceb	193	if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
	194	double bz = fMeasuredMap->GetBz(xyz);
	195	return (xyz[2]>fgkSol2DipZ \|\| fDipoleOFF) ? bzfFactorSol : bzfFactorDip;
db83d72f	196	}
9251fceb	197	else return 0.;
eeda4611	198	}
	199
	200	//_______________________________________________________________________
db83d72f	201	AliMagF& AliMagF::operator=(const AliMagF& src)
eeda4611	202	{
db83d72f	203	if (this != &src && src.fMeasuredMap) {
	204	if (fMeasuredMap) delete fMeasuredMap;
	205	fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
	206	SetName(src.GetName());
	207	fSolenoid = src.fSolenoid;
	208	fBeamType = src.fBeamType;
	209	fBeamEnergy = src.fBeamEnergy;
db83d72f	210	fInteg = src.fInteg;
	211	fPrecInteg = src.fPrecInteg;
	212	fFactorSol = src.fFactorSol;
	213	fFactorDip = src.fFactorDip;
	214	fMax = src.fMax;
	215	fDipoleOFF = src.fDipoleOFF;
	216	fParNames = src.fParNames;
	217	}
	218	return *this;
eeda4611	219	}
	220
	221	//_______________________________________________________________________
db83d72f	222	void AliMagF::InitMachineField(BeamType_t btype, Double_t benergy)
eeda4611	223	{
9251fceb	224	if (btype==kNoBeamField \|\| benergy<1.) {
db83d72f	225	fQuadGradient = fDipoleField = fCCorrField = fACorr1Field = fACorr2Field = 0.;
9251fceb	226	return;
db83d72f	227	}
db83d72f	228	//
9251fceb	229	double rigScale = benergy/7000.; // scale according to ratio of E/Enominal
	230	// for ions assume PbPb (with energy provided per nucleon) and account for A/Z
	231	if (btype == kBeamTypeAA) rigScale *= 208./82.;
	232	//
	233	fQuadGradient = 22.0002*rigScale;
	234	fDipoleField = 37.8781*rigScale;
	235	//
	236	// SIDE C
	237	fCCorrField = -9.6980;
	238	// SIDE A
	239	fACorr1Field = -13.2247;
	240	fACorr2Field = 11.7905;
db83d72f	241	//
eeda4611	242	}
eed8a1a2	243
db83d72f	244	//_______________________________________________________________________
db83d72f	245	void AliMagF::MachineField(const Double_t x, Double_t b) const
eed8a1a2	246	{
db83d72f	247	// ---- This is the ZDC part
9251fceb	248	// Compansators for Alice Muon Arm Dipole
	249	const Double_t kBComp1CZ = 1075., kBComp1hDZ = 260./2., kBComp1SqR = 4.0*4.0;
	250	const Double_t kBComp2CZ = 2049., kBComp2hDZ = 153./2., kBComp2SqR = 4.5*4.5;
	251	//
	252	const Double_t kTripQ1CZ = 2615., kTripQ1hDZ = 637./2., kTripQ1SqR = 3.5*3.5;
90ae20c9	253	const Double_t kTripQ2CZ = 3480., kTripQ2hDZ = 550./2., kTripQ2SqR = 3.5*3.5;
9251fceb	254	const Double_t kTripQ3CZ = 4130., kTripQ3hDZ = 550./2., kTripQ3SqR = 3.5*3.5;
9251fceb	255	const Double_t kTripQ4CZ = 5015., kTripQ4hDZ = 637./2., kTripQ4SqR = 3.5*3.5;
db83d72f	256	//
9251fceb	257	const Double_t kDip1CZ = 6310.8, kDip1hDZ = 945./2., kDip1SqRC = 4.54.5, kDip1SqRA = 3.3753.375;
	258	const Double_t kDip2CZ = 12640.3, kDip2hDZ = 945./2., kDip2SqRC = 4.54.5, kDip2SqRA = 3.753.75;
	259	const Double_t kDip2DXC = 9.7, kDip2DXA = 9.4;
db83d72f	260	//
	261	double rad2 = x[0] * x[0] + x[1] * x[1];
	262	//
9251fceb	263	b[0] = b[1] = b[2] = 0;
9251fceb	264	//
db83d72f	265	// SIDE C **************************************************
db83d72f	266	if(x[2]<0.){
9251fceb	267	if(TMath::Abs(x[2]+kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
9251fceb	268	b[0] = fCCorrField*fFactorDip;
db83d72f	269	}
9251fceb	270	else if(TMath::Abs(x[2]+kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
db83d72f	271	b[0] = fQuadGradient*x[1];
db83d72f	272	b[1] = fQuadGradient*x[0];
db83d72f	273	}
9251fceb	274	else if(TMath::Abs(x[2]+kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
db83d72f	275	b[0] = -fQuadGradient*x[1];
db83d72f	276	b[1] = -fQuadGradient*x[0];
db83d72f	277	}
9251fceb	278	else if(TMath::Abs(x[2]+kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
db83d72f	279	b[0] = -fQuadGradient*x[1];
db83d72f	280	b[1] = -fQuadGradient*x[0];
db83d72f	281	}
9251fceb	282	else if(TMath::Abs(x[2]+kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
db83d72f	283	b[0] = fQuadGradient*x[1];
db83d72f	284	b[1] = fQuadGradient*x[0];
db83d72f	285	}
9251fceb	286	else if(TMath::Abs(x[2]+kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRC){
db83d72f	287	b[1] = fDipoleField;
db83d72f	288	}
9251fceb	289	else if(TMath::Abs(x[2]+kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRC) {
	290	double dxabs = TMath::Abs(x[0])-kDip2DXC;
	291	if ( (dxabsdxabs + x[1]x[1])<kDip2SqRC) {
db83d72f	292	b[1] = -fDipoleField;
db83d72f	293	}
	294	}
	295	}
	296	//
	297	// SIDE A **************************************************
	298	else{
9251fceb	299	if(TMath::Abs(x[2]-kBComp1CZ)<kBComp1hDZ && rad2 < kBComp1SqR) {
db83d72f	300	// Compensator magnet at z = 1075 m
9251fceb	301	b[0] = fACorr1Field*fFactorDip;
db83d72f	302	}
db83d72f	303	//
9251fceb	304	if(TMath::Abs(x[2]-kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
	305	b[0] = fACorr2Field*fFactorDip;
	306	}
	307	else if(TMath::Abs(x[2]-kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
db83d72f	308	b[0] = -fQuadGradient*x[1];
db83d72f	309	b[1] = -fQuadGradient*x[0];
eed8a1a2	310	}
9251fceb	311	else if(TMath::Abs(x[2]-kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
	312	b[0] = fQuadGradient*x[1];
	313	b[1] = fQuadGradient*x[0];
eed8a1a2	314	}
9251fceb	315	else if(TMath::Abs(x[2]-kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
	316	b[0] = fQuadGradient*x[1];
	317	b[1] = fQuadGradient*x[0];
db83d72f	318	}
9251fceb	319	else if(TMath::Abs(x[2]-kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
db83d72f	320	b[0] = -fQuadGradient*x[1];
db83d72f	321	b[1] = -fQuadGradient*x[0];
db83d72f	322	}
9251fceb	323	else if(TMath::Abs(x[2]-kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRA){
db83d72f	324	b[1] = -fDipoleField;
db83d72f	325	}
9251fceb	326	else if(TMath::Abs(x[2]-kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRA) {
	327	double dxabs = TMath::Abs(x[0])-kDip2DXA;
	328	if ( (dxabsdxabs + x[1]x[1])<kDip2SqRA) {
db83d72f	329	b[1] = fDipoleField;
	330	}
	331	}
	332	}
9251fceb	333	//
db83d72f	334	}
	335
	336	//_______________________________________________________________________
	337	void AliMagF::GetTPCInt(const Double_t xyz, Double_t b) const
	338	{
	339	// Method to calculate the integral of magnetic integral from xyz to nearest cathode plane
	340	b[0]=b[1]=b[2]=0.0;
	341	if (fMeasuredMap) {
	342	fMeasuredMap->GetTPCInt(xyz,b);
	343	for (int i=3;i--;) b[i] *= fFactorSol;
	344	}
	345	}
	346
	347	//_______________________________________________________________________
	348	void AliMagF::GetTPCIntCyl(const Double_t rphiz, Double_t b) const
	349	{
	350	// Method to calculate the integral of magnetic integral from point to nearest cathode plane
	351	// in cylindrical coordiates ( -pi<phi<pi convention )
	352	b[0]=b[1]=b[2]=0.0;
	353	if (fMeasuredMap) {
	354	fMeasuredMap->GetTPCIntCyl(rphiz,b);
	355	for (int i=3;i--;) b[i] *= fFactorSol;
	356	}
eed8a1a2	357	}