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