[u/mrichter/AliRoot.git] / STEER / AliMagFCM.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 for Alice magnetic field with constant mesh
//  Used in the configuration macros (macros/Config.C, etc.)
//  Author:
//-----------------------------------------------------------------------
#include "TSystem.h"

#include "TVector.h"

#include "AliMagFCM.h"

ClassImp(AliMagFCM)

//_______________________________________________________________________
AliMagFCM::AliMagFCM():
  fXbeg(0),
  fYbeg(0),
  fZbeg(0),
  fXdel(0),
  fYdel(0),
  fZdel(0),
  fSolenoid(0),
  fXdeli(0),
  fYdeli(0),
  fZdeli(0),
  fXn(0),
  fYn(0),
  fZn(0),
  fB(0)
{
  //
  // Standard constructor
  //
  fType = kConMesh;
  fMap  = 2;
  SetSolenoidField();
}

//_______________________________________________________________________
AliMagFCM::AliMagFCM(const char *name, const char *title, Int_t integ, 
                     Float_t factor, Float_t fmax):
  AliMagFC(name,title,integ,factor,fmax),
  fXbeg(0),
  fYbeg(0),
  fZbeg(0),
  fXdel(0),
  fYdel(0),
  fZdel(0),
  fSolenoid(0),
  fXdeli(0),
  fYdeli(0),
  fZdeli(0),
  fXn(0),
  fYn(0),
  fZn(0),
  fB(0)
{
  //
  // Standard constructor
  //
  fType = kConMesh;
  fMap  = 2;
  SetSolenoidField();

  if(fDebug>-1) Info("ctor",
     "%s: Constant Mesh Field %s created: map= %d, factor= %f, file= %s\n",
	 ClassName(),fName.Data(), fMap, factor,fTitle.Data());
}

//_______________________________________________________________________
AliMagFCM::AliMagFCM(const AliMagFCM &magf):
  AliMagFC(magf),
  fXbeg(0),
  fYbeg(0),
  fZbeg(0),
  fXdel(0),
  fYdel(0),
  fZdel(0),
  fSolenoid(0),
  fXdeli(0),
  fYdeli(0),
  fZdeli(0),
  fXn(0),
  fYn(0),
  fZn(0),
  fB(0)
{
  //
  // Copy constructor
  //
  magf.Copy(*this);
}

//_______________________________________________________________________
void AliMagFCM::Field(Float_t *x, Float_t *b)
{
  //
  // Method to calculate the magnetic field
  //
  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;
  
  // --- find the position in the grid ---
  
  b[0]=b[1]=b[2]=0;

  
  if(-700 < -x[2] && -x[2] < fZbeg && x[0] * x[0] +(x[1]+30)*(x[1]+30) < 560*560) {
    b[2]=  fSolenoid;
  } else  {
      // The field map used here was calculated in a coordinate system where the muon arm is at z > 0
      // Transfom x -> -x and z -> -z 
      Float_t xm = - x[0];
      Float_t ym =   x[1];
      Float_t zm = - x[2];
 
      Bool_t infield=(fZbeg <= zm  && zm < fZbeg+fZdel*(fZn-1)
		      &&  ( fXbeg <= TMath::Abs(xm) && TMath::Abs(xm) < fXbeg+fXdel*(fXn-1) )
		      &&  ( fYbeg <= TMath::Abs(ym) && TMath::Abs(ym) < fYbeg+fYdel*(fYn-1) ));
    if(infield) {
      xl[0]=TMath::Abs(xm)-fXbeg;
      xl[1]=TMath::Abs(ym)-fYbeg;
      xl[2]=zm-fZbeg;
      
      // --- start with x
      
      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);
      
      if(fMap==2) {
	// ... simple interpolation
	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;
	//printf("ratx,raty,ratz,b[0],b[1],b[2] %f %f %f %f %f %f\n",
	//ratx,raty,ratz,b[0],b[1],b[2]);
	//
	// ... use the dipole symmetry
	if (xm*ym < 0) b[1]=-b[1];
	if (xm<0)      b[2]=-b[2];
	b[0] = -b[0];
	b[2] = -b[2];
	
      } else {
	printf("Invalid field map for constant mesh %d\n",fMap);
      }
    } else {
//This is the ZDC part
	ZDCField(x,b);
    }
    
    if(fFactor!=1) {
	b[0]*=fFactor;
	b[1]*=fFactor;
	b[2]*=fFactor;
    }
  }
}

