PeakCreator.cpp

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2004-2006 by Olivier Stezowski & Christian Finck        *
 *   stezow(AT)ipnl.in2p3.fr, cfinck(AT)ires.in2p3.fr                      *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#include <Riostream.h>


#ifndef ROOT_KeySymbols
#include "KeySymbols.h"
#endif

#ifndef ROOT_TCanvas
#include "TCanvas.h"
#endif

#ifndef ROOT_TH1
#include "TH1.h"
#endif

#ifndef ROOT_TList
#include "TList.h"
#endif

#ifndef ROOT_TMath
#include "TMath.h"
#endif

#ifndef ROOT_TROOT
#include "TROOT.h"
#endif

#include "TContextMenu.h"
#include "TFrame.h"

#ifndef Gw_PadManager
#include "PadManager.h"
#endif

#ifndef Gw_PeakCreator
#include "PeakCreator.h"
#endif

#include "TGraph.h"

using namespace Gw;

TList PeakCreator::gPeakCreatorList; TF1 *PeakCreator::gDefaultFWHM = 0x0;

//__________________________________________________________
PeakCreator::PeakCreator() 
    : TObject(),
      fEfficiencyGraph(0x0),
      fEfficiencyFunc(0x0),
      fRefArea(1),
      fRefAreaError(0),
      fDefaultPeak1D(0x0),
      fIsCollecting(false),
      fIsGate(false),
      fDefaultPeakFWHM(0x0),
      fCanvas(0x0),
      fLastX(0),
      fLastY(0),
      fLastMouseX(0),
      fLastMouseY(0),
      fIsHelpsPrintActive(true),
      fLog("PeakCreator")
{ 
    // default constructor
    fDefaultPeak1D   = new Peak1D();
    // to be sure it exists
    GetGlobalPeakFWHM();
    SetDefaultPeakFWHM("PCG_FWHM"); // PeakCreator Global _ FWHM
}

//__________________________________________________________
PeakCreator::PeakCreator(TCanvas *c) : TObject(),
    fEfficiencyGraph(0x0),
    fEfficiencyFunc(0x0),
    fRefArea(1),
    fRefAreaError(0),
    fDefaultPeak1D(0x0),
    fIsCollecting(true),
    fIsGate(false),
    fDefaultPeakFWHM(0x0),
    fCanvas(c),
    fLastX(0),
    fLastY(0),
    fLastMouseX(0),
    fLastMouseY(0),
    fIsHelpsPrintActive(true),
    fLog("PeakCreator")
{ 
    // default constructor
    fDefaultPeak1D   = new Peak1D();
    GetGlobalPeakFWHM();
    SetDefaultPeakFWHM("PCG_FWHM"); // PeakCreator Global _ FWHM

    // recursive delete
    SetBit(TObject::kMustCleanup);

    Connect();

    gPeakCreatorList.Add(this);
    //  gPeakCreatorList.ls();
}

//__________________________________________________________
PeakCreator::~PeakCreator()
{
    if ( fLog.IsDebug(5) )
        fLog << info << "Destroy PeakCreator " << nline;

    Disconnect();

    // default destructor
    if (fDefaultPeak1D)
    { delete fDefaultPeak1D; fDefaultPeak1D = 0x0; }
    if (fDefaultPeakFWHM)
    { delete fDefaultPeakFWHM; fDefaultPeakFWHM = 0x0; }

    if ( fLog.IsDebug(5) )
        fLog << dolog;
}

//__________________________________________________________

void PeakCreator::Connect()
{
    if(fCanvas == 0x0)
    {
        if ( TVirtualPad::Pad() )
            fCanvas = TVirtualPad::Pad()->GetCanvas();
        else {
            fLog.SetProcessMethod("Connect()");
            fLog << error << "gPad not existing" << dolog;
            return;
        }
    }

    fCanvas->Connect("ProcessedEvent(Int_t, Int_t, Int_t, TObject*)",
                     "Gw::PeakCreator",
                     this,
                     "HandleMovement(Int_t, Int_t, Int_t, TObject*)");

    fCanvas->Connect("Closed()",
                     "Gw::PeakCreator",
                     this,
                     "Delete()");
}

