GSI/EventPSAWatchers.C

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2010 by Olivier Stezowski                               *
 *   Copyright (C) 2010 by Christian Finck                                 *
 *   stezow(AT)ipnl.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.             *
 ***************************************************************************/

// Watcher definition
#ifndef _EventPSAWatchers
#include "EventPSAWatchers.h"
#endif

#ifndef ROOT_TBox
#include "TBox.h"
#endif

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

#ifndef ROOT_TGeoMatrix
#include "TGeoMatrix.h"
#endif

#ifndef ROOT_TGToolTip
#include "TGToolTip.h"
#endif

#ifndef ROOT_TLatex
#include "TLatex.h"
#endif

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

#ifndef ROOT_TObjArray
#include "TObjArray.h"
#endif

#ifndef ROOT_TPaveText
#include "TPaveText.h"
#endif

#ifndef ROOT_TPolyLine
#include "TPolyLine.h"
#endif

#ifndef ROOT_TRandom
#include "TRandom.h"
#endif

#ifndef ROOT_TString
#include "TString.h"
#endif

#include "TTimer.h"
#include "TSystem.h"

#include "DFAgent.h"
#include "AgataKeyFactory.h"
#include "MetaWatchers.h"

using namespace ADF;

ClassImp(EventPSAWatcher);

//__________________________________________________________
Bool_t EventPSAWatcher::SetTrigger( ADF::DFTrigger *trigger ) 
{
    if  ( !Watcher::GetFromTrigger(trigger, "event:data:psa", fFrame) )
        return false;

    // first extract from the current configuration the definition
    // of the data:psa frame on the data flow
    FactoryItem i_asked, k_defined, f_defined;

    i_asked.SetItem("Agata","data:psa",Version());
    k_defined = ConfAgent::theGlobalAgent()->GetDFAgent()->WhichKnownKey   (i_asked,i_asked);
    f_defined = ConfAgent::theGlobalAgent()->GetDFAgent()->WhichKnownFrame (i_asked,i_asked);

    Frame *f = MainFrameFactory::theMainFactory().New(k_defined,f_defined);
    if ( f ) { // it exists
        if ( fPSAFrame )
        { delete fPSAFrame; fPSAFrame = 0x0; }

        fPSAFrame = f;
    }
    else fFrame = 0x0;

    return fFrame != 0x0;
}

//__________________________________________________________
UInt_t EventPSAWatcher::GetNbFramePSA()
{
    if ( fFrame && fFrame->IsValid() ) {
        fFrame->GetFrame()->Scan();
        return fFrame->GetFrame()->GetNbSubFrame();
    }

    return 0u;
}  

//__________________________________________________________
PSAInterface *EventPSAWatcher::GetDataPSA(UInt_t which)
{
    // does not exist
    if ( !(which < GetNbFramePSA()) ) {
        if ( gDebug > 3 )
            printf("EventPSAWatcher::GetDataPSA %d %d \n",which,GetNbFramePSA());
        return 0x0;
    }
    // cannot load subframe #i in fPSAFrame
    if ( !fFrame->GetFrame()->LinkSubFrame(which,fPSAFrame) ) {
        if ( gDebug > 3 )
            printf("EventPSAWatcher::GetDataPSA cannot LinkSubFrame \n");
        return 0x0;
    }
    // cannot read the data
    if( fPSAFrame->Read() == 0 ) {
        if ( gDebug > 3 )
            printf("EventPSAWatcher::GetDataPSA cannot Read \n");
        return 0x0;
    }
    // return the data interface for that frame
    return GetDataPointer<PSAInterface>(fPSAFrame);
}

//**********************************************************


const Float_t gMetric = 10.;

ClassImp(HitsSpectra);

