[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 "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();

  if(fDebug>-1) Info("ctor",
     "%s: Constant Mesh Field %s created: map= %d, factor= %f, file= %s\n",
	 ClassName(),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)
{
  //
  // 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 {
	printf("Invalid field map for constant mesh %d\n",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;
  if(fDebug) printf("%s: Reading Magnetic Field %s from file %s\n",ClassName(),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);
    if(fDebug>1) printf("%s: fXn %d, fYn %d, fZn %d, fXdel %f, fYdel %f, fZdel %f, fXbeg %f, fYbeg %f, fZbeg %f\n",
			ClassName(),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 { 
    printf("%s: File %s not found !\n",ClassName(),fTitle.Data());
    exit(1);
  }
}

//_______________________________________________________________________
void AliMagFCM::Copy(TObject & /* magf */) const
{
  //
  // Copy *this onto magf -- Not implemented
  //
  Fatal("Copy","Not implemented!\n");
}
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 "AliMagFCM.h"
	29
	30	ClassImp(AliMagFCM)
	31
	32	//_______________________________________________________________________
	33	AliMagFCM::AliMagFCM():
	34	fXbeg(0),
	35	fYbeg(0),
	36	fZbeg(0),
	37	fXdel(0),
	38	fYdel(0),
	39	fZdel(0),
	40	fSolenoid(0),
	41	fXdeli(0),
	42	fYdeli(0),
	43	fZdeli(0),
	44	fXn(0),
	45	fYn(0),
	46	fZn(0),
	47	fB(0)
	48	{
	49	//
	50	// Standard constructor
	51	//
	52	fType = kConMesh;
	53	fMap = 2;
	54	SetSolenoidField();
	55	}
	56
	57	//_______________________________________________________________________
	58	AliMagFCM::AliMagFCM(const char name, const char title, Int_t integ,
	59	Float_t factor, Float_t fmax):
	60	AliMagFC(name,title,integ,factor,fmax),
	61	fXbeg(0),
	62	fYbeg(0),
	63	fZbeg(0),
	64	fXdel(0),
	65	fYdel(0),
	66	fZdel(0),
	67	fSolenoid(0),
	68	fXdeli(0),
	69	fYdeli(0),
	70	fZdeli(0),
	71	fXn(0),
	72	fYn(0),
	73	fZn(0),
	74	fB(0)
	75	{
	76	//
	77	// Standard constructor
	78	//
	79	fType = kConMesh;
	80	fMap = 2;
	81	SetSolenoidField();
	82
	83	if(fDebug>-1) Info("ctor",
	84	"%s: Constant Mesh Field %s created: map= %d, factor= %f, file= %s\n",
	85	ClassName(),fName.Data(), fMap, factor,fTitle.Data());
	86	}
	87
	88	//_______________________________________________________________________
	89	AliMagFCM::AliMagFCM(const AliMagFCM &magf):
	90	AliMagFC(magf),
	91	fXbeg(0),
	92	fYbeg(0),
	93	fZbeg(0),
	94	fXdel(0),
	95	fYdel(0),
	96	fZdel(0),
	97	fSolenoid(0),
	98	fXdeli(0),
	99	fYdeli(0),
	100	fZdeli(0),
	101	fXn(0),
	102	fYn(0),
	103	fZn(0),
	104	fB(0)
	105	{
	106	//
	107	// Copy constructor
	108	//
	109	magf.Copy(*this);
	110	}
	111
	112	//_______________________________________________________________________
	113	void AliMagFCM::Field(Float_t x, Float_t b)
	114	{
	115	//
	116	// Method to calculate the magnetic field
	117	//
	118	Double_t ratx, raty, ratz, hix, hiy, hiz, ratx1, raty1, ratz1,
	119	bhyhz, bhylz, blyhz, blylz, bhz, blz, xl[3];
	120	const Double_t kone=1;
	121	Int_t ix, iy, iz;
	122
	123	// --- find the position in the grid ---
	124
	125	b[0]=b[1]=b[2]=0;
	126
	127
	128	if(-700 < -x[2] && -x[2] < fZbeg && x[0] * x[0] +(x[1]+30)(x[1]+30) < 560560) {
	129	b[2]= fSolenoid;
	130	} else {
	131	// The field map used here was calculated in a coordinate system where the muon arm is at z > 0
	132	// Transfom x -> -x and z -> -z
	133	Float_t xm = - x[0];
	134	Float_t ym = x[1];
	135	Float_t zm = - x[2];
	136
	137	Bool_t infield=(fZbeg <= zm && zm < fZbeg+fZdel*(fZn-1)
	138	&& ( fXbeg <= TMath::Abs(xm) && TMath::Abs(xm) < fXbeg+fXdel*(fXn-1) )
	139	&& ( fYbeg <= TMath::Abs(ym) && TMath::Abs(ym) < fYbeg+fYdel*(fYn-1) ));
	140	if(infield) {
	141	xl[0]=TMath::Abs(xm)-fXbeg;
	142	xl[1]=TMath::Abs(ym)-fYbeg;
	143	xl[2]=zm-fZbeg;
	144
	145	// --- start with x
	146
	147	hix=xl[0]*fXdeli;
	148	ratx=hix-int(hix);
	149	ix=int(hix);
	150
	151	hiy=xl[1]*fYdeli;
	152	raty=hiy-int(hiy);
	153	iy=int(hiy);
	154
	155	hiz=xl[2]*fZdeli;
	156	ratz=hiz-int(hiz);
	157	iz=int(hiz);
	158
	159	if(fMap==2) {
	160	// ... simple interpolation
	161	ratx1=kone-ratx;
	162	raty1=kone-raty;
	163	ratz1=kone-ratz;
	164	bhyhz = Bx(ix ,iy+1,iz+1)ratx1+Bx(ix+1,iy+1,iz+1)ratx;
	165	bhylz = Bx(ix ,iy+1,iz )ratx1+Bx(ix+1,iy+1,iz )ratx;
	166	blyhz = Bx(ix ,iy ,iz+1)ratx1+Bx(ix+1,iy ,iz+1)ratx;
	167	blylz = Bx(ix ,iy ,iz )ratx1+Bx(ix+1,iy ,iz )ratx;
	168	bhz = blyhz raty1+bhyhz raty;
	169	blz = blylz raty1+bhylz raty;
	170	b[0] = blz ratz1+bhz ratz;
	171	//
	172	bhyhz = By(ix ,iy+1,iz+1)ratx1+By(ix+1,iy+1,iz+1)ratx;
	173	bhylz = By(ix ,iy+1,iz )ratx1+By(ix+1,iy+1,iz )ratx;
	174	blyhz = By(ix ,iy ,iz+1)ratx1+By(ix+1,iy ,iz+1)ratx;
	175	blylz = By(ix ,iy ,iz )ratx1+By(ix+1,iy ,iz )ratx;
	176	bhz = blyhz raty1+bhyhz raty;
	177	blz = blylz raty1+bhylz raty;
	178	b[1] = blz ratz1+bhz ratz;
	179	//
	180	bhyhz = Bz(ix ,iy+1,iz+1)ratx1+Bz(ix+1,iy+1,iz+1)ratx;
	181	bhylz = Bz(ix ,iy+1,iz )ratx1+Bz(ix+1,iy+1,iz )ratx;
	182	blyhz = Bz(ix ,iy ,iz+1)ratx1+Bz(ix+1,iy ,iz+1)ratx;
	183	blylz = Bz(ix ,iy ,iz )ratx1+Bz(ix+1,iy ,iz )ratx;
	184	bhz = blyhz raty1+bhyhz raty;
	185	blz = blylz raty1+bhylz raty;
	186	b[2] = blz ratz1+bhz ratz;
	187	//printf("ratx,raty,ratz,b[0],b[1],b[2] %f %f %f %f %f %f\n",
	188	//ratx,raty,ratz,b[0],b[1],b[2]);
	189	//
	190	// ... use the dipole symmetry
	191	if (xm*ym < 0) b[1]=-b[1];
	192	if (xm<0) b[2]=-b[2];
	193	b[0] = -b[0];
	194	b[2] = -b[2];
	195
	196	} else {
	197	printf("Invalid field map for constant mesh %d\n",fMap);
	198	}
	199	} else {
	200	//This is the ZDC part
	201	ZDCField(x,b);
	202	}
	203
	204	if(fFactor!=1) {
	205	b[0]*=fFactor;
	206	b[1]*=fFactor;
	207	b[2]*=fFactor;
	208	}
	209	}
	210	}
	211
	212	//_______________________________________________________________________
	213	void AliMagFCM::ReadField()
	214	{
	215	//
	216	// Method to read the magnetic field map from file
	217	//
	218	FILE *magfile;
	219	Int_t ix, iy, iz, ipx, ipy, ipz;
	220	Float_t bx, by, bz;
	221	char *fname;
	222	if(fDebug) printf("%s: Reading Magnetic Field %s from file %s\n",ClassName(),fName.Data(),fTitle.Data());
	223	fname = gSystem->ExpandPathName(fTitle.Data());
	224	magfile=fopen(fname,"r");
	225	delete [] fname;
	226	if (magfile) {
	227	fscanf(magfile,"%d %d %d %f %f %f %f %f %f",
	228	&fXn, &fYn, &fZn, &fXdel, &fYdel, &fZdel, &fXbeg, &fYbeg, &fZbeg);
	229	if(fDebug>1) printf("%s: fXn %d, fYn %d, fZn %d, fXdel %f, fYdel %f, fZdel %f, fXbeg %f, fYbeg %f, fZbeg %f\n",
	230	ClassName(),fXn, fYn, fZn, fXdel, fYdel, fZdel, fXbeg, fYbeg, fZbeg);
	231	fXdeli=1./fXdel;
	232	fYdeli=1./fYdel;
	233	fZdeli=1./fZdel;
	234	fB = new TVector(3fXnfYn*fZn);
	235	for (iz=0; iz<fZn; iz++) {
	236	ipz=iz3(fXn*fYn);
	237	for (iy=0; iy<fYn; iy++) {
	238	ipy=ipz+iy3fXn;
	239	for (ix=0; ix<fXn; ix++) {
	240	ipx=ipy+ix*3;
	241	fscanf(magfile,"%f %f %f",&bz,&by,&bx);
	242	(*fB)(ipx+2)=bz;
	243	(*fB)(ipx+1)=by;
	244	(*fB)(ipx )=bx;
	245	}
	246	}
	247	}
	248	} else {
	249	printf("%s: File %s not found !\n",ClassName(),fTitle.Data());
	250	exit(1);
	251	}
	252	}
	253
	254	//_______________________________________________________________________
	255	void AliMagFCM::Copy(TObject & /* magf */) const
	256	{
	257	//
	258	// Copy *this onto magf -- Not implemented
	259	//
	260	Fatal("Copy","Not implemented!\n");
	261	}
	262