ampIntegralMatrix.h

Go to the documentation of this file.

//
//    Copyright 2010
//
//    This file is part of rootpwa
//
//    rootpwa 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 3 of the License, or
//    (at your option) any later version.
//
//    rootpwa 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 rootpwa. If not, see <http://www.gnu.org/licenses/>.
//
//-------------------------------------------------------------------------
// File and Version Information:
// $Rev::                             $: revision of last commit
// $Author::                          $: author of last commit
// $Date::                            $: date of last commit
//
// Description:
//      container class for complex amplitude integral matrices
//
//
// Author List:
//      Boris Grube          TUM            (original author)
//
//
//-------------------------------------------------------------------------


#ifndef AMPINTEGRALMATRIX_H
#define AMPINTEGRALMATRIX_H


#include <vector>
#include <map>
#include <complex>
#include <string>
#include <iostream>
#include <fstream>


#ifndef __CINT__
#define BOOST_DISABLE_ASSERTS
#include <boost/multi_array.hpp>
#endif

#include "TObject.h"

#include "waveDescription.h"


class TTree;
namespace rpwa {
        class amplitudeTreeLeaf;
}


namespace rpwa {


        class ampIntegralMatrix : public TObject {


                typedef std::map<std::string, unsigned int>::const_iterator waveNameToIndexMapIterator;


        public:
        

#ifndef __CINT__
                typedef boost::multi_array<std::complex<double>, 2> integralMatrixType;
        private:
                typedef integralMatrixType::size_type               sizeType;
        public:
#endif


                ampIntegralMatrix();
                ampIntegralMatrix(const ampIntegralMatrix& integral);
                virtual ~ampIntegralMatrix();

                void clear();

                ampIntegralMatrix& operator =(const ampIntegralMatrix& integral);

                friend bool operator ==(const ampIntegralMatrix& lhsInt,
                                        const ampIntegralMatrix& rhsInt);

                // arithmetic operators for integrals
                ampIntegralMatrix& operator +=(const ampIntegralMatrix& integral);
                ampIntegralMatrix& operator -=(const ampIntegralMatrix& integral);
                ampIntegralMatrix& operator *=(const double             factor);
                ampIntegralMatrix& operator /=(const double             factor);
                
                // accessors
                unsigned int  nmbWaves () const { return _nmbWaves;  }  
                unsigned long nmbEvents() const { return _nmbEvents; }  

                void setNmbEvents(const unsigned long nmbEvents) { _nmbEvents = nmbEvents; }  
                
                bool               containsWave(const std::string& waveName ) const;  
                unsigned int       waveIndex   (const std::string& waveName ) const;  
                const std::string& waveName    (const unsigned int waveIndex) const;  

                const std::vector<rpwa::waveDescription>& waveDescriptions() const { return _waveDescriptions; }  
                const waveDescription* waveDesc(const unsigned int waveIndex) const;  
                const waveDescription* waveDesc(const std::string& waveName ) const { return waveDesc(waveIndex(waveName)); }  
                bool allWavesHaveDesc() const;  

#ifndef __CINT__
                integralMatrixType&       matrix()       { return _integrals; }  
                const integralMatrixType& matrix() const { return _integrals; }  
#endif

                std::complex<double> element(const unsigned int waveIndexI,
                                             const unsigned int waveIndexJ) const;  
                std::complex<double> element(const std::string& waveNameI,
                                             const std::string& waveNameJ)  const  
                { return element(waveIndex(waveNameI), waveIndex(waveNameJ)); }


                bool integrate(const std::vector<std::string>& binAmpFileNames,
                               const std::vector<std::string>& rootAmpFileNames,
                               const unsigned long             maxNmbEvents   = 0,
                               const std::string&              weightFileName = "");

                void renormalize(const unsigned long nmbEventsRenorm);

                bool writeAscii(std::ostream&      out = std::cout) const;
                bool readAscii (std::istream&      in  = std::cin );
                bool writeAscii(const std::string& outFileName) const;
                bool readAscii (const std::string& inFileName );

