StogsAFP.cpp

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2004 by Olivier Stezowski                               *
 *   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.             *
 ***************************************************************************/

#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <cstdio>

#include <list>
#include <utility>
#include <vector>

#include "TString.h"
#include "TObjArray.h"
#include "TObject.h"
#include "TSystem.h"
#include "TList.h"
#include "CompositeFrame.h"
#include "TFile.h"
#include "TChain.h"
#include "TEntryList.h"
#include "TKey.h"
#include "TChainElement.h"
#include "TError.h"
#include "TVector3.h"
#include "TRandom3.h"
#include "TROOT.h"

#include "PSAFrame.h"
#include "BashColor.h"
#include "AgataFrameFactory.h"
#include "AgataKeyFactory.h"
#include "BufferIO.h"
#include "SToGS_BaseROOTEvents.h"

#include "StogsAFP.h"


StogsAFP::StogsAFP():
    NarvalProducer(),
    fSourceOfFrames(aMByte,ConfAgent::kRead),
    fEndOfFrame("EndOfFrame"),
    fFrameEventPsa_Out(NULL),
    fSPSAFrame(NULL),
    fTrigger("event:data:psa"),
    fdebug(false),
    fBashColor(new BashColor()),
    fNEvts(0),
    fRejectedEvents(0),
    fEventNumber(0),
    fNEntries(0),
    fDataChain(new TChain("TreeMaster","TreeMaster")),
    fEntryList(0x0),
    fCurrentEntry(0),
    IsADHitTreeModeOn(false),
    fStogsEv(new SBREvent()),
    fProcessPacking(true),
    fSpatialResolution(5.), //default spatial resolution in mm
    fApplyEnergyResolution(true),
    fApplySmearing(true)
{       
    if(fdebug) std::cout<<"\E[31;1m"<<"Construct in"<<"\E[m"<<std::endl;

    Log.GetProcessName() = "StogsAFP"; GetFrameIO().SetName("StogsAFP");

    // to be sure unknown frames are written in the ouput
    GetFrameIO().SetModel(ConfAgent::kSafe);
    // to protect writings in the source of frames
    fSourceOfFrames.SetModeIO(ConfAgent::kRead);

    // the input datablock is owned by BasicAFP
    GetFrameIO().Attach(&fSourceOfFrames,0x0);

    gRandom->SetSeed();

    if(fdebug) std::cout<<"\E[31;1m"<<"Construct out"<<"\E[m"<<std::endl;
}


StogsAFP::~StogsAFP()
{  
    // in principle not needed ... just in case reset it has not been called by narval
    unsigned int error = 0u; process_reset(&error) ;
}


void StogsAFP::process_config (const char *directory_path, unsigned int *error_code)
{       
    // first init narval and ADF stuff (if required)
    NarvalInterface::process_config(directory_path,error_code);

    if ( (*error_code) == 0u ) { // now init your stuff
    }
}


unsigned int StogsAFP::ProcessBlock(ADF::FrameBlock &out)
{ 
    if(fdebug) std::cout<<"\E[31;1m"<<"ProcessBlock in"<<"\E[m"<<std::endl;

    // attach the input/output buffer to the FrameIO system

    GetFrameIO().Attach(NULL, &out);

    // start the processing

    while(true)
    {
        bool EventOk;
        if(IsADHitTreeModeOn) EventOk = ReadEventADHitTree();
        else EventOk = ReadEventStogs();

        if(!EventOk)
        {
            if(fCurrentEntry==fNEntries)
            {
                GetFrameIO().SetStatus(BaseFrameIO::kFinished);
                out.SetEoB(true);
                break;
            }
            else continue;
        }

        CompositeFrame *frame_out = (CompositeFrame*)fFrameEventPsa_Out->GetFrame();
        unsigned int EventLenght = frame_out->GetLength();

        if(fdebug)
        {
            fBashColor->SetErrorOut();
            cout<<"EventLenght = "<<EventLenght<<endl;
            cout<<out.GetBufferIO()->FreeSize()<<endl;
            fBashColor->ResetColor();
        }

        if(out.Reserve(EventLenght,true))
        {
            fNEvts++;
            fTrigger.Fired(true);
            GetFrameIO().Record();
            if(fCurrentEntry==fNEntries)
            {
                GetFrameIO().SetStatus(BaseFrameIO::kFinished);
                out.SetEoB(true);
                break;
            }
            else continue;
        }
        else
        {
            fCurrentEntry--;
            fEventNumber--;
            break;
        }
    }

    GetFrameIO().Detach(NULL, &out);

    if(fdebug) std::cout<<"\E[31;1m"<<"ProcessBlock out"<<"\E[m"<<std::endl;
    if ( GetFrameIO().GetStatus() == BaseFrameIO::kFinished) return 1u;
    else return 0u;
}       


