RADConverter.cpp

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2004-2006 by Olivier Stezowski & Christophe Theisen     *
 *   stezow(AT)ipnl.in2p3.fr, christophe.theisen(AT)cea.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.             *
 ***************************************************************************/
 
#ifndef GW_RADCONVERTER_H
#include "RADConverter.h"
#endif


#ifndef ROOT_TRegexp
#include <TRegexp.h>
#endif

#ifndef ROOT_TSystem
#include <TSystem.h>
#endif

#ifndef ROOT_TSystemDirectory
#include <TSystemDirectory.h>
#endif

#if ROOT_VERSION_CODE >= ROOT_VERSION(5,15,0)
#include "TList.h"
#endif

#include <iostream>
using namespace std;

namespace { RADConverter gRADConverter; } // prototype ONLY use to create a new DB system

//___________________________________________________________________
RADConverter::RADConverter(const char *name) : HistoConverter(name)
{ 
    fType    = "rad";
    CheckDirectory(name);
    fCycle = '_'; // because ; is not supported by all file systems 
    fLog.GetProcessName() = "RADConverter";
} 

//___________________________________________________________________
RADConverter::RADConverter() : HistoConverter()
{ 
    fType = "rad";
    fFileExt = ".spe";
    fOptions = "";
    // add this interface to the DB service 
    AddPrototype(GetType(),this);
} 

//___________________________________________________________________
RADConverter::~RADConverter()
{ 

} 

//___________________________________________________________________
void RADConverter::ls(Option_t *o) const 
{
        HistoConverter::ls(o);
        TSystemDirectory dir; dir.SetDirectory(fName.Data()); TList *list = dir.GetListOfFiles();
                
        TIter iter(list); TObject *entry; int count = 0;
        while ( (entry = iter()) ) {
                TString tmp = entry->GetName();
                //cout << " file name " << entry->GetName() << endl;            
                if ( tmp.EndsWith(fFileExt) ) {
                        cout << " [" << count++ <<  "] " << entry->GetName() ;
                        
                        if ( (count+1) % 5 == 0 ) cout << endl;
                        else cout << "\t";      
                }
        }
        cout << endl; if ( list ) delete list;
}

//___________________________________________________________________
HistoConverter *RADConverter::NewDB(const char *name) const 
{ 
  return new RADConverter(name);
}

//___________________________________________________________________
TH1 *RADConverter::Get(int which)
{ 
        // list all the file to find the entry name correponding to which
        TSystemDirectory dir; dir.SetDirectory(fName.Data()); TList *list = dir.GetListOfFiles();
                
        TIter iter(list); TObject *entry; int count = 0; TString result = "_EMPTY_";
        while ( (entry = iter()) ) {
                TString tmp = entry->GetName();
                if ( tmp.EndsWith(fFileExt) ) {
                        if ( which == count ) result = entry->GetName(); count++;
                }
        }
        if ( list ) delete list;

        if ( result != "_EMPTY_" ) return Get(result.Data()); // an entry has been found
        fLog.GetProcessMethod() = "Get(Int_t )";
        SetError(HistoConverter::kFail,
                Form("Cannot find entry # %d",which));
                        
        return NULL;
} 

//___________________________________________________________________
TH1 *RADConverter::Get(const char *name)
{ 
        TH1F *h = new TH1F(); 
        
        h->SetName(name);
        if ( !Read((*h)) ) { delete h; h = NULL; } 
        
        return h;       
} 

//___________________________________________________________________
bool RADConverter::Read(TH1 &h)
{
        TString filename = fName.Data();  // set file name as a full path name
        if ( filename.IsNull() )  filename = "."; // for mine types

        if ( ! filename.EndsWith("/") ) filename+= "/";  // set file name as a full path name
                
        filename += h.GetName();
        if ( filename.EndsWith(fFileExt) == false ) filename += fFileExt;
        
        if (fVerboseLevel)
                fLog << info << "Reading spectrum " << filename << nline;
        
        ifstream file;
        
        file.open(filename.Data());
        if ( ! file.is_open() ){
                fLog.GetProcessMethod() = "Read(TH1& )";
                SetError(HistoConverter::kFail,
                        Form("Cannot open file %s",filename.Data()));   
                return false;
        }
        
        // now read the spectrum
        char nom[24], nam[8]; int ndim[4], dum[2], i;
        float buffer[32768];
        
        file.read(nom,4);
        file.read(nam,8);
        if ( gDebug > 0 ) fLog << debug << "Spectrum name : " << nam << nline;
//      sprintf(name,nam);
        file.read((char*)ndim, 16);
        if (ndim[0] >= 65536) {// swapping is needeed
          fIsSwapping = true;
          for(i=0; i<4; i++) ndim[i] = Int(((char*)&(ndim[i])));
        }

        if ( gDebug > 0 ) fLog << debug << "Spectrum size : " << ndim[0] << nline;
//      printf("%d \n",ndim[1]);
//      printf("%d \n",ndim[2]);
//      printf("%d \n",ndim[3]);

        file.read((char*)dum, 8);
        file.read((char*)buffer, 4*ndim[0]);
//      printf("%d \n",n);
//      printf("Spectrum size : %d \n",ndim[0]);

         
        if ( file.good() ) { // change h only if everything was read correctly
                h.SetName(nam); 
                h.SetBins(ndim[0],0,ndim[0]); 
                for (i=0; i<ndim[0]; i++) 
                        h.Fill(i,Float((char*)&(buffer[i])));
        }

        return file.good();
} 

