/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without expressed
! * or implied warranty.
! *
!
!
!   Author(s): Robert Wagner 9/2002 <mailto:rw@rwagner.de>
!   Author(s): Robert Wagner 3/2003 <mailto:rw@rwagner.de>
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//                                                                          //
//  MHOnSubtraction                                                         //
//                                                                          //
//                                                                          //
//  extracts the gamma signal from a pure ON-signal given in an             //
//  ALPHA-Energy-Theta histogram. The class will take this histogram from   //
//  the parameter list and will provide result histograms in the            //
//  parameter list.                                                         //
//                                                                          //
//  This class still is work in progress.                                   //
//                                                                          //
//                                                                          //
//////////////////////////////////////////////////////////////////////////////

// This part of MARS is code still evolving. Please do not change the code 
// without prior feedback by the author.

#include "MHOnSubtraction.h"

#include <TPaveText.h>
#include <TPaveLabel.h>
#include <TF1.h>
#include <TLegend.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <TGraph.h>

#include "MBinning.h"
#include "MParList.h"
#include "MHArray.h"
#include "MH3.h"

#include "MHAlphaEnergyTheta.h"
#include "MHAlphaEnergyTime.h"

#include "MLog.h"
#include "MLogManip.h"

ClassImp(MHOnSubtraction);

using namespace std;

// --------------------------------------------------------------------------
//
// Default Constructor.
//
MHOnSubtraction::MHOnSubtraction(const char *name, const char *title) : fMaxSignif(0),	fMaxRedChiSq(0), fSlope(20.0)
{ 
  fName  = name  ? name  : "MHOnSubtraction";
  fTitle = title ? title : "Extracts Gamma signal from pure ON data";
  fHistogramType = "Theta";
  fChiSquareHisto=0;
  fSignificanceHisto=0;
  fSummedAlphaPlots=0;
  fThetaBin=0;
  
}


// --------------------------------------------------------------------------
//
// Destructor.
//
MHOnSubtraction::~MHOnSubtraction()
{ 
  if (fChiSquareHisto) delete fChiSquareHisto;
  if (fSignificanceHisto) delete fSignificanceHisto;
  if (fSummedAlphaPlots) delete fSummedAlphaPlots;
}

// -----------------------------------------------------------------------
//
// Calculate Significance according to Li and Ma, ApJ 272 (1983) 317, eq
// (17). n_{On}, n_{Off} and the ratio of On-times for On and Off 
// measurements have to be provided.
//
// This function underestimates the true significance for it does not take
// into account errors on the event numbers. A more accurate variation wil
// be provided soon
//
Double_t MHOnSubtraction::CalcSignificance(Double_t nOn, Double_t nOff, Double_t theta)
{ 
  if (nOn<=0)
    *fLog << warn << "Got " << nOn << " total events, " << flush;
  if (nOff<=0)
    *fLog << warn << "Got " << nOff << " background events, " << flush;
  if (nOn<=0 || nOff<=0.0) {
    *fLog << warn << "returning significance of 0.0" << endl;
    return 0.0;
  }
  return sqrt(2*(nOn*log((1+theta)*nOn/(theta*(nOn+nOff)))
		 +nOff*log((1+theta)*(nOff/(nOn+nOff)))));
}

// -----------------------------------------------------------------------
//
// This function takes a first look at a given ALPHA distribution
// and determines if a fit is applicable.
//
// Fits a given TH1 containing an ALPHA distribution with a combined
// gaussian plus polynomial of third grade and returns the fitted 
// function. By signalRegionFactor the width of the region in which
// signal extraction is to be performed can be specified in units of
// sigma of the gaussian which is fitted to the distribution (default: 
// 3.5). Accordingly, FitHistogram returns the range in which it suggests 
// signal extraction (lowerBin < ALPHA < upperBin). An estimated 
// significance of the signal is also provided. draw specifies whether
// FitHistogram should draw the distribution.
//
Bool_t MHOnSubtraction::FitHistogram(TH1 &alphaHisto, Double_t &sigLiMa,
			       Double_t &lowerBin, Double_t &upperBin,
			       Float_t signalRegionFactor, const Bool_t draw,
			       TString funcName)
{
  if (alphaHisto.GetEntries() == 0) {
    *fLog << warn << "Histogram contains no entries. Returning. " << endl;
    lowerBin=0;
    upperBin=0;
    sigLiMa=0;    
    if (draw) {
      TPaveLabel* lab = new TPaveLabel(0.1,0.2,0.9,0.8,"(no events)");
      lab->SetFillStyle(0);
      lab->SetBorderSize(0);
      lab->Draw();  
    }
    return kFALSE;
  }
  
  //cosmetics
  alphaHisto.GetXaxis()->SetTitle("ALPHA");
  alphaHisto.GetYaxis()->SetTitle("events");

  Double_t outerEvents = 0.0;
  Double_t innerEvents = 0.0;
  Int_t outerBins = 0;
  Int_t innerBins = 0;
  Int_t outerBinsZero = 0;
  Int_t innerBinsZero = 0;

  for (Int_t alphaBin = 1; alphaBin < alphaHisto.GetNbinsX(); alphaBin++) {
    if ((alphaHisto.GetBinLowEdge(alphaBin) >= -35.)&& //inner Region
  	(alphaHisto.GetBinLowEdge(alphaBin+1) <= 35.0)) {
      innerEvents+=alphaHisto.GetBinContent(alphaBin);
      innerBins++;
      if (alphaHisto.GetBinContent(alphaBin)==0.0)
  	innerBinsZero++;     
    } else {
      if ((alphaHisto.GetBinLowEdge(alphaBin) >= -89.5)&& //outer Region
	  (alphaHisto.GetBinLowEdge(alphaBin+1) <=89.5)) {
	outerEvents+=alphaHisto.GetBinContent(alphaBin);
	outerBins++;
	if (alphaHisto.GetBinContent(alphaBin)==0.0)
	  outerBinsZero++;     
      }
    }
  }

   *fLog << dbg << "Plot contains " << 
     outerEvents << " outer ev (" << outerEvents/outerBins << "/bin), " <<
     innerEvents << " inner ev (" << innerEvents/innerBins << "/bin) " << endl;

  if ((outerBinsZero!=0) || (innerBinsZero != 0))
    *fLog << warn << "There are ";
  if (outerBinsZero != 0)
    *fLog << dbg << outerBinsZero << " empty outer bins ";
  if (innerBinsZero != 0)
    *fLog << dbg <<innerBinsZero << " empty inner bins ";
  if ((outerBinsZero!=0) || (innerBinsZero != 0))
    *fLog << endl;

  if (outerEvents == 0.0) {
     *fLog << warn << "No events in OFF region. Assuming background = 0" 
	   << endl;
     TF1 *po = new TF1("pol0"+funcName,"pol0(0)",-89.5,89.5);
     po->SetLineWidth(2);
     po->SetLineColor(2);
     po->SetParNames("c");
     po->SetParameter(0,0.0);
     alphaHisto.GetListOfFunctions()->Add(po);
     alphaHisto.SetLineColor(97);
     if (draw) {
       alphaHisto.DrawCopy(); //rwagner
       po->Draw("SAME");
     }
     sigLiMa = 0;
     lowerBin = 1;
     upperBin = alphaHisto.GetNbinsX()-1;
     signalRegionFactor = 0;
    
     return kFALSE; //No gaus fit applied
  }

  if (outerEvents/outerBins < 2.5) {
    *fLog << warn << "Only " << /*setprecision(2) <<*/ outerEvents/outerBins 
	  << " events/bin in OFF region: "
	  << "Assuming this as background." << endl;
   
     TF1 *po = new TF1("pol0"+funcName,"pol0(0)",-89.5,89.5);
     po->SetLineWidth(2);
     po->SetLineColor(94);
     po->SetParNames("c");
     po->SetParameter(0,outerEvents/outerBins);
     alphaHisto.GetListOfFunctions()->Add(po);    
     if (draw) {
       alphaHisto.DrawCopy(); //rwagner
       po->Draw("SAME");
     }
    
     Int_t centerBin = alphaHisto.GetXaxis()->FindBin((Double_t)0.0);
     Int_t maxBin    = centerBin > alphaHisto.GetNbinsX() - centerBin ? 
       alphaHisto.GetNbinsX()- centerBin : centerBin;
     
     Int_t lowerSignalEdge = centerBin;
     for (Int_t i=3; i < maxBin; i++) {
       if ((Double_t)alphaHisto.GetBinContent(centerBin-i) 
	   < outerEvents/outerBins) {
	 lowerSignalEdge = centerBin-i+1;
	 break;
       }
     }
     if (centerBin<lowerSignalEdge) lowerSignalEdge=centerBin;

     Int_t upperSignalEdge = centerBin;
     for (Int_t i=3; i < maxBin; i++) {
       if ((Double_t)alphaHisto.GetBinContent(centerBin+i) 
	   <= outerEvents/outerBins) {
	 upperSignalEdge=centerBin+i-1;
	 break;
       } 
     }
     if (centerBin>upperSignalEdge) upperSignalEdge=centerBin;

     Double_t nOnInt = 0;
     for (Int_t i=1; i < alphaHisto.GetNbinsX(); i++)
       nOnInt += alphaHisto.GetBinContent(i) - outerEvents/outerBins;
     
     Double_t nOffInt = (upperSignalEdge - lowerSignalEdge + 1) 
       * (outerEvents/outerBins);
     
     sigLiMa = CalcSignificance(nOnInt, nOffInt, 1);

     if (sigLiMa < 3) 
       alphaHisto.SetLineColor(2);   

   
     *fLog << inf << "Estimated significance is " << sigLiMa 
	   << " sigma " << endl;   
     *fLog << inf << "Signal region is " 
	   << lowerBin << " < ALPHA < " << upperBin 
	   << " (Most likely you wanna ignore this)" << endl;

     return kFALSE; //No gaus fit applied    
  }
  
  //fit combined gaus+pol3 to data
  TF1 *gp = new TF1("gauspol3"+funcName,"gaus(0)+pol3(3)",-89.5,89.5);
  gp->SetLineWidth(2);
  gp->SetLineColor(2);
  gp->SetParNames("Excess","A","sigma","a","b","c","d");
  gp->SetParLimits(0, 200,2000); 
  gp->SetParLimits(1, -4.,4.); 
  gp->SetParLimits(2, 2.,20.);
  // gp->SetParameter(6,  0.0000); // include for slope(0)=0 constrain
  TString gpDrawOptions = draw ? "RIQ" : "RIQ0";
  gpDrawOptions = "RIQ0"; // rwagner
  alphaHisto.Fit("gauspol3"+funcName,gpDrawOptions);
  alphaHisto.DrawCopy(); //rwagner
  alphaHisto.SetDirectory(NULL); //rwagner
  
  Double_t gausMean  = gp->GetParameter(1);
  Double_t gausSigma = gp->GetParameter(2);

//   TF1 *p3 = new TF1("p3"+funcName,"pol3(0)",-signalRegionFactor*gausSigma+gausMean,
// 		    signalRegionFactor*gausSigma+gausMean);
  TF1 *p3 = new TF1("p3"+funcName,"pol3(0)",-100,
		    100);
  p3->SetParameters(gp->GetParameter(3),gp->GetParameter(4),
		    gp->GetParameter(5),gp->GetParameter(6));
  // p3->SetLineStyle(2);
  p3->SetLineColor(4);
  p3->SetLineWidth(1);
  if (draw) p3->Draw("SAME");   
   
  TF1 *ga = new TF1("ga"+funcName,"gaus(0)",-40,40);
  ga->SetParameters(gp->GetParameter(0),gp->GetParameter(1),gp->GetParameter(2));
  ga->SetLineColor(93);
  ga->SetLineWidth(2);
  //  if (draw) ga->Draw("SAME");

  // identify the signal region: signalRegionFactor*gausSigma
  // this allows us to
  // (1) determine n_{ON}, n_{OFF}
  Double_t scalingFactor = 
    (alphaHisto.GetXaxis()->GetBinLowEdge(alphaHisto.GetNbinsX()+1)
     - alphaHisto.GetXaxis()->GetBinLowEdge(1)) / alphaHisto.GetNbinsX();
  
  Double_t nOnInt  = (gp->Integral(-signalRegionFactor*gausSigma+gausMean, 
				 signalRegionFactor*gausSigma+gausMean) / scalingFactor);
  Double_t nOffInt = (p3->Integral(-signalRegionFactor*gausSigma+gausMean, 
				 signalRegionFactor*gausSigma+gausMean) / scalingFactor);

  // (2) determine the signal region from fit in degrees
  // we do it a bit more complicated: assuming that the binning in all 
  // histograms is the same, we want to be sure that summing up is always
  // done over the same bins.

  lowerBin = alphaHisto.GetXaxis()->FindBin(-signalRegionFactor*gausSigma+gausMean);
  upperBin = alphaHisto.GetXaxis()->FindBin( signalRegionFactor*gausSigma+gausMean);

  lowerBin = alphaHisto.GetBinLowEdge((Int_t)lowerBin);
  upperBin = alphaHisto.GetBinLowEdge((Int_t)upperBin)+alphaHisto.GetBinWidth((Int_t)upperBin);
  
  sigLiMa = CalcSignificance(nOnInt, nOffInt, 1);
//   if (sigLiMa < 3)
//     alphaHisto.SetLineColor(2);   


  *fLog << inf << "Fit estimates significance to be " 
	<< sigLiMa << " sigma " << endl; 

  *fLog << inf << "Fit yields signal region to be " 
	<< lowerBin << " < ALPHA < " << upperBin
	<< " (Chisquare/dof=" << gp->GetChisquare()/gp->GetNDF() 
	<< ", prob="  <<  gp->GetProb() << ")" << endl;

  return kTRUE; //returning gaus fit
}





