ZBinTemplateProcessor.h

#ifndef Z_BIN_TEMPLATE_PROCESSOR_H
#define Z_BIN_TEMPLATE_PROCESSOR_H 1

//LCIO
#include <lcio.h>

//MARLIN
#include "marlin/Processor.h"
//C++
#include <iostream>
#include <cmath>
#include <sstream>
//using namespace std;

//LCIO
#include <IMPL/LCCollectionVec.h>
#include <EVENT/Track.h>
#include <EVENT/TrackerHit.h>

// Aida
//#ifdef MARLIN_USE_AIDA // the whole processor does not make sense without aida
#include <AIDA/AIDA.h>
#include <marlin/AIDAProcessor.h>
#include <AIDA/IHistogramFactory.h>
#include <AIDA/IHistogram1D.h>
#include <AIDA/IHistogram2D.h>
// endif

#include <marlin/Global.h>

namespace AIDA{
  class IHistogram1D;
  class IHistogram2D;
}


namespace marlintpc
{

    template<unsigned int N1DHISTOS, unsigned int N2DHISTOS>
        class ZBinTemplateProcessor : public marlin::Processor
    {
      public:
            
        virtual Processor*  newProcessor() { throw Exception("You cannot create an instance of ZBinTemplateProcessor"); }
        
      virtual ~ZBinTemplateProcessor() {} 

      virtual void init() { 

        // usually a good idea to
        printParameters() ;
        
        // check parameter:
        // _nZBins >=0
        if (_nZBins < 0)
          {
            streamlog_out(WARNING3) << "invalid _nZBins < 0 ; setting to 0" 
                            << std::endl;
            _nZBins = 0;
          }

        // _minZ < _maxZ
        if ( _maxZ <= _minZ )
          {
            streamlog_out(WARNING3) << "invalid z range, setting to 0" 
                            << std::endl;
            
            _nZBins = 0;
          }

        // calculate width of a z bin
        if (_nZBins>0)
          _zBinWidth = (_maxZ - _minZ)/_nZBins;
        else
          _zBinWidth=1;  

        //resize histo arrays to correct number of z bins
        for (unsigned int i = 0 ; i < N1DHISTOS ; i++)
          {
            _1DHistos[i].resize(_nZBins +1);
          }

        for (unsigned int i = 0 ; i < N2DHISTOS ; i++)
          {
            _2DHistos[i].resize(_nZBins +1);
          }

        // get histograms from the AIDAProcessor
        
        // loop the z bins
        for (int zBin=0; zBin<=_nZBins ; zBin++)
          {
            // create unique names and titles for all histograms
            std::stringstream zRangeString;
            std::stringstream zIndexString;
            if (zBin < _nZBins)
              // no extention for the histo without z cut
              {
                zRangeString << " "<<_minZ + _zBinWidth*zBin
                             << " < z < "
                             << _minZ + _zBinWidth*(zBin+1);
                zIndexString << zBin ;
              }
            

            // now get all histograms in this z bin from the AIDAProcessor
            for (unsigned int histoIndex = 0 ; histoIndex < N1DHISTOS ; histoIndex ++)
              {
                streamlog_out(DEBUG3) << "_1DHistos["<<histoIndex<<"]["<<zBin<<"]" << std::endl;
                
                _1DHistos[histoIndex][zBin] =   marlin::AIDAProcessor::histogramFactory(this)->
                  createHistogram1D(  _histo1DBaseNames[histoIndex] + zIndexString.str(), 
                                      _histo1DTitles[histoIndex] + zRangeString.str()  , 
                                      _nBins1D[histoIndex], 
                                      _histoRanges1D[histoIndex].min,
                                      _histoRanges1D[histoIndex].max );
              }

            for (unsigned int histoIndex = 0 ; histoIndex < N2DHISTOS ; histoIndex ++)
              {
                
                _2DHistos[histoIndex][zBin] =   marlin::AIDAProcessor::histogramFactory(this)->
                  createHistogram2D(  _histo2DBaseNames[histoIndex] + zIndexString.str(), 
                                      _histo2DTitles[histoIndex] + zRangeString.str()  , 
                                      _nBins2D[histoIndex].nBinsX, 
                                      _histoRanges2D[histoIndex].xRange.min,
                                      _histoRanges2D[histoIndex].xRange.max,
                                      _nBins2D[histoIndex].nBinsY, 
                                      _histoRanges2D[histoIndex].yRange.min,
                                      _histoRanges2D[histoIndex].yRange.max );
              }

          }// for zBin

        _isInitialised=true;
      }


        
        struct HistoRange
        {
            float min;
            float max;
        };


        struct HistoRange2D
        {
            HistoRange xRange;
            HistoRange yRange;
        };

        struct HistoNBins2D
        {
            int nBinsX;
            int nBinsY;
        };

      protected:

      ZBinTemplateProcessor(const std::string & typeName )
        : marlin::Processor(typeName) , _isInitialised(false) {
        
        _description = "ZBinTemplateProcessor cannot be instantiated. Its a template class to derive from." ;
        
        registerProcessorParameter( "MaxZ" , 
                                    "MaxZ Maximum z value of highest bin in case of z binning (default: 260)",
                                    _maxZ,
                                    float(260) );
        registerProcessorParameter( "MinZ" , 
                                    "MinZ Minimum z value of lowest bin in case of z binning (default: 0)",
                                    _minZ,
                                    float(0) );
        registerProcessorParameter( "nZBins" , 
                                    "Number ob bins in z, 0 means one Histo without z cut (default: 10)",
                                    _nZBins,
                                    int(10) );
   
        if (_nZBins>0)
          _zBinWidth = (_maxZ - _minZ)/_nZBins;
        else
          _zBinWidth=1;  

        // initialize all vectors for the different histo types to the templated size
        // this should be done in the consturctor, it does not depend on processor parameters
        _1DHistos.resize(N1DHISTOS);
        _2DHistos.resize(N2DHISTOS);
        _histo1DBaseNames.resize(N1DHISTOS);
        _histo2DBaseNames.resize(N2DHISTOS);
        _histo1DTitles.resize(N1DHISTOS);
        _histo2DTitles.resize(N2DHISTOS);
        _histoRanges1D.resize(N1DHISTOS);
        _histoRanges2D.resize(N2DHISTOS);
        _nBins1D.resize(N1DHISTOS);
        _nBins2D.resize(N2DHISTOS);

        // now initialize the other variables for each histoIndex
        // for those parameters which don't depend on ProcessorParameters.
        // The zBin stuff has to be done in init()
        for (unsigned int i = 0 ; i < N1DHISTOS ; i++)
          {
            std::stringstream histoNameStream;
            histoNameStream << typeName<<"histo1D" << i << "_";
            //streamlog_out(DEBUG3) << histoNameStream.str() << std::endl;
      
            //      _1DHistos[i].resize(_nZBins +1);
            _histoRanges1D[i].min=0.;
            _histoRanges1D[i].max=1.;
            _nBins1D[i]=100;
            _histo1DTitles[i]=histoNameStream.str();
            _histo1DBaseNames[i]=histoNameStream.str();
          }
        for (unsigned int i = 0 ; i < N2DHISTOS ; i++)
          {
            std::stringstream histoNameStream;
            histoNameStream << typeName<<"histo2D" << i << "_";
 
            // the streamlog level is not set in the constructor yet, this is always printed
            //streamlog_out_T(DEBUG3) << histoNameStream.str() << std::endl;
 
            //      _2DHistos[i].resize(_nZBins +1);
            _histoRanges2D[i].xRange.min=0.;
            _histoRanges2D[i].xRange.max=1.;
            _histoRanges2D[i].yRange.min=0.;
            _histoRanges2D[i].yRange.max=1.;
            _nBins2D[i].nBinsX=100;
            _nBins2D[i].nBinsY=100;
            _histo2DTitles[i]=histoNameStream.str();
            _histo2DBaseNames[i]=histoNameStream.str();
          }

      }


        void fillHistogram(unsigned int histoIndex, float x, float z) { 
          if (!_isInitialised)
            throw Exception("ZBinTemplateProcessor is not initialised."
                            "Did you call ZBinTemplateProcessor::init() in your init function?");
          
          // don't confuse histoIndex with the zIndex
          //    unsigned int histoIndex = histoIndices1D[histo];
          
          if (_nZBins > 0)
            {
              // zBin of the track (from mean z if the track)
              int zBin = static_cast<int>( (z - _minZ) / _zBinWidth );
              
              if ( (zBin >= 0) && (zBin < _nZBins) )
                {
                  _1DHistos[histoIndex][zBin]->fill(x);
                } 
            } // in _nZBins >0
          
          // fill the histo without z cut
          _1DHistos[histoIndex][_nZBins]->fill(x);
        }
      

      void fillHistogram(unsigned int histoIndex, float x, float y, float z) { 

        if (!_isInitialised)
          throw Exception("ZBinTemplateProcessor is not initialised."
                          "Did you call ZBinTemplateProcessor::init() in your init function?");

        // don't confuse histoIndex with the zIndex
        //    unsigned int histoIndex = histoIndices2D[histo];

        if (_nZBins > 0)
          {
            // zBin of the track (from mean z if the track)
            int zBin = static_cast<int>( (z - _minZ) / _zBinWidth );
        
            if ( (zBin >= 0) && (zBin < _nZBins) )
              {
                _2DHistos[histoIndex][zBin]->fill(x,y);
              } 
          } // in _nZBins >0

        // fill the histo without z cut
        _1DHistos[histoIndex][_nZBins]->fill(x,y);
      }

            
        std::vector< std::vector<AIDA::IHistogram1D *> > _1DHistos; 
        std::vector< std::vector<AIDA::IHistogram2D *> > _2DHistos; 

        std::vector< std::string >  _histo1DBaseNames;
        std::vector< std::string >  _histo2DBaseNames;
        std::vector< std::string >  _histo1DTitles;
        std::vector< std::string >  _histo2DTitles;
        std::vector< HistoRange >   _histoRanges1D; 
        std::vector< HistoRange2D > _histoRanges2D; 
        std::vector< int > _nBins1D; 
        std::vector< HistoNBins2D > _nBins2D; 

        float _minZ; 
        float _maxZ; 
        
        int _nZBins; 
        
        float _zBinWidth; 

      private:
        bool _isInitialised;
 
    };
}

#endif // Z_BIN_TEMLATE_PROCESSOR_H