//___________________________________________________________________
bool RADConverter::Write(const TH1 &h)
{ 
        TString hname, filename;  // set file name as a full path name
        hname = h.GetName();
        filename = CheckFileVersion(hname);
        
       fLog.GetProcessMethod() = "Write(const TH1& )";
  
        filename.Prepend(fName.Data());
        if ( gDebug > 0 ) {
          fLog << error << "File to store the spectrum " 
               << filename << nline ;    
        }
        
        if (fVerboseLevel)
                fLog << info << "Writing spectrum " << filename << nline ;


        FILE *fprad;
    int /*dum,*/ i, j = 24, toto, tab[4] = {0,1,1,1}, size; char spname[8];
        
        ::strncpy(spname,h.GetName(),8);
        tab[0] = size = h.GetNbinsX(); toto = size * sizeof(float); 
        
        fprad = fopen(filename.Data(),"w");
        if (fprad == NULL){
          SetError(HistoConverter::kFail,
                   Form("[RADConverter::Write] Cannot open %s in writing mode",filename.Data()));
          return false;
        }
        float *spectre = new float[h.GetNbinsX()];
        for ( i = 0; i < tab[0]; i++ ) 
                spectre[i] = h.GetBinContent(i+1);
                
//      dum = fwrite(&j, sizeof(int), 1, fprad);
//      dum = fwrite(spname, 8, 1, fprad);
//      dum = fwrite(tab, sizeof(int), 4, fprad);
//      dum = fwrite(&j, sizeof(int), 1, fprad);
//      dum = fwrite(&toto, sizeof(int), 1, fprad);
//      dum = fwrite(spectre, sizeof(float), size, fprad);
//  dum = fwrite(&toto,sizeof(int),1,fprad);

    fwrite(&j, sizeof(int), 1, fprad);
    fwrite(spname, 8, 1, fprad);
    fwrite(tab, sizeof(int), 4, fprad);
    fwrite(&j, sizeof(int), 1, fprad);
    fwrite(&toto, sizeof(int), 1, fprad);
    fwrite(spectre, sizeof(float), size, fprad);
    fwrite(&toto,sizeof(int),1,fprad);

        fclose(fprad); delete spectre; return true;
} 

//___________________________________________________________________
bool RADConverter::Read(TH2 &h)
{       
        // set file name as a full path name
        TString filename = fName.Data(); filename += h.GetName();
        
  fLog.GetProcessMethod() = "Read(TH2& )";

        // CHECK OUT THE EXTENSION
        int whichext = 0;       
        if ( filename.EndsWith(".mat") ) whichext = 1;  
        if ( filename.EndsWith(".m4b") ) whichext = 2;  
        if ( filename.EndsWith(".spn") ) whichext = 3;  
        
        if ( whichext == 0 ) filename += ".mat"; //assumes mat as a standard
        
        if (fVerboseLevel)
                fLog << info << "Reading matrix " << filename << nline;
        
        // now open the file for reading
        ifstream file;
        
        file.open(filename.Data());
        if ( ! file.is_open() ){
          SetError(HistoConverter::kFail,
                   Form("Cannot open file %s",filename.Data()));
          return false;
        }
        // check out the file size (should be divided by 4096*4096) and set the item size
        int itemsize;
        file.seekg(0,ios::end); long filesize = file.tellg() ;
        //cout << filesize << endl;
        if ( (filesize != 2*4096*4096) && (filesize != 4*4096*4096) ) { 
          fLog.GetProcessMethod() = "Read(TH2& )";
          SetError(HistoConverter::kFail,
                   Form("File %s cannot be a Radware matrix",filename.Data()));
          file.close(); return false;           
        }
        else { itemsize = file.tellg() / (4096*4096) ; }
        file.seekg(0,ios::beg);
        bool ok = false;
        if ( itemsize == 2 ) { ok = DecodeMAT(h,file); }
        if ( itemsize == 4 ) { ok = DecodeSPN(h,file); }
        
        file.close(); return ok;
}

