]>
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 | // | |
f04e7f5f | 95 | if (fMapType == k2kG) parname = fDipoleOFF ? "Sol12_Dip0_Hole":"Sol12_Dip6_Hole"; |
96 | else if (fMapType == k5kG) parname = fDipoleOFF ? "Sol30_Dip0_Hole":"Sol30_Dip6_Hole"; | |
97 | else if (fMapType == k5kGUniform) parname = "Sol30_Dip6_Uniform"; | |
98 | else AliFatal(Form("Unknown field identifier %d is requested\n",fMapType)); | |
db83d72f | 99 | // |
100 | SetDataFileName(path); | |
101 | SetParamName(parname); | |
102 | // | |
db83d72f | 103 | LoadParameterization(); |
104 | InitMachineField(fBeamType,fBeamEnergy); | |
f04e7f5f | 105 | double xyz[3]={0.,0.,0.}; |
106 | fSolenoid = GetBz(xyz); | |
107 | SetFactorSol(factorSol); | |
108 | SetFactorDip(factorDip); | |
fe4da5cc | 109 | } |
110 | ||
eeda4611 | 111 | //_______________________________________________________________________ |
112 | AliMagF::AliMagF(const AliMagF &src): | |
db83d72f | 113 | TVirtualMagField(src), |
114 | fMeasuredMap(0), | |
115 | fMapType(src.fMapType), | |
116 | fSolenoid(src.fSolenoid), | |
117 | fBeamType(src.fBeamType), | |
118 | fBeamEnergy(src.fBeamEnergy), | |
eeda4611 | 119 | fInteg(src.fInteg), |
120 | fPrecInteg(src.fPrecInteg), | |
db83d72f | 121 | fFactorSol(src.fFactorSol), |
122 | fFactorDip(src.fFactorDip), | |
eeda4611 | 123 | fMax(src.fMax), |
db83d72f | 124 | fDipoleOFF(src.fDipoleOFF), |
125 | fQuadGradient(src.fQuadGradient), | |
126 | fDipoleField(src.fDipoleField), | |
127 | fCCorrField(src.fCCorrField), | |
128 | fACorr1Field(src.fACorr1Field), | |
129 | fACorr2Field(src.fACorr2Field), | |
130 | fParNames(src.fParNames) | |
eeda4611 | 131 | { |
db83d72f | 132 | if (src.fMeasuredMap) fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap); |
eeda4611 | 133 | } |
134 | ||
e2afb3b6 | 135 | //_______________________________________________________________________ |
db83d72f | 136 | AliMagF::~AliMagF() |
ff66b122 | 137 | { |
db83d72f | 138 | delete fMeasuredMap; |
139 | } | |
140 | ||
141 | //_______________________________________________________________________ | |
142 | Bool_t AliMagF::LoadParameterization() | |
143 | { | |
144 | if (fMeasuredMap) { | |
145 | AliError(Form("Field data %s are already loaded from %s\n",GetParamName(),GetDataFileName())); | |
146 | return kTRUE; | |
147 | } | |
ff66b122 | 148 | // |
db83d72f | 149 | char* fname = gSystem->ExpandPathName(GetDataFileName()); |
150 | TFile* file = TFile::Open(fname); | |
151 | if (!file) { | |
152 | AliError(Form("Failed to open magnetic field data file %s\n",fname)); | |
153 | return kFALSE; | |
154 | } | |
ff66b122 | 155 | // |
db83d72f | 156 | fMeasuredMap = dynamic_cast<AliMagWrapCheb*>(file->Get(GetParamName())); |
157 | if (!fMeasuredMap) { | |
158 | AliError(Form("Did not find field %s in %s\n",GetParamName(),fname)); | |
159 | return kFALSE; | |
160 | } | |
161 | file->Close(); | |
162 | delete file; | |
163 | return kTRUE; | |
ff66b122 | 164 | } |
165 | ||
db83d72f | 166 | |
ff66b122 | 167 | //_______________________________________________________________________ |
db83d72f | 168 | void AliMagF::Field(const Double_t *xyz, Double_t *b) |
fe4da5cc | 169 | { |
db83d72f | 170 | // Method to calculate the field at point xyz |
aee8290b | 171 | // |
9251fceb | 172 | // b[0]=b[1]=b[2]=0.0; |
173 | if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) { | |
db83d72f | 174 | fMeasuredMap->Field(xyz,b); |
175 | if (xyz[2]>fgkSol2DipZ || fDipoleOFF) for (int i=3;i--;) b[i] *= fFactorSol; | |
9251fceb | 176 | else for (int i=3;i--;) b[i] *= fFactorDip; |
db83d72f | 177 | } |
9251fceb | 178 | else MachineField(xyz, b); |
aee8290b | 179 | // |
fe4da5cc | 180 | } |
eeda4611 | 181 | |
182 | //_______________________________________________________________________ | |
db83d72f | 183 | Double_t AliMagF::GetBz(const Double_t *xyz) const |
eeda4611 | 184 | { |
db83d72f | 185 | // Method to calculate the field at point xyz |
186 | // | |
9251fceb | 187 | if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) { |
188 | double bz = fMeasuredMap->GetBz(xyz); | |
189 | return (xyz[2]>fgkSol2DipZ || fDipoleOFF) ? bz*fFactorSol : bz*fFactorDip; | |
db83d72f | 190 | } |
9251fceb | 191 | else return 0.; |
eeda4611 | 192 | } |
193 | ||
194 | //_______________________________________________________________________ | |
db83d72f | 195 | AliMagF& AliMagF::operator=(const AliMagF& src) |
eeda4611 | 196 | { |
db83d72f | 197 | if (this != &src && src.fMeasuredMap) { |
198 | if (fMeasuredMap) delete fMeasuredMap; | |
199 | fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap); | |
200 | SetName(src.GetName()); | |
201 | fSolenoid = src.fSolenoid; | |
202 | fBeamType = src.fBeamType; | |
203 | fBeamEnergy = src.fBeamEnergy; | |
db83d72f | 204 | fInteg = src.fInteg; |
205 | fPrecInteg = src.fPrecInteg; | |
206 | fFactorSol = src.fFactorSol; | |
207 | fFactorDip = src.fFactorDip; | |
208 | fMax = src.fMax; | |
209 | fDipoleOFF = src.fDipoleOFF; | |
210 | fParNames = src.fParNames; | |
211 | } | |
212 | return *this; | |
eeda4611 | 213 | } |
214 | ||
215 | //_______________________________________________________________________ | |
db83d72f | 216 | void AliMagF::InitMachineField(BeamType_t btype, Double_t benergy) |
eeda4611 | 217 | { |
9251fceb | 218 | if (btype==kNoBeamField || benergy<1.) { |
db83d72f | 219 | fQuadGradient = fDipoleField = fCCorrField = fACorr1Field = fACorr2Field = 0.; |
9251fceb | 220 | return; |
db83d72f | 221 | } |
222 | // | |
9251fceb | 223 | double rigScale = benergy/7000.; // scale according to ratio of E/Enominal |
224 | // for ions assume PbPb (with energy provided per nucleon) and account for A/Z | |
225 | if (btype == kBeamTypeAA) rigScale *= 208./82.; | |
226 | // | |
227 | fQuadGradient = 22.0002*rigScale; | |
228 | fDipoleField = 37.8781*rigScale; | |
229 | // | |
230 | // SIDE C | |
231 | fCCorrField = -9.6980; | |
232 | // SIDE A | |
233 | fACorr1Field = -13.2247; | |
234 | fACorr2Field = 11.7905; | |
db83d72f | 235 | // |
eeda4611 | 236 | } |
eed8a1a2 | 237 | |
db83d72f | 238 | //_______________________________________________________________________ |
239 | void AliMagF::MachineField(const Double_t *x, Double_t *b) const | |
eed8a1a2 | 240 | { |
db83d72f | 241 | // ---- This is the ZDC part |
9251fceb | 242 | // Compansators for Alice Muon Arm Dipole |
243 | const Double_t kBComp1CZ = 1075., kBComp1hDZ = 260./2., kBComp1SqR = 4.0*4.0; | |
244 | const Double_t kBComp2CZ = 2049., kBComp2hDZ = 153./2., kBComp2SqR = 4.5*4.5; | |
245 | // | |
246 | const Double_t kTripQ1CZ = 2615., kTripQ1hDZ = 637./2., kTripQ1SqR = 3.5*3.5; | |
90ae20c9 | 247 | const Double_t kTripQ2CZ = 3480., kTripQ2hDZ = 550./2., kTripQ2SqR = 3.5*3.5; |
9251fceb | 248 | const Double_t kTripQ3CZ = 4130., kTripQ3hDZ = 550./2., kTripQ3SqR = 3.5*3.5; |
249 | const Double_t kTripQ4CZ = 5015., kTripQ4hDZ = 637./2., kTripQ4SqR = 3.5*3.5; | |
db83d72f | 250 | // |
9251fceb | 251 | const Double_t kDip1CZ = 6310.8, kDip1hDZ = 945./2., kDip1SqRC = 4.5*4.5, kDip1SqRA = 3.375*3.375; |
252 | const Double_t kDip2CZ = 12640.3, kDip2hDZ = 945./2., kDip2SqRC = 4.5*4.5, kDip2SqRA = 3.75*3.75; | |
253 | const Double_t kDip2DXC = 9.7, kDip2DXA = 9.4; | |
db83d72f | 254 | // |
255 | double rad2 = x[0] * x[0] + x[1] * x[1]; | |
256 | // | |
9251fceb | 257 | b[0] = b[1] = b[2] = 0; |
258 | // | |
db83d72f | 259 | // SIDE C ************************************************** |
260 | if(x[2]<0.){ | |
9251fceb | 261 | if(TMath::Abs(x[2]+kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){ |
262 | b[0] = fCCorrField*fFactorDip; | |
db83d72f | 263 | } |
9251fceb | 264 | else if(TMath::Abs(x[2]+kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){ |
db83d72f | 265 | b[0] = fQuadGradient*x[1]; |
266 | b[1] = fQuadGradient*x[0]; | |
db83d72f | 267 | } |
9251fceb | 268 | else if(TMath::Abs(x[2]+kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){ |
db83d72f | 269 | b[0] = -fQuadGradient*x[1]; |
270 | b[1] = -fQuadGradient*x[0]; | |
db83d72f | 271 | } |
9251fceb | 272 | else if(TMath::Abs(x[2]+kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){ |
db83d72f | 273 | b[0] = -fQuadGradient*x[1]; |
274 | b[1] = -fQuadGradient*x[0]; | |
db83d72f | 275 | } |
9251fceb | 276 | else if(TMath::Abs(x[2]+kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){ |
db83d72f | 277 | b[0] = fQuadGradient*x[1]; |
278 | b[1] = fQuadGradient*x[0]; | |
db83d72f | 279 | } |
9251fceb | 280 | else if(TMath::Abs(x[2]+kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRC){ |
db83d72f | 281 | b[1] = fDipoleField; |
db83d72f | 282 | } |
9251fceb | 283 | else if(TMath::Abs(x[2]+kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRC) { |
284 | double dxabs = TMath::Abs(x[0])-kDip2DXC; | |
285 | if ( (dxabs*dxabs + x[1]*x[1])<kDip2SqRC) { | |
db83d72f | 286 | b[1] = -fDipoleField; |
db83d72f | 287 | } |
288 | } | |
289 | } | |
290 | // | |
291 | // SIDE A ************************************************** | |
292 | else{ | |
9251fceb | 293 | if(TMath::Abs(x[2]-kBComp1CZ)<kBComp1hDZ && rad2 < kBComp1SqR) { |
db83d72f | 294 | // Compensator magnet at z = 1075 m |
9251fceb | 295 | b[0] = fACorr1Field*fFactorDip; |
db83d72f | 296 | } |
297 | // | |
9251fceb | 298 | if(TMath::Abs(x[2]-kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){ |
299 | b[0] = fACorr2Field*fFactorDip; | |
300 | } | |
301 | else if(TMath::Abs(x[2]-kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){ | |
db83d72f | 302 | b[0] = -fQuadGradient*x[1]; |
303 | b[1] = -fQuadGradient*x[0]; | |
eed8a1a2 | 304 | } |
9251fceb | 305 | else if(TMath::Abs(x[2]-kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){ |
306 | b[0] = fQuadGradient*x[1]; | |
307 | b[1] = fQuadGradient*x[0]; | |
eed8a1a2 | 308 | } |
9251fceb | 309 | else if(TMath::Abs(x[2]-kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){ |
310 | b[0] = fQuadGradient*x[1]; | |
311 | b[1] = fQuadGradient*x[0]; | |
db83d72f | 312 | } |
9251fceb | 313 | else if(TMath::Abs(x[2]-kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){ |
db83d72f | 314 | b[0] = -fQuadGradient*x[1]; |
315 | b[1] = -fQuadGradient*x[0]; | |
db83d72f | 316 | } |
9251fceb | 317 | else if(TMath::Abs(x[2]-kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRA){ |
db83d72f | 318 | b[1] = -fDipoleField; |
db83d72f | 319 | } |
9251fceb | 320 | else if(TMath::Abs(x[2]-kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRA) { |
321 | double dxabs = TMath::Abs(x[0])-kDip2DXA; | |
322 | if ( (dxabs*dxabs + x[1]*x[1])<kDip2SqRA) { | |
db83d72f | 323 | b[1] = fDipoleField; |
324 | } | |
325 | } | |
326 | } | |
9251fceb | 327 | // |
db83d72f | 328 | } |
329 | ||
330 | //_______________________________________________________________________ | |
331 | void AliMagF::GetTPCInt(const Double_t *xyz, Double_t *b) const | |
332 | { | |
333 | // Method to calculate the integral of magnetic integral from xyz to nearest cathode plane | |
334 | b[0]=b[1]=b[2]=0.0; | |
335 | if (fMeasuredMap) { | |
336 | fMeasuredMap->GetTPCInt(xyz,b); | |
337 | for (int i=3;i--;) b[i] *= fFactorSol; | |
338 | } | |
339 | } | |
340 | ||
341 | //_______________________________________________________________________ | |
342 | void AliMagF::GetTPCIntCyl(const Double_t *rphiz, Double_t *b) const | |
343 | { | |
344 | // Method to calculate the integral of magnetic integral from point to nearest cathode plane | |
345 | // in cylindrical coordiates ( -pi<phi<pi convention ) | |
346 | b[0]=b[1]=b[2]=0.0; | |
347 | if (fMeasuredMap) { | |
348 | fMeasuredMap->GetTPCIntCyl(rphiz,b); | |
349 | for (int i=3;i--;) b[i] *= fFactorSol; | |
350 | } | |
eed8a1a2 | 351 | } |