PlotTime.cxx

Go to the documentation of this file.

// $Id: PlotTime.cxx,v 1.1 2008/03/31 17:33:25 rustem Exp $

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

// ROOT
#include "TDirectory.h"
#include "TH1.h"
#include "TH2.h"

// MINOS
#include "AstroUtil/Ast.h"
#include "AstroUtil/AstTime.h"
#include "Registry/Registry.h"
#include "Validity/VldTimeStamp.h"

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

// Local
#include "PlotTime.h"

REGISTER_ANP_OBJECT(AlgSnarl,PlotTime)

using namespace std;

//----------------------------------------------------------------------
Anp::HistTime::HistTime()
   :fPlot(false),
    fNMiss(0),
    fDay(0),
    fMonth(0),
    fYear(0),
    fDate(0),
    fSecs(0),
    fDir(0)
{
}

//----------------------------------------------------------------------
Anp::HistTime::HistTime(const VldTimeStamp &time)
   :fPlot(false),
    fNMiss(0),
    fDay(0),
    fMonth(0),
    fYear(0),
    fDate(0),
    fSecs(0),
    fDir(0)
{
   fDate = time.GetDate(true, 0, &fYear, &fMonth, &fDay);
   
   const VldTimeStamp btime(fYear, fMonth, fDay, 0, 0, 0);
   fSecs = btime.GetSec();
}

//----------------------------------------------------------------------
Anp::HistTime::~HistTime()
{
}

//----------------------------------------------------------------------
bool Anp::HistTime::Fill(const Record &record, const VldTimeStamp &time)
{
   //
   // Fill histograms if they have been created
   //
   if(!fPlot) return false;

   unsigned int hour = 0, min = 0, sec = 0;

   const int time_ = time.GetTime(true, 0, &hour, &min, &sec);

   Spill spill;
   if     (record.GetHeader().TorTgt() > 0.0) spill.protons = record.GetHeader().TorTgt();
   else if(record.GetHeader().TrTgtD() > 0.0) spill.protons = record.GetHeader().TrTgtD();
   else if(record.GetHeader().Tr101D() > 0.0) spill.protons = record.GetHeader().Tr101D();
   else if(record.GetHeader().Tor101() > 0.0) spill.protons = record.GetHeader().Tor101();
   
   spill.nevents = record.GetNEvents();
   spill.hour    = hour;
   spill.seconds = 3600*hour + 60*min + sec;

   return fSpill.insert(map<int, Spill>::value_type(time_, spill)).second;
}

//----------------------------------------------------------------------
bool Anp::HistTime::Make(TDirectory *dir)
{
   fPlot = false;

   if(!dir) return false;

   fDir = dir;

   if(fNMiss == 0)
   {
      fPlot = true;
   }
   else
   {
      cerr << "HistTime::Make - missed " << fNMiss << " histograms" << endl;
   }

   return fPlot;
}

//-----------------------------------------------------------------------------
TH1* Anp::HistTime::GetTH1(const std::string &key, const std::string &name)
{
   TH1 *h = HistMan::Instance().CreateTH1(key, "kinem");
   if(h)
   {
      Anp::SetDir(h, fDir, name);
   }
   else
   {
      ++fNMiss;
   }

   return h;
}

//----------------------------------------------------------------------
const pair<int, double> Anp::HistTime::GetSpills(const int hour)
{
   pair<int, double> spill(0, 0.0);

   for(map<int, Spill>::const_iterator sit = fSpill.begin(); sit != fSpill.end(); ++sit)
   {
      const Spill &spill_ = sit -> second;

      if(hour < 0 || hour == spill_.hour)
      {
         spill.first  += spill_.nevents;
         spill.second += spill_.protons;
      }
   }

   return spill;
}

//----------------------------------------------------------------------
void Anp::HistTime::FillSol(TH1 *h, const string &option)
{
   if(!h)
   {
      cerr << "HistTime::FillSol - invalid TH1 pointer" << endl;
      return;
   }

   const int tsec = 24*3600;
   for(map<int, Spill>::const_iterator sit = fSpill.begin(); sit != fSpill.end(); ++sit)
   {
      const Spill &spill = sit -> second;
      
      if(spill.seconds <= tsec)
      {
         if(option.find("protons") != string::npos)
         {
            h -> Fill(spill.seconds/double(tsec), spill.protons);
         }
         if(option.find("nevents") != string::npos)
         {
            h -> Fill(spill.seconds/double(tsec), spill.nevents);
         }
      }
      else
      {
         cerr << "HistTime::Fill - invalid time stamp: " << sit -> first << endl;
      }
   }
}