//___________________________________________________________________
bool RADConverter::DecodeMAT(TH2 &h, ifstream & file)
{
        // Read the matrix for the file. First scan to determine if the file can be read 
        // without any problem and if swapping is needed.
        // swapping will be determined by computing the value: Sum ( ABS(mat(i,j) )
        // The swapping that minimize this value will be considered as the good swapping
        
        int i,j; char buffer[2*4096]; 
        
        double swap = 0.0, noswap = 0.0;
        for(i=0; i<4096 ; i++){
                //printf("Reading row #%d \r",i);
                file.read(buffer,2*4096);
                if ( file.fail() ) {
                  fLog.GetProcessMethod() = "DecodeMAT(TH2&, ifstream& )";
                  SetError(HistoConverter::kFail,
                           "Cannot read correctly the file");
                  return false; 
                }
                else {
                        for (j=0; j<4096; j++) {
                                fIsSwapping = true;  swap   += TMath::Abs(Short(buffer+2*j));
                                fIsSwapping = false; noswap += TMath::Abs(Short(buffer+2*j));   
                        }
                }
        }       
        if ( swap > noswap ) fIsSwapping = false;
        else fIsSwapping = true;
        
        // rewing and re-start and fills the histograms
        file.seekg(0,ios::beg); 
        
        h.SetBins(4096,0,4096,4096,0,4096);
        for(i=0; i<4096 ; i++){
                printf("Reading row #%d \r",i);
                file.read(buffer,2*4096);
                for (j=0; j<4096; j++) h.Fill(i,j,Short(buffer+2*j));
        }

        return true;
}

//___________________________________________________________________
bool RADConverter::DecodeSPN(TH2 &h, ifstream & file)
{
        // Read the matrix for the file. First scan to determine if the file can be read 
        // without any problem and if swapping is needed.
        // swapping will be determined by computing the value: Sum ( ABS(mat(i,j) )
        // The swapping that minimize this value will be considered as the good swapping
        
        int i,j; char buffer[4*4096]; 
        
        double swap = 0.0, noswap = 0.0;
        for(i=0; i<4096 ; i++){
                // printf("Reading row #%d \r",i);
                file.read(buffer,4*4096);
                if ( file.fail() ) {
                  fLog.GetProcessMethod() = "DecodeSPN(TH2&, ifstream& )";
                  SetError(HistoConverter::kFail,
                           "Cannot read correctly the file");
                  return false; 
                }
                else {
                        for (j=0; j<4096; j++) {
                                fIsSwapping = true;  swap   += TMath::Abs(Int(buffer+4*j));
                                fIsSwapping = false; noswap += TMath::Abs(Int(buffer+4*j));     
                        }
                }
        }       
        if ( swap > noswap ) fIsSwapping = false;
        else fIsSwapping = true;
        
        // rewing and re-start and fills the histograms
        file.seekg(0,ios::beg);         
        
        h.SetBins(4096,0,4096,4096,0,4096);
        for(i=0; i<4096 ; i++){
                printf("Reading row #%d \r",i);
                file.read(buffer,4*4096);
                for (j=0; j<4096; j++){
                 h.Fill(i,j,Int(buffer+4*j));
                } 
        }

        return true;
}

