In This Package:

QMLFTool.cc

Go to the documentation of this file.

#include "QMLFTool.h"
#include "OpPara.h"
#include "GeomPara.h"

#include "Event/RecTrigger.h"

#include "Event/CalibReadout.h"
#include "Event/CalibReadoutPmtCrate.h"

#include "CLHEP/Vector/ThreeVector.h"
#include "CLHEP/Units/SystemOfUnits.h"
#include "DetHelpers/IPmtGeomInfoSvc.h"
#include "DetHelpers/ICoordSysSvc.h"
#include "Conventions/Detectors.h"
#include "Conventions/Reconstruction.h"
#include "DataSvc/ICableSvc.h"
#include "DataSvc/ICalibDataSvc.h"

#include "TMinuit.h"
#include "TMath.h"

#include "DetDesc/Material.h"
#include "DetDesc/IGeometryInfo.h"
#include "DetDesc/Surface.h"
#include "DetDesc/TabulatedProperty.h"

using namespace Gaudi;
using namespace DayaBay;
using namespace std;

ICableSvc            *QMLFTool::m_cableSvc = 0;
IPmtGeomInfoSvc      *QMLFTool::m_pmtGeomSvc = 0;
ICalibDataSvc        *QMLFTool::m_pmtCalibDataSvc = 0;

DayaBay::CalibReadoutPmtCrate *QMLFTool::m_calibCrate = 0;

QMLFTool::ExpChgMap QMLFTool::m_expChgMap;

OpPara               *QMLFTool::m_opPara = 0;
GeomPara             *QMLFTool::m_geomPara = 0;
QMLFTool::UsedFeeChannelMap        QMLFTool::m_timeLLFUsedChannels;
QMLFTool::UsedFeeChannelMap        QMLFTool::m_chgLLFUsedChannels;

int QMLFTool::m_LLFMode = 1;
bool QMLFTool::m_addTimeLLF = 0;
double QMLFTool::m_pmtChgCut = 0.0;
bool    QMLFTool::m_recExtInfo = 0;

// 
QMLFTool::QMLFTool(const std::string& type,
      const std::string& name,
      const IInterface* parent)
  :  GaudiTool(type, name, parent)
  //, m_recTrigExt(0)
{
  declareInterface<IReconTool>(this);

  declareProperty("CableSvcName", m_cableSvcName = "StaticCableSvc",
      "Name of service to map between detector, hardware, and electronic IDs");
  declareProperty("PmtGeomSvcName", m_pmtGeomSvcName = "PmtGeomInfoSvc", 
      "Name of Pmt Geometry Information Service");
  declareProperty("PmtCalibDataSvcName", m_pmtCalibDataSvcName = "StaticCalibDataSvc", 
      "Name of Pmt Calib Data Information Service");

  // declare property for those optical parameters
  declareProperty("LLFMode", m_LLFMode = 1, "Likelihood function mode");
  declareProperty("addTimeLLF", m_addTimeLLF = 0, "Add time likelihood or not");
  declareProperty("pmtChgCut", m_pmtChgCut = 12., "The charge cut for PMTs selection criteria");
  declareProperty("recExtInfo", m_recExtInfo = 0, "Extended information from reconstruction");
  declareProperty("opLocation", m_opLocation="DetDesc", "Optical parameters location");
  declareProperty("geomLocation", m_geomLocation="DetDesc", "Optical parameters location");
  declareProperty("AbsLength", m_ATTL=0., "Effective absorption length");
  declareProperty("TopRefZ", m_TopRefZpos=0., "Z position of top reflector");
  declareProperty("BotRefZ", m_BotRefZpos=0., "Z position of bottom reflector");
  declareProperty("TopReflectivity", m_TopRef=0., "Reflectivity of top reflector");
  declareProperty("BotReflectivity", m_BotRef=0., "Reflectivity of bottom reflector");
}

QMLFTool::~QMLFTool() {}