//__________________________________________________________

void PeakCreator::Disconnect()
{
    // in case this is connected to a canvas
    if ( fCanvas ) {

        fCanvas->Disconnect("ProcessedEvent(Int_t, Int_t, Int_t, TObject*)",
                            this,
                            "HandleMovement(Int_t, Int_t, Int_t, TObject*)");

        fCanvas->Disconnect("Closed()",
                            this,
                            "Delete()");

        gPeakCreatorList.Remove(this);
    }
}

//__________________________________________________________

TVirtualPad *PeakCreator::PadforAction(TVirtualPad *pad)
{ 
    if ( pad )
        return pad ;
    return TVirtualPad::Pad();
}

PeakCreator *PeakCreator::AddPeakCreator(TVirtualPad *pad)
{
    PeakCreator *pcreator = 0x0;

    // to get the pad on which the creator is added
    TVirtualPad *locpad = PeakCreator::PadforAction(pad);
    if ( locpad == 0x0 )
        return pcreator;

    TIter next(&gPeakCreatorList);
    while ( (pcreator = (PeakCreator *)next()) ) {
        if ( pcreator->GetCanvas() == locpad->GetCanvas() )
            return pcreator;
    }

    // create a PeakCreator and add it to the current canvas
    pcreator = new PeakCreator(locpad->GetCanvas());
    return pcreator;
}

PeakCreator *PeakCreator::IsPeakCreator(TVirtualPad *pad)
{
    PeakCreator *pcreator = 0x0;

    // to get the pad on which the creator is added
    TVirtualPad *locpad = PeakCreator::PadforAction(pad);
    if ( locpad == 0x0 )
        return pcreator;

    TIter next(&gPeakCreatorList);
    while ( (pcreator = (PeakCreator *)next()) ) {
        if ( pcreator->GetCanvas() == locpad->GetCanvas() ) {
            break;
        }
    }
    return pcreator;
}

void PeakCreator::Clear(Option_t * opt)
{
    TString o = opt;
    TObject *obj = 0x0;

    // collect peaks in the selected canvas
    TList listofpeaks;
    PadManager::Collect<BasePeak>(listofpeaks,fCanvas->GetSelectedPad());
    TIter nextp(&listofpeaks);
    while ( (obj = nextp()) ) {
        fCanvas->GetListOfPrimitives()->Remove(obj);
    }
    listofpeaks.Clear();

    PadManager::Collect<TPolyLine>(listofpeaks,fCanvas->GetSelectedPad());
    TIter nextpo(&listofpeaks);
    while ( (obj = nextpo()) ) {
        fCanvas->GetListOfPrimitives()->Remove(obj);
    }
    listofpeaks.Clear();
    TList listoffunctions;

    PadManager::Collect<TF1>(listoffunctions,fCanvas->GetSelectedPad());
    TIter nextf(&listoffunctions);
    while ( (obj = nextf()) ) {
        fCanvas->GetListOfPrimitives()->Remove(obj);
    }
    listoffunctions.Clear();
}

//__________________________________________________________
void PeakCreator::Collect(Bool_t do_collect)
{
    fIsCollecting = do_collect;
}

//__________________________________________________________
void PeakCreator::SetGateMode(Bool_t ngate)
{
    if ( ngate )
        fDefaultPeak1D->SetDrawAs(Peak1D::kGate);
    else
        fDefaultPeak1D->SetDrawAs(Peak1D::kPeak);
}

