/* ======================================================================== *\
!
! *
! * 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 express
! * or implied warranty.
! *
!
!
!   Author(s): Wolfgang Wittek, 06/2003 <mailto:wittek@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//  MPad
//
//  The aim of the padding is to make the noise level (NSB + electronic noise)
//  equal for the MC, ON and OFF data samples
//  
//  The target noise level is determined by 'merging' the (sigmabar vs. Theta)
//  distributions of MC, ON and OFF data.
//
//  The prescription on which fraction of events has to padded from
//  bin k of sigmabar to bin j is stored in the matrices 
//  fHgMC, fHgON and fHgOFF.
//
//  By the padding the number of photons, its error and the pedestal sigmas 
//  are altered. On average, the number of photons added is zero.
//
//  The formulas used can be found in Thomas Schweizer's Thesis,
//                                    Section 2.2.1
//
//  The padding is done as follows :
//
//     Increase the sigmabar according to the matrices fHg.... Then generate 
//     a pedestal sigma for each pixel using the 3D-histogram Theta, pixel no.,
//     Sigma^2-Sigmabar^2 (fHDiffPixTheta).
//
//  The padding has to be done before the image cleaning because the
//  image cleaning depends on the pedestal sigmas.
//
/////////////////////////////////////////////////////////////////////////////
#include "MPad.h"

#include <math.h>
#include <stdio.h>

#include <TH1.h>
#include <TH2.h>
#include <TH3.h>
#include <TRandom.h>
#include <TCanvas.h>
#include <TFile.h>
#include <TVirtualPad.h>

#include "MBinning.h"
#include "MSigmabar.h"
#include "MPointingPos.h"
#include "MLog.h"
#include "MLogManip.h"
#include "MParList.h"
#include "MGeomCam.h"

#include "MCerPhotPix.h"
#include "MCerPhotEvt.h"

#include "MH.h"

#include "MPedPhotCam.h"
#include "MPedPhotPix.h"

#include "MBlindPixels.h"

#include "MRead.h"
#include "MFilterList.h"
#include "MTaskList.h"
#include "MBlindPixelCalc.h"
#include "MHBlindPixels.h"
#include "MFillH.h"
#include "MHSigmaTheta.h"
#include "MEvtLoop.h"

ClassImp(MPad);

using namespace std;

// --------------------------------------------------------------------------
//
// Default constructor. 
//
MPad::MPad(const char *name, const char *title) 
{
  fName  = name  ? name  : "MPad";
  fTitle = title ? title : "Task for the padding";

  fType = "";

  fHSigmaThetaMC     = NULL;
  fHSigmaThetaON     = NULL;
  fHSigmaThetaOFF    = NULL;

  fHDiffPixThetaMC   = NULL;
  fHDiffPixThetaON   = NULL;
  fHDiffPixThetaOFF  = NULL;

  fHSigmaPixThetaMC   = NULL;
  fHSigmaPixThetaON   = NULL;
  fHSigmaPixThetaOFF  = NULL;

  fHgMC  = NULL;
  fHgON  = NULL;
  fHgOFF = NULL;

  fHBlindPixIdThetaMC  = NULL;
  fHBlindPixIdThetaON  = NULL;
  fHBlindPixIdThetaOFF = NULL;

  fHBlindPixNThetaMC   = NULL;
  fHBlindPixNThetaON   = NULL;
  fHBlindPixNThetaOFF  = NULL;

  fHSigmaPedestal = NULL;
  fHPhotons       = NULL;
  fHNSB           = NULL;
}

// --------------------------------------------------------------------------
//
// Destructor. 
//
MPad::~MPad()
{
  delete fHgMC;
  delete fHgON;
  delete fHgOFF;

  delete fHSigmaTheta;
  delete fHSigmaPixTheta;
  delete fHDiffPixTheta;
  delete fHBlindPixNTheta;
  delete fHBlindPixIdTheta;

  delete fHSigmaPedestal;
  delete fHPhotons;
  delete fHNSB;

  delete fInfile;
}
// --------------------------------------------------------------------------
//
// Merge the sigmabar distributions of 2 samples (samples A and B)
//
// input : the original distributions for samples A and B (hista, histb)
//
// output : the prescription which fraction of events should be padded
//          from the sigmabar bin k to bin j (fHgMC, fHgON, fHgOFF)
//

Bool_t MPad::Merge2Distributions(TH1D *hista, TH1D *histb, TH1D *histap,
                                      TH2D *fHga,  TH2D *fHgb, Int_t nbinssig )
{
  // hista is the normalized 1D histogram of sigmabar for sample A
  // histb is the normalized 1D histogram of sigmabar for sample B
  // histc is the difference A-B

  // at the beginning, histap is a copy of hista
  // at the end, it will be the 1D histogram for sample A after padding

  // at the beginning, histbp is a copy of histb
  // at the end, it will be the 1D histogram for sample B after padding

  // at the beginning, histcp is a copy of histc
  // at the end, it should be the difference histap-histbp,
  //             which should be zero

  // fHga[k][j]  tells which fraction of events from sample A should be padded
  //             from sigma_k to sigma_j


  // fHgb[k][j]  tells which fraction of events from sample B should be padded
  //             from sigma_k to sigma_j

  //*fLog << all << "MPad::Merge2Distributions(); Entry" << endl;

  Double_t eps = 1.e-10;

    TH1D *histbp  = new TH1D( (const TH1D&)*histb );

    TH1D *histc   = new TH1D( (const TH1D&)*hista );
    histc->Add(histb, -1.0);

    TH1D *histcp  = new TH1D( (const TH1D&)*histc );    


  //--------   start j loop   ------------------------------------------------
  // loop over bins in histc, 
  // starting from the end (i.e. at the largest sigmabar)
  Double_t v, s, w, t, x, u, a, b, arg;

  //*fLog << "Content of histcp : " << endl;

  for (Int_t j=nbinssig; j >= 1; j--)
  {
    //*fLog << "j, hista, histb , histap : " << j << ",  " 
    //                              <<  hista->GetBinContent(j) << ",  "
    //                              <<  histb->GetBinContent(j) << ",  "
    //                              <<  histap->GetBinContent(j) << endl;

    v = histcp->GetBinContent(j);
    //*fLog << "cp : j, v = " << j << ",  " << v << endl;

    if ( fabs(v) < eps ) continue;
    if (v >= 0.0) 
      s = 1.0;
    else
      s = -1.0;

    //................   start k loop   ......................................
    // look for a bin k which may compensate the content of bin j
    for (Int_t k=j-1; k >= 1; k--)
    {
      w = histcp->GetBinContent(k);
      if ( fabs(w) < eps ) continue;
      if (w >= 0.0) 
        t = 1.0;
      else
        t = -1.0;

      // if s==t bin k cannot compensate, go to next k bin
      if (s == t) continue;

      x = v + w;
      if (x >= 0.0) 
        u = 1.0;
      else
        u = -1.0;

      // if u==s bin k will partly compensate : pad the whole fraction
      // w from bin k to bin j
      if (u == s)
        arg = -w;

      // if u!=s bin k would overcompensate : pad only the fraction
      // v from bin k to bin j
      else
        arg = v;

      if (arg <=0.0)
      {
        fHga->SetBinContent(k, j, -arg);
        fHgb->SetBinContent(k, j,  0.0);
      }
      else
      {
        fHga->SetBinContent(k, j,  0.0);
        fHgb->SetBinContent(k, j,  arg);
      }

      //*fLog << all << "k, j, arg = " << k << ",  " << j 
      //      << ",  " << arg << endl;

      //......................................
      // this is for checking the procedure
      if (arg < 0.0)
      {
        a = histap->GetBinContent(k);
        histap->SetBinContent(k, a+arg);
        a = histap->GetBinContent(j);
        histap->SetBinContent(j, a-arg);
      }
      else
      {
        b = histbp->GetBinContent(k);
        histbp->SetBinContent(k, b-arg);
        b = histbp->GetBinContent(j);
        histbp->SetBinContent(j, b+arg);
      }
      //......................................

      if (u == s)
      {
        histcp->SetBinContent(k, 0.0);
        histcp->SetBinContent(j,   x);

        v = histcp->GetBinContent(j);
        if ( fabs(v) < eps ) break;

        // bin j was only partly compensated : go to next k bin
        continue;
      }
      else
      {
        histcp->SetBinContent(k,   x);
        histcp->SetBinContent(j, 0.0);

        // bin j was completely compensated : go to next j bin
        break;
      }

    }
    //................   end k loop   ......................................
  } 
  //--------   end j loop   ------------------------------------------------

  // check results 

  //*fLog << "Content of histap, histbp, histcp : " << endl;

  for (Int_t j=1; j<=nbinssig; j++)
  {
    Double_t a = histap->GetBinContent(j);
    Double_t b = histbp->GetBinContent(j);
    Double_t c = histcp->GetBinContent(j);

    //*fLog << "j, a, b, c = " << j << ":  " << a << ",  " << b << ",  "
    //      << c << endl;

    if( fabs(a-b)>3.0*eps  ||  fabs(c)>3.0*eps )
      *fLog << err << "MPad::Merge2Distributions(); inconsistency in results; j, a, b, c = "
            << j << ",  " << a << ",  " << b << ",  " << c << endl;
  }

  //---------------------------------------------------------------
  TCanvas *pad = MH::MakeDefCanvas("Merging", "", 600, 600); 
  gROOT->SetSelectedPad(NULL);

  pad->Divide(2, 2);

  pad->cd(1);
  gPad->SetBorderMode(0);
  hista->SetDirectory(NULL);
  hista->Draw();
  hista->SetBit(kCanDelete);    

  pad->cd(2);
  gPad->SetBorderMode(0);
  histb->SetDirectory(NULL);
  histb->Draw();
  histb->SetBit(kCanDelete);    

  pad->cd(3);
  gPad->SetBorderMode(0);
  histap->SetDirectory(NULL);
  histap->Draw();
  histap->SetBit(kCanDelete);    

  pad->DrawClone();

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

  delete histc;
  delete histbp;
  delete histcp;

  //*fLog << all << "MPad::Merge2Distributions(); Exit" << endl;

  return kTRUE;
}


// --------------------------------------------------------------------------
//
// Merge the distributions 
//   fHSigmaTheta      2D-histogram  (Theta, sigmabar)
//
// ===>   of ON, OFF and MC data   <==============
//
// and define the matrices fHgMC, fHgON and fHgOFF 
//
// These matrices tell which fraction of events should be padded 
// from bin k of sigmabar to bin j
//

Bool_t MPad::MergeONOFFMC(
   TH2D& sigthmc,  TH3D& diffpixthmc,  TH3D& sigmapixthmc,
   TH2D& blindidthmc,  TH2D& blindnthmc,
   TH2D& sigthon,  TH3D& diffpixthon,  TH3D& sigmapixthon,
   TH2D& blindidthon,  TH2D& blindnthon,
   TH2D& sigthoff, TH3D& diffpixthoff, TH3D& sigmapixthoff,
   TH2D& blindidthoff, TH2D& blindnthoff)
{
  //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  // for testing
  /*
  TAxis *xa = sigthon.GetXaxis();
  Int_t nbitheta = xa->GetNbins();

  TAxis *ya = sigthon.GetYaxis();
  Int_t nbisig = ya->GetNbins();
  for (Int_t i=1; i<=nbitheta; i++)
  {
    for (Int_t j=1; j<=nbisig; j++)
    {
      if (i>0)
      {
        sigthmc.SetBinContent( i, j, (Float_t) (625000.0+(nbisig*nbisig*nbisig-j*j*j)) );
        sigthon.SetBinContent( i, j, (Float_t)(1.0) );
        sigthoff.SetBinContent(  i, j, (Float_t)( (0.5*nbisig-j)*(0.5*nbisig-j)            ) );
      }
      else
      {
        sigthmc.SetBinContent( i, j, 0.0);
        sigthon.SetBinContent( i, j, 0.0);
        sigthoff.SetBinContent(i, j, 0.0);
      }
    }
  }
  */
  //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

  *fLog << all << "----------------------------------------------------------------------------------" << endl;
  *fLog << all << "MPad::MergeONOFFMC(); Merge the ON, OFF and MC histograms to obtain the target distributions for the padding"
        << endl;

  fHSigmaTheta = new TH2D( (const TH2D&) sigthon );
  fHSigmaTheta->SetNameTitle("2D-ThetaSigmabar", "2D-ThetaSigmabar (target)");

  *fLog << "after booking fHSigmaTheta" << endl;

  fHDiffPixTheta = new TH3D( (const TH3D&) diffpixthon );
  fHDiffPixTheta->SetNameTitle("3D-ThetaPixDiff", "3D-ThetaPixDiff (target)");

  *fLog << "after booking fHDiffPixTheta" << endl;

  fHSigmaPixTheta = new TH3D( (const TH3D&) sigmapixthon );
  fHSigmaPixTheta->SetNameTitle("3D-ThetaPixSigma", "3D-ThetaPixSigma (target)");
  *fLog << "after booking fHSigmaPixTheta" << endl;

  fHBlindPixNTheta = new TH2D( (const TH2D&) blindnthon );
  fHBlindPixNTheta->SetNameTitle("2D-ThetaBlindN", "2D-ThetaBlindN (target)");

  *fLog << "after booking fHBlindPixNTheta" << endl;

  fHBlindPixIdTheta = new TH2D( (const TH2D&) blindidthon );
  fHBlindPixIdTheta->SetNameTitle("2D-ThetaBlindId", "2D-ThetaBlindId (target)");

  *fLog << "after booking fHBlindPixIdTheta" << endl;
  *fLog << all << "Histograms for the merged padding plots were booked" 
        << endl;

  //--------------------------
  fHSigmaThetaMC = &sigthmc;
  fHSigmaThetaMC->SetNameTitle("2D-ThetaSigmabarMC", "2D-ThetaSigmabarMC (orig.)");
  fHSigmaThetaON = &sigthon;
  fHSigmaThetaON->SetNameTitle("2D-ThetaSigmabarON", "2D-ThetaSigmabarON (orig.)");
  fHSigmaThetaOFF = &sigthoff;
  fHSigmaThetaOFF->SetNameTitle("2D-ThetaSigmabarOFF", "2D-ThetaSigmabarOFF (orig.)");

  //*fLog << all << "addresses of SigmaTheta padding plots were copied" 
  //      << endl;

  //--------------------------
  fHBlindPixNThetaMC = &blindnthmc;
  fHBlindPixNThetaMC->SetNameTitle("2D-ThetaBlindNMC", "2D-ThetaBlindNMC (orig.)");
  fHBlindPixNThetaON = &blindnthon;
  fHBlindPixNThetaON->SetNameTitle("2D-ThetaBlindNON", "2D-ThetaBlindNON (orig.)");
  fHBlindPixNThetaOFF = &blindnthoff;
  fHBlindPixNThetaOFF->SetNameTitle("2D-ThetaBlindNOFF", "2D-ThetaBlindNOFF (orig.)");

  //*fLog << all << "addresses of BlindPixTheta padding plots were copied" 
  //      << endl;

  //--------------------------
  fHBlindPixIdThetaMC = &blindidthmc;
  fHBlindPixIdThetaMC->SetNameTitle("2D-ThetaBlindIdMC", "2D-ThetaBlindIdMC (orig.)");
  fHBlindPixIdThetaON = &blindidthon;
  fHBlindPixIdThetaON->SetNameTitle("2D-ThetaBlindIdON", "2D-ThetaBlindIdON (orig.)");
  fHBlindPixIdThetaOFF = &blindidthoff;
  fHBlindPixIdThetaOFF->SetNameTitle("2D-ThetaBlindIdOFF", "2D-ThetaBlindIdOFF (orig.)");

  //*fLog << all << "addresses of BlindPixIdTheta padding plots were copied" 
  //      << endl;

  //--------------------------
  fHDiffPixThetaMC = &diffpixthmc;
  fHDiffPixThetaMC->SetNameTitle("3D-ThetaPixDiffMC", "3D-ThetaPixDiffMC (orig.)");
  fHDiffPixThetaON = &diffpixthon;
  fHDiffPixThetaON->SetNameTitle("3D-ThetaPixDiffON", "3D-ThetaPixDiffON (orig.)");
  fHDiffPixThetaOFF = &diffpixthoff;
  fHDiffPixThetaOFF->SetNameTitle("3D-ThetaPixDiffOFF", "3D-ThetaPixDiffOFF (orig.)");

  //*fLog << all << "addresses of DiffPixTheta padding plots were copied" 
  //      << endl;

  //--------------------------
  fHSigmaPixThetaMC = &sigmapixthmc;
  fHSigmaPixThetaMC->SetNameTitle("3D-ThetaPixSigmaMC", "3D-ThetaPixSigmaMC (orig.)");
  fHSigmaPixThetaON = &sigmapixthon;
  fHSigmaPixThetaON->SetNameTitle("3D-ThetaPixSigmaON", "3D-ThetaPixSigmaON (orig.)");
  fHSigmaPixThetaOFF = &sigmapixthoff;
  fHSigmaPixThetaOFF->SetNameTitle("3D-ThetaPixSigmaOFF", "3D-ThetaPixSigmaOFF (orig.)");

  //*fLog << all << "addresses of SigmaPixTheta padding plots were copied" 
  //      << endl;

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


  // get binning for fHgON, fHgOFF and fHgMC  from sigthon
  // binning in Theta
  TAxis *ax = sigthon.GetXaxis();
  Int_t nbinstheta = ax->GetNbins();
  TArrayD edgesx;
  edgesx.Set(nbinstheta+1);
  for (Int_t i=0; i<=nbinstheta; i++)
  {
    edgesx[i] = ax->GetBinLowEdge(i+1);
    //*fLog << all << "i, theta low edge = " << i << ",  " << edgesx[i] 
    //      << endl; 
  }
  MBinning binth;
  binth.SetEdges(edgesx);

  // binning in sigmabar
  TAxis *ay = sigthon.GetYaxis();
  Int_t nbinssig = ay->GetNbins();
  TArrayD edgesy;
  edgesy.Set(nbinssig+1);
  for (Int_t i=0; i<=nbinssig; i++)
  {
    edgesy[i] = ay->GetBinLowEdge(i+1); 
    //*fLog << all << "i, sigmabar low edge = " << i << ",  " << edgesy[i] 
    //      << endl; 
  }
  MBinning binsg;
  binsg.SetEdges(edgesy);


  fHgMC = new TH3D;
  MH::SetBinning(fHgMC, &binth, &binsg, &binsg);
  fHgMC->SetNameTitle("3D-PaddingMatrixMC", "3D-PadMatrixThetaMC");

  fHgON = new TH3D;
  MH::SetBinning(fHgON, &binth, &binsg, &binsg);
  fHgON->SetNameTitle("3D-PaddingMatrixON", "3D-PadMatrixThetaON");

  fHgOFF = new TH3D;
  MH::SetBinning(fHgOFF, &binth, &binsg, &binsg);
  fHgOFF->SetNameTitle("3D-PaddingMatrixOFF", "3D-PadMatrixThetaOFF");

  //*fLog << all << "fHg histograms were booked" << endl;

  //............   loop over Theta bins   ...........................



  //*fLog << all << "MPad::MergeONOFFMC(); bins of Theta, Sigmabarold, Sigmabarnew, fraction of events to be padded" << endl;
  for (Int_t l=1; l<=nbinstheta; l++)
  {
    //-------------------------------------------
    // merge ON and OFF distributions
    // input : hista is the normalized 1D distr. of sigmabar for ON  data
    //         histb is the normalized 1D distr. of sigmabar for OFF data
    // output : histap    will be the merged distribution (ON-OFF)
    //          fHga(k,j) will tell which fraction of ON events should be 
    //                    padded from bin k to bin j
    //          fHgb(k,j) will tell which fraction of OFF events should be 
    //                    padded from bin k to bin j

    TH2D * fHga = new TH2D;
    MH::SetBinning(fHga, &binsg, &binsg);
    TH2D * fHgb = new TH2D;
    MH::SetBinning(fHgb, &binsg, &binsg);

    TH1D *hista  = sigthon.ProjectionY("sigon_y", l, l, "");
    hista->SetNameTitle("ON-orig", "ON (orig)"); 
    Stat_t suma = hista->Integral();
    if (suma != 0.0) hista->Scale(1./suma);

    TH1D *histap  = new TH1D( (const TH1D&)*hista );
    histap->SetNameTitle("ON-OFF-merged)", "ON-OFF (merged)"); 

    TH1D *histb  = sigthoff.ProjectionY("sigoff_y", l, l, "");
    histb->SetNameTitle("OFF-orig", "OFF (orig)"); 
    Stat_t sumb = histb->Integral();
    if (sumb != 0.0) histb->Scale(1./sumb);

    if (suma == 0.0  ||  sumb == 0.0)
    {
      delete hista;
      delete histb;
      delete fHga;
      delete fHgb;
      delete histap;

      *fLog << err 
            << "cannot call Merge2Distributions(ON, OFF)  for theta bin l = " 
            << l << ";  NevON, NevOFF = " << suma << ",  " << sumb << endl;

      // go to next theta bin (l)
      continue;
    }

    //*fLog << "call Merge2Distributions(ON, OFF)  for theta bin l = " 
    //      << l << endl;

    Merge2Distributions(hista, histb, histap, fHga, fHgb, nbinssig);

    // fill fHgON and fHgOFF
    for (Int_t k=1; k<=nbinssig; k++)
      for (Int_t j=1; j<=nbinssig; j++)
      {
        Double_t a = fHga->GetBinContent(k,j);
        fHgON->SetBinContent(l, k, j, a);

        Double_t b = fHgb->GetBinContent(k,j);
        fHgOFF->SetBinContent(l, k, j, b);
      }

    //-------------------------------------------
    // now merge ON-OFF and MC distributions
    // input : histe is the result of the merging of ON and OFF distributions  
    //         histf is the normalized 1D distr. of sigmabar for MC data
    // output : histep    will be the merged distribution (target distribution)
    //          fHge(k,j) will tell which fraction of ON, OFF events should be 
    //                    padded from bin k to bin j
    //          fHgf(k,j) will tell which fraction of MC events should be 
    //                    padded from bin k to bin j

    TH2D * fHge = new TH2D;
    MH::SetBinning(fHge, &binsg, &binsg);
    TH2D * fHgf = new TH2D;
    MH::SetBinning(fHgf, &binsg, &binsg);

    TH1D *histe  = new TH1D( (const TH1D&)*histap);
    histe->SetNameTitle("ON-OFF-merged", "ON-OFF (merged)"); 

    TH1D *histep  = new TH1D( (const TH1D&)*histe);
    histep->SetNameTitle("ON-OFF-MC-merged)", "ON-OFF-MC (merged)"); 

    TH1D *histf  = sigthmc.ProjectionY("sigmc_y", l, l, "");
    histf->SetNameTitle("MC-orig", "MC (orig)"); 
    Double_t sumf = histf->Integral();
    if (sumf != 0.0) histf->Scale(1./sumf);

    if (sumf == 0.0)
    {
      delete hista;
      delete histb;
      delete fHga;
      delete fHgb;
      delete histap;

      delete histe;
      delete histf;
      delete fHge;
      delete fHgf;
      delete histep;

      *fLog << err 
            << "cannot call Merge2Distributions(ON-OFF,MC)  for theta bin l = " 
            << l << ";  NevMC = " << sumf << endl;

      // go to next theta bin (l)
      continue;
    }

    //*fLog << "call Merge2Distributions(ON-OFF, MC)  for theta bin l = " 
    //      << l << endl;

    Merge2Distributions(histe, histf, histep, fHge, fHgf, nbinssig);

    // update fHgON and fHgOFF
    UpdateHg(fHga, histap, fHge, fHgON,  nbinssig, l);
    UpdateHg(fHgb, histap, fHge, fHgOFF, nbinssig, l);

    // fill fHgMC
    for (Int_t k=1; k<=nbinssig; k++)
      for (Int_t j=1; j<=nbinssig; j++)
      {
        Double_t f = fHgf->GetBinContent(k,j);
        fHgMC->SetBinContent(l, k, j, f);
      }

    // fill target distribution fHSigmaTheta 
    // (only for plotting, it will not  be used in the padding)
    for (Int_t j=1; j<=nbinssig; j++)
    {
      Double_t ep = histep->GetBinContent(j);
      fHSigmaTheta->SetBinContent(l, j, ep);
    }
    //-------------------------------------------

    delete hista;
    delete histb;
    delete fHga;
    delete fHgb;
    delete histap;

    delete histe;
    delete histf;
    delete fHge;
    delete fHgf;
    delete histep;
  }
  //............   end of loop over Theta bins   ....................

  
  // Define the target distributions 
  //        fHDiffPixTheta, fHSigmaPixTheta
  //               (they are calculated as 
  //               averages of the ON and OFF distributions)
  //        fHBlindPixNTheta, fHBlindPixIdTheta
  //               (they are calculated as 
  //               the sum of the ON and OFF distributions)
  // (fHDiffPixTheta will be used in the padding of MC events)

  //............   start of new loop over Theta bins   ....................
  for (Int_t j=1; j<=nbinstheta; j++)
  {
    // fraction of ON/OFF/MC events to be padded : fracON, fracOFF, fracMC

    Double_t fracON  = fHgON->Integral (j, j, 1, nbinssig, 1, nbinssig, "");
    Double_t fracOFF = fHgOFF->Integral(j, j, 1, nbinssig, 1, nbinssig, "");
    Double_t fracMC  = fHgMC->Integral(j, j, 1, nbinssig, 1, nbinssig, "");

    Double_t normON;
    Double_t normOFF;

    // set ranges for 2D-projections of 3D-histograms
    TAxis *ax = diffpixthon.GetXaxis();
    ax->SetRange(j, j);
    ax = diffpixthoff.GetXaxis();
    ax->SetRange(j, j);

    TH2D *hist;
    TH2D *histOFF;

    TH1D *hist1;
    TH1D *hist1OFF;

    // weights for ON and OFF distrubtions when
    // calculating the weighted average
    Double_t facON  = 0.5 * (1. - fracON + fracOFF);
    Double_t facOFF = 0.5 * (1. + fracON - fracOFF);
  
    *fLog << all << "Theta bin j : fracON, fracOFF, fracMC, facON, facOFF = " 
          << j << ",  " << fracON << ",  " << fracOFF << ",  "
          << fracMC << ",  " << facON << ",  " << facOFF << endl; 

    //------------------------------------------------------------------
    // define target distribution  'sigma^2-sigmavar^2 vs. pixel number'
    ay = diffpixthon.GetYaxis();
    Int_t nbinspixel = ay->GetNbins();

    TAxis *az = diffpixthon.GetZaxis();
    Int_t nbinsdiff = az->GetNbins();

    hist    = (TH2D*)diffpixthon.Project3D("zy");
    hist->SetName("dummy");
    histOFF = (TH2D*)diffpixthoff.Project3D("zy");

    normON  = hist->Integral();
    normOFF = histOFF->Integral();
    if (normON  != 0.0) hist->Scale(1.0/normON);
    if (normOFF != 0.0) histOFF->Scale(1.0/normOFF);

    // weighted average of ON and OFF distributions
    hist->Scale(facON);
    hist->Add(histOFF, facOFF); 

    for (Int_t k=1; k<=nbinspixel; k++)
      for (Int_t l=1; l<=nbinsdiff; l++)
      {
        Double_t cont = hist->GetBinContent(k,l);
        fHDiffPixTheta->SetBinContent(j, k, l, cont);  
      }

    delete hist;
    delete histOFF;

    //------------------------------------------------------------------
    // define target distribution  'sigma vs. pixel number'
    ay = sigmapixthon.GetYaxis();
    nbinspixel = ay->GetNbins();

    az = sigmapixthon.GetZaxis();
    Int_t nbinssigma = az->GetNbins();

    hist    = (TH2D*)sigmapixthon.Project3D("zy");
    hist->SetName("dummy");
    histOFF = (TH2D*)sigmapixthoff.Project3D("zy");

    normON  = hist->Integral();
    normOFF = histOFF->Integral();
    if (normON  != 0.0) hist->Scale(1.0/normON);
    if (normOFF != 0.0) histOFF->Scale(1.0/normOFF);

    // weighted average of ON and OFF distributions
    hist->Scale(facON);
    hist->Add(histOFF, facOFF); 

    for (Int_t k=1; k<=nbinspixel; k++)
      for (Int_t l=1; l<=nbinssigma; l++)
      {
        Double_t cont = hist->GetBinContent(k,l);
        fHSigmaPixTheta->SetBinContent(j, k, l, cont);  
      }

    delete hist;
    delete histOFF;


    //------------------------------------------------------------------
    // define target distribution  'number of blind pixels per event'
    ay = blindnthon.GetYaxis();
    Int_t nbinsn = ay->GetNbins();

    hist1    = fHBlindPixNThetaON->ProjectionY("", j, j, "");
    hist1->SetName("dummy");
    hist1OFF = fHBlindPixNThetaOFF->ProjectionY("", j, j, "");

    normON  = hist1->Integral();
    normOFF = hist1OFF->Integral();
    if (normON  != 0.0) hist1->Scale(1.0/normON);
    if (normOFF != 0.0) hist1OFF->Scale(1.0/normOFF);

    // sum of ON and OFF distributions
    hist1->Add(hist1OFF, 1.0); 
    Stat_t sum1 = hist1->Integral();
    if (sum1 != 0.0) hist1->Scale( 1.0/sum1 );

    for (Int_t k=1; k<=nbinsn; k++)
      {
        Double_t cont = hist1->GetBinContent(k);
        fHBlindPixNTheta->SetBinContent(j, k, cont);  
      }

    delete hist1;
    delete hist1OFF;

    //------------------------------------------------------------------
    // define target distribution  'id of blind pixel'
    ay = blindidthon.GetYaxis();
    Int_t nbinsid = ay->GetNbins();

    hist1    = fHBlindPixIdThetaON->ProjectionY("", j, j, "");
    hist1->SetName("dummy");
    hist1OFF = fHBlindPixIdThetaOFF->ProjectionY("", j, j, "");

    // divide by the number of events (from fHBlindPixNTheta)
    if (normON  != 0.0) hist1->Scale(1.0/normON);
    if (normOFF != 0.0) hist1OFF->Scale(1.0/normOFF);

    // sum of ON and OFF distributions
    hist1->Add(hist1OFF, 1.0); 

    for (Int_t k=1; k<=nbinsid; k++)
      {
        Double_t cont = hist1->GetBinContent(k);
        fHBlindPixIdTheta->SetBinContent(j, k, cont);  
      }

    delete hist1;
    delete hist1OFF;
  }
  //............   end of new loop over Theta bins   ....................

  *fLog << all 
        << "MPad::MergeONOFFMC(); The target distributions for the padding have been created" 
        << endl;
  *fLog << all << "----------------------------------------------------------" 
        << endl;
  //--------------------------------------------

  fHSigmaThetaMC->SetDirectory(NULL);
  fHSigmaThetaON->SetDirectory(NULL);
  fHSigmaThetaOFF->SetDirectory(NULL);

  fHDiffPixThetaMC->SetDirectory(NULL);
  fHDiffPixThetaON->SetDirectory(NULL);
  fHDiffPixThetaOFF->SetDirectory(NULL);

  fHSigmaPixThetaMC->SetDirectory(NULL);
  fHSigmaPixThetaON->SetDirectory(NULL);
  fHSigmaPixThetaOFF->SetDirectory(NULL);

  fHBlindPixIdThetaMC->SetDirectory(NULL);
  fHBlindPixIdThetaON->SetDirectory(NULL);
  fHBlindPixIdThetaOFF->SetDirectory(NULL);

  fHBlindPixNThetaMC->SetDirectory(NULL);
  fHBlindPixNThetaON->SetDirectory(NULL);
  fHBlindPixNThetaOFF->SetDirectory(NULL);

  fHgMC->SetDirectory(NULL);
  fHgON->SetDirectory(NULL);
  fHgOFF->SetDirectory(NULL);

  return kTRUE;
}

