]>
Commit | Line | Data |
---|---|---|
4c039060 | 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 | ||
4c039060 | 16 | |
db83d72f | 17 | #include <TClass.h> |
18 | #include <TFile.h> | |
19 | #include <TSystem.h> | |
fe4da5cc | 20 | |
21 | #include "AliMagF.h" | |
db83d72f | 22 | #include "AliMagWrapCheb.h" |
23 | #include "AliLog.h" | |
972ca52f | 24 | |
fe4da5cc | 25 | ClassImp(AliMagF) |
26 | ||
9251fceb | 27 | const Double_t AliMagF::fgkSol2DipZ = -700.; |
db83d72f | 28 | |
e2afb3b6 | 29 | //_______________________________________________________________________ |
30 | AliMagF::AliMagF(): | |
db83d72f | 31 | TVirtualMagField(), |
32 | fMeasuredMap(0), | |
33 | fMapType(k5kG), | |
34 | fSolenoid(0), | |
35 | fBeamType(kNoBeamField), | |
36 | fBeamEnergy(0), | |
db83d72f | 37 | // |
e2afb3b6 | 38 | fInteg(0), |
db83d72f | 39 | fPrecInteg(0), |
40 | fFactorSol(1.), | |
41 | fFactorDip(1.), | |
42 | fMax(15), | |
43 | fDipoleOFF(kFALSE), | |
e2afb3b6 | 44 | // |
db83d72f | 45 | fQuadGradient(0), |
46 | fDipoleField(0), | |
47 | fCCorrField(0), | |
48 | fACorr1Field(0), | |
49 | fACorr2Field(0), | |
50 | fParNames("","") | |
51 | { | |
e2afb3b6 | 52 | // Default constructor |
53 | // | |
54 | } | |
55 | ||
56 | //_______________________________________________________________________ | |
db83d72f | 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, | |
9251fceb | 60 | BeamType_t bt, Double_t be): |
db83d72f | 61 | TVirtualMagField(name), |
62 | fMeasuredMap(0), | |
63 | fMapType(maptype), | |
64 | fSolenoid(0), | |
65 | fBeamType(bt), | |
66 | fBeamEnergy(be), | |
db83d72f | 67 | // |
68 | fInteg(integ), | |
604e0531 | 69 | fPrecInteg(1), |
db83d72f | 70 | fFactorSol(1.), |
71 | fFactorDip(1.), | |
972ca52f | 72 | fMax(fmax), |
db83d72f | 73 | fDipoleOFF(factorDip==0.), |
74 | // | |
75 | fQuadGradient(0), | |
76 | fDipoleField(0), | |
77 | fCCorrField(0), | |
78 | fACorr1Field(0), | |
79 | fACorr2Field(0), | |
80 | fParNames("","") | |
fe4da5cc | 81 | { |
9251fceb | 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 | |
aee8290b | 85 | // |
db83d72f | 86 | SetTitle(title); |
87 | if(integ<0 || integ > 2) { | |
88 | AliWarning(Form("Invalid magnetic field flag: %5d; Helix tracking chosen instead",integ)); | |
89 | fInteg = 2; | |
90 | } | |
91 | if (fInteg == 0) fPrecInteg = 0; | |
aee8290b | 92 | // |
db83d72f | 93 | const char* parname = 0; |
94 | // | |
95 | if (fMapType == k2kG) { | |
96 | fSolenoid = 2.; | |
97 | parname = fDipoleOFF ? "Sol12_Dip0_Hole":"Sol12_Dip6_Hole"; | |
98 | } else if (fMapType == k5kG) { | |
99 | fSolenoid = 5.; | |
100 | parname = fDipoleOFF ? "Sol30_Dip0_Hole":"Sol30_Dip6_Hole"; | |
101 | } else if (fMapType == k5kGUniform) { | |
102 | fSolenoid = 5.; | |
103 | parname = "Sol30_Dip6_Uniform"; | |
104 | } else { | |
105 | AliFatal(Form("Unknown field identifier %d is requested\n",fMapType)); | |
106 | } | |
107 | // | |
108 | SetDataFileName(path); | |
109 | SetParamName(parname); | |
110 | // | |
111 | SetFactorSol(factorSol); | |
112 | SetFactorDip(factorDip); | |
113 | LoadParameterization(); | |
114 | InitMachineField(fBeamType,fBeamEnergy); | |
fe4da5cc | 115 | } |
116 | ||
eeda4611 | 117 | //_______________________________________________________________________ |
118 | AliMagF::AliMagF(const AliMagF &src): | |
db83d72f | 119 | TVirtualMagField(src), |
120 | fMeasuredMap(0), | |
121 | fMapType(src.fMapType), | |
122 | fSolenoid(src.fSolenoid), | |
123 | fBeamType(src.fBeamType), | |
124 | fBeamEnergy(src.fBeamEnergy), | |
eeda4611 | 125 | fInteg(src.fInteg), |
126 | fPrecInteg(src.fPrecInteg), | |
db83d72f | 127 | fFactorSol(src.fFactorSol), |
128 | fFactorDip(src.fFactorDip), | |
eeda4611 | 129 | fMax(src.fMax), |
db83d72f | 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) | |
eeda4611 | 137 | { |
db83d72f | 138 | if (src.fMeasuredMap) fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap); |
eeda4611 | 139 | } |
140 | ||
e2afb3b6 | 141 | //_______________________________________________________________________ |
db83d72f | 142 | AliMagF::~AliMagF() |
ff66b122 | 143 | { |
db83d72f | 144 | delete fMeasuredMap; |
145 | } | |
146 | ||
147 | //_______________________________________________________________________ | |
148 | Bool_t AliMagF::LoadParameterization() | |
149 | { | |
150 | if (fMeasuredMap) { | |
151 | AliError(Form("Field data %s are already loaded from %s\n",GetParamName(),GetDataFileName())); | |
152 | return kTRUE; | |
153 | } | |
ff66b122 | 154 | // |
db83d72f | 155 | char* fname = gSystem->ExpandPathName(GetDataFileName()); |
156 | TFile* file = TFile::Open(fname); | |
157 | if (!file) { | |
158 | AliError(Form("Failed to open magnetic field data file %s\n",fname)); | |
159 | return kFALSE; | |
160 | } | |
ff66b122 | 161 | // |
db83d72f | 162 | fMeasuredMap = dynamic_cast<AliMagWrapCheb*>(file->Get(GetParamName())); |
163 | if (!fMeasuredMap) { | |
164 | AliError(Form("Did not find field %s in %s\n",GetParamName(),fname)); | |
165 | return kFALSE; | |
166 | } | |
167 | file->Close(); | |
168 | delete file; | |
169 | return kTRUE; | |
ff66b122 | 170 | } |
171 | ||
db83d72f | 172 | |
ff66b122 | 173 | //_______________________________________________________________________ |
db83d72f | 174 | void AliMagF::Field(const Double_t *xyz, Double_t *b) |
fe4da5cc | 175 | { |
db83d72f | 176 | // Method to calculate the field at point xyz |
aee8290b | 177 | // |
9251fceb | 178 | // b[0]=b[1]=b[2]=0.0; |
179 | if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) { | |
db83d72f | 180 | fMeasuredMap->Field(xyz,b); |
181 | if (xyz[2]>fgkSol2DipZ || fDipoleOFF) for (int i=3;i--;) b[i] *= fFactorSol; | |
9251fceb | 182 | else for (int i=3;i--;) b[i] *= fFactorDip; |
db83d72f | 183 | } |
9251fceb | 184 | else MachineField(xyz, b); |
aee8290b | 185 | // |
fe4da5cc | 186 | } |
eeda4611 | 187 | |
188 | //_______________________________________________________________________ | |
db83d72f | 189 | Double_t AliMagF::GetBz(const Double_t *xyz) const |
eeda4611 | 190 | { |
db83d72f | 191 | // Method to calculate the field at point xyz |
192 | // | |
9251fceb | 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; | |
db83d72f | 196 | } |
9251fceb | 197 | else return 0.; |
eeda4611 | 198 | } |
199 | ||
200 | //_______________________________________________________________________ | |
db83d72f | 201 | AliMagF& AliMagF::operator=(const AliMagF& src) |
eeda4611 | 202 | { |
db83d72f | 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; | |
db83d72f | 210 | fInteg = src.fInteg; |
211 | fPrecInteg = src.fPrecInteg; | |
212 | fFactorSol = src.fFactorSol; | |
213 | fFactorDip = src.fFactorDip; | |
214 | fMax = src.fMax; | |
215 | fDipoleOFF = src.fDipoleOFF; | |
216 | fParNames = src.fParNames; | |
217 | } | |
218 | return *this; | |
eeda4611 | 219 | } |
220 | ||
221 | //_______________________________________________________________________ | |
db83d72f | 222 | void AliMagF::InitMachineField(BeamType_t btype, Double_t benergy) |
eeda4611 | 223 | { |
9251fceb | 224 | if (btype==kNoBeamField || benergy<1.) { |
db83d72f | 225 | fQuadGradient = fDipoleField = fCCorrField = fACorr1Field = fACorr2Field = 0.; |
9251fceb | 226 | return; |
db83d72f | 227 | } |
228 | // | |
9251fceb | 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.; | |
232 | // | |
233 | fQuadGradient = 22.0002*rigScale; | |
234 | fDipoleField = 37.8781*rigScale; | |
235 | // | |
236 | // SIDE C | |
237 | fCCorrField = -9.6980; | |
238 | // SIDE A | |
239 | fACorr1Field = -13.2247; | |
240 | fACorr2Field = 11.7905; | |
db83d72f | 241 | // |
eeda4611 | 242 | } |
eed8a1a2 | 243 | |
db83d72f | 244 | //_______________________________________________________________________ |
245 | void AliMagF::MachineField(const Double_t *x, Double_t *b) const | |
eed8a1a2 | 246 | { |
db83d72f | 247 | // ---- This is the ZDC part |
9251fceb | 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; | |
251 | // | |
252 | const Double_t kTripQ1CZ = 2615., kTripQ1hDZ = 637./2., kTripQ1SqR = 3.5*3.5; | |
90ae20c9 | 253 | const Double_t kTripQ2CZ = 3480., kTripQ2hDZ = 550./2., kTripQ2SqR = 3.5*3.5; |
9251fceb | 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; | |
db83d72f | 256 | // |
9251fceb | 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; | |
db83d72f | 260 | // |
261 | double rad2 = x[0] * x[0] + x[1] * x[1]; | |
262 | // | |
9251fceb | 263 | b[0] = b[1] = b[2] = 0; |
264 | // | |
db83d72f | 265 | // SIDE C ************************************************** |
266 | if(x[2]<0.){ | |
9251fceb | 267 | if(TMath::Abs(x[2]+kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){ |
268 | b[0] = fCCorrField*fFactorDip; | |
db83d72f | 269 | } |
9251fceb | 270 | else if(TMath::Abs(x[2]+kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){ |
db83d72f | 271 | b[0] = fQuadGradient*x[1]; |
272 | b[1] = fQuadGradient*x[0]; | |
db83d72f | 273 | } |
9251fceb | 274 | else if(TMath::Abs(x[2]+kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){ |
db83d72f | 275 | b[0] = -fQuadGradient*x[1]; |
276 | b[1] = -fQuadGradient*x[0]; | |
db83d72f | 277 | } |
9251fceb | 278 | else if(TMath::Abs(x[2]+kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){ |
db83d72f | 279 | b[0] = -fQuadGradient*x[1]; |
280 | b[1] = -fQuadGradient*x[0]; | |
db83d72f | 281 | } |
9251fceb | 282 | else if(TMath::Abs(x[2]+kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){ |
db83d72f | 283 | b[0] = fQuadGradient*x[1]; |
284 | b[1] = fQuadGradient*x[0]; | |
db83d72f | 285 | } |
9251fceb | 286 | else if(TMath::Abs(x[2]+kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRC){ |
db83d72f | 287 | b[1] = fDipoleField; |
db83d72f | 288 | } |
9251fceb | 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) { | |
db83d72f | 292 | b[1] = -fDipoleField; |
db83d72f | 293 | } |
294 | } | |
295 | } | |
296 | // | |
297 | // SIDE A ************************************************** | |
298 | else{ | |
9251fceb | 299 | if(TMath::Abs(x[2]-kBComp1CZ)<kBComp1hDZ && rad2 < kBComp1SqR) { |
db83d72f | 300 | // Compensator magnet at z = 1075 m |
9251fceb | 301 | b[0] = fACorr1Field*fFactorDip; |
db83d72f | 302 | } |
303 | // | |
9251fceb | 304 | if(TMath::Abs(x[2]-kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){ |
305 | b[0] = fACorr2Field*fFactorDip; | |
306 | } | |
307 | else if(TMath::Abs(x[2]-kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){ | |
db83d72f | 308 | b[0] = -fQuadGradient*x[1]; |
309 | b[1] = -fQuadGradient*x[0]; | |
eed8a1a2 | 310 | } |
9251fceb | 311 | else if(TMath::Abs(x[2]-kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){ |
312 | b[0] = fQuadGradient*x[1]; | |
313 | b[1] = fQuadGradient*x[0]; | |
eed8a1a2 | 314 | } |
9251fceb | 315 | else if(TMath::Abs(x[2]-kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){ |
316 | b[0] = fQuadGradient*x[1]; | |
317 | b[1] = fQuadGradient*x[0]; | |
db83d72f | 318 | } |
9251fceb | 319 | else if(TMath::Abs(x[2]-kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){ |
db83d72f | 320 | b[0] = -fQuadGradient*x[1]; |
321 | b[1] = -fQuadGradient*x[0]; | |
db83d72f | 322 | } |
9251fceb | 323 | else if(TMath::Abs(x[2]-kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRA){ |
db83d72f | 324 | b[1] = -fDipoleField; |
db83d72f | 325 | } |
9251fceb | 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) { | |
db83d72f | 329 | b[1] = fDipoleField; |
330 | } | |
331 | } | |
332 | } | |
9251fceb | 333 | // |
db83d72f | 334 | } |
335 | ||
336 | //_______________________________________________________________________ | |
337 | void AliMagF::GetTPCInt(const Double_t *xyz, Double_t *b) const | |
338 | { | |
339 | // Method to calculate the integral of magnetic integral from xyz to nearest cathode plane | |
340 | b[0]=b[1]=b[2]=0.0; | |
341 | if (fMeasuredMap) { | |
342 | fMeasuredMap->GetTPCInt(xyz,b); | |
343 | for (int i=3;i--;) b[i] *= fFactorSol; | |
344 | } | |
345 | } | |
346 | ||
347 | //_______________________________________________________________________ | |
348 | void AliMagF::GetTPCIntCyl(const Double_t *rphiz, Double_t *b) const | |
349 | { | |
350 | // Method to calculate the integral of magnetic integral from point to nearest cathode plane | |
351 | // in cylindrical coordiates ( -pi<phi<pi convention ) | |
352 | b[0]=b[1]=b[2]=0.0; | |
353 | if (fMeasuredMap) { | |
354 | fMeasuredMap->GetTPCIntCyl(rphiz,b); | |
355 | for (int i=3;i--;) b[i] *= fFactorSol; | |
356 | } | |
eed8a1a2 | 357 | } |