Data.cxx

Go to the documentation of this file.

// $Id: Data.cxx,v 1.11 2009/10/09 20:26:51 jyuko Exp $

// C/C++
#include <cmath>
#include <iostream>
#include <sstream>

// ROOT
#include "TDirectory.h"
#include "TH1D.h"
#include "TH2F.h"

// Local
#include "PhysicsNtuple/Default.h"
#include "PhysicsNtuple/Hist/HistMan.h"

// Local
#include "Data.h"

using namespace std;

//-------------------------------------------------------------------------------------
Anp::Plot::Data::Data()
   :fPlot1d(true),
    fPlot2d(false),
    fPlotAll(true),
    fCorrSize(0),
    fFillFrac(0),
    fFillNbin(0),
    fFillSize(0),
    fDir1d(0),
    fDir2d(0),
    fDirCorr(0),
    fMap1d(),
    fMap2d(),
    fMapCorr(),
    fCorr()
{
}

//-------------------------------------------------------------------------------------
Anp::Plot::Data::~Data()
{
}

//-------------------------------------------------------------------------------------
void Anp::Plot::Data::Fill(DataIter ibeg, DataIter iend, const double weight)
{
   for(DataIter xit = ibeg; xit != iend; ++xit)
   {
      const short xkey  = xit -> Key();
      const float xdata = xit -> Data();
 
      if(fDir1d)
      {
         Fill(xkey, xdata, weight);
      }

      if(!fDir2d)
      {
         continue;
      }

      for(DataIter yit = xit + 1; yit != iend; ++yit)
      {
         const short ykey  = yit -> Key();
         const float ydata = yit -> Data();
         
         TH2 *h2 = Get(fMap2d, fDir2d, xkey, ykey);
         if(h2)
         {
            h2 -> Fill(xdata, ydata, weight);
         }
      }
   }

   if(fDirCorr && fCorrSize > 1)
   {      
      if(fCorr.Fill(ibeg, iend) >= fCorrSize)
      {
         const Corr::CorrMap cmap = fCorr.GetCorrCoef();
         
         for(Corr::CorrMap::const_iterator cit = cmap.begin(); cit != cmap.end(); ++cit)
         {
            const Corr::Key &key = cit -> first;
            const double corr = cit -> second;
            
            TH1 *hc = Get(fMapCorr, fDirCorr, key.Key1(), key.Key2());
            if(hc)
            {
               hc -> Fill(corr);
            }
         }

         fCorr.Reset();
      }
   }
}

//-------------------------------------------------------------------------------------
void Anp::Plot::Data::Set(const string &key, const unsigned int value)
{
   //
   // Configure internal parameters
   //

   if(key == "Plot1d")
   {
      fPlot1d = static_cast<bool>(value);
   }
   else if(key == "Plot2d")
   {
      fPlot2d = static_cast<bool>(value);
   }
   else if(key == "PlotAll")
   {
      fPlotAll = static_cast<bool>(value);
   }
   else if(key == "CorrSize")
   {
      fCorrSize = value;
   }
   else if(key == "FillFrac")
   {
      fFillFrac = value;
   }
   else if(key == "FillNbin")
   {
      fFillNbin = value;
   }
   else if(key == "FillSize")
   {
      fFillSize = value;
   }
   else
   {
      cerr << "Data::Set - unknown key: " << key << endl;
   }
}

//-------------------------------------------------------------------------------------
void Anp::Plot::Data::Set(TDirectory *dir)
{
   //
   // Create TDirectory objects for internal use
   //

   if(!dir)
   {
      return;
   }

   if(fPlot1d)
   {
      fDir1d = dir;
   }
   if(fPlot2d)
   {
      fDir2d = Anp::GetDir("plot2d", dir);
   }      
   if(fCorrSize)
   {
      fDirCorr = Anp::GetDir("corr", dir);
   }
}

