Project
Loading...
Searching...
No Matches
MagneticWrapperChebyshev.cxx
Go to the documentation of this file.
1// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3// All rights not expressly granted are reserved.
4//
5// This software is distributed under the terms of the GNU General Public
6// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7//
8// In applying this license CERN does not waive the privileges and immunities
9// granted to it by virtue of its status as an Intergovernmental Organization
10// or submit itself to any jurisdiction.
11
15
16#include "Field/MagneticWrapperChebyshev.h"
17#include <TArrayF.h> // for TArrayF
18#include <TArrayI.h> // for TArrayI
19#include <TSystem.h> // for TSystem, gSystem
20#include <cstdio> // for printf, fprintf, fclose, fopen, FILE
21#include <cstring> // for memcpy
22#include <fairlogger/Logger.h> // for FairLogger
23#include "TMath.h" // for BinarySearch, Sort
24#include "TMathBase.h" // for Abs
25#include "TNamed.h" // for TNamed
26#include "TObjArray.h" // for TObjArray
27#include "TString.h" // for TString
28
29using namespace o2::field;
30using namespace o2::math_utils;
31
32ClassImp(MagneticWrapperChebyshev);
33
34MagneticWrapperChebyshev::MagneticWrapperChebyshev()
35 : mNumberOfParameterizationSolenoid(0),
36 mNumberOfDistinctZSegmentsSolenoid(0),
37 mNumberOfDistinctPSegmentsSolenoid(0),
38 mNumberOfDistinctRSegmentsSolenoid(0),
39 mCoordinatesSegmentsZSolenoid(nullptr),
40 mCoordinatesSegmentsPSolenoid(nullptr),
41 mCoordinatesSegmentsRSolenoid(nullptr),
42 mBeginningOfSegmentsPSolenoid(nullptr),
43 mNumberOfSegmentsPSolenoid(nullptr),
44 mBeginningOfSegmentsRSolenoid(nullptr),
45 mNumberOfRSegmentsSolenoid(nullptr),
46 mSegmentIdSolenoid(nullptr),
47 mMinZSolenoid(1.e6),
48 mMaxZSolenoid(-1.e6),
49 mParameterizationSolenoid(nullptr),
50 mMaxRadiusSolenoid(0),
51 mNumberOfParameterizationTPC(0),
52 mNumberOfDistinctZSegmentsTPC(0),
53 mNumberOfDistinctPSegmentsTPC(0),
54 mNumberOfDistinctRSegmentsTPC(0),
55 mCoordinatesSegmentsZTPC(nullptr),
56 mCoordinatesSegmentsPTPC(nullptr),
57 mCoordinatesSegmentsRTPC(nullptr),
58 mBeginningOfSegmentsPTPC(nullptr),
59 mNumberOfSegmentsPTPC(nullptr),
60 mBeginningOfSegmentsRTPC(nullptr),
61 mNumberOfRSegmentsTPC(nullptr),
62 mSegmentIdTPC(nullptr),
63 mMinZTPC(1.e6),
64 mMaxZTPC(-1.e6),
65 mParameterizationTPC(nullptr),
66 mMaxRadiusTPC(0),
67 mNumberOfParameterizationTPCRat(0),
68 mNumberOfDistinctZSegmentsTPCRat(0),
69 mNumberOfDistinctPSegmentsTPCRat(0),
70 mNumberOfDistinctRSegmentsTPCRat(0),
71 mCoordinatesSegmentsZTPCRat(nullptr),
72 mCoordinatesSegmentsPTPCRat(nullptr),
73 mCoordinatesSegmentsRTPCRat(nullptr),
74 mBeginningOfSegmentsPTPCRat(nullptr),
75 mNumberOfSegmentsPTPCRat(nullptr),
76 mBeginningOfSegmentsRTPCRat(nullptr),
77 mNumberOfRSegmentsTPCRat(nullptr),
78 mSegmentIdTPCRat(nullptr),
79 mMinZTPCRat(1.e6),
80 mMaxZTPCRat(-1.e6),
81 mParameterizationTPCRat(nullptr),
82 mMaxRadiusTPCRat(0),
83 mNumberOfParameterizationDipole(0),
84 mNumberOfDistinctZSegmentsDipole(0),
85 mNumberOfDistinctYSegmentsDipole(0),
86 mNumberOfDistinctXSegmentsDipole(0),
87 mCoordinatesSegmentsZDipole(nullptr),
88 mCoordinatesSegmentsYDipole(nullptr),
89 mCoordinatesSegmentsXDipole(nullptr),
90 mBeginningOfSegmentsYDipole(nullptr),
91 mNumberOfSegmentsYDipole(nullptr),
92 mBeginningOfSegmentsXDipole(nullptr),
93 mNumberOfSegmentsXDipole(nullptr),
94 mSegmentIdDipole(nullptr),
95 mMinDipoleZ(1.e6),
96 mMaxDipoleZ(-1.e6),
97 mParameterizationDipole(nullptr)
98{
99}
100
101MagneticWrapperChebyshev::MagneticWrapperChebyshev(const MagneticWrapperChebyshev& src)
102 : TNamed(src),
103 mNumberOfParameterizationSolenoid(0),
104 mNumberOfDistinctZSegmentsSolenoid(0),
105 mNumberOfDistinctPSegmentsSolenoid(0),
106 mNumberOfDistinctRSegmentsSolenoid(0),
107 mCoordinatesSegmentsZSolenoid(nullptr),
108 mCoordinatesSegmentsPSolenoid(nullptr),
109 mCoordinatesSegmentsRSolenoid(nullptr),
110 mBeginningOfSegmentsPSolenoid(nullptr),
111 mNumberOfSegmentsPSolenoid(nullptr),
112 mBeginningOfSegmentsRSolenoid(nullptr),
113 mNumberOfRSegmentsSolenoid(nullptr),
114 mSegmentIdSolenoid(nullptr),
115 mMinZSolenoid(1.e6),
116 mMaxZSolenoid(-1.e6),
117 mParameterizationSolenoid(nullptr),
118 mMaxRadiusSolenoid(0),
119 mNumberOfParameterizationTPC(0),
120 mNumberOfDistinctZSegmentsTPC(0),
121 mNumberOfDistinctPSegmentsTPC(0),
122 mNumberOfDistinctRSegmentsTPC(0),
123 mCoordinatesSegmentsZTPC(nullptr),
124 mCoordinatesSegmentsPTPC(nullptr),
125 mCoordinatesSegmentsRTPC(nullptr),
126 mBeginningOfSegmentsPTPC(nullptr),
127 mNumberOfSegmentsPTPC(nullptr),
128 mBeginningOfSegmentsRTPC(nullptr),
129 mNumberOfRSegmentsTPC(nullptr),
130 mSegmentIdTPC(nullptr),
131 mMinZTPC(1.e6),
132 mMaxZTPC(-1.e6),
133 mParameterizationTPC(nullptr),
134 mMaxRadiusTPC(0),
135 mNumberOfParameterizationTPCRat(0),
136 mNumberOfDistinctZSegmentsTPCRat(0),
137 mNumberOfDistinctPSegmentsTPCRat(0),
138 mNumberOfDistinctRSegmentsTPCRat(0),
139 mCoordinatesSegmentsZTPCRat(nullptr),
140 mCoordinatesSegmentsPTPCRat(nullptr),
141 mCoordinatesSegmentsRTPCRat(nullptr),
142 mBeginningOfSegmentsPTPCRat(nullptr),
143 mNumberOfSegmentsPTPCRat(nullptr),
144 mBeginningOfSegmentsRTPCRat(nullptr),
145 mNumberOfRSegmentsTPCRat(nullptr),
146 mSegmentIdTPCRat(nullptr),
147 mMinZTPCRat(1.e6),
148 mMaxZTPCRat(-1.e6),
149 mParameterizationTPCRat(nullptr),
150 mMaxRadiusTPCRat(0),
151 mNumberOfParameterizationDipole(0),
152 mNumberOfDistinctZSegmentsDipole(0),
153 mNumberOfDistinctYSegmentsDipole(0),
154 mNumberOfDistinctXSegmentsDipole(0),
155 mCoordinatesSegmentsZDipole(nullptr),
156 mCoordinatesSegmentsYDipole(nullptr),
157 mCoordinatesSegmentsXDipole(nullptr),
158 mBeginningOfSegmentsYDipole(nullptr),
159 mNumberOfSegmentsYDipole(nullptr),
160 mBeginningOfSegmentsXDipole(nullptr),
161 mNumberOfSegmentsXDipole(nullptr),
162 mSegmentIdDipole(nullptr),
163 mMinDipoleZ(1.e6),
164 mMaxDipoleZ(-1.e6),
165 mParameterizationDipole(nullptr)
166{
167 copyFrom(src);
168}
169
170void MagneticWrapperChebyshev::copyFrom(const MagneticWrapperChebyshev& src)
171{
172 Clear();
173 SetName(src.GetName());
174 SetTitle(src.GetTitle());
175
176 mNumberOfParameterizationSolenoid = src.mNumberOfParameterizationSolenoid;
177 mNumberOfDistinctZSegmentsSolenoid = src.mNumberOfDistinctZSegmentsSolenoid;
178 mNumberOfDistinctPSegmentsSolenoid = src.mNumberOfDistinctPSegmentsSolenoid;
179 mNumberOfDistinctRSegmentsSolenoid = src.mNumberOfDistinctRSegmentsSolenoid;
180 mMinZSolenoid = src.mMinZSolenoid;
181 mMaxZSolenoid = src.mMaxZSolenoid;
182 mMaxRadiusSolenoid = src.mMaxRadiusSolenoid;
183 if (src.mNumberOfParameterizationSolenoid) {
184 memcpy(mCoordinatesSegmentsZSolenoid = new Float_t[mNumberOfDistinctZSegmentsSolenoid],
185 src.mCoordinatesSegmentsZSolenoid, sizeof(Float_t) * mNumberOfDistinctZSegmentsSolenoid);
186 memcpy(mCoordinatesSegmentsPSolenoid = new Float_t[mNumberOfDistinctPSegmentsSolenoid],
187 src.mCoordinatesSegmentsPSolenoid, sizeof(Float_t) * mNumberOfDistinctPSegmentsSolenoid);
188 memcpy(mCoordinatesSegmentsRSolenoid = new Float_t[mNumberOfDistinctRSegmentsSolenoid],
189 src.mCoordinatesSegmentsRSolenoid, sizeof(Float_t) * mNumberOfDistinctRSegmentsSolenoid);
190 memcpy(mBeginningOfSegmentsPSolenoid = new Int_t[mNumberOfDistinctZSegmentsSolenoid],
191 src.mBeginningOfSegmentsPSolenoid, sizeof(Int_t) * mNumberOfDistinctZSegmentsSolenoid);
192 memcpy(mNumberOfSegmentsPSolenoid = new Int_t[mNumberOfDistinctZSegmentsSolenoid], src.mNumberOfSegmentsPSolenoid,
193 sizeof(Int_t) * mNumberOfDistinctZSegmentsSolenoid);
194 memcpy(mBeginningOfSegmentsRSolenoid = new Int_t[mNumberOfDistinctPSegmentsSolenoid],
195 src.mBeginningOfSegmentsRSolenoid, sizeof(Int_t) * mNumberOfDistinctPSegmentsSolenoid);
196 memcpy(mNumberOfRSegmentsSolenoid = new Int_t[mNumberOfDistinctPSegmentsSolenoid], src.mNumberOfRSegmentsSolenoid,
197 sizeof(Int_t) * mNumberOfDistinctPSegmentsSolenoid);
198 memcpy(mSegmentIdSolenoid = new Int_t[mNumberOfDistinctRSegmentsSolenoid], src.mSegmentIdSolenoid,
199 sizeof(Int_t) * mNumberOfDistinctRSegmentsSolenoid);
200 mParameterizationSolenoid = new TObjArray(mNumberOfParameterizationSolenoid);
201 for (int i = 0; i < mNumberOfParameterizationSolenoid; i++) {
202 mParameterizationSolenoid->AddAtAndExpand(new Chebyshev3D(*src.getParameterSolenoid(i)), i);
203 }
204 }
205
206 mNumberOfParameterizationTPC = src.mNumberOfParameterizationTPC;
207 mNumberOfDistinctZSegmentsTPC = src.mNumberOfDistinctZSegmentsTPC;
208 mNumberOfDistinctPSegmentsTPC = src.mNumberOfDistinctPSegmentsTPC;
209 mNumberOfDistinctRSegmentsTPC = src.mNumberOfDistinctRSegmentsTPC;
210 mMinZTPC = src.mMinZTPC;
211 mMaxZTPC = src.mMaxZTPC;
212 mMaxRadiusTPC = src.mMaxRadiusTPC;
213 if (src.mNumberOfParameterizationTPC) {
214 memcpy(mCoordinatesSegmentsZTPC = new Float_t[mNumberOfDistinctZSegmentsTPC], src.mCoordinatesSegmentsZTPC,
215 sizeof(Float_t) * mNumberOfDistinctZSegmentsTPC);
216 memcpy(mCoordinatesSegmentsPTPC = new Float_t[mNumberOfDistinctPSegmentsTPC], src.mCoordinatesSegmentsPTPC,
217 sizeof(Float_t) * mNumberOfDistinctPSegmentsTPC);
218 memcpy(mCoordinatesSegmentsRTPC = new Float_t[mNumberOfDistinctRSegmentsTPC], src.mCoordinatesSegmentsRTPC,
219 sizeof(Float_t) * mNumberOfDistinctRSegmentsTPC);
220 memcpy(mBeginningOfSegmentsPTPC = new Int_t[mNumberOfDistinctZSegmentsTPC], src.mBeginningOfSegmentsPTPC,
221 sizeof(Int_t) * mNumberOfDistinctZSegmentsTPC);
222 memcpy(mNumberOfSegmentsPTPC = new Int_t[mNumberOfDistinctZSegmentsTPC], src.mNumberOfSegmentsPTPC,
223 sizeof(Int_t) * mNumberOfDistinctZSegmentsTPC);
224 memcpy(mBeginningOfSegmentsRTPC = new Int_t[mNumberOfDistinctPSegmentsTPC], src.mBeginningOfSegmentsRTPC,
225 sizeof(Int_t) * mNumberOfDistinctPSegmentsTPC);
226 memcpy(mNumberOfRSegmentsTPC = new Int_t[mNumberOfDistinctPSegmentsTPC], src.mNumberOfRSegmentsTPC,
227 sizeof(Int_t) * mNumberOfDistinctPSegmentsTPC);
228 memcpy(mSegmentIdTPC = new Int_t[mNumberOfDistinctRSegmentsTPC], src.mSegmentIdTPC,
229 sizeof(Int_t) * mNumberOfDistinctRSegmentsTPC);
230 mParameterizationTPC = new TObjArray(mNumberOfParameterizationTPC);
231 for (int i = 0; i < mNumberOfParameterizationTPC; i++) {
232 mParameterizationTPC->AddAtAndExpand(new Chebyshev3D(*src.getParameterTPCIntegral(i)), i);
233 }
234 }
235
236 mNumberOfParameterizationTPCRat = src.mNumberOfParameterizationTPCRat;
237 mNumberOfDistinctZSegmentsTPCRat = src.mNumberOfDistinctZSegmentsTPCRat;
238 mNumberOfDistinctPSegmentsTPCRat = src.mNumberOfDistinctPSegmentsTPCRat;
239 mNumberOfDistinctRSegmentsTPCRat = src.mNumberOfDistinctRSegmentsTPCRat;
240 mMinZTPCRat = src.mMinZTPCRat;
241 mMaxZTPCRat = src.mMaxZTPCRat;
242 mMaxRadiusTPCRat = src.mMaxRadiusTPCRat;
243 if (src.mNumberOfParameterizationTPCRat) {
244 memcpy(mCoordinatesSegmentsZTPCRat = new Float_t[mNumberOfDistinctZSegmentsTPCRat], src.mCoordinatesSegmentsZTPCRat,
245 sizeof(Float_t) * mNumberOfDistinctZSegmentsTPCRat);
246 memcpy(mCoordinatesSegmentsPTPCRat = new Float_t[mNumberOfDistinctPSegmentsTPCRat], src.mCoordinatesSegmentsPTPCRat,
247 sizeof(Float_t) * mNumberOfDistinctPSegmentsTPCRat);
248 memcpy(mCoordinatesSegmentsRTPCRat = new Float_t[mNumberOfDistinctRSegmentsTPCRat], src.mCoordinatesSegmentsRTPCRat,
249 sizeof(Float_t) * mNumberOfDistinctRSegmentsTPCRat);
250 memcpy(mBeginningOfSegmentsPTPCRat = new Int_t[mNumberOfDistinctZSegmentsTPCRat], src.mBeginningOfSegmentsPTPCRat,
251 sizeof(Int_t) * mNumberOfDistinctZSegmentsTPCRat);
252 memcpy(mNumberOfSegmentsPTPCRat = new Int_t[mNumberOfDistinctZSegmentsTPCRat], src.mNumberOfSegmentsPTPCRat,
253 sizeof(Int_t) * mNumberOfDistinctZSegmentsTPCRat);
254 memcpy(mBeginningOfSegmentsRTPCRat = new Int_t[mNumberOfDistinctPSegmentsTPCRat], src.mBeginningOfSegmentsRTPCRat,
255 sizeof(Int_t) * mNumberOfDistinctPSegmentsTPCRat);
256 memcpy(mNumberOfRSegmentsTPCRat = new Int_t[mNumberOfDistinctPSegmentsTPCRat], src.mNumberOfRSegmentsTPCRat,
257 sizeof(Int_t) * mNumberOfDistinctPSegmentsTPCRat);
258 memcpy(mSegmentIdTPCRat = new Int_t[mNumberOfDistinctRSegmentsTPCRat], src.mSegmentIdTPCRat,
259 sizeof(Int_t) * mNumberOfDistinctRSegmentsTPCRat);
260 mParameterizationTPCRat = new TObjArray(mNumberOfParameterizationTPCRat);
261 for (int i = 0; i < mNumberOfParameterizationTPCRat; i++) {
262 mParameterizationTPCRat->AddAtAndExpand(new Chebyshev3D(*src.getParameterTPCRatIntegral(i)), i);
263 }
264 }
265
266 mNumberOfParameterizationDipole = src.mNumberOfParameterizationDipole;
267 mNumberOfDistinctZSegmentsDipole = src.mNumberOfDistinctZSegmentsDipole;
268 mNumberOfDistinctYSegmentsDipole = src.mNumberOfDistinctYSegmentsDipole;
269 mNumberOfDistinctXSegmentsDipole = src.mNumberOfDistinctXSegmentsDipole;
270 mMinDipoleZ = src.mMinDipoleZ;
271 mMaxDipoleZ = src.mMaxDipoleZ;
272 if (src.mNumberOfParameterizationDipole) {
273 memcpy(mCoordinatesSegmentsZDipole = new Float_t[mNumberOfDistinctZSegmentsDipole], src.mCoordinatesSegmentsZDipole,
274 sizeof(Float_t) * mNumberOfDistinctZSegmentsDipole);
275 memcpy(mCoordinatesSegmentsYDipole = new Float_t[mNumberOfDistinctYSegmentsDipole], src.mCoordinatesSegmentsYDipole,
276 sizeof(Float_t) * mNumberOfDistinctYSegmentsDipole);
277 memcpy(mCoordinatesSegmentsXDipole = new Float_t[mNumberOfDistinctXSegmentsDipole], src.mCoordinatesSegmentsXDipole,
278 sizeof(Float_t) * mNumberOfDistinctXSegmentsDipole);
279 memcpy(mBeginningOfSegmentsYDipole = new Int_t[mNumberOfDistinctZSegmentsDipole], src.mBeginningOfSegmentsYDipole,
280 sizeof(Int_t) * mNumberOfDistinctZSegmentsDipole);
281 memcpy(mNumberOfSegmentsYDipole = new Int_t[mNumberOfDistinctZSegmentsDipole], src.mNumberOfSegmentsYDipole,
282 sizeof(Int_t) * mNumberOfDistinctZSegmentsDipole);
283 memcpy(mBeginningOfSegmentsXDipole = new Int_t[mNumberOfDistinctYSegmentsDipole], src.mBeginningOfSegmentsXDipole,
284 sizeof(Int_t) * mNumberOfDistinctYSegmentsDipole);
285 memcpy(mNumberOfSegmentsXDipole = new Int_t[mNumberOfDistinctYSegmentsDipole], src.mNumberOfSegmentsXDipole,
286 sizeof(Int_t) * mNumberOfDistinctYSegmentsDipole);
287 memcpy(mSegmentIdDipole = new Int_t[mNumberOfDistinctXSegmentsDipole], src.mSegmentIdDipole,
288 sizeof(Int_t) * mNumberOfDistinctXSegmentsDipole);
289 mParameterizationDipole = new TObjArray(mNumberOfParameterizationDipole);
290 for (int i = 0; i < mNumberOfParameterizationDipole; i++) {
291 mParameterizationDipole->AddAtAndExpand(new Chebyshev3D(*src.getParameterDipole(i)), i);
292 }
293 }
294}
295
296MagneticWrapperChebyshev& MagneticWrapperChebyshev::operator=(const MagneticWrapperChebyshev& rhs)
297{
298 if (this != &rhs) {
299 Clear();
300 copyFrom(rhs);
301 }
302 return *this;
303}
304
305void MagneticWrapperChebyshev::Clear(const Option_t*)
306{
307 if (mNumberOfParameterizationSolenoid) {
308 mParameterizationSolenoid->SetOwner(kTRUE);
309 delete mParameterizationSolenoid;
310 mParameterizationSolenoid = nullptr;
311 delete[] mCoordinatesSegmentsZSolenoid;
312 mCoordinatesSegmentsZSolenoid = nullptr;
313 delete[] mCoordinatesSegmentsPSolenoid;
314 mCoordinatesSegmentsPSolenoid = nullptr;
315 delete[] mCoordinatesSegmentsRSolenoid;
316 mCoordinatesSegmentsRSolenoid = nullptr;
317 delete[] mBeginningOfSegmentsPSolenoid;
318 mBeginningOfSegmentsPSolenoid = nullptr;
319 delete[] mNumberOfSegmentsPSolenoid;
320 mNumberOfSegmentsPSolenoid = nullptr;
321 delete[] mBeginningOfSegmentsRSolenoid;
322 mBeginningOfSegmentsRSolenoid = nullptr;
323 delete[] mNumberOfRSegmentsSolenoid;
324 mNumberOfRSegmentsSolenoid = nullptr;
325 delete[] mSegmentIdSolenoid;
326 mSegmentIdSolenoid = nullptr;
327 }
328
329 mNumberOfParameterizationSolenoid = mNumberOfDistinctZSegmentsSolenoid = mNumberOfDistinctPSegmentsSolenoid =
330 mNumberOfDistinctRSegmentsSolenoid = 0;
331 mMinZSolenoid = 1e6;
332 mMaxZSolenoid = -1e6;
333 mMaxRadiusSolenoid = 0;
334
335 if (mNumberOfParameterizationTPC) {
336 mParameterizationTPC->SetOwner(kTRUE);
337 delete mParameterizationTPC;
338 mParameterizationTPC = nullptr;
339 delete[] mCoordinatesSegmentsZTPC;
340 mCoordinatesSegmentsZTPC = nullptr;
341 delete[] mCoordinatesSegmentsPTPC;
342 mCoordinatesSegmentsPTPC = nullptr;
343 delete[] mCoordinatesSegmentsRTPC;
344 mCoordinatesSegmentsRTPC = nullptr;
345 delete[] mBeginningOfSegmentsPTPC;
346 mBeginningOfSegmentsPTPC = nullptr;
347 delete[] mNumberOfSegmentsPTPC;
348 mNumberOfSegmentsPTPC = nullptr;
349 delete[] mBeginningOfSegmentsRTPC;
350 mBeginningOfSegmentsRTPC = nullptr;
351 delete[] mNumberOfRSegmentsTPC;
352 mNumberOfRSegmentsTPC = nullptr;
353 delete[] mSegmentIdTPC;
354 mSegmentIdTPC = nullptr;
355 }
356
357 mNumberOfParameterizationTPC = mNumberOfDistinctZSegmentsTPC = mNumberOfDistinctPSegmentsTPC =
358 mNumberOfDistinctRSegmentsTPC = 0;
359 mMinZTPC = 1e6;
360 mMaxZTPC = -1e6;
361 mMaxRadiusTPC = 0;
362
363 if (mNumberOfParameterizationTPCRat) {
364 mParameterizationTPCRat->SetOwner(kTRUE);
365 delete mParameterizationTPCRat;
366 mParameterizationTPCRat = nullptr;
367 delete[] mCoordinatesSegmentsZTPCRat;
368 mCoordinatesSegmentsZTPCRat = nullptr;
369 delete[] mCoordinatesSegmentsPTPCRat;
370 mCoordinatesSegmentsPTPCRat = nullptr;
371 delete[] mCoordinatesSegmentsRTPCRat;
372 mCoordinatesSegmentsRTPCRat = nullptr;
373 delete[] mBeginningOfSegmentsPTPCRat;
374 mBeginningOfSegmentsPTPCRat = nullptr;
375 delete[] mNumberOfSegmentsPTPCRat;
376 mNumberOfSegmentsPTPCRat = nullptr;
377 delete[] mBeginningOfSegmentsRTPCRat;
378 mBeginningOfSegmentsRTPCRat = nullptr;
379 delete[] mNumberOfRSegmentsTPCRat;
380 mNumberOfRSegmentsTPCRat = nullptr;
381 delete[] mSegmentIdTPCRat;
382 mSegmentIdTPCRat = nullptr;
383 }
384
385 mNumberOfParameterizationTPCRat = mNumberOfDistinctZSegmentsTPCRat = mNumberOfDistinctPSegmentsTPCRat =
386 mNumberOfDistinctRSegmentsTPCRat = 0;
387 mMinZTPCRat = 1e6;
388 mMaxZTPCRat = -1e6;
389 mMaxRadiusTPCRat = 0;
390
391 if (mNumberOfParameterizationDipole) {
392 mParameterizationDipole->SetOwner(kTRUE);
393 delete mParameterizationDipole;
394 mParameterizationDipole = nullptr;
395 delete[] mCoordinatesSegmentsZDipole;
396 mCoordinatesSegmentsZDipole = nullptr;
397 delete[] mCoordinatesSegmentsYDipole;
398 mCoordinatesSegmentsYDipole = nullptr;
399 delete[] mCoordinatesSegmentsXDipole;
400 mCoordinatesSegmentsXDipole = nullptr;
401 delete[] mBeginningOfSegmentsYDipole;
402 mBeginningOfSegmentsYDipole = nullptr;
403 delete[] mNumberOfSegmentsYDipole;
404 mNumberOfSegmentsYDipole = nullptr;
405 delete[] mBeginningOfSegmentsXDipole;
406 mBeginningOfSegmentsXDipole = nullptr;
407 delete[] mNumberOfSegmentsXDipole;
408 mNumberOfSegmentsXDipole = nullptr;
409 delete[] mSegmentIdDipole;
410 mSegmentIdDipole = nullptr;
411 }
412
413 mNumberOfParameterizationDipole = mNumberOfDistinctZSegmentsDipole = mNumberOfDistinctYSegmentsDipole =
414 mNumberOfDistinctXSegmentsDipole = 0;
415 mMinDipoleZ = 1e6;
416 mMaxDipoleZ = -1e6;
417}
418
419void MagneticWrapperChebyshev::Field(const Double_t* xyz, Double_t* b) const
420{
421 Double_t rphiz[3];
422
423#ifndef _BRING_TO_BOUNDARY_ // exact matching to fitted volume is requested
424 b[0] = b[1] = b[2] = 0;
425#endif
426
427 if (xyz[2] > mMinZSolenoid) {
428 cartesianToCylindrical(xyz, rphiz);
429 fieldCylindricalSolenoid(rphiz, b);
430 // convert field to cartesian system
431 cylindricalToCartesianCylB(rphiz, b, b);
432 return;
433 }
434
435 int iddip = findDipoleSegment(xyz);
436 if (iddip < 0) {
437 return;
438 }
439 Chebyshev3D* par = getParameterDipole(iddip);
440#ifndef _BRING_TO_BOUNDARY_
441 if (!par->isInside(xyz)) {
442 return;
443 }
444#endif
445 par->Eval(xyz, b);
446}
447
448Double_t MagneticWrapperChebyshev::getBz(const Double_t* xyz) const
449{
450 Double_t rphiz[3];
451
452 if (xyz[2] > mMinZSolenoid) {
453 cartesianToCylindrical(xyz, rphiz);
454 return fieldCylindricalSolenoidBz(rphiz);
455 }
456
457 int iddip = findDipoleSegment(xyz);
458 if (iddip < 0) {
459 return 0.;
460 }
461 Chebyshev3D* par = getParameterDipole(iddip);
462#ifndef _BRING_TO_BOUNDARY_
463 if (!par->isInside(xyz)) {
464 return 0.;
465 }
466#endif
467 return par->Eval(xyz, 2);
468}
469
470void MagneticWrapperChebyshev::Print(Option_t*) const
471{
472 printf("Alice magnetic field parameterized by Chebyshev polynomials\n");
473 printf("Segmentation for Solenoid (%+.2f<Z<%+.2f cm | R<%.2f cm)\n", mMinZSolenoid, mMaxZSolenoid,
474 mMaxRadiusSolenoid);
475
476 if (mParameterizationSolenoid) {
477 for (int i = 0; i < mNumberOfParameterizationSolenoid; i++) {
478 printf("SOL%4d ", i);
479 getParameterSolenoid(i)->Print();
480 }
481 }
482
483 printf("Segmentation for TPC field integral (%+.2f<Z<%+.2f cm | R<%.2f cm)\n", mMinZTPC, mMaxZTPC, mMaxRadiusTPC);
484
485 if (mParameterizationTPC) {
486 for (int i = 0; i < mNumberOfParameterizationTPC; i++) {
487 printf("TPC%4d ", i);
488 getParameterTPCIntegral(i)->Print();
489 }
490 }
491
492 printf("Segmentation for TPC field ratios integral (%+.2f<Z<%+.2f cm | R<%.2f cm)\n", mMinZTPCRat, mMaxZTPCRat,
493 mMaxRadiusTPCRat);
494
495 if (mParameterizationTPCRat) {
496 for (int i = 0; i < mNumberOfParameterizationTPCRat; i++) {
497 printf("TPCRat%4d ", i);
498 getParameterTPCRatIntegral(i)->Print();
499 }
500 }
501
502 printf("Segmentation for Dipole (%+.2f<Z<%+.2f cm)\n", mMinDipoleZ, mMaxDipoleZ);
503 if (mParameterizationDipole) {
504 for (int i = 0; i < mNumberOfParameterizationDipole; i++) {
505 printf("DIP%4d ", i);
506 getParameterDipole(i)->Print();
507 }
508 }
509}
510
511Int_t MagneticWrapperChebyshev::findDipoleSegment(const Double_t* xyz) const
512{
513 if (!mNumberOfParameterizationDipole) {
514 return -1;
515 }
516 int xid, yid, zid = TMath::BinarySearch(mNumberOfDistinctZSegmentsDipole, mCoordinatesSegmentsZDipole,
517 (Float_t)xyz[2]); // find zsegment
518
519 Bool_t reCheck = kFALSE;
520 while (true) {
521 int ysegBeg = mBeginningOfSegmentsYDipole[zid];
522
523 for (yid = 0; yid < mNumberOfSegmentsYDipole[zid]; yid++) {
524 if (xyz[1] < mCoordinatesSegmentsYDipole[ysegBeg + yid]) {
525 break;
526 }
527 }
528 if (--yid < 0) {
529 yid = 0;
530 }
531 yid += ysegBeg;
532
533 int xsegBeg = mBeginningOfSegmentsXDipole[yid];
534 for (xid = 0; xid < mNumberOfSegmentsXDipole[yid]; xid++) {
535 if (xyz[0] < mCoordinatesSegmentsXDipole[xsegBeg + xid]) {
536 break;
537 }
538 }
539
540 if (--xid < 0) {
541 xid = 0;
542 }
543 xid += xsegBeg;
544
545 // to make sure that due to the precision problems we did not pick the next Zbin
546 if (!reCheck && (xyz[2] - mCoordinatesSegmentsZDipole[zid] < 3.e-5) && zid &&
547 !getParameterDipole(mSegmentIdDipole[xid])->isInside(xyz)) { // check the previous Z bin
548 zid--;
549 reCheck = kTRUE;
550 continue;
551 }
552 break;
553 }
554 return mSegmentIdDipole[xid];
555}
556
557Int_t MagneticWrapperChebyshev::findSolenoidSegment(const Double_t* rpz) const
558{
559 if (!mNumberOfParameterizationSolenoid) {
560 return -1;
561 }
562 int rid, pid, zid = TMath::BinarySearch(mNumberOfDistinctZSegmentsSolenoid, mCoordinatesSegmentsZSolenoid,
563 (Float_t)rpz[2]); // find zsegment
564
565 Bool_t reCheck = kFALSE;
566 while (true) {
567 int psegBeg = mBeginningOfSegmentsPSolenoid[zid];
568 for (pid = 0; pid < mNumberOfSegmentsPSolenoid[zid]; pid++) {
569 if (rpz[1] < mCoordinatesSegmentsPSolenoid[psegBeg + pid]) {
570 break;
571 }
572 }
573 if (--pid < 0) {
574 pid = 0;
575 }
576 pid += psegBeg;
577
578 int rsegBeg = mBeginningOfSegmentsRSolenoid[pid];
579 for (rid = 0; rid < mNumberOfRSegmentsSolenoid[pid]; rid++) {
580 if (rpz[0] < mCoordinatesSegmentsRSolenoid[rsegBeg + rid]) {
581 break;
582 }
583 }
584 if (--rid < 0) {
585 rid = 0;
586 }
587 rid += rsegBeg;
588
589 // to make sure that due to the precision problems we did not pick the next Zbin
590 if (!reCheck && (rpz[2] - mCoordinatesSegmentsZSolenoid[zid] < 3.e-5) && zid &&
591 !getParameterSolenoid(mSegmentIdSolenoid[rid])->isInside(rpz)) { // check the previous Z bin
592 zid--;
593 reCheck = kTRUE;
594 continue;
595 }
596 break;
597 }
598 return mSegmentIdSolenoid[rid];
599}
600
601Int_t MagneticWrapperChebyshev::findTPCSegment(const Double_t* rpz) const
602{
603 if (!mNumberOfParameterizationTPC) {
604 return -1;
605 }
606 int rid, pid, zid = TMath::BinarySearch(mNumberOfDistinctZSegmentsTPC, mCoordinatesSegmentsZTPC,
607 (Float_t)rpz[2]); // find zsegment
608
609 Bool_t reCheck = kFALSE;
610 while (true) {
611 int psegBeg = mBeginningOfSegmentsPTPC[zid];
612
613 for (pid = 0; pid < mNumberOfSegmentsPTPC[zid]; pid++) {
614 if (rpz[1] < mCoordinatesSegmentsPTPC[psegBeg + pid]) {
615 break;
616 }
617 }
618 if (--pid < 0) {
619 pid = 0;
620 }
621 pid += psegBeg;
622
623 int rsegBeg = mBeginningOfSegmentsRTPC[pid];
624 for (rid = 0; rid < mNumberOfRSegmentsTPC[pid]; rid++) {
625 if (rpz[0] < mCoordinatesSegmentsRTPC[rsegBeg + rid]) {
626 break;
627 }
628 }
629 if (--rid < 0) {
630 rid = 0;
631 }
632 rid += rsegBeg;
633
634 // to make sure that due to the precision problems we did not pick the next Zbin
635 if (!reCheck && (rpz[2] - mCoordinatesSegmentsZTPC[zid] < 3.e-5) && zid &&
636 !getParameterTPCIntegral(mSegmentIdTPC[rid])->isInside(rpz)) { // check the previous Z bin
637 zid--;
638 reCheck = kTRUE;
639 continue;
640 }
641 break;
642 }
643 return mSegmentIdTPC[rid];
644}
645
646Int_t MagneticWrapperChebyshev::findTPCRatSegment(const Double_t* rpz) const
647{
648 if (!mNumberOfParameterizationTPCRat) {
649 return -1;
650 }
651 int rid, pid, zid = TMath::BinarySearch(mNumberOfDistinctZSegmentsTPCRat, mCoordinatesSegmentsZTPCRat,
652 (Float_t)rpz[2]); // find zsegment
653
654 Bool_t reCheck = kFALSE;
655 while (true) {
656 int psegBeg = mBeginningOfSegmentsPTPCRat[zid];
657
658 for (pid = 0; pid < mNumberOfSegmentsPTPCRat[zid]; pid++) {
659 if (rpz[1] < mCoordinatesSegmentsPTPCRat[psegBeg + pid]) {
660 break;
661 }
662 }
663 if (--pid < 0) {
664 pid = 0;
665 }
666 pid += psegBeg;
667
668 int rsegBeg = mBeginningOfSegmentsRTPCRat[pid];
669 for (rid = 0; rid < mNumberOfRSegmentsTPCRat[pid]; rid++) {
670 if (rpz[0] < mCoordinatesSegmentsRTPCRat[rsegBeg + rid]) {
671 break;
672 }
673 }
674 if (--rid < 0) {
675 rid = 0;
676 }
677 rid += rsegBeg;
678
679 // to make sure that due to the precision problems we did not pick the next Zbin
680 if (!reCheck && (rpz[2] - mCoordinatesSegmentsZTPCRat[zid] < 3.e-5) && zid &&
681 !getParameterTPCRatIntegral(mSegmentIdTPCRat[rid])->isInside(rpz)) { // check the previous Z bin
682 zid--;
683 reCheck = kTRUE;
684 continue;
685 }
686 break;
687 }
688 return mSegmentIdTPCRat[rid];
689}
690
691void MagneticWrapperChebyshev::getTPCIntegral(const Double_t* xyz, Double_t* b) const
692{
693 static Double_t rphiz[3];
694
695 // TPCInt region
696 // convert coordinates to cyl system
697 cartesianToCylindrical(xyz, rphiz);
698#ifndef _BRING_TO_BOUNDARY_
699 if ((rphiz[2] > getMaxZTPCIntegral() || rphiz[2] < getMinZTPCIntegral()) || rphiz[0] > getMaxRTPCIntegral()) {
700 for (int i = 3; i--;) {
701 b[i] = 0;
702 }
703 return;
704 }
705#endif
706
707 getTPCIntegralCylindrical(rphiz, b);
708
709 // convert field to cartesian system
710 cylindricalToCartesianCylB(rphiz, b, b);
711}
712
713void MagneticWrapperChebyshev::getTPCRatIntegral(const Double_t* xyz, Double_t* b) const
714{
715 static Double_t rphiz[3];
716
717 // TPCRatIntegral region
718 // convert coordinates to cylindrical system
719 cartesianToCylindrical(xyz, rphiz);
720#ifndef _BRING_TO_BOUNDARY_
721 if ((rphiz[2] > getMaxZTPCRatIntegral() || rphiz[2] < getMinZTPCRatIntegral()) ||
722 rphiz[0] > getMaxRTPCRatIntegral()) {
723 for (int i = 3; i--;) {
724 b[i] = 0;
725 }
726 return;
727 }
728#endif
729
730 getTPCRatIntegralCylindrical(rphiz, b);
731
732 // convert field to cartesian system
733 cylindricalToCartesianCylB(rphiz, b, b);
734}
735
736void MagneticWrapperChebyshev::fieldCylindricalSolenoid(const Double_t* rphiz, Double_t* b) const
737{
738 int id = findSolenoidSegment(rphiz);
739 if (id < 0) {
740 return;
741 }
742 Chebyshev3D* par = getParameterSolenoid(id);
743#ifndef _BRING_TO_BOUNDARY_ // exact matching to fitted volume is requested
744 if (!par->isInside(rphiz)) {
745 return;
746 }
747#endif
748 par->Eval(rphiz, b);
749 return;
750}
751
752Double_t MagneticWrapperChebyshev::fieldCylindricalSolenoidBz(const Double_t* rphiz) const
753{
754 int id = findSolenoidSegment(rphiz);
755 if (id < 0) {
756 return 0.;
757 }
758 Chebyshev3D* par = getParameterSolenoid(id);
759#ifndef _BRING_TO_BOUNDARY_
760 return par->isInside(rphiz) ? par->Eval(rphiz, 2) : 0;
761#else
762 return par->Eval(rphiz, 2);
763#endif
764}
765
766void MagneticWrapperChebyshev::getTPCIntegralCylindrical(const Double_t* rphiz, Double_t* b) const
767{
768 int id = findTPCSegment(rphiz);
769 if (id < 0) {
770 b[0] = b[1] = b[2] = 0;
771 return;
772 }
773 if (id >= mNumberOfParameterizationTPC) {
774 LOG(error) << "MagneticWrapperChebyshev::getTPCIntegralCylindrical: Wrong TPCParam segment " << id;
775 b[0] = b[1] = b[2] = 0;
776 return;
777 }
778 Chebyshev3D* par = getParameterTPCIntegral(id);
779 if (par->isInside(rphiz)) {
780 par->Eval(rphiz, b);
781 return;
782 }
783 b[0] = b[1] = b[2] = 0;
784 return;
785}
786
787void MagneticWrapperChebyshev::getTPCRatIntegralCylindrical(const Double_t* rphiz, Double_t* b) const
788{
789 int id = findTPCRatSegment(rphiz);
790 if (id < 0) {
791 b[0] = b[1] = b[2] = 0;
792 return;
793 }
794 if (id >= mNumberOfParameterizationTPCRat) {
795 LOG(error) << "MagneticWrapperChebyshev::getTPCRatIntegralCylindrical: Wrong TPCRatParam segment " << id;
796 b[0] = b[1] = b[2] = 0;
797 return;
798 }
799 Chebyshev3D* par = getParameterTPCRatIntegral(id);
800 if (par->isInside(rphiz)) {
801 par->Eval(rphiz, b);
802 return;
803 }
804 b[0] = b[1] = b[2] = 0;
805 return;
806}
807
808void checkExpected(char const* expected, TString& buffs)
809{
810 if (!buffs.BeginsWith(expected)) {
811 LOG(error) << R"(MagneticWrapperChebyshev::loadData: Expected: ")" << expected << R"( <name>", found ")" << buffs.Data() << "\"\nStop\n";
812 exit(1);
813 }
814}
815
816#ifdef _INC_CREATION_Chebyshev3D_
817
818void MagneticWrapperChebyshev::loadData(const char* inpfile)
819{
820 TString strf = inpfile;
821 gSystem->ExpandPathName(strf);
822 FILE* stream = fopen(strf, "r");
823
824 if (!stream) {
825 printf("Did not find input file %s\n", strf.Data());
826 return;
827 }
828
829 TString buffs;
830 Chebyshev3DCalc::readLine(buffs, stream);
831
832 checkExpected("START", buffs);
833
834 if (buffs.First(' ') > 0) {
835 SetName(buffs.Data() + buffs.First(' ') + 1);
836 }
837
838 // Solenoid part
839 Chebyshev3DCalc::readLine(buffs, stream);
840
841 checkExpected("START SOLENOID", buffs);
842
843 Chebyshev3DCalc::readLine(buffs, stream); // nparam
844 int nparSol = buffs.Atoi();
845
846 for (int ip = 0; ip < nparSol; ip++) {
847 auto* cheb = new Chebyshev3D();
848 cheb->loadData(stream);
849 addParameterSolenoid(cheb);
850 }
851
852 Chebyshev3DCalc::readLine(buffs, stream);
853 checkExpected("END SOLENOID", buffs);
854
855 // TPCInt part
856 Chebyshev3DCalc::readLine(buffs, stream);
857 checkExpected("START TPCINT", buffs);
858
859 Chebyshev3DCalc::readLine(buffs, stream); // nparam
860 int nparTPCInt = buffs.Atoi();
861
862 for (int ip = 0; ip < nparTPCInt; ip++) {
863 auto* cheb = new Chebyshev3D();
864 cheb->loadData(stream);
865 addParameterTPCIntegral(cheb);
866 }
867
868 Chebyshev3DCalc::readLine(buffs, stream);
869 checkExpected("END TPCINT", buffs);
870
871 // TPCRatInt part
872 Chebyshev3DCalc::readLine(buffs, stream);
873 checkExpected("START TPCRatINT", buffs);
874
875 Chebyshev3DCalc::readLine(buffs, stream); // nparam
876 int nparTPCRatInt = buffs.Atoi();
877
878 for (int ip = 0; ip < nparTPCRatInt; ip++) {
879 auto* cheb = new Chebyshev3D();
880 cheb->loadData(stream);
881 addParameterTPCRatIntegral(cheb);
882 }
883
884 Chebyshev3DCalc::readLine(buffs, stream);
885 checkExpected("END TPCRatINT", buffs);
886
887 // Dipole part
888 Chebyshev3DCalc::readLine(buffs, stream);
889 checkExpected("START DIPOLE", buffs);
890
891 Chebyshev3DCalc::readLine(buffs, stream); // nparam
892 int nparDip = buffs.Atoi();
893
894 for (int ip = 0; ip < nparDip; ip++) {
895 auto* cheb = new Chebyshev3D();
896 cheb->loadData(stream);
897 addParameterDipole(cheb);
898 }
899
900 Chebyshev3DCalc::readLine(buffs, stream);
901 checkExpected("END DIPOLE", buffs);
902
903 Chebyshev3DCalc::readLine(buffs, stream);
904
905 if (!buffs.BeginsWith("END ") && !buffs.Contains(GetName())) {
906 LOG(error) << R"(MagneticWrapperChebyshev::loadData: Expected: "END )"
907 << GetName() << R"( ", found ")" << buffs.Data() << "\"\nStop\n";
908 exit(1);
909 }
910
911 fclose(stream);
912 buildTableSolenoid();
913 buildTableDipole();
914 buildTableTPCIntegral();
915 buildTableTPCRatIntegral();
916
917 printf("Loaded magnetic field \"%s\" from %s\n", GetName(), strf.Data());
918}
919
920void MagneticWrapperChebyshev::buildTableSolenoid()
921{
922 buildTable(mNumberOfParameterizationSolenoid, mParameterizationSolenoid, mNumberOfDistinctZSegmentsSolenoid,
923 mNumberOfDistinctPSegmentsSolenoid, mNumberOfDistinctRSegmentsSolenoid, mMinZSolenoid, mMaxZSolenoid,
924 &mCoordinatesSegmentsZSolenoid, &mCoordinatesSegmentsPSolenoid, &mCoordinatesSegmentsRSolenoid,
925 &mBeginningOfSegmentsPSolenoid, &mNumberOfSegmentsPSolenoid, &mBeginningOfSegmentsRSolenoid,
926 &mNumberOfRSegmentsSolenoid, &mSegmentIdSolenoid);
927}
928
929void MagneticWrapperChebyshev::buildTableDipole()
930{
931 buildTable(mNumberOfParameterizationDipole, mParameterizationDipole, mNumberOfDistinctZSegmentsDipole,
932 mNumberOfDistinctYSegmentsDipole, mNumberOfDistinctXSegmentsDipole, mMinDipoleZ, mMaxDipoleZ,
933 &mCoordinatesSegmentsZDipole, &mCoordinatesSegmentsYDipole, &mCoordinatesSegmentsXDipole,
934 &mBeginningOfSegmentsYDipole, &mNumberOfSegmentsYDipole, &mBeginningOfSegmentsXDipole,
935 &mNumberOfSegmentsXDipole, &mSegmentIdDipole);
936}
937
938void MagneticWrapperChebyshev::buildTableTPCIntegral()
939{
940 buildTable(mNumberOfParameterizationTPC, mParameterizationTPC, mNumberOfDistinctZSegmentsTPC,
941 mNumberOfDistinctPSegmentsTPC, mNumberOfDistinctRSegmentsTPC, mMinZTPC, mMaxZTPC,
942 &mCoordinatesSegmentsZTPC, &mCoordinatesSegmentsPTPC, &mCoordinatesSegmentsRTPC, &mBeginningOfSegmentsPTPC,
943 &mNumberOfSegmentsPTPC, &mBeginningOfSegmentsRTPC, &mNumberOfRSegmentsTPC, &mSegmentIdTPC);
944}
945
946void MagneticWrapperChebyshev::buildTableTPCRatIntegral()
947{
948 buildTable(mNumberOfParameterizationTPCRat, mParameterizationTPCRat, mNumberOfDistinctZSegmentsTPCRat,
949 mNumberOfDistinctPSegmentsTPCRat, mNumberOfDistinctRSegmentsTPCRat, mMinZTPCRat, mMaxZTPCRat,
950 &mCoordinatesSegmentsZTPCRat, &mCoordinatesSegmentsPTPCRat, &mCoordinatesSegmentsRTPCRat,
951 &mBeginningOfSegmentsPTPCRat, &mNumberOfSegmentsPTPCRat, &mBeginningOfSegmentsRTPCRat,
952 &mNumberOfRSegmentsTPCRat, &mSegmentIdTPCRat);
953}
954
955#endif
956
957#ifdef _INC_CREATION_Chebyshev3D_
958MagneticWrapperChebyshev::MagneticWrapperChebyshev(const char* inputFile)
959 : mNumberOfParameterizationSolenoid(0),
960 mNumberOfDistinctZSegmentsSolenoid(0),
961 mNumberOfDistinctPSegmentsSolenoid(0),
962 mNumberOfDistinctRSegmentsSolenoid(0),
963 mCoordinatesSegmentsZSolenoid(nullptr),
964 mCoordinatesSegmentsPSolenoid(nullptr),
965 mCoordinatesSegmentsRSolenoid(nullptr),
966 mBeginningOfSegmentsPSolenoid(nullptr),
967 mNumberOfSegmentsPSolenoid(nullptr),
968 mBeginningOfSegmentsRSolenoid(nullptr),
969 mNumberOfRSegmentsSolenoid(nullptr),
970 mSegmentIdSolenoid(nullptr),
971 mMinZSolenoid(1.e6),
972 mMaxZSolenoid(-1.e6),
973 mParameterizationSolenoid(nullptr),
974 mMaxRadiusSolenoid(0),
975 mNumberOfParameterizationTPC(0),
976 mNumberOfDistinctZSegmentsTPC(0),
977 mNumberOfDistinctPSegmentsTPC(0),
978 mNumberOfDistinctRSegmentsTPC(0),
979 mCoordinatesSegmentsZTPC(nullptr),
980 mCoordinatesSegmentsPTPC(nullptr),
981 mCoordinatesSegmentsRTPC(nullptr),
982 mBeginningOfSegmentsPTPC(nullptr),
983 mNumberOfSegmentsPTPC(nullptr),
984 mBeginningOfSegmentsRTPC(nullptr),
985 mNumberOfRSegmentsTPC(nullptr),
986 mSegmentIdTPC(nullptr),
987 mMinZTPC(1.e6),
988 mMaxZTPC(-1.e6),
989 mParameterizationTPC(nullptr),
990 mMaxRadiusTPC(0),
991 mNumberOfParameterizationTPCRat(0),
992 mNumberOfDistinctZSegmentsTPCRat(0),
993 mNumberOfDistinctPSegmentsTPCRat(0),
994 mNumberOfDistinctRSegmentsTPCRat(0),
995 mCoordinatesSegmentsZTPCRat(nullptr),
996 mCoordinatesSegmentsPTPCRat(nullptr),
997 mCoordinatesSegmentsRTPCRat(nullptr),
998 mBeginningOfSegmentsPTPCRat(nullptr),
999 mNumberOfSegmentsPTPCRat(nullptr),
1000 mBeginningOfSegmentsRTPCRat(nullptr),
1001 mNumberOfRSegmentsTPCRat(nullptr),
1002 mSegmentIdTPCRat(nullptr),
1003 mMinZTPCRat(1.e6),
1004 mMaxZTPCRat(-1.e6),
1005 mParameterizationTPCRat(nullptr),
1006 mMaxRadiusTPCRat(0),
1007 mNumberOfParameterizationDipole(0),
1008 mNumberOfDistinctZSegmentsDipole(0),
1009 mNumberOfDistinctYSegmentsDipole(0),
1010 mNumberOfDistinctXSegmentsDipole(0),
1011 mCoordinatesSegmentsZDipole(nullptr),
1012 mCoordinatesSegmentsYDipole(nullptr),
1013 mCoordinatesSegmentsXDipole(nullptr),
1014 mBeginningOfSegmentsYDipole(nullptr),
1015 mNumberOfSegmentsYDipole(nullptr),
1016 mBeginningOfSegmentsXDipole(nullptr),
1017 mNumberOfSegmentsXDipole(nullptr),
1018 mSegmentIdDipole(nullptr),
1019 mMinDipoleZ(1.e6),
1020 mMaxDipoleZ(-1.e6),
1021 mParameterizationDipole(nullptr)
1022{
1023 loadData(inputFile);
1024}
1025
1026void MagneticWrapperChebyshev::addParameterSolenoid(const Chebyshev3D* param)
1027{
1028 if (!mParameterizationSolenoid) {
1029 mParameterizationSolenoid = new TObjArray();
1030 }
1031 mParameterizationSolenoid->Add((Chebyshev3D*)param);
1032 mNumberOfParameterizationSolenoid++;
1033 if (mMaxRadiusSolenoid < param->getBoundMax(0)) {
1034 mMaxRadiusSolenoid = param->getBoundMax(0);
1035 }
1036}
1037
1038void MagneticWrapperChebyshev::addParameterTPCIntegral(const Chebyshev3D* param)
1039{
1040 if (!mParameterizationTPC) {
1041 mParameterizationTPC = new TObjArray();
1042 }
1043 mParameterizationTPC->Add((Chebyshev3D*)param);
1044 mNumberOfParameterizationTPC++;
1045 if (mMaxRadiusTPC < param->getBoundMax(0)) {
1046 mMaxRadiusTPC = param->getBoundMax(0);
1047 }
1048}
1049
1050void MagneticWrapperChebyshev::addParameterTPCRatIntegral(const Chebyshev3D* param)
1051{
1052 if (!mParameterizationTPCRat) {
1053 mParameterizationTPCRat = new TObjArray();
1054 }
1055 mParameterizationTPCRat->Add((Chebyshev3D*)param);
1056 mNumberOfParameterizationTPCRat++;
1057 if (mMaxRadiusTPCRat < param->getBoundMax(0)) {
1058 mMaxRadiusTPCRat = param->getBoundMax(0);
1059 }
1060}
1061
1062void MagneticWrapperChebyshev::addParameterDipole(const Chebyshev3D* param)
1063{
1064 if (!mParameterizationDipole) {
1065 mParameterizationDipole = new TObjArray();
1066 }
1067 mParameterizationDipole->Add((Chebyshev3D*)param);
1068 mNumberOfParameterizationDipole++;
1069}
1070
1071void MagneticWrapperChebyshev::resetDipole()
1072{
1073 if (mNumberOfParameterizationDipole) {
1074 delete mParameterizationDipole;
1075 mParameterizationDipole = nullptr;
1076 delete[] mCoordinatesSegmentsZDipole;
1077 mCoordinatesSegmentsZDipole = nullptr;
1078 delete[] mCoordinatesSegmentsXDipole;
1079 mCoordinatesSegmentsXDipole = nullptr;
1080 delete[] mCoordinatesSegmentsYDipole;
1081 mCoordinatesSegmentsYDipole = nullptr;
1082 delete[] mBeginningOfSegmentsYDipole;
1083 mBeginningOfSegmentsYDipole = nullptr;
1084 delete[] mNumberOfSegmentsYDipole;
1085 mNumberOfSegmentsYDipole = nullptr;
1086 delete[] mBeginningOfSegmentsXDipole;
1087 mBeginningOfSegmentsXDipole = nullptr;
1088 delete[] mNumberOfSegmentsXDipole;
1089 mNumberOfSegmentsXDipole = nullptr;
1090 delete[] mSegmentIdDipole;
1091 mSegmentIdDipole = nullptr;
1092 }
1093 mNumberOfParameterizationDipole = mNumberOfDistinctZSegmentsDipole = mNumberOfDistinctXSegmentsDipole =
1094 mNumberOfDistinctYSegmentsDipole = 0;
1095 mMinDipoleZ = 1e6;
1096 mMaxDipoleZ = -1e6;
1097}
1098
1099void MagneticWrapperChebyshev::resetSolenoid()
1100{
1101 if (mNumberOfParameterizationSolenoid) {
1102 delete mParameterizationSolenoid;
1103 mParameterizationSolenoid = nullptr;
1104 delete[] mCoordinatesSegmentsZSolenoid;
1105 mCoordinatesSegmentsZSolenoid = nullptr;
1106 delete[] mCoordinatesSegmentsPSolenoid;
1107 mCoordinatesSegmentsPSolenoid = nullptr;
1108 delete[] mCoordinatesSegmentsRSolenoid;
1109 mCoordinatesSegmentsRSolenoid = nullptr;
1110 delete[] mBeginningOfSegmentsPSolenoid;
1111 mBeginningOfSegmentsPSolenoid = nullptr;
1112 delete[] mNumberOfSegmentsPSolenoid;
1113 mNumberOfSegmentsPSolenoid = nullptr;
1114 delete[] mBeginningOfSegmentsRSolenoid;
1115 mBeginningOfSegmentsRSolenoid = nullptr;
1116 delete[] mNumberOfRSegmentsSolenoid;
1117 mNumberOfRSegmentsSolenoid = nullptr;
1118 delete[] mSegmentIdSolenoid;
1119 mSegmentIdSolenoid = nullptr;
1120 }
1121 mNumberOfParameterizationSolenoid = mNumberOfDistinctZSegmentsSolenoid = mNumberOfDistinctPSegmentsSolenoid =
1122 mNumberOfDistinctRSegmentsSolenoid = 0;
1123 mMinZSolenoid = 1e6;
1124 mMaxZSolenoid = -1e6;
1125 mMaxRadiusSolenoid = 0;
1126}
1127
1128void MagneticWrapperChebyshev::resetTPCIntegral()
1129{
1130 if (mNumberOfParameterizationTPC) {
1131 delete mParameterizationTPC;
1132 mParameterizationTPC = nullptr;
1133 delete[] mCoordinatesSegmentsZTPC;
1134 mCoordinatesSegmentsZTPC = nullptr;
1135 delete[] mCoordinatesSegmentsPTPC;
1136 mCoordinatesSegmentsPTPC = nullptr;
1137 delete[] mCoordinatesSegmentsRTPC;
1138 mCoordinatesSegmentsRTPC = nullptr;
1139 delete[] mBeginningOfSegmentsPTPC;
1140 mBeginningOfSegmentsPTPC = nullptr;
1141 delete[] mNumberOfSegmentsPTPC;
1142 mNumberOfSegmentsPTPC = nullptr;
1143 delete[] mBeginningOfSegmentsRTPC;
1144 mBeginningOfSegmentsRTPC = nullptr;
1145 delete[] mNumberOfRSegmentsTPC;
1146 mNumberOfRSegmentsTPC = nullptr;
1147 delete[] mSegmentIdTPC;
1148 mSegmentIdTPC = nullptr;
1149 }
1150 mNumberOfParameterizationTPC = mNumberOfDistinctZSegmentsTPC = mNumberOfDistinctPSegmentsTPC =
1151 mNumberOfDistinctRSegmentsTPC = 0;
1152 mMinZTPC = 1e6;
1153 mMaxZTPC = -1e6;
1154 mMaxRadiusTPC = 0;
1155}
1156
1157void MagneticWrapperChebyshev::resetTPCRatIntegral()
1158{
1159 if (mNumberOfParameterizationTPCRat) {
1160 delete mParameterizationTPCRat;
1161 mParameterizationTPCRat = nullptr;
1162 delete[] mCoordinatesSegmentsZTPCRat;
1163 mCoordinatesSegmentsZTPCRat = nullptr;
1164 delete[] mCoordinatesSegmentsPTPCRat;
1165 mCoordinatesSegmentsPTPCRat = nullptr;
1166 delete[] mCoordinatesSegmentsRTPCRat;
1167 mCoordinatesSegmentsRTPCRat = nullptr;
1168 delete[] mBeginningOfSegmentsPTPCRat;
1169 mBeginningOfSegmentsPTPCRat = nullptr;
1170 delete[] mNumberOfSegmentsPTPCRat;
1171 mNumberOfSegmentsPTPCRat = nullptr;
1172 delete[] mBeginningOfSegmentsRTPCRat;
1173 mBeginningOfSegmentsRTPCRat = nullptr;
1174 delete[] mNumberOfRSegmentsTPCRat;
1175 mNumberOfRSegmentsTPCRat = nullptr;
1176 delete[] mSegmentIdTPCRat;
1177 mSegmentIdTPCRat = nullptr;
1178 }
1179 mNumberOfParameterizationTPCRat = mNumberOfDistinctZSegmentsTPCRat = mNumberOfDistinctPSegmentsTPCRat =
1180 mNumberOfDistinctRSegmentsTPCRat = 0;
1181 mMinZTPCRat = 1e6;
1182 mMaxZTPCRat = -1e6;
1183 mMaxRadiusTPCRat = 0;
1184}
1185
1186void MagneticWrapperChebyshev::buildTable(Int_t npar, TObjArray* parArr, Int_t& nZSeg, Int_t& nYSeg, Int_t& nXSeg,
1187 Float_t& minZ, Float_t& maxZ, Float_t** segZ, Float_t** segY, Float_t** segX,
1188 Int_t** begSegY, Int_t** nSegY, Int_t** begSegX, Int_t** nSegX, Int_t** segID)
1189{
1190 if (npar < 1) {
1191 return;
1192 }
1193 TArrayF segYArr, segXArr;
1194 TArrayI begSegYDipArr, begSegXDipArr;
1195 TArrayI nSegYDipArr, nSegXDipArr;
1196 TArrayI segIDArr;
1197 float *tmpSegZ, *tmpSegY, *tmpSegX;
1198
1199 // create segmentation in Z
1200 nZSeg = segmentDimension(&tmpSegZ, parArr, npar, 2, 1, -1, 1, -1, 1, -1) - 1;
1201 nYSeg = 0;
1202 nXSeg = 0;
1203
1204 // for each Z slice create segmentation in Y
1205 begSegYDipArr.Set(nZSeg);
1206 nSegYDipArr.Set(nZSeg);
1207 float xyz[3];
1208 for (int iz = 0; iz < nZSeg; iz++) {
1209 LOGF(debug, "\nZSegment#%d %+e : %+e\n", iz, tmpSegZ[iz], tmpSegZ[iz + 1]);
1210 int ny = segmentDimension(&tmpSegY, parArr, npar, 1, 1, -1, 1, -1, tmpSegZ[iz], tmpSegZ[iz + 1]) - 1;
1211 segYArr.Set(ny + nYSeg);
1212 for (int iy = 0; iy < ny; iy++) {
1213 segYArr[nYSeg + iy] = tmpSegY[iy];
1214 }
1215 begSegYDipArr[iz] = nYSeg;
1216 nSegYDipArr[iz] = ny;
1217 LOGF(debug, " Found %d YSegments, to start from %d\n", ny, begSegYDipArr[iz]);
1218
1219 // for each slice in Z and Y create segmentation in X
1220 begSegXDipArr.Set(nYSeg + ny);
1221 nSegXDipArr.Set(nYSeg + ny);
1222 xyz[2] = (tmpSegZ[iz] + tmpSegZ[iz + 1]) / 2.; // mean Z of this segment
1223
1224 for (int iy = 0; iy < ny; iy++) {
1225 int isg = nYSeg + iy;
1226 LOGF(debug, "\n YSegment#%d %+e : %+e\n", iy, tmpSegY[iy], tmpSegY[iy + 1]);
1227 int nx =
1228 segmentDimension(&tmpSegX, parArr, npar, 0, 1, -1, tmpSegY[iy], tmpSegY[iy + 1], tmpSegZ[iz], tmpSegZ[iz + 1]) -
1229 1;
1230
1231 segXArr.Set(nx + nXSeg);
1232 for (int ix = 0; ix < nx; ix++) {
1233 segXArr[nXSeg + ix] = tmpSegX[ix];
1234 }
1235 begSegXDipArr[isg] = nXSeg;
1236 nSegXDipArr[isg] = nx;
1237 LOGF(debug, " Found %d XSegments, to start from %d\n", nx, begSegXDipArr[isg]);
1238
1239 segIDArr.Set(nXSeg + nx);
1240
1241 // find corresponding params
1242 xyz[1] = (tmpSegY[iy] + tmpSegY[iy + 1]) / 2.; // mean Y of this segment
1243
1244 for (int ix = 0; ix < nx; ix++) {
1245 xyz[0] = (tmpSegX[ix] + tmpSegX[ix + 1]) / 2.; // mean X of this segment
1246 for (int ipar = 0; ipar < npar; ipar++) {
1247 Chebyshev3D* cheb = (Chebyshev3D*)parArr->At(ipar);
1248 if (!cheb->isInside(xyz)) {
1249 continue;
1250 }
1251 segIDArr[nXSeg + ix] = ipar;
1252 break;
1253 }
1254 }
1255 nXSeg += nx;
1256
1257 delete[] tmpSegX;
1258 }
1259 delete[] tmpSegY;
1260 nYSeg += ny;
1261 }
1262
1263 minZ = tmpSegZ[0];
1264 maxZ = tmpSegZ[nZSeg];
1265 (*segZ) = new Float_t[nZSeg];
1266 for (int i = nZSeg; i--;) {
1267 (*segZ)[i] = tmpSegZ[i];
1268 }
1269 delete[] tmpSegZ;
1270
1271 (*segY) = new Float_t[nYSeg];
1272 (*segX) = new Float_t[nXSeg];
1273 (*begSegY) = new Int_t[nZSeg];
1274 (*nSegY) = new Int_t[nZSeg];
1275 (*begSegX) = new Int_t[nYSeg];
1276 (*nSegX) = new Int_t[nYSeg];
1277 (*segID) = new Int_t[nXSeg];
1278
1279 for (int i = nYSeg; i--;) {
1280 (*segY)[i] = segYArr[i];
1281 }
1282 for (int i = nXSeg; i--;) {
1283 (*segX)[i] = segXArr[i];
1284 }
1285 for (int i = nZSeg; i--;) {
1286 (*begSegY)[i] = begSegYDipArr[i];
1287 (*nSegY)[i] = nSegYDipArr[i];
1288 }
1289 for (int i = nYSeg; i--;) {
1290 (*begSegX)[i] = begSegXDipArr[i];
1291 (*nSegX)[i] = nSegXDipArr[i];
1292 }
1293 for (int i = nXSeg; i--;) {
1294 (*segID)[i] = segIDArr[i];
1295 }
1296}
1297
1298// void MagneticWrapperChebyshev::BuildTableDip()
1299// {
1300// // build lookup table for dipole
1301//
1302// if (mNumberOfParameterizationDipole<1) return;
1303// TArrayF segY,segX;
1304// TArrayI begSegYDip,begSegXDip;
1305// TArrayI nsegYDip,nsegXDip;
1306// TArrayI segID;
1307// float *tmpSegZ,*tmpSegY,*tmpSegX;
1308//
1309// // create segmentation in Z
1310// mNumberOfDistinctZSegmentsDipole = segmentDimension(&tmpSegZ, mParameterizationDipole,
1311// mNumberOfParameterizationDipole, 2, 1,-1, 1,-1, 1,-1) - 1;
1312// mNumberOfDistinctYSegmentsDipole = 0;
1313// mNumberOfDistinctXSegmentsDipole = 0;
1314//
1315// // for each Z slice create segmentation in Y
1316// begSegYDip.Set(mNumberOfDistinctZSegmentsDipole);
1317// nsegYDip.Set(mNumberOfDistinctZSegmentsDipole);
1318// float xyz[3];
1319// for (int iz=0;iz<mNumberOfDistinctZSegmentsDipole;iz++) {
1320// printf("\nZSegment#%d %+e : %+e\n",iz,tmpSegZ[iz],tmpSegZ[iz+1]);
1321// int ny = segmentDimension(&tmpSegY, mParameterizationDipole, mNumberOfParameterizationDipole, 1,
1322// 1,-1, 1,-1, tmpSegZ[iz],tmpSegZ[iz+1]) - 1;
1323// segY.Set(ny + mNumberOfDistinctYSegmentsDipole);
1324// for (int iy=0;iy<ny;iy++) segY[mNumberOfDistinctYSegmentsDipole+iy] = tmpSegY[iy];
1325// begSegYDip[iz] = mNumberOfDistinctYSegmentsDipole;
1326// nsegYDip[iz] = ny;
1327// printf(" Found %d YSegments, to start from %d\n",ny, begSegYDip[iz]);
1328//
1329// // for each slice in Z and Y create segmentation in X
1330// begSegXDip.Set(mNumberOfDistinctYSegmentsDipole+ny);
1331// nsegXDip.Set(mNumberOfDistinctYSegmentsDipole+ny);
1332// xyz[2] = (tmpSegZ[iz]+tmpSegZ[iz+1])/2.; // mean Z of this segment
1333//
1334// for (int iy=0;iy<ny;iy++) {
1335// int isg = mNumberOfDistinctYSegmentsDipole+iy;
1336// printf("\n YSegment#%d %+e : %+e\n",iy, tmpSegY[iy],tmpSegY[iy+1]);
1337// int nx = segmentDimension(&tmpSegX, mParameterizationDipole, mNumberOfParameterizationDipole, 0,
1338// 1,-1, tmpSegY[iy],tmpSegY[iy+1], tmpSegZ[iz],tmpSegZ[iz+1]) - 1;
1339//
1340// segX.Set(nx + mNumberOfDistinctXSegmentsDipole);
1341// for (int ix=0;ix<nx;ix++) segX[mNumberOfDistinctXSegmentsDipole+ix] = tmpSegX[ix];
1342// begSegXDip[isg] = mNumberOfDistinctXSegmentsDipole;
1343// nsegXDip[isg] = nx;
1344// printf(" Found %d XSegments, to start from %d\n",nx, begSegXDip[isg]);
1345//
1346// segID.Set(mNumberOfDistinctXSegmentsDipole+nx);
1347//
1348// // find corresponding params
1349// xyz[1] = (tmpSegY[iy]+tmpSegY[iy+1])/2.; // mean Y of this segment
1350//
1351// for (int ix=0;ix<nx;ix++) {
1352// xyz[0] = (tmpSegX[ix]+tmpSegX[ix+1])/2.; // mean X of this segment
1353// for (int ipar=0;ipar<mNumberOfParameterizationDipole;ipar++) {
1354// Chebyshev3D* cheb = (Chebyshev3D*) mParameterizationDipole->At(ipar);
1355// if (!cheb->isInside(xyz)) continue;
1356// segID[mNumberOfDistinctXSegmentsDipole+ix] = ipar;
1357// break;
1358// }
1359// }
1360// mNumberOfDistinctXSegmentsDipole += nx;
1361//
1362// delete[] tmpSegX;
1363// }
1364// delete[] tmpSegY;
1365// mNumberOfDistinctYSegmentsDipole += ny;
1366// }
1367//
1368// mMinDipoleZ = tmpSegZ[0];
1369// mMaxDipoleZ = tmpSegZ[mNumberOfDistinctZSegmentsDipole];
1370// mCoordinatesSegmentsZDipole = new Float_t[mNumberOfDistinctZSegmentsDipole];
1371// for (int i=mNumberOfDistinctZSegmentsDipole;i--;) mCoordinatesSegmentsZDipole[i] = tmpSegZ[i];
1372// delete[] tmpSegZ;
1373//
1374// mCoordinatesSegmentsYDipole = new Float_t[mNumberOfDistinctYSegmentsDipole];
1375// mCoordinatesSegmentsXDipole = new Float_t[mNumberOfDistinctXSegmentsDipole];
1376// mBeginningOfSegmentsYDipole = new Int_t[mNumberOfDistinctZSegmentsDipole];
1377// mNumberOfSegmentsYDipole = new Int_t[mNumberOfDistinctZSegmentsDipole];
1378// mBeginningOfSegmentsXDipole = new Int_t[mNumberOfDistinctYSegmentsDipole];
1379// mNumberOfSegmentsXDipole = new Int_t[mNumberOfDistinctYSegmentsDipole];
1380// mSegmentIdDipole = new Int_t[mNumberOfDistinctXSegmentsDipole];
1381//
1382// for (int i=mNumberOfDistinctYSegmentsDipole;i--;) mCoordinatesSegmentsYDipole[i] = segY[i];
1383// for (int i=mNumberOfDistinctXSegmentsDipole;i--;) mCoordinatesSegmentsXDipole[i] = segX[i];
1384// for (int i=mNumberOfDistinctZSegmentsDipole;i--;) {mBeginningOfSegmentsYDipole[i] = begSegYDip[i];
1385// mNumberOfSegmentsYDipole[i] = nsegYDip[i];}
1386// for (int i=mNumberOfDistinctYSegmentsDipole;i--;) {mBeginningOfSegmentsXDipole[i] = begSegXDip[i];
1387// mNumberOfSegmentsXDipole[i] = nsegXDip[i];}
1388// for (int i=mNumberOfDistinctXSegmentsDipole;i--;) {mSegmentIdDipole[i] = segID[i];}
1389// }
1390
1391void MagneticWrapperChebyshev::saveData(const char* outfile) const
1392{
1393 TString strf = outfile;
1394 gSystem->ExpandPathName(strf);
1395 FILE* stream = fopen(strf, "w+");
1396
1397 // Solenoid part
1398 fprintf(stream, "# Set of Chebyshev parameterizations for ALICE magnetic field\nSTART %s\n", GetName());
1399 fprintf(stream, "START SOLENOID\n#Number of pieces\n%d\n", mNumberOfParameterizationSolenoid);
1400 for (int ip = 0; ip < mNumberOfParameterizationSolenoid; ip++) {
1401 getParameterSolenoid(ip)->saveData(stream);
1402 }
1403 fprintf(stream, "#\nEND SOLENOID\n");
1404
1405 // TPCIntegral part
1406 fprintf(stream, "START TPCINT\n#Number of pieces\n%d\n", mNumberOfParameterizationTPC);
1407 for (int ip = 0; ip < mNumberOfParameterizationTPC; ip++) {
1408 getParameterTPCIntegral(ip)->saveData(stream);
1409 }
1410 fprintf(stream, "#\nEND TPCINT\n");
1411
1412 // TPCRatIntegral part
1413 fprintf(stream, "START TPCRatINT\n#Number of pieces\n%d\n", mNumberOfParameterizationTPCRat);
1414 for (int ip = 0; ip < mNumberOfParameterizationTPCRat; ip++) {
1415 getParameterTPCRatIntegral(ip)->saveData(stream);
1416 }
1417 fprintf(stream, "#\nEND TPCRatINT\n");
1418
1419 // Dipole part
1420 fprintf(stream, "START DIPOLE\n#Number of pieces\n%d\n", mNumberOfParameterizationDipole);
1421 for (int ip = 0; ip < mNumberOfParameterizationDipole; ip++) {
1422 getParameterDipole(ip)->saveData(stream);
1423 }
1424 fprintf(stream, "#\nEND DIPOLE\n");
1425
1426 fprintf(stream, "#\nEND %s\n", GetName());
1427
1428 fclose(stream);
1429}
1430
1431Int_t MagneticWrapperChebyshev::segmentDimension(float** seg, const TObjArray* par, int npar, int dim, float xmn,
1432 float xmx, float ymn, float ymx, float zmn, float zmx)
1433{
1434 auto* tmpC = new float[2 * npar];
1435 auto* tmpInd = new int[2 * npar];
1436 int nseg0 = 0;
1437 for (int ip = 0; ip < npar; ip++) {
1438 Chebyshev3D* cheb = (Chebyshev3D*)par->At(ip);
1439 if (xmn < xmx && (cheb->getBoundMin(0) > (xmx + xmn) / 2 || cheb->getBoundMax(0) < (xmn + xmx) / 2)) {
1440 continue;
1441 }
1442 if (ymn < ymx && (cheb->getBoundMin(1) > (ymx + ymn) / 2 || cheb->getBoundMax(1) < (ymn + ymx) / 2)) {
1443 continue;
1444 }
1445 if (zmn < zmx && (cheb->getBoundMin(2) > (zmx + zmn) / 2 || cheb->getBoundMax(2) < (zmn + zmx) / 2)) {
1446 continue;
1447 }
1448
1449 tmpC[nseg0++] = cheb->getBoundMin(dim);
1450 tmpC[nseg0++] = cheb->getBoundMax(dim);
1451 }
1452 // range Dim's boundaries in increasing order
1453 TMath::Sort(nseg0, tmpC, tmpInd, kFALSE);
1454 // count number of really different Z's
1455 int nseg = 0;
1456 float cprev = -1e6;
1457 for (int ip = 0; ip < nseg0; ip++) {
1458 if (TMath::Abs(cprev - tmpC[tmpInd[ip]]) > 1e-4) {
1459 cprev = tmpC[tmpInd[ip]];
1460 nseg++;
1461 } else {
1462 tmpInd[ip] = -1; // supress redundant Z
1463 }
1464 }
1465
1466 *seg = new float[nseg]; // create final Z segmenations
1467 nseg = 0;
1468 for (int ip = 0; ip < nseg0; ip++) {
1469 if (tmpInd[ip] >= 0) {
1470 (*seg)[nseg++] = tmpC[tmpInd[ip]];
1471 }
1472 }
1473
1474 delete[] tmpC;
1475 delete[] tmpInd;
1476 return nseg;
1477}
1478
1479#endif
seg
o2::mch::mapping::CathodeSegmentation seg
Definition BenchCathodeSegmentation.cxx:92
i
int32_t i
Definition GPUCommonAlgorithm.h:443
checkExpected
void checkExpected(char const *expected, TString &buffs)
Definition MagneticWrapperChebyshev.cxx:808
ClassImp
ClassImp(MagneticWrapperChebyshev)
pid
uint16_t pid
Definition RawData.h:2
debug
std::ostringstream debug
Definition VariantJSONHelpers.h:307
Float_t
o2::field::MagneticWrapperChebyshev::fieldCylindricalSolenoid
void fieldCylindricalSolenoid(const Double_t *rphiz, Double_t *b) const
Definition MagneticWrapperChebyshev.cxx:736
o2::field::MagneticWrapperChebyshev::getParameterDipole
o2::math_utils::Chebyshev3D * getParameterDipole(Int_t ipar) const
Definition MagneticWrapperChebyshev.h:201
o2::field::MagneticWrapperChebyshev::operator=
MagneticWrapperChebyshev & operator=(const MagneticWrapperChebyshev &rhs)
Assignment operator.
Definition MagneticWrapperChebyshev.cxx:296
o2::field::MagneticWrapperChebyshev::Print
void Print(Option_t *="") const override
Prints info.
Definition MagneticWrapperChebyshev.cxx:470
o2::field::MagneticWrapperChebyshev::copyFrom
void copyFrom(const MagneticWrapperChebyshev &src)
Copy method.
Definition MagneticWrapperChebyshev.cxx:170
o2::field::MagneticWrapperChebyshev::getParameterTPCRatIntegral
o2::math_utils::Chebyshev3D * getParameterTPCRatIntegral(Int_t ipar) const
Definition MagneticWrapperChebyshev.h:191
o2::field::MagneticWrapperChebyshev::getBz
Double_t getBz(const Double_t *xyz) const
Definition MagneticWrapperChebyshev.cxx:448
o2::field::MagneticWrapperChebyshev::Clear
void Clear(const Option_t *="") override
Clears all dynamic parts.
Definition MagneticWrapperChebyshev.cxx:305
o2::field::MagneticWrapperChebyshev::fieldCylindricalSolenoidBz
Double_t fieldCylindricalSolenoidBz(const Double_t *rphiz) const
Definition MagneticWrapperChebyshev.cxx:752
o2::field::MagneticWrapperChebyshev::findDipoleSegment
Int_t findDipoleSegment(const Double_t *xyz) const
Finds the segment containing point xyz. If it is outside it finds the closest segment.
Definition MagneticWrapperChebyshev.cxx:511
o2::field::MagneticWrapperChebyshev::getTPCIntegralCylindrical
void getTPCIntegralCylindrical(const Double_t *rphiz, Double_t *b) const
Definition MagneticWrapperChebyshev.cxx:766
o2::field::MagneticWrapperChebyshev::cylindricalToCartesianCylB
static void cylindricalToCartesianCylB(const Double_t *rphiz, const Double_t *brphiz, Double_t *bxyz)
Converts field in cylindrical coordinates to cartesian system, point is in cyl.system.
Definition MagneticWrapperChebyshev.h:430
o2::field::MagneticWrapperChebyshev::getParameterSolenoid
o2::math_utils::Chebyshev3D * getParameterSolenoid(Int_t ipar) const
Definition MagneticWrapperChebyshev.h:186
o2::field::MagneticWrapperChebyshev::getTPCRatIntegral
void getTPCRatIntegral(const Double_t *xyz, Double_t *b) const
Definition MagneticWrapperChebyshev.cxx:713
o2::field::MagneticWrapperChebyshev::findTPCRatSegment
Int_t findTPCRatSegment(const Double_t *xyz) const
Finds the segment containing point xyz. If it is outside it finds the closest segment.
Definition MagneticWrapperChebyshev.cxx:646
o2::field::MagneticWrapperChebyshev::Field
virtual void Field(const Double_t *xyz, Double_t *b) const
Definition MagneticWrapperChebyshev.cxx:419
o2::field::MagneticWrapperChebyshev::getTPCIntegral
void getTPCIntegral(const Double_t *xyz, Double_t *b) const
Definition MagneticWrapperChebyshev.cxx:691
o2::field::MagneticWrapperChebyshev::findTPCSegment
Int_t findTPCSegment(const Double_t *xyz) const
Finds the segment containing point xyz. If it is outside it finds the closest segment.
Definition MagneticWrapperChebyshev.cxx:601
o2::field::MagneticWrapperChebyshev::getParameterTPCIntegral
o2::math_utils::Chebyshev3D * getParameterTPCIntegral(Int_t ipar) const
Definition MagneticWrapperChebyshev.h:196
o2::field::MagneticWrapperChebyshev::findSolenoidSegment
Int_t findSolenoidSegment(const Double_t *xyz) const
Finds the segment containing point xyz. If it is outside it finds the closest segment.
Definition MagneticWrapperChebyshev.cxx:557
o2::field::MagneticWrapperChebyshev::cartesianToCylindrical
static void cartesianToCylindrical(const Double_t *xyz, Double_t *rphiz)
Definition MagneticWrapperChebyshev.h:474
o2::field::MagneticWrapperChebyshev::getTPCRatIntegralCylindrical
void getTPCRatIntegralCylindrical(const Double_t *rphiz, Double_t *b) const
Definition MagneticWrapperChebyshev.cxx:787
o2::math_utils::Chebyshev3DCalc::readLine
static void readLine(TString &str, FILE *stream)
Reads single line from the stream, skipping empty and commented lines. EOF is not expected.
Definition Chebyshev3DCalc.cxx:385
o2::math_utils::Chebyshev3D
Definition Chebyshev3D.h:65
o2::math_utils::Chebyshev3D::isInside
Bool_t isInside(const Float_t *par) const
Checks if the point is inside of the fitted box.
Definition Chebyshev3D.h:271
o2::math_utils::Chebyshev3D::getBoundMin
Float_t getBoundMin(int i) const
Definition Chebyshev3D.h:162
o2::math_utils::Chebyshev3D::Print
void Print(const Option_t *opt="") const override
Definition Chebyshev3D.cxx:310
o2::math_utils::Chebyshev3D::Eval
void Eval(const Float_t *par, Float_t *res)
Evaluates Chebyshev parameterization for 3d->DimOut function.
Definition Chebyshev3D.h:293
o2::math_utils::Chebyshev3D::getBoundMax
Float_t getBoundMax(int i) const
Definition Chebyshev3D.h:167
src
GLenum src
Definition glcorearb.h:1767
b
GLboolean GLboolean GLboolean b
Definition glcorearb.h:1233
param
GLenum GLfloat param
Definition glcorearb.h:271
maxZ
GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat maxZ
Definition glcorearb.h:2910
stream
GLuint GLuint stream
Definition glcorearb.h:1806
id
GLuint id
Definition glcorearb.h:650
minZ
GLfloat GLfloat minZ
Definition glcorearb.h:2910
o2::math_utils
Definition CachingTF1.h:23
o2::track::iz
const int iz
Definition TrackUtils.h:69
o2::track::iy
const int iy
Definition TrackUtils.h:69
o2::track::ix
const int ix
Definition TrackUtils.h:69
loadData
std::tuple< TFile *, TTreeReader * > loadData(const std::string inFile)
Definition statusmap-to-rejectlist.cxx:36
expected
std::map< std::string, ID > expected
Definition testDCSAliases.cxx:166
LOG
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"