GANIL/ADFWatchers.C

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2010 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.             *
 ***************************************************************************/

// Watcher definition
#ifndef _ADFWatchers
#include "ADFWatchers.h"
#endif

#include "TSystem.h"
#include <stdexcept>

using namespace ADF;


Bool_t ADFWatcher::SetTrigger( ADF::DFTrigger *trigger ) 
{
    if ( ! Watcher::SetTrigger(trigger) ){
        fFrame = 0x0;
        return false;
    }

    // set the unique frame to be watched
    fFrame = GetTrigger()->GetInputSharedFP();

    return true;
}


void Dumper::Exec(Option_t */*option*/) 
{
    // Frame by Frame action
    // be sure the key (the frame) has been set properly
    if ( fFrame == 0x0 )
        return;

    // if new file is asked, next open erase the file
    if ( fIsNew )
    { fFrame->GetFrame()->Dump(fFilename.Data(), fIsNew); fIsNew = false; }
    else
    { fFrame->GetFrame()->Dump(fFilename.Data(), fIsNew); }
}


void Dumper::SetOutputFile(const Char_t *filename)
{
    fFilename = filename; fIsNew = true;
}


AgataKeyWatcher::AgataKeyWatcher(const Char_t *name, const Char_t *title, TDirectory *sp_dir, TDirectory *tag_dir):
    ADFWatcher(name,title,sp_dir,tag_dir),
    fKey(0x0),
    fWhichBin(-1),
    fKeyFilter("*",true),
    fSubKeyFilter("*",true),
    fLastEventNumber(0u),
    flastTimestamp(0UL),
    fShowGlobalConf(false),
    fLengthFrame(0x0),
    fKeyMessage(0x0),
    fEventNumber(0x0),
    fTimeStamp(0x0),
    fCEventNumber(0x0),
    fCTimeStamp(0x0),
    fCoincTimeStamp(0x0),
    fLengthSubFrame(0x0),
    fKeyMessageSubFrame(0x0)
{
    // allocation of the histograms you would like to fill
    fLengthFrame = MakeTH1<TH1F>("LengthFrame","Length of the Frames (bytes)", aKByte,0,20*aKByte);
    fKeyMessage  = MakeTH1<TH1F>("KeyMessage","Message type", 3,0,3) ;
    #if ROOT_VERSION_CODE >= ROOT_VERSION(6,0,0)
    fKeyMessage->GetXaxis()->SetCanExtend(true);
    #else
    fKeyMessage->SetBit(TH1::kCanRebin);
    #endif
    fEventNumber = MakeTH1<TH1D>("EventNumber","Last event numbers", kMaxCycle,0.,kMaxCycle) ;
    fTimeStamp   = MakeTH1<TH1D>("TimeStamp","Last time stamps", kMaxCycle,0.,kMaxCycle) ;
    fCEventNumber = MakeTH1<TH1F>("CEventNumber","Consecutive event numbers", 1000,0, 10000) ;
    fCTimeStamp = MakeTH1<TH1F>("CTimeStamp","Consecutive time stamps", 1000,0, 10000000) ;
    fCoincTimeStamp = MakeTH1<TH1F>("CoincTimeStamp","Time stamp distribution in composite frames", 500,0,1000 ) ;
    fLengthSubFrame = MakeTH1<TH1F>("LengthSubFrame","Length of the SubFrames (bytes)", aKByte,0,20*aKByte);
    fKeyMessageSubFrame = MakeTH1<TH1F>("KeyMessageSubFrame","Message type for SubFrames in Composite Frame", 1,0,1);
    #if ROOT_VERSION_CODE >= ROOT_VERSION(6,0,0)
    fKeyMessageSubFrame->GetXaxis()->SetCanExtend(true);
    #else
    fKeyMessageSubFrame->SetBit(TH1::kCanRebin);
    #endif
}


AgataKeyWatcher::~AgataKeyWatcher() 
{
    // destruction of the histograms done because they have been added to the pool
}


