In This Package:

RollingGain.cc

Go to the documentation of this file.

#include "RollingGain.h"
#include "Event/Readout.h"
#include "Event/CalibReadout.h"
#include "Event/CalibReadoutPmtCrate.h"
#include "Event/CalibReadoutPmtChannel.h"
#include "Event/RunData.h"
#include "DataSvc/ICableSvc.h"
#include "DataSvc/ICalibDataSvc.h"
#include "Event/RawEventHeader.h"
#include <math.h>
#include <fstream>
#include "TMath.h"

RollingGain::RollingGain(const std::string& name, ISvcLocator* pSvcLocator):
  GaudiAlgorithm(name, pSvcLocator)
{
  declareProperty("FileName", m_fileName = "AdcFit.root", "Root file name");
  declareProperty("FitPeriod", m_fitPeriod = 30*60, "Fit period in second, Default is 30 minutes");
  declareProperty("Fix", m_Fix = 0, "a temp fix for 40M noise");
}

RollingGain::~RollingGain()
{
}

double RollingGain::poisson(double u, int n){
  return pow(u,n)*exp(-u)/TMath::Gamma(n+1);
}

double RollingGain::gn(double x,double n, double q0, double sigma0,
          double q1, double sigma1){
  double sigman = sqrt(sigma0*sigma0+n*sigma1*sigma1);
  double xn = x-q0-n*q1;
  double tmp = pow(xn,2)/(2*sigman*sigman);
  return exp(-tmp)/(sigman*sqrt(2*TMath::Pi()));
}

double RollingGain::NIMmodel(double* xpar,double* par){
  double x =      xpar[0];
  double norm =   par[0];
  double q0 =     par[1];
  double sigma0 = par[2];
  double q1 =     par[3];  
  double sigma1 = par[4];
  double w =      par[5];
  double a =      par[6];
  double u =      par[7];
  int nmin =   1;
  int nmax =   6;

  double total = 0;
  for(int i=nmin; i<nmax; i++){
    int n = i;
    double tmp = poisson(u,n);
    double tmp2 = gn(x,n,q0,sigma0,q1,sigma1);
    double qn = q0+n*q1;
    double sigman = sqrt(sigma0*sigma0+n*sigma1*sigma1);

    double tmp3 = exp(-a*(x-qn-a*sigman*sigman/2.))*a/2.;
    double tmp4 = TMath::Erf(fabs(q0-qn-a*sigman*sigman)/(sigman*sqrt(2.)));
    double sign = 1.0;
    if((x-qn-sigman*sigman*a)<0) sign = -1.0;
    tmp4 += sign*TMath::Erf(fabs(x-qn-sigman*sigman*a)/(sigman*sqrt(2.)));
    tmp3 *= tmp4;
    
    double result = tmp*((1-w)*tmp2+w*tmp3);
    total += result;
  }
  return total*norm;
}

double RollingGain::Pedestal(double* xpar,double* par){
  double x =      xpar[0];
  double norm =   par[0];
  double q0 =     par[1];
  double sigma0 = par[2];
  double w =      par[3];
  double a =      par[4];

  double tmp = exp(-pow(x-q0,2)/(2*pow(sigma0,2)));    
  double result = (1-w)*tmp/(sigma0*sqrt(2*TMath::Pi()));
  if(x>q0) result += w*a*exp(-a*(x-q0));
  
  return result*norm;
}

