In This Package:

QMLFTool Class Reference

#include <QMLFTool.h>

Inheritance diagram for QMLFTool:

[legend]Collaboration diagram for QMLFTool:[legend]List of all members.


Public Types
	SUCCESS
	NO_INTERFACE
	VERSMISMATCH
	LAST_ERROR
	SUCCESS
	NO_INTERFACE
	VERSMISMATCH
	LAST_ERROR
enum	Status
Public Member Functions
	QMLFTool (const std::string &type, const std::string &name, const IInterface *parent)
virtual	~QMLFTool ()
virtual StatusCode	reconstruct (const DayaBay::CalibReadout &, DayaBay::RecTrigger &)
virtual StatusCode	initialize ()
virtual StatusCode	finalize ()
INTupleSvc *	ntupleSvc () const
INTupleSvc *	evtColSvc () const
IDataProviderSvc *	detSvc () const
IDataProviderSvc *	evtSvc () const
IIncidentSvc *	incSvc () const
IChronoStatSvc *	chronoSvc () const
IHistogramSvc *	histoSvc () const
IAlgContextSvc *	contextSvc () const
DataObject *	put (IDataProviderSvc svc, DataObject object, const std::string &address, const bool useRootInTES=true) const
DataObject *	put (DataObject *object, const std::string &address, const bool useRootInTES=true) const
Gaudi::Utils::GetData< TYPE >::return_type	get (IDataProviderSvc *svc, const std::string &location, const bool useRootInTES=true) const
Gaudi::Utils::GetData< TYPE >::return_type	get (const std::string &location, const bool useRootInTES=true) const
TYPE *	getDet (IDataProviderSvc *svc, const std::string &location) const
TYPE *	getDet (const std::string &location) const
bool	exist (IDataProviderSvc *svc, const std::string &location, const bool useRootInTES=true) const
bool	exist (const std::string &location, const bool useRootInTES=true) const
bool	existDet (IDataProviderSvc *svc, const std::string &location) const
bool	existDet (const std::string &location) const
TYPE *	getOrCreate (IDataProviderSvc *svc, const std::string &location, const bool useRootInTES=true) const
TYPE *	getOrCreate (const std::string &location, const bool useRootInTES=true) const
TOOL *	tool (const std::string &type, const std::string &name, const IInterface *parent=0, bool create=true) const
TOOL *	tool (const std::string &type, const IInterface *parent=0, bool create=true) const
SERVICE *	svc (const std::string &name, const bool create=true) const
IUpdateManagerSvc *	updMgrSvc () const
IDataProviderSvc *	fastContainersSvc () const
StatusCode	Error (const std::string &msg, const StatusCode st=StatusCode::FAILURE, const size_t mx=10) const
StatusCode	Warning (const std::string &msg, const StatusCode st=StatusCode::FAILURE, const size_t mx=10) const
StatusCode	Print (const std::string &msg, const StatusCode st=StatusCode::SUCCESS, const MSG::Level lev=MSG::INFO) const
StatusCode	Assert (const bool ok, const std::string &message="", const StatusCode sc=StatusCode(StatusCode::FAILURE, true)) const
StatusCode	Assert (const bool ok, const char *message, const StatusCode sc=StatusCode(StatusCode::FAILURE, true)) const
StatusCode	Exception (const std::string &msg, const GaudiException &exc, const StatusCode sc=StatusCode(StatusCode::FAILURE, true)) const
StatusCode	Exception (const std::string &msg, const std::exception &exc, const StatusCode sc=StatusCode(StatusCode::FAILURE, true)) const
StatusCode	Exception (const std::string &msg="no message", const StatusCode sc=StatusCode(StatusCode::FAILURE, true)) const
MsgStream &	msgStream (const MSG::Level level) const
MsgStream &	always () const
MsgStream &	fatal () const
MsgStream &	err () const
MsgStream &	error () const
MsgStream &	warning () const
MsgStream &	info () const
MsgStream &	debug () const
MsgStream &	verbose () const
MsgStream &	msg () const
const Statistics &	counters () const
StatEntity &	counter (const std::string &tag) const
MSG::Level	msgLevel () const
bool	msgLevel (const MSG::Level level) const
void	resetMsgStream () const
bool	typePrint () const
bool	propsPrint () const
bool	statPrint () const
bool	errorsPrint () const
long	printStat (const MSG::Level level=MSG::ALWAYS) const
long	printErrors (const MSG::Level level=MSG::ALWAYS) const
long	printProps (const MSG::Level level=MSG::ALWAYS) const
void	registerCondition (const std::string &condition, StatusCode(CallerClass::*mf)()=NULL)
void	registerCondition (const std::string &condition, CondType &condPtrDest, StatusCode(CallerClass::mf)()=NULL)
void	registerCondition (char condition, StatusCode(CallerClass::mf)()=NULL)
void	registerCondition (TargetClass condition, StatusCode(CallerClass::mf)()=NULL)
StatusCode	runUpdate ()
TransientFastContainer< T > *	getFastContainer (const std::string &location, typename TransientFastContainer< T >::size_type initial=0)
StatusCode	release (const IInterface *interface) const
virtual unsigned long	release ()
const std::string &	context () const
const std::string &	rootInTES () const
double	globalTimeOffset () const
virtual StatusCode	queryInterface (const InterfaceID &riid, void **ppvUnknown)
virtual unsigned long	addRef ()
virtual const std::string &	name () const
virtual const std::string &	type () const
virtual const IInterface *	parent () const
virtual StatusCode	configure ()
virtual StatusCode	start ()
virtual StatusCode	stop ()
virtual StatusCode	terminate ()
virtual StatusCode	reinitialize ()
virtual StatusCode	restart ()
virtual Gaudi::StateMachine::State	FSMState () const
virtual Gaudi::StateMachine::State	targetFSMState () const
virtual StatusCode	sysInitialize ()
virtual StatusCode	sysStart ()
virtual StatusCode	sysStop ()
virtual StatusCode	sysFinalize ()
virtual StatusCode	sysReinitialize ()
virtual StatusCode	sysRestart ()
virtual StatusCode	setProperty (const Property &p)
virtual StatusCode	setProperty (const std::string &s)
virtual StatusCode	setProperty (const std::string &n, const std::string &v)
StatusCode	setProperty (const std::string &name, const TYPE &value)
virtual StatusCode	getProperty (Property *p) const
virtual const Property &	getProperty (const std::string &name) const
virtual StatusCode	getProperty (const std::string &n, std::string &v) const
virtual const std::vector< Property * > &	getProperties () const
PropertyMgr *	getPropertyMgr ()
ISvcLocator *	serviceLocator () const
ISvcLocator *	svcLoc () const
IMessageSvc *	msgSvc () const
IToolSvc *	toolSvc () const
StatusCode	setProperties ()
StatusCode	service (const std::string &name, T *&svc, bool createIf=true) const
StatusCode	service (const std::string &type, const std::string &name, T *&svc) const
void	declInterface (const InterfaceID &, void *)
Property *	declareProperty (const std::string &name, T &property, const std::string &doc="none") const
Property *	declareRemoteProperty (const std::string &name, IProperty *rsvc, const std::string &rname="") const
IAuditorSvc *	auditorSvc () const
IMonitorSvc *	monitorSvc () const
void	declareInfo (const std::string &name, const T &var, const std::string &desc) const
void	declareInfo (const std::string &name, const std::string &format, const void *var, int size, const std::string &desc) const
virtual const std::string &	type () const =0
virtual const IInterface *	parent () const =0
virtual StatusCode	configure ()=0
virtual StatusCode	start ()=0
virtual StatusCode	stop ()=0
virtual StatusCode	terminate ()=0
virtual StatusCode	reinitialize ()=0
virtual StatusCode	restart ()=0
virtual Gaudi::StateMachine::State	FSMState () const =0
virtual StatusCode	sysInitialize ()=0
virtual StatusCode	sysStart ()=0
virtual StatusCode	sysStop ()=0
virtual StatusCode	sysFinalize ()=0
virtual StatusCode	sysReinitialize ()=0
virtual StatusCode	sysRestart ()=0
virtual unsigned long	refCount () const =0
virtual const std::string &	name () const =0
virtual StatusCode	queryInterface (const InterfaceID &riid, void **ppvInterface)=0
virtual unsigned long	addRef ()=0
virtual unsigned long	release ()=0
Static Public Member Functions
static const InterfaceID &	interfaceID ()
static const InterfaceID &	interfaceID ()
static const InterfaceID &	interfaceID ()
static const InterfaceID &	interfaceID ()
Public Attributes
std::string	m_cableSvcName
std::string	m_pmtGeomSvcName
std::string	m_pmtCalibDataSvcName
std::string	m_opLocation
std::string	m_geomLocation
double	m_ATTL
double	m_TopRef
double	m_BotRef
double	m_LightYield
double	m_PhotoCathodeArea
double	m_TopRefZpos
double	m_BotRefZpos
	SUCCESS
	NO_INTERFACE
	VERSMISMATCH
	LAST_ERROR
