CbmRoot
Loading...
Searching...
No Matches
CbmLitFieldFitter.cxx
Go to the documentation of this file.
1/* Copyright (C) 2009-2015 GSI/JINR-LIT, Darmstadt/Dubna
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Andrey Lebedev [committer] */
4
10#include "base/CbmLitFieldFitter.h"
11
12#include "CbmUtils.h"
13#include "FairField.h"
14#include "FairRunAna.h"
15
16//#include "TFitterMinuit.h"
17#include "Math/IFunction.h"
18#include "Minuit2/Minuit2Minimizer.h"
19#include "TMatrixD.h"
20#include "TVectorD.h"
21
22#include <cmath>
23
30class CbmLitPolynom0 : public CbmLitPolynom {
31 public:
35 CbmLitPolynom0() {}
36
40 virtual ~CbmLitPolynom0() {}
41
45 double Calculate(double x, double y, double c[]) const { return c[0]; }
46
50 unsigned int GetNofCoefficients() const { return 1; }
51};
52
59class CbmLitPolynom1 : public CbmLitPolynom {
60 public:
64 CbmLitPolynom1() {}
65
69 virtual ~CbmLitPolynom1() {}
70
74 double Calculate(double x, double y, double c[]) const { return c[0] + c[1] * x + c[2] * y; }
75
79 unsigned int GetNofCoefficients() const { return 3; }
80};
81
88class CbmLitPolynom2 : public CbmLitPolynom {
89 public:
93 CbmLitPolynom2() {}
94
98 virtual ~CbmLitPolynom2() {}
99
103 double Calculate(double x, double y, double c[]) const
104 {
105 return c[0] + c[1] * x + c[2] * y + c[3] * x * x + c[4] * x * y + c[5] * y * y;
106 }
107
111 unsigned int GetNofCoefficients() const { return 6; }
112};
113
120class CbmLitPolynom3 : public CbmLitPolynom {
121 public:
125 CbmLitPolynom3() {}
126
130 virtual ~CbmLitPolynom3() {}
131
135 double Calculate(double x, double y, double c[]) const
136 {
137 double x2 = x * x;
138 double y2 = y * y;
139 double x2y = x2 * y;
140 double xy2 = x * y2;
141 double x3 = x2 * x;
142 double y3 = y2 * y;
143
144 return c[0] + c[1] * x + c[2] * y + c[3] * x2 + c[4] * x * y + c[5] * y2 + c[6] * x3 + c[7] * x2y + c[8] * xy2
145 + c[9] * y3;
146 }
147
151 unsigned int GetNofCoefficients() const { return 10; }
152};
153
160class CbmLitPolynom4 : public CbmLitPolynom {
161 public:
165 CbmLitPolynom4() {}
166
170 virtual ~CbmLitPolynom4() {}
171
175 double Calculate(double x, double y, double c[]) const
176 {
177
178 double x2 = x * x;
179 double y2 = y * y;
180 double xy = x * y;
181
182 double x3 = x2 * x;
183 double y3 = y2 * y;
184 double xy2 = x * y2;
185 double x2y = x2 * y;
186
187 double x4 = x3 * x;
188 double y4 = y3 * y;
189 double xy3 = x * y3;
190 double x2y2 = x2 * y2;
191 double x3y = x3 * y;
192
193 return c[0] + c[1] * x + c[2] * y + c[3] * x2 + c[4] * xy + c[5] * y2 + c[6] * x3 + c[7] * x2y + c[8] * xy2
194 + c[9] * y3 + c[10] * x4 + c[11] * x3y + c[12] * x2y2 + c[13] * xy3 + c[14] * y4;
195 }
196
200 unsigned int GetNofCoefficients() const { return 15; }
201};
202
209class CbmLitPolynom5 : public CbmLitPolynom {
210 public:
214 CbmLitPolynom5() {}
215
219 virtual ~CbmLitPolynom5() {}
220
224 double Calculate(double x, double y, double c[]) const
225 {
226
227 double x2 = x * x;
228 double y2 = y * y;
229 double xy = x * y;
230
231 double x3 = x2 * x;
232 double y3 = y2 * y;
233 double xy2 = x * y2;
234 double x2y = x2 * y;
235
236 double x4 = x3 * x;
237 double y4 = y3 * y;
238 double xy3 = x * y3;
239 double x2y2 = x2 * y2;
240 double x3y = x3 * y;
241
242 double x5 = x4 * x;
243 double y5 = y4 * y;
244 double xy4 = x * y4;
245 double x2y3 = x2 * y3;
246 double x3y2 = x3 * y2;
247 double x4y = x4 * y;
248
249 return c[0] + c[1] * x + c[2] * y + c[3] * x2 + c[4] * xy + c[5] * y2 + c[6] * x3 + c[7] * x2y + c[8] * xy2
250 + c[9] * y3 + c[10] * x4 + c[11] * x3y + c[12] * x2y2 + c[13] * xy3 + c[14] * y4 + c[15] * x5 + c[16] * x4y
251 + c[17] * x3y2 + c[18] * x2y3 + c[19] * xy4 + c[20] * y5;
252 }
253
257 unsigned int GetNofCoefficients() const { return 21; }
258};
259
266class CbmLitPolynom6 : public CbmLitPolynom {
267 public:
271 CbmLitPolynom6() {}
272
276 virtual ~CbmLitPolynom6() {}
277
281 double Calculate(double x, double y, double c[]) const
282 {
283
284 double x2 = x * x;
285 double y2 = y * y;
286 double xy = x * y;
287
288 double x3 = x2 * x;
289 double y3 = y2 * y;
290 double xy2 = x * y2;
291 double x2y = x2 * y;
292
293 double x4 = x3 * x;
294 double y4 = y3 * y;
295 double xy3 = x * y3;
296 double x2y2 = x2 * y2;
297 double x3y = x3 * y;
298
299 double x5 = x4 * x;
300 double y5 = y4 * y;
301 double xy4 = x * y4;
302 double x2y3 = x2 * y3;
303 double x3y2 = x3 * y2;
304 double x4y = x4 * y;
305
306 double x6 = x5 * x;
307 double y6 = y5 * y;
308 double xy5 = x * y5;
309 double x2y4 = x2 * y4;
310 double x3y3 = x3 * y3;
311 double x4y2 = x4 * y2;
312 double x5y = x5 * y;
313
314 return c[0] + c[1] * x + c[2] * y + c[3] * x2 + c[4] * xy + c[5] * y2 + c[6] * x3 + c[7] * x2y + c[8] * xy2
315 + c[9] * y3 + c[10] * x4 + c[11] * x3y + c[12] * x2y2 + c[13] * xy3 + c[14] * y4 + c[15] * x5 + c[16] * x4y
316 + c[17] * x3y2 + c[18] * x2y3 + c[19] * xy4 + c[20] * y5 + c[21] * x6 + c[22] * x5y + c[23] * x4y2
317 + c[24] * x3y3 + c[25] * x2y4 + c[26] * xy5 + c[27] * y6;
318 }
319
323 unsigned int GetNofCoefficients() const { return 28; }
324};
325
332class CbmLitPolynom7 : public CbmLitPolynom {
333 public:
337 CbmLitPolynom7() {}
338
342 virtual ~CbmLitPolynom7() {}
343
347 double Calculate(double x, double y, double c[]) const
348 {
349
350 double x2 = x * x;
351 double y2 = y * y;
352 double xy = x * y;
353
354 double x3 = x2 * x;
355 double y3 = y2 * y;
356 double xy2 = x * y2;
357 double x2y = x2 * y;
358
359 double x4 = x3 * x;
360 double y4 = y3 * y;
361 double xy3 = x * y3;
362 double x2y2 = x2 * y2;
363 double x3y = x3 * y;
364
365 double x5 = x4 * x;
366 double y5 = y4 * y;
367 double xy4 = x * y4;
368 double x2y3 = x2 * y3;
369 double x3y2 = x3 * y2;
370 double x4y = x4 * y;
371
372 double x6 = x5 * x;
373 double y6 = y5 * y;
374 double xy5 = x * y5;
375 double x2y4 = x2 * y4;
376 double x3y3 = x3 * y3;
377 double x4y2 = x4 * y2;
378 double x5y = x5 * y;
379
380 double x7 = x6 * x;
381 double y7 = y6 * y;
382 double xy6 = x * y6;
383 double x2y5 = x2 * y5;
384 double x3y4 = x3 * y4;
385 double x4y3 = x4 * y3;
386 double x5y2 = x5 * y2;
387 double x6y = x6 * y;
388
389 return c[0] + c[1] * x + c[2] * y + c[3] * x2 + c[4] * xy + c[5] * y2 + c[6] * x3 + c[7] * x2y + c[8] * xy2
390 + c[9] * y3 + c[10] * x4 + c[11] * x3y + c[12] * x2y2 + c[13] * xy3 + c[14] * y4 + c[15] * x5 + c[16] * x4y
391 + c[17] * x3y2 + c[18] * x2y3 + c[19] * xy4 + c[20] * y5 + c[21] * x6 + c[22] * x5y + c[23] * x4y2
392 + c[24] * x3y3 + c[25] * x2y4 + c[26] * xy5 + c[27] * y6 + c[28] * x7 + c[29] * x6y + c[30] * x5y2
393 + c[31] * x4y3 + c[32] * x3y4 + c[33] * x2y5 + c[34] * xy6 + c[35] * y7;
394 }
395
399 unsigned int GetNofCoefficients() const { return 36; }
400};
401
408class CbmLitPolynom8 : public CbmLitPolynom {
409 public:
413 CbmLitPolynom8() {}
414
418 virtual ~CbmLitPolynom8() {}
419
423 double Calculate(double x, double y, double c[]) const
424 {
425
426 double x2 = x * x;
427 double y2 = y * y;
428 double xy = x * y;
429
430 double x3 = x2 * x;
431 double y3 = y2 * y;
432 double xy2 = x * y2;
433 double x2y = x2 * y;
434
435 double x4 = x3 * x;
436 double y4 = y3 * y;
437 double xy3 = x * y3;
438 double x2y2 = x2 * y2;
439 double x3y = x3 * y;
440
441 double x5 = x4 * x;
442 double y5 = y4 * y;
443 double xy4 = x * y4;
444 double x2y3 = x2 * y3;
445 double x3y2 = x3 * y2;
446 double x4y = x4 * y;
447
448 double x6 = x5 * x;
449 double y6 = y5 * y;
450 double xy5 = x * y5;
451 double x2y4 = x2 * y4;
452 double x3y3 = x3 * y3;
453 double x4y2 = x4 * y2;
454 double x5y = x5 * y;
455
456 double x7 = x6 * x;
457 double y7 = y6 * y;
458 double xy6 = x * y6;
459 double x2y5 = x2 * y5;
460 double x3y4 = x3 * y4;
461 double x4y3 = x4 * y3;
462 double x5y2 = x5 * y2;
463 double x6y = x6 * y;
464
465 double x8 = x7 * x;
466 double y8 = y7 * y;
467 double xy7 = x * y7;
468 double x2y6 = x2 * y6;
469 double x3y5 = x3 * y5;
470 double x4y4 = x4 * y4;
471 double x5y3 = x5 * y3;
472 double x6y2 = x6 * y2;
473 double x7y = x7 * y;
474
475 return c[0] + c[1] * x + c[2] * y + c[3] * x2 + c[4] * xy + c[5] * y2 + c[6] * x3 + c[7] * x2y + c[8] * xy2
476 + c[9] * y3 + c[10] * x4 + c[11] * x3y + c[12] * x2y2 + c[13] * xy3 + c[14] * y4 + c[15] * x5 + c[16] * x4y
477 + c[17] * x3y2 + c[18] * x2y3 + c[19] * xy4 + c[20] * y5 + c[21] * x6 + c[22] * x5y + c[23] * x4y2
478 + c[24] * x3y3 + c[25] * x2y4 + c[26] * xy5 + c[27] * y6 + c[28] * x7 + c[29] * x6y + c[30] * x5y2
479 + c[31] * x4y3 + c[32] * x3y4 + c[33] * x2y5 + c[34] * xy6 + c[35] * y7 + c[36] * x8 + c[37] * x7y
480 + c[38] * x6y2 + c[39] * x5y3 + c[40] * x4y4 + c[41] * x3y5 + c[42] * x2y6 + c[43] * xy7 + c[44] * y8;
481 }
482
486 unsigned int GetNofCoefficients() const { return 45; }
487};
488
495class CbmLitPolynom9 : public CbmLitPolynom {
496 public:
500 CbmLitPolynom9() {}
501
505 virtual ~CbmLitPolynom9() {}
506
510 double Calculate(double x, double y, double c[]) const
511 {
512
513 double x2 = x * x;
514 double y2 = y * y;
515 double xy = x * y;
516
517 double x3 = x2 * x;
518 double y3 = y2 * y;
519 double xy2 = x * y2;
520 double x2y = x2 * y;
521
522 double x4 = x3 * x;
523 double y4 = y3 * y;
524 double xy3 = x * y3;
525 double x2y2 = x2 * y2;
526 double x3y = x3 * y;
527
528 double x5 = x4 * x;
529 double y5 = y4 * y;
530 double xy4 = x * y4;
531 double x2y3 = x2 * y3;
532 double x3y2 = x3 * y2;
533 double x4y = x4 * y;
534
535 double x6 = x5 * x;
536 double y6 = y5 * y;
537 double xy5 = x * y5;
538 double x2y4 = x2 * y4;
539 double x3y3 = x3 * y3;
540 double x4y2 = x4 * y2;
541 double x5y = x5 * y;
542
543 double x7 = x6 * x;
544 double y7 = y6 * y;
545 double xy6 = x * y6;
546 double x2y5 = x2 * y5;
547 double x3y4 = x3 * y4;
548 double x4y3 = x4 * y3;
549 double x5y2 = x5 * y2;
550 double x6y = x6 * y;
551
552 double x8 = x7 * x;
553 double y8 = y7 * y;
554 double xy7 = x * y7;
555 double x2y6 = x2 * y6;
556 double x3y5 = x3 * y5;
557 double x4y4 = x4 * y4;
558 double x5y3 = x5 * y3;
559 double x6y2 = x6 * y2;
560 double x7y = x7 * y;
561
562 double x9 = x8 * x;
563 double y9 = y8 * y;
564 double xy8 = x * y8;
565 double x2y7 = x2 * y7;
566 double x3y6 = x3 * y6;
567 double x4y5 = x4 * y5;
568 double x5y4 = x5 * y4;
569 double x6y3 = x6 * y3;
570 double x7y2 = x7 * y2;
571 double x8y = x8 * y;
572
573 return c[0] + c[1] * x + c[2] * y + c[3] * x2 + c[4] * xy + c[5] * y2 + c[6] * x3 + c[7] * x2y + c[8] * xy2
574 + c[9] * y3 + c[10] * x4 + c[11] * x3y + c[12] * x2y2 + c[13] * xy3 + c[14] * y4 + c[15] * x5 + c[16] * x4y
575 + c[17] * x3y2 + c[18] * x2y3 + c[19] * xy4 + c[20] * y5 + c[21] * x6 + c[22] * x5y + c[23] * x4y2
576 + c[24] * x3y3 + c[25] * x2y4 + c[26] * xy5 + c[27] * y6 + c[28] * x7 + c[29] * x6y + c[30] * x5y2
577 + c[31] * x4y3 + c[32] * x3y4 + c[33] * x2y5 + c[34] * xy6 + c[35] * y7 + c[36] * x8 + c[37] * x7y
578 + c[38] * x6y2 + c[39] * x5y3 + c[40] * x4y4 + c[41] * x3y5 + c[42] * x2y6 + c[43] * xy7 + c[44] * y8
579 + c[45] * x9 + c[46] * x8y + c[47] * x7y2 + c[48] * x6y3 + c[49] * x5y4 + c[50] * x4y5 + c[51] * x3y6
580 + c[52] * x2y7 + c[53] * xy8 + c[54] * y9;
581 }
582
586 unsigned int GetNofCoefficients() const { return 55; }
587};
588
595class CbmLitPolynom10 : public CbmLitPolynom {
596 public:
600 CbmLitPolynom10() {}
601
605 virtual ~CbmLitPolynom10() {}
606
610 double Calculate(double x, double y, double c[]) const
611 {
612 double x2 = x * x;
613 double y2 = y * y;
614 double xy = x * y;
615
616 double x3 = x2 * x;
617 double y3 = y2 * y;
618 double xy2 = x * y2;
619 double x2y = x2 * y;
620
621 double x4 = x3 * x;
622 double y4 = y3 * y;
623 double xy3 = x * y3;
624 double x2y2 = x2 * y2;
625 double x3y = x3 * y;
626
627 double x5 = x4 * x;
628 double y5 = y4 * y;
629 double xy4 = x * y4;
630 double x2y3 = x2 * y3;
631 double x3y2 = x3 * y2;
632 double x4y = x4 * y;
633
634 double x6 = x5 * x;
635 double y6 = y5 * y;
636 double xy5 = x * y5;
637 double x2y4 = x2 * y4;
638 double x3y3 = x3 * y3;
639 double x4y2 = x4 * y2;
640 double x5y = x5 * y;
641
642 double x7 = x6 * x;
643 double y7 = y6 * y;
644 double xy6 = x * y6;
645 double x2y5 = x2 * y5;
646 double x3y4 = x3 * y4;
647 double x4y3 = x4 * y3;
648 double x5y2 = x5 * y2;
649 double x6y = x6 * y;
650
651 double x8 = x7 * x;
652 double y8 = y7 * y;
653 double xy7 = x * y7;
654 double x2y6 = x2 * y6;
655 double x3y5 = x3 * y5;
656 double x4y4 = x4 * y4;
657 double x5y3 = x5 * y3;
658 double x6y2 = x6 * y2;
659 double x7y = x7 * y;
660
661 double x9 = x8 * x;
662 double y9 = y8 * y;
663 double xy8 = x * y8;
664 double x2y7 = x2 * y7;
665 double x3y6 = x3 * y6;
666 double x4y5 = x4 * y5;
667 double x5y4 = x5 * y4;
668 double x6y3 = x6 * y3;
669 double x7y2 = x7 * y2;
670 double x8y = x8 * y;
671
672 double x10 = x9 * x;
673 double y10 = y9 * y;
674 double xy9 = x * y9;
675 double x2y8 = x2 * y8;
676 double x3y7 = x3 * y7;
677 double x4y6 = x4 * y6;
678 double x5y5 = x5 * y5;
679 double x6y4 = x6 * y4;
680 double x7y3 = x7 * y3;
681 double x8y2 = x8 * y2;
682 double x9y = x9 * y;
683
684 return c[0] + c[1] * x + c[2] * y + c[3] * x2 + c[4] * xy + c[5] * y2 + c[6] * x3 + c[7] * x2y + c[8] * xy2
685 + c[9] * y3 + c[10] * x4 + c[11] * x3y + c[12] * x2y2 + c[13] * xy3 + c[14] * y4 + c[15] * x5 + c[16] * x4y
686 + c[17] * x3y2 + c[18] * x2y3 + c[19] * xy4 + c[20] * y5 + c[21] * x6 + c[22] * x5y + c[23] * x4y2
687 + c[24] * x3y3 + c[25] * x2y4 + c[26] * xy5 + c[27] * y6 + c[28] * x7 + c[29] * x6y + c[30] * x5y2
688 + c[31] * x4y3 + c[32] * x3y4 + c[33] * x2y5 + c[34] * xy6 + c[35] * y7 + c[36] * x8 + c[37] * x7y
689 + c[38] * x6y2 + c[39] * x5y3 + c[40] * x4y4 + c[41] * x3y5 + c[42] * x2y6 + c[43] * xy7 + c[44] * y8
690 + c[45] * x9 + c[46] * x8y + c[47] * x7y2 + c[48] * x6y3 + c[49] * x5y4 + c[50] * x4y5 + c[51] * x3y6
691 + c[52] * x2y7 + c[53] * xy8 + c[54] * y9 + c[55] * x10 + c[56] * x9y + c[57] * x8y2 + c[58] * x7y3
692 + c[59] * x6y4 + c[60] * x5y5 + c[61] * x4y6 + c[62] * x3y7 + c[63] * x2y8 + c[64] * xy9 + c[65] * y10;
693 }
694
698 unsigned int GetNofCoefficients() const { return 66; }
699};
700
701
708class FCNPolynom : public ROOT::Math::IBaseFunctionMultiDim {
709 public:
713 FCNPolynom(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& z,
714 CbmLitPolynom* polynom)
715 : fX(x)
716 , fY(y)
717 , fZ(z)
718 , fPolynom(polynom)
719 {
720 }
721
725 ~FCNPolynom() {}
726
730 virtual double Up() const { return 1.; }
731
736 // virtual double operator()(
737 // const std::vector<double>& par) const
738 Double_t DoEval(const Double_t* x) const
739 {
740 // double* p = const_cast<double*>(&par[0]);
741 double* p = const_cast<double*>(x);
742 double r = 0.;
743 for (unsigned int i = 0; i < fX.size(); i++) {
744 // calculate weight
745 // double rad = std::sqrt(std::fabs(fX[i]*fX[i] + fY[i]*fY[i]));
746 // if( rad == 0. ) continue;
747 // Double_t w = 1./rad;//(r*r+1);
748 // Double_t T = 0.;
749 // Double_t w = (1 + T) / (1 + T * std::exp(rad));
750 double w = 1.;
751 double ri = w * (fZ[i] - fPolynom->Calculate(fX[i], fY[i], p));
752 r += ri * ri;
753 }
754 return r;
755 }
756
757 unsigned int NDim() const { return fPolynom->GetNofCoefficients(); }
758
759 ROOT::Math::IBaseFunctionMultiDim* Clone() const { return new FCNPolynom(fX, fY, fZ, fPolynom); }
760
765 const CbmLitPolynom* GetPolynom() const { return fPolynom; }
766
767 private:
768 std::vector<double> fX;
769 std::vector<double> fY;
770 std::vector<double> fZ;
771 CbmLitPolynom* fPolynom;
772};
773
774
775CbmLitFieldFitter::CbmLitFieldFitter(unsigned int polynomDegree)
776 : fXangle(27.)
777 , fYangle(27.)
778 , fNofBinsX(100)
779 , fNofBinsY(100)
780 , fUseEllipseAcc(true)
781 , fPolynomDegree(polynomDegree)
782{
783 fField = FairRunAna::Instance()->GetField();
784
785 switch (polynomDegree) {
786 case 0: fPolynom = new CbmLitPolynom0(); break;
787 case 1: fPolynom = new CbmLitPolynom1(); break;
788 case 2: fPolynom = new CbmLitPolynom2(); break;
789 case 3: fPolynom = new CbmLitPolynom3(); break;
790 case 4: fPolynom = new CbmLitPolynom4(); break;
791 case 5: fPolynom = new CbmLitPolynom5(); break;
792 case 6: fPolynom = new CbmLitPolynom6(); break;
793 case 7: fPolynom = new CbmLitPolynom7(); break;
794 case 8: fPolynom = new CbmLitPolynom8(); break;
795 case 9: fPolynom = new CbmLitPolynom9(); break;
796 default: fPolynom = NULL; break;
797 }
798}
799
800CbmLitFieldFitter::~CbmLitFieldFitter() {}
801
802void CbmLitFieldFitter::FitSliceScal(float Z, lit::parallel::LitFieldSlice<fscal>& slice) { FitSlice<fscal>(Z, slice); }
803
804void CbmLitFieldFitter::FitSliceVec(float Z, lit::parallel::LitFieldSlice<fvec>& slice) { FitSlice<fvec>(Z, slice); }
805
806void CbmLitFieldFitter::FitSlice(double Z, std::vector<double>& parBx, std::vector<double>& parBy,
807 std::vector<double>& parBz)
808{
809 double tanXangle = std::tan(fXangle * 3.14159265 / 180.); //
810 double tanYangle = std::tan(fYangle * 3.14159265 / 180.); //
811
812 double Xmax = Z * tanXangle;
813 double Ymax = Z * tanYangle;
814
815 std::cout << "Fit slice. Xmax=" << Xmax << " Ymax=" << Ymax << " Z=" << Z << " Xanlge=" << fXangle
816 << " Yangle=" << fYangle << std::endl;
817
818 double HX = 2 * Xmax / fNofBinsX; // step size for X position
819 double HY = 2 * Ymax / fNofBinsY; // step size for Y position
820 std::vector<double> x;
821 std::vector<double> y;
822 std::vector<double> Bx;
823 std::vector<double> By;
824 std::vector<double> Bz;
825 for (int j = 0; j < fNofBinsX; j++) { // loop over x position
826 double X = -Xmax + j * HX;
827 for (int k = 0; k < fNofBinsY; k++) { // loop over y position
828 double Y = -Ymax + k * HY;
829
830 //check the acceptance
831 if (fUseEllipseAcc) {
832 double el = (X * X) / (Xmax * Xmax) + (Y * Y) / (Ymax * Ymax);
833 if (el > 1.) {
834 continue;
835 }
836 }
837
838 // get field value
839 double pos[3] = {X, Y, Z};
840 double B[3];
841 fField->GetFieldValue(pos, B);
842
843 x.push_back(X);
844 y.push_back(Y);
845 Bx.push_back(B[0]);
846 By.push_back(B[1]);
847 Bz.push_back(B[2]);
848 }
849 }
850 FitSlice(x, y, Bx, parBx);
851 FitSlice(x, y, By, parBy);
852 FitSlice(x, y, Bz, parBz);
853}
854
855void CbmLitFieldFitter::FitSlice(const std::vector<double>& x, const std::vector<double>& y,
856 const std::vector<double>& z, std::vector<double>& par)
857{
858 FCNPolynom* theFCN = new FCNPolynom(x, y, z, fPolynom);
859
860 // TFitterMinuit theMinuit;
861 ROOT::Minuit2::Minuit2Minimizer theMinuit;
862
863 theMinuit.SetPrintLevel(-1);
864 theMinuit.SetFunction(*theFCN);
865 int nofParams = theFCN->GetPolynom()->GetNofCoefficients();
866
867 for (int i = 0; i < nofParams; i++) {
868 std::string ss = "c" + Cbm::ToString<int>(i);
869 theMinuit.SetVariable(i, ss.c_str(), 0., 0.1);
870 theMinuit.SetVariableLimits(i, 1., -1.);
871 }
872 // theMinuit.CreateMinimizer();
873 theMinuit.Minimize();
874
875 const double* fitResults = theMinuit.X();
876
877 par.clear();
878 for (int i = 0; i < nofParams; i++) {
879 par.push_back(fitResults[i]);
880 }
881 // delete theFCN;
882}
883
884
885void CbmLitFieldFitter::FitSliceMy(double Z, std::vector<double>& parBx, std::vector<double>& parBy,
886 std::vector<double>& parBz)
887{
888 const int M = fPolynomDegree; //5; // polinom order
889 const int N = (M + 1) * (M + 2) / 2;
890
891 double tanXangle = std::tan(fXangle * 3.14159265 / 180); // AL
892 double tanYangle = std::tan(fYangle * 3.14159265 / 180); // AL
893
894 double Xmax = Z * tanXangle; // AL
895 double Ymax = Z * tanYangle; // AL
896
897 double dx = 1.; // step for the field approximation
898 double dy = 1.;
899
900 if (dx > Xmax / N / 2) {
901 dx = Xmax / N / 4.;
902 }
903 if (dy > Ymax / N / 2) {
904 dy = Ymax / N / 4.;
905 }
906
907 // double C[3][N];
908 // for( int i=0; i<3; i++) {
909 // for( int k=0; k<N; k++) {
910 // C[i][k] = 0.;
911 // }
912 // }
913
914 TMatrixD A(N, N);
915 // A.Zero();
916 TVectorD b0(N), b1(N), b2(N);
917 // b0.Zero();
918 // b1.Zero();
919 // b2.Zero();
920 for (int i = 0; i < N; i++) {
921 for (int j = 0; j < N; j++) {
922 A(i, j) = 0.; // AL was ==
923 }
924 b0(i) = b1(i) = b2(i) = 0.;
925 }
926 for (double x = -Xmax; x <= Xmax; x += dx) {
927 for (double y = -Ymax; y <= Ymax; y += dy) {
928 double r = std::sqrt(std::fabs(x * x / Xmax / Xmax + y / Ymax * y / Ymax));
929 if (r > 1.) {
930 continue;
931 }
932 Double_t w = 1. / (r * r + 1);
933 Double_t p[3] = {x, y, Z};
934 Double_t B[3] = {0., 0., 0.};
935 FairRunAna::Instance()->GetField()->GetFieldValue(p, B);
936 TVectorD m(N);
937 m(0) = 1;
938 for (int i = 1; i <= M; i++) {
939 int k = (i - 1) * (i) / 2;
940 int l = i * (i + 1) / 2;
941 for (int j = 0; j < i; j++) {
942 m(l + j) = x * m(k + j);
943 }
944 m(l + i) = y * m(k + i - 1);
945 }
946
947 TVectorD mt = m;
948 for (int i = 0; i < N; i++) {
949 for (int j = 0; j < N; j++) {
950 A(i, j) += w * m(i) * m(j);
951 }
952 b0(i) += w * B[0] * m(i);
953 b1(i) += w * B[1] * m(i);
954 b2(i) += w * B[2] * m(i);
955 }
956 }
957 }
958 double det;
959 A.Invert(&det);
960 TVectorD c0 = A * b0, c1 = A * b1, c2 = A * b2;
961 // for(int i=0; i<N; i++){
962 // C[0][i] = c0(i);
963 // C[1][i] = c1(i);
964 // C[2][i] = c2(i);
965 // }
966 for (int i = 0; i < N; i++) {
967 parBx.push_back(c0(i));
968 parBy.push_back(c1(i));
969 parBz.push_back(c2(i));
970 }
971 // geo[ind++] = N;
972 // for( int k=0; k<3; k++ ){
973 // for( int j=0; j<N; j++) geo[ind++] = C[k][j];
974 // }
975}
CbmLitFieldFitter.h
Implementation of the polynomial field approximation.
CbmLitFieldFitter::FitSlice
void FitSlice(float Z, lit::parallel::LitFieldSlice< T > &slice)
Fits (X, Y) slice of the magnetic field at Z position.
Definition CbmLitFieldFitter.h:84
CbmLitFieldFitter::fPolynom
CbmLitPolynom * fPolynom
Definition CbmLitFieldFitter.h:162
CbmLitFieldFitter::FitSliceVec
void FitSliceVec(float Z, lit::parallel::LitFieldSlice< fvec > &slice)
FitSlice implementation using fvec data type.
Definition CbmLitFieldFitter.cxx:804
CbmLitFieldFitter::FitSliceMy
void FitSliceMy(double Z, std::vector< double > &parBx, std::vector< double > &parBy, std::vector< double > &parBz)
Fit (X, Y) slice of the magnetic field at Z position.
Definition CbmLitFieldFitter.cxx:885
CbmLitFieldFitter::~CbmLitFieldFitter
virtual ~CbmLitFieldFitter()
Destructor.
Definition CbmLitFieldFitter.cxx:800
CbmLitFieldFitter::FitSliceScal
void FitSliceScal(float Z, lit::parallel::LitFieldSlice< fscal > &slice)
FitSlice implementation using fscal data type.
Definition CbmLitFieldFitter.cxx:802
CbmLitFieldFitter::CbmLitFieldFitter
CbmLitFieldFitter(unsigned int polynomDegree)
Constructor.
Definition CbmLitFieldFitter.cxx:775
CbmLitFieldFitter::fPolynomDegree
unsigned int fPolynomDegree
Definition CbmLitFieldFitter.h:161
CbmLitPolynom0
0 degree polynomial with 1 coefficient.
Definition CbmLitFieldFitter.cxx:30
CbmLitPolynom0::~CbmLitPolynom0
virtual ~CbmLitPolynom0()
Destructor.
Definition CbmLitFieldFitter.cxx:40
CbmLitPolynom0::GetNofCoefficients
unsigned int GetNofCoefficients() const
Inherited from CbmLitPolynom.
Definition CbmLitFieldFitter.cxx:50
CbmLitPolynom0::Calculate
double Calculate(double x, double y, double c[]) const
Inherited from CbmLitPolynom.
Definition CbmLitFieldFitter.cxx:45
CbmLitPolynom0::CbmLitPolynom0
CbmLitPolynom0()
Constructor.
Definition CbmLitFieldFitter.cxx:35
CbmLitPolynom10
10th degree polynomial with 66 coefficients.
Definition CbmLitFieldFitter.cxx:595
CbmLitPolynom10::GetNofCoefficients
unsigned int GetNofCoefficients() const
Inherited from CbmLitPolynom.
Definition CbmLitFieldFitter.cxx:698
CbmLitPolynom10::CbmLitPolynom10
CbmLitPolynom10()
Constructor.
Definition CbmLitFieldFitter.cxx:600
CbmLitPolynom10::~CbmLitPolynom10
virtual ~CbmLitPolynom10()
Destructor.
Definition CbmLitFieldFitter.cxx:605
CbmLitPolynom10::Calculate
double Calculate(double x, double y, double c[]) const
Inherited from CbmLitPolynom.
Definition CbmLitFieldFitter.cxx:610
CbmLitPolynom1
1st degree polynomial with 3 coefficients.
Definition CbmLitFieldFitter.cxx:59
CbmLitPolynom1::GetNofCoefficients
unsigned int GetNofCoefficients() const
Inherited from CbmLitPolynom.
Definition CbmLitFieldFitter.cxx:79
CbmLitPolynom1::~CbmLitPolynom1
virtual ~CbmLitPolynom1()
Destructor.
Definition CbmLitFieldFitter.cxx:69
CbmLitPolynom1::Calculate
double Calculate(double x, double y, double c[]) const
Inherited from CbmLitPolynom.
Definition CbmLitFieldFitter.cxx:74
CbmLitPolynom1::CbmLitPolynom1
CbmLitPolynom1()
Constructor.
Definition CbmLitFieldFitter.cxx:64
CbmLitPolynom2
2nd degree polynomial with 6 coefficients.
Definition CbmLitFieldFitter.cxx:88
CbmLitPolynom2::GetNofCoefficients
unsigned int GetNofCoefficients() const
Inherited from CbmLitPolynom.
Definition CbmLitFieldFitter.cxx:111
CbmLitPolynom2::Calculate
double Calculate(double x, double y, double c[]) const
Inherited from CbmLitPolynom.
Definition CbmLitFieldFitter.cxx:103
CbmLitPolynom2::~CbmLitPolynom2
virtual ~CbmLitPolynom2()
Destructor.
Definition CbmLitFieldFitter.cxx:98
CbmLitPolynom2::CbmLitPolynom2
CbmLitPolynom2()
Constructor.
Definition CbmLitFieldFitter.cxx:93
CbmLitPolynom3
3rd degree polynomial with 10 coefficients.
Definition CbmLitFieldFitter.cxx:120
CbmLitPolynom3::GetNofCoefficients
unsigned int GetNofCoefficients() const
Inherited from CbmLitPolynom.
Definition CbmLitFieldFitter.cxx:151
CbmLitPolynom3::~CbmLitPolynom3
virtual ~CbmLitPolynom3()
Destructor.
Definition CbmLitFieldFitter.cxx:130
CbmLitPolynom3::CbmLitPolynom3
CbmLitPolynom3()
Constructor.
Definition CbmLitFieldFitter.cxx:125
CbmLitPolynom3::Calculate
double Calculate(double x, double y, double c[]) const
Inherited from CbmLitPolynom.
Definition CbmLitFieldFitter.cxx:135
CbmLitPolynom4
4th degree polynomial with 15 coefficients.
Definition CbmLitFieldFitter.cxx:160
CbmLitPolynom4::Calculate
double Calculate(double x, double y, double c[]) const
Inherited from CbmLitPolynom.
Definition CbmLitFieldFitter.cxx:175
CbmLitPolynom4::CbmLitPolynom4
CbmLitPolynom4()
Constructor.
Definition CbmLitFieldFitter.cxx:165
CbmLitPolynom4::GetNofCoefficients
unsigned int GetNofCoefficients() const
Inherited from CbmLitPolynom.
Definition CbmLitFieldFitter.cxx:200
CbmLitPolynom4::~CbmLitPolynom4
virtual ~CbmLitPolynom4()
Destructor.
Definition CbmLitFieldFitter.cxx:170
CbmLitPolynom5
5th degree polynomial with 21 coefficients.
Definition CbmLitFieldFitter.cxx:209
CbmLitPolynom5::CbmLitPolynom5
CbmLitPolynom5()
Constructor.
Definition CbmLitFieldFitter.cxx:214
CbmLitPolynom5::~CbmLitPolynom5
virtual ~CbmLitPolynom5()
Destructor.
Definition CbmLitFieldFitter.cxx:219
CbmLitPolynom5::Calculate
double Calculate(double x, double y, double c[]) const
Inherited from CbmLitPolynom.
Definition CbmLitFieldFitter.cxx:224
CbmLitPolynom5::GetNofCoefficients
unsigned int GetNofCoefficients() const
Inherited from CbmLitPolynom.
Definition CbmLitFieldFitter.cxx:257
CbmLitPolynom6
6th degree polynomial with 28 coefficients.
Definition CbmLitFieldFitter.cxx:266
CbmLitPolynom6::GetNofCoefficients
unsigned int GetNofCoefficients() const
Inherited from CbmLitPolynom.
Definition CbmLitFieldFitter.cxx:323
CbmLitPolynom6::CbmLitPolynom6
CbmLitPolynom6()
Constructor.
Definition CbmLitFieldFitter.cxx:271
CbmLitPolynom6::~CbmLitPolynom6
virtual ~CbmLitPolynom6()
Destructor.
Definition CbmLitFieldFitter.cxx:276
CbmLitPolynom6::Calculate
double Calculate(double x, double y, double c[]) const
Inherited from CbmLitPolynom.
Definition CbmLitFieldFitter.cxx:281
CbmLitPolynom7
7th degree polynomial with 36 coefficients.
Definition CbmLitFieldFitter.cxx:332
CbmLitPolynom7::~CbmLitPolynom7
virtual ~CbmLitPolynom7()
Destructor.
Definition CbmLitFieldFitter.cxx:342
CbmLitPolynom7::Calculate
double Calculate(double x, double y, double c[]) const
Inherited from CbmLitPolynom.
Definition CbmLitFieldFitter.cxx:347
CbmLitPolynom7::CbmLitPolynom7
CbmLitPolynom7()
Constructor.
Definition CbmLitFieldFitter.cxx:337
CbmLitPolynom7::GetNofCoefficients
unsigned int GetNofCoefficients() const
Inherited from CbmLitPolynom.
Definition CbmLitFieldFitter.cxx:399
CbmLitPolynom8
8th degree polynomial with 45 coefficients.
Definition CbmLitFieldFitter.cxx:408
CbmLitPolynom8::CbmLitPolynom8
CbmLitPolynom8()
Constructor.
Definition CbmLitFieldFitter.cxx:413
CbmLitPolynom8::Calculate
double Calculate(double x, double y, double c[]) const
Inherited from CbmLitPolynom.
Definition CbmLitFieldFitter.cxx:423
CbmLitPolynom8::~CbmLitPolynom8
virtual ~CbmLitPolynom8()
Destructor.
Definition CbmLitFieldFitter.cxx:418
CbmLitPolynom8::GetNofCoefficients
unsigned int GetNofCoefficients() const
Inherited from CbmLitPolynom.
Definition CbmLitFieldFitter.cxx:486
CbmLitPolynom9
9th degree polynomial with 55 coefficients.
Definition CbmLitFieldFitter.cxx:495
CbmLitPolynom9::Calculate
double Calculate(double x, double y, double c[]) const
Inherited from CbmLitPolynom.
Definition CbmLitFieldFitter.cxx:510
CbmLitPolynom9::GetNofCoefficients
unsigned int GetNofCoefficients() const
Inherited from CbmLitPolynom.
Definition CbmLitFieldFitter.cxx:586
CbmLitPolynom9::~CbmLitPolynom9
virtual ~CbmLitPolynom9()
Destructor.
Definition CbmLitFieldFitter.cxx:505
CbmLitPolynom9::CbmLitPolynom9
CbmLitPolynom9()
Constructor.
Definition CbmLitFieldFitter.cxx:500
CbmLitPolynom
Abstract class for polynomial function.
Definition CbmLitFieldFitter.h:26
CbmLitPolynom::GetNofCoefficients
virtual unsigned int GetNofCoefficients() const =0
Return number of coefficients for this polynomial function.
CbmLitPolynom::Calculate
virtual double Calculate(double x, double y, double c[]) const =0
Returns calculated value.
FCNPolynom
Implements FCNBase which is used for MINUIT minimization.
Definition CbmLitFieldFitter.cxx:708
FCNPolynom::NDim
unsigned int NDim() const
Definition CbmLitFieldFitter.cxx:757
FCNPolynom::GetPolynom
const CbmLitPolynom * GetPolynom() const
Return polynomial which is used for minimization.
Definition CbmLitFieldFitter.cxx:765
FCNPolynom::fPolynom
CbmLitPolynom * fPolynom
Definition CbmLitFieldFitter.cxx:771
FCNPolynom::DoEval
Double_t DoEval(const Double_t *x) const
Inherited from FCNBase.
Definition CbmLitFieldFitter.cxx:738
FCNPolynom::fZ
std::vector< double > fZ
Definition CbmLitFieldFitter.cxx:770
FCNPolynom::~FCNPolynom
~FCNPolynom()
Destructor.
Definition CbmLitFieldFitter.cxx:725
FCNPolynom::FCNPolynom
FCNPolynom(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &z, CbmLitPolynom *polynom)
Constructor.
Definition CbmLitFieldFitter.cxx:713
FCNPolynom::Clone
ROOT::Math::IBaseFunctionMultiDim * Clone() const
Definition CbmLitFieldFitter.cxx:759
FCNPolynom::Up
virtual double Up() const
Inherited from FCNBase.
Definition CbmLitFieldFitter.cxx:730
FCNPolynom::fX
std::vector< double > fX
Definition CbmLitFieldFitter.cxx:768
FCNPolynom::fY
std::vector< double > fY
Definition CbmLitFieldFitter.cxx:769
lit::parallel::LitFieldSlice
Approximated magnetic field slice in XY plane perpendicular to Z.
Definition LitFieldSlice.h:41
Cbm::ToString
std::string ToString(const T &value)
Definition CbmUtils.h:26