HitsSpectra::HitsSpectra(const char *name, const char *title, TDirectory *sp_dir, TDirectory *tag_dir) :
    EventPSAWatcher(name,title,sp_dir,tag_dir),
    fVertexBuilder(0x0),
    fTrans(0x0),
    fSpectra(180),
    fSpectraDoppler(180),
    fCoreCommon(0x0),
    fGxG(0x0),
    fGxT(0x0),
    fCoreCommonDoppler(0x0),
    fGxGDoppler(0x0),
    fDistri(0x0),
    fPosXYZ(0x0),
    fPosXYZ_Ta(0x0)
{
    fSpectra.SetOwner(false); fSpectraDoppler.SetOwner(false);

    // Vertex used for doppler correction, last registered
    fVertexBuilder = VertexBuilder::theCurrentVertexBuilder();
    if ( fVertexBuilder == 0x0 )
        std::cout << "Should not be NULL " << std::endl;

    fTrans = new AgataGeometryTransformer();

    if ( gSystem->AccessPathName("CrystalPositionLookUpTable") ) { // does not exist
        if ( gSystem->AccessPathName("CrystalPositionLookUpTable.dat") ) {
            fLog << error << " CrystalPositionLookUpTable required " << nline;
        }
        else {
            fLog << info << " Ge positions extracted from CrystalPositionLookUpTable.dat " << nline;
            fTrans->ReadTransformations("CrystalPositionLookUpTable.dat");
        }
    }
    else {
        fLog << info << " Ge positions extracted from CrystalPositionLookUpTable.dat " << nline;
        fTrans->ReadTransformations("CrystalPositionLookUpTable");
    }

    for (UInt_t i = 0u; i < gNbCrys; i++ ) {

        TString tmp; TH1F *h;

        tmp = Form("Crystal_%d",i);
        h =  MakeTH1<TH1F>(tmp.Data(),"Core energy of the crystal", 4096,0,4096);
        fSpectra.Add ( h ) ;

        tmp = Form("Crystal_track_%d",i);
        h =  MakeTH1<TH1F>(tmp.Data(),"Core energy of the crystal", 4096,0,4096);
        fSpectraDoppler.Add ( h ) ;

    }
    fCoreCommon = MakeTH1<TH1F>("CoreCommon","Core common spectrum",8192,0,4096);
    fSumSpectra = MakeTH1<TH1F>("SumSpectra","Sum of all spectra",8192,0,4096);
    fGxG    = MakeTH2<TH2F>("GxG","Gamma Gamma Matrix",4096,0,4096,4096,0,4096);
    fDistri = MakeTH1<TH1F>("CrystalDistri","Distribution of the crystal fired",180,0,180);
    fPosXYZ = MakeTH3<TH3F>("PosXYZ","Distribution of hits in Space",100,-200,200,100,-200,200,100,-400,0);
    
    // fPos3D =  MakeTH3<TH3F>("Pos3D","Positions of Tracked Gamma-rays ;x;y;z",200,-200,200,200,-200,200,200,-400,0);

    fPosXYZ_Ta = MakeTH3<TH3F>("PosXYZ_Ta","Distribution of hits in Space with Target hit",200,-400,400,200,-400,400,200,-400,400);
    

    
    fGxT = MakeTH2<TH2F>("GxT","Gamma Time Matrix",4096,0,4096,300,0,300);
    
    fCoreCommonDoppler = MakeTH1<TH1F>("DCoreCommon","Core common spectrum with doppler correction",8192,0,4096);
    fGxGDoppler    = MakeTH2<TH2F>("DGxG","Gamma Gamma Matrix with doppler correction",4096,0,4096,4096,0,4096);

    fFold        = MakeTH1<TH1F>("Fold","Crystal Fold distributions",10,0,10);
}

HitsSpectra::~HitsSpectra()
{
    if (fTrans) {
        delete fTrans; fTrans = 0x0;
    }
}

Double_t HitsSpectra::DoDopplerCorrection(const TrackedHit *hit, VertexInterface *vertex)
{
    Double_t e = hit->GetE(), beta = vertex->GetBeta(), gamma = vertex->GetGamma();

    Double_t px, py, pz, dx, dy, dz;
    vertex->GetPosition (px, py, pz);
    vertex->GetDirection(dx, dy, dz);
    
    //    std::cout << " " << px << " " << py << " " << pz <<std::endl;

    Double_t x = hit->GetX() - px;
    Double_t y = hit->GetY() - py;
    Double_t z = hit->GetZ() - pz;

    Double_t dd = x*x + y*y + z*z;
    if( !(dd == 0.0) ) {
        Double_t cosTheta = (x*dx + y*dy + z*dz)/::sqrt(dd);
        e = VertexInterface::DopplerCorrection(e,beta,gamma,cosTheta);
    }
    else e = -1.0;

    return e;
}

