/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

Public Member Functions | Static Public Member Functions | Public Attributes | Static Public Attributes | Private Types | Private Member Functions
LikelihoodTool Class Reference

#include <LikelihoodTool.h>

Inheritance diagram for LikelihoodTool:
Inheritance graph
[legend]
Collaboration diagram for LikelihoodTool:
Collaboration graph
[legend]

List of all members.

Public Member Functions

 LikelihoodTool (const std::string &type, const std::string &name, const IInterface *parent)
virtual ~LikelihoodTool ()
virtual StatusCode initialize ()
virtual StatusCode finalize ()
virtual double calculate (const DayaBay::CalibReadout &, const DayaBay::RecTrigger &)
 This is the extension. Sub classes must provide it.

Static Public Member Functions

static const InterfaceID & interfaceID ()
 Retrieve interface ID.

Public Attributes

std::string m_cableSvcName
std::string m_pmtGeomSvcName
std::string m_pmtCalibDataSvcName
std::string m_src
double m_ATTL
double m_TopRef
double m_BotRef
double m_LightYield
double m_PhotoCathodeArea
double recx
double recy
double recz
double rece
double charge [192]
double pmtX [200]
double pmtY [200]
double pmtZ [200]
double pmtNormX [200]
double pmtNormY [200]
double pmtNormZ [200]
double m_TopRefZpos
double m_BotRefZpos
float m_timeCutLow
float m_timeCutHigh

Static Public Attributes

static ICableSvcm_cableSvc
static IPmtGeomInfoSvcm_pmtGeomSvc
static ICalibDataSvcm_pmtCalibDataSvc
static
DayaBay::CalibReadoutPmtCrate
m_calibCrate = 0
static int m_LLFMode = 1
static bool m_addTimeLLF = 0
static double m_pmtChgCut = 0.0
static TGraph * m_gfangle = 0
static TGraph * m_gfdistance = 0
static std::vector< float > m_obsq
static std::vector< float > m_firstT

Private Types

typedef std::map
< DayaBay::DetectorSensor,
bool > 
UsedPMTsMap

Private Member Functions

StatusCode setUsedPMTs (DayaBay::CalibReadoutPmtCrate *)

Detailed Description

Definition at line 34 of file LikelihoodTool.h.


Member Typedef Documentation

typedef std::map<DayaBay::DetectorSensor, bool> LikelihoodTool::UsedPMTsMap [private]

Definition at line 65 of file LikelihoodTool.h.


Constructor & Destructor Documentation

LikelihoodTool::LikelihoodTool ( const std::string &  type,
const std::string &  name,
const IInterface *  parent 
)

Definition at line 55 of file LikelihoodTool.cc.