void StogsAFP::process_initialise (unsigned int *error_code)
{
    if(fdebug) std::cout<<"\E[31;1m"<<"process_initialise in"<<"\E[m"<<std::endl;

    *error_code = 0u; Log.ClearMessage(); Log.SetProcessMethod("process_initialise");

    // read an input file to get the base filename and the path
    std::string conffile = GetConfPath();
    conffile += "StogsAFP.conf";

    ReadConfFile(conffile.data());

    if(fEntryList == 0x0) fNEntries = fDataChain->GetEntries();
    else fNEntries = fEntryList->GetN();

    if(fNEntries)
    {
        fBashColor->SetInfoOut();
        std::cout<<"TChain created: "<<fDataChain->GetName()<<std::endl;
        std::cout<<"Number of entry to Treat: "<<fNEntries<<std::endl;
        fBashColor->ResetColor();

        if(IsADHitTreeModeOn)InitDataChainADHitTree();
        else InitDataChainStogs();
    }
    else
    {
        fBashColor->SetErrorOut();
        std::cout<<"Empty TChain: "<<fDataChain->GetName()<<std::endl;
        fBashColor->ResetColor();
        gROOT->ProcessLine(".qqqqqqqqqqqqqqqq");
    }

    //    FactoryItem keyItem("Agata", "data:psa",ADF::Version(4,0));
    //    FactoryItem fraItem("Agata", "data:psa",ADF::Version(0,0));
    //    fSPSAFrame = MainFrameFactory::theMainFactory().NewSharedFrame(keyItem, fraItem);

    fFrameEventPsa_Out = fTrigger.SetOutputFrame("Agata", "event:data:psa");

    fSPSAFrame = fTrigger.AddUtility("data:psa",ADF::ConfAgent::theGlobalAgent());

    GetFrameIO().Register(&fTrigger);

    Log << dolog;

    if(fdebug) std::cout<<"\E[31;1m"<<"process_initialise out"<<"\E[m"<<std::endl;
}


void StogsAFP::process_reset (unsigned int *error_code)
{
    *error_code = 0u;
    Log.ClearMessage(); Log.SetProcessMethod("process_reset");

    GetFrameIO().Print( Log() );

    Log << dolog;
}