//__________________________________________________________
void PeakCreator::HandleMovement(Int_t eventType, Int_t eventX, Int_t eventY, TObject* /*select*/)
{       
    // handle cursor mouvement
    TVirtualPad *pad = fCanvas->GetSelectedPad();
    if ( pad == 0x0 )
        return;
    if ( !pad->IsEditable() ) {
        return ;
    }
    if (eventType == kMouseMotion) {
        // to keep track of the last mouse position
        fLastMouseX = eventX;
        fLastMouseY = eventY;
    }
    if( eventType == kButton1Up ) {
        HandleRefresh(fCanvas);
    }
    if( eventType == kKeyPress) {

        bool CTRL = ( eventX == eventY-96);
        EKeySym Key = (EKeySym)eventY;

        bool LastCTRL = ( fLastX == fLastY-96);
        EKeySym LastKey = (EKeySym)fLastY;

        if(Key == kKey_0 && !CTRL && LastKey == kKey_c && !LastCTRL)
        {
            fLog << info << "Current color/style for new peak is now style 0" << nline;
            fDefaultPeak1D->SetLineColorPeak(BasePeak::gGateColor.k0);
            fDefaultPeak1D->SetFillColorPeak(BasePeak::gGateColor.k0);
        }
        if(Key == kKey_1 && !CTRL && LastKey == kKey_c && !LastCTRL)
        {
            fLog << info << "Current color/style for new peak is now style 1" << nline;
            fDefaultPeak1D->SetLineColorPeak(BasePeak::gGateColor.k1);
            fDefaultPeak1D->SetFillColorPeak(BasePeak::gGateColor.k1);
        }
        if(Key == kKey_2 && !CTRL && LastKey == kKey_c && !LastCTRL)
        {
            fLog << info << "Current color/style for new peak is now style 2" << nline;
            fDefaultPeak1D->SetLineColorPeak(BasePeak::gGateColor.k2);
            fDefaultPeak1D->SetFillColorPeak(BasePeak::gGateColor.k2);
        }
        if(Key == kKey_3 && !CTRL && LastKey == kKey_c && !LastCTRL)
        {
            fLog << info << "Current color/style for new peak is now style 3" << nline;
            fDefaultPeak1D->SetLineColorPeak(BasePeak::gGateColor.k3);
            fDefaultPeak1D->SetLineColorPeak(BasePeak::gGateColor.k3);
        }
        if(Key == kKey_p && !CTRL && LastKey == kKey_a && !LastCTRL)
        {
            if ( fIsCollecting ) {
                BasePeak *peak = CreatePeak(pad->AbsPixeltoX(fLastMouseX),"",pad);
                if ( peak && pad->GetFrame() )
                    peak->PaintFor(pad->GetFrame()->GetX1(),pad->GetFrame()->GetX2(),pad->GetFrame()->GetY1(),pad->GetFrame()->GetY2());
                pad->Update();
            }
        }
        if(Key == kKey_P && !CTRL && LastKey == kKey_A && !LastCTRL)
        {
            if ( fIsCollecting ) {
                BasePeak *peak = CreatePeak(pad->AbsPixeltoX(fLastMouseX),"bg",pad);
                if ( peak && pad->GetFrame() )
                    peak->PaintFor(pad->GetFrame()->GetX1(),pad->GetFrame()->GetX2(),pad->GetFrame()->GetY1(),pad->GetFrame()->GetY2());
                pad->Update();
            }
        }
        if(Key == kKey_p && CTRL && LastKey == kKey_a && LastCTRL)
        {
            if ( fIsCollecting ) {
                AddPointToPolyline(pad->AbsPixeltoX(fLastMouseX), pad->AbsPixeltoY(fLastMouseY),pad);
                pad->Update();
            }
        }
        if(Key == kKey_a && !CTRL && LastKey == kKey_c && !LastCTRL)
        {
            Clear();
            pad->Update();
        }
        if(Key == kKey_m && !CTRL && LastKey == kKey_o && !LastCTRL)
        {
            fCanvas->GetContextMenu()->Popup(fLastMouseX,fLastMouseY,this,fCanvas,pad);
        }
        if(Key == kKey_h && !CTRL && fIsHelpsPrintActive)
        {
            fLog << info;
            fLog << "\ta: Add a peak at mouse position" << nline;
            fLog << "\tA: Add a peak + background point at mouse position" << nline;
            fLog << "\tp: Add a point to the graphical polyline" << nline;
            fLog << "\t0: Switch Color to default #0" << nline;
            fLog << "\t1: Switch Color to default #1" << nline;
            fLog << "\t2: Switch Color to default #2" << nline;
            fLog << "\t3: Switch Color to default #3" << nline;
            fLog << "\tc: Clear all (peaks + graphical polyline)" << nline;

            fLog << "\t>: Popup PeakCreator menu" << nline;
            fLog << "\tM: Popup PeakCreator menu" << nline;
        }
    }

    fLastX = eventX;
    fLastY = eventY;

    fLog << dolog;
}