                std::ostream& print(std::ostream& out,
                                    const bool    printIntegralValues = false) const;  

                static bool debug() { return _debug; }                             
                static void setDebug(const bool debug = true) { _debug = debug; }  


        private:

                unsigned long openBinAmpFiles(std::vector<std::ifstream*>&    ampFiles,
                                              const std::vector<std::string>& ampFileNames,
                                              const unsigned int              waveIndexOffset = 0);  
                
                unsigned long openRootAmpFiles(std::vector<TTree*>&                   ampTrees,
                                               std::vector<rpwa::amplitudeTreeLeaf*>& ampTreeLeafs,
                                               const std::vector<std::string>&        ampFileNames,
                                               const unsigned int                     waveIndexOffset = 0,
                                               const std::string&                     ampLeafName     = "amplitude");  

                void rebuildWaveNameToIndexMap();  

                bool hasIdenticalWaveSet(const ampIntegralMatrix& integral) const;  


                static bool _debug;  

                unsigned int                        _nmbWaves;            
                std::map<std::string, unsigned int> _waveNameToIndexMap;  
                std::vector<std::string>            _waveNames;           
                unsigned long                       _nmbEvents;           

                void storeMultiArray();  
                void readMultiArray ();  
#ifndef __CINT__
                integralMatrixType        _integrals;  
#endif
                std::vector<unsigned int> _intStorageShape;        
                unsigned int              _intStorageNmbElements;  
                std::complex<double>*     _intStorageData;         //[_intStorageNmbElements]

                std::vector<rpwa::waveDescription> _waveDescriptions;  


#ifdef USE_STD_COMPLEX_TREE_LEAFS
                ClassDef(ampIntegralMatrix,1)
#endif

        };


        // comparison operators
        inline
        bool
        operator ==(const ampIntegralMatrix& lhsInt,
                    const ampIntegralMatrix& rhsInt)
        {
                if (not lhsInt.hasIdenticalWaveSet(rhsInt)
                    or (lhsInt.nmbEvents() != rhsInt.nmbEvents()))
                        return false;
                for (unsigned int i = 0; i < lhsInt.nmbWaves(); ++i) {
                        const std::string waveNameI = lhsInt.waveName(i);
                        for (unsigned int j = 0; j < lhsInt.nmbWaves(); ++j) {
                                const std::string waveNameJ = lhsInt.waveName(j);
                                if (lhsInt.matrix()[i][j]
                                    != rhsInt.matrix()[rhsInt.waveIndex(waveNameI)][rhsInt.waveIndex(waveNameJ)])
                                        return false;
                        }
                }
                return true;
        }

        inline
        bool
        operator !=(const ampIntegralMatrix& lhsInt,
                    const ampIntegralMatrix& rhsInt)
        {
                return not(lhsInt == rhsInt);
        }


        // arithmetic operators for integrals
        inline
        ampIntegralMatrix
        operator +(const ampIntegralMatrix& integralA,
                   const ampIntegralMatrix& integralB)
        {
                ampIntegralMatrix result = integralA;
                result += integralB;
                return result;
        }

        inline
        ampIntegralMatrix
        operator -(const ampIntegralMatrix& integralA,
                   const ampIntegralMatrix& integralB)
        {
                ampIntegralMatrix result = integralA;
                result -= integralB;
                return result;
        }

        inline
        ampIntegralMatrix
        operator *(const ampIntegralMatrix& integral,
                   const double             factor)
        {
                ampIntegralMatrix result = integral;
                result *= factor;
                return result;
        }

        inline
        ampIntegralMatrix
        operator *(const double             factor,
                   const ampIntegralMatrix& integral)
        {
                return integral * factor;
        }

        inline
        ampIntegralMatrix
        operator /(const ampIntegralMatrix& integral,
                   const double             factor)
        {
                ampIntegralMatrix result = integral;
                result /= factor;
                return result;
        }


        inline
        std::ostream&
        operator <<(std::ostream&            out,
                    const ampIntegralMatrix& integral)
        {
                return integral.print(out);
        }


}  // namespace rpwa


#endif  // AMPINTEGRALMATRIX_H