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   AliInfo(Form("Alice   B fields: Solenoid (%+.2f*)%.0f kG, Dipole %s (%+.2f) %s",
 128                factorSol,(fMapType==k5kG||fMapType==k5kGUniform)?5.:2.,
 129                fDipoleOFF ? "OFF":"ON",factorDip,fMapType==k5kGUniform?" |Constant Field!":""));
 130   AliInfo(Form("Machine B fields for %s beam (%.0f GeV): QGrad: %.4f Dipole: %.4f",
 131                bt==kBeamTypeAA ? "A-A":(bt==kBeamTypepp ? "p-p":"OFF"),be,fQuadGradient,fDipoleField));
 132 }
 133
 134 //_______________________________________________________________________
 135 AliMagF::AliMagF(const AliMagF &src):
 136   TVirtualMagField(src),
 137   fMeasuredMap(0),
 138   fMapType(src.fMapType),
 139   fSolenoid(src.fSolenoid),
 140   fBeamType(src.fBeamType),
 141   fBeamEnergy(src.fBeamEnergy),
 142   fInteg(src.fInteg),
 143   fPrecInteg(src.fPrecInteg),
 144   fFactorSol(src.fFactorSol),
 145   fFactorDip(src.fFactorDip),
 146   fMax(src.fMax),
 147   fDipoleOFF(src.fDipoleOFF),
 148   fQuadGradient(src.fQuadGradient),
 149   fDipoleField(src.fDipoleField),
 150   fCCorrField(src.fCCorrField),
 151   fACorr1Field(src.fACorr1Field),
 152   fACorr2Field(src.fACorr2Field),
 153   fParNames(src.fParNames)
 154 {
 155   if (src.fMeasuredMap) fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
 156 }
 157
 158 //_______________________________________________________________________
 159 AliMagF::~AliMagF()
 160 {
 161   delete fMeasuredMap;
 162 }
 163
 164 //_______________________________________________________________________
 165 Bool_t AliMagF::LoadParameterization()
 166 {
 167   if (fMeasuredMap) {
 168     AliError(Form("Field data %s are already loaded from %s\n",GetParamName(),GetDataFileName()));
 169     return kTRUE;
 170   }
 171   //
 172   char* fname = gSystem->ExpandPathName(GetDataFileName());
 173   TFile* file = TFile::Open(fname);
 174   if (!file) {
 175     AliError(Form("Failed to open magnetic field data file %s\n",fname));
 176     return kFALSE;
 177   }
 178   //
 179   fMeasuredMap = dynamic_cast<AliMagWrapCheb*>(file->Get(GetParamName()));
 180   if (!fMeasuredMap) {
 181     AliError(Form("Did not find field %s in %s\n",GetParamName(),fname));
 182     return kFALSE;
 183   }
 184   file->Close();
 185   delete file;
 186   return kTRUE;
 187 }
 188
 189
 190 //_______________________________________________________________________
 191 void AliMagF::Field(const Double_t *xyz, Double_t *b)
 192 {
 193   // Method to calculate the field at point  xyz
 194   //
 195   //  b[0]=b[1]=b[2]=0.0;
 196   if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
 197     fMeasuredMap->Field(xyz,b);
 198     if (xyz[2]>fgkSol2DipZ || fDipoleOFF) for (int i=3;i--;) b[i] *= fFactorSol;
 199     else                                  for (int i=3;i--;) b[i] *= fFactorDip;
 200   }
 201   else MachineField(xyz, b);
 202   //
 203 }
 204
 205 //_______________________________________________________________________
 206 Double_t AliMagF::GetBz(const Double_t *xyz) const
 207 {
 208   // Method to calculate the field at point  xyz
 209   //
 210   if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
 211     double bz = fMeasuredMap->GetBz(xyz);
 212     return (xyz[2]>fgkSol2DipZ || fDipoleOFF) ? bz*fFactorSol : bz*fFactorDip;
 213   }
 214   else return 0.;
 215 }
 216
 217 //_______________________________________________________________________
 218 AliMagF& AliMagF::operator=(const AliMagF& src)
 219 {
 220   if (this != &src && src.fMeasuredMap) {
 221     if (fMeasuredMap) delete fMeasuredMap;
 222     fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
 223     SetName(src.GetName());
 224     fSolenoid    = src.fSolenoid;
 225     fBeamType    = src.fBeamType;
 226     fBeamEnergy  = src.fBeamEnergy;
 227     fInteg       = src.fInteg;
 228     fPrecInteg   = src.fPrecInteg;
 229     fFactorSol   = src.fFactorSol;
 230     fFactorDip   = src.fFactorDip;
 231     fMax         = src.fMax;
 232     fDipoleOFF   = src.fDipoleOFF;
 233     fParNames    = src.fParNames;
 234   }
 235   return *this;
 236 }
 237
 238 //_______________________________________________________________________
 239 void AliMagF::InitMachineField(BeamType_t btype, Double_t benergy)
 240 {
 241   if (btype==kNoBeamField || benergy<1.) {
 242     fQuadGradient = fDipoleField = fCCorrField = fACorr1Field = fACorr2Field = 0.;
 243     return;
 244   }
 245   //
 246   double rigScale = benergy/7000.;   // scale according to ratio of E/Enominal
 247   // for ions assume PbPb (with energy provided per nucleon) and account for A/Z
 248   if (btype == kBeamTypeAA) rigScale *= 208./82.;
 249   //
 250   fQuadGradient = 22.0002*rigScale;
 251   fDipoleField  = 37.8781*rigScale;
 252   //
 253   // SIDE C
 254   fCCorrField   = -9.6980;
 255   // SIDE A
 256   fACorr1Field  = -13.2247;
 257   fACorr2Field  =  11.7905;
 258   //
 259 }
 260
 261 //_______________________________________________________________________
 262 void AliMagF::MachineField(const Double_t *x, Double_t *b) const
 263 {
 264   // ---- This is the ZDC part
 265   // Compansators for Alice Muon Arm Dipole
 266   const Double_t kBComp1CZ = 1075., kBComp1hDZ = 260./2., kBComp1SqR = 4.0*4.0;
 267   const Double_t kBComp2CZ = 2049., kBComp2hDZ = 153./2., kBComp2SqR = 4.5*4.5;
 268   //
 269   const Double_t kTripQ1CZ = 2615., kTripQ1hDZ = 637./2., kTripQ1SqR = 3.5*3.5;
 270   const Double_t kTripQ2CZ = 3480., kTripQ2hDZ = 550./2., kTripQ2SqR = 3.5*3.5;
 271   const Double_t kTripQ3CZ = 4130., kTripQ3hDZ = 550./2., kTripQ3SqR = 3.5*3.5;
 272   const Double_t kTripQ4CZ = 5015., kTripQ4hDZ = 637./2., kTripQ4SqR = 3.5*3.5;
 273   //
 274   const Double_t kDip1CZ = 6310.8,  kDip1hDZ = 945./2., kDip1SqRC = 4.5*4.5, kDip1SqRA = 3.375*3.375;
 275   const Double_t kDip2CZ = 12640.3, kDip2hDZ = 945./2., kDip2SqRC = 4.5*4.5, kDip2SqRA = 3.75*3.75;
 276   const Double_t kDip2DXC = 9.7, kDip2DXA = 9.4;
 277   //
 278   double rad2 = x[0] * x[0] + x[1] * x[1];
 279   //
 280   b[0] = b[1] = b[2] = 0;
 281   //
 282   // SIDE C **************************************************
 283   if(x[2]<0.){
 284     if(TMath::Abs(x[2]+kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
 285       b[0] = fCCorrField*fFactorDip;
 286     }
 287     else if(TMath::Abs(x[2]+kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
 288       b[0] = fQuadGradient*x[1];
 289       b[1] = fQuadGradient*x[0];
 290     }
 291     else if(TMath::Abs(x[2]+kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
 292       b[0] = -fQuadGradient*x[1];
 293       b[1] = -fQuadGradient*x[0];
 294     }
 295     else if(TMath::Abs(x[2]+kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
 296       b[0] = -fQuadGradient*x[1];
 297       b[1] = -fQuadGradient*x[0];
 298     }
 299     else if(TMath::Abs(x[2]+kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
 300       b[0] = fQuadGradient*x[1];
 301       b[1] = fQuadGradient*x[0];
 302     }
 303     else if(TMath::Abs(x[2]+kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRC){
 304       b[1] = fDipoleField;
 305     }
 306     else if(TMath::Abs(x[2]+kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRC) {
 307       double dxabs = TMath::Abs(x[0])-kDip2DXC;
 308       if ( (dxabs*dxabs + x[1]*x[1])<kDip2SqRC) {
 309         b[1] = -fDipoleField;
 310       }
 311     }
 312   }
 313   //
 314   // SIDE A **************************************************
 315   else{
 316     if(TMath::Abs(x[2]-kBComp1CZ)<kBComp1hDZ && rad2 < kBComp1SqR) {
 317       // Compensator magnet at z = 1075 m
 318       b[0] = fACorr1Field*fFactorDip;
 319     }
 320     //
 321     if(TMath::Abs(x[2]-kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
 322       b[0] = fACorr2Field*fFactorDip;
 323     }
 324     else if(TMath::Abs(x[2]-kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
 325       b[0] = -fQuadGradient*x[1];
 326       b[1] = -fQuadGradient*x[0];
 327     }
 328     else if(TMath::Abs(x[2]-kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
 329       b[0] =  fQuadGradient*x[1];
 330       b[1] =  fQuadGradient*x[0];
 331     }
 332     else if(TMath::Abs(x[2]-kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
 333       b[0] =  fQuadGradient*x[1];
 334       b[1] =  fQuadGradient*x[0];
 335     }
 336     else if(TMath::Abs(x[2]-kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
 337       b[0] = -fQuadGradient*x[1];
 338       b[1] = -fQuadGradient*x[0];
 339     }
 340     else if(TMath::Abs(x[2]-kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRA){
 341       b[1] = -fDipoleField;
 342     }
 343     else if(TMath::Abs(x[2]-kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRA) {
 344       double dxabs = TMath::Abs(x[0])-kDip2DXA;
 345       if ( (dxabs*dxabs + x[1]*x[1])<kDip2SqRA) {
 346         b[1] = fDipoleField;
 347       }
 348     }
 349   }
 350   //
 351 }
 352
 353 //_______________________________________________________________________
 354 void AliMagF::GetTPCInt(const Double_t *xyz, Double_t *b) const
 355 {
 356   // Method to calculate the integral of magnetic integral from xyz to nearest cathode plane
 357   b[0]=b[1]=b[2]=0.0;
 358   if (fMeasuredMap) {
 359     fMeasuredMap->GetTPCInt(xyz,b);
 360     for (int i=3;i--;) b[i] *= fFactorSol;
 361   }
 362 }
 363
 364 //_______________________________________________________________________
 365 void AliMagF::GetTPCIntCyl(const Double_t *rphiz, Double_t *b) const
 366 {
 367   // Method to calculate the integral of magnetic integral from point to nearest cathode plane
 368   // in cylindrical coordiates ( -pi<phi<pi convention )
 369   b[0]=b[1]=b[2]=0.0;
 370   if (fMeasuredMap) {
 371     fMeasuredMap->GetTPCIntCyl(rphiz,b);
 372     for (int i=3;i--;) b[i] *= fFactorSol;
 373   }
 374 }
 375
 376 //_______________________________________________________________________
 377 void AliMagF::SetFactorSol(Float_t fc)
 378 {
 379   // set the sign/scale of the current in the L3 according to fgkPolarityConvention
 380   switch (fgkPolarityConvention) {
 381   case kConvDCS2008: fFactorSol = -fc; break;
 382   case kConvLHC    : fFactorSol = -fc; break;
 383   default          : fFactorSol =  fc; break;  // case kConvMap2005: fFactorSol =  fc; break;
 384   }
 385 }
 386
 387 //_______________________________________________________________________
 388 void AliMagF::SetFactorDip(Float_t fc)
 389 {
 390   // set the sign*scale of the current in the Dipole according to fgkPolarityConvention
 391   switch (fgkPolarityConvention) {
 392   case kConvDCS2008: fFactorDip =  fc; break;
 393   case kConvLHC    : fFactorDip = -fc; break;
 394   default          : fFactorDip =  fc; break;  // case kConvMap2005: fFactorDip =  fc; break;
 395   }
 396 }
 397
 398 //_______________________________________________________________________
 399 Double_t AliMagF::GetFactorSol() const
 400 {
 401   // return the sign*scale of the current in the Dipole according to fgkPolarityConventionthe
 402   switch (fgkPolarityConvention) {
 403   case kConvDCS2008: return -fFactorSol;
 404   case kConvLHC    : return -fFactorSol;
 405   default          : return  fFactorSol;       //  case kConvMap2005: return  fFactorSol;
 406   }
 407 }
 408
 409 //_______________________________________________________________________
 410 Double_t AliMagF::GetFactorDip() const
 411 {
 412   // return the sign*scale of the current in the Dipole according to fgkPolarityConventionthe
 413   switch (fgkPolarityConvention) {
 414   case kConvDCS2008: return  fFactorDip;
 415   case kConvLHC    : return -fFactorDip;
 416   default          : return  fFactorDip;       //  case kConvMap2005: return  fFactorDip;
 417   }
 418 }