CalRawEnergyTool.cxx
/** @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; }