:  GaudiTool(type, name, parent)
{
  declareInterface< ILikelihoodTool >(this);
  //declareInterface<IReconTool>(this);

  declareProperty("CableSvcName", m_cableSvcName = "CableSvc",
      "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 = "CalibDataSvc", 
  //      "Name of Pmt Calib Data Information Service");

  // declare property for those optical parameters
  declareProperty("LLFMode", m_LLFMode = 1, "LikelihoodToolhood 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("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");
  declareProperty("m_timeCutLow", m_timeCutLow=-1700., "Lower limit of the TDC cut to get rid of noise");
  declareProperty("m_timeCutHigh", m_timeCutHigh=-1400., "Higher limit of the TDC cut to get rid of noise");
}
LikelihoodTool::~LikelihoodTool ( ) [virtual]

Definition at line 85 of file LikelihoodTool.cc.

{}

Member Function Documentation

StatusCode LikelihoodTool::initialize ( ) [virtual]

Definition at line 87 of file LikelihoodTool.cc.

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

  //  m_pmtGeomSvc = svc<IPmtGeomInfoSvc>(m_pmtGeomSvcName, true);

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

  //  std::string fname = getenv("LIKELYHOODALGROOT");
  //std::string fname = getenv("LIKELYHOODROOT");
  std::string fname = getenv("LIKELIHOODROOT");
  fname += "/share/PMTSet.txt";
  ifstream *inFile;
  inFile = new ifstream(fname.c_str());

  //ifstream *inFile;
  //inFile = new ifstream("/dybfs/rec/Temp/xiadm/PMTSet.txt");
  //  inFile = new ifstream("/publicfs/dyb/data/userdata/xiadm/xiadm/PMTSet.txt");
//  if(!inFile){
 //   cout << "Open the file: " << "PMTGeom.txt"<< " error!" << endl;
 //   exit(1);
 // } else {
  //  cout << "Successfully open the file " << "PMTGeom.txt" << "!" << endl;
 // }

  double pmtX[200], pmtY[200], pmtZ[200], pmtNormX[200], pmtNormY[200], pmtNormZ[200];
  int i=0;
 // while(1)
 // {
    //(*inFile)>>pmtX[i]>>pmtY[i]>>pmtZ[i]>>pmtNormX[i]>>pmtNormY[i]>>pmtNormZ[i];
    //if((*inFile).eof())
 //   inFile>>pmtX[i]>>pmtY[i]>>pmtZ[i]>>pmtNormX[i]>>pmtNormY[i]>>pmtNormZ[i];
   // if(inFile.eof())
    //  break;
    //i++;
  //}


  for(int i=0;i<192;i++)
  { 
    (*inFile)>>pmtX[i]>>pmtY[i]>>pmtZ[i]>>pmtNormX[i]>>pmtNormY[i]>>pmtNormZ[i];
    cout<<"PMTSet: "<<pmtX[i]<<" "<<pmtY[i]<<" "<<pmtZ[i]<<" "<<pmtNormX[i]<<" "<<pmtNormY[i]<<" "<<pmtNormZ[i]<<endl;
  }


  // Get Optical Parameters
  // getOpPara(m_opLocation);
  //createGraph();

  // Get Geometry Parameters
  //  getGeomPara(m_geomLocation);

  // Show Optical and Geom Parameters
  //printPara();

  return StatusCode::SUCCESS;
}
StatusCode LikelihoodTool::finalize ( ) [virtual]

Definition at line 150 of file LikelihoodTool.cc.

{
  return StatusCode::SUCCESS;
}
double LikelihoodTool::calculate ( const DayaBay::CalibReadout ,
const DayaBay::RecTrigger  
) [virtual]

This is the extension. Sub classes must provide it.

Implements ILikelihoodTool.

Definition at line 156 of file LikelihoodTool.cc.

{
  debug() << "running likelihood calculation " << endreq;

  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;
  }

  DayaBay::CalibReadout* calibTest 
    = const_cast<DayaBay::CalibReadout*> (&calibReadout);
  DayaBay::CalibReadoutPmtCrate* calibCrate 
    = dynamic_cast<DayaBay::CalibReadoutPmtCrate*> (calibTest);
  m_calibCrate = calibCrate;

  DetectorId::DetectorId_t  Detid = m_calibCrate->detector().detectorId();


  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() << "LikelihoodTool: test output" << endreq;

  // determine PMTs that used for reconstruction
  setUsedPMTs(calibCrate);//, &recTrigger);
  //const DayaBay::RecHeader* srcRecHeader = get<DayaBay::RecHeader>(m_src);
  //const DayaBay::RecTrigger& srcTrigger = srcRecHeader->recTrigger();
  const CLHEP::HepLorentzVector& srcPos = srcTrigger.position();
  //


  // pos.x();
  // Initial parameter   
  //double recx, recy, recz, rece;
  recx = srcPos.x();
  recy = srcPos.y();
  recz = srcPos.z();
  rece = srcTrigger.energy();


  double ECAL;
  // Global coefficient (QE, Yield...) 
  double m_photoCathodeArea=3.14*103*103;
  double m_yield=9000;
  double m_attl=15400;
  float m_botRefZ=-2027.5;
  float m_topRefZ=2104.55;
  double m_botRef=0.94;
  double m_topRef=0.94;
  ECAL = m_photoCathodeArea/4.0/TMath::Pi()*m_yield;
  ECAL = ECAL * 0.19; // Nominal 'absolute' QE: 0.20
  //ECAL = 1000.;

  double S_angle, dist;
  double expected_q, cosa_tmp, Theta;
  double AA, BB;
  double XD, YD, ZD, ED, ZF, dx, dy, dz, T0;
  double chg_likelihood, time_likelihood;

  XD = recx;
  YD = recy;
  ZD = recz;
  ED = rece*ECAL;

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

  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
  for(int local_id=0; local_id<192; local_id++) {

    //   unsigned int pmtid = m_chgLLFUsedPMTs[local_id];
    //  if ( pmtid == 0 ) continue;
    // IPmtGeomInfo* geoInfo = m_pmtGeomSvc->get(pmtid);
    //  const CLHEP::Hep3Vector& pos  = geoInfo->localPosition();
    //  const CLHEP::Hep3Vector& norm = geoInfo->localDirection();

    int ring = local_id/24 + 1;
    int col  = local_id%24 + 1;

    // determine whether to use this channel
    // m_chgLLFUsedPMTs[local_id] = true;
    // if(m_chgLLFUsedPMTs[local_id]==false) continue;

    AdPmtSensor adPmtId(ring, col, site, detid);
    // unsigned int pmtid = adPmtId.fullPackedData();

    // CLHEP::Hep3Vector pos =  m_pmtGeomSvc->get(pmtid)->localPosition();
    //  CLHEP::Hep3Vector norm = m_pmtGeomSvc->get(pmtid)->localDirection();

    //  const DayaBay::PmtCalibData* pmtData =
    //   m_pmtCalibDataSvc->pmtCalibData(adPmtId, svcMode);

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

    //  dx = XD - pos.x();
    //  dy = YD - pos.y();
    //  dz = ZD - pos.z();
    dx = XD - pmtX[local_id];
    dy = YD - pmtY[local_id];
    dz = ZD - pmtZ[local_id];

    //cout<<"set: "<<XD<<" "<<YD<<" "<<ZD<<" "<<dx<<" "<<dy<<" "<<dz<<endl;
    dist = dx*dx + dy*dy + dz*dz;
    if(dist<0.0001) dist = 0.0001;
    dist = sqrt(dist);

    //  cosa_tmp = ( AA*dx + BB*dy ) / dist;
    cosa_tmp = ( pmtNormX[local_id]*dx + pmtNormY[local_id]*dy ) / dist;
    if(cosa_tmp>1.0) cosa_tmp = 1.0;

    float test1;
    Theta = TMath::ACos(cosa_tmp);
    //  cout<<"Theta: "<<Theta<<endl;
    //    S_angle = m_gfangle->Eval(Theta,0,"");
    S_angle = 1.0164 - Theta*(0.1553 - Theta*(0.2859 - Theta*(0.7805 - 0.2926*Theta)));  
    //   if (Theta >1.) cout<<"S_angle: "<<S_angle<<endl;
    if(S_angle<1.0e-5) S_angle = 1.0e-5;

    expected_q = S_angle*TMath::Exp(-dist/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_botRefZ - recz;   // 1st-order reflection
      else if(ref_deg==1)
        ZF =  2* m_topRefZ - recz;   // 1st-order reflection
      else if(ref_deg==2)
        ZF = 2* m_botRefZ - 2* m_topRefZ + recz; // 2nd-order reflection
      else if(ref_deg==3)
        ZF = 2* m_topRefZ - 2* m_botRefZ + recz; // 2nd-order reflection
      else if(ref_deg==4)
        ZF = 4* m_botRefZ - 2* m_topRefZ - recz; // 3rd-order reflection
      else if(ref_deg==5)
        ZF = 4* m_topRefZ - 2* m_botRefZ - recz; // 3rd-order reflection
      else if(ref_deg==6)
        ZF = 4* m_botRefZ - 4* m_topRefZ + recz; // 4th-order reflection
      else if(ref_deg==7)
        ZF = 4* m_topRefZ - 4* m_botRefZ + recz; // 4th-order reflection
      else if(ref_deg==8)
        ZF = 6* m_botRefZ - 4* m_topRefZ - recz; // 5th-order reflection
      else
        ZF = 6* m_topRefZ - 4* m_botRefZ - recz; // 5th-order reflection

      //dz = ZF - pos.z();
      dz = ZF - pmtZ[local_id];

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

      //    cosa_tmp = ( AA*dx + BB*dy ) / dist;
      cosa_tmp = ( pmtNormX[local_id]*dx + pmtNormY[local_id]*dy ) / dist;
      if(cosa_tmp>1.0) cosa_tmp = 1.0;
      Theta = TMath::ACos(cosa_tmp);
      // S_angle = m_gfangle->Eval(Theta,0,"");

      S_angle = 1.0164 - Theta*(0.1553 - Theta*(0.2859 - Theta*(0.7805 - 0.2926*Theta)));
      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_topRef, index_top)*TMath::Power(m_botRef, index_bot)*S_angle*TMath::Exp(-dist/m_attl)/TMath::Power(dist, 2);
    }

    expected_q = ED * expected_q ;

    if(expected_q>0.)
    {
      if(m_obsq[local_id]>0.0) 
        chg_likelihood = chg_likelihood + 2*(expected_q-m_obsq[local_id]) 
          + 2*log(m_obsq[local_id]/expected_q)*m_obsq[local_id];
      else 
        chg_likelihood = chg_likelihood + 2*expected_q;
    }
    else {
      chg_likelihood = chg_likelihood + 2*m_obsq[local_id];
    }

  } // for(local_id....

  double fval;
  fval = chg_likelihood;
  return fval;
}
StatusCode LikelihoodTool::setUsedPMTs ( DayaBay::CalibReadoutPmtCrate calibCrate) [private]

