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

//-----------------------------------------------------------------------
//  Class for Alice magnetic field with constant mesh
//  Used in the configuration macros (macros/Config.C, etc.)
//  Author:
//-----------------------------------------------------------------------

#include "TSystem.h"

#include "TVector.h"

#include "AliLog.h"
#include "AliMagFCM.h"

ClassImp(AliMagFCM)

//_______________________________________________________________________
AliMagFCM::AliMagFCM():
  fXbeg(0),
  fYbeg(0),
  fZbeg(0),
  fXdel(0),
  fYdel(0),
  fZdel(0),
  fSolenoid(0),
  fXdeli(0),
  fYdeli(0),
  fZdeli(0),
  fXn(0),
  fYn(0),
  fZn(0),
  fB(0)
{
  //
  // Standard constructor
  //
  fType = kConMesh;
  fMap  = 2;
  SetSolenoidField();
}

//_______________________________________________________________________
AliMagFCM::AliMagFCM(const char *name, const char *title, Int_t integ, 
                     Float_t factor, Float_t fmax):
  AliMagFC(name,title,integ,factor,fmax),
  fXbeg(0),
  fYbeg(0),
  fZbeg(0),
  fXdel(0),
  fYdel(0),
  fZdel(0),
  fSolenoid(0),
  fXdeli(0),
  fYdeli(0),
  fZdeli(0),
  fXn(0),
  fYn(0),
  fZn(0),
  fB(0)
{
  //
  // Standard constructor
  //
  fType = kConMesh;
  fMap  = 2;
  SetSolenoidField();

  AliDebug(1, Form(
     "Constant Mesh Field %s created: map= %d, factor= %f, file= %s",
	 fName.Data(), fMap, factor,fTitle.Data()));
}

//_______________________________________________________________________
AliMagFCM::AliMagFCM(const AliMagFCM &magf):
  AliMagFC(magf),
  fXbeg(0),
  fYbeg(0),
  fZbeg(0),
  fXdel(0),
  fYdel(0),
  fZdel(0),
  fSolenoid(0),
  fXdeli(0),
  fYdeli(0),
  fZdeli(0),
  fXn(0),
  fYn(0),
  fZn(0),
  fB(0)
{
  //
  // Copy constructor
  //
  magf.Copy(*this);
}

