- Timestamp:
- 03/18/03 18:15:05 (22 years ago)
- Location:
- trunk/MagicSoft/Mars/mhist
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/mhist/MHMatrix.cc
r1826 r1828 688 688 // are not interested in this 689 689 // 690 Bool_t MHMatrix::DefRefMatrix(const TH1F &thsh, const Int_t refcolumn,690 Bool_t MHMatrix::DefRefMatrix(const UInt_t refcolumn, const TH1F &thsh, 691 691 const Int_t nmaxevts, TMatrix *rest) 692 692 { … … 694 694 { 695 695 *fLog << err << dbginf << "Matrix not initialized" << endl; 696 return kFALSE; 697 } 698 699 if (refcolumn==0) 700 { 701 *fLog << err << dbginf << "Reference column 0 unknown." << endl; 702 return kFALSE; 703 } 704 705 if (thsh.GetMinimum()<0) 706 { 707 *fLog << err << dbginf << "Renormalization not possible: Target Distribution has values < 0" << endl; 708 return kFALSE; 709 } 710 711 if (nmaxevts>fM.GetNrows()) 712 { 713 *fLog << err << dbginf << "Number of maximum events exceeds number of events" << endl; 714 return kFALSE; 715 } 716 717 if (nmaxevts<0) 718 { 719 *fLog << err << dbginf << "Number of maximum events < 0" << endl; 696 720 return kFALSE; 697 721 } … … 705 729 // the variable for which the target distribution is given 706 730 // 707 const Int_t refcol = refcolumn<0 ? -refcolumn - 1 : refcolumn - 1;708 731 709 732 // … … 722 745 TH1F hth("th", "data at input", nbins, frombin, tobin); 723 746 for (Int_t j=0; j<nrows; j++) 724 hth.Fill(fM(j, refcol)); 725 726 // 727 // ===== obtain correction factors (normalization factors) 728 // 747 hth.Fill(fM(j, refcolumn-1)); 748 729 749 TH1F hthd("thd", "correction factors", nbins, frombin, tobin); 730 731 Double_t cmax = 0.; 732 for (Int_t j=1; j<=nbins; j++) 733 { 734 const float a = hth.GetBinContent(j); 735 736 // if refcolumn < 0 set the correction factors equal to 1 737 const float b = refcolumn>=0 ? thsh.GetBinContent(j) : a; 738 739 if (a<0 || b<0) 740 { 741 *fLog << err << dbginf << "renormalization of input matrix is not possible" << endl; 742 *fLog << " a, b = " << a << ", " << b << endl; 743 return kFALSE; 744 } 745 746 const float c = a==0 ? 0 : b/a; 747 748 hthd.SetBinContent(j, c); 749 750 if (cmax < c) 751 cmax = c; 752 } 753 750 hthd.Divide((TH1F*)&thsh, &hth, 1, 1); 754 751 if (hthd.GetMaximum() <= 0) 755 752 { … … 758 755 } 759 756 757 // 758 // ===== obtain correction factors (normalization factors) 759 // 760 760 hthd.Scale(1/hthd.GetMaximum()); 761 762 //763 // ==== compress matrix fM into Mnew764 //765 761 766 762 // … … 781 777 Int_t evtcount1 = -1; 782 778 Int_t evtcount2 = 0; 779 783 780 TMatrix mnewtmp(nrows, ncols); 784 781 TMatrix mrest(nrows, ncols); … … 793 790 // select events (distribution after renormalization) 794 791 // 795 for (Int_t ir=0; ir<n rows; ir++)796 { 797 const Int_t indref = (Int_t)hindref(ind[ir], refcol );792 for (Int_t ir=0; ir<nmaxevts; ir++) 793 { 794 const Int_t indref = (Int_t)hindref(ind[ir], refcolumn); 798 795 799 796 cumulweight[indref] += hthd.GetBinContent(indref+1); … … 807 804 cumulweight[indref] -= 1.; 808 805 if (++evtcount1 >= nmaxevts) 809 {810 for (; ir<nrows; ir++)811 {812 TVector vold(fM.GetNcols());813 TMatrixRow(mrest, evtcount2++) = vold = TMatrixRow(fM, ind[ir]);814 }815 806 break; 816 }817 807 818 808 TVector vold(fM.GetNcols()); 819 809 TMatrixRow(mnewtmp, evtcount1) = vold = TMatrixRow(fM, ind[ir]); 810 } 811 812 for (Int_t ir=nmaxevts; ir<nrows; ir++) 813 { 814 TVector vold(fM.GetNcols()); 815 TMatrixRow(mrest, evtcount2++) = vold = TMatrixRow(fM, ind[ir]); 820 816 } 821 817 … … 827 823 // this is the matrix to be used in the g/h separation 828 824 // 829 830 825 fM.ResizeTo(evtcount1, ncols); 831 826 fNumRow = evtcount1; … … 907 902 } 908 903 904 // ------------------------------------------------------------------------ 905 // 906 // Define the reference matrix 907 // refcolumn number of the column (starting at 1)containing the variable, 908 // for which a target distribution may be given; 909 // if refcolumn is negative the target distribution will be set 910 // equal to the real distribution; the events in the reference 911 // matrix will then be simply a random selection of the events 912 // in the original matrix. 913 // thsh histogram containing the target distribution of the variable 914 // nmaxevts maximum number of events in the reference matrix 915 // rest a TMatrix conatining the resulting (not choosen) 916 // columns of the primary matrix. Maybe NULL if you 917 // are not interested in this 918 // 919 Bool_t MHMatrix::DefRefMatrix(const Int_t nmaxevts, TMatrix *rest) 920 { 921 if (!fM.IsValid()) 922 { 923 *fLog << err << dbginf << "Matrix not initialized" << endl; 924 return kFALSE; 925 } 926 927 const Int_t nrows = fM.GetNrows(); 928 const Int_t ncols = fM.GetNcols(); 929 930 // 931 // get random access 932 // 933 TArrayF ranx(nrows); 934 935 TRandom3 rnd(0); 936 for (Int_t i=0; i<nrows; i++) 937 ranx[i] = rnd.Rndm(i); 938 939 TArrayI ind(nrows); 940 TMath::Sort(nrows, ranx.GetArray(), ind.GetArray(), kTRUE); 941 942 // 943 // define new matrix 944 // 945 Int_t evtcount1 = 0; 946 Int_t evtcount2 = 0; 947 948 TMatrix mnewtmp(nrows, ncols); 949 TMatrix mrest(nrows, ncols); 950 951 // 952 // select events (distribution after renormalization) 953 // 954 for (Int_t ir=0; ir<nmaxevts; ir++) 955 { 956 TVector vold(fM.GetNcols()); 957 TMatrixRow(mnewtmp, evtcount1++) = vold = TMatrixRow(fM, ind[ir]); 958 } 959 for (Int_t ir=nmaxevts; ir<nrows; ir++) 960 { 961 TVector vold(fM.GetNcols()); 962 TMatrixRow(mrest, evtcount2++) = vold = TMatrixRow(fM, ind[ir]); 963 } 964 965 // 966 // reduce size 967 // 968 // matrix fM having the requested distribution 969 // and the requested number of rows; 970 // this is the matrix to be used in the g/h separation 971 // 972 fM.ResizeTo(evtcount1, ncols); 973 fNumRow = evtcount1; 974 for (Int_t ir=0; ir<evtcount1; ir++) 975 { 976 TVector vold(fM.GetNcols()); 977 TMatrixRow(fM, ir) = vold = TMatrixRow(mnewtmp, ir); 978 } 979 980 if (evtcount1 < nmaxevts) 981 *fLog << warn << "The reference sample contains less events (" << evtcount1 << ") than required (" << nmaxevts << ")" << endl; 982 983 if (!rest) 984 return kTRUE; 985 986 rest->ResizeTo(evtcount2, ncols); 987 for (Int_t ir=0; ir<evtcount1; ir++) 988 { 989 TVector vold(fM.GetNcols()); 990 TMatrixRow(*rest, ir) = vold = TMatrixRow(mrest, ir); 991 } 992 993 return kTRUE; 994 995 /* 996 TMatrix mnew(evtcount, nconl); 997 for (Int_t ir=0; ir<evtcount; ir++) 998 for (Int_t ic=0; ic<fNcols; ic++) 999 fM(ir,ic) = mnewtmp(ir,ic); 1000 1001 // 1002 // test: print new matrix (part) 1003 // 1004 *fLog << "DefRefMatrix: Event matrix (output) :" << endl; 1005 *fLog << "DefRefMatrix: Nrows, Ncols = " << mnew.GetNrows(); 1006 *fLog << " " << mnew.GetNcols() << endl; 1007 1008 for (Int_t ir=0;ir<10; ir++) 1009 { 1010 *fLog <<ir <<" "; 1011 for (Int_t ic=0; ic<mnew.GetNcols(); ic++) 1012 cout<<Mnew(ir,ic)<<" "; 1013 *fLog <<endl; 1014 } 1015 */ 1016 1017 /* 1018 // test print new bin contents 1019 *fLog << "MHMatrix::DefRefMatrix; new histogram: " << endl; 1020 for (Int_t j=1; j<=fnbins; j++) 1021 { 1022 float a = fHthaft->GetBinContent(j); 1023 *fLog << j << " "<< a << " "; 1024 } 1025 *fLog <<endl; 1026 */ 1027 1028 /* 1029 //--------------------------------------------------------- 1030 // ==== plot four histograms 1031 TCanvas *th1 = new TCanvas (fName, fName, 1); 1032 th1->Divide(2,2); 1033 1034 th1->cd(1); 1035 ((TH1F*)fHthsh)->DrawCopy(); // target 1036 1037 th1->cd(2); 1038 ((TH1F*)fHth)->DrawCopy(); // real histogram before 1039 1040 th1->cd(3); 1041 ((TH1F*)fHthd)->DrawCopy(); // correction factors 1042 1043 th1->cd(4); 1044 ((TH1F*)fHthaft)->DrawCopy(); // histogram after 1045 1046 //--------------------------------------------------------- 1047 */ 1048 //return kTRUE; 1049 } 1050 909 1051 // -------------------------------------------------------------------------- 910 1052 // -
trunk/MagicSoft/Mars/mhist/MHMatrix.h
r1826 r1828 99 99 Int_t Read(const char *name); 100 100 101 Bool_t DefRefMatrix(const TH1F &thsh, const Int_t refcolumn=-1,101 Bool_t DefRefMatrix(const UInt_t refcolumn, const TH1F &thsh, 102 102 const Int_t nmaxevts=0, TMatrix *mrest=NULL); 103 Bool_t DefRefMatrix(const Int_t nmaxevts=0, TMatrix *mrest=NULL); 103 104 104 105 ClassDef(MHMatrix, 1) // Multidimensional Matrix to store events
Note:
See TracChangeset
for help on using the changeset viewer.