// -----------------------------------------------------------------------
//
// Does the actual extraction of the gamma signal. For performance 
// reasons, fits already done by MHOnSubtraction::FitHistogram are used 
// and not redone.
// From it, a polynomial function for the background is evaluated and the
// gamma signal is extracted in the region given by lowerBin < ALPHA < 
// upperBin.
// Significance of the signal is also provided. draw specifies whether
// FitHistogram should draw the distribution.
//
Bool_t MHOnSubtraction::ExtractSignal(TH1 &alphaHisto, Double_t &sigLiMa,
   Double_t &lowerBin, Double_t &upperBin,
   Double_t &gammaSignal, Double_t &errorGammaSignal,
   Double_t &off, Double_t &errorOff,
   Float_t signalRegionFactor, const Bool_t draw, TString funcName,
   TPad *drawPad, Int_t drawBase)
{
  TF1 *gausPol = alphaHisto.GetFunction("gauspol3"+funcName);
  TF1 *pol = alphaHisto.GetFunction("pol0"+funcName);
  
  if (!gausPol && !pol) {
    *fLog << err << "Fatal: ALPHA histogram has no gauspol or pol "  
	  << " fit attached to it." << endl;
    TPaveLabel* lab = new TPaveLabel(0.1,0.2,0.9,0.8,"(no fit)");
    lab->SetFillStyle(3000);
    lab->SetBorderSize(0);
    lab->DrawClone();  
    lab->SetBit(kCanDelete);
    return kFALSE;
  } 

  TF1* po;
  po = NULL;

  if (gausPol) {
    Double_t gausMean  = gausPol->GetParameter(1);
    Double_t gausSigma = gausPol->GetParameter(2);
    po = new TF1("po"+funcName,"pol3(0)",
		  -signalRegionFactor*gausSigma+gausMean,
		  signalRegionFactor*gausSigma+gausMean);
    po->SetParameters(gausPol->GetParameter(3),gausPol->GetParameter(4),
		      gausPol->GetParameter(5),gausPol->GetParameter(6));

    TF1 *ga = new TF1("ga"+funcName,"gaus(0)",-40,40);
    ga->SetParameters(gausPol->GetParameter(0),gausPol->GetParameter(1),
		      gausPol->GetParameter(2));
    
    if (draw) {
      alphaHisto.Draw();
      gausPol->Draw("SAME"); //Maybe not even necessary?
      
      po->SetLineColor(4);
      po->SetLineWidth(2);
      po->Draw("SAME");   
      
      ga->SetLineColor(93);
      ga->SetLineWidth(2);    
      ga->Draw("SAME");   
      
      char legendTitle[80];
      sprintf(legendTitle, "Signal region: %2.1f < #alpha < %2.1f", 
	      lowerBin, upperBin);
      
      TLegend *legend = new TLegend(0.13, 0.52, 0.47, 0.72, legendTitle);
      
      legend->SetBorderSize(0);
      legend->SetFillColor(10);
      legend->SetFillStyle(0);
      legend->AddEntry(gausPol, "combined N_{on}","l");
      legend->AddEntry(po,"polynomial N_{bg} Signal region","l");
      legend->AddEntry(ga, "putative gaussian N_{S}","l");
      legend->Draw();
    }
  } // gausPol


  if (pol) {
    po = pol;
    
    if (draw) {
      alphaHisto.Draw();
      
      po->SetLineColor(6);
      po->SetLineWidth(2);
      po->Draw("SAME");   
            
      char legendTitle[80];
      sprintf(legendTitle, "Signal region: %2.1f < #alpha < %2.1f", 
	      lowerBin, upperBin);
      
      TLegend *legend = new TLegend(0.13, 0.52, 0.47, 0.72, legendTitle);
      
      legend->SetBorderSize(0);
      legend->SetFillColor(10);
      legend->SetFillStyle(0);
      legend->AddEntry(po,"polynomial N_{bg} Signal region","l");
      legend->Draw();
    }
  } // pol
    
  Double_t nOn = 0;
  Double_t eNOn = 0;

  Int_t binNumberLow = alphaHisto.GetXaxis()->FindBin(lowerBin);
  Int_t binNumberHigh = alphaHisto.GetXaxis()->FindBin(upperBin);
  
  for (Int_t bin=binNumberLow; bin<binNumberHigh+1; bin++) {
    nOn += alphaHisto.GetBinContent(bin);
    eNOn += alphaHisto.GetBinError(bin) * alphaHisto.GetBinError(bin);
  } //for bin
  eNOn = sqrt(eNOn);
     
  // Evaluate background
  
  Double_t nOff = 0;
  Double_t eNOff = 0;
  for (Int_t bin=binNumberLow; bin<binNumberHigh+1; bin++) {
    Double_t x = .5*(alphaHisto.GetBinLowEdge(bin)+alphaHisto.GetBinLowEdge(bin+1));
    Double_t offEvts = po->Eval(x);
    //cout << bin << ": " << offEvts << endl;
    nOff += offEvts;
  } //for bin
  eNOff = sqrt(fabs(nOff));
  
  if (nOn==0)   // there should not be a negative number of signal events
    nOff=0;
  
  if (nOff<0) { // there should not be a negative number of off events
    nOff=0;
    eNOff=0;
  }
    
  *fLog << inf << "nEvts = " << nOn << "+-" << eNOn << ", ";

  off = nOff;              errorOff = eNOff;
  gammaSignal = nOn-nOff;  errorGammaSignal = sqrt(eNOn*eNOn + eNOff*eNOff); 

  *fLog << inf << "nBack = " << nOff << "+-" << eNOff << ", ";
  *fLog << inf << "nSig = " << gammaSignal << "+-" << errorGammaSignal << endl;

  sigLiMa = CalcSignificance(nOn, nOff, 1);  
  //  Double_t sigFit=(ga->GetParameter(1))/(ga->GetParError(1)); //Mean / sigMean
  
  *fLog << inf << "Significance: "<<sigLiMa<<" sigma "<<endl;

  if (draw) { 
    TPaveText *pt = new TPaveText(0.11,.74,.57,.88,"NDC ");
    char tx[60];
    sprintf(tx, "Excess: %2.2f #pm %2.2f", nOn-nOff, sqrt(eNOn*eNOn + eNOff*eNOff));
    pt->AddText(tx);
    sprintf(tx, "Off:    %2.2f #pm %2.2f", nOff, eNOff);
    pt->AddText(tx);
    sprintf(tx, "Significance: %2.2f   #sigma", sigLiMa);
    pt->AddText(tx);
    pt->SetFillStyle(0);
    pt->SetBorderSize(0);
    pt->SetTextAlign(12);
    pt->Draw("");
  }  
  if (draw||drawPad) {

    // Plot significance vs alpha_cut
    
    Int_t centerBin = alphaHisto.GetXaxis()->FindBin((Double_t)0.0);
    Int_t maxBin    = centerBin > alphaHisto.GetNbinsX()- centerBin ? 
      alphaHisto.GetNbinsX()- centerBin : centerBin;
         
    const Int_t totalBins = centerBin;
    Float_t alpha[totalBins];
    Float_t signi[totalBins];
    Float_t maxSigni = 0;
    
    for (Int_t i=0; i < maxBin; i++) {
      Double_t nOn = 0;  Double_t eNOn = 0;
      Double_t nOff = 0; Double_t eNOff = 0;
      for (Int_t bin=centerBin-i; bin<centerBin+i+1; bin++) {
	nOn += alphaHisto.GetBinContent(bin);
	eNOn += alphaHisto.GetBinError(bin) * alphaHisto.GetBinError(bin);
	Double_t x = .5*(alphaHisto.GetBinLowEdge(bin)+alphaHisto.GetBinLowEdge(bin+1));
	Double_t offEvts = po->Eval(x);
	nOff += offEvts;
      } //for bin
      eNOn = sqrt(eNOn);
      eNOff = sqrt(nOff);  
      alpha[i] = 
	(alphaHisto.GetXaxis()->GetBinLowEdge(centerBin+i+1)-
	 alphaHisto.GetXaxis()->GetBinLowEdge(centerBin-i))/2;
      signi[i] = CalcSignificance(nOn, nOff, 1);  
      maxSigni = maxSigni > signi[i] ? maxSigni : signi[i];
    }

    if (!drawPad) {
      TCanvas *c3 = new TCanvas("c3"+funcName, "Significance vs ALPHA cut", 800,600);
      c3->cd();
    }

    if (drawPad) 
      drawPad->cd(drawBase-1);
    
    TGraph* gr = new TGraph(totalBins-2,alpha,signi);
    TString drawOpt = "L";
    
    if (draw || (drawPad && fSigniPlotIndex == 0))
      drawOpt += "A";
    
    gr->Draw(drawOpt);
    if (drawPad && fSigniPlotIndex == 0) {
      gr->GetXaxis()->SetTitle("|#alpha_{max}|");
      gr->GetYaxis()->SetTitle("Significance");
    }
    gr->SetMarkerStyle(2);
    gr->SetMarkerSize(1);
    gr->SetMarkerColor(4+fSigniPlotIndex);
    gr->SetLineColor(4+fSigniPlotIndex);
    gr->SetLineWidth(1);
    gr->SetTitle("Significance vs ALPHA integration range");
    gr->Draw("L");

    fSigniPlotIndex++;      
    
  } //draw

  return kTRUE;
}

