00001 #include <vector>
00002
00003 #include "GaudiKernel/MsgStream.h"
00004 #include "GaudiKernel/SmartDataPtr.h"
00005 #include "GaudiKernel/PropertyMgr.h"
00006 #include "GaudiKernel/Bootstrap.h"
00007
00008 #include "GaudiKernel/INTupleSvc.h"
00009 #include "GaudiKernel/NTuple.h"
00010 #include "GaudiKernel/ITHistSvc.h"
00011
00012 #include "CLHEP/Vector/ThreeVector.h"
00013 #include "CLHEP/Vector/LorentzVector.h"
00014
00015 #include "TH1F.h"
00016
00017 #include "EventModel/EventModel.h"
00018 #include "EventModel/Event.h"
00019
00020 #include "EvtRecEvent/EvtRecEvent.h"
00021 #include "EvtRecEvent/EvtRecTrack.h"
00022 #include "DstEvent/TofHitStatus.h"
00023 #include "EventModel/EventHeader.h"
00024
00025
00026 #include "ParticleID/ParticleID.h"
00027
00028 #include "DQAEvent/DQAEvent.h"
00029 #include "DQAFillEx/DQAFillEx.h"
00030
00031 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00032 typedef HepGeom::Point3D<double> HepPoint3D;
00033 #endif
00034 using CLHEP::HepLorentzVector;
00035
00036
00038
00039 DQAFillEx::DQAFillEx(const std::string& name, ISvcLocator* pSvcLocator) :
00040 Algorithm(name, pSvcLocator) {
00041
00042
00043 }
00044
00045
00046 StatusCode DQAFillEx::initialize(){
00047 MsgStream log(msgSvc(), name());
00048
00049 log << MSG::INFO << "in initialize()" << endmsg;
00050 StatusCode status;
00051
00052
00053
00054
00055
00056 NTuplePtr nt(ntupleSvc(), "DQAFILE/MDC");
00057 if ( nt ) m_tuple = nt;
00058 else {
00059 m_tuple = ntupleSvc()->book("DQAFILE/MDC", CLID_ColumnWiseTuple, "MDC ntuple");
00060 if( m_tuple ) {
00061 status = m_tuple->addItem("runNo", m_runNo);
00062 status = m_tuple->addItem("event", m_event);
00063 } else {
00064 log << MSG::ERROR << "Can not book N-tuple:" << long(m_tuple) << endreq;
00065 }
00066 }
00067
00068 if(service("THistSvc", m_thsvc).isFailure()) {
00069 log << MSG::ERROR << "Couldn't get THistSvc" << endreq;
00070 return StatusCode::FAILURE;
00071 }
00072
00073
00074
00075 TH1F* hrxy = new TH1F("Rxy", "Rxy distribution", 110, -1., 10.);
00076 if(m_thsvc->regHist("/DQAHist/MDC/hrxy", hrxy).isFailure()) {
00077 log << MSG::ERROR << "Couldn't register Rxy" << endreq;
00078 }
00079 TH1F* hz = new TH1F("z", "z distribution", 200, -100., 100.);
00080 if(m_thsvc->regHist("/DQAHist/MDC/hz", hz).isFailure()) {
00081 log << MSG::ERROR << "Couldn't register z" << endreq;
00082 }
00083
00084 log << MSG::INFO << "successfully return from initialize()" <<endmsg;
00085 return StatusCode::SUCCESS;
00086
00087
00088 }
00089
00090
00091 StatusCode DQAFillEx::execute() {
00092
00093 MsgStream log(msgSvc(), name());
00094 log << MSG::INFO << "in execute()" << endreq;
00095
00096 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00097 m_runNo=eventHeader->runNumber();
00098 m_event=eventHeader->eventNumber();
00099 log << MSG::DEBUG <<"run, evtnum = "
00100 << m_runNo << " , "
00101 << m_event <<endreq;
00102
00103
00104 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00105
00106 SmartDataPtr<DQAEvent::DQAEvent> dqaevt(eventSvc(), "/Event/DQATag");
00107 if ( dqaevt ) {
00108 log << MSG::INFO << "success get DQAEvent" << endreq;
00109 } else {
00110 log << MSG::ERROR << "Error accessing DQAEvent" << endreq;
00111 return StatusCode::FAILURE;
00112 }
00113
00114 log << MSG::DEBUG << "event tag = " << dqaevt->EventTag() << endreq;
00115
00116
00117 if ( dqaevt->Bhabha() ) {
00118 log << MSG::INFO << "Bhabha event" << endreq;
00119 SmartDataPtr<EvtRecTrackCol> evtRecTrkCol(eventSvc(), EventModel::EvtRec::EvtRecTrackCol);
00120 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
00121 EvtRecTrackIterator itTrk=evtRecTrkCol->begin() + i;
00122 log << MSG::DEBUG << i << " " << (*itTrk)->partId() << " "
00123 << (*itTrk)->quality() << endreq;
00124
00125
00126
00127
00128 if ( !(*itTrk)->isElectron() ) continue;
00129
00130
00131
00132
00133
00134
00135 int qual = (*itTrk)->quality();
00136 if ( qual != 0 && qual != 2) continue;
00137
00138 RecMdcTrack *mdcTrk = (*itTrk)->mdcTrack();
00139 if ( mdcTrk->charge() > 0 ) {
00140 log << MSG::DEBUG << "is electron" << endreq;
00141 } else {
00142 log << MSG::DEBUG << "is positron" << endreq;
00143 }
00144 double x0 =mdcTrk->x();
00145 double y0 =mdcTrk->y();
00146 double z0 =mdcTrk->z();
00147 double Rxy=sqrt(x0*x0+y0*y0);
00148
00149 TH1 *h(0);
00150 if (m_thsvc->getHist("/DQAHist/MDC/hrxy", h).isSuccess()) {
00151 h->Fill(Rxy);
00152 } else {
00153 log << MSG::ERROR << "Couldn't retrieve hrxy" << endreq;
00154 }
00155 if (m_thsvc->getHist("/DQAHist/MDC/hz", h).isSuccess()) {
00156 h->Fill(z0);
00157 } else {
00158 log << MSG::ERROR << "Couldn't retrieve hz" << endreq;
00159 }
00160 }
00161 }
00162
00163 m_tuple->write();
00164
00165 return StatusCode::SUCCESS;
00166
00167 }
00168
00169
00170 StatusCode DQAFillEx::finalize() {
00171
00172 MsgStream log(msgSvc(), name());
00173 log << MSG::INFO << "in finalize()" << endmsg;
00174 return StatusCode::SUCCESS;
00175 }
00176
00177