void RollingGain::Fit(){
  ostringstream tmp;
  tmp <<setw(2) <<setfill('0') << m_NStop;
  string filename = "pmtDataTable_"+tmp.str()+".txt";
  FILE* outfile = fopen(filename.c_str(), "w+");
  
  fprintf(outfile, "# First Trigger time: (TrigSecond, TrigNanoSec) = (%d,%d).\n", m_FirstTrigSecond, m_FirstTrigNanoSec);
  fprintf(outfile, "# Last Trigger time: (TrigSecond, TrigNanoSec) = (%d,%d).\n", m_LastTrigSecond, m_LastTrigNanoSec);
  fprintf(outfile, "# [pmtID] [description] [status] [nhits] [fitStat] [chi2ndf] [speHigh] [gainErr] [sigmaSpe] [speLow] [timeOffset] [timeSpread] [efficiency] [prePulse] [afterPulse] [darkRate] [darkRateErr] [preAdc] [preAdcErr]\n");

  m_rootfile->cd();
  string newdir = tmp.str();
  m_rootfile->mkdir(newdir.c_str());
  m_rootfile->cd(newdir.c_str());

  m_masterfile = fopen("Master.pmtCalibMap.txt", "a");
  TimeStamp lastT1(m_LastTrigSecond,m_LastTrigNanoSec+1);
  fprintf(m_masterfile, "0000000 0000000 %4.0d %4.0d %10.0d %10.0d %10.0d %10.0d %s\n",
          m_site, m_simFlag, m_FirstTrigSecond,m_FirstTrigNanoSec,(int)lastT1.GetSec(),lastT1.GetNanoSec(),filename.c_str());
  fclose(m_masterfile);

  m_StartTime = m_FirstTrigSecond + m_FirstTrigNanoSec*1e-9;
  m_EndTime = m_LastTrigSecond + m_LastTrigNanoSec*1e-9;
  
  map<int/*PmtId*/,PmtProp>::iterator it,idend=m_PmtPropMap.end();
  for(it=m_PmtPropMap.begin(); it!=idend; it++) {
    PmtProp& curr = it->second;

    //calculate dark rate
    double tlength = m_NTrigger*100*1.5625*1e-9;//1 tdc corresponding 1.5625ns
    curr.DarkRateErr = sqrt(curr.DarkRate)/tlength;
    curr.DarkRate = curr.DarkRate/tlength;

    if( curr.h_PreAdc->GetEntries() == 0 ) {          

      curr.Status=PmtCalibData::kUnknown;
    } else if ( (curr.h_PreAdc->GetEntries() > 0 &&           
                 curr.h_PreAdc->GetEntries() < 100) ||
                (curr.h_Adc->GetEntries() > 0 &&
                 curr.h_Adc->GetEntries() < 100) ) {
      curr.Status = PmtCalibData::kBadEfficiency;
      curr.NHits = curr.h_Adc->GetEntries();
    } else {                                          
      //fit pedestal
      ostringstream tmp3;
      tmp3<<setw(2)<<setfill('0')<< m_NStop;
      string histname = "PreAdc_"+curr.PmtLabel+"_"+tmp3.str();
      string histname2 = "Adc_"+curr.PmtLabel+"_"+tmp3.str();
      
      double pedarea = curr.h_PreAdc->GetEntries();
      double pedmean = curr.h_PreAdc->GetMean();
      double pedsigm = curr.h_PreAdc->GetRMS();
      f1->SetParameters(pedarea,pedmean,pedsigm);
      f1->SetRange(pedmean-3*pedsigm,pedmean+3*pedsigm);
      int status1 = curr.h_PreAdc->Fit("f1","R");
      
      double q0 = f1->GetParameter(1);
      double q0err = f1->GetParError(1);
      double sigma0 = f1->GetParameter(2);
      //double w = f1->GetParameter(3);
      //double a = f1->GetParameter(4);
      if(f1->GetChisquare()/f1->GetNDF()>10 || status1!=0){
        warning()<< "Batch "<<tmp3.str()<<" bad pedestal fit for PMT: "
                 << curr.PmtLabel <<endreq;;
        q0 = curr.h_PreAdc->GetMean();
        q0err = 0;
        sigma0 = 1.9;
        curr.Status |= PmtCalibData::kBadPedestal;
      }
      curr.Pedestal = q0;
      //curr.h_PreAdc->Write(histname.c_str());
      //----------------------------------------
      //Adc is pedestal subtracted, set q0=0, sigma0=0
      //why sigma0 set to 0? because preAdc and Adc are highly correlated
      q0 = 0;
      sigma0 = 0;
      //----------------------------------------
      double AreaGuess, Q1MeanGuess, Q1SigmaGuess;
      AreaGuess = curr.h_Adc->GetEntries();
      Q1MeanGuess = curr.h_Adc->GetMean();
      Q1SigmaGuess = curr.h_Adc->GetRMS();
      
      f2->SetParameters(AreaGuess*10, 0, 0, Q1MeanGuess, Q1SigmaGuess, 0.01, 0.03, 0.1);
      f2->SetParNames("N","Q_{0}","#sigma_{0}","Q_{1}","#sigma_{1}","w","#alpha","#mu");
      
      f2->FixParameter(1,q0);
      f2->FixParameter(2,sigma0);
      f2->FixParameter(5,0.);
      f2->FixParameter(6,0.);
      //f2->FixParameter(5,w);
      //f2->FixParameter(6,a);
      
      f2->SetRange(q0+sigma0, q0+Q1MeanGuess+5*Q1SigmaGuess); 
      int status2 = curr.h_Adc->Fit("f2","R");
      
      double q1 = f2->GetParameter(3);
      double sigma1 = f2->GetParameter(4);
      double q1err = f2->GetParError(3);
      double sigma1err = f2->GetParError(4);
      double avechi2 = f2->GetChisquare()/f2->GetNDF();
      if(sigma1<0) sigma1 = -sigma1;//sometimes fit gives a negative sigma
      if( status2!=0 || avechi2>10 || q1<10 || q1>80){
        warning()<<"Batch "<<tmp3.str()<<" bad gain fit for PMT: "<<
          curr.PmtLabel<<", try a simple Gaussian fit." <<endreq;
        
        double adcmean = curr.h_Adc->GetMean();
        f1->SetParameters(AreaGuess,adcmean,7.);
        f1->SetRange(adcmean-10, adcmean+12);
        status2 = curr.h_Adc->Fit("f1","R");
        
        q1 = f1->GetParameter(1)-q0;
        sigma1 = f1->GetParameter(2);
        q1err = f1->GetParError(1);
        sigma1err = f1->GetParError(2);
        avechi2 = f1->GetChisquare()/f1->GetNDF();
        if(status2!=0 || avechi2>10 || q1<10 || q1>80){
          warning()<< "Batch "<<tmp3.str()<<" gaussian fit for PMT: "
                   << curr.PmtLabel
                   << " is not good either. Set the gain to 0." <<endreq;
          q1 = 0; sigma1 = 0; q1err = 0; sigma1err = 0;
        }
      }
      //q0 is  always 0 here
      //curr.Pedestal = q0; 
      curr.PedestalErr = q0err;
      curr.Gain = q1;
      curr.Sigma = sigma1;
      curr.GainErr = q1err;
      curr.SigmaErr = sigma1err;
      curr.NHits = AreaGuess;
      curr.FitStatus = status2;
      curr.Chi2 = avechi2;
      if( curr.FitStatus!=0 || curr.Chi2>10 || curr.Gain<10 || curr.Gain>80 ) {
        curr.Status |= PmtCalibData::kBadGain;
      }
      if( curr.Sigma/curr.Gain > 1 ) {
        curr.Status |= PmtCalibData::kBadGainWidth;
      }
      if( curr.Status == PmtCalibData::kUnknown) {
        curr.Status = PmtCalibData::kGood;
      }
    }
    
    
    if(curr.h_Adc->GetEntries()>0){//skip 0s
      fprintf(outfile, "%d %s %2d %6d %2d %6.2f %8.3f %8.3f  %7.3f %.1f %.1f %.1f %.1f %.1f %.1f %9.1f %9.1f %8.3f %8.3f\n", 
              curr.PmtId, curr.PmtLabel.c_str(), curr.Status, curr.NHits, curr.FitStatus,
              curr.Chi2, curr.Gain, curr.GainErr, curr.Sigma, curr.Gain/19 ,0.,0.,0.,0.,0.,curr.DarkRate, curr.DarkRateErr, curr.Pedestal, curr.PedestalErr);
      
      curr.h_PreAdc->Write();
      curr.h_Adc->Write();
      curr.h_Tdc->Write();
    }
  }

  //m_tree->Fill();
  fclose(outfile);
}