//----------------------------------------------------------------------
void Anp::HistTime::FillSid(TH1 *h, const string &option)
{
   if(!h)
   {
      cerr << "HistTime::FillSid - invalid TH1 pointer" << endl;
      return;
   }
   
   for(map<int, Spill>::const_iterator sit = fSpill.begin(); sit != fSpill.end(); ++sit)
   {
      const Spill &spill = sit -> second;
      const double hour = spill.seconds/3600.0;

      double juld = 0.0, gmst = 0.0, lsid = 0.0;

      AstUtil::CalendarToJulian(fYear, fMonth, fDay, hour, juld);
      AstUtil::JulianToGAST(juld, gmst);
      AstUtil::GSTToLST(gmst, AstUtil::kNearDetLongitude, lsid);

      if(lsid < 0.0 || lsid > 24.0)
      {
         cerr << "HistTime::Fill - invalid time stamp: " << sit -> first << endl;
      }
      else
      {
         if(option.find("protons") != string::npos)
         {
            h -> Fill(lsid/24.0, spill.protons);
         }
         if(option.find("nevents") != string::npos)
         {
            h -> Fill(lsid/24.0, spill.nevents);
         }
      }
   }
}

//----------------------------------------------------------------------
Anp::PlotTime::PlotTime()
   :fDirName("time"),
    fDir(0),
    fPlot(false),
    fNMiss(0),
    fNFail(0),
    fPotCut(1.0e2),
    fPotDay(1.0e5),
    fPotHour(1.0e4)
{
}

//----------------------------------------------------------------------
Anp::PlotTime::~PlotTime()
{
}

//----------------------------------------------------------------------
bool Anp::PlotTime::Run(Record &record)
{

   const Header &header = record.GetHeader();
   const VldTimeStamp time(header.Sec(), header.NSec());   
   const int date = time.GetDate();

   PlotMap::iterator hit = fMap.find(date);
   if(hit == fMap.end())
   {
      Handle<HistTime> hist(new HistTime(time));

      //
      // Create histograms
      //
      if(!(hist -> Make(fDir)))
      {
         hist.release();
      }

      hit = fMap.insert(PlotMap::value_type(date, hist)).first;
   }
   
   if((hit -> second).valid())
   {
      if(!(hit -> second -> Fill(record, time)))
      {
         ++fNFail;
      }
   }

   return true;
}

//----------------------------------------------------------------------
void Anp::PlotTime::Config(const Registry &reg)
{
   const char *value_char = 0;
   if(reg.Get("PlotTimeDirName", value_char) && value_char)
   {
      fDirName = value_char;
   }

   reg.Get("PlotTimePotCut",  fPotCut);
   reg.Get("PlotTimePotDay",  fPotDay);
   reg.Get("PlotTimePotHour", fPotHour);

   if(reg.KeyExists("PrintConfig"))
   {
      cout << "PlotTime::Config" << endl
           << "   DirName = " << fDirName << endl
           << "   PotCut = " << fPotCut << endl
           << "   PotDay = " << fPotDay << endl
           << "   PotHour = " << fPotHour << endl;
   }
}

//----------------------------------------------------------------------
void Anp::PlotTime::Set(TDirectory *dir)
{
   fPlot = false;

   if(!dir) return;

   fDir = Anp::GetDir(fDirName, dir);

   fRatePerDay   = PlotTime::GetTH1("rate_per_day");
   fRatePerHour  = PlotTime::GetTH1("rate_per_hour");
   fProtonsDay   = PlotTime::GetTH1("protons_day");
   fProtonsHour  = PlotTime::GetTH1("protons_hour");
   fProtonsvsSol = PlotTime::GetTH1("protons_vs_time", "protons_vs_sol");
   fProtonsvsSid = PlotTime::GetTH1("protons_vs_time", "protons_vs_sid");
   fNEventsDay   = PlotTime::GetTH1("nevents_day");
   fNEventsHour  = PlotTime::GetTH1("nevents_hour");
   fNEventsvsSol = PlotTime::GetTH1("nevents_vs_time", "nevents_vs_sol");
   fNEventsvsSid = PlotTime::GetTH1("nevents_vs_time", "nevents_vs_sid");

   if(fNMiss == 0)
   {
      fPlot = true;
   }
   else
   {
      cerr << "PlotTime::Set - missed " << fNMiss << " histograms" << endl;
   }
}

