STEER/AliMagF.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 #include <TClass.h>
  18 #include <TFile.h>
  19 #include <TSystem.h>
  20
  21 #include "AliMagF.h"
  22 #include "AliMagWrapCheb.h"
  23 #include "AliLog.h"
  24
  25 ClassImp(AliMagF)
  26
  27 const Double_t AliMagF::fgkSol2DipZ    =  -700.;
  28
  29 //_______________________________________________________________________
  30 AliMagF::AliMagF():
  31   TVirtualMagField(),
  32   fMeasuredMap(0),
  33   fMapType(k5kG),
  34   fSolenoid(0),
  35   fBeamType(kNoBeamField),
  36   fBeamEnergy(0),
  37   //
  38   fInteg(0),
  39   fPrecInteg(0),
  40   fFactorSol(1.),
  41   fFactorDip(1.),
  42   fMax(15),
  43   fDipoleOFF(kFALSE),
  44   //
  45   fQuadGradient(0),
  46   fDipoleField(0),
  47   fCCorrField(0),
  48   fACorr1Field(0),
  49   fACorr2Field(0),
  50   fParNames("","")
  51 {
  52   // Default constructor
  53   //
  54 }
  55
  56 //_______________________________________________________________________
  57 AliMagF::AliMagF(const char *name, const char* title, Int_t integ,
  58                  Double_t factorSol, Double_t factorDip,
  59                  Double_t fmax, BMap_t maptype, const char* path,
  60                  BeamType_t bt, Double_t be):
  61   TVirtualMagField(name),
  62   fMeasuredMap(0),
  63   fMapType(maptype),
  64   fSolenoid(0),
  65   fBeamType(bt),
  66   fBeamEnergy(be),
  67   //
  68   fInteg(integ),
  69   fPrecInteg(1),
  70   fFactorSol(1.),
  71   fFactorDip(1.),
  72   fMax(fmax),
  73   fDipoleOFF(factorDip==0.),
  74   //
  75   fQuadGradient(0),
  76   fDipoleField(0),
  77   fCCorrField(0),
  78   fACorr1Field(0),
  79   fACorr2Field(0),
  80   fParNames("","")
  81 {
  82   // Initialize the field with Geant integration option "integ" and max field "fmax,
  83   // Impose scaling of parameterized L3 field by factorSol and of dipole by factorDip.
  84   // The "be" is the energy of the beam in GeV/nucleon
  85   //
  86   SetTitle(title);
  87   if(integ<0 || integ > 2) {
  88     AliWarning(Form("Invalid magnetic field flag: %5d; Helix tracking chosen instead",integ));
  89     fInteg = 2;
  90   }
  91   if (fInteg == 0) fPrecInteg = 0;
  92   //
  93   const char* parname = 0;
  94   //
  95   if (fMapType == k2kG) {
  96     fSolenoid = 2.;
  97     parname = fDipoleOFF ? "Sol12_Dip0_Hole":"Sol12_Dip6_Hole";
  98   } else if (fMapType == k5kG) {
  99     fSolenoid = 5.;
 100     parname = fDipoleOFF ? "Sol30_Dip0_Hole":"Sol30_Dip6_Hole";
 101   } else if (fMapType == k5kGUniform) {
 102     fSolenoid = 5.;
 103     parname = "Sol30_Dip6_Uniform";
 104   } else {
 105     AliFatal(Form("Unknown field identifier %d is requested\n",fMapType));
 106   }
 107   //
 108   SetDataFileName(path);
 109   SetParamName(parname);
 110   //
 111   SetFactorSol(factorSol);
 112   SetFactorDip(factorDip);
 113   LoadParameterization();
 114   InitMachineField(fBeamType,fBeamEnergy);
 115 }
 116
 117 //_______________________________________________________________________
 118 AliMagF::AliMagF(const AliMagF &src):
 119   TVirtualMagField(src),
 120   fMeasuredMap(0),
 121   fMapType(src.fMapType),
 122   fSolenoid(src.fSolenoid),
 123   fBeamType(src.fBeamType),
 124   fBeamEnergy(src.fBeamEnergy),
 125   fInteg(src.fInteg),
 126   fPrecInteg(src.fPrecInteg),
 127   fFactorSol(src.fFactorSol),
 128   fFactorDip(src.fFactorDip),
 129   fMax(src.fMax),
 130   fDipoleOFF(src.fDipoleOFF),
 131   fQuadGradient(src.fQuadGradient),
 132   fDipoleField(src.fDipoleField),
 133   fCCorrField(src.fCCorrField),
 134   fACorr1Field(src.fACorr1Field),
 135   fACorr2Field(src.fACorr2Field),
 136   fParNames(src.fParNames)
 137 {
 138   if (src.fMeasuredMap) fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
 139 }
 140
 141 //_______________________________________________________________________
 142 AliMagF::~AliMagF()
 143 {
 144   delete fMeasuredMap;
 145 }
 146
 147 //_______________________________________________________________________
 148 Bool_t AliMagF::LoadParameterization()
 149 {
 150   if (fMeasuredMap) {
 151     AliError(Form("Field data %s are already loaded from %s\n",GetParamName(),GetDataFileName()));
 152     return kTRUE;
 153   }
 154   //
 155   char* fname = gSystem->ExpandPathName(GetDataFileName());
 156   TFile* file = TFile::Open(fname);
 157   if (!file) {
 158     AliError(Form("Failed to open magnetic field data file %s\n",fname));
 159     return kFALSE;
 160   }
 161   //
 162   fMeasuredMap = dynamic_cast<AliMagWrapCheb*>(file->Get(GetParamName()));
 163   if (!fMeasuredMap) {
 164     AliError(Form("Did not find field %s in %s\n",GetParamName(),fname));
 165     return kFALSE;
 166   }
 167   file->Close();
 168   delete file;
 169   return kTRUE;
 170 }
 171
 172
 173 //_______________________________________________________________________
 174 void AliMagF::Field(const Double_t *xyz, Double_t *b)
 175 {
 176   // Method to calculate the field at point  xyz
 177   //
 178   //  b[0]=b[1]=b[2]=0.0;
 179   if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
 180     fMeasuredMap->Field(xyz,b);
 181     if (xyz[2]>fgkSol2DipZ || fDipoleOFF) for (int i=3;i--;) b[i] *= fFactorSol;
 182     else                                  for (int i=3;i--;) b[i] *= fFactorDip;
 183   }
 184   else MachineField(xyz, b);
 185   //
 186 }
 187
 188 //_______________________________________________________________________
 189 Double_t AliMagF::GetBz(const Double_t *xyz) const
 190 {
 191   // Method to calculate the field at point  xyz
 192   //
 193   if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
 194     double bz = fMeasuredMap->GetBz(xyz);
 195     return (xyz[2]>fgkSol2DipZ || fDipoleOFF) ? bz*fFactorSol : bz*fFactorDip;
 196   }
 197   else return 0.;
 198 }
 199
 200 //_______________________________________________________________________
 201 AliMagF& AliMagF::operator=(const AliMagF& src)
 202 {
 203   if (this != &src && src.fMeasuredMap) {
 204     if (fMeasuredMap) delete fMeasuredMap;
 205     fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
 206     SetName(src.GetName());
 207     fSolenoid    = src.fSolenoid;
 208     fBeamType    = src.fBeamType;
 209     fBeamEnergy  = src.fBeamEnergy;
 210     fInteg       = src.fInteg;
 211     fPrecInteg   = src.fPrecInteg;
 212     fFactorSol   = src.fFactorSol;
 213     fFactorDip   = src.fFactorDip;
 214     fMax         = src.fMax;
 215     fDipoleOFF   = src.fDipoleOFF;
 216     fParNames    = src.fParNames;
 217   }
 218   return *this;
 219 }
 220
 221 //_______________________________________________________________________
 222 void AliMagF::InitMachineField(BeamType_t btype, Double_t benergy)
 223 {
 224   if (btype==kNoBeamField || benergy<1.) {
 225     fQuadGradient = fDipoleField = fCCorrField = fACorr1Field = fACorr2Field = 0.;
 226     return;
 227   }
 228   //
 229   double rigScale = benergy/7000.;   // scale according to ratio of E/Enominal
 230   // for ions assume PbPb (with energy provided per nucleon) and account for A/Z
 231   if (btype == kBeamTypeAA) rigScale *= 208./82.;
 232   //
 233   fQuadGradient = 22.0002*rigScale;
 234   fDipoleField  = 37.8781*rigScale;
 235   //
 236   // SIDE C
 237   fCCorrField   = -9.6980;
 238   // SIDE A
 239   fACorr1Field  = -13.2247;
 240   fACorr2Field  =  11.7905;
 241   //
 242 }
 243
 244 //_______________________________________________________________________
 245 void AliMagF::MachineField(const Double_t *x, Double_t *b) const
 246 {
 247   // ---- This is the ZDC part
 248   // Compansators for Alice Muon Arm Dipole
 249   const Double_t kBComp1CZ = 1075., kBComp1hDZ = 260./2., kBComp1SqR = 4.0*4.0;
 250   const Double_t kBComp2CZ = 2049., kBComp2hDZ = 153./2., kBComp2SqR = 4.5*4.5;
 251   //
 252   const Double_t kTripQ1CZ = 2615., kTripQ1hDZ = 637./2., kTripQ1SqR = 3.5*3.5;
 253   const Double_t kTripQ2CZ = 3480., kTripQ2hDZ = 550./2., kTripQ2SqR = 3.5*3.5;
 254   const Double_t kTripQ3CZ = 4130., kTripQ3hDZ = 550./2., kTripQ3SqR = 3.5*3.5;
 255   const Double_t kTripQ4CZ = 5015., kTripQ4hDZ = 637./2., kTripQ4SqR = 3.5*3.5;
 256   //
 257   const Double_t kDip1CZ = 6310.8,  kDip1hDZ = 945./2., kDip1SqRC = 4.5*4.5, kDip1SqRA = 3.375*3.375;
 258   const Double_t kDip2CZ = 12640.3, kDip2hDZ = 945./2., kDip2SqRC = 4.5*4.5, kDip2SqRA = 3.75*3.75;
 259   const Double_t kDip2DXC = 9.7, kDip2DXA = 9.4;
 260   //
 261   double rad2 = x[0] * x[0] + x[1] * x[1];
 262   //
 263   b[0] = b[1] = b[2] = 0;
 264   //
 265   // SIDE C **************************************************
 266   if(x[2]<0.){
 267     if(TMath::Abs(x[2]+kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
 268       b[0] = fCCorrField*fFactorDip;
 269     }
 270     else if(TMath::Abs(x[2]+kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
 271       b[0] = fQuadGradient*x[1];
 272       b[1] = fQuadGradient*x[0];
 273     }
 274     else if(TMath::Abs(x[2]+kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
 275       b[0] = -fQuadGradient*x[1];
 276       b[1] = -fQuadGradient*x[0];
 277     }
 278     else if(TMath::Abs(x[2]+kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
 279       b[0] = -fQuadGradient*x[1];
 280       b[1] = -fQuadGradient*x[0];
 281     }
 282     else if(TMath::Abs(x[2]+kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
 283       b[0] = fQuadGradient*x[1];
 284       b[1] = fQuadGradient*x[0];
 285     }
 286     else if(TMath::Abs(x[2]+kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRC){
 287       b[1] = fDipoleField;
 288     }
 289     else if(TMath::Abs(x[2]+kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRC) {
 290       double dxabs = TMath::Abs(x[0])-kDip2DXC;
 291       if ( (dxabs*dxabs + x[1]*x[1])<kDip2SqRC) {
 292         b[1] = -fDipoleField;
 293       }
 294     }
 295   }
 296   //
 297   // SIDE A **************************************************
 298   else{
 299     if(TMath::Abs(x[2]-kBComp1CZ)<kBComp1hDZ && rad2 < kBComp1SqR) {
 300       // Compensator magnet at z = 1075 m
 301       b[0] = fACorr1Field*fFactorDip;
 302     }
 303     //
 304     if(TMath::Abs(x[2]-kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
 305       b[0] = fACorr2Field*fFactorDip;
 306     }
 307     else if(TMath::Abs(x[2]-kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
 308       b[0] = -fQuadGradient*x[1];
 309       b[1] = -fQuadGradient*x[0];
 310     }
 311     else if(TMath::Abs(x[2]-kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
 312       b[0] =  fQuadGradient*x[1];
 313       b[1] =  fQuadGradient*x[0];
 314     }
 315     else if(TMath::Abs(x[2]-kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
 316       b[0] =  fQuadGradient*x[1];
 317       b[1] =  fQuadGradient*x[0];
 318     }
 319     else if(TMath::Abs(x[2]-kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
 320       b[0] = -fQuadGradient*x[1];
 321       b[1] = -fQuadGradient*x[0];
 322     }
 323     else if(TMath::Abs(x[2]-kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRA){
 324       b[1] = -fDipoleField;
 325     }
 326     else if(TMath::Abs(x[2]-kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRA) {
 327       double dxabs = TMath::Abs(x[0])-kDip2DXA;
 328       if ( (dxabs*dxabs + x[1]*x[1])<kDip2SqRA) {
 329         b[1] = fDipoleField;
 330       }
 331     }
 332   }
 333   //
 334 }
 335
 336 //_______________________________________________________________________
 337 void AliMagF::GetTPCInt(const Double_t *xyz, Double_t *b) const
 338 {
 339   // Method to calculate the integral of magnetic integral from xyz to nearest cathode plane
 340   b[0]=b[1]=b[2]=0.0;
 341   if (fMeasuredMap) {
 342     fMeasuredMap->GetTPCInt(xyz,b);
 343     for (int i=3;i--;) b[i] *= fFactorSol;
 344   }
 345 }
 346
 347 //_______________________________________________________________________
 348 void AliMagF::GetTPCIntCyl(const Double_t *rphiz, Double_t *b) const
 349 {
 350   // Method to calculate the integral of magnetic integral from point to nearest cathode plane
 351   // in cylindrical coordiates ( -pi<phi<pi convention )
 352   b[0]=b[1]=b[2]=0.0;
 353   if (fMeasuredMap) {
 354     fMeasuredMap->GetTPCIntCyl(rphiz,b);
 355     for (int i=3;i--;) b[i] *= fFactorSol;
 356   }
 357 }