int RollingGain::InitPmtProp(const ServiceMode& svc)
{

  const vector<AdPmtSensor>& adPmts = m_cableSvc->adPmtSensors(svc);
  const vector<PoolPmtSensor>& poolPmts = m_cableSvc->poolPmtSensors(svc);
  
  PmtProp defPmtProp;
  defPmtProp.h_Adc=0;
  defPmtProp.h_Tdc=0;
  defPmtProp.h_PreAdc=0;
  defPmtProp.Pedestal=0;
  defPmtProp.PedestalErr=0;
  defPmtProp.DarkRate=0;
  defPmtProp.DarkRateErr=0;
  defPmtProp.Gain=0;
  defPmtProp.GainErr=0;
  defPmtProp.Sigma=0;
  defPmtProp.SigmaErr=0;
  defPmtProp.Status=PmtCalibData::kUnknown;
  defPmtProp.NHits=0;
  defPmtProp.FitStatus=-1;
  defPmtProp.Chi2=-1;
  defPmtProp.PmtLabel="";
  defPmtProp.PmtId=0;
  //define histograms
  m_PmtPropMap.clear();
  //  AD PMTs 
  for( unsigned int ct=0; ct<adPmts.size(); ct++ ) {
    AdPmtSensor sensor = adPmts[ct];
    int id = sensor.fullPackedData();
    m_PmtPropMap[id] = defPmtProp;
    // get label
    ostringstream sring, scolumn;
    sring<<setw(2)<<setfill('0')<<sensor.ring();
    scolumn<<setw(2)<<setfill('0')<<sensor.column();
    string label = sensor.detName()+"_R"+sring.str()+"_C"+scolumn.str();
    //
    string hname = "h_Adc_"+label;
    string hname2 = "h_PreAdc_"+label;
    string hname3 = "h_Tdc_"+label;
    m_PmtPropMap[id].h_Adc = new TH1F(hname.c_str(),hname.c_str(),350,-50,300);
    m_PmtPropMap[id].h_PreAdc = new TH1F(hname2.c_str(),hname2.c_str(),350,-50,300);
    m_PmtPropMap[id].h_Tdc = new TH1F(hname3.c_str(),hname3.c_str(),400,400,1200);
    m_PmtPropMap[id].PmtLabel=label;
    m_PmtPropMap[id].PmtId=id;
  }
  //  Water Poll PMTs
  for( unsigned int ct=0; ct<poolPmts.size(); ct++ ) {
    PoolPmtSensor sensor = poolPmts[ct];
    int id = sensor.fullPackedData();
    m_PmtPropMap[id] = defPmtProp;
    // get label
    ostringstream swall, sspot,sfacing;
    swall<<setw(2)<<setfill('0')<<sensor.wallNumber();
    sspot<<setw(2)<<setfill('0')<<sensor.wallSpot();
    sfacing<<setw(2)<<setfill('0')<<sensor.inwardFacing();
    string label = sensor.detName()+"_W"+swall.str()+"_S"+sspot.str()+"_F"+sfacing.str();
    // 
    string hname = "h_Adc_"+label;
    string hname2 = "h_PreAdc_"+label;
    string hname3 = "h_Tdc_"+label;
    m_PmtPropMap[id].h_Adc = new TH1F(hname.c_str(),hname.c_str(),350,-50,300);
    m_PmtPropMap[id].h_PreAdc = new TH1F(hname2.c_str(),hname2.c_str(),350,-50,300);
    m_PmtPropMap[id].h_Tdc = new TH1F(hname3.c_str(),hname3.c_str(),400,400,1200);
    m_PmtPropMap[id].PmtLabel=label;
    m_PmtPropMap[id].PmtId=id;
  }

  //map<int/*PmtId*/,PmtProp>::iterator it,idend=m_PmtPropMap.end();
  //for(it=m_PmtPropMap.begin(); it!=idend; it++) {
  //  PmtProp& curr = it->second;
  //}
  
  ResetPmtProp();//reset

  m_NStop = 1;
  m_startrun = -1;
  m_currentrun = -1;
  
  return 1;
}

