QsumEnergyTool.cc

Go to the documentation of this file.

#include "QsumEnergyTool.h"

#include "Event/CalibReadout.h"
#include "Event/CalibReadoutPmtCrate.h"
#include "Event/CalibReadoutPmtChannel.h"
#include "Event/RecTrigger.h"
#include "Event/UserDataHeader.h"

#include "Conventions/Reconstruction.h"
#include "Conventions/Detectors.h"
#include <iostream>
#include "DataSvc/ICableSvc.h"
#include "DataSvc/ICalibDataSvc.h"
#include "DybKernel/IDetCalibSvc.h"

#include "CLHEP/Units/SystemOfUnits.h"
#include "CLHEP/Vector/ThreeVector.h"

#define NODB_ERRORS 15

using namespace CLHEP;
using namespace DayaBay;
using namespace std;

QsumEnergyTool::QsumEnergyTool(const std::string& type,
                               const std::string& name, 
                               const IInterface* parent)
  : GaudiTool(type,name,parent)
  , m_cableSvc(0)
{

  declareInterface< IReconTool >(this) ;   
    declareProperty("CableSvcName",m_cableSvcName="CableSvc",
                    "Name of service to map between detector, hardware, and electronic IDs");
    declareProperty("RecDataSvcName", m_recDataSvcName="DybDetCalibSvc",
                    "Name of calibration data service");
    declareProperty("TaskForEnergyScale",m_task_es=1,
                    "Task number for DB energy scale lookup (0 for Co60; 1 for SN)");
    declareProperty("CalibStatsLocation",m_calibStatsLocation="/Event/Data/CalibStats",
                    "Location of the CalibStats");
    declareProperty("Charge",m_chargeLeaveName="NominalCharge",
                    "Name of the leave that hold the total charge in CalibStats");

    std::string Path = getenv ("DATASVCROOT");
    declareProperty("UniformityCorrFileName", m_uniformityCorrFileName = Path+"/share/uni-corr-p14a-v1.root",
                    "File name of the charge distribution template");

    declareProperty("UniformityCorrFileNameMC", m_uniformityCorrFileNameMC = Path+"/share/uni-corr-m14a-v1.root",
                    "File name of the charge distribution template for MC");

    // This root file include the new AD-by-AD non-uniformity correction.
    // Currently, this new non-uniformity correction works only for data.
    
}

QsumEnergyTool::~QsumEnergyTool(){}

StatusCode QsumEnergyTool::initialize()
{
  // Get Cable Service
  m_cableSvc = svc<ICableSvc>(m_cableSvcName,true);
  m_recDataSvc = svc<IDetCalibSvc>(m_recDataSvcName,true);

  // Initialize error counter
  nodb_count=0;
  
  info() << "Setting up non-uniformity correction using " << m_uniformityCorrFileName << endreq;
  m_corrector = new Corrector(m_uniformityCorrFileName.c_str());
  
  info() << "Setting up non-uniformity correction for MC using " << m_uniformityCorrFileNameMC << endreq;
  m_corrector_mc = new Corrector(m_uniformityCorrFileNameMC.c_str());
  
  return StatusCode::SUCCESS;

}

StatusCode QsumEnergyTool::finalize()
{

  if(nodb_count>0){
    warning() << "There were " << nodb_count << " warnings about not being able to find the energy scale constant in the database" << endreq;
  }

  return StatusCode::SUCCESS;
}

