////////////////////////////////////////////////////////////////////////////////////
//
//             _____ Plot Optimized cuts _____
//
//  Take as input an output from programs/optimizeCuts and plot the results
//  You only need to specify the input file name
//
//                 Jose Flix    <jflix@ifae.es>
//                 Javier Rico  <jrico@ifae.es>
////////////////////////////////////////////////////////////////////////////////////


void plotOptimal(TString fname="optimize.out")
{
  gROOT->Reset();
  gStyle->SetCanvasColor(0);
  gStyle->SetCanvasBorderMode(0);
  gStyle->SetPadBorderMode(0);
  gStyle->SetFrameBorderMode(0);
  gStyle->SetOptStat(000000);
  gStyle->SetPalette(1);

  ifstream ifun(fname.Data());
  Float_t width,length,excess,significance;

  Float_t wmin=9999999;
  Float_t lmin=9999999;
  Float_t wmax=-9999999;
  Float_t lmax=-9999999;
  Float_t emax=-9999999;
  Float_t smax=-9999999;

  Float_t dw=-1;
  Float_t dl=-1;
  Float_t prevw=-1;
  Float_t prevl=-1;

  // first loop to get parameters
  Int_t i=0;
  while(ifun >> width)
    {
      // read data
      ifun >> length;
      ifun >> excess;
      ifun >> significance;
      
      // keep max/min
      if(width<wmin)  wmin=width;
      if(length<lmin) lmin=length;
      if(width>wmax)  wmax=width;
      if(length>lmax) lmax=length;
      if(significance>smax) 
	{
	  smax=significance;
	  emax=excess;
	}

      if(width!=prevw)   dw=width-prevw;
      if(length!=prevl) dl=length-prevl;

      prevw=width;
      prevl=length;
      i++;
    }
  ifun.close();

  cout << "wmin: " << wmin << ", wmax: " << wmax+dw << ", lmin: " << lmin << ", lmax: " << lmax+dl << endl;

  // define histos
  Int_t nx = TMath::Nint((wmax-wmin)/dw)+1;
  Int_t ny = TMath::Nint((lmax-lmin)/dl)+1;
  TH2F *ex  = new TH2F("Excess","Excess",nx,wmin,wmax+dw,ny,lmin,lmax+dl);
  TH2F *sig = new TH2F("Significance","Significance",nx,wmin,wmax+dw,ny,lmin,lmax+dl);

  cout << "nx: " << nx << ", ny: " << ny << endl;
  cout << "dw: " << dw << ", dl: " << dl << endl;

  ifstream ifun2;
  ifun2.open(fname.Data());
  
  // second loop to fill the histograms
  Float_t wopt,lopt,eopt;
  while(ifun2 >> width)
    {
      // read data
      ifun2 >> length;
      ifun2 >> excess;
      ifun2 >> significance;

      cout << "Width: " << width << '\t' << "Length: " << length << '\t' << "Excess: " << excess << '\t' << "Significance: " << significance;
      if(significance==smax)
	{
	  cout << " *****";
	  wopt=width;
	  lopt=length;
	  eopt=excess;
	}
	    
      cout << endl;

      Int_t binx = TMath::Nint(((width-wmin)/dw));
      Int_t biny = TMath::Nint(((length-lmin)/dl));

      //      cout << "binnx: " << binx+1 << ", biny: " << biny+1 << endl;

      ex->SetBinContent(binx+1,biny+1,excess+.00001);
      sig->SetBinContent(binx+1,biny+1,significance+.00001);    
    }
  cout << "Maximal significance: " << smax << " (excess=" << eopt << ") for width="<<wopt<<", length=" << lopt << endl;

  // plot
  TCanvas *c3 = new TCanvas("c3","c3",0,0,500,750);
  c3->Divide(1,2);
  
  c3->cd(1);
  ex->GetXaxis()->SetTitle("Width (degrees)");
  ex->GetYaxis()->SetTitle("Length (degrees)");
  ex->Draw("colz");
  
  c3->cd(2);
  sig->GetXaxis()->SetTitle("Width (degrees)");
  sig->GetYaxis()->SetTitle("Length (degrees)");
  sig->Draw("colz");
}