Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MDispCalc.cc
===================================================================
--- trunk/MagicSoft/Mars/mtemp/mifae/library/MDispCalc.cc	(revision 5671)
+++ trunk/MagicSoft/Mars/mtemp/mifae/library/MDispCalc.cc	(revision 5671)
@@ -0,0 +1,294 @@
+/* ======================================================================== *\
+!
+! *
+! * 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;
+}
+//===========================================================================
+
+
+
+
Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MDispCalc.h
===================================================================
--- trunk/MagicSoft/Mars/mtemp/mifae/library/MDispCalc.h	(revision 5671)
+++ trunk/MagicSoft/Mars/mtemp/mifae/library/MDispCalc.h	(revision 5671)
@@ -0,0 +1,89 @@
+#ifndef MARS_MDispCalc
+#define MARS_MDispCalc
+
+#ifndef MARS_MTask
+#include "MTask.h"
+#endif
+
+#ifndef ROOT_TArrayD
+#include <TArrayD.h>
+#endif
+
+#ifndef ROOT_TArrayI
+#include <TArrayI.h>
+#endif
+
+class MSrcPosCam;
+class MHillas;
+class MHillasExt;
+class MNewImagePar;
+class MMcEvt;
+class MPointingPos;
+class MImageParDisp;
+class MDisp;
+class MHMatrix;
+class MParList;
+
+class MDispCalc : public MTask
+{
+private:
+
+    MSrcPosCam    *fSrcPos;
+    MHillas       *fHil;
+    MHillasExt    *fHilExt;
+    MNewImagePar  *fNewPar;
+    MMcEvt        *fMcEvt;
+    MPointingPos  *fPointing;      // NOTE: take (theta,phi) from this container
+                                   // when it is available (MC calibrated files 
+                                   // or data files)
+
+    MImageParDisp *fImageParDisp;  //! output container for Disp
+    MDisp         *fDisp;          //  container for Disp parameters
+
+    TString  fImageParDispName;    // name of container to store Disp
+    TString  fDispName;            // name of container for Disp parameters
+
+    Double_t fMm2Deg;              // conversion factor from mm to deg
+
+    MHMatrix *fMatrix;             // matrix defined to have as much columns as 
+                                   // variables wanted to be accessible during
+                                   // Disp optimization, and as much rows as 
+                                   // events used for the optimization
+                                   // (filled at MFindDisp)
+
+    TArrayI  fMap;                 // array to store the matrix mapping column
+                                   // numbers corresponding to each added variable
+
+
+    Int_t PreProcess(MParList *pList);
+    Int_t Process();
+
+    Double_t GetVal(Int_t i) const; // get value of variable entered in the matrix
+                                    // at column fMap[i]
+
+
+public:
+
+    MDispCalc(const char *imagepardispname = "MImageParDisp",
+	      const char *dispname         = "MDisp",
+              const char *name=NULL, const char *title=NULL);
+
+    void InitMapping(MHMatrix *mat);       // define the matrix of variables
+    MHMatrix* GetMatrixMap(TArrayI &map);  // get matrix and its mapping array
+
+    ClassDef(MDispCalc, 0) // Task to evaluate Disp
+};
+
+#endif
+
+
+
+
+
+
+
+
+
+
+
+
Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MFDisp.cc
===================================================================
--- trunk/MagicSoft/Mars/mtemp/mifae/library/MFDisp.cc	(revision 5671)
+++ trunk/MagicSoft/Mars/mtemp/mifae/library/MFDisp.cc	(revision 5671)
@@ -0,0 +1,291 @@
+/* ======================================================================== *\
+!
+! *
+! * 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
+!
+!
+\* ======================================================================== */
+
+//////////////////////////////////////////////////////////////////////////////
+//                                                                          //
+//  MFDisp                                                                  //
+//                                                                          //
+//  This is a class to set cuts to select an event sample                   //
+//  for the Disp optimization                                               //
+//                                                                          //
+//                                                                          //
+//////////////////////////////////////////////////////////////////////////////
+
+#include "MFDisp.h"
+
+#include "MGeomCam.h"
+#include "MHillas.h"
+#include "MHillasSrc.h"
+#include "MNewImagePar.h"
+
+#include "MParList.h"
+
+#include "MLog.h"
+#include "MLogManip.h"
+
+
+ClassImp(MFDisp);
+
+using namespace std;
+
+// --------------------------------------------------------------------------
+//
+// Default constructor.
+//
+MFDisp::MFDisp(const char *name, const char *title)
+{
+    fName  = name  ? name  : "MFDisp";
+    fTitle = title ? title : "Cuts for Disp optimization";
+
+    // default values of cuts
+    SetCuts(0, 1000,   0, 1000,   0, 1000, 
+            0.0, 1.e10,  0.0, 1.0,  0.0, 1.0,   0.0, 0.0);
+}
+
+
+// --------------------------------------------------------------------------
+//
+// Set the values for the cuts 
+// 
+void MFDisp::SetCuts(Int_t islandsmin,      Int_t islandsmax,
+                     Int_t usedpixelsmin,   Int_t usedpixelsmax, 
+                     Int_t corepixelsmin,   Int_t corepixelsmax,
+                     Float_t sizemin,       Float_t sizemax,
+                     Float_t leakage1min,   Float_t leakage1max,
+                     Float_t leakage2min,   Float_t leakage2max,
+                     Float_t lengthmin,     Float_t widthmin)
+{ 
+  fIslandsMin    = islandsmin;
+  fIslandsMax    = islandsmax;
+
+  fUsedPixelsMin = usedpixelsmin;
+  fUsedPixelsMax = usedpixelsmax;
+
+  fCorePixelsMin = corepixelsmin;
+  fCorePixelsMax = corepixelsmax;
+
+  fSizeMin       = sizemin;
+  fSizeMax       = sizemax;
+
+  fLeakage1Min   = leakage1min;
+  fLeakage1Max   = leakage1max;
+
+  fLeakage2Min   = leakage2min;
+  fLeakage2Max   = leakage2max;
+
+  fLengthMin     = lengthmin;
+  fWidthMin      = widthmin;
+}
+
+
+// --------------------------------------------------------------------------
+//
+Int_t MFDisp::PreProcess(MParList *pList)
+{
+    MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
+    if (!cam)
+    {
+        *fLog << err << "MGeomCam (Camera Geometry) not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    fMm2Deg = cam->GetConvMm2Deg();
+
+
+    fHil = (MHillas*)pList->FindObject("MHillas");
+    if (!fHil)
+    {
+        *fLog << err << "MHillas not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    fHilSrc = (MHillasSrc*)pList->FindObject("MHillasSrc");
+    if (!fHilSrc)
+    {
+        *fLog << err << "MHillasSrc not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    fNewImgPar = (MNewImagePar*)pList->FindObject("MNewImagePar");
+    if (!fNewImgPar)
+    {
+        *fLog << err << "MNewImagePar not found... aborting." << endl;
+        return kFALSE;
+    }
+
+
+    memset(fCut, 0, sizeof(fCut));
+
+
+    //--------------------
+    *fLog << inf << "MFDisp cut values :" << endl; 
+
+    *fLog << inf << "     fIslandsMin, fIslandsMax = "
+          << fIslandsMin << ",  " << fIslandsMax << endl;
+
+    *fLog << inf << "     fUsedPixelsMin, fUsedPixelsMax = "
+          << fUsedPixelsMin << ",  " << fUsedPixelsMax << endl;
+
+    *fLog << inf << "     fCorePixelsMin, fCorePixelsMax = "
+          << fCorePixelsMin << ",  " << fCorePixelsMax << endl;
+
+    *fLog << inf << "     fSizeMin, fSizeMax = " 
+          << fSizeMin << ",  " << fSizeMax << endl;
+
+    *fLog << inf << "     fLeakage1Min, fLeakage1Max = " 
+          << fLeakage1Min << ",  " << fLeakage1Max << endl;
+
+    *fLog << inf << "     fLeakage2Min, fLeakage2Max = " 
+          << fLeakage2Min << ",  " << fLeakage2Max << endl;
+
+    *fLog << inf << "     fLengthMin, fWidthMin = " 
+          << fLengthMin << ",  " << fWidthMin << endl;
+    //--------------------
+
+    return kTRUE;
+}
+
+
+// --------------------------------------------------------------------------
+//
+// Sum up the number of events passing each cut
+//
+Bool_t MFDisp::Set(Int_t rc)
+{
+    fResult = kTRUE;
+    fCut[rc]++;
+    return kTRUE;
+}
+
+
+// --------------------------------------------------------------------------
+//
+// Evaluate cuts for Disp optimization
+// if selections are fulfilled set fResult = kTRUE;
+// 
+Int_t MFDisp::Process()
+{
+    const Double_t length     = fHil->GetLength() * fMm2Deg;
+    const Double_t width      = fHil->GetWidth()  * fMm2Deg;
+    //const Double_t dist       = fHilSrc->GetDist()* fMm2Deg;
+    //const Double_t delta      = fHil->GetDelta()  * kRad2Deg;
+    const Double_t size       = fHil->GetSize();
+
+    const Double_t leakage1   = fNewImgPar->GetLeakage1();
+    const Double_t leakage2   = fNewImgPar->GetLeakage2();
+
+    const Int_t numusedpixels = fNewImgPar->GetNumUsedPixels();
+    const Int_t numcorepixels = fNewImgPar->GetNumCorePixels();
+    //const Int_t numislands    = fNewImgPar->GetNumIslands();
+
+    fResult = kFALSE;
+
+    //if  (numislands<fIslandsMin  || numislands>fIslandsMax )
+    //    return Set(1);
+
+    if  (numusedpixels<fUsedPixelsMin  || numusedpixels>fUsedPixelsMax )
+        return Set(2);
+
+    if  (numcorepixels<fCorePixelsMin  || numcorepixels>fCorePixelsMax )
+        return Set(3);
+
+    if (size<fSizeMin  ||  size>fSizeMax)
+        return Set(4);
+
+    if (leakage1<fLeakage1Min || leakage1>fLeakage1Max)
+        return Set(5);
+
+    if (leakage2<fLeakage2Min || leakage2>fLeakage2Max)
+        return Set(6);
+
+    if (length<fLengthMin || width<fWidthMin)
+        return Set(7);
+
+    fCut[0]++;
+    return kTRUE;
+}
+
+
+// --------------------------------------------------------------------------
+//
+//  Prints some statistics about the cuts selections
+//
+Int_t MFDisp::PostProcess()
+{
+    if (GetNumExecutions()==0)
+        return kTRUE;
+
+    *fLog << inf << endl;
+    *fLog << GetDescriptor() << " execution statistics:" << endl;
+
+    *fLog << dec << setfill(' ');
+    *fLog << " " << setw(7) << fCut[1] << " (" << setw(3);
+    *fLog << (int)(fCut[1]*100/GetNumExecutions()+0.5) ;
+    *fLog << "%) Evts skipped due to: No.of islands < " << fIslandsMin ;
+    *fLog << " or > " << fIslandsMax << endl;
+
+    *fLog << " " << setw(7) << fCut[2] << " (" << setw(3);
+    *fLog << (int)(fCut[2]*100/GetNumExecutions()+0.5) ;
+    *fLog << "%) Evts skipped due to: No.of used pixels < " << fUsedPixelsMin ;
+    *fLog << " or > " << fUsedPixelsMax << endl;
+
+    *fLog << " " << setw(7) << fCut[3] << " (" << setw(3);
+    *fLog << (int)(fCut[3]*100/GetNumExecutions()+0.5) ;
+    *fLog << "%) Evts skipped due to: No.of core pixels < " << fCorePixelsMin ;
+    *fLog << " or > " << fCorePixelsMax << endl;
+
+    *fLog << " " << setw(7) << fCut[4] << " (" << setw(3) ;
+    *fLog << (int)(fCut[4]*100/GetNumExecutions()+0.5) ;
+    *fLog << "%) Evts skipped due to: SIZE < " << fSizeMin;
+    *fLog << " or > " << fSizeMax << endl;
+
+    *fLog << " " << setw(7) << fCut[5] << " (" << setw(3) ;
+    *fLog << (int)(fCut[5]*100/GetNumExecutions()+0.5) ;
+    *fLog << "%) Evts skipped due to: Leakage1 < " << fLeakage1Min;
+    *fLog << " or > " << fLeakage1Max << endl;
+
+    *fLog << " " << setw(7) << fCut[6] << " (" << setw(3) ;
+    *fLog << (int)(fCut[6]*100/GetNumExecutions()+0.5) ;
+    *fLog << "%) Evts skipped due to: Leakage2 < " << fLeakage2Min;
+    *fLog << " or > " << fLeakage2Max << endl;
+
+    *fLog << " " << setw(7) << fCut[7] << " (" << setw(3) ;
+    *fLog << (int)(fCut[7]*100/GetNumExecutions()+0.5) ;
+    *fLog << "%) Evts skipped due to: LENGTH <= " << fLengthMin;
+    *fLog << " or WIDTH <= " << fWidthMin << endl;
+
+    *fLog << " " << fCut[0] << " (" ;
+    *fLog << (int)(fCut[0]*100/GetNumExecutions()+0.5) ;
+    *fLog << "%) Evts survived the Disp cuts!" << endl;
+    *fLog << endl;
+
+    return kTRUE;
+}
+
+
+
+
+
Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MFDisp.h
===================================================================
--- trunk/MagicSoft/Mars/mtemp/mifae/library/MFDisp.h	(revision 5671)
+++ trunk/MagicSoft/Mars/mtemp/mifae/library/MFDisp.h	(revision 5671)
@@ -0,0 +1,88 @@
+#ifndef MARS_MFDisp
+#define MARS_MFDisp
+
+#ifndef MARS_MFilter
+#include "MFilter.h"
+#endif
+
+class MHillas;
+class MHillasSrc;
+class MNewImagePar;
+
+class MFDisp : public MFilter
+{
+private:
+
+    MHillas      *fHil;
+    MHillasSrc   *fHilSrc;
+    MNewImagePar *fNewImgPar;
+
+    //--------------------------    
+
+    // cutting variables to select the sample
+    // for the Disp optimization
+
+    Int_t      fIslandsMin;
+    Int_t      fIslandsMax;
+
+    Int_t      fUsedPixelsMin;
+    Int_t      fUsedPixelsMax;
+
+    Int_t      fCorePixelsMin;
+    Int_t      fCorePixelsMax;
+
+    Float_t      fSizeMin;
+    Float_t      fSizeMax;
+
+    Float_t      fLeakage1Min;
+    Float_t      fLeakage1Max;
+
+    Float_t      fLeakage2Min;
+    Float_t      fLeakage2Max;
+
+    Float_t      fLengthMin;
+    Float_t      fWidthMin;
+
+    //--------------------------    
+
+    Double_t     fMm2Deg;     // conversion mm to degrees in camera
+
+    Int_t        fCut[8];     // array to save cuts statistics
+
+    Bool_t       fResult;
+
+    Int_t PreProcess(MParList *pList);
+    Int_t Process();
+    Int_t PostProcess();
+
+    Bool_t IsExpressionTrue() const  { return fResult; }
+
+    Bool_t Set(Int_t rc);
+
+public:
+
+    MFDisp(const char *name=NULL, const char *title=NULL);
+
+    void SetCuts(Int_t islandsmin,      Int_t islandsmax,
+                 Int_t usedpixelsmin,   Int_t usedpixelsmax,
+                 Int_t corepixelsmin,   Int_t corepixelsmax,
+                 Float_t sizemin,       Float_t sizemax, 
+                 Float_t leakage1min,   Float_t leakage1max,
+                 Float_t leakage2min,   Float_t leakage2max,
+                 Float_t lengthmin,     Float_t widthmin);
+
+    ClassDef(MFDisp, 0)   // Class to set cuts for Disp optimization
+};
+
+#endif
+
+
+
+
+
+
+
+
+
+
+
Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MFindDisp.cc
===================================================================
--- trunk/MagicSoft/Mars/mtemp/mifae/library/MFindDisp.cc	(revision 5671)
+++ trunk/MagicSoft/Mars/mtemp/mifae/library/MFindDisp.cc	(revision 5671)
@@ -0,0 +1,1073 @@
+/* ======================================================================== *\
+!
+! *
+! * 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
+!
+!
+\* ======================================================================== */
+
+/////////////////////////////////////////////////////////////////////////////
+//                                                                         //
+// MFindDisp                                                               //
+//                                                                         //
+// Class for otimizing the estimation of Disp                              //
+//                                                                         //
+//                                                                         //
+//                                                                         //
+/////////////////////////////////////////////////////////////////////////////
+#include "MFindDisp.h"
+
+#include <math.h>            // fabs 
+
+#include <TArrayD.h>
+#include <TCanvas.h>
+#include <TFile.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TMinuit.h>
+#include <TStopwatch.h>
+#include <TVirtualFitter.h>
+
+#include "MBinning.h"
+#include "MContinue.h"
+#include "MDisp.h"
+#include "MDispCalc.h"
+#include "MDataElement.h"
+#include "MDataMember.h"
+
+#include "MEvtLoop.h"
+#include "MFSelFinal.h"
+#include "MF.h"
+#include "MFDisp.h"
+#include "MFEventSelector.h"
+#include "MFEventSelector2.h"
+#include "MFillH.h"
+#include "MFEventSelector.h"
+#include "MGeomCamMagic.h"
+#include "MH3.h"
+#include "MHDisp.h"
+#include "MHFindSignificance.h"
+#include "MHMatrix.h"
+#include "MHOnSubtraction.h"
+
+#include "MLog.h"
+#include "MLogManip.h"
+#include "MMatrixLoop.h"
+#include "MMinuitInterface.h"
+#include "MParList.h"
+#include "MProgressBar.h"
+#include "MReadMarsFile.h"
+#include "MReadTree.h"
+#include "MStatusDisplay.h"
+#include "MTaskList.h"
+
+
+ClassImp(MFindDisp);
+
+using namespace std;
+
+
+//------------------------------------------------------------------------
+// fcnDisp 
+//------------------------------------------------------------------------
+//
+// - calculates the quantity to be minimized (by TMinuit)
+//
+// - the quantity to be minimized is defined in 
+//   MHDisp::Fill() and MHDisp::Finalize()
+//
+// - the parameters to be varied in the minimization are the parameters
+//   appearing in the parametrization of Disp (MDisp::Calc())
+//
+//------------------------------------------------------------------------
+
+static void fcnDisp(Int_t &npar, Double_t *gin, Double_t &f, 
+                    Double_t *par, Int_t iflag)
+{
+    cout <<  "entry fcnDisp" << endl;
+
+    // save pointer to the MINUIT oject (for optimizing Disp) for the case
+    // that it is overwritten by the pointer to another Minuit object
+    //    TMinuit *savePointer = gMinuit;
+    //-------------------------------------------------------------
+
+
+    MEvtLoop *evtloopfcn = (MEvtLoop*)gMinuit->GetObjectFit();
+
+    MParList *plistfcn   = (MParList*)evtloopfcn->GetParList();
+    MTaskList *tasklistfcn   = (MTaskList*)plistfcn->FindObject("MTaskList");
+
+    // get needed Disp containers from the existing parList
+    MDisp *disp = (MDisp*)plistfcn->FindObject("MDisp");
+    if (!disp)
+    {
+        gLog << "fcnDisp : MDisp object '" << "MDisp"
+            << "' not found... aborting" << endl;
+        return;
+    }
+
+    MHDisp *hdisp = (MHDisp*)plistfcn->FindObject("MHDisp");
+    if (!hdisp)
+    {
+        gLog << "fcnDisp : MHDisp object '" << "MHDisp"
+            << "' not found... aborting" << endl;
+        return;
+    }
+
+    //
+    // transfer current parameter values to MDisp
+    //
+    // Attention : npar is the number of variable parameters
+    //                  not the total number of parameters
+    //
+    Double_t fMin,   fEdm,   fErrdef;
+    Int_t    fNpari, fNparx, fIstat;
+    gMinuit->mnstat(fMin, fEdm, fErrdef, fNpari, fNparx, fIstat);
+
+    disp->SetParameters(TArrayD(fNparx, par));
+
+    // checking that parameters have been properly set
+    TArrayD checkparameters = disp->GetParameters();
+    gLog << "fcnDisp : fNpari, fNparx =" << fNpari << ",  " 
+         << fNparx  << endl;
+    gLog << "fcnDisp : i, par, checkparameters =" << endl;
+    for (Int_t i=0; i<fNparx; i++)
+    {
+      gLog << i << ",  " << par[i] << ",  " << checkparameters[i] << endl;
+    }
+
+    //
+    // loop over the events in the training matrix to compute the Chi^2
+    // for the current parameter values
+    evtloopfcn->Eventloop();
+
+    tasklistfcn->PrintStatistics(0, kTRUE);
+
+    //-------------------------------------------
+    // get the Chi^2 value for the current values of the parameters
+    f = hdisp->GetChi2();
+}
+
+
+// --------------------------------------------------------------------------
+//
+// Default constructor.
+//
+MFindDisp::MFindDisp(MFDisp *fdisp, const char *name, const char *title)
+  : fHowManyTrain(10000), fHowManyTest(10000)
+{
+    fName  = name  ? name  : "MFindDisp";
+    fTitle = title ? title : "Optimizer of Disp";
+
+    //---------------------------
+    // camera geometry is needed for conversion mm ==> degree
+    fCam = new MGeomCamMagic; 
+
+    // initialize MFDisp filter cuts to default (no cutting) values
+    // (done at MFDisp constructor)
+    fDispFilter = fdisp;
+    *fLog << inf << "MFindDisp::MFindDisp; address of MFDisp object = "
+          << fDispFilter << endl;    
+
+    // matrices to contain the training/test samples
+    fMatrixTrain = new MHMatrix("MatrixTrain");
+    fMatrixTest  = new MHMatrix("MatrixTest");
+
+    // objects of MDispCalc to which these matrices are attached
+    fDispCalcTrain = new MDispCalc("DispTrain");
+    fDispCalcTest  = new MDispCalc("DispTest");
+
+    // define columns of matrices
+    fDispCalcTrain->InitMapping(fMatrixTrain);
+    fDispCalcTest->InitMapping(fMatrixTest);
+}
+
+// --------------------------------------------------------------------------
+//
+// Default destructor.
+//
+MFindDisp::~MFindDisp()
+{
+    delete fCam;
+    delete fMatrixTrain;
+    delete fMatrixTest;
+    delete fDispCalcTrain;
+    delete fDispCalcTest;
+}
+
+
+// --------------------------------------------------------------------------
+//
+// Define the matrix 'fMatrixTrain' for the TRAINING sample
+//
+// alltogether 'howmanytrain' events are read from file 'nametrain';
+// the events are selected according to a target distribution 'hreftrain'
+//
+//
+Bool_t MFindDisp::DefineTrainMatrix(
+			  const TString &nametrain, MH3 &hreftrain,
+	                  const Int_t howmanytrain, const TString &filetrain,
+			  Int_t iflag)
+{
+    if (nametrain.IsNull() || howmanytrain <= 0)
+        return kFALSE;
+
+    *fLog   << "=============================================" << endl;
+    *fLog   << "Fill TRAINING Matrix from file '" << nametrain << endl;
+    *fLog   << "     select " << howmanytrain << " events " << endl;
+    if (!hreftrain.GetHist().GetEntries()==0)
+    {
+      *fLog << "     according to a distribution given by the MH3 object '"
+            << hreftrain.GetName() << "'" << endl;
+    }
+    else
+    {
+      *fLog << "     randomly" << endl;
+    }
+    *fLog   << "=============================================" << endl;
+
+
+    MParList  plist;
+    MTaskList tlist;
+
+    // initialize display to check the reference distribution
+    // for the event selection (just if iflag==1)
+    MStatusDisplay *display = NULL;
+    if (iflag)
+      display = new MStatusDisplay;
+
+    MReadMarsFile read("Events", nametrain);
+    read.DisableAutoScheme();
+
+    // apply cuts for event selection
+    MContinue contdisp(fDispFilter);
+    contdisp.SetName("ContFilterSelector2");
+
+    // choose a reference distribution
+    MFEventSelector2 seltrain(hreftrain);
+    seltrain.SetNumMax(howmanytrain);
+    seltrain.SetName("selectTrain");
+
+    MFillH filltrain(fMatrixTrain);
+    filltrain.SetFilter(&seltrain);
+    filltrain.SetName("fillMatrixTrain");
+
+    //******************************
+    // entries in MParList 
+    
+    plist.AddToList(&tlist);
+    plist.AddToList(fCam);
+    plist.AddToList(fMatrixTrain);
+
+    //******************************
+    // entries in MTaskList 
+
+    tlist.AddToList(&read);
+    if (fDispFilter != NULL)
+      tlist.AddToList(&contdisp);
+    tlist.AddToList(&seltrain);
+    tlist.AddToList(&filltrain);
+
+    //******************************
+
+    MProgressBar bar;
+    MEvtLoop evtloop;
+    if (display != NULL)
+      evtloop.SetDisplay(display);
+    evtloop.SetParList(&plist);
+    evtloop.SetName("EvtLoopMatrixTrain");
+    evtloop.SetProgressBar(&bar);
+
+    if (!evtloop.Eventloop())
+      return kFALSE;
+
+    tlist.PrintStatistics(0, kTRUE);
+
+
+    // print the filled Training Matrix
+    fMatrixTrain->Print("SizeCols");
+
+    // check that number of generated events is compatible with the resquested
+    Int_t howmanygenerated = fMatrixTrain->GetM().GetNrows();
+    if (TMath::Abs(howmanygenerated-howmanytrain) > TMath::Sqrt(9.*howmanytrain))
+    {
+      *fLog << "MFindDisp::DefineTrainMatrix; no.of generated events ("
+	    << howmanygenerated 
+            << ") is incompatible with the no.of requested events ("
+            << howmanytrain << ")" << endl;
+    }
+
+    *fLog << "TRAINING Matrix was filled" << endl;
+    *fLog << "=============================================" << endl;
+
+
+    //--------------------------
+    // write out training matrix
+
+    if (filetrain != "")
+    {
+      TFile filetr(filetrain, "RECREATE");
+      fMatrixTrain->Write();
+      filetr.Close();
+
+      *fLog << "MFindDisp::DefineTrainMatrix; Training matrix was written onto file '"
+            << filetrain << "'" << endl;
+    }
+
+    if (display != NULL)
+      delete display;
+
+    return kTRUE;
+}
+
+
+// --------------------------------------------------------------------------
+//
+// Define the matrix 'fMatrixTest' for the TEST sample
+//
+// alltogether 'howmanytest' events are read from file 'nametest'
+// the events are selected according to a target distribution 'hreftest'
+//
+//
+Bool_t MFindDisp::DefineTestMatrix(
+			  const TString &nametest, MH3 &hreftest,
+	                  const Int_t howmanytest, const TString &filetest,
+			  Int_t iflag)
+{
+    if (nametest.IsNull() || howmanytest<=0)
+        return kFALSE;
+
+    *fLog   << "=============================================" << endl;
+    *fLog   << "Fill TEST Matrix from file '" << nametest << endl;
+    *fLog   << "     select " << howmanytest << " events " << endl;
+    if (!hreftest.GetHist().GetEntries()==0)
+    {
+      *fLog << "     according to a distribution given by the MH3 object '"
+            << hreftest.GetName() << "'" << endl;
+    }
+    else
+    {
+      *fLog << "     randomly" << endl;
+    }
+
+
+    MParList  plist;
+    MTaskList tlist;
+
+    // initialize display to check the reference distribution
+    // for the event selection (just if iflag==1)
+    MStatusDisplay *display = NULL;
+    if (iflag)
+      display = new MStatusDisplay;
+
+    MReadMarsFile read("Events", nametest);
+    read.DisableAutoScheme();
+
+    // apply cuts for event selection
+    MContinue contdisp(fDispFilter);
+    contdisp.SetName("ContFilterSelector2");
+ 
+    // choose a reference distribution 
+    MFEventSelector2 seltest(hreftest);
+    seltest.SetNumMax(howmanytest);
+    seltest.SetName("selectTest");
+ 
+    MFillH filltest(fMatrixTest);
+    filltest.SetFilter(&seltest);
+    filltest.SetName("fillMatrixTest");
+
+    //******************************
+    // entries in MParList 
+    
+    plist.AddToList(&tlist);
+    plist.AddToList(fCam);
+    plist.AddToList(fMatrixTest);
+
+    //******************************
+    // entries in MTaskList 
+
+    tlist.AddToList(&read);
+    if (fDispFilter != NULL)
+      tlist.AddToList(&contdisp);
+    tlist.AddToList(&seltest);
+    tlist.AddToList(&filltest);
+
+    //******************************
+
+    MProgressBar bar;
+    MEvtLoop evtloop;
+    if (display != NULL)
+      evtloop.SetDisplay(display);
+    evtloop.SetParList(&plist);
+    evtloop.SetName("EvtLoopMatrixTest");
+    evtloop.SetProgressBar(&bar);
+
+    if (!evtloop.Eventloop())
+      return kFALSE;
+
+    tlist.PrintStatistics(0, kTRUE);
+
+
+    // print the filled Test Matrix
+    fMatrixTest->Print("SizeCols");
+
+    // check that number of generated events is compatible with the resquested
+    const Int_t howmanygenerated = fMatrixTest->GetM().GetNrows();
+    if (TMath::Abs(howmanygenerated-howmanytest) > TMath::Sqrt(9.*howmanytest))
+    {
+      *fLog << "MFindDisp::DefineTestMatrix; no.of generated events ("
+	    << howmanygenerated 
+            << ") is incompatible with the no.of requested events ("
+            << howmanytest << ")" << endl;
+    }
+
+    *fLog << "TEST Matrix was filled" << endl;
+    *fLog << "=============================================" << endl;
+
+
+    //--------------------------
+    // write out test matrix
+
+    if (filetest != "")
+    {
+      TFile filete(filetest, "RECREATE");
+      fMatrixTest->Write();
+      filete.Close();
+
+      *fLog << "MFindDisp::DefineTestMatrix; Test matrix was written onto file '"
+            << filetest << "'" << endl;
+    }
+
+    if (display != NULL)
+      delete display;
+
+    return kTRUE;
+}
+
+
+// --------------------------------------------------------------------------
+//
+// Define the matrices for the TRAINING and TEST sample respectively
+//
+//   - this is for the case that training and test samples are generated
+//     from the same input file (which has the name 'name')
+//
+Bool_t MFindDisp::DefineTrainTestMatrix(
+			  const TString &name, MH3 &href,
+	                  const Int_t howmanytrain, const Int_t howmanytest,
+                          const TString &filetrain, const TString &filetest,
+			  Int_t iflag)
+{
+    *fLog   << "=============================================" << endl;
+    *fLog   << "Fill TRAINING and TEST Matrix from file '" << name << endl;
+    *fLog   << "     select "   << howmanytrain 
+	    << " training and " << howmanytest << " test events " << endl;
+    if (!href.GetHist().GetEntries()==0)
+    {
+      *fLog << "     according to a distribution given by the MH3 object '"
+            << href.GetName() << "'" << endl;
+    }
+    else
+    {
+      *fLog << "     randomly" << endl;
+    }
+
+
+    MParList  plist;
+    MTaskList tlist;
+
+    // initialize display to check the reference distribution
+    // for the event selection (just if iflag==1)
+    MStatusDisplay *display = NULL;
+    if (iflag)
+      display = new MStatusDisplay;
+
+    MReadMarsFile read("Events", name);
+    read.DisableAutoScheme();
+
+    // apply cuts for event selection
+    MContinue contdisp(fDispFilter);
+    contdisp.SetName("ContFilterSelector2");
+
+    // choose a reference distribution 
+    MFEventSelector2 selector(href);
+    selector.SetNumMax(howmanytrain+howmanytest);
+    selector.SetName("selectTrainTest");
+    selector.SetInverted();
+    MContinue cont(&selector);
+    cont.SetName("ContTrainTest");
+
+    // choose randomly the events to fill the Training Matrix
+    Double_t prob =  ( (Double_t) howmanytrain )
+                   / ( (Double_t)(howmanytrain+howmanytest) );
+    MFEventSelector split;
+    split.SetSelectionRatio(prob);
+
+    MFillH filltrain(fMatrixTrain);
+    filltrain.SetFilter(&split);
+    filltrain.SetName("fillMatrixTrain");
+
+    // consider this event as candidate for a test event 
+    // only if event was not accepted as a training event
+    MContinue conttrain(&split);
+    conttrain.SetName("ContTrain");
+
+    MFillH filltest(fMatrixTest);
+    filltest.SetName("fillMatrixTest");
+
+
+    //******************************
+    // entries in MParList 
+    
+    plist.AddToList(&tlist);
+    plist.AddToList(fCam);
+    plist.AddToList(fMatrixTrain);
+    plist.AddToList(fMatrixTest);
+
+    //******************************
+    // entries in MTaskList 
+
+    tlist.AddToList(&read);
+    if (fDispFilter != NULL)
+      tlist.AddToList(&contdisp);
+    tlist.AddToList(&cont);
+    tlist.AddToList(&filltrain);
+    tlist.AddToList(&conttrain);
+    tlist.AddToList(&filltest);
+
+    //******************************
+
+    MProgressBar bar;
+    MEvtLoop evtloop;
+    if (display != NULL)
+      evtloop.SetDisplay(display);
+    evtloop.SetParList(&plist);
+    evtloop.SetName("EvtLoopMatrixTrainTest");
+    evtloop.SetProgressBar(&bar);
+
+    Int_t maxev = -1;
+    if (!evtloop.Eventloop(maxev))
+      return kFALSE;
+
+    tlist.PrintStatistics(0, kTRUE);
+
+
+    // print the filled Train Matrix
+    fMatrixTrain->Print("SizeCols");
+
+    // check that number of generated events is compatible with the resquested
+    const Int_t generatedtrain = fMatrixTrain->GetM().GetNrows();
+    if (TMath::Abs(generatedtrain-howmanytrain) > TMath::Sqrt(9.*howmanytrain))
+    {
+      *fLog << "MFindDisp::DefineTrainTestMatrix; no.of generated events ("
+	    << generatedtrain 
+            << ") is incompatible with the no.of requested events ("
+            << howmanytrain << ")" << endl;
+    }
+
+
+    // print the filled Test Matrix
+    fMatrixTest->Print("SizeCols");
+
+    // check that number of generated events is compatible with the resquested
+    const Int_t generatedtest = fMatrixTest->GetM().GetNrows();
+    if (TMath::Abs(generatedtest-howmanytest) > TMath::Sqrt(9.*howmanytest))
+    {
+      *fLog << "MFindDisp::DefineTrainTestMatrix; no.of generated events ("
+	    << generatedtest 
+            << ") is incompatible with the no.of requested events ("
+            << howmanytest << ")" << endl;
+    }
+
+
+    *fLog << "TRAINING and TEST Matrix were filled" << endl;
+    *fLog << "=============================================" << endl;
+
+
+    //--------------------------
+    // write out training matrix
+
+    if (filetrain != "")
+    {
+      TFile filetr(filetrain, "RECREATE");
+      fMatrixTrain->Write();
+      filetr.Close();
+
+      *fLog << "MFindDisp::DefineTrainTestMatrix; Training matrix was written onto file '"
+            << filetrain << "'" << endl;
+    }
+
+    //--------------------------
+    // write out test matrix
+
+    if (filetest != "")
+    {
+      TFile filete(filetest, "RECREATE");
+      fMatrixTest->Write();
+      filete.Close();
+
+      *fLog << "MFindDisp::DefineTrainTestMatrix; Test matrix was written onto file '"
+            << filetest << "'" << endl;
+    }
+
+    if (display != NULL)
+      delete display;
+
+    return kTRUE;
+}
+
+
+// --------------------------------------------------------------------------
+//
+// Read training and test matrices from files
+//
+//
+Bool_t MFindDisp::ReadMatrix(const TString &filetrain, const TString &filetest)
+{
+  //--------------------------
+  // read in training matrix
+
+  TFile filetr(filetrain);
+  fMatrixTrain->Read("MatrixTrain");
+  fMatrixTrain->Print("SizeCols");
+
+  *fLog << "MFindDisp::ReadMatrix; Training matrix was read in from file '"
+        << filetrain << "'" << endl;
+  filetr.Close();
+
+
+  //--------------------------
+  // read in test matrix
+
+  TFile filete(filetest);
+  fMatrixTest->Read("MatrixTest");
+  fMatrixTest->Print("SizeCols");
+
+  *fLog << "MFindDisp::ReadMatrix; Test matrix was read in from file '"
+        << filetest << "'" << endl;
+  filete.Close();
+
+  return kTRUE;  
+}
+
+
+//------------------------------------------------------------------------
+//
+// Steering program for optimizing Disp
+// ------------------------------------
+//
+//      the criterion for the 'optimum' is defined in 
+//      MHDisp::Fill()  and  MHDisp::Finalize(); 
+//      for example : smallest sum (over all events) of d^2, where d is the 
+//                    distance between the estimated source position 
+//                    (calculated using the current value of Disp) and
+//                    the true source position     
+//
+// The various steps are :
+//
+// - setup the event loop to be executed for each call to fcnDisp 
+//
+// - call TMinuit to do the minimization :
+//        the fcnDisp function calculates the Chi^2 for the current 
+//                                                parameter values;
+//        for this - Disp is calculated in the event loop by calling 
+//                   MDispCalc::Process() ==> MDisp::Calc() 
+//                 - the Chi^2 contributions are summed up in the event loop 
+//                   by calling MHDisp::Fill()
+//                 - after the event loop the Chi^2 is calculated by calling
+//                   MHDisp::Finalize()
+//
+// Needed as input : (to be set by the Set functions)
+//
+// - fFilenameParam      name of file to which optimum values of the 
+//                       parameters are written
+//
+// - for the minimization, the starting values of the parameters are taken  
+//     - from the file parDispInit (if it is != "")
+//     - or from the arrays params and/or steps 
+//     - or from the Disp constructor
+//
+//----------------------------------------------------------------------
+Bool_t MFindDisp::FindParams(TString parDispInit,
+                             TArrayD &params, TArrayD &steps)
+{
+    // Setup the event loop which will be executed in the 
+    //                 fcnDisp function  of MINUIT
+    //
+    // parDispInit is the name of the file containing the initial values 
+    // of the parameters; 
+    // if parDispInit = ""   'params' and 'steps' are taken as initial values
+    //
+
+    *fLog << "=============================================" << endl;
+    *fLog << "Setup event loop for fcnDisp" << endl;
+
+
+    if (fMatrixTrain == NULL)
+    {
+      *fLog << "MFindDisp::FindParams; training matrix is not defined... aborting"
+            << endl;
+      return kFALSE;
+    }
+
+    if (fMatrixTrain->GetM().GetNrows() <= 0)
+    {
+      *fLog << "MFindDisp::FindParams; training matrix has no entries"
+            << endl;
+      return kFALSE;
+    }
+
+    //---------------------------------------------------------
+    MParList  parlistfcn;
+    MTaskList tasklistfcn;
+
+    // loop over rows of matrix
+    MMatrixLoop loop(fMatrixTrain);
+
+    //--------------------------------
+    // create container for the Disp parameters
+    // and set them to their initial values
+    MDisp disp;
+
+    // take initial values from file parDispInit
+    if (parDispInit != "")
+    {
+      TFile inparam(parDispInit);
+      disp.Read("MDisp");
+      inparam.Close();
+      *fLog << "MFindDisp::FindParams; initial values of parameters are taken from file "
+            << parDispInit << endl;
+    }
+
+    // take initial values from 'params' and/or 'steps'
+    else if (params.GetSize() != 0  || steps.GetSize()  != 0 )
+    {
+      if (params.GetSize()  != 0)
+      {
+        *fLog << "MFindDisp::FindParams; initial values of parameters are taken from 'params'"
+              << endl;
+        disp.SetParameters(params);
+      }
+      if (steps.GetSize()  != 0)
+      {
+        *fLog << "MFindDisp::FindParams; initial step sizes are taken from 'steps'"
+              << endl;
+        disp.SetStepsizes(steps);
+      }
+    }
+    else
+    {
+        *fLog << "MFindDisp::FindParams; initial values and step sizes are taken "
+	      << "from the MDisp constructor" << endl;
+    }
+
+    // get the matrix pointer and mapping from MDispCalc object to set the same in MHDisp object
+    TArrayI map;
+    MHMatrix *tmpmatrix = NULL;
+    tmpmatrix = fDispCalcTrain->GetMatrixMap(map);
+
+    // attention: argument of MHDisp is the name of MImageParDisp container, that should
+    // be the same than the name given to it when creating MDispCalc object at the MFindDisp
+    // constructor:  fDispCalcTrain = new MDispCalc("DispTrain");
+    MHDisp hdisp("DispTrain");
+    hdisp.SetMatrixMap(tmpmatrix,map);
+    // fill the plots for Disp and sum up the Chi^2 contributions
+    MFillH filldispplots("MHDisp", "");
+
+    //******************************
+    // entries in MParList
+    
+    parlistfcn.AddToList(&tasklistfcn);
+    parlistfcn.AddToList(&disp);
+    parlistfcn.AddToList(&hdisp);
+    parlistfcn.AddToList(fCam);
+    parlistfcn.AddToList(fMatrixTrain);
+
+    //******************************
+    // entries in MTaskList
+
+    tasklistfcn.AddToList(&loop);
+    tasklistfcn.AddToList(fDispCalcTrain);
+    tasklistfcn.AddToList(&filldispplots);
+
+    //******************************
+
+    MEvtLoop evtloopfcn("EvtLoopFCN");
+    evtloopfcn.SetParList(&parlistfcn);
+    *fLog << "Event loop for fcnDisp has been setup" << endl;
+
+    // address of evtloopfcn is used in CallMinuit()
+
+
+    //-----------------------------------------------------------------------
+    //
+    //----------   Start of minimization part   --------------------
+    //
+    // Do the minimization with MINUIT
+    //
+    // Be careful: This is not thread safe
+    //
+    *fLog << "========================================================" << endl;
+    *fLog << "Start minimization for Disp" << endl;
+
+
+    // -------------------------------------------
+    // prepare call to MINUIT
+    //
+
+    // get initial values of parameters 
+    fVinit = disp.GetParameters();
+    fStep  = disp.GetStepsizes();
+
+    TString name[fVinit.GetSize()];
+    fStep.Set(fVinit.GetSize());
+    fLimlo.Set(fVinit.GetSize());
+    fLimup.Set(fVinit.GetSize());
+    fFix.Set(fVinit.GetSize());
+
+    fNpar = fVinit.GetSize();
+
+    // define names, step sizes, lower and upper limits of the parameters
+    // fFix[] = 0;  means the parameter is not fixed
+    for (UInt_t i=0; i<fNpar; i++)
+    {
+        name[i]   = "p";
+        name[i]  += i+1;
+        //fStep[i]  = TMath::Abs(fVinit[i]/10.0);
+        fLimlo[i] = -100.0;
+        fLimup[i] =  100.0;
+        fFix[i]   =     0;
+    }
+
+    // fix some parameters
+    //for (UInt_t i=0; i<fNpar; i++)
+    //{
+    //  if (i == 1  ||  i==3)
+    //	{
+    //      fStep[i] = 0.0;
+    //      fFix[i]  =   1;
+    //	}
+    //}
+ 
+
+    // -------------------------------------------
+    // call MINUIT
+
+    TStopwatch clock;
+    clock.Start();
+
+    *fLog << "before calling CallMinuit" << endl;
+
+    MMinuitInterface inter;               
+    Bool_t rc = inter.CallMinuit(fcnDisp, name,
+                                 fVinit, fStep, fLimlo, fLimup, fFix,
+                                 &evtloopfcn, "MIGRAD", kFALSE);
+ 
+    *fLog << "after calling CallMinuit" << endl;
+
+    *fLog << "Time spent for the minimization in MINUIT :   " << endl;;
+    clock.Stop();
+    clock.Print();
+
+
+    if (!rc)
+        return kFALSE;
+
+    *fLog << "Minimization for Disp finished" << endl;
+    *fLog << "========================================================" << endl;
+
+
+    // -----------------------------------------------------------------
+    // in 'fcnDisp' the optimum parameter values were put into MDisp
+
+    // write optimum parameter values onto file fFilenameParam
+    
+    TFile outparam(fFilenameParam, "RECREATE"); 
+    disp.Write();
+    outparam.Close();
+
+    *fLog << "Optimum parameter values for Disp were written onto file '"
+              << fFilenameParam << "' :" << endl;
+
+    const TArrayD &check = disp.GetParameters();
+    for (Int_t i=0; i<check.GetSize(); i++)
+        *fLog << check[i] << ",  ";
+    *fLog << endl;
+
+
+
+    //-------------------------------------------
+    // draw the plots
+    hdisp.Draw();
+
+    *fLog << "End of Disp optimization part" << endl;
+    *fLog << "======================================================" << endl;
+
+    return kTRUE;
+}
+
+
+// -----------------------------------------------------------------------
+//
+// Test the result of the Disp optimization on the test sample
+//
+
+Bool_t MFindDisp::TestParams()
+{
+    if (fMatrixTest == NULL)
+    {
+      *fLog << "MFindDisp::TestParams; test matrix is not defined... aborting"
+            << endl;
+      return kFALSE;
+    }
+
+    if (fMatrixTest->GetM().GetNrows() <= 0)
+    {
+        *fLog << "MFindDisp::TestParams; test matrix has no entries" 
+              << endl;
+        return kFALSE;
+    }
+
+    // -------------   TEST the Disp optimization    -------------------
+    //
+    *fLog << "Test the Disp optimization on the test sample" << endl;
+
+    // -----------------------------------------------------------------
+    // read optimum parameter values from file fFilenameParam
+    // into array 'dispPar'
+
+    TFile inparam(fFilenameParam);
+    MDisp dispin; 
+    dispin.Read("MDisp");
+    inparam.Close();
+
+    *fLog << "Optimum parameter values for Disp were read from file '";
+    *fLog << fFilenameParam << "' :" << endl;
+
+    const TArrayD &dispPar = dispin.GetParameters();
+    for (Int_t i=0; i<dispPar.GetSize(); i++)
+        *fLog << dispPar[i] << ",  ";
+    *fLog << endl;
+    //---------------------------
+
+
+    // -----------------------------------------------------------------
+    MParList  parlist2;
+    MTaskList tasklist2;
+
+    MDisp disp;
+    disp.SetParameters(dispPar);
+
+    MMatrixLoop loop(fMatrixTest);
+
+    // get the matrix pointer and mapping from MDispCalc object to set the same in MHDisp object
+    TArrayI map;
+    MHMatrix *tmpmatrix = NULL;
+    tmpmatrix = fDispCalcTest->GetMatrixMap(map);
+
+    // attention: argument of MHDisp is the name of MImageParDisp container, that should
+    // be the same than the name given to it when creating MDispCalc object at the MFindDisp
+    // constructor:  fDispCalcTrain = new MDispCalc("DispTest");
+    // fill the plots for Disp and sum up the Chi^2 contributions
+    MHDisp hdisp1("DispTest");
+    hdisp1.SetName("MHDispCorr");
+    hdisp1.SetSelectedPos(1);
+    hdisp1.SetMatrixMap(tmpmatrix,map);
+    MFillH filldispplots1("MHDispCorr[MHDisp]", "");
+    
+    MHDisp hdisp2("DispTest");
+    hdisp2.SetName("MHDispWrong");
+    hdisp2.SetSelectedPos(2);
+    hdisp2.SetMatrixMap(tmpmatrix,map);
+    MFillH filldispplots2("MHDispWrong[MHDisp]", "");
+    
+    MHDisp hdisp3("DispTest");
+    hdisp3.SetName("MHDispM3Long");
+    hdisp3.SetSelectedPos(3);
+    hdisp3.SetMatrixMap(tmpmatrix,map);
+    MFillH filldispplots3("MHDispM3Long[MHDisp]", "");
+    
+    MHDisp hdisp4("DispTest");
+    hdisp4.SetName("MHDispAsym");
+    hdisp4.SetSelectedPos(4);
+    hdisp4.SetMatrixMap(tmpmatrix,map);
+    MFillH filldispplots4("MHDispAsym[MHDisp]", "");
+    
+
+    //******************************
+    // entries in MParList
+
+    parlist2.AddToList(&tasklist2);
+    parlist2.AddToList(&disp);
+    parlist2.AddToList(&hdisp1);
+    parlist2.AddToList(&hdisp2);
+    parlist2.AddToList(&hdisp3);
+    parlist2.AddToList(&hdisp4);
+    parlist2.AddToList(fCam);
+    parlist2.AddToList(fMatrixTest);
+
+    //******************************
+    // entries in MTaskList
+
+    tasklist2.AddToList(&loop);
+    tasklist2.AddToList(fDispCalcTest);
+    tasklist2.AddToList(&filldispplots1);
+    tasklist2.AddToList(&filldispplots2);
+    tasklist2.AddToList(&filldispplots3);
+    tasklist2.AddToList(&filldispplots4);
+
+    //******************************
+
+    MProgressBar bar2;
+    MEvtLoop evtloop2;
+    evtloop2.SetParList(&parlist2);
+    evtloop2.SetName("EvtLoopTestParams");
+    evtloop2.SetProgressBar(&bar2);
+    Int_t maxevents2 = -1;
+    if (!evtloop2.Eventloop(maxevents2))
+        return kFALSE;
+
+    tasklist2.PrintStatistics(0, kTRUE);
+
+
+    //-------------------------------------------
+    // draw the plots
+
+    hdisp1.Draw();
+    hdisp2.Draw();
+    hdisp3.Draw();
+    hdisp4.Draw();
+
+    //-------------------------------------------
+
+
+    *fLog << "Test of Disp otimization finished" << endl;
+    *fLog << "======================================================" << endl;
+
+    return kTRUE;
+}
+
+
+
+
+
+
+
+
Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MFindDisp.h
===================================================================
--- trunk/MagicSoft/Mars/mtemp/mifae/library/MFindDisp.h	(revision 5671)
+++ trunk/MagicSoft/Mars/mtemp/mifae/library/MFindDisp.h	(revision 5671)
@@ -0,0 +1,126 @@
+#ifndef MARS_MFindDisp
+#define MARS_MFindDisp
+
+#ifndef MARS_MParContainer
+#include "MParContainer.h"
+#endif
+
+#ifndef ROOT_TArrayD
+#include <TArrayD.h>
+#endif
+
+#ifndef ROOT_TArrayI
+#include <TArrayI.h>
+#endif
+
+class MHMatrix;
+class MDispCalc;
+class MFDisp;
+class MFilter;
+class MH3;
+class MGeomCam;
+class MEvtLoop;
+
+class MFindDisp : public MParContainer
+{
+private:
+
+  TString   fFilenameParam;  // filename to store optimum Disp parameters
+
+  TString   fFilenameTrain;
+  TString   fFilenameTest;
+
+  Int_t     fHowManyTrain;   // number of events for training
+  Int_t     fHowManyTest;    // number of events for testing
+
+  MHMatrix  *fMatrixTrain;   // training matrix
+  MHMatrix  *fMatrixTest;    // testing matrix
+
+  MDispCalc *fDispCalcTrain;
+  MDispCalc *fDispCalcTest;
+
+  MFDisp    *fDispFilter;    // filter to select an events sample
+
+  MGeomCam  *fCam;
+
+  MEvtLoop  *fObjectFit;
+
+  //--------------------------------------------
+  // To comunicate with MINUIT
+  //--------------------------------------------
+  // attention : dimensions must agree with those in 
+  //             MMinuitInterface::CallMinuit()
+  //char  fParName [80][100];
+  TArrayD   fVinit;
+  TArrayD   fStep;
+  TArrayD   fLimlo;
+  TArrayD   fLimup;
+  TArrayI   fFix;
+
+  UInt_t    fNpar;
+
+  TString   fMethod;
+
+  Double_t  fMin,   fEdm,   fErrdef;
+  Int_t     fNpari, fNparx, fIstat;
+  Int_t     fErrMinimize;
+  //--------------------------------------------
+
+public:
+
+  MFindDisp(MFDisp *fdisp=NULL, const char *name=NULL, const char *title=NULL);
+  ~MFindDisp();
+
+  void SetFilenameParam(const TString &name)    {fFilenameParam  = name;}
+
+  void SetFilenameTraining(const TString &name, const Int_t howmany) 
+      {fFilenameTrain = name;  fHowManyTrain = howmany; }
+
+  void SetFilenameTest(const TString &name, const Int_t howmany)     
+      {fFilenameTest     = name;  fHowManyTest  = howmany; }
+
+  Bool_t DefineTrainMatrix(const TString &name, MH3 &href,
+                           const Int_t howmany, const TString &filetrain,
+			   Int_t iflag=0); 
+
+  Bool_t DefineTestMatrix(const TString &name, MH3 &href,
+                          const Int_t howmany, const TString &filetest,
+			  Int_t iflag=0); 
+
+  Bool_t DefineTrainTestMatrix(const TString &name, MH3 &href,
+			       const Int_t howmanytrain, const Int_t howmanytest, 
+			       const TString &filetrain, const TString &filetest,
+			       Int_t iflag=0); 
+
+  MHMatrix *GetMatrixTrain() { return fMatrixTrain; }
+  MHMatrix *GetMatrixTest()  { return fMatrixTest;  }
+
+  Bool_t ReadMatrix( const TString &filetrain, const TString &filetest);
+
+  Bool_t FindParams(TString parSCinit, TArrayD &params, TArrayD &steps);
+
+  Bool_t TestParams();
+
+  ClassDef(MFindDisp, 1) // Class for optimizing Disp
+};
+
+#endif
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MHDisp.cc
===================================================================
--- trunk/MagicSoft/Mars/mtemp/mifae/library/MHDisp.cc	(revision 5671)
+++ trunk/MagicSoft/Mars/mtemp/mifae/library/MHDisp.cc	(revision 5671)
@@ -0,0 +1,664 @@
+/* ======================================================================== *\
+!
+! *
+! * 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
+!
+!
+\* ======================================================================== */
+
+/////////////////////////////////////////////////////////////////////////////
+//                                                                         //
+//   MHDisp                                                                //
+//                                                                         //
+//   container holding the histograms for Disp                             //
+//   also computes the Chi^2 of the Disp optimization                      // 
+//                                                                         //
+/////////////////////////////////////////////////////////////////////////////
+#include "MHDisp.h"
+
+#include <math.h>
+
+#include <TH1.h>
+#include <TH2.h>
+#include <TProfile.h>
+#include <TArrayI.h>
+#include <TPad.h>
+#include <TCanvas.h>
+#include <TStyle.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 "MHMatrix.h"
+
+#include "MParList.h"
+
+#include "MLog.h"
+#include "MLogManip.h"
+
+ClassImp(MHDisp);
+
+using namespace std;
+
+// --------------------------------------------------------------------------
+//
+// Constructor 
+//
+MHDisp::MHDisp(const char *imagepardispname,
+	       const char *name, const char *title)
+  : fImageParDispName(imagepardispname)
+{
+    fName  = name  ? name  : "MHDisp";
+    fTitle = title ? title : "Histograms for Disp";
+
+    fSelectedPos = 1;  // default MHDisp flag (selects Correct Disp srcpos)
+
+    fMatrix = NULL;
+
+    //--------------------------------------------------
+    // creating the Disp related histograms
+
+    fHistEnergy   = new TH1F("fHistEnergy", 
+	 "log10(Energy)", 50, 1., 3.);
+    fHistEnergy->SetDirectory(NULL);
+    fHistEnergy->UseCurrentStyle();
+    fHistEnergy->SetXTitle("log10(Energy) [GeV]");
+    fHistEnergy->SetYTitle("# events");
+
+    fHistSize   = new TH1F("fHistSize", 
+	 "log10(Size)", 50, 2., 4.);
+    fHistSize->SetDirectory(NULL);
+    fHistSize->UseCurrentStyle();
+    fHistSize->SetXTitle("log10(Size) [#phot]");
+    fHistSize->SetYTitle("# events");
+
+    fHistcosZA   = new TH1F("fHistcosZA", 
+	 "cos(Zenith Angle)", 10, 0., 1.);
+    fHistcosZA->SetDirectory(NULL);
+    fHistcosZA->UseCurrentStyle();
+    fHistcosZA->SetXTitle("cos(Theta)");
+    fHistcosZA->SetYTitle("# events");
+
+    fSkymapXY = new TH2F("fSkymapXY", 
+         "Disp estimated source positions Skymap", 71, -2., 2., 71, -2., 2.);
+    fSkymapXY->SetDirectory(NULL);
+    fSkymapXY->UseCurrentStyle();
+    fSkymapXY->SetXTitle("X Disp source position [deg]");
+    fSkymapXY->SetYTitle("Y Disp source position [deg]");
+
+    fHistChi2   = new TH1F("fHistChi2"  , 
+         "Distance^2 between Disp and real srcpos", 100, 0., 2.);
+    fHistChi2->SetDirectory(NULL);
+    fHistChi2->UseCurrentStyle();
+    fHistChi2->SetXTitle("Chi^2 = distance^2 Disp to real srcpos [deg^2]");
+    fHistChi2->SetYTitle("# events");
+    
+    fHistDuDv   = new TH2F("fHistDuDv", 
+	 "Du vs. Dv (distances between Disp and real srcpos)", 
+	 100, -2., 2., 100, -2., 2.);
+    fHistDuDv->SetDirectory(NULL);
+    fHistDuDv->UseCurrentStyle();
+    fHistDuDv->SetXTitle("Dv = transveral distance [deg]");
+    fHistDuDv->SetYTitle("Du = longitudinal distance [deg]");
+
+    fHistChi2Energy   = new TH2F("fHistChi2Energy", 
+	 "Chi^2 vs. Energy", 50, 1., 3., 100, 0., 2.);
+    fHistChi2Energy->SetDirectory(NULL);
+    fHistChi2Energy->UseCurrentStyle();
+    fHistChi2Energy->SetXTitle("log10(Energy) [GeV]");
+    fHistChi2Energy->SetYTitle("Chi^2 = distance^2 Disp to real srcpos [deg^2]");
+
+    fHistChi2Size   = new TH2F("fHistChi2Size", 
+	 "Chi^2 vs. Size", 50, 2., 4., 100, 0., 2.);
+    fHistChi2Size->SetDirectory(NULL);
+    fHistChi2Size->UseCurrentStyle();
+    fHistChi2Size->SetXTitle("log10(Size) [#phot]");
+    fHistChi2Size->SetYTitle("Chi^2 = distance^2 Disp to real srcpos [deg^2]");
+
+    fHistDuEnergy   = new TH2F("fHistDuEnergy", 
+	 "Du vs. Energy", 50, 1., 3., 100, -2., 2.);
+    fHistDuEnergy->SetDirectory(NULL);
+    fHistDuEnergy->UseCurrentStyle();
+    fHistDuEnergy->SetXTitle("log10(Energy) [GeV]");
+    fHistDuEnergy->SetYTitle("Du = longitudinal distance Disp to real srcpos [deg]");
+
+    fHistDuSize   = new TH2F("fHistDuSize", 
+	 "Du vs. Size", 50, 2., 4., 100, -2., 2.);
+    fHistDuSize->SetDirectory(NULL);
+    fHistDuSize->UseCurrentStyle();
+    fHistDuSize->SetXTitle("log10(Size) [#phot]");
+    fHistDuSize->SetYTitle("Du = longitudinal distance Disp to real srcpos [deg]");
+
+    fHistDvEnergy   = new TH2F("fHistDvEnergy", 
+	 "Dv vs. Energy", 50, 1., 3., 100, -2., 2.);
+    fHistDvEnergy->SetDirectory(NULL);
+    fHistDvEnergy->UseCurrentStyle();
+    fHistDvEnergy->SetXTitle("log10(Energy) [GeV]");
+    fHistDvEnergy->SetYTitle("Dv = transveral distance Disp to real srcpos [deg]");
+
+    fHistDvSize   = new TH2F("fHistDvSize", 
+	 "Dv vs. Size", 50, 2., 4., 100, -2., 2.);
+    fHistDvSize->SetDirectory(NULL);
+    fHistDvSize->UseCurrentStyle();
+    fHistDvSize->SetXTitle("log10(Size) [#phot]");
+    fHistDvSize->SetYTitle("Dv = transveral distance Disp to real srcpos [deg]");
+
+    fEvCorrAssign   = new TProfile("fEvCorrAssign", 
+	 "Fraction events srcpos well assign vs. Size", 50, 2., 4., 0., 1.);
+    fEvCorrAssign->SetDirectory(NULL);
+    fEvCorrAssign->UseCurrentStyle();
+    fEvCorrAssign->SetXTitle("log10(Size) [#phot]");
+    fEvCorrAssign->SetYTitle("Fraction of events with srcpos well assign");
+}
+
+
+// --------------------------------------------------------------------------
+//
+// Destructor 
+//
+MHDisp::~MHDisp()
+{
+  delete fHistEnergy;
+  delete fHistSize;
+  delete fHistcosZA;
+  delete fSkymapXY;
+  delete fHistChi2;
+  delete fHistDuDv;
+  delete fHistChi2Energy;
+  delete fHistChi2Size;
+  delete fHistDuEnergy;
+  delete fHistDuSize;
+  delete fHistDvEnergy;
+  delete fHistDvSize;
+  delete fEvCorrAssign;
+}
+
+// --------------------------------------------------------------------------
+//
+// Set the pointers to the containers 
+// 
+//
+Bool_t MHDisp::SetupFill(const MParList *pList)
+{
+    // reset all histograms and Chi^2 computing variables
+    // before each new eventloop
+    fNumEv = 0;
+    fSumChi2  = 0.;
+    fChi2  = 0.;
+
+    fHistEnergy->Reset();
+    fHistSize->Reset();
+    fHistcosZA->Reset();
+    fSkymapXY->Reset();
+    fHistChi2->Reset();
+    fHistDuDv->Reset();
+    fHistChi2Energy->Reset();
+    fHistChi2Size->Reset();
+    fHistDuEnergy->Reset();
+    fHistDuSize->Reset();
+    fHistDvEnergy->Reset();
+    fHistDvSize->Reset();
+    fEvCorrAssign->Reset();
+
+
+    // 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->FindObject(fImageParDispName,
+						      "MImageParDisp");
+    if (!fImageParDisp)
+    {
+        *fLog << err << fImageParDispName 
+              << " [MImageParDisp] 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... aborting." << 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;
+}
+
+
+// --------------------------------------------------------------------------
+//
+// Set which selection algorithm for the Disp estimated source position 
+// we want to follow when filling the histograms:
+//  * iflag == 1 --> Correct source position, the closest to the real srcpos
+//                   (only applicable to MC or well known source real data)
+//  * iflag == 2 --> Wrong source position, the furthest to the real srcpos
+//                   (same applicability than case 1)
+//  * iflag == 3 --> Source position assigned as M3Long sign indicates
+//  * iflag == 4 --> Source position assigned as Asym sign indicates
+//
+void MHDisp::SetSelectedPos(Int_t iflag)
+{ 
+    fSelectedPos = iflag; 
+}    
+
+
+// --------------------------------------------------------------------------
+//
+// Sets the Matrix and the array of mapped values for each Matrix column
+// (see MDispCalc)
+//
+void MHDisp::SetMatrixMap(MHMatrix *matrix, TArrayI &map)
+{
+    fMatrix = matrix;
+    fMap = map;
+}
+
+
+// --------------------------------------------------------------------------
+//
+// Returns a mapped value from the Matrix
+//
+Double_t MHDisp::GetVal(Int_t i) const
+{
+    Double_t val = (*fMatrix)[fMap[i]];
+
+    //*fLog << "MHDisp::GetVal; i, fMatrix, fMap, val = "
+    //    << i << ",  " << fMatrix << ",  " << fMap[i] << ",  " << val << endl;
+
+    return val;
+}
+
+
+// --------------------------------------------------------------------------
+//
+// Fill the histograms and sum up the Chi^2 outcome of each processed event
+//
+Bool_t MHDisp::Fill(const MParContainer *par, const Stat_t w)
+{
+    Double_t energy = 0.;
+    Double_t impact = 0.;
+    Double_t xmax = 0.;
+    Double_t theta = 0.;
+    Double_t phi = 0.;
+
+    if ( fMatrix || (!fMatrix && fMcEvt) )
+    {  
+      energy   = fMatrix ? GetVal(13)  : fMcEvt->GetEnergy();
+      impact   = fMatrix ? GetVal(14)  : fMcEvt->GetImpact();
+      xmax     = fMatrix ? GetVal(15)  : fMcEvt->GetLongitmax();
+      
+      if (fPointing)
+      {
+	theta  = fMatrix ? GetVal(16)  : fPointing->GetZd();
+	phi    = fMatrix ? GetVal(17)  : fPointing->GetAz();
+      }
+      else
+      {
+	theta  = fMatrix ? GetVal(16)  : fMcEvt->GetTelescopeTheta()*kRad2Deg;
+	phi    = fMatrix ? GetVal(17)  : fMcEvt->GetTelescopePhi()*kRad2Deg;
+      }
+    }
+    else
+      *fLog << "No MMcEvt container available (no Energy,ImpactPar or Xmax)" << endl;
+
+    Double_t xsrcpos0   = fMatrix ? GetVal(0)   : fSrcPos->GetX();
+    Double_t ysrcpos0   = fMatrix ? GetVal(1)   : fSrcPos->GetY();
+    Double_t xmean0     = fMatrix ? GetVal(2)   : fHil->GetMeanX();  
+    Double_t ymean0     = fMatrix ? GetVal(3)   : fHil->GetMeanY();
+    Double_t delta      = fMatrix ? GetVal(4)   : fHil->GetDelta();
+    
+    Double_t size       = fMatrix ? GetVal(5)   : fHil->GetSize();
+    //    Double_t width0     = fMatrix ? GetVal(6)   : fHil->GetWidth();
+    //    Double_t length0    = fMatrix ? GetVal(7)   : fHil->GetLength();
+    
+    //    Double_t conc       = fMatrix ? GetVal(8)   : fNewPar->GetConc();
+    //    Double_t leakage1   = fMatrix ? GetVal(9)   : fNewPar->GetLeakage1();
+    //    Double_t leakage2   = fMatrix ? GetVal(10)  : fNewPar->GetLeakage2();
+    
+    Double_t m3long     = fMatrix ? GetVal(11)  : fHilExt->GetM3Long();
+    Double_t asym       = fMatrix ? GetVal(12)  : fHilExt->GetAsym();
+
+    //------------------------------------------
+    // convert to deg
+    //    Double_t width    = width0  * fMm2Deg;
+    //    Double_t length   = length0 * fMm2Deg;
+    Double_t xsrcpos  = xsrcpos0 * fMm2Deg;
+    Double_t ysrcpos  = ysrcpos0 * fMm2Deg;
+    Double_t xmean    = xmean0 * fMm2Deg;
+    Double_t ymean    = ymean0 * fMm2Deg;
+    
+    // calculate cosinus of the angle between d and a vectors
+    Double_t a = (xmean-xsrcpos)*cos(delta) + (ymean-ysrcpos)*sin(delta);
+    Double_t b = sqrt( (xmean-xsrcpos)*(xmean-xsrcpos) + (ymean-ysrcpos)*(ymean-ysrcpos) );
+    Double_t cosda = a/b;
+    
+    // sign of cosda
+    Int_t s0 = cosda>0 ? 1 : -1;   
+    
+    // get the sign of M3Long and Asym variables
+    Int_t sm3 = m3long>0 ? 1 : -1;
+    Int_t sa  =   asym>0 ? 1 : -1;
+    
+    // set the sign "s" that will select one Disp source position for each
+    // shower, according to each of the possible algorithms for solution selection:
+    //   SelectedPos = 1  means choose the right source position
+    //                 2                   wrong
+    //                 3               the position according to M3Long
+    //                 4               the position according to Asym
+    Int_t s = s0;
+    if (fSelectedPos == 1)    
+      s = s0;  
+    else if (fSelectedPos == 2)
+      s = -s0;
+    else if (fSelectedPos == 3)
+      s = sm3;
+    else if (fSelectedPos == 4)
+      s = sa;
+    else
+      *fLog << "Illegal value for Disp srcpos selection algorithm: " 
+	    << " fSelectedPos = " << fSelectedPos << endl;
+    
+    // count the number of events where source position has been correctly assigned
+    if (s == s0)
+      fEvCorrAssign->Fill(log10(size), 1.);
+    else
+      fEvCorrAssign->Fill(log10(size), 0.);
+    
+    // get estimated Disp value
+    Double_t disp = fImageParDisp->GetDisp();
+    
+    //------------------------------------------
+    // Disp estimated source position
+    Double_t xdisp = xmean - s*cos(delta)*disp;
+    Double_t ydisp = ymean - s*sin(delta)*disp;
+    fSkymapXY->Fill(xdisp,ydisp);
+    
+    // Distance between estimated Disp and real source position
+    Double_t d2 = (xdisp-xsrcpos)*(xdisp-xsrcpos) +  (ydisp-ysrcpos)*(ydisp-ysrcpos); 
+    fHistChi2->Fill(d2);
+    
+    // Longitudinal and transversal distances between Disp and real source positon
+    Double_t Du = -s*( (xdisp-xsrcpos)*cos(delta) + (ydisp-ysrcpos)*sin(delta));
+    Double_t Dv = -s*(-(xdisp-xsrcpos)*sin(delta) + (ydisp-ysrcpos)*cos(delta));
+    fHistDuDv->Fill(Dv,Du);
+    
+    // Fill Energy, Size and ZA distributions
+    fHistEnergy->Fill(log10(energy));
+    fHistSize->Fill(log10(size));
+    fHistcosZA->Fill(cos(theta/kRad2Deg));
+    
+    // to check the size dependence of the optimization
+    fHistChi2Energy->Fill(log10(energy),d2);
+    fHistChi2Size->Fill(log10(size),d2);
+    fHistDuEnergy->Fill(log10(energy),Du);
+    fHistDuSize->Fill(log10(size),Du);
+    fHistDvEnergy->Fill(log10(energy),Dv);
+    fHistDvSize->Fill(log10(size),Dv);
+    
+    // variables for Chi^2 calculation (= d^2)
+    fNumEv += 1;
+    fSumChi2 += d2;  
+    
+    return kTRUE;
+}
+
+
+// --------------------------------------------------------------------------
+//
+// Calculates the final Chi^2 of the Disp optimization
+//
+Bool_t MHDisp::Finalize()
+{
+    fChi2 = fSumChi2/fNumEv;
+    *fLog << endl;
+    *fLog << "MHDisp::Finalize: SumChi2, NumEv = " << fSumChi2 << ", " << fNumEv << endl;
+    *fLog << "MHDisp::Finalize: Final Chi2 = " << fChi2 << endl;
+    
+    fChi2 = fHistChi2->GetMean();
+    *fLog << "MHDisp::Finalize: Final Chi2 = " << fChi2 << endl;
+
+    SetReadyToSave();
+
+    return kTRUE;
+}
+
+
+// --------------------------------------------------------------------------
+//
+// Creates a new canvas and draws the Disp related histograms into it.
+// Be careful: The histograms belongs to this object and won't get deleted
+// together with the canvas.
+//
+void MHDisp::Draw(Option_t *opt)
+{
+    gStyle->SetPalette(1);
+
+    TString title = GetName();
+    title += ": ";
+    title += GetTitle();
+
+    TCanvas *pad = new TCanvas(GetName(),title,0,0,900,1500);    
+    pad->SetBorderMode(0);
+    pad->Divide(3,5);
+
+    pad->cd(1);
+    gPad->SetBorderMode(0);
+    fHistcosZA->SetTitleOffset(1.2,"Y");
+    fHistcosZA->DrawClone(opt);
+    fHistcosZA->SetBit(kCanDelete);
+    gPad->Modified();
+
+    pad->cd(2);
+    gPad->SetBorderMode(0);
+    fHistEnergy->SetTitleOffset(1.2,"Y");
+    fHistEnergy->DrawClone(opt);
+    fHistEnergy->SetBit(kCanDelete);
+    gPad->Modified();
+
+    pad->cd(3);
+    gPad->SetBorderMode(0);
+    fHistSize->SetTitleOffset(1.2,"Y");
+    fHistSize->DrawClone(opt);
+    fHistSize->SetBit(kCanDelete);
+    gPad->Modified();
+
+    pad->cd(4);
+    gPad->SetBorderMode(0);
+    fHistChi2->SetTitleOffset(1.2,"Y");
+    fHistChi2->DrawClone(opt);
+    fHistChi2->SetBit(kCanDelete);
+    gPad->Modified();
+
+    TProfile *profChi2Energy = fHistChi2Energy->ProfileX("Chi^2 vs. Energy",0,99999,"s");
+    profChi2Energy->SetXTitle("log10(Energy) [GeV]");
+    profChi2Energy->SetYTitle("Chi^2 = distance^2 Disp to real srcpos [deg^2]");
+    pad->cd(5);
+    gPad->SetBorderMode(0);
+    profChi2Energy->SetTitleOffset(1.2,"Y");
+    profChi2Energy->DrawClone(opt);
+    profChi2Energy->SetBit(kCanDelete);
+    gPad->Modified();
+
+    TProfile *profChi2Size = fHistChi2Size->ProfileX("Chi^2 vs. Size",0,99999,"s");
+    profChi2Size->SetXTitle("log10(Size) [#phot]");
+    profChi2Size->SetYTitle("Chi^2 = distance^2 Disp to real srcpos [deg^2]");
+    pad->cd(6);
+    gPad->SetBorderMode(0);
+    profChi2Size->SetTitleOffset(1.2,"Y");
+    profChi2Size->DrawClone(opt);
+    profChi2Size->SetBit(kCanDelete);
+    gPad->Modified();
+
+    pad->cd(7);
+    gPad->SetBorderMode(0);
+    fSkymapXY->SetTitleOffset(1.2,"Y");
+    fSkymapXY->DrawClone("COLZ");
+    fSkymapXY->SetBit(kCanDelete);
+    gPad->Modified();
+
+    pad->cd(8);
+    gPad->SetBorderMode(0);
+    fEvCorrAssign->SetTitleOffset(1.2,"Y");
+    fEvCorrAssign->DrawClone(opt);
+    fEvCorrAssign->SetBit(kCanDelete);
+    gPad->Modified();
+
+    pad->cd(9);
+    gPad->SetBorderMode(0);
+    fHistDuDv->SetTitleOffset(1.2,"Y");
+    fHistDuDv->DrawClone("COLZ");
+    fHistDuDv->SetBit(kCanDelete);
+    gPad->Modified();
+
+    TH1F *histDu = (TH1F*)fHistDuDv->ProjectionY("histDu");
+    histDu->SetTitle("Longitudinal distance Du");
+    histDu->SetXTitle("Du = longitudinal distance [deg]");
+    pad->cd(10);
+    gPad->SetBorderMode(0);
+    histDu->SetTitleOffset(1.2,"Y");
+    histDu->DrawClone(opt);
+    histDu->SetBit(kCanDelete);
+    gPad->Modified();
+
+    TProfile *profDuEnergy = fHistDuEnergy->ProfileX("Du vs. Energy",0,99999,"s");
+    profDuEnergy->SetXTitle("log10(Energy) [GeV]");
+    profDuEnergy->SetYTitle("Du = longitudinal distance [deg]");
+    pad->cd(11);
+    gPad->SetBorderMode(0);
+    profDuEnergy->SetTitleOffset(1.2,"Y");
+    profDuEnergy->DrawClone(opt);
+    profDuEnergy->SetBit(kCanDelete);
+    gPad->Modified();
+
+    TProfile *profDuSize = fHistDuSize->ProfileX("Du vs. Size",0,99999,"s");
+    profDuSize->SetXTitle("log10(Size) [#phot]");
+    profDuSize->SetYTitle("Du = longitudinal distance [deg]");
+    pad->cd(12);
+    gPad->SetBorderMode(0);
+    profDuSize->SetTitleOffset(1.2,"Y");
+    profDuSize->DrawClone(opt);
+    profDuSize->SetBit(kCanDelete);
+    gPad->Modified();
+
+    TH1F *histDv = (TH1F*)fHistDuDv->ProjectionX("histDv");
+    histDv->SetTitle("Transversal distance Dv");
+    histDv->SetXTitle("Dv = transversal distance [deg]");
+    pad->cd(13);
+    gPad->SetBorderMode(0);
+    histDv->SetTitleOffset(1.2,"Y");
+    histDv->DrawClone(opt);
+    histDv->SetBit(kCanDelete);
+    gPad->Modified();
+
+    TProfile *profDvEnergy = fHistDvEnergy->ProfileX("Dv vs. Energy",0,99999,"s");
+    profDvEnergy->SetXTitle("log10(Energy) [GeV]");
+    profDvEnergy->SetYTitle("Dv = transversal distance [deg]");
+    pad->cd(14);
+    gPad->SetBorderMode(0);
+    profDvEnergy->SetTitleOffset(1.2,"Y");
+    profDvEnergy->DrawClone(opt);
+    profDvEnergy->SetBit(kCanDelete);
+    gPad->Modified();
+
+    TProfile *profDvSize = fHistDvSize->ProfileX("Dv vs. Size",0,99999,"s");
+    profDvSize->SetXTitle("log10(Size) [#phot]");
+    profDvSize->SetYTitle("Dv = transversal distance [deg]");
+    pad->cd(15);
+    gPad->SetBorderMode(0);
+    profDvSize->SetTitleOffset(1.2,"Y");
+    profDvSize->DrawClone(opt);
+    profDvSize->SetBit(kCanDelete);
+    gPad->Modified();
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MHDisp.h
===================================================================
--- trunk/MagicSoft/Mars/mtemp/mifae/library/MHDisp.h	(revision 5671)
+++ trunk/MagicSoft/Mars/mtemp/mifae/library/MHDisp.h	(revision 5671)
@@ -0,0 +1,107 @@
+#ifndef MARS_MHDisp
+#define MARS_MHDisp
+
+#ifndef MARS_MH
+#include "MH.h"
+#endif
+
+#ifndef ROOT_TArrayI
+#include <TArrayI.h>
+#endif
+
+class TH1F;
+class TH2F;
+class TProfile;
+class MImageParDisp;
+class MSrcPosCam;
+class MHillas;
+class MHillasExt;
+class MNewImagePar;
+class MMcEvt;
+class MPointingPos;
+class MHMatrix;
+
+class MHDisp : public MH
+{
+private:
+
+    MImageParDisp *fImageParDisp;  // container with the estimated Disp
+    MSrcPosCam    *fSrcPos;
+    MHillas       *fHil;
+    MHillasExt    *fHilExt;
+    MNewImagePar  *fNewPar;
+    MMcEvt        *fMcEvt;
+    MPointingPos  *fPointing;
+
+    TString fImageParDispName;
+
+    Int_t fSelectedPos;      // flag to select which Disp estimated position
+                             // we want to fill the histograms (see Set function)
+
+    TH1F *fHistEnergy;       // Energy distribution of events
+    TH1F *fHistSize;         // Size distribution of events
+    TH1F *fHistcosZA;        // cosinus Zenith angle distribution of events
+    TH2F *fSkymapXY;         // 2D histogram for Disp estimated source positions
+    TH1F *fHistChi2;         // Histogram of the event Chi^2 (= d^2 = distance^2 
+                             // between Disp estimated and true source position)
+    TH2F *fHistDuDv;         // Distribution of longitudinal (Du) and transversal 
+                             // (Dv) distances between Disp and true source position    
+    TH2F *fHistChi2Energy;   // Chi^2 (= d^2) vs. Energy
+    TH2F *fHistChi2Size;     // Chi^2 (= d^2) vs. Size
+    TH2F *fHistDuEnergy;     // Du vs. Energy
+    TH2F *fHistDuSize;       // Du vs. Size
+    TH2F *fHistDvEnergy;     // Dv vs. Energy
+    TH2F *fHistDvSize;       // Dv vs. Size
+    TProfile *fEvCorrAssign; // % events with source position well assign vs. Size 
+
+
+    Double_t fMm2Deg;        // conversion factor from mm to deg
+
+    MHMatrix *fMatrix;       // matrix storing variables needed for Disp optimization
+    TArrayI fMap;            // array storing the matrix mapping of columns defined
+                             // in MDispCalc::InitMapping
+
+    Double_t GetVal(Int_t i) const;
+
+    Int_t fNumEv;            // total number of events
+    Double_t fSumChi2;       // current sum of Chi^2      
+    Double_t fChi2;          // final Chi^2 of the Disp optimization
+
+public:
+
+    MHDisp(const char *imagepardispname = "MImageParDisp",
+	   const char *name=NULL, const char *title=NULL);
+    ~MHDisp();
+
+    Bool_t SetupFill(const MParList *plist);
+    Bool_t Fill(const MParContainer *par, const Stat_t w=1);
+    Bool_t Finalize();  
+
+    void SetSelectedPos(Int_t iflag);
+
+    void SetMatrixMap(MHMatrix *matrix, TArrayI &map);    
+
+    Double_t GetChi2()     { return fChi2; }
+
+    TH1F *GetHistEnergy()       { return fHistEnergy; } 
+    TH1F *GetHistSize()         { return fHistSize; } 
+    TH1F *GetHistcosZA()        { return fHistcosZA; } 
+    TH2F *GetSkymapXY()         { return fSkymapXY; }
+    TH1F *GetHistChi2()         { return fHistChi2; }
+    TH2F *GetHistDuDv()         { return fHistDuDv; }
+    TH2F *GetHistChi2Energy()   { return fHistChi2Energy; }
+    TH2F *GetHistChi2Size()     { return fHistChi2Size; }
+    TH2F *GetHistDuEnergy()     { return fHistDuEnergy; }
+    TH2F *GetHistDuSize()       { return fHistDuSize; }
+    TH2F *GetHistDvEnergy()     { return fHistDvEnergy; }
+    TH2F *GetHistDvSize()       { return fHistDvSize; }
+    TProfile *GetEvCorrAssign() { return fEvCorrAssign; }
+
+    void Draw(Option_t *opt="");
+
+    ClassDef(MHDisp, 1) // Container holding the Histograms for Disp
+};
+
+#endif
+
+
