Ignore:
Timestamp:
03/18/03 20:41:39 (22 years ago)
Author:
tbretz
Message:
*** empty log message ***
File:
1 edited

Legend:

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

    r1830 r1831  
    591591
    592592    TMatrix m(n, fM.GetNcols());
     593    TVector vold(fM.GetNcols());
    593594    for (int i=0; i<n; i++)
    594     {
    595         TVector vold(fM.GetNcols());
    596595        TMatrixRow(m, i) = vold = TMatrixRow(fM, idx[i]);
    597     }
    598596
    599597    fM = m;
     
    666664    UInt_t newrow = 0;
    667665
     666    TVector vold(fM.GetNcols());
    668667    while (oldrow<rows)
    669668    {
     
    671670
    672671        if (newrow<=(unsigned int)sum)
    673         {
    674             TVector vold(fM.GetNcols());
    675672            TMatrixRow(fM, newrow++) = vold = TMatrixRow(fM, oldrow);
    676         }
    677673
    678674        oldrow++;
     
    743739    // Calculate normalization factors
    744740    //
    745     const int   nbins   = thsh.GetNbinsX();
    746     const float frombin = thsh.GetBinLowEdge(1);
    747     const float tobin   = thsh.GetBinLowEdge(nbins+1);
    748     const float dbin    = thsh.GetBinWidth(1);
    749     const Int_t nrows   = fM.GetNrows();
    750     const Int_t ncols   = fM.GetNcols();
     741    const int    nbins   = thsh.GetNbinsX();
     742    const double frombin = thsh.GetBinLowEdge(1);
     743    const double tobin   = thsh.GetBinLowEdge(nbins+1);
     744    const double dbin    = thsh.GetBinWidth(1);
     745    const Int_t  nrows   = fM.GetNrows();
     746    const Int_t  ncols   = fM.GetNcols();
    751747
    752748    //
     
    757753        hth.Fill(fM(j, refcolumn-1));
    758754
     755    hth.DrawCopy();
     756
    759757    TH1F hthd("thd", "correction factors", nbins, frombin, tobin);
    760758    hthd.Divide((TH1F*)&thsh, &hth, 1, 1);
     759
    761760    if (hthd.GetMaximum() <= 0)
    762761    {
     
    793792    TArrayF cumulweight(nrows); // keep track for each bin how many events
    794793
    795     TMatrix hindref(fM);
    796     hindref -= frombin;
    797     hindref *= 1/dbin;
     794    //
     795    // Project values in reference column into [0,1]
     796    //
     797    TVector v(fM.GetNrows());
     798    v = TMatrixColumn(fM, refcolumn-1);
     799    v += -frombin;
     800    v *= 1/dbin;
    798801
    799802    //
    800803    // select events (distribution after renormalization)
    801804    //
    802     for (Int_t ir=0; ir<nmaxevts; ir++)
    803     {
    804         const Int_t indref = (Int_t)hindref(ind[ir], refcolumn);
     805    Int_t ir;
     806    TVector vold(fM.GetNcols());
     807    for (ir=0; ir<nrows; ir++)
     808    {
     809        const Int_t indref = (Int_t)v(ind[ir]);
    805810
    806811        cumulweight[indref] += hthd.GetBinContent(indref+1);
    807812        if (cumulweight[indref]<=0.5)
    808813        {
    809             TVector vold(fM.GetNcols());
    810814            TMatrixRow(mrest, evtcount2++) = vold = TMatrixRow(fM, ind[ir]);
    811815            continue;
     
    813817
    814818        cumulweight[indref] -= 1.;
    815         if (++evtcount1 >= nmaxevts)
     819        if (++evtcount1 >=  nmaxevts)
    816820            break;
    817821
    818         TVector vold(fM.GetNcols());
    819822        TMatrixRow(mnewtmp, evtcount1) = vold = TMatrixRow(fM, ind[ir]);
    820823    }
    821824
    822     for (Int_t ir=nmaxevts; ir<nrows; ir++)
    823     {
    824         TVector vold(fM.GetNcols());
     825    for (/*empty*/; ir<nrows; ir++)
    825826        TMatrixRow(mrest, evtcount2++) = vold = TMatrixRow(fM, ind[ir]);
    826     }
    827827
    828828    //
     
    835835    fM.ResizeTo(evtcount1, ncols);
    836836    fNumRow = evtcount1;
    837     for (Int_t ir=0; ir<evtcount1; ir++)
    838     {
    839         TVector vold(fM.GetNcols());
     837    for (ir=0; ir<evtcount1; ir++)
    840838        TMatrixRow(fM, ir) = vold = TMatrixRow(mnewtmp, ir);
    841     }
    842839
    843840    if (evtcount1 < nmaxevts)
     
    848845
    849846    rest->ResizeTo(evtcount2, ncols);
    850     for (Int_t ir=0; ir<evtcount1; ir++)
     847    for (ir=0; ir<evtcount1; ir++)
    851848    {
    852849        TVector vold(fM.GetNcols());
     
    855852
    856853    return kTRUE;
    857 
    858     /*
    859      TMatrix mnew(evtcount, nconl);
    860      for (Int_t ir=0; ir<evtcount; ir++)
    861      for (Int_t ic=0; ic<fNcols; ic++)
    862      fM(ir,ic)   = mnewtmp(ir,ic);
    863 
    864      //
    865      //  test: print new matrix (part)
    866      //
    867      *fLog << "DefRefMatrix:  Event matrix (output) :" << endl;
    868      *fLog << "DefRefMatrix:  Nrows, Ncols = " << mnew.GetNrows();
    869      *fLog << " " << mnew.GetNcols() << endl;
    870 
    871      for (Int_t ir=0;ir<10; ir++)
    872      {
    873      *fLog <<ir <<" ";
    874      for (Int_t ic=0; ic<mnew.GetNcols(); ic++)
    875      cout<<Mnew(ir,ic)<<" ";
    876      *fLog <<endl;
    877      }
    878      */
    879 
    880     /*
    881      //  test print new bin contents
    882      *fLog << "MHMatrix::DefRefMatrix;  new histogram: " << endl;
    883      for (Int_t j=1; j<=fnbins; j++)
    884      {
    885      float a = fHthaft->GetBinContent(j);
    886      *fLog << j << "  "<< a << "   ";
    887      }
    888      *fLog <<endl;
    889      */
    890 
    891     /*
    892      //---------------------------------------------------------
    893      // ==== plot four histograms
    894      TCanvas *th1 = new TCanvas (fName, fName, 1);
    895      th1->Divide(2,2);
    896 
    897      th1->cd(1);
    898      ((TH1F*)fHthsh)->DrawCopy();      // target
    899 
    900      th1->cd(2);
    901      ((TH1F*)fHth)->DrawCopy();        // real histogram before
    902 
    903      th1->cd(3);
    904      ((TH1F*)fHthd)->DrawCopy();       // correction factors
    905 
    906      th1->cd(4);
    907      ((TH1F*)fHthaft)->DrawCopy();     // histogram after
    908 
    909      //---------------------------------------------------------
    910      */
    911     //return kTRUE;
    912854}
    913855
     
    977919    // select events (distribution after renormalization)
    978920    //
     921    TVector vold(fM.GetNcols());
    979922    for (Int_t ir=0; ir<nmaxevts; ir++)
    980     {
    981         TVector vold(fM.GetNcols());
    982923        TMatrixRow(mnewtmp, evtcount1++) = vold = TMatrixRow(fM, ind[ir]);
    983     }
     924
    984925    for (Int_t ir=nmaxevts; ir<nrows; ir++)
    985     {
    986         TVector vold(fM.GetNcols());
    987926        TMatrixRow(mrest, evtcount2++) = vold = TMatrixRow(fM, ind[ir]);
    988     }
    989927
    990928    //
Note: See TracChangeset for help on using the changeset viewer.