bool StogsAFP::ReadConfFile(TString ConfFilePath)
{
    TList *ListOfRootFiles;
    ListOfRootFiles = new TList();

    TString TreeName = "TreeMaster";

    fCurrentEntry = 0;

    std::ifstream filein(ConfFilePath);
    if ( filein.is_open())
    {
        std::string Line;
        TString Buffer;

        while ( filein.good() && !filein.eof() )
        {
            getline(filein,Line);
            Buffer = Line;

            if(Buffer.BeginsWith("#")) continue;
            if(Buffer.BeginsWith("b") || Buffer.BeginsWith("r")) //b /path/to PSAHits 0 //r /path/to PSAHits 0 99
            {
                TObjArray *TokenArray = Buffer.Tokenize(" ");

                if(Buffer.BeginsWith("b") && TokenArray->GetEntries() != 4) continue;
                if(Buffer.BeginsWith("r") && TokenArray->GetEntries() != 5) continue;

                TString Path = TokenArray->At(1)->GetName();
                if(!Path.EndsWith("/")) Path.Append("/");
                TString BaseName = TokenArray->At(2)->GetName();
                int FirstFileNbr = ((TString)TokenArray->At(3)->GetName()).Atoi();
                int LastFileNbr = FirstFileNbr+1;
                if(Buffer.BeginsWith("r")) LastFileNbr = ((TString)TokenArray->At(4)->GetName()).Atoi();

                TString FileName = Form("%s%s_%.4d.root",Path.Data(),BaseName.Data(),FirstFileNbr);

                int NFilesAdded = 0;
                while(gSystem->IsFileInIncludePath(FileName) && FirstFileNbr<=LastFileNbr)
                {
                    TFile *f = TFile::Open(FileName);
                    ListOfRootFiles->Add(f);
                    FileName = Form("%s%s_%.4d.root",Path.Data(),BaseName.Data(),++FirstFileNbr);

                    if(Buffer.BeginsWith("b")) LastFileNbr++;
                }
                if(!NFilesAdded)
                {
                    fBashColor->SetWarningOut();
                    cout<<"File "<<FileName<<" not found !"<<endl;
                    fBashColor->ResetColor();
                }
            }
            if(Buffer.BeginsWith("l"))
            {
                TObjArray *TokenArray = Buffer.Tokenize(" ");

                if(TokenArray->GetEntries() != 2) continue;

                TString FileName = TokenArray->Last()->GetName();

                TFile *f = TFile::Open(FileName);
                if(gSystem->IsFileInIncludePath(FileName)) ListOfRootFiles->Add(f);
                else
                {
                    fBashColor->SetWarningOut();
                    cout<<"File "<<FileName<<" not found !"<<endl;
                    fBashColor->ResetColor();
                }
            }
            if(Buffer.BeginsWith("el"))
            {
                TObjArray *TokenArray = Buffer.Tokenize(" ");
                int TokenSize = TokenArray->GetEntries();
                if(TokenSize != 2 && TokenSize != 3) continue;

                TString EntryListName = "";
                TEntryList *EntryList = 0x0;

                TString FileName = ((TString)TokenArray->At(1)->GetName());
                TFile *f = TFile::Open(FileName);
                if(!gSystem->IsFileInIncludePath(FileName))
                {
                    fBashColor->SetWarningOut();
                    cout<<"File "<<FileName<<" not found !"<<endl;
                    fBashColor->ResetColor();
                    continue;
                }

                if(TokenSize==3)
                {
                    EntryListName = TokenArray->Last()->GetName();
                    EntryList = (TEntryList*)f->Get(EntryListName)->Clone();
                    if(EntryList==0x0)
                    {
                        fBashColor->SetWarningOut();
                        cout<<"EntryList "<<EntryListName<<" not found in file "<<FileName<<"!"<<endl;
                        fBashColor->ResetColor();
                        continue;
                    }
                }
                else
                {
                    TList *ListOfKeys = f->GetListOfKeys();

                    for(int i=0 ; i<ListOfKeys->GetSize() ; i++)
                    {
                        TKey *key = (TKey*)ListOfKeys->At(i);

                        if((TString)key->GetClassName()=="TEntryList")
                        {
                            EntryList = (TEntryList*)f->Get(key->GetName())->Clone();
                        }
                    }
                    if(EntryList==0x0)
                    {
                        fBashColor->SetWarningOut();
                        cout<<"No EntryList found in file "<<FileName<<"!"<<endl;
                        fBashColor->ResetColor();
                        continue;
                    }
                }

                if(fEntryList == 0x0) fEntryList = EntryList;
                else fEntryList->Add(EntryList);

                EntryList->SetDirectory(0);
                f->Close();

                fBashColor->SetInfoOut();
                cout<<"EntryList "<<EntryList->GetName()<<" loaded for file "<<EntryList->GetFileName()<<endl;
                fBashColor->ResetColor();
            }
            if(Buffer.BeginsWith("n"))
            {
                TObjArray *TokenArray = Buffer.Tokenize(" ");
                if(TokenArray->GetEntries() != 2)
                {
                    fBashColor->SetWarningOut();
                    std::cout<<"The TTree name has not been undestood..."<<std::endl;
                    fBashColor->ResetColor();
                    continue;
                }
                TreeName = TokenArray->Last()->GetName();
                fDataChain->SetName(TreeName);
            }
        }
    }
    else
    {
        fBashColor->SetWarningOut();
        std::cout<<ConfFilePath<<" not found !"<<std::endl;
        gROOT->ProcessLine(".qqqqqqqqqqqqqqqqq");
        fBashColor->ResetColor();
    }
    filein.close();

    for(int i=0 ; i<ListOfRootFiles->GetSize() ; i++)
    {
        TFile *f = (TFile*)ListOfRootFiles->At(i);

        if(f!= 0x0) fDataChain->Add(f->GetName());
    }

    gErrorIgnoreLevel = kFatal;

    if(fEntryList != 0x0)
    {
        TObjArray *fFiles = fDataChain->GetListOfFiles();
        Int_t ne = fFiles->GetEntries();
        Int_t listfound=0;
        TString treename, filename;

        TEntryList *templist = 0;
        for (Int_t ie = 0; ie<ne; ie++){
           treename = gSystem->BaseName( ((TChainElement*)fFiles->UncheckedAt(ie))->GetName() );
           filename = ((TChainElement*)fFiles->UncheckedAt(ie))->GetTitle();
           templist = fEntryList->GetEntryList(treename.Data(), filename.Data());
           if(!templist)
           {
               TString BaseFilename = gSystem->BaseName(filename);
               templist = fEntryList->GetEntryList(treename.Data(), BaseFilename.Data());
               if(templist) templist->SetFileName(filename);
           }
           if (templist) listfound++;
        }

        if(!listfound)
        {
            fBashColor->SetWarningOut();
            cout<<"The selected EntryLists are not compatible with input TTrees, EntryList ignored."<<endl;
            fBashColor->ResetColor();

            delete fEntryList;
            fEntryList = 0x0;
        }
        else fDataChain->SetEntryList(fEntryList,"ne");


    }

    gErrorIgnoreLevel = kPrint;

    return true;
}