StatusCode QMLFTool::initialize()
{
  // Get Cable Service
  m_cableSvc = svc<ICableSvc>(m_cableSvcName, true);
  // Get Pmt Geometry Service
  m_pmtGeomSvc = svc<IPmtGeomInfoSvc>(m_pmtGeomSvcName, true);

  // Get Pmt Calib Data Service
  m_pmtCalibDataSvc = svc<ICalibDataSvc>(m_pmtCalibDataSvcName, true);

  // Get Optical Parameters
  getOpPara(m_opLocation);

  // Get Geometry Parameters
  getGeomPara(m_geomLocation);

  // Show Optical and Geom Parameters
  printPara();

  return StatusCode::SUCCESS;
}

StatusCode QMLFTool::finalize()
{
    return StatusCode::SUCCESS;
}

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

  //if(m_recExtInfo) {
  //   m_recTrigExt = new DayaBay::RecTriggerExt();
  //   m_expChgMap.clear();
  //}
   DayaBay::CalibReadout* calibTest 
     = const_cast<DayaBay::CalibReadout*> (&calibReadout);
   DayaBay::CalibReadoutPmtCrate* calibCrate 
     = dynamic_cast<DayaBay::CalibReadoutPmtCrate*> (calibTest);
   m_calibCrate = calibCrate;

   if(!m_calibCrate) {
     error() << "Incorrect type of readout crate for detector "
             << calibReadout.detector().detName() << endreq;
     recTrigger.setPositionStatus(ReconStatus::kBadReadout);
     recTrigger.setEnergyStatus(ReconStatus::kBadReadout);
     return StatusCode::FAILURE;
   }

   info() << "QMLFTool: test output" << endreq;

   // determine FeeChannels that used for reconstruction
   setUsedChannels(calibCrate, &recTrigger);
   
   //
   // Create TMinuit
   const int Npar = 4;
   TMinuit *recMinuit = new TMinuit(Npar);
   int ierflg;

   recMinuit->SetFCN(QMLFTool::NCLL_fcn);
   recMinuit->SetPrintLevel(1);

   // Initial parameter   
   double xini, yini, zini, eini, tini;
   xini = recTrigger.position().x();
   yini = recTrigger.position().y();
   zini = recTrigger.position().z();
   eini = recTrigger.energy();
   tini = recTrigger.triggerTime().GetNanoSec();

   if(eini==0.0) {
     error() << "No initial energy value was set!" << endreq;
   }
   if(xini==0.0&&yini==0.0&&zini==0.0) {
     warning() << "No initial vertex was set!" << endreq;
   }
   
   // Initial TMinuit Fitting Parameter
   recMinuit->mnparm(0, "xpos",   xini,  100.0, 0.0, 0.0, ierflg);
   recMinuit->mnparm(1, "ypos",   yini,  100.0, 0.0, 0.0, ierflg);
   recMinuit->mnparm(2, "zpos",   zini,  10.0, m_geomPara->m_botRefZ, m_geomPara->m_topRefZ, ierflg);
   recMinuit->mnparm(3, "energy", eini,   0.1, 0.001, 1000.0, ierflg);
   if(m_addTimeLLF==1) 
     recMinuit->mnparm(4, "T0", tini,  0.5, 0.0, 0.0, ierflg);
   // Set parameters for Minuit minimization
   int    err;
   double standin = 0.0;
   double maxCalls = 2000;
   double tolerance = 1.0e-3;
   double arglist[2];

   // Call FCN
   arglist[0] = 1;
   recMinuit->mnexcm("CALL FCN", arglist, 1, ierflg);

   // No Warnings
   recMinuit->mnexcm("SET NOW", &standin, 0, err);
   // Set error Definition
   // 1 for Chi square
   // 0.5 for negative log likelihood
   if(m_LLFMode==1)
     recMinuit->SetErrorDef(1);
   else 
     recMinuit->SetErrorDef(0.5);
   // Minimization strategy
   // 1 standard;  2 try to improve minimum (slower)
   arglist[0] = 2;
   recMinuit->mnexcm("SET STR", arglist, 1, ierflg); 
  
   //Do Minuit fit!  
   // set max iterations: recMinuit->SetMaxIteration(1000); recMinuit->Migrad();
   arglist[0] = maxCalls;
   arglist[1] = tolerance;
   recMinuit->mnexcm("MIGrad", arglist, 2, err);
   //recMinuit->mnexcm("MINIMIZE", arglist, 2, err);
   //recMinuit->mnexcm("SIMPlex", arglist, 2, err);

   //Scan is usefull to check FCN
   // arglist[0] = 0;
   // recMinuit->mnexcm("SCAN", arglist, 1, ierflg);

   //Collect the output
   //How to draw the n-sigma contour of a Minuit fit?
   //http://root.cern.ch/cgi-bin/print_hit_bold.pl/root/html508/examples/fitcont.C.html?minuit#first_hit

   // record fitter's status
   double ln0, edm, errdef;
   int nvpar, nparx, icstat;
   recMinuit->mnstat(ln0, edm, errdef, nvpar, nparx, icstat);
   
   // Get Fitted Parameter
   double x_fit, x_fit_err;
   double y_fit, y_fit_err;
   double z_fit, z_fit_err;
   double e_fit, e_fit_err;

   recMinuit->GetParameter(0, x_fit, x_fit_err);
   recMinuit->GetParameter(1, y_fit, y_fit_err);
   recMinuit->GetParameter(2, z_fit, z_fit_err);
   recMinuit->GetParameter(3, e_fit, e_fit_err);

   // error matrix; mnemat()
   double emat[Npar][Npar];
   CLHEP::HepMatrix matrix(Npar, Npar);
   recMinuit->mnemat(&emat[0][0], Npar);
   for(int row=0; row<Npar; row++) {
     for(int col=0; col<Npar; col++) {
       matrix[row][col] = emat[row][col];
     }
   }

   CLHEP::HepLorentzVector vertex(x_fit, y_fit, z_fit, 0.);
   recTrigger.setPosition(vertex);
   recTrigger.setEnergy(e_fit);
   recTrigger.setPositionQuality(ln0);
   recTrigger.setEnergyQuality(ln0);
   recTrigger.setErrorMatrix(matrix);

   if(m_recExtInfo) {
   arglist[0] = 3;
   recMinuit->mnexcm("CALL FCN", arglist, 1, ierflg);
   //m_recTrigExt->setExpChgMap(m_expChgMap);
   //recTrigger.setRecTrigExt(m_recTrigExt);
   }
   delete recMinuit;
       
   return StatusCode::SUCCESS;
}

