{
gROOT->Reset();
#include <iostream.h>
#include <string.h>
#include <iomanip.h>
//**********************************************************
 cout<< "\n\n******************************************************************"<< endl;
   cout<< "**Programma Mappe, output:intensity map  source+bkg"<< endl;
   cout<< "**legge il file int.dat e calcola 'Dim'= dimensione mappa"<< endl;
   cout<< "**ordina da solo l'istogramma se il vettore delle energie"<< endl;
   cout<< "**e' fatto con x variabile a y fissato"  << endl;
   cout<< "**CAMBIARE soltanto il range dell'istogramma" << endl;
cout<< "******************************************************************"<< endl;

//-------------------------------------


//-----------------------------------------
//  V A R I A B I L I 
   Int_t N=20000;   
   Float_t z,v,w,s1,s2,xmin,xmax,ymin,ymax;
   Double_t y,counts,lmin,lmax,bmin,bmax,av_b,av_l,binsize;
   Int_t n,i,j,k,s,u,t,Dim;//Dim=dim in pixel della mappa
   char f1[20], f2[20];
   char bkg[20]= "";   
   char file1[60],file2[60],file3[60],file4[60],file5[60];
   char path1[60]="";
   char path2[60]="../../Mapping/";
   
   
   Int_t *X = new Int_t[N];
   Double_t *Y = new Double_t[N]; 
   Double_t *Z = new Double_t[N]; 
     
//
   X[N]=0;
   Y[N]=0.0;
   Int_t nData =1;   
   s =i=j=1;   
  

   cout << "Nome file sorgente(es.'int1s100.dat')\n";
   cin  >> file1;
   cout << "Nome file sorgente + bkg(es.'int1s+bkg500.dat')\n";
   cin  >> file2; 

   strcat(path1,file1); 
   strcpy(file1,  path1);
   
   strcat(path2,file2); 
   strcpy(file2,  path2); 

   
cout << "\n\nNome sorgente e taglio in Mev in istogrammi\n";
cin  >> f1;
strcpy(f2,  f1); //copia la stringa
strcat(f2,bkg);     //concatena le due stringhe 

//----------------------------------------
// LEGGO FILE source in int.dat 
 ifstream in;
 in.open(file1, ios::in);

nData=0;
 while (1) {
   if(nData==0){
     in >> lmin >> lmax>> bmin >> bmax;
     nData=1;}
   in >> y >>u>>t;  // u e t sono i valori di x e y
      if (!in.good()) break;  
      Y[nData]=y;
      nData++;
 }
   printf("**Trovati in source %d points\n",nData-1);
   in.close();
     Dim=TMath::Power(nData-1,0.5);
 
//----------------------------------------
// LEGGO FILE source+bkg
 ifstream in;
 in.open(file2, ios::in);

nData=0;
 while (1) {
   if(nData==0){
     in >> lmin >> lmax>> bmin >> bmax;
     nData=1;}
   in >> y >>u>>t;  // u e t sono i valori di x e y
      if (!in.good()) break;  
      Z[nData]=y;
      nData++;
 }
   printf("**Trovati in source+bkg %d points\n",nData-1);
   in.close();
  lmin*=10.0;
  lmax*=10.0;
  bmin*=10.0;
  bmax*=10.0;
  lmin=ceil(lmin);
  lmax=floor(lmax);
  bmin=ceil(bmin);
  bmax=floor(bmax);
  lmin/=10.0;
  lmax/=10.0;
  bmin/=10.0;
  bmax/=10.0;
  av_l=(lmax+lmin)/2.0; 
  av_b=(bmax+bmin)/2.0;
  binsize=(lmax-lmin)/(Dim-1);
  double n_bins1 =(lmax-lmin)/binsize+1;
  int n_bins=n_bins1;

printf("\n %f, %f, %f, %f , %f\n",lmin,lmax,bmin,bmax,binsize);
//--------------------------
//**************************
// LISTA DEGLI ISTOGRAMMI
  
  
  TH2D *h2 = new TH2D("h2",f1,Dim,-lmax-binsize*0.5,-lmin+binsize*0.5,Dim,bmin-binsize*0.5,bmax+binsize*0.5);  
  TH2D *h3 = new TH2D("h3",f2,Dim,-lmax-binsize*0.5,-lmin+binsize*0.5,Dim,bmin-binsize*0.5,bmax+binsize*0.5);  
  
// Disegno la Nube e la PSR  
 double l1=78.0;  l1*=-1.0;//nube 
 double b1=2.3;   
 double largl=0.5; 
 double largb=0.5;
 double l2=6.1;  l2*=-1;//psr
 double b2=-0.1;   
 mc = new TEllipse(l1,b1,largl/2,largb/2);
 psr= new TMarker (l2,b2,3);

   for(int j=1;j<=Dim;j++) {     //bin(i,j) riempie prima x a y fisso
      for(int i=1;i<=Dim;i++){
	counts=counts+Y[s];
	s=Dim*(j-1)+i;           //
        k=(Dim+2)*j+i;              //relazione tra gli indici #bin(i,j)=k
          h2->SetBinContent(k,Y[s]);
          h3->SetBinContent(k,Z[s]);
	  }
     }    
//*************************//
   h2.GetXaxis()->SetTitle("l");
   h2.GetYaxis()->SetTitle("b");
   h3.GetXaxis()->SetTitle("l");
   h3.GetYaxis()->SetTitle("b");

   TAxis *xaxish2=h2->GetXaxis();
 //  xaxish2->SetLabelColor(0);
   TAxis *xaxish3=h3->GetXaxis();
   xaxish3->SetLabelColor(0);

   TF1 *f_1= new TF1("f_1","-x",lmin,lmax);
   TGaxis *xnew= new TGaxis(-lmax,bmin-binsize*0.,-lmin,bmin-binsize*0.,"f_1");

//     TGaxis *xnew= new TGaxis(-lmax,bmin,-lmin,bmin,lmax,lmin);

   
   
   gStyle->SetPalette(1);
   
   TPaveLabel pl;

   Float_t x1=0.37, y1=0.875, x2=0.57, y2=0.97;
//   h2.GetNent()->SetNent=Counts;
   Int_t cancolor = 17;

TCanvas c2h("c2h","2-d options",10,10,800,600);
  c2h.Divide(2,2);
  c2h.SetFillColor(cancolor);
   c2h.cd(1);
     h2.Draw("surf3");
   c2h.cd(2);   
     h2.Draw("lego1"); 
   c2h.cd(3);
     h3.Draw("surf3");   
   c2h.cd(4);
     h3.Draw("lego1"); 


TCanvas cont("contours","contours",100,100,800,600);
  cont->SetRightMargin(0.2);
  cont.Divide(2,2);
  gPad->SetGrid();
  cont.SetFillColor(cancolor);
   cont.cd(1);
     gPad->SetGrid();
     gPad->SetRightMargin(22);
     h2.Draw("contz"); pl.DrawPaveLabel(x1,y1,x2,y2,"ph/cm^2 s sr","blNDC");
     mc.Draw();
     psr.Draw();
     xnew.Draw();
   cont.cd(2);
     gPad->SetGrid();
     h2.Draw("colz"); pl.DrawPaveLabel(x1,y1,x2,y2,"ph/cm^2 s sr","brNDC");
     xnew.Draw();
   cont.cd(3);
     gPad->SetGrid();
     h3.Draw("contz"); pl.DrawPaveLabel(x1,y1,x2,y2,"ph/cm^2 s sr","brNDC");
     xnew.Draw();
     mc.Draw();
     psr.Draw();
   cont.cd(4);
     gPad->SetGrid();
     h3.Draw("cont1"); pl.DrawPaveLabel(x1,y1,x2,y2,"ph/cm^2 s sr","brNDC");
     xnew.Draw();

}