void PeakCreator::HandleRefresh(TCanvas *canvas)
{
    TCanvas *localCanvas = canvas;

    // in this case connect the current canvas
    if ( localCanvas == 0x0 ) {
        if ( TVirtualPad::Pad() )
          localCanvas = TVirtualPad::Pad()->GetCanvas();
    }
    if(localCanvas == 0x0) return;
    if(localCanvas->GetSelectedPad() == 0x0) return;

    TObject *obj; TFrame *frame; BasePeak *apeak;
    double xmin, xmax, ymin, ymax;

    // if no Frame, go out
    if ( (frame = localCanvas->GetSelectedPad()->GetFrame()) == 0x0 )
        return;

    xmin = frame->GetX1();
    xmax = frame->GetX2();
    ymin = frame->GetY1();
    ymax = frame->GetY2();

    if(localCanvas->GetLogy())
    {
        ymin = TMath::Exp(frame->GetY1()*TMath::Log(10));
        ymax = TMath::Exp(frame->GetY2()*TMath::Log(10));
    }
    if(localCanvas->GetLogx())
    {
        xmin = TMath::Exp(frame->GetX1()*TMath::Log(10));
        xmax = TMath::Exp(frame->GetX2()*TMath::Log(10));
    }

    // collect peaks in the selected canvas
    TList listofpeaks;
    PadManager::Collect<BasePeak>(listofpeaks,localCanvas->GetSelectedPad());

    //  printf(" %s TEST %d \n",fCanvas->GetSelectedPad()->GetFrame()->GetName(),listofpeaks.GetSize() );
    TIter next(&listofpeaks);
    while ( (obj = next()) ) {
        apeak = dynamic_cast<BasePeak *> (obj);

        if ( apeak )
            apeak->PaintFor(xmin,xmax,ymin,ymax);

        localCanvas->GetSelectedPad()->Modified();
    }
    // clear the list
    listofpeaks.Clear("nodelete");
}

TF1 *PeakCreator::GetGlobalPeakFWHM()
{
    if ( gROOT->GetListOfFunctions()->FindObject("PCG_FWHM") == 0x0 )
    {
        gDefaultFWHM = new TF1("PCG_FWHM","[0]");
        gDefaultFWHM->SetParameter(0,6.0);
    }

    if ( gROOT->GetListOfFunctions()->FindObject("AGATA_FWHM") == 0x0 )
    {
        //FWHM = 1.2 keV for Egamma = 60 keV
        double E1 = 60.;
        double FWHM1 = 1.2;
        //FWHM = 2.1 keV for Egamma = 1333 keV
        double E2 = 1333.;
        double FWHM2 = 2.1;

        new TF1("AGATA_FWHM",Form("(%lf-%lf)/(%lf-%lf)*(x-%lf)+%lf",FWHM2,FWHM1,E2,E1,E2,FWHM2));
    }

    return gDefaultFWHM;
}

const TF1 *PeakCreator::GetConstGlobalPeakFWHM()
{    
    return GetGlobalPeakFWHM();
}

