/************************************************************************** * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * * Author: The ALICE Off-line Project. * * Contributors are mentioned in the code where appropriate. * * * * Permission to use, copy, modify and distribute this software and its * * documentation strictly for non-commercial purposes is hereby granted * * without fee, provided that the above copyright notice appears in all * * copies and that both the copyright notice and this permission notice * * appear in the supporting documentation. The authors make no claims * * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ #include #include #include #include #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 -> negative 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(): TVirtualMagField(), fMeasuredMap(0), fMapType(k5kG), fSolenoid(0), fBeamType(kNoBeamField), fBeamEnergy(0), // fInteg(0), 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, 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), fFactorSol(1.), fFactorDip(1.), fMax(fmax), 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 // 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 (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): TVirtualMagField(src), fMeasuredMap(0), fMapType(src.fMapType), 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), 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() { 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.; } //_______________________________________________________________________ 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; // // SIDE C fCCorrField = -9.6980; // SIDE A fACorr1Field = -13.2247; fACorr2Field = 11.7905; // } //_______________________________________________________________________ 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]; // 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::SetFactorSol(Float_t fc) { // 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::SetFactorDip(Float_t fc) { // 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; } } //_______________________________________________________________________ Double_t AliMagF::GetFactorSol() const { // 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; } } //_______________________________________________________________________ Double_t AliMagF::GetFactorDip() const { // 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 sqrts, 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; } } // 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) { sclL3 = 0; map = k5kGUniform;} else { AliErrorGeneral("AliMagF",Form("Wrong L3 current (%f A)!",l3Cur)); return 0; } } // if (sclDip!=0 && (map==k5kG || map==k2kG) && ((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,sqrts/2,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; } }