// --------------------------------------------------------------------------
//
// Update the matrix fHgA
// which contains the final padding prescriptions
//

Bool_t MPad::UpdateHg(TH2D *fHga, TH1D *histap, TH2D *fHge, TH3D *fHgA,
                      Int_t nbinssig, Int_t l)
{
  // histap  target distribution after ON-OFF merging
  // fHga    padding prescription for ON (or OFF) data to get to histap
  // fHge    padding prescription to get from histap to the target
  //         distribution after the ON-OFF-MC merging
  // fHgA    updated padding prescription for ON (or OFF) data to get
  //         from the original ON (or OFF) histogram to the target
  //         distribution after the ON-OFF-MC merging
  // l       Theta bin

  Double_t eps = 1.e-10;

  for (Int_t k=1; k<=nbinssig; k++)
  {
    Double_t na  = fHga->Integral(1, nbinssig, k, k, " ");
    Double_t ne  = fHge->Integral(k, k, 1, nbinssig, " ");
    Double_t nap = histap->GetBinContent(k);

    if (ne <= eps) 
    {
      // go to next k
      continue;
    }

    Double_t r = nap - na;

    if (r >= ne-eps)
    {
      for (Int_t j=k+1; j<=nbinssig; j++)
      {
        Double_t e = fHge->GetBinContent(k,j);
        Double_t A = fHgA->GetBinContent(l,k,j);
        fHgA->SetBinContent(l, k, j, A + e);
      }
      // go to next k
      continue;
    }

    for (Int_t j=k+1; j<=nbinssig; j++)
    {
      Double_t e = fHge->GetBinContent(k,j);
      Double_t A = fHgA->GetBinContent(l,k,j);
      fHgA->SetBinContent(l, k, j, A + r*e/ne);
    }

    if (na <= eps) 
    {
      // go to next k
      continue;
    }

    for (Int_t i=1; i<k; i++)
    {
      Double_t a = fHga->GetBinContent(i,k);
      Double_t A = fHgA->GetBinContent(l,i,k);
      fHgA->SetBinContent(l, i, k, A - (ne-r)*a/na);
    }

    for (Int_t i=1; i<=nbinssig; i++)
    {
      if (i == k) continue;
      for (Int_t j=i+1; j<=nbinssig; j++)
      {
        if (j == k) continue;
        Double_t a = fHga->GetBinContent(i,k);
        Double_t e = fHge->GetBinContent(k,j);
        Double_t A = fHgA->GetBinContent(l,i,j);
        fHgA->SetBinContent(l, i, j, A + (ne-r)*a*e/(na*ne)  );
      }
    }
  }

  return kTRUE; 
}

