/////////////////////////////////////////////////////////////////////////////
//
//       This macro makes the two-dimension plot of the False Source Method:
//        root > .x signal.C(TString name, TString psname)
//
//                 Ester Aliu   <aliu@ifae.es>
//                 Oscar Blanch <blanch@ifae.es>
//                 Javier Rico  <jrico@ifae.es>
////////////////////////////////////////////////////////////////////////////
void signal(
TString fname ="hillasCrab/falseSourceCrab20040215.root",
TString psname ="kk.ps")
{
  Float_t excess;
  Float_t significance;
  Float_t array[2];
  
  Int_t nbin2d = 5;
  Int_t ntotal = 2*nbin2d+1;
  
  Int_t nx = 0;
  Int_t ny = 0;
  
  Float_t max = -10000;
  Int_t imax;
  Int_t jmax;
  Float_t A ;
  Float_t B;
  
  TH2F *hsignif = new TH2F("hsignif","Significance", ntotal, -1.0, 1.0, ntotal, -1.0, 1.0);
  TH2F *hexcess = new TH2F("hexcess","Excess", ntotal, -1.0, 1.0, ntotal, -1.0, 1.0);
  
  
  // loop on the histos and compute excess and significance
  for( Int_t i = -nbin2d; i <= nbin2d ; i++)  
    {    
      nx++;
      ny = 0;
      for( Int_t j = -nbin2d; j <= nbin2d ; j++)  
	{	
	  ny++;	  
	  step_signal(i,j,&array[0], fname);
	  
	  excess = array[0];
	  significance = array[1];

	  hexcess->SetBinContent(nx, ny, excess);
	  hsignif->SetBinContent(nx, ny, significance);

	  if ( significance > max)
	    {
	      max = significance;
	      imax = i;
	      jmax = j;
	      A = excess;
	      B = significance;
	    }
	}
    }
  
  cout << "The position of maximum significance is ( " << imax << " , " << jmax << " )" << endl;
  cout << "Excess: " << A << endl;
  cout << "Significance: " << B << endl;
    
  // Plot 
  gROOT->Reset();
  gStyle->SetCanvasColor(0);
  gStyle->SetCanvasBorderMode(0);
  gStyle->SetPadBorderMode(0);
  gStyle->SetFrameBorderMode(0);
  gStyle->SetOptTitle(0);
  gStyle->SetTitleOffset(1.7,"y");
  gStyle->SetPadLeftMargin(0.15);
  gStyle->SetOptStat(kFALSE);
  gStyle->SetStatColor(0);
  gStyle->SetStatBorderSize(1);
  gStyle->SetStatW(0.2);
  gStyle->SetStatH(0.1);
  gStyle->SetStatX(0.9);
  gStyle->SetStatY(0.9);
  

  TPostScript myps(psname,111);
  myps.Range(15,15);
  TCanvas *c2d = new TCanvas("c2d", "Matrices", 0, 0, 800,800);  
  c2d->Divide(2,2);
  gStyle->SetPalette(1); 
  
  c2d->cd(1);
  hexcess->Draw("colz");
  hexcess->SetXTitle("X position (deg)");
  hexcess->SetYTitle("Y position (deg)");
  gPad->Update();    
  
  c2d->cd(2);
  hsignif->Draw("colz");
  hsignif->SetXTitle("X position (deg)");
  hsignif->SetYTitle("Y position (deg)");
  gPad->Update();

  c2d->cd(3);
  hexcess->Draw("lego2");
  hexcess->SetXTitle("X position (deg)");
  hexcess->SetYTitle("Y position (deg)");
  gPad->Update();    
  
  c2d->cd(4);
  hsignif->Draw("lego2");
  hsignif->SetXTitle("X position (deg)");
  hsignif->SetYTitle("Y position (deg)");
  gPad->Update();    

  myps.Close();
}
  
  
void step_signal(Int_t i, Int_t j, Float_t *array, TString fname) 
{
  cout << "Bin (" << i << "," << j << ")" << endl;
  TH1F *onhisto;
  TH1F *offhisto;

  // open file and access to the histograms  
  TFile *f = new TFile(fname);
  Char_t name[20]; 
  Char_t title[50]; 

  // histo name
  sprintf(name,"hOnAlpha[%d][%d]", i, j);
  sprintf(title,"Alpha-Plot(On data) (%d ,%d)", i, j);
  onhisto = (TH1F*)gDirectory->Get(name);   
  
  sprintf(name,"hOffAlpha[%d][%d]", i, j);
  sprintf(title,"Alpha-Plot(Off data) (%d ,%d)", i, j);
  offhisto = (TH1F*)gDirectory->Get(name);   
  
  // line/normalization
  const Int_t nbins = 18;
  const Int_t inibin = 20./90.*nbins+1;
  Float_t level=0;
  Float_t leveloff=0;
  Float_t levelofferror=0;
  for(Int_t ibin = inibin; ibin<=nbins;ibin++)
    {
      level+=onhisto->GetBinContent(ibin);
      leveloff+=offhisto->GetBinContent(ibin);
    }
  level/=(nbins-inibin+1);
  leveloff/=(nbins-inibin+1);
  
  // normalize on/off
  offhisto->Sumw2(); // needed to compute correct errors after normalization
  const Float_t norm = level/leveloff;
  cout << "Normalizing by factor " << norm <<endl;
  offhisto->Scale(norm);    
  
  // significance:
  Float_t sig=0,bg=0,esig=0,ebg=0;
  Float_t significance=0;
  const Int_t signbins = inibin-1;
  for(Int_t ibin = 1; ibin<=signbins;ibin++)
    {
      sig  += onhisto->GetBinContent(ibin);
      esig += onhisto->GetBinError(ibin)*onhisto->GetBinError(ibin);
      bg   += offhisto->GetBinContent(ibin);
      ebg  += offhisto->GetBinError(ibin)*offhisto->GetBinError(ibin);
    }
  Float_t error= TMath::Sqrt(esig+ebg);
  Float_t excess = sig-bg;
  Float_t significance = excess/error;
  
  cout << "Excess:       " << excess << endl;
  cout << "Significance: " << significance << endl;
  
  array[0] = excess;
  array[1] = significance;
}