void HitsSpectra::DoCanvas(TCanvas *c, Option_t *o)
{
    TString opt = o;

    c->Divide(4,4,0.001,0.001);

    int PadNbr = 1;
    c->cd(PadNbr);
    fDistri->DrawCopy();

    for(int i=0 ; i<fDistri->GetNbinsX() ; i++)
    {
        if(fDistri->GetBinContent(i+1))
        {
            PadNbr++;
            c->cd(PadNbr%16);
            if(fSpectra[i+1]) fSpectra[fDistri->GetBinLowEdge(i+1)]->Draw();
        }
        if(PadNbr%16==0)
        {
            c = NewCanvas(Form("HitSpectra_%d_%d",PadNbr+1,PadNbr+16));
            c->Divide(4,4,0.001,0.001);
        }
    }


    TCanvas *OverImposedCore = NewCanvas("OverImposedCore");
    OverImposedCore->cd();
    bool FirstHist = true;
    for(int i=0 ; i<fDistri->GetNbinsX() ; i++)
    {
        if(fDistri->GetBinContent(i+1) && fSpectra[i+1])
        {
            if(FirstHist)
            {
                fSpectra[fDistri->GetBinLowEdge(i+1)]->Draw();
                FirstHist = false;
            }
            else fSpectra[fDistri->GetBinLowEdge(i+1)]->Draw("same");
        }
    }


    TCanvas *GlobalSpectra = NewCanvas("GlobalSpectra");
    GlobalSpectra->Divide(1,2,0.001,0.001);
    GlobalSpectra->cd(1);
    fCoreCommon->Draw();

    TPad *pad = new TPad("Crystal fold", "Pad for Crystal fold distribution",0.7,0.6,1,1);
    pad->Draw();
    pad->cd();
    fFold->Draw();

    GlobalSpectra->cd(2);
    fSumSpectra->Draw();
}

void HitsSpectra::Exec(Option_t * /*option*/)
{
    //  cout << "Exec Hits sectra" <<endl;
    Double_t xLocal, yLocal, zLocal, xGlobal, yGlobal, zGlobal; TrackedHitImp<EHit,EHit> gamma; Float_t gg[150], ggdc[150];

    // be sure the frame has been set properly
    if ( fFrame == 0x0  || !fFrame->IsValid() ) {
 //       cout << "non valid frame" << endl;
        return;
    }
    // the vertex used for doppler correction
    VertexInterface *vertex = 0x0;
    vertex =
            fVertexBuilder->GetVertex();
    TH1F *h;
    UInt_t nb_frames = GetNbFramePSA();

    fFold->Fill(nb_frames);

    double SumEnergy = 0.;

    for(UInt_t i = 0u; i <nb_frames; i++) {
        
        PSAInterface *data = GetDataPSA(i);
        gg[i] = data->GetE();
        h = (TH1F *)fSpectra.At(data->GetUID());
        if ( h ) {
            h->Fill(data->GetE());
        }
        fGxT->Fill(data->GetE(),data->GetT());

        SumEnergy += data->GetE();
        fCoreCommon->Fill(data->GetE());
        fDistri->Fill(data->GetUID());

        // Now look for the hit with the highest energy, taken for Doppler correction
        // Kind of tracking in crystal
        Int_t nHits = data->GetNbHits();
        Double_t max_e = 0.0;
        //test DR
        Double_t xt,yt,zt;
        vertex->GetPosition(xt,yt,zt);
        fPosXYZ_Ta->Fill(xt,yt,zt);

        for (UShort_t j = 0u; j < nHits; ++j) {

            const PSAHit* pHit = (PSAHit*)data->GetHit(j);
            //
            pHit->GetXYZ(xLocal, yLocal, zLocal);
            /*                  if ( gDebug > 3 ) {
             std::cout << "Global " << xGlobal << " " <<  yGlobal <<  " " << zGlobal << std::endl;
             std::cout << "Local " << xLocal << " " <<  yLocal <<  " " << zLocal << " " << pHit->GetID() << std::endl;
             } */
            fTrans->Local2Global(data->GetUID(), xLocal, yLocal, zLocal, xGlobal, yGlobal, zGlobal);//y,x,z
            fPosXYZ->Fill(xGlobal, yGlobal, zGlobal);
            fPosXYZ_Ta->Fill(xGlobal, yGlobal, zGlobal);
            
            if ( pHit->GetE() > max_e ) {
                gamma.SetXYZ(xGlobal, yGlobal, zGlobal);
                max_e = pHit->GetE();
            }
        } // j

        gamma.SetE(data->GetE()); h = (TH1F *)fSpectraDoppler.At(data->GetUID());
        if ( vertex ) {
            Double_t edc = DoDopplerCorrection(&gamma,vertex);
            h->Fill(edc);
            fCoreCommonDoppler->Fill(edc);
            ggdc[i] = edc;
        }
        else ggdc[i] = gg[i];
    }

    fSumSpectra->Fill(SumEnergy);

    for(UInt_t i = 0u; i <nb_frames; i++) {
        for(UInt_t j = i+1; j <nb_frames; j++) {

            fGxG->Fill(gg[i],gg[j]); fGxG->Fill(gg[j],gg[i]); fGxGDoppler->Fill(ggdc[i],ggdc[j]); fGxGDoppler->Fill(ggdc[j],ggdc[i]);
        }
    }

}


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

