In This Package:

CalibDataHistogram.cc

Go to the documentation of this file.

// ******************************* 
// ***** 

#include "CalibDataHistogram.h"
#include "Event/CalibReadoutHeader.h"
#include "Event/CalibReadoutPmtCrate.h"

#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/IHistogramSvc.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/DataObject.h"
#include "GaudiKernel/ITHistSvc.h"

#include "StatisticsSvc/IStatisticsSvc.h"

#include "TH1F.h"
#include "TH2F.h"
#include "TStyle.h"
using namespace std;
using namespace DayaBay;
CalibDataHistogram::CalibDataHistogram(const string& name, ISvcLocator* svcloc)
  : GaudiAlgorithm(name, svcloc),m_statsSvc(0)
{
  declareProperty("PrintFreq",        m_printFreq=0);
  declareProperty("MaxTotalCharge",   m_maxTotalCharge=10000.);
  declareProperty("MinChannelCharge", m_minCharge=0.);
  declareProperty("MaxChannelCharge", m_maxCharge=3000.);
  declareProperty("MinChannelTime",   m_minTime=-2000.);
  declareProperty("MaxChannelTime",   m_maxTime=-500.);


  for(int i=0;i!=DETNUMPLUS;i++)
    for(int j=0;j!=RINGPLUS;j++)
      for(int k=0;k!=COLUMNPLUS;k++)
        {
          //kh_pmt[i][j][k]   = 0;
          h_charge[i][j][k] = NULL;
          h_time[i][j][k]   = NULL;
        }

  for(int i=0;i!=DETNUMPLUS;i++)
    {
      //kh_sum[i]                = 0;
      h_sum_FEEQSUM[i]         = NULL;
      h_sum_Charge_PMT[i]      = NULL;
      h_sum_Mean_Charge_PMT[i] = NULL;
      h_sum_RMS_Charge_PMT[i]  = NULL;
      h_sum_Time_PMT[i]        = NULL;
      h_sum_Mean_Time_PMT[i]   = NULL;
      h_sum_RMS_Time_PMT[i]    = NULL;
      h_sum_Charge_Ratio[i]    = NULL;
      h_sum_Charge_Patton[i]   = NULL;
    }

}

CalibDataHistogram::~CalibDataHistogram()
{
}

StatusCode CalibDataHistogram::initialize()
{
  debug() << "initialize()" << endreq;
  
  
  StatusCode status;

  status = service("StatisticsSvc", m_statsSvc, true);
  //  status = service("THistSvc", m_thistsvc);
  if(status.isFailure() )
    {
      info() << "Unable to retrieve pointer to THistSvc" << endreq;
      return status;
    }

  return StatusCode::SUCCESS;
}

StatusCode CalibDataHistogram::finalize()
{
  //for(unsigned int i=0;i<MAXRING;i++) {
  //  for(unsigned int j=0;j<MAXCOLUMN;j++) {
  //    h_sum_Mean_Charge_PMT->Fill(i*24+j,h_charge[i][j]->GetMean());
  //    h_sum_RMS_Charge_PMT->Fill(i*24+j,h_charge[i][j]->GetRMS());
  //    h_sum_Mean_Time_PMT->Fill(i*24+j,h_time[i][j]->GetMean());
  //    h_sum_RMS_Time_PMT->Fill(i*24+j,h_time[i][j]->GetRMS());
  //  }
  //}

  debug() << MSG::DEBUG << "finalize()" << endreq;

  //  if( m_statsSvc ) {
  //    m_statsSvc->release();
  //  }

  return StatusCode::SUCCESS;
}

