treestream.cc

//----------------------------------------------------------------------------
// File: treestream.cc
//
// Description: The classes itreestream and otreestream provide a convenient 
//              interface to the kind of Root-trees typically used in 
//              analyses in Dzero and CMS. All the boilerplate code for 
//              getting at the tree, its branches and its leaves is hidden
//              from the user.
//
//              The design is based on the observation that, ultimately, data
//              reduce to a set of name/value pairs in which the value 
//              can be regarded as a homogeneous sequence of doubles, floats 
//              or ints. Moreover, an arbitrarily complex graph of the 
//              name/value pairs is readily achieved with, well, name/value 
//              pairs! Name/value pair is just another name for variable.
// 
//              Thinking in terms of objects, in the context of high energy 
//              physics analysis, has yielded great complexity. On the other 
//              hand, thinking in terms of name/value pairs reduces the 
//              complexity enormously. The irreducible complexity, which 
//              surely exists, properly resides in the graph generated by 
//              the set of name/value pairs, that is, variables.
//
//              Of course, it is often convenient to package name/value pairs 
//              into objects to avail oneself of their useful behaviour. This 
//              may be useful, for example, to model a particularly 
//              complicated sub-graph of name/value pairs. Sometimes this is
//              useful if an object defines functions on the name/value pairs.
//              A  nice example of packaging name/value pairs for this reason
//              is the packaging of the name/value pairs, Px, Py, Pz and E 
//              into TLorentzVector objects, which provides a host of useful 
//              mathematical operations on them.
//
// Created: 19-Feb-2005 Harrison B. Prosper, based on 12-Mar-2001 version of
//                      the same.
//          14-Jul-2005 HBP Make even simpler!
//          13-Aug-2005 HBP ..and simpler still!
//          30-Nov-2005 HBP fix counter loading bug
//          31-Oct-2009 HBP allow use of regexes in branch names
//                      fix looping bug so operator[] works for Python
//          02-Oct-2010 HBP minor change to itreestream to handle vector types
//                      directly.
//          22-Nov-2010 HBP allow reading of multiple trees
//          22-Nov-2011 HBP handle reading/storing of strings
//$Revision: 1.5 $
//----------------------------------------------------------------------------
#ifdef PROJECT_NAME
#include <boost/regex.hpp>
#endif
#include <glob.h>

#include <map>
#include <vector>
#include <string>
#include <iostream>
#include <sstream>
#include <algorithm>
#include <typeinfo>
#include <cctype>
#include <cassert>


#include "TList.h"
#include "TKey.h"
#include "TFile.h"
#include "TLeaf.h"
#include "TDirectory.h"
#include "TTree.h"
#include "TClass.h"
#include "TChain.h"
#include "TString.h"
#include "TList.h"
#include "TIterator.h"
#include "TFriendElement.h"

#ifdef PROJECT_NAME
#include "PhysicsTools/TheNtupleMaker/interface/treestream.h"
#else
#include "treestream.h"
#endif
//----------------------------------------------------------------------------

#ifdef __WITH_CINT__
ClassImp(itreestream)
  ClassImp(otreestream)
#endif

  using namespace std;

// Status codes

const int kSUCCESS=0;
const int kBADOPEN=1;
const int kBADTREE=2;
const int kBADBRANCH=3;
int DEBUGLEVEL=0;

// Make following visible only to this compilation unit.

namespace 
{
  void fatal(string message)
  {
    cout << "** Error ** " << message << endl;
    exit(1);
  }

  void warning(string message)
  {
    cout << "** Warning ** " << message << endl;
  }

  string blank("                                                            ");

  void DBUG(string message, int level=1)
  {
    if ( getenv("DBtreestream") > 0 )
      DEBUGLEVEL = atoi(getenv("DBtreestream"));
    else
      DEBUGLEVEL = 0;
    if ( DEBUGLEVEL >= level ) cout << message << endl;
  }

  void split(string str, vector<string>& vstr)
  {
    vstr.clear();
    istringstream stream(str);
    while ( stream )
      {
        string strg;
        stream >> strg;
        if ( stream ) vstr.push_back(strg);
      }
  }
 
  // For given leaf, find maximum number of items.
  int getmaxsize(TLeaf* leaf)
  {
    int count = 0;
    TLeaf* leafcounter = leaf->GetLeafCounter(count);
    if ( leafcounter != 0 )
      // Variable length array
      return leafcounter->GetMaximum();
    else
      // Either fixed length array or a simple variable
      return leaf->GetLen();
  }

  // ----------------------------------------------------------------------
  // return value of given internal buffer
  // ----------------------------------------------------------------------
  template <class T>
  inline
  T
  invalue(Field* field, int index=0, int which=0)
  {
    if ( which == 0 )
      return static_cast<T>(field->leaf->GetValue(index));
    else
      {
        assert(index >= 0);
        FieldBuffer<T>* d = dynamic_cast<FieldBuffer<T>*>(field);
        if ( d == 0 ) fatal("getinvalue - dynamic_cast failed " + 
                            string(field->branch->GetName()));
        assert(index < (int)d->value.size());
        return d->value[index];
      }
  }

  template <class T>
  inline
  int
  insize(Field* field)
  {
    FieldBuffer<T>* d = dynamic_cast<FieldBuffer<T>*>(field);
    if ( d == 0 ) fatal("getinvalue - dynamic_cast failed " + 
                        string(field->branch->GetName()));
    return (int)d->value.size();
  }


  // For debug only
  double
  getinvalue(Field* field, int index=0, int which=0)
  {
    if ( which == 0 )
      return field->leaf->GetValue(index);
    else
      {
        double val = 0;
        switch(field->iotype)
          {
          case 'D':
            val = invalue<double>(field, index, 1);
            break;
	    
          case 'F':
            val = static_cast<double>(invalue<float>(field, index, 1));
            break;
	      
          case 'L':
            val = static_cast<double>(invalue<long>(field, index, 1));
            break;

          case 'I':
            val = static_cast<double>(invalue<int>(field, index, 1));
            break;

          case 'S':
            val = static_cast<double>(invalue<short>(field, index, 1));
            break;
	    
          case 'B':
            val = static_cast<double>(invalue<char>(field, index, 1));
            break;

          case 'l':
            val = static_cast<double>(invalue<unsigned long>(field, index, 1));
            break;
	    
          case 'i':
            val = static_cast<double>(invalue<unsigned int>(field, index, 1));
            break;

          case 's':
            val=static_cast<double>(invalue<unsigned short>(field, index, 1));
            break;

          case 'b':
            val = static_cast<double>(invalue<unsigned char>(field, index, 1));
            break;

          default:
            val = invalue<double>(field, index, 1);
            break;
          }
        return val;
      }
  }

  int
  getinsize(Field* field)
  {
    int size = 0;
    switch(field->iotype)
      {
      case 'D':
        size = insize<double>(field);
        break;
        
      case 'F':
        size = insize<float>(field);
        break;
        
      case 'L':
        size = insize<long>(field);
        break;

      case 'I':
        size = insize<int>(field);
        break;

      case 'S':
        size = insize<short>(field);
        break;

      case 'B':
        size = insize<char>(field);
        break;

      case 'C':
        size = insize<string>(field);
        break;

      case 'l':
        size = insize<unsigned long>(field);
        break;
        
      case 'i':
        size = insize<unsigned int>(field);
        break;

      case 's':
        size = insize<unsigned short>(field);
        break;

      case 'b':
        size = insize<unsigned char>(field);
        break;

      default:
        size = insize<double>(field);
        break;
      }
    return size;
  }

