00001 #include <cmath>
00002 #include "GaudiKernel/MsgStream.h"
00003
00004 #include "GaudiKernel/SmartDataPtr.h"
00005
00006 #include "EventModel/Event.h"
00007 #include "EventModel/EventHeader.h"
00008 #include "Identifier/Identifier.h"
00009 #include "EmcRawEvent/EmcDigi.h"
00010 #include "Identifier/EmcID.h"
00011 #include "RawEvent/RawDataUtil.h"
00012 #include "HltAlgorithms/EFGlobalEnergy.h"
00013
00014 using namespace Event;
00015
00016 EFGlobalEnergy::EFGlobalEnergy(const std::string& name, ISvcLocator* pSvcLocator) :
00017 IEFAlgorithm(name, pSvcLocator) {
00018
00019 int output=(m_output%100000)/10000;
00020 MsgStream log(msgSvc(), name);
00021 msgSvc()->setOutputLevel(name,output);
00022 m_etot = new CriteriaItemValue;
00023 m_ebar = new CriteriaItemValue;
00024 m_eend = new CriteriaItemValue;
00025 m_ebal = new CriteriaItemValue;
00026 }
00027
00028 EFGlobalEnergy::~EFGlobalEnergy() {
00029 delete m_etot;
00030 delete m_ebar;
00031 delete m_eend;
00032 delete m_ebal;
00033 }
00034
00035 StatusCode EFGlobalEnergy::initialize(){
00036
00037 MsgStream log(msgSvc(), name());
00038 log << MSG::INFO << "in initialize()" << endreq;
00039
00040 IEFAlgorithm::initialize();
00041
00042 StatusCode sc;
00043 sc = m_HltStoreSvc->put("etot", m_etot);
00044 if ( sc.isFailure() ) {
00045 log << MSG::ERROR << "m_HltStoreSvc->put(etot) wrong" << endreq;
00046 return sc;
00047 }
00048 sc = m_HltStoreSvc->put("ebar", m_ebar);
00049 if ( sc.isFailure() ) {
00050 log << MSG::ERROR << "m_HltStoreSvc->put(ebar) wrong" << endreq;
00051 return sc;
00052 }
00053 sc = m_HltStoreSvc->put("eend", m_eend);
00054 if ( sc.isFailure() ) {
00055 log << MSG::ERROR << "m_HltStoreSvc->put(eend) wrong" << endreq;
00056 return sc;
00057 }
00058 sc = m_HltStoreSvc->put("ebal", m_ebal);
00059 if ( sc.isFailure() ) {
00060 log << MSG::ERROR << "m_HltStoreSvc->put(ebal) wrong" << endreq;
00061 return sc;
00062 }
00063
00064 sc = service("EmcCalibConstSvc", m_emcCalibConstSvc);
00065 if(sc != StatusCode::SUCCESS) {
00066 log << MSG::ERROR << "Can't get EmcCalibConstSvc." << endreq;
00067 m_emcCalibConstSvc=0;
00068 }
00069
00070 if(m_rawDigiSvc){
00071 if(m_rawDigiSvc->isOnlineMode()){
00072 char cbeam[20];
00073 sprintf(cbeam,"%f", m_beam);
00074 setenv("BEPCII_INFO.BER_PRB",cbeam,0);
00075 setenv("BEPCII_INFO.BPR_PRB",cbeam,0);
00076 }
00077 }
00078 return sc;
00079 }
00080
00081 StatusCode EFGlobalEnergy::execute() {
00082
00083
00084
00085
00086 MsgStream log(msgSvc(), name());
00087
00088 float e1=0,e2=0;
00089 char* electron=getenv("BEPCII_INFO.BER_PRB");
00090 if(electron){
00091 e1=atof(electron);
00092 }
00093 else{
00094 log << MSG::ERROR << "Cannot get beam energy (e-)! Please call online people" << endreq;
00095 }
00096 char* positron=getenv("BEPCII_INFO.BPR_PRB");
00097 if(positron){
00098 e2=atof(positron);
00099 }
00100 else{
00101 log << MSG::ERROR << "Cannot get beam energy (e+)! Please call online people" << endreq;
00102 }
00103 if(e1>0.1&&e2>0.1&&m_beam>0) {
00104 m_beam = 0.5*(e1+e2);
00105 }
00106 else{
00107
00108 }
00109 log << MSG::INFO << "beam energy = " << m_beam << endreq;
00110
00111
00112 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00113 if (!eventHeader) {
00114 log << MSG::FATAL << "Could not find Event Header" << endreq;
00115 return( StatusCode::FAILURE);
00116 }
00117
00118
00119
00120 bool calFlag=false;
00121 EmcDigiCol* emcDigiCol=0;
00122 if(m_rawDigiSvc){
00123 if(m_rawDigiSvc->isOnlineMode()){
00124 emcDigiCol=&(m_rawDigiSvc->getEmcDigiVec(EmcRawDataProvider::DropLow|EmcRawDataProvider::CutTime|EmcRawDataProvider::DropHot|EmcRawDataProvider::DropDead|EmcRawDataProvider::Redo));
00125 calFlag = true;
00126 }
00127 }
00128 if(!emcDigiCol){
00129 SmartDataPtr<EmcDigiCol> emcDigi(eventSvc(),"/Event/Digi/EmcDigiCol");
00130 if (!emcDigi) {
00131 log << MSG::FATAL << "Could not find Emc digi!!" << endreq;
00132 return( StatusCode::FAILURE);
00133 }
00134
00135 emcDigiCol=emcDigi;
00136 }
00137
00138 Identifier id;
00139 double adc,etot=0.,ebarrel=0.,eendcap=0.;
00140 double energyx=0.,energyy=0.,energyz=0.;
00141 double ebalance;
00142 unsigned int idBarrel_Endcap,itheta,iphi;
00143 double ewest=0.,eeast=0.;
00144
00145
00146 EmcDigiCol::iterator iterEMC=emcDigiCol->begin();
00147 while(iterEMC!= emcDigiCol->end()) {
00148 id=(*iterEMC)->identify();
00149 idBarrel_Endcap=EmcID::barrel_ec(id);
00150 itheta=EmcID::theta_module(id);
00151 iphi=EmcID::phi_module(id);
00152 adc=RawDataUtil::EmcCharge((*iterEMC)->getMeasure(),(*iterEMC)->getChargeChannel());
00153 adc /= 1000. ;
00154 if(!calFlag&&m_emcCalibConstSvc) {
00155 int index = m_emcCalibConstSvc->getIndex(idBarrel_Endcap,itheta,iphi);
00156 double adc2e = m_emcCalibConstSvc->getDigiCalibConst(index);
00157 log << MSG::DEBUG << "adc= " << adc << " and calibration constant: " << adc2e << " at " << idBarrel_Endcap << " " << itheta << " " << iphi << endreq;
00158 adc *= adc2e;
00159 }
00160
00161 etot += adc;
00162 double theta=0,phi=0;
00163
00164 if(idBarrel_Endcap==1){
00165 theta=(34.+112.*(itheta+0.5)/44.)/180.*3.1415926;
00166 phi=(iphi+0.5)/120.*6.2831852;
00167 ebarrel +=adc;
00168 if(itheta<22) eeast+=adc;
00169 else ewest+=adc;
00170 }
00171 else{
00172 theta=(90.+(90.-(itheta+0.5)/6.*(34.-21.56)+21.56)*(idBarrel_Endcap-1))
00173 /180.*3.1415926;
00174 if(itheta==0||itheta==1)phi=(iphi+0.5)/64.*6.2831852;
00175 if(itheta==2||itheta==3)phi=(iphi+0.5)/80.*6.2831852;
00176 if(itheta==4||itheta==5)phi=(iphi+0.5)/96.*6.2831852;
00177 eendcap +=adc;
00178 if(idBarrel_Endcap==0) {
00179 eeast+=adc;
00180 }
00181 else {
00182 ewest+=adc;
00183 }
00184 }
00185 energyx +=adc*sin(theta)*cos(phi);
00186 energyy +=adc*sin(theta)*sin(phi);
00187 energyz +=adc*cos(theta);
00188
00189
00190 iterEMC++;
00191 }
00192 if(etot>0)
00193 ebalance=sqrt(energyx*energyx+energyy*energyy+energyz*energyz)/etot;
00194 else
00195 ebalance=0;
00196
00197 log << MSG::INFO << "etot=" << etot << "(" << etot/2./abs(m_beam) << "); ebarrel=" <<ebarrel
00198 <<"; eendcap="<<eendcap <<"; ebalanece="<<ebalance<<endreq;
00199
00200
00201 m_etot->setValue(etot/2./abs(m_beam));
00202 m_ebar->setValue(ebarrel/2./abs(m_beam));
00203 m_eend->setValue(eendcap/2./abs(m_beam));
00204 m_ebal->setValue(ebalance);
00205
00206 m_ef->appToEFVec(etot, 31);
00207 m_ef->appToEFVec(ebalance,34);
00208 m_ef->appToEFVec(ebarrel,32);
00209 m_ef->appToEFVec(eendcap,33);
00210 m_ef->appToEFVec(eeast,35);
00211 m_ef->appToEFVec(ewest,36);
00212 m_ef->setVecBit(true, 0, 1);
00213 m_ef->addToEFVec(7<<16,1);
00214 m_run=1;
00215
00216 return StatusCode::SUCCESS;
00217 }
00218
00219 StatusCode EFGlobalEnergy::finalize() {
00220 MsgStream log(msgSvc(), name());
00221 log << MSG::INFO << "in finalize()" << endmsg;
00222 return StatusCode::SUCCESS;
00223 }
00224
00225 void EFGlobalEnergy::reset() {
00226
00227 if(m_run) {
00228 m_etot->reset();
00229 m_ebar->reset();
00230 m_eend->reset();
00231 m_ebal->reset();
00232 m_run=0;
00233 }
00234 return;
00235 }