[u/mrichter/AliRoot.git] / STEER / AliFieldMap.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 to handle the field
// I/O and interpolation
// of the field map to return the value
// of the magnetic field in an arbitrary position
// Author: Andreas Morsch <andreas.morsch@cern.ch>
//
//-----------------------------------------------------------------------

#include <TSystem.h>

#include "AliFieldMap.h"

ClassImp(AliFieldMap)

//_______________________________________________________________________
AliFieldMap::AliFieldMap():
  fXbeg(0),
  fYbeg(0),
  fZbeg(0),
  fXend(0),
  fYend(0),
  fZend(0),
  fXdel(0),
  fYdel(0),
  fZdel(0),
  fXdeli(0),
  fYdeli(0),
  fZdeli(0),
  fXn(0),
  fYn(0),
  fZn(0),
  fWriteEnable(0),
  fB(0)
{
  //
  // Standard constructor
  //
  SetWriteEnable();
}

//_______________________________________________________________________
AliFieldMap::AliFieldMap(const char *name, const char *title):
  TNamed(name,title),
  fXbeg(0),
  fYbeg(0),
  fZbeg(0),
  fXend(0),
  fYend(0),
  fZend(0),
  fXdel(0),
  fYdel(0),
  fZdel(0),
  fXdeli(0),
  fYdeli(0),
  fZdeli(0),
  fXn(0),
  fYn(0),
  fZn(0),
  fWriteEnable(0),
  fB(0)
{
  //
  // Standard constructor
  //
  ReadField();
  SetWriteEnable();
}

//_______________________________________________________________________
AliFieldMap::~AliFieldMap()
{
  //
  // Destructor
  //  
  delete fB;
}

//_______________________________________________________________________
AliFieldMap::AliFieldMap(const AliFieldMap &map):
  TNamed(map),
  fXbeg(0),
  fYbeg(0),
  fZbeg(0),
  fXend(0),
  fYend(0),
  fZend(0),
  fXdel(0),
  fYdel(0),
  fZdel(0),
  fXdeli(0),
  fYdeli(0),
  fZdeli(0),
  fXn(0),
  fYn(0),
  fZn(0),
  fWriteEnable(0),
  fB(0)
{
  //
  // Copy constructor
  //
  map.Copy(*this);
}

//_______________________________________________________________________
void AliFieldMap::ReadField()
{
  // 
  // Method to read the magnetic field map from file
  //
  FILE* magfile;
  //  FILE* endf = fopen("end.table", "r");
  //  FILE* out  = fopen("out", "w");
  
  Int_t   ix, iy, iz, ipx, ipy, ipz;
  Float_t bx, by, bz;
  char *fname = 0;
  printf("%s: Reading Magnetic Field Map %s from file %s\n",
	 ClassName(),fName.Data(),fTitle.Data());

  fname   = gSystem->ExpandPathName(fTitle.Data());
  magfile = fopen(fname,"r");
  delete [] fname;

  fscanf(magfile,"%d %d %d %f %f %f %f %f %f", 
	 &fXn, &fYn, &fZn, &fXbeg, &fYbeg, &fZbeg, &fXdel, &fYdel, &fZdel);
 
  fXdeli = 1./fXdel;
  fYdeli = 1./fYdel;	
  fZdeli = 1./fZdel;	
  fXend  = fXbeg + (fXn-1)*fXdel;
  fYend  = fYbeg + (fYn-1)*fYdel;
  fZend  = fZbeg + (fZn-1)*fZdel;
  
  Int_t nDim   = fXn*fYn*fZn;

  //  Float_t x,y,z,b;

  fB = new TVector(3*nDim);
  if (magfile) {
      for (ix = 0; ix < fXn; ix++) {
	  ipx=ix*3*(fZn*fYn);
	  for (iy = 0; iy < fYn; iy++) {
	      ipy=ipx+iy*3*fZn;
	      for (iz = 0; iz < fZn; iz++) {
		  ipz=ipy+iz*3;

		  if (iz == -1) {
//		      fscanf(endf,"%f %f %f", &bx,&by,&bz);
		  } else if (iz > -1) {
		      fscanf(magfile," %f %f %f", &bx, &by, &bz);
		  } else {
		      continue;
		  }
		  
//		  fscanf(magfile,"%f %f %f %f %f %f %f ",
//			 &x, &y, &z, &bx,&by,&bz, &b);
//		  fprintf(out, "%15.8e %15.8e %15.8e \n", bx, by, bz);

		  (*fB)(ipz+2) = 10.*bz;
		  (*fB)(ipz+1) = 10.*by;
		  (*fB)(ipz  ) = 10.*bx;
	      } //iz
	  } // iy
      } // ix
/*
//
// this part for interpolation between z = 700 and 720 cm to get
// z = 710 cm
//      
      for (ix = 0; ix < fXn; ix++) {
	  ipx=ix*3*(fZn*fYn);
	  for (iy = 0; iy < fYn; iy++) {
	      ipy=ipx+iy*3*fZn;
	      Float_t bxx = (Bx(ix,iy,0) + Bx(ix,iy,2))/2.;
	      Float_t byy = (By(ix,iy,0) + By(ix,iy,2))/2.;
	      Float_t bzz = (Bz(ix,iy,0) + Bz(ix,iy,2))/2.;	      
	      ipz=ipy+3;
	      (*fB)(ipz+2) = bzz;
	      (*fB)(ipz+1) = byy;
	      (*fB)(ipz  ) = bxx;
	  } // iy
      } // ix
*/      
  } else { 
      printf("%s: File %s not found !\n",ClassName(),fTitle.Data());
      exit(1);
  } // if mafile
}