void AgataKeyWatcher::SetKeyFilter(const Char_t *keypattern, const Char_t *subkeypattern)
{
    TString stmp1(keypattern), stmp2(subkeypattern);
    
    Bool_t wildcard1 = false, wildcard2 = false;
    if ( stmp1.Contains("*") )
        wildcard1 = true;
    if ( stmp2.Contains("*") )
        wildcard2 = true;
    
    TRegexp tmp1(keypattern,wildcard1), tmp2(subkeypattern,wildcard2);

    // change the filters only if it is a valid expression
    if ( tmp1.Status() == TRegexp::kOK ) {
        fKeyFilter = tmp1;
    }
    if ( tmp2.Status() == TRegexp::kOK ) {
        fSubKeyFilter = tmp2;
    }
}


void AgataKeyWatcher::DoCanvas(TCanvas *c, Option_t *o)
{
    TString option = o;
    if ( option == "" ) {
        c->Divide(3,2,0.001,0.001);
        c->cd(1);
        fLengthFrame->Draw();
        c->cd(2);
        fLengthSubFrame->Draw();
        c->cd(3);
        fEventNumber->Draw();
        c->cd(4);
        fKeyMessage->Draw();
        c->cd(5);
        fKeyMessageSubFrame->Draw();
        c->cd(6);
        fTimeStamp->Draw();
    }
    if ( option == "q" ) {
        c->Divide(1,3,0.001,0.001);
        c->cd(1);
        fCEventNumber->Draw();
        c->cd(2);
        fCTimeStamp->Draw();
        c->cd(3);
        fCoincTimeStamp->Draw();
    }
}


