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