//_______________________________________________________________________
void AliFieldMap::Field(Float_t *x, Float_t *b) const
{
  //
  // Use simple interpolation to obtain field at point x
  //
    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;
    b[0]=b[1]=b[2]=0;
    //
    
    xl[0]=TMath::Abs(x[0])-fXbeg;
    xl[1]=TMath::Abs(x[1])-fYbeg;
    xl[2]=x[2]-fZbeg;
    
    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);

    ratx1=kone-ratx;
    raty1=kone-raty;
    ratz1=kone-ratz;

    if (!fB) return;

    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;
}

//_______________________________________________________________________
void AliFieldMap::Copy(TObject & /* magf */) const
{
  //
  // Copy *this onto magf -- Not implemented
  //
  Fatal("Copy","Not implemented!\n");
}

//_______________________________________________________________________
AliFieldMap & AliFieldMap::operator =(const AliFieldMap &magf)
{
  magf.Copy(*this);
  return *this;
}

//_______________________________________________________________________
void AliFieldMap::Streamer(TBuffer &R__b)
{
   // Stream an object of class AliFieldMap.
    TVector* save = 0;
    
    if (R__b.IsReading()) {
	AliFieldMap::Class()->ReadBuffer(R__b, this);
    } else {
	if (!fWriteEnable) {
	    save = fB;
	    fB = 0;
	}
	AliFieldMap::Class()->WriteBuffer(R__b, this);
	if (!fWriteEnable) fB = save;
    }
}
Commit	Line	Data
84737f5e	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$ */
84737f5e	17
116cbefd	18	//-----------------------------------------------------------------------
90dd3681	19	//
af7ba10c	20	// Class to handle the field
	21	// I/O and interpolation
	22	// of the field map to return the value
	23	// of the magnetic field in an arbitrary position
84737f5e	24	// Author: Andreas Morsch <andreas.morsch@cern.ch>
90dd3681	25	//
116cbefd	26	//-----------------------------------------------------------------------
84737f5e	27
116cbefd	28	#include <TSystem.h>
116cbefd	29
84737f5e	30	#include "AliFieldMap.h"
84737f5e	31
	32	ClassImp(AliFieldMap)
	33
e2afb3b6	34	//_______________________________________________________________________
	35	AliFieldMap::AliFieldMap():
	36	fXbeg(0),
	37	fYbeg(0),
	38	fZbeg(0),
	39	fXend(0),
	40	fYend(0),
	41	fZend(0),
	42	fXdel(0),
	43	fYdel(0),
	44	fZdel(0),
	45	fXdeli(0),
	46	fYdeli(0),
	47	fZdeli(0),
	48	fXn(0),
	49	fYn(0),
	50	fZn(0),
	51	fWriteEnable(0),
	52	fB(0)
84737f5e	53	{
	54	//
	55	// Standard constructor
	56	//
7b6cddfa	57	SetWriteEnable();
84737f5e	58	}
84737f5e	59
e2afb3b6	60	//_______________________________________________________________________
	61	AliFieldMap::AliFieldMap(const char name, const char title):
	62	TNamed(name,title),
	63	fXbeg(0),
	64	fYbeg(0),
	65	fZbeg(0),
	66	fXend(0),
	67	fYend(0),
	68	fZend(0),
	69	fXdel(0),
	70	fYdel(0),
	71	fZdel(0),
	72	fXdeli(0),
	73	fYdeli(0),
	74	fZdeli(0),
	75	fXn(0),
	76	fYn(0),
	77	fZn(0),
	78	fWriteEnable(0),
	79	fB(0)
84737f5e	80	{
	81	//
	82	// Standard constructor
	83	//
84737f5e	84	ReadField();
7b6cddfa	85	SetWriteEnable();
84737f5e	86	}
84737f5e	87
e2afb3b6	88	//_______________________________________________________________________
84737f5e	89	AliFieldMap::~AliFieldMap()
84737f5e	90	{
e2afb3b6	91	//
	92	// Destructor
	93	//
84737f5e	94	delete fB;
	95	}
	96