Static Public Attributes
static ICableSvc *	m_cableSvc = 0
static IPmtGeomInfoSvc *	m_pmtGeomSvc = 0
static ICalibDataSvc *	m_pmtCalibDataSvc = 0
static DayaBay::CalibReadoutPmtCrate *	m_calibCrate = 0
static ExpChgMap	m_expChgMap
static int	m_LLFMode = 1
static bool	m_addTimeLLF = 0
static double	m_pmtChgCut = 0.0
static bool	m_recExtInfo = 0
static OpPara *	m_opPara = 0
static GeomPara *	m_geomPara = 0
static UsedFeeChannelMap	m_timeLLFUsedChannels
static UsedFeeChannelMap	m_chgLLFUsedChannels
Protected Types
typedef std::map< std::string, StatEntity >	Statistics
typedef std::map< std::string, unsigned int >	Counter
typedef std::vector< IAlgTool * >	AlgTools
typedef std::pair< IInterface *, std::string >	ServiceEntry
typedef std::vector< ServiceEntry >	Services
Protected Member Functions
StatusCode	releaseTool (const IAlgTool *tool) const
StatusCode	releaseSvc (const IInterface *svc) const
int	outputLevel () const
virtual unsigned long	refCount () const
IntegerProperty &	outputLevelProperty ()
void	initOutputLevel (Property &prop)
Static Protected Attributes
static const bool	IgnoreRootInTES
static const bool	UseRootInTES
Private Types
typedef std::map< DayaBay::FeeChannelId, bool >	UsedFeeChannelMap
typedef std::map< DayaBay::FeeChannelId, double >	ExpChgMap
Private Member Functions
StatusCode	getOpPara (std::string &data_src)
StatusCode	getGeomPara (std::string &data_src)
StatusCode	printPara ()
StatusCode	setUsedChannels (const DayaBay::CalibReadoutPmtCrate , DayaBay::RecTrigger )
Static Private Member Functions
static void	NCLL_fcn (int &npar, double grad, double &fval, double xval, int iflag)

Detailed Description

Definition at line 40 of file QMLFTool.h.

Member Typedef Documentation

typedef std::map<DayaBay::FeeChannelId, bool> QMLFTool::UsedFeeChannelMap [private]

Definition at line 66 of file QMLFTool.h.

typedef std::map<DayaBay::FeeChannelId, double> QMLFTool::ExpChgMap [private]

Definition at line 67 of file QMLFTool.h.

Constructor & Destructor Documentation

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

Definition at line 50 of file QMLFTool.cc.

