In This Package:

DarkNoiseHistogram.cc

Go to the documentation of this file.

#include "DarkNoiseHistogram.h"
#include "Event/ReadoutHeader.h"
#include "Event/ReadoutPmtCrate.h"

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

#include <iomanip>
#include "TLine.h"
#include "TF1.h"
#include "TMath.h"
using namespace std;
using namespace DayaBay;

DarkNoiseHistogram::DarkNoiseHistogram(const string& name, ISvcLocator* svcloc)
  : GaudiAlgorithm(name, svcloc)
  , m_log(msgSvc(), name)
{
  declareProperty("Strategy", m_Strategy=0);
}

DarkNoiseHistogram::~DarkNoiseHistogram()
{
}

StatusCode DarkNoiseHistogram::initialize()
{
  m_log << MSG::DEBUG << "initialize()" << endreq;
  
  StatusCode status;
  status = service("THistSvc", m_thistsvc);
  if(status.isFailure() ){
    m_log << MSG::INFO << "Unable to retrieve pointer to THistSvc" << endreq;
    return status;
  }

  string path;
  ostringstream   oss1[NBOARD][NCHANNEL],oss2[NBOARD][NCHANNEL];
  for(int i=0; i<NBOARD; i++){
    if(i<5||i>17) continue;
    for(int j=0; j<NCHANNEL; j++){
      path="/DQMHist/Raw/";
      oss1[i][j]<<setw(2)<<setfill('0')<<i;
      oss2[i][j]<<setw(2)<<setfill('0')<<j;
      string hname = "h_DarkNoise_"+oss1[i][j].str()+"_"+oss2[i][j].str();
      h_DarkNoise[i][j] = new TH1F(hname.c_str(),hname.c_str(),250,-50,200);
      status = m_thistsvc->regHist((path+hname).c_str(),h_DarkNoise[i][j]);
      
      m_DarkRate[i][j] = 0;
    }
  }
  m_NTrigger = 0;

  //h_sum_Gain = new TH1F("h_sum_Gain","h_Gain",300,0,300);
  h_sum_DarkRate = new TH1F("h_sum_DarkRate","h_DarkRate",300,0,300);
  //h_sum_DarkRate = new TProfile("h_sum_DarkRate","h_DarkRate",300,0,300,0,30);
  h_sum_DarkNoise = new TProfile("h_sum_DarkNoise","h_sum_DarkNoise",300,0,300,0,300);
  //status = m_thistsvc->regHist("/DQMHist/Raw/h_sum_Gain", h_sum_Gain);
  status = m_thistsvc->regHist("/DQMHist/Raw/h_sum_DarkRate", h_sum_DarkRate);
  status = m_thistsvc->regHist("/DQMHist/Raw/h_sum_DarkNoise", h_sum_DarkNoise);
  if(status.isFailure()){
    m_log << MSG::INFO << "Unable to regist histograms" << endreq;
    return status;
  }
  return StatusCode::SUCCESS;
}

StatusCode DarkNoiseHistogram::execute()
{
  m_log << MSG::DEBUG << "execute()" << endreq;

  ReadoutHeader* roh = get<ReadoutHeader>("/Event/Readout/ReadoutHeader");
  if (roh == 0) {
    m_log << MSG::ERROR << "Cannot find readouth header." << endreq;
    return StatusCode::FAILURE;
  }

  m_log << MSG::DEBUG << "Readout: " << (Readout)(*(roh->readout())) << endreq;

  Readout* readout = const_cast<Readout*>(roh->readout());
  ReadoutPmtCrate* crate = dynamic_cast<ReadoutPmtCrate*>(readout);
  ReadoutPmtCrate::PmtChannelReadouts channels = crate->channelReadout();
  m_log << MSG::DEBUG << "#channels=" << crate->channelReadout().size()  << endreq;
  ReadoutPmtCrate::PmtChannelReadouts::iterator it;
  for (it=channels.begin(); it != channels.end(); it++) {
    m_log << MSG::DEBUG << it->first << " => " << it->second  << endreq;
    
    int board=it->second.channelId().board();
    int channel=it->second.channelId().connector()-1;
    //cout << "board:" << board << " channel:" << channel << endl;

    int firsthit[NBOARD][NCHANNEL];
    for(int i=0; i<NBOARD; i++){
      if(i<5||i>17) continue;
      for(int j=0; j<NCHANNEL; j++){
        firsthit[i][j] = 1;
      }
    }

    long fee_charge;
    unsigned long fee_time;
    for(int i=0;i<it->second.size();i++){
      fee_charge=it->second.adc(i);           
      fee_time=it->second.tdc(i);

      if(fee_time>1070 && fee_time<1170){//select dark noise before trigger
        m_DarkRate[board][channel]++;
      }

      if(it->second.peakAdcRange(i)==FeeGain::kLow) {
        fee_charge = (it->second.adc(i)-it->second.pedestal(i))*20;
      } else {
        fee_charge = it->second.adc(i)-it->second.pedestal(i);
      }
      if(fee_time>1050){//select dark noise before trigger
        h_DarkNoise[board][channel]->Fill(fee_charge);          
  h_sum_DarkNoise->Fill(board*16+channel,fee_charge);
      }
      if(m_Strategy==1){//include first hits after trigger
        if(firsthit[board][channel]&&fee_time<800){
          h_DarkNoise[board][channel]->Fill(fee_charge);
        }       
      }

      //after first hit
      firsthit[board][channel] = 0;
    }

  }
  m_NTrigger++;

  return StatusCode::SUCCESS;
}