//----------------------------------------------------------------------
void Anp::PlotTime::End(const DataBlock &)
{
   if(!fDir || !fPlot)
   {
      return;
   }

   if(fNFail > 0)
   {
      cout << "PlotTime::End - failed to add " << fNFail << " snarls" << endl;
   }

   unsigned secs_beg = 0, secs_end = 0;
   for(PlotMap::iterator pit = fMap.begin(); pit != fMap.end(); ++pit)
   {
      Handle<HistTime> hist = pit -> second;

      if(pit == fMap.begin())
      {
         secs_beg = hist -> GetSecs();
         secs_end = hist -> GetSecs();
      }
      else
      {
         secs_beg = std::min<int>(secs_beg, hist -> GetSecs());
         secs_end = std::max<int>(secs_end, hist -> GetSecs());
      }
   }

   const int nhours = 24 + (secs_end - secs_beg)/3600;
   const int ndays  =  1 + (secs_end - secs_beg)/(24*3600);
   
   TH1 *hday  = 0;
   TH1 *hhour = 0;

   if(nhours > 1)
   {
      hhour = new TH1D("rate_per_hour_vs_time", "rate_per_hour_vs_time", nhours, 0, nhours);
      Anp::SetDir(hhour, fDir);
   }
   if(ndays > 1)
   {
      hday = new TH1D("rate_per_day_vs_time", "rate_per_day_vs_time", ndays, 0, ndays);
      Anp::SetDir(hday, fDir);
   }

   for(PlotMap::iterator pit = fMap.begin(); pit != fMap.end(); ++pit)
   {
      Handle<HistTime> hist = pit -> second;

      const double hour = double(hist -> GetSecs() - secs_beg)/3600.0;
      const double day  = double(hist -> GetSecs() - secs_beg)/(24.0*3600.0);

      const pair<int, double> aspill = hist -> GetSpills();

      if(aspill.first > 0 && aspill.second > fPotCut)
      {
         const double drate  = aspill.first*fPotDay/aspill.second;
         const double derror = drate*(0.01 + std::pow(aspill.first, -0.5));

         fNEventsDay -> Fill(aspill.first);
         fProtonsDay -> Fill(aspill.second);
         fRatePerDay -> Fill(drate);
         
         hist -> FillSol(fProtonsvsSol, "protons");
         hist -> FillSol(fNEventsvsSol, "nevents");

         hist -> FillSid(fProtonsvsSid, "protons");
         hist -> FillSid(fNEventsvsSid, "nevents");

         if(hday)
         {
            const int ibin = hday -> FindBin(day);
            hday -> SetBinContent(ibin, drate);
            hday -> SetBinError(ibin, derror);
         }
      }
      
      for(int ihour = 0; ihour < 24; ++ihour)
      {
         const pair<int, double> hspill = hist -> GetSpills(ihour);
         
         if(hspill.first > 0 && hspill.second > fPotCut)
         {
            const double hrate  = hspill.first*fPotHour/hspill.second;
            const double herror = hrate*(0.01 + std::pow(hspill.first, -0.5));

            fNEventsHour -> Fill(hspill.first);
            fProtonsHour -> Fill(hspill.second);
            fRatePerHour -> Fill(hrate);

            if(hhour)
            {
               const int ibin = hhour -> FindBin(hour + ihour);
               hhour -> SetBinContent(ibin, hrate);
               hhour -> SetBinError(ibin, herror);
            }
         }
      }
   }
}

//-----------------------------------------------------------------------------
TH1* Anp::PlotTime::GetTH1(const std::string &key, const std::string &name)
{
   TH1 *h = HistMan::Instance().CreateTH1(key, "time");
   if(h)
   {
      Anp::SetDir(h, fDir, name);
   }
   else
   {
      ++fNMiss;
   }

   return h;
}

---