//_______________________________________________________________________
void AliMagFCM::ReadField()
{
  // 
  // Method to read the magnetic field map from file
  //
  FILE *magfile;
  Int_t ix, iy, iz, ipx, ipy, ipz;
  Float_t bx, by, bz;
  char *fname;
  if(fDebug) printf("%s: Reading Magnetic Field %s from file %s\n",ClassName(),fName.Data(),fTitle.Data());
  fname = gSystem->ExpandPathName(fTitle.Data());
  magfile=fopen(fname,"r");
  delete [] fname;
  if (magfile) {
    fscanf(magfile,"%d %d %d %f %f %f %f %f %f",
    	   &fXn, &fYn, &fZn, &fXdel, &fYdel, &fZdel, &fXbeg, &fYbeg, &fZbeg);
    if(fDebug>1) printf("%s: fXn %d, fYn %d, fZn %d, fXdel %f, fYdel %f, fZdel %f, fXbeg %f, fYbeg %f, fZbeg %f\n",
			ClassName(),fXn, fYn, fZn, fXdel, fYdel, fZdel, fXbeg, fYbeg, fZbeg);
    fXdeli=1./fXdel;
    fYdeli=1./fYdel;
    fZdeli=1./fZdel;
    fB = new TVector(3*fXn*fYn*fZn);
    for (iz=0; iz<fZn; iz++) {
      ipz=iz*3*(fXn*fYn);
      for (iy=0; iy<fYn; iy++) {
	ipy=ipz+iy*3*fXn;
	for (ix=0; ix<fXn; ix++) {
	  ipx=ipy+ix*3;
	  fscanf(magfile,"%f %f %f",&bz,&by,&bx);
	  (*fB)(ipx+2)=bz;
	  (*fB)(ipx+1)=by;
	  (*fB)(ipx  )=bx;
	}
      }
    }
  } else { 
    printf("%s: File %s not found !\n",ClassName(),fTitle.Data());
    exit(1);
  }
}

//_______________________________________________________________________
void AliMagFCM::Copy(TObject & /* magf */) const
{
  //
  // Copy *this onto magf -- Not implemented
  //
  Fatal("Copy","Not implemented!\n");
}
Commit	Line	Data
aee8290b	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$ */
e2afb3b6	17
5d8718b8	18	//-----------------------------------------------------------------------
	19	// Class for Alice magnetic field with constant mesh
	20	// Used in the configuration macros (macros/Config.C, etc.)
	21	// Author:
	22	//-----------------------------------------------------------------------
fb17acd4	23	#include "TSystem.h"
88cb7938	24
65fb704d	25	#include "TVector.h"
aee8290b	26
aee8290b	27	#include "AliMagFCM.h"
aee8290b	28
	29	ClassImp(AliMagFCM)
	30
e2afb3b6	31	//_______________________________________________________________________
	32	AliMagFCM::AliMagFCM():
	33	fXbeg(0),
	34	fYbeg(0),
	35	fZbeg(0),
	36	fXdel(0),
	37	fYdel(0),
	38	fZdel(0),
	39	fSolenoid(0),
	40	fXdeli(0),
	41	fYdeli(0),
	42	fZdeli(0),
	43	fXn(0),
	44	fYn(0),
	45	fZn(0),
	46	fB(0)
	47	{
	48	//
	49	// Standard constructor
	50	//
	51	fType = kConMesh;
	52	fMap = 2;
	53	SetSolenoidField();
	54	}
	55
	56	//_______________________________________________________________________
d0f1ee3b	57	AliMagFCM::AliMagFCM(const char name, const char title, Int_t integ,
d0f1ee3b	58	Float_t factor, Float_t fmax):
57754f18	59	AliMagFC(name,title,integ,factor,fmax),
e2afb3b6	60	fXbeg(0),
	61	fYbeg(0),
	62	fZbeg(0),
	63	fXdel(0),
	64	fYdel(0),
	65	fZdel(0),
	66	fSolenoid(0),
	67	fXdeli(0),
	68	fYdeli(0),
	69	fZdeli(0),
	70	fXn(0),
	71	fYn(0),
	72	fZn(0),
	73	fB(0)