00053   :  GaudiTool(type, name, parent)
00054   //, m_recTrigExt(0)
00055 {
00056   declareInterface<IReconTool>(this);
00057 
00058   declareProperty("CableSvcName", m_cableSvcName = "StaticCableSvc",
00059       "Name of service to map between detector, hardware, and electronic IDs");
00060   declareProperty("PmtGeomSvcName", m_pmtGeomSvcName = "PmtGeomInfoSvc", 
00061       "Name of Pmt Geometry Information Service");
00062   declareProperty("PmtCalibDataSvcName", m_pmtCalibDataSvcName = "StaticCalibDataSvc", 
00063       "Name of Pmt Calib Data Information Service");
00064 
00065   // declare property for those optical parameters
00066   declareProperty("LLFMode", m_LLFMode = 1, "Likelihood function mode");
00067   declareProperty("addTimeLLF", m_addTimeLLF = 0, "Add time likelihood or not");
00068   declareProperty("pmtChgCut", m_pmtChgCut = 12., "The charge cut for PMTs selection criteria");
00069   declareProperty("recExtInfo", m_recExtInfo = 0, "Extended information from reconstruction");
00070   declareProperty("opLocation", m_opLocation="DetDesc", "Optical parameters location");
00071   declareProperty("geomLocation", m_geomLocation="DetDesc", "Optical parameters location");
00072   declareProperty("AbsLength", m_ATTL=0., "Effective absorption length");
00073   declareProperty("TopRefZ", m_TopRefZpos=0., "Z position of top reflector");
00074   declareProperty("BotRefZ", m_BotRefZpos=0., "Z position of bottom reflector");
00075   declareProperty("TopReflectivity", m_TopRef=0., "Reflectivity of top reflector");
00076   declareProperty("BotReflectivity", m_BotRef=0., "Reflectivity of bottom reflector");
00077 }

QMLFTool::~QMLFTool ( ) [virtual]

Definition at line 79 of file QMLFTool.cc.

00079 {}

Member Function Documentation

StatusCode QMLFTool::reconstruct	(	const DayaBay::CalibReadout &	,
		DayaBay::RecTrigger &
	)			`[virtual]`

Implements IReconTool.

Definition at line 108 of file QMLFTool.cc.

00110 {
00111   if( !calibReadout.detector().isAD() ) {
00112     debug() << "Not an AD readout; ignoring detector " 
00113           << calibReadout.detector().detName() << endreq;
00114     recTrigger.setPositionStatus( ReconStatus::kNotProcessed);
00115     recTrigger.setEnergyStatus( ReconStatus::kNotProcessed);
00116     return StatusCode::SUCCESS;
00117   }
00118 
00119   //if(m_recExtInfo) {
00120   //   m_recTrigExt = new DayaBay::RecTriggerExt();
00121   //   m_expChgMap.clear();
00122   //}
00123    DayaBay::CalibReadout* calibTest 
00124      = const_cast<DayaBay::CalibReadout*> (&calibReadout);
00125    DayaBay::CalibReadoutPmtCrate* calibCrate 
00126      = dynamic_cast<DayaBay::CalibReadoutPmtCrate*> (calibTest);
00127    m_calibCrate = calibCrate;
00128 
00129    if(!m_calibCrate) {
00130      error() << "Incorrect type of readout crate for detector "
00131              << calibReadout.detector().detName() << endreq;
00132      recTrigger.setPositionStatus(ReconStatus::kBadReadout);
00133      recTrigger.setEnergyStatus(ReconStatus::kBadReadout);
00134      return StatusCode::FAILURE;
00135    }
00136 
00137    info() << "QMLFTool: test output" << endreq;
00138 
00139    // determine FeeChannels that used for reconstruction
00140    setUsedChannels(calibCrate, &recTrigger);
00141    
00142    //
00143    // Create TMinuit
00144    const int Npar = 4;
00145    TMinuit *recMinuit = new TMinuit(Npar);
00146    int ierflg;
00147 
00148    recMinuit->SetFCN(QMLFTool::NCLL_fcn);
00149    recMinuit->SetPrintLevel(1);
00150 
00151    // Initial parameter   
00152    double xini, yini, zini, eini, tini;
00153    xini = recTrigger.position().x();
00154    yini = recTrigger.position().y();
00155    zini = recTrigger.position().z();
00156    eini = recTrigger.energy();
00157    tini = recTrigger.triggerTime().GetNanoSec();
00158 
00159    if(eini==0.0) {
00160      error() << "No initial energy value was set!" << endreq;
00161    }
00162    if(xini==0.0&&yini==0.0&&zini==0.0) {
00163      warning() << "No initial vertex was set!" << endreq;
00164    }
00165    
00166    // Initial TMinuit Fitting Parameter
00167    recMinuit->mnparm(0, "xpos",   xini,  100.0, 0.0, 0.0, ierflg);
00168    recMinuit->mnparm(1, "ypos",   yini,  100.0, 0.0, 0.0, ierflg);
00169    recMinuit->mnparm(2, "zpos",   zini,  10.0, m_geomPara->m_botRefZ, m_geomPara->m_topRefZ, ierflg);
00170    recMinuit->mnparm(3, "energy", eini,   0.1, 0.001, 1000.0, ierflg);
00171    if(m_addTimeLLF==1) 
00172      recMinuit->mnparm(4, "T0", tini,  0.5, 0.0, 0.0, ierflg);
00173    // Set parameters for Minuit minimization
00174    int    err;
00175    double standin = 0.0;
00176    double maxCalls = 2000;
00177    double tolerance = 1.0e-3;
00178    double arglist[2];
00179 
00180    // Call FCN
00181    arglist[0] = 1;
00182    recMinuit->mnexcm("CALL FCN", arglist, 1, ierflg);
00183 
00184    // No Warnings
00185    recMinuit->mnexcm("SET NOW", &standin, 0, err);
00186    // Set error Definition
00187    // 1 for Chi square
00188    // 0.5 for negative log likelihood
00189    if(m_LLFMode==1)
00190      recMinuit->SetErrorDef(1);
00191    else 
00192      recMinuit->SetErrorDef(0.5);
00193    // Minimization strategy
00194    // 1 standard;  2 try to improve minimum (slower)
00195    arglist[0] = 2;
00196    recMinuit->mnexcm("SET STR", arglist, 1, ierflg); 
00197   
00198    //Do Minuit fit!  
00199    // set max iterations: recMinuit->SetMaxIteration(1000); recMinuit->Migrad();
00200    arglist[0] = maxCalls;
00201    arglist[1] = tolerance;
00202    recMinuit->mnexcm("MIGrad", arglist, 2, err);
00203    //recMinuit->mnexcm("MINIMIZE", arglist, 2, err);
00204    //recMinuit->mnexcm("SIMPlex", arglist, 2, err);
00205 
00206    //Scan is usefull to check FCN
00207    // arglist[0] = 0;
00208    // recMinuit->mnexcm("SCAN", arglist, 1, ierflg);
00209 
00210    //Collect the output
00211    //How to draw the n-sigma contour of a Minuit fit?
00212    //http://root.cern.ch/cgi-bin/print_hit_bold.pl/root/html508/examples/fitcont.C.html?minuit#first_hit
00213 
00214    // record fitter's status
00215    double ln0, edm, errdef;
00216    int nvpar, nparx, icstat;
00217    recMinuit->mnstat(ln0, edm, errdef, nvpar, nparx, icstat);
00218    
00219    // Get Fitted Parameter
00220    double x_fit, x_fit_err;
00221    double y_fit, y_fit_err;
00222    double z_fit, z_fit_err;
00223    double e_fit, e_fit_err;
00224 
00225    recMinuit->GetParameter(0, x_fit, x_fit_err);
00226    recMinuit->GetParameter(1, y_fit, y_fit_err);
00227    recMinuit->GetParameter(2, z_fit, z_fit_err);
00228    recMinuit->GetParameter(3, e_fit, e_fit_err);
00229 
00230    // error matrix; mnemat()
00231    double emat[Npar][Npar];
00232    CLHEP::HepMatrix matrix(Npar, Npar);
00233    recMinuit->mnemat(&emat[0][0], Npar);
00234    for(int row=0; row<Npar; row++) {
00235      for(int col=0; col<Npar; col++) {
00236        matrix[row][col] = emat[row][col];
00237      }
00238    }
00239 
00240    CLHEP::HepLorentzVector vertex(x_fit, y_fit, z_fit, 0.);
00241    recTrigger.setPosition(vertex);
00242    recTrigger.setEnergy(e_fit);
00243    recTrigger.setPositionQuality(ln0);
00244    recTrigger.setEnergyQuality(ln0);
00245    recTrigger.setErrorMatrix(matrix);
00246 
00247    if(m_recExtInfo) {
00248    arglist[0] = 3;
00249    recMinuit->mnexcm("CALL FCN", arglist, 1, ierflg);
00250    //m_recTrigExt->setExpChgMap(m_expChgMap);
00251    //recTrigger.setRecTrigExt(m_recTrigExt);
00252    }
00253    delete recMinuit;
00254        
00255    return StatusCode::SUCCESS;
00256 }