void StogsAFP::InitDataChainADHitTree()
{
    fDataChain->SetBranchAddress("nbHits", &NbHits);
    fDataChain->SetBranchAddress("hitE", HitE);
    fDataChain->SetBranchAddress("hitX", HitX);
    fDataChain->SetBranchAddress("hitY", HitY);
    fDataChain->SetBranchAddress("hitZ", HitZ);

    fDataChain->SetBranchAddress("hitId", HitId);
    fDataChain->SetBranchAddress("hitSg", HitSg);


    fDataChain->SetBranchAddress("nbCores", &NbCores);
    fDataChain->SetBranchAddress("coreE0", CoreE0);
    fDataChain->SetBranchAddress("coreE1", CoreE1);
    fDataChain->SetBranchAddress("coreT0", CoreT0);
    fDataChain->SetBranchAddress("coreT1", CoreT1);
    fDataChain->SetBranchAddress("hit_per_cryst", NbHitsPerCrystal);

    fDataChain->SetBranchAddress("TSHit", &fTimestamp);

    NbHits=0;
    NbCores=0;
    for(int i=0 ; i<MaxHits ; i++)
    {
        HitE[i] = 0;
        HitX[i] = 0;
        HitY[i] = 0;
        HitZ[i] = 0;
        HitId[i] = 0;
        HitSg[i] = 0;

        if(i<MaxCore)
        {
            CoreE0[i]=0;
            CoreE1[i]=0;
            CoreT0[i]=0;
            CoreT1[i]=0;
            NbHitsPerCrystal[i]=0;
        }

    }
}


void StogsAFP::InitDataChainStogs()
{
    fDataChain->SetBranchAddress("Ev.", &fStogsEv);

    fEventNumber = 0;

    fStogsEv->Clear("");
}