// -----------------------------------------------------------------------
//
// Extraction of Gamma signal is performed on a MHAlphaEnergyTheta or 
// MHAlphaEnergyTime object expected in the ParList.
//
Bool_t MHOnSubtraction::Calc(MParList *parlist, const Bool_t Draw)
{
  //Find source histograms
  fHistogramType = "Theta";
  MHAlphaEnergyTheta* mhisttheta = (MHAlphaEnergyTheta*)parlist->FindObject("MHAlphaEnergyTheta"); 
  if (mhisttheta) return Calc((TH3D*)mhisttheta->GetHist(), parlist, Draw);

  fHistogramType = "Time";
  MHAlphaEnergyTime* mhisttime = (MHAlphaEnergyTime*)parlist->FindObject("MHAlphaEnergyTime");  
  if (mhisttime) return Calc((TH3D*)mhisttime->GetHist(), parlist, Draw);

  *fLog << err << "No MHAlphaEnergyTheta type object found in the parameter list. Aborting." << endl;
  return kFALSE; 
}

// -----------------------------------------------------------------------
//
// Extraction of Gamma signal is performed on a TH3D histogram in
// ALPHA, Energy and Theta.
//
Bool_t MHOnSubtraction::CalcAET(TH3D *aetHisto, MParList *parlist, const Bool_t Draw)
{
  // Analyze aetHisto
  // -------------------------------------
  Int_t alphaBins =  aetHisto->GetNbinsX();  // # of alpha bins 
  Int_t energyBins = aetHisto->GetNbinsY();  
  Int_t thetaBins =  aetHisto->GetNbinsZ();  
  
  // We output an array of histograms to the parlist.
  // Create a template for such a histogram, needed by MHArray
  // -------------------------------------
  MBinning *binsfft = new MBinning("BinningMH3X");
  binsfft->SetEdgesLog(energyBins,aetHisto->GetYaxis()->GetBinLowEdge(1),
		aetHisto->GetYaxis()->GetBinLowEdge(energyBins+1));
  parlist->AddToList(binsfft);  
  MH3 *energyHistogram = new MH3(1); // The template histogram in energy
                                     // for a given theta value
  energyHistogram->SetName("MHOnSubtractionEnergyHistogramTemplate");
  parlist->AddToList(energyHistogram);

  fThetaHistoArray = new MHArray("MHOnSubtractionEnergyHistogramTemplate", kTRUE, 
				 "MHOnSubtractionGammaSignalArray", 
				 "Array of histograms in energy bins");
  fThetaHistoArray->SetupFill(parlist);
  parlist->AddToList(fThetaHistoArray);
  
  // Canvases---direct debug output for the time being
  // -------------------------------------
  TCanvas *c3 = new TCanvas("c3", "Plots by MHOnSubtraction::ExtractSignal", 800,600);
  cout << thetaBins << " x " << energyBins << endl;
  c3->Divide(thetaBins,energyBins);

  TCanvas *c4a = new TCanvas("c4a", "Energy distributions for different ZA", 800,600);
  
  TH1D* histalphaon[energyBins*thetaBins]; // ALPHA histograms

  fChiSquareHisto = new TH1D("fChiSquareHisto", "#chi^{2}/d.o.f. of fits", 50, 0, 5);
  fSignificanceHisto = new TH1D("fSignificanceHisto", "Significances", 41, -0.5, 40.5);
  fSummedAlphaPlots = new TH1D("fSummedAlphaPlots", "Cumulative Alpha", 
			       alphaBins,aetHisto->GetXaxis()->GetBinLowEdge(1),
			       aetHisto->GetXaxis()->GetBinLowEdge(alphaBins+1));
  
  // -------------------------------------

  fThetaLegend = new TLegend(0.83, 0.07, 0.98, 0.42, "Energy Distributions");  
  fThetaLegend->SetBorderSize(1);
  
  Double_t overallSigLiMa = 0;

  for (Int_t thetaBin = 1; thetaBin < thetaBins+1; thetaBin++) {
    
    char hname[80];  
    sprintf(hname, "Energy distribution for %s bin %d", fHistogramType.Data(), 
	    thetaBin);
    *fLog << all << "Calculating " << hname << endl;

    Double_t minLowerBin = -10; // minimum integration range
    Double_t maxUpperBin =  10; // minimum integration range
    Double_t maxAlpha =  70; // maximum integration range

    Double_t sumSigLiMa = 0;
    
    // This loop just fixes the integration range
    // And alerts when no significant excess is found.

    for (Int_t energyBin = 1; energyBin < energyBins+1; energyBin++) {

      sprintf(hname, "histalphaon%d", (thetaBin-1)*(energyBins)+energyBin);
      histalphaon[(thetaBin-1)*(energyBins)+energyBin] = 
	new TH1D(hname, "ALPHA histogram",
		 alphaBins,aetHisto->GetXaxis()->GetBinLowEdge(1),
		 aetHisto->GetXaxis()->GetBinLowEdge(alphaBins+1));
      histalphaon[(thetaBin-1)*(energyBins)+energyBin]->Sumw2();    
      histalphaon[(thetaBin-1)*(energyBins)+energyBin]->SetTitle(hname);
     
      sprintf(hname,"ALPHA distribution for E bin %d, theta bin %d",
	      energyBin, thetaBin);
      *fLog << inf << hname << endl;

      // now we fill one histogram for one particular set
      // of Energy and Theta...
      
      if (aetHisto->InheritsFrom("TH3D")||aetHisto->InheritsFrom("TH2D")) {
	aetHisto->GetYaxis()->SetRange(energyBin,energyBin); // E
	if (aetHisto->InheritsFrom("TH3D"))
	  aetHisto->GetZaxis()->SetRange(thetaBin,thetaBin); // theta
	histalphaon[(thetaBin-1)*(energyBins)+energyBin] = (TH1D*)aetHisto->Project3D("xe");
      } else {
	histalphaon[(thetaBin-1)*(energyBins)+energyBin] = (TH1D*)aetHisto;
      }

      *fLog << inf << "This histogram has " 
	    << histalphaon[(thetaBin-1)*(energyBins)+energyBin]->GetEntries()
	    << " entries" << endl;

      sprintf(hname, "histalphaon%d", (thetaBin-1)*(energyBins)+energyBin);
      histalphaon[(thetaBin-1)*(energyBins)+energyBin]->SetName(hname);
      
      // Test: Find signal region and significance in histogram
      Double_t lowerBin, upperBin, sSigLiMa;
      FitHistogram(*histalphaon[(thetaBin-1)*(energyBins)+energyBin],
		   sSigLiMa, lowerBin, upperBin, (Float_t)3.5, kFALSE);

      Double_t redChiSq = histalphaon[(thetaBin-1)*(energyBins)+energyBin]->GetFunction("gauspol3")->GetChisquare()/histalphaon[(thetaBin-1)*(energyBins)+energyBin]->GetFunction("gauspol3")->GetNDF();

      fChiSquareHisto->Fill(redChiSq);
      fMaxRedChiSq = redChiSq > fMaxRedChiSq ? redChiSq : fMaxRedChiSq;

      // histalphaon[(thetaBin-1)*(energyBins)+energyBin]->GetFunction("gauspol3")->GetProb() <<endl;

      fSummedAlphaPlots->Add(histalphaon[(thetaBin-1)*(energyBins)+energyBin], 1);
      
      if (sSigLiMa < 3)
	*fLog << inf << "No significant excess for this bin (sigma="<< sSigLiMa <<")!"<<endl;
      sumSigLiMa+=sSigLiMa;
   
      // Evaluate integration range

      lowerBin = (lowerBin < -maxAlpha) ? -maxAlpha : lowerBin;
      minLowerBin = (lowerBin < minLowerBin) ? lowerBin : minLowerBin;

      upperBin = (upperBin > maxAlpha) ? maxAlpha : upperBin;
      maxUpperBin = (upperBin > maxUpperBin) ? upperBin : maxUpperBin;

    } //energyBin
    
    *fLog << inf << "Integration range for this " << fHistogramType 
	  << " value will be " << minLowerBin << " < ALPHA < "
	  << maxUpperBin << endl;
    
    *fLog << inf << "Summed estd. significance for this " << fHistogramType 
	  << " value is " << sumSigLiMa << endl;

    // Create Histogram in Energy for this Theta value
    cout << "STARTIGN REAL CALC1 "<< endl;
    fThetaHistoArray->AddHistogram();
    cout << "STARTIGN REAL CALC1a "<< endl;
    fThetaHistoArray->SetIndex(thetaBin-1);
    
    MH3 &mThetaHisto = *(MH3*)(fThetaHistoArray->GetH());
    mThetaHisto.Print();
   
    TH1D &thetaHisto = (TH1D&)(mThetaHisto.GetHist());
    sprintf(hname,"Energy distribution for theta bin %d", thetaBin);
    thetaHisto.SetTitle(hname);
    thetaHisto.SetEntries(0);
      
    // we have a rough idea of what is going on in the ALPHA plot...
    // So now we can indeed extract the signal.
         
//     cout << "STARTIGN REAL CALC "<< endl;

    sumSigLiMa = 0;
    for (Int_t energyBin = 1; energyBin < energyBins+1; energyBin++) {

      sprintf(hname,"ALPHA distribution for E bin %d, theta bin %d",
	      energyBin, thetaBin);
      *fLog << inf << hname << endl;
  
      Double_t gammaSignal, errorGammaSignal, off, errorOff, sigLiMa;

      c3->cd((thetaBin-1)*(energyBins)+energyBin);
  
      ExtractSignal(*histalphaon[(thetaBin-1)*(energyBins)+energyBin],
		    sigLiMa, minLowerBin, maxUpperBin, 
		    gammaSignal, errorGammaSignal, 
		    off, errorOff, (Float_t)3, kTRUE);
      
      fSignificanceHisto->Fill(sigLiMa);
      fMaxSignif = sigLiMa > fMaxSignif ? sigLiMa : fMaxSignif;

      thetaHisto.SetBinContent(energyBin, gammaSignal*TMath::Exp(7-energyBin));
      thetaHisto.SetBinError(energyBin, errorGammaSignal);
      thetaHisto.SetEntries(thetaHisto.GetEntries()+gammaSignal);
     
      sumSigLiMa += sigLiMa;
      
    }//energyBin

    *fLog << inf << "Summed significance for this " << fHistogramType 
	  << " value is " << sumSigLiMa << endl;

    //fit exponential function to data
    TF1* e = new TF1("expF","expo",0,5);
    e->SetLineWidth(3);
    e->SetLineColor(thetaBin);
    // e->SetParLimits(1, -0.5,3); //limits on slope
    
     if (fSlope!=20.0) {
       e->SetParameter(1, fSlope);
       *fLog << inf << "Fixing slope to " << e->GetParameter(1) << endl;
     }

    TString eDrawOptions = Draw ? "RIQ" : "RIQ0";
    cout << "doing the fit" << endl;
    thetaHisto.Fit("expF");
    
    Double_t expoSlope = e->GetParameter(1);
    
    *fLog << inf << "Determined slope for energy distribution is " 
	  << expoSlope << endl;

    Double_t integralEvts =
      e->Integral(aetHisto->GetYaxis()->GetBinLowEdge(1), 
		  aetHisto->GetYaxis()->GetBinLowEdge(energyBins+1));
    
    *fLog << inf << "Integral in E range [" <<
      aetHisto->GetYaxis()->GetBinLowEdge(1) << ";" <<
      aetHisto->GetYaxis()->GetBinLowEdge(energyBins+1) << "] is " << 
      integralEvts << endl;
      
    // Plot Energy histogram
  
    c4a->cd(0);
    thetaHisto.SetLineColor(thetaBin);
    if (thetaBin==1) {      
      thetaHisto.Draw();
    } else {
      thetaHisto.Draw("SAME");
    }

    sprintf(hname,"Theta bin %d",thetaBin);
    fThetaLegend->AddEntry(&thetaHisto, hname, "l");

    overallSigLiMa += sumSigLiMa;

  } //thetaBin

  fThetaLegend->Draw();

  Double_t sigLiMa, lowerBin, upperBin;    
  FitHistogram(*fSummedAlphaPlots, sigLiMa, lowerBin, upperBin);
  fSummedAlphaPlots->SetTitle("Cumulative ALPHA distribution");

  *fLog << inf << "Summed overall significance is " << overallSigLiMa << endl;

  *fLog << inf << "Setting range for Significance histo " << fMaxSignif << endl;
  // fSignificanceHisto->GetXaxis()->SetRange(0,fMaxSignif+1);

  // *fLog << inf << "Setting range for chisq histo " << fMaxRedChiSq << endl;
  // fChiSquareHisto->GetXaxis()->SetRange(0,fMaxRedChiSq+5);
  
  return kTRUE;
}








