X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;ds=sidebyside;f=STEER%2FAliMagF.cxx;h=f7960fece85c5cf9392810c8ab378abb62a3721d;hb=8f6e7d1078fc857b9751f03ce74b3a45173dc2fc;hp=a3cf3e2353525d283721b5928b05b69054b6df63;hpb=ff66b122d1a51ad36bcdadf001d1c1cf619bdfda;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliMagF.cxx b/STEER/AliMagF.cxx index a3cf3e23535..f7960fece85 100644 --- a/STEER/AliMagF.cxx +++ b/STEER/AliMagF.cxx @@ -13,142 +13,507 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* $Id$ */ -//---------------------------------------------------------------------- -// Basic magnetic field class -// Used in all the detectors, and also in the traking classes -// Author: -//---------------------------------------------------------------------- +#include +#include +#include +#include -#include "AliLog.h" #include "AliMagF.h" +#include "AliMagWrapCheb.h" +#include "AliLog.h" ClassImp(AliMagF) +const Double_t AliMagF::fgkSol2DipZ = -700.; +const UShort_t AliMagF::fgkPolarityConvention = kConvLHC; + +/* + Explanation for polarity conventions: these are the mapping between the + current signs and main field components in L3 (Bz) and Dipole (Bx) (in Alice frame) + 1) kConvMap2005: used for the field mapping in 2005 + positive L3 current -> negative Bz + positive Dip current -> positive Bx + 2) kConvMapDCS2008: defined by the microswitches/cabling of power converters as of 2008 - 1st half 2009 + positive L3 current -> positive Bz + positive Dip current -> positive Bx + 3) kConvLHC : defined by LHC + positive L3 current -> positive Bz + positive Dip current -> negative Bx + + Note: only "negative Bz(L3) with postive Bx(Dipole)" and its inverse was mapped in 2005. Hence + the GRP Manager will reject the runs with the current combinations (in the convention defined by the + static Int_t AliMagF::GetPolarityConvention()) which do not lead to such field polarities. +*/ //_______________________________________________________________________ AliMagF::AliMagF(): - fMap(0), - fType(0), + TVirtualMagField(), + fMeasuredMap(0), + fMapType(k5kG), + fSolenoid(0), + fBeamType(kNoBeamField), + fBeamEnergy(0), + // fInteg(0), - fPrecInteg(1), - fFactor(0), - fMax(0), - fReadField(1) -{ + fPrecInteg(0), + fFactorSol(1.), + fFactorDip(1.), + fMax(15), + fDipoleOFF(kFALSE), // + fQuadGradient(0), + fDipoleField(0), + fCCorrField(0), + fACorr1Field(0), + fACorr2Field(0), + fParNames("","") +{ // Default constructor // } //_______________________________________________________________________ -AliMagF::AliMagF(const char *name, const char *title, Int_t integ, - Float_t factor, Float_t fmax): - TNamed(name,title), - fMap(0), - fType(0), - fInteg(0), +AliMagF::AliMagF(const char *name, const char* title, Double_t factorSol, Double_t factorDip, + BMap_t maptype, BeamType_t bt, Double_t be,Int_t integ, Double_t fmax, const char* path): + TVirtualMagField(name), + fMeasuredMap(0), + fMapType(maptype), + fSolenoid(0), + fBeamType(bt), + fBeamEnergy(be), + // + fInteg(integ), fPrecInteg(1), - fFactor(factor), + fFactorSol(1.), + fFactorDip(1.), fMax(fmax), - fReadField(1) + fDipoleOFF(factorDip==0.), + // + fQuadGradient(0), + fDipoleField(0), + fCCorrField(0), + fACorr1Field(0), + fACorr2Field(0), + fParNames("","") { + // Initialize the field with Geant integration option "integ" and max field "fmax, + // Impose scaling of parameterized L3 field by factorSol and of dipole by factorDip. + // The "be" is the energy of the beam in GeV/nucleon // - // Standard constructor + SetTitle(title); + if(integ<0 || integ > 2) { + AliWarning(Form("Invalid magnetic field flag: %5d; Helix tracking chosen instead",integ)); + fInteg = 2; + } + if (fInteg == 0) fPrecInteg = 0; // - if(integ<0 || integ > 2) { - AliWarning(Form( - "Invalid magnetic field flag: %5d; Helix tracking chosen instead" - ,integ)); - fInteg = 2; - } else { - fInteg = integ; - } - - if (fInteg == 0) fPrecInteg = 0; - - fType = kUndef; - // + if (fBeamEnergy<=0 && fBeamType!=kNoBeamField) { + if (fBeamType == kBeamTypepp) fBeamEnergy = 7000.; // max proton energy + else if (fBeamType == kBeamTypeAA) fBeamEnergy = 2750; // max PbPb energy + AliInfo("Maximim possible beam energy for requested beam is assumed"); + } + const char* parname = 0; + // + if (fMapType == k2kG) parname = fDipoleOFF ? "Sol12_Dip0_Hole":"Sol12_Dip6_Hole"; + else if (fMapType == k5kG) parname = fDipoleOFF ? "Sol30_Dip0_Hole":"Sol30_Dip6_Hole"; + else if (fMapType == k5kGUniform) parname = "Sol30_Dip6_Uniform"; + else AliFatal(Form("Unknown field identifier %d is requested\n",fMapType)); + // + SetDataFileName(path); + SetParamName(parname); + // + LoadParameterization(); + InitMachineField(fBeamType,fBeamEnergy); + double xyz[3]={0.,0.,0.}; + fSolenoid = GetBz(xyz); + SetFactorSol(factorSol); + SetFactorDip(factorDip); + AliInfo(Form("Alice B fields: Solenoid (%+.2f*)%.0f kG, Dipole %s (%+.2f) %s", + factorSol,(fMapType==k5kG||fMapType==k5kGUniform)?5.:2., + fDipoleOFF ? "OFF":"ON",factorDip,fMapType==k5kGUniform?" |Constant Field!":"")); + AliInfo(Form("Machine B fields for %s beam (%.0f GeV): QGrad: %.4f Dipole: %.4f", + bt==kBeamTypeAA ? "A-A":(bt==kBeamTypepp ? "p-p":"OFF"),be,fQuadGradient,fDipoleField)); } //_______________________________________________________________________ AliMagF::AliMagF(const AliMagF &src): - TNamed(src), - fMap(src.fMap), - fType(src.fType), + TVirtualMagField(src), + fMeasuredMap(0), + fMapType(src.fMapType), + fSolenoid(src.fSolenoid), + fBeamType(src.fBeamType), + fBeamEnergy(src.fBeamEnergy), fInteg(src.fInteg), fPrecInteg(src.fPrecInteg), - fFactor(src.fFactor), + fFactorSol(src.fFactorSol), + fFactorDip(src.fFactorDip), fMax(src.fMax), - fReadField(src.fReadField) + fDipoleOFF(src.fDipoleOFF), + fQuadGradient(src.fQuadGradient), + fDipoleField(src.fDipoleField), + fCCorrField(src.fCCorrField), + fACorr1Field(src.fACorr1Field), + fACorr2Field(src.fACorr2Field), + fParNames(src.fParNames) +{ + if (src.fMeasuredMap) fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap); +} + +//_______________________________________________________________________ +AliMagF::~AliMagF() +{ + delete fMeasuredMap; +} + +//_______________________________________________________________________ +Bool_t AliMagF::LoadParameterization() { - // Copy constructor + if (fMeasuredMap) { + AliError(Form("Field data %s are already loaded from %s\n",GetParamName(),GetDataFileName())); + return kTRUE; + } + // + char* fname = gSystem->ExpandPathName(GetDataFileName()); + TFile* file = TFile::Open(fname); + if (!file) { + AliError(Form("Failed to open magnetic field data file %s\n",fname)); + return kFALSE; + } + // + fMeasuredMap = dynamic_cast(file->Get(GetParamName())); + if (!fMeasuredMap) { + AliError(Form("Did not find field %s in %s\n",GetParamName(),fname)); + return kFALSE; + } + file->Close(); + delete file; + return kTRUE; +} + + +//_______________________________________________________________________ +void AliMagF::Field(const Double_t *xyz, Double_t *b) +{ + // Method to calculate the field at point xyz + // + // b[0]=b[1]=b[2]=0.0; + if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]GetMaxZ()) { + fMeasuredMap->Field(xyz,b); + if (xyz[2]>fgkSol2DipZ || fDipoleOFF) for (int i=3;i--;) b[i] *= fFactorSol; + else for (int i=3;i--;) b[i] *= fFactorDip; + } + else MachineField(xyz, b); + // +} + +//_______________________________________________________________________ +Double_t AliMagF::GetBz(const Double_t *xyz) const +{ + // Method to calculate the field at point xyz + // + if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]GetMaxZ()) { + double bz = fMeasuredMap->GetBz(xyz); + return (xyz[2]>fgkSol2DipZ || fDipoleOFF) ? bz*fFactorSol : bz*fFactorDip; + } + else return 0.; } //_______________________________________________________________________ -void AliMagF::Field(float*, float *b) const +AliMagF& AliMagF::operator=(const AliMagF& src) { + if (this != &src && src.fMeasuredMap) { + if (fMeasuredMap) delete fMeasuredMap; + fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap); + SetName(src.GetName()); + fSolenoid = src.fSolenoid; + fBeamType = src.fBeamType; + fBeamEnergy = src.fBeamEnergy; + fInteg = src.fInteg; + fPrecInteg = src.fPrecInteg; + fFactorSol = src.fFactorSol; + fFactorDip = src.fFactorDip; + fMax = src.fMax; + fDipoleOFF = src.fDipoleOFF; + fParNames = src.fParNames; + } + return *this; +} + +//_______________________________________________________________________ +void AliMagF::InitMachineField(BeamType_t btype, Double_t benergy) +{ + if (btype==kNoBeamField) { + fQuadGradient = fDipoleField = fCCorrField = fACorr1Field = fACorr2Field = 0.; + return; + } + // + double rigScale = benergy/7000.; // scale according to ratio of E/Enominal + // for ions assume PbPb (with energy provided per nucleon) and account for A/Z + if (btype == kBeamTypeAA) rigScale *= 208./82.; + // + fQuadGradient = 22.0002*rigScale; + fDipoleField = 37.8781*rigScale; // - // Method to return the field in one point -- dummy in this case + // SIDE C + fCCorrField = -9.6980; + // SIDE A + fACorr1Field = -13.2247; + fACorr2Field = 11.7905; // - AliWarning("Undefined MagF Field called, returning 0"); - b[0]=b[1]=b[2]=0; } //_______________________________________________________________________ -void AliMagF::Field(double*, double *b) const +void AliMagF::MachineField(const Double_t *x, Double_t *b) const { + // ---- This is the ZDC part + // Compansators for Alice Muon Arm Dipole + const Double_t kBComp1CZ = 1075., kBComp1hDZ = 260./2., kBComp1SqR = 4.0*4.0; + const Double_t kBComp2CZ = 2049., kBComp2hDZ = 153./2., kBComp2SqR = 4.5*4.5; + // + const Double_t kTripQ1CZ = 2615., kTripQ1hDZ = 637./2., kTripQ1SqR = 3.5*3.5; + const Double_t kTripQ2CZ = 3480., kTripQ2hDZ = 550./2., kTripQ2SqR = 3.5*3.5; + const Double_t kTripQ3CZ = 4130., kTripQ3hDZ = 550./2., kTripQ3SqR = 3.5*3.5; + const Double_t kTripQ4CZ = 5015., kTripQ4hDZ = 637./2., kTripQ4SqR = 3.5*3.5; + // + const Double_t kDip1CZ = 6310.8, kDip1hDZ = 945./2., kDip1SqRC = 4.5*4.5, kDip1SqRA = 3.375*3.375; + const Double_t kDip2CZ = 12640.3, kDip2hDZ = 945./2., kDip2SqRC = 4.5*4.5, kDip2SqRA = 3.75*3.75; + const Double_t kDip2DXC = 9.7, kDip2DXA = 9.4; + // + double rad2 = x[0] * x[0] + x[1] * x[1]; // - // Method to return the field in one point -- dummy in this case + b[0] = b[1] = b[2] = 0; // - AliWarning("Undefined MagF Field called, returning 0"); - b[0]=b[1]=b[2]=0; + // SIDE C ************************************************** + if(x[2]<0.){ + if(TMath::Abs(x[2]+kBComp2CZ)GetTPCInt(xyz,b); + for (int i=3;i--;) b[i] *= fFactorSol; + } +} + +//_______________________________________________________________________ +void AliMagF::GetTPCIntCyl(const Double_t *rphiz, Double_t *b) const +{ + // Method to calculate the integral of magnetic integral from point to nearest cathode plane + // in cylindrical coordiates ( -piGetTPCIntCyl(rphiz,b); + for (int i=3;i--;) b[i] *= fFactorSol; + } } //_______________________________________________________________________ -void AliMagF::GetTPCInt(Float_t *, Float_t *b) const +void AliMagF::SetFactorSol(Float_t fc) { -// -// Obtain the integral of the field components in the TPC from given point -// to the closest cathod plane -// - AliWarning("Undefined MagF TPCIntegral called, returning 0"); - b[0]=b[1]=b[2]=0; + // set the sign/scale of the current in the L3 according to fgkPolarityConvention + switch (fgkPolarityConvention) { + case kConvDCS2008: fFactorSol = -fc; break; + case kConvLHC : fFactorSol = -fc; break; + default : fFactorSol = fc; break; // case kConvMap2005: fFactorSol = fc; break; + } } //_______________________________________________________________________ -void AliMagF::GetTPCIntCyl(Float_t *, Float_t *b) const +void AliMagF::SetFactorDip(Float_t fc) { -// -// Obtain the integral of the field components in the TPC from given point -// to the closest cathod plane -// - AliWarning("Undefined MagF TPCIntegral called, returning 0"); - b[0]=b[1]=b[2]=0; + // set the sign*scale of the current in the Dipole according to fgkPolarityConvention + switch (fgkPolarityConvention) { + case kConvDCS2008: fFactorDip = fc; break; + case kConvLHC : fFactorDip = -fc; break; + default : fFactorDip = fc; break; // case kConvMap2005: fFactorDip = fc; break; + } } //_______________________________________________________________________ -AliMagF& AliMagF::operator=(const AliMagF& rhs) +Double_t AliMagF::GetFactorSol() const { - // Asignment operator - fMap = rhs.fMap; - fType = rhs.fType; - fInteg = rhs.fInteg; - fPrecInteg = rhs.fPrecInteg; - fFactor = rhs.fFactor; - fMax = rhs.fMax; - fReadField = rhs.fReadField; - return *this; + // return the sign*scale of the current in the Dipole according to fgkPolarityConventionthe + switch (fgkPolarityConvention) { + case kConvDCS2008: return -fFactorSol; + case kConvLHC : return -fFactorSol; + default : return fFactorSol; // case kConvMap2005: return fFactorSol; + } } -void AliMagF::SetPrecInteg(Int_t integ) +//_______________________________________________________________________ +Double_t AliMagF::GetFactorDip() const { - if (fInteg > 0) { - fPrecInteg = integ; + // return the sign*scale of the current in the Dipole according to fgkPolarityConventionthe + switch (fgkPolarityConvention) { + case kConvDCS2008: return fFactorDip; + case kConvLHC : return -fFactorDip; + default : return fFactorDip; // case kConvMap2005: return fFactorDip; + } +} + +//_____________________________________________________________________________ +AliMagF* AliMagF::CreateFieldMap(Float_t l3Cur, Float_t diCur, Int_t convention, Bool_t uniform, + Float_t beamenergy, const Char_t *beamtype, const Char_t *path) +{ + //------------------------------------------------ + // The magnetic field map, defined externally... + // L3 current 30000 A -> 0.5 T + // L3 current 12000 A -> 0.2 T + // dipole current 6000 A + // The polarities must match the convention (LHC or DCS2008) + // unless the special uniform map was used for MC + //------------------------------------------------ + const Float_t l3NominalCurrent1=30000.; // (A) + const Float_t l3NominalCurrent2=12000.; // (A) + const Float_t diNominalCurrent =6000. ; // (A) + + const Float_t tolerance=0.03; // relative current tolerance + const Float_t zero=77.; // "zero" current (A) + // + BMap_t map; + double sclL3,sclDip; + // + Float_t l3Pol = l3Cur > 0 ? 1:-1; + Float_t diPol = diCur > 0 ? 1:-1; + + l3Cur = TMath::Abs(l3Cur); + diCur = TMath::Abs(diCur); + // + if (TMath::Abs((sclDip=diCur/diNominalCurrent)-1.) > tolerance && !uniform) { + if (diCur <= zero) sclDip = 0.; // some small current.. -> Dipole OFF + else { + AliErrorGeneral("AliMagF",Form("Wrong dipole current (%f A)!",diCur)); + return 0; } - else if (integ != 0) - { - AliWarning("Precision integration flag set to 0 \n"); - fPrecInteg = 0; + } + // + if (uniform) { + // special treatment of special MC with uniform mag field (normalized to 0.5 T) + // no check for scaling/polarities are done + map = k5kGUniform; + sclL3 = l3Cur/l3NominalCurrent1; + } + else { + if (TMath::Abs((sclL3=l3Cur/l3NominalCurrent1)-1.) < tolerance) map = k5kG; + else if (TMath::Abs((sclL3=l3Cur/l3NominalCurrent2)-1.) < tolerance) map = k2kG; + else if (l3Cur <= zero && diCur<=zero) { sclL3=0; sclDip=0; map = k5kGUniform;} + else { + AliErrorGeneral("AliMagF",Form("Wrong L3 current (%f A)!",l3Cur)); + return 0; + } + } + // + if (sclDip!=0 && map!=k5kGUniform) { + if ( (l3Cur<=zero) || ((convention==kConvLHC && l3Pol!=diPol) || (convention==kConvDCS2008 && l3Pol==diPol)) ) { + AliErrorGeneral("AliMagF",Form("Wrong combination for L3/Dipole polarities (%c/%c) for convention %d", + l3Pol>0?'+':'-',diPol>0?'+':'-',GetPolarityConvention())); + return 0; } + } + // + if (l3Pol<0) sclL3 = -sclL3; + if (diPol<0) sclDip = -sclDip; + // + BeamType_t btype = kNoBeamField; + TString btypestr = beamtype; + btypestr.ToLower(); + TPRegexp protonBeam("(proton|p)\\s*-?\\s*\\1"); + TPRegexp ionBeam("(lead|pb|ion|a)\\s*-?\\s*\\1"); + if (btypestr.Contains(ionBeam)) btype = kBeamTypeAA; + else if (btypestr.Contains(protonBeam)) btype = kBeamTypepp; + else AliInfoGeneral("AliMagF",Form("Assume no LHC magnet field for the beam type %s, ",beamtype)); + char ttl[80]; + sprintf(ttl,"L3: %+5d Dip: %+4d kA; %s | Polarities in %s convention",(int)TMath::Sign(l3Cur,float(sclL3)), + (int)TMath::Sign(diCur,float(sclDip)),uniform ? " Constant":"", + convention==kConvLHC ? "LHC":"DCS2008"); + // LHC and DCS08 conventions have opposite dipole polarities + if ( GetPolarityConvention() != convention) sclDip = -sclDip; + // + return new AliMagF("MagneticFieldMap", ttl,sclL3,sclDip,map,btype,beamenergy,2,10.,path); + // } + +//_____________________________________________________________________________ +const char* AliMagF::GetBeamTypeText() const +{ + const char *beamNA = "No Beam"; + const char *beamPP = "p-p"; + const char *beamPbPb= "Pb-Pb"; + switch ( fBeamType ) { + case kBeamTypepp : return beamPP; + case kBeamTypeAA : return beamPbPb; + case kNoBeamField: + default: return beamNA; + } +} +