StatusCode QsumEnergyTool::reconstruct(const CalibReadout& readout,
                                       RecTrigger& recTrigger)
{
  if( !readout.detector().isAD() ){
    debug() << "Not an AD readout; ignoring detector " 
            << readout.detector().detName() << endreq;
    recTrigger.setEnergyStatus( ReconStatus::kNotProcessed );
    return StatusCode::SUCCESS;
  }

  const CalibReadoutPmtCrate* pmtReadout 
    = dynamic_cast<const CalibReadoutPmtCrate*>(&readout);
  if(!pmtReadout){
    error() << "Incorrect type of readout crate for detector "
            << readout.detector().detName() << endreq;
    recTrigger.setEnergyStatus( ReconStatus::kBadReadout );
    return StatusCode::FAILURE;
  }


  // Context for this data
  // Note: when looking up pmt status (right below) need task of zero. When looking up energy scale the task may change, according to the value of m_task_es (jpochoa)
  ServiceMode svcMode(readout.header()->context(), 0);
  ServiceMode svcMode_es(readout.header()->context(), m_task_es);


  // Loop over hits and add up charge
  double qSum = 0;
  double qSquareSum = 0;
  int nHits = 0;

  
  DayaBay::UserDataHeader* calibStatsHdr = 0;
  if(exist<UserDataHeader>(m_calibStatsLocation)){
    calibStatsHdr = get<UserDataHeader>(m_calibStatsLocation);
  }else{
    error() << "Cannot get CalibStats information" << endreq;
    error() << "QsumEnergyTool now require to run CalibStats before this!!!!" << endreq;
    return StatusCode::FAILURE;

    // followings are obsolete code
    
    CalibReadoutPmtCrate::PmtChannelReadouts::const_iterator chanIter, 
      chanEnd = pmtReadout->channelReadout().end();
    for(chanIter=pmtReadout->channelReadout().begin(); chanIter != chanEnd; 
        chanIter++){
      const CalibReadoutPmtChannel& channel = *chanIter;
      AdPmtSensor pmtId(channel.pmtSensorId().fullPackedData());
      //AdPmtSensor pmtId = m_cableSvc->adPmtSensor(channel.channelId(), svcMode);
      if( pmtId.bogus() || pmtId.ring() < 1 ){ continue; }  // Ignore 2" pmts and non-pmt channels  (timing reference, noise... etc)
         
      //if(pmtId.ring()==1&&pmtId.column()==1) cout<<readout.detector().detName()<<' '<<pmtCalib->m_speHigh<<endl;
      //    if( pmtCalib->m_status != PmtCalibData::kGood ) { continue; } // this is obsolete
      
      double peakAdc = channel.earliestCharge(-1650,-1250);
      qSum += peakAdc;
      qSquareSum += peakAdc*peakAdc;
      if(peakAdc>0){
        nHits++;
      }
    }

  }

  
  qSum = calibStatsHdr->getFloat(m_chargeLeaveName);
  nHits = (int)calibStatsHdr->getInt("nHit");
  float nActivePMTs = calibStatsHdr->getFloat("nActivePMTs");
  
  
  if(nHits<1){
    recTrigger.setEnergyStatus( ReconStatus::kNoHits );
    return StatusCode::SUCCESS;
  }

  if(recTrigger.positionStatus()==ReconStatus::kNotConverged){
    return StatusCode::SUCCESS;
  }

  
  //Get position-dependent correction
  double posDepCorr=1.;
  if(recTrigger.positionStatus()==ReconStatus::kGood){
    const CLHEP::Hep3Vector pos = recTrigger.position();
    double position_unit = CLHEP::meter;
    double posx = pos.x()/position_unit;
    double posy = pos.y()/position_unit;
    double posz = pos.z()/position_unit;
    if ((readout.header()->context()).GetSimFlag()== SimFlag::kMC){
      int site = (readout.header()->context()).GetSite();
      int detid = (readout.header()->context()).GetDetId();
      posDepCorr = m_corrector_mc->get_corr(site,detid,posx,posy,posz);
    }else{
      // posDepCorr = uniformCorr(posx,posy,posz);//<-- uniformCorr takes reco pos in meters
      int site = (readout.header()->context()).GetSite();
      int detid = (readout.header()->context()).GetDetId();
      posDepCorr = m_corrector->get_corr(site,detid,posx,posy,posz);

    }
  }  

  // Perform scaling 
  //-->Get number from database
  double scaleConstant=170.; //default value to use for now (PEtoEvis)
  // Note: change task number when looking-up energy scale number 
  const DayaBay::DetEnergyReconData* srcEn = m_recDataSvc->detEnergyReconData(svcMode_es);
  if(!srcEn){
    nodb_count+=1;
    if(nodb_count < NODB_ERRORS){
      warning() << "No energy scale found in DB for context " << readout.header()->context().AsString() << " and task " << m_task_es << ". Will use a default value of " << scaleConstant << endreq;
    }
    if(nodb_count == NODB_ERRORS){
      warning() << "No more errors about not finding energy scale in DB will be reported. A summary count of these errors will appear at the end." << endreq;
    }
  } else {
    scaleConstant = srcEn->m_PeEvis;
  }

  // Apply scale constant based on the number of live PMTs. Aug. 2012, Yasu.
  if (nActivePMTs > 0)
    scaleConstant = scaleConstant * nActivePMTs/192.;
  
  //-->do the math
  double energy = 0;
  double rawEvis = 0;
  if(scaleConstant>0) {
    rawEvis = (qSum *1./ scaleConstant);
    energy = rawEvis*posDepCorr;
  }else {
    error() << "scaleConstant found to be negative or zero! Something is seriously wrong, so will abort" << endreq;
    return StatusCode::FAILURE;
  }
  //--> storage
  recTrigger.setRawEvis( rawEvis );
  recTrigger.setEnergy( energy );  
  // (set quality to standard deviation of hits)
  recTrigger.setEnergyQuality(sqrt((qSquareSum - qSum*qSum/nHits)/(nHits-1)));
  recTrigger.setEnergyStatus( ReconStatus::kGood );

  return StatusCode::SUCCESS;
}