aee8290b	74	{
	75	//
	76	// Standard constructor
	77	//
	78	fType = kConMesh;
d8408e76	79	fMap = 2;
4cc8933f	80	SetSolenoidField();
4cc8933f	81
e2afb3b6	82	if(fDebug>-1) Info("ctor",
e2afb3b6	83	"%s: Constant Mesh Field %s created: map= %d, factor= %f, file= %s\n",
9e1a0ddb	84	ClassName(),fName.Data(), fMap, factor,fTitle.Data());
aee8290b	85	}
aee8290b	86
e2afb3b6	87	//_______________________________________________________________________
e2afb3b6	88	AliMagFCM::AliMagFCM(const AliMagFCM &magf):
57754f18	89	AliMagFC(magf),
e2afb3b6	90	fXbeg(0),
	91	fYbeg(0),
	92	fZbeg(0),
	93	fXdel(0),
	94	fYdel(0),
	95	fZdel(0),
	96	fSolenoid(0),
	97	fXdeli(0),
	98	fYdeli(0),
	99	fZdeli(0),
	100	fXn(0),
	101	fYn(0),
	102	fZn(0),
	103	fB(0)
aee8290b	104	{
	105	//
	106	// Copy constructor
	107	//
	108	magf.Copy(*this);
	109	}
	110