void PeakCreator::SetDefaultPeakFWHM( const char *fname )
{
    TObject *o = gROOT->GetListOfFunctions()->FindObject(fname); TF1 *f;
    if ( o == 0x0 )
        return;
    else
        f = dynamic_cast<TF1 *> (o);

    if ( fDefaultPeakFWHM )
        delete fDefaultPeakFWHM;

    fDefaultPeakFWHM = (TF1 *)(f->Clone());
    gROOT->GetListOfFunctions()->Remove(fDefaultPeakFWHM);
}

//__________________________________________________________
BasePeak *PeakCreator::CreatePeak(Double_t x, Option_t *opt, TVirtualPad *pad)
{
    // create peak from current position in pad
    BasePeak* peak = 0x0; TH1* h;

    // to get the pad on which the creator is added
    TVirtualPad *locpad = PeakCreator::PadforAction(pad);
    if ( locpad == 0x0 )
        return peak;

    h = PadManager::GetHisto(locpad);
    if ( h ) {
        // look in the current pad for the list of peak that already exist
        std::vector<Gw::Peak1D *> peaks_in_pad;
        // collect
        PadManager::Collect<Gw::Peak1D>(peaks_in_pad,locpad);
        // set the name of the peak used as a cloner of the new peak to be created
        TString pname; Int_t i = 0;
        while ( 1 ) {
            pname = Form("Peak_%d",i++); Bool_t is_in = false;
            for (size_t i = 0; i < peaks_in_pad.size(); i++) {
                if ( pname == peaks_in_pad[i]->GetName() ) {
                    is_in = true;
                    break;
                }
            }
            if ( is_in == false ) {
                fDefaultPeak1D->SetName( pname.Data() );
                break;
            }
        }

        // now creates the peak
        peak = CreatePeak(h,x,opt);
        if ( peak ) {
            //
            peak->SetBit(TObject::kCanDelete);
            peak->SetBit(TObject::kMustCleanup);

            // add peak to pad
            locpad->GetListOfPrimitives()->Add(peak,"lm");
            locpad->Modified();
        }
    }
    return peak;
}