void QMLFTool::NCLL_fcn(int& npar, double* grad, double& fval, double* xval, int iflag) {

  if(iflag==1){
    // Initialization
    // read input data, calculate any necessary constants, etc
  }
  if(iflag==2){
    // Compute derivatives (GRAD), the first derivatives of FVAL
    // (this is optional)
  }

  // Always calculate the value of the function, FVAL,
  // which is usually a chisquare or log likelihood.
  
    double ECAL;
    // Global coefficient (QE, Yield...) 
    ECAL = m_opPara->m_photoCathodeArea/4.0/TMath::Pi()*m_opPara->m_yield;
    ECAL = ECAL * 0.20; // Nominal 'absolute' QE: 0.20
    //ECAL = 1000.;

    double S_angle, dist;
    double expected_q, qhit, cosa_tmp;
    double sum_exp, AA, BB, LLF;
    double XD, YD, ZD, ED, ZF, dx, dy, dz, T0;
    double chg_likelihood, time_likelihood;
    int ipmt, npe;

    XD = xval[0];
    YD = xval[1];
    ZD = xval[2];
    ED = xval[3]*ECAL;
    if(m_addTimeLLF==1) T0 = xval[4];

    double sum_act = 0.0;
    sum_exp = 0.0;
    chg_likelihood = 0.0;
    time_likelihood = 0.0;
    double c_n = 200.; // light speed in liquid, mm/s

    DayaBay::CalibReadoutPmtCrate::PmtChannelReadouts::const_iterator cpcrIter,
      done = (m_calibCrate->channelReadout()).end();
    
    Context context = m_calibCrate->header()->context();
    int task = 0;
    ServiceMode svcMode(context, task);

    Site::Site_t site = m_calibCrate->detector().site();
    DetectorId::DetectorId_t  detid = m_calibCrate->detector().detectorId();

    // Loop over calib readout channels
    // Remember: only non-zero readout channels are stored the map
    double sumQ = 0.0, Xc, Yc, Zc;
    for(int local_id=0; local_id<192; local_id++) {
      int ring = local_id/24 + 1;
      int col  = local_id%24 + 1;

      AdPmtSensor adPmtId(ring, col, site, detid);
      FeeChannelId channelId = m_cableSvc->feeChannelId(adPmtId, svcMode);

      // determine whether to use this channel
      if(iflag!=3&&m_chgLLFUsedChannels[channelId]==false) continue;

      DayaBay::CalibReadoutPmtChannel* cpcrPtr = m_calibCrate->sensor(adPmtId);
      if(cpcrPtr) { // map has the record of this channelId
        qhit = cpcrPtr->maxCharge();
      }else { // map has no record for this channel, 
              // yet it still has contribution to likelihood
        qhit = 0.0;
      }

      unsigned int pmtid = adPmtId.fullPackedData();
      if(adPmtId.site() == Site::kSAB) {
        pmtid = DayaBay::AdPmtSensor(adPmtId.ring(),adPmtId.column(),Site::kDayaBay,adPmtId.detectorId()).fullPackedData();
      }
      CLHEP::Hep3Vector pos =  m_pmtGeomSvc->get(pmtid)->localPosition();
      CLHEP::Hep3Vector norm = m_pmtGeomSvc->get(pmtid)->localDirection();

      // Get pmt data from ICalibDataSvc
      const DayaBay::DetectorSensor& sensDetId 
        = m_cableSvc->sensor(channelId, svcMode);
      const DayaBay::PmtCalibData* pmtData =
        m_pmtCalibDataSvc->pmtCalibData(sensDetId, svcMode);

      AA = norm.x();
      BB = norm.y();

      dx = XD - pos.x();
      dy = YD - pos.y();
      dz = ZD - pos.z();

      dist = dx*dx + dy*dy + dz*dz;
      if(dist<0.0001) dist = 0.0001;
      dist = sqrt(dist);

      cosa_tmp = ( AA*dx + BB*dy ) / dist;
      if(cosa_tmp>1.0) cosa_tmp = 1.0;

      // whether to use time information of this channel
      if(m_addTimeLLF==1) {
        if(m_timeLLFUsedChannels[channelId]==true){
          // get the earliest TDC time as the 'first hit' time, 
          // calculate the time likelihood
          double hittime = cpcrIter->earliestTime();
          time_likelihood = time_likelihood 
                  - 2*log(1./TMath::Sqrt(2*TMath::Pi())/pmtData->m_timeSpread)
                  + TMath::Power((hittime-pmtData->m_timeOffset-dist/c_n-T0)/pmtData->m_timeSpread, 2);
        }
      }

      S_angle = 0.378 + 0.5099*cosa_tmp + 0.1131*pow(cosa_tmp,2);
      if(S_angle<1.0e-5) S_angle = 1.0e-5;

      expected_q = S_angle*exp(-dist/m_opPara->m_attl)/TMath::Power(dist, 2);


      int ref_deg;
      for(ref_deg=0; ref_deg<10; ref_deg++){
     
        if(ref_deg==0)
          ZF =  2*m_geomPara->m_botRefZ - xval[2];   // 1st-order reflection
        else if(ref_deg==1)
          ZF =  2*m_geomPara->m_topRefZ - xval[2];   // 1st-order reflection
        else if(ref_deg==2)
          ZF = 2*m_geomPara->m_botRefZ - 2*m_geomPara->m_topRefZ + xval[2]; // 2nd-order reflection
        else if(ref_deg==3)
          ZF = 2*m_geomPara->m_topRefZ - 2*m_geomPara->m_botRefZ + xval[2]; // 2nd-order reflection
        else if(ref_deg==4)
          ZF = 4*m_geomPara->m_botRefZ - 2*m_geomPara->m_topRefZ - xval[2]; // 3rd-order reflection
        else if(ref_deg==5)
          ZF = 4*m_geomPara->m_topRefZ - 2*m_geomPara->m_botRefZ - xval[2]; // 3rd-order reflection
        else if(ref_deg==6)
          ZF = 4*m_geomPara->m_botRefZ - 4*m_geomPara->m_topRefZ + xval[2]; // 4th-order reflection
        else if(ref_deg==7)
          ZF = 4*m_geomPara->m_topRefZ - 4*m_geomPara->m_botRefZ + xval[2]; // 4th-order reflection
        else if(ref_deg==8)
          ZF = 6*m_geomPara->m_botRefZ - 4*m_geomPara->m_topRefZ - xval[2]; // 5th-order reflection
        else
          ZF = 6*m_geomPara->m_topRefZ - 4*m_geomPara->m_botRefZ - xval[2]; // 5th-order reflection

        dz = ZF - pos.z();

        dist = dx*dx + dy*dy + dz*dz;
        if(dist<0.0001) dist = 0.0001;
        dist = sqrt(dist);

        cosa_tmp = ( AA*dx + BB*dy ) / dist;
        if(cosa_tmp>1.0) cosa_tmp = 1.0;

        S_angle = 0.378 + 0.5099*cosa_tmp + 0.1131*pow(cosa_tmp,2);
        if(S_angle<1.0e-5) S_angle = 1.0e-5;

        int index = (int)(ref_deg/2.0) + 1;
        int index_top, index_bot;
        if(ref_deg == 0) 
          index_top = 0;
        else 
          index_top = (int)((ref_deg-1)/4.0) + 1;

        index_bot = index - index_top;

        expected_q += TMath::Power(m_opPara->m_topRef, index_top)*TMath::Power(m_opPara->m_botRef, index_bot)*S_angle*TMath::Exp(-dist/m_opPara->m_attl)/TMath::Power(dist, 2);
      }


      double eff = pmtData->m_efficiency;

      expected_q = ED * expected_q * eff;
      sum_exp = sum_exp + expected_q;
      sum_act = sum_act + qhit;

      if(iflag==3) {
        //after the fit is finished, get Exp_Q for each PMT
        if(m_recExtInfo) m_expChgMap[channelId] = expected_q;

      }
      else {
        if(m_LLFMode!=1){           


        }
        else
        {
          if(qhit!=0.0) 
            chg_likelihood = chg_likelihood + 2*(expected_q-qhit) + 2*log(qhit/expected_q)*qhit;
          else 
            chg_likelihood = chg_likelihood + 2*expected_q;
        }
      }
    
    } // for(local_id....

    fval = chg_likelihood;
    if(m_addTimeLLF==1) fval += time_likelihood;
}

