// Ulli Schwanke, Louise Oakes

void CountPhotons(TH1D* h){  
  const int N = h->GetNbinsX();
  for(int i=0;i<3;i++){
    double thresh = pow(10,i)+0.00001;
    int k = h->FindBin(thresh);  
    double photons = h->Integral(k,N);
    double sum = 0;
    double sumw = 0;
    for(int j=k;j<=N;j++){
      sum += h->GetBinCenter(j)*h->GetBinContent(j);
      sumw += h->GetBinContent(j);
    }
    std::cout << "Above " << thresh << " GeV: " << photons << " photons, " << sum/sumw << " GeV" << std::endl;
  }
}


Double_t fun(Double_t* E,Double_t* par){
  double x = E[0]/par[0];
  if(x<=1) return 1./par[0]*pow(x,-1.5)*exp(0.047-8.7*x+9.14*x*x-10.3*x*x*x);
  return 0;
}

void plot_spectrum(char * filename = "bm.log.model.1.dat"){
 
  TH1D * spec1 = new TH1D("spec1","#gamma spectrum",500,0,500);
  TH1D * spec2 = new TH1D("spec2","#gamma spectrum",500,0,500);
  TH1D * spec3 = new TH1D("spec3","#gamma spectrum",500,0,500);
  TH2D * frame = new TH2D("frame","frame",499,1,500,500,1e-9,20);
  frame->GetXaxis()->SetTitle("E_{#gamma} [GeV]");
  frame->GetYaxis()->SetTitle("dN/dE");
  frame->SetStats(0);
  frame->SetTitle(0);
  TF1 * func2 = new TF1("func2",fun,0,500.0,1);
  func2->SetParameter(0,500); //500 GeV
  TF1 * func3 = new TF1("func3",fun,0,500.0,1);
  ifstream fin(filename);
  int i=0;
  double dnde = 0;
  double mass = 0;
  while (!fin.eof()){   
    if (i>0){
      fin >> dnde;      
      spec1->SetBinContent(i,dnde);   
    } else {
      fin >> mass; 
      std::cout << "WIMP mass is " << mass << " GeV" << std::endl;
      func3->SetParameter(0,mass);
      fin >> dnde; //ignore
    }
    i++;    
  }
  for(int i=1;i<=spec2->GetNbinsX();i++){
    double e1 = spec2->GetBinLowEdge(i);
    double e2 = e1 + spec2->GetBinWidth(i);
    double s = func2->Integral(e1,e2);   
    spec2->SetBinContent(i,s);   
    s = func3->Integral(e1,e2);   
    spec3->SetBinContent(i,s);
  }
 
  gStyle->SetOptLogx(1);
  gStyle->SetOptLogy(1);

  TCanvas * can = new TCanvas("can","HAP iDM Workshop",900,750);
  spec1->SetStats(0);
  spec1->SetLineColor(kBlue);
  spec1->SetLineWidth(3);
  spec1->GetXaxis()->SetTitle("E_{#gamma} [GeV]");
  spec1->GetYaxis()->SetTitle("dN/dE");
  spec3->SetLineColor(kRed);

  frame->Draw();
  spec1->Draw("same");
  spec2->Draw("same");
  spec3->Draw("same");
  can->Print("plot_spectrum.png");

  CountPhotons(spec1);
  CountPhotons(spec2);
  CountPhotons(spec3);
}