1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
22 #include "AliMagWrapCheb.h"
27 const Double_t AliMagF::fgkSol2DipZ = -700.;
29 //_______________________________________________________________________
35 fBeamType(kNoBeamField),
52 // Default constructor
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),
73 fDipoleOFF(factorDip==0.),
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
87 if(integ<0 || integ > 2) {
88 AliWarning(Form("Invalid magnetic field flag: %5d; Helix tracking chosen instead",integ));
91 if (fInteg == 0) fPrecInteg = 0;
93 const char* parname = 0;
95 if (fMapType == k2kG) {
97 parname = fDipoleOFF ? "Sol12_Dip0_Hole":"Sol12_Dip6_Hole";
98 } else if (fMapType == k5kG) {
100 parname = fDipoleOFF ? "Sol30_Dip0_Hole":"Sol30_Dip6_Hole";
101 } else if (fMapType == k5kGUniform) {
103 parname = "Sol30_Dip6_Uniform";
105 AliFatal(Form("Unknown field identifier %d is requested\n",fMapType));
108 SetDataFileName(path);
109 SetParamName(parname);
111 SetFactorSol(factorSol);
112 SetFactorDip(factorDip);
113 LoadParameterization();
114 InitMachineField(fBeamType,fBeamEnergy);
117 //_______________________________________________________________________
118 AliMagF::AliMagF(const AliMagF &src):
119 TVirtualMagField(src),
121 fMapType(src.fMapType),
122 fSolenoid(src.fSolenoid),
123 fBeamType(src.fBeamType),
124 fBeamEnergy(src.fBeamEnergy),
126 fPrecInteg(src.fPrecInteg),
127 fFactorSol(src.fFactorSol),
128 fFactorDip(src.fFactorDip),
130 fDipoleOFF(src.fDipoleOFF),
131 fQuadGradient(src.fQuadGradient),
132 fDipoleField(src.fDipoleField),
133 fCCorrField(src.fCCorrField),
134 fACorr1Field(src.fACorr1Field),
135 fACorr2Field(src.fACorr2Field),
136 fParNames(src.fParNames)
138 if (src.fMeasuredMap) fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
141 //_______________________________________________________________________
147 //_______________________________________________________________________
148 Bool_t AliMagF::LoadParameterization()
151 AliError(Form("Field data %s are already loaded from %s\n",GetParamName(),GetDataFileName()));
155 char* fname = gSystem->ExpandPathName(GetDataFileName());
156 TFile* file = TFile::Open(fname);
158 AliError(Form("Failed to open magnetic field data file %s\n",fname));
162 fMeasuredMap = dynamic_cast<AliMagWrapCheb*>(file->Get(GetParamName()));
164 AliError(Form("Did not find field %s in %s\n",GetParamName(),fname));
173 //_______________________________________________________________________
174 void AliMagF::Field(const Double_t *xyz, Double_t *b)
176 // Method to calculate the field at point xyz
178 // b[0]=b[1]=b[2]=0.0;
179 if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
180 fMeasuredMap->Field(xyz,b);
181 if (xyz[2]>fgkSol2DipZ || fDipoleOFF) for (int i=3;i--;) b[i] *= fFactorSol;
182 else for (int i=3;i--;) b[i] *= fFactorDip;
184 else MachineField(xyz, b);
188 //_______________________________________________________________________
189 Double_t AliMagF::GetBz(const Double_t *xyz) const
191 // Method to calculate the field at point xyz
193 if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
194 double bz = fMeasuredMap->GetBz(xyz);
195 return (xyz[2]>fgkSol2DipZ || fDipoleOFF) ? bz*fFactorSol : bz*fFactorDip;
200 //_______________________________________________________________________
201 AliMagF& AliMagF::operator=(const AliMagF& src)
203 if (this != &src && src.fMeasuredMap) {
204 if (fMeasuredMap) delete fMeasuredMap;
205 fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
206 SetName(src.GetName());
207 fSolenoid = src.fSolenoid;
208 fBeamType = src.fBeamType;
209 fBeamEnergy = src.fBeamEnergy;
211 fPrecInteg = src.fPrecInteg;
212 fFactorSol = src.fFactorSol;
213 fFactorDip = src.fFactorDip;
215 fDipoleOFF = src.fDipoleOFF;
216 fParNames = src.fParNames;
221 //_______________________________________________________________________
222 void AliMagF::InitMachineField(BeamType_t btype, Double_t benergy)
224 if (btype==kNoBeamField || benergy<1.) {
225 fQuadGradient = fDipoleField = fCCorrField = fACorr1Field = fACorr2Field = 0.;
229 double rigScale = benergy/7000.; // scale according to ratio of E/Enominal
230 // for ions assume PbPb (with energy provided per nucleon) and account for A/Z
231 if (btype == kBeamTypeAA) rigScale *= 208./82.;
233 fQuadGradient = 22.0002*rigScale;
234 fDipoleField = 37.8781*rigScale;
237 fCCorrField = -9.6980;
239 fACorr1Field = -13.2247;
240 fACorr2Field = 11.7905;
244 //_______________________________________________________________________
245 void AliMagF::MachineField(const Double_t *x, Double_t *b) const
247 // ---- This is the ZDC part
248 // Compansators for Alice Muon Arm Dipole
249 const Double_t kBComp1CZ = 1075., kBComp1hDZ = 260./2., kBComp1SqR = 4.0*4.0;
250 const Double_t kBComp2CZ = 2049., kBComp2hDZ = 153./2., kBComp2SqR = 4.5*4.5;
252 const Double_t kTripQ1CZ = 2615., kTripQ1hDZ = 637./2., kTripQ1SqR = 3.5*3.5;
253 const Double_t kTripQ2CZ = 3480., kTripQ2hDZ = 550./2., kTripQ2SqR = 3.5*3.5;
254 const Double_t kTripQ3CZ = 4130., kTripQ3hDZ = 550./2., kTripQ3SqR = 3.5*3.5;
255 const Double_t kTripQ4CZ = 5015., kTripQ4hDZ = 637./2., kTripQ4SqR = 3.5*3.5;
257 const Double_t kDip1CZ = 6310.8, kDip1hDZ = 945./2., kDip1SqRC = 4.5*4.5, kDip1SqRA = 3.375*3.375;
258 const Double_t kDip2CZ = 12640.3, kDip2hDZ = 945./2., kDip2SqRC = 4.5*4.5, kDip2SqRA = 3.75*3.75;
259 const Double_t kDip2DXC = 9.7, kDip2DXA = 9.4;
261 double rad2 = x[0] * x[0] + x[1] * x[1];
263 b[0] = b[1] = b[2] = 0;
265 // SIDE C **************************************************
267 if(TMath::Abs(x[2]+kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
268 b[0] = fCCorrField*fFactorDip;
270 else if(TMath::Abs(x[2]+kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
271 b[0] = fQuadGradient*x[1];
272 b[1] = fQuadGradient*x[0];
274 else if(TMath::Abs(x[2]+kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
275 b[0] = -fQuadGradient*x[1];
276 b[1] = -fQuadGradient*x[0];
278 else if(TMath::Abs(x[2]+kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
279 b[0] = -fQuadGradient*x[1];
280 b[1] = -fQuadGradient*x[0];
282 else if(TMath::Abs(x[2]+kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
283 b[0] = fQuadGradient*x[1];
284 b[1] = fQuadGradient*x[0];
286 else if(TMath::Abs(x[2]+kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRC){
289 else if(TMath::Abs(x[2]+kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRC) {
290 double dxabs = TMath::Abs(x[0])-kDip2DXC;
291 if ( (dxabs*dxabs + x[1]*x[1])<kDip2SqRC) {
292 b[1] = -fDipoleField;
297 // SIDE A **************************************************
299 if(TMath::Abs(x[2]-kBComp1CZ)<kBComp1hDZ && rad2 < kBComp1SqR) {
300 // Compensator magnet at z = 1075 m
301 b[0] = fACorr1Field*fFactorDip;
304 if(TMath::Abs(x[2]-kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
305 b[0] = fACorr2Field*fFactorDip;
307 else if(TMath::Abs(x[2]-kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
308 b[0] = -fQuadGradient*x[1];
309 b[1] = -fQuadGradient*x[0];
311 else if(TMath::Abs(x[2]-kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
312 b[0] = fQuadGradient*x[1];
313 b[1] = fQuadGradient*x[0];
315 else if(TMath::Abs(x[2]-kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
316 b[0] = fQuadGradient*x[1];
317 b[1] = fQuadGradient*x[0];
319 else if(TMath::Abs(x[2]-kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
320 b[0] = -fQuadGradient*x[1];
321 b[1] = -fQuadGradient*x[0];
323 else if(TMath::Abs(x[2]-kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRA){
324 b[1] = -fDipoleField;
326 else if(TMath::Abs(x[2]-kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRA) {
327 double dxabs = TMath::Abs(x[0])-kDip2DXA;
328 if ( (dxabs*dxabs + x[1]*x[1])<kDip2SqRA) {
336 //_______________________________________________________________________
337 void AliMagF::GetTPCInt(const Double_t *xyz, Double_t *b) const
339 // Method to calculate the integral of magnetic integral from xyz to nearest cathode plane
342 fMeasuredMap->GetTPCInt(xyz,b);
343 for (int i=3;i--;) b[i] *= fFactorSol;
347 //_______________________________________________________________________
348 void AliMagF::GetTPCIntCyl(const Double_t *rphiz, Double_t *b) const
350 // Method to calculate the integral of magnetic integral from point to nearest cathode plane
351 // in cylindrical coordiates ( -pi<phi<pi convention )
354 fMeasuredMap->GetTPCIntCyl(rphiz,b);
355 for (int i=3;i--;) b[i] *= fFactorSol;