/* ======================================================================== *\
!
! *
! * 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): Oscar Blanch    12/2001 <mailto:blanch@ifae.es>
!   Author(s): Thomas Bretz    08/2002 <mailto:tbretz@astro.uni.wuerzburg.de>
!   Author(s): Nadia Tonello   02/2004 <mailto:tonello@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//  MBlindPixelsCalc2
//
//  NT:: WARNING: This class will change soon to fit into mbadpixels!!
//  It does all MBlindPixelCalc did, plus check the pedestal Rms 
//  to define bad pixels.
// 
//  This is the specific image cleaning for a list of pixels. This task
//  sets the pixels listed in fPixelsIdx to unused so they should not be
//  used for analysis (eg calculation of hillas parameters).
//
//  NT::You can use MBlindPixelsCalc2::SetCheckPedestalRms to see which
//  pixels to blind. If the PedestalRms in one pixel is more than 
//  3 times the mean or less than 1/3 the mean, they are set as blind.
//
//  You can use MBlindPixelsCalc2::SetUseInterpolation to replace the
//  blind pixels by the average of its neighbors instead of unmapping
//  them. If you want to include the central pixel use
//  MBlindPixelsCalc2::SetUseCentralPixel.
//  NT:: The interpolation will be done only if there are at least
//  3 non-blind pixels. 
//
//  Example:
//     MBadPixelCalc2 blind;
//     blind.SetCheckPedestalRms();
//     blind.SetUseInterpolation();
//
//  If you do not use SetUseInterpolation(), the pixels selected with 
//  SetCheckPedestalRms will be Unmapped and not used for the image
//  analysis.
//
//  Input Containers:
//   MCerPhotEvt
//   [MBlindPixels]
//
//  Output Containers:
//   MCerPhotEvt
//   MBlindPixels
//
/////////////////////////////////////////////////////////////////////////////
#include "MBlindPixelsCalc2.h"

#include <fstream>

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

#include "MParList.h"

#include "MGeomPix.h"
#include "MGeomCam.h"
#include "MCerPhotPix.h"
#include "MCerPhotEvt.h"
#include "MPedPhotPix.h"
#include "MPedPhotCam.h"
#include "MPedestalPix.h"
#include "MPedestalCam.h"
#include "MBlindPixels.h"

#include "MMcRunHeader.hxx"

ClassImp(MBlindPixelsCalc2);

using namespace std;

static const TString gsDefName  = "MBlindPixelsCalc2";
static const TString gsDefTitle = "Find hot spots (star, broken pixels, etc)";

// --------------------------------------------------------------------------
//
// Default constructor.
//
MBlindPixelsCalc2::MBlindPixelsCalc2(const char *name, const char *title)
  : fFlags(0), fErrors(3)
{
    fName  = name  ? name  : gsDefName.Data();
    fTitle = title ? title : gsDefTitle.Data();
}

// --------------------------------------------------------------------------
//
//  - Try to find or create MBlindPixels in parameter list.
//  - get the MCerPhotEvt from the parlist (abort if missing)
//
Int_t MBlindPixelsCalc2::PreProcess (MParList *pList)
{
    if (TESTBIT(fFlags, kUseBlindPixels))
      fPixels = (MBlindPixels*)pList->FindObject(AddSerialNumber("MBlindPixels"));
    else
        fPixels = (MBlindPixels*)pList->FindCreateObj(AddSerialNumber("MBlindPixels"));

    if (!fPixels)
        {*fLog << err << dbginf << "Try to remove  -SetUseBlindPixels-  from the task list to go on" << endl; 
    return kFALSE;}

    fEvt = (MCerPhotEvt*)pList->FindObject(AddSerialNumber("MCerPhotEvt"));
    if (!fEvt)
    {
        *fLog << err << dbginf << "MCerPhotEvt not found... aborting." << endl;
        return kFALSE;
    }

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

    fGeomCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
    if (!fGeomCam)
        *fLog << warn << dbginf << "No camera geometry available... can't use interpolation." << endl;

    const UShort_t size = fPixelsIdx.GetSize();
    if (size == 0)
    {
        if (!pList->FindObject("MMcRunHeader") && !fPixels)
        {
            *fLog << warn << "Warning - Neither blind pixels are given nor a MMcRunHeader was found... removing MBlindPixelsCalc2 from list." << endl;
            return kSKIP;
        }
        return kTRUE;
    }

    // Set as blind pixels the global blind pixels, which are given
    // through the macros

    UShort_t numids = fPixelsIdx.GetSize();

    for(Int_t i = 0; i<numids; i++)
	fPixels->SetPixelBlind(fPixelsIdx[i]);


   memset(fErrors.GetArray(), 0, sizeof(Char_t)*fErrors.GetSize());


    fErrors[0]= 0;
    fErrors[1]= 0;
    fErrors[2]= 0;

    return kTRUE;
}

// --------------------------------------------------------------------------
//
//  Replaces each pixel (signal, signal error, pedestal, pedestal rms)
//  by the average of its surrounding pixels.
//  If TESTBIT(fFlags, kUseCentralPixel) is set the central pixel is also
//  included.
//
void MBlindPixelsCalc2::Interpolate() const
{
    const UShort_t entries = fGeomCam->GetNumPixels();

    //
    // Create arrays
    //
    Double_t *nphot  = new Double_t[entries];
    Double_t *perr   = new Double_t[entries];
 
    //
    // Loop over all pixels
    //
    for (UShort_t i=0; i<entries; i++)
    {
        //
        // Check whether pixel with idx i is blind
        //
        if (!fPixels->IsBlind(i))
            continue;

        //
        // Get a pointer to this pixel. If it is not yet existing
        // create a new entry for this pixel in MCerPhotEvt
        //

	MCerPhotPix *pix = fEvt->GetPixById(i);
        if (!pix)
            pix = fEvt->AddPixel(i, 0, 0);

        //
        // Get the corresponding geometry and pedestal
        //
        const MGeomPix &gpix = (*fGeomCam)[i];

        // Do Not-Use-Central-Pixel
        const Bool_t nucp = !TESTBIT(fFlags, kUseCentralPixel);

        Int_t num = nucp ? 0 : 1;

        nphot[i]  = nucp ? 0 : pix->GetNumPhotons();
        perr[i]   = nucp ? 0 : Pow2(pix->GetErrorPhot());

        //
	// The values are rescaled to the small pixels area for the right comparison
        //
        const Double_t ratio = fGeomCam->GetPixRatio(i);

        nphot[i]  *= ratio;
        perr[i]   *= ratio;

        //
        // Loop over all its neighbors
        //
        const Int_t n = gpix.GetNumNeighbors();
        for (int j=0; j<n; j++)
        {
            const UShort_t nidx = gpix.GetNeighbor(j);

            //
            // Do not use blind neighbors
            //
            if (fPixels->IsBlind(nidx))
                continue;

            //
            // Check whether the neighbor has a signal stored
            //
            const MCerPhotPix *evtpix = fEvt->GetPixById(nidx);
            if (!evtpix)
                continue;

            //
            // Get the geometry for the neighbor
            //
            const Double_t nratio = fGeomCam->GetPixRatio(nidx);

            //
	    //The error is calculated as the quadratic sum of the errors
	    //
            nphot[i]  += nratio*evtpix->GetNumPhotons();
            perr[i]   += nratio* Pow2(evtpix->GetErrorPhot());

            num++;
        }

        //
	// Now the mean is calculated and the values rescaled back to the pixel area
        //
	nphot[i] /= (num*ratio);
        perr[i]   = TMath::Sqrt(perr[i]/(num*ratio));
 

	// Check if there are enough neighbors to calculate the mean
	// If not, unmap the pixel. The maximum number of blind neighbors should be 2 

	if (num < 3) 
	  { pix->SetPixelUnmapped();
	  //	    pix->Set(nphot[i], perr[i]);
 	    continue;
	  }

	    pix->Set(nphot[i], perr[i]);
    }
    //
    // Delete the intermediat arrays
    //
    delete nphot;
    delete perr;
}

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

void MBlindPixelsCalc2::InterpolatePedestals() const
{
  const Int_t entries = fPed->GetSize();
  Double_t *ped    = new Double_t[entries];
  Double_t *pedrms = new Double_t[entries];
  
  //
  // Loop over all pixels
  //
  for (UShort_t i=0; i<entries; i++)
    {
      //
      // Check whether pixel with idx i is blind
      //
      if (!fPixels->IsBlind(i))
	continue;
      
      //
      // Get the corresponding geometry and pedestal
      //
      const MGeomPix &gpix = (*fGeomCam)[i];
      const MPedPhotPix &ppix = (*fPed)[i];

      // Do Not-Use-Central-Pixel
      const Bool_t nucp = !TESTBIT(fFlags, kUseCentralPixel);
      
      Int_t num = nucp ? 0 : 1;
      
      ped[i]   = nucp ? 0 : ppix.GetMean();
      pedrms[i] = nucp ? 0 : Pow2(ppix.GetRms());
      
      //
      // The values are rescaled to the small pixels area for the right comparison
      //
      const Double_t ratio = fGeomCam->GetPixRatio(i);
      
      ped[i]    *= ratio;
      pedrms[i] *= ratio;
      //
      // Loop over all its neighbors
      //
      const Int_t n = gpix.GetNumNeighbors();
      for (int j=0; j<n; j++)
	{
	  const UShort_t nidx = gpix.GetNeighbor(j);
	  
	  //
	  // Do not use blind neighbors
	  //
	  if (fPixels->IsBlind(nidx))
	    continue;
	  
	  //
	  // Get the geometry for the neighbor
	  //
	  const Double_t nratio = fGeomCam->GetPixRatio(nidx);
	  MPedPhotPix &nppix = (*fPed)[nidx];
	  
	  //
	  //The error is calculated as the quadratic sum of the errors
	  //
	  ped[i]    += (nratio*nppix.GetMean());
	  pedrms[i] += nratio*Pow2(nppix.GetRms());
	  
	  num++;
	}
      
      // Check if there are enough neighbors to calculate the mean
      // If not, unmap the pixel. The minimum number of good neighbors 
      // should be 3 
      
      if (num < 3) 
	{ 	
	  MCerPhotPix *pix = fEvt->GetPixById(i);
	  if (!pix)
	    pix = fEvt->AddPixel(i, 0, 0);
	  pix->SetPixelUnmapped();
	  continue;
	}

      //
      // Now the mean is calculated and the values rescaled back to the pixel area
      //
      ped[i]   /=  (num*ratio);
      pedrms[i] = TMath::Sqrt(pedrms[i]/(num*ratio));
      
     (*fPed)[i].Set(ped[i], pedrms[i]);
    }
  
  //
  // Delete the arrays
  //
  delete ped;
  delete pedrms;
  
  return;
}

// --------------------------------------------------------------------------
//
//  Removes all blind pixels from the analysis by setting their state
//  to unmapped.
//
void MBlindPixelsCalc2::Unmap() const
{
    const UShort_t entries = fEvt->GetNumPixels();

    //
    // remove the pixels in fPixelsIdx if they are set to be used,
    // (set them to 'unused' state)
    //
    for (UShort_t i=0; i<entries; i++)
    {
        MCerPhotPix &pix = (*fEvt)[i];

        if (fPixels->IsBlind(pix.GetPixId()))
	      pix.SetPixelUnmapped();
    }
}

// --------------------------------------------------------------------------
// Check the pedestal Rms of the pixels: compute with 2 iterations the mean 
// for inner and outer pixels. Set as blind the pixels with too small or 
// too high pedestal Rms with respect to the mean.
// 
Bool_t MBlindPixelsCalc2::CheckPedestalRms() 

{
  const Int_t entries = fPed->GetSize();
  Int_t pixtoblind = 0;
  Int_t pixin = 0;
  Int_t pixout = 0;

  Float_t meanPedRms =0.;
  Float_t meanPedRmsout = 0.;
  Float_t meanPedRms2 = 0.;
  Float_t meanPedRmsout2 = 0.;

 for (Int_t i=0; i<entries; i++)
    {
      //
      // get pixel as entry from list
      // 
      const MPedPhotPix &ppix = (*fPed)[i];
      const Double_t nratio = fGeomCam->GetPixRatioSqrt(i);
      const Double_t pixPedRms = ppix.GetRms();
     
      //Calculate the corrected means:
      
      if (nratio == 1 )
	{
	  if(  pixPedRms >0. &&  pixPedRms <200. )
	  {meanPedRms += pixPedRms;
	  pixin++;}
	}
      else 
	{
	  if(  pixPedRms >0. &&  pixPedRms <500. )
	  {meanPedRmsout += pixPedRms;
	  pixout++;}
	}
    }
 
      //if no pixel has a minimum signal, return
 if( pixin==0 || pixout==0 )
   {fErrors[1]++;          //no valid Pedestals Rms
   return kFALSE;
   } 
 
 meanPedRms =  meanPedRms / pixin ;
 meanPedRmsout =  meanPedRmsout / pixout ;

 if (  meanPedRms == 0. ||   meanPedRmsout ==0.)
   {fErrors[1]++;         //no valid Pedestals Rms
   return kFALSE;
   }
 
 pixin  = 0;
 pixout = 0;

 for (Int_t i=0; i<entries; i++)
   {
     //
     // get pixel as entry from list
     // 
     const MPedPhotPix &ppix = (*fPed)[i];
     const Double_t nratio = fGeomCam->GetPixRatioSqrt(i);
     const Double_t pixPedRms = ppix.GetRms();
     
     //Calculate the corrected means:
     
     if (nratio == 1)
       {if (pixPedRms > 0.5 * meanPedRms  && pixPedRms < 1.5 * meanPedRms)  
	 {meanPedRms2 += pixPedRms;
	 pixin++;}
       }
     else 
       {if (pixPedRms > 0.5 * meanPedRmsout  && pixPedRms < 1.5 * meanPedRmsout)  
	 {meanPedRmsout2 += pixPedRms;
	 pixout++;}
       }
   }
     
      //if no pixel has a minimum signal, return
 if( pixin==0 || pixout==0 )
   {fErrors[1]++;            //no valid Pedestals Rms
   return kFALSE;} 
 
 meanPedRms =  meanPedRms2 / pixin ;
 meanPedRmsout =  meanPedRmsout2 / pixout ;

 //Blind the Bad Pixels
 for (Int_t i=0; i<entries; i++)
   {
     const MPedPhotPix &ppix = (*fPed)[i];
     const Float_t nratio = fGeomCam->GetPixRatioSqrt(i);
     const Double_t pixPedRms = ppix.GetRms();
     
     if (nratio == 1)	
       { if (pixPedRms > 3.* meanPedRms || pixPedRms <= meanPedRms/3.)
	 { fPixels->SetPixelBlind(i); 
	 pixtoblind++;}
       }
     else
       { if (pixPedRms > 3.* meanPedRmsout || pixPedRms <= meanPedRmsout/3.)
	 { fPixels->SetPixelBlind(i); 
	 pixtoblind++; }
       }
   }
 
 // Check if the number of pixels to blind is > 60% of total number of pixels
 // 
 if (pixtoblind>0.6*entries)
   { fErrors[2]++;
   return kFALSE; }
 
 
 return kTRUE;
 
}

// --------------------------------------------------------------------------
//
// Treat the blind pixels
//
Int_t MBlindPixelsCalc2::Process()
{
   if (TESTBIT(fFlags, kCheckPedestalRms))
     {
       // if the number of blind pixels is too high, do not interpolate
       if (CheckPedestalRms()== kFALSE) 
	   return kTRUE;
       else
	 {if (TESTBIT(fFlags, kUseInterpolation))
	   InterpolatePedestals();	         
	 } 
     }

   if (TESTBIT(fFlags, kUseInterpolation) && fGeomCam)
     {   Interpolate();
     } 
   else
     { Unmap();
     } 
     fErrors[0]++;
   return kTRUE;
}


// --------------------------------------------------------------------------
//
//  - Check whether pixels to disable are available. If pixels are
//    given by the user nothing more is done.
//
Bool_t MBlindPixelsCalc2::ReInit(MParList *pList)
{
    if (TESTBIT(fFlags, kUseBlindPixels))
        return kTRUE;

    //
    // If pixels are given by the user, we are already done
    //
    if (fPixelsIdx.GetSize() > 0)
        return kTRUE;

    //
    // Delete the old array holding the blind pixels for the last file
    //
    fPixels->Clear();

    if (!fGeomCam->InheritsFrom("MGeomCamMagic"))
    {
        *fLog << warn << "MBlindPixelsCalc2::ReInit: Warning - Starfield only implemented for Magic standard Camera... no action." << endl;
        return kTRUE;
    }
   return kTRUE;
}

void MBlindPixelsCalc2::StreamPrimitive(ofstream &out) const
{
    out << "   MBlindPixelsCalc2 " << GetUniqueName();
    if (fName!=gsDefName || fTitle!=gsDefTitle)
    {
        out << "(\"" << fName << "\"";
        if (fTitle!=gsDefTitle)
            out << ", \"" << fTitle << "\"";
        out <<")";
    }
    out << ";" << endl;

    if (TESTBIT(fFlags, kUseInterpolation))
        out << "   " << GetUniqueName() << ".SetUseInterpolation();" << endl;
    if (TESTBIT(fFlags, kUseCentralPixel))
        out << "   " << GetUniqueName() << ".SetUseCentralPixel();" << endl;
    if (TESTBIT(fFlags, kCheckPedestalRms))
        out << "   " << GetUniqueName() << ".SetCheckPedestalRms();" << endl;

    if (fPixelsIdx.GetSize()==0)
        return;

    out << "   {" << endl;
    out << "      TArrayS dummy;" << endl;
    for (int i=0; i<fPixelsIdx.GetSize(); i++)
        out << "      dummy[" << i << "]=" << ((TArrayS)fPixelsIdx)[i] << ";" << endl;
    out << "      " << GetUniqueName() << ".SetPixels(dummy);" << endl;
    out << "   }" << endl;
}


// --------------------------------------------------------------------------
//
//  This is used to print the output in the PostProcess. Later (having
//  millions of events) we may want to improve the output.
//
void MBlindPixelsCalc2::PrintSkipped(int i, const char *str) const
{
    *fLog << " " << setw(7) << fErrors[i] << " (";
    *fLog << setw(3) << (int)(100.*fErrors[i]/GetNumExecutions());
    *fLog << "%) Evts skipped due to: " << str << endl;
}

// --------------------------------------------------------------------------
//
//  Prints some statistics about the calculation. The percentage
//  is calculated with respect to the number of executions of this task.
//
Int_t MBlindPixelsCalc2::PostProcess()
{
    if (GetNumExecutions()==0)
        return kTRUE;

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

    // Number of events where the computation of the mean pedestal Rms failed 
    // On these events no blind pixels have been selected

    if (TESTBIT(fFlags, kCheckPedestalRms))
    PrintSkipped(1, "Big fluctuations in Pedestal Rms, no CHECK");

    // Number of events with too many bad pixels have been found
    // On these events, no interpolation has been applied

    if (TESTBIT(fFlags, kUseInterpolation))
    PrintSkipped(2, " >60% of pixels are bad, no INTERPOLATION");

    //Number of events that have been  succesfully checked and corrected

    *fLog << " " << (int)fErrors[0] << " (" << (int)(100.*fErrors[0]/GetNumExecutions()) << "%) Evts well treated!" << endl;
    *fLog << endl;

    return kTRUE;
}