StatusCode QMLFTool::initialize ( ) [virtual]

Reimplemented from GaudiTool.

Definition at line 81 of file QMLFTool.cc.

00082 {
00083   // Get Cable Service
00084   m_cableSvc = svc<ICableSvc>(m_cableSvcName, true);
00085   // Get Pmt Geometry Service
00086   m_pmtGeomSvc = svc<IPmtGeomInfoSvc>(m_pmtGeomSvcName, true);
00087 
00088   // Get Pmt Calib Data Service
00089   m_pmtCalibDataSvc = svc<ICalibDataSvc>(m_pmtCalibDataSvcName, true);
00090 
00091   // Get Optical Parameters
00092   getOpPara(m_opLocation);
00093 
00094   // Get Geometry Parameters
00095   getGeomPara(m_geomLocation);
00096 
00097   // Show Optical and Geom Parameters
00098   printPara();
00099 
00100   return StatusCode::SUCCESS;
00101 }

StatusCode QMLFTool::finalize ( ) [virtual]

Reimplemented from GaudiTool.

Definition at line 103 of file QMLFTool.cc.

00104 {
00105     return StatusCode::SUCCESS;
00106 }

void QMLFTool::NCLL_fcn	(	int &	npar,
		double *	grad,
		double &	fval,
		double *	xval,
		int	iflag
	)			`[static, private]`

Definition at line 258 of file QMLFTool.cc.

