[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.;
const Double_t AliMagF::fgkBMachineZ1 =   919.;
const Double_t AliMagF::fgkBMachineZ2 = -1972.;

//_______________________________________________________________________
AliMagF::AliMagF():
  TVirtualMagField(),
  fMeasuredMap(0),
  fMapType(k5kG),
  fSolenoid(0),
  fBeamType(kNoBeamField),
  fBeamEnergy(0),
  fCompensator(kFALSE),
  //
  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, Bool_t compensator):
  TVirtualMagField(name),
  fMeasuredMap(0),
  fMapType(maptype),
  fSolenoid(0),
  fBeamType(bt),
  fBeamEnergy(be),
  fCompensator(compensator),
  //
  fInteg(integ),
  fPrecInteg(1),
  fFactorSol(1.),
  fFactorDip(1.),
  fMax(fmax),
  fDipoleOFF(factorDip==0.),
  //
  fQuadGradient(0),
  fDipoleField(0),
  fCCorrField(0), 
  fACorr1Field(0),
  fACorr2Field(0),
  fParNames("","")
{
  //
  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),
  fCompensator(src.fCompensator),
  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 (xyz[2] > fgkBMachineZ1 || xyz[2] < fgkBMachineZ2) MachineField(xyz, b);
  else if (fMeasuredMap) {
    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;
  }
  //
}

//_______________________________________________________________________
Double_t AliMagF::GetBz(const Double_t *xyz) const
{
  // Method to calculate the field at point  xyz
  //
  if (xyz[2] <= fgkBMachineZ1 && xyz[2] >= fgkBMachineZ2) return fMeasuredMap->GetBz(xyz);
  else {
    double b[3] = {0,0,0};
    MachineField(xyz, b);
    return b[2];
  }
  //
}

//_______________________________________________________________________
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;
    fCompensator = src.fCompensator;
    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)
{
  const double kToler = 0.1;
  if (btype==kNoBeamField) {
    fQuadGradient = fDipoleField = fCCorrField = fACorr1Field = fACorr2Field = 0.;
  }
  //
  else if (btype==kBeamTypepp && TMath::Abs(1.-benergy/5000.)<kToler ){
    // p-p @ 5+5 TeV
    fQuadGradient = 15.7145;
    fDipoleField  = 27.0558;
    // SIDE C
    fCCorrField   = 9.7017;
    // SIDE A
    fACorr1Field  = -13.2143;
    fACorr2Field  = -11.9909;
  } else if (btype == kBeamTypepp && TMath::Abs(1.-benergy/450.)<kToler) {
    // p-p 0.45+0.45 TeV
    Double_t const kEnergyRatio = benergy / 7000.;
    fQuadGradient = 22.0002 * kEnergyRatio;
    fDipoleField  = 37.8781 * kEnergyRatio;
    // SIDE C
    fCCorrField   =  9.6908;
    // SIDE A
    fACorr1Field  = -13.2014;
    fACorr2Field  = -9.6908;
  } else if ( (btype == kBeamTypepp && TMath::Abs(1.-benergy/7000.)<kToler) || 
	      (fBeamType == kBeamTypeAA) ) {
    // Pb-Pb @ 2.7+2.7 TeV or p-p @ 7+7 TeV
    fQuadGradient = 22.0002;
    fDipoleField  = 37.8781;
    // SIDE C
    fCCorrField   = 9.6908;
    // SIDE A
    fACorr1Field  = -13.2014;
    fACorr2Field  = -9.6908;
  }
  //
}