ClassImp(HitDisplay3D);

Float_t HitDisplay3D::fgMetric = 10.;

HitDisplay3D::HitDisplay3D(const char *name, const char *title, TDirectory *sp_dir, TDirectory *tag_dir): 
    EventPSAWatcher(name,title,sp_dir,tag_dir),
    fMaxCycle(100),
    fPathFile(Gw::AgataEventDisplay::GetDefaultAgataPath())
{
    fTrans = new AgataGeometryTransformer();
    fTrans->ReadTransformations("CrystalPositionLookUpTable.dat");
}

HitDisplay3D::~HitDisplay3D()
{
    delete fTrans;
}

void HitDisplay3D::SetActive(Bool_t active) 
{
    TTask::SetActive(active);
    AgataEventDisplay::Instance()->Reset();
} 

void HitDisplay3D::DoCanvas(TCanvas* /*c*/, Option_t* /*o*/)
{
}

void HitDisplay3D::ShowEve(const Char_t* agataPathFile)
{
    AgataEventDisplay::SetDefaultAgataPath(agataPathFile);
    AgataEventDisplay::Instance()->BuildDefaultGeometry();
    AgataEventDisplay::Instance()->ShowDisplay();
}

void HitDisplay3D::Exec(Option_t */*option*/) 
{
    // be sure the frame has been set properly
    if ( fFrame == 0x0  || !fFrame->IsValid() ) {
        return;
    }

    if ( gDebug > 3 )
        printf("HitDisplay3D::::Exec is called %d \n",GetNbFramePSA());

    AgataEventContainer* agataContainer = AgataEventDisplay::Instance()->GetEventContainer();

    if (agataContainer->GetNofEvents() > fMaxCycle) {
        TTask::SetActive(false);
        return;
    }

    static UInt_t evtNbOld = 0;
    Frame *f = GetTrigger()->GetInputFrame(0);
    const AgataKey* key = dynamic_cast<const AgataKey *>(f->GetKey());
    UInt_t evtNb = key->GetEventNumber();
    if (evtNb == evtNbOld)
        return;

    for(UInt_t i = 0u; i < GetNbFramePSA(); i++) {

        PSAInterface *data = GetDataPSA(i);

        UInt_t nHits = data->GetNbHits();
        //cout << "Fold " << nHits << endl;

        for (UShort_t j = 0u; j < nHits; ++j) {
            const PSAHit* pHit = (PSAHit*)data->GetHit(j);
            StdHit* hit = agataContainer->NewHit();
            Float_t xLocal = pHit->GetX();
            Float_t yLocal = pHit->GetY();
            Float_t zLocal = pHit->GetZ();

            Float_t xGlobal, yGlobal, zGlobal;
            fTrans->Local2Global(pHit->GetID(), xLocal, yLocal, zLocal, xGlobal, yGlobal, zGlobal);

            if ( gDebug > 3 ) {
                std::cout << "Global " << xGlobal << " " <<  yGlobal <<  " " << zGlobal << std::endl;
                std::cout << "Local " << xLocal << " " <<  yLocal <<  " " << zLocal << " " << pHit->GetE() << std::endl;
            }
            hit->SetX(xGlobal/fgMetric);
            hit->SetY(yGlobal/fgMetric);
            hit->SetZ(zGlobal/fgMetric);
            hit->SetE(pHit->GetE());
        }
    }
    agataContainer->FillHits();
    evtNbOld = evtNb;
}


//**********************************************************

ClassImp(RateDisplay);

const Double_t RateDisplay::fgkDRho   = 3.;
const Double_t RateDisplay::fgkRho1   = 5.;
const Double_t RateDisplay::fgkRho2   = 80.;
const Double_t RateDisplay::fgkDRange = 400.;
const Int_t    RateDisplay::fgkNofPoints   = 5;
const Int_t    RateDisplay::fgkNofSectors  = 6;
const Int_t    RateDisplay::fgkNofSlices   = 6;
const Int_t    RateDisplay::fgkNofCrystals = 60;
const Int_t    RateDisplay::fgkNDivX       = 5;
const Int_t    RateDisplay::fgkNDivY       = 3;

