source: trunk/MagicSoft/Mars/mimage/M2dimFunction.cc@ 2521

/* ======================================================================== *\
!
! *
! * 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 10/2003 <mailto:wittek@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// M2dimFunction
//
// Storage Container for the parameters of the 2-dim function describing
// the shower image
//
//
/////////////////////////////////////////////////////////////////////////////
#include "M2dimFunction.h"

#include <fstream>>

#include "TMath.h"
#include "TVectorD.h"
#include "TMatrixD.h"

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

#include "MHillas.h"

#include "MGeomCam.h"
#include "MGeomPix.h"

#include "MCerPhotEvt.h"
#include "MCerPhotPix.h"

ClassImp(M2dimFunction);

using namespace std;

// --------------------------------------------------------------------------
//
// TwodimFunction
//
// this function calculates for a given set of parameters
//
//   - the log likelihood function L to be minimized
//   - beta_k  = -1/2 * dL/da_k          (a kind of gradient of L)
//   - alfa_kl =  1/2 * dL/(da_k da_l)   (a kind of second derivative of L)
//
// (see Numerical recipes (2nd ed.), W.H.Press et al., p. 687)
//
// Note : this is not a member function of M2dimFunction;
//        it will be called by the minimization class MMarquardt;
//        the address of this function is passed to MMarquardt 
//        in M2dimFunction::Fit() by the call MMarquardt::Loop()
//

Bool_t TwodimFunction(TVectord &a, TMatrixD &alfa, TVectord &beta, Double_t &L)
{
  Int_t npar    = a.GetNrows();
  Int_t npixels =;

  TMatrixD dmuda  (npar, npixels);
  TVector  dLdmu  (npixels);
  TVector  d2Ldmu2(npixels);

  //------------------------------------------
  // these are the parameters for which alfa, beta and L are calculated
 
  Double_t xbar  = a(0);
  Double_t ybar  = a(1);
  Double_t delta = a(2);
  Double_t amp   = a(3);
  Double_t leng  = a(4);
  Double_t wid   = a(5);
  Double_t asy   = a(6);


  for (Int_t m=0; m<npar; m++)
  {
    beta(m) = 0.0;
    for (Int_t n=0; n<=m; n++)
    {
      alfa(m,n) = 0.0;
    }
  }

  //------------------------------------------
  // loop over the pixels
  //
  // quantities for each pixel :
  //
  // ar        pixel area
  // b         quantum efficiency * geometrical acceptance of 1st dynode
  // c         = ar * b
  // q         measured number of photo electrons
  //           (after subtraction of the pedestal)  
  // S         'measured' number of photo electrons (including the NSB)
  // mu        fitted number of  Cherenkov photons per area
  // lambda    average number of NSB photons per area 
  //           (obtained from the pedestal fluctuations)
  // sigma_el  sigma of electronic noise
  // F         the probability of measuring S


  L = 0.0;
  Double_t Fexcessnoise2 = 1.2 * 1.2;
  for (Int_t i=0; i<npixels; i++)
  {
    Double_t lambda   =
    Double_t sigma_el =

    Double_t S        = 
    Double_t mu = 2-dim function evaluated for pixel i, at the position (x,y);

    Double_t F      = 0.0;
    Double_t dFdmu  = 0.0;
    Double_t dFdmu2 = 0.0;

    // range of n for which Poisson and Gaus are not too small
    Int_t nmin =
    Int_t nmax =

    for (Int_t n=nmin; n<=nmax; n++)
    {
      Double_t sigma_n = sqrt(   n*(Fecessnoise2-1.0) + sigma_el*sigma_el );
      Double_t probn =   TMath::Poisson(n, c*(mu+lambda)) 
                       * TMath::Gaus(S, n, sigma_n, kTRUE);
      Double_t brack = n/(mu+lambda)-c;

      F      += probn;
      dFdmu  += probn * brack;
      dFdmu2 += probn * (brack * brack - n/( (mu+lambda)*(mu+lambda) );
    }

    // log-likelihood function 
    L -= 2.0 * log(F);

    // derivatives of log-likelihood function, up to factors of 2) 
    dLdmu(i)   = dFdmu / F;
    d2Ldmu2(i) = dFdmu*dFdmu / (F*F) - dFdmu2 / F;


    // derivatives of 2-dim function mu
    // - with respect to xbar, ybar and delta :
    dmudu =
    dmudv =
    dmuda(0,i) = -dmudu * cos(delta)  +  dmudv * sin(delta);
    dmuda(1,i) = -dmudu * sin(delta)  -  dmudv * cos(delta);
    dmuda(2,i) =  dmudu * ( -(x-xbar)*sin(delta) +(y-ybar)*cos(delta) )
                + dmudv * ( -(x-xbar)*cos(delta) -(y-ybar)*sin(delta) );
            
    // - with respect to the other variables :
    dmuda(3,i) = 
    dmuda(4,i) = 
    dmuda(5,i) = 
    dmuda(6,i) = 
  }


  //------------------------------------------
  // calculate alfa and beta
  // 
  // sum over all pixels
  for (Int_t i=0; i<npixels; i++)
  {
    for (Int_t m=0; m<npar; m++)
    {
      Double_t g = dmuda(m,i);
      for (Int_t n=0; n<=m; n++)
      {
        alfa(m,n) += d2Ldmu2(i) * g * dmuda(n,i);
      }
      beta(m) += dLdmu(i) * g;
    }
  }

  return kTRUE;
}

// --------------------------------------------------------------------------
//
// Default constructor.
//
M2dimFunction::M2dimFunction(const char *name, const char *title)
{
    fName  = name  ? name  : "M2dimFunction";
    fTitle = title ? title : "Parameters of 2-dim function";
}


// --------------------------------------------------------------------------
//
//  SetVinit
//
//  set the initial values of the parameters
//
void M2dimFunction::SetVinit(MHillas *fhillas)
{
  // get some initial values from the Hillas class

  if (fhillas)
  {
    fVinit(0) = fhillas->GetMeanX();
    fVinit(1) = fhillas->GetMeanY();
    fVinit(2) = fhillas->GetDelta();
  }
  else
  {
    fVinit(0) = ;
    fVinit(1) = ;
    fVinit(2) = ;
  }

  fVinit(3) =
  fVinit(4) =
  fVinit(5) =
  fVinit(6) =

  return;
} 



// --------------------------------------------------------------------------
//
//  Fit of the 2-dim function to the shower image
//
void M2dimFunction::Fit()
{
  fMarquardt.Loop(TwodimFunction, fVinit);

  SetReadyToSave();
} 

// --------------------------------------------------------------------------
//
void M2dimFunction::Print(Option_t *) const
{
    *fLog << all;
    *fLog << "Parameters of 2-dim function (" << GetName() << ")" << endl;
    *fLog << " - fXbar           = " << fXbar  << endl;
    *fLog << " - fYbar           = " << fYbar  << endl;
    *fLog << " - fAmp            = " << fAmp   << endl;
    *fLog << " - fMajor          = " << fMajor << endl;
    *fLog << " - fMinor          = " << fMinor << endl;
    *fLog << " - fAsym           = " << fAsym  << endl;
}
//============================================================================