107 {
108
109
110
111
112
113
114
115 MsgStream log(
msgSvc(), name() );
116 log << MSG::INFO << "in execute()" << endmsg;
117
118 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
119 int run = eventHeader->runNumber();
120
121 m_run = run;
122
123 int event = eventHeader->eventNumber();
124 if ( m_event % 1000 == 0 )
125 {
126 std::cout << m_event << "-------------TofEnergyCalib run,event: " << run << "," << event
127 << std::endl;
128 }
129 m_event++;
130
133
134 log << MSG::DEBUG << "ncharg, nneu, tottks = " << evtRecEvent->totalCharged() << " , "
135 << evtRecEvent->totalNeutral() << " , " << evtRecEvent->totalTracks() << endmsg;
136
137 if ( evtRecEvent->totalCharged() != 2 ) { return StatusCode::SUCCESS; }
138
139
141 StatusCode sc = service(
"RawDataProviderSvc",
tofDigiSvc );
142 if ( sc != StatusCode::SUCCESS )
143 {
144 log << MSG::ERROR << "TofRec Get Tof Raw Data Service Failed !! " << endmsg;
145 return StatusCode::SUCCESS;
146 }
147
148
150 StatusCode scc = service(
"TofCaliSvc",
tofCaliSvc );
151 if ( scc != StatusCode::SUCCESS )
152 {
153 log << MSG::ERROR << "TofRec Get Calibration Service Failed !! " << endmsg;
154 return StatusCode::SUCCESS;
155 }
156
157
160 if ( sc != StatusCode::SUCCESS )
161 {
162 log << MSG::ERROR << "TofRec Get QCorr Service Failed !! " << endmsg;
163 return StatusCode::SUCCESS;
164 }
165
168 if ( sc != StatusCode::SUCCESS )
169 {
170 log << MSG::ERROR << "TofRec Get QElec Service Failed !! " << endmsg;
171 return StatusCode::SUCCESS;
172 }
173
174
175 std::vector<TofData*> tofDataVec;
177 const unsigned int ADC_MASK = 0x00001FFF;
178
179 for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
180 {
182
183 if ( !( *itTrk )->isExtTrackValid() ) continue;
184 RecExtTrack* extTrk = ( *itTrk )->extTrack();
185
186 if ( !( *itTrk )->isTofTrackValid() ) continue;
187
188
189 int npart, tof1, tof2;
190 Hep3Vector pos1, pos2;
192
195
196
197
198 if ( tof1 >= 0 && tof1 < 88 && tof2 >= 88 && tof2 < 176 )
199 {
200 npart = 1;
201 }
202 else if ( tof1 >= 176 && tof1 < 224 )
203 {
204 tof1 -= 176;
205 npart = 2;
206 }
207 else if ( tof1 >= 224 && tof1 < 272 )
208 {
209 tof1 -= 224;
210 npart = 0;
211 }
212 else { continue; }
213
216
219
220 double zpos1 = 0, zpos2 = 0, l1 = 0, l2 = 0, lext = 0;
221 double z1, xy1, z2, xy2;
222
223 if ( npart == 1 )
224 {
225
226 const double tofDPhi = CLHEP::twopi / 88.;
227 double tofPhi1 = tof1 * tofDPhi + tofDPhi / 2;
228 double tofPhi2 = ( tof2 - 88 ) * tofDPhi;
229
230
231 double extTheta1, extTheta2, extPhi1, extPhi2;
232 extTheta1 = pos1.theta();
233 extTheta2 = pos2.theta();
234 extPhi1 = pos1.phi();
235 extPhi2 = pos2.phi();
236
237
238
239 z1 = 5 /
tan( extTheta1 );
240 xy1 = 5 /
cos( extPhi1 - tofPhi1 );
241 l1 = sqrt( z1 * z1 + xy1 * xy1 );
242 z2 = 5 /
tan( extTheta2 );
243 xy2 = 5 /
cos( extPhi2 - tofPhi2 );
244 l2 = sqrt( z2 * z2 + xy2 * xy2 );
245 zpos1 = ( pos1.z() + z1 / 2 ) / 100;
246 zpos2 = ( pos2.z() + z2 / 2 ) / 100;
247
248
250 }
251
252 else
253 {
254
255 double extTheta1 = pos1.theta();
256
257
258
259 l1 = 5. / fabs(
cos( extTheta1 ) );
260 zpos1 =
261 ( sqrt( pos1.x() * pos1.x() + pos1.y() * pos1.y() ) + 5. *
tan( extTheta1 ) / 2 ) /
262 100;
263 }
264
265 vector<TofData*>::iterator it;
266
267
268 for ( it = tofDataVec.begin(); it != tofDataVec.end(); it++ )
269 {
270
271 Identifier idd( ( *it )->identify() );
274 if ( im + layer * 88 == tof1 - 1 || im + layer * 88 == tof1 + 1 ||
275 im + layer * 88 == tof1 - 2 || im + layer * 88 == tof1 + 2 )
276 {
277
278 return StatusCode::SUCCESS;
279 }
280 }
281
282 for ( it = tofDataVec.begin(); it != tofDataVec.end(); it++ )
283 {
284
285 Identifier idd( ( *it )->identify() );
288
289 double adc1=0, adc2=0, tdc1=0, tdc2=0, ztdc=0, tofq=0;
290
291 if ( ( *it )->barrel() )
292 {
293 TofData* bTofData = *it;
294 tdc1 = bTofData->
tdc1();
295 tdc2 = bTofData->
tdc2();
296 adc1 = bTofData->
adc1();
297 adc2 = bTofData->
adc2();
298
301 if ( tofq < 100 || tofq > 10000 ) continue;
302
303 npart = 1;
304
305 }
306 else
307 {
308 TofData* eTofData = *it;
309
310 adc1 = eTofData->
adc();
311 tdc1 = eTofData->
tdc();
312 npart = 2;
313 }
314
315
316 if ( im + layer * 88 == tof1 )
317 {
318 m_adc1 = adc1;
319 m_adc2 = adc2;
320 m_tdc1 = tdc1;
321 m_tdc2 = tdc2;
322 m_ztdc = ztdc;
323 m_q = tofq;
324 m_npart = npart;
325 m_number = tof1;
326 m_zpos = zpos1;
327 m_length = l1;
328 m_length_ext = lext;
329
330 m_energy = l1 * 10.25 / 5.;
331 m_energy_ext = lext * 10.25 / 5.;
332
333
334 if (
abs( m_zpos ) < 0.6 && m_q > 500 && m_q < 4000 )
335 {
336 m_num++;
337 m_EvsQ = m_EvsQ + m_energy / m_q;
338 }
339 m_tuple->write();
340 }
341 else if ( ( *it )->barrel() && im + layer * 88 == tof2 )
342 {
343 m_adc1 = adc1;
344 m_adc2 = adc2;
345 m_tdc1 = tdc1;
346 m_tdc2 = tdc2;
347 m_ztdc = ztdc;
348 m_q = tofq;
349 m_npart = npart;
350 m_number = tof2;
351 m_zpos = zpos2;
352 m_length = l2;
353 m_length_ext = lext;
354
355 m_energy = l2 * 10.25 / 5.;
356 m_energy_ext = lext * 10.25 / 5.;
357
358 m_tuple->write();
359
360 if (
abs( m_zpos ) < 0.6 && m_q > 500 && m_q < 4000 )
361 {
362 m_num++;
363 m_EvsQ = m_EvsQ + m_energy / m_q;
364 }
365 }
366 }
367 }
368
369 return sc;
370}
EvtRecTrackCol::iterator EvtRecTrackIterator
double tan(const BesAngle a)
double cos(const BesAngle a)
ITofQElecSvc * tofQElecSvc
ITofQCorrSvc * tofQCorrSvc
IRawDataProviderSvc * tofDigiSvc
const double tof1Path() const
const Hep3Vector tof1Position() const
const int tof1VolumeNumber() const
const Hep3Vector tof2Momentum() const
const Hep3Vector tof1Momentum() const
const double tof2Path() const
const int tof2VolumeNumber() const
const Hep3Vector tof2Position() const
virtual TofDataVector & tofDataVectorTof(double estime=0)=0
virtual const double BPh(double ADC1, double ADC2, double zHit, unsigned int id)=0
virtual const double ZTDC(double tleft, double tright, unsigned id)=0
static int phi_module(const Identifier &id)
static int layer(const Identifier &id)
_EXTERN_ std::string EvtRecEvent
_EXTERN_ std::string EvtRecTrackCol