Index: trunk/Mars/hawc/atmprof.C
===================================================================
--- trunk/Mars/hawc/atmprof.C	(revision 19770)
+++ trunk/Mars/hawc/atmprof.C	(revision 19770)
@@ -0,0 +1,220 @@
+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;
+}