//__________________________________________________________
BasePeak* PeakCreator::CreatePeak(const TH1 *histo, Double_t x, Option_t *opt)
{
    // to check options
    TString o = opt;

    // the peak to be created
    BasePeak* peak = 0x0;

    if (histo->GetDimension() == 0 || histo->GetDimension() > 1) {
        fLog.SetProcessMethod("CreatePeak(TH1*, Double_t)");
        fLog << error << "Dimension not handled" << histo->GetDimension() << dolog;
        return peak;
    }
    switch ( histo->GetDimension() ) {
    case 1:
    { // 1D :
        // scan the histogram 6 sigma around the given position
        Double_t width =
                fDefaultPeakFWHM->Eval(x);
        // in case with is greater that the full range of the spectra, reduce width to 1/10 of the range
        if ( histo->GetXaxis()->GetXmax() < histo->GetXaxis()->GetXmin() + width )
            width = ( histo->GetXaxis()->GetXmax() - histo->GetXaxis()->GetXmin() ) / 10.;

        Int_t binmin = histo->GetXaxis()->FindFixBin(x-3*width);
        Int_t binmax = histo->GetXaxis()->FindFixBin(x+3*width);

        // to avoid including under/over flow
        binmin = TMath::Max(1,binmin);
        binmax = TMath::Min(histo->GetNbinsX(),binmax);

        Double_t yminleft, yminright, ymin, ymax, ymax_abs; Int_t bin_ymax, bin_ymax_abs;

        bin_ymax = bin_ymax_abs = histo->GetXaxis()->FindFixBin(x); ymax = ymax_abs = histo->GetBinContent(histo->GetXaxis()->FindFixBin(x));

        yminleft  =  histo->GetBinContent(binmin);
        yminright =  histo->GetBinContent(binmax);

        // first compute some characteritic values
        for (Int_t i = binmin ; i < binmax; i++) {

            if ( histo->GetBinCenter(i) < x-width && histo->GetBinContent(i) < yminleft ) {
                yminleft = histo->GetBinContent(i);
            }
            if ( histo->GetBinCenter(i) > x+width && histo->GetBinContent(i) < yminright ) {
                yminright = histo->GetBinContent(i);
            }
            // fine maximum
            if ( histo->GetBinContent(i) > ymax && histo->GetBinCenter(i) > (x-0.4*width) && histo->GetBinCenter(i) < (x+0.4*width) )
            { ymax = histo->GetBinContent(i); bin_ymax = i; }

            if ( histo->GetBinContent(i) > ymax_abs )
            { ymax_abs = histo->GetBinContent(i); bin_ymax_abs = i; }
        }

        // set maximum ... avoid having the max in a hole
        if ( (ymax < yminleft || ymax < yminright) && !o.Contains("ath") )
            ymax = ymax_abs ;

        // in case of empty histogram
        if ( ymax == 0.0 )
            ymax = 1.0;

        /*                              Double_t DeltaY = TMath::Max ( TMath::Abs(ymax-yminleft) , TMath::Abs(ymax-yminright) ) ;
                for (Int_t i = binmin; i < binmax; i++) {


                }                               */
        // yminleft /= width; yminright /= width;
        ymin = (yminleft+yminright)/2;

        // Set peak properties
        Peak1D* rpeak =
                dynamic_cast<Gw::Peak1D*>(fDefaultPeak1D->Clone());

        rpeak->SetBackground(x-5*width,
                             x-4*width,
                             x+4*width,
                             x+5*width, yminleft, yminleft, yminright, yminright );
        /*
                 histo->GetBinContent(histo->GetXaxis()->FindFixBin(x-3*width)),
                 histo->GetBinContent(histo->GetXaxis()->FindFixBin(x-2*width)),
                 histo->GetBinContent(histo->GetXaxis()->FindFixBin(x+2*width)),
                 histo->GetBinContent(histo->GetXaxis()->FindFixBin(x+3*width)) ); */

        if ( o.Contains("atp") ) // does not change the position
            rpeak->SetPeak(x, ymax-ymin, width);
        else
            rpeak->SetPeak(histo->GetBinCenter(bin_ymax), ymax-ymin, width);

        if ( o.Contains("bg") )
            rpeak->WithBkg(true);
        else
            rpeak->WithBkg(false);

        if(fEfficiencyFunc != 0x0)
            rpeak->SetEfficiencyFunc(fEfficiencyFunc,fEffRelError);
        else if(fEfficiencyGraph != 0x0)
            rpeak->SetEfficiencyGraph(fEfficiencyGraph,fEffRelError);

        rpeak->SetRefArea(fRefArea,fRefAreaError);

        peak = rpeak;
    }
        break;
    default:
        break;
    }

    return peak;
}


//__________________________________________________________
void PeakCreator::AddPointToPolyline(Double_t x, Double_t y, TVirtualPad *pad)
{       
    // to get the pad on which the creator is added
    TVirtualPad *locpad = PeakCreator::PadforAction(pad);
    if ( locpad == 0x0 )
        return ;

    TPolyLine* line = 0x0;
    TList* list = locpad->GetListOfPrimitives();
    for (Int_t i = 0; i < list->GetEntries(); ++i) {
        TString tmp(list->At(i)->ClassName());
        if ( tmp.CompareTo("TPolyLine") == 0 )
            line = static_cast<TPolyLine*>(list->At(i));
    }

    if (!line) {
        line = new TPolyLine(0);

        line->SetBit(kCanDelete, true); line->SetBit(kMustCleanup, true);

        line->SetLineWidth(2);
        line->SetLineColor(kRed-5);
        line->SetFillColor(kRed-5);
        line->SetFillStyle(3010);
        line->SetNextPoint(x,y);

        locpad->GetListOfPrimitives()->Add(line,"f");

    }
    else line->SetNextPoint(x,y);

    locpad->Modified();
}


//
ClassImp(PeakCreator)