Algorithm execution.
281 {
282 MsgStream log(
msgSvc(),name());
283
284 log<<MSG::DEBUG<< "in execute()" <<endreq;
285
286 int event, run;
287 ifpass = false;
288
289
290
291
292 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
293 if (!eventHeader) {
294 log << MSG::FATAL << "Could not find Event Header" << endreq;
295 return( StatusCode::FAILURE);
296 }
297 run = eventHeader->runNumber();
298 event = eventHeader->eventNumber();
299
300 if(mTrigRootFlag) {
301
302 m_RunId = run;
303 m_EventId = event;
304 }
305
306
307
308
309 if(writeFile==1 && event == 0)
310 {
311 readin.open(input.c_str(),ios_base::in);
312 if(readin) cout<<"Input File is ok "<<input<<endl;
313 readout.open(output.c_str(),ios_base::out);
314 if(readout) cout<<"Output File is ok "<<output<<endl;
316 readin >> version;
317 readout << version;
318 }
319
320
321
322
323 multimap<int,uint32_t,less<int> > mdc_hitmap;
324 mdc_hitmap.clear();
325
326
327
328
329 multimap<int,int,less<int> > tof_hitmap;
330 tof_hitmap.clear();
331
332
333
334
335 SmartDataPtr<MdcDigiCol> mdcDigiCol(eventSvc(),"/Event/Digi/MdcDigiCol");
336 if (!mdcDigiCol) {
337 log << MSG::FATAL << "Could not find MDC digi" << endreq;
338 return( StatusCode::FAILURE);
339 }
340 for(MdcDigiCol::iterator iter3=mdcDigiCol->begin();iter3!=mdcDigiCol->end();iter3++)
341 {
345 int tdc = (*iter3)->getTimeChannel();
346 if(tdc < 0x7FFFFFFF && tdc > 0) {
347 if(layer<=19) {
348 typedef pair<int, uint32_t > vpair;
349 uint32_t mdc_Id = (layer & 0xFFFF ) << 16;
350 mdc_Id = mdc_Id | (wire & 0xFFFF);
351 mdc_hitmap.insert(vpair(tdc,mdc_Id));
352 }
353 if(layer>=36&&layer<=39)
354 {
355 typedef pair<int, uint32_t > vpair;
356 uint32_t mdc_Id = (layer & 0xFFFF ) << 16;
357 mdc_Id = mdc_Id | (wire & 0xFFFF);
358 mdc_hitmap.insert(vpair(tdc,mdc_Id));
359 }
360 }
361 }
362
363
364
365
367 for(TofDataMap::iterator iter3 = tofDigiMap.begin();iter3 != tofDigiMap.end(); iter3++) {
369 TofData * p_tofDigi = iter3->second;
373 if(barrel_ec == 1) {
374 if(((p_tofDigi->
quality()) & 0x5) != 0x5)
continue;
375 double tdc1 = p_tofDigi->
tdc1();
376 double tdc2 = p_tofDigi->
tdc2();
377 int tdc;
378 if(tdc1 > tdc2) tdc = (int) tdc1;
379 else tdc = (int) tdc2;
380 typedef pair<int, int > vpair;
381 tdc = tdc;
382 tof_hitmap.insert(vpair(tdc,10000*barrel_ec+1000*layer+im*10));
383 }
384 else {
385 int tdc1 = (int)p_tofDigi->
tdc1();
386 typedef pair<int, int > vpair;
387 tdc1 = tdc1;
388 tof_hitmap.insert(vpair(tdc1,10000*barrel_ec+1000*layer+im*10));
389 }
390 }
391
392
393
394
395 SmartDataPtr<EmcDigiCol> emcDigiCol(eventSvc(),"/Event/Digi/EmcDigiCol");
396 if (!emcDigiCol) {
397 log << MSG::FATAL << "Could not find EMC digi" << endreq;
398 return( StatusCode::FAILURE);
399 }
400
401
402
403
405
406
407
408
409 multimap<int,uint32_t,less<int> > emc_TC;
410 emc_TC.clear();
411
412
413
414
416
417
418
419
420 double peak_time = 0;
422
423
424
425
426 SmartDataPtr<MucDigiCol> mucDigiCol(eventSvc(),"/Event/Digi/MucDigiCol");
427 if (!mucDigiCol) {
428 log << MSG::FATAL << "Could not find MUC digi" << endreq;
429 return( StatusCode::FAILURE);
430 }
431 if(m_mucDigi) m_mucDigi->
getMucDigi(mucDigiCol);
432
433
434
435
436 if(event%10000 == 0) std::cout << "---> EventNo is " << event << std::endl;
437 nEvent++;
438
439
440
441
442
443
444
445
446 double stime = -1, etime = -1;
447 findSETime(mdc_hitmap,tof_hitmap,emc_TC,stime,etime);
448
449
450 int nclock = 0;
451 if(stime >= 0) {
452 nclock = int ((etime - stime)/24) + 1;
453 }
454 else {
455 nclock = 0;
456 }
457
458
459
460
461 int** trgcond = new int*[48];
462 for(int condId = 0; condId < 48; condId++) {
463 trgcond[condId] = new int[nclock];
464 }
465
466
467 int idle_status = -1;
468
469 for(int iclock = 0; iclock < nclock; iclock++) {
470
471
472
474
475
476
477
479
480
481
482
484
485
486
487
488
489
490
491
492
494 if(status!=StatusCode::SUCCESS) {
495 log<<MSG::FATAL<< "Could not set trigger condition index" <<endreq;
496 return StatusCode::FAILURE;
497 }
498
499
500
501
502 for(int condId = 0; condId < 48; condId++) {
503 trgcond[condId][iclock] = m_pIBGT->
getTrigCond(condId);
504 }
505 }
506
507
508
509
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532 bool ifClus1 = false;
533 for(int i = 0; i < nclock; i++) {
534 if(trgcond[0][i] > 0) ifClus1 = true;
535 }
536
537 if(ifClus1 == false) {
538 for(int i = 0; i < nclock; i++) {
539 for(int j = 0; j < 16; j++) {
540 trgcond[j][i] = 0;
541 }
542 }
543 }
544
545
546
547
548 for(int i = 0; i < nclock; i++) {
549 for(int j = 0; j < 48; j++) {
551 }
552 }
553
554
555
556
558
559
560
561
563 if(ifpass)
564 {
565 passNo++;
566 log<<MSG::INFO<<"pass event number is "<<passNo<<endl;
567 }
568
569
570
571
572 if(writeFile == 2) {
573 if(ifpass)
574 {
575 setFilterPassed(true);
576 }
577 else
578 {
579 setFilterPassed(false);
580 }
581 }
582
583 if(mTrigRootFlag) {
584
585
586
587 for(int i = 0; i < 48; i++) {
588 bool edge = false;
589 int NOne = 0;
590 m_condNOne[i] = -9;
591 m_condNZero[i] = -9;
592 for(int j = 0; j < nclock; j++) {
593 m_mc_cond[i] += trgcond[i][j];
594 if(trgcond[i][j] != 0) {
595 if (NOne == 0) {
596 m_condNZero[i] = j;
597 m_trigCondi_MC->Fill(i);
598
599 }
600 edge = true;
601 NOne++;
602 }
603 else {
604 edge = false;
605 }
606 if(edge == false && NOne != 0) break;
607 }
608 m_condNOne[i] = NOne;
609 }
610 m_cond_id = 48;
611
612
613
614
615 if(m_runMode == 0) {
616 SmartDataPtr<TrigData> trigData(eventSvc(), "/Event/Trig/TrigData");
617 if (!trigData) {
618 log << MSG::FATAL << "Could not find Trigger Data for physics analysis" << endreq;
619 return StatusCode::FAILURE;
620 }
621
622 for(int id = 0; id < 48; id++) {
623 if(trigData->getTrigCondition(id) != 0) { m_trigCondi_Data->Fill(id); }
624 m_data_cond[id] = trigData->getTrigCondition(id);
625 }
626 m_cond_id = 48;
627 }
628 }
629
630
631
632
633 for(int condId = 0; condId < 48; condId++) {
634 delete trgcond[condId];
635 }
636 delete trgcond;
637
638
639 if(mTrigRootFlag) {
640 m_evtId = event;
641 m_tuple3->write();
642
646
647 map<int,vector<complex<int> >, greater<int> > mymap;
649 log<<MSG::INFO<<"EMC test: "<<endreq;
650 int emc_btc_id = 0;
651 for(map<
int,vector<complex<int> >, greater<int> >::iterator
iter=mymap.begin();
iter!=mymap.end();
iter++) {
652 if((
iter->first)==1) {
653 for(
unsigned int i=0; i<(
iter->second).size(); i++) {
654 log<<MSG::INFO<<
"barrel theta is "<<(
iter->second[i]).
real()<<
" phi is "<<(
iter->second[i]).
imag()<<endreq;
655 emc_btc_id++;
656 }
657 }
658 if((
iter->first)==0) {
659 for(
unsigned int i=0; i<(
iter->second).size(); i++)
660 log<<MSG::INFO<<
"east theta is "<<(
iter->second[i]).real()<<
" phi is "<<(
iter->second[i]).
imag()<<endreq;
661 }
662 if((
iter->first)==2) {
663 for(
unsigned int i=0; i<(
iter->second).size(); i++)
664 log<<MSG::INFO<<
"west theta is "<<(
iter->second[i]).real()<<
" phi is "<<(
iter->second[i]).
imag()<<endreq;
665 }
666 }
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742 m_tuple1->write();
743
744
745
746
747 vector<int> vstrkId;
748 vector<int> vltrkId;
751 log<<MSG::INFO<<"Mdc test: "<<endreq;
752 for(unsigned int i=0; i<vstrkId.size(); i++) log<<MSG::INFO<<"short is "<<vstrkId[i]<<endreq;
753 for(unsigned int j=0; j<vltrkId.size(); j++) { log<<MSG::INFO<<"long is "<<vltrkId[j]<<endreq; }
754
755 map<int,vector<int>,greater<int> > tofmap;
757 log<<MSG::INFO<<"TOF test: "<<endreq;
758 for(map<
int,vector<int>,greater<int> >::iterator
iter=tofmap.begin();
iter!=tofmap.end();
iter++) {
759 if(
iter->first == 0) {
760 for(
unsigned int i=0; i<
iter->second.size(); i++) {
761 log<<MSG::INFO<<
"east tof Id is "<<
iter->second[i]<<endreq;
762 }
763 }
764 if(
iter->first == 1) {
765 for(
unsigned int i=0; i<
iter->second.size(); i++) { log<<MSG::INFO<<
" barrel tof Id is "<<
iter->second[i]<<endreq; }
766 }
767 if(
iter->first == 2) {
768 for(
unsigned int i=0; i<
iter->second.size(); i++) { log<<MSG::INFO<<
"west tof Id is "<<
iter->second[i]<<endreq; }
769 }
770 }
771
772
773 std::vector<int> vtmp;
774
776 m_index2 = 0;
777 for(std::vector<int>::iterator
iter = vtmp.begin();
iter != vtmp.end();
iter++) {
778 m_fireLayer[m_index2] = (long) *
iter;
779 m_index2++;
780 if(m_index2 > m_index2->range().distance()) { break; cerr<<"*********** too many index ************"<<endl; }
781 }
782
783 long trackb3=0, tracke3=0, trackb2=0, tracke2=0, trackb1=0, tracke1=0;
784 int trackwe = 0, trackee = 0;
785 for(unsigned int i=0; i<vtmp.size(); i++) {
786 if(0<=vtmp[i]&&vtmp[i]<100) {
787 if((vtmp[i]%10)>=3) { tracke3++; trackee++; }
788 }
789 if(200<=vtmp[i]) {
790 if(((vtmp[i]-200)%10)>=3) { tracke3++; trackwe++; }
791 }
792 if(100<=vtmp[i]&&vtmp[i]<200) {
793 if(((vtmp[i]-100)%10)>=3) trackb3++;
794 }
795 }
796 m_ntrack3 = trackb3 + tracke3;
797
798 for(unsigned int i=0; i<vtmp.size(); i++) {
799 if(0<=vtmp[i]&&vtmp[i]<100) {
800 if((vtmp[i]%10)>=2) tracke2++;
801 }
802 if(200<=vtmp[i]) {
803 if(((vtmp[i]-200)%10)>=2) tracke2++;
804 }
805 if(100<=vtmp[i]&&vtmp[i]<200) {
806 if(((vtmp[i]-100)%10)>=2) trackb2++;
807 }
808 }
809 m_ntrack2 = trackb2 + tracke2;
810
811 for(unsigned int i=0; i<vtmp.size(); i++) {
812 if(0<=vtmp[i]&&vtmp[i]<100) {
813 if((vtmp[i]%10)>=1) tracke1++;
814 }
815 if(200<=vtmp[i]) {
816 if(((vtmp[i]-200)%10)>=1) tracke1++;
817 }
818 if(100<=vtmp[i]&&vtmp[i]<200) {
819 if(((vtmp[i]-100)%10)>=1) trackb1++;
820 }
821 }
822 m_ntrack1 = trackb1 + tracke1;
823
824
826 m_index3 = 0;
827 for(std::vector<int>::iterator
iter = vtmp.begin();
iter != vtmp.end();
iter++) {
828 m_hitLayer[m_index3] = (long) *
iter;
829 m_index3++;
830 if(m_index3 > m_index3->range().distance()) { break; cerr<<"*********** too many index ************"<<endl; }
831 }
832
834 m_index4 = 0;
835 for(std::vector<int>::iterator
iter = vtmp.begin();
iter != vtmp.end();
iter++) {
836 m_hitSeg[m_index4] = *(
iter);
837 m_index4++;
838 if(m_index4 > m_index4->range().distance()) { break; cerr<<"*********** too many index ************"<<endl; }
839 }
840 }
841
842
843
844
845 if(ifoutEvtId==1)
846 {
847 ofstream eventnum(outEvtId.c_str(),ios_base::app);
848 if(!ifpass)
849 eventnum<<event<<endl;
850 eventnum.close();
851 }
852
853
854
855
856 if(ifoutEvtId==2)
857 {
858 ofstream eventnum(outEvtId.c_str(),ios_base::app);
859 if(ifpass)
860 eventnum<<event<<endl;
861 eventnum.close();
862 }
863
864
865
866
867 if(writeFile==1)
868 {
870 readin >> asciiEvt;
872 {
873 if(ifpass==true)
874 readout<<asciiEvt<<endl;
875 }
876 else
877 cout<<"********* Event No. from AsciiFile do not equal Event No. from TDS "
879 }
880
881
882
883
884 if(m_runMode == 1) {
887 int window = 0;
888 int timing = 0;
889 bool preScale = false;
890
891 StatusCode sc = StatusCode::SUCCESS ;
893 sc = eventSvc()->registerObject("/Event/Trig",aTrigEvent);
894 if(sc!=StatusCode::SUCCESS) {
895 log<<MSG::DEBUG<< "Could not register TrigEvent, you can ignore." <<endreq;
896 }
897
898 TrigData* aTrigData =
new TrigData(window, timing, trigcond, trigchan, preScale);
899 sc = eventSvc()->registerObject("/Event/Trig/TrigData",aTrigData);
900 if(sc!=StatusCode::SUCCESS) {
901 log<<MSG::ERROR<< "Could not register TrigData!!!!!" <<endreq;
902 }
903 }
904
905 return StatusCode::SUCCESS;
906}
std::multimap< unsigned int, TofData * > TofDataMap
double imag(const EvtComplex &c)
double real(const EvtComplex &c)
void setTrigCond(int i, bool j)
map< int, vector< complex< int > >, greater< int > > getEmcClusId()
std::vector< int > getMdcLtrkId()
const int getTrigCond(int i)
std::vector< int > getMdcStrkId()
std::vector< int > getMuchitLayer()
std::vector< int > getMuchitSeg()
const int getTrigChan(int i)
map< int, vector< int >, greater< int > > getTofHitPos()
StatusCode setTrigCondition()
std::vector< int > getMuclayerSeg()
void runAclock_emc(int iclock, double stime, std::multimap< int, uint32_t, less< int > > emc_TC, EmcWaveform *blockWave)
void findEmcPeakTime(double &peak_time, EmcWaveform *blockWave)
void runAclock_tof(int iclock, double stime, int &idle_status, std::multimap< int, int, less< int > > tof_hitmap)
void getEmcAnalogSig(EmcDigiCol *emcDigiCol, EmcWaveform(&blockWave)[16], multimap< int, uint32_t, less< int > > &emc_TC)
void findSETime(multimap< int, uint32_t, less< int > > mdc_hitmap, multimap< int, int, less< int > > tof_hitmap, multimap< int, uint32_t, less< int > > emc_TC, double &stime, double &etime)
void stretchTrgCond(int nclock, int **&trgcond)
void runAclock_mdc(int iclock, double stime, multimap< int, uint32_t, less< int > > mdc_hitmap)
virtual TofDataMap & tofDataMapEstime()=0
static int layer(const Identifier &id)
Values of different levels (failure returns 0)
static int wire(const Identifier &id)
void getMucDigi(MucDigiCol *mucDigiCol)
unsigned int quality() const
static int phi_module(const Identifier &id)
static int barrel_ec(const Identifier &id)
Values of different levels (failure returns 0)
static int layer(const Identifier &id)