  // ----------------------------------------------------------------------
  // return value of given external buffer
  // T is type of external buffer
  // ----------------------------------------------------------------------
  template <class T>
  inline
  T
  exvalue(Field* field, int index=0)
  {
    T val;
    if ( field->isvector )
      {
        vector<T>* d = reinterpret_cast<vector<T>*>(field->address);
        if ( d == 0 ) fatal("exvalue - reinterpret_cast failed" + 
                            string(field->branch->GetName()));
        val = (*d)[index];
      }
    else
      {
        T* d = reinterpret_cast<T*>(field->address);
        if ( d == 0 ) fatal("exvalue - reinterpret_cast failed" + 
                            string(field->branch->GetName()));
        val = *d;
      }
    if ( DEBUGLEVEL > 0 )
      cout << "\texternal value(" << val << ")" << endl;
    return val;
  }

  template <class T>
  inline
  int
  exsize(Field* field)
  {
    int val = 1;
    if ( field->isvector )
      {
        vector<T>* d = reinterpret_cast<vector<T>*>(field->address);
        if ( d == 0 ) fatal("exvalue - reinterpret_cast failed" + 
                            string(field->branch->GetName()));
        val = d->size();
        //
        if ( DEBUGLEVEL > 1 )
          cout << "EXSIZE(" << val << ") " 
               << string(field->branch->GetName()) << endl;
      }
    return val;
  }


  string 
  getiotype(Field* field)
  {
    char iotype[2]; iotype[0] = field->iotype;  iotype[1] = 0;
    return string(iotype);
  }

  string 
  getsrctype(Field* field)
  {
    char srctype[2]; srctype[0] = field->srctype; srctype[1] = 0;
    return string(srctype);
  }

  // For debug only

  double
  getexvalue(Field* field, int index=0)
  {
    double d = 0;
    switch(field->srctype)
      {
      case 'D':
        d = exvalue<double>(field, index);
        break;
        
      case 'F':
        d = static_cast<double>(exvalue<float>(field, index));
        break;
        
      case 'L':
        d = static_cast<double>(exvalue<long>(field, index));
        break;
 
      case 'I':
        d = static_cast<double>(exvalue<int>(field, index));
        break;

      case 'S':
        d = static_cast<double>(exvalue<short>(field, index));
        break;
       
      case 'B':
        d = static_cast<double>(exvalue<char>(field, index));
        break;
       
      case 'l':
        d = static_cast<double>(exvalue<unsigned long>(field, index));
        break;
        
      case 'i':
        d = static_cast<double>(exvalue<unsigned int>(field, index));
        break;

      case 's':
        d = static_cast<double>(exvalue<unsigned short>(field, index));
        break;

      case 'b':
        d = static_cast<double>(exvalue<unsigned char>(field, index));
        break;
      }
    return d;
  }

  int
  getexsize(Field* field)
  {
    int size = 0;
    switch(field->srctype)
      {
      case 'D':
        size = exsize<double>(field);
        break;
        
      case 'F':
        size = exsize<float>(field);
        break;

      case 'L':
        size = exsize<long>(field);
        break;
              
      case 'I':
        size = exsize<int>(field);           

      case 'S':
        size = exsize<short>(field);
     
      case 'B':
        size = exsize<char>(field);

      case 'C':
        size = exsize<string>(field);
     
      case 'l':
        size = exsize<unsigned long>(field);
        break;
        
      case 'i':
        size = exsize<unsigned int>(field);
        break;
   
      case 's':
        size = exsize<unsigned short>(field);
        break;

      case 'b':
        size = exsize<unsigned char>(field);

      default:
        size = exsize<double>(field);
        break;
      }
    return size;
  }

  // ----------------------------------------------------------------------
  // Copy from external buffer to internal buffer.
  // Note: T is the type of the internal buffer
  // ----------------------------------------------------------------------

  template <class T>
  inline
  void
  fromexternal(Field* field, int count)
  {
    if ( DEBUGLEVEL > 0 )
      {
        char preamble[160];
        sprintf(preamble,
                "BEGIN fromexternal\n\t%-16s (%s <- %s) count = %d",
                field->leaf->GetName(), 
                getiotype (field).c_str(),
                getsrctype(field).c_str(),
                count);
        cout << preamble << endl;
      }        
    
    FieldBuffer<T>* d = dynamic_cast<FieldBuffer<T>*>(field);
    if ( d == 0 ) fatal("fromexternal - dynamic_cast failed " + 
                        string(field->branch->GetName()));
    
    switch(field->srctype)
      {
      case 'D':
        for(int i=0; i < count; i++)
          d->value[i] = static_cast<T>(exvalue<double>(field, i));
        break;
        
      case 'F':
        for(int i=0; i < count; i++)
          d->value[i] = static_cast<T>(exvalue<float>(field, i));
        break;
        
      case 'L':
        for(int i=0; i < count; i++)
          d->value[i] = static_cast<T>(exvalue<long>(field, i));
        break;

      case 'I':
        for(int i=0; i < count; i++)
          d->value[i] = static_cast<T>(exvalue<int>(field, i));
        break;

      case 'S':
        for(int i=0; i < count; i++)
          d->value[i] = static_cast<T>(exvalue<short>(field, i));
        break;
        
      case 'B':
        for(int i=0; i < count; i++)
          d->value[i] = static_cast<T>(exvalue<char>(field, i));
        break;

      case 'l':
        for(int i=0; i < count; i++)
          d->value[i] = static_cast<T>(exvalue<unsigned long>(field, i));
        break;
        
      case 'i':
        for(int i=0; i < count; i++)
          d->value[i] = static_cast<T>(exvalue<unsigned int>(field, i));
        break;

      case 's':
        for(int i=0; i < count; i++)
          d->value[i] = static_cast<T>(exvalue<unsigned short>(field, i));
        break;

      case 'b':
        for(int i=0; i < count; i++)
          d->value[i] = static_cast<T>(exvalue<unsigned char>(field, i));
        break;
        
      }
    
    if ( DEBUGLEVEL > 0 )
      {
        cout << "\tinternal value = " << d->value[0] << endl;
        cout << "\taddress        = " << &d->value[0] << endl;
        cout << "END fromexternal" << endl;
      }
  }


  template <>
  inline
  void
  fromexternal<string>(Field* field, int count)
  {
    if ( DEBUGLEVEL > 0 )
      {
        char preamble[160];
        sprintf(preamble,
                "BEGIN fromexternal\n\t%-16s (%s <- %s) count = %d",
                field->leaf->GetName(), 
                getiotype (field).c_str(),
                getsrctype(field).c_str(),
                count);
        cout << preamble << endl;
      }        
    
    FieldBuffer<string>* d = dynamic_cast<FieldBuffer<string>*>(field);

    if ( d == 0 ) fatal("fromexternal - dynamic_cast failed " + 
                        string(field->branch->GetName()));

    // IMPORTANT: use c_str() so that bytes are copied into existing
    // string.
    for(int i=0; i < count; i++)
      d->value[i] = exvalue<string>(field, i).c_str();
    
 
    if ( DEBUGLEVEL > 0 )
      {
        cout << "\tinternal value(" << d->value[0] << ")" << endl;
        void* address = (void*)(d->value[0].c_str());
        cout << "\tAddr(d->value[0].c_str()): " << address << endl;
        cout << "END fromexternal" << endl;
      }
  }