bool StogsAFP::ReadEventADHitTree()
{

    int EntrySize = fDataChain->GetEntry(fCurrentEntry++);

    if(EntrySize==0)
    {
        fBashColor->SetWarningOut();
        cout<<"Warning ! Entry "<<fCurrentEntry<<" Is empty, is it normal ?"<<endl;
        fBashColor->ResetColor();
        return false;
    }

    CompositeFrame *frame_out = (CompositeFrame*)fFrameEventPsa_Out->GetFrame();
    frame_out->Reset();

    ((AgataKey *)frame_out->GetKey())->SetEventNumber(fCurrentEntry);
    ((AgataKey *)frame_out->GetKey())->SetTimeStamp(fTimestamp);

    if(fdebug)
    {
        TString GlobMessType = AgataKeyFactory::theFactory()->GetMessage(((AgataKey*)frame_out->GetKey())->GetMessage());

        fBashColor->SetInfoOut();
        cout<<"New Main frame: "<<GlobMessType<<endl;
        cout<<"  --timestamp = "<<((AgataKey*)frame_out->GetKey())->GetTimeStamp()<<endl;
        cout<<"  --EventNumber = "<<((AgataKey*)frame_out->GetKey())->GetEventNumber()<<endl;
        cout<<"  --Number of total hits: "<<NbHits<<" in "<<NbCores<<" core"<<endl;
        fBashColor->ResetColor();
    }

    int hit = 0;
    for(int icore=0 ; icore< NbCores ; icore++)
    {
        Frame *PSAFrame = fSPSAFrame->GetFrame();
        PSAFrame->Reset();

        ((AgataKey *)PSAFrame->GetKey())->SetEventNumber(fCurrentEntry);
        ((AgataKey *)PSAFrame->GetKey())->SetTimeStamp(fTimestamp);

        PSAInterface *data_out  = GetDataPointer<PSAInterface>(PSAFrame);
        data_out->SetE(CoreE0[icore],0);
        data_out->SetE(CoreE1[icore],1);
        data_out->SetT(CoreT0[icore],0);
        data_out->SetT(CoreT1[icore],1);
        data_out->SetUID(HitId[hit]);

        if(fdebug)
        {
            TString MessageType = AgataKeyFactory::theFactory()->GetMessage(((AgataKey*)PSAFrame->GetKey())->GetMessage());
            fBashColor->SetInfoOut();
            cout<<"    New subframe: "<<MessageType<<endl;
            cout<<"      --CoreID: "<<HitId[hit]<<" ==> E0 = "<<CoreE0[icore]<<", T0 = "<<CoreT0[icore]<<endl;
            cout<<"      --Number of crystal hits: "<<NbHitsPerCrystal[icore]<<endl;
            fBashColor->ResetColor();
        }

        ADF::PSAHit *pOut = NULL;
        for(int iseg=0 ; iseg<NbHitsPerCrystal[icore] ; iseg++)
        {
            pOut = (PSAHit*)data_out->NewHit();
            pOut->Reset();
            pOut->SetXYZ(HitX[hit],HitY[hit],HitZ[hit]);
            pOut->SetE(HitE[hit]);
            pOut->SetID(HitSg[hit],0);
            pOut->SetID(HitId[hit],1);

            if(fdebug)
            {
                fBashColor->SetInfoOut();
                cout<<"        hit nb: "<<iseg<<" : E = "<<HitE[hit]<<", SegID = "<<HitId[hit]<<", SegNbr = "<<HitSg[hit]<<", (X,Y,Z): ("<<HitX[hit]<<","<<HitY[hit]<<","<<HitZ[hit]<<")"<<endl;
                fBashColor->ResetColor();
            }
            hit++;
        }

        PSAFrame->Write();

        frame_out->AddSubFrame((*PSAFrame));
    }

    if(fdebug)
    {
        fBashColor->SetInfoOut();
        cout<<"  ==> Number of subframes in the global frame: "<<frame_out->GetNbSubFrame()<<endl;
        cout<<"  ==> Event Lenght = "<<frame_out->GetLength()<<endl;
        cout<<"_____________________________________________"<<endl<<endl;
        fBashColor->ResetColor();
    }

    return true;
}