// --------------------------------------------------------------------------
//
// Merge the distributions 
//   fHSigmaTheta      2D-histogram  (Theta, sigmabar)
//
// ===>   of ON and MC data   <==============
//
// and define the matrices fHgMC and fHgON  
//
// These matrices tell which fraction of events should be padded 
// from bin k of sigmabar to bin j
//

Bool_t MPad::MergeONMC(
   TH2D& sigthmc,  TH3D& diffpixthmc,  TH3D& sigmapixthmc,
   TH2D& blindidthmc,  TH2D& blindnthmc,
   TH2D& sigthon,  TH3D& diffpixthon,  TH3D& sigmapixthon,
   TH2D& blindidthon,  TH2D& blindnthon)
{
  *fLog << all << "----------------------------------------------------------------------------------" << endl;
  *fLog << all << "MPad::MergeONMC(); Merge the ON  and MC histograms to obtain the target distributions for the padding"
        << endl;

  fHSigmaTheta = new TH2D( (TH2D&) sigthon );
  fHSigmaTheta->SetNameTitle("2D-ThetaSigmabar", "2D-ThetaSigmabar (target)");

  fHDiffPixTheta = new TH3D( (TH3D&) diffpixthon );
  fHDiffPixTheta->SetNameTitle("3D-ThetaPixDiff", "3D-ThetaPixDiff (target)");

  fHSigmaPixTheta = new TH3D( (TH3D&) sigmapixthon );
  fHSigmaPixTheta->SetNameTitle("3D-ThetaPixSigma", "3D-ThetaPixSigma (target)");

  fHBlindPixNTheta = new TH2D( (TH2D&) blindnthon );
  fHBlindPixNTheta->SetNameTitle("2D-ThetaBlindN", "2D-ThetaBlindN (target)");

  fHBlindPixIdTheta = new TH2D( (TH2D&) blindidthon );
  fHBlindPixIdTheta->SetNameTitle("2D-ThetaBlindId", "2D-ThetaBlindId (target)");

  //--------------------------
  fHSigmaThetaMC = &sigthmc;
  fHSigmaThetaMC->SetNameTitle("2D-ThetaSigmabarMC", "2D-ThetaSigmabarMC (orig.)");
  fHSigmaThetaON = &sigthon;
  fHSigmaThetaON->SetNameTitle("2D-ThetaSigmabarON", "2D-ThetaSigmabarON (orig.)");

  //--------------------------
  fHBlindPixNThetaMC = &blindnthmc;
  fHBlindPixNThetaMC->SetNameTitle("2D-ThetaBlindNMC", "2D-ThetaBlindNMC (orig.)");
  fHBlindPixNThetaON = &blindnthon;
  fHBlindPixNThetaON->SetNameTitle("2D-ThetaBlindNON", "2D-ThetaBlindNON (orig.)");

  //--------------------------
  fHBlindPixIdThetaMC = &blindidthmc;
  fHBlindPixIdThetaMC->SetNameTitle("2D-ThetaBlindIdMC", "2D-ThetaBlindIdMC (orig.)");
  fHBlindPixIdThetaON = &blindidthon;
  fHBlindPixIdThetaON->SetNameTitle("2D-ThetaBlindIdON", "2D-ThetaBlindIdON (orig.)");

  //--------------------------
  fHDiffPixThetaMC = &diffpixthmc;
  fHDiffPixThetaMC->SetNameTitle("3D-ThetaPixDiffMC", "3D-ThetaPixDiffMC (orig.)");
  fHDiffPixThetaON = &diffpixthon;
  fHDiffPixThetaON->SetNameTitle("3D-ThetaPixDiffON", "3D-ThetaPixDiffON (orig.)");

  //--------------------------
  fHSigmaPixThetaMC = &sigmapixthmc;
  fHSigmaPixThetaMC->SetNameTitle("3D-ThetaPixSigmaMC", "3D-ThetaPixSigmaMC (orig.)");
  fHSigmaPixThetaON = &sigmapixthon;
  fHSigmaPixThetaON->SetNameTitle("3D-ThetaPixSigmaON", "3D-ThetaPixSigmaON (orig.)");

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


  // get binning for fHgON and fHgMC  from sigthon
  // binning in Theta
  TAxis *ax = sigthon.GetXaxis();
  Int_t nbinstheta = ax->GetNbins();
  TArrayD edgesx;
  edgesx.Set(nbinstheta+1);
  for (Int_t i=0; i<=nbinstheta; i++)
  {
    edgesx[i] = ax->GetBinLowEdge(i+1);
    *fLog << all << "i, theta low edge = " << i << ",  " << edgesx[i] 
          << endl; 
  }
  MBinning binth;
  binth.SetEdges(edgesx);

  // binning in sigmabar
  TAxis *ay = sigthon.GetYaxis();
  Int_t nbinssig = ay->GetNbins();
  TArrayD edgesy;
  edgesy.Set(nbinssig+1);
  for (Int_t i=0; i<=nbinssig; i++)
  {
    edgesy[i] = ay->GetBinLowEdge(i+1); 
    *fLog << all << "i, sigmabar low edge = " << i << ",  " << edgesy[i] 
          << endl; 
  }
  MBinning binsg;
  binsg.SetEdges(edgesy);


  fHgMC = new TH3D;
  MH::SetBinning(fHgMC, &binth, &binsg, &binsg);
  fHgMC->SetNameTitle("3D-PaddingMatrixMC", "3D-PadMatrixThetaMC");

  fHgON = new TH3D;
  MH::SetBinning(fHgON, &binth, &binsg, &binsg);
  fHgON->SetNameTitle("3D-PaddingMatrixON", "3D-PadMatrixThetaON");


  //............   loop over Theta bins   ...........................



  *fLog << all << "MPad::MergeONMC(); bins of Theta, Sigmabarold, Sigmabarnew, fraction of events to be padded" << endl;
  for (Int_t l=1; l<=nbinstheta; l++)
  {
    TH2D * fHga = new TH2D;
    MH::SetBinning(fHga, &binsg, &binsg);
    TH2D * fHgb = new TH2D;
    MH::SetBinning(fHgb, &binsg, &binsg);

    //-------------------------------------------
    // merge ON and MC distributions
    // input : hista is the normalized 1D distr. of sigmabar for ON data
    //         histb is the normalized 1D distr. of sigmabar for MC data
    // output : histap    will be the merged distribution (ON-MC)
    //          fHga(k,j) will tell which fraction of ON events should be 
    //                    padded from bin k to bin j
    //          fHgb(k,j) will tell which fraction of MC events should be 
    //                    padded from bin k to bin j

    TH1D *hista  = sigthon.ProjectionY("sigon_y", l, l, "");
    Stat_t suma = hista->Integral();
    if (suma != 0.0) hista->Scale(1./suma);

    TH1D *histap  = new TH1D( (const TH1D&)*hista );

    TH1D *histb  = sigthmc.ProjectionY("sigmc_y", l, l, "");
    Stat_t sumb = histb->Integral();
    if (sumb != 0.0) histb->Scale(1./sumb);

    Merge2Distributions(hista, histb, histap, fHga, fHgb, nbinssig);
    delete hista;
    delete histb;

    // fill fHgON and fHgMC
    for (Int_t k=1; k<=nbinssig; k++)
      for (Int_t j=1; j<=nbinssig; j++)
      {
        Double_t a = fHga->GetBinContent(k,j);
        fHga->SetBinContent (k, j, 0.0);
        fHgON->SetBinContent(l, k, j, a);

        Double_t b = fHgb->GetBinContent(k,j);
        fHgb->SetBinContent (k, j, 0.0);
        fHgMC->SetBinContent(l, k, j, b);
      }


    // fill target distribution fHSigmaTheta 
    // (only for plotting, it will not  be used in the padding)
    for (Int_t j=1; j<=nbinssig; j++)
    {
      Double_t a = histap->GetBinContent(j);
      fHSigmaTheta->SetBinContent(l, j, a);
    }

    delete fHga;
    delete fHgb;

    delete histap;
  }
  //............   end of loop over Theta bins   ....................

  
  // Define the target distributions 
  //        fHDiffPixTheta, fHSigmaPixTheta
  //               (they are calculated as 
  //               averages of the ON and MCdistributions)
  //        fHBlindPixNTheta, fHBlindPixIdTheta
  //               (they are calculated as 
  //               the sum of the ON and MC distributions)
  // (fHDiffPixTheta will be used in the padding of MC events)

  //............   start of new loop over Theta bins   ....................
  for (Int_t j=1; j<=nbinstheta; j++)
  {
    // fraction of ON/MC events to be padded : fracON, fracMC

    Double_t fracON  = fHgON->Integral (j, j, 1, nbinssig, 1, nbinssig, "");
    Double_t fracMC  = fHgMC->Integral(j, j, 1, nbinssig, 1, nbinssig, "");

    Double_t normON;
    Double_t normMC;

    // set ranges for 2D-projections of 3D-histograms
    TAxis *ax = diffpixthon.GetXaxis();
    ax->SetRange(j, j);
    ax = diffpixthmc.GetXaxis();
    ax->SetRange(j, j);

    TH2D *hist;
    TH2D *histMC;

    TH1D *hist1;
    TH1D *hist1MC;

    // weights for ON and MC distrubtions when
    // calculating the weighted average
    Double_t facON = 0.5 * (1. - fracON + fracMC);
    Double_t facMC = 0.5 * (1. + fracON - fracMC);
  
    *fLog << all << "Theta bin j : fracON, fracMC, facON, facMC = " 
          << j << ",  " << fracON << ",  " << fracMC << ",  "
          << facON << ",  " << facMC << endl; 

    //------------------------------------------------------------------
    // define target distribution  'sigma^2-sigmavar^2 vs. pixel number'
    ay = diffpixthon.GetYaxis();
    Int_t nbinspixel = ay->GetNbins();

    TAxis *az = diffpixthon.GetZaxis();
    Int_t nbinsdiff = az->GetNbins();

    hist    = (TH2D*)diffpixthon.Project3D("zy");
    hist->SetName("dummy");
    histMC = (TH2D*)diffpixthmc.Project3D("zy");

    normON  = hist->Integral();
    normMC = histMC->Integral();
    if (normON  != 0.0) hist->Scale(1.0/normON);
    if (normMC != 0.0)  histMC->Scale(1.0/normMC);

    // weighted average of ON and MC distributions
    hist->Scale(facON);
    hist->Add(histMC, facMC); 

    for (Int_t k=1; k<=nbinspixel; k++)
      for (Int_t l=1; l<=nbinsdiff; l++)
      {
        Double_t cont = hist->GetBinContent(k,l);
        fHDiffPixTheta->SetBinContent(j, k, l, cont);  
      }

    delete hist;
    delete histMC;

    //------------------------------------------------------------------
    // define target distribution  'sigma vs. pixel number'
    ay = sigmapixthon.GetYaxis();
    nbinspixel = ay->GetNbins();

    az = sigmapixthon.GetZaxis();
    Int_t nbinssigma = az->GetNbins();

    hist    = (TH2D*)sigmapixthon.Project3D("zy");
    hist->SetName("dummy");
    histMC = (TH2D*)sigmapixthmc.Project3D("zy");

    normON  = hist->Integral();
    normMC = histMC->Integral();
    if (normON  != 0.0) hist->Scale(1.0/normON);
    if (normMC  != 0.0) histMC->Scale(1.0/normMC);

    // weighted average of ON and MC distributions
    hist->Scale(facON);
    hist->Add(histMC, facMC); 

    for (Int_t k=1; k<=nbinspixel; k++)
      for (Int_t l=1; l<=nbinssigma; l++)
      {
        Double_t cont = hist->GetBinContent(k,l);
        fHSigmaPixTheta->SetBinContent(j, k, l, cont);  
      }

    delete hist;
    delete histMC;


    //------------------------------------------------------------------
    // define target distribution  'number of blind pixels per event'
    ay = blindnthon.GetYaxis();
    Int_t nbinsn = ay->GetNbins();

    hist1    = fHBlindPixNThetaON->ProjectionY("", j, j, "");
    hist1->SetName("dummy");
    hist1MC  = fHBlindPixNThetaMC->ProjectionY("", j, j, "");

    normON  = hist1->Integral();
    normMC  = hist1MC->Integral();
    if (normON  != 0.0) hist1->Scale(1.0/normON);
    if (normMC  != 0.0) hist1MC->Scale(1.0/normMC);

    // sum of ON and MC distributions
    hist1->Add(hist1MC, 1.0); 
    Stat_t sum1 = hist1->Integral();
    if (sum1 != 0.0) hist1->Scale( 1.0/sum1 );

    for (Int_t k=1; k<=nbinsn; k++)
      {
        Double_t cont = hist1->GetBinContent(k);
        fHBlindPixNTheta->SetBinContent(j, k, cont);  
      }

    delete hist1;
    delete hist1MC;

    //------------------------------------------------------------------
    // define target distribution  'id of blind pixel'
    ay = blindidthon.GetYaxis();
    Int_t nbinsid = ay->GetNbins();

    hist1    = fHBlindPixIdThetaON->ProjectionY("", j, j, "");
    hist1->SetName("dummy");
    hist1MC = fHBlindPixIdThetaMC->ProjectionY("", j, j, "");

    // divide by the number of events (from fHBlindPixNTheta)
    if (normON  != 0.0) hist1->Scale(1.0/normON);
    if (normMC  != 0.0) hist1MC->Scale(1.0/normMC);

    // sum of ON and MC distributions
    hist1->Add(hist1MC, 1.0); 

    for (Int_t k=1; k<=nbinsid; k++)
      {
        Double_t cont = hist1->GetBinContent(k);
        fHBlindPixIdTheta->SetBinContent(j, k, cont);  
      }

    delete hist1;
    delete hist1MC;
  }
  //............   end of new loop over Theta bins   ....................

  *fLog << all 
        << "MPad::MergeONMC(); The target distributions for the padding have been created" 
        << endl;
  *fLog << all << "----------------------------------------------------------" 
        << endl;
  //--------------------------------------------

  fHSigmaThetaMC->SetDirectory(NULL);
  fHSigmaThetaON->SetDirectory(NULL);

  fHDiffPixThetaMC->SetDirectory(NULL);
  fHDiffPixThetaON->SetDirectory(NULL);

  fHSigmaPixThetaMC->SetDirectory(NULL);
  fHSigmaPixThetaON->SetDirectory(NULL);

  fHBlindPixIdThetaMC->SetDirectory(NULL);
  fHBlindPixIdThetaON->SetDirectory(NULL);

  fHBlindPixNThetaMC->SetDirectory(NULL);
  fHBlindPixNThetaON->SetDirectory(NULL);

  fHgMC->SetDirectory(NULL);
  fHgON->SetDirectory(NULL);

  return kTRUE;
}