void AgataKeyWatcher::Exec(Option_t */*option*/) 
{
    if ( gDebug > 2 )
        printf("AgataKeyWatcher::Exec is called \n");

    UInt_t length, current, eventnumber, message; ULong64_t timestamp; TString mess, mess_sec, mess_tmp;

    if ( fFrame != 0x0 ) {
        if ( fFrame->IsValid() )
            fKey = dynamic_cast<const AgataKey *>(fFrame->GetFrame()->GetKey());
        else
            fKey = 0x0;
    }
    else
        fKey = 0x0;

    // be sure the key (the frame) has been set properly
    if ( fKey == 0x0 )
        return;

    // first test if a filter is applied
    mess = AgataKeyFactory::theFactory()->GetMessage(fKey->GetMessage());
    if ( ! mess.Contains(fKeyFilter) )
        return;

    // now fill histograms
    fLengthFrame->Fill(fKey->GetFrameLength()); fKeyMessage->Fill(mess.Data(),1);
    if ( mess.Contains("UnknownKey") ) { // an unknown key
        std::cout << " WARNING, unknown key " << std::endl;
        fKey->Print();
    }
    // try and read the possible event# and timestamp (may not exists in case of low version number of agata key
    const BufferIO *buf = fFrame->GetFrame()->GetRealBuffer();
    buf->SetOffset();
    if ( fKey->GetFrameLength() < 12u ) {
        eventnumber = 0u;
        timestamp   = 0UL;
    }
    else {
        (*buf) >> eventnumber;
        if ( fKey->GetFrameLength() < 20u )
            timestamp = 0UL;
        else (*buf) >> timestamp;
    }

    // which is negative if this method has never been called
    if ( fWhichBin == -1 ) {
        fLastEventNumber = eventnumber;
        flastTimestamp   = timestamp;
        fWhichBin = 0;
    }

    // fill timestamp and event number related
    fEventNumber->Fill(fWhichBin,eventnumber);
    fTimeStamp->Fill(fWhichBin,timestamp);

    fCEventNumber->Fill(eventnumber-fLastEventNumber);
    fCTimeStamp->Fill(timestamp-flastTimestamp);

    /*
    if ( eventnumber < fLastEventNumber ) {
        fEventNumberQuality->Fill(-1);
    //  cout << " Current event # " << eventnumber << " Previous event # " << fLastEventNumber << endl;
    }
    else
        eventnumber == fLastEventNumber ? fEventNumberQuality->Fill(0) : fEventNumberQuality->Fill(1);

    if ( timestamp < flastTimestamp ) {
        fTimeStampQuality->Fill(-1);
    //  cout << " Current timestamp # " << timestamp << " Previous timestamp # " << flastTimestamp << endl;
    }
    else
        timestamp == flastTimestamp ? fTimeStampQuality->Fill(0) : fTimeStampQuality->Fill(1); */

    fLastEventNumber = eventnumber;
    flastTimestamp       = timestamp;

    // check if message contains conf. if yes dump the content in a file
    if ( mess.Contains("conf") && fShowGlobalConf ) {
        FILE *fp; Short_t which = 0;
        while ( which < 20 ) {
            if ( gSystem->AccessPathName(Form("./%s%d.conf",mess.Data(),which)) )
                break;
            which++;
        }
        fp = ::fopen(Form("./%s%d.conf",mess.Data(),which),"w"); fFrame->GetFrame()->Dump(fp); ::fclose(fp);
    }
    if ( mess.Contains("event") ) {

        fKeyMessageSubFrame->Fill(mess.Data(),1);

        // try to look for a true message type at the beginning of the buffer
        Bool_t is_a_start = false; Short_t max_loop = 0u;
        buf->SetOffset();
        while ( buf->FreeSize() > 20u ) {

            // up to ten embedded frames in one frame
            if ( max_loop++ == 10u )
            { fLog << " Max has been reached in scan of composite frame " << nline; break; }

            current = buf->Offset();

            (*buf) >> length >> message >> eventnumber >> timestamp ;

            mess_tmp = AgataKeyFactory::theFactory()->GetMessage(message);

            //                  cout << "[" << mess_tmp.Data() << "]" << endl;
            if ( mess_tmp.Contains("UnknownKey") ) { // an unknown key
                is_a_start = false;
            }
            else is_a_start = true;

            if ( is_a_start )  { // a known key has been found

                buf->SetOffset( current + length );

                mess_sec  = "|";
                mess_sec += mess_tmp;
                mess_sec += "|";

                mess_sec = Form("|%s|",mess_tmp.Data());

                if ( mess_tmp.Contains(fSubKeyFilter ) ) {

                    fLengthSubFrame->Fill(length);
                    fKeyMessageSubFrame->Fill(mess_sec.Data(),1);

                    fCoincTimeStamp->Fill(timestamp-flastTimestamp);
                }
            }
            else buf->SetOffset( current + 4u );
        }
    }

    if ( fWhichBin == kMaxCycle ) {
        fWhichBin = 0u;
        fEventNumber->Reset("ICE"); fTimeStamp->Reset("ICE");
    }
    else fWhichBin++;
}


TSRate::TSRate(const Char_t *name, const Char_t *title, TDirectory *sp_dir, TDirectory *tag_dir) :
    ADFWatcher(name,title,sp_dir,tag_dir),
    fTSRates(0)
{ 
}


void TSRate::Exec(Option_t * /*option*/)
{       
    // fLog << debug << " In TSRate::Exec() " << nline;

    if ( fFrame && fFrame->IsValid() ) {
        AgataKey *key = dynamic_cast<AgataKey *>(fFrame->GetFrame()->GetKey());
        if ( key ) {
            ULong64_t ts = key->GetTimeStamp() ;
            if ( fTSRates.size() > 0 ) {
                if ( fTSRates.back() < ts )
                    fTSRates.push_back( ts ) ;
                //                              else
                //                                      std::cout << " Two consecutive TS in a bunch inversed !" << std::endl;
            }
            else fTSRates.push_back( ts ) ;
        }
    }

    // fLog << debug << " Out TSRate::Exec() " << dolog;

}


#include "TMath.h"