//_______________________________________________________________________
void AliMagFCM::Field(const float *x, float *b) const
{
  //
  // Method to calculate the magnetic field
  //
  Float_t ratx, raty, ratz, hix, hiy, hiz, ratx1, raty1, ratz1, 
    bhyhz, bhylz, blyhz, blylz, bhz, blz, xl[3];
  const Float_t kone=1;
  Int_t ix, iy, iz;
  
  // --- find the position in the grid ---
  
  b[0]=b[1]=b[2]=0;

  
  if(-700 < -x[2] && -x[2] < fZbeg && x[0] * x[0] +(x[1]+30)*(x[1]+30) < 560*560) {
    b[2]=  fSolenoid;
  } else  {
      // The field map used here was calculated in a coordinate system where the muon arm is at z > 0
      // Transfom x -> -x and z -> -z 
      Float_t xm = - x[0];
      Float_t ym =   x[1];
      Float_t zm = - x[2];
 
      Bool_t infield=(fZbeg <= zm  && zm < fZbeg+fZdel*(fZn-1)
		      &&  ( fXbeg <= TMath::Abs(xm) && TMath::Abs(xm) < fXbeg+fXdel*(fXn-1) )
		      &&  ( fYbeg <= TMath::Abs(ym) && TMath::Abs(ym) < fYbeg+fYdel*(fYn-1) ));
    if(infield) {
      xl[0]=TMath::Abs(xm)-fXbeg;
      xl[1]=TMath::Abs(ym)-fYbeg;
      xl[2]=zm-fZbeg;
      
      // --- start with x
      
      hix=xl[0]*fXdeli;
      ratx=hix-int(hix);
      ix=int(hix);
      
      hiy=xl[1]*fYdeli;
      raty=hiy-int(hiy);
      iy=int(hiy);
      
      hiz=xl[2]*fZdeli;
      ratz=hiz-int(hiz);
      iz=int(hiz);
      
      if(fMap==2) {
	// ... simple interpolation
	ratx1=kone-ratx;
	raty1=kone-raty;
	ratz1=kone-ratz;
	bhyhz = Bx(ix  ,iy+1,iz+1)*ratx1+Bx(ix+1,iy+1,iz+1)*ratx;
	bhylz = Bx(ix  ,iy+1,iz  )*ratx1+Bx(ix+1,iy+1,iz  )*ratx;
	blyhz = Bx(ix  ,iy  ,iz+1)*ratx1+Bx(ix+1,iy  ,iz+1)*ratx;
	blylz = Bx(ix  ,iy  ,iz  )*ratx1+Bx(ix+1,iy  ,iz  )*ratx;
	bhz   = blyhz             *raty1+bhyhz             *raty;
	blz   = blylz             *raty1+bhylz             *raty;
	b[0]  = blz               *ratz1+bhz               *ratz;
	//
	bhyhz = By(ix  ,iy+1,iz+1)*ratx1+By(ix+1,iy+1,iz+1)*ratx;
	bhylz = By(ix  ,iy+1,iz  )*ratx1+By(ix+1,iy+1,iz  )*ratx;
	blyhz = By(ix  ,iy  ,iz+1)*ratx1+By(ix+1,iy  ,iz+1)*ratx;
	blylz = By(ix  ,iy  ,iz  )*ratx1+By(ix+1,iy  ,iz  )*ratx;
	bhz   = blyhz             *raty1+bhyhz             *raty;
	blz   = blylz             *raty1+bhylz             *raty;
	b[1]  = blz               *ratz1+bhz               *ratz;
	//
	bhyhz = Bz(ix  ,iy+1,iz+1)*ratx1+Bz(ix+1,iy+1,iz+1)*ratx;
	bhylz = Bz(ix  ,iy+1,iz  )*ratx1+Bz(ix+1,iy+1,iz  )*ratx;
	blyhz = Bz(ix  ,iy  ,iz+1)*ratx1+Bz(ix+1,iy  ,iz+1)*ratx;
	blylz = Bz(ix  ,iy  ,iz  )*ratx1+Bz(ix+1,iy  ,iz  )*ratx;
	bhz   = blyhz             *raty1+bhyhz             *raty;
	blz   = blylz             *raty1+bhylz             *raty;
	b[2]  = blz               *ratz1+bhz               *ratz;
	//printf("ratx,raty,ratz,b[0],b[1],b[2] %f %f %f %f %f %f\n",
	//ratx,raty,ratz,b[0],b[1],b[2]);
	//
	// ... use the dipole symmetry
	if (xm*ym < 0) b[1]=-b[1];
	if (xm<0)      b[2]=-b[2];
	b[0] = -b[0];
	b[2] = -b[2];
	
      } else {
	AliError(Form("Invalid field map for constant mesh %d",fMap));
      }
    } else {
//This is the ZDC part
	ZDCField(x,b);
    }
    
    if(fFactor!=1) {
	b[0]*=fFactor;
	b[1]*=fFactor;
	b[2]*=fFactor;
    }
  }
}