e2afb3b6	111	//_______________________________________________________________________
aee8290b	112	void AliMagFCM::Field(Float_t x, Float_t b)
	113	{
	114	//
	115	// Method to calculate the magnetic field
	116	//
	117	Double_t ratx, raty, ratz, hix, hiy, hiz, ratx1, raty1, ratz1,
	118	bhyhz, bhylz, blyhz, blylz, bhz, blz, xl[3];
	119	const Double_t kone=1;
	120	Int_t ix, iy, iz;
	121
	122	// --- find the position in the grid ---
	123
	124	b[0]=b[1]=b[2]=0;
57754f18	125
	126
	127	if(-700 < -x[2] && -x[2] < fZbeg && x[0] * x[0] +(x[1]+30)(x[1]+30) < 560560) {
	128	b[2]= fSolenoid;
aee8290b	129	} else {
57754f18	130	// The field map used here was calculated in a coordinate system where the muon arm is at z > 0
	131	// Transfom x -> -x and z -> -z
	132	Float_t xm = - x[0];
	133	Float_t ym = x[1];
	134	Float_t zm = - x[2];
	135
	136	Bool_t infield=(fZbeg <= zm && zm < fZbeg+fZdel*(fZn-1)
	137	&& ( fXbeg <= TMath::Abs(xm) && TMath::Abs(xm) < fXbeg+fXdel*(fXn-1) )
	138	&& ( fYbeg <= TMath::Abs(ym) && TMath::Abs(ym) < fYbeg+fYdel*(fYn-1) ));
aee8290b	139	if(infield) {
57754f18	140	xl[0]=TMath::Abs(xm)-fXbeg;
	141	xl[1]=TMath::Abs(ym)-fYbeg;
	142	xl[2]=zm-fZbeg;
aee8290b	143
	144	// --- start with x
	145
	146	hix=xl[0]*fXdeli;
	147	ratx=hix-int(hix);
	148	ix=int(hix);
	149
	150	hiy=xl[1]*fYdeli;
	151	raty=hiy-int(hiy);
	152	iy=int(hiy);
	153
	154	hiz=xl[2]*fZdeli;
	155	ratz=hiz-int(hiz);
	156	iz=int(hiz);
	157
	158	if(fMap==2) {
	159	// ... simple interpolation
	160	ratx1=kone-ratx;
	161	raty1=kone-raty;
	162	ratz1=kone-ratz;
	163	bhyhz = Bx(ix ,iy+1,iz+1)ratx1+Bx(ix+1,iy+1,iz+1)ratx;
	164	bhylz = Bx(ix ,iy+1,iz )ratx1+Bx(ix+1,iy+1,iz )ratx;
	165	blyhz = Bx(ix ,iy ,iz+1)ratx1+Bx(ix+1,iy ,iz+1)ratx;
	166	blylz = Bx(ix ,iy ,iz )ratx1+Bx(ix+1,iy ,iz )ratx;
	167	bhz = blyhz raty1+bhyhz raty;
	168	blz = blylz raty1+bhylz raty;
	169	b[0] = blz ratz1+bhz ratz;
	170	//
	171	bhyhz = By(ix ,iy+1,iz+1)ratx1+By(ix+1,iy+1,iz+1)ratx;
	172	bhylz = By(ix ,iy+1,iz )ratx1+By(ix+1,iy+1,iz )ratx;
	173	blyhz = By(ix ,iy ,iz+1)ratx1+By(ix+1,iy ,iz+1)ratx;
	174	blylz = By(ix ,iy ,iz )ratx1+By(ix+1,iy ,iz )ratx;
	175	bhz = blyhz raty1+bhyhz raty;
	176	blz = blylz raty1+bhylz raty;
	177	b[1] = blz ratz1+bhz ratz;
	178	//
	179	bhyhz = Bz(ix ,iy+1,iz+1)ratx1+Bz(ix+1,iy+1,iz+1)ratx;
	180	bhylz = Bz(ix ,iy+1,iz )ratx1+Bz(ix+1,iy+1,iz )ratx;
	181	blyhz = Bz(ix ,iy ,iz+1)ratx1+Bz(ix+1,iy ,iz+1)ratx;
	182	blylz = Bz(ix ,iy ,iz )ratx1+Bz(ix+1,iy ,iz )ratx;
	183	bhz = blyhz raty1+bhyhz raty;
	184	blz = blylz raty1+bhylz raty;
	185	b[2] = blz ratz1+bhz ratz;
	186	//printf("ratx,raty,ratz,b[0],b[1],b[2] %f %f %f %f %f %f\n",
	187	//ratx,raty,ratz,b[0],b[1],b[2]);
	188	//
	189	// ... use the dipole symmetry
57754f18	190	if (xm*ym < 0) b[1]=-b[1];
	191	if (xm<0) b[2]=-b[2];
	192	b[0] = -b[0];
	193	b[2] = -b[2];
	194
aee8290b	195	} else {
	196	printf("Invalid field map for constant mesh %d\n",fMap);
	197	}
	198	} else {
d95ea23f	199	//This is the ZDC part
57754f18	200	ZDCField(x,b);
aee8290b	201	}
d95ea23f	202
57754f18	203	if(fFactor!=1) {
	204	b[0]*=fFactor;
	205	b[1]*=fFactor;
	206	b[2]*=fFactor;
d95ea23f	207	}
aee8290b	208	}
aee8290b	209	}
aee8290b	210
e2afb3b6	211	//_______________________________________________________________________
aee8290b	212	void AliMagFCM::ReadField()
	213	{
	214	//
	215	// Method to read the magnetic field map from file
	216	//
	217	FILE *magfile;
	218	Int_t ix, iy, iz, ipx, ipy, ipz;
	219	Float_t bx, by, bz;
	220	char *fname;
9e1a0ddb	221	if(fDebug) printf("%s: Reading Magnetic Field %s from file %s\n",ClassName(),fName.Data(),fTitle.Data());
aee8290b	222	fname = gSystem->ExpandPathName(fTitle.Data());
	223	magfile=fopen(fname,"r");
	224	delete [] fname;
	225	if (magfile) {
	226	fscanf(magfile,"%d %d %d %f %f %f %f %f %f",
	227	&fXn, &fYn, &fZn, &fXdel, &fYdel, &fZdel, &fXbeg, &fYbeg, &fZbeg);
9e1a0ddb	228	if(fDebug>1) printf("%s: fXn %d, fYn %d, fZn %d, fXdel %f, fYdel %f, fZdel %f, fXbeg %f, fYbeg %f, fZbeg %f\n",
9e1a0ddb	229	ClassName(),fXn, fYn, fZn, fXdel, fYdel, fZdel, fXbeg, fYbeg, fZbeg);
aee8290b	230	fXdeli=1./fXdel;
	231	fYdeli=1./fYdel;
	232	fZdeli=1./fZdel;
	233	fB = new TVector(3fXnfYn*fZn);
	234	for (iz=0; iz<fZn; iz++) {
	235	ipz=iz3(fXn*fYn);
	236	for (iy=0; iy<fYn; iy++) {
	237	ipy=ipz+iy3fXn;
	238	for (ix=0; ix<fXn; ix++) {
	239	ipx=ipy+ix*3;
	240	fscanf(magfile,"%f %f %f",&bz,&by,&bx);
	241	(*fB)(ipx+2)=bz;
	242	(*fB)(ipx+1)=by;
	243	(*fB)(ipx )=bx;
	244	}
	245	}
	246	}
	247	} else {
9e1a0ddb	248	printf("%s: File %s not found !\n",ClassName(),fTitle.Data());
aee8290b	249	exit(1);
	250	}
	251	}
	252
e2afb3b6	253	//_______________________________________________________________________
6c4904c2	254	void AliMagFCM::Copy(TObject & /* magf */) const
aee8290b	255	{
aee8290b	256	//
8918e700	257	// Copy *this onto magf -- Not implemented
aee8290b	258	//
	259	Fatal("Copy","Not implemented!\n");
	260	}
	261