source: trunk/MagicSoft/Mars/mhist/MHFalseSource.cc@ 3552

/* ======================================================================== *\
!
! *
! * 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): Thomas Bretz, 3/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
// MHFalseSource
//
// Create a false source plot. For the Binning in x,y the object MBinning
// "BinningFalseSource" is taken from the paremeter list.
//
// The binning in alpha is currently fixed as 18bins from 0 to 90deg.
//
// If MTime, MObservatory and MPointingPos is available in the paremeter
// list each image is derotated.
//
// MHFalseSource fills a 3D histogram with alpha distribution for
// each (derotated) x and y position.
//
// Only a radius of 1.2deg around the camera center is taken into account.
//
// The displayed histogram is the projection of alpha<fAlphaCut into
// the x,y plain and the projection of alpha>90-fAlphaCut
//
// By using the context menu (find it in a small region between the single
// pads) you can change the AlphaCut 'online'
//
// Here is a slightly simplified version of the algorithm:
// ------------------------------------------------------
//    MHillas hil; // Taken as second argument in MFillH
//
//    MSrcPosCam src;
//    MHillasSrc hsrc;
//    hsrc.SetSrcPos(&src);
//
//    for (int ix=0; ix<nx; ix++)
//        for (int iy=0; iy<ny; iy++)
//        {
//            TVector2 v(cx[ix], cy[iy]); //cx center of x-bin ix
//            if (rho!=0)                 //cy center of y-bin iy
//                v.Rotate(-rho);         //image rotation angle
//
//            src.SetXY(v);               //source position in the camera
//
//            if (!hsrc.Calc(hil))        //calc source dependant hillas
//                return;
//
//            //fill absolute alpha into histogram
//            const Double_t alpha = hsrc.GetAlpha();
//            fHist.Fill(cx[ix], cy[iy], TMath::Abs(alpha), w);
//        }
//    }
//
// Use MHFalse Source like this:
// -----------------------------
//    MFillH fill("MHFalseSource", "MHillas");
// or
//    MHFalseSource hist;
//    hist.SetAlphaCut(12.5);  // The default value
//    hist.SetBgmean(55);      // The default value
//    MFillH fill(&hist, "MHillas");
//
// To be implemented:
// ------------------
//  - a switch to use alpha or |alpha|
//  - taking the binning for alpha from the parlist (binning is
//    currently fixed)
//  - a possibility to change the fit-function (eg using a different TF1)
//  - a possibility to change the fit algorithm (eg which paremeters
//    are fixed at which time)
//  - use a different varaible than alpha
//  - a possiblity to change the algorithm which is used to calculate
//    alpha (or another parameter) is missing (this could be done using
//    a tasklist somewhere...)
//  - a more clever (and faster) algorithm to fill the histogram, eg by
//    calculating alpha once and fill the two coils around the mean
//
//////////////////////////////////////////////////////////////////////////////
#include "MHFalseSource.h"

#include <TF1.h>
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>

#include "MGeomCam.h"
#include "MSrcPosCam.h"
#include "MHillasSrc.h"
#include "MTime.h"
#include "MObservatory.h"
#include "MPointingPos.h"
#include "MAstroSky2Local.h"
#include "MStatusDisplay.h"

#include "MBinning.h"
#include "MParList.h"

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

ClassImp(MHFalseSource);

using namespace std;

// --------------------------------------------------------------------------
//
// Default Constructor
//
MHFalseSource::MHFalseSource(const char *name, const char *title)
    : fMm2Deg(-1), fUseMmScale(kTRUE), fAlphaCut(12.5), fBgMean(55)
{
    //
    //   set the name and title of this object
    //
    fName  = name  ? name  : "MHFalseSource";
    fTitle = title ? title : "3D-plot of Alpha vs x, y";

    fHist.SetDirectory(NULL);

    fHist.SetName("Alpha");
    fHist.SetTitle("3D-plot of Alpha vs x, y");
    fHist.SetXTitle("x [\\circ]");
    fHist.SetYTitle("y [\\circ]");
    fHist.SetZTitle("\\alpha [\\circ]");
}

// --------------------------------------------------------------------------
//
// Set the alpha cut (|alpha|<fAlphaCut) which is use to estimate the
// number of excess events
//
void MHFalseSource::SetAlphaCut(Float_t alpha)
{
    if (alpha<0)
        *fLog << warn << "Alpha<0... taking absolute value." << endl;

    fAlphaCut = TMath::Abs(alpha);

    Modified();
}

// --------------------------------------------------------------------------
//
// Set mean alpha around which the off sample is determined
// (fBgMean-fAlphaCut/2<|fAlpha|<fBgMean+fAlphaCut/2) which is use
// to estimate the number of off events
//
void MHFalseSource::SetBgMean(Float_t alpha)
{
    if (alpha<0)
        *fLog << warn << "Alpha<0... taking absolute value." << endl;

    fBgMean = TMath::Abs(alpha);

    Modified();
}

// --------------------------------------------------------------------------
//
// Use this function to setup your own conversion factor between degrees
// and millimeters. The conversion factor should be the one calculated in
// MGeomCam. Use this function with Caution: You could create wrong values
// by setting up your own scale factor.
//
void MHFalseSource::SetMm2Deg(Float_t mmdeg)
{
    if (mmdeg<0)
    {
        *fLog << warn << dbginf << "Warning - Conversion factor < 0 - nonsense. Ignored." << endl;
        return;
    }

    if (fMm2Deg>=0)
        *fLog << warn << dbginf << "Warning - Conversion factor already set. Overwriting" << endl;

    fMm2Deg = mmdeg;
}

// --------------------------------------------------------------------------
//
// With this function you can convert the histogram ('on the fly') between
// degrees and millimeters.
//
void MHFalseSource::SetMmScale(Bool_t mmscale)
{
    if (fUseMmScale == mmscale)
        return;

    if (fMm2Deg<0)
    {
        *fLog << warn << dbginf << "Warning - Sorry, no conversion factor for conversion available." << endl;
        return;
    }

    if (fUseMmScale)
    {
        fHist.SetXTitle("x [mm]");
        fHist.SetYTitle("y [mm]");

        fHist.Scale(1./fMm2Deg);
    }
    else
    {
        fHist.SetXTitle("x [\\circ]");
        fHist.SetYTitle("y [\\circ]");

        fHist.Scale(1./fMm2Deg);
    }

    fUseMmScale = mmscale;
}

// --------------------------------------------------------------------------
//
// Calculate Significance as
// significance = (s-b)/sqrt(s+k*k*b) mit k=s/b
// 
Double_t MHFalseSource::Significance(Double_t s, Double_t b)
{
    const Double_t k = b==0 ? 0 : s/b;
    const Double_t f = s+k*k*b;

    return f==0 ? 0 : (s-b)/TMath::Sqrt(f);
}

// --------------------------------------------------------------------------
//
// Set binnings (takes BinningFalseSource) and prepare filling of the
// histogram.
//
// Also search for MTime, MObservatory and MPointingPos
// 
Bool_t MHFalseSource::SetupFill(const MParList *plist)
{
    const MGeomCam *geom = (MGeomCam*)plist->FindObject("MGeomCam");
    if (geom)
    {
        fMm2Deg = geom->GetConvMm2Deg();
        fUseMmScale = kFALSE;

        fHist.SetXTitle("x [\\circ]");
        fHist.SetYTitle("y [\\circ]");
    }

    MBinning binsa;
    binsa.SetEdges(18, 0, 90);

    const MBinning *bins = (MBinning*)plist->FindObject("BinningFalseSource");
    if (!bins)
    {
        Float_t r = geom ? geom->GetMaxRadius() : 600;
        r /= 3;
        if (!fUseMmScale)
            r *= fMm2Deg;

        MBinning b;
        b.SetEdges(100, -r, r);
        SetBinning(&fHist, &b, &b, &binsa);
    }
    else
        SetBinning(&fHist, bins, bins, &binsa);

    fPointPos = (MPointingPos*)plist->FindObject(AddSerialNumber("MPointingPos"));
    if (!fPointPos)
        *fLog << warn << "MPointingPos not found... no derotation." << endl;

    fTime = (MTime*)plist->FindObject(AddSerialNumber("MTime"));
    if (!fTime)
        *fLog << warn << "MTime not found... no derotation." << endl;

    fObservatory = (MObservatory*)plist->FindObject(AddSerialNumber("MObservatory"));
    if (!fObservatory)
        *fLog << err << "MObservatory not found...  no derotation." << endl;


    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Fill the histogram. For details see the code or the class description
// 
Bool_t MHFalseSource::Fill(const MParContainer *par, const Stat_t w)
{
    MHillas *hil = (MHillas*)par;

    const Double_t maxr = 0.98*TMath::Abs(fHist.GetBinCenter(1));

    // Get camera rotation angle
    Double_t rho = 0;
    if (fTime && fObservatory && fPointPos)
    {
        const Double_t ra  = fPointPos->GetRa();
        const Double_t dec = fPointPos->GetDec();

        rho = MAstroSky2Local(*fTime, *fObservatory).RotationAngle(ra, dec);
    }

    MSrcPosCam src;
    MHillasSrc hsrc;
    hsrc.SetSrcPos(&src);

    const Int_t nx = fHist.GetNbinsX();
    const Int_t ny = fHist.GetNbinsY();

    Axis_t cx[nx];
    Axis_t cy[ny];
    fHist.GetXaxis()->GetCenter(cx);
    fHist.GetYaxis()->GetCenter(cy);

    for (int ix=0; ix<nx; ix++)
    {
        for (int iy=0; iy<ny; iy++)
        {
            if (TMath::Hypot(cx[ix], cy[iy])>maxr)
                continue;

            TVector2 v(cx[ix], cy[iy]);
            if (rho!=0)
                v.Rotate(-rho);

            if (!fUseMmScale)
                v *= 1./fMm2Deg;

            src.SetXY(v);

            if (!hsrc.Calc(hil))
            {
                *fLog << warn << "Calculation of MHillasSrc failed for x=" << cx[ix] << " y=" << cy[iy] << endl;
                return kFALSE;
            }

            const Double_t alpha = hsrc.GetAlpha();
            fHist.Fill(cx[ix], cy[iy], TMath::Abs(alpha), w);
        }
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Create projection for off data, taking sum of bin contents of
// range (fBgMean-fAlphaCut/2, fBgMean+fAlphaCut) Making sure to take
// the same number of bins than for on-data
//
void MHFalseSource::ProjectOff(TH2D *h2)
{
    TAxis &axe = *fHist.GetZaxis();

    // Find range to cut (left edge and width)
    const Int_t f = axe.FindBin(fBgMean-fAlphaCut/2);
    const Int_t l = axe.FindBin(fAlphaCut)+f-1;

    axe.SetRange(f, l);
    const Float_t cut1 = axe.GetBinLowEdge(f);
    const Float_t cut2 = axe.GetBinUpEdge(l);
    h2->SetTitle(Form("Distribution of %.1f\\circ<|\\alpha|<%.1f\\circ in x,y", cut1, cut2));

    // Get projection for range
    TH2D *p = (TH2D*)fHist.Project3D("xy_off");

    // Reset range
    axe.SetRange(0,9999);

    // Move contents from projection to h2
    h2->Reset();
    h2->Add(p);

    // Delete p
    delete p;

    // Set Minimum as minimum value Greater Than 0
    h2->SetMinimum(GetMinimumGT(*h2));
}

// --------------------------------------------------------------------------
//
// Create projection for on data, taking sum of bin contents of
// range (0, fAlphaCut)
//
void MHFalseSource::ProjectOn(TH2D *h3)
{
    TAxis &axe = *fHist.GetZaxis();

    // Find range to cut
    axe.SetRangeUser(0, fAlphaCut);
    const Float_t cut = axe.GetBinUpEdge(axe.GetLast());
    h3->SetTitle(Form("Distribution of |\\alpha|<%.1f\\circ in x,y", cut));

    // Get projection for range
    TH2D *p = (TH2D*)fHist.Project3D("xy_on");

    // Reset range
    axe.SetRange(0,9999);

    // Move contents from projection to h3
    h3->Reset();
    h3->Add(p);
    delete p;

    // Set Minimum as minimum value Greater Than 0
    h3->SetMinimum(GetMinimumGT(*h3));
}

// --------------------------------------------------------------------------
//
// Update the projections
//
void MHFalseSource::Paint(Option_t *opt)
{
    // sigma = (s-b)/sqrt(s+k*k*b) mit k=s/b

    gStyle->SetPalette(1, 0);

    TVirtualPad *padsave = gPad;

    TH1D* h1;
    TH2D* h2;
    TH2D* h3;
    TH2D* h4;

    padsave->cd(3);
    if ((h3 = (TH2D*)gPad->FindObject("Alpha_xy_on")))
        ProjectOn(h3);

    padsave->cd(4);
    if ((h2 = (TH2D*)gPad->FindObject("Alpha_xy_off")))
        ProjectOff(h2);

    padsave->cd(2);
    if (h2 && h3 && (h4 = (TH2D*)gPad->FindObject("Alpha_xy_sig")))
    {
        const Int_t nx = h4->GetXaxis()->GetNbins();
        const Int_t ny = h4->GetYaxis()->GetNbins();

        Int_t maxx=nx/2;
        Int_t maxy=ny/2;

        Int_t max = h4->GetBin(maxx, maxy);

        for (int ix=0; ix<nx; ix++)
            for (int iy=0; iy<ny; iy++)
            {
                const Int_t n = h4->GetBin(ix+1, iy+1);

                const Double_t s = h3->GetBinContent(n);
                const Double_t b = h2->GetBinContent(n);

                const Double_t sig = Significance(s, b);

                h4->SetBinContent(n, sig);

                if (sig>h4->GetBinContent(max) && sig!=0)
                {
                    max = n;
                    maxx=ix;
                    maxy=iy;
                }
            }

        padsave->cd(1);
        if ((h1 = (TH1D*)gPad->FindObject("Alpha")))
        {
            //h1->Reset();

            const Double_t x = fHist.GetXaxis()->GetBinCenter(maxx);
            const Double_t y = fHist.GetYaxis()->GetBinCenter(maxy);
            const Double_t s = h4->GetBinContent(max);

            TH1 *h = fHist.ProjectionZ("Alpha", maxx, maxx, maxy, maxy);
            h->SetTitle(Form("Distribution of \\alpha for x=%.2f y=%.2f (\\sigma=%.1f)", x, y, s));
        }
    }

    gPad = padsave;
}

// --------------------------------------------------------------------------
//
// Draw the histogram
//
void MHFalseSource::Draw(Option_t *opt)
{
    TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
    pad->SetBorderMode(0);

    AppendPad("");

    pad->Divide(2, 2);

    // draw the 2D histogram Sigmabar versus Theta
    pad->cd(1);
    gPad->SetBorderMode(0);
    TH1 *h1 = fHist.ProjectionZ("Alpha");
    h1->SetDirectory(NULL);
    h1->SetTitle("Distribution of \\alpha");
    h1->SetXTitle(fHist.GetZaxis()->GetTitle());
    h1->SetYTitle("Counts");
    h1->Draw(opt);
    h1->SetBit(kCanDelete);

    pad->cd(4);
    gPad->SetBorderMode(0);
    fHist.GetZaxis()->SetRangeUser(fBgMean-fAlphaCut/2, fBgMean+fAlphaCut/2);
    TH1 *h2 = fHist.Project3D("xy_off");
    h2->SetDirectory(NULL);
    h2->SetXTitle(fHist.GetXaxis()->GetTitle());
    h2->SetYTitle(fHist.GetYaxis()->GetTitle());
    h2->Draw("colz");
    h2->SetBit(kCanDelete);

    pad->cd(3);
    gPad->SetBorderMode(0);
    fHist.GetZaxis()->SetRangeUser(0,fAlphaCut);
    TH1 *h3 = fHist.Project3D("xy_on");
    fHist.GetZaxis()->SetRange(0,9999);
    h3->SetDirectory(NULL);
    h3->SetXTitle(fHist.GetXaxis()->GetTitle());
    h3->SetYTitle(fHist.GetYaxis()->GetTitle());
    h3->Draw("colz");
    h3->SetBit(kCanDelete);

    pad->cd(2);
    gPad->SetBorderMode(0);
    fHist.GetZaxis()->SetRange(0,0);
    TH1 *h4 = fHist.Project3D("xy_sig"); // Do this to get the correct binning....
    fHist.GetZaxis()->SetRange(0,9999);
    h4->SetTitle("Significance");
    h4->SetDirectory(NULL);
    h4->SetXTitle(fHist.GetXaxis()->GetTitle());
    h4->SetYTitle(fHist.GetYaxis()->GetTitle());
    h4->Draw("colz");
    h4->SetBit(kCanDelete);

}

// --------------------------------------------------------------------------
//
// Everything which is in the main pad belongs to this class!
//
Int_t MHFalseSource::DistancetoPrimitive(Int_t px, Int_t py)
{
    return 0;
}

// --------------------------------------------------------------------------
//
// Set all sub-pads to Modified()
//
void MHFalseSource::Modified()
{
    if (!gPad)
        return;

    TVirtualPad *padsave = gPad;
    padsave->Modified();
    padsave->cd(1);
    gPad->Modified();
    padsave->cd(2);
    gPad->Modified();
    padsave->cd(3);
    gPad->Modified();
    padsave->cd(4);
    gPad->Modified();
    gPad->cd();
}

Double_t fcn(Double_t *arg, Double_t *p)
{
    const Double_t x = arg[0];

    const Double_t dx = (x-p[1])/p[2];

    const Double_t f1 = p[0]*TMath::Exp(-0.5*dx*dx);
    const Double_t f2 = p[3]*x*x + p[4];

    return f1 + f2;
}

Double_t FcnI1(Double_t x, Double_t *p)
{
    return (p[3]*x*x/3+p[4])*x;
}
Double_t FcnI2(Double_t x, Double_t *p)
{
    static const Double_t sqrt2   = TMath::Sqrt(2.);
    static const Double_t sqrt2pi = TMath::Sqrt(TMath::TwoPi());

    const Double_t dx = (x-p[1])/p[2];

    const Double_t f2 = p[0]*p[2]*sqrt2pi*TMath::Erf(dx/sqrt2)/2;

    return f2;
}


void MHFalseSource::FitSignificance()
{
    TH1D h0a("A",     "Parameter A",      50,  0, 4000);
    TH1D h0b("a",     "Parameter a",      50,  0, 4000);
    TH1D h1("mu",    "Parameter mu",      50, -1, 1);
    TH1D h2("sigma", "Parameter sigma",   50,  0, 90);
    TH1D h3("b",     "Parameter b",       50, -0.1, 0.1);
    TH1D h4a("chisq1", "\\chi^{2}",       50,  0, 35);
    TH1D h5a("prob1",  "Fit probability", 50,  0, 1.1);
    TH1D h4b("chisq2", "\\chi^{2}",       50,  0, 35);
    TH1D h5b("prob2",  "Fit probability", 50,  0, 1.1);
    TH1D h6("Signif",  "Significance", 50,  -20, 20);
    h0a.SetDirectory(NULL);
    h0b.SetDirectory(NULL);
    h1.SetDirectory(NULL);
    h2.SetDirectory(NULL);
    h3.SetDirectory(NULL);
    h4a.SetDirectory(NULL);
    h5b.SetDirectory(NULL);
    h4a.SetDirectory(NULL);
    h5b.SetDirectory(NULL);
    h6.SetDirectory(NULL);

    h0a.SetLineColor(kBlue);
    h4a.SetLineColor(kBlue);
    h5a.SetLineColor(kBlue);
    h0b.SetLineColor(kRed);
    h4b.SetLineColor(kRed);
    h5b.SetLineColor(kRed);

    TH1 *hist = fHist.Project3D("xy_fit");

    //                      xmin, xmax, npar
    TF1 func("MyFunc", fcn, 0,    90,   5);

    TArrayD samx(50);
    TArrayD samw(50);

    // func.CalcGaussLegendreSamplingPoints(50, samx.GetArray(), samw.GetArray());
    // func.IntegralFast(50, samx.GetArray(), samw.GetArray(), 0, 15);

    func.SetParName(0, "A");
    func.SetParName(1, "mu");
    func.SetParName(2, "sigma");
    func.SetParName(3, "a");
    func.SetParName(4, "b");

    func.SetParLimits(3, -1, 1);

    const Int_t nx = fHist.GetXaxis()->GetNbins();
    const Int_t ny = fHist.GetYaxis()->GetNbins();

    Double_t maxs=3;
    TH1 *result=0;
    TF1 *fres=0;


    TH1 *h=0;
    for (int ix=0; ix<nx; ix++)
        for (int iy=0; iy<ny; iy++)
        {
            h = fHist.ProjectionZ("AlphaFit", ix+1, ix+1, iy+1, iy+1);

            // Check for the regios which is not filled...
            if (h->GetBinContent(1)==0)
                continue;

            // First fit a polynom in the off region
            func.FixParameter(0, 0);
            func.FixParameter(1, 0);
            func.FixParameter(2, 1);
            func.ReleaseParameter(3);
            func.ReleaseParameter(4);

            h->Fit(&func, "N0Q", "", 35, 75);
            //*fLog << dbg << ix << "/" << iy << ":  " << func.GetParameter(3) << "    " << func.GetParameter(4) << endl;

            h4a.Fill(func.GetChisquare());
            h5a.Fill(func.GetProb());

            // Now fit a gaus in the on region on top of the polynom
            func.SetParameter(0, h->GetBinContent(1)-func.GetParameter(4));
            func.SetParameter(2, 80);

            //func.SetParLimits(0, 0.5*h->GetBinContent(1), 1.5*h->GetBinContent(1));
            //func.SetParLimits(2, 5, 90);

            func.ReleaseParameter(0);
            //func.ReleaseParameter(1);
            func.ReleaseParameter(2);
            func.FixParameter(3, func.GetParameter(3));
            func.FixParameter(4, func.GetParameter(4));

            func.SetParLimits(2, 0, 80);
            h->Fit(&func, "N0Q", "", 0, 35);
            //*fLog << dbg << "     " << func.GetParameter(0) << "    " << func.GetParameter(1) << "    " << func.GetParameter(2) << endl;

            // Fill results into some histograms
            h0a.Fill(func.GetParameter(0));
            h0b.Fill(func.GetParameter(4));
            h1.Fill(func.GetParameter(1));
            h2.Fill(func.GetParameter(2));
            h3.Fill(func.GetParameter(3));
            h4b.Fill(func.GetChisquare());
            h5b.Fill(func.GetProb());

            const Int_t n = hist->GetBin(ix+1, iy+1);
            if (func.GetParameter(0)>h->GetBinContent(1)*2 ||
                func.GetParameter(2)<2.5 || func.GetParameter(2)>70
                /*func.GetProb()<0.005 ||*/)
            {
                hist->SetBinContent(n, 0);
                continue;
            }

            //*fLog << dbg << "     Chisq=" << func.GetChisquare() << " NDF=" << func.GetNDF()
            //    << " Prob=" << func.GetProb() << " FitP=" << func.GetNumberFitPoints() << endl;

            Double_t b = FcnI1(15, func.GetParameters());
            Double_t s = FcnI2(15, func.GetParameters());

            //*fLog << dbg << func.GetParameter(3) << ": " << b << "    ";
            //*fLog << dbg << func.GetParameter(0) << ": " << s;
            //*fLog << endl;

            if (b<0)
                continue;

            const Double_t sig = Significance(s+b, b);

            hist->SetBinContent(n, sig);
            h6.Fill(sig);

            if (sig>maxs)
            {
                maxs = sig;
                if (result)
                {
                    delete result;
                    delete fres;
                }
                result=(TH1*)h->Clone("Result \\alpha");
                fres  =(TF1*)func.Clone("Result Func");
            }
        }

    hist->SetMinimum(-10);
    hist->SetMaximum(maxs);

    TCanvas *c=new TCanvas;

    c->Divide(3,2);
    c->cd(1);
    h0b.DrawCopy();
    h0a.DrawCopy("same");
    c->cd(2);
    hist->Draw("colz");
    hist->SetDirectory(NULL);
    hist->SetBit(kCanDelete);
    c->cd(3);
    h2.DrawCopy();
    c->cd(4);
    h3.DrawCopy();
    c->cd(5);
    h4a.DrawCopy();
    h4b.DrawCopy("same");
    h6.DrawCopy("same");
    c->cd(6);
    h5a.DrawCopy();
    h5b.DrawCopy("same");

    if (result)
    {
        c=new TCanvas;
        result->SetBit(kCanDelete);
        result->Draw();

        TF1 f1("MyFunc1", fcn, 0, 90, 5);
        TF1 f2("MyFunc2", fcn, 0, 90, 5);
        f1.SetParameters(fres->GetParameters());
        f2.SetParameters(fres->GetParameters());
        f2.FixParameter(0, 0);
        f2.FixParameter(1, 0);
        f2.FixParameter(2, 1);
        f1.SetLineColor(kGreen);
        f2.SetLineColor(kRed);

        f2.DrawCopy("same");
        f1.DrawCopy("same");

        delete fres;
    }
}