// -----------------------------------------------------------------------
//
// Extraction of Gamma signal is performed on a TH3D histogram in
// ALPHA, Energy and Theta. Output to parameter list: TH2 histograms
// in energy and Theta.
//
Bool_t MHOnSubtraction::Calc(TH3 *aetHisto, MParList *parlist, const Bool_t Draw)
{
  // Analyze aetHisto
  // -------------------------------------
  Int_t energyBins = aetHisto->GetNbinsY();  
  Int_t thetaBins =  aetHisto->GetNbinsZ();  

  *fLog << "Total events:     " << aetHisto->GetEntries() << endl;
  *fLog << "Total energy bins: " << energyBins << endl;
  *fLog << "Total theta bins:  " << thetaBins << endl;
  
  // We output results in an array of histograms to the parameter list.
  // Create a template for such a histogram, needed by MH3
  // -------------------------------------
  MBinning *binsmh3x = new MBinning("BinningMH3X");
  // dummy binning, real one follows some lines down
  binsmh3x->SetEdgesLog(energyBins,aetHisto->GetYaxis()->GetBinLowEdge(1),
			aetHisto->GetYaxis()->GetBinLowEdge(energyBins+1));
  parlist->AddToList(binsmh3x);  

  MBinning *binsmh3y = new MBinning("BinningMH3Y");
  // dummy binning, real one follows some lines down
  binsmh3y->SetEdges(thetaBins,aetHisto->GetZaxis()->GetBinLowEdge(1),
		     aetHisto->GetZaxis()->GetBinLowEdge(thetaBins+1));
  parlist->AddToList(binsmh3y);  
  
  MH3 *signalHisto = new MH3(2); // Signal(E,t)
  signalHisto->SetupFill(parlist);
  parlist->AddToList(signalHisto);
  signalHisto->SetName("MHOnSubtractionSignal");
  signalHisto->SetTitle("Gamma Events");

  MH3 *offHisto = new MH3(2); // Off(E,t)
  offHisto->SetupFill(parlist);
  parlist->AddToList(offHisto);
  offHisto->SetName("MHOnSubtractionOff");
  offHisto->SetTitle("Background Events");

  MH3 *signifHisto = new MH3(2); // Significance(E,t)
  signifHisto->SetupFill(parlist);
  parlist->AddToList(signifHisto);
  signifHisto->SetName("MHOnSubtractionSignif");
  signifHisto->SetTitle("Significance");

//   if (!fChiSquareHisto) 
//     fChiSquareHisto = new TH1D("fChiSquareHisto", "#chi^{2}/d.o.f. of fits", 50, 0, 5);
//   if (!fSignificanceHisto) 
//     fSignificanceHisto = new TH1D("fSignificanceHisto", "Significances", 40.5, -0.5, 41);
//   if (!fSummedAlphaPlots)
//     fSummedAlphaPlots = new TH1D("fSummedAlphaPlots", "Cumulative Alpha", 
// 				 alphaBins,aetHisto->GetXaxis()->GetBinLowEdge(1),
// 				 aetHisto->GetXaxis()->GetBinLowEdge(alphaBins+1));
  
  TCanvas *c4 = new TCanvas("c4", "MHOnSubtraction Detailed Result Plots", 800,600);
  c4->Divide(energyBins+3,thetaBins,0,0,0);
  
  TH2D& signalTH2DHist = (TH2D&)signalHisto->GetHist();
  TH2D& offTH2DHist =    (TH2D&)offHisto->GetHist();
  TH2D& signifTH2DHist = (TH2D&)signifHisto->GetHist();
  
  signalTH2DHist.Reset();
  signalTH2DHist.SetTitle("Gamma Signal");
  signalTH2DHist.SetXTitle("E [GeV]");
  signalTH2DHist.SetYTitle("theta [deg]");
  signalHisto->SetBinning(&signalTH2DHist, aetHisto->GetYaxis(), aetHisto->GetZaxis());
  signalTH2DHist.Sumw2();
  
  offTH2DHist.Reset();
  offTH2DHist.SetTitle("Off Signal");
  offTH2DHist.SetXTitle("E [GeV]");
  offTH2DHist.SetYTitle("theta [deg]");
  offHisto->SetBinning(&offTH2DHist, aetHisto->GetYaxis(), aetHisto->GetZaxis());
  offTH2DHist.Sumw2();  

  signifTH2DHist.Reset();
  signifTH2DHist.SetTitle("Significance");
  signifTH2DHist.SetXTitle("E [GeV]");
  signifTH2DHist.SetYTitle("theta [deg]");
  signifHisto->SetBinning(&signifTH2DHist, aetHisto->GetYaxis(), aetHisto->GetZaxis());
  signifTH2DHist.Sumw2();
 
  for (Int_t thetaBin = 1; thetaBin < thetaBins+1; thetaBin++) {  

    *fLog << dbg << "--------------------------------------" << endl;
    *fLog << dbg << "Processing ALPHA plots for theta bin " << thetaBin << endl;           
    *fLog << dbg << "--------------------------------------" << endl;
    aetHisto->GetZaxis()->SetRange(thetaBin, thetaBin);
    TH2F* aeHisto = (TH2F*)aetHisto->Project3D("yxe");
    aeHisto->SetDirectory(NULL);
    char hname[80];
    sprintf(hname, "%2.1f < #theta < %2.1f", 
	    aetHisto->GetZaxis()->GetBinLowEdge(thetaBin),
	    aetHisto->GetZaxis()->GetBinLowEdge(thetaBin+1));
    *fLog << inf << "There are " << aeHisto->GetEntries() 
 	  << " entries in " << hname << endl;
    aeHisto->SetTitle(hname);
    sprintf(hname, "MHOnSubtractionThetaBin%d", thetaBin);
    aeHisto->SetName(hname);

    c4->cd((energyBins+3)*(thetaBin-1)+1);
    aeHisto->SetMarkerColor(3);
    aeHisto->DrawCopy();

    c4->Update();

    // We hand over a nonstandard parameter list, which
    // contails lower-dimensional result histograms 
    // signalHisto, offHisto and signifHisto

//     cout << fLog->GetDebugLevel() << endl;
    fLog->SetDebugLevel(1);

    MParList *parlistth2 = new MParList();
    //    parlistth2->SetNullOutput();
    parlistth2->AddToList(binsmh3x);  

    MH3* signalHisto2 = new MH3(1); // Signal(E)
    signalHisto2->SetupFill(parlistth2);
    parlistth2->AddToList(signalHisto2);
    signalHisto2->SetName("MHOnSubtractionSignal");
    signalHisto2->SetTitle("Gamma Events");

    MH3* offHisto2 = new MH3(1); // Off(E)
    offHisto2->SetupFill(parlistth2);
    parlistth2->AddToList(offHisto2);
    offHisto2->SetName("MHOnSubtractionOff");
    offHisto2->SetTitle("Background Events");

    MH3* signifHisto2 = new MH3(1); // Significance(E)
    signifHisto2->SetupFill(parlistth2);
    parlistth2->AddToList(signifHisto2);
    signifHisto2->SetName("MHOnSubtractionSignif");
    signifHisto2->SetTitle("Significance");

    fLog->SetDebugLevel(-1);

    TH1D& signalTH1DHist = (TH1D&)signalHisto2->GetHist();
    TH1D& offTH1DHist =    (TH1D&)offHisto2->GetHist();
    TH1D& signifTH1DHist = (TH1D&)signifHisto2->GetHist();

    signalTH1DHist.Sumw2();
    offTH1DHist.Sumw2();
    signifTH1DHist.Sumw2();

    TH2Calc(aeHisto, parlistth2, kFALSE, c4,   
	    (energyBins+3)*(thetaBin-1)+4);


    for (Int_t energyBin = 1; energyBin < energyBins+1; energyBin++) {  
    
   //    cout << "filling " << thetaBin << " " << energyBin << ": " 
// 		      << signalTH1DHist.GetBinContent(energyBin) << "+-" 
// 	   << signalTH1DHist.GetBinError(energyBin)  <<  " " 
// 		      << offTH1DHist.GetBinContent(energyBin) << "+-" 
// 		      << offTH1DHist.GetBinError(energyBin) << endl;
        
      if (signalTH1DHist.GetBinContent(energyBin)>=0.0) {
	signalTH2DHist.SetBinContent(energyBin, thetaBin,
				     signalTH1DHist.GetBinContent(energyBin));
	signalTH2DHist.SetBinError(energyBin, thetaBin,
				   signalTH1DHist.GetBinError(energyBin));
      }
      
      if (offTH1DHist.GetBinContent(energyBin)>=0.0) {
	offTH2DHist.SetBinContent(energyBin, thetaBin, 
				  offTH1DHist.GetBinContent(energyBin));
	offTH2DHist.SetBinError(energyBin, thetaBin,
				offTH1DHist.GetBinError(energyBin));
      }
      
      signifTH2DHist.SetBinContent(energyBin, thetaBin,
 				   signifTH1DHist.GetBinContent(energyBin));
      signifTH2DHist.SetBinError(energyBin, thetaBin,
 				 signifTH1DHist.GetBinError(energyBin));           
    } //energyBin

  
    c4->cd((energyBins+3)*(thetaBin-1)+2);
    sprintf(hname,"c4_%d",(energyBins+3)*(thetaBin-1)+2);     
    TPad* tp = (TPad*)(gROOT->FindObject(hname));
    tp->SetLogx();     
    signalTH1DHist.SetLineColor(2);
    signalTH1DHist.DrawCopy();     
    offTH1DHist.SetLineColor(4);
    offTH1DHist.DrawCopy("SAME");          
    c4->Update();
           
    signalTH2DHist.SetEntries(signalTH2DHist.GetEntries()+signalTH1DHist.GetEntries());
    offTH2DHist.SetEntries(offTH2DHist.GetEntries()+offTH1DHist.GetEntries());
    signifTH2DHist.SetEntries(signifTH2DHist.GetEntries()+signifTH1DHist.GetEntries());
    
    delete signifHisto2;
    delete offHisto2;
    delete signalHisto2;
    delete parlistth2;
  
  }


  c4->Update();

  return kTRUE;
}



















































