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) parname = fDipoleOFF ? "Sol12_Dip0_Hole":"Sol12_Dip6_Hole";
  96   else if (fMapType == k5kG) parname = fDipoleOFF ? "Sol30_Dip0_Hole":"Sol30_Dip6_Hole";
  97   else if (fMapType == k5kGUniform) parname = "Sol30_Dip6_Uniform";
  98   else AliFatal(Form("Unknown field identifier %d is requested\n",fMapType));
  99   //
 100   SetDataFileName(path);
 101   SetParamName(parname);
 102   //
 103   LoadParameterization();
 104   InitMachineField(fBeamType,fBeamEnergy);
 105   double xyz[3]={0.,0.,0.};
 106   fSolenoid = GetBz(xyz);
 107   SetFactorSol(factorSol);
 108   SetFactorDip(factorDip);
 109 }
 110
 111 //_______________________________________________________________________
 112 AliMagF::AliMagF(const AliMagF &src):
 113   TVirtualMagField(src),
 114   fMeasuredMap(0),
 115   fMapType(src.fMapType),
 116   fSolenoid(src.fSolenoid),
 117   fBeamType(src.fBeamType),
 118   fBeamEnergy(src.fBeamEnergy),
 119   fInteg(src.fInteg),
 120   fPrecInteg(src.fPrecInteg),
 121   fFactorSol(src.fFactorSol),
 122   fFactorDip(src.fFactorDip),
 123   fMax(src.fMax),
 124   fDipoleOFF(src.fDipoleOFF),
 125   fQuadGradient(src.fQuadGradient),
 126   fDipoleField(src.fDipoleField),
 127   fCCorrField(src.fCCorrField),
 128   fACorr1Field(src.fACorr1Field),
 129   fACorr2Field(src.fACorr2Field),
 130   fParNames(src.fParNames)
 131 {
 132   if (src.fMeasuredMap) fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
 133 }
 134
 135 //_______________________________________________________________________
 136 AliMagF::~AliMagF()
 137 {
 138   delete fMeasuredMap;
 139 }
 140
 141 //_______________________________________________________________________
 142 Bool_t AliMagF::LoadParameterization()
 143 {
 144   if (fMeasuredMap) {
 145     AliError(Form("Field data %s are already loaded from %s\n",GetParamName(),GetDataFileName()));
 146     return kTRUE;
 147   }
 148   //
 149   char* fname = gSystem->ExpandPathName(GetDataFileName());
 150   TFile* file = TFile::Open(fname);
 151   if (!file) {
 152     AliError(Form("Failed to open magnetic field data file %s\n",fname));
 153     return kFALSE;
 154   }
 155   //
 156   fMeasuredMap = dynamic_cast<AliMagWrapCheb*>(file->Get(GetParamName()));
 157   if (!fMeasuredMap) {
 158     AliError(Form("Did not find field %s in %s\n",GetParamName(),fname));
 159     return kFALSE;
 160   }
 161   file->Close();
 162   delete file;
 163   return kTRUE;
 164 }
 165
 166
 167 //_______________________________________________________________________
 168 void AliMagF::Field(const Double_t *xyz, Double_t *b)
 169 {
 170   // Method to calculate the field at point  xyz
 171   //
 172   //  b[0]=b[1]=b[2]=0.0;
 173   if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
 174     fMeasuredMap->Field(xyz,b);
 175     if (xyz[2]>fgkSol2DipZ || fDipoleOFF) for (int i=3;i--;) b[i] *= fFactorSol;
 176     else                                  for (int i=3;i--;) b[i] *= fFactorDip;
 177   }
 178   else MachineField(xyz, b);
 179   //
 180 }
 181
 182 //_______________________________________________________________________
 183 Double_t AliMagF::GetBz(const Double_t *xyz) const
 184 {
 185   // Method to calculate the field at point  xyz
 186   //
 187   if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
 188     double bz = fMeasuredMap->GetBz(xyz);
 189     return (xyz[2]>fgkSol2DipZ || fDipoleOFF) ? bz*fFactorSol : bz*fFactorDip;
 190   }
 191   else return 0.;
 192 }
 193
 194 //_______________________________________________________________________
 195 AliMagF& AliMagF::operator=(const AliMagF& src)
 196 {
 197   if (this != &src && src.fMeasuredMap) {
 198     if (fMeasuredMap) delete fMeasuredMap;
 199     fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
 200     SetName(src.GetName());
 201     fSolenoid    = src.fSolenoid;
 202     fBeamType    = src.fBeamType;
 203     fBeamEnergy  = src.fBeamEnergy;
 204     fInteg       = src.fInteg;
 205     fPrecInteg   = src.fPrecInteg;
 206     fFactorSol   = src.fFactorSol;
 207     fFactorDip   = src.fFactorDip;
 208     fMax         = src.fMax;
 209     fDipoleOFF   = src.fDipoleOFF;
 210     fParNames    = src.fParNames;
 211   }
 212   return *this;
 213 }
 214
 215 //_______________________________________________________________________
 216 void AliMagF::InitMachineField(BeamType_t btype, Double_t benergy)
 217 {
 218   if (btype==kNoBeamField || benergy<1.) {
 219     fQuadGradient = fDipoleField = fCCorrField = fACorr1Field = fACorr2Field = 0.;
 220     return;
 221   }
 222   //
 223   double rigScale = benergy/7000.;   // scale according to ratio of E/Enominal
 224   // for ions assume PbPb (with energy provided per nucleon) and account for A/Z
 225   if (btype == kBeamTypeAA) rigScale *= 208./82.;
 226   //
 227   fQuadGradient = 22.0002*rigScale;
 228   fDipoleField  = 37.8781*rigScale;
 229   //
 230   // SIDE C
 231   fCCorrField   = -9.6980;
 232   // SIDE A
 233   fACorr1Field  = -13.2247;
 234   fACorr2Field  =  11.7905;
 235   //
 236 }
 237
 238 //_______________________________________________________________________
 239 void AliMagF::MachineField(const Double_t *x, Double_t *b) const
 240 {
 241   // ---- This is the ZDC part
 242   // Compansators for Alice Muon Arm Dipole
 243   const Double_t kBComp1CZ = 1075., kBComp1hDZ = 260./2., kBComp1SqR = 4.0*4.0;
 244   const Double_t kBComp2CZ = 2049., kBComp2hDZ = 153./2., kBComp2SqR = 4.5*4.5;
 245   //
 246   const Double_t kTripQ1CZ = 2615., kTripQ1hDZ = 637./2., kTripQ1SqR = 3.5*3.5;
 247   const Double_t kTripQ2CZ = 3480., kTripQ2hDZ = 550./2., kTripQ2SqR = 3.5*3.5;
 248   const Double_t kTripQ3CZ = 4130., kTripQ3hDZ = 550./2., kTripQ3SqR = 3.5*3.5;
 249   const Double_t kTripQ4CZ = 5015., kTripQ4hDZ = 637./2., kTripQ4SqR = 3.5*3.5;
 250   //
 251   const Double_t kDip1CZ = 6310.8,  kDip1hDZ = 945./2., kDip1SqRC = 4.5*4.5, kDip1SqRA = 3.375*3.375;
 252   const Double_t kDip2CZ = 12640.3, kDip2hDZ = 945./2., kDip2SqRC = 4.5*4.5, kDip2SqRA = 3.75*3.75;
 253   const Double_t kDip2DXC = 9.7, kDip2DXA = 9.4;
 254   //
 255   double rad2 = x[0] * x[0] + x[1] * x[1];
 256   //
 257   b[0] = b[1] = b[2] = 0;
 258   //
 259   // SIDE C **************************************************
 260   if(x[2]<0.){
 261     if(TMath::Abs(x[2]+kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
 262       b[0] = fCCorrField*fFactorDip;
 263     }
 264     else if(TMath::Abs(x[2]+kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
 265       b[0] = fQuadGradient*x[1];
 266       b[1] = fQuadGradient*x[0];
 267     }
 268     else if(TMath::Abs(x[2]+kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
 269       b[0] = -fQuadGradient*x[1];
 270       b[1] = -fQuadGradient*x[0];
 271     }
 272     else if(TMath::Abs(x[2]+kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
 273       b[0] = -fQuadGradient*x[1];
 274       b[1] = -fQuadGradient*x[0];
 275     }
 276     else if(TMath::Abs(x[2]+kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
 277       b[0] = fQuadGradient*x[1];
 278       b[1] = fQuadGradient*x[0];
 279     }
 280     else if(TMath::Abs(x[2]+kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRC){
 281       b[1] = fDipoleField;
 282     }
 283     else if(TMath::Abs(x[2]+kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRC) {
 284       double dxabs = TMath::Abs(x[0])-kDip2DXC;
 285       if ( (dxabs*dxabs + x[1]*x[1])<kDip2SqRC) {
 286         b[1] = -fDipoleField;
 287       }
 288     }
 289   }
 290   //
 291   // SIDE A **************************************************
 292   else{
 293     if(TMath::Abs(x[2]-kBComp1CZ)<kBComp1hDZ && rad2 < kBComp1SqR) {
 294       // Compensator magnet at z = 1075 m
 295       b[0] = fACorr1Field*fFactorDip;
 296     }
 297     //
 298     if(TMath::Abs(x[2]-kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
 299       b[0] = fACorr2Field*fFactorDip;
 300     }
 301     else if(TMath::Abs(x[2]-kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
 302       b[0] = -fQuadGradient*x[1];
 303       b[1] = -fQuadGradient*x[0];
 304     }
 305     else if(TMath::Abs(x[2]-kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
 306       b[0] =  fQuadGradient*x[1];
 307       b[1] =  fQuadGradient*x[0];
 308     }
 309     else if(TMath::Abs(x[2]-kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
 310       b[0] =  fQuadGradient*x[1];
 311       b[1] =  fQuadGradient*x[0];
 312     }
 313     else if(TMath::Abs(x[2]-kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
 314       b[0] = -fQuadGradient*x[1];
 315       b[1] = -fQuadGradient*x[0];
 316     }
 317     else if(TMath::Abs(x[2]-kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRA){
 318       b[1] = -fDipoleField;
 319     }
 320     else if(TMath::Abs(x[2]-kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRA) {
 321       double dxabs = TMath::Abs(x[0])-kDip2DXA;
 322       if ( (dxabs*dxabs + x[1]*x[1])<kDip2SqRA) {
 323         b[1] = fDipoleField;
 324       }
 325     }
 326   }
 327   //
 328 }
 329
 330 //_______________________________________________________________________
 331 void AliMagF::GetTPCInt(const Double_t *xyz, Double_t *b) const
 332 {
 333   // Method to calculate the integral of magnetic integral from xyz to nearest cathode plane
 334   b[0]=b[1]=b[2]=0.0;
 335   if (fMeasuredMap) {
 336     fMeasuredMap->GetTPCInt(xyz,b);
 337     for (int i=3;i--;) b[i] *= fFactorSol;
 338   }
 339 }
 340
 341 //_______________________________________________________________________
 342 void AliMagF::GetTPCIntCyl(const Double_t *rphiz, Double_t *b) const
 343 {
 344   // Method to calculate the integral of magnetic integral from point to nearest cathode plane
 345   // in cylindrical coordiates ( -pi<phi<pi convention )
 346   b[0]=b[1]=b[2]=0.0;
 347   if (fMeasuredMap) {
 348     fMeasuredMap->GetTPCIntCyl(rphiz,b);
 349     for (int i=3;i--;) b[i] *= fFactorSol;
 350   }
 351 }