//_______________________________________________________________________
void AliMagFCM::Field(const double *x, double *b) const
{
  //
  // Method to calculate the magnetic field
  //
  Double_t ratx, raty, ratz, hix, hiy, hiz, ratx1, raty1, ratz1, 
    bhyhz, bhylz, blyhz, blylz, bhz, blz, xl[3];
  const Double_t kone=1;
  Int_t ix, iy, iz;
  
  // --- find the position in the grid ---
  
  b[0]=b[1]=b[2]=0;

  
  if(-700 < -x[2] && -x[2] < fZbeg && x[0] * x[0] +(x[1]+30)*(x[1]+30) < 560*560) {
    b[2]=  fSolenoid;
  } else  {
      // The field map used here was calculated in a coordinate system where the muon arm is at z > 0
      // Transfom x -> -x and z -> -z 
      Float_t xm = - x[0];
      Float_t ym =   x[1];
      Float_t zm = - x[2];
 
      Bool_t infield=(fZbeg <= zm  && zm < fZbeg+fZdel*(fZn-1)
		      &&  ( fXbeg <= TMath::Abs(xm) && TMath::Abs(xm) < fXbeg+fXdel*(fXn-1) )
		      &&  ( fYbeg <= TMath::Abs(ym) && TMath::Abs(ym) < fYbeg+fYdel*(fYn-1) ));
    if(infield) {
      xl[0]=TMath::Abs(xm)-fXbeg;
      xl[1]=TMath::Abs(ym)-fYbeg;
      xl[2]=zm-fZbeg;
      
      // --- start with x
      
      hix=xl[0]*fXdeli;
      ratx=hix-int(hix);
      ix=int(hix);
      
      hiy=xl[1]*fYdeli;
      raty=hiy-int(hiy);
      iy=int(hiy);
      
      hiz=xl[2]*fZdeli;
      ratz=hiz-int(hiz);
      iz=int(hiz);
      
      if(fMap==2) {
	// ... simple interpolation
	ratx1=kone-ratx;
	raty1=kone-raty;
	ratz1=kone-ratz;
	bhyhz = Bx(ix  ,iy+1,iz+1)*ratx1+Bx(ix+1,iy+1,iz+1)*ratx;
	bhylz = Bx(ix  ,iy+1,iz  )*ratx1+Bx(ix+1,iy+1,iz  )*ratx;
	blyhz = Bx(ix  ,iy  ,iz+1)*ratx1+Bx(ix+1,iy  ,iz+1)*ratx;
	blylz = Bx(ix  ,iy  ,iz  )*ratx1+Bx(ix+1,iy  ,iz  )*ratx;
	bhz   = blyhz             *raty1+bhyhz             *raty;
	blz   = blylz             *raty1+bhylz             *raty;
	b[0]  = blz               *ratz1+bhz               *ratz;
	//
	bhyhz = By(ix  ,iy+1,iz+1)*ratx1+By(ix+1,iy+1,iz+1)*ratx;
	bhylz = By(ix  ,iy+1,iz  )*ratx1+By(ix+1,iy+1,iz  )*ratx;
	blyhz = By(ix  ,iy  ,iz+1)*ratx1+By(ix+1,iy  ,iz+1)*ratx;
	blylz = By(ix  ,iy  ,iz  )*ratx1+By(ix+1,iy  ,iz  )*ratx;
	bhz   = blyhz             *raty1+bhyhz             *raty;
	blz   = blylz             *raty1+bhylz             *raty;
	b[1]  = blz               *ratz1+bhz               *ratz;
	//
	bhyhz = Bz(ix  ,iy+1,iz+1)*ratx1+Bz(ix+1,iy+1,iz+1)*ratx;
	bhylz = Bz(ix  ,iy+1,iz  )*ratx1+Bz(ix+1,iy+1,iz  )*ratx;
	blyhz = Bz(ix  ,iy  ,iz+1)*ratx1+Bz(ix+1,iy  ,iz+1)*ratx;
	blylz = Bz(ix  ,iy  ,iz  )*ratx1+Bz(ix+1,iy  ,iz  )*ratx;
	bhz   = blyhz             *raty1+bhyhz             *raty;
	blz   = blylz             *raty1+bhylz             *raty;
	b[2]  = blz               *ratz1+bhz               *ratz;
	//printf("ratx,raty,ratz,b[0],b[1],b[2] %f %f %f %f %f %f\n",
	//ratx,raty,ratz,b[0],b[1],b[2]);
	//
	// ... use the dipole symmetry
	if (xm*ym < 0) b[1]=-b[1];
	if (xm<0)      b[2]=-b[2];
	b[0] = -b[0];
	b[2] = -b[2];
	
      } else {
	AliError(Form("Invalid field map for constant mesh %d",fMap));
      }
    } else {
//This is the ZDC part
	ZDCField(x,b);
    }
    
    if(fFactor!=1) {
	b[0]*=fFactor;
	b[1]*=fFactor;
	b[2]*=fFactor;
    }
  }
}

