340 {
341 MsgStream log(
msgSvc(), name());
342 log << MSG::INFO << "in execute()" << endreq;
343
344 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
345 int runNo=eventHeader->runNumber();
346
347
350 if (runNum>=runFrom&&runNum<=runTo) {
351 m_inFlag=true;
352 } else {
353 m_inFlag=false;
354 }
355
356 if (m_inFlag==false){
357
358 runFrom=m_EmcShEnCalibSvc -> getRunFrom();
359 runTo=m_EmcShEnCalibSvc -> getRunTo();
360
361
362
363 cout <<
"current run=" <<
runNo<<endl;
364 cout << "in AbsCor open getPi0CalibFile()= " << m_EmcShEnCalibSvc -> getPi0CalibFile()<<endl;
365 cout << "open getSingleGammaCalibFile()= " << m_EmcShEnCalibSvc -> getSingleGammaCalibFile()<<endl;
366 cout <<"getRunFrom()="<< runFrom<<endl;
367 cout <<"getRunTo()="<< runTo<<endl;
368
369
370 string CorFunparaPath;
371 if(MCuseTof==1){
372 if (m_ifReadDB){
373 CorFunparaPath = m_EmcShEnCalibSvc -> getSingleGammaCalibFile();
374 } else {
375 CorFunparaPath = m_CorFunparaPath;
376 cout<<"no read database, but using the set:"<<CorFunparaPath<<endl;
377 }
378 }
379 if(MCuseTof==0){
380 CorFunparaPath = getenv("ABSCORROOT");
381 CorFunparaPath += "/share/evsetCorFunctionPar.txt";
382 }
383
384 ifstream in2corfun;
385 in2corfun.open(CorFunparaPath.c_str());
386 assert(in2corfun);
387
388 for(int i=0;i<28;i++){
389 in2corfun>>m_corFunPar[i][0];
390 in2corfun>>m_corFunPar[i][1];
391 in2corfun>>m_corFunPar[i][2];
392 in2corfun>>m_corFunPar[i][3];
393 in2corfun>>m_corFunPar[i][4];
394 in2corfun>>m_corFunPar[i][5];
395 }
396 in2corfun.close();
397
398
399 string DataPathc3p;
400 DataPathc3p = getenv("ABSCORROOT");
401 DataPathc3p += "/share/c3p.txt";
402
403 string DataPathc3ptof;
404 if (m_ifReadDB){
405 DataPathc3ptof = m_EmcShEnCalibSvc -> getPi0CalibFile();
406 } else {
407 DataPathc3ptof = m_DataPathc3ptof;
408 cout<<"no read database, but using the set:"<<DataPathc3ptof<<endl;
409 }
410
411 ifstream inc3p;
412 inc3p.exceptions( ifstream::failbit | ifstream::badbit );
413 if(usetof==0)inc3p.open(DataPathc3p.c_str(),ios::in);
414 if(usetof==1)inc3p.open(DataPathc3ptof.c_str(),ios::in);
415 for(int i=0; i<4; i++){
416 double am,ame;
417 inc3p>>am;
418 inc3p>>ame;
420 }
421 }
422
423
426 if(evtRecEvent->totalTracks()>evtRecTrkCol->size())return SUCCESS;
427 if(evtRecEvent->totalTracks()>50)return SUCCESS;
428
430 ISvcLocator* svcLocator = Gaudi::svcLocator();
431 StatusCode sc = svcLocator->service("EmcRecGeoSvc",iGeoSvc);
432 if(sc!=StatusCode::SUCCESS) {
433 cout<<"Error: Can't get EmcRecGeoSvc"<<endl;
434 }
435
436 for(int i = 0; i< evtRecEvent->totalTracks(); i++) {
438 if(!(*itTrk)->isEmcShowerValid())continue;
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
475
476 unsigned int module, ntheta, nphi;
477 module = EmcID::barrel_ec(id);
480
481
482
485
486
487 double e5x5=emcTrk->
e5x5();
488 double etof=0;
489
490 if(usetof==1 && (*itTrk)->isTofTrackValid()){
491 SmartRefVector<RecTofTrack> recTofTrackVec=(*itTrk)->tofTrack();
492 if(!recTofTrackVec.empty())etof=recTofTrackVec[0]->energy();
493 if(etof>100.)etof=0;
494 }
495
496 double energyC;
497
498 int thetaId;
499 if (module==0||module==2) thetaId = thetaModule;
500 if (module==1 && thetaModule<=21) thetaId = thetaModule + 6;
501 if (module==1 && thetaModule>21) thetaId = 43 - thetaModule + 6;
502 double DthetaId=double(thetaId);
503
504 if(MCuseTof==1){
505 if (thetaId<6) etof=0.0;
506 if (MCCorUseFunction==1){
507 energyC=ECorrFunctionMC(e5x5+etof,DthetaId);
508 } else if (MCCorUseFunction==0){
509 energyC=ECorrMC(e5x5+etof,DthetaId);
510 }
511 }
512 if (MCuseTof==0) {
513 if (MCCorUseFunction==1){
514 energyC=ECorrFunctionMC(e5x5,DthetaId);
515 } else if (MCCorUseFunction==0){
516 energyC=ECorrMC(e5x5,DthetaId);
517 }
518 }
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540 double lnE = std::log(energyC);
541
542 if (energyC>1.0) lnE=std::log(1.0);
543 if (energyC<0.05) lnE=std::log(0.05);
544
545 double lnEcor=1.0;
546 if(dopi0Cor==1){
547 if(
runNo>0 && dodatacor==1) {
548 lnEcor =
ai[0]+
ai[1]*lnE+
ai[2]*lnE*lnE+
ai[3]*lnE*lnE*lnE;
549 }
550 }
551 if(lnEcor<0.5)continue;
552
553
554
555
557
558
559 double enecor=1.;
560 if(hotcellmask==1){
561 for(int ih=0;ih<10;ih++){
562 if(hrunstart[ih]==-1||hrunend[ih]==-1||hcell[ih]==-1)continue;
565 }
566 }
567
568 if(edgecor==1){
569
570 if(module==1) {
571 unsigned int theModule;
572 if(thetaModule>21) {
573 theModule = 43 - thetaModule;
575 pos.setZ(-pos.z());
576 } else {
577 theModule = thetaModule;
578 }
579
580 double IRShift;
581 if(theModule==21) {
582 IRShift = 0;
583
584 } else if(theModule==20) {
585 IRShift = 2.5;
586
587 } else {
588 IRShift = 5.;
589
590 }
591
592
597
598 double theta01,theta23,length,d;
599 theta01 = (center01-IR).theta();
600 theta23 = (center23-IR).theta();
601 length = (center01-IR).mag();
602 d = length*
tan(theta01-theta23);
603
604 double xIn;
605 xIn = length*
tan(theta01-(pos-IR).theta())-d/2;
606 if (xIn < (-d/2) && theModule!=21){
607
608 theModule=theModule+1 ;
609
611 if(theModule==21) {
612 IRShift = 0;
613
614 } else if(theModule==20) {
615 IRShift = 2.5;
616
617 } else {
618 IRShift = 5.;
619
620 }
622 IR=IR1;
625 theta01 = (center01-IR).theta();
626 theta23 = (center23-IR).theta();
627 length = (center01-IR).mag();
628 d = length*
tan(theta01-theta23);
629
630 xIn = length*
tan(theta01-(pos-IR).theta())-d/2;
631
632 } else if (xIn > (d/2) && theModule!=0 ){
633
634 theModule=theModule-1 ;
636 if(theModule==21) {
637 IRShift = 0;
638
639 } else if(theModule==20) {
640 IRShift = 2.5;
641
642 } else {
643 IRShift = 5.;
644
645 }
646
648 IR=IR1;
651 theta01 = (center01-IR).theta();
652 theta23 = (center23-IR).theta();
653 length = (center01-IR).mag();
654 d = length*
tan(theta01-theta23);
655
656 xIn = length*
tan(theta01-(pos-IR).theta())-d/2;
657 }
658
659 double yIn;
660
661
662
663 yIn = pos.phi() < 0 ? pos.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843
664 : pos.phi()*180./CLHEP::pi-phiModule*3.-1.843;
665
666 if(nphi==0&&yIn>100) yIn=yIn-360;
667 if(nphi==119&&yIn<-100) yIn=yIn+360;
668
669 if(yIn<-1.5-0.15){
670
671 if(nphi!=0) phiModule=phiModule-1;
672 if(nphi==0) phiModule=119;
673 yIn = pos.phi() < 0 ? pos.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843
674 : pos.phi()*180./CLHEP::pi-phiModule*3.-1.843;
675 if(phiModule==0&&yIn>100) yIn=yIn-360;
676 if(phiModule==119&&yIn<-100) yIn=yIn+360;
677
678 }
679
680 if(yIn>1.5-0.15){
681
682
683 if(nphi!=119) phiModule=phiModule+1;
684 if(nphi==119) phiModule=0;
685
686 yIn = pos.phi() < 0 ? pos.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843
687 : pos.phi()*180./CLHEP::pi-phiModule*3.-1.843;
688 if(phiModule==0&&yIn>100) yIn=yIn-360;
689 if(phiModule==119&&yIn<-100) yIn=yIn+360;
690 }
691
692 enecor=m_par[theModule][0]
693 +m_par[theModule][1]*xIn
694 +m_par[theModule][2]*xIn*xIn
695 +m_par[theModule][3]*xIn*xIn*xIn
696 +m_par[theModule][4]*xIn*xIn*xIn*xIn
697 +m_par[theModule][5]*xIn*xIn*xIn*xIn*xIn
698 +m_par[theModule][6]*xIn*xIn*xIn*xIn*xIn*xIn
699
700
701 +m_par[theModule][7]*yIn
702 +m_par[theModule][8]*yIn*yIn
703 +m_par[theModule][9]*yIn*yIn*yIn
704 +m_par[theModule][10]*yIn*yIn*yIn*yIn;
705
706 } else{
707 enecor=1;
708 }
709 }
710
711
712 double energyCC= energyC/(lnEcor*enecor);
713
715
716 if(ntOut==true){
717 m_ef=energyCC;
719 m_ec=energyC;
720 m_ct=lnEcor ;
721 m_cedge=enecor;
722 m_tuple1->write();
723 }
724
725
726 }
727
729 if(!recEmcShowerCol){
730 return( StatusCode::SUCCESS);
731 }
733 if(!dstEmcShowerCol){
734 return( StatusCode::SUCCESS);
735 }
736
737
738 if(recEmcShowerCol->size()!=dstEmcShowerCol->size())return SUCCESS;
739 for(int i=0;i<recEmcShowerCol->size();i++){
740 RecEmcShowerCol::iterator iter_rec = recEmcShowerCol->begin()+i;
741 DstEmcShowerCol::iterator iter_dst = dstEmcShowerCol->begin()+i;
742
743 (*iter_dst)->setEnergy((*iter_rec)->energy());
744 (*iter_dst)->setStatus((*iter_rec)->status());
745
746
747 }
748 return( StatusCode::SUCCESS);
749}
double tan(const BesAngle a)
EvtRecTrackCol::iterator EvtRecTrackIterator
HepPoint3D position() 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 theta_module(const Identifier &id)
static unsigned int phi_module(const Identifier &id)
virtual HepPoint3D GetCrystalPoint(const Identifier &id, const int i) const =0
_EXTERN_ std::string DstEmcShowerCol
_EXTERN_ std::string EvtRecEvent
_EXTERN_ std::string EvtRecTrackCol
_EXTERN_ std::string RecEmcShowerCol