|
33 | 33 | #include <Rtypes.h> |
34 | 34 | #include <RtypesCore.h> |
35 | 35 |
|
| 36 | +#include <algorithm> |
36 | 37 | #include <array> |
37 | 38 | #include <cmath> |
38 | 39 | #include <cstddef> |
@@ -75,8 +76,9 @@ int VarManager::fgCalibrationType = 0; // 0 - no calibration, 1 - |
75 | 76 | bool VarManager::fgUseInterpolatedCalibration = true; // use interpolated calibration histograms (default: true) |
76 | 77 | int VarManager::fgEfficiencyType = 0; // type of efficiency to be applied, default is no efficiency |
77 | 78 | TObject* VarManager::fgEfficiencyHist = nullptr; // histogram for efficiency |
78 | | -TH3D* VarManager::fgPosiPhiMap = nullptr; |
79 | | -TH3D* VarManager::fgNegaPhiMap = nullptr; |
| 79 | +TObject* VarManager::fgPosiPhiMap = nullptr; |
| 80 | +TObject* VarManager::fgNegaPhiMap = nullptr; |
| 81 | +bool VarManager::fgUsePhiCorrection = false; |
80 | 82 | //__________________________________________________________________ |
81 | 83 | VarManager::VarManager() : TObject() |
82 | 84 | { |
@@ -454,127 +456,74 @@ void VarManager::FillEfficiency(float* values) |
454 | 456 | } |
455 | 457 | } |
456 | 458 |
|
457 | | -void VarManager::SetPhiMap(TH3D* hposi, TH3D* hnega) |
| 459 | +void VarManager::SetPhiMap(TObject* hposi, TObject* hnega, bool option) |
458 | 460 | { |
459 | 461 | fgPosiPhiMap = hposi; |
460 | 462 | fgNegaPhiMap = hnega; |
| 463 | + fgUsePhiCorrection = option; |
461 | 464 | } |
462 | 465 |
|
463 | 466 | double VarManager::SampleRotationPhi(double pt, double eta, int charge) |
464 | 467 | { |
465 | 468 | // each type only alarm once |
466 | | - static bool warnedInvalidCharge = false; |
467 | | - static bool warnedMissingPositiveMap = false; |
468 | | - static bool warnedMissingNegativeMap = false; |
469 | | - static bool warnedOutOfRange = false; |
470 | 469 | static bool warnedEmptyPhi = false; |
471 | 470 |
|
472 | | - TH3D* hMap = nullptr; |
473 | | - |
474 | | - if (charge > 0) { |
475 | | - hMap = fgPosiPhiMap; |
476 | | - |
477 | | - if (!hMap) { |
478 | | - if (!warnedMissingPositiveMap) { |
479 | | - LOGF(warn, |
480 | | - "Positive phi correction map is not available. " |
481 | | - "Falling back to uniform phi sampling."); |
482 | | - warnedMissingPositiveMap = true; |
483 | | - } |
| 471 | + if (!fgUsePhiCorrection) { |
| 472 | + return gRandom->Uniform(0., o2::constants::math::TwoPI); |
| 473 | + } else { |
484 | 474 |
|
485 | | - return gRandom->Uniform( |
486 | | - 0., o2::constants::math::TwoPI); |
| 475 | + TH3D* hMap = nullptr; |
| 476 | + if (charge > 0) { |
| 477 | + hMap = dynamic_cast<TH3D*>(fgPosiPhiMap); |
| 478 | + } else { |
| 479 | + hMap = dynamic_cast<TH3D*>(fgNegaPhiMap); |
487 | 480 | } |
488 | | - } else if (charge < 0) { |
489 | | - hMap = fgNegaPhiMap; |
490 | 481 |
|
491 | 482 | if (!hMap) { |
492 | | - if (!warnedMissingNegativeMap) { |
493 | | - LOGF(warn, |
494 | | - "Negative phi correction map is not available. " |
495 | | - "Falling back to uniform phi sampling."); |
496 | | - warnedMissingNegativeMap = true; |
497 | | - } |
498 | | - |
499 | | - return gRandom->Uniform( |
500 | | - 0., o2::constants::math::TwoPI); |
501 | | - } |
502 | | - } else { |
503 | | - if (!warnedInvalidCharge) { |
504 | | - LOGF(warn, |
505 | | - "Track with charge=0 passed to SampleRotationPhi. " |
506 | | - "Falling back to uniform phi sampling."); |
507 | | - warnedInvalidCharge = true; |
508 | | - } |
509 | | - |
510 | | - return gRandom->Uniform( |
511 | | - 0., o2::constants::math::TwoPI); |
512 | | - } |
513 | | - |
514 | | - // TH3 axes: X=pT, Y=phi, Z=eta |
515 | | - const int ptBin = hMap->GetXaxis()->FindBin(pt); |
516 | | - const int etaBin = hMap->GetZaxis()->FindBin(eta); |
517 | | - |
518 | | - if (ptBin < 1 || ptBin > hMap->GetNbinsX() || |
519 | | - etaBin < 1 || etaBin > hMap->GetNbinsZ()) { |
520 | | - if (!warnedOutOfRange) { |
521 | | - LOGF(warn, |
522 | | - "Track outside phi correction map: " |
523 | | - "pt=%f, eta=%f, charge=%d, ptBin=%d, etaBin=%d. " |
524 | | - "Map ranges: pt=[%f,%f], eta=[%f,%f]. " |
525 | | - "Falling back to uniform phi sampling.", |
526 | | - pt, |
527 | | - eta, |
528 | | - charge, |
529 | | - ptBin, |
530 | | - etaBin, |
531 | | - hMap->GetXaxis()->GetXmin(), |
532 | | - hMap->GetXaxis()->GetXmax(), |
533 | | - hMap->GetZaxis()->GetXmin(), |
534 | | - hMap->GetZaxis()->GetXmax()); |
535 | | - |
536 | | - warnedOutOfRange = true; |
| 483 | + LOGF(fatal, "Phi map is not a TH3D"); |
537 | 484 | } |
| 485 | + // TH3 axes: X=pT, Y=phi, Z=eta |
| 486 | + int ptBin = hMap->GetXaxis()->FindBin(pt); |
| 487 | + int etaBin = hMap->GetZaxis()->FindBin(eta); |
538 | 488 |
|
539 | | - return gRandom->Uniform( |
540 | | - 0., o2::constants::math::TwoPI); |
541 | | - } |
| 489 | + ptBin = std::clamp(ptBin, 1, hMap->GetNbinsX()); |
| 490 | + etaBin = std::clamp(etaBin, 1, hMap->GetNbinsZ()); |
542 | 491 |
|
543 | | - TH1D* hPhi = hMap->ProjectionY( |
| 492 | + TH1D* hPhi = hMap->ProjectionY( |
544 | 493 | Form("hTRPhi_tmp_charge%d_ptbin%d_etabin%d", |
545 | 494 | charge, ptBin, etaBin), |
546 | 495 | ptBin, |
547 | 496 | ptBin, |
548 | 497 | etaBin, |
549 | 498 | etaBin); |
550 | 499 |
|
551 | | - if (!hPhi || |
552 | | - hPhi->Integral(1, hPhi->GetNbinsX()) <= 0.) { |
553 | | - if (!warnedEmptyPhi) { |
554 | | - LOGF(warn, |
555 | | - "Empty phi distribution for " |
556 | | - "pt=%f, eta=%f, charge=%d, ptBin=%d, etaBin=%d. " |
557 | | - "Falling back to uniform phi sampling.", |
558 | | - pt, |
559 | | - eta, |
560 | | - charge, |
561 | | - ptBin, |
562 | | - etaBin); |
563 | | - |
564 | | - warnedEmptyPhi = true; |
565 | | - } |
| 500 | + if (!hPhi || hPhi->Integral(1, hPhi->GetNbinsX()) <= 0.) { |
| 501 | + if (!warnedEmptyPhi) { |
| 502 | + LOGF(warn, |
| 503 | + "Empty phi distribution for " |
| 504 | + "pt=%f, eta=%f, charge=%d, ptBin=%d, etaBin=%d. " |
| 505 | + "Falling back to uniform phi sampling.", |
| 506 | + pt, |
| 507 | + eta, |
| 508 | + charge, |
| 509 | + ptBin, |
| 510 | + etaBin); |
| 511 | + |
| 512 | + warnedEmptyPhi = true; |
| 513 | + } |
566 | 514 |
|
567 | | - delete hPhi; |
| 515 | + delete hPhi; |
568 | 516 |
|
569 | | - return gRandom->Uniform( |
| 517 | + return gRandom->Uniform( |
570 | 518 | 0., o2::constants::math::TwoPI); |
571 | | - } |
| 519 | + } |
572 | 520 |
|
573 | | - const double phi = |
| 521 | + const double phi = |
574 | 522 | RecoDecay::constrainAngle(hPhi->GetRandom()); |
575 | 523 |
|
576 | | - delete hPhi; |
577 | | - return phi; |
| 524 | + delete hPhi; |
| 525 | + return phi; |
| 526 | + } |
578 | 527 | } |
579 | 528 |
|
580 | 529 | //__________________________________________________________________ |
|
0 commit comments