/* ======================================================================== *\
!
! *
! * 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
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
//   MDisp                                                                 //
//                                                                         //
//   this is the container for the Disp parameterization                   //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////
#include "MDisp.h"

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

#include <TArrayD.h>

#include "MImageParDisp.h"

#include "MParList.h"

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


ClassImp(MDisp);

using namespace std;

// --------------------------------------------------------------------------
//
// constructor
//
MDisp::MDisp(const char *name, const char *title)
  : fLogSize0(3.), fLength0(0.1), fWidth0(0.05), fConc0(0.35),
    fLeakage10(0.05),fLeakage20(0.1)
{
    fName  = name  ? name  : "MDisp";
    fTitle = title ? title : "Container for the Disp parameters";

    // set arrays size to the number of parameters used
    fParameters.Set(5);
    fStepsizes.Set(5);

    // set Disp parameters to their default values
    InitParameters();
}


// --------------------------------------------------------------------------
//
// set default values for the Disp parameters
//
void MDisp::InitParameters()
{
    //---------------------------------------------------
    //  these are the default values

    fParameters[0] =  1.0;
    fParameters[1] =  0.6;
    fParameters[2] = -0.8;
    fParameters[3] = -0.8;
    fParameters[4] = -1.2;


    //---------------------------------------------------
    // fStepsizes 
    // if == 0.0    the parameter will be fixed in the minimization
    //    != 0.0    initial step sizes for the parameters

    fStepsizes[0] = 0.010;
    fStepsizes[1] = 0.006;
    fStepsizes[2] = 0.008;
    fStepsizes[3] = 0.008;
    fStepsizes[4] = 0.012;
}


// --------------------------------------------------------------------------
//
// Set the parameter values from the array 'd'
//
//
Bool_t MDisp::SetParameters(const TArrayD &d)
{
    if (d.GetSize() != fParameters.GetSize())
    {
        *fLog << err << "Sizes of d and of fParameters are different : "
              << d.GetSize() << ",  " << fParameters.GetSize() << endl;
        return kFALSE;
    }

    fParameters = d;

    return kTRUE;
}


// --------------------------------------------------------------------------
//
// Set the step sizes from the array 'd'
//
//
Bool_t MDisp::SetStepsizes(const TArrayD &d)
{
    if (d.GetSize() != fStepsizes.GetSize())
    {
        *fLog << err << "Sizes of d and of fStepsizes are different : "
              << d.GetSize() << ",  " << fStepsizes.GetSize() << endl;
        return kFALSE;
    }

    fStepsizes = d;

    return kTRUE;
}


// --------------------------------------------------------------------------
//
// Set the Disp parameterization and Calculate Disp
//
//
Double_t MDisp::Calc(TArrayD &imagevar)
{
    // get parameters
    Int_t numpar = fParameters.GetSize();
    TArrayD p(numpar);
    for (Int_t i=0; i<numpar; i++)
      p[i] = fParameters[i];

    // get image variables to be used in the Disp parameterization
    Double_t logsize  = imagevar[0];
    Double_t width    = imagevar[1];
    Double_t length   = imagevar[2];
    //    Double_t conc     = imagevar[3];
    //    Double_t leakage1 = imagevar[4];
    //    Double_t leakage2 = imagevar[5];
    
    // Disp parameterization to be optimized
    //  Note: fLogSize0, fLength0... variables are introduced to uncorrelate as much
    //        as possible the parameters in the Disp expression, with the purpose of
    //        helping the minimization algorithm to converge. They are set approx.
    //        to their distribution mean value in the MDisp constructor, but can be 
    //        changed using the corresponding set function.
    //

    //    Double_t disp = p[0] + p[1]*(logsize-fLogSize0) 
    //      + (p[2] + p[3]*(logsize-fLogSize0))*width/length;

    Double_t disp = p[0] + p[1]*(logsize-fLogSize0) + p[4]*(length-fLength0) 
      + (p[2] + p[3]*(logsize-fLogSize0))*width/length;

    //    Double_t disp = p[0] + p[1]*(logsize-fLogSize0) + p[4]*(length-fLength0) 
    //      + (p[2] + p[3]*(logsize-fLogSize0))*length/width;

    //    Double_t disp = p[0] + p[1]*(logsize-fLogSize0) + p[4]*(length-fLength0) 
    //      + (p[2] + p[3]*(logsize-fLogSize0) + p[5]*(length-fLength0))*width/length;

    //    Double_t disp = p[0] + p[1]*(logsize-fLogSize0) + p[4]*(width-fWidth0) 
    //      + (p[2] + p[3]*(logsize-fLogSize0))*width/length;   // + p[5]*(width-fWidth0))*width/length;

    //    Double_t disp = p[0] + p[1]*(logsize-fLogSize0) + p[4]*(conc-fConc0) 
    //      + (p[2] + p[3]*(logsize-fLogSize0))*width/length;   // + p[5]*(conc-fConc0))*width/length;

    //    Double_t disp = ( p[0] + p[1]*(logsize-fLogSize0) 
    //		      + p[2]*pow(logsize-fLogSize0,2)
    //		      + p[3]*pow(logsize-fLogSize0,3) 
    //		      + p[4]*pow(logsize-fLogSize0,4) )
    //                    *( 1/(1+width/length) );

    //    Double_t disp = ( p[0] + p[1]*(logsize-fLogSize0) 
    //		      + p[2]*pow(logsize-fLogSize0,2)
    //		      + p[3]*pow(logsize-fLogSize0,3) 
    //		      + p[4]*pow(logsize-fLogSize0,4) )
    //                    *(1-width/length);

    return disp;
}


