00001
00012
00013 #include "Trigger/EmcTCFinder.h"
00014
00015 #include "GaudiKernel/ISvcLocator.h"
00016 #include "GaudiKernel/Bootstrap.h"
00017 #include "GaudiKernel/IDataProviderSvc.h"
00018 #include "CLHEP/Random/RanecuEngine.h"
00019 #include "CLHEP/Random/RandGauss.h"
00020
00021 #include "Identifier/Identifier.h"
00022 #include "Identifier/EmcID.h"
00023 #include "RawEvent/RawDataUtil.h"
00024 #include "RawEvent/DigiEvent.h"
00025 #include "EmcCalibConstSvc/EmcCalibConstSvc.h"
00026 #include <math.h>
00027 using namespace std;
00028 using namespace CLHEP;
00029
00030 EmcTCFinder* EmcTCFinder::emc_Pointer = 0;
00031
00032 EmcTCFinder* EmcTCFinder::get_Emc(void) {
00033 if(!emc_Pointer) emc_Pointer = new EmcTCFinder();
00034 return emc_Pointer;
00035 }
00036
00037 EmcTCFinder::EmcTCFinder()
00038 {
00039 ISvcLocator* svcLocator = Gaudi::svcLocator();
00040 IRealizationSvc *tmpReal;
00041 StatusCode status = svcLocator->service("RealizationSvc",tmpReal);
00042 if (!status.isSuccess())
00043 {
00044 cout << "FATAL: Could not initialize Realization Service" << endl;
00045 } else {
00046 m_RealizationSvc=dynamic_cast<RealizationSvc*>(tmpReal);
00047 }
00048
00049
00050 status = svcLocator->service("EmcCalibConstSvc", emcCalibConstSvc);
00051 if(status != StatusCode::SUCCESS) {
00052 cout << "EmcRecDigit2Hit Error: Can't get EmcCalibConstSvc." << endl;
00053 }
00054 }
00055 EmcTCFinder::~EmcTCFinder()
00056 {
00057 }
00058 void EmcTCFinder::setEmcDigi(EmcDigiCol* emcDigiCol)
00059 {
00060
00061 double tot_adc = 0.;
00062 double adc = 0., adc1 = 0., tdc = 0.;
00063 unsigned int measure;
00064
00065 for(int i=0;i<TrigConf::TCTHETANO_B;i++)
00066 for(int j=0;j<TrigConf::TCPHINO_B;j++) {
00067 BTCEnergy[i][j] = 0;
00068 BTCEnergy_adc[i][j] = 0;
00069 }
00070 for(int i=0;i<TrigConf::TCTHETANO_E;i++)
00071 for(int j=0;j<TrigConf::TCPHINO_E;j++)
00072 {
00073 EETCEnergy[i][j] = 0;
00074 WETCEnergy[i][j] = 0;
00075 EETCEnergy_adc[i][j] = 0;
00076 WETCEnergy_adc[i][j] = 0;
00077 }
00078 EmcDigiCol::iterator iter3;
00079 Identifier id;
00080 for (iter3=emcDigiCol->begin();iter3!= emcDigiCol->end();iter3++) {
00081 id=(*iter3)->identify();
00082
00083 unsigned int module;
00084 unsigned int theta;
00085 unsigned int phi;
00086 module = EmcID::barrel_ec(id);
00087 theta = EmcID::theta_module(id);
00088 phi = EmcID::phi_module(id);
00089 adc = double ((*iter3)->getChargeChannel());
00090 adc1 = double ((*iter3)->getChargeChannel());
00091 measure = (*iter3)->getMeasure();
00092 tdc = (*iter3)->getTimeChannel();
00093
00094 int index = emcCalibConstSvc->getIndex(module,theta,phi);
00095
00096
00097 double trgGain = m_RealizationSvc->getTrgGain(index);
00098 std::cout <<"partId, thetaId, phiId, trgGain: " << module << ", " << theta << ", " << phi << ", " << trgGain << std::endl;
00099
00100
00101
00102 if((*iter3)->getMeasure()==0) adc = adc*2*800.*2/65535.*(trgGain);
00103 else if((*iter3)->getMeasure()==1) adc = adc*16*800.*2/65535*(trgGain);
00104 else adc = adc*64*800.*2/65535*(trgGain);
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116 int TCThetaId = getTCThetaId(module,theta,phi);
00117 int TCPhiId = getTCPhiId(module,theta,phi);
00118
00119 if(module==1) BTCEnergy[TCThetaId][TCPhiId] += adc;
00120 if(module==0) EETCEnergy[TCThetaId][TCPhiId] += adc;
00121 if(module==2) WETCEnergy[TCThetaId][TCPhiId] += adc;
00122 if(module==1) BTCEnergy_adc[TCThetaId][TCPhiId] += adc;
00123 if(module==0) EETCEnergy_adc[TCThetaId][TCPhiId] += adc;
00124 if(module==2) WETCEnergy_adc[TCThetaId][TCPhiId] += adc;
00125 }
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142 }
00143
00144 void EmcTCFinder::setEmcTC(std::vector<uint32_t> vTC) {
00145
00146 for(int i=0;i<TrigConf::TCTHETANO_B;i++)
00147 for(int j=0;j<TrigConf::TCPHINO_B;j++) {
00148 BTC[i][j] = 0;
00149 }
00150 for(int i=0;i<TrigConf::TCTHETANO_E;i++)
00151 for(int j=0;j<TrigConf::TCPHINO_E;j++)
00152 {
00153 EETC[i][j] = 0;
00154 WETC[i][j] = 0;
00155 }
00156 for(std::vector<uint32_t>::iterator iter = vTC.begin(); iter != vTC.end(); iter++) {
00157 int par_TC = (*iter & 0xFF0000) >> 16;
00158 int the_TC = (*iter & 0xFF00) >> 8;
00159 int phi_TC = (*iter & 0xFF);
00160 if(par_TC == 0) EETC[the_TC][phi_TC] = 1;
00161 if(par_TC == 1) BTC[the_TC][phi_TC] = 1;
00162 if(par_TC == 2) WETC[the_TC][phi_TC] = 1;
00163 }
00164 }
00165
00166 void EmcTCFinder::setEmcBE(std::vector<double> vBE) {
00167
00168 for(int i = 0; i < 16; i++) {
00169 BlkE[i] = 0;
00170 }
00171 if(vBE.size() != 0) {
00172 if(vBE.size() != 16 ) std::cerr << "The number of block is not equal 16, please check it (in EmcTCFinder::setEmcBE() )" << std::endl;
00173 for(int i = 0; i < vBE.size(); i++) {
00174 BlkE[i] = vBE[i];
00175 }
00176 }
00177 }
00178
00179 int EmcTCFinder::getTCThetaId(int partId,int ThetaNb,int PhiNb)
00180 {
00181
00182
00183 if(partId==1)
00184 {
00185
00186
00187
00188
00189
00190
00191 TCThetaNb = (int)ThetaNb/4;
00192 }
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218 if(partId==0) TCThetaNb = 0;
00219 if(partId==2) TCThetaNb = 0;
00220 return TCThetaNb;
00221 }
00222 int EmcTCFinder::getTCPhiId(int partId,int ThetaNb,int PhiNb)
00223 {
00224 if(partId==1)
00225 TCPhiNb = int(PhiNb/4);
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242 if(partId==0)
00243 {
00244 if(ThetaNb<2) TCPhiNb = int(PhiNb/2);
00245 if(ThetaNb==2)
00246 {
00247 int quot = int(PhiNb/5);
00248 int rema = PhiNb%5;
00249 if(rema <= 1) TCPhiNb = 2*quot;
00250 if(rema > 1) TCPhiNb = 2*quot + 1;
00251 }
00252 if(ThetaNb==3)
00253 {
00254 int quot = int(PhiNb/5);
00255 int rema = PhiNb%5;
00256 if(rema <= 2) TCPhiNb = 2*quot;
00257 if(rema > 2) TCPhiNb = 2*quot + 1;
00258 }
00259 if(ThetaNb>=4) TCPhiNb = int(PhiNb/3);
00260 }
00261 if(partId==2)
00262 {
00263 if(ThetaNb<2) TCPhiNb = int(PhiNb/2);
00264 if(ThetaNb==2)
00265 {
00266 int quot = int(PhiNb/5);
00267 int rema = PhiNb%5;
00268 if(rema <= 1) TCPhiNb = 2*quot;
00269 if(rema > 1) TCPhiNb = 2*quot + 1;
00270 }
00271 if(ThetaNb==3)
00272 {
00273 int quot = int(PhiNb/5);
00274 int rema = PhiNb%5;
00275 if(rema <= 2) TCPhiNb = 2*quot;
00276 if(rema > 2) TCPhiNb = 2*quot + 1;
00277 }
00278 if(ThetaNb>=4) TCPhiNb = int(PhiNb/3);
00279 }
00280
00281 return TCPhiNb;
00282 }
00283 int EmcTCFinder::getBLKId(int TCTheta,int TCPhi) const
00284 {
00285 int id,parity;
00286 parity = (int)TCPhi/5;
00287 if(parity%2==0)
00288 {
00289 if(TCTheta<6) id = parity;
00290 if(TCTheta>=6) id = parity+6;
00291 }
00292 if(parity%2!=0)
00293 {
00294 if(TCTheta<5) id = parity;
00295 if(TCTheta>=5) id = parity+6;
00296 }
00297
00298 return id;
00299 }