{
		 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 x,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]="../../Mapping/";
   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]; 
   Int_t *binx = new Int_t[N]; 
   Int_t *biny = new Int_t[N]; 
     
//
   X[N]=0;
   Y[N]=0.0;
   Int_t nData =1;   
   s =i=j=1;   
   x=0;
	 /*
   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); 

	 */
//----------------------------------------
// LEGGO FILE source in int.dat 
   ifstream in;
   in.open("w66prova.dat", 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;
		 binx[nData]=u;
		 biny[nData]=t;			
		 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
  
		TCanvas *cont=new TCanvas("contours","contours",10,10,800,600);
    cont->SetRightMargin(0.2);
    cont->SetFillColor(kWhite);
    gPad->SetGrid();

		
		TH2D *h2 = new TH2D("h2","",Dim,lmin-binsize*0.5,lmax+binsize*0.5,Dim,bmin-binsize*0.5,bmax+binsize*0.5);
			TH2D *h3 = new TH2D("h3","",Dim,lmin-binsize*0.5,lmax+binsize*0.5,Dim,bmin-binsize*0.5,bmax+binsize*0.5); 
			/*		TH2D *h3 = new TH2D("h3","",Dim,lmax-binsize*0.5,-lmin+binsize*0.5,Dim,bmin-binsize*0.5,bmax+binsize*0.5); */


		gStyle->SetOptStat(0000000);
		gStyle->SetOptFit(0000);

		// Disegno la Nube e la PSR  
		double l1=78.0;  //l1*=-1.0;//nube 
		double b1=2.3;   
		double largl=0.1; 
		double largb=0.2;
		double l2=78.2;  //l2*=-1;//psr
		double b2=2.1;   
		TEllipse *mc = new TEllipse(l1,b1,largl/2,largb/2);
		TMarker *psr= new TMarker (l2,b2,3);
		mc->SetLineColor(kWhite);
		mc->Draw();
		
		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]);
				h2->SetCellContent(Dim-i,j,Y[s]);
				//			h3->SetBinContent(k,Z[s]);
			}
		}
		
//------------------------//

		h2.GetXaxis()->SetTitle("l [longituge]");
		h2.GetYaxis()->SetTitle("b [latitude]");
		h3.GetXaxis()->SetTitle("l []");
		h3.GetYaxis()->SetTitle("b");
		
		TAxis *xaxish2=h2->GetXaxis();
		xaxish2->SetTickLength(0);
		xaxish2->SetLabelOffset(0.02);
		xaxish2->SetLabelSize(0.03);
		
		TAxis *yaxish2=h2->GetYaxis();
		yaxish2->SetTickLength(0);
		yaxish2->SetLabelOffset(0.02);
		yaxish2->SetLabelSize(0.03);
		
		TAxis *zaxish2=h2->GetZaxis();
		zaxish2->SetNoExponent(kTRUE);
		zaxish2->SetLabelSize(0.03);
		
		gStyle->SetPalette(1);
		h2->SetContour(80);
		TPaveLabel pl,pr;

		Float_t x1=0.37, y1=0.875, x2=0.57, y2=0.97;

			gPad->SetTheta(90); gPad->SetPhi(0.005);

			h2->Draw("surf3z");
		// h2->Draw("cont4z");

     
		TPaveText lom(0.7,0.,0.90,0.06); 
		lom.SetFillColor(kWhite);
		lom.SetBorderSize(0);
		TText *lo3=lom.AddText("ph/cm^2/s/sr");
     
		lo3.SetTextAngle(90);     
		lo3.SetTextAlign(22);
		lom.Draw();
	 

		//cont->Update();

		h3->Draw("same");
		mc->Draw("same");
		psr->Draw("same");
}