00258                                                                                       {
00259 
00260   if(iflag==1){
00261     // Initialization
00262     // read input data, calculate any necessary constants, etc
00263   }
00264   if(iflag==2){
00265     // Compute derivatives (GRAD), the first derivatives of FVAL
00266     // (this is optional)
00267   }
00268 
00269   // Always calculate the value of the function, FVAL,
00270   // which is usually a chisquare or log likelihood.
00271   
00272     double ECAL;
00273     // Global coefficient (QE, Yield...) 
00274     ECAL = m_opPara->m_photoCathodeArea/4.0/TMath::Pi()*m_opPara->m_yield;
00275     ECAL = ECAL * 0.20; // Nominal 'absolute' QE: 0.20
00276     //ECAL = 1000.;
00277 
00279     double S_angle, dist;
00280     double expected_q, qhit, cosa_tmp;
00281     double sum_exp, AA, BB, LLF;
00282     double XD, YD, ZD, ED, ZF, dx, dy, dz, T0;
00283     double chg_likelihood, time_likelihood;
00284     int ipmt, npe;
00285 
00287     XD = xval[0];
00288     YD = xval[1];
00289     ZD = xval[2];
00290     ED = xval[3]*ECAL;
00291     if(m_addTimeLLF==1) T0 = xval[4];
00292 
00293     double sum_act = 0.0;
00294     sum_exp = 0.0;
00295     chg_likelihood = 0.0;
00296     time_likelihood = 0.0;
00297     double c_n = 200.; // light speed in liquid, mm/s
00298 
00299     DayaBay::CalibReadoutPmtCrate::PmtChannelReadouts::const_iterator cpcrIter,
00300       done = (m_calibCrate->channelReadout()).end();
00301     
00302     Context context = m_calibCrate->header()->context();
00303     int task = 0;
00304     ServiceMode svcMode(context, task);
00305 
00306     Site::Site_t site = m_calibCrate->detector().site();
00307     DetectorId::DetectorId_t  detid = m_calibCrate->detector().detectorId();
00308 
00309     // Loop over calib readout channels
00310     // Remember: only non-zero readout channels are stored the map
00311     double sumQ = 0.0, Xc, Yc, Zc;
00312     for(int local_id=0; local_id<192; local_id++) {
00313       int ring = local_id/24 + 1;
00314       int col  = local_id%24 + 1;
00315 
00316       AdPmtSensor adPmtId(ring, col, site, detid);
00317       FeeChannelId channelId = m_cableSvc->feeChannelId(adPmtId, svcMode);
00318 
00319       // determine whether to use this channel
00320       if(iflag!=3&&m_chgLLFUsedChannels[channelId]==false) continue;
00321 
00322       DayaBay::CalibReadoutPmtChannel* cpcrPtr = m_calibCrate->sensor(adPmtId);
00323       if(cpcrPtr) { // map has the record of this channelId
00324         qhit = cpcrPtr->maxCharge();
00325       }else { // map has no record for this channel, 
00326               // yet it still has contribution to likelihood
00327         qhit = 0.0;
00328       }
00329 
00330       unsigned int pmtid = adPmtId.fullPackedData();
00331       if(adPmtId.site() == Site::kSAB) {
00332         pmtid = DayaBay::AdPmtSensor(adPmtId.ring(),adPmtId.column(),Site::kDayaBay,adPmtId.detectorId()).fullPackedData();
00333       }
00334       CLHEP::Hep3Vector pos =  m_pmtGeomSvc->get(pmtid)->localPosition();
00335       CLHEP::Hep3Vector norm = m_pmtGeomSvc->get(pmtid)->localDirection();
00336 
00337       // Get pmt data from ICalibDataSvc
00338       const DayaBay::DetectorSensor& sensDetId 
00339         = m_cableSvc->sensor(channelId, svcMode);
00340       const DayaBay::PmtCalibData* pmtData =
00341         m_pmtCalibDataSvc->pmtCalibData(sensDetId, svcMode);
00342 
00343       AA = norm.x();
00344       BB = norm.y();
00345 
00346       dx = XD - pos.x();
00347       dy = YD - pos.y();
00348       dz = ZD - pos.z();
00349 
00350       dist = dx*dx + dy*dy + dz*dz;
00351       if(dist<0.0001) dist = 0.0001;
00352       dist = sqrt(dist);
00353 
00354       cosa_tmp = ( AA*dx + BB*dy ) / dist;
00355       if(cosa_tmp>1.0) cosa_tmp = 1.0;
00356 
00357       // whether to use time information of this channel
00358       if(m_addTimeLLF==1) {
00359         if(m_timeLLFUsedChannels[channelId]==true){
00360           // get the earliest TDC time as the 'first hit' time, 
00361           // calculate the time likelihood
00362           double hittime = cpcrIter->earliestTime();
00363           time_likelihood = time_likelihood 
00364                   - 2*log(1./TMath::Sqrt(2*TMath::Pi())/pmtData->m_timeSpread)
00365                   + TMath::Power((hittime-pmtData->m_timeOffset-dist/c_n-T0)/pmtData->m_timeSpread, 2);
00366         }
00367       }
00368 
00369       S_angle = 0.378 + 0.5099*cosa_tmp + 0.1131*pow(cosa_tmp,2);
00370       if(S_angle<1.0e-5) S_angle = 1.0e-5;
00371 
00372       expected_q = S_angle*exp(-dist/m_opPara->m_attl)/TMath::Power(dist, 2);
00373 
00375 
00376       int ref_deg;
00377       for(ref_deg=0; ref_deg<10; ref_deg++){
00378      
00379         if(ref_deg==0)
00380           ZF =  2*m_geomPara->m_botRefZ - xval[2];   // 1st-order reflection
00381         else if(ref_deg==1)
00382           ZF =  2*m_geomPara->m_topRefZ - xval[2];   // 1st-order reflection
00383         else if(ref_deg==2)
00384           ZF = 2*m_geomPara->m_botRefZ - 2*m_geomPara->m_topRefZ + xval[2]; // 2nd-order reflection
00385         else if(ref_deg==3)
00386           ZF = 2*m_geomPara->m_topRefZ - 2*m_geomPara->m_botRefZ + xval[2]; // 2nd-order reflection
00387         else if(ref_deg==4)
00388           ZF = 4*m_geomPara->m_botRefZ - 2*m_geomPara->m_topRefZ - xval[2]; // 3rd-order reflection
00389         else if(ref_deg==5)
00390           ZF = 4*m_geomPara->m_topRefZ - 2*m_geomPara->m_botRefZ - xval[2]; // 3rd-order reflection
00391         else if(ref_deg==6)
00392           ZF = 4*m_geomPara->m_botRefZ - 4*m_geomPara->m_topRefZ + xval[2]; // 4th-order reflection
00393         else if(ref_deg==7)
00394           ZF = 4*m_geomPara->m_topRefZ - 4*m_geomPara->m_botRefZ + xval[2]; // 4th-order reflection
00395         else if(ref_deg==8)
00396           ZF = 6*m_geomPara->m_botRefZ - 4*m_geomPara->m_topRefZ - xval[2]; // 5th-order reflection
00397         else
00398           ZF = 6*m_geomPara->m_topRefZ - 4*m_geomPara->m_botRefZ - xval[2]; // 5th-order reflection
00399 
00400         dz = ZF - pos.z();
00401 
00402         dist = dx*dx + dy*dy + dz*dz;
00403         if(dist<0.0001) dist = 0.0001;
00404         dist = sqrt(dist);
00405 
00406         cosa_tmp = ( AA*dx + BB*dy ) / dist;
00407         if(cosa_tmp>1.0) cosa_tmp = 1.0;
00408 
00409         S_angle = 0.378 + 0.5099*cosa_tmp + 0.1131*pow(cosa_tmp,2);
00410         if(S_angle<1.0e-5) S_angle = 1.0e-5;
00411 
00412         int index = (int)(ref_deg/2.0) + 1;
00413         int index_top, index_bot;
00414         if(ref_deg == 0) 
00415           index_top = 0;
00416         else 
00417           index_top = (int)((ref_deg-1)/4.0) + 1;
00418 
00419         index_bot = index - index_top;
00420 
00421         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);
00422       }
00423 
00425 
00426       double eff = pmtData->m_efficiency;
00427 
00428       expected_q = ED * expected_q * eff;
00429       sum_exp = sum_exp + expected_q;
00430       sum_act = sum_act + qhit;
00431 
00432       if(iflag==3) {
00433         //after the fit is finished, get Exp_Q for each PMT
00434         if(m_recExtInfo) m_expChgMap[channelId] = expected_q;
00435 
00436       }
00437       else {
00438         if(m_LLFMode!=1){           
00439 
00440 
00441         }
00442         else
00443         {
00445           if(qhit!=0.0) 
00446             chg_likelihood = chg_likelihood + 2*(expected_q-qhit) + 2*log(qhit/expected_q)*qhit;
00447           else 
00448             chg_likelihood = chg_likelihood + 2*expected_q;
00449         }
00450       }
00451     
00452     } // for(local_id....
00453 
00454     fval = chg_likelihood;
00455     if(m_addTimeLLF==1) fval += time_likelihood;
00456 }

