[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 <cstdlib>
#include <TClass.h>
#include <TSystem.h>

#include "AliLog.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;
  AliInfo(Form("Reading Magnetic Field Map %s from file %s",
	       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] = x[0] - fXbeg;
    xl[1] = x[1] - fYbeg;
    xl[2] = x[2] - fZbeg;
    
    hix=TMath::Max(0.,TMath::Min(xl[0]*fXdeli,fXn-1.0001));
    ratx=hix-int(hix);
    ix=int(hix);
    
    hiy=TMath::Max(0.,TMath::Min(xl[1]*fYdeli,fYn-1.0001));
    raty=hiy-int(hiy);
    iy=int(hiy);
    
    hiz=TMath::Max(0.,TMath::Min(xl[2]*fZdeli,fZn-1.0001));
    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
  //
  AliFatal("Not implemented!");
}

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