STEER/AliMagFCM.cxx

   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.9  2001/09/04 15:05:13  hristov
  19 Pointer fB is initialised to 0 in the default constructor
  20
  21 Revision 1.8  2001/05/28 14:10:35  morsch
  22 SetSolenoidField method to set the L3 field strength. 2 kG is default.
  23
  24 Revision 1.7  2001/05/16 14:57:22  alibrary
  25 New files for folders and Stack
  26
  27 Revision 1.6  2000/12/18 10:44:01  morsch
  28 Possibility to set field map by passing pointer to objet of type AliMagF via
  29 SetField().
  30 Example:
  31 gAlice->SetField(new AliMagFCM("Map2", "$(ALICE_ROOT)/data/field01.dat",2,1.,10.));
  32
  33 Revision 1.5  2000/12/01 11:20:27  alibrary
  34 Corrector dipole removed from ZDC
  35
  36 Revision 1.4  2000/11/30 07:12:49  alibrary
  37 Introducing new Rndm and QA classes
  38
  39 Revision 1.3  2000/11/10 18:09:55  fca
  40 New field map for the ZDC
  41
  42 Revision 1.2  2000/07/12 08:56:25  fca
  43 Coding convention correction and warning removal
  44
  45 Revision 1.1  2000/07/11 18:24:59  fca
  46 Coding convention corrections + few minor bug fixes
  47
  48 */
  49 #include "TVector.h"
  50
  51 #include "AliMagFCM.h"
  52 #include "TSystem.h"
  53
  54 ClassImp(AliMagFCM)
  55
  56 //_______________________________________________________________________
  57 AliMagFCM::AliMagFCM():
  58   fXbeg(0),
  59   fYbeg(0),
  60   fZbeg(0),
  61   fXdel(0),
  62   fYdel(0),
  63   fZdel(0),
  64   fSolenoid(0),
  65   fXdeli(0),
  66   fYdeli(0),
  67   fZdeli(0),
  68   fXn(0),
  69   fYn(0),
  70   fZn(0),
  71   fB(0)
  72 {
  73   //
  74   // Standard constructor
  75   //
  76   fType = kConMesh;
  77   fMap  = 2;
  78   SetSolenoidField();
  79 }
  80
  81 //_______________________________________________________________________
  82 AliMagFCM::AliMagFCM(const char *name, const char *title, const Int_t integ,
  83                      const Float_t factor, const Float_t fmax):
  84   AliMagF(name,title,integ,factor,fmax),
  85   fXbeg(0),
  86   fYbeg(0),
  87   fZbeg(0),
  88   fXdel(0),
  89   fYdel(0),
  90   fZdel(0),
  91   fSolenoid(0),
  92   fXdeli(0),
  93   fYdeli(0),
  94   fZdeli(0),
  95   fXn(0),
  96   fYn(0),
  97   fZn(0),
  98   fB(0)
  99 {
 100   //
 101   // Standard constructor
 102   //
 103   fType = kConMesh;
 104   fMap  = 2;
 105   SetSolenoidField();
 106
 107   if(fDebug>-1) Info("ctor",
 108      "%s: Constant Mesh Field %s created: map= %d, factor= %f, file= %s\n",
 109          ClassName(),fName.Data(), fMap, factor,fTitle.Data());
 110 }
 111
 112 //_______________________________________________________________________
 113 AliMagFCM::AliMagFCM(const AliMagFCM &magf):
 114   AliMagF(magf),
 115   fXbeg(0),
 116   fYbeg(0),
 117   fZbeg(0),
 118   fXdel(0),
 119   fYdel(0),
 120   fZdel(0),
 121   fSolenoid(0),
 122   fXdeli(0),
 123   fYdeli(0),
 124   fZdeli(0),
 125   fXn(0),
 126   fYn(0),
 127   fZn(0),
 128   fB(0)
 129 {
 130   //
 131   // Copy constructor
 132   //
 133   magf.Copy(*this);
 134 }
 135
 136 //_______________________________________________________________________
 137 void AliMagFCM::Field(Float_t *x, Float_t *b)
 138 {
 139   //
 140   // Method to calculate the magnetic field
 141   //
 142   Double_t ratx, raty, ratz, hix, hiy, hiz, ratx1, raty1, ratz1,
 143     bhyhz, bhylz, blyhz, blylz, bhz, blz, xl[3];
 144   const Double_t kone=1;
 145   Int_t ix, iy, iz;
 146
 147   // --- find the position in the grid ---
 148
 149   b[0]=b[1]=b[2]=0;
 150   if(-700<x[2] && x[2]<fZbeg && x[0]*x[0]+(x[1]+30)*(x[1]+30) < 560*560) {
 151     b[2]= fSolenoid;
 152   } else  {
 153     Bool_t infield=(fZbeg<=x[2] && x[2]<fZbeg+fZdel*(fZn-1)
 154                     &&  ( fXbeg <= TMath::Abs(x[0]) && TMath::Abs(x[0]) < fXbeg+fXdel*(fXn-1) )
 155                     &&  ( fYbeg <= TMath::Abs(x[1]) && TMath::Abs(x[1]) < fYbeg+fYdel*(fYn-1) ));
 156     if(infield) {
 157       xl[0]=TMath::Abs(x[0])-fXbeg;
 158       xl[1]=TMath::Abs(x[1])-fYbeg;
 159       xl[2]=x[2]-fZbeg;
 160
 161       // --- start with x
 162
 163       hix=xl[0]*fXdeli;
 164       ratx=hix-int(hix);
 165       ix=int(hix);
 166
 167       hiy=xl[1]*fYdeli;
 168       raty=hiy-int(hiy);
 169       iy=int(hiy);
 170
 171       hiz=xl[2]*fZdeli;
 172       ratz=hiz-int(hiz);
 173       iz=int(hiz);
 174
 175       if(fMap==2) {
 176         // ... simple interpolation
 177         ratx1=kone-ratx;
 178         raty1=kone-raty;
 179         ratz1=kone-ratz;
 180         bhyhz = Bx(ix  ,iy+1,iz+1)*ratx1+Bx(ix+1,iy+1,iz+1)*ratx;
 181         bhylz = Bx(ix  ,iy+1,iz  )*ratx1+Bx(ix+1,iy+1,iz  )*ratx;
 182         blyhz = Bx(ix  ,iy  ,iz+1)*ratx1+Bx(ix+1,iy  ,iz+1)*ratx;
 183         blylz = Bx(ix  ,iy  ,iz  )*ratx1+Bx(ix+1,iy  ,iz  )*ratx;
 184         bhz   = blyhz             *raty1+bhyhz             *raty;
 185         blz   = blylz             *raty1+bhylz             *raty;
 186         b[0]  = blz               *ratz1+bhz               *ratz;
 187         //
 188         bhyhz = By(ix  ,iy+1,iz+1)*ratx1+By(ix+1,iy+1,iz+1)*ratx;
 189         bhylz = By(ix  ,iy+1,iz  )*ratx1+By(ix+1,iy+1,iz  )*ratx;
 190         blyhz = By(ix  ,iy  ,iz+1)*ratx1+By(ix+1,iy  ,iz+1)*ratx;
 191         blylz = By(ix  ,iy  ,iz  )*ratx1+By(ix+1,iy  ,iz  )*ratx;
 192         bhz   = blyhz             *raty1+bhyhz             *raty;
 193         blz   = blylz             *raty1+bhylz             *raty;
 194         b[1]  = blz               *ratz1+bhz               *ratz;
 195         //
 196         bhyhz = Bz(ix  ,iy+1,iz+1)*ratx1+Bz(ix+1,iy+1,iz+1)*ratx;
 197         bhylz = Bz(ix  ,iy+1,iz  )*ratx1+Bz(ix+1,iy+1,iz  )*ratx;
 198         blyhz = Bz(ix  ,iy  ,iz+1)*ratx1+Bz(ix+1,iy  ,iz+1)*ratx;
 199         blylz = Bz(ix  ,iy  ,iz  )*ratx1+Bz(ix+1,iy  ,iz  )*ratx;
 200         bhz   = blyhz             *raty1+bhyhz             *raty;
 201         blz   = blylz             *raty1+bhylz             *raty;
 202         b[2]  = blz               *ratz1+bhz               *ratz;
 203         //printf("ratx,raty,ratz,b[0],b[1],b[2] %f %f %f %f %f %f\n",
 204         //ratx,raty,ratz,b[0],b[1],b[2]);
 205         //
 206         // ... use the dipole symmetry
 207         if (x[0]*x[1] < 0) b[1]=-b[1];
 208         if (x[0]<0) b[2]=-b[2];
 209       } else {
 210         printf("Invalid field map for constant mesh %d\n",fMap);
 211       }
 212     } else {
 213 //This is the ZDC part
 214     Float_t rad2=x[0]*x[0]+x[1]*x[1];
 215     if(x[2]>kCORBEG2 && x[2]<kCOREND2){
 216       if(rad2<kCOR2RA2){
 217         b[0] = kFCORN2;
 218       }
 219     }
 220     else if(x[2]>kZ1BEG && x[2]<kZ1END){
 221       if(rad2<kZ1RA2){
 222         b[0] = -kG1*x[1];
 223         b[1] = -kG1*x[0];
 224       }
 225     }
 226     else if(x[2]>kZ2BEG && x[2]<kZ2END){
 227       if(rad2<kZ2RA2){
 228         b[0] = kG1*x[1];
 229         b[1] = kG1*x[0];
 230       }
 231     }
 232     else if(x[2]>kZ3BEG && x[2]<kZ3END){
 233       if(rad2<kZ3RA2){
 234         b[0] = kG1*x[1];
 235         b[1] = kG1*x[0];
 236       }
 237     }
 238     else if(x[2]>kZ4BEG && x[2]<kZ4END){
 239       if(rad2<kZ4RA2){
 240         b[0] = -kG1*x[1];
 241         b[1] = -kG1*x[0];
 242       }
 243     }
 244     else if(x[2]>kD1BEG && x[2]<kD1END){
 245       if(rad2<kD1RA2){
 246         b[1] = -kFDIP;
 247       }
 248     }
 249     else if(x[2]>kD2BEG && x[2]<kD2END){
 250       if(((x[0]-kXCEN1D2)*(x[0]-kXCEN1D2)+(x[1]-kYCEN1D2)*(x[1]-kYCEN1D2))<kD2RA2
 251         || ((x[0]-kXCEN2D2)*(x[0]-kXCEN2D2)+(x[1]-kYCEN2D2)*(x[1]-kYCEN2D2))<kD2RA2){
 252         b[1] = kFDIP;
 253       }
 254     }
 255
 256     }
 257   }
 258   if(fFactor!=1) {
 259     b[0]*=fFactor;
 260     b[1]*=fFactor;
 261     b[2]*=fFactor;
 262   }
 263 }
 264
 265 //_______________________________________________________________________
 266 void AliMagFCM::ReadField()
 267 {
 268   //
 269   // Method to read the magnetic field map from file
 270   //
 271   FILE *magfile;
 272   Int_t ix, iy, iz, ipx, ipy, ipz;
 273   Float_t bx, by, bz;
 274   char *fname;
 275   if(fDebug) printf("%s: Reading Magnetic Field %s from file %s\n",ClassName(),fName.Data(),fTitle.Data());
 276   fname = gSystem->ExpandPathName(fTitle.Data());
 277   magfile=fopen(fname,"r");
 278   delete [] fname;
 279   if (magfile) {
 280     fscanf(magfile,"%d %d %d %f %f %f %f %f %f",
 281            &fXn, &fYn, &fZn, &fXdel, &fYdel, &fZdel, &fXbeg, &fYbeg, &fZbeg);
 282     if(fDebug>1) printf("%s: fXn %d, fYn %d, fZn %d, fXdel %f, fYdel %f, fZdel %f, fXbeg %f, fYbeg %f, fZbeg %f\n",
 283                         ClassName(),fXn, fYn, fZn, fXdel, fYdel, fZdel, fXbeg, fYbeg, fZbeg);
 284     fXdeli=1./fXdel;
 285     fYdeli=1./fYdel;
 286     fZdeli=1./fZdel;
 287     fB = new TVector(3*fXn*fYn*fZn);
 288     for (iz=0; iz<fZn; iz++) {
 289       ipz=iz*3*(fXn*fYn);
 290       for (iy=0; iy<fYn; iy++) {
 291         ipy=ipz+iy*3*fXn;
 292         for (ix=0; ix<fXn; ix++) {
 293           ipx=ipy+ix*3;
 294           fscanf(magfile,"%f %f %f",&bz,&by,&bx);
 295           (*fB)(ipx+2)=bz;
 296           (*fB)(ipx+1)=by;
 297           (*fB)(ipx  )=bx;
 298         }
 299       }
 300     }
 301   } else {
 302     printf("%s: File %s not found !\n",ClassName(),fTitle.Data());
 303     exit(1);
 304   }
 305 }
 306
 307 //_______________________________________________________________________
 308 void AliMagFCM::Copy(AliMagFCM & /* magf */) const
 309 {
 310   //
 311   // Copy *this onto magf -- Not implemented
 312   //
 313   Fatal("Copy","Not implemented!\n");
 314 }
 315