#include #include #include #include #include void optPad() { Int_t N = 10; Double_t eps = 1.e-10; TMatrix g(N, N); g.Zero(); // hista is the normalized 1D histogram of sigmabar for ON data // histb is the normalized 1D histogram of sigmabar for OFF data // histc is the difference ON-OFF TH1D *hista = new TH1D("ON", "ON data before padding", N, 0., 5.); TH1D *histb = new TH1D("OFF", "OFF data before padding", N, 0., 5.); TH1D *histc = new TH1D("ON-OFF", "ON-OFF before padding", N, 0., 5.); hista->SetMaximum( 5.0/(Double_t)N ); hista->SetMinimum(-1.0/(Double_t)N ); histb->SetMaximum( 5.0/(Double_t)N ); histb->SetMinimum(-1.0/(Double_t)N ); histc->SetMaximum( 5.0/(Double_t)N ); histc->SetMinimum(-1.0/(Double_t)N ); // at the beginning, histap is a copy of hista // at the end, it will be the 1D histogram for ON data after padding // at the beginning, histbp is a copy of histb // at the end, it will be the 1D histogram for OFF data after padding // at the beginning, histcp is a copy of histc // at the end, it should be zero TH1D *histap = new TH1D("ONp", "ON data after padding", N, 0., 5.); TH1D *histbp = new TH1D("OFFp", "OFF data after padding", N, 0., 5.); TH1D *histcp = new TH1D("ON-OFFp", "ON-OFF after padding", N, 0., 5.); histap->SetMaximum( 5.0/(Double_t)N ); histap->SetMinimum(-1.0/(Double_t)N ); histbp->SetMaximum( 5.0/(Double_t)N ); histbp->SetMinimum(-1.0/(Double_t)N ); histcp->SetMaximum( 1.0/(Double_t)N); histcp->SetMinimum(-1.0/(Double_t)N); Double_t shoulda[] = {1.0, 2.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 2.0, 1.0}; Double_t shouldb[] = {0.0, 0.0, 1.0, 2.0, 3.0, 2.0, 1.0, 0.0, 1.0, 2.0}; for (Int_t i=1; i<=N; i++) { hista->SetBinContent(i, shoulda[i-1]); histb->SetBinContent(i, shouldb[i-1]); } // normalize histograms Double_t suma, sumb, conta, contb; suma = 0.0; sumb = 0.0; for (Int_t i=1; i<=N; i++) { conta = hista->GetBinContent(i); contb = histb->GetBinContent(i); suma += conta; sumb += contb; } for (Int_t i=1; i<=N; i++) { conta = hista->GetBinContent(i); contb = histb->GetBinContent(i); hista->SetBinContent(i,conta/suma); histb->SetBinContent(i,contb/sumb); histc->SetBinContent(i, conta/suma - contb/sumb); histap->SetBinContent(i,conta/suma); histbp->SetBinContent(i,contb/sumb); histcp->SetBinContent(i, conta/suma - contb/sumb); } // g[k-1][j-1] (if <0.0) tells how many ON events in bin k should be padded // from sigma_k to sigma_j // g[k-1][j-1] (if >0.0) tells how many OFF events in bin k should be padded // from sigma_k to sigma_j //-------- start j loop ------------------------------------------------ // loop over bins in histc, starting from the end Double_t v, s, w, t, x, u, a, b, arg; for (Int_t j=N; j >= 1; j--) { v = histcp->GetBinContent(j); if ( fabs(v) < eps ) continue; if (v >= 0.0) s = 1.0; else s = -1.0; //................ start k loop ...................................... // look for a bin k which may compensate the content of bin j for (Int_t k=j-1; k >= 1; k--) { w = histcp->GetBinContent(k); if ( fabs(w) < eps ) continue; if (w >= 0.0) t = 1.0; else t = -1.0; if (s==t) continue; x = v + w; if (x >= 0.0) u = 1.0; else u = -1.0; if (u == s) { arg = -w; g(k-1, j-1) = arg; cout << "k-1, j-1, arg = " << k-1 << ", " << j-1 << ", " << arg << endl; //...................................... // this is for checking the procedure if (arg < 0.0) { a = histap->GetBinContent(k); histap->SetBinContent(k, a+arg); a = histap->GetBinContent(j); histap->SetBinContent(j, a-arg); } else { b = histbp->GetBinContent(k); histbp->SetBinContent(k, b-arg); b = histbp->GetBinContent(j); histbp->SetBinContent(j, b+arg); } //...................................... histcp->SetBinContent(k, 0.0); histcp->SetBinContent(j, x); //...................................... // redefine v v = histcp->GetBinContent(j); if ( fabs(v) < eps ) break; if (v >= 0.0) s = 1.0; else s = -1.0; //...................................... continue; } arg = v; g(k-1, j-1) = arg; cout << "k-1, j-1, arg = " << k-1 << ", " << j-1 << ", " << arg << endl; //...................................... // this is for checking the procedure if (arg < 0.0) { a = histap->GetBinContent(k); histap->SetBinContent(k, a+arg); a = histap->GetBinContent(j); histap->SetBinContent(j, a-arg); } else { b = histbp->GetBinContent(k); histbp->SetBinContent(k, b-arg); b = histbp->GetBinContent(j); histbp->SetBinContent(j, b+arg); } //...................................... histcp->SetBinContent(k, x); histcp->SetBinContent(j, 0.0); break; } //................ end k loop ...................................... } //-------- end j loop ------------------------------------------------ TVector index(N); index.Zero(); TVector map(N); Int_t indexold = 0; cout << "=========================================================" << endl; for (Int_t k=0; k eps) { map(indexold) = g(k, j); cout << "k, j, indexold, map = " << k << ", " << j << ", " << indexold << ", " << map(indexold) << endl;; indexold += 1; add += 1; } } if (add == 0) index(k) = 0; cout << endl; cout << "index(k), add = " << index(k) << ", " << add << endl; } cout << "=========================================================" << endl; cout << " " << endl; //------------------------------------------------------------------------ TCanvas *c1 = new TCanvas("c1","", 600, 900); c1->Divide(2,3); c1->cd(1); hista->Draw(); c1->cd(2); histap->Draw(); c1->cd(3); histb->Draw(); c1->cd(4); histbp->Draw(); c1->cd(5); histc->Draw(); c1->cd(6); histcp->Draw(); } //=========================================================================