StatusCode QMLFTool::getOpPara ( std::string & data_src ) [private]

Definition at line 458 of file QMLFTool.cc.

00458                                                   {
00459 
00460   // construct OpPara, default value
00461   m_opPara = new OpPara(1.25E+04, 0.98, 0.98, 9000, 3.14*103*103);
00462 
00463   if(data_src=="DetDesc") {
00464 
00465     // Get Optical Parameters from Det/DetDesc
00466     // A tabular was defined in DetDesc, but only one value was get here
00467     // FIXME
00468     // Get the Detector Data Service to access TDS:
00469     // This should be done in your tools initialize()
00470     // and save m_DDS as a data member
00471     IDataProviderSvc* dds = 0;
00472     StatusCode sc = service( "DetectorDataSvc", dds, true );
00473     if (sc.isFailure()) {
00474         // ... handle error ...
00475     }
00476     // Later, look up some material.  If this material is always used you
00477     // should look it up once in initialize() and save it for later use
00478 
00479     // absorption length
00480     std::string location = "/dd/Materials/GdDopedLS";
00481     Material* gdls = GaudiCommon<AlgTool>::get<Material>(dds,location);
00482     
00483     Material::Tables& gdls_tab = gdls->tabulatedProperties();
00484     Material::Tables::const_iterator mtIter = gdls_tab.begin(), 
00485       done = gdls_tab.end();
00486     TabulatedProperty::Table attl_vec;
00487     for(; mtIter!=done; mtIter++) {
00488       std::string type = (*mtIter)->type();
00489       if(type=="ABSLENGTH") {
00490         attl_vec = (*mtIter)->table(); 
00491       }
00492     }
00493     TabulatedProperty::Table::const_iterator tb_iter;
00494     int count = 0;
00495     double attl_tmp = 0.0, ref_tmp = 0.0;
00496     for(tb_iter=attl_vec.begin(); tb_iter!=attl_vec.end(); tb_iter++){
00497       if(tb_iter->first<3.18e-06 && tb_iter->first>2.07e-06 ){ // 390nm ~ 600nm
00498         count++; attl_tmp += tb_iter->second;
00499       }
00500     }
00501     attl_tmp = attl_tmp/count; // average value between 390nm~600nm
00502     m_opPara->m_attl = attl_tmp;
00503       
00504     // ESR surface location
00505     std::string topESR_location = "/dd/Geometry/AdDetails/AdSurfacesAll/ESRAirSurfaceTop";
00506     std::string botESR_location = "/dd/Geometry/AdDetails/AdSurfacesAll/ESRAirSurfaceBot";
00507 
00508     Surface* esrtop = GaudiCommon<AlgTool>::get<Surface>(dds, topESR_location);
00509     Surface* esrbot = GaudiCommon<AlgTool>::get<Surface>(dds, botESR_location);
00510     
00511     Surface::Tables& esrtop_tab = esrtop->tabulatedProperties();
00512     Surface::Tables& esrbot_tab = esrbot->tabulatedProperties();
00513 
00514     Surface::Tables::const_iterator sfIter; 
00515     TabulatedProperty::Table ref_vec;
00516 
00517     // look up Tables, get an average top ESR reflectivity
00518     for(sfIter=esrtop_tab.begin(); sfIter!=esrtop_tab.end(); sfIter++) {
00519       std::string type = (*sfIter)->type();
00520       if(type=="REFLECTIVITY") {
00521         ref_vec = (*sfIter)->table(); 
00522       }
00523     }
00524     count = 0;
00525     for(tb_iter=ref_vec.begin(); tb_iter!=ref_vec.end(); tb_iter++){
00526       if(tb_iter->first<3.18e-06 && tb_iter->first>2.07e-06 ){ // 390nm ~ 600nm
00527         count++; ref_tmp += tb_iter->second;
00528       }
00529     }
00530     ref_tmp = ref_tmp/count; // average value between 390nm~600nm
00531     m_opPara->m_topRef = ref_tmp;
00532  
00533     // look up Tables, get an average bot ESR reflectivity
00534     for(sfIter=esrbot_tab.begin(); sfIter!=esrbot_tab.end(); sfIter++) {
00535       std::string type = (*sfIter)->type();
00536       if(type=="REFLECTIVITY") {
00537         ref_vec = (*sfIter)->table(); 
00538       }
00539     }
00540     count = 0; ref_tmp = 0.0;
00541     for(tb_iter=ref_vec.begin(); tb_iter!=ref_vec.end(); tb_iter++){
00542       if(tb_iter->first<3.18e-06 && tb_iter->first>2.07e-06 ){ // 390nm ~ 600nm
00543         count++; ref_tmp += tb_iter->second;
00544       }
00545     }
00546     ref_tmp = ref_tmp/count; // average value between 390nm~600nm
00547     m_opPara->m_botRef = ref_tmp;
00548 
00549   } else if(data_src=="OpCalibSvc"){
00550 
00551   } else {
00552 
00553     warning() <<"No external optical data source, use default value."<< endreq;
00554   }
00555 
00556   // if set new values from declareProperty
00557   if(m_ATTL!=0.) m_opPara->m_attl = m_ATTL;
00558   if(m_TopRef!=0.) m_opPara->m_topRef = m_TopRef;
00559   if(m_BotRef!=0.) m_opPara->m_botRef = m_BotRef;
00560 
00561   return StatusCode::SUCCESS;
00562 }

