//Double_t kRad2Deg = 180./TMath::Pi();

void analp()
{
    TChain chain("Photons");
    chain.Add("cascade_500kpc_0*.root");
    chain.Add("cascade_500kpc_2*.root");
//    chain.Add("cascade_100kpc_0*.root");
//    chain.Add("cascade_100kpc_1*.root");

    MPhoton *p = new MPhoton;
    chain.SetBranchAddress("MPhoton.", &p);

    Int_t n = chain.GetEntries();

    cout << "Found " << n << " entries." << endl;

    if (n==0)
        return;

    MBinning binsx;
    binsx.SetEdgesLog(14, 1e4, 1e11);

    MBinning binspolx;
    MBinning binspoly1;
    MBinning binspoly2;
    binspolx.SetEdges(16, -180, 180);

    const Double_t ls = 299792458;        // [m/s]
    const Double_t ly = 3600.*24.*365.*ls; // [m/ly]
    const Double_t pc = 1./3.258;         // [pc/ly]

    binspoly1.SetEdges(20, 0, 2e-6*3600);
    binspoly2.SetEdges(20, 0, 1e-5*1e3/pc*365);

    TH1D h;
    h.SetName("Photons");
    h.SetTitle(" E^{2} Spectra ");
    h.SetXTitle("E\\gamma [GeV]");
    h.SetYTitle("E^{2} * Counts");
    MH::SetBinning(&h, &binsx);

    TH1D prim;
    MH::SetBinning(&prim, &binsx);

    TH2D r;
    TH2D a;
    MH::SetBinning(&a, &binspolx, &binspoly1);
    MH::SetBinning(&r, &binspolx, &binspoly2);

    Double_t weight = 0;
    Double_t alpha = -2;
    for (int i=0; i<n; i++)
    {
        chain.GetEntry(i);

        Double_t Ep = p->GetEnergy();

        if (i==0)
            weight = pow(Ep, alpha);

        if (p->IsPrimary())
        {
            weight = pow(Ep, alpha);
            prim.Fill(Ep, Ep*Ep * weight);
            continue;
        }

        h.Fill(Ep, Ep*Ep * weight);
        r.Fill(p->GetPhi()*kRad2Deg, p->GetR()*1e3/pc*365, weight);
        a.Fill(p->GetPsi()*kRad2Deg, p->GetTheta()*kRad2Deg*3600, weight);
    }

    delete p;

    gStyle->SetOptStat(10);

//    delete gROOT->FindObject("Analysis Arrival");
    c = MH::MakeDefCanvas("Analysis Arrival", "");
    c->Divide(2,2);

    c->cd(1);
    gPad->SetLogy();
    gPad->SetLogz();
    gPad->SetTheta(70);
    gPad->SetPhi(90);
    r.SetTitle(" Distance from Observer ");
    r.GetXaxis()->SetLabelOffset(-0.015);
    r.GetXaxis()->SetTitleOffset(1.1);
    r.GetXaxis()->SetRangeUser(1e4, 1e9);
    r.SetXTitle("\\Phi [\\circ]");
    r.SetYTitle("R [light days]");
    r.DrawCopy("surf1pol");

    c->cd(2);
    gPad->SetLogy();
    gPad->SetLogz();
    gPad->SetTheta(70);
    gPad->SetPhi(90);
    a.SetTitle(" Hit Angle for Observer ");
    a.GetXaxis()->SetLabelOffset(-0.015);
    a.GetXaxis()->SetTitleOffset(1.1);
    a.GetXaxis()->SetRangeUser(1e4, 1e9);
    a.SetXTitle("\\Psi [\\circ]");
    a.SetYTitle("\\Theta [\"]");
    a.DrawCopy("surf1pol");

    c->cd(3);
    gPad->SetLogy();
    TH1 *g=r.ProjectionY("p1");
    g->SetBit(kCanDelete);
    g->SetTitle(" Hit Observers Plain ");
    g->GetXaxis()->SetLabelOffset(0);
    g->GetXaxis()->SetTitleOffset(1.1);
    g->GetYaxis()->SetTitleOffset(1.3);
    g->SetXTitle("R [light days]");
    g->SetYTitle("Counts");
    g->Draw();

    c->cd(4);
    gPad->SetLogy();
    g=a.ProjectionY("p2");
    g->SetBit(kCanDelete);
    g->SetTitle("Hit Angle");
    g->GetXaxis()->SetLabelOffset(0);
    g->GetXaxis()->SetTitleOffset(1.1);
    g->GetYaxis()->SetTitleOffset(1.3);
    g->SetXTitle("\\Phi [\"]");
    g->SetYTitle("Counts");
    g->Draw();

  //  delete gROOT->FindObject("Analysis Photons");
    TCanvas *c = MH::MakeDefCanvas("Analysis Photons", "", 580, 870);
    c->Divide(1,2);

    c->cd(1);
    gPad->SetLogx();
    gPad->SetLogy();
    prim.SetFillStyle(0);
    h.SetFillStyle(0);
    h.GetXaxis()->SetLabelOffset(-0.015);
    h.GetXaxis()->SetTitleOffset(1.1);
//    h.GetXaxis()->SetRangeUser(1e4, 1e9);
    prim.SetMarkerStyle(kPlus);
    h.SetMarkerStyle(kMultiply);
    prim.SetMarkerColor(kRed);
    prim.SetLineColor(kRed);

    h.DrawCopy("P");
    prim.DrawCopy("Psame");
    prim.DrawCopy("Csame");
    h.DrawCopy("Csame");
 
    c->cd(2);
    gPad->SetLogx();
    TH1D div;
    MH::SetBinning(&div, &binsx);
    div.Divide(&h, &prim);
    div.SetTitle("Spectrum / Primara Spectrum");
    div.GetXaxis()->SetLabelOffset(-0.015);
    div.GetXaxis()->SetTitleOffset(1.1);
    div.SetXTitle("E\\gamma [GeV]");
    div.SetYTitle("Ratio");
    div.DrawCopy();
}