source: trunk/Cosy/macros/testtrack.C

#include <iomanip.h>

void testtrack()
{
    TGraph g;

    TH2F h("Hist", "dX/dY",  77, -768/2-.5,  768/2+.5, 58, -576/2-.5,  576/2+.5); // 3
    TH1F hmag("HistMag", "Mag",  19,  0, 100);

    TArrayF x;
    TArrayF y;
    TArrayF mag;

    ifstream fin0("data/tracking_702.999039.txt");

    int i=0;

    while (!fin0.eof())
    {
        double d, p, m;
        fin0 >> d >> p >> m;
        if (!fin0)
            break;

        y.Set(i+1);
        x.Set(i+1);
        mag.Set(i+1);

        y.AddAt(d, i);
        x.AddAt(p, i);
        mag.AddAt(m, i);

        i++;

        hmag.Fill(m);
    }

    TCanvas *c = new TCanvas;
    c->Divide(2,2);
    c->cd(1);
    hmag.DrawCopy();

    for (i=0; i<mag.GetSize(); i++)
    {
        if (mag[i]<48+15 && mag[i]>48-15)
        {
            h.Fill(x[i], y[i]);
        }
    }

    c->cd(2);
    h.DrawCopy("surf4");

    c->cd(3);
    TH1D *hist = h.ProjectionX();
    hist->SetLineColor(kGreen);
    hist->Draw();
    hist->SetBit(kCanDelete);
    hist = h.ProjectionY();
    hist->SetLineColor(kRed);
    hist->Draw("same");
    hist->SetBit(kCanDelete);

    Int_t y1, y2, y3;
    h.GetMaximumBin(y1, y2, y3);

    double ymax = h.GetYaxis()->GetBinCenter(y2);
    double dy   = h.GetYaxis()->GetBinWidth(y2);

    double xmax = h.GetXaxis()->GetBinCenter(y1);
    double dx   = h.GetXaxis()->GetBinWidth(y1);


    cout << "*** " << h.GetEntries() << " " << y1 << " " << y2 << endl;

    cout << "Cut-X: " << xmax << " " << dx*1 << endl;
    cout << "Cut-Y: " << ymax << " " << dy*1 << endl;

    TGraph g2;
    Int_t pnt=0;
    for (i=0; i<mag.GetSize(); i++)
    {
        if (y[i]>ymax-1*dy && y[i]<ymax+1*dy &&
            x[i]>xmax-1*dx && x[i]<xmax+1*dx &&
            mag[i]>48-15 && mag[i]<48+15)
            g2.SetPoint(g2.GetN(), x[i], y[i]);
    }

    c->cd(4);
    g2.DrawClone("A*");

    cout << setprecision(3);
    cout << "X:  " << g2.GetMean(2) << "' +- " << g2.GetRMS(2) << "'" << endl;
    cout << "Y: " << g2.GetMean(1) << "' +- " << g2.GetRMS(1) << "'" << endl;
}