//__________________________________________________________
RateDisplay::RateDisplay(const Char_t *name, const Char_t *title, TDirectory *sp_dir, TDirectory *tag_dir) : 
    EventPSAWatcher(name,title,sp_dir,tag_dir),
    fSectorList(new TObjArray(fgkNofSectors*fgkNofCrystals)),
    fLabelList(new TObjArray(fgkNofCrystals)),
    fRateList(new TObjArray(fgkNofSectors*fgkNofCrystals)),
    fTransf(new AgataGeometryTransformer()),
    fCountCoreMax(1),
    fCountMaxAll(1),
    fCurrentCycle(0),
    fMaxCycle(200),
    fDepth(1),
    fNormByAll(true),
    fTimer(new TTimer()),
    fRefresh(3),
    fDepthText(new TLatex()),
    fBox(new TBox(0, fgkDRange, 0, fgkDRange)),
    fInfoTip(new TGToolTip(fBox, "", 350))
{

    // Set owner for delete
    fSectorList->SetOwner(true);
    fLabelList->SetOwner(true);
    fRateList->SetOwner(true);

    // count arrays
    for (Int_t i = 0; i < fgkNofSlices+1; ++i) {
        fCountArray[i] = new TArrayI(fgkNofSectors*fgkNofCrystals);
        fCountArray[i]->Reset();
    }

    fCountMax = new TArrayI(fgkNofCrystals);
    fCountMax->Reset(1);

    // create cluster & rate display
    for (Int_t j = 0; j < fgkNDivY; ++j)
        for (Int_t i = 0; i < fgkNDivX; ++i) {
            AddCluster(fgkDRange*i,fgkDRange*j);
        }

    AddRates();

    // histos
    fHistoCntCrystal  = MakeTH1<TH1I>("HistoCntCrystal","Number of counts per crystal", fgkNofCrystals, 0, fgkNofCrystals) ;

    for (Int_t i = 0; i < fgkNofCrystals; ++i) {
        for (Int_t j = 0; j < fgkNofSectors; ++j) {
            for (Int_t k = 0; k < fgkNofSlices; ++k) {
                Int_t id = i*fgkNofSectors*fgkNofSlices+j*fgkNofSlices + k;
                fHistoCntSlice[id]  = MakeTH1<TH1I>(Form("%s%d_%d_%d","fHistoCntSlice", i,j, k),
                                                    "Number of counts per sector and per slice", fMaxCycle+1, 0, fMaxCycle);
                fHistoCntSlice[id]->SetLineStyle(k+1);
                fHistoCntSlice[id]->SetLineColor(j+1);
                fHistoCntSlice[id]->SetLineWidth(2);
            }
        }
    }

    fTimer->Connect("Timeout()", "RateDisplay", this, "UpdateCanvas()");
}

//__________________________________________________________
RateDisplay::~RateDisplay()
{
    delete fSectorList;
    delete fLabelList;
    delete fRateList;
    delete fTransf;
    delete fDepthText;
    delete fTimer;
    delete fInfoTip;
    delete fBox;

    for (Int_t i = 0; i < fgkNofSlices+1; ++i)
        delete fCountArray[i];
}

//__________________________________________________________
Bool_t RateDisplay::SetTrigger( ADF::DFTrigger *trigger ) 
{
    // trigger on event:data:psa
    if ( EventPSAWatcher::SetTrigger(trigger) ) {

        // the crystal id is needed
        GObject *glob =
                GetDataPointer<PSAInterface>(fPSAFrame)->Global();
        fCrysID = glob->Get<UShort_t>("CrystalID");
        fCoreE1 = glob->Get<Float_t>("CoreE1");
        return true;
    } else return false;
}

//__________________________________________________________
void RateDisplay::Zero(Option_t *opt1, Option_t *opt2)
{
    fCurrentCycle = 0;
    return Watcher::Zero(opt1, opt2);
}

//__________________________________________________________
void RateDisplay::HandleMovement(Int_t eventType, Int_t eventX, Int_t eventY, TObject* select)
{

    TVirtualPad* pad = TVirtualPad::Pad();
    pad->cd();

    if (eventType == kMouseMotion) {
        fLastX = eventX; fLastY = eventY;

        // to keep track of the last mouse position
        fLastX = eventX; fLastY = eventY;
        if ( select->InheritsFrom("TPolyLine") && select->GetUniqueID() != 0) {
            fInfoTip->SetText(Form("Counts: %d", select->GetUniqueID()));
            fInfoTip->Show(fLastX, fLastY);
        } else
            fInfoTip->Hide();
    }

}

