Changeset 5199 for trunk/MagicSoft/Mars
- Timestamp:
- 10/07/04 15:34:32 (20 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/mcalib/MCalibrationChargePix.cc
r5198 r5199 413 413 Float_t MCalibrationChargePix::GetConvertedRSigma() const 414 414 { 415 if (fRSigmaSquare < 0) 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) 416 423 return -1; 417 424 418 const Float_t rsigma = TMath::Sqrt(fRSigmaSquare); 419 420 return IsHiGainSaturation() ? rsigma*fConversionHiLo : rsigma ; 425 return rsigma ; 421 426 } 422 427 … … 435 440 { 436 441 437 if (fRSigmaSquareVar < 0) 438 return -1; 439 442 Float_t rsigmavar = 0.; 440 443 // 441 444 // SigmaSquareVar = 4. * Sigma * Sigma * Var(sigma) 442 445 // ==> Var(sigma) = 0.25 * SigmaSquareVar / (Sigma * Sigma) 443 446 // 444 const Float_t rsigmaVar = 0.25 * fRSigmaSquareVar / fRSigmaSquare;445 446 447 if (IsHiGainSaturation()) 447 return TMath::Sqrt(rsigmaVar/fRSigmaSquare + GetConversionHiLoRelVar()) * GetRSigma(); 448 rsigmavar = 0.25 * fLoGainRSigmaSquareVar / fLoGainRSigmaSquare; 449 else 450 rsigmavar = 0.25 * fHiGainRSigmaSquareVar / fHiGainRSigmaSquare; 451 452 if (rsigmavar < 0) 453 return -1; 454 455 if (IsHiGainSaturation()) 456 return TMath::Sqrt(rsigmavar/fLoGainRSigmaSquare + GetConversionHiLoRelVar()) * GetLoGainRSigma(); 448 457 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) 461 return -1; 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 458 return TMath::Sqrt(rsigmavar); 459 460 } 461 462 537 463 // -------------------------------------------------------------------------- 538 464 // 539 465 // Get the reduced Sigma Square: 540 // - If f RSigmaSquare is smaller than 0 (i.e. has not yet been set), return -1.466 // - If fLoGainRSigmaSquare is smaller than 0 (i.e. has not yet been set), return -1. 541 467 // - Test bit kHiGainSaturation: 542 468 // If yes, return fRSigmaSquare, multiplied with fConversionHiLo^2, … … 545 471 Float_t MCalibrationChargePix::GetConvertedRSigmaSquare() const 546 472 { 547 if (fRSigmaSquare < 0) 473 474 Float_t rsigmasquare = 0.; 475 if (IsHiGainSaturation()) 476 rsigmasquare = fLoGainRSigmaSquare*fConversionHiLo*fConversionHiLo; 477 else 478 rsigmasquare = fHiGainRSigmaSquare; 479 480 if (rsigmasquare < 0) 548 481 return -1; 549 482 550 return IsHiGainSaturation() ? fRSigmaSquare*fConversionHiLo*fConversionHiLo : fRSigmaSquare ;483 return rsigmasquare ; 551 484 } 552 485 … … 561 494 // Else returns the relative variance of the reduced sigma 562 495 // 563 Float_t MCalibrationChargePix::Get RSigmaRelVar() const564 { 565 566 if (f RSigmaSquareVar < 0)496 Float_t MCalibrationChargePix::GetHiGainRSigmaRelVar() const 497 { 498 499 if (fHiGainRSigmaSquareVar < 0) 567 500 return -1; 568 501 … … 571 504 // ==> Var(sigma) = 0.25 * SigmaSquareVar / (Sigma * Sigma) 572 505 // 573 return 0.25 * fRSigmaSquareVar / ( fRSigmaSquare * fRSigmaSquare ); 574 575 } 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 } 576 699 577 700 // -------------------------------------------------------------------------- … … 581 704 // - Else returns the square root of fPheFFactorMethodVar 582 705 // 583 Float_t MCalibrationChargePix::GetPheFFactorMethodErr() const 584 { 585 if (fPheFFactorMethodVar < 0.) 586 return -1.; 587 return TMath::Sqrt(fPheFFactorMethodVar); 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); 588 724 } 589 725 … … 608 744 // - Else returns fPheFFactorMethodVar / fPheFFactorMethod^2 609 745 // 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); 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); 618 771 } 619 772 … … 667 820 // in GetRSigma() and GetRSigmaErr() 668 821 // 669 Bool_t MCalibrationChargePix::Calc ReducedSigma()670 { 671 672 if (Get Sigma() < 0.)822 Bool_t MCalibrationChargePix::CalcHiGainReducedSigma() 823 { 824 825 if (GetHiGainSigma() < 0.) 673 826 return kFALSE; 674 827 … … 676 829 return kFALSE; 677 830 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;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; 682 835 683 836 if (IsDebug()) … … 698 851 // Calculate the reduced sigmas 699 852 // 700 f RSigmaSquare = sigmaSquare - pedRmsSquare;701 702 if (IsDebug()) 703 { 704 *fLog << dbginf << "ID: " << GetPixId() 705 << " Red.Sigma Square: " << f RSigmaSquare706 << endl; 707 } 708 709 if (f RSigmaSquare <= 0.)853 fHiGainRSigmaSquare = sigmaSquare - pedRmsSquare; 854 855 if (IsDebug()) 856 { 857 *fLog << dbginf << "ID: " << GetPixId() 858 << " Red.Sigma Square: " << fHiGainRSigmaSquare 859 << endl; 860 } 861 862 if (fHiGainRSigmaSquare <= 0.) 710 863 { 711 864 if (IsDebug()) … … 717 870 718 871 719 fRSigmaSquareVar = 4. * (sigmaSquareVar + pedRmsSquareVar); 720 721 if (IsDebug()) 722 { 723 *fLog << dbginf << "ID: " << GetPixId() 724 << " Var.Red.Sigma Square: " << fRSigmaSquareVar 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 725 965 << endl; 726 966 } … … 755 995 // set kFFactorMethodValid to kFALSE and return kFALSE 756 996 // 757 Bool_t MCalibrationChargePix::Calc FFactor()758 { 759 760 if (f RSigmaSquare < 0.)997 Bool_t MCalibrationChargePix::CalcHiGainFFactor() 998 { 999 1000 if (fHiGainRSigmaSquare < 0.) 761 1001 return kFALSE; 762 1002 … … 764 1004 // Square all variables in order to avoid applications of square root 765 1005 // 766 const Float_t meanSquare = Get Mean() * GetMean();767 const Float_t meanSquareRelVar = 4.* Get MeanRelVar();1006 const Float_t meanSquare = GetHiGainMean() * GetHiGainMean(); 1007 const Float_t meanSquareRelVar = 4.* GetHiGainMeanRelVar(); 768 1008 769 1009 const Float_t ffactorsquare = gkFFactor * gkFFactor; 770 1010 const Float_t ffactorsquareRelVar = 4.* GetFFactorRelVar(); 771 1011 772 const Float_t rsigmaSquareRelVar = f RSigmaSquareVar / fRSigmaSquare / fRSigmaSquare;1012 const Float_t rsigmaSquareRelVar = fHiGainRSigmaSquareVar / fHiGainRSigmaSquare / fHiGainRSigmaSquare; 773 1013 // 774 1014 // Calculate the number of phe's from the F-Factor method 775 1015 // (independent on Hi Gain or Lo Gain) 776 1016 // 777 f PheFFactorMethod = ffactorsquare * meanSquare / fRSigmaSquare;1017 fHiGainPheFFactorMethod = ffactorsquare * meanSquare / fHiGainRSigmaSquare; 778 1018 779 1019 if (IsDebug()) … … 782 1022 << " F-Factor Square: " << ffactorsquare 783 1023 << " Mean Square: " << meanSquare 784 << " Red.Sigma Square: " << f RSigmaSquare785 << " Photo-electrons: " << f PheFFactorMethod786 << endl; 787 } 788 789 if (f PheFFactorMethod < fPheFFactorMethodLimit)1024 << " Red.Sigma Square: " << fHiGainRSigmaSquare 1025 << " Photo-electrons: " << fHiGainPheFFactorMethod 1026 << endl; 1027 } 1028 1029 if (fHiGainPheFFactorMethod < fPheFFactorMethodLimit) 790 1030 return kFALSE; 791 1031 … … 793 1033 // Calculate the Error of Nphe 794 1034 // 795 const Float_t pheRelVar = ffactorsquareRelVar + meanSquareRelVar + rsigmaSquareRelVar;796 f PheFFactorMethodVar = pheRelVar * fPheFFactorMethod * fPheFFactorMethod;1035 const Float_t pheRelVar = ffactorsquareRelVar + meanSquareRelVar + rsigmaSquareRelVar; 1036 fHiGainPheFFactorMethodVar = pheRelVar * fHiGainPheFFactorMethod * fHiGainPheFFactorMethod; 797 1037 798 1038 if (IsDebug()) … … 806 1046 } 807 1047 808 if (fPheFFactorMethodVar < 0. ) 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. ) 809 1132 return kFALSE; 810 1133 … … 841 1164 // return kTRUE 842 1165 // 843 Bool_t MCalibrationChargePix::Calc ConvFFactor()844 { 845 846 if (f PheFFactorMethod <= 0.)1166 Bool_t MCalibrationChargePix::CalcHiGainConvFFactor() 1167 { 1168 1169 if (fHiGainPheFFactorMethod <= 0.) 847 1170 return kFALSE; 848 1171 849 const Float_t convmean = Get ConvertedMean();1172 const Float_t convmean = GetHiGainMean(); 850 1173 851 1174 if (convmean <= 0.) 852 1175 return kFALSE; 853 1176 854 f MeanConvFADC2Phe = fPheFFactorMethod / convmean;1177 fHiGainMeanConvFADC2Phe = fHiGainPheFFactorMethod / convmean; 855 1178 856 1179 if (IsDebug()) … … 858 1181 *fLog << dbginf << "ID: " << GetPixId() 859 1182 << " Converted Mean: " << convmean 860 << " Conversion FADC2Phe: " << f MeanConvFADC2Phe1183 << " Conversion FADC2Phe: " << fHiGainMeanConvFADC2Phe 861 1184 << endl; 862 1185 } 863 1186 864 1187 const Float_t ffactorsquareRelVar = 4.* GetFFactorRelVar(); 865 const Float_t rsigmaSquareRelVar = f RSigmaSquareVar / fRSigmaSquare / fRSigmaSquare;1188 const Float_t rsigmaSquareRelVar = fHiGainRSigmaSquareVar / fHiGainRSigmaSquare / fHiGainRSigmaSquare; 866 1189 // 867 1190 // In the calculation of the number of phe's one mean square has already been used. … … 869 1192 // the errors, but have to take account of this cancellation: 870 1193 // 871 Float_t convrelvar = ffactorsquareRelVar + Get MeanRelVar() + rsigmaSquareRelVar;872 const Float_t limit = IsHiGainSaturation() ? fConvFFactorRelVarLimit * 4. :fConvFFactorRelVarLimit;1194 Float_t convrelvar = ffactorsquareRelVar + GetHiGainMeanRelVar() + rsigmaSquareRelVar; 1195 const Float_t limit = fConvFFactorRelVarLimit; 873 1196 874 1197 // … … 876 1199 // have a fRSigmaSquareVar, calculate their error directly! 877 1200 // 878 if (f RSigmaSquareVar < 0.)879 convrelvar = Get MeanRelVar() + GetPheFFactorMethodRelVar();1201 if (fHiGainRSigmaSquareVar < 0.) 1202 convrelvar = GetHiGainMeanRelVar() + GetHiGainPheFFactorMethodRelVar(); 880 1203 881 1204 if (IsDebug()) … … 883 1206 *fLog << dbginf << "ID: " << GetPixId() 884 1207 << " Rel.Var.Red.Sigma: " << rsigmaSquareRelVar 885 << " Rel.Var.Mean: " << Get MeanRelVar()1208 << " Rel.Var.Mean: " << GetHiGainMeanRelVar() 886 1209 << " Rel.Var.F-Factor: " << ffactorsquareRelVar 887 1210 << " Rel.Var.Conversion FADC2Phe: " << convrelvar … … 896 1219 } 897 1220 898 f MeanConvFADC2PheVar = convrelvar * fMeanConvFADC2Phe * fMeanConvFADC2Phe;1221 fHiGainMeanConvFADC2PheVar = convrelvar * fHiGainMeanConvFADC2Phe * fHiGainMeanConvFADC2Phe; 899 1222 900 1223 SetFFactorMethodValid(kTRUE); … … 902 1225 } 903 1226 1227 // ------------------------------------------------------------------ 1228 // 1229 // If fPheFFactorMethod is smaller than 0 (i.e. has not yet been set), 1230 // return kFALSE 1231 // 1232 // If GetCovertedMean() is smaller than 0 (i.e. has not yet been set), 1233 // return kFALSE 1234 // 1235 // Calculate fMeanConvFADC2Phe with the following formula: 1236 // 1237 // fMeanConvFADC2Phe = fPheFFactorMethod / GetConvMean(); 1238 // 1239 // Calculate the rel. variance of fMeanConvFADC2Phe, taking into account that 1240 // 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 propagate 1242 // 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 kFALSE 1248 // 1249 // Calculate the variance of fMeanConvFADC2Phe with the formula: 1250 // 1251 // fMeanConvFADC2PheVar = convrelvar * fMeanConvFADC2Phe * fMeanConvFADC2Phe; 1252 // 1253 // Set kFFactorMethodValid to kTRUE and 1254 // return kTRUE 1255 // 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: " << convmean 1273 << " Conversion FADC2Phe: " << fLoGainMeanConvFADC2Phe 1274 << 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 propagate 1282 // 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 not 1289 // 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: " << rsigmaSquareRelVar 1298 << " Rel.Var.Mean: " << GetLoGainMeanRelVar() 1299 << " Rel.Var.F-Factor: " << ffactorsquareRelVar 1300 << " Rel.Var.Conversion FADC2Phe: " << convrelvar 1301 << 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 904 1317 // ---------------------------------------------------------------------------------- 905 1318 // … … 911 1324 // Calculate the error of the total F-Factor 912 1325 // 913 Bool_t MCalibrationChargePix::Calc MeanFFactor( const Float_t nphotons, const Float_t nphotonsrelvar )1326 Bool_t MCalibrationChargePix::CalcHiGainMeanFFactor( const Float_t nphotons, const Float_t nphotonsrelvar ) 914 1327 { 915 1328 … … 920 1333 << " Number photons: " << nphotons 921 1334 << " Rel.Var.Number photons: " << nphotonsrelvar 922 << " Red.Sigma Square: " << f RSigmaSquare923 << " Mean: " << Get Mean()1335 << " Red.Sigma Square: " << fHiGainRSigmaSquare 1336 << " Mean: " << GetHiGainMean() 924 1337 << endl; 925 1338 } … … 938 1351 } 939 1352 940 f MeanFFactorFADC2Phot = TMath::Sqrt(fRSigmaSquare * nphotons) / GetMean() ;941 942 if (IsDebug()) 943 { 944 *fLog << dbginf << "ID: " << GetPixId() 945 << " F-Factor FADC2Phot: " << f MeanFFactorFADC2Phot946 << endl; 947 } 948 949 if (f MeanFFactorFADC2Phot < 0.)1353 fHiGainMeanFFactorFADC2Phot = TMath::Sqrt(fHiGainRSigmaSquare * nphotons) / GetHiGainMean() ; 1354 1355 if (IsDebug()) 1356 { 1357 *fLog << dbginf << "ID: " << GetPixId() 1358 << " F-Factor FADC2Phot: " << fHiGainMeanFFactorFADC2Phot 1359 << endl; 1360 } 1361 1362 if (fHiGainMeanFFactorFADC2Phot < 0.) 950 1363 { 951 1364 *fLog << warn << GetDescriptor() << ": F-Factor photons to FADC counts smaller than 0." << endl; … … 953 1366 } 954 1367 955 const Float_t ffactorrelvar = 0.25 * f RSigmaSquareVar / ( fRSigmaSquare * fRSigmaSquare)956 + Get MeanRelVar()1368 const Float_t ffactorrelvar = 0.25 * fHiGainRSigmaSquareVar / ( fHiGainRSigmaSquare * fHiGainRSigmaSquare) 1369 + GetHiGainMeanRelVar() 957 1370 + 0.25 * nphotonsrelvar; 958 1371 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() 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() 966 1449 << " Rel.Var.photons: " << 0.25 * nphotonsrelvar 967 1450 << " Rel.Var.F-Factor FADC2Phot: " << ffactorrelvar
Note:
See TracChangeset
for help on using the changeset viewer.