StatusCode QMLFTool::getGeomPara ( std::string & data_src ) [private]

Definition at line 564 of file QMLFTool.cc.

00564                                                     {
00565 
00566   // construct GeomPara, default value
00567   m_geomPara = new GeomPara(2010., -2010.);
00568 
00569   if(data_src=="DetDesc") {
00570 
00571     IDataProviderSvc* dds = 0;
00572     StatusCode sc = service( "DetectorDataSvc", dds, true );
00573     if (sc.isFailure()) {
00574        error() << "No DetectorDataSvc available!" << endreq;
00575     }
00576     
00577     // ESR structure location
00578     std::string topEsrPath = "/dd/Structure/AdReflectorStructure/db-ad/db-ad1-topesr";
00579     std::string botEsrPath = "/dd/Structure/AdReflectorStructure/db-ad/db-ad1-botesr";
00580     SmartDataPtr<IDetectorElement> topEsr(dds,topEsrPath);
00581     SmartDataPtr<IDetectorElement> botEsr(dds,botEsrPath);
00582 
00583     // define a pointer to ICoordSysSvc to
00584     // use to cache the service during the job.
00585     ICoordSysSvc* css;
00586     sc = service( "CoordSysSvc", css, true );
00587 
00588     Gaudi::XYZPoint zero(0,0,0);
00589     // Set top ESR position
00590     Gaudi::XYZPoint globalPoint = topEsr->geometry()->toGlobal(zero);
00591     IDetectorElement* theAdDe = css->coordSysDE(globalPoint);
00592     Gaudi::XYZPoint adLocalPoint = theAdDe->geometry()->toLocal(globalPoint);
00593     m_geomPara->m_topRefZ = adLocalPoint.z();
00594 
00595     // Set bot ESR position
00596     globalPoint = botEsr->geometry()->toGlobal(zero);
00597     theAdDe = css->coordSysDE(globalPoint);
00598     adLocalPoint = theAdDe->geometry()->toLocal(globalPoint);
00599     m_geomPara->m_botRefZ = adLocalPoint.z();
00600         
00601   } else if(data_src=="OnsiteSuveyData"){
00602     // Get geometry info from onsite suvey data
00603 
00604   } else {
00605 
00606     warning() <<"Unknown geometry data source, use default value."<< endreq;
00607   }
00608 
00609   // if set new values from declareProperty
00610   if(m_TopRefZpos!=0.) m_geomPara->m_topRefZ = m_TopRefZpos;
00611   if(m_BotRefZpos!=0.) m_geomPara->m_botRefZ = m_BotRefZpos;
00612 
00613   return StatusCode::SUCCESS;
00614 }

StatusCode QMLFTool::printPara ( ) [private]

Definition at line 616 of file QMLFTool.cc.

00616                               {
00617   info() << "Op and Geom parameters used in reconstruction." << endreq;
00618   info() << "Attenuation length : " << m_opPara->m_attl <<"mm"<< endreq;
00619   info() << "Top reflectivity : " << m_opPara->m_topRef << endreq;
00620   info() << "Bottom reflectivity : " << m_opPara->m_botRef << endreq;
00621   info() << "Top reflector Z pos :" << m_geomPara->m_topRefZ <<"mm"<< endreq;
00622   info() << "Bot reflector Z pos :" << m_geomPara->m_botRefZ <<"mm"<< endreq;
00623 
00624   return StatusCode::SUCCESS;
00625 }

StatusCode QMLFTool::setUsedChannels	(	const DayaBay::CalibReadoutPmtCrate *	,
		DayaBay::RecTrigger *
	)			`[private]`

Definition at line 627 of file QMLFTool.cc.