StatusCode QMLFTool::getOpPara(std::string &data_src) {

  // construct OpPara, default value
  m_opPara = new OpPara(1.25E+04, 0.98, 0.98, 9000, 3.14*103*103);

  if(data_src=="DetDesc") {

    // Get Optical Parameters from Det/DetDesc
    // A tabular was defined in DetDesc, but only one value was get here
    // FIXME
    // Get the Detector Data Service to access TDS:
    // This should be done in your tools initialize()
    // and save m_DDS as a data member
    IDataProviderSvc* dds = 0;
    StatusCode sc = service( "DetectorDataSvc", dds, true );
    if (sc.isFailure()) {
        // ... handle error ...
    }
    // Later, look up some material.  If this material is always used you
    // should look it up once in initialize() and save it for later use

    // absorption length
    std::string location = "/dd/Materials/GdDopedLS";
    Material* gdls = GaudiCommon<AlgTool>::get<Material>(dds,location);
    
    Material::Tables& gdls_tab = gdls->tabulatedProperties();
    Material::Tables::const_iterator mtIter = gdls_tab.begin(), 
      done = gdls_tab.end();
    TabulatedProperty::Table attl_vec;
    for(; mtIter!=done; mtIter++) {
      std::string type = (*mtIter)->type();
      if(type=="ABSLENGTH") {
        attl_vec = (*mtIter)->table(); 
      }
    }
    TabulatedProperty::Table::const_iterator tb_iter;
    int count = 0;
    double attl_tmp = 0.0, ref_tmp = 0.0;
    for(tb_iter=attl_vec.begin(); tb_iter!=attl_vec.end(); tb_iter++){
      if(tb_iter->first<3.18e-06 && tb_iter->first>2.07e-06 ){ // 390nm ~ 600nm
        count++; attl_tmp += tb_iter->second;
      }
    }
    attl_tmp = attl_tmp/count; // average value between 390nm~600nm
    m_opPara->m_attl = attl_tmp;
      
    // ESR surface location
    std::string topESR_location = "/dd/Geometry/AdDetails/AdSurfacesAll/ESRAirSurfaceTop";
    std::string botESR_location = "/dd/Geometry/AdDetails/AdSurfacesAll/ESRAirSurfaceBot";

    Surface* esrtop = GaudiCommon<AlgTool>::get<Surface>(dds, topESR_location);
    Surface* esrbot = GaudiCommon<AlgTool>::get<Surface>(dds, botESR_location);
    
    Surface::Tables& esrtop_tab = esrtop->tabulatedProperties();
    Surface::Tables& esrbot_tab = esrbot->tabulatedProperties();

    Surface::Tables::const_iterator sfIter; 
    TabulatedProperty::Table ref_vec;

    // look up Tables, get an average top ESR reflectivity
    for(sfIter=esrtop_tab.begin(); sfIter!=esrtop_tab.end(); sfIter++) {
      std::string type = (*sfIter)->type();
      if(type=="REFLECTIVITY") {
        ref_vec = (*sfIter)->table(); 
      }
    }
    count = 0;
    for(tb_iter=ref_vec.begin(); tb_iter!=ref_vec.end(); tb_iter++){
      if(tb_iter->first<3.18e-06 && tb_iter->first>2.07e-06 ){ // 390nm ~ 600nm
        count++; ref_tmp += tb_iter->second;
      }
    }
    ref_tmp = ref_tmp/count; // average value between 390nm~600nm
    m_opPara->m_topRef = ref_tmp;
 
    // look up Tables, get an average bot ESR reflectivity
    for(sfIter=esrbot_tab.begin(); sfIter!=esrbot_tab.end(); sfIter++) {
      std::string type = (*sfIter)->type();
      if(type=="REFLECTIVITY") {
        ref_vec = (*sfIter)->table(); 
      }
    }
    count = 0; ref_tmp = 0.0;
    for(tb_iter=ref_vec.begin(); tb_iter!=ref_vec.end(); tb_iter++){
      if(tb_iter->first<3.18e-06 && tb_iter->first>2.07e-06 ){ // 390nm ~ 600nm
        count++; ref_tmp += tb_iter->second;
      }
    }
    ref_tmp = ref_tmp/count; // average value between 390nm~600nm
    m_opPara->m_botRef = ref_tmp;

  } else if(data_src=="OpCalibSvc"){

  } else {

    warning() <<"No external optical data source, use default value."<< endreq;
  }

  // if set new values from declareProperty
  if(m_ATTL!=0.) m_opPara->m_attl = m_ATTL;
  if(m_TopRef!=0.) m_opPara->m_topRef = m_TopRef;
  if(m_BotRef!=0.) m_opPara->m_botRef = m_BotRef;

  return StatusCode::SUCCESS;
}

