source: trunk/MagicSoft/Mars/mtemp/mmpi/MApplySupercuts.cc@ 6985

/* 
======================================================================== *\
!
! *
! * 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): Nicola Galante  02/2005 <mailto:nicola.galante@pi.infn.it>
!
!   Copyright: MAGIC Software Development, 2004
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//   MApplySupercuts
//
// Apply supercuts to a file containing Hillas parameters
// (usually a star file).
//
// Input:
//   rootfile (starfile) of Hillas parameters
//
// Output:
//   MHadronness estimated on event
//
/////////////////////////////////////////////////////////////////////////////
#include <TMath.h>
#include <TFile.h>
#include <TH1F.h>
#include <TCanvas.h>

#include "MGeomCamMagic.h"
#include "MLog.h"
#include "MLogManip.h"
#include "MParList.h"
#include "MRawEvtHeader.h"
#include "MRawRunHeader.h"
#include "MSupercuts.h"
#include "MHadronness.h"
#include "MHillas.h"
#include "MHillasSrc.h"
#include "MHFindSignificance.h"
#include "MEnergyEst.h"
#include "MApplySupercuts.h"

ClassImp(MApplySupercuts);

using namespace std;

// --------------------------------------------------------------------------
//
// Default constructor
//
MApplySupercuts::MApplySupercuts(const char *name, const char *title)
  : fRunHeader(0), fEvtHeader(0), fSupercuts(0)
{
    fName  = name  ? name  : "MApplySupercuts";
    fTitle = title ? title : "Task to apply supercuts to star files";

    fPhe = kTRUE;
    fPlot = kFALSE;
}


// --------------------------------------------------------------------------
//
Int_t MApplySupercuts::PreProcess(MParList *pList)
{
  fGeom = (MGeomCamMagic *)pList->FindCreateObj("MGeomCamMagic");
  if (!fGeom)
    {
      *fLog << err << "MGeomCamMagic not found... abort." << endl;
      return kFALSE;
    }
  
  fRunHeader = (MRawRunHeader*)pList->FindCreateObj("MRawRunHeader");
  if (!fRunHeader)
    {
      *fLog << err << "MRawRunHeader not found... abort." << endl;
      return kFALSE;
    }
  
  fEvtHeader = (MRawEvtHeader *)pList->FindObject("MRawEvtHeader");
  if (!fEvtHeader)
    {
      *fLog << err << "MRawEvtHeader not found... but ok." << endl;
    }
  
  fHillas = (MHillas *)pList->FindObject("MHillas");
  if (!fHillas)
    {
      *fLog << err << "MHillas not found... abort." << endl;
      return kFALSE;
    }
  
  fHillasSrc = (MHillasSrc *)pList->FindObject("MHillasSrc");
  if (!fHillasSrc)
    {
      *fLog << err << "MHillasSrc not found... abort." << endl;
      return kFALSE;
    }
  
  fHadronness = (MHadronness *)pList->FindCreateObj("MHadronness",AddSerialNumber("MHadronness"));
  if (!fHadronness)
    {
      *fLog << err << "MHadronness not found... abort." << endl;
      return kFALSE;
    }

  fEnergyEst = (MEnergyEst *)pList->FindCreateObj("MEnergyEst");
  if (!fEnergyEst)
    {
      *fLog << err << "MEnergyEst not found... abort." << endl;
      return kFALSE;
    }

  if(fPlot){
    fAlpha = new TH1F("fHistAlpha","Alpha Plot",30,0,90);
      if (!fAlpha)
        {
          *fLog << err << "Impossible to create Alpha histogram... abort." << endl;
          return kFALSE;
        }
  }

  fHFindSigma = (MHFindSignificance *)pList->FindCreateObj("MHFindSignificance");
  if (!fHFindSigma)
    {
      *fLog << err << "MHFindSignificance not found... abort." << endl;
      return kFALSE;
    }

  // Start reading supercuts
  fSupercuts = (MSupercuts *)pList->FindCreateObj("MSupercuts",AddSerialNumber("MSupercuts"));
  if (!fSupercuts)
    {
      *fLog << err << "MSupercuts not found... abort." << endl;
      return kFALSE;
    }
  // read the cuts coefficients
  TFile inparam(fSCFilename);
  fSupercuts->Read("MSupercuts");
  inparam.Close();
  *fLog << "MFindSupercutsONOFF::FindParams; initial values of parameters are taken from file "
        << fSCFilename << endl;
  
  TArrayD params = fSupercuts->GetParameters();
  
  for (Int_t i=0; i<8; i++)
    {
      fLengthUp[i] = params[i];
      fLengthLo[i] = params[i+8];
      fWidthUp[i]  = params[i+16];
      fWidthLo[i]  = params[i+24];
      fDistUp[i]   = params[i+32];
      fDistLo[i]   = params[i+40];
    }
    
  // end reading supercuts
    return kTRUE;
}

// --------------------------------------------------------------------------
//
Int_t MApplySupercuts::Process()
{
  fHadronness->SetHadronness(kHadronHadronness);

  Float_t fMm2Deg = fGeom->GetConvMm2Deg();
  Double_t log3000 = TMath::Log(fSizeOffset);
  Float_t fsize = fHillas->GetSize();
  if(fPhe)
    fsize /= kPhe2Ph;
  Float_t fdist = fHillasSrc->GetDist()*fMm2Deg;
  Double_t logsize = TMath::Log((Double_t)fsize);
  Double_t lgsize = logsize-log3000;
  Double_t lgsize2 = lgsize*lgsize;
  Double_t dist2 =  (Double_t)fdist*fdist - fDistOffset*fDistOffset;
  //Double_t dist2 = (fDist-fDistOffset)*(fDist-fDistOffset);

  // parameters:
  Float_t flength = (fHillas->GetLength()) * fMm2Deg;
  Float_t fwidth = (fHillas->GetWidth())*fMm2Deg;
  //Float_t fsize = fHillas.GetSize()/0.18;
  //Float_t fmeanx = (fHillas.GetMeanX())*fMm2Deg;
  //Float_t fmeany = (fHillas.GetMeanY())*fMm2Deg;
  //Float_t falpha = fHillasSrc.GetAlpha();
  //Float_t fDCAdelta = fHillasSrc.GetDCADelta();



  if ( flength < CalcLimit(fLengthUp, lgsize, lgsize2, dist2) &&
       flength > CalcLimit(fLengthLo, lgsize, lgsize2, dist2) &&
       fwidth  < CalcLimit(fWidthUp, lgsize, lgsize2, dist2)  &&
       fwidth  > CalcLimit(fWidthLo, lgsize, lgsize2, dist2)  &&
       fdist   < CalcLimit(fDistUp, lgsize, lgsize2, dist2)   &&
       fdist   > CalcLimit(fDistLo, lgsize, lgsize2, dist2) )
    {
      // gamma candidates!
      fHadronness->SetHadronness(kGammaHadronness);
      if( fPlot && fHillas->GetSize()>fSizeLow && fHillas->GetSize()<fSizeUp )
        fAlpha->Fill(TMath::Floor(TMath::Abs(fHillasSrc->GetAlpha())));
    }

  return kTRUE;
}


Int_t MApplySupercuts::PostProcess()
{
  if(fPlot){
    TFile f("shit_file.root","RECREATE");
    fAlpha->Write();
    f.Close();
    TCanvas c;
    if(fHFindSigma!=NULL){
      cout << "W el cogno " << fAlphaSignal <<endl;
      fHFindSigma->FindSigma((TH1 *)fAlpha,fAlphaMin,fAlphaMax,fDegree,fAlphaSignal,kTRUE,kTRUE,kTRUE);
      *fLog << "SIGNIFICANCE = " << fHFindSigma->GetSignificance() << endl;
      c.SaveAs("Alpha.gif");
    }
    else
      *fLog << err << "ERROR: fHFindSigma already deleted" << endl;
  }
  return kTRUE;
}


Double_t MApplySupercuts::CalcLimit(Double_t *a, Double_t ls, Double_t ls2, Double_t dd2)
{
  
  Double_t  limit = a[0] + a[1] * dd2 +
    ls  * (a[3] + a[4] * dd2) +
    ls2 * (a[6] + a[7] * dd2);
  
  return limit;
}


void MApplySupercuts::SetPrintOutSignificance(Double_t bgalphamin=50.0, Double_t bgalphamax=90.0, 
                                              Double_t alphasig=15.0, Int_t degree=4)
{
  fAlphaMin = bgalphamin;
  fAlphaMax = bgalphamax;
  fDegree = degree;
  fAlphaSignal = alphasig;

  fPlot = kTRUE;
}