/*

// -----------------------------------------------------------------------
//
// Extraction of Gamma signal is performed on a TH3D histogram in
// ALPHA, Energy and Theta. Output to parameter list: TH2 histograms
// in energy and Theta.
//
Bool_t MHOnSubtraction::CalcLightCurve(TH3 *aetHisto, MParList *parlist, const Bool_t draw)
{
  // Analyze aetHisto
  // -------------------------------------
  Int_t alphaBins = aetHisto->GetNbinsX(); 
  Int_t energyBins = aetHisto->GetNbinsY();  
  Int_t timeBins =  aetHisto->GetNbinsZ();  

  *fLog << "Total events:      " << aetHisto->GetEntries() << endl;
  *fLog << "Total energy bins: " << energyBins << endl;
  *fLog << "Total time bins:   " << timeBins << endl;

  // We output results in an array of histograms to the parameter list.
  // Create a template for such a histogram, needed by MH3
  // -------------------------------------
  MBinning *binsmh3x = new MBinning("BinningMH3X");
  // dummy binning, real one follows some lines down
  binsmh3x->SetEdges(timeBins,aetHisto->GetZaxis()->GetBinLowEdge(1),
		     aetHisto->GetZaxis()->GetBinLowEdge(timeBins+1));
  parlist->AddToList(binsmh3x); 
  
  MH3 *signalHisto = new MH3(1); // Signal(t)
  signalHisto->SetupFill(parlist);
  parlist->AddToList(signalHisto);
  signalHisto->SetName("MHOnSubtractionSignal");
  signalHisto->SetTitle("Gamma Events");

  MH3 *offHisto = new MH3(1); // Off(t)
  offHisto->SetupFill(parlist);
  parlist->AddToList(offHisto);
  offHisto->SetName("MHOnSubtractionOff");
  offHisto->SetTitle("Background Events");

  MH3 *signifHisto = new MH3(1); // Significance(t)
  signifHisto->SetupFill(parlist);
  parlist->AddToList(signifHisto);
  signifHisto->SetName("MHOnSubtractionSignif");
  signifHisto->SetTitle("Significance");

  TH1D& signalTH1DHist = (TH1D&)signalHisto->GetHist();
  TH1D& offTH1DHist =    (TH1D&)offHisto->GetHist();
  TH1D& signifTH1DHist = (TH1D&)signifHisto->GetHist();
  
  signalTH1DHist.Reset();
  signalTH1DHist.SetTitle("Gamma Signal");
  signalTH1DHist.SetXTitle("Time [MJD]");
  signalTH1DHist.SetYTitle("Events");
  signalHisto->SetBinning(&signalTH1DHist, aetHisto->GetZaxis());
  signalTH1DHist.Sumw2();
  
  offTH1DHist.Reset();
  offTH1DHist.SetTitle("Background Signal");
  offTH1DHist.SetXTitle("Time [MJD]");
  offTH1DHist.SetYTitle("Events");
  offHisto->SetBinning(&offTH1DHist,  aetHisto->GetZaxis());
  offTH1DHist.Sumw2();  

  signifTH1DHist.Reset();
  signifTH1DHist.SetTitle("Significance");
  signifTH1DHist.SetXTitle("Time [MJD]");
  signifTH1DHist.SetYTitle("Significance #sigma");
  signifHisto->SetBinning(&signifTH1DHist, aetHisto->GetZaxis());
  signifTH1DHist.Sumw2();
 
  //The following histogram is an additional histogram needed for
  //the light curve

  TCanvas *c4 = new TCanvas("c4", "MHOnSubtraction Detailed Result Plots", 800,600);
  c4->Divide(8,7,0,0,0);

 

  // The following loop fixes a common integration region for each
  // energy bin and alerts when no significant excess is found.
  
  Double_t minLowerBin = -10; // minimum integration range
  Double_t maxUpperBin =  10; // minimum integration range
  Double_t maxAlpha =  70; // maximum integration range
  Double_t sumSigLiMa = 0;

  TH1F* alphaHisto[timeBins];

  for (Int_t timeBin = 1; timeBin < timeBins+1; timeBin++) {  

    *fLog << dbg << "--------------------------------------" << endl;
    *fLog << dbg << "Processing ALPHA plots for time bin " << timeBin << endl;           
    *fLog << dbg << "--------------------------------------" << endl;
   
    aetHisto->GetZaxis()->SetRange(timeBin, timeBin);

    char hname[80];
    sprintf(hname, "MHOnSubtractionTimeBin%d", timeBin);    
    alphaHisto[timeBin-1] =
      new TH1F(hname, "ALPHA histogram",
	       alphaBins,aetHisto->GetXaxis()->GetBinLowEdge(1),
	       aetHisto->GetXaxis()->GetBinLowEdge(alphaBins+1));
    alphaHisto[timeBin-1]->SetDirectory(NULL);
    alphaHisto[timeBin-1]->Sumw2();
    alphaHisto[timeBin-1] = (TH1F*)aetHisto->Project3D("xe");
    alphaHisto[timeBin-1]->SetName(hname);
    alphaHisto[timeBin-1]->SetDirectory(NULL);

    sprintf(hname, "%6.0f < t < %6.0f", 
	    aetHisto->GetZaxis()->GetBinLowEdge(timeBin),
	    aetHisto->GetZaxis()->GetBinLowEdge(timeBin+1));
    *fLog << inf << "There are " << alphaHisto[timeBin-1]->GetEntries() 
 	  << " entries in " << hname << endl;
    alphaHisto[timeBin-1]->SetTitle(hname);
      
    // Find signal region and significance in histogram
    Double_t lowerBin, upperBin, sSigLiMa;
    char funcName[20];
    sprintf(funcName, "%2d", timeBin);

    Bool_t drawFit = kTRUE;

    c4->cd(timeBin);
//     alphaHisto[timeBin-1]->SetMarkerColor(3);
    alphaHisto[timeBin-1]->DrawCopy();

    c4->Update();

    fSigniPlotColor = 0;
 ;
    Bool_t fit = FitHistogram(*alphaHisto[timeBin-1], sSigLiMa, 
			      lowerBin, upperBin, (Float_t)3.0, drawFit, 
			      funcName);

    if (fit) { 
      if (sSigLiMa < 3)
	*fLog << warn << "No significant excess for this bin (sigma="<< sSigLiMa <<")!"<<endl;
      sumSigLiMa+=sSigLiMa;
   
      // Evaluate integration range
      lowerBin = (lowerBin < -maxAlpha) ? -maxAlpha : lowerBin;
      minLowerBin = (lowerBin < minLowerBin) ? lowerBin : minLowerBin;    
      upperBin = (upperBin > maxAlpha) ? maxAlpha : upperBin;
      maxUpperBin = (upperBin > maxUpperBin) ? upperBin : maxUpperBin;
    }


    



  } //timeBin
  
  *fLog << inf << "=> Integration range will be " << minLowerBin << " < ALPHA < "
	<< maxUpperBin << "," << endl;    
  *fLog << inf << "   Summed estimated significance is " << sumSigLiMa << endl;
      
  // we have an idea of what is going on in the ALPHA plot...
  // So now we can indeed extract the signal.
         
  sumSigLiMa = 0;
  for (Int_t timeBin = 1; timeBin < timeBins+1; timeBin++) {
    *fLog << inf << "Processing ALPHA distribution for time bin " << timeBin << endl;
    if (alphaHisto[timeBin-1]->GetEntries() == 0) {
       *fLog << warn << "No attempt to extract a signal from 0 events." << endl;       
       if (draw) {
	 TPaveLabel* lab = new TPaveLabel(0.2,0.4,0.8,0.6,"-(nothing to extract)-");
	 lab->SetBorderSize(0);
	 lab->Draw();  
       }
    } else {
      char funcName[20];
      sprintf(funcName, "%2d", timeBin);
      
      Double_t gammaSignal, errorGammaSignal, off, errorOff, sigLiMa;    

      Bool_t drawFit = kFALSE;     

      ExtractSignal(*alphaHisto[timeBin-1],
		    sigLiMa, minLowerBin, maxUpperBin, 
		    gammaSignal, errorGammaSignal, off, errorOff, (Float_t)3, drawFit,
		    funcName, NULL, 1);

      sumSigLiMa += sigLiMa;      

      fMaxSignif = sigLiMa > fMaxSignif ? sigLiMa : fMaxSignif;

      // Fill Signal 
      TH1D &h = (TH1D&)(signalHisto->GetHist());
      h.SetBinContent(timeBin, gammaSignal);
      h.SetBinError(timeBin, errorGammaSignal);
      h.SetEntries(h.GetEntries()+gammaSignal);
      
      // Fill Off Events
      TH1D &ho = (TH1D&)(offHisto->GetHist());
      ho.SetBinContent(timeBin, off);
      ho.SetBinError(timeBin, errorOff);
      ho.SetEntries(ho.GetEntries()+off);
      
      // Fill Significance
      TH1D &hg = (TH1D&)(signifHisto->GetHist());
      hg.SetBinContent(timeBin, sigLiMa);
      hg.SetEntries(hg.GetEntries()+sigLiMa);           
    }
    
  }//timeBin

  *fLog << inf << "Summed significance is " << sumSigLiMa << endl;


  if (draw) {
    Double_t sigLiMa, lowerBin, upperBin;    
    FitHistogram(*fSummedAlphaPlots, sigLiMa, lowerBin, upperBin);
    fSummedAlphaPlots->SetTitle("Cumulative ALPHA distribution");

    *fLog << inf << "Setting range for Significance histo " << fMaxSignif << endl;
    // fSignificanceHisto->GetXaxis()->SetRange(0,fMaxSignif+1);
    
    // *fLog << inf << "Setting range for chisq histo " << fMaxRedChiSq << endl;
    // fChiSquareHisto->GetXaxis()->SetRange(0,fMaxRedChiSq+5);
    
    fChiSquareHisto->DrawClone();
    fSignificanceHisto->DrawClone();
    fSummedAlphaPlots->DrawClone();
  }
  
  return kTRUE;
}





*/












































































