BOSS 7.0.8
BESIII Offline Software System
Loading...
Searching...
No Matches
EmcTCFinder.cxx
Go to the documentation of this file.
1//---------------------------------------------------------------------------//
2//// BOOST --- BESIII Object_Oriented Simulation Tool /
3////
4////---------------------------------------------------------------------------/
5////
6////Description:
7////Author: Caogf
8////Created: Feb, 2006
9////Modified:
10////Comment:
11////
12//
13#include "Trigger/EmcTCFinder.h"
14
15#include "GaudiKernel/ISvcLocator.h"
16#include "GaudiKernel/Bootstrap.h"
17#include "GaudiKernel/IDataProviderSvc.h"
18#include "CLHEP/Random/RanecuEngine.h"
19#include "CLHEP/Random/RandGauss.h"
20
22#include "Identifier/EmcID.h"
24#include "RawEvent/DigiEvent.h"
26#include <math.h>
27using namespace std;
28using namespace CLHEP;
29
30EmcTCFinder* EmcTCFinder::emc_Pointer = 0;
31
33 if(!emc_Pointer) emc_Pointer = new EmcTCFinder();
34 return emc_Pointer;
35}
36
38{
39 ISvcLocator* svcLocator = Gaudi::svcLocator();
40 IRealizationSvc *tmpReal;
41 StatusCode status = svcLocator->service("RealizationSvc",tmpReal);
42 if (!status.isSuccess())
43 {
44 cout << "FATAL: Could not initialize Realization Service" << endl;
45 } else {
46 m_RealizationSvc=dynamic_cast<RealizationSvc*>(tmpReal);
47 }
48
49 // Get EmcCalibConstSvc.
50 status = svcLocator->service("EmcCalibConstSvc", emcCalibConstSvc);
51 if(status != StatusCode::SUCCESS) {
52 cout << "EmcRecDigit2Hit Error: Can't get EmcCalibConstSvc." << endl;
53 }
54}
56{
57}
59{
60 //double energy;
61 double tot_adc = 0.;
62 double adc = 0., adc1 = 0., tdc = 0.;
63 unsigned int measure;
64 //reset tc energy of barrel and endcap
65 for(int i=0;i<TrigConf::TCTHETANO_B;i++)
66 for(int j=0;j<TrigConf::TCPHINO_B;j++) {
67 BTCEnergy[i][j] = 0;
68 BTCEnergy_adc[i][j] = 0;
69 }
70 for(int i=0;i<TrigConf::TCTHETANO_E;i++)
71 for(int j=0;j<TrigConf::TCPHINO_E;j++)
72 {
73 EETCEnergy[i][j] = 0;
74 WETCEnergy[i][j] = 0;
75 EETCEnergy_adc[i][j] = 0;
76 WETCEnergy_adc[i][j] = 0;
77 }
78 EmcDigiCol::iterator iter3;
79 Identifier id;
80 for (iter3=emcDigiCol->begin();iter3!= emcDigiCol->end();iter3++) {
81 id=(*iter3)->identify();
82
83 unsigned int module;
84 unsigned int theta;
85 unsigned int phi;
86 module = EmcID::barrel_ec(id);
87 theta = EmcID::theta_module(id);
88 phi = EmcID::phi_module(id);
89 adc = double ((*iter3)->getChargeChannel());
90 adc1 = double ((*iter3)->getChargeChannel());
91 measure = (*iter3)->getMeasure();
92 tdc = (*iter3)->getTimeChannel();
93
94 int index = emcCalibConstSvc->getIndex(module,theta,phi);
95 //divided by electronics calibration constants and times trigger gain
96 //double trgGain = ((50./255.)*(m_RealizationSvc->getTrgGain(index))/60.)*5*0.5;
97 double trgGain = m_RealizationSvc->getTrgGain(index);
98 std::cout <<"partId, thetaId, phiId, trgGain: " << module << ", " << theta << ", " << phi << ", " << trgGain << std::endl;
99 //if((*iter3)->getMeasure()==0) adc = adc*(m_RealizationSvc->getEleCalib(index,7))*(trgGain)*800.*2/65535;
100 //else if((*iter3)->getMeasure()==1) adc = adc*(m_RealizationSvc->getEleCalib(index,4))*(trgGain)*800.*2/65535;
101 //else adc = adc*(m_RealizationSvc->getEleCalib(index,1))*(trgGain)*800.*2/65535;
102 if((*iter3)->getMeasure()==0) adc = adc*2*800.*2/65535.*(trgGain);
103 else if((*iter3)->getMeasure()==1) adc = adc*16*800.*2/65535*(trgGain);
104 else adc = adc*64*800.*2/65535*(trgGain);
105
106 //double trgConst = m_RealizationSvc->getTrgConst(index);
107 //if((*iter3)->getMeasure()==0) adc1 = adc1*2*800.*2/65535.*(trgConst)*0.333;
108 //else if((*iter3)->getMeasure()==1) adc1 = adc1*16*800.*2/65535*(trgConst)*0.333;
109 //else adc1 = adc1*64*800.*2/65535*(trgConst)*0.333;
110
111 //std::cout << "eCalib: " << m_RealizationSvc->getEleCalib(index,7) << "," << m_RealizationSvc->getEleCalib(index,4) << "," << m_RealizationSvc->getEleCalib(index,1) << " TrigGain: " << m_RealizationSvc->getTrgGain(index) << std::endl;
112 //std::cout << " TrgGain: " << trgGain << " adc: " << adc << std::endl;
113 //adc = RawDataUtil::EmcCharge(measure,adc); //MeV
114 //energy = RandGauss::shoot(energy,0.1*energy); //add 10% smear in each ADC channel
115
116 int TCThetaId = getTCThetaId(module,theta,phi);
117 int TCPhiId = getTCPhiId(module,theta,phi);
118 //Get energy of each trigger cell
119 if(module==1) BTCEnergy[TCThetaId][TCPhiId] += adc;
120 if(module==0) EETCEnergy[TCThetaId][TCPhiId] += adc;
121 if(module==2) WETCEnergy[TCThetaId][TCPhiId] += adc;
122 if(module==1) BTCEnergy_adc[TCThetaId][TCPhiId] += adc;
123 if(module==0) EETCEnergy_adc[TCThetaId][TCPhiId] += adc;
124 if(module==2) WETCEnergy_adc[TCThetaId][TCPhiId] += adc;
125 }
126
127 //subtract ped
128/*
129 for(int i=0;i<TrigConf::TCTHETANO_B;i++)
130 for(int j=0;j<TrigConf::TCPHINO_B;j++) {
131 if((BTCEnergy[i][j] - 0x0a) < 0) BTCEnergy[i][j] = 0;
132 else BTCEnergy[i][j] = BTCEnergy[i][j] - 0x0a;
133 }
134 for(int i=0;i<TrigConf::TCTHETANO_E;i++)
135 for(int j=0;j<TrigConf::TCPHINO_E;j++)
136 {
137 if((EETCEnergy[i][j] - 0x0a) < 0) EETCEnergy[i][j] = 0;
138 else EETCEnergy[i][j] = EETCEnergy[i][j] - 0x0a;
139 if((WETCEnergy[i][j] - 0x0a) < 0) WETCEnergy[i][j] = 0;
140 else WETCEnergy[i][j] = WETCEnergy[i][j] - 0x0a;
141 }*/
142}
143
144void EmcTCFinder::setEmcTC(std::vector<uint32_t> vTC) {
145 //reset tc energy of barrel and endcap
146 for(int i=0;i<TrigConf::TCTHETANO_B;i++)
147 for(int j=0;j<TrigConf::TCPHINO_B;j++) {
148 BTC[i][j] = 0;
149 }
150 for(int i=0;i<TrigConf::TCTHETANO_E;i++)
151 for(int j=0;j<TrigConf::TCPHINO_E;j++)
152 {
153 EETC[i][j] = 0;
154 WETC[i][j] = 0;
155 }
156 for(std::vector<uint32_t>::iterator iter = vTC.begin(); iter != vTC.end(); iter++) {
157 int par_TC = (*iter & 0xFF0000) >> 16;
158 int the_TC = (*iter & 0xFF00) >> 8;
159 int phi_TC = (*iter & 0xFF);
160 if(par_TC == 0) EETC[the_TC][phi_TC] = 1;
161 if(par_TC == 1) BTC[the_TC][phi_TC] = 1;
162 if(par_TC == 2) WETC[the_TC][phi_TC] = 1;
163 }
164}
165
166void EmcTCFinder::setEmcBE(std::vector<double> vBE) {
167 //reset block energy
168 for(int i = 0; i < 16; i++) {
169 BlkE[i] = 0;
170 }
171 if(vBE.size() != 0) {
172 if(vBE.size() != 16 ) std::cerr << "The number of block is not equal 16, please check it (in EmcTCFinder::setEmcBE() )" << std::endl;
173 for(int i = 0; i < vBE.size(); i++) {
174 BlkE[i] = vBE[i];
175 }
176 }
177}
178
179int EmcTCFinder::getTCThetaId(int partId,int ThetaNb,int PhiNb)
180{
181 //Note: There are 44*120(theta*phi) crystals and they are divided 11*30 trigger cell in
182 //barrel emc. So each trigger cell includes 4*4 crystals
183 if(partId==1)
184 {
185 //if(ThetaNb<2)
186 //TCThetaNb = 0;
187 //else if(ThetaNb>41)
188 // TCThetaNb = 11;
189 //else
190 //TCThetaNb = int((ThetaNb+2)/4);
191 TCThetaNb = (int)ThetaNb/4;
192 }
193/* if(partId==0)
194 {
195 if(ThetaNb<3) TCThetaNb = 0;
196 if(ThetaNb>3) TCThetaNb = 1;
197 if(ThetaNb==3)
198 {
199 int n = PhiNb%10;
200 if(n>2&&n<7)
201 TCThetaNb = 0;
202 else TCThetaNb = 1;
203 }
204 }
205 if(partId==2)
206 {
207 if(ThetaNb<3) TCThetaNb = 0;
208 if(ThetaNb>3) TCThetaNb = 1;
209 if(ThetaNb==3)
210 {
211 int n = PhiNb%10;
212 if(n>2&&n<7)
213 TCThetaNb = 0;
214 else TCThetaNb = 1;
215 }
216 }
217*/
218 if(partId==0) TCThetaNb = 0;
219 if(partId==2) TCThetaNb = 0;
220 return TCThetaNb;
221}
222int EmcTCFinder::getTCPhiId(int partId,int ThetaNb,int PhiNb)
223{
224 if(partId==1)
225 TCPhiNb = int(PhiNb/4);
226/*
227 //trigger cell in endcaps theta*phi = 2*16
228 if(partId==0)
229 {
230 if(ThetaNb<2) TCPhiNb = int(PhiNb/4);
231 if(ThetaNb>=2&&ThetaNb<4) TCPhiNb = int(PhiNb/5);
232 if(ThetaNb>=4) TCPhiNb = int(PhiNb/6);
233 }
234 if(partId==2)
235 {
236 if(ThetaNb<2) TCPhiNb = int(PhiNb/4);
237 if(ThetaNb>=2&&ThetaNb<4) TCPhiNb = int(PhiNb/5);
238 if(ThetaNb>=4) TCPhiNb = int(PhiNb/6);
239 }
240*/
241 // trigger cell in endcaps theta*phi = 1*32
242 if(partId==0)
243 {
244 if(ThetaNb<2) TCPhiNb = int(PhiNb/2);
245 if(ThetaNb==2)
246 {
247 int quot = int(PhiNb/5);
248 int rema = PhiNb%5;
249 if(rema <= 1) TCPhiNb = 2*quot;
250 if(rema > 1) TCPhiNb = 2*quot + 1;
251 }
252 if(ThetaNb==3)
253 {
254 int quot = int(PhiNb/5);
255 int rema = PhiNb%5;
256 if(rema <= 2) TCPhiNb = 2*quot;
257 if(rema > 2) TCPhiNb = 2*quot + 1;
258 }
259 if(ThetaNb>=4) TCPhiNb = int(PhiNb/3);
260 }
261 if(partId==2)
262 {
263 if(ThetaNb<2) TCPhiNb = int(PhiNb/2);
264 if(ThetaNb==2)
265 {
266 int quot = int(PhiNb/5);
267 int rema = PhiNb%5;
268 if(rema <= 1) TCPhiNb = 2*quot;
269 if(rema > 1) TCPhiNb = 2*quot + 1;
270 }
271 if(ThetaNb==3)
272 {
273 int quot = int(PhiNb/5);
274 int rema = PhiNb%5;
275 if(rema <= 2) TCPhiNb = 2*quot;
276 if(rema > 2) TCPhiNb = 2*quot + 1;
277 }
278 if(ThetaNb>=4) TCPhiNb = int(PhiNb/3);
279 }
280
281 return TCPhiNb;
282}
283int EmcTCFinder::getBLKId(int TCTheta,int TCPhi) const
284{
285 int id,parity;
286 parity = (int)TCPhi/5;
287 if(parity%2==0)
288 {
289 if(TCTheta<6) id = parity;
290 if(TCTheta>=6) id = parity+6;
291 }
292 if(parity%2!=0)
293 {
294 if(TCTheta<5) id = parity;
295 if(TCTheta>=5) id = parity+6;
296 }
297
298 return id;
299}
ObjectVector< EmcDigi > EmcDigiCol
Definition: EmcDigi.h:43
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
static unsigned int barrel_ec(const Identifier &id)
Values of different levels (failure returns 0)
Definition: EmcID.cxx:38
static unsigned int theta_module(const Identifier &id)
Definition: EmcID.cxx:43
static unsigned int phi_module(const Identifier &id)
Definition: EmcID.cxx:48
int getTCThetaId(int partId, int ThetaNb, int PhiNb)
static EmcTCFinder * get_Emc(void)
Definition: EmcTCFinder.cxx:32
void setEmcBE(std::vector< double > vBE)
int getTCPhiId(int partId, int ThetaNb, int PhiNb)
void setEmcTC(std::vector< uint32_t > vTC)
void setEmcDigi(EmcDigiCol *emcDigiCol)
Definition: EmcTCFinder.cxx:58
int getBLKId(int TCTheta, int TCPhi) const
virtual int getIndex(unsigned int PartId, unsigned int ThetaIndex, unsigned int PhiIndex) const =0
double getTrgGain(int cry_id)