source: trunk/Mars/hawc/atmabs.C@ 19859

void atmabs()
{
    // =============================================================================
    // measured flux (lambda) = emitted flux (lambda) * trans(H, lambda)
    // =============================================================================
    TH2F hist1 ("TransCeres",         "Ceres Transmission",            700-250, 250, 700, 50, 0, 50);
    TH2F hist1a("TransCeresRayleigh", "Ceres Transmission (Rayleigh)", 700-250, 250, 700, 50, 0, 50);
    TH2F hist1b("TransCeresOzone",    "Ceres Transmission (Ozone)",    700-250, 250, 700, 50, 0, 50);
    TH2F hist1c("TransCeresAerosol",  "Ceres Transmission (Aerosol)",  700-250, 250, 700, 50, 0, 50);

    const char *ozone    = "resmc/atmosphere-ozone.txt";
    const char *aerosols = "resmc/atmosphere-aerosols.txt";

    MAtmosphere atm(ozone, aerosols);

    cout << atm.HasValidAerosol() << " " << atm.HasValidOzone() << endl;

    for (int h=0; h<50; h++)
    {
        for (int l=250; l<700; l++)
        {
            MPhotonData ph;
            ph.SetWavelength(l);
            ph.SetProductionHeight(h*1e3*1e2);

            hist1a.Fill(l, h, atm.CalcTransmission(h*1e5, l, 0));
            hist1b.Fill(l, h, atm.CalcOzoneAbsorption(h*1e5, l, 0));
            hist1c.Fill(l, h, atm.CalcAerosolAbsorption(h*1e5, l, 0));
            hist1.Fill(l, h, atm.GetTransmission(ph));
        }
    }

    // =============================================================================
    // The CCD Photometric Calibration Cookbook
    // 8 Atmospheric Extinction and Air Mass
    // http://star-www.rl.ac.uk/star/docs/sc6.htx/sc6se8.html
    //
    // measured magnitude (lambda) = emitted magnitude (lambda) + extco(lambda) * X(z)
    //
    //  ==> measured flux = emitted flux * exp(extco * X(h))
    //
    // =============================================================================

    TH2F hist2("AtmAbs", "Corsika Atmosphere atmabs.txt", (700-180)/5+1, 178.5, 702.5, 51, -0.5, 50.5);

    ifstream fin("/home/tbretz/SW/corsika-76900/run/atmabs.dat");

    TString buf;
    buf.ReadLine(fin);
    cout << buf << endl;

    while (1)
    {
        double wl;
        fin >> wl;
        if (!fin)
            break;

        for (int h=0; h<51; h++)
        {
            double cc;
            fin >> cc;
            hist2.Fill(wl, h, cc);
        }
    }

    // =============================================================================
    hist1.SetMinimum(0);
    hist1.SetMaximum(1);
    hist1.SetContour(500);

    hist1a.SetMinimum(0);
    hist1a.SetMaximum(1);
    hist1a.SetContour(500);

    hist1b.SetMinimum(0);
    hist1b.SetMaximum(1);
    hist1b.SetContour(500);

    hist1c.SetMinimum(0);
    hist1c.SetMaximum(1);
    hist1c.SetContour(500);

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

    gStyle->SetOptStat(11);

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

    c->cd(1);
    hist1a.DrawCopy("colz");

    c->cd(2);
    hist1b.DrawCopy("colz");

    c->cd(3);
    hist1c.DrawCopy("colz");

    c->cd(4);
    hist1.DrawCopy("colz");

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

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

    c->cd(1);
    hist1.DrawCopy("colz");

    c->cd(2);
    hist2.SetMinimum(0);
    //hist2.SetMaximum(1);
    hist2.SetContour(500);
    hist2.DrawCopy("colz");

    // =============================================================================
    //
    //  see CORSIKA ATABSO()
    //
    //  CORSIKA calculates a coefficient COATEX which is the value interpolated
    //  in wavelength and height from the table.
    //
    //  From this value, the interpolated value at observation level is subtraced.
    //
    //  The probability to survive (=transmission) is then exp(-COATEX/WEMIS)
    //
    //  WEMIS defined as the particles z-direction cosine (-> correction for
    //  zenith angle (flat atmosphere))
    //
    // =============================================================================

    TH2F hist4("TransCorsika", "Corsika Transmission", (700-180)/5+1, 178.5, 702.5, 51, -0.5, 50.5);
    for (int x=1; x<=hist2.GetNbinsX(); x++)
    {
        double sum = 0;
        for (int y=1; y<=hist2.GetNbinsY(); y++)
        {
            const double cosz = 1;
            hist4.SetBinContent(x, y, exp(-hist2.GetBinContent(x, y)/cosz));
        }
    }

    c->cd(4);
    hist4.SetMinimum(0);
    hist4.SetMaximum(1);
    hist4.SetContour(500);
    hist4.DrawCopy("colz");

    c->cd(3);
    TH1D *p1 = hist1.ProjectionX("TransCeres50km",   hist1.GetNbinsY(), hist1.GetNbinsY());
    TH1D *p4 = hist4.ProjectionX("TransCorsika50km", hist4.GetNbinsY(), hist4.GetNbinsY());
    TH1D *p1b = hist1.ProjectionX("TransCeres5km",   5, 5);
    TH1D *p4b = hist4.ProjectionX("TransCorsika5km", 5, 5);

    p1->SetNameTitle("Trans", "Transmission");

    p4->SetLineColor(kBlue);
    p4b->SetLineColor(kBlue);
    p1->SetMinimum(0);
    p1->SetMaximum(1);
    p1->Draw();
    p4->Draw("same");
    p1b->Draw("same");
    p4b->Draw("same");

    TLatex text;
    text.DrawLatex(275, 0.850, "Ceres");
    text.DrawLatex(340, 0.250, "50km");
    text.DrawLatex(320, 0.650, "5km");
    text.SetTextColor(kBlue);
    text.DrawLatex(275, 0.925, "Corsika");
}