[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(Float_t *x, Float_t *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
	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/* $Id$ */
	17
	18	//-----------------------------------------------------------------------
	19	// Class for Alice magnetic field with constant mesh
	20	// Used in the configuration macros (macros/Config.C, etc.)
	21	// Author:
	22	//-----------------------------------------------------------------------
	23
	24	#include "TSystem.h"
	25
	26	#include "TVector.h"
	27
	28	#include "AliLog.h"
	29	#include "AliMagFCM.h"
	30
	31	ClassImp(AliMagFCM)
	32
	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	//_______________________________________________________________________
	59	AliMagFCM::AliMagFCM(const char name, const char title, Int_t integ,
	60	Float_t factor, Float_t fmax):
	61	AliMagFC(name,title,integ,factor,fmax),
	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)
	76	{
	77	//
	78	// Standard constructor
	79	//
	80	fType = kConMesh;
	81	fMap = 2;
	82	SetSolenoidField();
	83
	84	AliDebug(1, Form(
	85	"Constant Mesh Field %s created: map= %d, factor= %f, file= %s",
	86	fName.Data(), fMap, factor,fTitle.Data()));
	87	}
	88
	89	//_______________________________________________________________________
	90	AliMagFCM::AliMagFCM(const AliMagFCM &magf):
	91	AliMagFC(magf),
	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)
	106	{
	107	//
	108	// Copy constructor
	109	//
	110	magf.Copy(*this);
	111	}
	112
	113	//_______________________________________________________________________
	114	void AliMagFCM::Field(Float_t x, Float_t b) const
	115	{
	116	//
	117	// Method to calculate the magnetic field
	118	//
	119	Double_t ratx, raty, ratz, hix, hiy, hiz, ratx1, raty1, ratz1,
	120	bhyhz, bhylz, blyhz, blylz, bhz, blz, xl[3];
	121	const Double_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	//_______________________________________________________________________
	214	void AliMagFCM::ReadField()
	215	{
	216	//
	217	// Method to read the magnetic field map from file
	218	//
	219	FILE *magfile;
	220	Int_t ix, iy, iz, ipx, ipy, ipz;
	221	Float_t bx, by, bz;
	222	char *fname;
	223	AliDebug(1,Form("Reading Magnetic Field %s from file %s",fName.Data(),fTitle.Data()));
	224	fname = gSystem->ExpandPathName(fTitle.Data());
	225	magfile=fopen(fname,"r");
	226	delete [] fname;
	227	if (magfile) {
	228	fscanf(magfile,"%d %d %d %f %f %f %f %f %f",
	229	&fXn, &fYn, &fZn, &fXdel, &fYdel, &fZdel, &fXbeg, &fYbeg, &fZbeg);
	230	AliDebug(2,Form("fXn %d, fYn %d, fZn %d, fXdel %f, fYdel %f, fZdel %f, fXbeg %f, fYbeg %f, fZbeg %f",
	231	fXn, fYn, fZn, fXdel, fYdel, fZdel, fXbeg, fYbeg, fZbeg));
	232	fXdeli=1./fXdel;
	233	fYdeli=1./fYdel;
	234	fZdeli=1./fZdel;
	235	fB = new TVector(3fXnfYn*fZn);
	236	for (iz=0; iz<fZn; iz++) {
	237	ipz=iz3(fXn*fYn);
	238	for (iy=0; iy<fYn; iy++) {
	239	ipy=ipz+iy3fXn;
	240	for (ix=0; ix<fXn; ix++) {
	241	ipx=ipy+ix*3;
	242	fscanf(magfile,"%f %f %f",&bz,&by,&bx);
	243	(*fB)(ipx+2)=bz;
	244	(*fB)(ipx+1)=by;
	245	(*fB)(ipx )=bx;
	246	}
	247	}
	248	}
	249	} else {
	250	AliFatal(Form("File %s not found !",fTitle.Data()));
	251	}
	252	}
	253
	254	//_______________________________________________________________________
	255	void AliMagFCM::Copy(TObject & /* magf */) const
	256	{
	257	//
	258	// Copy *this onto magf -- Not implemented
	259	//
	260	AliFatal("Not implemented!");
	261	}
	262