Changeset 5200 for trunk/MagicSoft/Mars/mcalib
- Timestamp:
- 10/07/04 15:50:03 (20 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/mcalib/MCalibrationChargePix.cc
r5199 r5200 413 413 Float_t MCalibrationChargePix::GetConvertedRSigma() const 414 414 { 415 416 Float_t rsigma = 0.; 417 if (IsHiGainSaturation()) 418 rsigma = TMath::Sqrt(fLoGainRSigmaSquare)*fConversionHiLo; 419 else 420 rsigma = TMath::Sqrt(fHiGainRSigmaSquare); 421 422 if (rsigma < 0) 415 if (fRSigmaSquare < 0) 423 416 return -1; 424 417 425 return rsigma ; 418 const Float_t rsigma = TMath::Sqrt(fRSigmaSquare); 419 420 return IsHiGainSaturation() ? rsigma*fConversionHiLo : rsigma ; 426 421 } 427 422 … … 440 435 { 441 436 442 Float_t rsigmavar = 0.; 437 if (fRSigmaSquareVar < 0) 438 return -1; 439 443 440 // 444 441 // SigmaSquareVar = 4. * Sigma * Sigma * Var(sigma) 445 442 // ==> Var(sigma) = 0.25 * SigmaSquareVar / (Sigma * Sigma) 446 443 // 444 const Float_t rsigmaVar = 0.25 * fRSigmaSquareVar / fRSigmaSquare; 445 447 446 if (IsHiGainSaturation()) 448 rsigmavar = 0.25 * fLoGainRSigmaSquareVar / fLoGainRSigmaSquare; 449 else 450 rsigmavar = 0.25 * fHiGainRSigmaSquareVar / fHiGainRSigmaSquare; 451 452 if (rsigmavar < 0) 447 return TMath::Sqrt(rsigmaVar/fRSigmaSquare + GetConversionHiLoRelVar()) * GetRSigma(); 448 else 449 return TMath::Sqrt(rsigmaVar); 450 451 } 452 453 // -------------------------------------------------------------------------- 454 // 455 // Get the reduced Sigma: 456 // - If fRSigmaSquare is smaller than 0 (i.e. has not yet been set), return -1. 457 // 458 Float_t MCalibrationChargePix::GetRSigma() const 459 { 460 if (fRSigmaSquare < 0) 453 461 return -1; 454 462 455 if (IsHiGainSaturation()) 456 return TMath::Sqrt(rsigmavar/fLoGainRSigmaSquare + GetConversionHiLoRelVar()) * GetLoGainRSigma(); 457 else 458 return TMath::Sqrt(rsigmavar); 459 460 } 461 462 463 return TMath::Sqrt(fRSigmaSquare); 464 465 } 466 467 // -------------------------------------------------------------------------- 468 // 469 // Get the error of the reduced Sigma: 470 // - If fRSigmaSquareVar is smaller than 0 (i.e. has not yet been set), return -1. 471 // - Calculate the absolute variance of the reduced sigma with the formula: 472 // reduced sigma variance = 0.25 * fRSigmaSquareVar / fRSigmaSquare 473 // 474 Float_t MCalibrationChargePix::GetRSigmaErr() const 475 { 476 477 if (fRSigmaSquareVar < 0) 478 return -1; 479 480 // 481 // SigmaSquareVar = 4. * Sigma * Sigma * Var(sigma) 482 // ==> Var(sigma) = 0.25 * SigmaSquareVar / (Sigma * Sigma) 483 // 484 return TMath::Sqrt(0.25 * fRSigmaSquareVar / fRSigmaSquare); 485 486 } 487 488 // -------------------------------------------------------------------------- 489 // 490 // Get the reduced Sigma per Charge: 491 // - If GetRSigma() is smaller or equal 0. (i.e. has not yet been set), return -1. 492 // - If GetMean() is 0. or -1. (i.e. has not yet been set), return -1. 493 // - Return GetRSigma() / GetMean() 494 // 495 Float_t MCalibrationChargePix::GetRSigmaPerCharge() const 496 { 497 498 const Float_t rsigma = GetRSigma(); 499 500 if (rsigma <= 0) 501 return -1.; 502 503 504 const Float_t mean = GetMean(); 505 506 if (mean == 0. || mean == -1.) 507 return -1.; 508 509 return rsigma / mean; 510 } 511 512 513 // -------------------------------------------------------------------------- 514 // 515 // Get the error of the reduced Sigma per Charge: 516 // - If GetRSigmaRelVar() is smaller or equal 0. (i.e. has not yet been set), return -1. 517 // - If GetMeanRelVar() is smaller or equal 0. (i.e. has not yet been set), return -1. 518 // - Return the propagated error of GetRSigmaPerCharge() 519 // 520 Float_t MCalibrationChargePix::GetRSigmaPerChargeErr() const 521 { 522 523 const Float_t rsigmarelvar = GetRSigmaRelVar(); 524 525 if (rsigmarelvar <= 0) 526 return -1.; 527 528 529 const Float_t meanrelvar = GetMeanRelVar(); 530 531 if (meanrelvar <= 0.) 532 return -1.; 533 534 return TMath::Sqrt(rsigmarelvar + meanrelvar) * GetRSigmaPerCharge(); 535 } 536 463 537 // -------------------------------------------------------------------------- 464 538 // 465 539 // Get the reduced Sigma Square: 466 // - If f LoGainRSigmaSquare is smaller than 0 (i.e. has not yet been set), return -1.540 // - If fRSigmaSquare is smaller than 0 (i.e. has not yet been set), return -1. 467 541 // - Test bit kHiGainSaturation: 468 542 // If yes, return fRSigmaSquare, multiplied with fConversionHiLo^2, … … 471 545 Float_t MCalibrationChargePix::GetConvertedRSigmaSquare() const 472 546 { 473 474 Float_t rsigmasquare = 0.; 475 if (IsHiGainSaturation()) 476 rsigmasquare = fLoGainRSigmaSquare*fConversionHiLo*fConversionHiLo; 477 else 478 rsigmasquare = fHiGainRSigmaSquare; 479 480 if (rsigmasquare < 0) 547 if (fRSigmaSquare < 0) 481 548 return -1; 482 549 483 return rsigmasquare ;550 return IsHiGainSaturation() ? fRSigmaSquare*fConversionHiLo*fConversionHiLo : fRSigmaSquare ; 484 551 } 485 552 … … 494 561 // Else returns the relative variance of the reduced sigma 495 562 // 496 Float_t MCalibrationChargePix::Get HiGainRSigmaRelVar() const497 { 498 499 if (f HiGainRSigmaSquareVar < 0)563 Float_t MCalibrationChargePix::GetRSigmaRelVar() const 564 { 565 566 if (fRSigmaSquareVar < 0) 500 567 return -1; 501 568 … … 504 571 // ==> Var(sigma) = 0.25 * SigmaSquareVar / (Sigma * Sigma) 505 572 // 506 return 0.25 * fHiGainRSigmaSquareVar / ( fHiGainRSigmaSquare * fHiGainRSigmaSquare ); 507 508 } 509 510 // -------------------------------------------------------------------------- 511 // 512 // Get the relative variance of the reduced Sigma: 513 // - If fRSigmaSquareVar is smaller than 0 (i.e. has not yet been set), return -1. 514 // - Calculate the relative variance of the reduced sigma squares with the formula: 515 // reduced sigma rel. variance = 0.25 * fRSigmaSquareVar / fRSigmaSquare / fRSigmaSquare 516 // - Test bit kHiGainSaturation: 517 // If yes, returns the sum of the relative variances of the reduced sigma and fConversionHiLo 518 // Else returns the relative variance of the reduced sigma 519 // 520 Float_t MCalibrationChargePix::GetLoGainRSigmaRelVar() const 521 { 522 523 if (fLoGainRSigmaSquareVar < 0) 524 return -1; 525 526 // 527 // SigmaSquareVar = 4. * Sigma * Sigma * Var(sigma) 528 // ==> Var(sigma) = 0.25 * SigmaSquareVar / (Sigma * Sigma) 529 // 530 return 0.25 * fLoGainRSigmaSquareVar / ( fLoGainRSigmaSquare * fLoGainRSigmaSquare ); 531 532 } 533 534 // -------------------------------------------------------------------------- 535 // 536 // Get the reduced Sigma: 537 // - If fHiGainRSigmaSquare is smaller than 0 (i.e. has not yet been set), return -1. 538 // 539 Float_t MCalibrationChargePix::GetHiGainRSigma() const 540 { 541 if (fHiGainRSigmaSquare < 0) 542 return -1; 543 544 return TMath::Sqrt(fHiGainRSigmaSquare); 545 546 } 547 548 // -------------------------------------------------------------------------- 549 // 550 // Get the error of the reduced Sigma: 551 // - If fHiGainRSigmaSquareVar is smaller than 0 (i.e. has not yet been set), return -1. 552 // - Calculate the absolute variance of the reduced sigma with the formula: 553 // reduced sigma variance = 0.25 * fHiGainRSigmaSquareVar / fHiGainRSigmaSquare 554 // 555 Float_t MCalibrationChargePix::GetHiGainRSigmaErr() const 556 { 557 558 if (fHiGainRSigmaSquareVar < 0) 559 return -1; 560 561 // 562 // SigmaSquareVar = 4. * Sigma * Sigma * Var(sigma) 563 // ==> Var(sigma) = 0.25 * SigmaSquareVar / (Sigma * Sigma) 564 // 565 return TMath::Sqrt(0.25 * fHiGainRSigmaSquareVar / fHiGainRSigmaSquare); 566 567 } 568 569 // -------------------------------------------------------------------------- 570 // 571 // Get the reduced Sigma per Charge: 572 // - If GetHiGainRSigma() is smaller or equal 0. (i.e. has not yet been set), return -1. 573 // - If GetMean() is 0. or -1. (i.e. has not yet been set), return -1. 574 // - Return GetHiGainRSigma() / GetMean() 575 // 576 Float_t MCalibrationChargePix::GetHiGainRSigmaPerCharge() const 577 { 578 579 const Float_t rsigma = GetHiGainRSigma(); 580 581 if (rsigma <= 0) 582 return -1.; 583 584 585 const Float_t mean = GetHiGainMean(); 586 587 if (mean == 0. || mean == -1.) 588 return -1.; 589 590 return rsigma / mean; 591 } 592 593 594 // -------------------------------------------------------------------------- 595 // 596 // Get the error of the reduced Sigma per Charge: 597 // - If GetHiGainRSigmaRelVar() is smaller or equal 0. (i.e. has not yet been set), return -1. 598 // - If GetMeanRelVar() is smaller or equal 0. (i.e. has not yet been set), return -1. 599 // - Return the propagated error of GetHiGainRSigmaPerCharge() 600 // 601 Float_t MCalibrationChargePix::GetHiGainRSigmaPerChargeErr() const 602 { 603 604 const Float_t rsigmarelvar = GetHiGainRSigmaRelVar(); 605 606 if (rsigmarelvar <= 0) 607 return -1.; 608 609 const Float_t meanrelvar = GetHiGainMeanRelVar(); 610 611 if (meanrelvar <= 0.) 612 return -1.; 613 614 return TMath::Sqrt(rsigmarelvar + meanrelvar) * GetHiGainRSigmaPerCharge(); 615 } 616 617 // -------------------------------------------------------------------------- 618 // 619 // Get the reduced Sigma: 620 // - If fLoGainRSigmaSquare is smaller than 0 (i.e. has not yet been set), return -1. 621 // 622 Float_t MCalibrationChargePix::GetLoGainRSigma() const 623 { 624 if (fLoGainRSigmaSquare < 0) 625 return -1; 626 627 return TMath::Sqrt(fLoGainRSigmaSquare); 628 629 } 630 631 // -------------------------------------------------------------------------- 632 // 633 // Get the error of the reduced Sigma: 634 // - If fLoGainRSigmaSquareVar is smaller than 0 (i.e. has not yet been set), return -1. 635 // - Calculate the absolute variance of the reduced sigma with the formula: 636 // reduced sigma variance = 0.25 * fLoGainRSigmaSquareVar / fLoGainRSigmaSquare 637 // 638 Float_t MCalibrationChargePix::GetLoGainRSigmaErr() const 639 { 640 641 if (fLoGainRSigmaSquareVar < 0) 642 return -1; 643 644 // 645 // SigmaSquareVar = 4. * Sigma * Sigma * Var(sigma) 646 // ==> Var(sigma) = 0.25 * SigmaSquareVar / (Sigma * Sigma) 647 // 648 return TMath::Sqrt(0.25 * fLoGainRSigmaSquareVar / fLoGainRSigmaSquare); 649 650 } 651 652 // -------------------------------------------------------------------------- 653 // 654 // Get the reduced Sigma per Charge: 655 // - If GetLoGainRSigma() is smaller or equal 0. (i.e. has not yet been set), return -1. 656 // - If GetMean() is 0. or -1. (i.e. has not yet been set), return -1. 657 // - Return GetLoGainRSigma() / GetMean() 658 // 659 Float_t MCalibrationChargePix::GetLoGainRSigmaPerCharge() const 660 { 661 662 const Float_t rsigma = GetLoGainRSigma(); 663 664 if (rsigma <= 0) 665 return -1.; 666 667 668 const Float_t mean = GetLoGainMean(); 669 670 if (mean == 0. || mean == -1.) 671 return -1.; 672 673 return rsigma / mean; 674 } 675 676 677 // -------------------------------------------------------------------------- 678 // 679 // Get the error of the reduced Sigma per Charge: 680 // - If GetLoGainRSigmaRelVar() is smaller or equal 0. (i.e. has not yet been set), return -1. 681 // - If GetMeanRelVar() is smaller or equal 0. (i.e. has not yet been set), return -1. 682 // - Return the propagated error of GetLoGainRSigmaPerCharge() 683 // 684 Float_t MCalibrationChargePix::GetLoGainRSigmaPerChargeErr() const 685 { 686 687 const Float_t rsigmarelvar = GetLoGainRSigmaRelVar(); 688 689 if (rsigmarelvar <= 0) 690 return -1.; 691 692 const Float_t meanrelvar = GetLoGainMeanRelVar(); 693 694 if (meanrelvar <= 0.) 695 return -1.; 696 697 return TMath::Sqrt(rsigmarelvar + meanrelvar) * GetLoGainRSigmaPerCharge(); 698 } 573 return 0.25 * fRSigmaSquareVar / ( fRSigmaSquare * fRSigmaSquare ); 574 575 } 699 576 700 577 // -------------------------------------------------------------------------- … … 704 581 // - Else returns the square root of fPheFFactorMethodVar 705 582 // 706 Float_t MCalibrationChargePix::GetHiGainPheFFactorMethodErr() const 707 { 708 if (fHiGainPheFFactorMethodVar < 0.) 709 return -1.; 710 return TMath::Sqrt(fHiGainPheFFactorMethodVar); 711 } 712 713 // -------------------------------------------------------------------------- 714 // 715 // Get the error on the number of photo-electrons (F-Factor Method): 716 // - If fPheFFactorMethodVar is smaller than 0 (i.e. has not yet been set), return -1. 717 // - Else returns the square root of fPheFFactorMethodVar 718 // 719 Float_t MCalibrationChargePix::GetLoGainPheFFactorMethodErr() const 720 { 721 if (fLoGainPheFFactorMethodVar < 0.) 722 return -1.; 723 return TMath::Sqrt(fLoGainPheFFactorMethodVar); 583 Float_t MCalibrationChargePix::GetPheFFactorMethodErr() const 584 { 585 if (fPheFFactorMethodVar < 0.) 586 return -1.; 587 return TMath::Sqrt(fPheFFactorMethodVar); 724 588 } 725 589 … … 744 608 // - Else returns fPheFFactorMethodVar / fPheFFactorMethod^2 745 609 // 746 Float_t MCalibrationChargePix::GetHiGainPheFFactorMethodRelVar() const 747 { 748 if (fHiGainPheFFactorMethodVar < 0.) 749 return -1.; 750 if (fHiGainPheFFactorMethod == 0.) 751 return -1.; 752 753 return fHiGainPheFFactorMethodVar / (fHiGainPheFFactorMethod * fHiGainPheFFactorMethod); 754 } 755 756 // -------------------------------------------------------------------------- 757 // 758 // Get the relative variance on the number of photo-electrons (F-Factor Method): 759 // - If fPheFFactorMethodVar is smaller than 0 (i.e. has not yet been set), return -1. 760 // - If fPheFFactorMethod is 0, return -1. 761 // - Else returns fPheFFactorMethodVar / fPheFFactorMethod^2 762 // 763 Float_t MCalibrationChargePix::GetLoGainPheFFactorMethodRelVar() const 764 { 765 if (fLoGainPheFFactorMethodVar < 0.) 766 return -1.; 767 if (fLoGainPheFFactorMethod == 0.) 768 return -1.; 769 770 return fLoGainPheFFactorMethodVar / (fLoGainPheFFactorMethod * fLoGainPheFFactorMethod); 610 Float_t MCalibrationChargePix::GetPheFFactorMethodRelVar() const 611 { 612 if (fPheFFactorMethodVar < 0.) 613 return -1.; 614 if (fPheFFactorMethod == 0.) 615 return -1.; 616 617 return fPheFFactorMethodVar / (fPheFFactorMethod * fPheFFactorMethod); 771 618 } 772 619 … … 820 667 // in GetRSigma() and GetRSigmaErr() 821 668 // 822 Bool_t MCalibrationChargePix::Calc HiGainReducedSigma()823 { 824 825 if (Get HiGainSigma() < 0.)669 Bool_t MCalibrationChargePix::CalcReducedSigma() 670 { 671 672 if (GetSigma() < 0.) 826 673 return kFALSE; 827 674 … … 829 676 return kFALSE; 830 677 831 const Float_t sigma = fHiGainSigma ;832 const Float_t sigmavar = fHiGainSigmaVar;833 const Float_t pedRmsSquare = fPedRms*fPedRms;834 const Float_t pedRmsSquareVar = 0.25*fPedVar*pedRmsSquare;678 const Float_t sigma = IsHiGainSaturation() ? fLoGainSigma : fHiGainSigma ; 679 const Float_t sigmavar = IsHiGainSaturation() ? fLoGainSigmaVar : fHiGainSigmaVar; 680 const Float_t pedRmsSquare = IsHiGainSaturation() ? fLoGainPedRmsSquare : fPedRms*fPedRms; 681 const Float_t pedRmsSquareVar = IsHiGainSaturation() ? fLoGainPedRmsSquareVar : 0.25*fPedVar*pedRmsSquare; 835 682 836 683 if (IsDebug()) … … 851 698 // Calculate the reduced sigmas 852 699 // 853 f HiGainRSigmaSquare = sigmaSquare - pedRmsSquare;854 855 if (IsDebug()) 856 { 857 *fLog << dbginf << "ID: " << GetPixId() 858 << " Red.Sigma Square: " << f HiGainRSigmaSquare859 << endl; 860 } 861 862 if (f HiGainRSigmaSquare <= 0.)700 fRSigmaSquare = sigmaSquare - pedRmsSquare; 701 702 if (IsDebug()) 703 { 704 *fLog << dbginf << "ID: " << GetPixId() 705 << " Red.Sigma Square: " << fRSigmaSquare 706 << endl; 707 } 708 709 if (fRSigmaSquare <= 0.) 863 710 { 864 711 if (IsDebug()) … … 870 717 871 718 872 fHiGainRSigmaSquareVar = 4. * (sigmaSquareVar + pedRmsSquareVar); 873 874 if (IsDebug()) 875 { 876 *fLog << dbginf << "ID: " << GetPixId() 877 << " Var.Red.Sigma Square: " << fHiGainRSigmaSquareVar 878 << endl; 879 } 880 881 return kTRUE; 882 } 883 884 // ---------------------------------------------------------------------------- 885 // 886 // - If fSigma is smaller than 0 (i.e. has not yet been set), return kFALSE 887 // - If fPedRms is smaller than 0 (i.e. has not yet been set), return kFALSE 888 // 889 // Calculate the reduced sigma of the low-Gain FADC slices: 890 // - Test bit IsHiGainSaturation() for the Sigma: 891 // If yes, take fLoGainSigma and fLoGainSigmaVar 892 // If no , take fHiGainSigma and fHiGainSigmaVar 893 // 894 // - Test bit IsHiGainSaturation() for the pedRMS: 895 // If yes, take fLoGainPedRmsSquare and fLoGainPedRmsSquareVar 896 // If no , take fPedRms and fPedVar 897 // 898 // - Calculate the reduced sigma with the formula: 899 // fRSigmaSquare = Sigma*Sigma - pedRMS*pedRMS 900 // 901 // - If fRSigmaSquare is smaller than 0, give a warning and return kFALSE 902 // 903 // - Calculate the variance of the reduced sigma with the formula: 904 // fRSigmaSquareVar = 4.* (sigmaVar*Sigma*Sigma + pedRmsVar*pedRMS*pedRMS) 905 // 906 // A back-transformation to the corr. amplification factor of the High-Gain is done 907 // in GetRSigma() and GetRSigmaErr() 908 // 909 Bool_t MCalibrationChargePix::CalcLoGainReducedSigma() 910 { 911 912 if (GetLoGainSigma() < 0.) 913 return kFALSE; 914 915 if (fLoGainPedRmsSquare < 0.) 916 return kFALSE; 917 918 const Float_t sigma = fLoGainSigma ; 919 const Float_t sigmavar = fLoGainSigmaVar ; 920 const Float_t pedRmsSquare = fLoGainPedRmsSquare ; 921 const Float_t pedRmsSquareVar = fLoGainPedRmsSquareVar ; 922 923 if (IsDebug()) 924 { 925 *fLog << dbginf << "ID: " << GetPixId() 926 << " HiGainSaturation: " << IsHiGainSaturation() 927 << " Sigma: " << sigma 928 << " Var.Sigma: " << sigmavar 929 << " PedRmsSquare: " << pedRmsSquare 930 << " pedRmsSquareVar: " << pedRmsSquareVar 931 << endl; 932 } 933 934 const Float_t sigmaSquare = sigma * sigma; 935 const Float_t sigmaSquareVar = 4. * sigmavar * sigmaSquare; 936 937 // 938 // Calculate the reduced sigmas 939 // 940 fLoGainRSigmaSquare = sigmaSquare - pedRmsSquare; 941 942 if (IsDebug()) 943 { 944 *fLog << dbginf << "ID: " << GetPixId() 945 << " Red.Sigma Square: " << fLoGainRSigmaSquare 946 << endl; 947 } 948 949 if (fLoGainRSigmaSquare <= 0.) 950 { 951 if (IsDebug()) 952 *fLog << warn 953 << "WARNING: Cannot calculate the reduced sigma: smaller than 0 in pixel " 954 << fPixId << endl; 955 return kFALSE; 956 } 957 958 959 fLoGainRSigmaSquareVar = 4. * (sigmaSquareVar + pedRmsSquareVar); 960 961 if (IsDebug()) 962 { 963 *fLog << dbginf << "ID: " << GetPixId() 964 << " Var.Red.Sigma Square: " << fLoGainRSigmaSquareVar 719 fRSigmaSquareVar = 4. * (sigmaSquareVar + pedRmsSquareVar); 720 721 if (IsDebug()) 722 { 723 *fLog << dbginf << "ID: " << GetPixId() 724 << " Var.Red.Sigma Square: " << fRSigmaSquareVar 965 725 << endl; 966 726 } … … 995 755 // set kFFactorMethodValid to kFALSE and return kFALSE 996 756 // 997 Bool_t MCalibrationChargePix::Calc HiGainFFactor()998 { 999 1000 if (f HiGainRSigmaSquare < 0.)757 Bool_t MCalibrationChargePix::CalcFFactor() 758 { 759 760 if (fRSigmaSquare < 0.) 1001 761 return kFALSE; 1002 762 … … 1004 764 // Square all variables in order to avoid applications of square root 1005 765 // 1006 const Float_t meanSquare = Get HiGainMean() * GetHiGainMean();1007 const Float_t meanSquareRelVar = 4.* Get HiGainMeanRelVar();766 const Float_t meanSquare = GetMean() * GetMean(); 767 const Float_t meanSquareRelVar = 4.* GetMeanRelVar(); 1008 768 1009 769 const Float_t ffactorsquare = gkFFactor * gkFFactor; 1010 770 const Float_t ffactorsquareRelVar = 4.* GetFFactorRelVar(); 1011 771 1012 const Float_t rsigmaSquareRelVar = f HiGainRSigmaSquareVar / fHiGainRSigmaSquare / fHiGainRSigmaSquare;772 const Float_t rsigmaSquareRelVar = fRSigmaSquareVar / fRSigmaSquare / fRSigmaSquare; 1013 773 // 1014 774 // Calculate the number of phe's from the F-Factor method 1015 775 // (independent on Hi Gain or Lo Gain) 1016 776 // 1017 f HiGainPheFFactorMethod = ffactorsquare * meanSquare / fHiGainRSigmaSquare;777 fPheFFactorMethod = ffactorsquare * meanSquare / fRSigmaSquare; 1018 778 1019 779 if (IsDebug()) … … 1022 782 << " F-Factor Square: " << ffactorsquare 1023 783 << " Mean Square: " << meanSquare 1024 << " Red.Sigma Square: " << f HiGainRSigmaSquare1025 << " Photo-electrons: " << f HiGainPheFFactorMethod1026 << endl; 1027 } 1028 1029 if (f HiGainPheFFactorMethod < fPheFFactorMethodLimit)784 << " Red.Sigma Square: " << fRSigmaSquare 785 << " Photo-electrons: " << fPheFFactorMethod 786 << endl; 787 } 788 789 if (fPheFFactorMethod < fPheFFactorMethodLimit) 1030 790 return kFALSE; 1031 791 … … 1033 793 // Calculate the Error of Nphe 1034 794 // 1035 const Float_t pheRelVar 1036 f HiGainPheFFactorMethodVar = pheRelVar * fHiGainPheFFactorMethod * fHiGainPheFFactorMethod;795 const Float_t pheRelVar = ffactorsquareRelVar + meanSquareRelVar + rsigmaSquareRelVar; 796 fPheFFactorMethodVar = pheRelVar * fPheFFactorMethod * fPheFFactorMethod; 1037 797 1038 798 if (IsDebug()) … … 1046 806 } 1047 807 1048 if (fHiGainPheFFactorMethodVar < 0. ) 1049 return kFALSE; 1050 1051 return kTRUE; 1052 } 1053 1054 // ------------------------------------------------------------------ 1055 // 1056 // If fRSigmaSquare is smaller than 0 (i.e. has not yet been set), 1057 // return kFALSE 1058 // 1059 // Calculate the number of photo-electrons with the F-Factor method: 1060 // - Test bit IsLoGainSaturation() for the Mean Sum of FADC slices: 1061 // If yes, take fLoGainMean and fLoGainMeanVar 1062 // If no , take fHiGainMean and fHiGainMeanVar 1063 // 1064 // - Test bit IsHiGainSaturation() for the pedRMS: 1065 // If yes, take fLoGainPedRmsSquare and fLoGainPedRmsSquareVar 1066 // If no , take fPedRms and fPedVar 1067 // 1068 // - Calculate the number of photo-electrons with the formula: 1069 // fPheFFactorMethod = gkFFactor*gkFFactor * Mean * Mean / fRSigmaSquare 1070 // 1071 // - Calculate the Variance on the photo-electrons with the formula: 1072 // fPheFFactorMethodVar = ( 4. * gkFFactorErr * gkFFactorErr / ( gkFFactor * gkFFactor ) 1073 // + 4. * Mean Var. / ( Mean * Mean ) 1074 // + fRSigmaSquareVar / fRSigmaSquare 1075 // ) * fPheFFactorMethod * fPheFFactorMethod 1076 // 1077 // - If fPheFFactorMethod is less than fPheFFactorMethodLimit, 1078 // set kFFactorMethodValid to kFALSE and return kFALSE 1079 // 1080 Bool_t MCalibrationChargePix::CalcLoGainFFactor() 1081 { 1082 1083 if (fLoGainRSigmaSquare < 0.) 1084 return kFALSE; 1085 1086 // 1087 // Square all variables in order to avoid applications of square root 1088 // 1089 const Float_t meanSquare = GetLoGainMean() * GetLoGainMean(); 1090 const Float_t meanSquareRelVar = 4.* GetLoGainMeanRelVar(); 1091 1092 const Float_t ffactorsquare = gkFFactor * gkFFactor; 1093 const Float_t ffactorsquareRelVar = 4.* GetFFactorRelVar(); 1094 1095 const Float_t rsigmaSquareRelVar = fLoGainRSigmaSquareVar / fLoGainRSigmaSquare / fLoGainRSigmaSquare; 1096 // 1097 // Calculate the number of phe's from the F-Factor method 1098 // (independent on Lo Gain or Lo Gain) 1099 // 1100 fLoGainPheFFactorMethod = ffactorsquare * meanSquare / fLoGainRSigmaSquare; 1101 1102 if (IsDebug()) 1103 { 1104 *fLog << dbginf << "ID: " << GetPixId() 1105 << " F-Factor Square: " << ffactorsquare 1106 << " Mean Square: " << meanSquare 1107 << " Red.Sigma Square: " << fLoGainRSigmaSquare 1108 << " Photo-electrons: " << fLoGainPheFFactorMethod 1109 << endl; 1110 } 1111 1112 if (fLoGainPheFFactorMethod < fPheFFactorMethodLimit) 1113 return kFALSE; 1114 1115 // 1116 // Calculate the Error of Nphe 1117 // 1118 const Float_t pheRelVar = ffactorsquareRelVar + meanSquareRelVar + rsigmaSquareRelVar; 1119 fLoGainPheFFactorMethodVar = pheRelVar * fLoGainPheFFactorMethod * fLoGainPheFFactorMethod; 1120 1121 if (IsDebug()) 1122 { 1123 *fLog << dbginf << "ID: " << GetPixId() 1124 << " Rel.Var.F-Factor Square: " << ffactorsquareRelVar 1125 << " Rel.Var. Mean Square: " << meanSquareRelVar 1126 << " Rel.Var. Red.Sigma Square: " << rsigmaSquareRelVar 1127 << " Rel.Var. Photo-electrons: " << pheRelVar 1128 << endl; 1129 } 1130 1131 if (fLoGainPheFFactorMethodVar < 0. ) 808 if (fPheFFactorMethodVar < 0. ) 1132 809 return kFALSE; 1133 810 … … 1164 841 // return kTRUE 1165 842 // 1166 Bool_t MCalibrationChargePix::Calc HiGainConvFFactor()1167 { 1168 1169 if (f HiGainPheFFactorMethod <= 0.)843 Bool_t MCalibrationChargePix::CalcConvFFactor() 844 { 845 846 if (fPheFFactorMethod <= 0.) 1170 847 return kFALSE; 1171 848 1172 const Float_t convmean = Get HiGainMean();849 const Float_t convmean = GetConvertedMean(); 1173 850 1174 851 if (convmean <= 0.) 1175 852 return kFALSE; 1176 853 1177 f HiGainMeanConvFADC2Phe = fHiGainPheFFactorMethod / convmean;854 fMeanConvFADC2Phe = fPheFFactorMethod / convmean; 1178 855 1179 856 if (IsDebug()) … … 1181 858 *fLog << dbginf << "ID: " << GetPixId() 1182 859 << " Converted Mean: " << convmean 1183 << " Conversion FADC2Phe: " << f HiGainMeanConvFADC2Phe860 << " Conversion FADC2Phe: " << fMeanConvFADC2Phe 1184 861 << endl; 1185 862 } 1186 863 1187 864 const Float_t ffactorsquareRelVar = 4.* GetFFactorRelVar(); 1188 const Float_t rsigmaSquareRelVar = f HiGainRSigmaSquareVar / fHiGainRSigmaSquare / fHiGainRSigmaSquare;865 const Float_t rsigmaSquareRelVar = fRSigmaSquareVar / fRSigmaSquare / fRSigmaSquare; 1189 866 // 1190 867 // In the calculation of the number of phe's one mean square has already been used. … … 1192 869 // the errors, but have to take account of this cancellation: 1193 870 // 1194 Float_t convrelvar = ffactorsquareRelVar + Get HiGainMeanRelVar() + rsigmaSquareRelVar;1195 const Float_t limit = fConvFFactorRelVarLimit;871 Float_t convrelvar = ffactorsquareRelVar + GetMeanRelVar() + rsigmaSquareRelVar; 872 const Float_t limit = IsHiGainSaturation() ? fConvFFactorRelVarLimit * 4. : fConvFFactorRelVarLimit; 1196 873 1197 874 // … … 1199 876 // have a fRSigmaSquareVar, calculate their error directly! 1200 877 // 1201 if (f HiGainRSigmaSquareVar < 0.)1202 convrelvar = Get HiGainMeanRelVar() + GetHiGainPheFFactorMethodRelVar();878 if (fRSigmaSquareVar < 0.) 879 convrelvar = GetMeanRelVar() + GetPheFFactorMethodRelVar(); 1203 880 1204 881 if (IsDebug()) … … 1206 883 *fLog << dbginf << "ID: " << GetPixId() 1207 884 << " Rel.Var.Red.Sigma: " << rsigmaSquareRelVar 1208 << " Rel.Var.Mean: " << Get HiGainMeanRelVar()885 << " Rel.Var.Mean: " << GetMeanRelVar() 1209 886 << " Rel.Var.F-Factor: " << ffactorsquareRelVar 1210 887 << " Rel.Var.Conversion FADC2Phe: " << convrelvar … … 1219 896 } 1220 897 1221 f HiGainMeanConvFADC2PheVar = convrelvar * fHiGainMeanConvFADC2Phe * fHiGainMeanConvFADC2Phe;898 fMeanConvFADC2PheVar = convrelvar * fMeanConvFADC2Phe * fMeanConvFADC2Phe; 1222 899 1223 900 SetFFactorMethodValid(kTRUE); … … 1225 902 } 1226 903 1227 // ------------------------------------------------------------------1228 //1229 // If fPheFFactorMethod is smaller than 0 (i.e. has not yet been set),1230 // return kFALSE1231 //1232 // If GetCovertedMean() is smaller than 0 (i.e. has not yet been set),1233 // return kFALSE1234 //1235 // Calculate fMeanConvFADC2Phe with the following formula:1236 //1237 // fMeanConvFADC2Phe = fPheFFactorMethod / GetConvMean();1238 //1239 // Calculate the rel. variance of fMeanConvFADC2Phe, taking into account that1240 // in the calculation of the number of phe's one mean square has already been used.1241 // Now, dividing by another mean, one mean calcels out, one cannot directly propagate1242 // the errors, but instead havs to take into account this cancellation:1243 //1244 // convrelvar = ffactorsquareRelVar + GetMeanRelVar() + rsigmaSquareRelVar;1245 //1246 // If confrelvar is smaller than 0. or greater than fConvFFactorRelVarLimit,1247 // return kFALSE1248 //1249 // Calculate the variance of fMeanConvFADC2Phe with the formula:1250 //1251 // fMeanConvFADC2PheVar = convrelvar * fMeanConvFADC2Phe * fMeanConvFADC2Phe;1252 //1253 // Set kFFactorMethodValid to kTRUE and1254 // return kTRUE1255 //1256 Bool_t MCalibrationChargePix::CalcLoGainConvFFactor()1257 {1258 1259 if (fLoGainPheFFactorMethod <= 0.)1260 return kFALSE;1261 1262 const Float_t convmean = GetConvertedLoGainMean();1263 1264 if (convmean <= 0.)1265 return kFALSE;1266 1267 fLoGainMeanConvFADC2Phe = fLoGainPheFFactorMethod / convmean;1268 1269 if (IsDebug())1270 {1271 *fLog << dbginf << "ID: " << GetPixId()1272 << " Converted Mean: " << convmean1273 << " Conversion FADC2Phe: " << fLoGainMeanConvFADC2Phe1274 << endl;1275 }1276 1277 const Float_t ffactorsquareRelVar = 4.* GetFFactorRelVar();1278 const Float_t rsigmaSquareRelVar = fLoGainRSigmaSquareVar / fLoGainRSigmaSquare / fLoGainRSigmaSquare;1279 //1280 // In the calculation of the number of phe's one mean square has already been used.1281 // Now, we divide by another mean, so one mean calcels out, we cannot directly propagate1282 // the errors, but have to take account of this cancellation:1283 //1284 Float_t convrelvar = ffactorsquareRelVar + GetConvertedLoGainMeanRelVar() + rsigmaSquareRelVar;1285 const Float_t limit = fConvFFactorRelVarLimit;1286 1287 //1288 // Also have to take into account the pixels labelled MBadPixelsPix::kChargeSigmaNotValid which do not1289 // have a fRSigmaSquareVar, calculate their error directly!1290 //1291 if (fLoGainRSigmaSquareVar < 0.)1292 convrelvar = GetConvertedLoGainMeanRelVar() + GetLoGainPheFFactorMethodRelVar();1293 1294 if (IsDebug())1295 {1296 *fLog << dbginf << "ID: " << GetPixId()1297 << " Rel.Var.Red.Sigma: " << rsigmaSquareRelVar1298 << " Rel.Var.Mean: " << GetLoGainMeanRelVar()1299 << " Rel.Var.F-Factor: " << ffactorsquareRelVar1300 << " Rel.Var.Conversion FADC2Phe: " << convrelvar1301 << endl;1302 }1303 1304 if (convrelvar > limit || convrelvar < 0.)1305 {1306 *fLog << warn << GetDescriptor() << ": Conv. F-Factor Method Rel. Var.: "1307 << Form("%4.3f out of limits: [0,%3.2f] in pixel:%4i",convrelvar,limit,fPixId) << endl;1308 return kFALSE;1309 }1310 1311 fLoGainMeanConvFADC2PheVar = convrelvar * fLoGainMeanConvFADC2Phe * fLoGainMeanConvFADC2Phe;1312 1313 SetFFactorMethodValid(kTRUE);1314 return kTRUE;1315 }1316 1317 904 // ---------------------------------------------------------------------------------- 1318 905 // … … 1324 911 // Calculate the error of the total F-Factor 1325 912 // 1326 Bool_t MCalibrationChargePix::Calc HiGainMeanFFactor( const Float_t nphotons, const Float_t nphotonsrelvar )913 Bool_t MCalibrationChargePix::CalcMeanFFactor( const Float_t nphotons, const Float_t nphotonsrelvar ) 1327 914 { 1328 915 … … 1333 920 << " Number photons: " << nphotons 1334 921 << " Rel.Var.Number photons: " << nphotonsrelvar 1335 << " Red.Sigma Square: " << f HiGainRSigmaSquare1336 << " Mean: " << Get HiGainMean()922 << " Red.Sigma Square: " << fRSigmaSquare 923 << " Mean: " << GetMean() 1337 924 << endl; 1338 925 } … … 1351 938 } 1352 939 1353 f HiGainMeanFFactorFADC2Phot = TMath::Sqrt(fHiGainRSigmaSquare * nphotons) / GetHiGainMean() ;1354 1355 if (IsDebug()) 1356 { 1357 *fLog << dbginf << "ID: " << GetPixId() 1358 << " F-Factor FADC2Phot: " << f HiGainMeanFFactorFADC2Phot1359 << endl; 1360 } 1361 1362 if (f HiGainMeanFFactorFADC2Phot < 0.)940 fMeanFFactorFADC2Phot = TMath::Sqrt(fRSigmaSquare * nphotons) / GetMean() ; 941 942 if (IsDebug()) 943 { 944 *fLog << dbginf << "ID: " << GetPixId() 945 << " F-Factor FADC2Phot: " << fMeanFFactorFADC2Phot 946 << endl; 947 } 948 949 if (fMeanFFactorFADC2Phot < 0.) 1363 950 { 1364 951 *fLog << warn << GetDescriptor() << ": F-Factor photons to FADC counts smaller than 0." << endl; … … 1366 953 } 1367 954 1368 const Float_t ffactorrelvar = 0.25 * f HiGainRSigmaSquareVar / ( fHiGainRSigmaSquare * fHiGainRSigmaSquare)1369 + Get HiGainMeanRelVar()955 const Float_t ffactorrelvar = 0.25 * fRSigmaSquareVar / ( fRSigmaSquare * fRSigmaSquare) 956 + GetMeanRelVar() 1370 957 + 0.25 * nphotonsrelvar; 1371 958 1372 fHiGainMeanFFactorFADC2PhotVar = ffactorrelvar * fHiGainMeanFFactorFADC2Phot * fHiGainMeanFFactorFADC2Phot; 1373 1374 if (IsDebug()) 1375 { 1376 *fLog << dbginf << "ID: " << GetPixId() 1377 << " Rel.Var.Red.Sigma: " << 0.25 * fHiGainRSigmaSquareVar / ( fHiGainRSigmaSquare * fHiGainRSigmaSquare) 1378 << " Rel.Var.Mean: " << GetHiGainMeanRelVar() 1379 << " Rel.Var.photons: " << 0.25 * nphotonsrelvar 1380 << " Rel.Var.F-Factor FADC2Phot: " << ffactorrelvar 1381 << endl; 1382 } 1383 1384 return kTRUE; 1385 } 1386 1387 // ---------------------------------------------------------------------------------- 1388 // 1389 // If photflux is smaller or equal 0, return kFALSE 1390 // 1391 // Calculate the total F-Factor with the formula: 1392 // fMeanFFactorFADC2Phot = Sqrt ( fRSigmaSquare ) / GetMean() * sqrt(nphotons) 1393 // 1394 // Calculate the error of the total F-Factor 1395 // 1396 Bool_t MCalibrationChargePix::CalcLoGainMeanFFactor( const Float_t nphotons, const Float_t nphotonsrelvar ) 1397 { 1398 1399 1400 if (IsDebug()) 1401 { 1402 *fLog << dbginf << "ID: " << GetPixId() 1403 << " Number photons: " << nphotons 1404 << " Rel.Var.Number photons: " << nphotonsrelvar 1405 << " Red.Sigma Square: " << fLoGainRSigmaSquare 1406 << " Mean: " << GetLoGainMean() 1407 << endl; 1408 } 1409 1410 1411 if (nphotons <= 0.) 1412 { 1413 *fLog << warn << GetDescriptor() << ": Assumed photon flux is smaller or equal 0." << endl; 1414 return kFALSE; 1415 } 1416 1417 if (nphotonsrelvar < 0.) 1418 { 1419 *fLog << warn << GetDescriptor() << ": Assumed photon flux variance is smaller than 0." << endl; 1420 return kFALSE; 1421 } 1422 1423 fLoGainMeanFFactorFADC2Phot = TMath::Sqrt(fLoGainRSigmaSquare * nphotons) / GetLoGainMean() ; 1424 1425 if (IsDebug()) 1426 { 1427 *fLog << dbginf << "ID: " << GetPixId() 1428 << " F-Factor FADC2Phot: " << fLoGainMeanFFactorFADC2Phot 1429 << endl; 1430 } 1431 1432 if (fLoGainMeanFFactorFADC2Phot < 0.) 1433 { 1434 *fLog << warn << GetDescriptor() << ": F-Factor photons to FADC counts smaller than 0." << endl; 1435 return kFALSE; 1436 } 1437 1438 const Float_t ffactorrelvar = 0.25 * fLoGainRSigmaSquareVar / ( fLoGainRSigmaSquare * fLoGainRSigmaSquare) 1439 + GetLoGainMeanRelVar() 1440 + 0.25 * nphotonsrelvar; 1441 1442 fLoGainMeanFFactorFADC2PhotVar = ffactorrelvar * fLoGainMeanFFactorFADC2Phot * fLoGainMeanFFactorFADC2Phot; 1443 1444 if (IsDebug()) 1445 { 1446 *fLog << dbginf << "ID: " << GetPixId() 1447 << " Rel.Var.Red.Sigma: " << 0.25 * fLoGainRSigmaSquareVar / ( fLoGainRSigmaSquare * fLoGainRSigmaSquare) 1448 << " Rel.Var.Mean: " << GetLoGainMeanRelVar() 959 fMeanFFactorFADC2PhotVar = ffactorrelvar * fMeanFFactorFADC2Phot * fMeanFFactorFADC2Phot; 960 961 if (IsDebug()) 962 { 963 *fLog << dbginf << "ID: " << GetPixId() 964 << " Rel.Var.Red.Sigma: " << 0.25 * fRSigmaSquareVar / ( fRSigmaSquare * fRSigmaSquare) 965 << " Rel.Var.Mean: " << GetMeanRelVar() 1449 966 << " Rel.Var.photons: " << 0.25 * nphotonsrelvar 1450 967 << " Rel.Var.F-Factor FADC2Phot: " << ffactorrelvar
Note:
See TracChangeset
for help on using the changeset viewer.