// --------------------------------------------------------------------------
//
// Read padding distributions from a file
//
//
Bool_t MPad::ReadPaddingDist(const char* namefilein)
{
  *fLog << all << "MPad : Read padding histograms from file " 
        << namefilein << endl;

  fInfile = new TFile(namefilein);
  fInfile->ls();

    //------------------------------------
  
      /*
      fHSigmaThetaMC = 
      (TH2D*) gROOT->FindObject("2D-ThetaSigmabarMC");
      if (!fHSigmaThetaMC)
	{
          *fLog << all 
                << "MPad : Object '2D-ThetaSigmabarMC' not found on root file" 
                << endl;
          return kFALSE;
	}
      */  
      fHSigmaThetaMC = new TH2D; 
      fHSigmaThetaMC->Read("2D-ThetaSigmabarMC");
      *fLog << all 
            << "MPad : Object '2D-ThetaSigmabarMC' was read in" << endl;

      /*
      fHSigmaThetaON = 
      (TH2D*) gROOT->FindObject("2D-ThetaSigmabarON");
      if (!fHSigmaThetaON)
	{
          *fLog << all 
                << "MPad : Object '2D-ThetaSigmabarON' not found on root file" 
                << endl;
          return kFALSE;
	}
      */
      fHSigmaThetaON = new TH2D;
      fHSigmaThetaON->Read("2D-ThetaSigmabarON");
      *fLog << all 
            << "MPad : Object '2D-ThetaSigmabarON' was read in" << endl;

      /*
      fHSigmaThetaOFF = 
      (TH2D*) gROOT->FindObject("2D-ThetaSigmabarOFF");
      if (!fHSigmaThetaOFF)
	{
          *fLog << all 
                << "MPad : Object '2D-ThetaSigmabarOFF' not found on root file" 
                << endl;
          return kFALSE;
	}
      */
      fHSigmaThetaOFF = new TH2D;
      fHSigmaThetaOFF->Read("2D-ThetaSigmabarOFF");
      *fLog << all 
            << "MPad:Object '2D-ThetaSigmabarOFF' was read in" << endl;

    //------------------------------------
      /*
      fHDiffPixTheta = 
      (TH3D*) gROOT->FindObject("3D-ThetaPixDiff");
      if (!fHDiffPixTheta)
	{
          *fLog << all 
                << "MPad : Object '3D-ThetaPixDiff' not found on root file" 
                << endl;
          return kFALSE;
	}
      */
      fHDiffPixTheta = new TH3D;
      fHDiffPixTheta->Read("3D-ThetaPixDiff");
      *fLog << all 
            << "MPad : Object '3D-ThetaPixDiff' was read in" << endl;

      /*
      fHDiffPixThetaMC = 
      (TH3D*) gROOT->FindObject("3D-ThetaPixDiffMC");
      if (!fHDiffPixThetaMC)
	{
          *fLog << all 
                << "MPad : Object '3D-ThetaPixDiffMC' not found on root file" 
                << endl;
          return kFALSE;
	}
      */
      fHDiffPixThetaMC = new TH3D;
      fHDiffPixThetaMC->Read("3D-ThetaPixDiffMC");
      *fLog << all 
            << "MPad : Object '3D-ThetaPixDiffMC' was read in" << endl;

      /*
      fHDiffPixThetaON = 
      (TH3D*) gROOT->FindObject("3D-ThetaPixDiffON");
      if (!fHDiffPixThetaON)
	{
          *fLog << all 
                << "MPad : Object '3D-ThetaPixDiffON' not found on root file" 
                << endl;
          return kFALSE;
	}
      */
      fHDiffPixThetaON = new TH3D;
      fHDiffPixThetaON->Read("3D-ThetaPixDiffON");
      *fLog << all 
            << "MPad : Object '3D-ThetaPixDiffON' was read in" << endl;

      /*
      fHDiffPixThetaOFF = 
      (TH3D*) gROOT->FindObject("3D-ThetaPixDiffOFF");
      if (!fHDiffPixThetaOFF)
	{
          *fLog << all 
                << "MPad : Object '3D-ThetaPixDiffOFF' not found on root file" 
                << endl;
          return kFALSE;
	}
      */
      fHDiffPixThetaOFF = new TH3D;
      fHDiffPixThetaOFF->Read("3D-ThetaPixDiffOFF");
      *fLog << all 
            << "MPad : Object '3D-ThetaPixDiffOFF' was read in" << endl;


    //------------------------------------
      /*
       fHSigmaPixTheta = 
      (TH3D*) gROOT->FindObject("3D-ThetaPixSigma");
      if (!fHSigmaPixTheta)
	{
          *fLog << all 
                << "MPad : Object '3D-ThetaPixSigma' not found on root file" 
                << endl;
          return kFALSE;
	}
      */
      fHSigmaPixTheta = new TH3D;
      fHSigmaPixTheta->Read("3D-ThetaPixSigma");
      *fLog << all 
            << "MPad : Object '3D-ThetaPixSigma' was read in" << endl;

      /*
       fHSigmaPixThetaMC = 
      (TH3D*) gROOT->FindObject("3D-ThetaPixSigmaMC");
      if (!fHSigmaPixThetaMC)
	{
          *fLog << all 
                << "MPad : Object '3D-ThetaPixSigmaMC' not found on root file" 
                << endl;
          return kFALSE;
	}
      */
      fHSigmaPixThetaMC = new TH3D;
      fHSigmaPixThetaMC->Read("3D-ThetaPixSigmaMC");
      *fLog << all 
            << "MPad : Object '3D-ThetaPixSigmaMC' was read in" << endl;

      /*
      fHSigmaPixThetaON = 
      (TH3D*) gROOT->FindObject("3D-ThetaPixSigmaON");
      if (!fHSigmaPixThetaON)
	{
          *fLog << all 
                << "MPad : Object '3D-ThetaPixSigmaON' not found on root file" 
                << endl;
          return kFALSE;
	}
      */
      fHSigmaPixThetaON = new TH3D;
      fHSigmaPixThetaON->Read("3D-ThetaPixSigmaON");
      *fLog << all 
            << "MPad : Object '3D-ThetaPixSigmaON' was read in" << endl;

      /*
      fHSigmaPixThetaOFF = 
      (TH3D*) gROOT->FindObject("3D-ThetaPixSigmaOFF");
      if (!fHSigmaPixThetaOFF)
	{
          *fLog << all 
                << "MPad : Object '3D-ThetaPixSigmaOFF' not found on root file" 
                << endl;
          return kFALSE;
	}
      */
      fHSigmaPixThetaOFF = new TH3D;
      fHSigmaPixThetaOFF->Read("3D-ThetaPixSigmaOFF");
      *fLog << all 
            << "MPad : Object '3D-ThetaPixSigmaOFF' was read in" << endl;

    //------------------------------------
      /*
      fHgMC = 
      (TH3D*) gROOT->FindObject("3D-PaddingMatrixMC");
      if (!fHgMC)
	{
          *fLog << all 
                << "MPad : Object '3D-PaddingMatrixMC' not found on root file" 
                << endl;
          return kFALSE;
	}
      */
      fHgMC = new TH3D;
      fHgMC->Read("3D-PaddingMatrixMC");
      *fLog << all 
            << "MPad : Object '3D-PaddingMatrixMC' was read in" << endl;

      /*
      fHgON = 
      (TH3D*) gROOT->FindObject("3D-PaddingMatrixON");
      if (!fHgON)
	{
          *fLog << all 
                << "MPad : Object '3D-PaddingMatrixON' not found on root file" 
                << endl;
          return kFALSE;
	}
      */
      fHgON = new TH3D;
      fHgON->Read("3D-PaddingMatrixON");
      *fLog << all 
            << "MPad : Object '3D-PaddingMatrixON' was read in" << endl;

      /*
      fHgOFF = 
      (TH3D*) gROOT->FindObject("3D-PaddingMatrixOFF");
      if (!fHgOFF)
	{
          *fLog << all 
                << "MPad : Object '3D-PaddingMatrixOFF' not found on root file"
                << endl;
          return kFALSE;
	}
      */
      fHgOFF = new TH3D;
      fHgOFF->Read("3D-PaddingMatrixOFF");
      *fLog << all 
            << "MPad : Object '3D-PaddingMatrixOFF' was read in" << endl;

    //------------------------------------
      /*
      fHBlindPixIdTheta = 
      (TH2D*) gROOT->FindObject("2D-ThetaBlindId");
      if (!fHBlindPixIdTheta)
	{
          *fLog << all 
                << "MPad : Object '2D-ThetaBlindId' not found on root file" 
                << endl;
          return kFALSE;
	}
      */
      fHBlindPixIdTheta = new TH2D;
      fHBlindPixIdTheta->Read("2D-ThetaBlindId");
      *fLog << all 
            << "MPad : Object '2D-ThetaBlindId' was read in" << endl;

      /*
      fHBlindPixIdThetaMC = 
      (TH2D*) gROOT->FindObject("2D-ThetaBlindIdMC");
      if (!fHBlindPixIdThetaMC)
	{
          *fLog << all 
                << "MPad : Object '2D-ThetaBlindIdMC' not found on root file" 
                << endl;
          return kFALSE;
	}
      */
      fHBlindPixIdThetaMC = new TH2D;
      fHBlindPixIdThetaMC->Read("2D-ThetaBlindIdMC");
      *fLog << all 
            << "MPad : Object '2D-ThetaBlindIdMC' was read in" << endl;

      /*
      fHBlindPixIdThetaON = 
      (TH2D*) gROOT->FindObject("2D-ThetaBlindIdON");
      if (!fHBlindPixIdThetaON)
	{
          *fLog << all 
                << "MPad : Object '2D-ThetaBlindIdON' not found on root file" 
                << endl;
          return kFALSE;
	}
      */
      fHBlindPixIdThetaON = new TH2D;
      fHBlindPixIdThetaON->Read("2D-ThetaBlindIdON");
      *fLog << all 
            << "MPad : Object '2D-ThetaBlindIdON' was read in" << endl;

      /*
      fHBlindPixIdThetaOFF = 
      (TH2D*) gROOT->FindObject("2D-ThetaBlindIdOFF");
      if (!fHBlindPixIdThetaOFF)
	{
          *fLog << all 
                << "MPad : Object '2D-ThetaBlindIdOFF' not found on root file" 
                << endl;
          return kFALSE;
	}
      */
      fHBlindPixIdThetaOFF = new TH2D;
      fHBlindPixIdThetaOFF->Read("2D-ThetaBlindIdOFF");
      *fLog << all 
            << "MPad : Object '2D-ThetaBlindIdOFF' was read in" << endl;

    //------------------------------------
      /*
      fHBlindPixNTheta = 
      (TH2D*) gROOT->FindObject("2D-ThetaBlindN");
      if (!fHBlindPixNTheta)
	{
          *fLog << all 
                << "MPad : Object '2D-ThetaBlindN' not found on root file" 
                << endl;
          return kFALSE;
	}
      */
      fHBlindPixNTheta = new TH2D;
      fHBlindPixNTheta->Read("2D-ThetaBlindN");
      *fLog << all 
            << "MPad : Object '2D-ThetaBlindN' was read in" << endl;

      /*
      fHBlindPixNThetaMC = 
      (TH2D*) gROOT->FindObject("2D-ThetaBlindNMC");
      if (!fHBlindPixNThetaMC)
	{
          *fLog << all 
                << "MPad : Object '2D-ThetaBlindNMC' not found on root file" 
                << endl;
          return kFALSE;
	}
      */
      fHBlindPixNThetaMC = new TH2D;
      fHBlindPixNThetaMC->Read("2D-ThetaBlindNMC");
      *fLog << all 
            << "MPad : Object '2D-ThetaBlindNMC' was read in" << endl;

      /*
      fHBlindPixNThetaON = 
      (TH2D*) gROOT->FindObject("2D-ThetaBlindNON");
      if (!fHBlindPixNThetaON)
	{
          *fLog << all 
                << "MPad : Object '2D-ThetaBlindNON' not found on root file" 
                << endl;
          return kFALSE;
	}
      */
      fHBlindPixNThetaON = new TH2D;
      fHBlindPixNThetaON->Read("2D-ThetaBlindNON");
      *fLog << all 
            << "MPad : Object '2D-ThetaBlindNON' was read in" << endl;

      /*
      fHBlindPixNThetaOFF = 
      (TH2D*) gROOT->FindObject("2D-ThetaBlindNOFF");
      if (!fHBlindPixNThetaOFF)
	{
          *fLog << all 
                << "MPad : Object '2D-ThetaBlindNOFF' not found on root file" 
                << endl;
          return kFALSE;
	}
      */
      fHBlindPixNThetaOFF = new TH2D;
      fHBlindPixNThetaOFF->Read("2D-ThetaBlindNOFF");
      *fLog << all 
            << "MPad : Object '2D-ThetaBlindNOFF' was read in" << endl;

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

  return kTRUE;
}


