Moved old AliMagFCheb to AliMagF, small fixes/optimizations in field classes
[u/mrichter/AliRoot.git] / STEER / AliMagFC.cxx
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 /* $Id$ */
17
18 //-------------------------------------------------------------------------
19 //     Constant magnetic field class
20 //     Used by AliRun class
21 //     Author:
22 //-------------------------------------------------------------------------
23
24 #include <stdlib.h>
25
26 #include "AliLog.h"
27 #include "AliMagFC.h"
28
29 ClassImp(AliMagFC)
30
31 //________________________________________
32 AliMagFC::AliMagFC()
33     :AliMagF(),
34     fCompensator(kFALSE),
35     fBeamType(kBeamTypepp),
36     fBeamEnergy(0),
37     fQuadGradient(0),
38     fDipoleField(0),
39     fCCorrField(0), 
40     fACorr1Field(0),
41     fACorr2Field(0)
42 {
43   // 
44   // Default constructor
45   //
46 }
47
48 //________________________________________
49 AliMagFC::AliMagFC(const char *name, const char *title, Int_t integ, 
50                    Float_t factor, Float_t fmax)
51     : AliMagF(name,title,integ,factor,fmax),
52     fCompensator(kFALSE),
53     fBeamType(kBeamTypepp), 
54     fBeamEnergy(7000.),
55     fQuadGradient(0),
56     fDipoleField(0),
57     fCCorrField(0), 
58     fACorr1Field(0),
59     fACorr2Field(0)
60
61 {
62   // 
63   // Standard constructor
64   //
65   fType = kConst;
66   fMap  = 1;
67   
68 }
69
70 //________________________________________
71 void AliMagFC::Field(const float *x, float *b) const
72 {
73   //
74   // Method to return the field in a point
75   //
76   b[0]=b[1]=b[2]=0;
77   if(fMap==1) {
78     if(TMath::Abs(x[2])<700 && x[0]*x[0]+(x[1]+30)*(x[1]+30) < 560*560) {
79       b[2]=2;
80     } 
81     else {
82       if(-725 >= x[2] && x[2] >= -1225 ){
83         Float_t dz = TMath::Abs(-975-x[2])*0.01;
84         b[0] = - (1-0.1*dz*dz)*7;
85         if(fFactor!=1) {
86             b[0]*=fFactor;
87             b[1]*=fFactor;
88             b[2]*=fFactor;
89         }
90       }
91       else {
92           ZDCField(x, b);
93       }
94     }
95
96   } 
97   else {
98       AliFatal(Form("Invalid field map for constant field %d",fMap));
99   }
100 }
101 //________________________________________
102 void AliMagFC::Field(const double *x, double *b) const
103 {
104   //
105   // Method to return the field in a point
106   //
107   b[0]=b[1]=b[2]=0;
108   if(fMap==1) {
109     if(TMath::Abs(x[2])<700 && x[0]*x[0]+(x[1]+30)*(x[1]+30) < 560*560) {
110       b[2]=2;
111     } 
112     else {
113       if(-725 >= x[2] && x[2] >= -1225 ){
114         Float_t dz = TMath::Abs(-975-x[2])*0.01;
115         b[0] = - (1-0.1*dz*dz)*7;
116         if(fFactor!=1) {
117             b[0]*=fFactor;
118             b[1]*=fFactor;
119             b[2]*=fFactor;
120         }
121       }
122       else {
123           ZDCField(x, b);
124       }
125     }
126
127   } 
128   else {
129       AliFatal(Form("Invalid field map for constant field %d",fMap));
130   }
131 }
132
133 //___________________________________________________
134 void AliMagFC::ZDCField(const float *x, float *b) const
135 {
136   // ---- This is the ZDC part
137   
138   float rad2 = x[0] * x[0] + x[1] * x[1];
139   static Bool_t init = kFALSE;
140
141   if (! init) {
142       init = kTRUE;
143       //////////////////////////////////////////////////////////////////////
144       // ---- Magnetic field values (according to beam type and energy) ----
145       if(fBeamType==kBeamTypepp && fBeamEnergy == 5000.){
146           // p-p @ 5+5 TeV
147           fQuadGradient = 15.7145;
148           fDipoleField  = 27.0558;
149           // SIDE C
150           fCCorrField   = 9.7017;
151           // SIDE A
152           fACorr1Field  = -13.2143;
153           fACorr2Field  = -11.9909;
154       } else if (fBeamType == kBeamTypepp && fBeamEnergy == 450.) {
155           // p-p 0.45+0.45 TeV
156           Float_t const kEnergyRatio = fBeamEnergy / 7000.;
157           
158           fQuadGradient = 22.0002 * kEnergyRatio;
159           fDipoleField  = 37.8781 * kEnergyRatio;
160           // SIDE C
161           fCCorrField   =  9.6908;
162           // SIDE A
163           fACorr1Field  = -13.2014;
164           fACorr2Field  = -9.6908;
165       } else if ((fBeamType == kBeamTypepp && fBeamEnergy == 7000.) ||
166                  (fBeamType == kBeamTypeAA))
167       {
168           // Pb-Pb @ 2.7+2.7 TeV or p-p @ 7+7 TeV
169           fQuadGradient = 22.0002;
170           fDipoleField  = 37.8781;
171           // SIDE C
172           fCCorrField   = 9.6908;
173           // SIDE A
174           fACorr1Field  = -13.2014;
175           fACorr2Field  = -9.6908;
176       }
177   }
178   
179   
180   // SIDE C **************************************************
181   if(x[2]<0.){  
182     if(x[2] < kCCorrBegin && x[2] > kCCorrEnd && rad2 < kCCorrSqRadius){
183         if (fFactor != 0.) {
184             b[0] = fCCorrField;
185             b[1] = 0.;
186             b[2] = 0.;
187         } 
188     }
189     else if(x[2] < kCQ1Begin && x[2] > kCQ1End && rad2 < kCQ1SqRadius){
190         b[0] = fQuadGradient*x[1];
191         b[1] = fQuadGradient*x[0];
192         b[2] = 0.;
193     }
194     else if(x[2] < kCQ2Begin && x[2] > kCQ2End && rad2 < kCQ2SqRadius){
195         b[0] = -fQuadGradient*x[1];
196         b[1] = -fQuadGradient*x[0];
197         b[2] = 0.;
198     }
199     else if(x[2] < kCQ3Begin && x[2] > kCQ3End && rad2 < kCQ3SqRadius){
200         b[0] = -fQuadGradient*x[1];
201         b[1] = -fQuadGradient*x[0];
202         b[2] = 0.;
203     }
204     else if(x[2] < kCQ4Begin && x[2] > kCQ4End && rad2 < kCQ4SqRadius){
205         b[0] = fQuadGradient*x[1];
206         b[1] = fQuadGradient*x[0];
207         b[2] = 0.;
208     }
209     else if(x[2] < kCD1Begin && x[2] > kCD1End && rad2 < kCD1SqRadius){
210         b[1] = fDipoleField;
211         b[2] = 0.;
212         b[2] = 0.;
213     }
214     else if(x[2] < kCD2Begin && x[2] > kCD2End){
215         if(((x[0]-kCD2XCentre1)*(x[0]-kCD2XCentre1)+(x[1]*x[1]))<kCD2SqRadius
216            || ((x[0]-kCD2XCentre2)*(x[0]-kCD2XCentre2)+(x[1]*x[1]))<kCD2SqRadius){
217           b[1] = -fDipoleField;
218           b[2] = 0.;
219           b[2] = 0.;
220         }
221     }
222   }
223   
224   // SIDE A **************************************************
225   else{        
226     if(fCompensator && (x[2] > kACorr1Begin && x[2] < kACorr1End) && rad2 < kCCorr1SqRadius) {
227       // Compensator magnet at z = 1075 m 
228         if (fFactor != 0.) {
229             b[0] = fACorr1Field;
230             b[1] = 0.;
231             b[2] = 0.;
232         }
233         return;
234     }
235     
236     if(x[2] > kACorr2Begin && x[2] < kACorr2End && rad2 < kCCorr2SqRadius){
237         if (fFactor != 0.) {
238             b[0] = fACorr2Field;
239             b[1] = 0.;
240             b[2] = 0.;
241         }
242     }          
243     else if(x[2] > kAQ1Begin && x[2] < kAQ1End && rad2 < kAQ1SqRadius){
244         // First quadrupole of inner triplet de-focussing in x-direction
245         b[0] = -fQuadGradient*x[1];
246         b[1] = -fQuadGradient*x[0];
247         b[2] = 0.;
248     }
249     else if(x[2] > kAQ2Begin && x[2] < kAQ2End && rad2 < kAQ2SqRadius){
250         b[0] = fQuadGradient*x[1];
251         b[1] = fQuadGradient*x[0];
252         b[2] = 0.;
253     }
254     else if(x[2] > kAQ3Begin && x[2] < kAQ3End && rad2 < kAQ3SqRadius){
255         b[0] = fQuadGradient*x[1];
256         b[1] = fQuadGradient*x[0];
257         b[2] = 0.;
258     }
259     else if(x[2] > kAQ4Begin && x[2] < kAQ4End && rad2 < kAQ4SqRadius){
260         b[0] = -fQuadGradient*x[1];
261         b[1] = -fQuadGradient*x[0];
262         b[2] = 0.;
263     }
264     else if(x[2] > kAD1Begin && x[2] < kAD1End && rad2 < kAD1SqRadius){
265         b[0] = 0.;
266         b[1] = -fDipoleField;
267         b[2] = 0.;
268     }
269     else if(x[2] > kAD2Begin && x[2] < kAD2End){
270         if(((x[0]-kAD2XCentre1)*(x[0]-kAD2XCentre1)+(x[1]*x[1])) < kAD2SqRadius
271            || ((x[0]-kAD2XCentre2)*(x[0]-kAD2XCentre2)+(x[1]*x[1])) < kAD2SqRadius){
272             b[1] = fDipoleField;
273         }
274     }
275   }
276 }
277
278 void AliMagFC::ZDCField(const double *x, double *b) const
279 {
280   // ---- This is the ZDC part
281   
282   double rad2 = x[0] * x[0] + x[1] * x[1];
283   static Bool_t init = kFALSE;
284
285   if (! init) {
286       init = kTRUE;
287       //////////////////////////////////////////////////////////////////////
288       // ---- Magnetic field values (according to beam type and energy) ----
289       if(fBeamType==kBeamTypepp && fBeamEnergy == 5000.){
290           // p-p @ 5+5 TeV
291           fQuadGradient = 15.7145;
292           fDipoleField  = 27.0558;
293           // SIDE C
294           fCCorrField   = 9.7017;
295           // SIDE A
296           fACorr1Field  = -13.2143;
297           fACorr2Field  = -11.9909;
298       } else if (fBeamType == kBeamTypepp && fBeamEnergy == 450.) {
299           // p-p 0.45+0.45 TeV
300           Float_t const kEnergyRatio = fBeamEnergy / 7000.;
301           
302           fQuadGradient = 22.0002 * kEnergyRatio;
303           fDipoleField  = 37.8781 * kEnergyRatio;
304           // SIDE C
305           fCCorrField   =  9.6908;
306           // SIDE A
307           fACorr1Field  = -13.2014;
308           fACorr2Field  = -9.6908;
309       } else if ((fBeamType == kBeamTypepp && fBeamEnergy == 7000.) ||
310                  (fBeamType == kBeamTypeAA))
311       {
312           // Pb-Pb @ 2.7+2.7 TeV or p-p @ 7+7 TeV
313           fQuadGradient = 22.0002;
314           fDipoleField  = 37.8781;
315           // SIDE C
316           fCCorrField   = 9.6908;
317           // SIDE A
318           fACorr1Field  = -13.2014;
319           fACorr2Field  = -9.6908;
320       }
321   }
322   
323   
324   // SIDE C **************************************************
325   if(x[2]<0.){  
326     if(x[2] < kCCorrBegin && x[2] > kCCorrEnd && rad2 < kCCorrSqRadius){
327         if (fFactor != 0.) {
328             b[0] = fCCorrField;
329             b[1] = 0.;
330             b[2] = 0.;
331         } 
332     }
333     else if(x[2] < kCQ1Begin && x[2] > kCQ1End && rad2 < kCQ1SqRadius){
334         b[0] = fQuadGradient*x[1];
335         b[1] = fQuadGradient*x[0];
336         b[2] = 0.;
337     }
338     else if(x[2] < kCQ2Begin && x[2] > kCQ2End && rad2 < kCQ2SqRadius){
339         b[0] = -fQuadGradient*x[1];
340         b[1] = -fQuadGradient*x[0];
341         b[2] = 0.;
342     }
343     else if(x[2] < kCQ3Begin && x[2] > kCQ3End && rad2 < kCQ3SqRadius){
344         b[0] = -fQuadGradient*x[1];
345         b[1] = -fQuadGradient*x[0];
346         b[2] = 0.;
347     }
348     else if(x[2] < kCQ4Begin && x[2] > kCQ4End && rad2 < kCQ4SqRadius){
349         b[0] = fQuadGradient*x[1];
350         b[1] = fQuadGradient*x[0];
351         b[2] = 0.;
352     }
353     else if(x[2] < kCD1Begin && x[2] > kCD1End && rad2 < kCD1SqRadius){
354         b[1] = fDipoleField;
355         b[2] = 0.;
356         b[2] = 0.;
357     }
358     else if(x[2] < kCD2Begin && x[2] > kCD2End){
359         if(((x[0]-kCD2XCentre1)*(x[0]-kCD2XCentre1)+(x[1]*x[1]))<kCD2SqRadius
360            || ((x[0]-kCD2XCentre2)*(x[0]-kCD2XCentre2)+(x[1]*x[1]))<kCD2SqRadius){
361           b[1] = -fDipoleField;
362           b[2] = 0.;
363           b[2] = 0.;
364         }
365     }
366   }
367   
368   // SIDE A **************************************************
369   else{        
370     if(fCompensator && (x[2] > kACorr1Begin && x[2] < kACorr1End) && rad2 < kCCorr1SqRadius) {
371       // Compensator magnet at z = 1075 m 
372         if (fFactor != 0.) {
373             b[0] = fACorr1Field;
374             b[1] = 0.;
375             b[2] = 0.;
376         }
377         return;
378     }
379     
380     if(x[2] > kACorr2Begin && x[2] < kACorr2End && rad2 < kCCorr2SqRadius){
381         if (fFactor != 0.) {
382             b[0] = fACorr2Field;
383             b[1] = 0.;
384             b[2] = 0.;
385         }
386     }          
387     else if(x[2] > kAQ1Begin && x[2] < kAQ1End && rad2 < kAQ1SqRadius){
388         // First quadrupole of inner triplet de-focussing in x-direction
389         b[0] = -fQuadGradient*x[1];
390         b[1] = -fQuadGradient*x[0];
391         b[2] = 0.;
392     }
393     else if(x[2] > kAQ2Begin && x[2] < kAQ2End && rad2 < kAQ2SqRadius){
394         b[0] = fQuadGradient*x[1];
395         b[1] = fQuadGradient*x[0];
396         b[2] = 0.;
397     }
398     else if(x[2] > kAQ3Begin && x[2] < kAQ3End && rad2 < kAQ3SqRadius){
399         b[0] = fQuadGradient*x[1];
400         b[1] = fQuadGradient*x[0];
401         b[2] = 0.;
402     }
403     else if(x[2] > kAQ4Begin && x[2] < kAQ4End && rad2 < kAQ4SqRadius){
404         b[0] = -fQuadGradient*x[1];
405         b[1] = -fQuadGradient*x[0];
406         b[2] = 0.;
407     }
408     else if(x[2] > kAD1Begin && x[2] < kAD1End && rad2 < kAD1SqRadius){
409         b[0] = 0.;
410         b[1] = -fDipoleField;
411         b[2] = 0.;
412     }
413     else if(x[2] > kAD2Begin && x[2] < kAD2End){
414         if(((x[0]-kAD2XCentre1)*(x[0]-kAD2XCentre1)+(x[1]*x[1])) < kAD2SqRadius
415            || ((x[0]-kAD2XCentre2)*(x[0]-kAD2XCentre2)+(x[1]*x[1])) < kAD2SqRadius){
416             b[1] = fDipoleField;
417         }
418     }
419   }
420 }
421