Definition at line 372 of file LikelihoodTool.cc.

{
  m_obsq.clear();

  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
  int nhit=0;
  for(int local_id=0; local_id<192; local_id++) {
    int ring = local_id/24 + 1;
    int col  = local_id%24 + 1;
    nhit++;
    AdPmtSensor adPmtId(ring, col, site, detid);
    DayaBay::CalibReadoutPmtChannel* calibChannel = m_calibCrate->sensor(adPmtId);
    float charge = 0.;
    double hittime = 1.e4;
    if(calibChannel) {
      for(unsigned int i=0; i<calibChannel->time().size(); i++){
        nhit++;
        if(calibChannel->time(i)>m_timeCutLow && calibChannel->time(i)<m_timeCutHigh){
          if(charge<calibChannel->charge(i)) charge = calibChannel->charge(i);
          if(hittime>calibChannel->time(i)) hittime = calibChannel->time(i);
        }
      } 
    }
 //   cout<<"nhit: "<<nhit<<" "<<ring<<" "<<col<<" "<<charge<<endl;
    m_obsq.push_back(charge);
  }
  return StatusCode::SUCCESS;
}
const InterfaceID & ILikelihoodTool::interfaceID ( ) [static, inherited]

Retrieve interface ID.

Definition at line 8 of file ILikelihoodTool.cc.