int RollingGain::ResetPmtProp()
{
  map<int/*PmtId*/,PmtProp>::iterator it,idend=m_PmtPropMap.end();
  for(it=m_PmtPropMap.begin(); it!=idend; it++) {
    PmtProp& curr = it->second;    
    curr.h_Adc->Reset();
    curr.h_PreAdc->Reset();
    curr.h_Tdc->Reset();
    curr.Pedestal=0;
    curr.PedestalErr=0;
    curr.DarkRate=0;
    curr.DarkRateErr=0;
    curr.Gain=0;
    curr.GainErr=0;
    curr.Sigma=0;
    curr.SigmaErr=0;
    curr.Status=PmtCalibData::kUnknown;
    curr.NHits=0;
    curr.FitStatus=-1;
    curr.Chi2=-1;
  }

  m_NTrigger = 0;
  m_startrun = -1;

  return 1;
}

StatusCode RollingGain::initialize()
{
  StatusCode sc;
  sc = this->GaudiAlgorithm::initialize();
  if( sc.isFailure() ) {
    error() << "Base class initialization error" << endreq;
    return sc;
  }

  sc = service("RunDataSvc",m_runDataSvc);
  if( sc.isFailure() ) {
    error() << "Can't get RunDataSvc" << endreq;
    return sc;
  }

  sc = service("StaticCableSvc",m_cableSvc);
  if( sc.isFailure() ) {
    error() << "Can't get StaticCableSvc" << endreq;
    return sc;
  }

  f1 = new TF1("f1","gaus",20,250);
  //f1 = new TF1("f1",Pedestal,0,150,5);
  f2 = new TF1("f2",NIMmodel,20,300,8);

  m_rootfile = new TFile(m_fileName.c_str(),"recreate");
  //m_tree = new TTree("t", "My Tree");
  //m_tree->SetMaxTreeSize(200000000);  // 200M

  m_masterfile = fopen("Master.pmtCalibMap.txt", "a");
  fprintf(m_masterfile, "# [StartRun] [EndRun] [Sites] [SimFlag] [StartTimeSec] [StartTimeNanoSec] [EndTimeSec] [EndTimeNanoSec] [map file]\n");
  fclose(m_masterfile);

  return sc;
}

