153 {
154
155 MsgStream log(
msgSvc(), name() );
156 log << MSG::DEBUG << "in execute()" << endmsg;
157
158
159 int event, run;
160
161 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
162 if ( !eventHeader )
163 {
164 log << MSG::FATAL << name() << " Could not find Event Header" << endmsg;
165 return ( StatusCode::FAILURE );
166 }
167 run = eventHeader->runNumber();
168 event = eventHeader->eventNumber();
170
173
174 if ( fEventNb != 0 && m_event % fEventNb == 0 )
175 { log << MSG::FATAL << m_event << "-------: " << run << "," << event << endmsg; }
176 m_event++;
177 if ( fOutput >= 2 )
178 {
179 log << MSG::DEBUG << "====================================" << endmsg;
180 log << MSG::DEBUG << "run= " << run << "; event= " << event << endmsg;
181 }
182
183 Hep3Vector posG;
184#ifndef OnlineMode
185 if ( !( m_rawDataProviderSvc->isOnlineMode() ) )
186 {
187 if ( fOutput >= 1 && run < 0 )
188 {
189
190 SmartDataPtr<McParticleCol> mcParticleCol( eventSvc(), "/Event/MC/McParticleCol" );
191 if ( !mcParticleCol ) { log << MSG::WARNING << "Could not find McParticle" << endmsg; }
192 else
193 {
194 HepLorentzVector pG;
195 McParticleCol::iterator iter_mc = mcParticleCol->begin();
196 for ( ; iter_mc != mcParticleCol->end(); iter_mc++ )
197 {
198 log << MSG::INFO << " particleId = " << ( *iter_mc )->particleProperty() << endmsg;
199 pG = ( *iter_mc )->initialFourMomentum();
200 posG = ( *iter_mc )->initialPosition().v();
201 }
202 ttheta = pG.theta();
203 tphi = pG.phi();
204 if ( tphi < 0 ) { tphi = twopi + tphi; }
205 tp = pG.rho();
206
207
208 SmartDataPtr<EmcMcHitCol> emcMcHitCol( eventSvc(), "/Event/MC/EmcMcHitCol" );
209 if ( !emcMcHitCol ) { log << MSG::WARNING << "Could not find EMC truth" << endmsg; }
210
212 unsigned int mcTrackIndex;
213 double mcPosX = 0, mcPosY = 0, mcPosZ = 0;
214 double mcPx = 0, mcPy = 0, mcPz = 0;
215 double mcEnergy = 0;
216
217 EmcMcHitCol::iterator iterMc;
218 for ( iterMc = emcMcHitCol->begin(); iterMc != emcMcHitCol->end(); iterMc++ )
219 {
220 mcId = ( *iterMc )->identify();
221 mcTrackIndex = ( *iterMc )->getTrackIndex();
222 mcPosX = ( *iterMc )->getPositionX();
223 mcPosY = ( *iterMc )->getPositionY();
224 mcPosZ = ( *iterMc )->getPositionZ();
225 mcPx = ( *iterMc )->getPx();
226 mcPy = ( *iterMc )->getPy();
227 mcPz = ( *iterMc )->getPz();
228 mcEnergy = ( *iterMc )->getDepositEnergy();
229
230 if ( fOutput >= 2 )
231 {
232
233
234
235
236 }
237 }
238 }
239 }
240 }
241
242#endif
243
244
245 fDigitMap.clear();
246 fHitMap.clear();
247 fClusterMap.clear();
248 fShowerMap.clear();
249
250
251 IEmcCalibConstSvc* emcCalibConstSvc = 0;
252 StatusCode sc = service( "EmcCalibConstSvc", emcCalibConstSvc );
253 if ( sc != StatusCode::SUCCESS )
254 {
255 ;
256
257 }
258
259 SmartDataPtr<EmcDigiCol> emcDigiCol( eventSvc(), "/Event/Digi/EmcDigiCol" );
260 if ( !emcDigiCol )
261 {
262 log << MSG::FATAL << "Could not find EMC digi" << endmsg;
263 return ( StatusCode::FAILURE );
264 }
265
266 EmcDigiCol::iterator iter3;
267 for ( iter3 = emcDigiCol->begin(); iter3 != emcDigiCol->end(); iter3++ )
268 {
269 RecEmcID id( ( *iter3 )->identify() );
270 double adc =
273
274
278
279 int index = emcCalibConstSvc->
getIndex( nnpart, nnthe, nnphi );
281
282 if ( run > 0 && ixtalNumber > 0 )
283 {
284 unsigned int npartNew = emcCalibConstSvc->
getPartID( ixtalNumber );
285 unsigned int ntheNew = emcCalibConstSvc->
getThetaIndex( ixtalNumber );
286 unsigned int nphiNew = emcCalibConstSvc->
getPhiIndex( ixtalNumber );
288 }
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345 if ( run >= 52940 && run <= 52995 )
346 {
347 unsigned int ntheNew, nphiNew;
348 ntheNew = nnthe;
349 nphiNew = nnphi;
350 if ( nnpart == 2 )
351 {
352 if ( nnthe == 0 )
353 {
354 if ( nnphi == 58 ) nphiNew = 60;
355 if ( nnphi == 59 ) nphiNew = 61;
356 if ( nnphi == 60 ) nphiNew = 58;
357 if ( nnphi == 61 ) nphiNew = 59;
358 ntheNew = 0;
359 }
360 if ( nnthe == 1 )
361 {
362 if ( nnphi == 58 ) nphiNew = 60;
363 if ( nnphi == 59 ) nphiNew = 61;
364 if ( nnphi == 60 ) nphiNew = 58;
365 if ( nnphi == 61 ) nphiNew = 59;
366 ntheNew = 1;
367 }
368 if ( nnthe == 2 )
369 {
370 if ( nnphi == 73 ) nphiNew = 75;
371 if ( nnphi == 74 ) nphiNew = 76;
372 if ( nnphi == 75 ) nphiNew = 73;
373 if ( nnphi == 76 ) nphiNew = 74;
374 ntheNew = 2;
375
376 }
377 if ( nnthe == 3 )
378 {
379 if ( nnphi == 73 ) nphiNew = 75;
380 if ( nnphi == 74 ) nphiNew = 76;
381 if ( nnphi == 75 ) nphiNew = 73;
382 if ( nnphi == 76 ) nphiNew = 74;
383 ntheNew = 3;
384 }
385 if ( nnthe == 3 && nnphi == 72 )
386 {
387 ntheNew = 2;
388 nphiNew = 77;
389 }
390 if ( nnthe == 2 && nnphi == 77 )
391 {
392 ntheNew = 3;
393 nphiNew = 72;
394 }
395
396 if ( nnthe == 4 )
397 {
398 if ( nnphi == 87 ) nphiNew = 90;
399 if ( nnphi == 88 ) nphiNew = 91;
400 if ( nnphi == 89 ) nphiNew = 92;
401 if ( nnphi == 90 ) nphiNew = 87;
402 if ( nnphi == 91 ) nphiNew = 88;
403 if ( nnphi == 92 ) nphiNew = 89;
404 ntheNew = 4;
405 }
406 if ( nnthe == 5 )
407 {
408 if ( nnphi == 87 ) nphiNew = 90;
409 if ( nnphi == 88 ) nphiNew = 91;
410 if ( nnphi == 89 ) nphiNew = 92;
411 if ( nnphi == 90 ) nphiNew = 87;
412 if ( nnphi == 91 ) nphiNew = 88;
413 if ( nnphi == 92 ) nphiNew = 89;
414 ntheNew = 5;
415 }
416
418 }
419 }
420
421
422
423 if ( ixtalNumber == -9 ) { adc = 0.0; }
424
425
426 if ( ixtalNumber == -99 ) { adc = 0.0; }
427
428 RecEmcDigit aDigit;
429 aDigit.
Assign(
id, adc, tdc );
430 fDigitMap[id] = aDigit;
431 }
432
433 if ( fOutput >= 2 )
434 {
435 RecEmcDigitMap::iterator iDigitMap;
436 for ( iDigitMap = fDigitMap.begin(); iDigitMap != fDigitMap.end(); iDigitMap++ )
437 {
438
439 }
440 }
441
442
443
444
445 fDigit2Hit.Convert( fDigitMap, fHitMap );
446 if ( fOutput >= 2 )
447 {
448 RecEmcHitMap::iterator iHitMap;
449 for ( iHitMap = fHitMap.begin(); iHitMap != fHitMap.end(); iHitMap++ )
450 {
451
452 }
453 fDigit2Hit.Output( fHitMap );
454 }
455
456
457
458 fHit2Cluster.Convert( fHitMap, fClusterMap );
459
460 if ( fOutput >= 2 )
461 {
462 RecEmcClusterMap::iterator iClusterMap;
463 for ( iClusterMap = fClusterMap.begin(); iClusterMap != fClusterMap.end(); iClusterMap++ )
464 {
465
466 }
467 }
468
469
470 fCluster2Shower->Convert( fClusterMap, fShowerMap );
471
472 if ( fOutput >= 2 )
473 {
474 RecEmcShowerMap::iterator iShowerMap;
475 for ( iShowerMap = fShowerMap.begin(); iShowerMap != fShowerMap.end(); iShowerMap++ )
476 {
477
478 }
479 }
480
481
482
483 EmcRecTDS tds;
486
487
488#ifndef OnlineMode
489 if ( !( m_rawDataProviderSvc->isOnlineMode() ) )
490 {
491 if ( fOutput >= 1 )
492 {
493 nrun = run;
494 nrec = event;
495
496
497 ndigi = fDigitMap.size();
498 nhit = fHitMap.size();
499 ncluster = fClusterMap.size();
501 RecEmcShowerMap::iterator iShowerMap;
502 for ( iShowerMap = fShowerMap.begin(); iShowerMap != fShowerMap.end(); iShowerMap++ )
503 { fShowerVec.push_back( iShowerMap->second ); }
504 sort( fShowerVec.begin(), fShowerVec.end(), greater<RecEmcShower>() );
505 nneu = fShowerMap.size();
506
507 RecEmcShower aShower;
508
509 RecEmcShowerVec::iterator iShowerVec;
510 iShowerVec = fShowerVec.begin();
511 npart = -99;
512 ntheta = -99;
513 nphi = -99;
514 if ( iShowerVec != fShowerVec.end() )
515 {
516 aShower = *iShowerVec;
524 pp1[3] = aShower.
energy();
525 theta1 = ( aShower.
position() - posG ).theta();
526 phi1 = ( aShower.
position() - posG ).phi();
527 if ( phi1 < 0 ) { phi1 = twopi + phi1; }
529 phigap1 = aShower.
PhiGap();
533 eseed = aShower.
eSeed();
534 e3x3 = aShower.
e3x3();
535 e5x5 = aShower.
e5x5();
544 if ( nseed.
is_valid() ) { enseed = fHitMap[nseed].getEnergy(); }
545 else { enseed = 0; }
546
547 dphi1 = phi1 - tphi;
548 if ( dphi1 < -
pi ) dphi1 = dphi1 + twopi;
549 if ( dphi1 >
pi ) dphi1 = dphi1 - twopi;
550 }
551
552 if ( iShowerVec != fShowerVec.end() )
553 {
554 iShowerVec++;
555 if ( iShowerVec != fShowerVec.end() )
556 {
557 aShower = *iShowerVec;
561 pp2[3] = aShower.
energy();
562 }
563 }
564
565
566 if ( fShowerVec.size() >= 2 )
567 {
568 RecEmcShowerVec::iterator iShowerVec1, iShowerVec2;
569 iShowerVec1 = fShowerVec.begin();
570 iShowerVec2 = fShowerVec.begin() + 1;
571 double e1 = ( *iShowerVec1 ).energy();
572 double e2 = ( *iShowerVec2 ).energy();
573 double angle = ( *iShowerVec1 ).position().angle( ( *iShowerVec2 ).position() );
574 mpi0 = sqrt( 2 *
e1 *
e2 * ( 1 -
cos( angle ) ) );
575 }
576 m_tuple->write();
577 }
578 }
579#endif
580
581 return StatusCode::SUCCESS;
582}
vector< RecEmcShower > RecEmcShowerVec
double cos(const BesAngle a)
HepPoint3D position() const
double secondMoment() const
static Identifier crystal_id(const unsigned int barrel_ec, const unsigned int theta_module, const unsigned int phi_module)
For a single crystal.
static unsigned int barrel_ec(const Identifier &id)
Values of different levels (failure returns 0).
static unsigned int theta_module(const Identifier &id)
static unsigned int phi_module(const Identifier &id)
static EmcRecParameter & GetInstance()
void SetDataMode(double en)
StatusCode RegisterToTDS(RecEmcHitMap &aHitMap, RecEmcClusterMap &aClusterMap, RecEmcShowerMap &aShowerMap)
StatusCode CheckRegister()
virtual unsigned int getPartID(int Index) const =0
virtual int getIxtalNumber(int No) const =0
virtual unsigned int getPhiIndex(int Index) const =0
virtual unsigned int getThetaIndex(int Index) const =0
virtual int getIndex(unsigned int PartId, unsigned int ThetaIndex, unsigned int PhiIndex) const =0
bool is_valid() const
Check if id is in a valid state.
static double EmcTime(int timeChannel)
static double EmcCharge(int measure, int chargeChannel)
double getSecondMoment() const
void Assign(const RecEmcID &CellId, const RecEmcADC &ADC, const RecEmcTDC &TDC)
RecEmcID getShowerId() const
RecEmcCluster * getCluster() const
RecEmcEnergy getETof2x1() const
RecEmcID NearestSeed() const
RecEmcEnergy getETof2x3() const