  // ----------------------------------------------------------------------
  // Copy to external buffer from internal buffer.
  // Note: T is the type of the external buffer
  //       Make sure external buffer is of the correct size
  // ----------------------------------------------------------------------
  template <class T>
  inline
  void
  toexternal(Field* field)
  {
    int count = min(field->leaf->GetLen(), field->maxsize);
    //int count = getmaxsize(field->leaf);

    if ( DEBUGLEVEL > 0 )
      {
        char preamble[160];
        sprintf(preamble,
                "BEGIN toexternal\n\t%-16s (%s -> %s) count = %d",
                field->leaf->GetName(), 
                getiotype (field).c_str(),
                getsrctype(field).c_str(),
                count);
        cout << preamble << endl;
      }        

    // Make sure external buffer is of the correct size

    if ( field->isvector )
      {
        vector<T>* d = reinterpret_cast<vector<T>*>(field->address);
        if ( d == 0 ) fatal("toexternal - reinterpret_cast failed" + 
                            string(field->fullname));
        
        if ( (int)d->size() != count ) d->resize(count, 0);
        
        if ( (int)d->size() != count ) 
          fatal("toexternal - unable to resize external buffer for " + 
                field->leafname);
        
        for(int index=0; index < count; index++)
          (*d)[index] = static_cast<T>(field->leaf->GetValue(index));

        if ( DEBUGLEVEL > 0 )
          cout << "\texternal value = " << (*d)[0] << endl;
      }
    else
      {
        T* d = reinterpret_cast<T*>(field->address);
        if ( d == 0 ) fatal("toexternal - reinterpret_cast failed" + 
                            string(field->fullname));
        
        *d = static_cast<T>(field->leaf->GetValue());

        if ( DEBUGLEVEL > 0 )
          {
            cout << "\texternal value = " << *d;
            if ( field->iscounter ) cout << " (counter)";
            cout << endl;
          }
      }
    if ( DEBUGLEVEL > 0 )
      cout << "END toexternal" << endl;
  }

  template <>
  inline
  void
  toexternal<string>(Field* field)
  {
    int count = field->leaf->GetLen();

    if ( DEBUGLEVEL > 0 )
      {
        char preamble[160];
        sprintf(preamble,
                "BEGIN toexternal\n\t%-16s (%s -> %s) maxsize = %d",
                field->leaf->GetName(), 
                getiotype (field).c_str(),
                getsrctype(field).c_str(),
                count);
        cout << preamble << endl;
      }        

    // Make sure external buffer is of the correct size

    string* d = reinterpret_cast<string*>(field->address);
    if ( d == 0 ) fatal("toexternal - reinterpret_cast failed" + 
                        string(field->fullname));
        
    d->assign(count, 0); // assign new content to string
    int size=0;
    for(int index=0; index < count; index++)
      {
        (*d)[index] = static_cast<char>(field->leaf->GetValue(index));
        if ( (*d)[index] == 0 ) break;
        size++;
      }
    d->resize(size);

    if ( DEBUGLEVEL > 0 )
      cout << "\texternal value = " << *d << endl;

    if ( DEBUGLEVEL > 0 )
      cout << "END toexternal" << endl;
  }

  // ----------------------------------------------------------------------
  // IMPORTANT: In tree->Branch(..), it is necessary to use the correct 
  // type when giving the address of the variable from which data are to 
  // be read, otherwise tree->Fill() gets really confused and can even hang
  // your computer!
  // ----------------------------------------------------------------------


  template <class T>
  inline
  void
  createbranch(TTree* tree, Field* field, 
               const char* format, 
               SelectedData& selecteddata)
  {
    DBUG("\tcreatebranch: BEGIN");

    if ( tree == 0 ) fatal("createbranch - tree pointer is zero ");
    if ( field == 0 ) fatal("createbranch - field pointer is zero ");

    FieldBuffer<T>* v = new FieldBuffer<T>();
    assert(v != 0);

    v->iotype  = field->iotype;
    v->srctype = field->srctype;
    v->address = field->address;
    v->isvector= field->isvector;
    v->maxsize = field->maxsize;
    v->branchname= field->branchname;
    v->leafname= field->branchname;
    v->iscounter = false;
    v->branch  = 0;
    v->leaf    = 0;
    v->leafname= ""; 

    v->value.clear();
    v->value.reserve(v->maxsize);
    v->value.resize(v->maxsize, 0);

    if ( DEBUGLEVEL > 0 )
      {
        cout << "\tcreatebranch: maxsize: " << v->value.size() << endl;
        cout << "\tcreatebranch: ADDRESS: " << v->address << endl;
      }

    // Store in map
    selecteddata[v->branchname] = v;

    DBUG("\tcreatebranch: " + v->branchname + "\t" + format);

    void* address = &(v->value[0]);
    tree->Branch(v->branchname.c_str(), address, format);


    TBranch* branch = tree->GetBranch(v->branchname.c_str());
    if ( branch == 0 ) fatal("createbranch - unable create branch " + 
                             v->branchname);

    v->branch = branch;
    v->leaf   = branch->GetLeaf(v->branchname.c_str());
    assert(v->leaf);

    DBUG("\tcreatebranch: END");
  }
 

  template <>
  inline
  void
  createbranch<string>(TTree* tree, Field* field, 
                       const char* format, 
                       SelectedData& selecteddata)
  {
    DBUG("\tcreatebranch: BEGIN");

    if ( tree == 0 ) fatal("createbranch - tree pointer is zero ");
    if ( field == 0 ) fatal("createbranch - field pointer is zero ");

    FieldBuffer<string>* v = new FieldBuffer<string>();
    assert(v != 0);

    v->iotype  = field->iotype;
    v->srctype = field->srctype;
    v->address = field->address;
    v->isvector= field->isvector;
    v->maxsize = field->maxsize;
    v->branchname= field->branchname;
    v->leafname= field->branchname;
    v->iscounter = false;
    v->branch  = 0;
    v->leaf    = 0;
    v->leafname= ""; 

    string* d = reinterpret_cast<string*>(field->address);
    if ( d == 0 ) fatal("createbranch - reinterpret_cast failed" + 
                        string(field->fullname));
    v->value.clear();
    v->value.reserve(v->maxsize);
    v->value.resize(v->maxsize, *d);

    if ( DEBUGLEVEL > 0 )
      {
        cout << "\tcreatebranch: init value(" << v->value[0] << ")" << endl;
        cout << "\tcreatebranch: maxsize:   " << v->value.size() << endl;
      }
    // Store in map
    selecteddata[v->branchname] = v;

    DBUG("\tcreatebranch: " + v->branchname + "\t" + format);

    void* address = (void*)(v->value[0].c_str());
    tree->Branch(v->branchname.c_str(), address, format);

    TBranch* branch = tree->GetBranch(v->branchname.c_str());
    if ( branch == 0 ) fatal("createbranch - unable create branch " + 
                             v->branchname);

    v->branch = branch;
    v->leaf   = branch->GetLeaf(v->branchname.c_str());
    assert(v->leaf);

    DBUG("\tcreatebranch: END");
  }


  void 
  readbranch(Field* field, int entry)
  {
    assert(field != 0);
    if ( field->branch ==0 ) return;

    assert(field->branch != 0);
    assert(field->leaf != 0);

    // Read entry for current branch

    field->branch->GetEntry(entry);

    // If address field is zero, this signals that 
    // the caller has not provided a location into which 
    // the value of the current variable is to be written.
    // This should happen only for leaf counter variables.
    // In this case, just return

    if ( field->address == 0) return;

    // If this is intrinsically a vector type, we let Root handle it
    // directly
    if ( field->iotype == 'v' ) return;

    // Copy data from internal to external buffers
    // iotype -> srctype

    switch(field->srctype)
      {
      case 'D':
        toexternal<double>(field);
        break;
	
      case 'F':
        toexternal<float>(field);
        break;
        
      case 'L':
        toexternal<long>(field);
        break;
        
      case 'I':
        toexternal<int>(field);
        break;
        
      case 'S':
        toexternal<short>(field);
        break;

      case 'B':
        toexternal<char>(field);
        break;

      case 'C':
        toexternal<string>(field);
        break;
	
      case 'l':
        toexternal<unsigned long>(field);
        break;
	
      case 'i':
        toexternal<unsigned int>(field);
        break;

      case 's':
        toexternal<unsigned short>(field);
        break;

      case 'b':
        toexternal<unsigned char>(field);
        break;
	
      default:
        toexternal<double>(field);
        break;
      }
  }
}


// Default constructor

itreestream::itreestream()
  : _tree(0),
    _chain(0),
    _statuscode(kSUCCESS),
    _current(-1),
    _entries(0),
    _entry(0),
    _index(0),
    _buffer(vector<double>(1000)),
    data(Data()),
    selecteddata(SelectedData()),
    _delete(true)
{}

