[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.                  *
 **************************************************************************/

/*
$Log$
Revision 1.3  2002/02/26 16:13:07  morsch
Correction in streamer.

Revision 1.2  2002/02/22 14:00:20  morsch
Protection against replication of fieldmap data in gAlice.

Revision 1.1  2002/02/14 11:41:28  morsch
Magnetic field map for ALICE for L3+muon spectrometer stored in 3 seperate
root files.

*/

//
// Author: Andreas Morsch <andreas.morsch@cern.ch>
//

#include <TVector.h>
#include "AliFieldMap.h"
#include "TSystem.h"

ClassImp(AliFieldMap)

//________________________________________
AliFieldMap::AliFieldMap()
{
  //
  // Standard constructor
  //
  fB = 0;
  SetWriteEnable();
}

AliFieldMap::AliFieldMap(const char *name, const char *title)
    : TNamed(name,title)
{
  //
  // Standard constructor
  //
  fB = 0;
  ReadField();
  SetWriteEnable();
}

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

//________________________________________
AliFieldMap::AliFieldMap(const AliFieldMap &map)
{
  //
  // 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)
{
//
// 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;

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