/* ======================================================================== *\
!
! *
! * 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  02/2003 <wittek@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2002
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
//  MSelBasic                                                              //
//                                                                         //
//  This is a task to evaluate basic cuts                                  //
//                                                                         //
//  to be called after the calibration (when the number of photons is      //
//               available for all pixels)                                  //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////

#include "MSelBasic.h"

#include "MParList.h"

#include "MHillas.h"
#include "MCerPhotEvt.h"
#include "MMcEvt.hxx"
#include "MRawRunHeader.h"
#include "MGeomCam.h"
#include "MPedestalCam.h"
#include "MPedestalPix.h"
#include "MGeomPix.h"

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

ClassImp(MSelBasic);

// --------------------------------------------------------------------------
//
// Default constructor.
//
MSelBasic::MSelBasic(const char *name, const char *title)
{
    fName  = name  ? name  : "MSelBasic";
    fTitle = title ? title : "Task to evaluate basic cuts";

    ThetaMin =  0.0;
    ThetaMax = 60.0;
}

// --------------------------------------------------------------------------
//
// 
// 
// 
//
Bool_t MSelBasic::PreProcess(MParList *pList)
{
    fRawRun = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
    if (!fRawRun)
    {
        *fLog << dbginf << "MRawRunHeader not found... aborting." << endl;
        return kFALSE;
    }

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

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


    fCam = (MGeomCam*)pList->FindObject("MGeomCam");
    if (!fCam)
    {
        *fLog << dbginf << "MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl;
        return kFALSE;
    }

    fPed = (MPedestalCam*)pList->FindObject("MPedestalCam");
    if (!fPed)
    {
        *fLog << dbginf << "MPedestalCam missing in Parameter List... aborting." << endl;
        return kFALSE;
    }

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

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Evaluate basic cuts
// 
// if cuts are fulfilled      : return 0
// if they are not fullfilled : skip remaining tasks for this event
//
Bool_t MSelBasic::Process()
{
  /*
  //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  *fLog << "=========================================================" << endl;
  *fLog << "" << endl;
  *fLog << "MmcEvt data : " << endl;
  *fLog << "" << endl;
  fMcEvt->Print();
  *fLog << "" << endl;
  *fLog << "PartId() = " << fMcEvt->GetPartId() << endl;
  *fLog << "Energy() = " << fMcEvt->GetEnergy() << endl;
  *fLog << "Theta()  = " << fMcEvt->GetTheta() << endl;

  *fLog << "Phi()    = " << fMcEvt->GetPhi() << endl;
  //*fLog << "CoreD()  = " << fMcEvt->GetCoreD() << endl;
  //*fLog << "CoreX()  = " << fMcEvt->GetCoreX() << endl;

  //*fLog << "CoreY()  = " << fMcEvt->GetCoreY() << endl;
  *fLog << "Impact() = " << fMcEvt->GetImpact() << endl;
  //*fLog << "PhotIni()= " << fMcEvt->GetPhotIni() << endl;

  //*fLog << "PassPhotAtm() = " << fMcEvt->GetPassPhotAtm() << endl;
  //*fLog << "PassPhotRef() = " << fMcEvt->GetPassPhotRef() << endl;
  //*fLog << "PassPhotCone()  = " << fMcEvt->GetPassPhotCone() << endl;

  //*fLog << "PhotElfromShower()    = " << fMcEvt->GetPhotElfromShower() << endl;
  //*fLog << "PhotElinCamera()  = " << fMcEvt->GetPhotElinCamera() << endl;
  *fLog << "TelescopePhi()  = " << fMcEvt->GetTelescopePhi() << endl;

  *fLog << "TelescopeTheta()  = " << fMcEvt->GetTelescopeTheta() << endl;
  *fLog << "Impact() = " << fMcEvt->GetImpact() << endl;
  //*fLog << "PhotIni()= " << fMcEvt->GetPhotIni() << endl;
  *fLog << "" << endl;
  *fLog << "=========================================================" << endl;

  *fLog << "=========================================================" << endl;
  *fLog << "" << endl;
  *fLog << "MPedestalPix data : " << endl;
  *fLog << "" << endl;

  Int_t ntot;
  ntot = fPed->GetSize();
    *fLog << "MeanRms() :" << endl;    
  for (Int_t i=0; i<ntot; i++)
  {
    MPedestalPix &pix = (*fPed)[i];
    //*fLog << "Mean()    = " << i << ",  " << pix.GetMean() << endl;    
    //*fLog << "Sigma()   = " << i << ",  " << pix.GetSigma() << endl;    
    *fLog << pix.GetMeanRms() << " ";    
    //*fLog << "SigmaRms()= " << i << ",  " << pix.GetSigmaRms() << endl;    
  }
  *fLog << "" << endl;

  *fLog << "" << endl;
  ntot = fEvt->GetNumPixels();
    *fLog << "MCerPhotPix :  pix.GetNumPhotons()" << endl;    
  for (Int_t i=0; i<ntot; i++)
  {
    MCerPhotPix &pix = (*fEvt)[i];
    *fLog << pix.GetNumPhotons() << " ";    
  }
  *fLog << "" << endl;

  *fLog << "" << endl;
  ntot = fEvt->GetNumPixels();
    *fLog << "MCerPhotPix :  pix.GetErrorPhot()" << endl;    
  for (Int_t i=0; i<ntot; i++)
  {
    MCerPhotPix &pix = (*fEvt)[i];
    *fLog << pix.GetErrorPhot() << " ";    
  }
  *fLog << "" << endl;

  //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
  */


    Int_t rc = 0;

    //if ( fRawRun->GetRunNumber() < 16279 ) 
    //{
    //   rc = 1;
    //   return kCONTINUE;
    //}

    Double_t Theta = kRad2Deg*fMcEvt->GetTelescopeTheta();
    if ( Theta < ThetaMin )
    {
      *fLog << "MSelBasic::Process; Run, Event, Theta = " 
            << fRawRun->GetRunNumber()<< ",  " 
            << fMcEvt->GetEvtNumber() << ",  " << Theta << endl;
      rc = 1;
    }    

    else if ( Theta > ThetaMax )
    {
      rc = 2;
    }    

    else if ( !SwTrigger() )
    {
      //*fLog << "MSelBasic::Process; SwTrigger = " << SwTrigger << endl;
      rc = 3;
    }    

    fErrors[rc]++;

    return rc==0 ? kTRUE : kCONTINUE;
}
// --------------------------------------------------------------------------
//
// Software trigger
// 
// require 2 neighboring pixels (which are not in the outermost ring), 
//                       each having at least 'minphotons' photons
// 
// 
Bool_t MSelBasic::SwTrigger()
{
    // minimum number of photons required
    Double_t minphotons = 13.0;

    const Int_t entries = fEvt->GetNumPixels();
 
    //$$$$$$$$$$$$$$$$$$
    //const Int_t nall = fPed->GetSize();  
    //*fLog << "nall = " << nall << endl;
    //for (Int_t id=0; id<nall; id++)
    //{
    //  MGeomPix &gpix = (*fCam)[id];
    //  if ( gpix.IsInOutermostRing() ) 
    //  {
    //    *fLog << "IsInOutermostRing : pixel no. = " << id << endl;
    //  }

    //  if ( gpix.IsInOuterRing() ) 
    //  {
    //    *fLog << "IsInOuterRing : pixel no. = " << id << endl;
    //  }
    //}
    //$$$$$$$$$$$$$$$$$$


    for (Int_t i=0; i<entries; i++)
    {
      MCerPhotPix &pix = (*fEvt)[i];
      Int_t id = pix.GetPixId();
      if (!pix.IsPixelUsed()) continue;

      Double_t photons = pix.GetNumPhotons();
      if (photons < minphotons) continue;

      // this pixel is used and has the required no.of photons
      // check whether this is also true for a neigboring pixel

      MGeomPix &gpix = (*fCam)[id];
      if ( gpix.IsInOutermostRing() ) continue;

      const Int_t nneighbors = gpix.GetNumNeighbors();
      for (Int_t n=0; n<nneighbors; n++)
      {
        const Int_t id1 =  gpix.GetNeighbor(n);
        if ( !fEvt->IsPixelUsed(id1) ) continue;

        MGeomPix &gpix1 = (*fCam)[id1];
        if ( gpix1.IsInOutermostRing() ) continue;

        //MCerPhotPix &pix1 = (*fEvt)[id1];      
        MCerPhotPix &pix1 = *(fEvt->GetPixById(id1));      
        if ( &pix1 == NULL ) 
	{
          *fLog << "MSelBasic::SwTrigger; &pix1 is NULL" << endl;
          continue;
	}

        Double_t photons1 = pix1.GetNumPhotons();
        if (photons1 >= minphotons) return kTRUE;
      }
    }
    return kFALSE;
}

// --------------------------------------------------------------------------
//
//  Prints some statistics about the Basic selections.
//
Bool_t MSelBasic::PostProcess()
{
    if (GetNumExecutions()==0)
        return kTRUE;

    *fLog << inf << endl;
    *fLog << GetDescriptor() << " execution statistics:" << endl;
    *fLog << dec << setfill(' ');
    *fLog << " " << setw(7) << fErrors[1] << " (" << setw(3) 
          << (int)(fErrors[1]*100/GetNumExecutions()) 
          << "%) Evts skipped due to: Zenith angle < " << ThetaMin << endl;

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

    *fLog << " " << setw(7) << fErrors[3] << " (" << setw(3) 
          << (int)(fErrors[3]*100/GetNumExecutions()) 
          << "%) Evts skipped due to: Software trigger not fullfilled" << endl;

    *fLog << " " << fErrors[0] << " (" << (int)(fErrors[0]*100/GetNumExecutions()) 
          << "%) Evts survived Basic selections!" << endl;
    *fLog << endl;

    return kTRUE;
}