itreestream::itreestream(string filename_, int bufsize)
  : _tree(0),
    _chain(0),
    _statuscode(kSUCCESS),
    _current(-1),
    _entries(0),
    _entry(0),
    _index(0),
    _buffer(vector<double>(bufsize)),
    data(Data()),
    selecteddata(SelectedData()),
    _delete(true)
{
  vector<string> fname;
  split(filename_, fname);
  vector<string> tname;
  _open(fname, tname);
}

itreestream::itreestream(vector<string>& fname, int bufsize)
  : _tree(0),
    _chain(0),
    _statuscode(kSUCCESS),
    _current(-1),
    _entries(0),
    _entry(0),
    _index(0),
    _buffer(vector<double>(bufsize)),
    data(Data()),
    selecteddata(SelectedData()),
    _delete(true)
{
  vector<string> tname;
  _open(fname, tname);
}

itreestream::itreestream(string filename_, string treename, int bufsize)
  : _tree(0),
    _chain(0),
    _statuscode(kSUCCESS),
    _current(-1),
    _entries(0),
    _entry(0),
    _index(0),
    _buffer(vector<double>(bufsize)),
    data(Data()),
    selecteddata(SelectedData()),
    _delete(true)
{
  vector<string> fname;
  split(filename_, fname);
  vector<string> tname;
  split(treename, tname);
  _open(fname, tname);
}

itreestream::itreestream(vector<string>& fname, string treename, int bufsize)
  : _tree(0),
    _chain(0),
    _statuscode(kSUCCESS),
    _current(-1),
    _entries(0),
    _entry(0),
    _index(0),
    _buffer(vector<double>(bufsize)),
    data(Data()),
    selecteddata(SelectedData()),
    _delete(true)
{
  vector<string> tname;
  split(treename, tname);
  _open(fname, tname);
}

void
itreestream::init(TTree* tree_)
{
  _delete = false;
  _tree = tree_;
  vector<string> fname;
  vector<string> tname;
  _open(fname, tname);
}

// ------------------------------------------------------------------------
// R E A D
// ------------------------------------------------------------------------
// Handle wildcarded filenames.
// 1. Use glob to get first filename
// 2. Get tree from file
// 3. Create a chain
// ------------------------------------------------------------------------
void
itreestream::_open(vector<string>& fname, vector<string>& tname)
{
  DBUG("itreestream::ctor - BEGIN", 1);

  // Clear internal buffer

  _buffer.clear();

  string treename("");
  if ( tname.size() > 0 ) treename = tname[0];

  // If tree pointer is zero, get tree from file
  if ( _tree == 0 )
    {
      // ----------------------------------------
      // Get list of files
      // Format of filename:
      // file1 file2 ..
      // ----------------------------------------
      filepath.clear();
#ifdef PROJECT_NAME
      for(int i=0; i < (int)fname.size(); ++i)
        {
          glob_t g;
          glob(fname[i].c_str(), GLOB_ERR | GLOB_NOCHECK, NULL, &g);
          for (int j=0; j < (int)g.gl_pathc; ++j)
            filepath.push_back(g.gl_pathv[j]);
          globfree(&g);
        }
#else
     for (unsigned int j=0; j < fname.size(); ++j){
       filepath.push_back(fname[j]);
     }
     //filepath.push_back(fname[0]);
      std::cout << "filepath size: " << filepath.size() << endl;  
#endif
      //DBUG("itreestream::ctor - new TFile ", 2);
      DBUG("itreestream::ctor - TFile::Open ", 2);
  
      // ----------------------------------------
      // Open first file
      // ----------------------------------------
      TFile* file_ = TFile::Open(filepath[0].c_str());
      //TFile* file_ = new TFile(filepath[0].c_str());
      if ( ! file_ || (file_ != 0 && ! file_->IsOpen()) )
        fatal("itreestream - unable to open file " + filepath[0]);
      file_->cd();
      
      if ( treename == "" )
        {      
          // ----------------------------------------
          // No tree name was given. Here is the default
          // action: If one of the trees is called 
          // Events then use it and warn user. If not,
          // use the first tree and warn user.
          // ----------------------------------------
          _tree = 0; // make sure to zero
          
          TIter nextkey(file_->GetListOfKeys());
          
          while ( TKey* key = (TKey*)nextkey() )
            {
              TObject* o = key->ReadObj();
              if ( o->IsA()->InheritsFrom("TTree") )
                {
                  if ( _tree == 0 ) _tree = (TTree*)o; // Record first tree
                  
                  if ( string(_tree->GetName()) == "Events" )
                    {
                      // Found a tree called Events, so use it
                      _tree = (TTree*)o;
                      break;
                    }
                }
            }
          if ( ! _tree )
            fatal("itreestream - NO tree found in file " + filepath[0]);
          
          treename = string(_tree->GetName());
          
          cout << endl << "** NB. itreestream - using tree: " 
               << treename << endl << endl;
        }
      else
        {
          _tree = (TTree*)file_->Get(treename.c_str());
          if ( ! _tree )
            fatal("itreestream - NO tree found in file " + filepath[0]);
        }
      
      string message("itreestream::ctor - treename: " + treename);
      DBUG(message, 2);

      // Remember to close file. It will be re-opened as part of a
      // chain.
      
      file_->Close();
      
      DBUG("itreestream::ctor - after file->Close", 2);

      // ----------------------------------------
      // Create a chain of files
      // ----------------------------------------
      // WARNING: This might be slow for large chains. 
      
      DBUG("itreestream::ctor - new TChain", 2);
      _chain = new TChain(treename.c_str());
      if ( ! _chain ) fatal("itreestream - Unable to create chain");

      _chainlist.push_back(_chain);
      // ----------------------------------------
      // Add possible friends
      // ----------------------------------------
      for(unsigned int i=1; i < tname.size(); i++)
        {
          DBUG("itreestream::ctor - AddFriend " + tname[i], 2);
          _chainlist.push_back(new TChain(tname[i].c_str()));
          _chain->AddFriend(_chainlist.back());
        }

      for(int i=0; i < (int)filepath.size(); i++)
        {
          for(unsigned int k=0; k < _chainlist.size(); k++)
            {
              _chainlist[k]->Add(filepath[i].c_str());
            }
        }

      DBUG("itreestream::ctor - GetEntries", 2);
      _entries = _chain->GetEntries();
      
      // ----------------------------------------
      // Update tree pointer
      // ----------------------------------------
      _tree = _chain; 

    }

  // ----------------------------------------
  // Get all branches and create a data
  // structure (Field) for every leaf. For 
  // leaves whose names differ from the branch, 
  // use the syntax:
  //    <branchname>.<leafname> 
  // ----------------------------------------
  
  for(unsigned int k=0; k < _chainlist.size(); k++)
    {
      TChain* chain = _chainlist[k];
      TObjArray* array = chain->GetListOfBranches();
      if ( !array ) fatal("itreestream::ctor - "
                          "Unable to GetListOfBranches for " +
                          string(chain->GetName()));

      int nitems = array->GetEntries();

      if ( DEBUGLEVEL > 1 )
        {
          char message[80];
          sprintf(message, "Number of branches in tree %s: %d", 
                  string(chain->GetName()).c_str(), nitems);
          cout << message << endl;
        }

      for (int i = 0; i < nitems; i++)
        {
          TBranch* branch = (TBranch*)((*array)[i]);      
          _getbranches(branch, 0);
        }

    }

  // ----------------------------------------
  // Flag branches that are leaf counters
  // ----------------------------------------
  Data::iterator it;
  for(it=data.begin(); it != data.end(); it++)
    {
      int count = 0;
      TLeaf* leafcounter = it->second.leaf->GetLeafCounter(count);
      if ( leafcounter != 0 )
        if ( data.find(leafcounter->GetBranch()->GetName()) != data.end() )
          data[leafcounter->GetBranch()->GetName()].iscounter = true;
    }

  if ( DEBUGLEVEL > 0 ) 
    cout << "itreestream::ctor - DATA.COUNT(" << data.size() << ")" << endl;

  DBUG("itreestream::ctor - exit OK", 1);
}