StatusCode QMLFTool::getGeomPara(std::string &data_src) {

  // construct GeomPara, default value
  m_geomPara = new GeomPara(2010., -2010.);

  if(data_src=="DetDesc") {

    IDataProviderSvc* dds = 0;
    StatusCode sc = service( "DetectorDataSvc", dds, true );
    if (sc.isFailure()) {
       error() << "No DetectorDataSvc available!" << endreq;
    }
    
    // ESR structure location
    std::string topEsrPath = "/dd/Structure/AdReflectorStructure/db-ad/db-ad1-topesr";
    std::string botEsrPath = "/dd/Structure/AdReflectorStructure/db-ad/db-ad1-botesr";
    SmartDataPtr<IDetectorElement> topEsr(dds,topEsrPath);
    SmartDataPtr<IDetectorElement> botEsr(dds,botEsrPath);

    // define a pointer to ICoordSysSvc to
    // use to cache the service during the job.
    ICoordSysSvc* css;
    sc = service( "CoordSysSvc", css, true );

    Gaudi::XYZPoint zero(0,0,0);
    // Set top ESR position
    Gaudi::XYZPoint globalPoint = topEsr->geometry()->toGlobal(zero);
    IDetectorElement* theAdDe = css->coordSysDE(globalPoint);
    Gaudi::XYZPoint adLocalPoint = theAdDe->geometry()->toLocal(globalPoint);
    m_geomPara->m_topRefZ = adLocalPoint.z();

    // Set bot ESR position
    globalPoint = botEsr->geometry()->toGlobal(zero);
    theAdDe = css->coordSysDE(globalPoint);
    adLocalPoint = theAdDe->geometry()->toLocal(globalPoint);
    m_geomPara->m_botRefZ = adLocalPoint.z();
        
  } else if(data_src=="OnsiteSuveyData"){
    // Get geometry info from onsite suvey data

  } else {

    warning() <<"Unknown geometry data source, use default value."<< endreq;
  }

  // if set new values from declareProperty
  if(m_TopRefZpos!=0.) m_geomPara->m_topRefZ = m_TopRefZpos;
  if(m_BotRefZpos!=0.) m_geomPara->m_botRefZ = m_BotRefZpos;

  return StatusCode::SUCCESS;
}

