Double_t Log10(Double_t lg)
{
    // Avoid numerical limits
    return lg<=1e-6 ? -6 : log10(lg);
}

// The mass overburden in each layer is given by
// T = AATM + BATM * exp(-h/CATM)
// Note that the AATM are correlated with the layer border
Double_t fcn(Double_t *xx, Double_t *par)
{
    Double_t *A = par;
    Double_t *B = par+5;
    Double_t *C = par+10;
    Double_t *H = par+15;

    Double_t x = xx[0]*1e5;

    double Ai = A[0];

    for (int i=0; i<4; i++)
    {
        if (x>=H[i] && x<H[i+1])
            return Log10(Ai + B[i]*exp(-x/C[i]));

        Ai += B[i]  *exp(-H[i+1]/C[i]);    // y-value of current bin function [i]   at transition border H[i+1] to next layer
        Ai -= B[i+1]*exp(-H[i+1]/C[i+1]);  // y-value of next    bin function [i+1] at transition border H[i+1] to next layer
    }

    return Log10(Ai + B[4]*exp(-x/C[4]));
}

void atmprof()
{
    // As in CORSIKA the atmosphere for atmabs.dat is required up to 50 km
    // the best is to do a fit which is bext up to at least 50 km
    double MAX = 50.5;

    // MAGIC Winter (MAGIC Summer: atmprof9)
    TGraph atmprof("/home/tbretz/SW/corsika-76900/bernlohr/atmprof8.dat", "%lg %*lg %lg, %*lg");

    // The last point has a thickness of 0 (log!)
    atmprof.RemovePoint(atmprof.GetN()-1);

    // Calculate the logarithm (better fit)
    double miny = log10(atmprof.GetY()[0]);
    for (int i=0; i<atmprof.GetN(); i++)
    {
        atmprof.GetY()[i] = log10(atmprof.GetY()[i]);
        if (atmprof.GetX()[i]<MAX)
            miny = TMath::Min(miny, atmprof.GetY()[i]);
    }

    // =============================================================================

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

    atmprof.SetMarkerStyle(kFullDotMedium);
    atmprof.GetHistogram()->GetXaxis()->SetRangeUser(0, MAX);
    atmprof.GetHistogram()->GetYaxis()->SetRangeUser(miny, atmprof.GetY()[0]*1.05);
    atmprof.DrawClone("ALP");

    TMarker m;
    m.SetMarkerStyle(kFullDotLarge);

    // =============================================================================

    // Define the fit function betwene 0 and 50km with 15 free parameters
    TF1 f("thickness", fcn, 0, MAX, 4*5);
    f.SetNpx(10000);

    // This is the Keilhauer U.S. Standard Atmosphere (as starting values)
    Double_t params[] =
    {
        -149.801663, -57.932486, 0.63631894, 4.35453690e-4, 0.01128292, // A
        1183.6071,   1143.0425,  1322.9748,  655.67307, 1,              // B
        954248.34,   800005.34,  629568.93,  737521.77, 1e9,            // C
        0, 7.0e5, 11.4e5, 37.0e5, 100e5,                                // H

    };
    f.SetParameters(params);

    // =============================================================================
    // Calc residuals

    TGraph g1;
    for (int i=0; i<atmprof.GetN(); i++)
    {
        double x = atmprof.GetX()[i];
        if (x>MAX)
            break;

        double y = atmprof.GetY()[i];
        g1.SetPoint(g1.GetN(), x, 1-pow(10, f.Eval(x))/pow(10, y));
    }


    // =============================================================================
    // Plot and print

    m.SetMarkerColor(kRed);
    f.SetLineColor(kRed);

    f.DrawCopy("same");

    cout << endl;
    for (int i=0; i<5; i++)
    {
        double A[5] = { f.GetParameters()[0], 0, 0, 0, 0 };
        double *B = f.GetParameters()+5;
        double *C = f.GetParameters()+10;
        double *H = f.GetParameters()+15;

        A[1] = A[0] + B[0]*exp(-H[1]/C[0]) - B[1]*exp(-H[1]/C[1]);
        A[2] = A[1] + B[1]*exp(-H[2]/C[1]) - B[2]*exp(-H[2]/C[2]);
        A[3] = A[2] + B[2]*exp(-H[3]/C[2]) - B[3]*exp(-H[3]/C[3]);
        A[4] = A[3] + B[3]*exp(-H[4]/C[3]) - B[4]*exp(-H[4]/C[4]);

        cout << Form("A[%d]=%13f   B[%d]=%13f   C[%d]=%13f", i, A[i], i, B[i], i, C[i]);
        if (i<4)
        {
            cout << Form("   H=%13f", H[i+1]/1e5);
            m.DrawMarker(H[i+1]/1e5, f.Eval(H[i+1]/1e5));
        }
        cout << endl;
    }

    // =============================================================================

    // The parameters above 50km should not be left free
    f.FixParameter( 1, params[1]);
    f.FixParameter( 2, params[2]);
    f.FixParameter( 3, params[3]);
    f.FixParameter( 4, params[4]);
    f.FixParameter( 9, params[9]);
    f.FixParameter(14, params[14]);
    f.FixParameter(15, 0);

    // Fit
    float max = 5;
    while (max<MAX)
    {
        atmprof.Fit(&f, "", "",  0, max);
        max += 5;
    }

    // =============================================================================
    // Plot and print

    m.SetMarkerColor(kBlue);
    f.SetLineColor(kBlue);

    f.DrawCopy("same");

    cout << endl;
    for (int i=0; i<5; i++)
    {
        double A[5] = { f.GetParameters()[0], 0, 0, 0, 0 };
        double *B = f.GetParameters()+5;
        double *C = f.GetParameters()+10;
        double *H = f.GetParameters()+15;

        A[1] = A[0] + B[0]*exp(-H[1]/C[0]) - B[1]*exp(-H[1]/C[1]);
        A[2] = A[1] + B[1]*exp(-H[2]/C[1]) - B[2]*exp(-H[2]/C[2]);
        A[3] = A[2] + B[2]*exp(-H[3]/C[2]) - B[3]*exp(-H[3]/C[3]);
        A[4] = A[3] + B[3]*exp(-H[4]/C[3]) - B[4]*exp(-H[4]/C[4]);

        cout << Form("A[%d]=%13f   B[%d]=%13f   C[%d]=%13f", i, A[i], i, B[i], i, C[i]);
        if (i<4)
        {
            cout << Form("   H=%13f", H[i+1]/1e5);
            m.DrawMarker(H[i+1]/1e5, f.Eval(H[i+1]/1e5));
        }
        cout << endl;
    }

    // =============================================================================
    // Plot residuals

    c->cd(2);

    TGraph g2;
    for (int i=0; i<atmprof.GetN(); i++)
    {
        double x = atmprof.GetX()[i];
        if (x>MAX)
            break;

        double y = atmprof.GetY()[i];
        g2.SetPoint(g2.GetN(), x, 1-pow(10, f.Eval(x))/pow(10, y));
    }

    g1.SetLineColor(kRed);
    g1.SetMarkerColor(kRed);
    g1.SetMarkerStyle(kFullDotMedium);

    g2.SetLineColor(kBlue);
    g2.SetMarkerColor(kBlue);
    g2.SetMarkerStyle(kFullDotMedium);

    TH1C h("", "Residuals", 100, 0, MAX);
    Double_t xmin, ymin, xmax, ymax;
    g1.ComputeRange(xmin, ymin, xmax, ymax);
    Double_t maxy = TMath::Max(fabs(ymin)*1.05, fabs(ymax)*1.05);
    h.SetMinimum(-maxy);
    h.SetMaximum( maxy);
    h.SetStats(kFALSE);
    h.DrawCopy();

    g1.GetHistogram()->GetXaxis()->SetRangeUser(0, MAX);
    g1.DrawClone("LP");
    g2.DrawClone("LP");

    g2.ComputeRange(xmin, ymin, xmax, ymax);
    maxy = TMath::Max(fabs(ymin)*1.05, fabs(ymax)*1.05);

    cout << "\nMax residual after fit: " << maxy*100 << " %\n" << endl;
}