STEER/AliMagFC.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 /* $Id$ */
  17
  18 //-------------------------------------------------------------------------
  19 //     Constant magnetic field class
  20 //     Used by AliRun class
  21 //     Author:
  22 //-------------------------------------------------------------------------
  23
  24 #include <stdlib.h>
  25
  26 #include "AliLog.h"
  27 #include "AliMagFC.h"
  28
  29 ClassImp(AliMagFC)
  30
  31 //________________________________________
  32 AliMagFC::AliMagFC()
  33     :AliMagF(),
  34     fCompensator(kFALSE),
  35     fBeamType(kBeamTypepp),
  36     fBeamEnergy(0),
  37     fQuadGradient(0),
  38     fDipoleField(0),
  39     fCCorrField(0),
  40     fACorr1Field(0),
  41     fACorr2Field(0)
  42 {
  43   //
  44   // Default constructor
  45   //
  46 }
  47
  48 //________________________________________
  49 AliMagFC::AliMagFC(const char *name, const char *title, Int_t integ,
  50                    Float_t factor, Float_t fmax)
  51     : AliMagF(name,title,integ,factor,fmax),
  52     fCompensator(kFALSE),
  53     fBeamType(kBeamTypepp),
  54     fBeamEnergy(7000.),
  55     fQuadGradient(0),
  56     fDipoleField(0),
  57     fCCorrField(0),
  58     fACorr1Field(0),
  59     fACorr2Field(0)
  60
  61 {
  62   //
  63   // Standard constructor
  64   //
  65   fType = kConst;
  66   fMap  = 1;
  67
  68 }
  69
  70 //________________________________________
  71 void AliMagFC::Field(const float *x, float *b) const
  72 {
  73   //
  74   // Method to return the field in a point
  75   //
  76   b[0]=b[1]=b[2]=0;
  77   if(fMap==1) {
  78     if(TMath::Abs(x[2])<700 && x[0]*x[0]+(x[1]+30)*(x[1]+30) < 560*560) {
  79       b[2]=2;
  80     }
  81     else {
  82       if(-725 >= x[2] && x[2] >= -1225 ){
  83         Float_t dz = TMath::Abs(-975-x[2])*0.01;
  84         b[0] = - (1-0.1*dz*dz)*7;
  85         if(fFactor!=1) {
  86             b[0]*=fFactor;
  87             b[1]*=fFactor;
  88             b[2]*=fFactor;
  89         }
  90       }
  91       else {
  92           ZDCField(x, b);
  93       }
  94     }
  95
  96   }
  97   else {
  98       AliFatal(Form("Invalid field map for constant field %d",fMap));
  99   }
 100 }
 101 //________________________________________
 102 void AliMagFC::Field(const double *x, double *b) const
 103 {
 104   //
 105   // Method to return the field in a point
 106   //
 107   b[0]=b[1]=b[2]=0;
 108   if(fMap==1) {
 109     if(TMath::Abs(x[2])<700 && x[0]*x[0]+(x[1]+30)*(x[1]+30) < 560*560) {
 110       b[2]=2;
 111     }
 112     else {
 113       if(-725 >= x[2] && x[2] >= -1225 ){
 114         Float_t dz = TMath::Abs(-975-x[2])*0.01;
 115         b[0] = - (1-0.1*dz*dz)*7;
 116         if(fFactor!=1) {
 117             b[0]*=fFactor;
 118             b[1]*=fFactor;
 119             b[2]*=fFactor;
 120         }
 121       }
 122       else {
 123           ZDCField(x, b);
 124       }
 125     }
 126
 127   }
 128   else {
 129       AliFatal(Form("Invalid field map for constant field %d",fMap));
 130   }
 131 }
 132
 133 //___________________________________________________
 134 void AliMagFC::ZDCField(const float *x, float *b) const
 135 {
 136   // ---- This is the ZDC part
 137
 138   float rad2 = x[0] * x[0] + x[1] * x[1];
 139   static Bool_t init = kFALSE;
 140
 141   if (! init) {
 142       init = kTRUE;
 143       //////////////////////////////////////////////////////////////////////
 144       // ---- Magnetic field values (according to beam type and energy) ----
 145       if(fBeamType==kBeamTypepp && fBeamEnergy == 5000.){
 146           // p-p @ 5+5 TeV
 147           fQuadGradient = 15.7145;
 148           fDipoleField  = 27.0558;
 149           // SIDE C
 150           fCCorrField   = 9.7017;
 151           // SIDE A
 152           fACorr1Field  = -13.2143;
 153           fACorr2Field  = -11.9909;
 154       } else if (fBeamType == kBeamTypepp && fBeamEnergy == 450.) {
 155           // p-p 0.45+0.45 TeV
 156           Float_t const kEnergyRatio = fBeamEnergy / 7000.;
 157
 158           fQuadGradient = 22.0002 * kEnergyRatio;
 159           fDipoleField  = 37.8781 * kEnergyRatio;
 160           // SIDE C
 161           fCCorrField   =  9.6908;
 162           // SIDE A
 163           fACorr1Field  = -13.2014;
 164           fACorr2Field  = -9.6908;
 165       } else if ((fBeamType == kBeamTypepp && fBeamEnergy == 7000.) ||
 166                  (fBeamType == kBeamTypeAA))
 167       {
 168           // Pb-Pb @ 2.7+2.7 TeV or p-p @ 7+7 TeV
 169           fQuadGradient = 22.0002;
 170           fDipoleField  = 37.8781;
 171           // SIDE C
 172           fCCorrField   = 9.6908;
 173           // SIDE A
 174           fACorr1Field  = -13.2014;
 175           fACorr2Field  = -9.6908;
 176       }
 177   }
 178
 179
 180   // SIDE C **************************************************
 181   if(x[2]<0.){
 182     if(x[2] < kCCorrBegin && x[2] > kCCorrEnd && rad2 < kCCorrSqRadius){
 183         if (fFactor != 0.) {
 184             b[0] = fCCorrField;
 185             b[1] = 0.;
 186             b[2] = 0.;
 187         }
 188     }
 189     else if(x[2] < kCQ1Begin && x[2] > kCQ1End && rad2 < kCQ1SqRadius){
 190         b[0] = fQuadGradient*x[1];
 191         b[1] = fQuadGradient*x[0];
 192         b[2] = 0.;
 193     }
 194     else if(x[2] < kCQ2Begin && x[2] > kCQ2End && rad2 < kCQ2SqRadius){
 195         b[0] = -fQuadGradient*x[1];
 196         b[1] = -fQuadGradient*x[0];
 197         b[2] = 0.;
 198     }
 199     else if(x[2] < kCQ3Begin && x[2] > kCQ3End && rad2 < kCQ3SqRadius){
 200         b[0] = -fQuadGradient*x[1];
 201         b[1] = -fQuadGradient*x[0];
 202         b[2] = 0.;
 203     }
 204     else if(x[2] < kCQ4Begin && x[2] > kCQ4End && rad2 < kCQ4SqRadius){
 205         b[0] = fQuadGradient*x[1];
 206         b[1] = fQuadGradient*x[0];
 207         b[2] = 0.;
 208     }
 209     else if(x[2] < kCD1Begin && x[2] > kCD1End && rad2 < kCD1SqRadius){
 210         b[1] = fDipoleField;
 211         b[2] = 0.;
 212         b[2] = 0.;
 213     }
 214     else if(x[2] < kCD2Begin && x[2] > kCD2End){
 215         if(((x[0]-kCD2XCentre1)*(x[0]-kCD2XCentre1)+(x[1]*x[1]))<kCD2SqRadius
 216            || ((x[0]-kCD2XCentre2)*(x[0]-kCD2XCentre2)+(x[1]*x[1]))<kCD2SqRadius){
 217           b[1] = -fDipoleField;
 218           b[2] = 0.;
 219           b[2] = 0.;
 220         }
 221     }
 222   }
 223
 224   // SIDE A **************************************************
 225   else{
 226     if(fCompensator && (x[2] > kACorr1Begin && x[2] < kACorr1End) && rad2 < kCCorr1SqRadius) {
 227       // Compensator magnet at z = 1075 m
 228         if (fFactor != 0.) {
 229             b[0] = fACorr1Field;
 230             b[1] = 0.;
 231             b[2] = 0.;
 232         }
 233         return;
 234     }
 235
 236     if(x[2] > kACorr2Begin && x[2] < kACorr2End && rad2 < kCCorr2SqRadius){
 237         if (fFactor != 0.) {
 238             b[0] = fACorr2Field;
 239             b[1] = 0.;
 240             b[2] = 0.;
 241         }
 242     }
 243     else if(x[2] > kAQ1Begin && x[2] < kAQ1End && rad2 < kAQ1SqRadius){
 244         // First quadrupole of inner triplet de-focussing in x-direction
 245         b[0] = -fQuadGradient*x[1];
 246         b[1] = -fQuadGradient*x[0];
 247         b[2] = 0.;
 248     }
 249     else if(x[2] > kAQ2Begin && x[2] < kAQ2End && rad2 < kAQ2SqRadius){
 250         b[0] = fQuadGradient*x[1];
 251         b[1] = fQuadGradient*x[0];
 252         b[2] = 0.;
 253     }
 254     else if(x[2] > kAQ3Begin && x[2] < kAQ3End && rad2 < kAQ3SqRadius){
 255         b[0] = fQuadGradient*x[1];
 256         b[1] = fQuadGradient*x[0];
 257         b[2] = 0.;
 258     }
 259     else if(x[2] > kAQ4Begin && x[2] < kAQ4End && rad2 < kAQ4SqRadius){
 260         b[0] = -fQuadGradient*x[1];
 261         b[1] = -fQuadGradient*x[0];
 262         b[2] = 0.;
 263     }
 264     else if(x[2] > kAD1Begin && x[2] < kAD1End && rad2 < kAD1SqRadius){
 265         b[0] = 0.;
 266         b[1] = -fDipoleField;
 267         b[2] = 0.;
 268     }
 269     else if(x[2] > kAD2Begin && x[2] < kAD2End){
 270         if(((x[0]-kAD2XCentre1)*(x[0]-kAD2XCentre1)+(x[1]*x[1])) < kAD2SqRadius
 271            || ((x[0]-kAD2XCentre2)*(x[0]-kAD2XCentre2)+(x[1]*x[1])) < kAD2SqRadius){
 272             b[1] = fDipoleField;
 273         }
 274     }
 275   }
 276 }
 277
 278 void AliMagFC::ZDCField(const double *x, double *b) const
 279 {
 280   // ---- This is the ZDC part
 281
 282   double rad2 = x[0] * x[0] + x[1] * x[1];
 283   static Bool_t init = kFALSE;
 284
 285   if (! init) {
 286       init = kTRUE;
 287       //////////////////////////////////////////////////////////////////////
 288       // ---- Magnetic field values (according to beam type and energy) ----
 289       if(fBeamType==kBeamTypepp && fBeamEnergy == 5000.){
 290           // p-p @ 5+5 TeV
 291           fQuadGradient = 15.7145;
 292           fDipoleField  = 27.0558;
 293           // SIDE C
 294           fCCorrField   = 9.7017;
 295           // SIDE A
 296           fACorr1Field  = -13.2143;
 297           fACorr2Field  = -11.9909;
 298       } else if (fBeamType == kBeamTypepp && fBeamEnergy == 450.) {
 299           // p-p 0.45+0.45 TeV
 300           Float_t const kEnergyRatio = fBeamEnergy / 7000.;
 301
 302           fQuadGradient = 22.0002 * kEnergyRatio;
 303           fDipoleField  = 37.8781 * kEnergyRatio;
 304           // SIDE C
 305           fCCorrField   =  9.6908;
 306           // SIDE A
 307           fACorr1Field  = -13.2014;
 308           fACorr2Field  = -9.6908;
 309       } else if ((fBeamType == kBeamTypepp && fBeamEnergy == 7000.) ||
 310                  (fBeamType == kBeamTypeAA))
 311       {
 312           // Pb-Pb @ 2.7+2.7 TeV or p-p @ 7+7 TeV
 313           fQuadGradient = 22.0002;
 314           fDipoleField  = 37.8781;
 315           // SIDE C
 316           fCCorrField   = 9.6908;
 317           // SIDE A
 318           fACorr1Field  = -13.2014;
 319           fACorr2Field  = -9.6908;
 320       }
 321   }
 322
 323
 324   // SIDE C **************************************************
 325   if(x[2]<0.){
 326     if(x[2] < kCCorrBegin && x[2] > kCCorrEnd && rad2 < kCCorrSqRadius){
 327         if (fFactor != 0.) {
 328             b[0] = fCCorrField;
 329             b[1] = 0.;
 330             b[2] = 0.;
 331         }
 332     }
 333     else if(x[2] < kCQ1Begin && x[2] > kCQ1End && rad2 < kCQ1SqRadius){
 334         b[0] = fQuadGradient*x[1];
 335         b[1] = fQuadGradient*x[0];
 336         b[2] = 0.;
 337     }
 338     else if(x[2] < kCQ2Begin && x[2] > kCQ2End && rad2 < kCQ2SqRadius){
 339         b[0] = -fQuadGradient*x[1];
 340         b[1] = -fQuadGradient*x[0];
 341         b[2] = 0.;
 342     }
 343     else if(x[2] < kCQ3Begin && x[2] > kCQ3End && rad2 < kCQ3SqRadius){
 344         b[0] = -fQuadGradient*x[1];
 345         b[1] = -fQuadGradient*x[0];
 346         b[2] = 0.;
 347     }
 348     else if(x[2] < kCQ4Begin && x[2] > kCQ4End && rad2 < kCQ4SqRadius){
 349         b[0] = fQuadGradient*x[1];
 350         b[1] = fQuadGradient*x[0];
 351         b[2] = 0.;
 352     }
 353     else if(x[2] < kCD1Begin && x[2] > kCD1End && rad2 < kCD1SqRadius){
 354         b[1] = fDipoleField;
 355         b[2] = 0.;
 356         b[2] = 0.;
 357     }
 358     else if(x[2] < kCD2Begin && x[2] > kCD2End){
 359         if(((x[0]-kCD2XCentre1)*(x[0]-kCD2XCentre1)+(x[1]*x[1]))<kCD2SqRadius
 360            || ((x[0]-kCD2XCentre2)*(x[0]-kCD2XCentre2)+(x[1]*x[1]))<kCD2SqRadius){
 361           b[1] = -fDipoleField;
 362           b[2] = 0.;
 363           b[2] = 0.;
 364         }
 365     }
 366   }
 367
 368   // SIDE A **************************************************
 369   else{
 370     if(fCompensator && (x[2] > kACorr1Begin && x[2] < kACorr1End) && rad2 < kCCorr1SqRadius) {
 371       // Compensator magnet at z = 1075 m
 372         if (fFactor != 0.) {
 373             b[0] = fACorr1Field;
 374             b[1] = 0.;
 375             b[2] = 0.;
 376         }
 377         return;
 378     }
 379
 380     if(x[2] > kACorr2Begin && x[2] < kACorr2End && rad2 < kCCorr2SqRadius){
 381         if (fFactor != 0.) {
 382             b[0] = fACorr2Field;
 383             b[1] = 0.;
 384             b[2] = 0.;
 385         }
 386     }
 387     else if(x[2] > kAQ1Begin && x[2] < kAQ1End && rad2 < kAQ1SqRadius){
 388         // First quadrupole of inner triplet de-focussing in x-direction
 389         b[0] = -fQuadGradient*x[1];
 390         b[1] = -fQuadGradient*x[0];
 391         b[2] = 0.;
 392     }
 393     else if(x[2] > kAQ2Begin && x[2] < kAQ2End && rad2 < kAQ2SqRadius){
 394         b[0] = fQuadGradient*x[1];
 395         b[1] = fQuadGradient*x[0];
 396         b[2] = 0.;
 397     }
 398     else if(x[2] > kAQ3Begin && x[2] < kAQ3End && rad2 < kAQ3SqRadius){
 399         b[0] = fQuadGradient*x[1];
 400         b[1] = fQuadGradient*x[0];
 401         b[2] = 0.;
 402     }
 403     else if(x[2] > kAQ4Begin && x[2] < kAQ4End && rad2 < kAQ4SqRadius){
 404         b[0] = -fQuadGradient*x[1];
 405         b[1] = -fQuadGradient*x[0];
 406         b[2] = 0.;
 407     }
 408     else if(x[2] > kAD1Begin && x[2] < kAD1End && rad2 < kAD1SqRadius){
 409         b[0] = 0.;
 410         b[1] = -fDipoleField;
 411         b[2] = 0.;
 412     }
 413     else if(x[2] > kAD2Begin && x[2] < kAD2End){
 414         if(((x[0]-kAD2XCentre1)*(x[0]-kAD2XCentre1)+(x[1]*x[1])) < kAD2SqRadius
 415            || ((x[0]-kAD2XCentre2)*(x[0]-kAD2XCentre2)+(x[1]*x[1])) < kAD2SqRadius){
 416             b[1] = fDipoleField;
 417         }
 418     }
 419   }
 420 }
 421