{ 
    return IID_ILikelihoodTool; 
}

Member Data Documentation

Definition at line 70 of file LikelihoodTool.h.

Definition at line 73 of file LikelihoodTool.h.

Definition at line 76 of file LikelihoodTool.h.

Definition at line 79 of file LikelihoodTool.h.

Definition at line 83 of file LikelihoodTool.h.

Definition at line 86 of file LikelihoodTool.h.

Definition at line 88 of file LikelihoodTool.h.

int LikelihoodTool::m_LLFMode = 1 [static]

Definition at line 93 of file LikelihoodTool.h.

bool LikelihoodTool::m_addTimeLLF = 0 [static]

Definition at line 97 of file LikelihoodTool.h.

double LikelihoodTool::m_pmtChgCut = 0.0 [static]

Definition at line 101 of file LikelihoodTool.h.

TGraph * LikelihoodTool::m_gfangle = 0 [static]

Definition at line 102 of file LikelihoodTool.h.

TGraph * LikelihoodTool::m_gfdistance = 0 [static]

Definition at line 103 of file LikelihoodTool.h.

std::string LikelihoodTool::m_src

Definition at line 113 of file LikelihoodTool.h.

Definition at line 114 of file LikelihoodTool.h.

Definition at line 115 of file LikelihoodTool.h.

Definition at line 116 of file LikelihoodTool.h.

Definition at line 117 of file LikelihoodTool.h.

Definition at line 118 of file LikelihoodTool.h.

Definition at line 119 of file LikelihoodTool.h.

Definition at line 120 of file LikelihoodTool.h.

Definition at line 121 of file LikelihoodTool.h.

Definition at line 122 of file LikelihoodTool.h.

Definition at line 123 of file LikelihoodTool.h.

double LikelihoodTool::pmtX[200]

Definition at line 124 of file LikelihoodTool.h.

double LikelihoodTool::pmtY[200]

Definition at line 125 of file LikelihoodTool.h.

double LikelihoodTool::pmtZ[200]

Definition at line 126 of file LikelihoodTool.h.

Definition at line 127 of file LikelihoodTool.h.

Definition at line 128 of file LikelihoodTool.h.

Definition at line 129 of file LikelihoodTool.h.

Definition at line 133 of file LikelihoodTool.h.

Definition at line 134 of file LikelihoodTool.h.

Definition at line 136 of file LikelihoodTool.h.

Definition at line 137 of file LikelihoodTool.h.

std::vector< float > LikelihoodTool::m_obsq [static]

Definition at line 139 of file LikelihoodTool.h.

std::vector< float > LikelihoodTool::m_firstT [static]

Definition at line 140 of file LikelihoodTool.h.


The documentation for this class was generated from the following files:
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:13:37 for Likelihood by doxygen 1.7.4