TSRateCollector::TSRateCollector(const Char_t *name, const Char_t *title, TDirectory *sp_dir, TDirectory *tag_dir) :
    Watcher(name,title,sp_dir,tag_dir),
    fDefinitionRate(),
    fCurrentRate(),
    fTSRates(),
    fRates(),
    fPacketWidth(),
    fDeadTime(),
    fT0(kMaxULong64)
{
    fTSRates.SetOwner(false); fPacketWidth.SetOwner(false); fDeadTime.SetOwner(false);

    // rate every 30 second for one hour
    fDefinitionRate.Add(30*fgkOneSecond,120);
    // rate every 10 minutes for 12 hour
    fDefinitionRate.Add(600*fgkOneSecond,72);
    // rate every hour for 48 hours
    fDefinitionRate.Add(3600*fgkOneSecond,48);
}


Bool_t TSRateCollector::SetT0()
{
    ULong64_t val, min = kMaxULong64;

    // try and see if first time it is called or if timestamp has been reset
    TIter next(&fTSRates);
    TSRate *ts_rate;
    while ( (ts_rate = (TSRate *)next()) ) {
        if ( ts_rate->LastValues().size() ) {
            val = ts_rate->LastValues().at(0);
            if ( val < min )
                min = val;
        }
    }
    if ( min < fT0  ) {
        // change histograms title to include starting time
        TDatime dt;
        for (Int_t i = 0; i < fRates.GetHists()->GetSize(); i++) {

            // reset stat
            fCurrentRate[i].Reset(fDefinitionRate.fNbScales);

            if ( fRates.GetHists()->At(i) == 0x0 )
                continue;
            TString tmp = fRates.GetHists()->At(i)->GetName();
            tmp += " since ";
            tmp += dt.AsString();

            ((TH1 *)fRates.GetHists()->At(i))->Reset();
            ((TH1 *)fRates.GetHists()->At(i))->SetTitle(tmp.Data());

        }
        fT0 = min;
        std::cout << " New TIMESTAMP in TSRateCollector" << fT0 << std::endl;
    }

    return (fT0 != kMaxULong64);
}


void TSRateCollector::AddHook(Watcher *ts_rate, const char *fd_name)
{
    const Short_t wheel_size = 5;
    Color_t color_wheel[wheel_size] = { kBlue, kGreen+2, kOrange, kRed, kViolet };
    Style_t style_wheel[wheel_size] = { 3001, 3002, 3003, 3004, 3005 };

    if ( !ts_rate->InheritsFrom("TSRate") )
        return;

    if ( ts_rate ) {
        // keep an hook on it and creates the associated struct to compute stats
        fTSRates.Add ( ts_rate );
        C_Rate c_rate(fDefinitionRate.fNbScales); fCurrentRate.push_back( c_rate );

        // now add histograms
        TString hname = fd_name;
        hname += "_Rate";

        Int_t hsize = 0;
        for (Int_t i = 0; i < fDefinitionRate.fNbScales; i++)
            hsize += fDefinitionRate.fNbSamples[i] + gEmptyBins;

        hsize += 10;
        TH1F *h = MakeTH1<TH1F>(hname.Data(),"Rate",hsize,0,hsize);
        fRates.Add(h);
        h->SetFillColor(color_wheel[(fRates.GetHists()->GetSize()-1)%wheel_size]);
        h->SetFillStyle(style_wheel[(fRates.GetHists()->GetSize()-1)%wheel_size]);
        h->SetXTitle("[0][30s/bin = 1hour] -- [170][10m/bin = 12hours] -- [292][1h/bin = 2days]");
        h->SetYTitle("k/s");

        hname = fd_name; hname += "_Width";
        h = MakeTH1<TH1F>(hname.Data(),"Width distribution of packets (s)",200,0,20);
        fPacketWidth.Add(h);
        h->SetFillColor(color_wheel[(fRates.GetHists()->GetSize()-1)%wheel_size]);
        h->SetFillStyle(style_wheel[(fRates.GetHists()->GetSize()-1)%wheel_size]);

        hname = fd_name; hname += "_Deadtime";
        h = MakeTH1<TH1F>(hname.Data(),"Dead time distribution between consecutive packets (s)",100,0,100);
        fDeadTime.Add(h);
        h->SetFillColor(color_wheel[(fRates.GetHists()->GetSize()-1)%wheel_size]);
        h->SetFillStyle(style_wheel[(fRates.GetHists()->GetSize()-1)%wheel_size]);
    }
}