// --------------------------------------------------------------------------
//
// Write padding distributions onto a file
//
//
Bool_t MPad::WritePaddingDist(const char* namefileout)
{
  *fLog << all << "MPad : Write padding histograms onto file " 
        << namefileout << endl;

  TFile outfile(namefileout, "RECREATE");

  fHBlindPixNTheta->Write();
  fHBlindPixNThetaMC->Write();
  fHBlindPixNThetaON->Write();
  fHBlindPixNThetaOFF->Write();

  fHBlindPixIdTheta->Write();
  fHBlindPixIdThetaMC->Write();
  fHBlindPixIdThetaON->Write();
  fHBlindPixIdThetaOFF->Write();

  fHSigmaThetaMC->Write();
  fHSigmaThetaON->Write();
  fHSigmaThetaOFF->Write();

  fHDiffPixTheta->Write();
  fHDiffPixThetaMC->Write();
  fHDiffPixThetaON->Write();
  fHDiffPixThetaOFF->Write();

  fHSigmaPixTheta->Write();
  fHSigmaPixThetaMC->Write();
  fHSigmaPixThetaON->Write();
  fHSigmaPixThetaOFF->Write();

  fHgMC->Write();
  fHgON->Write();
  fHgOFF->Write();

  *fLog << all 
        << "MPad::WriteTargetDist(); padding histograms written onto file "
        << namefileout << endl;

  return kTRUE;
}