itreestream::~itreestream()
{
  _statuscode = kSUCCESS;
  close();
}

// ------------------------------------------------------------------------
// Get all branches from the root-file, recursively. Create a
// Field for each branch/leaf. By doing this recursively we do not have 
// to know the structure of the tree.
// ------------------------------------------------------------------------
void
itreestream::_getbranches(TBranch* branch, int depth)
{
  depth++;
  if ( depth > 10 ) fatal("_getbranches is lost in tree! Go boil your head!");

  /*
    branchname.push_back(branch->GetName());
    branchtab.push_back(depth);
    _getleaf(branch);            // Create leaf structure for current branch
  */

  // Check current branch for leaves

  TObjArray* array = 0;
  array = branch->GetListOfLeaves();
  if ( array != 0 )
    {
      int nleaves = array->GetEntries();
      for(int i=0; i < nleaves; i++)
        {
          TLeaf* leaf = (TLeaf*)((*array)[i]);
          if ( leaf != 0 )
            {
              _getleaf(branch, leaf); // Create leaf structure
              branchname.push_back(leaf->GetName());
              branchtab.push_back(depth);
            }
        }
    }
  
  // Check current branch for branches. For now, we assume branch names are
  // unique.

  array = 0;
  array = branch->GetListOfBranches();
  if ( array != 0 )
    {
      int nbranches = array->GetEntries();
      for (int i = 0; i < nbranches; i++)
        {
          TBranch* branch_ = (TBranch*)((*array)[i]);
          _getbranches(branch_, depth);
        }
    }
}

// ------------------------------------------------------------------------
// Create a leaf structure
// ------------------------------------------------------------------------
void
itreestream::_getleaf(TBranch* branch, TLeaf* leaf)
{
  if ( branch == 0 ) fatal("_getleaf - branch pointer zero");

  //bool storeBranchname = leaf == 0;

  if ( leaf == 0 )
    {
      TObjArray* array = branch->GetListOfLeaves();
      if ( array != 0 ) leaf = (TLeaf*)((*array)[0]);
    }
  if ( leaf == 0 )
    fatal("getleaf - leaf pointer zero for\n"
          "\tbranch " + string(branch->GetName()));

  Field v;
  v.address = 0;      // Set by caller
  v.srctype = 0;      // Set by caller
  v.maxsize = 0;      // Set by caller
  v.isvector= false;  // Set by caller (true if external variable is a vector)
  v.branch  = branch;
  v.leaf    = leaf;
  v.branchname= branch->GetName();
  v.leafname  = leaf->GetName();
  v.iotype    = leaf->GetTypeName()[0];
  v.iscounter = false;

  // If the branch contains a single leaf, it is sufficient
  // to identify the leaf with the branch. Otherwise, create
  // a name for each leaf within the branch.
  if ( branch->GetListOfLeaves()->GetEntries() == 1 )
    {
      v.fullname  = v.branchname;
      data[v.branchname] = v;
    }
  else
    {
      v.fullname  = v.branchname + "." + v.leafname;
      data[v.fullname] = v;
    }

  //   if ( storeBranchname )
  //     data[v.branchname] = v;
  //   else
  //     data[v.fullname] = v;

  if ( DEBUGLEVEL > 0 ) 
    cout << "getleaf(" << v.fullname << ")" << endl;
}

// Shutdown stream

void
itreestream::close()
{
  _statuscode = kSUCCESS;
  
  if ( _tree == 0 ) return;
  DBUG("itreestream::close file",3);
  if ( _delete ) delete  _tree;
  _tree = 0;
}

bool
itreestream::good() { return _statuscode == kSUCCESS; }

int
itreestream::status() { return _statuscode; }


void 
itreestream::select(string namen, double& d)
{
  _select(namen, &d, 1, 'D');
}

void 
itreestream::select(string namen, float& d)
{
  _select(namen, &d, 1, 'F');
}

void 
itreestream::select(string namen, long& d)
{
  _select(namen, &d, 1, 'L');
}

void 
itreestream::select(string namen, int& d)
{
  _select(namen, &d, 1, 'I');
}

void 
itreestream::select(string namen, short& d)
{
  _select(namen, &d, 1, 'S');
}

void 
itreestream::select(string namen, string& d)
{
  _select(namen, &d, 1, 'C');
}

void 
itreestream::select(string namen, unsigned long& d)
{
  _select(namen, &d, 1, 'l');
}

void 
itreestream::select(string namen, unsigned int& d)
{
  _select(namen, &d, 1, 'i');
}

void 
itreestream::select(string namen, unsigned short& d)
{
  _select(namen, &d, 1, 's');
}


// Now vectors

void 
itreestream::select(string namen, vector<double>& d)
{
  _select(namen, &d, d.size(), 'D', true);
}

void 
itreestream::select(string namen, vector<float>& d)
{
  _select(namen, &d, d.size(), 'F', true);
}

void 
itreestream::select(string namen, vector<long>& d)
{
  _select(namen, &d, d.size(), 'L', true);
}

void 
itreestream::select(string namen, vector<int>& d)
{
  _select(namen, &d, d.size(), 'I', true);
}

void 
itreestream::select(string namen, vector<short>& d)
{
  _select(namen, &d, d.size(), 'S', true);
}

void 
itreestream::select(string namen, vector<char>& d)
{
  _select(namen, &d, d.size(), 'B', true);
}

void 
itreestream::select(string namen, vector<unsigned long>& d)
{
  _select(namen, &d, d.size(), 'l', true);
}

void 
itreestream::select(string namen, vector<unsigned int>& d)
{
  _select(namen, &d, d.size(), 'i', true);
}

void 
itreestream::select(string namen, vector<unsigned short>& d)
{
  _select(namen, &d, d.size(), 's', true);
}

void 
itreestream::select(string namen)
{
  _bufmap[namen] = _index;
  _select(namen, &_buffer[_index], 1, 'D');
  _index++;
}

void 
itreestream::select(vector<string>& namen)
{
  _bufoffset =_index;
  _bufcount  = namen.size();
  for(unsigned i=0; i < namen.size(); i++) select(namen[i]);
}

vector<double>
itreestream::vget() 
{ 
  vector<double> v(_bufcount);
  for(int i=0; i < _bufcount; i++) v[i] = _buffer[i +_bufoffset];
  return v; 
}

double
itreestream::get(string namen) 
{ 
  if (_bufmap.find(namen) != _bufmap.end())
    return _buffer[ _bufmap[namen] ];
  else
    return 0; 
}

// ------------------------------------------------------------------------
// Read tree with ordinal value entry.
// ------------------------------------------------------------------------
int 
itreestream::read(int entry)
{
  _statuscode = kSUCCESS;
  int localentry = 0;

  // If entry is negative, we assume that the tree is already in
  // memory, in which case we do nothing.
  if ( entry > -1 )
    {
      _entry = entry;
      if ( _chain == 0 ) fatal("chain pointer is zero");

      // Load tree into memory

      localentry = _chain->LoadTree(entry);
      
      if (localentry < 0) return localentry;

      if ( DEBUGLEVEL > 0 ) 
        cout << "entry(" << entry << ")"
             << "localentry(" << localentry << ")" << endl;

      // Update pointers to tree, branches and leaves.
      
      if ( _chain->GetTreeNumber() != _current) _update();
    }
  else
    {
      localentry = _entry;

      // Update pointers to tree, branches and leaves.
      
      if ( _tree->GetTreeNumber() != _current) _update();

      _entry++;
    }
     
  // Copy data into external buffers

  SelectedData::iterator it;

  // IMPORTANT: Read leaf counters first...

  for(it=selecteddata.begin(); it != selecteddata.end(); it++)
    {
      Field* field = it->second;
      assert(field != 0);
      if ( ! field->iscounter ) continue;
      readbranch(field, localentry);
    }

  // ..then other variables

  for(it=selecteddata.begin(); it != selecteddata.end(); it++)
    {
      Field* field = it->second;
      assert(field != 0);
      if ( field->iscounter ) continue;
      readbranch(field, localentry);
    }

  return localentry; // Return ordinal value within current tree.
}

