333 {
334 MsgStream log(
msgSvc(), name());
335 log << MSG::INFO << "in execute()" << endreq;
336
337 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
338 int runNo=eventHeader->runNumber();
339
340
343 unsigned int runFrom=m_EmcShEnCalibSvc -> getRunFrom();
344 unsigned int runTo=m_EmcShEnCalibSvc -> getRunTo();
345 if (!m_ReadFile) {
346 runFrom0 = runFrom;
347 runTo0 = runTo;
348 }
349
350 if (runNum>=runFrom0&&runNum<=runTo0&&m_ReadFile==true){
351
352
353
354
355 } else {
356 if (runFrom0!=runFrom) runFrom0 = runFrom;
357 if (runTo0!=runTo) runTo0 = runTo;
358
359
360 cout <<
"current run=" <<
runNo<<endl;
361 cout << "in AbsCor open getPi0CalibFile()= " << m_EmcShEnCalibSvc -> getPi0CalibFile()<<endl;
362 cout << "open getSingleGammaCalibFile()= " << m_EmcShEnCalibSvc -> getSingleGammaCalibFile()<<endl;
363 cout <<"getRunFrom()="<< runFrom<<endl;
364 cout <<"getRunTo()="<< runTo<<endl;
365
366
367 string CorFunparaPath;
368 if(MCuseTof==1){
369 CorFunparaPath = m_EmcShEnCalibSvc -> getSingleGammaCalibFile();
370 }
371 if(MCuseTof==0){
372 CorFunparaPath = getenv("ABSCORROOT");
373 CorFunparaPath += "/share/evsetCorFunctionPar.txt";
374 }
375
377 in2corfun.open(CorFunparaPath.c_str());
378 assert(in2corfun);
379
380 for(int i=0;i<28;i++){
381 in2corfun>>m_corFunPar[i][0];
382 in2corfun>>m_corFunPar[i][1];
383 in2corfun>>m_corFunPar[i][2];
384 in2corfun>>m_corFunPar[i][3];
385 in2corfun>>m_corFunPar[i][4];
386 in2corfun>>m_corFunPar[i][5];
387 }
388 in2corfun.close();
389
390
391 string DataPathc3p;
392 DataPathc3p = getenv("ABSCORROOT");
393 DataPathc3p += "/share/c3p.txt";
394
395 string DataPathc3ptof;
396 DataPathc3ptof = m_EmcShEnCalibSvc -> getPi0CalibFile();
397
399 inc3p.exceptions( ifstream::failbit | ifstream::badbit );
400 if(usetof==0)inc3p.open(DataPathc3p.c_str(),ios::in);
401 if(usetof==1)inc3p.open(DataPathc3ptof.c_str(),ios::in);
402 for(int i=0; i<4; i++){
403 double am,ame;
404 inc3p>>am;
405 inc3p>>ame;
407 }
408
409
410 m_ReadFile=true;
411 cout<<"-----------read parameter file-------"<<endl;
412 }
413
414
417 if(evtRecEvent->totalTracks()>evtRecTrkCol->size())return SUCCESS;
418 if(evtRecEvent->totalTracks()>50)return SUCCESS;
419
421 ISvcLocator* svcLocator = Gaudi::svcLocator();
422 StatusCode sc = svcLocator->service("EmcRecGeoSvc",iGeoSvc);
423 if(sc!=StatusCode::SUCCESS) {
424 cout<<"Error: Can't get EmcRecGeoSvc"<<endl;
425 }
426
427 for(int i = 0; i< evtRecEvent->totalTracks(); i++) {
429 if(!(*itTrk)->isEmcShowerValid())continue;
431
432
433
434
435
436
437
438
439
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
466
467 unsigned int module, ntheta, nphi;
471
472
473
476
477
478 double e5x5=emcTrk->
e5x5();
479 double etof=0;
480
481 if(usetof==1 && (*itTrk)->isTofTrackValid()){
482 SmartRefVector<RecTofTrack> recTofTrackVec=(*itTrk)->tofTrack();
483 if(!recTofTrackVec.empty())etof=recTofTrackVec[0]->energy();
484 if(etof>100.)etof=0;
485 }
486
487 double energyC;
488
489 int thetaId;
490 if (module==0||module==2) thetaId = thetaModule;
491 if (module==1 && thetaModule<=21) thetaId = thetaModule + 6;
492 if (module==1 && thetaModule>21) thetaId = 43 - thetaModule + 6;
493 double DthetaId=double(thetaId);
494
495 if(MCuseTof==1){
496 if (thetaId<6) etof=0.0;
497 if (MCCorUseFunction==1){
498 energyC=ECorrFunctionMC(e5x5+etof,DthetaId);
499 } else if (MCCorUseFunction==0){
500 energyC=ECorrMC(e5x5+etof,DthetaId);
501 }
502 }
503 if (MCuseTof==0) {
504 if (MCCorUseFunction==1){
505 energyC=ECorrFunctionMC(e5x5,DthetaId);
506 } else if (MCCorUseFunction==0){
507 energyC=ECorrMC(e5x5,DthetaId);
508 }
509 }
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531 double lnE = std::log(energyC);
532
533 if (energyC>1.0) lnE=std::log(1.0);
534 if (energyC<0.05) lnE=std::log(0.05);
535
536 double lnEcor=1.0;
537 if(dopi0Cor==1){
538 if(
runNo>0 && dodatacor==1) {
539 lnEcor =
ai[0]+
ai[1]*lnE+
ai[2]*lnE*lnE+
ai[3]*lnE*lnE*lnE;
540 }
541 }
542 if(lnEcor<0.5)continue;
543
544
545
546
548
549
550 double enecor=1.;
551 if(hotcellmask==1){
552 for(int ih=0;ih<10;ih++){
553 if(hrunstart[ih]==-1||hrunend[ih]==-1||hcell[ih]==-1)continue;
556 }
557 }
558
559 if(edgecor==1){
560
561 if(module==1) {
562 unsigned int theModule;
563 if(thetaModule>21) {
564 theModule = 43 - thetaModule;
566 pos.setZ(-pos.z());
567 } else {
568 theModule = thetaModule;
569 }
570
571 double IRShift;
572 if(theModule==21) {
573 IRShift = 0;
574
575 } else if(theModule==20) {
576 IRShift = 2.5;
577
578 } else {
579 IRShift = 5.;
580
581 }
582
583
588
589 double theta01,theta23,length,d;
590 theta01 = (center01-IR).theta();
591 theta23 = (center23-IR).theta();
592 length = (center01-IR).mag();
593 d = length*
tan(theta01-theta23);
594
595 double xIn;
596 xIn = length*
tan(theta01-(pos-IR).theta())-d/2;
597 if (xIn < (-d/2) && theModule!=21){
598
599 theModule=theModule+1 ;
600
602 if(theModule==21) {
603 IRShift = 0;
604
605 } else if(theModule==20) {
606 IRShift = 2.5;
607
608 } else {
609 IRShift = 5.;
610
611 }
613 IR=IR1;
616 theta01 = (center01-IR).theta();
617 theta23 = (center23-IR).theta();
618 length = (center01-IR).mag();
619 d = length*
tan(theta01-theta23);
620
621 xIn = length*
tan(theta01-(pos-IR).theta())-d/2;
622
623 } else if (xIn > (d/2) && theModule!=0 ){
624
625 theModule=theModule-1 ;
627 if(theModule==21) {
628 IRShift = 0;
629
630 } else if(theModule==20) {
631 IRShift = 2.5;
632
633 } else {
634 IRShift = 5.;
635
636 }
637
639 IR=IR1;
642 theta01 = (center01-IR).theta();
643 theta23 = (center23-IR).theta();
644 length = (center01-IR).mag();
645 d = length*
tan(theta01-theta23);
646
647 xIn = length*
tan(theta01-(pos-IR).theta())-d/2;
648 }
649
650 double yIn;
651
652
653
654 yIn = pos.phi() < 0 ? pos.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843
655 : pos.phi()*180./CLHEP::pi-phiModule*3.-1.843;
656
657 if(nphi==0&&yIn>100) yIn=yIn-360;
658 if(nphi==119&&yIn<-100) yIn=yIn+360;
659
660 if(yIn<-1.5-0.15){
661
662 if(nphi!=0) phiModule=phiModule-1;
663 if(nphi==0) phiModule=119;
664 yIn = pos.phi() < 0 ? pos.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843
665 : pos.phi()*180./CLHEP::pi-phiModule*3.-1.843;
666 if(phiModule==0&&yIn>100) yIn=yIn-360;
667 if(phiModule==119&&yIn<-100) yIn=yIn+360;
668
669 }
670
671 if(yIn>1.5-0.15){
672
673
674 if(nphi!=119) phiModule=phiModule+1;
675 if(nphi==119) phiModule=0;
676
677 yIn = pos.phi() < 0 ? pos.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843
678 : pos.phi()*180./CLHEP::pi-phiModule*3.-1.843;
679 if(phiModule==0&&yIn>100) yIn=yIn-360;
680 if(phiModule==119&&yIn<-100) yIn=yIn+360;
681 }
682
683 enecor=m_par[theModule][0]
684 +m_par[theModule][1]*xIn
685 +m_par[theModule][2]*xIn*xIn
686 +m_par[theModule][3]*xIn*xIn*xIn
687 +m_par[theModule][4]*xIn*xIn*xIn*xIn
688 +m_par[theModule][5]*xIn*xIn*xIn*xIn*xIn
689 +m_par[theModule][6]*xIn*xIn*xIn*xIn*xIn*xIn
690
691
692 +m_par[theModule][7]*yIn
693 +m_par[theModule][8]*yIn*yIn
694 +m_par[theModule][9]*yIn*yIn*yIn
695 +m_par[theModule][10]*yIn*yIn*yIn*yIn;
696
697 } else{
698 enecor=1;
699 }
700 }
701
702
703 double energyCC= energyC/(lnEcor*enecor);
704
706
707 if(ntOut==true){
708 m_ef=energyCC;
710 m_ec=energyC;
711 m_ct=lnEcor ;
712 m_cedge=enecor;
713 m_tuple1->write();
714 }
715
716
717 }
718
720 if(!recEmcShowerCol){
721 return( StatusCode::SUCCESS);
722 }
724 if(!dstEmcShowerCol){
725 return( StatusCode::SUCCESS);
726 }
727
728
729 if(recEmcShowerCol->size()!=dstEmcShowerCol->size())return SUCCESS;
730 for(int i=0;i<recEmcShowerCol->size();i++){
731 RecEmcShowerCol::iterator iter_rec = recEmcShowerCol->begin()+i;
732 DstEmcShowerCol::iterator iter_dst = dstEmcShowerCol->begin()+i;
733
734 (*iter_dst)->setEnergy((*iter_rec)->energy());
735 (*iter_dst)->setStatus((*iter_rec)->status());
736
737
738 }
739 return( StatusCode::SUCCESS);
740}
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 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)
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