// --------------------------------------------------------------------------
//
// Set type of data to be padded
//
//     this is not necessary if the type of the data can be recognized
//     directly from the input
//
//
void MPad::SetDataType(const char* type)
{
  fType = type;
  *fLog << all << "MPad::SetDataType(); type of data to be padded : " 
        << fType << endl; 

  return;
}

// --------------------------------------------------------------------------
//
//  Set the pointers and prepare the histograms
//
Int_t MPad::PreProcess(MParList *pList)
{
  if ( !fHSigmaThetaMC     || !fHSigmaThetaON     || !fHSigmaThetaOFF      ||  
       !fHDiffPixThetaMC   || !fHDiffPixThetaON   || !fHDiffPixThetaOFF    ||
       !fHBlindPixIdThetaMC|| !fHBlindPixIdThetaON|| !fHBlindPixIdThetaOFF ||
       !fHBlindPixNThetaMC || !fHBlindPixNThetaON || !fHBlindPixNThetaOFF  ||
       !fHgMC              || !fHgON              || !fHgOFF                 )
  { 
       *fLog << err 
             << "MPad : At least one of the histograms needed for the padding is not defined ... aborting." 
             << endl;
       return kFALSE;
  }

  fPointPos = (MPointingPos*)pList->FindObject("MPointingPos");
  if (!fPointPos)
    {
       *fLog << err << dbginf << "MPad : MPointingPos not found... aborting." 
             << endl;
       return kFALSE;
    }
  
   fPed = (MPedPhotCam*)pList->FindObject("MPedPhotCam");
   if (!fPed)
     {
       *fLog << err << "MPad : MPedPhotCam not found... aborting." 
             << endl;
       return kFALSE;
     }

   fCam = (MGeomCam*)pList->FindObject("MGeomCam");
   if (!fCam)
     {
       *fLog << err 
             << "MPad : MGeomCam not found (no geometry information available)... aborting." 
             << endl;
       return kFALSE;
     }
  
   fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt");
   if (!fEvt)
     {
       *fLog << err << "MPad : MCerPhotEvt not found... aborting." 
             << endl;
       return kFALSE;
     }

   fSigmabar = (MSigmabar*)pList->FindCreateObj("MSigmabar");
   if (!fSigmabar)
     {
       *fLog << err << "MPad : MSigmabar not found... aborting." << endl;
       return kFALSE;
     }

   fBlindPix = (MBlindPixels*)pList->FindCreateObj("MBlindPixels");
   if (!fBlindPix)
     {
       *fLog << err << "MPad : MBlindPixels not found... aborting." 
             << endl;
       return kFALSE;
     }
   
   if (fType !="ON"  &&  fType !="OFF"  &&  fType !="MC")
     {
       *fLog << err 
             << "MPad : Type of data to be padded not defined... aborting." 
             << endl;
       return kFALSE;
     }


   //--------------------------------------------------------------------
   // histograms for checking the padding
   //
   // plot pedestal sigmas
   fHSigmaPedestal = new TH2D("SigPed","Sigma: after vs. before padding", 
                     100, 0.0, 75.0, 100, 0.0, 75.0);
   fHSigmaPedestal->SetXTitle("Pedestal sigma before padding");
   fHSigmaPedestal->SetYTitle("Pedestal sigma after padding");

   // plot no.of photons (before vs. after padding) 
   fHPhotons = new TH2D("Photons","Photons: after vs.before padding", 
                        100, -100.0, 300.0, 100, -100, 300);
   fHPhotons->SetXTitle("no.of photons before padding");
   fHPhotons->SetYTitle("no.of photons after padding");

   // plot of added NSB
   fHNSB = new TH1D("NSB","Distribution of added NSB", 100, -100.0, 200.0);
   fHNSB->SetXTitle("no.of added NSB photons");
   fHNSB->SetYTitle("no.of pixels");


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

   fIter = 10;

   memset(fErrors,   0, sizeof(fErrors));
   memset(fWarnings, 0, sizeof(fWarnings));

   return kTRUE;
}

// --------------------------------------------------------------------------
//
// Do the Padding (noise adjustment)
//
// input for the padding : 
//  - the matrices fHgON, fHgOFF, fHgMC
//  - the original distributions fHBlindPixNTheta, fHBlindPixIdTheta
//                               fHSigmaTheta, fHDiffPixTheta
//

