[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.6.1  2002/07/24 10:08:13  alibrary
Updating VirtualMC

Revision 1.4  2002/07/10 15:06:22  morsch
Initialize TVector* save.

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