StatusCode QMLFTool::printPara(){
  info() << "Op and Geom parameters used in reconstruction." << endreq;
  info() << "Attenuation length : " << m_opPara->m_attl <<"mm"<< endreq;
  info() << "Top reflectivity : " << m_opPara->m_topRef << endreq;
  info() << "Bottom reflectivity : " << m_opPara->m_botRef << endreq;
  info() << "Top reflector Z pos :" << m_geomPara->m_topRefZ <<"mm"<< endreq;
  info() << "Bot reflector Z pos :" << m_geomPara->m_botRefZ <<"mm"<< endreq;

  return StatusCode::SUCCESS;
}

StatusCode QMLFTool::setUsedChannels(
                               const DayaBay::CalibReadoutPmtCrate* calibCrate,
                               DayaBay::RecTrigger* recTrigger)
{
  m_chgLLFUsedChannels.clear();
  m_timeLLFUsedChannels.clear();

  DayaBay::CalibReadoutPmtCrate::PmtChannelReadouts::const_iterator cpcrIter,
    done = (calibCrate->channelReadout()).end();
    
  Context context = calibCrate->header()->context();
  int task = 0;
  ServiceMode svcMode(context, task);

  Site::Site_t site = calibCrate->detector().site();
  DetectorId::DetectorId_t  detid = calibCrate->detector().detectorId();

  // default is use pmts whose status are Good
  for(int local_id=0; local_id<192; local_id++) {
    int ring = local_id/24 + 1;
    int col  = local_id%24 + 1;

    AdPmtSensor adPmtId(ring, col, site, detid);
    FeeChannelId channelId = m_cableSvc->feeChannelId(adPmtId, svcMode);
    
    const DayaBay::DetectorSensor& sensDetId = m_cableSvc->sensor(channelId, svcMode);
    const DayaBay::PmtCalibData* pmtData =
      m_pmtCalibDataSvc->pmtCalibData(sensDetId, svcMode);
    
    if(pmtData->m_status==DayaBay::PmtCalibData::kGood) { 
      m_chgLLFUsedChannels[channelId] = true;
      m_timeLLFUsedChannels[channelId] = false;
    }
    else {
      m_chgLLFUsedChannels[channelId] = false;
      m_timeLLFUsedChannels[channelId] = false;
    }

    if(m_addTimeLLF==1) { 
      DayaBay::CalibReadoutPmtChannel* cpcrPtr = m_calibCrate->sensor(adPmtId);
      if(cpcrPtr) { // map has the record of this channelId
        double qhit = cpcrPtr->maxCharge();
        if(qhit>m_pmtChgCut)
          m_timeLLFUsedChannels[channelId] = true;
      }
    }
    
  }
  // Other options, depending on CalibPmtCrate data
  //recTrigger->setUsedChannels(m_chgLLFUsedChannels);

  return StatusCode::SUCCESS;
}

---