00001
00002
00003
00004
00005
00006 #include "TruthExamples/TruthDemo.h"
00007 #include "GeneratorObject/McEventCollection.h"
00008
00009 #include "GaudiKernel/MsgStream.h"
00010 #include "GaudiKernel/AlgFactory.h"
00011
00012 #include "GaudiKernel/SmartDataPtr.h"
00013 #include "GaudiKernel/DataSvc.h"
00014
00015 #include "GaudiKernel/IHistogramSvc.h"
00016
00017 #include "TruthHelper/IsGenStable.h"
00018
00019 #include "AIDA/IHistogram1D.h"
00020 #include "AIDA/IHistogram2D.h"
00021
00022 #include "GaudiKernel/PropertyMgr.h"
00023 #include "GaudiKernel/INTupleSvc.h"
00024 #include "GaudiKernel/ISvcLocator.h"
00025 #include "GaudiKernel/IDataProviderSvc.h"
00026 #include "GaudiKernel/MsgStream.h"
00027 #include "GaudiKernel/ObjectList.h"
00028
00029 #include "HepMC/GenEvent.h"
00030 #include "HepMC/GenParticle.h"
00031
00032
00033
00034 #include "HepMC/GenVertex.h"
00035
00036 typedef std::vector<const HepMC::GenParticle*> MCparticleCollection ;
00037
00038
00039 static const AlgFactory<TruthDemo> Factory;
00040 const IAlgFactory& TruthDemoFactory = Factory;
00041
00042 TruthDemo::TruthDemo(const std::string& name, ISvcLocator* pSvcLocator) :
00043 Algorithm(name, pSvcLocator)
00044 {
00045
00046 declareProperty("HistogramFlag", m_produceHistogram = true );
00047
00048 }
00049
00050 StatusCode TruthDemo::initialize(){
00051 StatusCode result = StatusCode::SUCCESS;
00052 MsgStream msglog(messageService(), name());
00053 msglog << MSG::INFO << ">>> Truthdemo from Initialize" << endreq;
00054 m_hgenerated = histoSvc()->book("/stat/1Dhist/1","Generated",100,0,1200);
00055 if (0 == m_hgenerated) {
00056 msglog << MSG::ERROR << " ERROR booking histogram" << endreq;
00057 result = StatusCode::FAILURE;
00058 }
00059 m_pxBalance = histoSvc()->book("/stat/1Dhist/25","px balance",50,-10.,10.);
00060 m_pyBalance = histoSvc()->book("/stat/1Dhist/26","py balance",50,-10.,10.);
00061 m_totEnergy = histoSvc()->book("/stat/1Dhist/27","total energy",50,10000.,20000.);
00062 m_tesIO = new GenAccessIO();
00063
00064
00065 return result;
00066 }
00067 StatusCode TruthDemo::execute() {
00068
00069
00070
00071 MsgStream msglog(messageService(), name());
00072 msglog << MSG::INFO << ">>> TruthDemo from execute" << endreq;
00073
00074
00075
00076 float totenergy = 0.;
00077 float pxbalance = 0.;
00078 float pybalance = 0.;
00079
00080 IsGenStable ifs;
00081 std::vector<const HepMC::GenParticle*> particles;
00082 StatusCode stat = m_tesIO->getMC(particles, &ifs);
00083 for (std::vector<const HepMC::GenParticle*>::iterator pitr = particles.begin();
00084 pitr != particles.end(); pitr++) {
00085 pxbalance += (*pitr)->momentum().x();
00086 pybalance += (*pitr)->momentum().y();
00087 totenergy += (*pitr)->momentum().e();
00088 }
00089 m_pxBalance->fill(pxbalance, 1.);
00090 m_pyBalance->fill(pybalance, 1.);
00091 m_totEnergy->fill(totenergy, 1.);
00092
00093 return StatusCode::SUCCESS;
00094 }
00095
00096 StatusCode TruthDemo::finalize() {
00097 MsgStream msglog(messageService(), name());
00098 msglog << MSG::INFO << ">>> TruthDemo from finalize" << endreq;
00099 return StatusCode::SUCCESS;
00100 }
00101
00102