//__________________________________________________________
void RateDisplay::DoCanvas(TCanvas* c, Option_t* o)
{

    TString tmp(o);
    tmp.ToLower();

    // rate display
    if (tmp.IsNull() ) {

        if ( !c->HasConnection("ProcessedEvent(Int_t, Int_t, Int_t, TObject*)") ) {
            c->Connect("ProcessedEvent(Int_t, Int_t, Int_t, TObject*)",
                       "RateDisplay",
                       this,
                       "HandleMovement(Int_t, Int_t, Int_t, TObject*)");
        }

        gPad->Range(-fgkDRange/2, -fgkDRange/2, fgkDRange/2+(fgkNDivX-1)*fgkDRange, fgkDRange/2+(fgkNDivY-1)*fgkDRange);

        fSectorList->Draw();
        fLabelList->Draw();
        fRateList->Draw();
        fDepthText->Draw();

        gPad->SetEditable(false);
    }

    // histo display
    if (tmp.BeginsWith("h") ) {
        c->SetWindowSize(fgkDRange*fgkNDivX,fgkDRange*fgkNDivY);
        c->Divide(fgkNDivX, fgkNDivY*3);

        for (Int_t i = 0; i < fgkNofCrystals; ++i) {
            c->cd(i+1);
            for (Int_t j = 0; j < fgkNofSectors; ++j) {
                for (Int_t k = 0; k < fgkNofSlices; ++k) {
                    Int_t id = i*fgkNofSectors*fgkNofSlices+j*fgkNofSlices + k;
                    if (j == 0 && k == 0)
                        fHistoCntSlice[id]->Draw();
                    else
                        fHistoCntSlice[id]->Draw("SAME");
                }
            }
        }
    }
    SetDepth(0, false);
    fTimer->Start(fRefresh*1000);
}

//__________________________________________________________
void RateDisplay::UpdateCanvas()
{
    static Int_t depth = 0;

    SetDepth(depth++, false);
    if (depth == fgkNofSlices+1) {
        depth = 0;
        Reset();
    }

    //  if (gPad) { // TO BE DONE
    //    gPad->Update();
    //  }
}

//__________________________________________________________
void RateDisplay::Exec(Option_t */*option*/) 
{
    // Frame by Frame action

    // be sure the frame has been set properly
    if ( fFrame == 0x0 || !fFrame->IsValid() )
        return;

    if ( gDebug > 3 )
        printf("RateDisplay::Exec is called %d \n",GetNbFramePSA());


    Int_t tmpID = 0;
    Int_t tmp   = 0;

    // depth 0 for core energy
    for (Int_t s = 0; s < fgkNofSectors; ++s) {
        tmpID = fCrysID->Get()*fgkNofSectors + s;
        tmp   = fCountArray[0]->At(tmpID)+1;
        if (fCoreE1->Get() > 0) {
            fCountArray[0]->AddAt(tmp, tmpID);
            if (tmp > fCountCoreMax)
                fCountCoreMax = tmp;
        }
    }
    fHistoCntCrystal->Fill(fCrysID->Get());

    // loop on PSAFrame and extract the information you need
    for(UInt_t i = 0u; i < GetNbFramePSA(); ++i) {
        PSAInterface *data = GetDataPSA(i);

        if ( data ) {
            for (Int_t j = 0; j < data->GetNbHits(); ++j) {
                const PSAHit* pHit = (PSAHit*)data->GetHit(j);
                Int_t depthID  = fTransf->GetDepthId(pHit->GetZ());
                Int_t sectorID = fTransf->GetSectorId(pHit->GetX(), pHit->GetY());

                tmpID = fCrysID->Get()*fgkNofSectors + sectorID;
                tmp   = fCountArray[depthID]->At(tmpID)+1;

                Int_t histoID = fCrysID->Get()*fgkNofSectors*fgkNofSlices + sectorID*fgkNofSlices + (depthID-1);
                fHistoCntSlice[histoID]->Fill(fCurrentCycle, tmp);

                fCountArray[depthID]->AddAt(tmp, tmpID);

                if ( tmp > fCountMax->At(fCrysID->Get()) )
                    fCountMax->AddAt(tmp, fCrysID->Get());

                if (tmp > fCountMaxAll)
                    fCountMaxAll = tmp;


                if ( gDebug > 3 ) {
                    printf("hit %d with pos (%f, %f, %f)\n", j, pHit->GetX(), pHit->GetY(), pHit->GetZ());
                    printf("Count %d for crystal %d depth %d and sector %d max %d\n",
                           fCountArray[depthID]->At(tmpID), fCrysID->Get(), depthID, sectorID, fCountMax->At(fCrysID->Get()));
                }
            }
        }
    }
}