bool StogsAFP::ReadEventStogs()
{
    fStogsEv->Clear("");

    int EntrySize;

    if(fEntryList == 0x0) EntrySize = fDataChain->GetEntry(fCurrentEntry++);
    else if(fDataChain->InheritsFrom("TChain")) EntrySize = fDataChain->GetEntry(fEntryList->GetEntry(fCurrentEntry++));
    else EntrySize = fDataChain->GetEntry(fDataChain->GetEntryNumber(fCurrentEntry++));


    if(EntrySize==0)
    {
        fBashColor->SetWarningOut();
        cout<<"Warning ! Entry "<<fCurrentEntry -1 <<" Is empty, is it normal ?"<<endl;
        fBashColor->ResetColor();
        TString test;
        cin>>test;
        if(!test.EqualTo("y",TString::kIgnoreCase)) return false;
    }

    SBREvent *FilteredEvent = FilterEvent(fStogsEv);
    SBREvent *PackingEvt = FilteredEvent;

    if(fProcessPacking) PackingEvt = ProcessPacking(FilteredEvent);
    if(fApplyEnergyResolution) ApplyEnergyResolution(PackingEvt);
    if(fApplySmearing) ApplySmearing(PackingEvt);

    CompositeFrame *frame_out = (CompositeFrame*)fFrameEventPsa_Out->GetFrame();
    frame_out->Reset();

    fTimestamp = fCurrentEntry*1e8/500.; //pour avoir un rate de 500 coups/s
    ((AgataKey *)frame_out->GetKey())->SetEventNumber(fEventNumber);
    ((AgataKey *)frame_out->GetKey())->SetTimeStamp(fTimestamp);

    int NbHits = PackingEvt->GetNbHits();

    int nbCores=0;

    if(fdebug)
    {
        fBashColor->SetInfoOut();
        cout<<"*** New Event : NbHits = "<<NbHits<<endl;
        fBashColor->ResetColor();
    }

    int hit=0;

    while(hit<NbHits)
    {
        Frame *PSAFrame = fSPSAFrame->GetFrame();
        PSAFrame->Reset();

        ((AgataKey *)PSAFrame->GetKey())->SetEventNumber(fEventNumber);
        ((AgataKey *)PSAFrame->GetKey())->SetTimeStamp(fTimestamp);

        PSAInterface *data_out  = GetDataPointer<PSAInterface>(PSAFrame);

        SBRHit *a_hit = PackingEvt->GetHit(hit);

        int CoreId = a_hit->fUID;
        if(CoreId>180)
        {
            hit++;
            continue;
        }
        double CoreE = 0.;

        int nHitsPerCrystal = 0;

        ADF::PSAHit *pOut = NULL;
        while(( hit<NbHits ) && ( a_hit->fUID == CoreId ))
        {
            if(nHitsPerCrystal ==100) break;
            pOut = (PSAHit*)data_out->NewHit();
            pOut->Reset();
            pOut->SetXYZ(a_hit->fX,a_hit->fY,a_hit->fZ); // Warning !!! for the moment, the space variables are given in the global frame, not in the crystal frame
            double ahitE = a_hit->fE;
            pOut->SetE(ahitE);
            CoreE += ahitE;
            pOut->SetID(0,0); // Warning, without access to the lookuptable, the seg Id is cannot be deduced from the global coordinates
            pOut->SetID(CoreId,1);

            if(fdebug || ahitE<=0)
            {
                fBashColor->SetInfoOut();
                if(ahitE<=0) fBashColor->SetWarningOut();
                cout<<"        hit nb: "<<nHitsPerCrystal<<" : E = "<<a_hit->fE<<", SegID = "<<a_hit->fUID<<", (X,Y,Z): ("<<a_hit->fX<<","<<a_hit->fY<<","<<a_hit->fZ<<")"<<endl;
                fBashColor->ResetColor();
            }

            hit++;
            nHitsPerCrystal++;

            a_hit = PackingEvt->GetHit(hit);
        }

        data_out->SetE(CoreE,0);
        data_out->SetE(0,1);
        data_out->SetT(0,0);
        data_out->SetT(0,1);
        data_out->SetUID(CoreId);

        nbCores++;

        int FrameSize = PSAFrame->Write();
        if(FrameSize == 0)
        {
            fRejectedEvents++;
            return false;
        }

        frame_out->AddSubFrame((*PSAFrame));

        if(fdebug)
        {
            TString MessageType = AgataKeyFactory::theFactory()->GetMessage(((AgataKey*)PSAFrame->GetKey())->GetMessage());
            fBashColor->SetInfoOut();
            cout<<"    New subframe: "<<MessageType<<endl;
            cout<<"      --CoreID: "<<CoreId<<" ==> E0 = "<<CoreE<<endl;
            cout<<"      --Number of crystal hits: "<<nHitsPerCrystal<<endl;
            fBashColor->ResetColor();
        }
    }

    if(fProcessPacking) delete PackingEvt;
    delete FilteredEvent;

    if(fdebug)
    {
        TString GlobMessType = AgataKeyFactory::theFactory()->GetMessage(((AgataKey*)frame_out->GetKey())->GetMessage());

        fBashColor->SetInfoOut();
        cout<<"New Main frame: "<<GlobMessType<<endl;
        cout<<"  --timestamp = "<<((AgataKey*)frame_out->GetKey())->GetTimeStamp()<<endl;
        cout<<"  --EventNumber = "<<((AgataKey*)frame_out->GetKey())->GetEventNumber()<<endl;
        cout<<"  --Number of total hits: "<<NbHits<<" in "<<nbCores<<" core"<<endl;

        cout<<"  ==> Number of subframes in the global frame: "<<frame_out->GetNbSubFrame()<<endl;
        cout<<"  ==> Event Lenght = "<<frame_out->GetLength()<<endl;
        cout<<"_____________________________________________"<<endl<<endl;
        fBashColor->ResetColor();
    }

    if(NbHits==0) return false;

    fEventNumber++;

    return true;
}

