/** @file CalRawEnergyTool.cxx
@brief implementation of the class CalRawEnergyTool
$Header: /nfs/slac/g/glast/ground/cvs/CalRecon/src/CalRawEnergyTool.cxx,v 1.12 2005/04/12 16:42:37 chamont Exp $
*/
#include "IEnergyCorr.h"
#include "GaudiKernel/AlgTool.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/GaudiException.h"
#include "Event/TopLevel/EventModel.h"
#include "Event/Recon/CalRecon/CalCluster.h"
#include "Event/Recon/CalRecon/CalEventEnergy.h"
#include "CLHEP/Matrix/Vector.h"
/**
* @class CalRawEnergyTool
* @author CalRecon Rewrite Group
*
* This sets the "raw" energy for an event (after clustering)
*
* $Header: /nfs/slac/g/glast/ground/cvs/CalRecon/src/CalRawEnergyTool.h,v 0.1 2005/04/11 13:28:50 chamont Exp $
*/
class CalRawEnergyTool : public AlgTool, virtual public IEnergyCorr
{
public:
//! Constructor
CalRawEnergyTool( const std::string& type, const std::string& name, const IInterface* parent);
//! destructor
virtual ~CalRawEnergyTool() {};
StatusCode initialize();
// worker function to get the corrected energy
StatusCode doEnergyCorr( Event::CalCluster * ) ;
private:
/// Pointer to the Gaudi data provider service
IDataProviderSvc* m_dataSvc;
};
#include "GaudiKernel/DeclareFactoryEntries.h"
DECLARE_TOOL_FACTORY(CalRawEnergyTool) ;
CalRawEnergyTool::CalRawEnergyTool( const std::string & type,
const std::string & name,
const IInterface * parent )
: AlgTool(type,name,parent)
{
declareInterface<IEnergyCorr>(this) ;
}
StatusCode CalRawEnergyTool::initialize()
{
// This function does following initialization actions:
// - Sets up a pointer to the data service as a convenience for the tool
MsgStream log(msgSvc(), name());
StatusCode sc = StatusCode::SUCCESS;
log << MSG::INFO << "Initializing CalRawEnergyTool" <<endreq;
//Locate and store a pointer to the geometry service
IService* iService = 0;
//Locate and store a pointer to the data service which allows access to the TDS
if ((sc = serviceLocator()->getService("EventDataSvc", iService)).isFailure())
{
throw GaudiException("Service [EventDataSvc] not found", name(), sc);
}
m_dataSvc = dynamic_cast<IDataProviderSvc*>(iService);
return sc;
}
StatusCode CalRawEnergyTool::doEnergyCorr( Event::CalCluster* )
{
//Purpose and method:
//
// This function determines the overal "raw" energy for an event
// Note: This sums the energy contribution from all clusters and
// does a weighted average for centroid position and axis
//
// TDS input: CalClusters
// TDS output: CalEventEnergy
MsgStream lm(msgSvc(), name());
StatusCode sc = StatusCode::SUCCESS;
// Recover the collection of all clusters
SmartDataPtr<Event::CalClusterCol> calClusters(m_dataSvc, EventModel::CalRecon::CalClusterCol);
// Set up to loop over all clusters to get total raw energy
double rawEnergy = 0.;
double rawEneError = 0.;
HepVector posSum(3);
HepMatrix posWghtSum(3,3,0.);
HepVector axisSum(3);
HepMatrix axisWghtSum(3,3,0.);
// Do the loop and accumulate information
for(Event::CalClusterCol::iterator clusIter = calClusters->begin(); clusIter != calClusters->end(); clusIter++)
{
Event::CalCluster* cluster = *clusIter;
const Event::CalParams& params = cluster->getCalParams();
HepVector centroid(3);
centroid[0] = params.getCentroid().x();
centroid[1] = params.getCentroid().y();
centroid[2] = params.getCentroid().z();
HepVector axis(3);
axis[0] = params.getAxis().x();
axis[1] = params.getAxis().y();
axis[2] = params.getAxis().z();
rawEnergy += params.getEnergy();
rawEneError += params.getEnergyErr() * params.getEnergyErr();
HepMatrix posCovInv = params.getCentroidErrs();
int matInvErr = 0;
posCovInv.invert(matInvErr);
posWghtSum += posCovInv;
posSum += posCovInv * centroid;
HepMatrix axisCovInv = params.getAxisErrs();
axisCovInv.invert(matInvErr);
axisWghtSum += axisCovInv;
axisSum += axisCovInv * axis;
}
// Get new errors and weighted average centroid
int matInvErr = 0;
posWghtSum.invert(matInvErr);
posSum = posWghtSum * posSum;
// Get new errors and weighted average axis
axisWghtSum.invert(matInvErr);
axisSum = axisWghtSum * axisSum;
// New estimate of energy error
rawEneError = sqrt(rawEneError);
// Create a CalParams object to contain the results
Point centroid(posSum[0], posSum[1], posSum[2]);
Vector axis(axisSum[0], axisSum[1], axisSum[3]);
Event::CalParams params(rawEnergy, rawEneError, centroid, posWghtSum, axis, axisWghtSum);
// Create a CalCorToolResult object to hold the information
Event::CalCorToolResult* corResult = new Event::CalCorToolResult();
corResult->setStatusBit(Event::CalCorToolResult::RAWENERGY);
corResult->setParams(params);
corResult->setChiSquare(1.);
// Note that this tool **cannot** be called without an CalEventEnergy object already
// existing in the TDS. No checking done to insure this though...
SmartDataPtr<Event::CalEventEnergy> eventEnergy(m_dataSvc,EventModel::CalRecon::CalEventEnergy);
// Add the corrected result to the CalEventEnergy object
eventEnergy->push_back(corResult);
return sc;
}
|