source: branches/Mars_McMismatchStudy/manalysis/MEnergyTable.cc@ 18303

/* ======================================================================== *\
!
! *
! * 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): Thomas Bretz  9/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
!              Wolfgang Wittek  7/2003 <mailto:wittek@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
// MEnergyTable                                                            //
// Energy Tables generated at discrete Zenith and threshold                //
// We then interpolate between them to estimate the energy                 //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////
#include "MEnergyTable.h"

#include "MParList.h"

#include "MGeomCam.h"
#include "MMcEvt.hxx"
#include "MPointingPos.h"
#include "MHMatrix.h"
#include "MHillas.h"
#include "MHillasSrc.h"
#include "MNewImagePar.h"
#include "MEnergyEstimate.h"
#include "MParameters.h"

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

#include "TFile.h"

ClassImp(MEnergyTable);

// --------------------------------------------------------------------------
//
// Default constructor. Give the name of the parameter container (MHillas)
// storing width, length and size (Default="MHillas").
// For the Zenith Angle MMcEvt.fTheta is used.
//
MEnergyTable::MEnergyTable(const char *name, const char *title)
{

    fName  = name  ? name  : "MEnergyTable";
    fTitle = title ? title : "Task to estimate the energy using MC Tables";

    fNameHillasSrc = "MHillasSrc";
    fNameHillas = "MHillas";
    fNameMC = "MMcEvt";
    fNamePoint = "MPointingPos";

    //fNameParameter = "MParameterD";

    fNameParameter  = "MEnergyTable";
    fNameParameterL = "ScaledLength";
    fNameParameterW = "ScaledWidth";

    SetNameParameter("MEnergyTable");
}


// --------------------------------------------------------------------------
//
// Destructor.
//
MEnergyTable::~MEnergyTable()
{

    delete fHillasSrc;
    delete fHillas;
    delete fParameter;
    delete fParameterL;
    delete fParameterW;
    delete fPointing;

    fTableFile->Close();

}

Int_t MEnergyTable::PreProcess(MParList *pList)
{

    fHillasSrc = (MHillasSrc*)pList->FindCreateObj("MHillasSrc", AddSerialNumber(fNameHillasSrc));
    if (!fHillasSrc)
        return kFALSE;

    fHillas = (MHillas*)pList->FindCreateObj("MHillas", AddSerialNumber(fNameHillas));
    if (!fHillas)
        return kFALSE;

    fPointing = (MPointingPos*)pList->FindCreateObj("MPointingPos", AddSerialNumber(fNamePoint));
    if (!fPointing)
        return kFALSE;

    fParameter = (MParameterD*)pList->FindCreateObj("MParameterD", fNameParameter);
    if (!fParameter)
        return kFALSE;

    fParameterW = (MParameterD*)pList->FindCreateObj("MParameterD", fNameParameterW);
    if (!fParameterW)
        return kFALSE;

    fParameterL = (MParameterD*)pList->FindCreateObj("MParameterD", fNameParameterL);
    if (!fParameterL)
        return kFALSE;


    return kTRUE;

}

Int_t MEnergyTable::Process()
{

  Double_t fDist = fHillasSrc->GetDist();
  Double_t fSize = fHillas->GetSize();
  fZenith = fPointing->GetZd();
  fThreshold = 305.; //TODO

  iZenith = 0;
  iThreshold = 0;

  for( int i = 0; i < 2; i++ )
    if( fZenith >= iZen[i] && fZenith < iZen[i+1] ) iZenith = i;

  for( int i = 0; i < 2; i++ )
    if( fThreshold >= iThr[i] && fThreshold < iThr[i+1] ) iThreshold = i;

  if( iZenith != -999 && iThreshold != -999 )
  {

// Cannot get Threshold Yet so doing 1D interpolation on Zenith angle
// Between i and i+1 bin

    iX = hEnergy[iZenith][0]->GetXaxis()->FindBin(log10(fSize));
    iY = hEnergy[iZenith][0]->GetYaxis()->FindBin(fDist);
    fEnergy1 = hEnergy[iZenith][0]->GetBinContent(iX,iY);
    fLength1 = hLength[iZenith][0]->GetBinContent(iX,iY);
    fLengthSigma1 = hLength[iZenith][0]->GetBinContent(iX,iY);
    fWidth1 = hLength[iZenith][0]->GetBinContent(iX,iY);
    fWidthSigma1 = hLength[iZenith][0]->GetBinContent(iX,iY);


    iX = hEnergy[iZenith+1][0]->GetXaxis()->FindBin(log10(fSize));
    iY = hEnergy[iZenith+1][0]->GetYaxis()->FindBin(fDist);
    fEnergy2 = hEnergy[iZenith+1][0]->GetBinContent(iX,iY);
    fLength2 = hLength[iZenith+1][0]->GetBinContent(iX,iY);
    fLengthSigma2 = hLength[iZenith+1][0]->GetBinContent(iX,iY);
    fWidth2 = hLength[iZenith+1][0]->GetBinContent(iX,iY);
    fWidthSigma2 = hLength[iZenith+1][0]->GetBinContent(iX,iY);

    fEnergyRecon = fEnergy2-fEnergy1;
    fEnergyRecon /= cos(D2R*iZen[iZenith+1]) - cos(D2R*iZen[iZenith]);
    fEnergyRecon *= cos(D2R*fZenith) - cos(D2R*iZen[iZenith]);
    fEnergyRecon += fEnergy1;

    fLength = fLength2-fLength1;
    fLength /= cos(D2R*iZen[iZenith+1]) - cos(D2R*iZen[iZenith]);
    fLength *= cos(D2R*fZenith) - cos(D2R*iZen[iZenith]);
    fLength += fLength1;

    fLengthSigma = fLengthSigma2-fLengthSigma1;
    fLengthSigma /= cos(D2R*iZen[iZenith+1]) - cos(D2R*iZen[iZenith]);
    fLengthSigma *= cos(D2R*fZenith) - cos(D2R*iZen[iZenith]);
    fLengthSigma += fLengthSigma1;

    fWidth = fWidth2-fWidth1;
    fWidth /= cos(D2R*iZen[iZenith+1]) - cos(D2R*iZen[iZenith]);
    fWidth *= cos(D2R*fZenith) - cos(D2R*iZen[iZenith]);
    fWidth += fWidth1;

    fWidthSigma = fWidthSigma2-fWidthSigma1;
    fWidthSigma /= cos(D2R*iZen[iZenith+1]) - cos(D2R*iZen[iZenith]);
    fWidthSigma *= cos(D2R*fZenith) - cos(D2R*iZen[iZenith]);
    fWidthSigma += fWidthSigma1;


    RSW = (fHillas->GetWidth()-fWidth)/fWidthSigma;
    RSL = (fHillas->GetLength()-fLength)/fLengthSigma;

  } 

   fParameter->SetVal(fEnergyRecon);
   fParameter->SetReadyToSave();

   fParameterW->SetVal(RSW);
   fParameterW->SetReadyToSave();
   fParameterL->SetVal(RSL);
   fParameterL->SetReadyToSave();

   return kTRUE;

}

Int_t MEnergyTable::SetTableFile( const TString sTableFile )
{


   fTableFile = new TFile(sTableFile,"OPEN");
   if( fTableFile->IsZombie() )
   {
     cout << "Cannot open Table File: " << sTableFile << endl;
     return kFALSE;
   }

// TODO
   iZen[0] = 5;
   iZen[1] = 30;
   iZen[2] = 45;
   iThr[0] = 300;
   iThr[1] = 405;
   iThr[2] = 500;

   for( int i = 0; i < 3; i++ )
   {
     for( int j = 0; j < 3; j++ )
     {

        sprintf(sHistDir,"z%.2d/t%.2d",iZen[i],iThr[j]);
        sprintf(sHistName,"z%.2d/t%.2d/hEnergyAverage_%.2d_%.2d",iZen[i],iThr[j],iZen[i],iThr[j]);
        sprintf(sHistNameL,"z%.2d/t%.2d/hLengthAverage_%.2d_%.2d",iZen[i],iThr[j],iZen[i],iThr[j]);
        sprintf(sHistNameLS,"z%.2d/t%.2d/hLengthSigma_%.2d_%.2d",iZen[i],iThr[j],iZen[i],iThr[j]);
        sprintf(sHistNameW,"z%.2d/t%.2d/hWidthAverage_%.2d_%.2d",iZen[i],iThr[j],iZen[i],iThr[j]);
        sprintf(sHistNameWS,"z%.2d/t%.2d/hWidthSigma_%.2d_%.2d",iZen[i],iThr[j],iZen[i],iThr[j]);

        if( fTableFile->GetDirectory(sHistDir) ) 
        {
          fTableFile->cd(sHistDir);
          cout << "Adding: " << i << " " << j << " " << sHistDir << endl;

          htempE.push_back(  (TH2F*)fTableFile->Get(sHistName) );
          htempL.push_back(  (TH2F*)fTableFile->Get(sHistNameL) );
          htempLS.push_back( (TH2F*)fTableFile->Get(sHistNameLS) );
          htempW.push_back(  (TH2F*)fTableFile->Get(sHistNameW) );
          htempWS.push_back( (TH2F*)fTableFile->Get(sHistNameWS) );

        } else
        {
          cout << "Error Table File not complete!" << endl;
          return kFALSE;
        }
     }

     hEnergy.push_back(htempE);
     hLength.push_back(htempL);
     hLengthSigma.push_back(htempLS);
     hWidth.push_back(htempW);
     hWidthSigma.push_back(htempWS);

     htempE.clear();
     htempL.clear();
     htempLS.clear();
     htempW.clear();
     htempWS.clear();

   }

   return kTRUE;

}