void StogsAFP::Print()
{
    fBashColor->SetInfoOut();

    std::cout<<"********  StogsAFP producer stats:  ********"<<std::endl<<std::endl;

    std::cout<<"Number of entries in the input tree : "<<fNEntries<<std::endl;
    std::cout<<"Number of good events               : "<<fNEvts<<std::endl;
    std::cout<<"Number of empty events              : "<<fNEntries-fNEvts-fRejectedEvents<<std::endl;
    std::cout<<"Number of rejected events           : "<<fRejectedEvents<<std::endl;
    std::cout<<std::endl;
    std::cout<<"********************************************"<<std::endl;

    fBashColor->ResetColor();
}

SBREvent* StogsAFP::ProcessPacking(SBREvent *StogsEvt)
{
    TList *ListOfHits;
    ListOfHits = new TList();
    TList *ListOfPackedHits;
    ListOfPackedHits = new TList();

    int NbHits = StogsEvt->GetNbHits();

    for(int i=0 ; i<NbHits ; i++)
    {
        SBRHit *a_hit = StogsEvt->GetHit(i);

        //*10 to set distances in mm
        a_hit->fX = a_hit->fX*10.;
        a_hit->fY = a_hit->fY*10.;
        a_hit->fZ = a_hit->fZ*10.;

        //*1e3 to set energies in keV
        a_hit->fE = a_hit->fE*1e3;

        fBashColor->SetWarningOut();
        if(a_hit->fE > 1e5) cout<<"Warning !!! E = "<<a_hit->fE<<" keV, initial energies should be already in keV ?!"<<endl;
        if(TMath::Abs(a_hit->fZ) > 1000.)  cout<<"Warning !!! Z = "<<a_hit->fZ<<" mm, initial distances should be already in mm ?!"<<endl;
        fBashColor->ResetColor();

        ListOfHits->Add(a_hit);
    }

    while(ListOfHits->GetSize())
    {
        SBRHit *PackedHit = 0x0;

        for(int i=0 ; i<ListOfHits->GetSize() ; i++)
        {
            SBRHit *a_hit = (SBRHit*)ListOfHits->At(i);
            if(PackedHit == 0x0) PackedHit = a_hit;
            else if(a_hit->fE > PackedHit->fE) PackedHit = a_hit;
        }

        ListOfHits->Remove(PackedHit);
        ListOfPackedHits->Add(PackedHit);

        TVector3 PackedHitPosPond(PackedHit->fX*PackedHit->fE,PackedHit->fY*PackedHit->fE,PackedHit->fZ*PackedHit->fE);
        double PondFactor = PackedHit->fE;

        TVector3 PackedHitPos(PackedHit->fX,PackedHit->fY,PackedHit->fZ);

        for(int i=0 ; i<ListOfHits->GetSize() ; i++)
        {
            SBRHit *a_hit = (SBRHit*)ListOfHits->At(i);
            TVector3 HitPos(a_hit->fX,a_hit->fY,a_hit->fZ);

            TVector3 Diff = HitPos-PackedHitPos;

            double distance = Diff.Mag();

            if(distance<fSpatialResolution)
            {
                TVector3 HitPosPond(a_hit->fX*a_hit->fE,a_hit->fY*a_hit->fE,a_hit->fZ*a_hit->fE);
                PackedHitPosPond += HitPosPond;
                PackedHit->fE = PackedHit->fE + a_hit->fE;
                PondFactor += a_hit->fE;

                ListOfHits->Remove(a_hit);
                i--;
            }
        }

        PackedHitPosPond = PackedHitPosPond*(1./PondFactor);
        PackedHit->fX = PackedHitPosPond.x();
        PackedHit->fY = PackedHitPosPond.y();
        PackedHit->fZ = PackedHitPosPond.z();
    }

    SBREvent *PackedEvt;
    PackedEvt = new SBREvent();

    for(int i=0 ; i<ListOfPackedHits->GetSize() ; i++)
    {
        SBRHit *a_hit = (SBRHit*)ListOfPackedHits->At(i);
        PackedEvt->AddHit(a_hit);
    }

    delete ListOfHits;
    delete ListOfPackedHits;

    return PackedEvt;
}