//_______________________________________________________________________
void AliMagFCM::ReadField()
{
  // 
  // Method to read the magnetic field map from file
  //
  FILE *magfile;
  Int_t ix, iy, iz, ipx, ipy, ipz;
  Float_t bx, by, bz;
  char *fname;
  AliDebug(1,Form("Reading Magnetic Field %s from file %s",fName.Data(),fTitle.Data()));
  fname = gSystem->ExpandPathName(fTitle.Data());
  magfile=fopen(fname,"r");
  delete [] fname;
  if (magfile) {
    fscanf(magfile,"%d %d %d %f %f %f %f %f %f",
    	   &fXn, &fYn, &fZn, &fXdel, &fYdel, &fZdel, &fXbeg, &fYbeg, &fZbeg);
    AliDebug(2,Form("fXn %d, fYn %d, fZn %d, fXdel %f, fYdel %f, fZdel %f, fXbeg %f, fYbeg %f, fZbeg %f",
			fXn, fYn, fZn, fXdel, fYdel, fZdel, fXbeg, fYbeg, fZbeg));
    fXdeli=1./fXdel;
    fYdeli=1./fYdel;
    fZdeli=1./fZdel;
    fB = new TVector(3*fXn*fYn*fZn);
    for (iz=0; iz<fZn; iz++) {
      ipz=iz*3*(fXn*fYn);
      for (iy=0; iy<fYn; iy++) {
	ipy=ipz+iy*3*fXn;
	for (ix=0; ix<fXn; ix++) {
	  ipx=ipy+ix*3;
	  fscanf(magfile,"%f %f %f",&bz,&by,&bx);
	  (*fB)(ipx+2)=bz;
	  (*fB)(ipx+1)=by;
	  (*fB)(ipx  )=bx;
	}
      }
    }
  } else { 
    AliFatal(Form("File %s not found !",fTitle.Data()));
  }
}

//_______________________________________________________________________
void AliMagFCM::Copy(TObject & /* magf */) const
{
  //
  // Copy *this onto magf -- Not implemented
  //
  AliFatal("Not implemented!");
}
Commit	Line	Data
aee8290b	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
acd84897	16	/* $Id$ */
e2afb3b6	17
5d8718b8	18	//-----------------------------------------------------------------------
	19	// Class for Alice magnetic field with constant mesh
	20	// Used in the configuration macros (macros/Config.C, etc.)
	21	// Author:
	22	//-----------------------------------------------------------------------
0742d588	23
fb17acd4	24	#include "TSystem.h"
88cb7938	25
65fb704d	26	#include "TVector.h"
aee8290b	27
594d8990	28	#include "AliLog.h"
aee8290b	29	#include "AliMagFCM.h"
aee8290b	30
	31	ClassImp(AliMagFCM)
	32
e2afb3b6	33	//_______________________________________________________________________
	34	AliMagFCM::AliMagFCM():
	35	fXbeg(0),
	36	fYbeg(0),
	37	fZbeg(0),
	38	fXdel(0),
	39	fYdel(0),
	40	fZdel(0),
	41	fSolenoid(0),
	42	fXdeli(0),
	43	fYdeli(0),
	44	fZdeli(0),
	45	fXn(0),
	46	fYn(0),
	47	fZn(0),
	48	fB(0)
	49	{
	50	//
	51	// Standard constructor
	52	//
	53	fType = kConMesh;
	54	fMap = 2;
	55	SetSolenoidField();
	56	}
	57
	58	//_______________________________________________________________________
d0f1ee3b	59	AliMagFCM::AliMagFCM(const char name, const char title, Int_t integ,
d0f1ee3b	60	Float_t factor, Float_t fmax):
57754f18	61	AliMagFC(name,title,integ,factor,fmax),
e2afb3b6	62	fXbeg(0),
	63	fYbeg(0),
	64	fZbeg(0),
	65	fXdel(0),
	66	fYdel(0),
	67	fZdel(0),
	68	fSolenoid(0),
	69	fXdeli(0),
	70	fYdeli(0),
	71	fZdeli(0),
	72	fXn(0),
	73	fYn(0),
	74	fZn(0),
	75	fB(0)
aee8290b	76	{
	77	//
	78	// Standard constructor
	79	//
	80	fType = kConMesh;
d8408e76	81	fMap = 2;
4cc8933f	82	SetSolenoidField();
4cc8933f	83
594d8990	84	AliDebug(1, Form(
	85	"Constant Mesh Field %s created: map= %d, factor= %f, file= %s",
	86	fName.Data(), fMap, factor,fTitle.Data()));
aee8290b	87	}
aee8290b	88
e2afb3b6	89	//_______________________________________________________________________
e2afb3b6	90	AliMagFCM::AliMagFCM(const AliMagFCM &magf):
57754f18	91	AliMagFC(magf),
e2afb3b6	92	fXbeg(0),
	93	fYbeg(0),
	94	fZbeg(0),
	95	fXdel(0),
	96	fYdel(0),
	97	fZdel(0),
	98	fSolenoid(0),
	99	fXdeli(0),
	100	fYdeli(0),
	101	fZdeli(0),
	102	fXn(0),
	103	fYn(0),
	104	fZn(0),
	105	fB(0)
aee8290b	106	{
	107	//
	108	// Copy constructor
	109	//
	110	magf.Copy(*this);
	111	}
	112