StatusCode CalibDataHistogram::execute()
{
  debug() << "execute() ______________________________ begin" << endreq;
  
  
  CalibReadoutHeader* croh = get<CalibReadoutHeader>("/Event/CalibReadout/CalibReadoutHeader");
  if (croh == 0) 
    {
      error() << " =======> Requested Object can not be accessable." << endreq;
      return StatusCode::FAILURE;
    }
  
  CalibReadout* calibReadout = const_cast<CalibReadout*>(croh->calibReadout());
  CalibReadoutPmtCrate* calibCrate = dynamic_cast<CalibReadoutPmtCrate*>(calibReadout);
  CalibReadoutPmtCrate::PmtChannelReadouts calibChannels = calibCrate->channelReadout();
  
  //******************************* which detector == DetID
  const DayaBay::Detector detector = calibReadout->detector();
  
  //*******************************
  
  double qsum  = 0;
  double qmax  = 0;

  int ring     = 0;
  int column   = 0;

  double pmtid = 0;

  double sumCharge    = 0;
  double earliestTime = 0;

  CalibReadoutPmtCrate::PmtChannelReadouts::iterator it;

  unsigned int DetId   = detector.detectorId();

  for (it=calibChannels.begin(); it != calibChannels.end(); it++) 
    {
      AdPmtSensor sensorId(it->pmtSensorId().fullPackedData());
      ring = sensorId.ring();
      column = sensorId.column();
      
      stringstream path;
      path<<"/file1/";

      ostringstream ossSiteDetId;
      ossSiteDetId<<Site::AsString(detector.site())<<"/"<<DetectorId::AsString(detector.detectorId())<<"/";
      
      ostringstream ossPmtRing,ossPmtColumn;
      ossPmtRing<<setw(2)<<setfill('0')<<ring;
      ossPmtColumn<<setw(2)<<setfill('0')<<column;

      string pathDir = path.str()+ossSiteDetId.str();

      string pathPmtRingColumn = "Ring"+ossPmtRing.str()+"_"+"Column"+ossPmtColumn.str()+"/";

      string pathALL    = pathDir+pathPmtRingColumn;
      string pathSumDir = pathDir+"Calib/";

      if(h_charge[DetId][ring][column]==NULL)
        {
          h_charge[DetId][ring][column] = new TH1F("Charge","Charge",1000, 0, m_maxCharge);
          m_statsSvc->put((pathALL+"Charge").c_str(), h_charge[DetId][ring][column]);

          h_time[DetId][ring][column] = new TH1F("Time","Time",1000, m_minTime, m_maxTime);
          m_statsSvc->put((pathALL+"Time").c_str(), h_time[DetId][ring][column]);

          //kh_pmt[DetId][ring][column] = 1;
        }

      if(h_sum_FEEQSUM[DetId]==NULL)
        {
          h_sum_FEEQSUM[DetId]         = new TH1F("h_sum_FEEQSUM", "h_FEEQSUM", 5000 , 0, m_maxTotalCharge);
          m_statsSvc->put((pathSumDir+"h_sum_FEEQSUM").c_str(), h_sum_FEEQSUM[DetId]);


          h_sum_Charge_PMT[DetId]      = new TH2F("h_sum_ChargePMT", "h_sum_ChargePMT", 
                                                  MAXRING*MAXCOLUMN, 0.5, MAXRING*MAXCOLUMN+0.5, 3000, 0, m_maxCharge);
          m_statsSvc->put((pathSumDir+"h_sum_ChargePMT").c_str(), h_sum_Charge_PMT[DetId]);


          h_sum_Mean_Charge_PMT[DetId] = new TH1F("h_sum_MeanChargePMT", "h_sum_MeanChargePMT",
                                                  MAXRING*MAXCOLUMN, 0.5, MAXRING*MAXCOLUMN+0.5);
          m_statsSvc->put((pathSumDir+"h_sum_MeanChargePMT").c_str(), h_sum_Mean_Charge_PMT[DetId]);


          h_sum_RMS_Charge_PMT[DetId]  = new TH1F("h_sum_RMSChargePMT", "h_sum_RMSChargePMT",
                                                  MAXRING*MAXCOLUMN, 0.5, MAXRING*MAXCOLUMN+0.5);
          m_statsSvc->put((pathSumDir+"h_sum_RMSChargePMT").c_str(), h_sum_RMS_Charge_PMT[DetId]);


          h_sum_Time_PMT[DetId]        = new TH2F("h_sum_TimePMT", "h_sum_TimePMT", 
                                                  MAXRING*MAXCOLUMN, 0.5, MAXRING*MAXCOLUMN+0.5, 1000, m_minTime, m_maxTime);
          m_statsSvc->put((pathSumDir+"h_sum_TimePMT").c_str(), h_sum_Time_PMT[DetId]);


          h_sum_Mean_Time_PMT[DetId]   =new TH1F("h_sum_MeanTimePMT", "h_sum_MeanTimePMT",
                                                 MAXRING*MAXCOLUMN, 0.5, MAXRING*MAXCOLUMN+0.5) ;
          m_statsSvc->put((pathSumDir+"h_sum_MeanTimePMT").c_str(), h_sum_Mean_Time_PMT[DetId]);


          h_sum_RMS_Time_PMT[DetId]    = new TH1F("h_sum_RMSTimePMT", "h_sum_RMSTimePMT",
                                                  MAXRING*MAXCOLUMN, 0.5, MAXRING*MAXCOLUMN+0.5);
          m_statsSvc->put((pathSumDir+"h_sum_RMSTimePMT").c_str(), h_sum_RMS_Time_PMT[DetId]);


          h_sum_Charge_Ratio[DetId]    = new TH2F("h_sum_ChargeRatio", "h_sum_ChargeRatio",
                                                  5000 , 0, m_maxTotalCharge, 1000, 0, 1);
          m_statsSvc->put((pathSumDir+"h_sum_ChargeRatio").c_str(), h_sum_Charge_Ratio[DetId]);


          h_sum_Charge_Patton[DetId]   = new TH2F("h_sum_ChargePatton","h_sum_ChargePatton",49,0.25,24.75,19,-0.75,8.75);
          m_statsSvc->put((pathSumDir+"h_sum_ChargePatton").c_str(), h_sum_Charge_Patton[DetId]);
          

          //kh_sum[DetId] = 1;
        }
      
      
      if(ring==0) 
        {
          pmtid = 192+column-1;
        } 
      else 
        {
          pmtid = (ring-1)*24+column-1;
        }
      
      sumCharge = it->sumCharge();
      earliestTime = it->earliestTime();
      
      if(sumCharge > qmax) qmax = sumCharge;
      qsum += sumCharge;
      
      for(unsigned int i=0; i<it->size(); i++) 
        {

          h_charge[DetId][ring][column]->Fill(it->charge(i));
          h_time[DetId][ring][column]->Fill(it->time(i));

          h_sum_Charge_PMT[DetId]->Fill(pmtid,it->charge(i));
          h_sum_Time_PMT[DetId]->Fill(pmtid,it->time(i));

          //      h_charge->Fill(it->charge(i));
          //      h_time->Fill(it->time(i));
          
          //      h_sum_Charge_PMT_Det1st->Fill(pmtid,it->charge(i));
          //      h_sum_Time_PMT_Det1st->Fill(pmtid,it->time(i));
        }
    }

  
  h_sum_FEEQSUM[DetId]->Fill(qsum);
  if(qsum>0) h_sum_Charge_Ratio[DetId]->Fill(qsum,qmax/qsum);

  debug() << "execute() ______________________________ end" << endreq;
  return StatusCode::SUCCESS;
}

---