//_______________________________________________________________________
void AliMagF::MachineField(const Double_t *x, Double_t *b) const
{
  // ---- This is the ZDC part
  const Double_t kCCorrBegin = fgkBMachineZ2-0.5,kCCorrEnd = kCCorrBegin - 153., kCCorrSqRadius = 4.5*4.5;
  //
  const Double_t kCTripletBegin  = -2296.5;
  const Double_t kCQ1Begin = kCTripletBegin,        kCQ1End = kCQ1Begin-637., kCQ1SqRadius = 3.5*3.5;
  const Double_t kCQ2Begin = kCTripletBegin-908.5,  kCQ2End = kCQ2Begin-550., kCQ2SqRadius = 3.5*3.5;
  const Double_t kCQ3Begin = kCTripletBegin-1558.5, kCQ3End = kCQ3Begin-550., kCQ3SqRadius = 3.5*3.5;
  const Double_t kCQ4Begin = kCTripletBegin-2400.,  kCQ4End = kCQ4Begin-637., kCQ4SqRadius = 3.5*3.5;
  //
  const Double_t kCD1Begin = -5838.3,  kCD1End = kCD1Begin-945., kCD1SqRadius = 4.5*4.5;
  const Double_t kCD2Begin = -12167.8, kCD2End = kCD2Begin-945., kCD2SqRadius = 4.5*4.5;
  const Double_t kCD2XCentre1 = -9.7;
  const Double_t kCD2XCentre2 =  9.7;
  //
  // -> SIDE A
  // NB -> kACorr1Begin = 919. to be checked
  const Double_t kACorr1Begin = fgkBMachineZ1, kACorr1End = kACorr1Begin+260., kCCorr1SqRadius = 4.*4.;
  const Double_t kACorr2Begin = -fgkBMachineZ2 + 0.5, kACorr2End = kACorr2Begin+153., kCCorr2SqRadius = 4.5*4.5;
  const Double_t kATripletBegin  = 2296.5;
  const Double_t kAQ1Begin = kATripletBegin,	kAQ1End = kAQ1Begin+637., kAQ1SqRadius = 3.5*3.5;
  const Double_t kAQ2Begin = kATripletBegin+908.5,  kAQ2End = kAQ2Begin+550., kAQ2SqRadius = 3.5*3.5;
  const Double_t kAQ3Begin = kATripletBegin+1558.5, kAQ3End = kAQ3Begin+550., kAQ3SqRadius = 3.5*3.5;
  const Double_t kAQ4Begin = kATripletBegin+2400.,  kAQ4End = kAQ4Begin+637., kAQ4SqRadius = 3.5*3.5;
  //
  const Double_t kAD1Begin = 5838.3,  kAD1End = kAD1Begin+945., kAD1SqRadius = 3.375*3.375;
  const Double_t kAD2Begin = 12167.8, kAD2End = kAD2Begin+945., kAD2SqRadius = 3.75*3.75;
  const Double_t kAD2XCentre1 = -9.4;
  const Double_t kAD2XCentre2 =  9.4;
  //
  double rad2 = x[0] * x[0] + x[1] * x[1];
  //
  // SIDE C **************************************************
  if(x[2]<0.){  
    if(x[2] < kCCorrBegin && x[2] > kCCorrEnd && rad2 < kCCorrSqRadius){
      b[0] = fCCorrField;
      b[1] = 0.;
      b[2] = 0.;
    } 
    else if(x[2] < kCQ1Begin && x[2] > kCQ1End && rad2 < kCQ1SqRadius){
      b[0] = fQuadGradient*x[1];
      b[1] = fQuadGradient*x[0];
      b[2] = 0.;
    }
    else if(x[2] < kCQ2Begin && x[2] > kCQ2End && rad2 < kCQ2SqRadius){
      b[0] = -fQuadGradient*x[1];
      b[1] = -fQuadGradient*x[0];
      b[2] = 0.;
    }
    else if(x[2] < kCQ3Begin && x[2] > kCQ3End && rad2 < kCQ3SqRadius){
      b[0] = -fQuadGradient*x[1];
      b[1] = -fQuadGradient*x[0];
      b[2] = 0.;
    }
    else if(x[2] < kCQ4Begin && x[2] > kCQ4End && rad2 < kCQ4SqRadius){
      b[0] = fQuadGradient*x[1];
      b[1] = fQuadGradient*x[0];
      b[2] = 0.;
    }
    else if(x[2] < kCD1Begin && x[2] > kCD1End && rad2 < kCD1SqRadius){
      b[1] = fDipoleField;
      b[2] = 0.;
      b[2] = 0.;
    }
    else if(x[2] < kCD2Begin && x[2] > kCD2End){
      if(((x[0]-kCD2XCentre1)*(x[0]-kCD2XCentre1)+(x[1]*x[1]))<kCD2SqRadius
	 || ((x[0]-kCD2XCentre2)*(x[0]-kCD2XCentre2)+(x[1]*x[1]))<kCD2SqRadius){
	b[1] = -fDipoleField;
	b[2] = 0.;
	b[2] = 0.;
      }
    }
  }
  //
  // SIDE A **************************************************
  else{        
    if(fCompensator && (x[2] > kACorr1Begin && x[2] < kACorr1End) && rad2 < kCCorr1SqRadius) {
      // Compensator magnet at z = 1075 m 
      b[0] = fACorr1Field;
      b[1] = 0.;
      b[2] = 0.;
      return;
    }
    //
    if(x[2] > kACorr2Begin && x[2] < kACorr2End && rad2 < kCCorr2SqRadius){
      b[0] = fACorr2Field;
      b[1] = 0.;
      b[2] = 0.;
    }          
    else if(x[2] > kAQ1Begin && x[2] < kAQ1End && rad2 < kAQ1SqRadius){
      // First quadrupole of inner triplet de-focussing in x-direction
      b[0] = -fQuadGradient*x[1];
      b[1] = -fQuadGradient*x[0];
      b[2] = 0.;
    }
    else if(x[2] > kAQ2Begin && x[2] < kAQ2End && rad2 < kAQ2SqRadius){
      b[0] = fQuadGradient*x[1];
      b[1] = fQuadGradient*x[0];
      b[2] = 0.;
    }
    else if(x[2] > kAQ3Begin && x[2] < kAQ3End && rad2 < kAQ3SqRadius){
      b[0] = fQuadGradient*x[1];
      b[1] = fQuadGradient*x[0];
      b[2] = 0.;
    }
    else if(x[2] > kAQ4Begin && x[2] < kAQ4End && rad2 < kAQ4SqRadius){
      b[0] = -fQuadGradient*x[1];
      b[1] = -fQuadGradient*x[0];
      b[2] = 0.;
    }
    else if(x[2] > kAD1Begin && x[2] < kAD1End && rad2 < kAD1SqRadius){
      b[0] = 0.;
      b[1] = -fDipoleField;
      b[2] = 0.;
    }
    else if(x[2] > kAD2Begin && x[2] < kAD2End){
      if(((x[0]-kAD2XCentre1)*(x[0]-kAD2XCentre1)+(x[1]*x[1])) < kAD2SqRadius
	 || ((x[0]-kAD2XCentre2)*(x[0]-kAD2XCentre2)+(x[1]*x[1])) < kAD2SqRadius){
	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
db83d72f	27	const Double_t AliMagF::fgkSol2DipZ = -700.;
	28	const Double_t AliMagF::fgkBMachineZ1 = 919.;
	29	const Double_t AliMagF::fgkBMachineZ2 = -1972.;
	30
e2afb3b6	31	//_______________________________________________________________________
e2afb3b6	32	AliMagF::AliMagF():
db83d72f	33	TVirtualMagField(),
	34	fMeasuredMap(0),
	35	fMapType(k5kG),
	36	fSolenoid(0),
	37	fBeamType(kNoBeamField),
	38	fBeamEnergy(0),
	39	fCompensator(kFALSE),
	40	//
e2afb3b6	41	fInteg(0),
db83d72f	42	fPrecInteg(0),
	43	fFactorSol(1.),
	44	fFactorDip(1.),
	45	fMax(15),
	46	fDipoleOFF(kFALSE),
e2afb3b6	47	//
db83d72f	48	fQuadGradient(0),
	49	fDipoleField(0),
	50	fCCorrField(0),
	51	fACorr1Field(0),
	52	fACorr2Field(0),
	53	fParNames("","")
	54	{
e2afb3b6	55	// Default constructor
	56	//
	57	}
	58
	59	//_______________________________________________________________________
db83d72f	60	AliMagF::AliMagF(const char name, const char title, Int_t integ,
	61	Double_t factorSol, Double_t factorDip,
	62	Double_t fmax, BMap_t maptype, const char* path,
	63	BeamType_t bt, Double_t be, Bool_t compensator):
	64	TVirtualMagField(name),
	65	fMeasuredMap(0),
	66	fMapType(maptype),
	67	fSolenoid(0),
	68	fBeamType(bt),
	69	fBeamEnergy(be),
	70	fCompensator(compensator),
	71	//
	72	fInteg(integ),
604e0531	73	fPrecInteg(1),
db83d72f	74	fFactorSol(1.),
db83d72f	75	fFactorDip(1.),
972ca52f	76	fMax(fmax),
db83d72f	77	fDipoleOFF(factorDip==0.),
	78	//
	79	fQuadGradient(0),
	80	fDipoleField(0),
	81	fCCorrField(0),
	82	fACorr1Field(0),
	83	fACorr2Field(0),
	84	fParNames("","")
fe4da5cc	85	{
aee8290b	86	//
db83d72f	87	SetTitle(title);
	88	if(integ<0 \|\| integ > 2) {
	89	AliWarning(Form("Invalid magnetic field flag: %5d; Helix tracking chosen instead",integ));
	90	fInteg = 2;
	91	}
	92	if (fInteg == 0) fPrecInteg = 0;
aee8290b	93	//
db83d72f	94	const char* parname = 0;
	95	//
	96	if (fMapType == k2kG) {
	97	fSolenoid = 2.;
	98	parname = fDipoleOFF ? "Sol12_Dip0_Hole":"Sol12_Dip6_Hole";
	99	} else if (fMapType == k5kG) {
	100	fSolenoid = 5.;
	101	parname = fDipoleOFF ? "Sol30_Dip0_Hole":"Sol30_Dip6_Hole";
	102	} else if (fMapType == k5kGUniform) {
	103	fSolenoid = 5.;
	104	parname = "Sol30_Dip6_Uniform";
	105	} else {
	106	AliFatal(Form("Unknown field identifier %d is requested\n",fMapType));
	107	}
	108	//
	109	SetDataFileName(path);
	110	SetParamName(parname);
	111	//
	112	SetFactorSol(factorSol);
	113	SetFactorDip(factorDip);
	114	LoadParameterization();
	115	InitMachineField(fBeamType,fBeamEnergy);
fe4da5cc	116	}
fe4da5cc	117
eeda4611	118	//_______________________________________________________________________
eeda4611	119	AliMagF::AliMagF(const AliMagF &src):
db83d72f	120	TVirtualMagField(src),
	121	fMeasuredMap(0),
	122	fMapType(src.fMapType),
	123	fSolenoid(src.fSolenoid),
	124	fBeamType(src.fBeamType),
	125	fBeamEnergy(src.fBeamEnergy),
	126	fCompensator(src.fCompensator),
eeda4611	127	fInteg(src.fInteg),
eeda4611	128	fPrecInteg(src.fPrecInteg),
db83d72f	129	fFactorSol(src.fFactorSol),
db83d72f	130	fFactorDip(src.fFactorDip),
eeda4611	131	fMax(src.fMax),
db83d72f	132	fDipoleOFF(src.fDipoleOFF),
	133	fQuadGradient(src.fQuadGradient),
	134	fDipoleField(src.fDipoleField),
	135	fCCorrField(src.fCCorrField),
	136	fACorr1Field(src.fACorr1Field),
	137	fACorr2Field(src.fACorr2Field),
	138	fParNames(src.fParNames)
eeda4611	139	{
db83d72f	140	if (src.fMeasuredMap) fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
eeda4611	141	}
eeda4611	142
e2afb3b6	143	//_______________________________________________________________________
db83d72f	144	AliMagF::~AliMagF()
ff66b122	145	{
db83d72f	146	delete fMeasuredMap;
	147	}
	148
	149	//_______________________________________________________________________
	150	Bool_t AliMagF::LoadParameterization()
	151	{
	152	if (fMeasuredMap) {
	153	AliError(Form("Field data %s are already loaded from %s\n",GetParamName(),GetDataFileName()));
	154	return kTRUE;
	155	}
ff66b122	156	//
db83d72f	157	char* fname = gSystem->ExpandPathName(GetDataFileName());
	158	TFile* file = TFile::Open(fname);
	159	if (!file) {
	160	AliError(Form("Failed to open magnetic field data file %s\n",fname));
	161	return kFALSE;
	162	}
ff66b122	163	//
db83d72f	164	fMeasuredMap = dynamic_cast<AliMagWrapCheb*>(file->Get(GetParamName()));
	165	if (!fMeasuredMap) {
	166	AliError(Form("Did not find field %s in %s\n",GetParamName(),fname));
	167	return kFALSE;
	168	}
	169	file->Close();
	170	delete file;
	171	return kTRUE;
ff66b122	172	}
ff66b122	173
db83d72f	174
ff66b122	175	//_______________________________________________________________________
db83d72f	176	void AliMagF::Field(const Double_t xyz, Double_t b)
fe4da5cc	177	{
db83d72f	178	// Method to calculate the field at point xyz
aee8290b	179	//
db83d72f	180	b[0]=b[1]=b[2]=0.0;
	181	if (xyz[2] > fgkBMachineZ1 \|\| xyz[2] < fgkBMachineZ2) MachineField(xyz, b);
	182	else if (fMeasuredMap) {
	183	fMeasuredMap->Field(xyz,b);
	184	if (xyz[2]>fgkSol2DipZ \|\| fDipoleOFF) for (int i=3;i--;) b[i] *= fFactorSol;
	185	else for (int i=3;i--;) b[i] *= fFactorDip;
	186	}
aee8290b	187	//
fe4da5cc	188	}
eeda4611	189
eeda4611	190	//_______________________________________________________________________
db83d72f	191	Double_t AliMagF::GetBz(const Double_t *xyz) const
eeda4611	192	{
db83d72f	193	// Method to calculate the field at point xyz
	194	//
	195	if (xyz[2] <= fgkBMachineZ1 && xyz[2] >= fgkBMachineZ2) return fMeasuredMap->GetBz(xyz);
	196	else {
	197	double b[3] = {0,0,0};
	198	MachineField(xyz, b);
	199	return b[2];
	200	}
	201	//
eeda4611	202	}
	203
	204	//_______________________________________________________________________
db83d72f	205	AliMagF& AliMagF::operator=(const AliMagF& src)
eeda4611	206	{
db83d72f	207	if (this != &src && src.fMeasuredMap) {
	208	if (fMeasuredMap) delete fMeasuredMap;
	209	fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
	210	SetName(src.GetName());
	211	fSolenoid = src.fSolenoid;
	212	fBeamType = src.fBeamType;
	213	fBeamEnergy = src.fBeamEnergy;
	214	fCompensator = src.fCompensator;
	215	fInteg = src.fInteg;
	216	fPrecInteg = src.fPrecInteg;
	217	fFactorSol = src.fFactorSol;
	218	fFactorDip = src.fFactorDip;
	219	fMax = src.fMax;
	220	fDipoleOFF = src.fDipoleOFF;
	221	fParNames = src.fParNames;
	222	}
	223	return *this;
eeda4611	224	}
	225
	226	//_______________________________________________________________________
db83d72f	227	void AliMagF::InitMachineField(BeamType_t btype, Double_t benergy)
eeda4611	228	{
db83d72f	229	const double kToler = 0.1;
	230	if (btype==kNoBeamField) {
	231	fQuadGradient = fDipoleField = fCCorrField = fACorr1Field = fACorr2Field = 0.;
	232	}
	233	//
	234	else if (btype==kBeamTypepp && TMath::Abs(1.-benergy/5000.)<kToler ){
	235	// p-p @ 5+5 TeV
	236	fQuadGradient = 15.7145;
	237	fDipoleField = 27.0558;
	238	// SIDE C
	239	fCCorrField = 9.7017;
	240	// SIDE A
	241	fACorr1Field = -13.2143;
	242	fACorr2Field = -11.9909;
	243	} else if (btype == kBeamTypepp && TMath::Abs(1.-benergy/450.)<kToler) {
	244	// p-p 0.45+0.45 TeV
	245	Double_t const kEnergyRatio = benergy / 7000.;
	246	fQuadGradient = 22.0002 * kEnergyRatio;
	247	fDipoleField = 37.8781 * kEnergyRatio;
	248	// SIDE C
	249	fCCorrField = 9.6908;
	250	// SIDE A
	251	fACorr1Field = -13.2014;
	252	fACorr2Field = -9.6908;
	253	} else if ( (btype == kBeamTypepp && TMath::Abs(1.-benergy/7000.)<kToler) \|\|
	254	(fBeamType == kBeamTypeAA) ) {
	255	// Pb-Pb @ 2.7+2.7 TeV or p-p @ 7+7 TeV
	256	fQuadGradient = 22.0002;
	257	fDipoleField = 37.8781;
	258	// SIDE C
	259	fCCorrField = 9.6908;
	260	// SIDE A
	261	fACorr1Field = -13.2014;
	262	fACorr2Field = -9.6908;
	263	}
	264	//
eeda4611	265	}
eed8a1a2	266
db83d72f	267	//_______________________________________________________________________
db83d72f	268	void AliMagF::MachineField(const Double_t x, Double_t b) const
eed8a1a2	269	{
db83d72f	270	// ---- This is the ZDC part
	271	const Double_t kCCorrBegin = fgkBMachineZ2-0.5,kCCorrEnd = kCCorrBegin - 153., kCCorrSqRadius = 4.5*4.5;
	272	//
	273	const Double_t kCTripletBegin = -2296.5;
	274	const Double_t kCQ1Begin = kCTripletBegin, kCQ1End = kCQ1Begin-637., kCQ1SqRadius = 3.5*3.5;
	275	const Double_t kCQ2Begin = kCTripletBegin-908.5, kCQ2End = kCQ2Begin-550., kCQ2SqRadius = 3.5*3.5;
	276	const Double_t kCQ3Begin = kCTripletBegin-1558.5, kCQ3End = kCQ3Begin-550., kCQ3SqRadius = 3.5*3.5;
	277	const Double_t kCQ4Begin = kCTripletBegin-2400., kCQ4End = kCQ4Begin-637., kCQ4SqRadius = 3.5*3.5;
	278	//
	279	const Double_t kCD1Begin = -5838.3, kCD1End = kCD1Begin-945., kCD1SqRadius = 4.5*4.5;
	280	const Double_t kCD2Begin = -12167.8, kCD2End = kCD2Begin-945., kCD2SqRadius = 4.5*4.5;
	281	const Double_t kCD2XCentre1 = -9.7;
	282	const Double_t kCD2XCentre2 = 9.7;
	283	//
	284	// -> SIDE A
	285	// NB -> kACorr1Begin = 919. to be checked
	286	const Double_t kACorr1Begin = fgkBMachineZ1, kACorr1End = kACorr1Begin+260., kCCorr1SqRadius = 4.*4.;
	287	const Double_t kACorr2Begin = -fgkBMachineZ2 + 0.5, kACorr2End = kACorr2Begin+153., kCCorr2SqRadius = 4.5*4.5;
	288	const Double_t kATripletBegin = 2296.5;
	289	const Double_t kAQ1Begin = kATripletBegin, kAQ1End = kAQ1Begin+637., kAQ1SqRadius = 3.5*3.5;
	290	const Double_t kAQ2Begin = kATripletBegin+908.5, kAQ2End = kAQ2Begin+550., kAQ2SqRadius = 3.5*3.5;
	291	const Double_t kAQ3Begin = kATripletBegin+1558.5, kAQ3End = kAQ3Begin+550., kAQ3SqRadius = 3.5*3.5;
	292	const Double_t kAQ4Begin = kATripletBegin+2400., kAQ4End = kAQ4Begin+637., kAQ4SqRadius = 3.5*3.5;
	293	//
	294	const Double_t kAD1Begin = 5838.3, kAD1End = kAD1Begin+945., kAD1SqRadius = 3.375*3.375;
	295	const Double_t kAD2Begin = 12167.8, kAD2End = kAD2Begin+945., kAD2SqRadius = 3.75*3.75;
	296	const Double_t kAD2XCentre1 = -9.4;
	297	const Double_t kAD2XCentre2 = 9.4;
	298	//
	299	double rad2 = x[0] * x[0] + x[1] * x[1];
	300	//
	301	// SIDE C **************************************************
	302	if(x[2]<0.){
	303	if(x[2] < kCCorrBegin && x[2] > kCCorrEnd && rad2 < kCCorrSqRadius){
	304	b[0] = fCCorrField;
	305	b[1] = 0.;
	306	b[2] = 0.;
	307	}
	308	else if(x[2] < kCQ1Begin && x[2] > kCQ1End && rad2 < kCQ1SqRadius){
	309	b[0] = fQuadGradient*x[1];
	310	b[1] = fQuadGradient*x[0];
	311	b[2] = 0.;
	312	}
	313	else if(x[2] < kCQ2Begin && x[2] > kCQ2End && rad2 < kCQ2SqRadius){
	314	b[0] = -fQuadGradient*x[1];
	315	b[1] = -fQuadGradient*x[0];
	316	b[2] = 0.;
	317	}
	318	else if(x[2] < kCQ3Begin && x[2] > kCQ3End && rad2 < kCQ3SqRadius){
	319	b[0] = -fQuadGradient*x[1];
	320	b[1] = -fQuadGradient*x[0];
	321	b[2] = 0.;
	322	}
	323	else if(x[2] < kCQ4Begin && x[2] > kCQ4End && rad2 < kCQ4SqRadius){
	324	b[0] = fQuadGradient*x[1];
	325	b[1] = fQuadGradient*x[0];
	326	b[2] = 0.;
	327	}
	328	else if(x[2] < kCD1Begin && x[2] > kCD1End && rad2 < kCD1SqRadius){
	329	b[1] = fDipoleField;
	330	b[2] = 0.;
	331	b[2] = 0.;
	332	}
	333	else if(x[2] < kCD2Begin && x[2] > kCD2End){
334	if(((x[0]-kCD2XCentre1)(x[0]-kCD2XCentre1)+(x[1]x[1]))<kCD2SqRadius
335	\|\| ((x[0]-kCD2XCentre2)(x[0]-kCD2XCentre2)+(x[1]x[1]))<kCD2SqRadius){
336	b[1] = -fDipoleField;
337	b[2] = 0.;
338	b[2] = 0.;
339	}
340	}
341	}
342	//
343	// SIDE A **************************************************
344	else{
345	if(fCompensator && (x[2] > kACorr1Begin && x[2] < kACorr1End) && rad2 < kCCorr1SqRadius) {
346	// Compensator magnet at z = 1075 m
347	b[0] = fACorr1Field;
348	b[1] = 0.;
349	b[2] = 0.;
350	return;
351	}
352	//
353	if(x[2] > kACorr2Begin && x[2] < kACorr2End && rad2 < kCCorr2SqRadius){
354	b[0] = fACorr2Field;
355	b[1] = 0.;
356	b[2] = 0.;
357	}
358	else if(x[2] > kAQ1Begin && x[2] < kAQ1End && rad2 < kAQ1SqRadius){
359	// First quadrupole of inner triplet de-focussing in x-direction
360	b[0] = -fQuadGradient*x[1];
361	b[1] = -fQuadGradient*x[0];
362	b[2] = 0.;
eed8a1a2	363	}
db83d72f	364	else if(x[2] > kAQ2Begin && x[2] < kAQ2End && rad2 < kAQ2SqRadius){
	365	b[0] = fQuadGradient*x[1];
	366	b[1] = fQuadGradient*x[0];
	367	b[2] = 0.;
eed8a1a2	368	}
db83d72f	369	else if(x[2] > kAQ3Begin && x[2] < kAQ3End && rad2 < kAQ3SqRadius){
	370	b[0] = fQuadGradient*x[1];
	371	b[1] = fQuadGradient*x[0];
	372	b[2] = 0.;
	373	}
	374	else if(x[2] > kAQ4Begin && x[2] < kAQ4End && rad2 < kAQ4SqRadius){
	375	b[0] = -fQuadGradient*x[1];
	376	b[1] = -fQuadGradient*x[0];
	377	b[2] = 0.;
	378	}
	379	else if(x[2] > kAD1Begin && x[2] < kAD1End && rad2 < kAD1SqRadius){
	380	b[0] = 0.;
	381	b[1] = -fDipoleField;
	382	b[2] = 0.;
	383	}
	384	else if(x[2] > kAD2Begin && x[2] < kAD2End){
	385	if(((x[0]-kAD2XCentre1)(x[0]-kAD2XCentre1)+(x[1]x[1])) < kAD2SqRadius
	386	\|\| ((x[0]-kAD2XCentre2)(x[0]-kAD2XCentre2)+(x[1]x[1])) < kAD2SqRadius){
	387	b[1] = fDipoleField;
	388	}
	389	}
	390	}
	391	}
	392
	393	//_______________________________________________________________________
	394	void AliMagF::GetTPCInt(const Double_t xyz, Double_t b) const
	395	{
	396	// Method to calculate the integral of magnetic integral from xyz to nearest cathode plane
	397	b[0]=b[1]=b[2]=0.0;
	398	if (fMeasuredMap) {
	399	fMeasuredMap->GetTPCInt(xyz,b);
	400	for (int i=3;i--;) b[i] *= fFactorSol;
	401	}
	402	}
	403
	404	//_______________________________________________________________________
	405	void AliMagF::GetTPCIntCyl(const Double_t rphiz, Double_t b) const
	406	{
	407	// Method to calculate the integral of magnetic integral from point to nearest cathode plane
	408	// in cylindrical coordiates ( -pi<phi<pi convention )
	409	b[0]=b[1]=b[2]=0.0;
	410	if (fMeasuredMap) {
	411	fMeasuredMap->GetTPCIntCyl(rphiz,b);
	412	for (int i=3;i--;) b[i] *= fFactorSol;
	413	}
eed8a1a2	414	}