source: trunk/MagicSoft/Mars/mtemp/mifae/library/MDispCalc.cc@ 5775

/* ======================================================================== *\
!
! *
! * 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): Eva Domingo    , 12/2004 <mailto:domingo@ifae.es>
!              Wolfgang Wittek, 12/2004 <mailto:wittek@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2005
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
//   MDispCalc                                                             //
//                                                                         //
//   This task calculates Disp                                             //
//   (parameters are taken from the container MDisp)                       //
//   Also the matrix of variables to be used in the Disp optimization      //
//   is defined                                                            //
//                                                                         //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////
#include "MDispCalc.h"

#include <math.h>
#include <TArray.h>

#include "MGeomCam.h"
#include "MSrcPosCam.h"
#include "MHillas.h"
#include "MHillasExt.h"
#include "MNewImagePar.h"
#include "MMcEvt.hxx"
#include "MPointingPos.h"
#include "MImageParDisp.h"
#include "MDisp.h"

#include "MHMatrix.h"

#include "MParList.h"

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


ClassImp(MDispCalc);

using namespace std;

// --------------------------------------------------------------------------
//
// constructor
//
MDispCalc::MDispCalc(const char *imagepardispname, const char *dispname, 
                     const char *name, const char *title)
  : fImageParDispName(imagepardispname), fDispName(dispname)
{
    fName  = name  ? name  : "MDispCalc";
    fTitle = title ? title : "Class to calculate Disp";

    fMatrix = NULL;

    // initialize mapping array dimension to the number of columns of fMatrix
    fMap.Set(18);
}


// --------------------------------------------------------------------------
//
Int_t MDispCalc::PreProcess(MParList *pList)
{
    // look for the defined camera geometry to get mm to deg conversion factor
    MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
    if (!cam)
    {
        *fLog << err << "MGeomCam (Camera Geometry) not found... aborting." 
              << endl;
        return kFALSE;
    }

    fMm2Deg = cam->GetConvMm2Deg();


    // look for Disp related containers
    fImageParDisp = (MImageParDisp*)pList->FindCreateObj("MImageParDisp", 
                                                         fImageParDispName);
    if (!fImageParDisp)
    {
        *fLog << err << fImageParDispName 
              << " [MImageParDisp] not found... aborting." << endl;
        return kFALSE;
    }

    fDisp = (MDisp*)pList->FindObject(fDispName, "MDisp");
    if (!fDisp)
    {
        *fLog << err << fDispName << " [MDisp] not found... aborting." << endl;
        return kFALSE;
    }


    //-----------------------------------------------------------
    // if fMatrix exists, skip preprocessing and just read events from matrix;
    // if doesn't exist, read variables values from containers, so look for them
    if (fMatrix)
        return kTRUE;

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

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

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

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

    fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
    if (!fMcEvt)
    {
        *fLog << err << "MMcEvt not found... This is not a MC file,"
              << " you are not trying to optimize Disp, just calculating it."
              << endl;
        //        return kFALSE;
    }

    fPointing = (MPointingPos*)pList->FindObject("MPointingPos");
    if (!fPointing)
    {
        *fLog << err << "MPointingPos not found... " 
              << "MMcEvt is going to be used to get Theta and Phi." 
              << endl;
        //        return kFALSE;
    }

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

    return kTRUE;
}


// --------------------------------------------------------------------------
//
// You can use this function if you want to use a MHMatrix instead of the
// given containers for the Disp optimization. This function adds all 
// necessary columns to the given matrix. Afterwards, you should fill 
// the matrix with the corresponding data (eg from a file by using 
// MHMatrix::Fill). Then, if you loop through the matrix (eg using 
// MMatrixLoop), MDispCalc::Process will take the values from the matrix 
// instead of the containers.
//
void MDispCalc::InitMapping(MHMatrix *mat)
{
    if (fMatrix)
      return;

    fMatrix = mat;

    fMap[0]  = fMatrix->AddColumn("MSrcPosCam.fX");
    fMap[1]  = fMatrix->AddColumn("MSrcPosCam.fY");
    fMap[2]  = fMatrix->AddColumn("MHillas.fMeanX");
    fMap[3]  = fMatrix->AddColumn("MHillas.fMeanY");
    fMap[4]  = fMatrix->AddColumn("MHillas.fDelta");

    fMap[5]  = fMatrix->AddColumn("MHillas.fSize");
    fMap[6]  = fMatrix->AddColumn("MHillas.fWidth");
    fMap[7]  = fMatrix->AddColumn("MHillas.fLength");

    fMap[8]  = fMatrix->AddColumn("MNewImagePar.fConc");
    fMap[9]  = fMatrix->AddColumn("MNewImagePar.fLeakage1");
    fMap[10] = fMatrix->AddColumn("MNewImagePar.fLeakage2");

    fMap[11] = fMatrix->AddColumn("MHillasExt.fM3Long");
    fMap[12] = fMatrix->AddColumn("MHillasExt.fAsym");

    fMap[13]  = fMatrix->AddColumn("MMcEvt.fEnergy");
    fMap[14]  = fMatrix->AddColumn("MMcEvt.fImpact");
    fMap[15]  = fMatrix->AddColumn("MMcEvt.fLongitmax");
    
    if (fPointing)
    {
      fMap[16]  = fMatrix->AddColumn("MPointingPos.fZd");
      fMap[17]  = fMatrix->AddColumn("MPointingPos.fAz");
    }
    else
    {
      fMap[16]  = fMatrix->AddColumn("MMcEvt.fTelescopeTheta*kRad2Deg");
      fMap[17]  = fMatrix->AddColumn("MMcEvt.fTelescopePhi*kRad2Deg");
    }
}


// --------------------------------------------------------------------------
//
// Returns the Matrix and the mapped value for each Matrix column
//
MHMatrix* MDispCalc::GetMatrixMap(TArrayI &map)
{
    map = fMap;

    return fMatrix;
}


// --------------------------------------------------------------------------
//
// Returns a mapped value from the Matrix
//
Double_t MDispCalc::GetVal(Int_t i) const
{
    Double_t val = (*fMatrix)[fMap[i]];

    //    *fLog << "MDispCalc::GetVal; i, fMatrix, fMap, val = "
    //          << i << ",  " << fMatrix << ",  " << fMap[i] << ",  " << val << endl;

    return val;
}


// ---------------------------------------------------------------------------
//
// Calculate Disp 
// 
//
Int_t MDispCalc::Process()
{
    // get variables needed to compute disp from the matrix or the containers
    const Double_t size     = fMatrix ? GetVal(5)  : fHil->GetSize();
    const Double_t width0   = fMatrix ? GetVal(6)  : fHil->GetWidth();
    const Double_t length0  = fMatrix ? GetVal(7)  : fHil->GetLength();
    const Double_t conc     = fMatrix ? GetVal(8)  : fNewPar->GetConc();
    const Double_t leakage1 = fMatrix ? GetVal(9)  : fNewPar->GetLeakage1();
    const Double_t leakage2 = fMatrix ? GetVal(10) : fNewPar->GetLeakage2();

    // convert to deg
    const Double_t width   = width0  * fMm2Deg;
    const Double_t length  = length0 * fMm2Deg;

    // create an array to pass the variables to MDisp::Calc
    Int_t numimagevars = 6;  
    TArrayD imagevars(numimagevars);
    imagevars[0]  = log10(size);
    imagevars[1]  = width;
    imagevars[2]  = length;
    imagevars[3]  = conc;
    imagevars[4]  = leakage1;
    imagevars[5]  = leakage2;

    // compute Disp
    Double_t disp = fDisp->Calc(imagevars); 

    // save value into MImageParDisp container
    fImageParDisp->SetDisp(disp);
    fImageParDisp->SetReadyToSave();

    return kTRUE;
}
//===========================================================================