StatusCode DarkNoiseHistogram::finalize()
{
  //each Tdc corresponding 1.5625ns
  double timeinterval = m_NTrigger*100*1.5625*1e-9;
  cout << "Time interval " << timeinterval << endl;

  m_log << MSG::DEBUG << "finalize()" << endreq;
  m_log << MSG::INFO << "# of triggers: " << m_NTrigger << endreq;

  TF1 *f1 = new TF1("f1", "gaus", 0, 200);
  for(int i=0; i<NBOARD;i++){       
    if(i<5||i>17) continue;
    for(int j=0; j<NCHANNEL; j++){       
      double darkrate = m_DarkRate[i][j]/timeinterval;
      double darkrate_err = TMath::Sqrt(m_DarkRate[i][j])/timeinterval;
      int id = i*16+j;
      if(darkrate>30000){//set maximum 30000
        //h_sum_DarkRate->SetBinContent(id,30000); 
        //h_sum_DarkRate->SetBinError(id,1000);
  //h_sum_DarkRate->Fill(id,30);
        h_sum_DarkRate->SetBinContent(id,30); 
        h_sum_DarkRate->SetBinError(id,1);

        m_log << MSG::INFO << "Dark Rate of (board, channel):(" 
              << i << "," << j 
              << ") is larger than 30K Hz! Set to 30K Hz in histogram." 
              << endreq; 
      }
      else{
  //h_sum_DarkRate->Fill(id,darkrate/1000.);
        h_sum_DarkRate->SetBinContent(id,darkrate/1000.); 
        h_sum_DarkRate->SetBinError(id,darkrate_err/1000.); 
  //cout << darkrate << "," << darkrate_err << ","
  //  << h_sum_DarkRate->GetBinContent(id) << "," 
  //  << h_sum_DarkRate->GetBinError(id) << endl;

        if(darkrate>10000){
          m_log << MSG::INFO << "Dark Rate of (board, channel):(" 
                << i << "," << j << ") is " << darkrate << " Hz."
                << endreq; 
        }
      }
      //------fit the gain-----------------
//      double gain = 0;
//      double gain_err = 0;
//      double nentry = h_DarkNoise[i][j]->GetEntries();
//      if(nentry<100){
//      m_log << MSG::INFO << "Only " << nentry 
//            << " entries, skip the gain fit."
//            << endreq;        
//      }
//      else{
//        double maxx = -50+h_DarkNoise[i][j]->GetMaximumBin();
//      //f1->SetRange(maxx-10, maxx+12);
//        f1->SetParameters(20, 7);
//        f1->SetLineColor(2);
//        h_DarkNoise[i][j]->Fit("f1", "R");
//
//      gain = f1->GetParameter(1);
//      gain_err = f1->GetParError(1);
//      }
//      if(gain>50||gain<0){//something wrong
//      h_sum_Gain->SetBinContent(id,0);
//      h_sum_Gain->SetBinError(id,0);
//      }
//      else{
//      h_sum_Gain->SetBinContent(id,gain);
//      h_sum_Gain->SetBinError(id,gain_err);
//      }
    }
  }
//  h_sum_Gain->SetXTitle("board*16+channel");
//  h_sum_Gain->SetYTitle("Adc/PE");
//  h_sum_Gain->SetMaximum(50);
//  h_sum_Gain->SetMinimum(0);

  h_sum_DarkRate->SetXTitle("board*16+channel");
  h_sum_DarkRate->SetYTitle("Dark Rate (kHz)");
  h_sum_DarkRate->SetMaximum(31);
  h_sum_DarkRate->SetMinimum(0);
  h_sum_DarkRate->GetXaxis()->SetRangeUser(80,280);
  
  h_sum_DarkNoise->SetXTitle("board*16+channel");
  h_sum_DarkNoise->SetYTitle("Dark Noise (ADC)");
  h_sum_DarkNoise->SetMinimum(0);
  h_sum_DarkNoise->GetXaxis()->SetRangeUser(80,280);
  
  delete f1;
  return StatusCode::SUCCESS;
}

---