TF1* flatdist;
TF1* flatdist91;
char gtitle[200];


void DrawLDF() {
  gStyle->SetOptStat(0);
  e_g->Draw();
  e_g->GetXaxis()->SetRange(0,80);
  gPad->Update();
  gPad->SetLogy();
  gPad->SetGrid();
  m_g->Draw("same");
  TH1F* hsum = (TH1F*)e_g->Clone();
  hsum->Add(m_g);
  hsum->Draw("same");
  hsum->SetLineColor(kBlack);
  headerInfo();
  Double_t energy = TMath::Power(10, 9+header->fLogE);
  ldist(energy, header->fTheta);
  ldist91(energy, header->fTheta);
  flatdist->Draw("same");
  flatdist91->Draw("same");
  TLegend* leg = new TLegend(0.5, 0.6, 0.95, 0.85);
  leg->AddEntry(e_g,"Electrons from Aires", "L");
  leg->AddEntry(m_g,"Muons from Aires", "L");
  leg->AddEntry(hsum,"Electrons + Muons", "L");
  leg->AddEntry(flatdist,"Lat Dist from CTShower","L");
  leg->AddEntry(flatdist91,"Lat Dist from CTShower (M_{R} = 91.6 m)","L");
  leg->Draw();
  e_g->SetTitle(gtitle);
  e_g->SetYTitle("Intensity /m^{2}");
  gPad->Update();
  

}

void ldist(Double_t aEnergy, Double_t aThetaDegrees) {
  flatdist = new TF1("latdist",cjjLateralDistribution,1,10000,2);
  Double_t theta = TMath::DegToRad()*aThetaDegrees;
  Double_t logC = EnergytoLogC(theta, aEnergy);
  Double_t C = TMath::Exp(logC);
  flatdist->SetParameter(0, TMath::Pi() - theta);
  flatdist->SetParameter(1, C);
  flatdist->SetLineColor(8);

}

void ldist91(Double_t aEnergy, Double_t aThetaDegrees) {
  flatdist91 = new TF1("latdist91",cjjLateralDistribution_Moliere91,1,10000,2);
  Double_t theta = TMath::DegToRad()*aThetaDegrees;
  Double_t logC = EnergytoLogC_Moliere91(theta, aEnergy);
  Double_t C = TMath::Exp(logC);
  flatdist91->SetParameter(0, TMath::Pi() - theta);
  flatdist91->SetParameter(1, C);
  flatdist91->SetLineColor(5);
  flatdist91->SetLineStyle(2);
}

void headerInfo() {
  char primary[64];
  if( header->fPrimaryParticleCode ==14 ) sprintf(primary,"Primary: proton");
  else sprintf(primary,"Primary: %d",  header->fPrimaryParticleCode);
  char za[64];
  sprintf(za,"ZA = %3.0f deg", header->fTheta);
  char energy[64];
  sprintf(energy,"Energy = 10^{%4.2f} eV",9+header->fLogE);
  sprintf(gtitle,"%s, %s, %s",primary,energy,za);
}