00630 {
00631   m_chgLLFUsedChannels.clear();
00632   m_timeLLFUsedChannels.clear();
00633 
00634   DayaBay::CalibReadoutPmtCrate::PmtChannelReadouts::const_iterator cpcrIter,
00635     done = (calibCrate->channelReadout()).end();
00636     
00637   Context context = calibCrate->header()->context();
00638   int task = 0;
00639   ServiceMode svcMode(context, task);
00640 
00641   Site::Site_t site = calibCrate->detector().site();
00642   DetectorId::DetectorId_t  detid = calibCrate->detector().detectorId();
00643 
00644   // default is use pmts whose status are Good
00645   for(int local_id=0; local_id<192; local_id++) {
00646     int ring = local_id/24 + 1;
00647     int col  = local_id%24 + 1;
00648 
00649     AdPmtSensor adPmtId(ring, col, site, detid);
00650     FeeChannelId channelId = m_cableSvc->feeChannelId(adPmtId, svcMode);
00651     
00652     const DayaBay::DetectorSensor& sensDetId = m_cableSvc->sensor(channelId, svcMode);
00653     const DayaBay::PmtCalibData* pmtData =
00654       m_pmtCalibDataSvc->pmtCalibData(sensDetId, svcMode);
00655     
00656     if(pmtData->m_status==DayaBay::PmtCalibData::kGood) { 
00657       m_chgLLFUsedChannels[channelId] = true;
00658       m_timeLLFUsedChannels[channelId] = false;
00659     }
00660     else {
00661       m_chgLLFUsedChannels[channelId] = false;
00662       m_timeLLFUsedChannels[channelId] = false;
00663     }
00664 
00665     if(m_addTimeLLF==1) { 
00666       DayaBay::CalibReadoutPmtChannel* cpcrPtr = m_calibCrate->sensor(adPmtId);
00667       if(cpcrPtr) { // map has the record of this channelId
00668         double qhit = cpcrPtr->maxCharge();
00669         if(qhit>m_pmtChgCut)
00670           m_timeLLFUsedChannels[channelId] = true;
00671       }
00672     }
00673     
00674   }
00675   // Other options, depending on CalibPmtCrate data
00676   //recTrigger->setUsedChannels(m_chgLLFUsedChannels);
00677 
00678   return StatusCode::SUCCESS;
00679 }

Member Data Documentation

std::string QMLFTool::m_cableSvcName

Definition at line 72 of file QMLFTool.h.

ICableSvc * QMLFTool::m_cableSvc = 0 [static]

Definition at line 75 of file QMLFTool.h.

std::string QMLFTool::m_pmtGeomSvcName

Definition at line 78 of file QMLFTool.h.

IPmtGeomInfoSvc * QMLFTool::m_pmtGeomSvc = 0 [static]

Definition at line 81 of file QMLFTool.h.

std::string QMLFTool::m_pmtCalibDataSvcName

Definition at line 84 of file QMLFTool.h.

ICalibDataSvc * QMLFTool::m_pmtCalibDataSvc = 0 [static]

Definition at line 87 of file QMLFTool.h.

DayaBay::CalibReadoutPmtCrate * QMLFTool::m_calibCrate = 0 [static]

Definition at line 89 of file QMLFTool.h.

QMLFTool::ExpChgMap QMLFTool::m_expChgMap [static]

Definition at line 92 of file QMLFTool.h.

int QMLFTool::m_LLFMode = 1 [static]

Definition at line 97 of file QMLFTool.h.

bool QMLFTool::m_addTimeLLF = 0 [static]

Definition at line 101 of file QMLFTool.h.

double QMLFTool::m_pmtChgCut = 0.0 [static]

Definition at line 105 of file QMLFTool.h.

bool QMLFTool::m_recExtInfo = 0 [static]

Definition at line 109 of file QMLFTool.h.

OpPara * QMLFTool::m_opPara = 0 [static]

Definition at line 112 of file QMLFTool.h.

std::string QMLFTool::m_opLocation

Definition at line 113 of file QMLFTool.h.

GeomPara * QMLFTool::m_geomPara = 0 [static]

Definition at line 116 of file QMLFTool.h.

std::string QMLFTool::m_geomLocation

Definition at line 117 of file QMLFTool.h.

double QMLFTool::m_ATTL

Definition at line 119 of file QMLFTool.h.

double QMLFTool::m_TopRef

Definition at line 120 of file QMLFTool.h.

double QMLFTool::m_BotRef

Definition at line 121 of file QMLFTool.h.

double QMLFTool::m_LightYield

Definition at line 122 of file QMLFTool.h.

double QMLFTool::m_PhotoCathodeArea

Definition at line 123 of file QMLFTool.h.

double QMLFTool::m_TopRefZpos

Definition at line 126 of file QMLFTool.h.

double QMLFTool::m_BotRefZpos

Definition at line 127 of file QMLFTool.h.

QMLFTool::UsedFeeChannelMap QMLFTool::m_timeLLFUsedChannels [static]

Definition at line 130 of file QMLFTool.h.

QMLFTool::UsedFeeChannelMap QMLFTool::m_chgLLFUsedChannels [static]

Definition at line 131 of file QMLFTool.h.

The documentation for this class was generated from the following files:

/home/dayabay/installs/trunk_2011_04_11_opt/NuWa-trunk/dybgaudi/Reconstruction/AdRec/src/components/QMLFTool.h
/home/dayabay/installs/trunk_2011_04_11_opt/NuWa-trunk/dybgaudi/Reconstruction/AdRec/src/components/QMLFTool.cc

Generated on Mon Apr 11 20:42:46 2011 for AdRec by 1.4.7

In This Package:

QMLFTool Class Reference

Public Types

Public Member Functions

Static Public Member Functions

Public Attributes

Static Public Attributes

Protected Types

Protected Member Functions

Static Protected Attributes

Private Types

Private Member Functions

Static Private Member Functions

Detailed Description

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation