StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StarMagField.cxx
1 /***********************************************************************
2  *
3  * $Id: StarMagField.cxx,v 1.34 2020/01/15 02:26:57 perev Exp $
4  *
5  * Author: Jim Thomas 11/1/2000
6  *
7  ***********************************************************************
8  *
9  * Description: the STAR Magnetic Field
10  *
11  ***********************************************************************
12  *
13  * $Log: StarMagField.cxx,v $
14  * Revision 1.34 2020/01/15 02:26:57 perev
15  * Print++
16  *
17  * Revision 1.33 2018/06/29 21:46:25 smirnovd
18  * Revert iTPC-related changes committed on 2018-06-20 through 2018-06-28
19  *
20  * Revert "NoDead option added"
21  * Revert "Fill mag field more carefully"
22  * Revert "Assert commented out"
23  * Revert "Merging with TPC group code"
24  * Revert "Remove too strong assert"
25  * Revert "Restore removed by mistake line"
26  * Revert "Remove not used anymore file"
27  * Revert "iTPCheckIn"
28  *
29  * Revision 1.31 2017/04/28 19:44:35 perev
30  * Fix wrong default. Non const field is default
31  *
32  * Revision 1.30 2017/04/26 21:11:48 perev
33  * Add setConstBz
34  *
35  * Revision 1.29 2014/07/27 13:21:41 fisyak
36  * Add cast for c++11 option
37  *
38  * Revision 1.28 2014/06/27 14:27:11 fisyak
39  * Add switch between new and old schema
40  *
41  * Revision 1.27 2014/06/26 21:50:17 fisyak
42  * New Tpc Alignment, v632
43  *
44  * Revision 1.26 2013/03/26 13:38:18 fisyak
45  * restore back modififcations as not related to drop in no. of reconstructed tracks
46  *
47  * Revision 1.24 2013/02/22 18:40:03 perev
48  * Remove gufld definition
49  *
50  * Revision 1.23 2013/02/22 17:20:30 fisyak
51  * gufld => agufld
52  *
53  * Revision 1.22 2013/01/17 15:11:33 fisyak
54  * More clear handling ROOT and non ROOT versions
55  *
56  * Revision 1.20 2013/01/15 23:45:02 fisyak
57  * Account ROOT version with TVirtualMagField
58  *
59  * Revision 1.19 2013/01/15 17:35:23 fisyak
60  * Create clean versions of ROOT and non ROOT StarMagField
61  *
62  * Revision 1.18 2011/07/21 16:52:10 fisyak
63  * Comment out Lijuan correction which broke B3DField
64  *
65  * Revision 1.17 2010/08/10 19:46:18 fisyak
66  * Lock mag.field if it was initialized from GEANT
67  *
68  * Revision 1.16 2010/05/27 14:52:02 fisyak
69  * Clean up, set assert when mag. field is not initialized
70  *
71  * Revision 1.15 2009/12/07 23:38:15 fisyak
72  * Move size definition from #define to enumerations
73  *
74  * Revision 1.14 2009/11/10 21:18:53 fisyak
75  * use local fMap variable
76  *
77  * Revision 1.13 2009/02/03 15:53:30 fisyak
78  * Clean up
79  *
80  * Revision 1.12 2009/01/26 15:15:56 fisyak
81  * Add missing (ROOT Version >= 5.22)include
82  *
83  * Revision 1.11 2009/01/13 03:19:43 perev
84  * Mag field nou controlled from starsim. BugFix
85  *
86  * Revision 1.10 2007/09/21 21:07:08 fisyak
87  * Remove ClassDef and ClassImp
88  *
89  * Revision 1.9 2007/09/21 17:30:33 perev
90  * Root dependecies removed
91  *
92  * Revision 1.8 2007/09/13 00:00:27 fisyak
93  * add mag.field in steel, from Lijuan Ruan
94  *
95  * Revision 1.7 2005/10/10 19:26:59 fisyak
96  * Ignore default request for StarMagField from mfldgeo
97  *
98  * Revision 1.6 2005/09/12 13:56:27 fisyak
99  * Fix B[0] = B[1] = 0 at r = 0
100  *
101  * Revision 1.5 2005/08/31 15:45:43 fisyak
102  * Make agufld a function to have comis happy (FPE)
103  *
104  * Revision 1.4 2005/08/29 22:59:56 fisyak
105  * Fix printouts
106  *
107  * Revision 1.3 2005/07/28 19:46:01 fisyak
108  * Add:
109  * - frindge magnetic field from P.Nevski extrapolation,
110  * - lock mechanism for Magnetic Field parameters
111  * Remove:
112  * - Electric Field (Still part of StDbUtilitie/StMagUtilities)
113  *
114  * Comparision between mfldgeo and StarMagField calculation for Z and R components of mag.field
115  * is at http://www.star.bnl.gov/~fisyak/star/MagField/
116  *
117  * Revision 1.2 2005/07/15 22:17:37 fisyak
118  * fix misleading print out
119  *
120  * Revision 1.1.1.1 2005/07/07 14:13:55 fisyak
121  * The version of STAR mag. field extracted from StDbUtilities/StMagUtilities to be used in Simulation and Reconstruction instead of agufld
122  *
123  * Revision 1.2 2005/07/07 14:07:55 fisyak
124  * Final version before moving to official repository
125  *
126  * Revision 1.1 2004/03/12 13:26:24 fisyak
127  * Singleton for STAR magnetic field
128  *
129  ***********************************************************************/
131 // //
132 // StarMagField Class //
133 // //
135 
171 #include <string.h>
172 #include <assert.h>
173 #include "StarMagField.h"
174 #include "StarCallf77.h"
175 #include <string>
176 #include "TMath.h"
177 StarMagField *StarMagField::fgInstance = 0;
178 //________________________________________________________________________________
179 
180 #define agufld F77_NAME(agufld,AGUFLD)
181 #define mfldgeo F77_NAME(mfldgeo,MFLDGEO)
182 #ifdef __ROOT__
183 #include "TString.h"
184 #include "TSystem.h"
185 #include "TFile.h"
186 #include "TError.h"
187 #include "TEnv.h"
188 ClassImp(StarMagField);
189 #endif
190 //________________________________________________________________________________
191 StarMagField* StarMagField::Instance() {return fgInstance;}
192 //________________________________________________________________________________
193 R__EXTERN "C" {
194 
195  Float_t type_of_call agufld(Float_t *x, Float_t *bf) {
196  bf[0] = bf[1] = bf[2] = 0;
197  if (StarMagField::Instance())
198  StarMagField::Instance()->BField(x,bf);
199  else {
200  printf("agufld:: request for non initialized mag.field, return 0\n");
201  assert(StarMagField::Instance());
202  }
203  return 0;
204  }
205 //________________________________________________________________________________
206  void type_of_call mfldgeo(Float_t &factor) {
207  if (StarMagField::Instance()) {
208  printf("StarMagField mfldgeo: The field has been already instantiated.\n");
209  } else {
210  printf("StarMagField instantiate starsim field=%g\n",factor);
211  (new StarMagField(StarMagField::kMapped,factor/5.))->SetLock();
212  }
213  Float_t x[3]={0},b[3];
214  agufld(x,b);
215  printf("StarMagField:mfldgeo(%g) Bz=%g\n",factor,b[2]);
216  }
217 }
218 //________________________________________________________________________________
219 struct BFLD_t {
220  Int_t version;
221  const Char_t *code;
222  Float_t date; Int_t kz; Float_t rmaxx, zmaxx, rrm, zz1, zz2;
223  Float_t RmaxInn, ZmaxInn;
224  Int_t nrp, nzp;
225 };
226 static const BFLD_t BFLD = {// real field
227  3 , // version
228  "opt1" , // code: fit version code
229  22.10 , // date: fit date
230  22 , // kz: number of z lines
231  270 , // rmaxx: maximum radius of extrapolated measurements
232  290 , // zmaxx: maximum length of extrapolated measurements
233  400 , // rrm: maximum radius of all fields
234  270 , // zz1: length of measured field interpolation
235  800 , // zz2: max length of all fields
236  264.265 , // RmaxInn: Inner field volume radius
237  312.500 , // ZmaxInn: Inner field volume length
238  200 , // nrp: number of R nodes in the map
239  800 // nzp: number of Z nodes in the map
240 };
241 
242 struct BDAT_t {
243  Int_t N;
244  Float_t Zi, Ri[20], Bzi[20], Bri[20];
245 };
246 
247 static const Int_t nZext = 23;
248 static const BDAT_t BDAT[nZext] = { // calculated STAR field
249  { 15 , // Number of field points
250  0.0 , // distance to Z=0
251  { 0.0, 25.0, 50.0, 75.0, 100.0, 125.0, 150.0, 175.0,
252  200.0, 225.0, 250.0, 275.0, 300.0, 375.0, 400.0 } , // Radius
253  {4704.6,4768.0,4946.6,5148.3,5128.6,4927.3,4844.9,4830.1,
254  4823.3,4858.0,5110.9,3402.4, -18.1, -13.6, -25.0 } , // Axial field
255  { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
256  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 } },// Radial field == 0 @ z = 0
257  // { 0.0, 131.3, 188.4, 74.1,-148.6,-164.8, -53.2, 23.8,
258  // 97.4, 213.7, 329.3, 75.3, 18.2, -44.3, -36.5 } },// Radial field
259  { 15 , // Number of field points
260  270.0 , // distance to Z=0
261  { 0.0, 25.0, 50.0, 75.0, 100.0, 125.0, 150.0, 175.0,
262  200.0, 225.0, 250.0, 275.0, 300.0, 375.0, 400.0 } , // Radius
263  {4704.6,4768.0,4946.6,5148.3,5128.6,4927.3,4844.9,4830.1,
264  4823.3,4858.0,5110.9,3402.4, -18.1, -13.6, -25.0 } , // Axial field
265  { 0.0, 131.3, 188.4, 74.1,-148.6,-164.8, -53.2, 23.8,
266  97.4, 213.7, 329.3, 75.3, 18.2, -44.3, -36.5 } },// Radial field
267  { 15 , // Number of field points
268  280.0 , // distance to Z=0
269  { 0.0, 25.0, 50.0, 75.0, 100.0, 125.0, 150.0, 175.0,
270  200.0, 225.0, 250.0, 275.0, 300.0, 375.0, 400.0 } , // Radius
271  {4568.8,4649.8,4898.6,5241.5,5234.5,4883.9,4806.5,4802.4,
272  4781.7,4771.8,5057.8,3504.2,-144.8, -15.0, -28.0 } , // Axial field
273  { 0.0, 188.6, 297.5, 151.9,-241.2,-242.5, -60.1, 19.5,
274  92.3, 244.3, 541.5, 396.8, 83.6, -49.9, -40.6 } },// Radial field
275  { 15 , // Number of field points
276  290.0 , // distance to Z=0
277  { 0.0, 25.0, 50.0, 75.0, 100.0, 125.0, 150.0, 175.0,
278  200.0, 225.0, 250.0, 275.0, 300.0, 375.0, 400.0 } , // Radius
279  {4383.8,4478.4,4801.1,5378.7,5431.1,4771.4,4765.5,4778.2,
280  4741.4,4651.9,4852.9,3684.4, 6.9, -16.9, -31.6 } , // Axial field
281  { 0.0, 260.2, 456.5, 312.9,-414.5,-349.8, -51.7, 14.4,
282  74.7, 234.0, 858.0, 726.3, 355.0, -56.5, -45.0 } },// Radial field
283  { 15 , // Number of field points
284  300.0 , // distance to Z=0
285  { 0.0, 25.0, 50.0, 75.0, 100.0, 125.0, 150.0, 175.0,
286  200.0, 225.0, 250.0, 275.0, 300.0, 375.0, 400.0 } , // Radius
287  {4142.1,4240.8,4614.8,5546.4,5829.1,4450.0,4737.6,4761.4,
288  4711.3,4534.1,4231.0,4067.5,-880.0, -19.3, -36.2 } , // Axial field
289  { 0.0, 341.1, 669.5, 661.0,-766.7,-480.9, -24.5, 8.8,
290  43.5, 149.9,1333.6, 999.3, 53.6, -64.2, -49.8 } },// Radial field
291  { 15 , // Number of field points
292  310.0 , // distance to Z=0
293  { 0.0, 25.0, 50.0, 75.0, 100.0, 125.0, 150.0, 175.0,
294  200.0, 225.0, 250.0, 275.0, 300.0, 375.0, 400.0 } , // Radius
295  {3842.7,3930.2,4292.5,5589.1,6643.0,3236.8,4733.0,4755.1,
296  4699.4,4485.0,1931.8,4782.0, 50.2, -22.8, -42.0 } , // Axial field
297  { 0.0, 421.2, 915.6,1382.6,-1482.8,-1019.7, 1.2, 2.0,
298  1.9, -2.3,2069.4, 791.7, 240.6, -73.6, -54.9 } },// Radial field
299  { 8 , // Number of field points
300  320.0 , // distance to Z=0
301  { 0.0, 25.0, 50.0, 75.0, 100.0, 125.0, 375.0, 400.0 } , // Radius
302  {3491.2,3552.1,3807.3,4923.7,7889.6,1983.9, -28.0, -49.4 } , // Axial field
303  { 0.0, 485.7,1133.5,2502.8, -38.8,-174.8, -85.1, -60.0 } },// Radial field
304  { 7 , // Number of field points
305  330.0 , // distance to Z=0
306  { 0.0, 25.0, 50.0, 75.0, 250.0, 375.0, 400.0 } , // Radius
307  {3105.3,3127.0,3200.4,3268.9, -3.5, -36.6, -59.0 } , // Axial field
308  { 0.0, 521.1,1246.1,3029.5,9199.2, -99.4, -64.5 } },// Radial field
309  { 6 , // Number of field points
310  340.0 , // distance to Z=0
311  { 0.0, 25.0, 50.0, 75.0, 375.0, 400.0 } , // Radius
312  {2706.4,2686.8,2574.5,1826.7, -51.8, -71.0 } , // Axial field
313  { 0.0, 520.6,1218.1,2485.3,-116.9, -67.3 } },// Radial field
314  { 6 , // Number of field points
315  350.0 , // distance to Z=0
316  { 0.0, 25.0, 50.0, 75.0, 375.0, 400.0 } , // Radius
317  {2317.7,2264.6,2026.3,1142.6, -80.8, -85.1 } , // Axial field
318  { 0.0, 487.6,1082.3,1787.2,-133.8, -67.0 } },// Radial field
319  { 8 , // Number of field points
320  360.0 , // distance to Z=0
321  { 0.0, 25.0, 50.0, 75.0, 100.0, 250.0, 375.0, 400.0 } , // Radius
322  {1958.5,1885.6,1595.6, 829.2,-563.7,4895.8,-127.6, -99.8 } , // Axial field
323  { 0.0, 432.4, 901.7,1265.8, 788.0,9507.4,-134.0, -62.2 } },// Radial field
324  { 17 , // Number of field points
325  370.0 , // distance to Z=0
326  { 0.0, 25.0, 50.0, 75.0, 100.0, 125.0, 150.0, 175.0, 200.0,
327  225.0, 250.0, 275.0, 300.0, 325.0, 350.0, 375.0, 400.0 } , // Radius
328  {1637.8,1562.2,1276.7, 678.9, 15.7, 251.6, 384.9, 503.7, 683.3,
329  1087.1,1868.1,-1320.5,-593.9,-391.5,-345.9,-168.2,-112.9 } , // Axial field
330  { 0.0, 367.3, 720.6, 900.1, 421.6, 60.4, 37.1, 44.5, 79.7,
331  229.6,2339.4, 654.6, 114.6, 35.9, -30.0,-101.8, -52.4 } },// Radial field
332  { 17 , // Number of field points
333  380.0 , // distance to Z=0
334  { 0.0, 25.0, 50.0, 75.0, 100.0, 125.0, 150.0, 175.0, 200.0,
335  225.0, 250.0, 275.0, 300.0, 325.0, 350.0, 375.0, 400.0 } , // Radius
336  {1373.5,1296.6,1045.5, 603.2, 221.4, 278.6, 382.7, 488.2, 638.7,
337  892.4, 708.6,-709.9,-515.0,-364.7,-293.1,-181.5,-122.1 } , // Axial field
338  { 0.0, 302.3, 563.3, 650.3, 369.7, 120.0, 79.6, 96.2, 169.1,
339  430.1,1454.7, 860.7, 228.6, 77.5, -10.8, -60.2, -39.4 } },// Radial field
340  { 17 , // Number of field points
341  390.0 , // distance to Z=0
342  { 0.0, 25.0, 50.0, 75.0, 100.0, 125.0, 150.0, 175.0, 200.0,
343  225.0, 250.0, 275.0, 300.0, 325.0, 350.0, 375.0, 400.0 } , // Radius
344  {1151.2,1083.6, 877.2, 557.6, 308.9, 305.2, 377.6, 463.3, 573.3,
345  684.5, 377.5,-376.2,-415.2,-326.0,-258.2,-179.8,-126.9 } , // Axial field
346  { 0.0, 243.7, 437.7, 486.1, 319.9, 155.1, 115.2, 139.4, 232.4,
347  494.6,1019.1, 751.4, 289.6, 112.2, 19.4, -26.7, -25.0 } },// Radial field
348  { 17 , // Number of field points
349  400.0 , // distance to Z=0
350  { 0.0, 25.0, 50.0, 75.0, 100.0, 125.0, 150.0, 175.0, 200.0,
351  225.0, 250.0, 275.0, 300.0, 325.0, 350.0, 375.0, 400.0 } , // Radius
352  {971.6, 914.8, 751.6, 520.8, 348.0, 323.7, 369.1, 432.0, 500.9,
353  520.4, 251.2,-214.3,-320.8,-282.3,-230.2,-171.7,-127.7 } , // Axial field
354  { 0.0, 194.5, 341.1, 375.8, 277.5, 171.7, 142.1, 172.1, 269.4,
355  486.6, 769.0, 624.1, 308.0, 137.2, 44.9, -1.4, -11.2 } },// Radial field
356  { 9 , // Number of field points
357  450.0 , // distance to Z=0
358  { 0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0, 350.0, 400.0 } , // Radius
359  {481.5, 423.2, 325.1, 283.8, 242.6, 88.2, -85.2,-123.1,-103.9 } , // Axial
360  { 0.0, 119.3, 157.5, 174.6, 248.4, 314.8, 220.8, 95.4, 32.3 } },// Radial
361  { 9 , // Number of field points
362  500.0 , // distance to Z=0
363  { 0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0, 350.0, 400.0 } , // Radius
364  {291.3, 273.2, 234.4, 192.6, 136.7, 53.6, -25.6, -64.6, -72.5 } , // Axial
365  { 0.0, 60.7, 103.2, 135.9, 168.1, 177.4, 140.6, 84.7, 41.9 } },// Radial
366  { 9 , // Number of field points
367  550.0 , // distance to Z=0
368  { 0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0, 350.0, 400.0 } , // Radius
369  {190.4, 181.9, 159.6, 127.7, 86.0, 37.1, -7.3, -35.9, -48.9 } , // Axial
370  { 0.0, 37.8, 69.7, 94.3, 110.3, 110.8, 92.6, 64.2, 37.3 } },// Radial
371  { 9 , // Number of field points
372  600.0 , // distance to Z=0
373  { 0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0, 350.0, 400.0 } , // Radius
374  {126.9, 121.7, 107.1, 84.9, 56.8, 26.4, -1.2, -21.2, -32.9 } , // Axial
375  { 0.0, 25.0, 46.9, 63.4, 72.6, 72.2, 62.3, 46.3, 29.2 } },// Radial
376  { 9 , // Number of field points
377  650.0 , // distance to Z=0
378  { 0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0, 350.0, 400.0 } , // Radius
379  { 84.8, 81.4, 71.7, 56.8, 38.3, 18.7, 0.8, -13.1, -22.2 } , // Axial
380  { 0.0, 16.7, 31.4, 42.4, 48.2, 48.1, 42.4, 32.8, 21.6 } },// Radial
381  { 9 , // Number of field points
382  700.0 , // distance to Z=0
383  { 0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0, 350.0, 400.0 } , // Radius
384  { 56.7, 54.4, 47.9, 38.0, 25.9, 13.1, 1.3, -8.3, -15.0 } , // Axial
385  { 0.0, 11.2, 21.1, 28.4, 32.4, 32.5, 29.1, 23.0, 15.5 } },// Radial
386  { 9 , // Number of field points
387  750.0 , // distance to Z=0
388  { 0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0, 350.0, 400.0 } , // Radius
389  { 37.7, 36.2, 31.9, 25.4, 17.5, 9.1, 1.2, -5.4, -10.1 } , // Axial
390  { 0.0, 7.6, 14.3, 19.3, 22.0, 22.2, 20.1, 16.1, 11.0 } },// Radial
391  { 9 , // Number of field points
392  800.0 , // distance to Z=0
393  { 0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0, 350.0, 400.0 } , // Radius
394  { 24.8, 23.8, 21.0, 16.8, 11.6, 6.1, 0.9, -3.5, -6.7 } , // Axial
395  { 0.0, 5.2, 9.8, 13.3, 15.2, 15.4, 14.1, 11.4, 7.9 } },// Radial
396 };
397 #ifdef __ROOT__
398 //________________________________________________________________________________
399 void StarMagField::SetStarMagFieldRotation(TGeoRotation &rot) {
400  fStarMagFieldRotation = rot;
401  fStarMagFieldRotation.SetName("StarMagFieldRotation");
402  fStarMagFieldRotation.Print();
403 }
404 //________________________________________________________________________________
405 void StarMagField::SetStarMagFieldRotation(Double_t *r) {
406  TGeoRotation rot;
407  rot.SetMatrix(r);
408  SetStarMagFieldRotation(rot);
409 }
410 #endif
411 
412 bool StarMagField::mConstBz = false;
413 
414 //________________________________________________________________________________
415 StarMagField::StarMagField ( EBField map, Float_t factor,
416  Bool_t lock, Float_t rescale,
417  Float_t BDipole, Float_t RmaxDip,
418  Float_t ZminDip, Float_t ZmaxDip) :
419 #ifdef __ROOT__
420 #if ROOT_VERSION_CODE >= 335360 /* ROOT_VERSION(5,30,0) */
421  TVirtualMagField("StarMagField"),
422 #endif
423  fBzdZCorrection(0),
424  fBrdZCorrection(0),
425 #endif
426  fMap(map),
427  fFactor(factor), fRescale(rescale),
428  fBDipole(BDipole), fRmaxDip(RmaxDip),
429  fZminDip(ZminDip), fZmaxDip(ZmaxDip),
430  fLock(lock)
431 {
432  if (fgInstance) {
433  printf("Cannot initialise twice StarMagField class\n");
434  assert(0);
435  }
436  fgInstance = this;
437  if (fMap == kUndefined) {
438  printf("StarMagField is instantiated with predefined factor %f and map %i\n",fFactor, fMap);
439  } else {
440  if (fLock) printf("StarMagField is locked, no modification from DB will be accepted\n");
441  }
442  ReadField() ; // Read the Magnetic
443  float myX[3]={0},myB[3];
444  BField(myX,myB);
445  printf ("StarMagField(0,0,0) = %g",myB[2]);
446 
447 #ifdef __ROOT__
448  fStarMagFieldRotation = TGeoRotation("StarMagFieldRotation");
449 #endif /* __ROOT__ */
450 }
451 //________________________________________
453 void StarMagField::BField( const Double_t x[], Double_t B[] ) {
454  Float_t xx[3] = {(Float_t) x[0], (Float_t) x[1], (Float_t) x[2]};
455  Float_t bb[3];
456  BField(xx,bb);
457  B[0] = bb[0]; B[1] = bb[1]; B[2] = bb[2];
458 }
459 //________________________________________________________________________________
460 void StarMagField::BField( const Float_t x[], Float_t B[] )
461 
462 {
463 
464 
465  Float_t r, z, Br_value, Bz_value ;
466  Float_t phi, Bphi_value,phi1;
467  Bphi_value=0;
468  Br_value = Bz_value = 0;
469  B[0] = B[1] = B[2] = 0;
470  z = x[2] ;
471  r = sqrt( x[0]*x[0] + x[1]*x[1] ) ;
472  phi = atan2( x[1], x[0] ) ;
473  if ( phi < 0 ) phi += 2*TMath::Pi() ; // Table uses phi from 0 to 2*Pi
474 
475 
476 
477  if ( mConstBz )
478  {
479  B[0] = B[1] = B[2] = 0.;
480  if ( abs(z) < 380.0 && r < 300.0 ) B[2] = +5.0;
481  return;
482  }
483 
484 
485 
486  Float_t za = fabs(z);
487  if (za > fZminDip && za < fZmaxDip && r < fRmaxDip) {// beam Dipole
488  B[1] = TMath::Sign(fBDipole, z);
489  B[2] = fabs(B[1]/1000.);
490  return;
491  }
492  if (z >= ZList[0] && z <= ZList[nZ-1] && r <= Radius[nR-1]) { // within Map
493  Interpolate2DBfield( r, z, Br_value, Bz_value ) ;
494  Double_t BL[3] = {0, 0, Bz_value};
495  if ( r != 0.0 ) {
496  BL[0] = Br_value * (x[0]/r) ;
497  BL[1] = Br_value * (x[1]/r) ;
498  }
499 #ifdef __ROOT__
500  Double_t BG[3];
501  fStarMagFieldRotation.LocalToMaster(BL,BG);
502  for (Int_t i = 0; i < 3; i++) B[i] = BG[i];
503 #else /* ! __ROOT__ */
504  for (Int_t i = 0; i < 3; i++) B[i] = BL[i];
505 #endif /* __ROOT__ */
506  return;
507  }
508 
509 // //added by Lijuan within the steel
510 
511 
512  if (za <=342.20 && r>=303.29 && r <= 364.25) { // within Map
513 
514  phi1=phi*TMath::RadToDeg();
515  if(phi1>12) phi1=phi1-int(phi1/12)*12;
516 
517  Interpolate3DBSteelfield( r, za, phi1, Br_value, Bz_value, Bphi_value ) ;
518  B[0] = Br_value * (x[0]/r) - Bphi_value * (x[1]/r) ;
519  B[1] = Br_value * (x[1]/r) + Bphi_value * (x[0]/r) ;
520  B[2] = Bz_value ;
521 
522  if(z<0) {
523  B[0]=-B[0];
524  B[1]=-B[1];
525  }
526 
527  //cout<<B[0]<<" --- "<<B[1]<<" --- "<<B[2]<<" --- "<<endl;
528 
529  return;
530  }
531 
532 //end added by Lijuan within the steel
533 
534 
535 // cout<<" <<<<<<<<<<<<<<<<<<<<<<<<<<<debug--- "<<endl;
536 
537  Interpolate2ExtDBfield( r, z, Br_value, Bz_value ) ;
538  if (za <= BFLD.zmaxx && r <= BFLD.rmaxx) {
539  static const Float_t zero = 0;
540  static const Float_t one = 1;
541  Float_t wz = (za - ZList[nZ-1] )/(BFLD.zmaxx - ZList[nZ-1]);
542  Float_t wr = (r - Radius[nR-1])/(BFLD.rmaxx - Radius[nR-1]);
543  Float_t w = TMath::Min(TMath::Max(zero,TMath::Max(wz,wr)),one);
544  Float_t rm = TMath::Min(r,Radius[nR-1]);
545  Float_t zm = TMath::Sign(TMath::Min(za,ZList[nZ-1]),z);
546  Float_t BrI, BzI;
547  Interpolate2DBfield( rm, zm, BrI, BzI ) ;
548  Br_value = (1-w)*BrI + w*Br_value;
549  Bz_value = (1-w)*BzI + w*Bz_value;
550  }
551  B[2] = Bz_value ;
552  if ( r != 0.0 ) {
553  B[0] = Br_value * (x[0]/r) ;
554  B[1] = Br_value * (x[1]/r) ;
555  }
556 
557  // cout<<"r=== "<<r<<" z=== "<<z<<" phi=== "<<phi<<endl;
558  return;
559 }
560 //________________________________________________________________________________
562 void StarMagField::B3DField( const Float_t x[], Float_t B[] )
563 {
564  Float_t r, z, phi, Br_value, Bz_value, Bphi_value ;
565  Bphi_value=0;
566  Br_value = Bz_value = 0;
567  B[0] = B[1] = B[2] = 0;
568 #if 0
569  Float_t phi1;
570 #endif
571  z = x[2] ;
572  r = sqrt( x[0]*x[0] + x[1]*x[1] ) ;
573 
574  if ( r != 0.0 )
575  {
576  phi = TMath::ATan2( x[1], x[0] ) ;
577  if ( phi < 0 ) phi += 2*TMath::Pi() ; // Table uses phi from 0 to 2*Pi
578 #if 0
579  //added by Lijuan
580  phi1=phi*TMath::RadToDeg();
581  //added by Lijuan
582 
583  Interpolate3DBfield( r, z, phi1, Br_value, Bz_value, Bphi_value ) ;
584 #else
585  Interpolate3DBfield( r, z, phi, Br_value, Bz_value, Bphi_value ) ;
586 #endif
587  B[0] = Br_value * (x[0]/r) - Bphi_value * (x[1]/r) ;
588  B[1] = Br_value * (x[1]/r) + Bphi_value * (x[0]/r) ;
589  B[2] = Bz_value ;
590  }
591  else
592  {
593  phi = 0 ;
594  Interpolate3DBfield( r, z, phi, Br_value, Bz_value, Bphi_value ) ;
595  B[0] = Br_value ;
596  B[1] = Bphi_value ;
597  B[2] = Bz_value ;
598  }
599  Double_t BL[3] = {B[0], B[1], B[2]};
600 #ifdef __ROOT__
601  Double_t BG[3];
602  fStarMagFieldRotation.LocalToMaster(BL,BG);
603  for (Int_t i = 0; i < 3; i++) B[i] = BG[i];
604 #else
605  for (Int_t i = 0; i < 3; i++) B[i] = BL[i];
606 #endif
607  return ;
608 
609 }
610 void StarMagField::B3DField( const Double_t x[], Double_t B[] ) {
611  Float_t xx[3] = {(Float_t) x[0], (Float_t) x[1], (Float_t) x[2]};
612  Float_t bb[3];
613  B3DField(xx,bb);
614  B[0] = bb[0]; B[1] = bb[1]; B[2] = bb[2];
615 }
616 
618 
619 void StarMagField::BrBzField( const Float_t r, const Float_t z, Float_t &Br_value, Float_t &Bz_value )
620 
621 {
622 
623 
624  Br_value = Bz_value = 0;
625  if(r>0) Interpolate2DBfield( r, z, Br_value, Bz_value ) ;
626  return;
627 
628 }
629 
630 
631 // /// B field in Radial coordinates - 3D field
632 
633 void StarMagField::BrBz3DField( const Float_t r, const Float_t z, const Float_t phi,
634  Float_t &Br_value, Float_t &Bz_value, Float_t &Bphi_value )
635 
636 {
637 
638  Bphi_value=0;
639  Br_value = Bz_value = 0;
640 
641 #if 0
642 
643  Float_t phiprime ;
644 
645  phiprime = phi ;
646  if ( phiprime < 0 ) phiprime += 2*TMath::Pi() ; // Table uses phi from 0 to 2*Pi
647  //added by Lijuan
648  phiprime=phiprime*TMath::RadToDeg();
649  //added by Lijuan
650 #endif
651 
652  if(r>0) {
653 #if 0
654  Interpolate3DBfield( r, z, phiprime, Br_value, Bz_value, Bphi_value ) ;
655 #else
656  Interpolate3DBfield( r, z, phi, Br_value, Bz_value, Bphi_value ) ;
657 #endif
658  }
659 
660  return;
661 
662 }
663 
664 
665 //________________________________________
666 
668 
669 void StarMagField::ReadField( )
670 
671 {
672  FILE *magfile, *b3Dfile ;
673  std::string comment, filename, filename3D ;
674  std::string MapLocation ;
675  std::string BaseLocation = getenv("STAR") ; // Base Directory for Maps
676  BaseLocation += "/StarDb/StMagF/" ; // Base Directory for Maps
677 #ifdef __ROOT__
678  if (gEnv->GetValue("NewTpcAlignment",0) != 0) {
679  TString rootf("StarFieldZ.root");
680  TString path(".:./StarDb/StMagF:$STAR/StarDb/StMagF");
681  Char_t *file = gSystem->Which(path,rootf,kReadPermission);
682  if (! file) {
683  Error("StarMagField::ReadField","File %s has not been found in path %s",rootf.Data(),path.Data());
684  } else {
685  Warning("StarMagField::ReadField","File %s has been found",rootf.Data());
686  TFile *pFile = new TFile(file);
687  TH2F *Br0 = (TH2F *) pFile->Get("Br0");
688  TH2F *Bz0 = (TH2F *) pFile->Get("Bz0");
689  if (Br0 && Bz0) {
690  TH2F *Br5cm = (TH2F *) pFile->Get("Br5cm");
691  TH2F *Bz5cm = (TH2F *) pFile->Get("Bz5cm");
692  assert(Br5cm && Bz5cm);
693  TH2F *Br10cm = (TH2F *) pFile->Get("Br10cm");
694  TH2F *Bz10cm = (TH2F *) pFile->Get("Bz10cm");
695  assert(Br10cm && Bz10cm);
696  fBzdZCorrection = new TH2F(*Bz5cm); fBzdZCorrection->SetDirectory(0);
697  fBzdZCorrection->Scale(0.5);
698  fBzdZCorrection->Add(Bz10cm,0.5);
699  fBzdZCorrection->Add(Bz0,-1.0);
700  fBrdZCorrection = new TH2F(*Br5cm); fBrdZCorrection->SetDirectory(0);
701  fBrdZCorrection->Scale(0.5);
702  fBrdZCorrection->Add(Br10cm,0.5);
703  fBrdZCorrection->Add(Br0,-1.0);
704  Warning("StarMagField::ReadField","Use effective PMT box dZ = 7.5 cm");
705  }
706  delete pFile;
707  }
708  delete [] file;
709  }
710 #endif
711  if ( fMap == kMapped ) // Mapped field values
712  {
713  if ( fabs(fFactor) > 0.8 ) // Scale from full field data
714  {
715  if ( fFactor > 0 )
716  {
717  filename = "bfield_full_positive_2D.dat" ;
718  filename3D = "bfield_full_positive_3D.dat" ;
719  comment = "Measured Full Field" ;
720  fRescale = 1 ; // Normal field
721  }
722  else
723  {
724  filename = "bfield_full_negative_2D.dat" ;
725  filename3D = "bfield_full_negative_3D.dat" ;
726  comment = "Measured Full Field Reversed" ;
727  fRescale = -1 ; // Reversed field
728  }
729  }
730  else // Scale from half field data
731  {
732  filename = "bfield_half_positive_2D.dat" ;
733  filename3D = "bfield_half_positive_3D.dat" ;
734  comment = "Measured Half Field" ;
735  fRescale = 2 ; // Adjust scale factor to use half field data
736  }
737  }
738  else if ( fMap == kConstant ) // Constant field values
739  {
740  filename = "const_full_positive_2D.dat" ;
741  comment = "Constant Full Field" ;
742  fRescale = 1 ; // Normal field
743  }
744  else
745  {
746  fprintf(stderr,"StarMagField::ReadField No map available - you must choose a mapped field or a constant field\n");
747  exit(1) ;
748  }
749 
750  printf("StarMagField::ReadField Reading Magnetic Field %s, Scale factor = %f \n",comment.c_str(),fFactor);
751  printf("StarMagField::ReadField Filename is %s, Adjusted Scale factor = %f \n",filename.c_str(),fFactor*fRescale);
752 
753  MapLocation = BaseLocation + filename ;
754  magfile = fopen(MapLocation.c_str(),"r") ;
755  printf("StarMagField::ReadField Reading 2D Magnetic Field file: %s \n",filename.c_str());
756 
757  if (magfile)
758 
759  {
760  Char_t cname[128] ;
761  fgets ( cname, sizeof(cname) , magfile ) ; // Read comment lines at begining of file
762  fgets ( cname, sizeof(cname) , magfile ) ;
763  fgets ( cname, sizeof(cname) , magfile ) ;
764  fgets ( cname, sizeof(cname) , magfile ) ;
765  fgets ( cname, sizeof(cname) , magfile ) ;
766 
767  for ( Int_t j=0 ; j < nZ ; j++ )
768  {
769  for ( Int_t k=0 ; k < nR ; k++ )
770  {
771  fgets ( cname, sizeof(cname) , magfile ) ;
772  sscanf ( cname, " %f %f %f %f ", &Radius[k], &ZList[j], &Br[j][k], &Bz[j][k] ) ;
773 #if defined(__ROOT__)
774  if (fBzdZCorrection && fBrdZCorrection) {
775  Br[j][k] += fFactor*fBrdZCorrection->Interpolate(ZList[j],Radius[k]);
776  Bz[j][k] += fFactor*fBzdZCorrection->Interpolate(ZList[j],Radius[k]);
777  }
778 #endif
779  }
780  }
781  }
782 
783  else
784 
785  {
786  fprintf(stderr,"StarMagField::ReadField File %s not found !\n",MapLocation.c_str());
787  exit(1);
788  }
789 
790  fclose(magfile) ;
791 
792  MapLocation = BaseLocation + filename3D ;
793  b3Dfile = fopen(MapLocation.c_str(),"r") ;
794  printf("StarMagField::ReadField Reading 3D Magnetic Field file: %s \n",filename3D.c_str());
795 
796  if (b3Dfile)
797 
798  {
799  Char_t cname[128] ;
800  fgets ( cname, sizeof(cname) , b3Dfile ) ; // Read comment lines at begining of file
801  fgets ( cname, sizeof(cname) , b3Dfile ) ; // Read comment lines at begining of file
802  fgets ( cname, sizeof(cname) , b3Dfile ) ; // Read comment lines at begining of file
803  fgets ( cname, sizeof(cname) , b3Dfile ) ; // Read comment lines at begining of file
804  fgets ( cname, sizeof(cname) , b3Dfile ) ; // Read comment lines at begining of file
805  fgets ( cname, sizeof(cname) , b3Dfile ) ; // Read comment lines at begining of file
806 
807  for ( Int_t i=0 ; i < nPhi ; i++ )
808  {
809  for ( Int_t j=0 ; j < nZ ; j++ )
810  {
811  for ( Int_t k=0 ; k < nR ; k++ )
812  {
813  fgets ( cname, sizeof(cname) , b3Dfile ) ;
814  sscanf ( cname, " %f %f %f %f %f %f ",
815  &R3D[k], &Z3D[j], &Phi3D[i], &Br3D[i][j][k], &Bz3D[i][j][k], &Bphi3D[i][j][k] ) ;
816  Phi3D[i] *= TMath::Pi() / 180. ; // Convert to Radians phi = 0 to 2*Pi
817 #if defined(__ROOT__)
818  if (fBzdZCorrection && fBrdZCorrection) {
819  Br3D[i][j][k] += fFactor*fBrdZCorrection->Interpolate(Z3D[j],R3D[k]);
820  Bz3D[i][j][k] += fFactor*fBzdZCorrection->Interpolate(Z3D[j],R3D[k]);
821  }
822 #endif
823  }
824  }
825  }
826  }
827 
828  else if ( fMap == kConstant ) // Constant field values
829 
830  {
831  for ( Int_t i=0 ; i < nPhi ; i++ )
832  {
833  for ( Int_t j=0 ; j < nZ ; j++ )
834  {
835  for ( Int_t k=0 ; k < nR ; k++ )
836  {
837  Br3D[i][j][k] = Br[j][k] ;
838  Bz3D[i][j][k] = Bz[j][k] ;
839  Bphi3D[i][j][k] = 0 ;
840  }
841  }
842  }
843  }
844 
845  else
846 
847  {
848  fprintf(stderr,"StarMagField::ReadField File %s not found !\n",MapLocation.c_str());
849  exit(1);
850  }
851 
852  fclose(b3Dfile) ;
853  #if 1
854  //cout<<"---------------"<<endl;
855 // memset(R3DSteel, 0, nRSteel*sizeof(Float_t));
856 // memset(Z3DSteel, 0, nZSteel*sizeof(Float_t));
857 // memset(Phi3DSteel, 0, nPhiSteel*sizeof(Float_t));
858 // memset(Bx3DSteel, 0, nPhiSteel*nZSteel*nRSteel*sizeof(Float_t));
859 // memset(By3DSteel, 0, nPhiSteel*nZSteel*nRSteel*sizeof(Float_t));
860 // memset(Bz3DSteel, 0, nPhiSteel*nZSteel*nRSteel*sizeof(Float_t));
861  MapLocation = BaseLocation + "steel_magfieldmap.dat";
862  magfile = fopen(MapLocation.c_str(),"r") ;
863  if (magfile) {
864  printf("StarMagField::ReadField Reading 3D Magnetic Field file: %s \n",filename.c_str());
865  Char_t cname[128] ;
866  for (;;) {
867  fgets ( cname, sizeof(cname) , magfile ) ; // Read comment lines at begining of file
868  if (cname[0] == '#') continue;
869  break;
870  }
871  for ( Int_t i=0 ; i < nPhiSteel ; i++ )
872  {
873 
874  for ( Int_t k=0 ; k < nRSteel ; k++ )
875  {
876  for ( Int_t j=0 ; j < nZSteel ; j++ )
877  {
878 
879 
880  fgets ( cname, sizeof(cname) , magfile ) ;
881  sscanf ( cname, " %f %f %f %f %f %f ",
882  &R3DSteel[k], &Z3DSteel[j], &Phi3DSteel[i], &Bx3DSteel[i][j][k], &Bz3DSteel[i][j][k], &By3DSteel[i][j][k] ) ;
883 
884  //added by Lijuan
885  Br3DSteel[i][j][k]=cos(Phi3DSteel[i]*TMath::DegToRad())*Bx3DSteel[i][j][k]+sin(Phi3DSteel[i]*TMath::DegToRad())*By3DSteel[i][j][k];
886 
887  Bphi3DSteel[i][j][k]=0-sin(Phi3DSteel[i]*TMath::DegToRad())*Bx3DSteel[i][j][k]+cos(Phi3DSteel[i]*TMath::DegToRad())*By3DSteel[i][j][k];
888 
889 
890  //cout<<R3DSteel[k]<<" "<<Z3DSteel[j]<<" "<<Phi3DSteel[i]<<" "<<Bx3DSteel[i][j][k]<<" "<<Bz3DSteel[i][j][k]<<" "<<By3DSteel[i][j][k]<<endl;
891 
892  //end added by Lijuan
893 
894  //cout<<Br3DSteel[i][j][k]<<"--------------------"<<Bphi3DSteel[i][j][k]<<endl;
895 
896  }
897  }
898  }
899  fclose(magfile);
900  }
901  #endif
902 #if 1
903 #endif
904  return ;
905 
906 }
907 
908 
909 //________________________________________
910 
912 
913 void StarMagField::Interpolate2DBfield( const Float_t r, const Float_t z, Float_t &Br_value, Float_t &Bz_value )
914 
915 {
916 
917  Float_t fscale ;
918 
919  fscale = 0.001*fFactor*fRescale ; // Scale STAR maps to work in kGauss, cm
920 
921 
922  const Int_t ORDER = 1 ; // Linear interpolation = 1, Quadratic = 2
923  static Int_t jlow=0, klow=0 ;
924  Float_t save_Br[ORDER+1] ;
925  Float_t save_Bz[ORDER+1] ;
926 
927  Search ( nZ, ZList, z, jlow ) ;
928  Search ( nR, Radius, r, klow ) ;
929  if ( jlow < 0 ) jlow = 0 ; // artifact of Root's binsearch, returns -1 if out of range
930  if ( klow < 0 ) klow = 0 ;
931  if ( jlow + ORDER >= nZ - 1 ) jlow = nZ - 1 - ORDER ;
932  if ( klow + ORDER >= nR - 1 ) klow = nR - 1 - ORDER ;
933 
934  for ( Int_t j = jlow ; j < jlow + ORDER + 1 ; j++ )
935  {
936  save_Br[j-jlow] = Interpolate( &Radius[klow], &Br[j][klow], ORDER, r ) ;
937  save_Bz[j-jlow] = Interpolate( &Radius[klow], &Bz[j][klow], ORDER, r ) ;
938  }
939  Br_value = fscale * Interpolate( &ZList[jlow], save_Br, ORDER, z ) ;
940  Bz_value = fscale * Interpolate( &ZList[jlow], save_Bz, ORDER, z ) ;
941 
942 }
943 //________________________________________________________________________________
944 void StarMagField::Interpolate2ExtDBfield( const Float_t r, const Float_t z, Float_t &Br_value, Float_t &Bz_value ) {
945  static Float_t ZExtList[nZext];
946  static Bool_t first = kTRUE;
947  if (first) {
948  for (Int_t j = 0; j < nZext; j++) ZExtList[j] = BDAT[j].Zi;
949  first = kFALSE;
950  }
951  Float_t za = fabs(z);
952  if (za > BFLD.zz2 || r > BFLD.rrm) return;
953  if (za < ZList[nZ-1] && r < Radius[nR-1]) return;
954 
955  //added by Lijuan
956  if (za <=342.20 && r>=303.29 && r <= 363.29) return;
957  //end added by Lijuan
958 
959 
960  Float_t fscale = 0.001*fFactor;// Scale STAR maps to work in kGauss, cm. Table only for Full Field, no Rescale !
961 
962  const Int_t ORDER = 1 ; // Linear interpolation = 1, Quadratic = 2
963  static Int_t jlow=0, klow=0 ;
964  Float_t save_Br[ORDER+1] ;
965  Float_t save_Bz[ORDER+1] ;
966  Search ( nZext, ZExtList, za, jlow ) ;
967  if ( jlow < 0 ) jlow = 0 ; // artifact of Root's binsearch, returns -1 if out of range
968  if ( jlow + ORDER >= nZext - 1 ) jlow = nZext - 1 - ORDER ;
969 
970  for ( Int_t j = jlow ; j < jlow + ORDER + 1 ; j++ ) {
971  Int_t N = BDAT[j].N;
972  Search ( N, (Float_t *) (&BDAT[j].Ri[0]), r, klow ) ;
973  if ( klow < 0 ) klow = 0 ;
974  if ( klow + ORDER >= BDAT[j].N - 1 ) klow = BDAT[j].N - 1 - ORDER ;
975  save_Br[j-jlow] = Interpolate( &BDAT[j].Ri[klow], &BDAT[j].Bri[klow], ORDER, r ) ;
976  save_Bz[j-jlow] = Interpolate( &BDAT[j].Ri[klow], &BDAT[j].Bzi[klow], ORDER, r ) ;
977  }
978  Br_value = fscale * Interpolate( &ZExtList[jlow], save_Br, ORDER, za ) ;
979  Bz_value = fscale * Interpolate( &ZExtList[jlow], save_Bz, ORDER, za ) ;
980  if (z < 0) Br_value = - Br_value;
981 }
982 
984 
985 void StarMagField::Interpolate3DBfield( const Float_t r, const Float_t z, const Float_t phi,
986  Float_t &Br_value, Float_t &Bz_value, Float_t &Bphi_value )
987 {
988 
989  Float_t fscale ;
990 
991  fscale = 0.001*fFactor*fRescale ; // Scale STAR maps to work in kGauss, cm
992 
993  const Int_t ORDER = 1 ; // Linear interpolation = 1, Quadratic = 2
994  static Int_t ilow=0, jlow=0, klow=0 ;
995  Float_t save_Br[ORDER+1], saved_Br[ORDER+1] ;
996  Float_t save_Bz[ORDER+1], saved_Bz[ORDER+1] ;
997  Float_t save_Bphi[ORDER+1], saved_Bphi[ORDER+1] ;
998 
999 
1000 
1001  //cout<<"r=== "<<r<<" z=== "<<z<<" phi=== "<<phi<<endl;
1002  if(r<0) return;
1003 
1004  Search( nPhi, Phi3D, phi, ilow ) ;
1005  Search( nZ, Z3D, z, jlow ) ;
1006  Search( nR, R3D, r, klow ) ;
1007  if ( ilow < 0 ) ilow = 0 ; // artifact of Root's binsearch, returns -1 if out of range
1008  if ( jlow < 0 ) jlow = 0 ;
1009  if ( klow < 0 ) klow = 0 ;
1010 
1011  if ( ilow + ORDER >= nPhi - 1 ) ilow = nPhi - 1 - ORDER ;
1012  if ( jlow + ORDER >= nZ - 1 ) jlow = nZ - 1 - ORDER ;
1013  if ( klow + ORDER >= nR - 1 ) klow = nR - 1 - ORDER ;
1014 
1015  for ( Int_t i = ilow ; i < ilow + ORDER + 1 ; i++ )
1016  {
1017  for ( Int_t j = jlow ; j < jlow + ORDER + 1 ; j++ )
1018  {
1019  save_Br[j-jlow] = Interpolate( &R3D[klow], &Br3D[i][j][klow], ORDER, r ) ;
1020  save_Bz[j-jlow] = Interpolate( &R3D[klow], &Bz3D[i][j][klow], ORDER, r ) ;
1021  save_Bphi[j-jlow] = Interpolate( &R3D[klow], &Bphi3D[i][j][klow], ORDER, r ) ;
1022  }
1023  saved_Br[i-ilow] = Interpolate( &Z3D[jlow], save_Br, ORDER, z ) ;
1024  saved_Bz[i-ilow] = Interpolate( &Z3D[jlow], save_Bz, ORDER, z ) ;
1025  saved_Bphi[i-ilow] = Interpolate( &Z3D[jlow], save_Bphi, ORDER, z ) ;
1026  }
1027  Br_value = fscale * Interpolate( &Phi3D[ilow], saved_Br, ORDER, phi ) ;
1028  Bz_value = fscale * Interpolate( &Phi3D[ilow], saved_Bz, ORDER, phi ) ;
1029  Bphi_value = fscale * Interpolate( &Phi3D[ilow], saved_Bphi, ORDER, phi ) ;
1030 
1031 }
1032 
1033 
1034 
1035 
1036 
1037 
1038 
1039 
1040 //added by Lijuan for the magnetic field in steel.
1041 
1042 
1044 
1045 void StarMagField::Interpolate3DBSteelfield( const Float_t r, const Float_t z, const Float_t phi,
1046  Float_t &Br_value, Float_t &Bz_value, Float_t &Bphi_value )
1047 {
1048 
1049  Float_t fscale ;
1050 
1051  //This is different from the usual bfield map, changed by Lijuan
1052 
1053  // fscale = 0.001*fFactor*fRescale ; // Scale STAR maps to work in kGauss, cm
1054  fscale = 0.001*fFactor; // Scale STAR maps to work in kGauss, cm
1055 
1056  const Int_t ORDER = 1 ; // Linear interpolation = 1, Quadratic = 2
1057  static Int_t ilow=0, jlow=0, klow=0 ;
1058  Float_t save_Br[ORDER+1], saved_Br[ORDER+1] ;
1059  Float_t save_Bz[ORDER+1], saved_Bz[ORDER+1] ;
1060  Float_t save_Bphi[ORDER+1], saved_Bphi[ORDER+1] ;
1061  // phi=phi+1;
1062 
1063  Search( nPhiSteel, Phi3DSteel, phi, ilow ) ;
1064  Search( nZSteel, Z3DSteel, z, jlow ) ;
1065  Search( nRSteel, R3DSteel, r, klow ) ;
1066  if ( ilow < 0 ) ilow = 0 ; // artifact of Root's binsearch, returns -1 if out of range
1067  if ( jlow < 0 ) jlow = 0 ;
1068  if ( klow < 0 ) klow = 0 ;
1069 
1070  if ( ilow + ORDER >= nPhiSteel - 1 ) ilow = nPhiSteel - 1 - ORDER ;
1071  if ( jlow + ORDER >= nZSteel - 1 ) jlow = nZSteel - 1 - ORDER ;
1072  if ( klow + ORDER >= nRSteel - 1 ) klow = nRSteel - 1 - ORDER ;
1073 
1074  for ( Int_t i = ilow ; i < ilow + ORDER + 1 ; i++ )
1075  {
1076  for ( Int_t j = jlow ; j < jlow + ORDER + 1 ; j++ )
1077  {
1078  save_Br[j-jlow] = Interpolate( &R3DSteel[klow], &Br3DSteel[i][j][klow], ORDER, r ) ;
1079  save_Bz[j-jlow] = Interpolate( &R3DSteel[klow], &Bz3DSteel[i][j][klow], ORDER, r ) ;
1080  save_Bphi[j-jlow] = Interpolate( &R3DSteel[klow], &Bphi3DSteel[i][j][klow], ORDER, r ) ;
1081  }
1082  saved_Br[i-ilow] = Interpolate( &Z3DSteel[jlow], save_Br, ORDER, z ) ;
1083  saved_Bz[i-ilow] = Interpolate( &Z3DSteel[jlow], save_Bz, ORDER, z ) ;
1084  saved_Bphi[i-ilow] = Interpolate( &Z3DSteel[jlow], save_Bphi, ORDER, z ) ;
1085  }
1086  Br_value = fscale * Interpolate( &Phi3DSteel[ilow], saved_Br, ORDER, phi ) ;
1087  Bz_value = fscale * Interpolate( &Phi3DSteel[ilow], saved_Bz, ORDER, phi ) ;
1088  Bphi_value = fscale * Interpolate( &Phi3DSteel[ilow], saved_Bphi, ORDER, phi ) ;
1089 
1090 }
1091 
1092 
1093 
1094 
1095 
1096 //end added by Lijuan
1097 
1098 //________________________________________
1099 #if 0
1100 
1102 void StarMagField::Interpolate2DEdistortion( const Float_t r, const Float_t z,
1103  const Float_t Er[neZ][neR], Float_t &Er_value )
1104 
1105 {
1106 
1107  const Int_t ORDER = 1 ; // Linear interpolation = 1, Quadratic = 2
1108  static Int_t jlow=0, klow=0 ;
1109  Float_t save_Er[ORDER+1] ;
1110 
1111  Search( neZ, eZList, z, jlow ) ;
1112  Search( neR, eRadius, r, klow ) ;
1113  if ( jlow < 0 ) jlow = 0 ; // artifact of Root's binsearch, returns -1 if out of range
1114  if ( klow < 0 ) klow = 0 ;
1115  if ( jlow + ORDER >= neZ - 1 ) jlow = neZ - 1 - ORDER ;
1116  if ( klow + ORDER >= neR - 1 ) klow = neR - 1 - ORDER ;
1117 
1118  for ( Int_t j = jlow ; j < jlow + ORDER + 1 ; j++ )
1119  {
1120  save_Er[j-jlow] = Interpolate( &eRadius[klow], &Er[j][klow], ORDER, r ) ;
1121  }
1122  Er_value = Interpolate( &eZList[jlow], save_Er, ORDER, z ) ;
1123 
1124 }
1125 
1127 
1128 void StarMagField::Interpolate3DEdistortion( const Float_t r, const Float_t phi, const Float_t z,
1129  const Float_t Er[neZ][nePhi][neR], const Float_t Ephi[neZ][nePhi][neR],
1130  Float_t &Er_value, Float_t &Ephi_value )
1131 
1132 {
1133 
1134  const Int_t ORDER = 1 ; // Linear interpolation = 1, Quadratic = 2
1135  static Int_t ilow=0, jlow=0, klow=0 ;
1136  Float_t save_Er[ORDER+1], saved_Er[ORDER+1] ;
1137  Float_t save_Ephi[ORDER+1], saved_Ephi[ORDER+1] ;
1138 
1139  Search( neZ, eZList, z, ilow ) ;
1140  Search( nePhi, ePhiList, phi, jlow ) ;
1141  Search( neR, eRadius, r, klow ) ;
1142  if ( ilow < 0 ) ilow = 0 ; // artifact of Root's binsearch, returns -1 if out of range
1143  if ( jlow < 0 ) jlow = 0 ;
1144  if ( klow < 0 ) klow = 0 ;
1145 
1146  if ( ilow + ORDER >= neZ - 1 ) ilow = neZ - 1 - ORDER ;
1147  if ( jlow + ORDER >= nePhi - 1 ) jlow = nePhi - 1 - ORDER ;
1148  if ( klow + ORDER >= neR - 1 ) klow = neR - 1 - ORDER ;
1149 
1150  for ( Int_t i = ilow ; i < ilow + ORDER + 1 ; i++ )
1151  {
1152  for ( Int_t j = jlow ; j < jlow + ORDER + 1 ; j++ )
1153  {
1154  save_Er[j-jlow] = Interpolate( &eRadius[klow], &Er[i][j][klow], ORDER, r ) ;
1155  save_Ephi[j-jlow] = Interpolate( &eRadius[klow], &Ephi[i][j][klow], ORDER, r ) ;
1156  }
1157  saved_Er[i-ilow] = Interpolate( &ePhiList[jlow], save_Er, ORDER, phi ) ;
1158  saved_Ephi[i-ilow] = Interpolate( &ePhiList[jlow], save_Ephi, ORDER, phi ) ;
1159  }
1160  Er_value = Interpolate( &eZList[ilow], saved_Er, ORDER, z ) ;
1161  Ephi_value = Interpolate( &eZList[ilow], saved_Ephi, ORDER, z ) ;
1162 
1163 }
1164 #endif
1165 
1166 //________________________________________
1167 
1169 
1170 Float_t StarMagField::Interpolate( const Float_t Xarray[], const Float_t Yarray[],
1171  const Int_t ORDER, const Float_t x )
1172 
1173 {
1174 
1175  Float_t y ;
1176 
1177 
1178  if ( ORDER == 2 ) // Quadratic Interpolation = 2
1179 
1180  {
1181  y = (x-Xarray[1]) * (x-Xarray[2]) * Yarray[0] / ( (Xarray[0]-Xarray[1]) * (Xarray[0]-Xarray[2]) ) ;
1182  y += (x-Xarray[2]) * (x-Xarray[0]) * Yarray[1] / ( (Xarray[1]-Xarray[2]) * (Xarray[1]-Xarray[0]) ) ;
1183  y += (x-Xarray[0]) * (x-Xarray[1]) * Yarray[2] / ( (Xarray[2]-Xarray[0]) * (Xarray[2]-Xarray[1]) ) ;
1184 
1185  }
1186 
1187  else // Linear Interpolation = 1
1188 
1189  {
1190  y = Yarray[0] + ( Yarray[1]-Yarray[0] ) * ( x-Xarray[0] ) / ( Xarray[1] - Xarray[0] ) ;
1191  }
1192 
1193  return (y) ;
1194 
1195 }
1196 
1197 
1198 //________________________________________
1199 
1201 
1202 void StarMagField::Search( Int_t N, const Float_t Xarray[], Float_t x, Int_t &low )
1203 
1204 {
1205  assert(! TMath::IsNaN(x));
1206  Long_t middle, high ;
1207  Int_t ascend = 0, increment = 1 ;
1208 
1209  if ( Xarray[N-1] >= Xarray[0] ) ascend = 1 ; // Ascending ordered table if true
1210 
1211  if ( low < 0 || low > N-1 ) { low = -1 ; high = N ; }
1212 
1213  else // Ordered Search phase
1214  {
1215  if ( (Int_t)( x >= Xarray[low] ) == ascend )
1216  {
1217  if ( low == N-1 ) return ;
1218  high = low + 1 ;
1219  while ( (Int_t)( x >= Xarray[high] ) == ascend )
1220  {
1221  low = high ;
1222  increment *= 2 ;
1223  high = low + increment ;
1224  if ( high > N-1 ) { high = N ; break ; }
1225  }
1226  }
1227  else
1228  {
1229  if ( low == 0 ) { low = -1 ; return ; }
1230  high = low - 1 ;
1231  while ( (Int_t)( x < Xarray[low] ) == ascend )
1232  {
1233  high = low ;
1234  increment *= 2 ;
1235  if ( increment >= high ) { low = -1 ; break ; }
1236  else low = high - increment ;
1237  }
1238  }
1239  }
1240 
1241  while ( (high-low) != 1 ) // Binary Search Phase
1242  {
1243  middle = ( high + low ) / 2 ;
1244  if ( (Int_t)( x >= Xarray[middle] ) == ascend )
1245  low = middle ;
1246  else
1247  high = middle ;
1248  }
1249 
1250  if ( x == Xarray[N-1] ) low = N-2 ;
1251  if ( x == Xarray[0] ) low = 0 ;
1252 
1253  }
1254 //________________________________________________________________________________
1255 #define PPLOCK(A) \
1256  void StarMagField::Set ## A (Float_t m) { \
1257  if (!fLock) f ## A = m; \
1258  else printf("StarMagField::Set"#A"() "#A" is locked at %f; Set to %f is ignored\n", f ## A ,m); \
1259  }
1260 PPLOCK(Factor)
1261 PPLOCK(Rescale)
1262 PPLOCK(BDipole)
1263 PPLOCK(RmaxDip)
1264 PPLOCK(ZminDip)
1265 PPLOCK(ZmaxDip)
1266 #undef PPLOCK
1267 //________________________________________________________________________________
1268 void StarMagField::SetLock () {
1269  if (! fLock) {
1270  fLock = kTRUE;
1271  printf("StarMagField::SetLock lock StarMagField parameters\n");
1272  Print();
1273  }
1274 }
1275 //________________________________________________________________________________
1276 #define PrintPar(A) printf("StarMagField:: "#A"\t%f\n",f ## A)
1277 void StarMagField::Print (Option_t*) const {
1278  if (fLock) printf("StarMagField parameters are locked\n");
1279  printf("StarMagField:: Map\t%i\n",fMap );
1280  PrintPar(Factor );
1281  PrintPar(Rescale);
1282  PrintPar(BDipole);
1283  PrintPar(RmaxDip);
1284  PrintPar(ZminDip);
1285  PrintPar(ZmaxDip);
1286 }
1287 #undef PrintPar
StarMagField::Interpolate3DBSteelfield
virtual void Interpolate3DBSteelfield(const Float_t r, const Float_t z, const Float_t phi, Float_t &Br_value, Float_t &Bz_value, Float_t &Bphi_value)
Interpolate the B field map - 3D interpolation.
Definition: StarMagField.cxx:1045
Error
Definition: FJcore.h:367
StarMagField::Interpolate2DBfield
virtual void Interpolate2DBfield(const Float_t r, const Float_t z, Float_t &Br_value, Float_t &Bz_value)
Interpolate the B field map - 2D interpolation.
Definition: StarMagField.cxx:913
StarMagField::B3DField
virtual void B3DField(const Float_t x[], Float_t B[])
Bfield in Cartesian coordinates - 3D field.
Definition: StarMagField.cxx:562
BDAT_t
Definition: StarMagField.cxx:242
BFLD_t
Definition: StarMagField.cxx:219
StarMagField::Search
static void Search(Int_t N, const Float_t Xarray[], Float_t x, Int_t &low)
Search an ordered table by starting at the most recently used point.
Definition: StarMagField.cxx:1202
StarMagField::Interpolate
virtual Float_t Interpolate(const Float_t Xarray[], const Float_t Yarray[], const Int_t ORDER, const Float_t x)
Interpolate a 3x2 table (quadratic) or a 2x2 table (linear)
Definition: StarMagField.cxx:1170
StarMagField::Interpolate3DBfield
virtual void Interpolate3DBfield(const Float_t r, const Float_t z, const Float_t phi, Float_t &Br_value, Float_t &Bz_value, Float_t &Bphi_value)
Interpolate the B field map - 3D interpolation.
Definition: StarMagField.cxx:985
StarMagField::BrBzField
virtual void BrBzField(const Float_t r, const Float_t z, Float_t &Br_value, Float_t &Bz_value)
B field in Radial coordinates - 2D field (ie Phi symmetric)
Definition: StarMagField.cxx:619
Generated by   doxygen 1.8.5