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 Double_t AliMagF::fgkBMachineZ1 =   919.;
  29 const Double_t AliMagF::fgkBMachineZ2 = -1972.;
  30
  31 //_______________________________________________________________________
  32 AliMagF::AliMagF():
  33   TVirtualMagField(),
  34   fMeasuredMap(0),
  35   fMapType(k5kG),
  36   fSolenoid(0),
  37   fBeamType(kNoBeamField),
  38   fBeamEnergy(0),
  39   fCompensator(kFALSE),
  40   //
  41   fInteg(0),
  42   fPrecInteg(0),
  43   fFactorSol(1.),
  44   fFactorDip(1.),
  45   fMax(15),
  46   fDipoleOFF(kFALSE),
  47   //
  48   fQuadGradient(0),
  49   fDipoleField(0),
  50   fCCorrField(0),
  51   fACorr1Field(0),
  52   fACorr2Field(0),
  53   fParNames("","")
  54 {
  55   // Default constructor
  56   //
  57 }
  58
  59 //_______________________________________________________________________
  60 AliMagF::AliMagF(const char *name, const char* title, Int_t integ,
  61                  Double_t factorSol, Double_t factorDip,
  62                  Double_t fmax, BMap_t maptype, const char* path,
  63                  BeamType_t bt, Double_t be, Bool_t compensator):
  64   TVirtualMagField(name),
  65   fMeasuredMap(0),
  66   fMapType(maptype),
  67   fSolenoid(0),
  68   fBeamType(bt),
  69   fBeamEnergy(be),
  70   fCompensator(compensator),
  71   //
  72   fInteg(integ),
  73   fPrecInteg(1),
  74   fFactorSol(1.),
  75   fFactorDip(1.),
  76   fMax(fmax),
  77   fDipoleOFF(factorDip==0.),
  78   //
  79   fQuadGradient(0),
  80   fDipoleField(0),
  81   fCCorrField(0),
  82   fACorr1Field(0),
  83   fACorr2Field(0),
  84   fParNames("","")
  85 {
  86   //
  87   SetTitle(title);
  88   if(integ<0 || integ > 2) {
  89     AliWarning(Form("Invalid magnetic field flag: %5d; Helix tracking chosen instead",integ));
  90     fInteg = 2;
  91   }
  92   if (fInteg == 0) fPrecInteg = 0;
  93   //
  94   const char* parname = 0;
  95   //
  96   if (fMapType == k2kG) {
  97     fSolenoid = 2.;
  98     parname = fDipoleOFF ? "Sol12_Dip0_Hole":"Sol12_Dip6_Hole";
  99   } else if (fMapType == k5kG) {
 100     fSolenoid = 5.;
 101     parname = fDipoleOFF ? "Sol30_Dip0_Hole":"Sol30_Dip6_Hole";
 102   } else if (fMapType == k5kGUniform) {
 103     fSolenoid = 5.;
 104     parname = "Sol30_Dip6_Uniform";
 105   } else {
 106     AliFatal(Form("Unknown field identifier %d is requested\n",fMapType));
 107   }
 108   //
 109   SetDataFileName(path);
 110   SetParamName(parname);
 111   //
 112   SetFactorSol(factorSol);
 113   SetFactorDip(factorDip);
 114   LoadParameterization();
 115   InitMachineField(fBeamType,fBeamEnergy);
 116 }
 117
 118 //_______________________________________________________________________
 119 AliMagF::AliMagF(const AliMagF &src):
 120   TVirtualMagField(src),
 121   fMeasuredMap(0),
 122   fMapType(src.fMapType),
 123   fSolenoid(src.fSolenoid),
 124   fBeamType(src.fBeamType),
 125   fBeamEnergy(src.fBeamEnergy),
 126   fCompensator(src.fCompensator),
 127   fInteg(src.fInteg),
 128   fPrecInteg(src.fPrecInteg),
 129   fFactorSol(src.fFactorSol),
 130   fFactorDip(src.fFactorDip),
 131   fMax(src.fMax),
 132   fDipoleOFF(src.fDipoleOFF),
 133   fQuadGradient(src.fQuadGradient),
 134   fDipoleField(src.fDipoleField),
 135   fCCorrField(src.fCCorrField),
 136   fACorr1Field(src.fACorr1Field),
 137   fACorr2Field(src.fACorr2Field),
 138   fParNames(src.fParNames)
 139 {
 140   if (src.fMeasuredMap) fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
 141 }
 142
 143 //_______________________________________________________________________
 144 AliMagF::~AliMagF()
 145 {
 146   delete fMeasuredMap;
 147 }
 148
 149 //_______________________________________________________________________
 150 Bool_t AliMagF::LoadParameterization()
 151 {
 152   if (fMeasuredMap) {
 153     AliError(Form("Field data %s are already loaded from %s\n",GetParamName(),GetDataFileName()));
 154     return kTRUE;
 155   }
 156   //
 157   char* fname = gSystem->ExpandPathName(GetDataFileName());
 158   TFile* file = TFile::Open(fname);
 159   if (!file) {
 160     AliError(Form("Failed to open magnetic field data file %s\n",fname));
 161     return kFALSE;
 162   }
 163   //
 164   fMeasuredMap = dynamic_cast<AliMagWrapCheb*>(file->Get(GetParamName()));
 165   if (!fMeasuredMap) {
 166     AliError(Form("Did not find field %s in %s\n",GetParamName(),fname));
 167     return kFALSE;
 168   }
 169   file->Close();
 170   delete file;
 171   return kTRUE;
 172 }
 173
 174
 175 //_______________________________________________________________________
 176 void AliMagF::Field(const Double_t *xyz, Double_t *b)
 177 {
 178   // Method to calculate the field at point  xyz
 179   //
 180   b[0]=b[1]=b[2]=0.0;
 181   if (xyz[2] > fgkBMachineZ1 || xyz[2] < fgkBMachineZ2) MachineField(xyz, b);
 182   else if (fMeasuredMap) {
 183     fMeasuredMap->Field(xyz,b);
 184     if (xyz[2]>fgkSol2DipZ || fDipoleOFF) for (int i=3;i--;) b[i] *= fFactorSol;
 185     else                                  for (int i=3;i--;) b[i] *= fFactorDip;
 186   }
 187   //
 188 }
 189
 190 //_______________________________________________________________________
 191 Double_t AliMagF::GetBz(const Double_t *xyz) const
 192 {
 193   // Method to calculate the field at point  xyz
 194   //
 195   if (xyz[2] <= fgkBMachineZ1 && xyz[2] >= fgkBMachineZ2) return fMeasuredMap->GetBz(xyz);
 196   else {
 197     double b[3] = {0,0,0};
 198     MachineField(xyz, b);
 199     return b[2];
 200   }
 201   //
 202 }
 203
 204 //_______________________________________________________________________
 205 AliMagF& AliMagF::operator=(const AliMagF& src)
 206 {
 207   if (this != &src && src.fMeasuredMap) {
 208     if (fMeasuredMap) delete fMeasuredMap;
 209     fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
 210     SetName(src.GetName());
 211     fSolenoid    = src.fSolenoid;
 212     fBeamType    = src.fBeamType;
 213     fBeamEnergy  = src.fBeamEnergy;
 214     fCompensator = src.fCompensator;
 215     fInteg       = src.fInteg;
 216     fPrecInteg   = src.fPrecInteg;
 217     fFactorSol   = src.fFactorSol;
 218     fFactorDip   = src.fFactorDip;
 219     fMax         = src.fMax;
 220     fDipoleOFF   = src.fDipoleOFF;
 221     fParNames    = src.fParNames;
 222   }
 223   return *this;
 224 }
 225
 226 //_______________________________________________________________________
 227 void AliMagF::InitMachineField(BeamType_t btype, Double_t benergy)
 228 {
 229   const double kToler = 0.1;
 230   if (btype==kNoBeamField) {
 231     fQuadGradient = fDipoleField = fCCorrField = fACorr1Field = fACorr2Field = 0.;
 232   }
 233   //
 234   else if (btype==kBeamTypepp && TMath::Abs(1.-benergy/5000.)<kToler ){
 235     // p-p @ 5+5 TeV
 236     fQuadGradient = 15.7145;
 237     fDipoleField  = 27.0558;
 238     // SIDE C
 239     fCCorrField   = 9.7017;
 240     // SIDE A
 241     fACorr1Field  = -13.2143;
 242     fACorr2Field  = -11.9909;
 243   } else if (btype == kBeamTypepp && TMath::Abs(1.-benergy/450.)<kToler) {
 244     // p-p 0.45+0.45 TeV
 245     Double_t const kEnergyRatio = benergy / 7000.;
 246     fQuadGradient = 22.0002 * kEnergyRatio;
 247     fDipoleField  = 37.8781 * kEnergyRatio;
 248     // SIDE C
 249     fCCorrField   =  9.6908;
 250     // SIDE A
 251     fACorr1Field  = -13.2014;
 252     fACorr2Field  = -9.6908;
 253   } else if ( (btype == kBeamTypepp && TMath::Abs(1.-benergy/7000.)<kToler) ||
 254               (fBeamType == kBeamTypeAA) ) {
 255     // Pb-Pb @ 2.7+2.7 TeV or p-p @ 7+7 TeV
 256     fQuadGradient = 22.0002;
 257     fDipoleField  = 37.8781;
 258     // SIDE C
 259     fCCorrField   = 9.6908;
 260     // SIDE A
 261     fACorr1Field  = -13.2014;
 262     fACorr2Field  = -9.6908;
 263   }
 264   //
 265 }
 266
 267 //_______________________________________________________________________
 268 void AliMagF::MachineField(const Double_t *x, Double_t *b) const
 269 {
 270   // ---- This is the ZDC part
 271   const Double_t kCCorrBegin = fgkBMachineZ2-0.5,kCCorrEnd = kCCorrBegin - 153., kCCorrSqRadius = 4.5*4.5;
 272   //
 273   const Double_t kCTripletBegin  = -2296.5;
 274   const Double_t kCQ1Begin = kCTripletBegin,        kCQ1End = kCQ1Begin-637., kCQ1SqRadius = 3.5*3.5;
 275   const Double_t kCQ2Begin = kCTripletBegin-908.5,  kCQ2End = kCQ2Begin-550., kCQ2SqRadius = 3.5*3.5;
 276   const Double_t kCQ3Begin = kCTripletBegin-1558.5, kCQ3End = kCQ3Begin-550., kCQ3SqRadius = 3.5*3.5;
 277   const Double_t kCQ4Begin = kCTripletBegin-2400.,  kCQ4End = kCQ4Begin-637., kCQ4SqRadius = 3.5*3.5;
 278   //
 279   const Double_t kCD1Begin = -5838.3,  kCD1End = kCD1Begin-945., kCD1SqRadius = 4.5*4.5;
 280   const Double_t kCD2Begin = -12167.8, kCD2End = kCD2Begin-945., kCD2SqRadius = 4.5*4.5;
 281   const Double_t kCD2XCentre1 = -9.7;
 282   const Double_t kCD2XCentre2 =  9.7;
 283   //
 284   // -> SIDE A
 285   // NB -> kACorr1Begin = 919. to be checked
 286   const Double_t kACorr1Begin = fgkBMachineZ1, kACorr1End = kACorr1Begin+260., kCCorr1SqRadius = 4.*4.;
 287   const Double_t kACorr2Begin = -fgkBMachineZ2 + 0.5, kACorr2End = kACorr2Begin+153., kCCorr2SqRadius = 4.5*4.5;
 288   const Double_t kATripletBegin  = 2296.5;
 289   const Double_t kAQ1Begin = kATripletBegin,    kAQ1End = kAQ1Begin+637., kAQ1SqRadius = 3.5*3.5;
 290   const Double_t kAQ2Begin = kATripletBegin+908.5,  kAQ2End = kAQ2Begin+550., kAQ2SqRadius = 3.5*3.5;
 291   const Double_t kAQ3Begin = kATripletBegin+1558.5, kAQ3End = kAQ3Begin+550., kAQ3SqRadius = 3.5*3.5;
 292   const Double_t kAQ4Begin = kATripletBegin+2400.,  kAQ4End = kAQ4Begin+637., kAQ4SqRadius = 3.5*3.5;
 293   //
 294   const Double_t kAD1Begin = 5838.3,  kAD1End = kAD1Begin+945., kAD1SqRadius = 3.375*3.375;
 295   const Double_t kAD2Begin = 12167.8, kAD2End = kAD2Begin+945., kAD2SqRadius = 3.75*3.75;
 296   const Double_t kAD2XCentre1 = -9.4;
 297   const Double_t kAD2XCentre2 =  9.4;
 298   //
 299   double rad2 = x[0] * x[0] + x[1] * x[1];
 300   //
 301   // SIDE C **************************************************
 302   if(x[2]<0.){
 303     if(x[2] < kCCorrBegin && x[2] > kCCorrEnd && rad2 < kCCorrSqRadius){
 304       b[0] = fCCorrField;
 305       b[1] = 0.;
 306       b[2] = 0.;
 307     }
 308     else if(x[2] < kCQ1Begin && x[2] > kCQ1End && rad2 < kCQ1SqRadius){
 309       b[0] = fQuadGradient*x[1];
 310       b[1] = fQuadGradient*x[0];
 311       b[2] = 0.;
 312     }
 313     else if(x[2] < kCQ2Begin && x[2] > kCQ2End && rad2 < kCQ2SqRadius){
 314       b[0] = -fQuadGradient*x[1];
 315       b[1] = -fQuadGradient*x[0];
 316       b[2] = 0.;
 317     }
 318     else if(x[2] < kCQ3Begin && x[2] > kCQ3End && rad2 < kCQ3SqRadius){
 319       b[0] = -fQuadGradient*x[1];
 320       b[1] = -fQuadGradient*x[0];
 321       b[2] = 0.;
 322     }
 323     else if(x[2] < kCQ4Begin && x[2] > kCQ4End && rad2 < kCQ4SqRadius){
 324       b[0] = fQuadGradient*x[1];
 325       b[1] = fQuadGradient*x[0];
 326       b[2] = 0.;
 327     }
 328     else if(x[2] < kCD1Begin && x[2] > kCD1End && rad2 < kCD1SqRadius){
 329       b[1] = fDipoleField;
 330       b[2] = 0.;
 331       b[2] = 0.;
 332     }
 333     else if(x[2] < kCD2Begin && x[2] > kCD2End){
 334       if(((x[0]-kCD2XCentre1)*(x[0]-kCD2XCentre1)+(x[1]*x[1]))<kCD2SqRadius
 335          || ((x[0]-kCD2XCentre2)*(x[0]-kCD2XCentre2)+(x[1]*x[1]))<kCD2SqRadius){
 336         b[1] = -fDipoleField;
 337         b[2] = 0.;
 338         b[2] = 0.;
 339       }
 340     }
 341   }
 342   //
 343   // SIDE A **************************************************
 344   else{
 345     if(fCompensator && (x[2] > kACorr1Begin && x[2] < kACorr1End) && rad2 < kCCorr1SqRadius) {
 346       // Compensator magnet at z = 1075 m
 347       b[0] = fACorr1Field;
 348       b[1] = 0.;
 349       b[2] = 0.;
 350       return;
 351     }
 352     //
 353     if(x[2] > kACorr2Begin && x[2] < kACorr2End && rad2 < kCCorr2SqRadius){
 354       b[0] = fACorr2Field;
 355       b[1] = 0.;
 356       b[2] = 0.;
 357     }
 358     else if(x[2] > kAQ1Begin && x[2] < kAQ1End && rad2 < kAQ1SqRadius){
 359       // First quadrupole of inner triplet de-focussing in x-direction
 360       b[0] = -fQuadGradient*x[1];
 361       b[1] = -fQuadGradient*x[0];
 362       b[2] = 0.;
 363     }
 364     else if(x[2] > kAQ2Begin && x[2] < kAQ2End && rad2 < kAQ2SqRadius){
 365       b[0] = fQuadGradient*x[1];
 366       b[1] = fQuadGradient*x[0];
 367       b[2] = 0.;
 368     }
 369     else if(x[2] > kAQ3Begin && x[2] < kAQ3End && rad2 < kAQ3SqRadius){
 370       b[0] = fQuadGradient*x[1];
 371       b[1] = fQuadGradient*x[0];
 372       b[2] = 0.;
 373     }
 374     else if(x[2] > kAQ4Begin && x[2] < kAQ4End && rad2 < kAQ4SqRadius){
 375       b[0] = -fQuadGradient*x[1];
 376       b[1] = -fQuadGradient*x[0];
 377       b[2] = 0.;
 378     }
 379     else if(x[2] > kAD1Begin && x[2] < kAD1End && rad2 < kAD1SqRadius){
 380       b[0] = 0.;
 381       b[1] = -fDipoleField;
 382       b[2] = 0.;
 383     }
 384     else if(x[2] > kAD2Begin && x[2] < kAD2End){
 385       if(((x[0]-kAD2XCentre1)*(x[0]-kAD2XCentre1)+(x[1]*x[1])) < kAD2SqRadius
 386          || ((x[0]-kAD2XCentre2)*(x[0]-kAD2XCentre2)+(x[1]*x[1])) < kAD2SqRadius){
 387         b[1] = fDipoleField;
 388       }
 389     }
 390   }
 391 }
 392
 393 //_______________________________________________________________________
 394 void AliMagF::GetTPCInt(const Double_t *xyz, Double_t *b) const
 395 {
 396   // Method to calculate the integral of magnetic integral from xyz to nearest cathode plane
 397   b[0]=b[1]=b[2]=0.0;
 398   if (fMeasuredMap) {
 399     fMeasuredMap->GetTPCInt(xyz,b);
 400     for (int i=3;i--;) b[i] *= fFactorSol;
 401   }
 402 }
 403
 404 //_______________________________________________________________________
 405 void AliMagF::GetTPCIntCyl(const Double_t *rphiz, Double_t *b) const
 406 {
 407   // Method to calculate the integral of magnetic integral from point to nearest cathode plane
 408   // in cylindrical coordiates ( -pi<phi<pi convention )
 409   b[0]=b[1]=b[2]=0.0;
 410   if (fMeasuredMap) {
 411     fMeasuredMap->GetTPCIntCyl(rphiz,b);
 412     for (int i=3;i--;) b[i] *= fFactorSol;
 413   }
 414 }