e2afb3b6	97	//_______________________________________________________________________
	98	AliFieldMap::AliFieldMap(const AliFieldMap &map):
	99	TNamed(map),
	100	fXbeg(0),
	101	fYbeg(0),
	102	fZbeg(0),
	103	fXend(0),
	104	fYend(0),
	105	fZend(0),
	106	fXdel(0),
	107	fYdel(0),
	108	fZdel(0),
	109	fXdeli(0),
	110	fYdeli(0),
	111	fZdeli(0),
	112	fXn(0),
	113	fYn(0),
	114	fZn(0),
	115	fWriteEnable(0),
	116	fB(0)
84737f5e	117	{
	118	//
	119	// Copy constructor
	120	//
	121	map.Copy(*this);
	122	}
	123
e2afb3b6	124	//_______________________________________________________________________
84737f5e	125	void AliFieldMap::ReadField()
	126	{
	127	//
	128	// Method to read the magnetic field map from file
	129	//
	130	FILE* magfile;
e2afb3b6	131	// FILE* endf = fopen("end.table", "r");
e2afb3b6	132	// FILE* out = fopen("out", "w");
84737f5e	133
	134	Int_t ix, iy, iz, ipx, ipy, ipz;
	135	Float_t bx, by, bz;
	136	char *fname = 0;
	137	printf("%s: Reading Magnetic Field Map %s from file %s\n",
	138	ClassName(),fName.Data(),fTitle.Data());
	139
	140	fname = gSystem->ExpandPathName(fTitle.Data());
	141	magfile = fopen(fname,"r");
	142	delete [] fname;
	143
	144	fscanf(magfile,"%d %d %d %f %f %f %f %f %f",
	145	&fXn, &fYn, &fZn, &fXbeg, &fYbeg, &fZbeg, &fXdel, &fYdel, &fZdel);
	146
	147	fXdeli = 1./fXdel;
	148	fYdeli = 1./fYdel;
	149	fZdeli = 1./fZdel;
	150	fXend = fXbeg + (fXn-1)*fXdel;
	151	fYend = fYbeg + (fYn-1)*fYdel;
	152	fZend = fZbeg + (fZn-1)*fZdel;
	153
	154	Int_t nDim = fXnfYnfZn;
	155
e2afb3b6	156	// Float_t x,y,z,b;
84737f5e	157
	158	fB = new TVector(3*nDim);
	159	if (magfile) {
	160	for (ix = 0; ix < fXn; ix++) {
	161	ipx=ix3(fZn*fYn);
	162	for (iy = 0; iy < fYn; iy++) {
	163	ipy=ipx+iy3fZn;
	164	for (iz = 0; iz < fZn; iz++) {
	165	ipz=ipy+iz*3;
	166
7b6cddfa	167	if (iz == -1) {
	168	// fscanf(endf,"%f %f %f", &bx,&by,&bz);
	169	} else if (iz > -1) {
84737f5e	170	fscanf(magfile," %f %f %f", &bx, &by, &bz);
7b6cddfa	171	} else {
84737f5e	172	continue;
7b6cddfa	173	}
7b6cddfa	174
84737f5e	175	// fscanf(magfile,"%f %f %f %f %f %f %f ",
	176	// &x, &y, &z, &bx,&by,&bz, &b);
	177	// fprintf(out, "%15.8e %15.8e %15.8e \n", bx, by, bz);
	178
	179	(fB)(ipz+2) = 10.bz;
	180	(fB)(ipz+1) = 10.by;
	181	(fB)(ipz ) = 10.bx;
	182	} //iz
	183	} // iy
	184	} // ix
	185	/*
	186	//
	187	// this part for interpolation between z = 700 and 720 cm to get
	188	// z = 710 cm
	189	//
	190	for (ix = 0; ix < fXn; ix++) {
	191	ipx=ix3(fZn*fYn);
	192	for (iy = 0; iy < fYn; iy++) {
	193	ipy=ipx+iy3fZn;
	194	Float_t bxx = (Bx(ix,iy,0) + Bx(ix,iy,2))/2.;
	195	Float_t byy = (By(ix,iy,0) + By(ix,iy,2))/2.;
	196	Float_t bzz = (Bz(ix,iy,0) + Bz(ix,iy,2))/2.;
	197	ipz=ipy+3;
	198	(*fB)(ipz+2) = bzz;
	199	(*fB)(ipz+1) = byy;
	200	(*fB)(ipz ) = bxx;
	201	} // iy
	202	} // ix
	203	*/
	204	} else {
	205	printf("%s: File %s not found !\n",ClassName(),fTitle.Data());
	206	exit(1);
	207	} // if mafile
	208	}
	209