//-------------------------------------------------------------------------------------
TH1* Anp::Plot::Data::Fill(const short key, const float data, const double weight)
{
   //
   // Find (create for a first time), fill and return TH1 histogram
   //

   TH1 *h = 0;
   
   Short1Map::iterator sit = fMap1d.find(key);
   if(sit == fMap1d.end())
   {
      //
      // Attempt to read histogram from histman, if it fails then make one here
      //
      
      h = HistMan::Instance().CreateTH1(key, ".");
      if(h)
      {
         Anp::SetDir(h, fDir1d);
         sit = fMap1d.insert(Short1Map::value_type(key, h)).first;
      }
      else if(fPlotAll && fFillSize > 9 && fFillNbin > 1)
      {
         //
         // Create TH1 histogram, store first fFillSize events 
         // and use them to determine histogram range.
         //      
         vector<DataItem<float, double> > &fvec = fFloat[key];

         if(fvec.size() < fFillSize)
         {
            fvec.push_back(Anp::DataItem<float, double>(data, weight));
         }
         else if(!fvec.empty())
         {
            stringstream hname;
            hname << "put_" << std::setw(5) << std::setfill('0') << key;
            
            //
            // Sort vector, then exclude lowest and highest 5% of events
            //
            unsigned int lpos =  (1*fvec.size())/100;
            unsigned int hpos = (99*fvec.size())/100;
            if(fFillFrac < 100)
            {
               lpos = (fFillFrac*fvec.size())/100;
               hpos = (static_cast<unsigned int>(100 - int(fFillFrac))*fvec.size())/100;
            }

            assert(lpos < fvec.size() && hpos < fvec.size());

            std::sort(fvec.begin(), fvec.end());
            const double xmin = fvec[lpos].Key();
            const double xmax = fvec[hpos].Key();

            //
            // Create TH1
            //
            if(!hname.str().empty())
            {
               Anp::Lock<Anp::Mutex> lock(Anp::GetROOTMutex());
               h = new TH1D(hname.str().c_str(), hname.str().c_str(), fFillNbin, xmin, xmax);
            }
            Anp::SetDir(h, fDir1d);

            //
            // Insert TH1 into the map, now we are done with this key
            //
            sit = fMap1d.insert(Short1Map::value_type(key, h)).first;
            
            for(vector<DataItem<float, double> >::iterator fit = fvec.begin(); fit != fvec.end(); ++fit)
            {
               h -> Fill(fit -> Key(), fit -> Data());
            }
            
            fvec.clear();
         }
      }
      else
      {
         sit = fMap1d.insert(Short1Map::value_type(key, 0)).first;
      }
   }
   else
   {
      //
      // Histogram might already exists, get TH1 pointer 
      //
      h = sit -> second;
   }

   //
   // For valid TH1 pointer fill histogram
   //
   if(h)
   {
      h -> Fill(data, weight);
   }

   return h;
}

//-------------------------------------------------------------------------------------
TH2* Anp::Plot::Data::Get(Short2Map &smap, TDirectory *dir, const short xkey, const short ykey)
{
   //
   // Find (create for a first time) and return TH2 histogram
   //
   const HistKey key(xkey, ykey);

   Short2Map::iterator sit = smap.find(key);
   if(sit == smap.end())
   {
      TH2 *h = HistMan::Instance().GetTH2(xkey, ykey);
      if(h)
      {
         Anp::SetDir(h, dir);
      }
      
      sit = smap.insert(Short2Map::value_type(key, h)).first;
   }      

   return sit -> second;
}

//-------------------------------------------------------------------------------------
TH1* Anp::Plot::Data::Get(CorrMap &cmap, TDirectory *dir, const short xkey, const short ykey)
{
   const HistKey key(xkey, ykey);

   CorrMap::iterator sit = cmap.find(key);
   if(sit == cmap.end())
   {
      stringstream hname, aname;
      hname << "key_" << xkey << "_and_key_" << ykey;
      aname << xkey << " and " << ykey << " correlation coefficient";

      TH1 *h = 0;
      if(!hname.str().empty())
      {
         Anp::Lock<Anp::Mutex> lock(Anp::GetROOTMutex());
         TH1 *h = new TH1D(hname.str().c_str(), hname.str().c_str(), 
                           100, -1.00001, 1.00001);      
         h -> GetXaxis() -> SetTitle(aname.str().c_str());
         h -> GetXaxis() -> CenterTitle();
      }
      Anp::SetDir(h, dir);
      
      sit = cmap.insert(CorrMap::value_type(key, h)).first;
   }      

   return sit -> second;
}

---