// -----------------------------------------------------------------------
//
// Extraction of Gamma signal is performed on a TH2 histogram in
// ALPHA and Energy. Output to parameter list is TH1 histogram in
// energy
//
Bool_t MHOnSubtraction::TH2Calc(TH2 *aeHisto, MParList *parlist, const Bool_t draw, TPad *drawPad, Int_t drawBase)
{

  fSigniPlotColor = 0;

  // Analyze aeHisto
  // -------------------------------------
  const Int_t alphaBins =  aeHisto->GetNbinsX();  // # of alpha bins 
  const Int_t energyBins = aeHisto->GetNbinsY();  

  // Check availability of result histograms
  // -------------------------------------

  MH3* signalHisto = (MH3*)parlist->FindObject("MHOnSubtractionSignal","MH3");
  MH3* offHisto =    (MH3*)parlist->FindObject("MHOnSubtractionOff","MH3");
  MH3* signifHisto = (MH3*)parlist->FindObject("MHOnSubtractionSignif","MH3");

  if (signalHisto && offHisto && signifHisto) {  
//     *fLog << dbg << "Result histograms are present in parameter list " <<
//       "and will be used further." << endl;
  } else {

    *fLog << warn << "Result histograms are not present in parameter " <<
      "list and thus are going to be created now." << endl;

    MBinning *binsmh3x = new MBinning("BinningMH3X");
    binsmh3x->SetEdgesLog(energyBins,aeHisto->GetYaxis()->GetBinLowEdge(1),
		       aeHisto->GetYaxis()->GetBinLowEdge(energyBins+1));
    parlist->AddToList(binsmh3x);  
  
    signalHisto = new MH3(1); // Signal(E)
    signalHisto->SetName("MHOnSubtractionSignal");
    signalHisto->SetTitle("Extracted Signal");
    parlist->AddToList(signalHisto);
    signalHisto->SetBinning(&((TH1D&)signalHisto->GetHist()), 
			    aeHisto->GetYaxis());

    offHisto = new MH3(1); // Off(E)
    offHisto->SetName("MHOnSubtractionOff");
    offHisto->SetTitle("Off Signal");
    parlist->AddToList(offHisto);
    offHisto->SetBinning(&((TH1D&)offHisto->GetHist()), 
			    aeHisto->GetYaxis());
    
    signifHisto = new MH3(1); // Significance(E)
    signifHisto->SetName("MHOnSubtractionSignif");
    signifHisto->SetTitle("Significance");
    parlist->AddToList(signifHisto);
    signifHisto->SetBinning(&((TH1D&)signifHisto->GetHist()), 
			    aeHisto->GetYaxis());

  }
  if (fChiSquareHisto==0x0)
    fChiSquareHisto = new TH1D("fChiSquareHisto", "#chi^{2}/d.o.f. of fits", 50, 0, 5);
  if (fSignificanceHisto==0x0)
    fSignificanceHisto = new TH1D("fSignificanceHisto", "Significances", 41, -0.5, 40.5);
  if (fSummedAlphaPlots==0x0)
    fSummedAlphaPlots = new TH1D("fSummedAlphaPlots", "Cumulative Alpha", 
				 alphaBins,aeHisto->GetXaxis()->GetBinLowEdge(1),
				 aeHisto->GetXaxis()->GetBinLowEdge(alphaBins+1));  
         
  // The following loop fixes a common integration region for each
  // energy bin and alerts when no significant excess is found.
  
  Double_t minLowerBin = -10; // minimum integration range
  Double_t maxUpperBin =  10; // minimum integration range
  Double_t maxAlpha =  70; // maximum integration range
  Double_t sumSigLiMa = 0;

  TH1F* alphaHisto[energyBins];

  fSigniPlotIndex = 0; // cf. ExtractSignal
 
  for (Int_t energyBin = 1; energyBin < energyBins+1; energyBin++) {
    *fLog << inf << "Preprocessing ALPHA distribution for E bin " << energyBin << endl;
    char hname[80];
    sprintf(hname, "MHOnSubtractionAlpha%d", energyBin);    
    alphaHisto[energyBin-1] =
      new TH1F(hname, "ALPHA histogram",
	       alphaBins,aeHisto->GetXaxis()->GetBinLowEdge(1),
	       aeHisto->GetXaxis()->GetBinLowEdge(alphaBins+1));
    alphaHisto[energyBin-1]->SetDirectory(NULL);
    alphaHisto[energyBin-1]->Sumw2();
    alphaHisto[energyBin-1] = (TH1F*)aeHisto->ProjectionX("xe", energyBin, energyBin);
    alphaHisto[energyBin-1]->SetName(hname);
    alphaHisto[energyBin-1]->SetDirectory(NULL);

    sprintf(hname, "%2.1f < E < %2.1f", 
	    aeHisto->GetYaxis()->GetBinLowEdge(energyBin),
	    aeHisto->GetYaxis()->GetBinLowEdge(energyBin+1));
    alphaHisto[energyBin-1]->SetTitle(hname);
    //    alphaHisto[energyBin-1]->DrawCopy();
    alphaHisto[energyBin-1]->SetDirectory(NULL);
      
    // Find signal region and significance in histogram
    Double_t lowerBin, upperBin, sSigLiMa;
    char funcName[20];
    sprintf(funcName, "%2d", energyBin);

    Bool_t drawFit = kFALSE;

    if (drawPad) {
      drawPad->cd(drawBase+energyBin-1);
      drawFit = kTRUE;
    }

    Bool_t fit = FitHistogram(*alphaHisto[energyBin-1], sSigLiMa, 
			      lowerBin, upperBin, (Float_t)3.0, drawFit, 
			      funcName);

    if (fit) { 
      Double_t redChiSq = alphaHisto[energyBin-1]->GetFunction("gauspol3"+(TString)funcName)->GetChisquare()/
	alphaHisto[energyBin-1]->GetFunction("gauspol3"+(TString)funcName)->GetNDF();
      fChiSquareHisto->Fill(redChiSq);
      fMaxRedChiSq = redChiSq > fMaxRedChiSq ? redChiSq : fMaxRedChiSq;
      fSummedAlphaPlots->Add(alphaHisto[energyBin-1], 1);
      if (sSigLiMa < 3)
	*fLog << warn << "No significant excess for this bin (sigma="<< sSigLiMa <<")!"<<endl;
      sumSigLiMa+=sSigLiMa;
   
      // Evaluate integration range
      lowerBin = (lowerBin < -maxAlpha) ? -maxAlpha : lowerBin;
      minLowerBin = (lowerBin < minLowerBin) ? lowerBin : minLowerBin;    
      upperBin = (upperBin > maxAlpha) ? maxAlpha : upperBin;
      maxUpperBin = (upperBin > maxUpperBin) ? upperBin : maxUpperBin;
    } else {
      fChiSquareHisto->Fill(0);
    }

//     debugOut->Update();

  } //energyBin
  
  *fLog << inf << "=> Integration range for this theta bin will be " << minLowerBin << " < ALPHA < "
	<< maxUpperBin << "," << endl;    
  *fLog << inf << "   Summed estimated significance is " << sumSigLiMa << endl;
      
  // we have an idea of what is going on in the ALPHA plot...
  // So now we can indeed extract the signal.
         
  sumSigLiMa = 0;
  for (Int_t energyBin = 1; energyBin < energyBins+1; energyBin++) {
    *fLog << inf << "Processing ALPHA distribution for E bin " << energyBin << endl;
    if (alphaHisto[energyBin-1]->GetEntries() == 0) {
       *fLog << warn << "No attempt to extract a signal from 0 events." << endl;       
       if (draw) {
	 TPaveLabel* lab = new TPaveLabel(0.2,0.4,0.8,0.6,"-(nothing to extract)-");
	 lab->SetBorderSize(0);
	 lab->Draw();  
       }
    } else {
      char funcName[20];
      sprintf(funcName, "%2d", energyBin);
      
      Double_t gammaSignal, errorGammaSignal, off, errorOff, sigLiMa;    

      Bool_t drawFit = kFALSE;
      
//       if (drawPad) {
// 	cout << "Going to use pad " <<drawBase+energyBin-1 << endl;
// 	drawPad->cd(drawBase+energyBin-1);
// 	drawFit = kTRUE;
//       }


      ExtractSignal(*alphaHisto[energyBin-1],
		    sigLiMa, minLowerBin, maxUpperBin, 
		    gammaSignal, errorGammaSignal, off, errorOff, (Float_t)3, drawFit,
		    funcName, drawPad, drawBase);

      sumSigLiMa += sigLiMa;      
      fSignificanceHisto->Fill(sigLiMa);
      fMaxSignif = sigLiMa > fMaxSignif ? sigLiMa : fMaxSignif;

      // Fill Signal 
      TH1D &h = (TH1D&)(signalHisto->GetHist());
      h.SetBinContent(energyBin, gammaSignal);
      h.SetBinError(energyBin, errorGammaSignal);
      h.SetEntries(h.GetEntries()+gammaSignal);
      
      // Fill Off Events
      TH1D &ho = (TH1D&)(offHisto->GetHist());
      ho.SetBinContent(energyBin, off);
      ho.SetBinError(energyBin, errorOff);
      ho.SetEntries(ho.GetEntries()+off);
      
      // Fill Significance
      TH1D &hg = (TH1D&)(signifHisto->GetHist());
      hg.SetBinContent(energyBin, sigLiMa);
      hg.SetEntries(hg.GetEntries()+sigLiMa);           
    }
    
  }//energyBin

  *fLog << inf << "Summed significance is " << sumSigLiMa << endl;


  if (draw) {
    Double_t sigLiMa, lowerBin, upperBin;    
    FitHistogram(*fSummedAlphaPlots, sigLiMa, lowerBin, upperBin);
    fSummedAlphaPlots->SetTitle("Cumulative ALPHA distribution");

    *fLog << inf << "Setting range for Significance histo " << fMaxSignif << endl;
    // fSignificanceHisto->GetXaxis()->SetRange(0,fMaxSignif+1);
    
    // *fLog << inf << "Setting range for chisq histo " << fMaxRedChiSq << endl;
    // fChiSquareHisto->GetXaxis()->SetRange(0,fMaxRedChiSq+5);
    
    fChiSquareHisto->DrawClone();
    fSignificanceHisto->DrawClone();
    fSummedAlphaPlots->DrawClone();
  }
  
  return kTRUE;
}
















