Changeset 1828 for trunk


Ignore:
Timestamp:
03/18/03 18:15:05 (22 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mhist
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/mhist/MHMatrix.cc

    r1826 r1828  
    688688//              are not interested in this
    689689//
    690 Bool_t MHMatrix::DefRefMatrix(const TH1F &thsh, const Int_t refcolumn,
     690Bool_t MHMatrix::DefRefMatrix(const UInt_t refcolumn, const TH1F &thsh,
    691691                              const Int_t nmaxevts, TMatrix *rest)
    692692{
     
    694694    {
    695695        *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;
    696720        return kFALSE;
    697721    }
     
    705729    // the variable for which the target distribution is given
    706730    //
    707     const Int_t refcol = refcolumn<0 ? -refcolumn - 1 : refcolumn - 1;
    708731
    709732    //
     
    722745    TH1F hth("th", "data at input", nbins, frombin, tobin);
    723746    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
    729749    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);
    754751    if (hthd.GetMaximum() <= 0)
    755752    {
     
    758755    }
    759756
     757    //
     758    // ===== obtain correction factors (normalization factors)
     759    //
    760760    hthd.Scale(1/hthd.GetMaximum());
    761 
    762     //
    763     // ==== compress matrix fM into Mnew
    764     //
    765761
    766762    //
     
    781777    Int_t evtcount1 = -1;
    782778    Int_t evtcount2 =  0;
     779
    783780    TMatrix mnewtmp(nrows, ncols);
    784781    TMatrix mrest(nrows, ncols);
     
    793790    // select events (distribution after renormalization)
    794791    //
    795     for (Int_t ir=0; ir<nrows; 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);
    798795
    799796        cumulweight[indref] += hthd.GetBinContent(indref+1);
     
    807804        cumulweight[indref] -= 1.;
    808805        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             }
    815806            break;
    816         }
    817807
    818808        TVector vold(fM.GetNcols());
    819809        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]);
    820816    }
    821817
     
    827823    // this is the matrix to be used in the g/h separation
    828824    //
    829 
    830825    fM.ResizeTo(evtcount1, ncols);
    831826    fNumRow = evtcount1;
     
    907902}
    908903
     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//
     919Bool_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
    9091051// --------------------------------------------------------------------------
    9101052//
  • trunk/MagicSoft/Mars/mhist/MHMatrix.h

    r1826 r1828  
    9999    Int_t Read(const char *name);
    100100
    101     Bool_t DefRefMatrix(const TH1F &thsh, const Int_t refcolumn=-1,
     101    Bool_t DefRefMatrix(const UInt_t refcolumn, const TH1F &thsh,
    102102                        const Int_t nmaxevts=0, TMatrix *mrest=NULL);
     103    Bool_t DefRefMatrix(const Int_t nmaxevts=0, TMatrix *mrest=NULL);
    103104
    104105    ClassDef(MHMatrix, 1) // Multidimensional Matrix to store events
Note: See TracChangeset for help on using the changeset viewer.