void StogsAFP::ApplyEnergyResolution(SBREvent *StogsEvt)
{
    int NbHits = StogsEvt->GetNbHits();

    for(int i=0 ; i<NbHits ; i++)
    {
        SBRHit *a_hit = StogsEvt->GetHit(i);

        double Ener = a_hit->fE;

        double WFHM = (2.1-1.2)/(1333.-60.)*(Ener-1333.)+2.1;
        double sigma = WFHM/sqrt(8.*log(2.));

        Ener = TMath::Abs(gRandom->Gaus(Ener,sigma));

        a_hit->fE = Ener;
    }
}

void StogsAFP::ApplySmearing(SBREvent *StogsEvt)
{
    int NbHits = StogsEvt->GetNbHits();

    for(int i=0 ; i<NbHits ; i++)
    {
        SBRHit *a_hit = StogsEvt->GetHit(i);

        double x,y,z;
        gRandom->Sphere(x,y,z,fSpatialResolution);

        a_hit->fX += x;
        a_hit->fY += y;
        a_hit->fZ += z;
    }
}

SBREvent *StogsAFP::FilterEvent(SBREvent *GlobEvt)
{
    TList *ListOfHits;
    ListOfHits = new TList();

    TList *ListOfOrderedHits;
    ListOfOrderedHits = new TList();

    int NbHits = GlobEvt->GetNbHits();

    for(int i=0 ; i<NbHits ; i++)
    {
        SBRHit *a_hit = GlobEvt->GetHit(i);
        if(a_hit->fUID < 180) ListOfHits->Add(a_hit);
    }


    while(ListOfHits->GetSize())
    {
        SBRHit *FirstHit = (SBRHit*)ListOfHits->First();
        int firstUID = FirstHit->fUID;
        ListOfOrderedHits->Add(FirstHit);
        ListOfHits->Remove(FirstHit);

        for(int i=0 ; i<ListOfHits->GetSize() ; i++)
        {
            SBRHit *a_hit = (SBRHit*)ListOfHits->At(i);

            if(a_hit->fUID == firstUID)
            {
                ListOfOrderedHits->Add(a_hit);
                ListOfHits->Remove(a_hit);
                i--;
            }
        }
    }

    SBREvent *FilterdEvent;
    FilterdEvent = new SBREvent();

    for(int i=0 ; i<ListOfOrderedHits->GetSize() ; i++)
    {
        SBRHit *a_hit = (SBRHit*)ListOfOrderedHits->At(i);
        FilterdEvent->AddHit(a_hit);
    }

    delete ListOfHits;
    delete ListOfOrderedHits;

    return FilterdEvent;
}