int 
itreestream::entries() { return _entries; }

int 
itreestream::size()    { return _entries; }

string
itreestream::name() { return _tree ? _tree->GetName() : ""; }

string
itreestream::title() { return _tree ? _tree->GetTitle() : ""; }

int 
itreestream::number() { return _current; }

vector<string>
itreestream::names() { return branchname; }

vector<string>
itreestream::filenames() { return filepath; }

string
itreestream::filename() { return filepath[_current]; }

TTree*
itreestream::tree() { return _tree; }

TFile*
itreestream::file() { return ((TChain*)(_tree))->GetFile(); }

string
itreestream::str() const
{
  if ( ! _tree ) return string("");

  ostringstream out;
  out << endl;
  out << "Tree Name          " << _tree->GetName()      << endl;
  out << "Number of entries  " << _tree->GetEntries()   << endl;
  out << endl;

  int count = 0;
  Data::const_iterator it;
  for(it=data.begin(); it != data.end(); it++)
    {
      const Field& field = it->second;
      count++;

      int maxcount=0;
      TLeaf* leafcounter = field.leaf->GetLeafCounter(maxcount);

      string lfsym("");
	  if (field.iscounter) lfsym = string(" *");

      char record[512];
      if ( leafcounter != 0 )
        {
          // This variable has a leaf counter
 //          sprintf(record, "%5d %s \t/ %s (%d) [%s]\n",
//                   count,
//                   field.fullname.c_str(),
//                   field.leaf->GetTypeName(),
//                   leafcounter->GetMaximum(),
//                   leafcounter->GetName());
          sprintf(record, "%5d %s \t/ %s (%d)\n",
                  count,
                  field.fullname.c_str(),
                  field.leaf->GetTypeName(),
                  leafcounter->GetMaximum());
        }
      else
        {
          // This variable does not have a leaf counter
          maxcount = field.leaf->GetLen();
          if ( maxcount > 1 )
            sprintf(record, "%5d %s \t/ %s [%d]\n",
                    count,
                    field.fullname.c_str(),
                    field.leaf->GetTypeName(), 
                    maxcount);
          else
            sprintf(record, "%5d %s\t/ %s%s\n",
                    count,
                    field.fullname.c_str(),
                    field.leaf->GetTypeName(),
                    lfsym.c_str());
        }
      out << record;
    }
  return out.str();
}

// ------------------------------------------------------------------------
// The following routine is used to specify the addresses into which data
// are to be read, when reading a Root-file
// Note 1: When itreestream is created a Field object is created for each 
//         branch/leaf in the root-tuple. However, a branch is activated 
//         only if the _select method has been called for that branch and
//         leaf and, of course, only if the named branch/leaf (that is, 
//         variable) exists.
//
// Note 2: The _select/read sequence does not use the internal buffers
//         field.value, but loads data directly into the user's
//         variables.
// ------------------------------------------------------------------------
void 
itreestream::_select(string namen, void* address, int maxsize, char srctype,
                     bool isvector)
{
  _statuscode = kSUCCESS;

  // If variable has already been selected, just update its address and
  // source type, otherwise get the branch and leaf.

  if ( DEBUGLEVEL > 0 )
    cout << "_select - " << namen << "\t" 
         << " address: " << address
         << " maxsize: " << maxsize
         << " srctype: " << srctype
         << " isvector: " << isvector
         << endl;

#ifdef PROJECT_NAME

  // If this is a regular expression, expand to full name of branch
  
  boost::regex e("[+&*?<>=()]|[[][^0-9]+[]]");
  boost::smatch w;
  
  if ( boost::regex_search(namen, w, e) )
    {
      // This seems like a regular expression

      TString tnom(namen.c_str());
      namen = tnom.ReplaceAll("*",".*") + "\\b"; // Add a boundary

      if ( DEBUGLEVEL > 0 )
        cout << "_select REGEX( " << namen << " )" << endl;

      vector<string> mnames;
      string nom = namen;
      
      if ( DEBUGLEVEL > 0 )
        cout << "_select DATA.COUNT( " << data.size() << " )" << endl;

      Data::iterator it;
      for(it=data.begin(); it != data.end(); it++)
        {
          Field& field = it->second;
          string name = field.fullname;

          if ( DEBUGLEVEL > 0 )
            cout << "_select(" << nom << ")[" << name << "]" << endl;

          boost::regex e2(nom);
          boost::smatch w2;
          if ( boost::regex_search(name, w2, e2) ) 
            {
              namen = name;
              mnames.push_back(namen);
              if ( DEBUGLEVEL > 0 )
                cout << "_select - EXPAND: " << nom << " --> " 
                     << namen << endl;
            }
        }

      // Fall on sword if regex selection is not unique

      if ( mnames.size() > 1 )
        {
          string rec("");
          for(int i=0; i < (int)mnames.size(); i++) rec += "\n\t" + mnames[i];
          fatal("_select - regex selection " + nom + 
                " is not unique:" + rec); 
        }
    }
#endif
 

  if ( DEBUGLEVEL > 0 )
    cout << "_select - NAMEN(" << namen << ")" << endl;

  // Check for new name
  if ( selecteddata.find(namen) != selecteddata.end() )
    {
      selecteddata[namen]->address = address;
      selecteddata[namen]->srctype = srctype;
    }
  else if ( data.find(namen) != data.end() )
    {
      // Known branch - activate it

      Field& field  = data[namen]; // NB: Get a reference, not a copy!
      field.srctype = srctype;
      field.maxsize = maxsize;
      field.address = address;     // source address
      field.isvector= isvector;
      selecteddata[namen] = &field;

      // If this branch has a branch counter, select its branch unless
      // already selected by caller.

      int count = 0;
      TLeaf* leafcounter = field.leaf->GetLeafCounter(count);
      if ( leafcounter != 0 )
        {
          // Current field has a leaf counter
          
          TBranch* branch = leafcounter->GetBranch();
          string name_(branch->GetName());
          
          if ( DEBUGLEVEL > 1 )
            cout << "_select - " << namen << " needs counter " 
                 << name_ << endl;
          
          // If leaf counter field not yet created, create it
          if ( selecteddata.find(name_) == selecteddata.end() )
            {
              if ( DEBUGLEVEL > 0 )
                cout << "_select - " << namen << " get counter " 
                     << name_ << endl;
              
              _getleaf(branch);

              // Do not set the address. If needed, it will be set by the 
              // caller in a subsequent call to the select method.
              // If the caller does not need it, then the address field 
              // should remain zero.
              
              Field& f = data[name_];  // NB: Get a reference, not a copy!
              f.srctype = 'I';         // Default source type
              f.maxsize = 1;
              f.iscounter = true;
              selecteddata[name_] = &f;
            }
        }
    }
  else
    {
      warning("itreestream - branch " + namen + " not found");
      _statuscode = kBADBRANCH;
    }
}  