e2afb3b6	113	//_______________________________________________________________________
611fa94a	114	void AliMagFCM::Field(const float x, float b) const
ff66b122	115	{
	116	//
	117	// Method to calculate the magnetic field
	118	//
	119	Float_t ratx, raty, ratz, hix, hiy, hiz, ratx1, raty1, ratz1,
	120	bhyhz, bhylz, blyhz, blylz, bhz, blz, xl[3];
	121	const Float_t kone=1;
	122	Int_t ix, iy, iz;
	123
	124	// --- find the position in the grid ---
	125
	126	b[0]=b[1]=b[2]=0;
	127
	128
	129	if(-700 < -x[2] && -x[2] < fZbeg && x[0] * x[0] +(x[1]+30)(x[1]+30) < 560560) {
	130	b[2]= fSolenoid;
	131	} else {
	132	// The field map used here was calculated in a coordinate system where the muon arm is at z > 0
	133	// Transfom x -> -x and z -> -z
	134	Float_t xm = - x[0];
	135	Float_t ym = x[1];
	136	Float_t zm = - x[2];
	137
	138	Bool_t infield=(fZbeg <= zm && zm < fZbeg+fZdel*(fZn-1)
	139	&& ( fXbeg <= TMath::Abs(xm) && TMath::Abs(xm) < fXbeg+fXdel*(fXn-1) )
	140	&& ( fYbeg <= TMath::Abs(ym) && TMath::Abs(ym) < fYbeg+fYdel*(fYn-1) ));
	141	if(infield) {
	142	xl[0]=TMath::Abs(xm)-fXbeg;
	143	xl[1]=TMath::Abs(ym)-fYbeg;
	144	xl[2]=zm-fZbeg;
	145
	146	// --- start with x
	147
	148	hix=xl[0]*fXdeli;
	149	ratx=hix-int(hix);
	150	ix=int(hix);
	151
	152	hiy=xl[1]*fYdeli;
	153	raty=hiy-int(hiy);
	154	iy=int(hiy);
	155
	156	hiz=xl[2]*fZdeli;
	157	ratz=hiz-int(hiz);
	158	iz=int(hiz);
	159
	160	if(fMap==2) {
	161	// ... simple interpolation
	162	ratx1=kone-ratx;
	163	raty1=kone-raty;
	164	ratz1=kone-ratz;
	165	bhyhz = Bx(ix ,iy+1,iz+1)ratx1+Bx(ix+1,iy+1,iz+1)ratx;
	166	bhylz = Bx(ix ,iy+1,iz )ratx1+Bx(ix+1,iy+1,iz )ratx;
	167	blyhz = Bx(ix ,iy ,iz+1)ratx1+Bx(ix+1,iy ,iz+1)ratx;
	168	blylz = Bx(ix ,iy ,iz )ratx1+Bx(ix+1,iy ,iz )ratx;
	169	bhz = blyhz raty1+bhyhz raty;
	170	blz = blylz raty1+bhylz raty;
	171	b[0] = blz ratz1+bhz ratz;
	172	//
	173	bhyhz = By(ix ,iy+1,iz+1)ratx1+By(ix+1,iy+1,iz+1)ratx;
	174	bhylz = By(ix ,iy+1,iz )ratx1+By(ix+1,iy+1,iz )ratx;
	175	blyhz = By(ix ,iy ,iz+1)ratx1+By(ix+1,iy ,iz+1)ratx;
	176	blylz = By(ix ,iy ,iz )ratx1+By(ix+1,iy ,iz )ratx;
	177	bhz = blyhz raty1+bhyhz raty;
	178	blz = blylz raty1+bhylz raty;
179	b[1] = blz ratz1+bhz ratz;
180	//
181	bhyhz = Bz(ix ,iy+1,iz+1)ratx1+Bz(ix+1,iy+1,iz+1)ratx;
182	bhylz = Bz(ix ,iy+1,iz )ratx1+Bz(ix+1,iy+1,iz )ratx;
183	blyhz = Bz(ix ,iy ,iz+1)ratx1+Bz(ix+1,iy ,iz+1)ratx;
184	blylz = Bz(ix ,iy ,iz )ratx1+Bz(ix+1,iy ,iz )ratx;
185	bhz = blyhz raty1+bhyhz raty;
186	blz = blylz raty1+bhylz raty;
187	b[2] = blz ratz1+bhz ratz;
188	//printf("ratx,raty,ratz,b[0],b[1],b[2] %f %f %f %f %f %f\n",
189	//ratx,raty,ratz,b[0],b[1],b[2]);
190	//
191	// ... use the dipole symmetry
192	if (xm*ym < 0) b[1]=-b[1];
193	if (xm<0) b[2]=-b[2];
194	b[0] = -b[0];
195	b[2] = -b[2];
196
197	} else {
198	AliError(Form("Invalid field map for constant mesh %d",fMap));
199	}
200	} else {
201	//This is the ZDC part
202	ZDCField(x,b);
203	}
204
205	if(fFactor!=1) {
206	b[0]*=fFactor;
207	b[1]*=fFactor;
208	b[2]*=fFactor;
209	}
210	}
211	}
212
213	//_______________________________________________________________________
611fa94a	214	void AliMagFCM::Field(const double x, double b) const
aee8290b	215	{
	216	//
	217	// Method to calculate the magnetic field
	218	//
	219	Double_t ratx, raty, ratz, hix, hiy, hiz, ratx1, raty1, ratz1,
	220	bhyhz, bhylz, blyhz, blylz, bhz, blz, xl[3];
	221	const Double_t kone=1;
	222	Int_t ix, iy, iz;
	223
	224	// --- find the position in the grid ---
	225
	226	b[0]=b[1]=b[2]=0;
57754f18	227
	228
	229	if(-700 < -x[2] && -x[2] < fZbeg && x[0] * x[0] +(x[1]+30)(x[1]+30) < 560560) {
	230	b[2]= fSolenoid;
aee8290b	231	} else {
57754f18	232	// The field map used here was calculated in a coordinate system where the muon arm is at z > 0
	233	// Transfom x -> -x and z -> -z
	234	Float_t xm = - x[0];
	235	Float_t ym = x[1];
	236	Float_t zm = - x[2];
	237
	238	Bool_t infield=(fZbeg <= zm && zm < fZbeg+fZdel*(fZn-1)
	239	&& ( fXbeg <= TMath::Abs(xm) && TMath::Abs(xm) < fXbeg+fXdel*(fXn-1) )
	240	&& ( fYbeg <= TMath::Abs(ym) && TMath::Abs(ym) < fYbeg+fYdel*(fYn-1) ));
aee8290b	241	if(infield) {
57754f18	242	xl[0]=TMath::Abs(xm)-fXbeg;
	243	xl[1]=TMath::Abs(ym)-fYbeg;
	244	xl[2]=zm-fZbeg;
aee8290b	245
	246	// --- start with x
	247
	248	hix=xl[0]*fXdeli;
	249	ratx=hix-int(hix);
	250	ix=int(hix);
	251
	252	hiy=xl[1]*fYdeli;
	253	raty=hiy-int(hiy);
	254	iy=int(hiy);
	255
	256	hiz=xl[2]*fZdeli;
	257	ratz=hiz-int(hiz);
	258	iz=int(hiz);
	259
	260	if(fMap==2) {
	261	// ... simple interpolation
	262	ratx1=kone-ratx;
	263	raty1=kone-raty;
	264	ratz1=kone-ratz;
	265	bhyhz = Bx(ix ,iy+1,iz+1)ratx1+Bx(ix+1,iy+1,iz+1)ratx;
	266	bhylz = Bx(ix ,iy+1,iz )ratx1+Bx(ix+1,iy+1,iz )ratx;
	267	blyhz = Bx(ix ,iy ,iz+1)ratx1+Bx(ix+1,iy ,iz+1)ratx;
	268	blylz = Bx(ix ,iy ,iz )ratx1+Bx(ix+1,iy ,iz )ratx;
	269	bhz = blyhz raty1+bhyhz raty;
	270	blz = blylz raty1+bhylz raty;
	271	b[0] = blz ratz1+bhz ratz;
	272	//
	273	bhyhz = By(ix ,iy+1,iz+1)ratx1+By(ix+1,iy+1,iz+1)ratx;
	274	bhylz = By(ix ,iy+1,iz )ratx1+By(ix+1,iy+1,iz )ratx;
	275	blyhz = By(ix ,iy ,iz+1)ratx1+By(ix+1,iy ,iz+1)ratx;
	276	blylz = By(ix ,iy ,iz )ratx1+By(ix+1,iy ,iz )ratx;
	277	bhz = blyhz raty1+bhyhz raty;
	278	blz = blylz raty1+bhylz raty;
	279	b[1] = blz ratz1+bhz ratz;
	280	//
	281	bhyhz = Bz(ix ,iy+1,iz+1)ratx1+Bz(ix+1,iy+1,iz+1)ratx;
	282	bhylz = Bz(ix ,iy+1,iz )ratx1+Bz(ix+1,iy+1,iz )ratx;
	283	blyhz = Bz(ix ,iy ,iz+1)ratx1+Bz(ix+1,iy ,iz+1)ratx;
	284	blylz = Bz(ix ,iy ,iz )ratx1+Bz(ix+1,iy ,iz )ratx;
	285	bhz = blyhz raty1+bhyhz raty;
	286	blz = blylz raty1+bhylz raty;
	287	b[2] = blz ratz1+bhz ratz;
	288	//printf("ratx,raty,ratz,b[0],b[1],b[2] %f %f %f %f %f %f\n",
	289	//ratx,raty,ratz,b[0],b[1],b[2]);
	290	//
	291	// ... use the dipole symmetry
57754f18	292	if (xm*ym < 0) b[1]=-b[1];
	293	if (xm<0) b[2]=-b[2];
	294	b[0] = -b[0];
	295	b[2] = -b[2];
	296
aee8290b	297	} else {
594d8990	298	AliError(Form("Invalid field map for constant mesh %d",fMap));
aee8290b	299	}
aee8290b	300	} else {
d95ea23f	301	//This is the ZDC part
57754f18	302	ZDCField(x,b);
aee8290b	303	}
d95ea23f	304
57754f18	305	if(fFactor!=1) {
	306	b[0]*=fFactor;
	307	b[1]*=fFactor;
	308	b[2]*=fFactor;
d95ea23f	309	}
aee8290b	310	}
aee8290b	311	}
aee8290b	312
e2afb3b6	313	//_______________________________________________________________________
aee8290b	314	void AliMagFCM::ReadField()
	315	{
	316	//
	317	// Method to read the magnetic field map from file
	318	//
	319	FILE *magfile;
	320	Int_t ix, iy, iz, ipx, ipy, ipz;
	321	Float_t bx, by, bz;
	322	char *fname;
594d8990	323	AliDebug(1,Form("Reading Magnetic Field %s from file %s",fName.Data(),fTitle.Data()));
aee8290b	324	fname = gSystem->ExpandPathName(fTitle.Data());
	325	magfile=fopen(fname,"r");
	326	delete [] fname;
	327	if (magfile) {
	328	fscanf(magfile,"%d %d %d %f %f %f %f %f %f",
	329	&fXn, &fYn, &fZn, &fXdel, &fYdel, &fZdel, &fXbeg, &fYbeg, &fZbeg);
594d8990	330	AliDebug(2,Form("fXn %d, fYn %d, fZn %d, fXdel %f, fYdel %f, fZdel %f, fXbeg %f, fYbeg %f, fZbeg %f",
594d8990	331	fXn, fYn, fZn, fXdel, fYdel, fZdel, fXbeg, fYbeg, fZbeg));
aee8290b	332	fXdeli=1./fXdel;
	333	fYdeli=1./fYdel;
	334	fZdeli=1./fZdel;
	335	fB = new TVector(3fXnfYn*fZn);
	336	for (iz=0; iz<fZn; iz++) {
	337	ipz=iz3(fXn*fYn);
	338	for (iy=0; iy<fYn; iy++) {
	339	ipy=ipz+iy3fXn;
	340	for (ix=0; ix<fXn; ix++) {
	341	ipx=ipy+ix*3;
	342	fscanf(magfile,"%f %f %f",&bz,&by,&bx);
	343	(*fB)(ipx+2)=bz;
	344	(*fB)(ipx+1)=by;
	345	(*fB)(ipx )=bx;
	346	}
	347	}
	348	}
	349	} else {
594d8990	350	AliFatal(Form("File %s not found !",fTitle.Data()));
aee8290b	351	}
	352	}
	353
e2afb3b6	354	//_______________________________________________________________________
6c4904c2	355	void AliMagFCM::Copy(TObject & /* magf */) const
aee8290b	356	{
aee8290b	357	//
8918e700	358	// Copy *this onto magf -- Not implemented
aee8290b	359	//
594d8990	360	AliFatal("Not implemented!");
aee8290b	361	}
aee8290b	362