#include "TFile.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TGaxis.h"
#include "TRandom.h"
#include "TLegend.h"
#include "TPaveStats.h"
#include "TGraph.h"

void denseprim()
{
   
   char datafilegamma1[200]= "G:/cors/gamma/gammatot/distribuzioni-gammatot.txt";
   char datafileele1[200]= "G:/cors/gamma/gammatot/distribuzioni-eletot.txt";
   char datafilemuo1[200]= "G:/cors/gamma/gammatot/distribuzioni-muotot.txt";
   char datafilehad1[200]= "G:/cors/gamma/gammatot/distribuzioni-hadtot.txt";
   char dataout[200]="G:/cors/gamma/gamma100m_e_primari-dens.png";
   
   gStyle->SetOptFit();
   TCanvas *c36 = new TCanvas("c36","multigraph",1280,1024);
   float offx=1.1;
   float offy=1.7;
   float margr=0.18;
   float margl=0.12;
   gPad->SetLeftMargin(margl);
   gPad->SetRightMargin(margr);
   //c36->SetGrid();
    gStyle->SetOptFit(0001);
	Double_t a1=1e-5;
	Double_t b1=1e-5;
	Double_t c1=-9e-4;
	
   TMultiGraph *mg = new TMultiGraph();
   
    //gamma secondari
   TGraphErrors *gamma10_1 = new TGraphErrors(datafilegamma1,"%lg %*lg %lg %*lg");
   gamma10_1->SetMarkerColor(kBlack);
   gamma10_1->SetMarkerStyle(2);
   gamma10_1->SetTitle("#gamma");
   TF1 *fitg10_1 = new TF1("fitg10_1", "[0]*x*x+[1]*x+[2]", 1, 5e3);
   fitg10_1->SetParameters(a1, b1, c1);
   fitg10_1->SetParName(0,"a");
   fitg10_1->SetParName(1,"b");
   fitg10_1->SetParName(2,"c");
   fitg10_1->SetLineColorAlpha(kBlack,0.35);
   gamma10_1->Fit("fitg10_1","q");
   //gr1->Fit("pol1","q");
   mg->Add(gamma10_1);
   
   //elettroni secondari
   
   TGraphErrors *ele10_1 = new TGraphErrors(datafileele1,"%lg %*lg %lg %*lg");
   ele10_1->SetMarkerColor(kRed);
   ele10_1->SetMarkerStyle(20);
   TF1 *fite10_1 = new TF1("fite10_1", "[0]*x*x+[1]*x+[2]", 1, 5e3);
   fite10_1->SetParameters(a1, b1, c1);
   fite10_1->SetParName(0,"a");
   fite10_1->SetParName(1,"b");
   fite10_1->SetParName(2,"c");
   fite10_1->SetLineColorAlpha(kRed,0.35);
   ele10_1->Fit("fite10_1","q");
   //gr1->Fit("pol1","q");
   mg->Add(ele10_1);

    //muoni secondari
   
   TGraphErrors *muo10_1 = new TGraphErrors(datafilemuo1,"%lg %*lg %lg %*lg");
   muo10_1->SetMarkerColor(kGreen);
   muo10_1->SetMarkerStyle(25);
   TF1 *fitm10_1 = new TF1("fitm10_1", "[0]*x*x+[1]*x+[2]", 1, 5e3);
   fitm10_1->SetParameters(a1, b1, c1);
   fitm10_1->SetParName(0,"a");
   fitm10_1->SetParName(1,"b");
   fitm10_1->SetParName(2,"c");
   fitm10_1->SetLineColorAlpha(kGreen,0.35);
   muo10_1->Fit("fitm10_1","q");
   //gr1->Fit("pol1","q");
   mg->Add(muo10_1);
   
    //adroni secondari
   
   TGraphErrors *had10_1 = new TGraphErrors(datafilehad1,"%lg %*lg %lg %*lg");
   had10_1->SetMarkerColor(kBlue);
   had10_1->SetMarkerStyle(30);
   TF1 *fith10_1 = new TF1("fith10_1", "[0]*x*x+[1]*x+[2]", 1, 5e3);
   fith10_1->SetParameters(a1, b1, c1);
   fith10_1->SetParName(0,"a");
   fith10_1->SetParName(1,"b");
   fith10_1->SetParName(2,"c");
   fith10_1->SetLineColorAlpha(kBlue,0.35);
   had10_1->Fit("fith10_1","q");
   //gr1->Fit("pol1","q");
   mg->Add(had10_1);
   
   mg->Draw("ap");
   //force drawing of canvas to generate the fit TPaveStats
   //mg->SetTitle("Densita' di #gamma ed e^{#pm} a 600 m dal core degli EAS indotti da #font[12]{p}");
  mg->SetTitle("Densita' dei secondari a 100 m dal core degli EAS indotti da #font[12]{#gamma}");
   mg->GetXaxis()->SetTitle("Energia_{#font[12]{#gamma}} (GeV)");
   mg->GetYaxis()->SetTitle("Densita'(Particelle/m^{2})");
   mg->GetYaxis()->SetTitleOffset(offy);
   c36->Modified();
   //force drawing of canvas to generate the fit TPaveStats
   c36->Update();
    TLegend* leg10 = new TLegend(0.55, 0.7, .35, .8);
    leg10->SetNColumns(1);
    leg10->AddEntry(gamma10_1, "#gamma", "lp");
    leg10->AddEntry(ele10_1, "e^{#pm}", "lp");
    
    leg10->AddEntry(muo10_1, "#mu^{#pm}", "lp");
    leg10->AddEntry(had10_1, "Adroni", "lp");
    
    leg10->Draw("same");
    TPaveStats *stats1 = (TPaveStats*)gamma10_1->GetListOfFunctions()->FindObject("stats");
    TPaveStats *stats2 = (TPaveStats*)ele10_1->GetListOfFunctions()->FindObject("stats");
    TPaveStats *stats3 = (TPaveStats*)muo10_1->GetListOfFunctions()->FindObject("stats");
    TPaveStats *stats4 = (TPaveStats*)had10_1->GetListOfFunctions()->FindObject("stats");
   stats1->SetTextColor(kBlack);
   stats2->SetTextColor(kRed);
   stats3->SetTextColor(kGreen);
   stats4->SetTextColor(kBlue);
   
   stats1->SetX1NDC(0.80); stats1->SetX2NDC(0.98);
   stats1->SetY1NDC(0.77); stats1->SetY2NDC(0.92);
   stats2->SetX1NDC(0.80); stats2->SetX2NDC(0.98);
   stats2->SetY1NDC(0.60); stats2->SetY2NDC(0.75);
   stats3->SetX1NDC(0.80); stats3->SetX2NDC(0.98);
   stats3->SetY1NDC(0.43); stats3->SetY2NDC(0.58);
   stats4->SetX1NDC(0.80); stats4->SetX2NDC(0.98);
   stats4->SetY1NDC(0.26); stats4->SetY2NDC(0.41);
   c36->Modified();
   gROOT->ProcessLine("gROOT->SetBatch()");
   c36.Print(dataout);delete c36;
   gROOT->ProcessLine("gROOT->SetBatch(kFALSE)");
   exit();
  
   }