source: trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.cc@ 1917

/* ======================================================================== *\
!
! *
! * 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  04/2003 <mailto:wittek@mppmu.mpg.de>

!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
//   MCT1SupercutsCalc                                                     //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////
#include "MCT1SupercutsCalc.h"

#include <math.h>
#include <fstream.h>

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

#include "MParList.h"
#include "MHillasExt.h"
#include "MHillasSrc.h"
#include "MMcEvt.hxx"
#include "MCerPhotEvt.h"
#include "MGeomCam.h"
#include "MHadronness.h"

ClassImp(MCT1SupercutsCalc);


// --------------------------------------------------------------------------
//
MCT1SupercutsCalc::MCT1SupercutsCalc(const char *hilname, 
               const char *hilsrcname, const char *name, const char *title) 
{
    fName  = name  ? name  : "MCT1SupercutsCalc";
    fTitle = title ? title : "Class to evaluate the Supercuts";

    fHilName    = hilname;
    fHilSrcName = hilsrcname;
}

// --------------------------------------------------------------------------
//
MCT1SupercutsCalc::~MCT1SupercutsCalc() 
{
}


// --------------------------------------------------------------------------
//
Bool_t MCT1SupercutsCalc::PreProcess(MParList *pList)
{
    fHadronness = (MHadronness*)pList->FindCreateObj("MHadronness", fHadronnessName);
    if (!fHadronness)
    {
      *fLog << dbginf << "MHadronness object " << fHadronnessName << " not found... aborting." << endl;
      return kFALSE;
    }


    fHil    = (MHillasExt*)pList->FindObject(fHilName, "MHillasExt");
    if (!fHil)
    {
      *fLog << dbginf << "MHillasExt object " << fHilName << " not found... aborting." << endl;
      return kFALSE;
    }

    fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc");
    if (!fHilSrc)
    {
      *fLog << dbginf << "MHillasSrc object " << fHilSrcName << " not found... aborting." << endl;
      return kFALSE;
    }

    fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
    if (!fMcEvt)
    {
        *fLog << dbginf << "MMcEvt 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;
    }

    fMm2Deg = fCam->GetConvMm2Deg();

    //*fLog << "fMm2Deg = " << fMm2Deg << endl;


    //---------------------------------
    // cut parameters
       
    lengthup[0] = 0.315585; 
    lengthup[1] = 0.001455; 
    lengthup[2] = 0.203198; 
    lengthup[3] = 0.005532; 
    lengthup[4] =-0.001670;
    lengthup[5] =-0.020362; 
    lengthup[6] = 0.007388; 
    lengthup[7] =-0.013463;

    widthup[0] = 0.145412; 
    widthup[1] =-0.001771; 
    widthup[2] = 0.054462; 
    widthup[3] = 0.022280; 
    widthup[4] =-0.009893;
    widthup[5] = 0.056353; 
    widthup[6] = 0.020711; 
    widthup[7] =-0.016703;

    distup[0] = 1.787943; 
    distup[1] = 0.; 
    distup[2] = 2.942310; 
    distup[3] = 0.199815; 
    distup[4] = 0.; 
    distup[5] = 0.249909;
    distup[6] = 0.189697; 
    distup[7] = 0.;

    lengthlow[0] = 0.151530; 
    lengthlow[1] = 0.028323; 
    lengthlow[2] = 0.510707; 
    lengthlow[3] = 0.053089; 
    lengthlow[4] = 0.013708;
    lengthlow[5] = 2.357993; 
    lengthlow[6] = 0.000080; 
    lengthlow[7] =-0.007157;

    widthlow[0] = 0.089187; 
    widthlow[1] =-0.006430; 
    widthlow[2] = 0.074442; 
    widthlow[3] = 0.003738;
    widthlow[4] =-0.004256; 
    widthlow[5] =-0.014101; 
    widthlow[6] = 0.006126; 
    widthlow[7] =-0.002849;

    distlow[0] = 0.589406; 
    distlow[1] = 0.; 
    distlow[2] =-0.083964; 
    distlow[3] =-0.007975; 
    distlow[4] = 0.; 
    distlow[5] = 0.045374;
    distlow[6] =-0.001750; 
    distlow[7] = 0.;

    asymup[0] = 0.061267; 
    asymup[1] = 0.014462; 
    asymup[2] = 0.014327; 
    asymup[3] = 0.014540; 
    asymup[4] = 0.013391;
    asymup[5] = 0.012319; 
    asymup[6] = 0.010444; 
    asymup[7] = 0.008328;

    asymlow[0] =-0.012055; 
    asymlow[1] = 0.009157; 
    asymlow[2] = 0.005441; 
    asymlow[3] = 0.000399; 
    asymlow[4] = 0.001433;
    asymlow[5] =-0.002050; 
    asymlow[6] =-0.000104; 
    asymlow[7] =-0.001188;

    alphaup[0] = 13.123440; 
    alphaup[1] = 0.; 
    alphaup[2] = 0.; 
    alphaup[3] = 0.; 
    alphaup[4] = 0.; 
    alphaup[5] = 0.; 
    alphaup[6] = 0.; 
    alphaup[7] = 0.;
    //---------------------------------

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Calculation of upper and lower limits
//
Double_t MCT1SupercutsCalc::CtsMCut(Double_t *a,  Double_t ls, Double_t ct,
                                 Double_t ls2, Double_t dd2)
{
  // define cut-function
  //
  //    dNOMLOGSIZE = 4.1 (=log(60.0)
  //    dNOMCOSZA   = 1.0
  //
  //      a: array of cut parameters
  //     ls: log(SIZE) - dNOMLOGSIZE
  //    ls2: ls^2
  //     ct: Cos(ZA.) - dNOMCOSZA
  //    dd2: DIST^2

  Double_t limit = *a + *(a+1) * dd2  + *(a+2) * ct 
                      + ls  * (*(a+3) + *(a+4) * dd2 + *(a+5) * ct)
                      + ls2 * (*(a+6) + *(a+7) * dd2);

  //*fLog << "MCT1SupercutsCalc::CtsMCut; *a = "
  //      << *a     << ",  " << *(a+1) << ",  " << *(a+2) << ",  " 
  //      << *(a+3) << ",  " << *(a+4) << ",  " << *(a+5) << ",  " 
  //      << *(a+6) << ",  " << *(a+7) << endl;

  //*fLog << "MCT1SupercutsCalc::CtsMCut; ls, ls2, ct, dd2, limit = " << ls 
  //      << ",  " << ls2 << ",  " << ct << ",  " << dd2 << ",  "
  //      << limit << endl; 
    
  return limit;
}

// ---------------------------------------------------------------------------
//
// Evaluate supercuts for CT1 Mkn421 2001 (Daniel Kranich)
//
Bool_t MCT1SupercutsCalc::Process()
{
  // apply dynamical Scuts 
  // optimized for mkn 421 2001 data

  Double_t dNOMLOGSIZE = 4.1;
  Double_t dNOMCOSZA   = 1.0;


  Double_t newdist = fHilSrc->GetDist()          * fMm2Deg;

  Double_t xbar    = fHil->GetMeanX();
  Double_t ybar    = fHil->GetMeanY();
  Double_t dist    = sqrt(xbar*xbar + ybar*ybar) * fMm2Deg;
  Double_t dd2     = dist * dist;

  Double_t dmls    = log(fHil->GetSize())           - dNOMLOGSIZE;
  Double_t dmls2   = dmls * dmls;

  Double_t dmcza = cos(fMcEvt->GetTelescopeTheta()) - dNOMCOSZA;

  Double_t length  = fHil->GetLength()           *fMm2Deg;
  Double_t width   = fHil->GetWidth()            *fMm2Deg;

  //*fLog << "MCT1SupercutsCalc::Process; dmls, dmcza, dmls2, dd2 = "
  //      << dmls << ",  " << dmcza << ",  " << dmls2 << ",  "
  //      << dd2 << endl;

  //*fLog << "MCT1SupercutsCalc::Process; newdist, dist, length, width = "
  //      << newdist << ",  " << dist << ",  " << length << ",  "
  //      << width << endl;

  if (    newdist  < 1.05
       && newdist  < CtsMCut (distup,    dmls, dmcza, dmls2, dd2)
       && newdist  > CtsMCut (distlow,   dmls, dmcza, dmls2, dd2)
       && dist     < 1.05
       && length   < CtsMCut (lengthup,  dmls, dmcza, dmls2, dd2)
       && length   > CtsMCut (lengthlow, dmls, dmcza, dmls2, dd2)
       && width    < CtsMCut (widthup,   dmls, dmcza, dmls2, dd2)
       && width    > CtsMCut (widthlow,  dmls, dmcza, dmls2, dd2)
       //&& asym   < CtsMCut (asymup,    dmls, dmcza, dmls2, dd2)
       //&& asym   > CtsMCut (asymlow,   dmls, dmcza, dmls2, dd2)
       && dist     < CtsMCut (distup,    dmls, dmcza, dmls2, dd2)
       && dist     > CtsMCut (distlow,   dmls, dmcza, dmls2, dd2)  )
    fHadronness->SetHadronness(1.e-10);
  else
    fHadronness->SetHadronness(0.9999999999);

  return kTRUE;
}
// ---------------------------------------------------------------------------