//___________________________________________________________________
bool RADConverter::Write(const TH2 &h)
{ 
        // only TH2S and TH2I can be exported
        fLog.GetProcessMethod() = "Write(TH2& )";
  
        if ( !h.InheritsFrom("TH2S") && !h.InheritsFrom("TH2I") ) {
          fLog.GetProcessMethod() = "Write(const TH2& )";
          SetError(HistoConverter::kFail,
                   "Only TH2S or TH2I can be exported as Radware matrix");  
          return false; 
        }
        
        TString hname, filename, ext;  // set file name as a full path name
        
        // check out the histogram name to avoid blank in a filename
        hname = h.GetName(); hname.ReplaceAll(" ",""); 
        
        // CHECK OUT THE EXTENSION
        if ( hname.EndsWith(".mat") ) { hname.Remove(hname.Last('.'));} 
        if ( hname.EndsWith(".m4b") ) { hname.Remove(hname.Last('.'));} 
        if ( hname.EndsWith(".spn") ) { hname.Remove(hname.Last('.'));} 
        
        filename  = hname;  
        if ( h.InheritsFrom("TH2S") ) 
                { filename += ".mat"; ext = "mat"; }
        else 
                { filename += ".m4b"; ext = "m4b"; }
                

        
        // determine the name of the file taking into account the versionning system
        TSystemDirectory dir; 
        dir.SetDirectory(fName.Data()); TList *list = dir.GetListOfFiles();     

        if ( list ) {
                // look for the next cycle 
                int cycle;
                for ( int i = 0; i <= MAXCYCLE; i++ ) { // up to MAXCYCLE files
                
                        TString tmp;
                        if ( i == 0 ) tmp.Form("%s.%s",hname.Data(),ext.Data()); 
                        else tmp.Form("%s%c%d.%s",hname.Data(),fCycle,i,ext.Data());    
                        
                        cycle = i;
                        if ( !list->FindObject(tmp.Data()) ) { break; }
                } 
                if ( cycle == 0) { // the file doesn't exist yet 
                        filename.Form("%s.%s",hname.Data(),ext.Data());
                }
                if ( cycle == MAXCYCLE ) { // force overwriting even if in versionning mode
                        filename.Form("%s%c%d.%s",hname.Data(),fCycle,MAXCYCLE,ext.Data());     
                }
                if ( 0 < cycle && cycle < MAXCYCLE ) { 
                        if ( fMode == HistoConverter::kVersion ) // ready for a new file
                          filename.Form("%s%c%d.%s",hname.Data(),fCycle,cycle,ext.Data());
                        else { // comes back to the previous one
                          cycle--;
                          if ( cycle == 0 ) 
                                filename.Form("%s.%s",hname.Data(),ext.Data());
                          else 
                                filename.Form("%s%c%d.%s",hname.Data(),fCycle,cycle,ext.Data());          
                        }
                }       
                if ( gDebug > 0 ) 
                   fLog << debug << "cycle is " << cycle << nline;    
                
                delete list;
        } 
        filename.Prepend(fName.Data());
        if ( gDebug > 0 ) {
                fLog << error << "File to store the matrix " 
                        << filename << nline;    
        }
        
        if (fVerboseLevel)
                fLog << info << "Writing matrix " << filename << nline; 
                
        // now the matrix could be store
        ofstream file;
        
        file.open(filename.Data());
        if ( ! file.is_open() ) {
                  SetError(HistoConverter::kFail,
                           Form("Cannot open file %s",filename.Data()));
                return false;
        }
        
        bool ok;
        for(int i=0; i<4096 ; i++) {
                for (int j=0; j<4096; j++) {    
                        if ( h.InheritsFrom("TH2S") ) {
                          short value = 0;
                          if ( i < h.GetNbinsX() &&  j < h.GetNbinsY() ) 
                                value = (short)h.GetBinContent(i+1,j+1);        
                          file.write((char *)&value,2);
                        }
                        else {
                          int value;
                          if ( i < h.GetNbinsX() &&  j < h.GetNbinsY() ) 
                                value = (int)h.GetBinContent(i+1,j+1);  
                          file.write((char *)&value,4);
                        }
                }
        }
        ok = !file.fail(); file.close(); return ok;
} 

//___________________________________________________________________
bool RADConverter::Read(THStack &hs)
{       
        TString hsname = hs.GetName();
        
        // to avoid blanck in name and to get the correct extention .spe in pattern
        hsname.ReplaceAll(" ",""); 
        if ( !hsname.EndsWith(fFileExt) )
        hsname.Append(fFileExt);
    
        // pattern expression
    TString reg = "*" + fFileExt; Bool_t wildcard = false;
    if ( hsname.Contains("*") )
        wildcard = true;
        TRegexp all(reg,kTRUE), pattern(hsname.Data(),wildcard);
        
        if ( pattern.Status() != TRegexp::kOK ) { 
                if ( gDebug > 0 ) cerr << "pattern [*] intead of " << hsname << endl ;
                pattern = all; 
        }
        else { if ( gDebug > 0 ) cerr << "pattern " << hsname << endl ; }
        
        // loop on files to load the good ones.
        TSystemDirectory dir; dir.SetDirectory(fName.Data()); TList *list = dir.GetListOfFiles();
                
        TIter iter(list); TObject *entry; 
        while ( (entry = iter()) ) {
        
                TString tmp = entry->GetName();
                
                if ( tmp.Contains(pattern) ) { // should be added to the collection
                  if (fVerboseLevel) {
                    fLog.GetProcessMethod() = "Read(THStack& )";
                    fLog << info << "Reading spectrum " << entry->GetName() << nline;
                  }
                        TH1 *h = Get(entry->GetName()); if ( h ) { hs.Add(h); }
                }
        }
        if ( list ) delete list;
        
        return true;
} 

//___________________________________________________________________
bool RADConverter::Write(const THStack &/*hs*/)
{ 
  fLog.GetProcessMethod() = "Write(const THStack& )";
  SetError(HistoConverter::kFail,"NO IMPLEMENTATION !!");
  return false;
}