// ------------------------------------------------------------------------
// Update the branch and leaf pointers. We do this when we switch from one 
// tree to another in a chain of Root-files.
// ------------------------------------------------------------------------
void 
itreestream::_update()
{
  _statuscode = kSUCCESS; 
  if ( _chain != 0 ) 
    _current = _chain->GetTreeNumber();
  else
    _current = _tree->GetTreeNumber();

  SelectedData::iterator it;
  for(it=selecteddata.begin(); it != selecteddata.end(); it++)
    {
      Field* field = it->second;
      if ( field == 0 )  fatal("update - zero field pointer");

      TBranch* branch = _tree->GetBranch(field->branchname.c_str());
      if ( branch == 0 )
        { 
          warning("update - pointer is zero for branch " 
                  + field->branchname);
          field->branch = 0;
          field->leaf   = 0;
          continue;
        }
      TLeaf* leaf = branch->GetLeaf(field->leafname.c_str());
      if ( leaf == 0 ) fatal("update - pointer is zero for leaf "
                             + field->leafname);

      field->branch = branch;
      field->leaf   = leaf;

      // We let Root handle vector types directly
      if ( field->iotype == 'v')
        {
          _chain->SetBranchAddress(field->branchname.c_str(), 
                                   &field->address, 
                                   &field->branch);
        }
      else if ( field->maxsize < 1 )
        fatal("_update - external buffer for " 
              + field->fullname + " is of zero length!");
    }
  DBUG("\tdone updating branch pointers", 1);
}

int 
itreestream::maximum(string name_)
{
  if ( data.find(name_) != data.end() )
    return getmaxsize(data[name_].leaf);
  else
    return 1;
}

bool
itreestream::present(string name_) { return data.find(name_) != data.end(); }

// int 
// itreestream::operator[](int entry) 
// {
//   int ev = read(entry);
//   if ( ev < 0 )
//     {
//       PyErr_SetString(PyExc_IndexError, "index out of range");
//       boost::python::throw_error_already_set(); // bail out
//     } 
//   return ev; 
// }

std::ostream& operator<<(std::ostream& os, const itreestream& tuple)
{
  os << tuple.str();
  return os;
}

// ------------------------------------------------------------------------
// W R I T E
// ------------------------------------------------------------------------
otreestream::otreestream()
  : _file(0),
    _tree(0),
    _statuscode(kSUCCESS),
    _entries(0)
{}

otreestream::otreestream(std::string filename, 
                         std::string treename, 
                         std::string title_,
                         int complevel,
                         int bufsize)
  : _file(0),
    _tree(0),
    _statuscode(kSUCCESS),
    _entries(0),
    _idatabuf(0),
    _databuf(vector<double>(bufsize)),
    _autosavecount(-1)
{
  DBUG("create file "+filename,1);

  _file = TFile::Open(filename.c_str(), "RECREATE");
  //_file = new TFile(filename.c_str(), "RECREATE");
  if ( ! _file )
    {
      cerr << "itreestream **Error** unable to open " << filename << endl;
      _statuscode = kBADOPEN;
      return;
    }
  _file->SetCompressionLevel(complevel);
  _file->cd();

  _tree = new TTree(treename.c_str(), title_.c_str());
  if ( ! _tree )
    {
      cerr << "itreestream **Error** can't make tree " << treename << endl;
      _statuscode = kBADTREE;
      return;
    }
}

otreestream::otreestream(TFile* file_, 
                         std::string treename, 
                         std::string title_,
                         int complevel,
                         int bufsize)
  : _file(file_),
    _tree(0),
    _statuscode(kSUCCESS),
    _entries(0),
    _idatabuf(0),
    _databuf(vector<double>(bufsize)),
    _autosavecount(-1)
{
  if ( ! _file )
    {
      cerr << "otreestream **Error** TFile pointer is zero" << endl;
      _statuscode = kBADOPEN;
      assert(_file);
      return;
    }
  _file->SetCompressionLevel(complevel);
  _file->cd();

  _tree = new TTree(treename.c_str(), title_.c_str());
  if ( ! _tree )
    {
      cerr << "otreestream **Error** can't make tree " << treename << endl;
      _statuscode = kBADTREE;
      return;
    }
}

otreestream::~otreestream()
{
  _statuscode = kSUCCESS;
  close();
}

void
otreestream::close(bool closefile)
{
  _statuscode = kSUCCESS;

  if ( _file == 0 ) return;
  
  DBUG("otreestream::close file (before)", 1);

  SelectedData::iterator it = selecteddata.begin();
  while ( it != selecteddata.end() ) 
    {
      delete it->second;
      it->second = 0;
      it++;
    }
  DBUG("otreestream::close file (after)", 1);

  if ( closefile )
    {
      _file->Write("", TObject::kOverwrite);
      _file->Close();

      delete _file;

      _tree = 0;
      _file = 0;
    }
}

bool
otreestream::good() { return _statuscode == kSUCCESS; }

int
otreestream::status() { return _statuscode; }

string
otreestream::name() { return _tree ? _tree->GetName() : ""; }

string
otreestream::title() { return _tree ? _tree->GetTitle() : ""; }

void
otreestream::add(string namen, double& datum, char iotype)
{

  // The last two variables are the source type and I/O type of variable.
  // We allow for the possibility that they may differ. For instance,
  // we may work with doubles but wish only to store floats.
  
  _add(namen, &datum, 1, 'D', iotype);
}

void
otreestream::add(string namen, float& datum)
{
  _add(namen, &datum, 1, 'F', 'F');
}

void
otreestream::add(string namen, long& datum)
{
  _add(namen, &datum, 1, 'L', 'L');
}

void
otreestream::add(string namen, int& datum)
{
  _add(namen, &datum, 1, 'I', 'I');
}

void
otreestream::add(string namen, short& datum)
{
  _add(namen, &datum, 1, 'S', 'S');
}

void
otreestream::add(string namen, string& datum)
{
  _add(namen, &datum, 1, 'C', 'C');
}

void
otreestream::add(string namen, unsigned long& datum)
{
  _add(namen, &datum, 1, 'l', 'l');
}

void
otreestream::add(string namen, unsigned int& datum)
{
  _add(namen, &datum, 1, 'i', 'i');
}

void
otreestream::add(string namen, unsigned short& datum)
{
  _add(namen, &datum, 1, 'l', 'l');
}

// Vectors

void 
otreestream::add(string namen, vector<double>& d, char iotype)
{
  _add(namen, &d, d.size(), 'D', iotype, true);
}

void 
otreestream::add(string namen, vector<float>& d)
{
  _add(namen, &d, d.size(), 'F', 'F', true);
}

void 
otreestream::add(string namen, vector<long>& d)
{
  _add(namen, &d, d.size(), 'L', 'L', true);
}

void 
otreestream::add(string namen, vector<int>& d)
{
  _add(namen, &d, d.size(), 'I', 'I', true);
}

void 
otreestream::add(string namen, vector<short>& d)
{
  _add(namen, &d, d.size(), 'S', 'S', true);
}

void 
otreestream::add(string namen, vector<char>& d)
{
  _add(namen, &d, d.size(), 'B', 'B', true);
}

void 
otreestream::add(string namen, vector<unsigned long>& d)
{
  _add(namen, &d, d.size(), 'l', 'l', true);
}

void 
otreestream::add(string namen, vector<unsigned int>& d)
{
  _add(namen, &d, d.size(), 'i', 'i', true);
}

void 
otreestream::add(string namen, vector<unsigned short>& d)
{
  _add(namen, &d, d.size(), 's', 's', true);
}

void 
otreestream::add(string namen)
{
  if ( (unsigned)_idatabuf >= _databuf.size() )
    {
      cout << "You need to increase internal buffer size" << endl;
      cout << "It is currently set to " << _databuf.size() << endl;
      exit(0);
    }
  _add(namen, &_databuf[_idatabuf], 1, 'D', 'D');
  _idatabuf++;
}

// ------------------------------------------------------------------------
// Copy each variable from external buffer to internal buffer, taking 
// care of type conversions. Then store data away.
// ------------------------------------------------------------------------
void
otreestream::insert(vector<double>& data) 
{ 
  copy(data.begin(), data.end(), _databuf.begin());
}