StatusCode RollingGain::execute()
{
  debug()<<"RollingGain executing ... "<<endreq;

  StatusCode sc;
  
//  RawEventHeader* reh = get<RawEventHeader>(RawEventHeader::defaultLocation());
//  if(!reh) {
//    warning()<<"Failed to get raw event header"<<endreq;
//  }
//  else{
//    if(m_startrun==-1) m_startrun = reh->runNum();
//    m_currentrun = reh->runNum();
//  }

  ReadoutHeader* readoutHdr = get<ReadoutHeader>( ReadoutHeader::defaultLocation() );
  if(!readoutHdr) {
    error()<<"Failed to get readout header"<<endreq;
    return StatusCode::FAILURE;
  }

  const DaqCrate* daqCrate = readoutHdr->daqCrate();
  if(!daqCrate) {
    error()<<"Failed to get daqCrate from header"<<endreq;
    return StatusCode::FAILURE;
  }

  TimeStamp trigTime = daqCrate->triggerTime();
  m_currentTime = trigTime.GetSeconds();

  static bool first = true;
  if( first ) {
    first = false;
    Context ctx = readoutHdr->context();
    m_site=ctx.GetSite();
    m_simFlag=ctx.GetSimFlag();
    ctx.SetSite(Site::kAll);
    ctx.SetDetId(DetectorId::kAll);
    ServiceMode svc(ctx,0);
    //
    InitPmtProp(svc);
    //
    m_prevTrigTime = trigTime.GetSeconds();
    m_lastFitTime = trigTime.GetSeconds();
  }

  if( !(daqCrate->detector().isAD() ||
        daqCrate->detector().isWaterShield()) ) {
    return StatusCode::SUCCESS;
  }

  const DaqPmtCrate* pmtCrate = daqCrate->asPmtCrate();
  if (!pmtCrate) {
    error()<<"Failed to get DaqPmtCrate"<<endreq;
    return StatusCode::FAILURE;
  }
  m_NTrigger++;

  // Keep tracking the first trigger time and 
  // the laster trigger for this group of fit
  m_LastTrigSecond = trigTime.GetSec();
  m_LastTrigNanoSec = trigTime.GetNanoSec();
  if(m_NTrigger==1){
    m_FirstTrigSecond = trigTime.GetSec();
    m_FirstTrigNanoSec = trigTime.GetNanoSec();
  }

  Detector det = pmtCrate->detector();
  ServiceMode svcMode( readoutHdr->context(),0 );
  //  cout << "run number:" << m_runDataSvc->runData(svcMode)->runNumber() << endl;
  
  //int TrigType = pmtCrate->triggerType();
  //if( TrigType != 7 ) {

  const DaqPmtCrate::PmtChannelPtrList& channels = pmtCrate->pmtChannelReadouts();
  DaqPmtCrate::PmtChannelPtrList::const_iterator ci, ci_end = channels.end();

  for(ci = channels.begin(); ci!=ci_end; ci++) {
    const DaqPmtChannel* channel = *ci;
    FeeChannelId channelId = channel->channelId();
    
    // number of hits in this readout
    unsigned int multi = channel->hitCount();
    // ring and column
    AdPmtSensor pmtId = m_cableSvc->adPmtSensor( channelId, svcMode );
    int id = pmtId.fullPackedData();
    PmtProp& curPmt = m_PmtPropMap[id];

    for(unsigned int nr = 0; nr<multi; ++nr) {
      int tdc = channel->tdc(nr);
      int adc = channel->adc(nr);
      float preAdc = channel->preAdcAvg(nr);
      int gain = channel->isHighGainAdc(nr);

      curPmt.h_Tdc->Fill(tdc);
      
      if(tdc>1070 && tdc<1170){//select dark noise before trigger
        curPmt.DarkRate++;
      }

      if( multi==1 &&                               
          tdc>1070 &&                               
          gain==FeeGain::kHigh &&                  
          m_currentTime - m_prevTrigTime > 20e-6     
          ) {
        if(m_Fix){                //  special treat for run 5773, lose one third data.
          int residue = (tdc-15)%16;
          if(residue>4){
            curPmt.h_Adc->Fill(adc-preAdc);
            curPmt.h_PreAdc->Fill(preAdc);
          }
        }
        else{
          curPmt.h_Adc->Fill(adc-preAdc);
          curPmt.h_PreAdc->Fill(preAdc);
        }
      }
    }
  } // end of channels

  //-------------------------------------
  //fit the data
  if(m_currentTime - m_lastFitTime > m_fitPeriod) {
    //Fit
    Fit();
    //after fit
    ResetPmtProp();
    m_NStop++;
    m_lastFitTime = m_currentTime;
  }
  
  m_prevTrigTime = m_currentTime;

  return StatusCode::SUCCESS;
}

StatusCode RollingGain::finalize()
{
  //Fit
  Fit();
  
  map<int/*PmtId*/,PmtProp>::iterator it,idend=m_PmtPropMap.end();
  for(it=m_PmtPropMap.begin(); it!=idend; it++) {
    PmtProp& curr = it->second;
    delete curr.h_Adc;
    delete curr.h_PreAdc;
    delete curr.h_Tdc;
  }

  m_rootfile->Write();
  m_rootfile->Close();

  delete f1;
  delete f2;
  //delete m_tree;
  delete m_rootfile;

  StatusCode sc;
  sc = this->GaudiAlgorithm::finalize();
  return sc;
}

---