/* ======================================================================== *\
!
! *
! * 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 minimization parameter 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 "MParameters.h"

#include "MParList.h"

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

ClassImp(MHDisp);

using namespace std;

enum dispVar_t {kX,kY,kMeanX,kMeanY,kDelta,kSize,kM3Long,kAsym,
		kEnergy,kImpact,kLongitmax,kZd,kAz,kTotVar};    // enum variables for the
                                                                // matrix columns mapping

// --------------------------------------------------------------------------
//
// 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 source position solution)

    fMatrix = NULL;

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

    //--------------------------------------------------
    // 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]");

    fHistMinPar   = new TH1F("fHistMinPar"  , 
         "Distance^2 between Disp and real srcpos", 100, 0., 2.);
    fHistMinPar->SetDirectory(NULL);
    fHistMinPar->UseCurrentStyle();
    fHistMinPar->SetXTitle("Minimization parameter = d^2 Disp to real srcpos [deg^2]");
    fHistMinPar->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]");

    fHistMinParEnergy   = new TH2F("fHistMinParEnergy", 
	 "Minimization parameter vs. Energy", 50, 1., 3., 100, 0., 2.);
    fHistMinParEnergy->SetDirectory(NULL);
    fHistMinParEnergy->UseCurrentStyle();
    fHistMinParEnergy->SetXTitle("log10(Energy) [GeV]");
    fHistMinParEnergy->SetYTitle("Minimization parameter = d^2 Disp to real srcpos [deg^2]");

    fHistMinParSize   = new TH2F("fHistMinParSize", 
	 "Minimization parameter vs. Size", 50, 2., 4., 100, 0., 2.);
    fHistMinParSize->SetDirectory(NULL);
    fHistMinParSize->UseCurrentStyle();
    fHistMinParSize->SetXTitle("log10(Size) [#phot]");
    fHistMinParSize->SetYTitle("Minimization parameter = d^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->SetStats(0);
    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 fHistMinPar;
  delete fHistDuDv;
  delete fHistMinParEnergy;
  delete fHistMinParSize;
  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 Minimization parameter computing variables
    // before each new eventloop
    fNumEv = 0;
    fSumMinPar  = 0.;
    fMinPar = (MParameterD*)const_cast<MParList*>(pList)->FindCreateObj("MParameterD", "MinimizationParameter");
    if (!fMinPar)
    {
      *fLog << err << "MParameterD (MinimizationParameter) not found and could not be created... aborting." 
	    << endl;
      return kFALSE;
    }
    fMinPar->SetVal(0);

    fHistEnergy->Reset();
    fHistSize->Reset();
    fHistcosZA->Reset();
    fSkymapXY->Reset();
    fHistMinPar->Reset();
    fHistDuDv->Reset();
    fHistMinParEnergy->Reset();
    fHistMinParSize->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... 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... aborting." << endl;
	return kFALSE;
    }

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

    return kTRUE;
}


// --------------------------------------------------------------------------
//
// Set which selection algorithm for the Disp estimated source position 
// solutions 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
//
void MHDisp::SetMatrixMap(MHMatrix *matrix, TArrayI &map)
{
    fMatrix = matrix;
    fMap = map;
}


// --------------------------------------------------------------------------
//
// Returns the Matrix and the mapped value for each Matrix column
//
void MHDisp::GetMatrixMap(MHMatrix* &matrix, TArrayI &map)
{
    map = fMap;
    matrix = fMatrix;
}


// --------------------------------------------------------------------------
//
// 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), MHDisp::Fill will take the values from the matrix 
// instead of the containers.
//
void MHDisp::InitMapping(MHMatrix *mat)
{
    if (fMatrix)
      return;

    fMatrix = mat;

    fMap[kX]          = fMatrix->AddColumn("MSrcPosCam.fX");
    fMap[kY]          = fMatrix->AddColumn("MSrcPosCam.fY");
    fMap[kMeanX]      = fMatrix->AddColumn("MHillas.fMeanX");
    fMap[kMeanY]      = fMatrix->AddColumn("MHillas.fMeanY");
    fMap[kDelta]      = fMatrix->AddColumn("MHillas.fDelta");

    fMap[kSize]       = fMatrix->AddColumn("MHillas.fSize");
    
    fMap[kM3Long]     = fMatrix->AddColumn("MHillasExt.fM3Long");
    fMap[kAsym]       = fMatrix->AddColumn("MHillasExt.fAsym");
    
    fMap[kEnergy]     = fMatrix->AddColumn("MMcEvt.fEnergy");
    fMap[kImpact]     = fMatrix->AddColumn("MMcEvt.fImpact");
    fMap[kLongitmax]  = fMatrix->AddColumn("MMcEvt.fLongitmax");
    
    fMap[kZd]         = fMatrix->AddColumn("MPointingPos.fZd");
    fMap[kAz]         = fMatrix->AddColumn("MPointingPos.fAz");
}


// --------------------------------------------------------------------------
//
// 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 Minimization paramter 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.;

    if ( fMatrix || (!fMatrix && fMcEvt) )
    {  
      energy   = fMatrix ? GetVal(kEnergy)     : fMcEvt->GetEnergy();
      impact   = fMatrix ? GetVal(kImpact)     : fMcEvt->GetImpact();
      xmax     = fMatrix ? GetVal(kLongitmax)  : fMcEvt->GetLongitmax();      
    }

    Double_t theta      = fMatrix ? GetVal(kZd)  : fPointing->GetZd();
    //    Double_t phi        = fMatrix ? GetVal(kAz)  : fPointing->GetAz();

    Double_t xsrcpos0   = fMatrix ? GetVal(kX)       : fSrcPos->GetX();
    Double_t ysrcpos0   = fMatrix ? GetVal(kY)       : fSrcPos->GetY();
    Double_t xmean0     = fMatrix ? GetVal(kMeanX)   : fHil->GetMeanX();  
    Double_t ymean0     = fMatrix ? GetVal(kMeanY)   : fHil->GetMeanY();
    Double_t delta      = fMatrix ? GetVal(kDelta)   : fHil->GetDelta();
    
    Double_t size       = fMatrix ? GetVal(kSize)    : fHil->GetSize();
    
    Double_t m3long     = fMatrix ? GetVal(kM3Long)  : fHilExt->GetM3Long();
    Double_t asym       = fMatrix ? GetVal(kAsym)    : fHilExt->GetAsym();

    //------------------------------------------
    // convert to deg
    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); 
    fHistMinPar->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
    if (fMatrix || (!fMatrix && fMcEvt))
      fHistEnergy->Fill(log10(energy));
    fHistSize->Fill(log10(size));
    fHistcosZA->Fill(cos(theta/kRad2Deg));
    
    // to check the size and energy dependence of the optimization
    if (fMatrix || (!fMatrix && fMcEvt))
    {
      fHistMinParEnergy->Fill(log10(energy),d2);
      fHistDuEnergy->Fill(log10(energy),Du);
      fHistDvEnergy->Fill(log10(energy),Dv);
    }
    fHistMinParSize->Fill(log10(size),d2);
    fHistDuSize->Fill(log10(size),Du);
    fHistDvSize->Fill(log10(size),Dv);
    
    // variables for the Minimization parameter calculation (= d^2)
    fNumEv += 1;
    fSumMinPar += d2;  
    
    return kTRUE;
}


// --------------------------------------------------------------------------
//
// Calculates the final Minimization parameter of the Disp optimization
//
Bool_t MHDisp::Finalize()
{
    fMinPar->SetVal(fSumMinPar/fNumEv);
    *fLog << endl;
    *fLog << "MHDisp::Finalize: SumMinPar, NumEv = " << fSumMinPar << ", " << fNumEv << endl;
    *fLog << "MHDisp::Finalize: Final MinPar = " << fMinPar->GetVal() << endl;
    
    fMinPar->SetVal(fHistMinPar->GetMean());
    *fLog << "MHDisp::Finalize: Final MinPar = " << fMinPar->GetVal() << 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);
    fHistMinPar->SetTitleOffset(1.2,"Y");
    fHistMinPar->DrawClone(opt);
    fHistMinPar->SetBit(kCanDelete);
    gPad->Modified();

    TProfile *profMinParEnergy = fHistMinParEnergy->ProfileX("Minimization parameter vs. Energy",0,99999,"s");
    profMinParEnergy->SetXTitle("log10(Energy) [GeV]");
    profMinParEnergy->SetYTitle("Minimization parameter = d^2 Disp to real srcpos [deg^2]");
    pad->cd(5);
    gPad->SetBorderMode(0);
    profMinParEnergy->SetTitleOffset(1.2,"Y");
    profMinParEnergy->SetStats(0);
    profMinParEnergy->DrawClone(opt);
    profMinParEnergy->SetBit(kCanDelete);
    gPad->Modified();

    TProfile *profMinParSize = fHistMinParSize->ProfileX("Minimization parameter vs. Size",0,99999,"s");
    profMinParSize->SetXTitle("log10(Size) [#phot]");
    profMinParSize->SetYTitle("Minimization parameter = d^2 Disp to real srcpos [deg^2]");
    pad->cd(6);
    gPad->SetBorderMode(0);
    profMinParSize->SetTitleOffset(1.2,"Y");
    profMinParSize->SetStats(0);
    profMinParSize->DrawClone(opt);
    profMinParSize->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->SetStats(0);
    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->SetStats(0);
    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->SetStats(0);
    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->SetStats(0);
    profDvSize->DrawClone(opt);
    profDvSize->SetBit(kCanDelete);
    gPad->Modified();
}