void
otreestream::store() 
{ 
  DBUG("BEGIN store");
  _statuscode = kSUCCESS;
  _entry = _entries;

  map<string, int> counter;

  SelectedData::iterator it;

  for(it = selecteddata.begin(); it != selecteddata.end(); it++)
    {
      Field* field  = it->second;
      assert(field);
      assert(field->leaf);

      DBUG("\tcommit - variable: " + string(field->leaf->GetName()));

      // Ok, get the count. If the variable has a leafcounter that we
      // know about, get the count set by the caller and use it if it
      // is smaller than the external buffer size.

      int size1 = getexsize(field);
      if ( DEBUGLEVEL > 0 )
        cout << "\tcommit - getexsize: " << size1 
             << "\tfield->maxsize: " << field->maxsize << endl;

      int count = min(size1, field->maxsize);
      int flag = 0;
      TLeaf* leafcounter = field->leaf->GetLeafCounter(flag);
      if ( leafcounter != 0 )
        {
          string name_(leafcounter->GetName());
          DBUG("\t\t- leaf counter: " + name_);

          if ( selecteddata.find(name_) != selecteddata.end() )
            //count = min(count, exvalue<int>(selecteddata[name_]));
            count = min(count, 
                        static_cast<int>(getexvalue(selecteddata[name_])));

          if ( DEBUGLEVEL > 1 )
            cout << "\t\t - with count: " << count << endl;
        }

      //if ( count == 0 ) fatal("exvalue: count is zero");

      // Copy data to internal from external buffers
      // iotype <- srctype

      switch(field->iotype)
        {
        case 'D':
          fromexternal<double>(field, count);
          break;
          
        case 'F':
          fromexternal<float> (field, count);
          break;
          
        case 'L':
          fromexternal<long>  (field, count);
          break;
          
        case 'I':
          fromexternal<int>   (field, count);
          break;

        case 'S':
          fromexternal<short> (field, count);
          break;

        case 'B':
          fromexternal<char> (field, count);
          break;

        case 'C':
          fromexternal<string> (field, count);
          break;

        case 'l':
          fromexternal<unsigned long>  (field, count);
          break;
          
        case 'i':
          fromexternal<unsigned int>   (field, count);
          break;

        case 's':
          fromexternal<unsigned short> (field, count);
          break;
          
        default:
          fromexternal<double>(field, count);
          break;
        }
    }
  DBUG("END store");
  DBUG("");
}

void
otreestream::save()
{
  // ..and store away.

  _file->cd();
  _tree->Fill();
  _entries++;

  //   // Save header every _autosavecount events
  //   if ( _entries % _autosavecount == 0 )
  //     {
  //       cout << _entries << "\totreestream::commit is saving header" << endl;
  //       _tree->AutoSave("SaveSelf");
  //     }
}

void
otreestream::commit()
{
  store();
  save();
}

void
otreestream::autosave(int count) 
{ 
  // count in Mbytes
  _autosavecount = count;
  if ( _autosavecount > 0 ) _tree->SetAutoSave(_autosavecount * 1000000);
}

int 
otreestream::entries() { return _entries; }

int 
otreestream::size()    { return _entries; }

vector<string>
otreestream::names() { return branchname; }

TFile*
otreestream::file() { return _file; }

TTree*
otreestream::tree() { return _tree; }

string 
otreestream::str() const
{
  ostringstream out;
  out << "Tree Name          " << _tree->GetName()      << endl;
  out << "Number of entries  " << _tree->GetEntries()   << endl;
  out << "Number of branches " << _tree->GetNbranches() << endl;
  return out.str();
}

void
otreestream::ls(ostream& out)
{
  out << str();
}

// ------------------------------------------------------------------------
// Create branch, of given type - int, float or double. The method can be 
// called repeatedly for the same branch name. On the first call the branch 
// is created and the address of the source variable is set. On subsequent 
// calls the address is reset.
// ------------------------------------------------------------------------
void 
otreestream::_add(string namen, void* address, int maxsize,
                  char srctype, char iotype, bool isvector)
{
  DBUG("_add: BEGIN");

  _statuscode = kSUCCESS;
  if ( maxsize < 1 )
    fatal("add - external buffer for " 
          + namen + " is of zero length!");

  int k = namen.find("/");
  if ( k > -1 ) namen = namen.substr(0, k);

  // Get format

  char format[256];
  char IOtype[2]; IOtype[0] = iotype; IOtype[1]=0;
  sprintf(format, "%s/%s", namen.c_str(), IOtype);


  // This could be a variable length array with name: <name>[countername].
  // if so, strip away [<countername>]
  k = namen.rfind("[");
  if ( k >= 0 )
    {
      namen = namen.substr(0,k);
    }

  // If branch is not present, create it; otherwise 
  // just set the source address.

  if ( selecteddata.find(namen) == selecteddata.end() )
    {      
      Field field;
      field.srctype = srctype;
      field.iotype  = iotype;
      field.address = address; // Source address
      field.maxsize = maxsize;
      field.isvector= isvector;
      field.branchname = namen;

      // Allocate a buffer of appropriate type and size
      // and create a corresponding branch

      switch(iotype)
        {
        case 'D':
          createbranch<double>(_tree, &field, format, selecteddata);
          break;

        case 'F':
          createbranch<float> (_tree, &field, format, selecteddata);
          break;
          
        case 'L':
          createbranch<long>  (_tree, &field, format, selecteddata);
          break;
          
        case 'I':
          createbranch<int>   (_tree, &field, format, selecteddata);
          break;
          
        case 'S':
          createbranch<short> (_tree, &field, format, selecteddata);
          break;

        case 'B':
          createbranch<char> (_tree, &field, format, selecteddata);
          break;

        case 'C':
          createbranch<string> (_tree, &field, format, selecteddata);
          break;

        case 'l':
          createbranch<unsigned long>  (_tree, &field, format, selecteddata);
          break;
          
        case 'i':
          createbranch<unsigned int>   (_tree, &field, format, selecteddata);
          break;
          
        case 's':
          createbranch<unsigned short> (_tree, &field, format, selecteddata);
          break;
    
        case 'b':
          createbranch<unsigned char> (_tree, &field, format, selecteddata);
          break;

        default:
          fatal(string("add - unsupported type ") + iotype);
        }
    }
  else
    // Update source address, because it may have changed.
    selecteddata[namen]->address = address;

  DBUG("_add: END");
}

std::ostream& operator<<(std::ostream& os, const otreestream& tuple)
{
  os << tuple.str();;
  return os;
}


//////////////////////////////////////////////////////////////////////////////

#ifdef __TEST__

#include <iomanip>
#include "TRandom3.h"

int main()
{
  cout << endl << endl << "Read/Write Test" << endl << endl;

  otreestream oustream("test.root", "Events", "Test");

  int njet = 0;
  oustream.add("njet", njet);

  vector<double> jetet(20);
  oustream.add("jetEt[njet]",jetet);

  string str(80, 0);
  oustream.add("str", str);

  TRandom3 rand;
  int entries = 10;
  int step    = 1;
  for(int entry=0; entry < entries; entry++)
    {
      njet = rand.Integer(10);
      jetet.clear();
      for(int i=0; i < njet; i++) 
        jetet.push_back(rand.Exp(10));

      char rec[80];
      string delim(entry+1, '-');
      sprintf(rec, "event: %d %s njet: %d", entry + 1, delim.c_str(), njet);
      str = string(rec);
 
      oustream.commit();

      if ( entry % step == 0 )
        cout << setw(5) << entry 
             << setw(5) << jetet.size()  
             << setw(10)<< jetet[0] 
             << " (" << str << ")"
             << endl;
    }
  oustream.close();


  itreestream stream("test.root", "Events");
  
  int nentries = stream.entries();
  //cout << endl << "Number of entries " << nentries << endl;
  cout << stream << endl;

  vector<float> JETET(20, -1);
  string STR;

  stream.select("jetEt", JETET);
  stream.select("str", STR);
  
  for(int entry=0; entry < nentries; entry++)
    {
      stream.read(entry);

      if ( entry % step == 0 )
        cout << setw(5) << entry 
             << setw(5) << JETET.size()  
             << setw(10)<< JETET[0]
             << " (" << STR << ")"
             << endl;
    }

  stream.close();
  return 0;
}

#endif