void TSRateCollector::DoCanvas(TCanvas *c, Option_t * /*o*/)
{
    Int_t nbcolumn = 3; c->Divide(nbcolumn,fRates.GetHists()->GetSize()/nbcolumn,0.001,0.001);

    for (Int_t i = 0; i < fRates.GetHists()->GetSize(); i++) {
        c->cd(i+1);
        TH1 *h = (TH1F *)fRates.GetHists()->At(i);
        if ( h )
            h->Draw();
    }
    TString cmd = Form("RunLoupeOnPad( (TCanvas *)gROOT->GetListOfCanvases()->FindObject(\"%s\") )", c->GetName());
    gROOT->ProcessLine(cmd.Data());
}


void TSRateCollector::Exec(Option_t * /*option*/)
{       
    const static Double_t to_rate = fgkOneSecond/(1000.);

    TSRate *ts_rate; TH1 *h; ULong64_t dt0, dt, slot, cycle;

    if ( gDebug > 0 )
        printf("TSRateCollector::Exec is called \n");

    if ( fTSRates.GetSize() == 0 )
        return; // to avoid doing things in case there are no input

    // first time this is called, computes T0
    if ( ! SetT0() ) {
        printf(" WARNING in TSRateCollector::Exec Cannot get First timestamp \n");
        return;
    }

    // loop on all registered TSRate

    for (Int_t i = 0; i < fTSRates.GetSize(); i++) {

        // get TSRate
        ts_rate = (TSRate *)fTSRates.At(i);

        //              std::cout << " TSRate name : " << ts_rate->GetName()
        //                              << " , registered @ position " << i << " , has view # TS : " << ts_rate->LastValues().size() << std::endl;

        // nothing to do
        if ( ts_rate->LastValues().size() < 2 )
            continue;

        // compute some value to have an idea of the TS structure
        h = (TH1F *)fPacketWidth.At(i);
        Double_t width = ts_rate->LastValues().at(ts_rate->LastValues().size()-1) - ts_rate->LastValues().at(0);
        width /= fgkOneSecond;
        h->Fill(width);

        h = (TH1F *)fDeadTime.At(i);
        if ( ts_rate->LastValues().at(0) > fCurrentRate[i].fLastTS ) {
            h->Fill((ts_rate->LastValues().at(0)-fCurrentRate[i].fLastTS)/fgkOneSecond);
        }
        else { // ste underflow
            h->Fill(-(fCurrentRate[i].fLastTS-ts_rate->LastValues().at(0))/fgkOneSecond);
            //                      continue;
        }
        // h->Fill(ts_rate->LastValues().size());

        //              std::cout << " WIDTH " << width << std::endl;

        // keep last TS of the packet for next packet
        fCurrentRate[i].fLastTS = ts_rate->LastValues().at(ts_rate->LastValues().size()-1);

        /*              // only one frame ... very low counting rate ...
        if ( ts_rate->LastValues().size() < 2 ) {
//                      ts_rate->LastValues().resize(0);
            // do not resize so that the value is kept in the TSRate and a new one is going to be added
            // not really a good way to proceed but should give a quite good idea in case of very low counting rates
            continue;
        } */

        // statistics
        h = (TH1F *)fRates.GetHists()->At(i);

        for (size_t j = 1; j < ts_rate->LastValues().size(); j++) {

            if ( ts_rate->LastValues().at(j) < fT0 ) {
                std::cout << " ERROR Timestamp in array smaller that T0 " << std::endl;
                continue;
            }

            // compute time since T0 and since first timestamp
            dt0 = ts_rate->LastValues().at(j) - fT0;
            // compute time between two consecutive TS and the associated rate
            dt  = ts_rate->LastValues().at(j) - ts_rate->LastValues().at(j-1);

            // loop on the different scales to compute statistics
            Int_t h0 = 0;
            for (Int_t k = 0; k < fDefinitionRate.fNbScales; k++) {
                // cycle since T0
                cycle = dt0 / fDefinitionRate.fPrecision[k];
                //
                if ( cycle == fCurrentRate[i].fCycle[k] ) { // add to current for integration
                    fCurrentRate[i].fSum[k] += dt;
                    fCurrentRate[i].fNbCounts[k]++;
                }
                else { // a new input in the histogram

                    // bin number to be changed
                    slot = fCurrentRate[i].fCycle[k] % fDefinitionRate.fNbSamples[k];

                    // fill h
                    // it has to be with Sum to take into account the real time counting time and remove empty cycle
                    /*                                  if ( slot == 0 ) {
                        for (Int_t bin = h0; bin < h0+fDefinitionRate.fNbSamples[k]; bin++) {
                            h->SetBinContent(bin+1,0);
                        }
                    } */
                    if ( fCurrentRate[i].fSum[k] > 0
                         && ((slot+h0+1) < UInt_t(h->GetNbinsX()))
                         && ((slot+h0+1) > 0) ) {
                        h->SetBinContent(slot+h0+1,fCurrentRate[i].fNbCounts[k]*to_rate/fCurrentRate[i].fSum[k]);
                    }
                    else {
                        std::cout << " SUM is O !! "
                                  << fCurrentRate[i].fSum[k] << " j " << j << " " << ts_rate->LastValues().size()
                                  << " " << fCurrentRate[i].fNbCounts[k] << " " << dt << std::endl;
                    }

                    /*                                  else
                     std::cout << slot+h0+1 << std::endl;
                     std::cout << fCurrentRate[i].fNbCounts[k]*to_rate/fCurrentRate[i].fSum[k] << std::endl; */

                    // keep the current situation
                    fCurrentRate[i].fSum[k] = dt;
                    fCurrentRate[i].fNbCounts[k]  = 1;

                    fCurrentRate[i].fCycle[k] = cycle;
                }
                h0 += (fDefinitionRate.fNbSamples[k] + gEmptyBins);
            }
        }

        ts_rate->LastValues().resize(0);
    } // i

    if ( 0 )
        printf("TSRateCollector::Exec is done \n");
}