Int_t MPad::Process()
{
  //*fLog << all << "Entry MPad::Process();" << endl;

  // this is the conversion factor from the number of photons
  //                                 to the number of photo electrons
  // later to be taken from MCalibration
  // fPEperPhoton = PW * LG * QE * 1D
  Double_t fPEperPhoton = 0.92 * 0.94 * 0.23 * 0.9;

  //------------------------------------------------
  // Add pixels to MCerPhotEvt which are not yet in;
  // set their number of photons equal to zero. 
  // Purpose : by the padding the number of photons is changed
  // so that cleaning and cuts may be affected.
  // However, no big effect is expectec unless the padding is strong
  // (strong increase of the noise level)
  // At present comment out this code

  /*
  UInt_t imaxnumpix = fCam->GetNumPixels();

  for (UInt_t i=0; i<imaxnumpix; i++)
  {
    Bool_t alreadythere = kFALSE;
    UInt_t npix = fEvt->GetNumPixels();
    for (UInt_t j=0; j<npix; j++)
    {
      MCerPhotPix &pix = (*fEvt)[j];
      UInt_t id = pix.GetPixId();
      if (i==id)
      {
        alreadythere = kTRUE;
        break;
      }
    }
    if (alreadythere)
      continue;

    fEvt->AddPixel(i, 0.0, (*fPed)[i].GetRms());
  }
  */


  //-----------------------------------------
  Int_t rc=0;
  Int_t rw=0;

  const UInt_t npix = fEvt->GetNumPixels();

  //-------------------------------------------
  // Calculate average sigma of the event
  //

  fSigmabar->Calc(*fCam, *fPed, *fEvt);
  //*fLog << all << "MPad::Process(); before padding : " << endl;
  //fSigmabar->Print("");

  // inner sigmabar/sqrt(area)
  Double_t sigbarInnerold_phot = fSigmabar->GetSigmabarInner();
  Double_t sigbarInnerold  = sigbarInnerold_phot * fPEperPhoton;
  Double_t sigbarInnerold2 = sigbarInnerold*sigbarInnerold;

  // outer sigmabar/sqrt(area)
  Double_t sigbarOuterold_phot = fSigmabar->GetSigmabarOuter();
  Double_t sigbarOuterold  = sigbarOuterold_phot * fPEperPhoton;
  Double_t sigbarOuterold2 = sigbarOuterold*sigbarOuterold;

  const Double_t theta = fPointPos->GetZd();

  //*fLog << all << "theta = " << theta << endl;

  Int_t binTheta = fHBlindPixNTheta->GetXaxis()->FindBin(theta);
  if ( binTheta < 1  ||  binTheta > fHBlindPixNTheta->GetNbinsX() )
  {
    *fLog << warn 
          << "MPad::Process(); binNumber out of range : theta, binTheta = "
          << theta << ",  " << binTheta << ";  aborting " << endl;

    rc = 2;
    fErrors[rc]++;
    return kCONTINUE;
  }

  //*fLog << all << "binTheta = " << binTheta << endl;
  

  //-------------------------------------------
  // get number of events in this theta bin
  // and number of events in this sigbarInnerold_phot bin

  Double_t nTheta;
  Double_t nSigma;

  TH2D *st = NULL;
  if      (fType == "MC")  st = fHSigmaThetaMC;
  else if (fType == "ON")  st = fHSigmaThetaON;
  else if (fType == "OFF") st = fHSigmaThetaOFF;
  else
  {
    *fLog << err << "MPad : illegal data type '" << fType 
          <<  "'" << endl;
    return kERROR;  
  }

  TH1D *hn;
  hn = st->ProjectionY("", binTheta, binTheta, "");

  //*fLog << "Title of histogram : " << st->GetTitle() << endl;
  //for (Int_t i=1; i<=hn->GetNbinsX(); i++)
  //{
  //  *fLog << "bin " << i << " : no.of entries = " << hn->GetBinContent(i)
  //        << endl;
  //}

  nTheta = hn->Integral();
  Int_t binSigma = hn->FindBin(sigbarInnerold_phot);
  nSigma = hn->GetBinContent(binSigma);

  //*fLog << all           
  //       << "Theta, sigbarold_phot, binTheta, binSigma, nTheta, nSigma = "
  //       << theta << ",  " << sigbarInnerold_phot << ",  " << binTheta 
  //       << ",  "
  //       << binSigma << ",  " << nTheta << ",  " << nSigma   << endl;      

  delete hn;
 
  //-------------------------------------------
  // for the current theta :
  //     generate blind pixels according to the histograms 
  //     fHBlindPixNTheta and fHBlindPixIDTheta
  //
  // ON  : add the blind pixels from the OFF data
  // OFF : add the blind pixels from the ON data
  // MC  : add the blind pixels from the ON and OFF data

  //-----------------------------------
  if (fType == "ON"  ||  fType == "MC")
  {
    // numBlind is the number of blind pixels in an event
    TH1D *nblind;
    UInt_t numBlind;

    nblind = fHBlindPixNThetaOFF->ProjectionY("", binTheta, binTheta, "");
    if ( nblind->Integral() == 0.0 )
    {
      *fLog << warn << "MPad::Process(); projection of '"
            << fHBlindPixNThetaOFF->GetName() << "' for Theta bin " 
            << binTheta << " has no entries; aborting " << endl;
      // event cannot be padded; skip event
      delete nblind;

      rc = 7;
      fErrors[rc]++;
      return kCONTINUE;
    }
    else
    {
      numBlind = (Int_t) (nblind->GetRandom()+0.5);
    }
    delete nblind;


    // throw the Id of numBlind different pixels in this event
    if ( numBlind > 0)
    {
      TH1D *hblind;
      UInt_t idBlind;
      UInt_t listId[npix];
      UInt_t nlist = 0;
      Bool_t equal;

      hblind = fHBlindPixIdThetaOFF->ProjectionY("", binTheta, binTheta, "");
      if ( hblind->Integral() > 0.0 )
        for (UInt_t i=0; i<numBlind; i++)
        {
          while(1)
          {
            idBlind = (Int_t) (hblind->GetRandom()+0.5);
            equal = kFALSE;
            for (UInt_t j=0; j<nlist; j++)
              if (idBlind == listId[j])
  	      { 
                equal = kTRUE;
                break;
	      }
            if (!equal) break;
          }
          listId[nlist] = idBlind;
          nlist++;

          fBlindPix->SetPixelBlind(idBlind);
          //*fLog << all << "idBlind = " << idBlind << endl;
        }
      fBlindPix->SetReadyToSave();

      delete hblind;
    }
  }

  //------------------------------------
  if (fType == "OFF"  ||  fType == "MC")
  {
    // throw numBlind;
    // numBlind is the number of blind pixels in an event
    TH1D *nblind;
    UInt_t numBlind;

    nblind = fHBlindPixNThetaON->ProjectionY("", binTheta, binTheta, "");
    if ( nblind->Integral() == 0.0 )
    {
      *fLog << warn << "MPad::Process(); projection of '" 
            << fHBlindPixNThetaON->GetName() << "' for Theta bin " 
            << binTheta << " has no entries; Skip event " << endl;
      delete nblind;

      rc = 7;
      fErrors[rc]++;
      return kCONTINUE;
    }
    else
    {
      numBlind = (Int_t) (nblind->GetRandom()+0.5);
    }
    delete nblind;


    // throw the Id of numBlind different pixels in this event
    if ( numBlind > 0)
    {
      TH1D *hblind;
      UInt_t idBlind;
      UInt_t listId[npix];
      UInt_t nlist = 0;
      Bool_t equal;

      hblind = fHBlindPixIdThetaON->ProjectionY("", binTheta, binTheta, "");
      if ( hblind->Integral() > 0.0 )
        for (UInt_t i=0; i<numBlind; i++)
        {
          while(1)
          {
            idBlind = (Int_t) (hblind->GetRandom()+0.5);
            equal = kFALSE;
            for (UInt_t j=0; j<nlist; j++)
              if (idBlind == listId[j])
  	      { 
                equal = kTRUE;
                break;
	      }
            if (!equal) break;
          }
          listId[nlist] = idBlind;
          nlist++;

          fBlindPix->SetPixelBlind(idBlind);
          //*fLog << all << "idBlind = " << idBlind << endl;
        }
      fBlindPix->SetReadyToSave();

      delete hblind;
    }
  }

  //******************************************************************
  // has the event to be padded ?
  // if yes : to which sigmabar should it be padded ?
  //

  Int_t binSig = fHgON->GetYaxis()->FindBin(sigbarInnerold_phot);
  //*fLog << all << "binSig, sigbarInnerold_phot = " << binSig << ",  " 
  //        << sigbarInnerold_phot << endl;

  Double_t prob;
  TH1D *hpad = NULL;

  TH3D *Hg = NULL;
  if      (fType == "MC")  Hg = fHgMC;
  else if (fType == "ON")  Hg = fHgON;
  else if (fType == "OFF") Hg = fHgOFF;
  else
  {
    *fLog << err << "MPad : illegal type of data '" << fType << "'"
          << endl;
    return kERROR;
  }

  hpad = Hg->ProjectionZ("zON", binTheta, binTheta, binSig, binSig, "");

  //Int_t nb = hpad->GetNbinsX();
  //for (Int_t i=1; i<=nb; i++)
  //{
  //  if (hpad->GetBinContent(i) != 0.0)
  //    *fLog << all << "i, fHgON = " << i << ",  " 
  //          << hpad->GetBinContent(i) << endl;
  //}

  if (nSigma != 0.0)
    prob = hpad->Integral() * nTheta/nSigma;
  else
    prob = 0.0;

  //*fLog << all << "nTheta, nSigma, prob = " << nTheta << ",  " << nSigma 
  //       << ",  " << prob << endl;

  if ( prob <= 0.0  ||  gRandom->Uniform() > prob )
  {
    delete hpad;
    // event should not be padded
    //*fLog << all << "event is not padded" << endl;

    rc = 8;
    fErrors[rc]++;
    return kTRUE;
  }
  // event should be padded
  //*fLog << all << "event will be padded" << endl;  


  //-------------------------------------------
  // for the current theta, generate a sigmabar :
  //     for MC/ON/OFF data according to the matrix fHgMC/ON/OFF
  //
  Double_t sigmabarInner_phot = 0;
  Double_t sigmabarInner      = 0;

  
  //Int_t nbinsx = hpad->GetNbinsX();
  //for (Int_t i=1; i<=nbinsx; i++)
  //{
  //  if (hpad->GetBinContent(i) != 0.0)
  //    *fLog << all << "i, fHg = " << i << ",  " << hpad->GetBinContent(i) 
  //          << endl;
  //}

  sigmabarInner_phot = hpad->GetRandom();
  sigmabarInner = sigmabarInner_phot * fPEperPhoton;

  //*fLog << all << "sigmabarInner_phot = " << sigmabarInner_phot << endl;

  delete hpad;
  
  // new inner sigmabar2/area
  const Double_t sigmabarInner2 = sigmabarInner*sigmabarInner;

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

  //*fLog << all << "MPad::Process(); sigbarInnerold, sigmabarInner = " 
  //      << sigbarInnerold << ",  "<< sigmabarInner << endl;

  // Stop if target sigmabar is <= sigbarold
  if (sigmabarInner <= sigbarInnerold)
  {
    *fLog << all << "MPad::Process(); target sigmabarInner is less than sigbarInnerold : "
          << sigmabarInner << ",  " << sigbarInnerold << ",   aborting" 
          << endl;

    rc = 4;
    fErrors[rc]++;
    return kCONTINUE;     
  }


  //-------------------------------------------
  //
  // Calculate average number of NSB photo electrons to be added (lambdabar)
  // from the value of sigmabar, 
  //  - making assumptions about the average electronic noise (elNoise2) and
  //  - using a fixed value (F2excess)  for the excess noise factor
  
  Double_t elNoise2;         // [photo electrons]  
  Double_t F2excess  = 1.3;
  Double_t lambdabar;        // [photo electrons]


  Int_t bincheck = fHDiffPixTheta->GetXaxis()->FindBin(theta);
  if (binTheta != bincheck)
  {
    *fLog << err 
          << "MPad::Process(); The binnings of the 2 histograms are not identical; aborting"
          << endl;
    return kERROR;
  }

  // In this Theta bin, get RMS of (Sigma^2-sigmabar^2) for inner pixels
  // The average electronic noise (to be added) has to be in the order of
  // this RMS, otherwise the electronic noise of an individual pixel 
  // (elNoise2Pix) may become negative

  // Attention : maximum ID of inner pixels hard coded !!!
  Int_t idmaxpixinner = 396;
  Int_t binmaxpixinner = 
        fHDiffPixTheta->GetYaxis()->FindBin( (Double_t)idmaxpixinner );

  TH1D *hnoise;
  if (fType == "MC")
    hnoise = fHDiffPixTheta->ProjectionZ("", binTheta, binTheta, 0, binmaxpixinner, "");
  else if (fType == "ON")
    hnoise = fHDiffPixThetaOFF->ProjectionZ("", binTheta, binTheta, 0, binmaxpixinner, "");
  else if (fType == "OFF")
    hnoise = fHDiffPixThetaON->ProjectionZ("", binTheta, binTheta, 0, binmaxpixinner, "");
  else
  {
    *fLog << err << "MPad::Process(); illegal data type... aborting" 
          << endl;
    return kERROR;
  }

  Double_t RMS_phot = hnoise->GetRMS(1);  
  Double_t RMS = RMS_phot * fPEperPhoton * fPEperPhoton;
  delete hnoise;

  elNoise2 = TMath::Min(2.0*RMS,  sigmabarInner2 - sigbarInnerold2);
  //*fLog << all << "elNoise2 = " << elNoise2 << endl; 

  lambdabar = (sigmabarInner2 - sigbarInnerold2 - elNoise2) / F2excess;  

  // This value of lambdabar is the same for all pixels;
  // note that lambdabar is normalized to the area of pixel 0

  //----------   start loop over pixels   ---------------------------------
  // do the padding for each pixel
  //
  // pad only pixels   - which are used and not blind (before image cleaning)
  //

  Double_t sigma2      = 0;

  Double_t diff_phot   = 0;
  Double_t diff        = 0;

  Double_t addSig2_phot= 0;
  Double_t addSig2     = 0;

  Double_t elNoise2Pix = 0;


  // throw the Sigma for the pixels from the distribution fHDiffPixTheta
  // MC  : from fHDiffPixTheta
  // ON  : from fHDiffPixThetaOFF
  // OFF : from fHDiffPixThetaON

  TH3D *sp = NULL;

  if      (fType == "MC")  sp = fHDiffPixTheta;
  else if (fType == "ON")  sp = fHDiffPixThetaOFF;
  else if (fType == "OFF") sp = fHDiffPixThetaON;
  else
  {
    *fLog << err << "MPad::Process(); illegal data type... aborting" 
          << endl;
    return kERROR;
  }

  Double_t sigbarold2;
  Double_t sigmabar2;
  Double_t sigmabarOuter2;


  for (UInt_t i=0; i<npix; i++) 
  {   

    MCerPhotPix &pix = (*fEvt)[i];
    if ( !pix.IsPixelUsed() )
      continue;

    //if ( pix.GetNumPhotons() == 0.0)
    //{
    //  *fLog << warn 
    //        << "MPad::Process(); no.of photons is 0 for used pixel" 
    //        << endl;
    //  continue;
    //}

    Int_t j = pix.GetPixId();

    // GetPixRatio returns (area of pixel 0 / area of current pixel)
    Double_t ratioArea = 1.0 / fCam->GetPixRatio(j);

    if (ratioArea < 1.5)
    {
      sigbarold2 = sigbarInnerold2;
      sigmabar2  = sigmabarInner2;
    }
    else
    {
      sigbarold2    = sigbarOuterold2;
      //sigbarOuter2  = sigmabarInner2 * sigbarOuterold2 / sigbarInnerold2;
      sigmabarOuter2  = sigmabarInner2 + (sigbarOuterold2 - sigbarInnerold2);
      sigmabar2     = sigmabarOuter2;
    }

    MPedPhotPix &ppix = (*fPed)[j];

    if ( fBlindPix != NULL  &&  fBlindPix->IsBlind(j) )
    {
       // this should never occur, because blind pixels should have
       // been set unused by MBlindPixelsCalc2::UnMap()
       //*fLog << all << "MPad::Process; blind pixel found which is used, j = "
       //       << j << "... go to next pixel." << endl;
       continue;
    }	

    // count number of pixels treated
    fWarnings[0]++;


    Double_t oldsigma_phot = ppix.GetRms();
    Double_t oldsigma = oldsigma_phot * fPEperPhoton;
    Double_t oldsigma2 = oldsigma*oldsigma / ratioArea;

    //---------------------------------
    // throw the Sigma for this pixel 
    //
    Int_t binPixel = fHDiffPixTheta->GetYaxis()->FindBin( (Double_t)j );

    Int_t count;
    Bool_t ok;

    TH1D *hdiff;


     hdiff = sp->ProjectionZ("", binTheta, binTheta,
                                 binPixel, binPixel, "");
     Double_t integral =  hdiff->Integral();
     // if there are no entries in hdiff, diff cannot be thrown
     // in this case diff will be set to zero
     if ( integral == 0 )
     {
       //*fLog << warn << "MPad::Process(); fType = " << fType 
       //      << ", projection of '"
       //      << sp->GetName() << "' for Theta bin " 
       //      << binTheta << "  and pixel " << j  
       //      << " has no entries; set diff equal to the old difference  " 
       //      << endl;

       diff = TMath::Max(oldsigma2 - sigbarold2,
                         lambdabar*F2excess + oldsigma2 - sigmabar2);

       rc = 2;
       fWarnings[rc]++;
     }

     // start of else -------------------
     else
     {
       count = 0;
       ok = kFALSE;
       for (Int_t m=0; m<fIter; m++)
       {
         count += 1;
         diff_phot = hdiff->GetRandom();

         //*fLog << "after GetRandom : j, m, diff_phot = " << j << " : "
         //      << m << ",  " << diff_phot << endl;

         diff = diff_phot * fPEperPhoton * fPEperPhoton;
 
         // the following condition ensures that elNoise2Pix > 0.0 
         if ( (diff + sigmabar2 - oldsigma2
                                - lambdabar*F2excess) > 0.0 )
         {
           ok = kTRUE;
           break;
         }
       }

       if (!ok)
       {
         //*fLog << all << "theta, j, count, sigmabar2, diff, oldsigma2, ratioArea, lambdabar = " 
         //      << theta << ",  " 
         //      << j << ",  " << count << ",  " << sigmabar2 << ",  " 
         //      << diff << ",  " << oldsigma2 << ",  " << ratioArea << ",  "
         //      << lambdabar << endl;
         diff = lambdabar*F2excess + oldsigma2 - sigmabar2; 

         rw = 1;
         fWarnings[rw]++;
       }
     }
     // end of else --------------------

    delete hdiff;
    sigma2 = diff + sigmabar2;


    //---------------------------------
    // get the additional sigma^2 for this pixel (due to the padding)

    addSig2 = (sigma2 - oldsigma2) * ratioArea;
    addSig2_phot = addSig2 / (fPEperPhoton * fPEperPhoton);

    //---------------------------------
    // get the additional electronic noise for this pixel

    elNoise2Pix = addSig2 - lambdabar*F2excess*ratioArea;


    //---------------------------------
    // throw actual number of additional NSB photo electrons (NSB)
    //       and its RMS (sigmaNSB) 

    Double_t NSB0 = gRandom->Poisson(lambdabar*ratioArea);
    Double_t arg  = NSB0*(F2excess-1.0) + elNoise2Pix;
    Double_t sigmaNSB0;

    if (arg >= 0.0)
    {
      sigmaNSB0 = sqrt( arg );
    }
    else
    {
      sigmaNSB0 = 0.0000001;      
      if (arg < -1.e-10)
      {
        *fLog << warn << "MPad::Process(); argument of sqrt < 0 : "
              << arg << endl;
      }
    }


    //---------------------------------
    // smear NSB0 according to sigmaNSB0
    // and subtract lambdabar because of AC coupling

    Double_t NSB = gRandom->Gaus(NSB0, sigmaNSB0) - lambdabar*ratioArea;
    Double_t NSB_phot = NSB / fPEperPhoton;

    //---------------------------------
 
    // add additional NSB to the number of photons
    Double_t oldphotons_phot = pix.GetNumPhotons();
    Double_t oldphotons = oldphotons_phot * fPEperPhoton;
    Double_t newphotons = oldphotons + NSB;
    Double_t newphotons_phot = newphotons / fPEperPhoton;    
    pix.SetNumPhotons(	newphotons_phot );


    fHNSB->Fill( NSB_phot/ratioArea );
    fHPhotons->Fill( oldphotons_phot/ratioArea, 
                     newphotons_phot/ratioArea );


    // error: add sigma of padded noise quadratically
    Double_t olderror_phot = pix.GetErrorPhot();
    Double_t newerror_phot = 
                           sqrt( olderror_phot*olderror_phot + addSig2_phot );
    pix.SetErrorPhot( newerror_phot );


    Double_t newsigma = sqrt( sigma2 * ratioArea ); 
    Double_t newsigma_phot = newsigma / fPEperPhoton; 
    ppix.SetRms( newsigma_phot );

    fHSigmaPedestal->Fill( oldsigma_phot, newsigma_phot );
  } 
  //----------   end of loop over pixels   ---------------------------------

  // Calculate sigmabar again and crosscheck

  fSigmabar->Calc(*fCam, *fPed, *fEvt);
  //*fLog << all << "MPad::Process(); after padding : " << endl;
  //fSigmabar->Print("");

  //*fLog << all << "Exit MPad::Process();" << endl;

  rc = 0;
  fErrors[rc]++;
  return kTRUE;
  //******************************************************************
}