// Below are correction used for P13A, M13B and before, and no longer used.
// P14A, M14A and later uses new correction function "Corrector()"



// // Position-dependent correction 
// // Important: reconstructed position must be provided in meters 
// double QsumEnergyTool::uniformCorr(double reco_x, double reco_y, double reco_z){
  

//   double reco_r = sqrt(pow(reco_x,2)+pow(reco_y,2));

//   // Set maximum r and z position to avoid over-correcting
  
//   double max_reco_r = 2.5;
//   double max_reco_z = 2.5;
//   double min_reco_z = -2.5;
  
//   if (reco_r > max_reco_r) reco_r = max_reco_r;
//   if (reco_z > max_reco_z) reco_z = max_reco_z;
//   if (reco_z < min_reco_z) reco_z = min_reco_z;
  
//   double corr = 8.05 / (7.84628 * (1 + 3.41294e-02*reco_r*reco_r) * (1 - 1.21750e-02*reco_z  - 1.64275e-02*reco_z*reco_z + 7.33006e-04*pow(reco_z,3)));

  
//   return corr;
// }

// double QsumEnergyTool::uniformCorrMC(double reco_x, double reco_y, double reco_z){
  
//   double reco_r = sqrt(pow(reco_x,2)+pow(reco_y,2));

//   // Set maximum r and z position to avoid over-correcting
  
//   double max_reco_r = 2.5;
//   double max_reco_z = 2.5;
//   double min_reco_z = -2.5;
  
//   if (reco_r > max_reco_r) reco_r = max_reco_r;
//   if (reco_z > max_reco_z) reco_z = max_reco_z;
//   if (reco_z < min_reco_z) reco_z = min_reco_z;
  
//   double corr = 7.937 / (7.66412 * (1 + 0.0399542*pow(reco_r,2) + 0.00116836*pow(reco_r,4)) * (1 - 0.00994061*reco_z -  0.0100378*reco_z*reco_z + 0.000827075*pow(reco_z,3)));

//   return corr;
// }


Corrector::Corrector(const char* fname){
  infile = new TFile(fname);
  char hname[100];
  for (int detidx=0; detidx<8; ++detidx){
    sprintf(hname, "resp_r2z_%d", detidx+1);
    resp_r2z[detidx] = (TH2D*)infile->Get(hname);

    sprintf(hname, "resp_phi_%d", detidx+1);
    resp_phi[detidx] = (TF1*)infile->Get(hname);
  }
}

Corrector::~Corrector(){
  if (infile != 0 && infile->IsOpen())
    infile->Close();
  if (infile != 0)
    delete infile;
}

int Corrector::get_ad(int site, int detid){
  switch(site) {
  case 1:
    return detid;
    break;
  case 2:
    if (detid == 1)
      return 3;
    else if (detid == 2)
      return 8;
    break;
  case 4:
    return detid + 3;
    break;
  default:
    return 0;
  }
  return 0;
}

// site: 1,2,4 for EH1,2,3 respectively
// detid: 1,2,3,4
// x, y, z vertex in meter
double Corrector::get_corr(int site, int detid,
                           double x, double y, double z) {
  int detidx = get_ad(site, detid) - 1;
  double r2 = x*x + y*y;
  double phi = atan2(y, x);

  if (abs(z) >= 2)
    z = z > 0 ? 1.999 : -1.999;
  if (r2 >= 4)
    r2 = 3.999;

  return 1. / resp_r2z[detidx]->Interpolate(r2, z)
    / resp_phi[detidx]->Eval(phi);
}