#include <iostream>
#include <fstream>
#include <string>
#include <sstream>
#include <cmath>

void Graph_Coherent_Detect(){

   ifstream input_no_trig("Detect_0_Alp_No_Trig.txt");
   ifstream input_trig("Detect_0_Alp_Trig.txt");

   vector<double> xmed_no_trig, xmed_trig;
   vector<double> y_no_trig, y_trig;
   vector<double> ex_no_trig, ex_trig;
   vector<double> ey_no_trig, ey_trig;

   Double_t energy1_trig=0;
   Double_t energy2_trig=0;
   Double_t energymed_trig=0;
   Double_t flux_trig=0;
   Double_t flux_err_trig=0;
   Double_t flux_err_rel_trig=0;

   Double_t energy1_no_trig=0;
   Double_t energy2_no_trig=0;
   Double_t energymed_no_trig=0;
   Double_t flux_no_trig=0;
   Double_t flux_err_no_trig=0;
   Double_t flux_err_rel_no_trig=0;

   if(!input_trig){cout<<"errore apertura file all_sim"<<endl;}
   if(!input_no_trig){cout<<"errore apertura file depth"<<endl;}

   while(input_no_trig>>energy1_no_trig>>energy2_no_trig>>flux_no_trig>>flux_err_no_trig){
      energymed_no_trig=((energy1_no_trig+energy2_no_trig)/2)*1000;
      xmed_no_trig.push_back(energymed_no_trig);
      flux_err_no_trig=((flux_err_no_trig)/100)*flux_no_trig;

      if( flux_no_trig==0) {
         flux_err_rel_no_trig=0;
      } else{
         flux_err_rel_no_trig=flux_err_no_trig/flux_no_trig;
      }

      flux_err_no_trig=flux_err_rel_no_trig*flux_no_trig;
      y_no_trig.push_back(flux_no_trig);
      ex_no_trig.push_back(0);
      ey_no_trig.push_back(flux_err_no_trig);
   }

   Int_t j_no_trig=xmed_no_trig.size();

   while(input_trig>>energy1_trig>>energy2_trig>>flux_trig>>flux_err_trig){
      energymed_trig=((energy1_trig+energy2_trig)/2)*1000;
      xmed_trig.push_back(energymed_trig);
      flux_err_trig=((flux_err_trig)/100)*flux_trig;

      if(flux_trig==0){
         flux_err_rel_trig=0;
      } else{
         flux_err_rel_trig=flux_err_trig/flux_trig;
      }
      flux_err_trig=flux_err_rel_trig*flux_trig;
      y_trig.push_back(flux_trig);
      ex_trig.push_back(0);
      ey_trig.push_back(flux_err_trig);
   }

   Int_t j_trig=xmed_trig.size();

   TCanvas *spectrum=new TCanvas("Spectrum","Spectrum");

   TMultiGraph *g = new TMultiGraph();


   TGraphErrors *g_trig=new TGraphErrors(j_trig,&xmed_trig[0],&y_trig[0],&ex_trig[0],&ey_trig[0]);
   g_trig->SetLineWidth(2);
   g_trig->SetLineColor(2);

   g->Add(g_trig);
   g->Draw("APC");

   TLegend *legend = new TLegend(0.6,0.7,0.89,0.88);

   legend->SetFillColor(0);
   legend->SetLineColor(1);
   legend->SetTextSize(0.035);
   legend->AddEntry(g_trig,"Neutron Trigger","lep");
   legend->Draw();

   g->GetXaxis()->SetTitle("Energy (MeV)");
   g->GetYaxis()->SetTitle("Energy Deposited");

   xmed_no_trig.clear();
   y_no_trig.clear();
   ex_no_trig.clear();
   ey_no_trig.clear();

   xmed_trig.clear();
   y_trig.clear();
   ex_trig.clear();
   ey_trig.clear();
}