// --------------------------------------------------------------------------
//
//
Int_t MPad::PostProcess()
{
    if (GetNumExecutions() != 0)
    {

    *fLog << inf << endl;
    *fLog << GetDescriptor() << " execution statistics:" << endl;
    *fLog << dec << setfill(' ');

    if (fWarnings[0] == 0 ) fWarnings[0] = 1;

    *fLog << " " << setw(7) << fWarnings[1] << " (" << setw(3) 
          << (int)(fWarnings[1]*100/fWarnings[0])
          << "%) Pixel: iteration to find acceptable sigma failed" 
          << ", fIter = " << fIter << endl;

    *fLog << " " << setw(7) << fWarnings[2] << " (" << setw(3) 
          << (int)(fWarnings[2]*100/fWarnings[0]) 
          << "%) Pixel: No data for generating Sigma^2-Sigmabar^2" << endl;


    *fLog << " " << setw(7) << fErrors[2] << " (" << setw(3) 
          << (int)(fErrors[2]*100/GetNumExecutions()) 
          << "%) Evts skipped due to: Zenith angle out of range" << endl;

    *fLog << " " << setw(7) << fErrors[4] << " (" << setw(3) 
          << (int)(fErrors[4]*100/GetNumExecutions()) 
          << "%) Evts skipped due to: Target sigma <= Sigmabar_old" << endl;

    *fLog << " " << setw(7) << fErrors[7] << " (" << setw(3) 
          << (int)(fErrors[7]*100/GetNumExecutions()) 
          << "%) Evts skipped due to: No data for generating Blind pixels" << endl;

    *fLog << " " << setw(7) << fErrors[8] << " (" << setw(3) 
          << (int)(fErrors[8]*100/GetNumExecutions()) 
          << "%) Evts didn't have to be padded" << endl;

    *fLog << " " << fErrors[0] << " (" 
          << (int)(fErrors[0]*100/GetNumExecutions()) 
          << "%) Evts were successfully padded!" << endl;
    *fLog << endl;

    }

    //---------------------------------------------------------------
    TCanvas &c = *(MH::MakeDefCanvas("Pad", "", 900, 600)); 
    c.Divide(3, 2);
    gROOT->SetSelectedPad(NULL);

    c.cd(1);
    fHNSB->SetDirectory(NULL);
    fHNSB->DrawCopy();
    fHNSB->SetBit(kCanDelete);    

    c.cd(2);
    fHSigmaPedestal->SetDirectory(NULL);
    fHSigmaPedestal->DrawCopy();
    fHSigmaPedestal->SetBit(kCanDelete);    

    c.cd(3);
    fHPhotons->SetDirectory(NULL);
    fHPhotons->DrawCopy();
    fHPhotons->SetBit(kCanDelete);    

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

    /*
    c.cd(4);
    fHSigmaTheta->SetDirectory(NULL);
    fHSigmaTheta->SetTitle("(Target) 2D : Sigmabar, \\Theta");
    fHSigmaTheta->DrawCopy();     
    fHSigmaTheta->SetBit(kCanDelete);     
    */

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

    /*
    c.cd(7);
    fHBlindPixNTheta->SetDirectory(NULL);
    fHBlindPixNTheta->SetTitle("(Target) 2D : no.of blind pixels, \\Theta");
    fHBlindPixNTheta->DrawCopy();     
    fHBlindPixNTheta->SetBit(kCanDelete);     
    */

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

    /*
    c.cd(10);
    fHBlindPixIdTheta->SetDirectory(NULL);
    fHBlindPixIdTheta->SetTitle("(Target) 2D : blind pixel Id, \\Theta");
    fHBlindPixIdTheta->DrawCopy();     
    fHBlindPixIdTheta->SetBit(kCanDelete);     
    */

    //--------------------------------------------------------------------
    // draw the 3D histogram (target): Theta, pixel, Sigma^2-Sigmabar^2

    /*
    c.cd(5);
    TH2D *l1 = NULL;
    l1 = (TH2D*) ((TH3*)fHDiffPixTheta)->Project3D("zx");
    l1->SetDirectory(NULL);
    l1->SetTitle("(Target) Sigma^2-Sigmabar^2 vs. \\Theta (all pixels)");
    l1->SetXTitle("\\Theta [\\circ]");
    l1->SetYTitle("Sigma^2-Sigmabar^2");

    l1->DrawCopy("box");
    l1->SetBit(kCanDelete);;

    c.cd(8);
    TH2D *l2 = NULL;
    l2 = (TH2D*) ((TH3*)fHDiffPixTheta)->Project3D("zy");
    l2->SetDirectory(NULL);
    l2->SetTitle("(Target) Sigma^2-Sigmabar^2 vs. pixel number (all \\Theta)");
    l2->SetXTitle("pixel");
    l2->SetYTitle("Sigma^2-Sigmabar^2");

    l2->DrawCopy("box");
    l2->SetBit(kCanDelete);;
    */

    //--------------------------------------------------------------------
    // draw the 3D histogram (target): Theta, pixel, Sigma

    /*
    c.cd(6);
    TH2D *k1 = NULL;
    k1 = (TH2D*) ((TH3*)fHSigmaPixTheta)->Project3D("zx");
    k1->SetDirectory(NULL);
    k1->SetTitle("(Target) Sigma vs. \\Theta (all pixels)");
    k1->SetXTitle("\\Theta [\\circ]");
    k1->SetYTitle("Sigma");

    k1->DrawCopy("box");
    k1->SetBit(kCanDelete);

    c.cd(9);
    TH2D *k2 = NULL;
    k2 = (TH2D*) ((TH3*)fHSigmaPixTheta)->Project3D("zy");
    k2->SetDirectory(NULL);
    k2->SetTitle("(Target) Sigma vs. pixel number (all \\Theta)");
    k2->SetXTitle("pixel");
    k2->SetYTitle("Sigma");

    k2->DrawCopy("box");
    k2->SetBit(kCanDelete);;
    */

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

    
    c.cd(4);
    TH2D *m1;
    m1 = (TH2D*) ((TH3*)fHgON)->Project3D("zy");
    m1->SetDirectory(NULL);
    m1->SetTitle("(fHgON) Sigmabar-new vs. Sigmabar-old (ON, all  \\Theta)");
    m1->SetXTitle("Sigmabar-old");
    m1->SetYTitle("Sigmabar-new");

    m1->DrawCopy("box");
    m1->SetBit(kCanDelete);;

    c.cd(5);
    TH2D *m2;
    m2 = (TH2D*) ((TH3*)fHgOFF)->Project3D("zy");
    m2->SetDirectory(NULL);
    m2->SetTitle("(fHgOFF) Sigmabar-new vs. Sigmabar-old (OFF, all  \\Theta)");
    m2->SetXTitle("Sigmabar-old");
    m2->SetYTitle("Sigmabar-new");

    m2->DrawCopy("box");
    m2->SetBit(kCanDelete);;

    c.cd(6);
    TH2D *m3;
    m3 = (TH2D*) ((TH3*)fHgMC)->Project3D("zy");
    m3->SetDirectory(NULL);
    m3->SetTitle("(fHgMC) Sigmabar-new vs. Sigmabar-old (MC, all  \\Theta)");
    m3->SetXTitle("Sigmabar-old");
    m3->SetYTitle("Sigmabar-new");

    m3->DrawCopy("box");
    m3->SetBit(kCanDelete);;


  return kTRUE;
}

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