// -----------------------------------------------------------------------
//
// Extraction of Gamma signal is performed on a TH1 histogram in ALPHA
// 
//
Bool_t MHOnSubtraction::Calc(TH1 *aHisto, MParList *parlist, 
			     const Bool_t draw, Float_t signalRegion)
{ 
  // Find signal region and estimate significance
  Double_t lowerBin, upperBin, sigLiMa;
  FitHistogram(*aHisto,	sigLiMa, lowerBin, upperBin, (Float_t)3.5, kFALSE);

  //if (sigLiMa < 3)
  //  *fLog << err << "No significant excess (sigma=" << sigLiMa <<")!"<<endl;

  if (signalRegion!=0) {
    lowerBin = -signalRegion;
    upperBin =  signalRegion;
  }
       
  Double_t gammaSignal, errorGammaSignal, off, errorOff;

  ExtractSignal(*aHisto, sigLiMa, lowerBin, upperBin, 
		gammaSignal, errorGammaSignal, off, errorOff, 
		(Float_t)3.0, draw);

  fSignificance = sigLiMa;
  
  *fLog << inf << "Signal is " 
	<< gammaSignal << "+-" << errorGammaSignal << endl
	<< inf << "Significance is " << sigLiMa << endl;

  return kTRUE;
}

// -------------------------------------------------------------------------
//
// Draw a copy of the histogram
//
TObject *MHOnSubtraction::DrawClone(Option_t *opt) const
{  
    TCanvas &c = *MakeDefCanvas("OnSubtraction", "Results from OnSubtraction");
    c.Divide(2,2);

    gROOT->SetSelectedPad(NULL);
    gStyle->SetOptStat(0);

    c.cd(1);
    fSummedAlphaPlots->DrawCopy(opt);

    c.cd(2);
    fSignificanceHisto->DrawCopy(opt);

    c.cd(3);  
    TString drawopt="";
    for (fThetaHistoArray->SetIndex(0); 
 	 MH3* h=(MH3*)(fThetaHistoArray->GetH()); 
 	 fThetaHistoArray->IncIndex())
      {
	// Get the current histogram
	TH1D& hist = (TH1D&)h->GetHist();
	fSummedAlphaPlots->SetTitle("Energy distributions for different theta bins");
	hist.Draw(drawopt);
	drawopt="SAME";
      }
    fThetaLegend->Draw();
    
    c.cd(4);
    fChiSquareHisto->DrawCopy(opt);
    
    c.Modified();
    c.Update();

    return &c;
}

// -------------------------------------------------------------------------
//
// Draw the histogram
//
void MHOnSubtraction::Draw(Option_t *opt) 
{ 
    if (!gPad)
      MakeDefCanvas("OnSubtraction", "Results from OnSubtraction");

    gPad->Divide(2,3);

// Ok, at some point this will be implemented

    gPad->Modified();
    gPad->Update();
}





