Line data Source code
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 :
16 :
17 : #include <TClass.h>
18 : #include <TFile.h>
19 : #include <TSystem.h>
20 : #include <TPRegexp.h>
21 :
22 : #include "AliMagF.h"
23 : #include "AliMagWrapCheb.h"
24 : #include "AliLog.h"
25 :
26 176 : ClassImp(AliMagF)
27 :
28 : const Double_t AliMagF::fgkSol2DipZ = -700.;
29 : const UShort_t AliMagF::fgkPolarityConvention = AliMagF::kConvLHC;
30 : /*
31 : Explanation for polarity conventions: these are the mapping between the
32 : current signs and main field components in L3 (Bz) and Dipole (Bx) (in Alice frame)
33 : 1) kConvMap2005: used for the field mapping in 2005
34 : positive L3 current -> negative Bz
35 : positive Dip current -> positive Bx
36 : 2) kConvMapDCS2008: defined by the microswitches/cabling of power converters as of 2008 - 1st half 2009
37 : positive L3 current -> positive Bz
38 : positive Dip current -> positive Bx
39 : 3) kConvLHC : defined by LHC
40 : positive L3 current -> positive Bz
41 : positive Dip current -> negative Bx
42 :
43 : Note: only "negative Bz(L3) with postive Bx(Dipole)" and its inverse was mapped in 2005. Hence
44 : the GRP Manager will reject the runs with the current combinations (in the convention defined by the
45 : static Int_t AliMagF::GetPolarityConvention()) which do not lead to such field polarities.
46 :
47 : -----------------------------------------------
48 :
49 : Explanation on integrals in the TPC region
50 : GetTPCInt(xyz,b) and GetTPCRatInt(xyz,b) give integrals from point (x,y,z) to point (x,y,0)
51 : (irrespectively of the z sign) of the following:
52 : TPCInt: b contains int{bx}, int{by}, int{bz}
53 : TPCRatInt: b contains int{bx/bz}, int{by/bz}, int{(bx/bz)^2+(by/bz)^2}
54 :
55 : The same applies to integral in cylindrical coordinates:
56 : GetTPCIntCyl(rphiz,b)
57 : GetTPCIntRatCyl(rphiz,b)
58 : They accept the R,Phi,Z coordinate (-pi<phi<pi) and return the field
59 : integrals in cyl. coordinates.
60 :
61 : Thus, to compute the integral from arbitrary xy_z1 to xy_z2, one should take
62 : b = b1-b2 with b1 and b2 coming from GetTPCInt(xy_z1,b1) and GetTPCInt(xy_z2,b2)
63 :
64 : Note: the integrals are defined for the range -300<Z<300 and 0<R<300
65 : */
66 : //_______________________________________________________________________
67 : AliMagF::AliMagF():
68 3 : TVirtualMagField(),
69 3 : fMeasuredMap(0),
70 3 : fMapType(k5kG),
71 3 : fSolenoid(0),
72 3 : fBeamType(kNoBeamField),
73 3 : fBeamEnergy(0),
74 : //
75 3 : fInteg(0),
76 3 : fPrecInteg(0),
77 3 : fFactorSol(1.),
78 3 : fFactorDip(1.),
79 3 : fMax(15),
80 3 : fDipoleOFF(kFALSE),
81 : //
82 3 : fQuadGradient(0),
83 3 : fDipoleField(0),
84 3 : fCCorrField(0),
85 3 : fACorr1Field(0),
86 3 : fACorr2Field(0),
87 3 : fParNames("","")
88 15 : {
89 : // Default constructor
90 : //
91 6 : }
92 :
93 : //_______________________________________________________________________
94 : AliMagF::AliMagF(const char *name, const char* title, Double_t factorSol, Double_t factorDip,
95 : BMap_t maptype, BeamType_t bt, Double_t be,Int_t integ, Double_t fmax, const char* path):
96 5 : TVirtualMagField(name),
97 5 : fMeasuredMap(0),
98 5 : fMapType(maptype),
99 5 : fSolenoid(0),
100 5 : fBeamType(bt),
101 5 : fBeamEnergy(be),
102 : //
103 5 : fInteg(integ),
104 5 : fPrecInteg(1),
105 5 : fFactorSol(1.),
106 5 : fFactorDip(1.),
107 5 : fMax(fmax),
108 5 : fDipoleOFF(factorDip==0.),
109 : //
110 5 : fQuadGradient(0),
111 5 : fDipoleField(0),
112 5 : fCCorrField(0),
113 5 : fACorr1Field(0),
114 5 : fACorr2Field(0),
115 5 : fParNames("","")
116 25 : {
117 : // Initialize the field with Geant integration option "integ" and max field "fmax,
118 : // Impose scaling of parameterized L3 field by factorSol and of dipole by factorDip.
119 : // The "be" is the energy of the beam in GeV/nucleon
120 : //
121 5 : SetTitle(title);
122 5 : if(integ<0 || integ > 2) {
123 0 : AliWarning(Form("Invalid magnetic field flag: %5d; Helix tracking chosen instead",integ));
124 0 : fInteg = 2;
125 0 : }
126 5 : if (fInteg == 0) fPrecInteg = 0;
127 : //
128 5 : if (fBeamEnergy<=0 && fBeamType!=kNoBeamField) {
129 0 : if (fBeamType == kBeamTypepp) fBeamEnergy = 7000.; // max proton energy
130 0 : else if (fBeamType == kBeamTypeAA) fBeamEnergy = 2760; // max PbPb energy
131 0 : else if (fBeamType == kBeamTypepA || fBeamType == kBeamTypeAp) fBeamEnergy = 2760; // same rigitiy max PbPb energy
132 0 : AliInfo("Maximim possible beam energy for requested beam is assumed");
133 : }
134 : const char* parname = 0;
135 : //
136 5 : if (fMapType == k2kG) parname = fDipoleOFF ? "Sol12_Dip0_Hole":"Sol12_Dip6_Hole";
137 10 : else if (fMapType == k5kG) parname = fDipoleOFF ? "Sol30_Dip0_Hole":"Sol30_Dip6_Hole";
138 0 : else if (fMapType == k5kGUniform) parname = "Sol30_Dip6_Uniform";
139 0 : else AliFatal(Form("Unknown field identifier %d is requested\n",fMapType));
140 : //
141 5 : SetDataFileName(path);
142 5 : SetParamName(parname);
143 : //
144 5 : LoadParameterization();
145 5 : InitMachineField(fBeamType,fBeamEnergy);
146 5 : double xyz[3]={0.,0.,0.};
147 10 : fSolenoid = GetBz(xyz);
148 5 : SetFactorSol(factorSol);
149 5 : SetFactorDip(factorDip);
150 5 : Print("a");
151 10 : }
152 :
153 : //_______________________________________________________________________
154 : AliMagF::AliMagF(const AliMagF &src):
155 0 : TVirtualMagField(src),
156 0 : fMeasuredMap(0),
157 0 : fMapType(src.fMapType),
158 0 : fSolenoid(src.fSolenoid),
159 0 : fBeamType(src.fBeamType),
160 0 : fBeamEnergy(src.fBeamEnergy),
161 0 : fInteg(src.fInteg),
162 0 : fPrecInteg(src.fPrecInteg),
163 0 : fFactorSol(src.fFactorSol),
164 0 : fFactorDip(src.fFactorDip),
165 0 : fMax(src.fMax),
166 0 : fDipoleOFF(src.fDipoleOFF),
167 0 : fQuadGradient(src.fQuadGradient),
168 0 : fDipoleField(src.fDipoleField),
169 0 : fCCorrField(src.fCCorrField),
170 0 : fACorr1Field(src.fACorr1Field),
171 0 : fACorr2Field(src.fACorr2Field),
172 0 : fParNames(src.fParNames)
173 0 : {
174 0 : if (src.fMeasuredMap) fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
175 0 : }
176 :
177 : //_______________________________________________________________________
178 : AliMagF::~AliMagF()
179 30 : {
180 10 : delete fMeasuredMap;
181 15 : }
182 :
183 : //_______________________________________________________________________
184 : Bool_t AliMagF::LoadParameterization()
185 : {
186 10 : if (fMeasuredMap) {
187 0 : AliFatal(Form("Field data %s are already loaded from %s\n",GetParamName(),GetDataFileName()));
188 0 : }
189 : //
190 5 : char* fname = gSystem->ExpandPathName(GetDataFileName());
191 5 : TFile* file = TFile::Open(fname);
192 5 : if (!file) {
193 0 : AliFatal(Form("Failed to open magnetic field data file %s\n",fname));
194 0 : }
195 : //
196 15 : fMeasuredMap = dynamic_cast<AliMagWrapCheb*>(file->Get(GetParamName()));
197 5 : if (!fMeasuredMap) {
198 0 : AliFatal(Form("Did not find field %s in %s\n",GetParamName(),fname));
199 0 : }
200 5 : file->Close();
201 10 : delete file;
202 5 : return kTRUE;
203 : }
204 :
205 :
206 : //_______________________________________________________________________
207 : void AliMagF::Field(const Double_t *xyz, Double_t *b)
208 : {
209 : // Method to calculate the field at point xyz
210 : //
211 : // b[0]=b[1]=b[2]=0.0;
212 15867394 : if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
213 2876176 : fMeasuredMap->Field(xyz,b);
214 28090408 : if (xyz[2]>fgkSol2DipZ || fDipoleOFF) for (int i=3;i--;) b[i] *= fFactorSol;
215 671352 : else for (int i=3;i--;) b[i] *= fFactorDip;
216 : }
217 1090874 : else MachineField(xyz, b);
218 : //
219 3967050 : }
220 :
221 : //_______________________________________________________________________
222 : Double_t AliMagF::GetBz(const Double_t *xyz) const
223 : {
224 : // Method to calculate the field at point xyz
225 : //
226 1189696 : if (fMeasuredMap && xyz[2]>fMeasuredMap->GetMinZ() && xyz[2]<fMeasuredMap->GetMaxZ()) {
227 297424 : double bz = fMeasuredMap->GetBz(xyz);
228 892272 : return (xyz[2]>fgkSol2DipZ || fDipoleOFF) ? bz*fFactorSol : bz*fFactorDip;
229 : }
230 0 : else return 0.;
231 297424 : }
232 :
233 : //_______________________________________________________________________
234 : AliMagF& AliMagF::operator=(const AliMagF& src)
235 : {
236 0 : if (this != &src) {
237 0 : if (src.fMeasuredMap) {
238 0 : if (fMeasuredMap) delete fMeasuredMap;
239 0 : fMeasuredMap = new AliMagWrapCheb(*src.fMeasuredMap);
240 0 : }
241 0 : SetName(src.GetName());
242 0 : fSolenoid = src.fSolenoid;
243 0 : fBeamType = src.fBeamType;
244 0 : fBeamEnergy = src.fBeamEnergy;
245 0 : fInteg = src.fInteg;
246 0 : fPrecInteg = src.fPrecInteg;
247 0 : fFactorSol = src.fFactorSol;
248 0 : fFactorDip = src.fFactorDip;
249 0 : fMax = src.fMax;
250 0 : fDipoleOFF = src.fDipoleOFF;
251 0 : fParNames = src.fParNames;
252 0 : }
253 0 : return *this;
254 0 : }
255 :
256 : //_______________________________________________________________________
257 : void AliMagF::InitMachineField(BeamType_t btype, Double_t benergy)
258 : {
259 10 : if (btype==kNoBeamField) {
260 4 : fQuadGradient = fDipoleField = fCCorrField = fACorr1Field = fACorr2Field = 0.;
261 4 : return;
262 : }
263 : //
264 1 : double rigScale = benergy/7000.; // scale according to ratio of E/Enominal
265 : // for ions assume PbPb (with energy provided per nucleon) and account for A/Z
266 2 : if (btype==kBeamTypeAA/* || btype==kBeamTypepA || btype==kBeamTypeAp */) rigScale *= 208./82.;
267 : // Attention: in p-Pb the energy recorded in the GRP is the PROTON energy, no rigidity
268 : // rescaling is needed
269 : //
270 1 : fQuadGradient = 22.0002*rigScale;
271 1 : fDipoleField = 37.8781*rigScale;
272 : //
273 : // SIDE C
274 1 : fCCorrField = -9.6980;
275 : // SIDE A
276 1 : fACorr1Field = -13.2247;
277 1 : fACorr2Field = 11.7905;
278 : //
279 6 : }
280 :
281 : //_______________________________________________________________________
282 : void AliMagF::MachineField(const Double_t *x, Double_t *b) const
283 : {
284 : // ---- This is the ZDC part
285 : // Compansators for Alice Muon Arm Dipole
286 : const Double_t kBComp1CZ = 1075., kBComp1hDZ = 260./2., kBComp1SqR = 4.0*4.0;
287 : const Double_t kBComp2CZ = 2049., kBComp2hDZ = 153./2., kBComp2SqR = 4.5*4.5;
288 : //
289 : const Double_t kTripQ1CZ = 2615., kTripQ1hDZ = 637./2., kTripQ1SqR = 3.5*3.5;
290 : const Double_t kTripQ2CZ = 3480., kTripQ2hDZ = 550./2., kTripQ2SqR = 3.5*3.5;
291 : const Double_t kTripQ3CZ = 4130., kTripQ3hDZ = 550./2., kTripQ3SqR = 3.5*3.5;
292 : const Double_t kTripQ4CZ = 5015., kTripQ4hDZ = 637./2., kTripQ4SqR = 3.5*3.5;
293 : //
294 : const Double_t kDip1CZ = 6310.8, kDip1hDZ = 945./2., kDip1SqRC = 4.5*4.5, kDip1SqRA = 3.375*3.375;
295 : const Double_t kDip2CZ = 12640.3, kDip2hDZ = 945./2., kDip2SqRC = 4.5*4.5, kDip2SqRA = 3.75*3.75;
296 : const Double_t kDip2DXC = 9.7, kDip2DXA = 9.4;
297 : //
298 2181748 : double rad2 = x[0] * x[0] + x[1] * x[1];
299 : //
300 1090874 : b[0] = b[1] = b[2] = 0;
301 : //
302 : // SIDE C **************************************************
303 1090874 : if(x[2]<0.){
304 806 : if(TMath::Abs(x[2]+kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
305 0 : b[0] = fCCorrField*fFactorDip;
306 0 : }
307 806 : else if(TMath::Abs(x[2]+kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
308 0 : b[0] = fQuadGradient*x[1];
309 0 : b[1] = fQuadGradient*x[0];
310 0 : }
311 806 : else if(TMath::Abs(x[2]+kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
312 0 : b[0] = -fQuadGradient*x[1];
313 0 : b[1] = -fQuadGradient*x[0];
314 0 : }
315 806 : else if(TMath::Abs(x[2]+kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
316 0 : b[0] = -fQuadGradient*x[1];
317 0 : b[1] = -fQuadGradient*x[0];
318 0 : }
319 806 : else if(TMath::Abs(x[2]+kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
320 0 : b[0] = fQuadGradient*x[1];
321 0 : b[1] = fQuadGradient*x[0];
322 0 : }
323 806 : else if(TMath::Abs(x[2]+kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRC){
324 0 : b[1] = fDipoleField;
325 0 : }
326 806 : else if(TMath::Abs(x[2]+kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRC) {
327 0 : double dxabs = TMath::Abs(x[0])-kDip2DXC;
328 0 : if ( (dxabs*dxabs + x[1]*x[1])<kDip2SqRC) {
329 0 : b[1] = -fDipoleField;
330 0 : }
331 0 : }
332 : }
333 : //
334 : // SIDE A **************************************************
335 : else{
336 1090068 : if(TMath::Abs(x[2]-kBComp1CZ)<kBComp1hDZ && rad2 < kBComp1SqR) {
337 : // Compensator magnet at z = 1075 m
338 7228 : b[0] = fACorr1Field*fFactorDip;
339 7228 : }
340 : //
341 1090068 : if(TMath::Abs(x[2]-kBComp2CZ)<kBComp2hDZ && rad2 < kBComp2SqR){
342 61 : b[0] = fACorr2Field*fFactorDip;
343 61 : }
344 1090007 : else if(TMath::Abs(x[2]-kTripQ1CZ)<kTripQ1hDZ && rad2 < kTripQ1SqR){
345 0 : b[0] = -fQuadGradient*x[1];
346 0 : b[1] = -fQuadGradient*x[0];
347 0 : }
348 1090007 : else if(TMath::Abs(x[2]-kTripQ2CZ)<kTripQ2hDZ && rad2 < kTripQ2SqR){
349 0 : b[0] = fQuadGradient*x[1];
350 0 : b[1] = fQuadGradient*x[0];
351 0 : }
352 1090007 : else if(TMath::Abs(x[2]-kTripQ3CZ)<kTripQ3hDZ && rad2 < kTripQ3SqR){
353 0 : b[0] = fQuadGradient*x[1];
354 0 : b[1] = fQuadGradient*x[0];
355 0 : }
356 1090007 : else if(TMath::Abs(x[2]-kTripQ4CZ)<kTripQ4hDZ && rad2 < kTripQ4SqR){
357 0 : b[0] = -fQuadGradient*x[1];
358 0 : b[1] = -fQuadGradient*x[0];
359 0 : }
360 1090007 : else if(TMath::Abs(x[2]-kDip1CZ)<kDip1hDZ && rad2 < kDip1SqRA){
361 0 : b[1] = -fDipoleField;
362 0 : }
363 1090007 : else if(TMath::Abs(x[2]-kDip2CZ)<kDip2hDZ && rad2 < kDip2SqRA) {
364 0 : double dxabs = TMath::Abs(x[0])-kDip2DXA;
365 0 : if ( (dxabs*dxabs + x[1]*x[1])<kDip2SqRA) {
366 0 : b[1] = fDipoleField;
367 0 : }
368 0 : }
369 : }
370 : //
371 1090874 : }
372 :
373 : //_______________________________________________________________________
374 : void AliMagF::GetTPCInt(const Double_t *xyz, Double_t *b) const
375 : {
376 : // Method to calculate the integral_0^z of br,bt,bz
377 0 : b[0]=b[1]=b[2]=0.0;
378 0 : if (fMeasuredMap) {
379 0 : fMeasuredMap->GetTPCInt(xyz,b);
380 0 : for (int i=3;i--;) b[i] *= fFactorSol;
381 0 : }
382 0 : }
383 :
384 : //_______________________________________________________________________
385 : void AliMagF::GetTPCRatInt(const Double_t *xyz, Double_t *b) const
386 : {
387 : // Method to calculate the integral_0^z of bx/bz,by/bz and (bx/bz)^2+(by/bz)^2
388 0 : b[0]=b[1]=b[2]=0.0;
389 0 : if (fMeasuredMap) {
390 0 : fMeasuredMap->GetTPCRatInt(xyz,b);
391 0 : b[2] /= 100;
392 0 : }
393 0 : }
394 :
395 : //_______________________________________________________________________
396 : void AliMagF::GetTPCIntCyl(const Double_t *rphiz, Double_t *b) const
397 : {
398 : // Method to calculate the integral_0^z of br,bt,bz
399 : // in cylindrical coordiates ( -pi<phi<pi convention )
400 0 : b[0]=b[1]=b[2]=0.0;
401 0 : if (fMeasuredMap) {
402 0 : fMeasuredMap->GetTPCIntCyl(rphiz,b);
403 0 : for (int i=3;i--;) b[i] *= fFactorSol;
404 0 : }
405 0 : }
406 :
407 : //_______________________________________________________________________
408 : void AliMagF::GetTPCRatIntCyl(const Double_t *rphiz, Double_t *b) const
409 : {
410 : // Method to calculate the integral_0^z of bx/bz,by/bz and (bx/bz)^2+(by/bz)^2
411 : // in cylindrical coordiates ( -pi<phi<pi convention )
412 0 : b[0]=b[1]=b[2]=0.0;
413 0 : if (fMeasuredMap) {
414 0 : fMeasuredMap->GetTPCRatIntCyl(rphiz,b);
415 0 : b[2] /= 100;
416 0 : }
417 0 : }
418 :
419 : //_______________________________________________________________________
420 : void AliMagF::SetFactorSol(Float_t fc)
421 : {
422 : // set the sign/scale of the current in the L3 according to fgkPolarityConvention
423 : switch (fgkPolarityConvention) {
424 : case kConvDCS2008: fFactorSol = -fc; break;
425 10 : case kConvLHC : fFactorSol = -fc; break;
426 : default : fFactorSol = fc; break; // case kConvMap2005: fFactorSol = fc; break;
427 : }
428 5 : }
429 :
430 : //_______________________________________________________________________
431 : void AliMagF::SetFactorDip(Float_t fc)
432 : {
433 : // set the sign*scale of the current in the Dipole according to fgkPolarityConvention
434 : switch (fgkPolarityConvention) {
435 : case kConvDCS2008: fFactorDip = fc; break;
436 10 : case kConvLHC : fFactorDip = -fc; break;
437 : default : fFactorDip = fc; break; // case kConvMap2005: fFactorDip = fc; break;
438 : }
439 5 : }
440 :
441 : //_______________________________________________________________________
442 : Double_t AliMagF::GetFactorSol() const
443 : {
444 : // return the sign*scale of the current in the Dipole according to fgkPolarityConventionthe
445 : switch (fgkPolarityConvention) {
446 : case kConvDCS2008: return -fFactorSol;
447 44 : case kConvLHC : return -fFactorSol;
448 : default : return fFactorSol; // case kConvMap2005: return fFactorSol;
449 : }
450 : }
451 :
452 : //_______________________________________________________________________
453 : Double_t AliMagF::GetFactorDip() const
454 : {
455 : // return the sign*scale of the current in the Dipole according to fgkPolarityConventionthe
456 : switch (fgkPolarityConvention) {
457 : case kConvDCS2008: return fFactorDip;
458 44 : case kConvLHC : return -fFactorDip;
459 : default : return fFactorDip; // case kConvMap2005: return fFactorDip;
460 : }
461 : }
462 :
463 : //_____________________________________________________________________________
464 : AliMagF* AliMagF::CreateFieldMap(Float_t l3Cur, Float_t diCur, Int_t convention, Bool_t uniform,
465 : Float_t beamenergy, const Char_t *beamtype, const Char_t *path,
466 : Bool_t returnNULLOnInvalidCurrent)
467 : {
468 : //------------------------------------------------
469 : // The magnetic field map, defined externally...
470 : // L3 current 30000 A -> 0.5 T
471 : // L3 current 12000 A -> 0.2 T
472 : // dipole current 6000 A
473 : // The polarities must match the convention (LHC or DCS2008)
474 : // unless the special uniform map was used for MC
475 : //------------------------------------------------
476 4 : AliLog::EType_t error_type = returnNULLOnInvalidCurrent ? AliLog::kError : AliLog::kFatal;
477 :
478 : const Float_t l3NominalCurrent1=30000.; // (A)
479 : const Float_t l3NominalCurrent2=12000.; // (A)
480 : const Float_t diNominalCurrent =6000. ; // (A)
481 :
482 : const Float_t tolerance=0.03; // relative current tolerance
483 : const Float_t zero=77.; // "zero" current (A)
484 : //
485 : BMap_t map = k5kG;
486 : double sclL3,sclDip;
487 : //
488 4 : Float_t l3Pol = l3Cur > 0 ? 1:-1;
489 4 : Float_t diPol = diCur > 0 ? 1:-1;
490 :
491 4 : l3Cur = TMath::Abs(l3Cur);
492 4 : diCur = TMath::Abs(diCur);
493 : //
494 4 : if (TMath::Abs((sclDip=diCur/diNominalCurrent)-1.) > tolerance && !uniform) {
495 0 : if (diCur <= zero) sclDip = 0.; // some small current.. -> Dipole OFF
496 : else {
497 0 : AliMessageGeneral("AliMagF",error_type,Form("Wrong dipole current (%f A)!",diCur));
498 0 : return NULL;
499 : }
500 0 : }
501 : //
502 4 : if (uniform) {
503 : // special treatment of special MC with uniform mag field (normalized to 0.5 T)
504 : // no check for scaling/polarities are done
505 : map = k5kGUniform;
506 0 : sclL3 = l3Cur/l3NominalCurrent1;
507 0 : }
508 : else {
509 8 : if (TMath::Abs((sclL3=l3Cur/l3NominalCurrent1)-1.) < tolerance) map = k5kG;
510 0 : else if (TMath::Abs((sclL3=l3Cur/l3NominalCurrent2)-1.) < tolerance) map = k2kG;
511 0 : else if (l3Cur <= zero && diCur<=zero) { sclL3=0; sclDip=0; map = k5kGUniform;}
512 : else {
513 0 : AliMessageGeneral("AliMagF",error_type,Form("Wrong L3 current (%f A)!",l3Cur));
514 0 : return NULL;
515 : }
516 : }
517 : //
518 4 : if (sclDip!=0 && map!=k5kGUniform) {
519 16 : if ( (l3Cur<=zero) || ((convention==kConvLHC && l3Pol!=diPol) || (convention==kConvDCS2008 && l3Pol==diPol)) ) {
520 0 : AliMessageGeneral("AliMagF",error_type,Form("Wrong combination for L3/Dipole polarities (%c/%c) for convention %d",
521 : l3Pol>0?'+':'-',diPol>0?'+':'-',GetPolarityConvention()));
522 0 : return NULL;
523 : }
524 : }
525 : //
526 8 : if (l3Pol<0) sclL3 = -sclL3;
527 8 : if (diPol<0) sclDip = -sclDip;
528 : //
529 : BeamType_t btype = kNoBeamField;
530 4 : TString btypestr = beamtype;
531 4 : btypestr.ToLower();
532 12 : TPRegexp protonBeam("(proton|p)\\s*-?\\s*\\1");
533 12 : TPRegexp ionBeam("(lead|pb|ion|a|A)\\s*-?\\s*\\1");
534 12 : TPRegexp protonionBeam("(proton|p)\\s*-?\\s*(lead|pb|ion|a|A)");
535 12 : TPRegexp ionprotonBeam("(lead|pb|ion|a|A)\\s*-?\\s*(proton|p)");
536 8 : if (btypestr.Contains(ionBeam)) btype = kBeamTypeAA;
537 8 : else if (btypestr.Contains(protonBeam)) btype = kBeamTypepp;
538 8 : else if (btypestr.Contains(protonionBeam)) btype = kBeamTypepA;
539 8 : else if (btypestr.Contains(ionprotonBeam)) btype = kBeamTypeAp;
540 8 : else AliInfoGeneral("AliMagF",Form("Assume no LHC magnet field for the beam type %s, ",beamtype));
541 4 : char ttl[80];
542 8 : snprintf(ttl,79,"L3: %+5d Dip: %+4d kA; %s | Polarities in %s convention",(int)TMath::Sign(l3Cur,float(sclL3)),
543 4 : (int)TMath::Sign(diCur,float(sclDip)),uniform ? " Constant":"",
544 4 : convention==kConvLHC ? "LHC":"DCS2008");
545 : // LHC and DCS08 conventions have opposite dipole polarities
546 4 : if ( GetPolarityConvention() != convention) sclDip = -sclDip;
547 : //
548 8 : return new AliMagF("MagneticFieldMap", ttl,sclL3,sclDip,map,btype,beamenergy,2,10.,path);
549 : //
550 8 : }
551 :
552 : //_____________________________________________________________________________
553 : const char* AliMagF::GetBeamTypeText() const
554 : {
555 : // beam type in text form
556 : const char *beamNA = "No Beam";
557 : const char *beamPP = "p-p";
558 : const char *beamPbPb= "A-A";
559 : const char *beamPPb = "p-A";
560 : const char *beamPbP = "A-p";
561 10 : switch ( fBeamType ) {
562 0 : case kBeamTypepp : return beamPP;
563 1 : case kBeamTypeAA : return beamPbPb;
564 0 : case kBeamTypepA : return beamPPb;
565 0 : case kBeamTypeAp : return beamPbP;
566 : case kNoBeamField:
567 4 : default: return beamNA;
568 : }
569 5 : }
570 :
571 : //_____________________________________________________________________________
572 : void AliMagF::Print(Option_t *opt) const
573 : {
574 : // print short or long info
575 10 : TString opts = opt; opts.ToLower();
576 25 : AliInfo(Form("%s:%s",GetName(),GetTitle()));
577 20 : AliInfo(Form("Solenoid (%+.2f*)%.0f kG, Dipole %s (%+.2f) %s",
578 : GetFactorSol(),(fMapType==k5kG||fMapType==k5kGUniform)?5.:2.,
579 : fDipoleOFF ? "OFF":"ON",GetFactorDip(),fMapType==k5kGUniform?" |Constant Field!":""));
580 10 : if (opts.Contains("a")) {
581 15 : AliInfo(Form("Machine B fields for %s beam (%.0f GeV): QGrad: %.4f Dipole: %.4f",
582 : GetBeamTypeText(),
583 : fBeamEnergy,fQuadGradient,fDipoleField));
584 25 : AliInfo(Form("Uses %s of %s",GetParamName(),GetDataFileName()));
585 : }
586 5 : }
|