//__________________________________________________________
void RateDisplay::AddCluster(Double_t xOffset, Double_t yOffset)
{
    TPaveText* pave = 0x0;
    Double_t xPave  = 0;
    Double_t yPave  = 0;
    Double_t dxPave = 10;
    Double_t dyPave = 10;

    Double_t xPos = 0;
    Double_t yPos = 2*TMath::Sin(TMath::Pi()/6)*(fgkRho2+3*fgkDRho);
    AddCristals(xPos+xOffset,yPos+yOffset);

    xPave = xPos+xOffset-0.8*fgkRho2;
    yPave = yPos+yOffset+fgkRho2;
    pave  = new TPaveText(xPave-dxPave,yPave-dyPave, xPave+dxPave, yPave+dyPave);
    pave->SetOption("");
    pave->AddText(Form("%d", (fSectorList->GetEntries()-1)/fgkNofSectors));
    fLabelList->Add(pave);

    xPos = -TMath::Cos(TMath::Pi()/6)*(fgkRho2+3*fgkDRho);
    yPos = -TMath::Sin(TMath::Pi()/6)*(fgkRho2+3*fgkDRho);
    AddCristals(xPos+xOffset, yPos+yOffset);

    xPave = xPos+xOffset-0.8*fgkRho2;
    yPave = yPos+yOffset+fgkRho2;
    pave  = new TPaveText(xPave-dxPave,yPave-dyPave, xPave+dxPave, yPave+dyPave);
    pave->SetOption("l");
    pave->AddText(Form("%d", (fSectorList->GetEntries()-1)/fgkNofSectors));
    fLabelList->Add(pave);

    xPos =  TMath::Cos(TMath::Pi()/6)*(fgkRho2+3*fgkDRho);
    yPos = -TMath::Sin(TMath::Pi()/6)*(fgkRho2+3*fgkDRho);
    AddCristals(xPos+xOffset, yPos+yOffset);

    xPave = xPos+xOffset+0.8*fgkRho2;
    yPave = yPos+yOffset+fgkRho2;
    pave  = new TPaveText(xPave-dxPave,yPave-dyPave, xPave+dxPave, yPave+dyPave);
    pave->SetOption("r");
    pave->AddText(Form("%d", (fSectorList->GetEntries()-1)/fgkNofSectors));
    fLabelList->Add(pave);
}

//__________________________________________________________
void RateDisplay::AddCristals(Double_t xPos, Double_t yPos)
{  
    TPolyLine* detectorLine = 0x0;

    Double_t x[fgkNofPoints];
    Double_t y[fgkNofPoints];

    Double_t offset = TMath::Pi()/6.;
    Double_t dPhi   = TMath::Pi()/3.;

    for (Int_t i = 0; i < fgkNofSectors; ++i) {
        Double_t angle = dPhi*i + offset;
        Double_t shiftX = fgkDRho*TMath::Cos(angle+dPhi/2.);
        Double_t shiftY = fgkDRho*TMath::Sin(angle+dPhi/2.);

        x[0] = x[4] = fgkRho1*TMath::Cos(angle) + shiftX;
        y[0] = y[4] = fgkRho1*TMath::Sin(angle) + shiftY;
        x[1] = fgkRho2*TMath::Cos(angle) + shiftX;
        y[1] = fgkRho2*TMath::Sin(angle) + shiftY;

        angle = dPhi*(i+1) + offset;
        x[2] = fgkRho2*TMath::Cos(angle) + shiftX;
        y[2] = fgkRho2*TMath::Sin(angle) + shiftY;

        x[3] = fgkRho1*TMath::Cos(angle) + shiftX;
        y[3] = fgkRho1*TMath::Sin(angle) + shiftY;

        for (Int_t j = 0; j  < fgkNofPoints; ++j) {
            x[j] += xPos;
            y[j] += yPos;
        }

        detectorLine = new TPolyLine(fgkNofPoints,x,y);
        detectorLine->SetLineWidth(1);
        detectorLine->SetLineStyle(2);
        fSectorList->Add(detectorLine);
    }
}