5d8718b8	210	//_______________________________________________________________________
af7ba10c	211	void AliFieldMap::Field(Float_t x, Float_t b) const
5d8718b8	212	{
	213	//
	214	// Use simple interpolation to obtain field at point x
90dd3681	215	//
	216	Double_t ratx, raty, ratz, hix, hiy, hiz, ratx1, raty1, ratz1,
	217	bhyhz, bhylz, blyhz, blylz, bhz, blz, xl[3];
	218	const Double_t kone=1;
	219	Int_t ix, iy, iz;
	220	b[0]=b[1]=b[2]=0;
	221	//
	222
	223	xl[0]=TMath::Abs(x[0])-fXbeg;
	224	xl[1]=TMath::Abs(x[1])-fYbeg;
	225	xl[2]=x[2]-fZbeg;
	226
	227	hix=xl[0]*fXdeli;
	228	ratx=hix-int(hix);
	229	ix=int(hix);
	230
	231	hiy=xl[1]*fYdeli;
	232	raty=hiy-int(hiy);
	233	iy=int(hiy);
	234
	235	hiz=xl[2]*fZdeli;
	236	ratz=hiz-int(hiz);
	237	iz=int(hiz);
	238
	239	ratx1=kone-ratx;
	240	raty1=kone-raty;
	241	ratz1=kone-ratz;
	242
39b41d10	243	if (!fB) return;
39b41d10	244
90dd3681	245	bhyhz = Bx(ix ,iy+1,iz+1)ratx1+Bx(ix+1,iy+1,iz+1)ratx;
	246	bhylz = Bx(ix ,iy+1,iz )ratx1+Bx(ix+1,iy+1,iz )ratx;
	247	blyhz = Bx(ix ,iy ,iz+1)ratx1+Bx(ix+1,iy ,iz+1)ratx;
	248	blylz = Bx(ix ,iy ,iz )ratx1+Bx(ix+1,iy ,iz )ratx;
	249	bhz = blyhz raty1+bhyhz raty;
	250	blz = blylz raty1+bhylz raty;
	251	b[0] = blz ratz1+bhz ratz;
	252	//
	253	bhyhz = By(ix ,iy+1,iz+1)ratx1+By(ix+1,iy+1,iz+1)ratx;
	254	bhylz = By(ix ,iy+1,iz )ratx1+By(ix+1,iy+1,iz )ratx;
	255	blyhz = By(ix ,iy ,iz+1)ratx1+By(ix+1,iy ,iz+1)ratx;
	256	blylz = By(ix ,iy ,iz )ratx1+By(ix+1,iy ,iz )ratx;
	257	bhz = blyhz raty1+bhyhz raty;
	258	blz = blylz raty1+bhylz raty;
	259	b[1] = blz ratz1+bhz ratz;
	260	//
	261	bhyhz = Bz(ix ,iy+1,iz+1)ratx1+Bz(ix+1,iy+1,iz+1)ratx;
	262	bhylz = Bz(ix ,iy+1,iz )ratx1+Bz(ix+1,iy+1,iz )ratx;
	263	blyhz = Bz(ix ,iy ,iz+1)ratx1+Bz(ix+1,iy ,iz+1)ratx;
	264	blylz = Bz(ix ,iy ,iz )ratx1+Bz(ix+1,iy ,iz )ratx;
	265	bhz = blyhz raty1+bhyhz raty;
	266	blz = blylz raty1+bhylz raty;
	267	b[2] = blz ratz1+bhz ratz;
84737f5e	268	}
84737f5e	269
e2afb3b6	270	//_______________________________________________________________________
6c4904c2	271	void AliFieldMap::Copy(TObject & /* magf */) const
84737f5e	272	{
	273	//
	274	// Copy *this onto magf -- Not implemented
	275	//
	276	Fatal("Copy","Not implemented!\n");
	277	}
	278
e2afb3b6	279	//_______________________________________________________________________
84737f5e	280	AliFieldMap & AliFieldMap::operator =(const AliFieldMap &magf)
	281	{
	282	magf.Copy(*this);
	283	return *this;
	284	}
7b6cddfa	285
e2afb3b6	286	//_______________________________________________________________________
7b6cddfa	287	void AliFieldMap::Streamer(TBuffer &R__b)
	288	{
	289	// Stream an object of class AliFieldMap.
f9aab227	290	TVector* save = 0;
b4cfe01e	291
	292	if (R__b.IsReading()) {
	293	AliFieldMap::Class()->ReadBuffer(R__b, this);
	294	} else {
	295	if (!fWriteEnable) {
	296	save = fB;
	297	fB = 0;
	298	}
	299	AliFieldMap::Class()->WriteBuffer(R__b, this);
	300	if (!fWriteEnable) fB = save;
	301	}
7b6cddfa	302	}