00001 #include "ZddReconAlg/ZddReconAlg.h"
00002 #include "ReconEvent/ReconEvent.h"
00003 #include "ZddEvent/ZddBoard.h"
00004 #include "ZddEvent/RecZddChannel.h"
00005 #include "EvTimeEvent/RecEsTime.h"
00006 #include "EventModel/EventHeader.h"
00007 #include "GaudiKernel/MsgStream.h"
00008 #include "GaudiKernel/SmartDataPtr.h"
00009
00010 using namespace std;
00011
00013 ZddReconAlg::ZddReconAlg(const std::string& name, ISvcLocator* pSvcLocator) :
00014 Algorithm(name, pSvcLocator),
00015 m_errStat(false)
00016 {
00017 declareProperty("BreakWithError", m_errQuit = true);
00018 }
00019
00020
00021 StatusCode ZddReconAlg::initialize()
00022 {
00023 MsgStream log(msgSvc(), name());
00024 log << MSG::INFO << "in initialize()" << endreq;
00025
00026 return StatusCode::SUCCESS;
00027 }
00028
00029
00030 StatusCode ZddReconAlg::execute()
00031 {
00032 MsgStream log(msgSvc(), name());
00033 log << MSG::INFO << "in execute()" << endreq;
00034
00035
00036 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00037 if (!eventHeader) {
00038 log << MSG::FATAL << "Could not find Event Header" << endreq;
00039 return StatusCode::FAILURE;
00040 }
00041 log << MSG::DEBUG << "Retrieved event: " << eventHeader->eventNumber()
00042 << " run: " << eventHeader->runNumber() << endreq;
00043
00044
00045 DataObject* aRecEvent = 0;
00046 eventSvc()->findObject("/Event/Recon", aRecEvent);
00047 if ( aRecEvent == 0 ) {
00048 aRecEvent = new ReconEvent();
00049 StatusCode sc = eventSvc()->registerObject("/Event/Recon", aRecEvent);
00050 if ( sc.isFailure() ) {
00051 log << MSG::FATAL << "Could not register ReconEvent" << endreq;
00052 return StatusCode::FAILURE;
00053 }
00054 }
00055
00056 RecZddChannelCol* recZddCol = new RecZddChannelCol;
00057 StatusCode sc = eventSvc()->registerObject(EventModel::Recon::RecZddChannelCol, recZddCol);
00058 if ( sc.isFailure() ) {
00059 log << MSG::FATAL << "Could not register RecZddChannelCol" << endreq;
00060 return StatusCode::FAILURE;
00061 }
00062
00063
00064
00065 if ( m_errStat ) return StatusCode::SUCCESS;
00066
00067
00068 SmartDataPtr<Event::ZddEvent> zddEvt(eventSvc(),"/Event/ZddEvent");
00069 int zddCheck = zddDataStat(zddEvt.ptr(), eventHeader->eventNumber()+1);
00070
00071 if ( zddCheck != 0 ) {
00072 if ( zddCheck < 0 ) {
00073
00074 m_errStat = true;
00075
00076 if ( m_errQuit ) {
00077 return StatusCode::FAILURE;
00078 }
00079 }
00080
00081
00082 return StatusCode::SUCCESS;
00083 }
00084
00085
00086 double bes3_t0 = -10000.0;
00087 SmartDataPtr<RecEsTimeCol> evTimeCol(eventSvc(), "/Event/Recon/RecEsTimeCol");
00088 if ( !evTimeCol || evTimeCol->size() == 0 ) {
00089 log << MSG::WARNING << " Could not find RecEsTimeCol" << endreq;
00090
00091 }
00092 else {
00093 RecEsTimeCol::iterator iter_evt = evTimeCol->begin();
00094 if (iter_evt != evTimeCol->end()) {
00095 bes3_t0 = (*iter_evt)->getTest();
00096
00097 bes3_t0 = 6400 - bes3_t0;
00098 }
00099 }
00100
00101
00102 const Event::ZddEvent::Channels& chs = zddEvt->channels();
00103
00104
00105 Event::ZddEvent::Channels::const_iterator end_ch = chs.end();
00106 for ( Event::ZddEvent::Channels::const_iterator it = chs.begin(); it != end_ch; ++it ) {
00107 const ZddChannel* pch = *it;
00108 const ZddChannel::Fragments& frags = pch->fragments();
00109
00110
00111 RecZddChannel* recZddCh = new RecZddChannel;
00112 recZddCh->setChannelId(pch->getChId());
00113 recZddCh->setScanCode(pch->getScanCode());
00114 double e_K = getEK(pch->getChId());
00115
00116
00117 int maxSamples = 800;
00118 unsigned char waveform[800];
00119 memset(waveform, 0, maxSamples);
00120 ZddChannel::Fragments::const_iterator end_fg = frags.end();
00121 bool quit_event = false;
00122 int start = 0;
00123 for ( ZddChannel::Fragments::const_iterator jt = frags.begin(); jt != end_fg; ++jt) {
00124 const ZddFragment& frag = *jt;
00125
00126
00127
00128
00129
00130 start = frag.start_index;
00131
00132
00133 if ( start+frag.length > maxSamples ) {
00134 recZddCh->setScanCode( 2*10 + pch->getScanCode() );
00135 MsgStream log(msgSvc(), name());
00136 log << MSG::ERROR << "ZDD BAD DATA: CAEN corruption problem" << endreq;
00137 quit_event = true;
00138 break;
00139 }
00140
00141 for ( int i = 0; i < frag.length; ++i ) waveform[start++] = frag.sample[i];
00142 }
00143
00144
00145 unsigned char threshold = 20;
00146 unsigned char rephaseThreshold = 40;
00147 unsigned char minSample = 255, maxSample = -1;
00148 int maxTime = -1;
00149 bool closed = true;
00150 int phases[4] = {-1,-1,-1,-1};
00151 for ( int pt=0; pt<maxSamples; pt++ ) {
00152 bool notZero = waveform[pt]>0;
00153 bool smaller = waveform[pt] < minSample;
00154 if ( notZero && smaller ) minSample = waveform[pt];
00155 if ( waveform[pt] > threshold ) {
00156 if ( closed ) {
00157 closed = false;
00158 maxSample = waveform[pt];
00159 maxTime = pt;
00160 } else {
00161 if ( waveform[pt] > maxSample ) {
00162 maxSample = waveform[pt];
00163 maxTime = pt;
00164 }
00165 }
00166 } else {
00167 if ( ! closed ) {
00168 closed = true;
00169 double tNsec = 2.*maxTime;
00170 recZddCh->addFragment(tNsec, maxSample*e_K);
00171
00172 if ( maxSample > rephaseThreshold ) {
00173 int phase = maxTime%4;
00174 phases[phase]++;
00175 }
00176 }
00177 }
00178 }
00179 if ( ! closed ) {
00180 closed = true;
00181 double tNsec = 2.*maxTime;
00182 recZddCh->addFragment(tNsec, maxSample*e_K);
00183
00184 if ( maxSample > rephaseThreshold ) {
00185 int phase = maxTime%4;
00186 phases[phase]++;
00187 }
00188 }
00189
00190
00191 int mostProb = -1;
00192 int chPhase = -1;
00193 for (int ph=0; ph<4; ph++) {
00194 if ( phases[ph] > mostProb ) {
00195 mostProb = phases[ph];
00196 chPhase = ph;
00197 }
00198 }
00199
00200 if ( chPhase==-1 ) {
00201 recZddCh->setScanCode( 3*10 + pch->getScanCode() );
00202 chPhase = -2;
00203 }
00204
00205 recZddCh->setBaseLine(minSample);
00206 recZddCh->setPhase(chPhase);
00207 recZddCol->push_back(recZddCh);
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220 if ( quit_event ) break;
00221
00222 }
00223 return StatusCode::SUCCESS;
00224 }
00225
00226
00227 StatusCode ZddReconAlg::finalize()
00228 {
00229 MsgStream log(msgSvc(), name());
00230 log << MSG::INFO << "in finalize()" << endreq;
00231
00232 return StatusCode::SUCCESS;
00233 }
00234
00235 double ZddReconAlg::getEK(int chId)
00236 {
00237 return 1.0/16.0;
00238 }
00239
00240 int ZddReconAlg::zddDataStat(const Event::ZddEvent* zddEvt, int evtId)
00241 {
00242 MsgStream log(msgSvc(), name());
00243
00244 if ( !zddEvt ) {
00245 log << MSG::FATAL << "Could not find ZddEvent" << endreq;
00246 return 1;
00247 }
00248
00249
00250 const std::vector<ZddBoard*>& boards = zddEvt->boards();
00251 if ( boards.size() != 2 ) {
00252 log << MSG::FATAL << "incomplete ZDD data, no more ZDD data this run!" << endreq;
00253 return -2;
00254 }
00255
00256 if ( boards[0]->getCounter() != evtId || boards[1]->getCounter() != evtId )
00257 {
00258 log << MSG::FATAL << "missaligned ZDD triggers, no more ZDD data this run!" << endreq;
00259 return -3;
00260 }
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280 return 0;
00281 }
00282
00283
00284 int ZddReconAlg::calFragEandT(const ZddFragment& frag, int& efrag, int& tfrag)
00285 {
00286 int minIndex = 0;
00287 unsigned char min = 255, max = 0;
00288 for ( int i = 0; i < frag.length; ++i ) {
00289 unsigned char& sample = frag.sample[i];
00290 if ( sample < min ) {
00291 min = sample;
00292 minIndex = i;
00293 }
00294 if ( sample > max ) {
00295 max = sample;
00296 }
00297 }
00298
00299 efrag = max - min;
00300 tfrag = 8160 - 2 * (minIndex + frag.start_index);
00301
00302 return 0;
00303 }