TimeStampWatcher::TimeStampWatcher( const char * name, const char * title, TDirectory *sp_dir, TDirectory *tag_dir):
    ADFWatcher(name,title,sp_dir,tag_dir)
{
    int TMax = 3600;
    TimeStampHist = MakeTH1<TH1F>("TimeStampHist", "Time stamp distribution",TMax*10, 0, TMax);
    FirstTimeStamp = 0;
}


TimeStampWatcher::~TimeStampWatcher()
{
}


void TimeStampWatcher::DoCanvas(TCanvas * /*c*/, Option_t *o)
{
    TimeStampHist->Draw(o);
}


void TimeStampWatcher::Exec(Option_t * /*option*/)
{

    ULong64_t fTimestamp = ((AgataKey *)fFrame->GetFrame()->GetKey())->GetTimeStamp();

    if(fTimestamp != 0 && FirstTimeStamp==0) FirstTimeStamp = fTimestamp;

    double TSinSec = (fTimestamp-FirstTimeStamp)*10*1e-9;

    TimeStampHist->Fill(TSinSec);
    LastTimeStamp = fTimestamp;

}


void TimeStampWatcher::Zero(Option_t * /*opt1*/, Option_t * /*opt2*/)
{
    TimeStampHist->Reset();
    FirstTimeStamp = LastTimeStamp;
}


ClassImp(ADFWatcher)
ClassImp(Dumper)
ClassImp(AgataKeyWatcher)
ClassImp(TSRate)
ClassImp(TSRateCollector)
ClassImp(TimeStampWatcher)