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