//__________________________________________________________
void RateDisplay::AddRates()
{
    TPolyLine* rateLine = 0x0;

    for (Int_t iCrystal = 0; iCrystal < fSectorList->GetEntries(); ++iCrystal) {
        Color_t col = 0;
        if (iCrystal/fgkNofSectors % 3 == 0)
            col = kRed;
        if (iCrystal/fgkNofSectors % 3 == 1)
            col = kGreen;
        if (iCrystal/fgkNofSectors % 3 == 2)
            col = kBlue;

        rateLine = new TPolyLine(fgkNofPoints);
        rateLine->SetLineColor(2);
        rateLine->SetLineWidth(1);
        rateLine->SetFillColor(col);
        rateLine->SetOption("f");
        fRateList->Add(rateLine);

    }
}

//__________________________________________________________
void RateDisplay::SetRefresh(Long_t r) 
{ 
    fRefresh = r;
    fTimer->Start(r*1000);
}

//__________________________________________________________
void RateDisplay::SetDepth(Int_t d, Bool_t stop) 
{ 
    if (stop)
        fTimer->Stop();

    if (d > 6 || d < 0) return;

    fDepth = d;
    UpdateRates();
} 

//__________________________________________________________
void RateDisplay::UpdateRates()
{
    TPolyLine* rateLine     = 0x0;
    TPolyLine* detectorLine = 0x0;

    Double_t xm;
    Double_t ym;

    Double_t x[fgkNofPoints];
    Double_t y[fgkNofPoints];

    for (Int_t iCrystal = 0; iCrystal < fSectorList->GetEntries()/fgkNofSectors; ++iCrystal) {

        Int_t max = 1;

        if (IsNormByAll()) {
            max = fCountMaxAll;
        } else {
            max = fCountMax->At(iCrystal);
        }

        if (max == 0) continue;

        for (Int_t iSector = 0; iSector < fgkNofSectors; ++iSector) {
            detectorLine = static_cast<TPolyLine*> ( fSectorList->At(iCrystal*fgkNofSectors + iSector) );
            rateLine     = static_cast<TPolyLine*> ( fRateList->At(iCrystal*fgkNofSectors + iSector) );

            Double_t x1 = (detectorLine->GetX()[1] + detectorLine->GetX()[2])/2.;
            Double_t y1 = (detectorLine->GetY()[1] + detectorLine->GetY()[2])/2.;

            Double_t x2 = (detectorLine->GetX()[0] + detectorLine->GetX()[3])/2.;
            Double_t y2 = (detectorLine->GetY()[0] + detectorLine->GetY()[3])/2.;

            xm = (x1+x2)/2.;
            ym = (y1+y2)/2.;

            Double_t rateRatio = 1.;

            if (fDepth != 0)
                rateRatio = (Double_t)fCountArray[fDepth]->At(iCrystal*fgkNofSectors + iSector)/max*0.85;
            else
                rateRatio = (Double_t)fCountArray[fDepth]->At(iCrystal*fgkNofSectors + iSector)/fCountCoreMax*0.85;

            if ( gDebug > 3 ) {
                printf("count %d ratio %lf for sector %d and crystal %d\n",
                       fCountArray[fDepth]->At(iCrystal*fgkNofSectors + iSector), rateRatio,
                       iSector, iCrystal);
            }

            for (Int_t j = 0; j < fgkNofPoints; ++j) {
                x[j] = (detectorLine->GetX()[j] - xm)* rateRatio + xm;
                y[j] = (detectorLine->GetY()[j] - ym)* rateRatio + ym;
            }

            rateLine->SetPolyLine(fgkNofPoints, x, y, "f");
            rateLine->SetUniqueID(fCountArray[fDepth]->At(iCrystal*fgkNofSectors + iSector));
        }
    }

    // show depth
    fDepthText->Clear();
    fDepthText->SetTextSize(0.04);
    fDepthText->SetText(-fgkDRange/2+20, -fgkDRange/2+20, Form("Depth: %d", fDepth));

    fCurrentCycle++;
}

//__________________________________________________________
void RateDisplay::Reset()
{
    // reset counters + cycle
    fCountMax->Reset();
    fCountCoreMax = 1;
    for (Int_t i = 0; i < fgkNofSlices+1; ++i)
        fCountArray[i]->Reset();

    fCountMaxAll = 0;
}