STEER/AliMagF.cxx

   1
   2 #include "AliMagF.h"
   3 #include "TSystem.h"
   4 #include <stdlib.h>
   5 #include <stdio.h>
   6
   7 //ZDC part -------------------------------------------------------------------
   8
   9   static const Float_t G1=20.03;
  10   static const Float_t FDIP=-37.34;
  11   static const Float_t FDIMU=6.;
  12   static const Float_t FCORN=11.72;
  13 //
  14 // ZBEG       Beginning of the inner triplet
  15 // D1BEG      Beginning of separator dipole 1
  16 // D2BEG      Beginning of separator dipole 2
  17 // CORBEG     Corrector dipole beginning (because of dimuon arm)
  18 //
  19   static const Float_t CORBEG=1920,COREND=CORBEG+190, CORRA2=4.5*4.5;
  20 //
  21   static const Float_t ZBEG=2300;
  22   static const Float_t Z1BEG=ZBEG+   0,Z1END=Z1BEG+630,Z1RA2=3.5*3.5;
  23   static const Float_t Z2BEG=ZBEG+ 880,Z2END=Z2BEG+550,Z2RA2=3.5*3.5;
  24   static const Float_t Z3BEG=ZBEG+1530,Z3END=Z3BEG+550,Z3RA2=3.5*3.5;
  25   static const Float_t Z4BEG=ZBEG+2430,Z4END=Z4BEG+630,Z4RA2=3.5*3.5;
  26   static const Float_t D1BEG=5843.5   ,D1END=D1BEG+945,D1RA2=4.5*4.5;
  27   static const Float_t D2BEG=12113.2  ,D2END=D2BEG+945,D2RA2=4.5*.5;
  28
  29 //ZDC part -------------------------------------------------------------------
  30
  31 ClassImp(AliMagF)
  32
  33 //________________________________________
  34 AliMagF::AliMagF(const char *name, const char *title, const Int_t integ, const Int_t map,
  35                  const Float_t factor, const Float_t fmax)
  36   : TNamed(name,title)
  37 {
  38   fMap = map;
  39   fType = Undef;
  40   fInteg = integ;
  41   fFactor = factor;
  42   fMax = fmax;
  43 }
  44
  45 //________________________________________
  46 void AliMagF::Field(Float_t*, Float_t *b)
  47 {
  48   printf("Undefined MagF Field called, returning 0\n");
  49   b[0]=b[1]=b[2]=0;
  50 }
  51
  52 ClassImp(AliMagFC)
  53
  54 //________________________________________
  55 AliMagFC::AliMagFC(const char *name, const char *title, const Int_t integ, const Int_t map,
  56                  const Float_t factor, const Float_t fmax)
  57   : AliMagF(name,title,integ,map,factor,fmax)
  58 {
  59   printf("Constant Field %s created: map= %d, factor= %f\n",fName.Data(),map,factor);
  60   fType = Const;
  61 }
  62
  63 //________________________________________
  64 void AliMagFC::Field(Float_t *x, Float_t *b)
  65 {
  66   b[0]=b[1]=b[2]=0;
  67   if(fMap==1) {
  68     if(TMath::Abs(x[2])<700 && x[0]*x[0]+(x[1]+30)*(x[1]+30) < 560*560) {
  69       b[2]=2;
  70     } else {
  71       if ( 725 <= x[2] && x[2] <= 1225 ) {
  72         Float_t dz = TMath::Abs(975-x[2])*0.01;
  73         b[0]=(1-0.1*dz*dz)*7;
  74       }
  75       else {
  76 //This is the ZDC part
  77         Float_t rad2=x[0]*x[0]+x[1]*x[1];
  78         if(rad2<D2RA2) {
  79           if(x[2]>D2BEG) {
  80
  81 //    Separator Dipole D2
  82             if(x[2]<D2END) b[1]=FDIP;
  83           } else if(x[2]>D1BEG) {
  84
  85 //    Separator Dipole D1
  86             if(x[2]<D1END) b[1]=-FDIP;
  87           }
  88           if(rad2<CORRA2) {
  89
  90 //    First quadrupole of inner triplet de-focussing in x-direction
  91 //    Inner triplet
  92             if(x[2]>Z4BEG) {
  93               if(x[2]<Z4END) {
  94
  95 //    2430 <-> 3060
  96                 b[0]=-G1*x[1];
  97                 b[1]=-G1*x[0];
  98               }
  99             } else if(x[2]>Z3BEG) {
 100               if(x[2]<Z3END) {
 101
 102 //    1530 <-> 2080
 103                 b[0]=G1*x[1];
 104                 b[1]=G1*x[0];
 105               }
 106             } else if(x[2]>Z2BEG) {
 107               if(x[2]<Z2END) {
 108
 109 //    890 <-> 1430
 110                 b[0]=G1*x[1];
 111                 b[1]=G1*x[0];
 112               }
 113             } else if(x[2]>Z1BEG) {
 114               if(x[2]<Z1END) {
 115
 116 //    0 <->  630
 117                 b[0]=-G1*x[1];
 118                 b[1]=-G1*x[0];
 119               }
 120             } else if(x[2]>CORBEG) {
 121               if(x[2]<COREND) {
 122 //    Corrector dipole (because of dimuon arm)
 123                 b[0]=FCORN;
 124               }
 125             }
 126           }
 127         }
 128       }
 129     }
 130   } else {
 131     printf("Invalid field map for constant field %d\n",fMap);
 132     exit(1);
 133   }
 134 }
 135
 136 ClassImp(AliMagFCM)
 137
 138 //________________________________________
 139 AliMagFCM::AliMagFCM(const char *name, const char *title, const Int_t integ, const Int_t map,
 140                  const Float_t factor, const Float_t fmax)
 141   : AliMagF(name,title,integ,map,factor,fmax)
 142 {
 143   fType = ConMesh;
 144   printf("Constant Mesh Field %s created: map= %d, factor= %f, file= %s\n",fName.Data(),map,factor,fTitle.Data());
 145 }
 146
 147 //________________________________________
 148 void AliMagFCM::Field(Float_t *x, Float_t *b)
 149 {
 150   Double_t ratx, raty, ratz, hix, hiy, hiz, ratx1, raty1, ratz1,
 151     bhyhz, bhylz, blyhz, blylz, bhz, blz, xl[3];
 152   const Double_t one=1;
 153   Int_t ix, iy, iz;
 154
 155   // --- find the position in the grid ---
 156
 157   b[0]=b[1]=b[2]=0;
 158   if(-700<x[2] && x[2]<fZbeg && x[0]*x[0]+(x[1]+30)*(x[1]+30) < 560*560) {
 159     b[2]=2;
 160   } else  {
 161     if(fZbeg<=x[2] && x[2]<fZbeg+fZdel*(fZn-1)
 162        &&  ( fXbeg <= TMath::Abs(x[0]) && TMath::Abs(x[0]) < fXbeg+fXdel*(fXn-1) )
 163        &&  ( fYbeg <= TMath::Abs(x[1]) && TMath::Abs(x[1]) < fYbeg+fYdel*(fYn-1) )) {
 164       xl[0]=TMath::Abs(x[0])-fXbeg;
 165       xl[1]=TMath::Abs(x[1])-fYbeg;
 166       xl[2]=x[2]-fZbeg;
 167
 168     // --- start with x
 169
 170       hix=xl[0]*fXdeli;
 171       ratx=hix-int(hix);
 172       ix=int(hix);
 173
 174       hiy=xl[1]*fYdeli;
 175       raty=hiy-int(hiy);
 176       iy=int(hiy);
 177
 178       hiz=xl[2]*fZdeli;
 179       ratz=hiz-int(hiz);
 180       iz=int(hiz);
 181
 182       if(fMap==2) {
 183       // ... simple interpolation
 184         ratx1=one-ratx;
 185         raty1=one-raty;
 186         ratz1=one-ratz;
 187         bhyhz = Bx(ix  ,iy+1,iz+1)*ratx1+Bx(ix+1,iy+1,iz+1)*ratx;
 188         bhylz = Bx(ix  ,iy+1,iz  )*ratx1+Bx(ix+1,iy+1,iz  )*ratx;
 189         blyhz = Bx(ix  ,iy  ,iz+1)*ratx1+Bx(ix+1,iy  ,iz+1)*ratx;
 190         blylz = Bx(ix  ,iy  ,iz  )*ratx1+Bx(ix+1,iy  ,iz  )*ratx;
 191         bhz   = blyhz             *raty1+bhyhz             *raty;
 192         blz   = blylz             *raty1+bhylz             *raty;
 193         b[0]  = blz               *ratz1+bhz               *ratz;
 194         //
 195         bhyhz = By(ix  ,iy+1,iz+1)*ratx1+By(ix+1,iy+1,iz+1)*ratx;
 196         bhylz = By(ix  ,iy+1,iz  )*ratx1+By(ix+1,iy+1,iz  )*ratx;
 197         blyhz = By(ix  ,iy  ,iz+1)*ratx1+By(ix+1,iy  ,iz+1)*ratx;
 198         blylz = By(ix  ,iy  ,iz  )*ratx1+By(ix+1,iy  ,iz  )*ratx;
 199         bhz   = blyhz             *raty1+bhyhz             *raty;
 200         blz   = blylz             *raty1+bhylz             *raty;
 201         b[1]  = blz               *ratz1+bhz               *ratz;
 202         //
 203         bhyhz = Bz(ix  ,iy+1,iz+1)*ratx1+Bz(ix+1,iy+1,iz+1)*ratx;
 204         bhylz = Bz(ix  ,iy+1,iz  )*ratx1+Bz(ix+1,iy+1,iz  )*ratx;
 205         blyhz = Bz(ix  ,iy  ,iz+1)*ratx1+Bz(ix+1,iy  ,iz+1)*ratx;
 206         blylz = Bz(ix  ,iy  ,iz  )*ratx1+Bz(ix+1,iy  ,iz  )*ratx;
 207         bhz   = blyhz             *raty1+bhyhz             *raty;
 208         blz   = blylz             *raty1+bhylz             *raty;
 209         b[2]  = blz               *ratz1+bhz               *ratz;
 210       //printf("ratx,raty,ratz,b[0],b[1],b[2] %f %f %f %f %f %f\n",
 211       //ratx,raty,ratz,b[0],b[1],b[2]);
 212       //
 213     // ... use the dipole symmetry
 214         if (x[0]*x[1] < 0) b[1]=-b[1];
 215         if (x[0]<0) b[2]=-b[2];
 216       } else {
 217         printf("Invalid field map for constant mesh %d\n",fMap);
 218       }
 219     } else {
 220 //This is the ZDC part
 221       Float_t rad2=x[0]*x[0]+x[1]*x[1];
 222       if(rad2<D2RA2) {
 223         if(x[2]>D2BEG) {
 224
 225 //    Separator Dipole D2
 226           if(x[2]<D2END) b[1]=FDIP;
 227         } else if(x[2]>D1BEG) {
 228
 229 //    Separator Dipole D1
 230           if(x[2]<D1END) b[1]=-FDIP;
 231         }
 232         if(rad2<CORRA2) {
 233
 234 //    First quadrupole of inner triplet de-focussing in x-direction
 235 //    Inner triplet
 236           if(x[2]>Z4BEG) {
 237             if(x[2]<Z4END) {
 238
 239 //    2430 <-> 3060
 240               b[0]=-G1*x[1];
 241               b[1]=-G1*x[0];
 242             }
 243           } else if(x[2]>Z3BEG) {
 244             if(x[2]<Z3END) {
 245
 246 //    1530 <-> 2080
 247               b[0]=G1*x[1];
 248               b[1]=G1*x[0];
 249             }
 250           } else if(x[2]>Z2BEG) {
 251             if(x[2]<Z2END) {
 252
 253 //    890 <-> 1430
 254               b[0]=G1*x[1];
 255               b[1]=G1*x[0];
 256             }
 257           } else if(x[2]>Z1BEG) {
 258             if(x[2]<Z1END) {
 259
 260 //    0 <->  630
 261               b[0]=-G1*x[1];
 262               b[1]=-G1*x[0];
 263             }
 264           } else if(x[2]>CORBEG) {
 265             if(x[2]<COREND) {
 266 //    Corrector dipole (because of dimuon arm)
 267               b[0]=FCORN;
 268             }
 269           }
 270         }
 271       }
 272     }
 273   }
 274 }
 275
 276 //________________________________________
 277 void AliMagFCM::ReadField()
 278 {
 279   FILE *magfile;
 280   Int_t ix, iy, iz, ipx, ipy, ipz;
 281   Float_t bx, by, bz;
 282   char *fname;
 283   printf("Reading Magnetic Field %s from file %s\n",fName.Data(),fTitle.Data());
 284   fname = gSystem->ExpandPathName(fTitle.Data());
 285   magfile=fopen(fname,"r");
 286   delete [] fname;
 287   if (magfile) {
 288     fscanf(magfile,"%d %d %d %f %f %f %f %f %f",
 289            &fXn, &fYn, &fZn, &fXdel, &fYdel, &fZdel, &fXbeg, &fYbeg, &fZbeg);
 290     printf("fXn %d, fYn %d, fZn %d, fXdel %f, fYdel %f, fZdel %f, fXbeg %f, fYbeg %f, fZbeg %f\n",
 291            fXn, fYn, fZn, fXdel, fYdel, fZdel, fXbeg, fYbeg, fZbeg);
 292     fXdeli=1./fXdel;
 293     fYdeli=1./fYdel;
 294     fZdeli=1./fZdel;
 295     fB = new TVector(3*fXn*fYn*fZn);
 296     for (iz=0; iz<fZn; iz++) {
 297       ipz=iz*3*(fXn*fYn);
 298       for (iy=0; iy<fYn; iy++) {
 299         ipy=ipz+iy*3*fXn;
 300         for (ix=0; ix<fXn; ix++) {
 301           ipx=ipy+ix*3;
 302           fscanf(magfile,"%f %f %f",&bz,&by,&bx);
 303           (*fB)(ipx+2)=bz;
 304           (*fB)(ipx+1)=by;
 305           (*fB)(ipx  )=bx;
 306         }
 307       }
 308     }
 309   } else {
 310     printf("File %s not found !\n",fTitle.Data());
 311     exit(1);
 312   }
 313 }
 314
 315
 316