BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EmcRec.cxx
Go to the documentation of this file.
1#include <fstream>
2#include <iostream>
3
4#include "EmcRawEvent/EmcDigi.h"
5#include "EmcRecEventModel/RecEmcEventModel.h"
6#include "EventModel/EventHeader.h"
7#include "GaudiKernel/ISvcLocator.h"
8#include "GaudiKernel/MsgStream.h"
9#include "GaudiKernel/PropertyMgr.h"
10#include "GaudiKernel/SmartDataPtr.h"
11#include "McTruth/EmcMcHit.h"
12#include "McTruth/McParticle.h"
13#include "RawEvent/RawDataUtil.h"
14
15#include "EmcRec.h"
16#include "EmcRecGeoSvc/IEmcRecGeoSvc.h"
17
18#include "EmcRec/EmcRecCluster2Shower.h"
19#include "EmcRec/EmcRecFastCluster2Shower.h"
20#include "EmcRec/EmcRecParameter.h"
21#include "EmcRec/EmcRecTDS.h"
22#include "EmcRec/EmcRecTofDigitCalib.h"
23#include "EmcRec/EmcRecTofMatch.h"
24
25#include "RawDataProviderSvc/EmcRawDataProvider.h"
26#include "RawDataProviderSvc/IRawDataProviderSvc.h"
27// tianhl for mt
28#include "GaudiKernel/Service.h"
29// #include "GaudiKernel/ThreadGaudi.h"
30// tianhl for mt
31
32using namespace std;
33using namespace Event;
34DECLARE_COMPONENT( EmcRec )
35/////////////////////////////////////////////////////////////////////////////
36EmcRec::EmcRec( const std::string& name, ISvcLocator* pSvcLocator )
37 : Algorithm( name, pSvcLocator ), fCluster2Shower( 0 ) {
38 m_event = 0;
39 fPositionMode.push_back( "log" );
40 fPositionMode.push_back( "5x5" );
41 // Declare the properties
42 declareProperty( "Output", fOutput = 0 );
43 declareProperty( "EventNb", fEventNb = 0 );
44 declareProperty( "DigiCalib", fDigiCalib = false );
45 declareProperty( "TofEnergy", fTofEnergy = false );
46 declareProperty( "OnlineMode", fOnlineMode = false );
47 declareProperty( "TimeMin", fTimeMin = 0 );
48 declareProperty( "TimeMax", fTimeMax = 60 );
49 declareProperty( "PositionMode", fPositionMode );
50
51 // Method == 0 use old correction parameter
52 // Method == 1 use new correction parameter
53 declareProperty( "MethodMode", fMethodMode = 1 );
54 // PosCorr=0, no position correction, and PosCorr=1 with position correction
55 declareProperty( "PosCorr", fPosCorr = 1 );
56 declareProperty( "ElecSaturation", fElecSaturation = 1 );
57}
58
59// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
60StatusCode EmcRec::initialize() {
61
62 MsgStream log( msgSvc(), name() );
63 log << MSG::INFO << "in initialize()" << endmsg;
64
65 // Get EmcRecParameter
66 EmcRecParameter::lock();
67 EmcRecParameter& Para = EmcRecParameter::GetInstance();
68 Para.SetDigiCalib( fDigiCalib );
69 Para.SetTimeMin( fTimeMin );
70 Para.SetTimeMax( fTimeMax );
71 Para.SetMethodMode( fMethodMode );
72 Para.SetPosCorr( fPosCorr );
73 Para.SetPositionMode( fPositionMode );
74 Para.SetElecSaturation( fElecSaturation );
75
76 EmcRecParameter::unlock();
77
78 // Get RawDataProviderSvc
79 // tianhl for mt
80 std::string rawDataProviderSvc_name( "RawDataProviderSvc" );
81 StatusCode sc = service( rawDataProviderSvc_name.c_str(), m_rawDataProviderSvc );
82 if ( sc != StatusCode::SUCCESS )
83 { log << MSG::ERROR << "EmcRec Error: Can't get RawDataProviderSvc." << endmsg; }
84
85 fDigit2Hit.SetAlgName( name() );
86
87#ifndef OnlineMode
88 if ( !( m_rawDataProviderSvc->isOnlineMode() ) )
89 {
90 m_tuple = 0;
91 StatusCode status;
92
93 if ( fOutput >= 1 )
94 {
95 NTuplePtr nt( ntupleSvc(), "FILE301/n1" );
96 if ( nt ) m_tuple = nt;
97 else
98 {
99 m_tuple = ntupleSvc()->book( "FILE301/n1", CLID_ColumnWiseTuple, "EmcRec" );
100 if ( m_tuple )
101 {
102 status = m_tuple->addItem( "pid", pid );
103 status = m_tuple->addItem( "tp", tp );
104 status = m_tuple->addItem( "ttheta", ttheta );
105 status = m_tuple->addItem( "tphi", tphi );
106 status = m_tuple->addItem( "nrun", nrun );
107 status = m_tuple->addItem( "nrec", nrec );
108 status = m_tuple->addItem( "nneu", nneu );
109 status = m_tuple->addItem( "ncluster", ncluster );
110 status = m_tuple->addItem( "npart", npart );
111 status = m_tuple->addItem( "ntheta", ntheta );
112 status = m_tuple->addItem( "nphi", nphi );
113 status = m_tuple->addItem( "ndigi", ndigi );
114 status = m_tuple->addItem( "nhit", nhit );
115 status = m_tuple->addItem( "pp1", 4, pp1 );
116 status = m_tuple->addItem( "theta1", theta1 );
117 status = m_tuple->addItem( "phi1", phi1 );
118 status = m_tuple->addItem( "dphi1", dphi1 );
119 status = m_tuple->addItem( "eseed", eseed );
120 status = m_tuple->addItem( "e3x3", e3x3 );
121 status = m_tuple->addItem( "e5x5", e5x5 );
122 status = m_tuple->addItem( "enseed", enseed );
123 status = m_tuple->addItem( "etof2x1", etof2x1 );
124 status = m_tuple->addItem( "etof2x3", etof2x3 );
125 status = m_tuple->addItem( "cluster2ndMoment", cluster2ndMoment );
126 status = m_tuple->addItem( "secondMoment", secondMoment );
127 status = m_tuple->addItem( "latMoment", latMoment );
128 status = m_tuple->addItem( "a20Moment", a20Moment );
129 status = m_tuple->addItem( "a42Moment", a42Moment );
130 status = m_tuple->addItem( "mpi0", mpi0 );
131 status = m_tuple->addItem( "thtgap1", thtgap1 );
132 status = m_tuple->addItem( "phigap1", phigap1 );
133 status = m_tuple->addItem( "pp2", 4, pp2 );
134 }
135 else
136 { // did not manage to book the N tuple....
137 log << MSG::ERROR << "Cannot book N-tuple:" << long( m_tuple ) << endmsg;
138 return StatusCode::FAILURE;
139 }
140 }
141 }
142 fCluster2Shower = new EmcRecCluster2Shower;
143 }
144 else { fCluster2Shower = new EmcRecFastCluster2Shower; }
145#else
146 fCluster2Shower = new EmcRecFastCluster2Shower;
147#endif
148
149 return StatusCode::SUCCESS;
150}
151
152// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
153StatusCode EmcRec::execute() {
154
155 MsgStream log( msgSvc(), name() );
156 log << MSG::DEBUG << "in execute()" << endmsg;
157 /// Reconstruction begins.
158 /// State 1: Initialize digit map
159 int event, run;
160
161 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
162 if ( !eventHeader )
163 {
164 log << MSG::FATAL << name() << " Could not find Event Header" << endmsg;
165 return ( StatusCode::FAILURE );
166 }
167 run = eventHeader->runNumber();
168 event = eventHeader->eventNumber();
169 EmcRecParameter& Para = EmcRecParameter::GetInstance();
170
171 if ( run > 0 ) Para.SetDataMode( 1 );
172 else Para.SetDataMode( 0 );
173
174 if ( fEventNb != 0 && m_event % fEventNb == 0 )
175 { log << MSG::FATAL << m_event << "-------: " << run << "," << event << endmsg; }
176 m_event++;
177 if ( fOutput >= 2 )
178 {
179 log << MSG::DEBUG << "====================================" << endmsg;
180 log << MSG::DEBUG << "run= " << run << "; event= " << event << endmsg;
181 }
182
183 Hep3Vector posG; // initial photon position
184#ifndef OnlineMode
185 if ( !( m_rawDataProviderSvc->isOnlineMode() ) )
186 {
187 if ( fOutput >= 1 && run < 0 )
188 { // run<0 means MC
189 // Retrieve mc truth
190 SmartDataPtr<McParticleCol> mcParticleCol( eventSvc(), "/Event/MC/McParticleCol" );
191 if ( !mcParticleCol ) { log << MSG::WARNING << "Could not find McParticle" << endmsg; }
192 else
193 {
194 HepLorentzVector pG;
195 McParticleCol::iterator iter_mc = mcParticleCol->begin();
196 for ( ; iter_mc != mcParticleCol->end(); iter_mc++ )
197 {
198 log << MSG::INFO << " particleId = " << ( *iter_mc )->particleProperty() << endmsg;
199 pG = ( *iter_mc )->initialFourMomentum();
200 posG = ( *iter_mc )->initialPosition().v();
201 }
202 ttheta = pG.theta();
203 tphi = pG.phi();
204 if ( tphi < 0 ) { tphi = twopi + tphi; }
205 tp = pG.rho();
206
207 // Retrieve EMC truth
208 SmartDataPtr<EmcMcHitCol> emcMcHitCol( eventSvc(), "/Event/MC/EmcMcHitCol" );
209 if ( !emcMcHitCol ) { log << MSG::WARNING << "Could not find EMC truth" << endmsg; }
210
211 RecEmcID mcId;
212 unsigned int mcTrackIndex;
213 double mcPosX = 0, mcPosY = 0, mcPosZ = 0;
214 double mcPx = 0, mcPy = 0, mcPz = 0;
215 double mcEnergy = 0;
216
217 EmcMcHitCol::iterator iterMc;
218 for ( iterMc = emcMcHitCol->begin(); iterMc != emcMcHitCol->end(); iterMc++ )
219 {
220 mcId = ( *iterMc )->identify();
221 mcTrackIndex = ( *iterMc )->getTrackIndex();
222 mcPosX = ( *iterMc )->getPositionX();
223 mcPosY = ( *iterMc )->getPositionY();
224 mcPosZ = ( *iterMc )->getPositionZ();
225 mcPx = ( *iterMc )->getPx();
226 mcPy = ( *iterMc )->getPy();
227 mcPz = ( *iterMc )->getPz();
228 mcEnergy = ( *iterMc )->getDepositEnergy();
229
230 if ( fOutput >= 2 )
231 {
232 // cout<<"mcId="<<mcId<<"\t"<<"mcTrackIndex="<<mcTrackIndex<<endl;
233 // cout<<"mcPosition:\t"<<mcPosX<<"\t"<<mcPosY<<"\t"<<mcPosZ<<endl;
234 // cout<<"mcP:\t"<<mcPx<<"\t"<<mcPy<<"\t"<<mcPz<<endl;
235 // cout<<"mcDepositEnergy=\t"<<mcEnergy<<endl;
236 }
237 }
238 }
239 } // fOutput>=1
240 } // isOnlineMode
241
242#endif
243 // cout<<"EmcRec1--1"<<endl;
244 /// Retrieve EMC digi
245 fDigitMap.clear();
246 fHitMap.clear();
247 fClusterMap.clear();
248 fShowerMap.clear();
249
250 // Get EmcCalibConstSvc.
251 IEmcCalibConstSvc* emcCalibConstSvc = 0;
252 StatusCode sc = service( "EmcCalibConstSvc", emcCalibConstSvc );
253 if ( sc != StatusCode::SUCCESS )
254 {
255 ;
256 // cout << "EmcRec Error: Can't get EmcCalibConstSvc." << endl;
257 }
258
259 SmartDataPtr<EmcDigiCol> emcDigiCol( eventSvc(), "/Event/Digi/EmcDigiCol" );
260 if ( !emcDigiCol )
261 {
262 log << MSG::FATAL << "Could not find EMC digi" << endmsg;
263 return ( StatusCode::FAILURE );
264 }
265 // cout<<"EmcRec1--2"<<endl;
266 EmcDigiCol::iterator iter3;
267 for ( iter3 = emcDigiCol->begin(); iter3 != emcDigiCol->end(); iter3++ )
268 {
269 RecEmcID id( ( *iter3 )->identify() );
270 double adc =
271 RawDataUtil::EmcCharge( ( *iter3 )->getMeasure(), ( *iter3 )->getChargeChannel() );
272 double tdc = RawDataUtil::EmcTime( ( *iter3 )->getTimeChannel() );
273
274 // for cable connection correction
275 unsigned int nnpart = EmcID::barrel_ec( id );
276 unsigned int nnthe = EmcID::theta_module( id );
277 unsigned int nnphi = EmcID::phi_module( id );
278
279 int index = emcCalibConstSvc->getIndex( nnpart, nnthe, nnphi );
280 int ixtalNumber = emcCalibConstSvc->getIxtalNumber( index );
281
282 if ( run > 0 && ixtalNumber > 0 )
283 {
284 unsigned int npartNew = emcCalibConstSvc->getPartID( ixtalNumber );
285 unsigned int ntheNew = emcCalibConstSvc->getThetaIndex( ixtalNumber );
286 unsigned int nphiNew = emcCalibConstSvc->getPhiIndex( ixtalNumber );
287 id = EmcID::crystal_id( npartNew, ntheNew, nphiNew );
288 } //-------
289
290 //*2016-12
291
292 /*
293 unsigned int nphiNew;
294 nphiNew=nnphi;
295 if(nnpart==0) {
296 if (nnthe==0){
297 if (nnphi==36) nphiNew=40;
298 if (nnphi==40) nphiNew=36;
299 if (nnphi==37) nphiNew=41;
300 if (nnphi==41) nphiNew=37;
301 }//----
302 if (nnthe==1){
303 if (nnphi==36) nphiNew=40;
304 if (nnphi==40) nphiNew=36;
305 if (nnphi==37) nphiNew=41;
306 if (nnphi==41) nphiNew=37;
307 }//----
308 if (nnthe==2){
309 if (nnphi==45) nphiNew=50;
310 if (nnphi==50) nphiNew=45;
311 if (nnphi==46) nphiNew=51;
312 if (nnphi==51) nphiNew=46;
313 if (nnphi==47) nphiNew=52;
314 if (nnphi==52) nphiNew=47;
315 }//----
316 if (nnthe==3){
317 if (nnphi==45) nphiNew=50;
318 if (nnphi==50) nphiNew=45;
319 if (nnphi==46) nphiNew=51;
320 if (nnphi==51) nphiNew=46;
321
322 }//----
323 if (nnthe==4){
324 if (nnphi==54) nphiNew=60;
325 if (nnphi==60) nphiNew=54;
326 if (nnphi==55) nphiNew=61;
327 if (nnphi==61) nphiNew=55;
328 if (nnphi==56) nphiNew=62;
329 if (nnphi==62) nphiNew=56;
330 }//----
331 if (nnthe==5){
332 if (nnphi==54) nphiNew=60;
333 if (nnphi==60) nphiNew=54;
334 if (nnphi==55) nphiNew=61;
335 if (nnphi==61) nphiNew=55;
336 if (nnphi==56) nphiNew=62;
337 if (nnphi==62) nphiNew=56;
338 }//----
339
340 id=EmcID::crystal_id(nnpart,nnthe,nphiNew);
341 }
342 */
343
344 // 2017-12
345 if ( run >= 52940 && run <= 52995 )
346 {
347 unsigned int ntheNew, nphiNew;
348 ntheNew = nnthe;
349 nphiNew = nnphi;
350 if ( nnpart == 2 )
351 {
352 if ( nnthe == 0 )
353 {
354 if ( nnphi == 58 ) nphiNew = 60;
355 if ( nnphi == 59 ) nphiNew = 61;
356 if ( nnphi == 60 ) nphiNew = 58;
357 if ( nnphi == 61 ) nphiNew = 59;
358 ntheNew = 0;
359 } //----
360 if ( nnthe == 1 )
361 {
362 if ( nnphi == 58 ) nphiNew = 60;
363 if ( nnphi == 59 ) nphiNew = 61;
364 if ( nnphi == 60 ) nphiNew = 58;
365 if ( nnphi == 61 ) nphiNew = 59;
366 ntheNew = 1;
367 } //----
368 if ( nnthe == 2 )
369 {
370 if ( nnphi == 73 ) nphiNew = 75;
371 if ( nnphi == 74 ) nphiNew = 76;
372 if ( nnphi == 75 ) nphiNew = 73;
373 if ( nnphi == 76 ) nphiNew = 74;
374 ntheNew = 2;
375
376 } //----
377 if ( nnthe == 3 )
378 {
379 if ( nnphi == 73 ) nphiNew = 75;
380 if ( nnphi == 74 ) nphiNew = 76;
381 if ( nnphi == 75 ) nphiNew = 73;
382 if ( nnphi == 76 ) nphiNew = 74;
383 ntheNew = 3;
384 } //----
385 if ( nnthe == 3 && nnphi == 72 )
386 {
387 ntheNew = 2;
388 nphiNew = 77;
389 }
390 if ( nnthe == 2 && nnphi == 77 )
391 {
392 ntheNew = 3;
393 nphiNew = 72;
394 }
395
396 if ( nnthe == 4 )
397 {
398 if ( nnphi == 87 ) nphiNew = 90;
399 if ( nnphi == 88 ) nphiNew = 91;
400 if ( nnphi == 89 ) nphiNew = 92;
401 if ( nnphi == 90 ) nphiNew = 87;
402 if ( nnphi == 91 ) nphiNew = 88;
403 if ( nnphi == 92 ) nphiNew = 89;
404 ntheNew = 4;
405 } //----
406 if ( nnthe == 5 )
407 {
408 if ( nnphi == 87 ) nphiNew = 90;
409 if ( nnphi == 88 ) nphiNew = 91;
410 if ( nnphi == 89 ) nphiNew = 92;
411 if ( nnphi == 90 ) nphiNew = 87;
412 if ( nnphi == 91 ) nphiNew = 88;
413 if ( nnphi == 92 ) nphiNew = 89;
414 ntheNew = 5;
415 } //----
416
417 id = EmcID::crystal_id( nnpart, ntheNew, nphiNew );
418 }
419 }
420 ///==========
421
422 // ixtalNumber=-9 tag dead channel
423 if ( ixtalNumber == -9 ) { adc = 0.0; }
424
425 // ixtalNumber=-99 tag hot channel
426 if ( ixtalNumber == -99 ) { adc = 0.0; }
427
428 RecEmcDigit aDigit;
429 aDigit.Assign( id, adc, tdc );
430 fDigitMap[id] = aDigit;
431 }
432 // cout<<"EmcRec1--3"<<endl;
433 if ( fOutput >= 2 )
434 {
435 RecEmcDigitMap::iterator iDigitMap;
436 for ( iDigitMap = fDigitMap.begin(); iDigitMap != fDigitMap.end(); iDigitMap++ )
437 {
438 // cout<<iDigitMap->second;
439 }
440 }
441 // DigitMap ok
442
443 // oooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
444 /// State 2: DigitMap --> HitMap
445 fDigit2Hit.Convert( fDigitMap, fHitMap );
446 if ( fOutput >= 2 )
447 {
448 RecEmcHitMap::iterator iHitMap;
449 for ( iHitMap = fHitMap.begin(); iHitMap != fHitMap.end(); iHitMap++ )
450 {
451 // cout<<iHitMap->second;
452 }
453 fDigit2Hit.Output( fHitMap );
454 }
455 // cout<<"EmcRec1--4"<<endl;
456
457 /// State 3: HitMap --> ClusterMap
458 fHit2Cluster.Convert( fHitMap, fClusterMap );
459 //
460 if ( fOutput >= 2 )
461 {
462 RecEmcClusterMap::iterator iClusterMap;
463 for ( iClusterMap = fClusterMap.begin(); iClusterMap != fClusterMap.end(); iClusterMap++ )
464 {
465 // cout<<iClusterMap->second;
466 }
467 }
468 // cout<<"EmcRec1--5"<<endl;
469 /// State 4: ClusterMap --> ShowerMap
470 fCluster2Shower->Convert( fClusterMap, fShowerMap );
471 //
472 if ( fOutput >= 2 )
473 {
474 RecEmcShowerMap::iterator iShowerMap;
475 for ( iShowerMap = fShowerMap.begin(); iShowerMap != fShowerMap.end(); iShowerMap++ )
476 {
477 // cout<<iShowerMap->second;
478 }
479 }
480 // cout<<"EmcRec1--6"<<endl;
481
482 /// State 6: Register to TDS
483 EmcRecTDS tds;
484 tds.RegisterToTDS( fHitMap, fClusterMap, fShowerMap );
485 if ( fOutput >= 2 ) { tds.CheckRegister(); }
486 // cout<<"EmcRec1--7"<<endl;
487 // oooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
488#ifndef OnlineMode
489 if ( !( m_rawDataProviderSvc->isOnlineMode() ) )
490 {
491 if ( fOutput >= 1 )
492 {
493 nrun = run;
494 nrec = event;
495
496 // cout<<"cm="<<cm<<"\tmm="<<mm<<"\tGeV="<<GeV<<"\tMeV="<<MeV<<endl;
497 ndigi = fDigitMap.size();
498 nhit = fHitMap.size();
499 ncluster = fClusterMap.size();
500 RecEmcShowerVec fShowerVec;
501 RecEmcShowerMap::iterator iShowerMap;
502 for ( iShowerMap = fShowerMap.begin(); iShowerMap != fShowerMap.end(); iShowerMap++ )
503 { fShowerVec.push_back( iShowerMap->second ); }
504 sort( fShowerVec.begin(), fShowerVec.end(), greater<RecEmcShower>() );
505 nneu = fShowerMap.size();
506
507 RecEmcShower aShower;
508 // aShower.e5x5(-99.);
509 RecEmcShowerVec::iterator iShowerVec;
510 iShowerVec = fShowerVec.begin();
511 npart = -99;
512 ntheta = -99;
513 nphi = -99;
514 if ( iShowerVec != fShowerVec.end() )
515 {
516 aShower = *iShowerVec;
517 RecEmcID id = aShower.getShowerId();
518 npart = EmcID::barrel_ec( id );
519 ntheta = EmcID::theta_module( id );
520 nphi = EmcID::phi_module( id );
521 pp1[0] = aShower.energy() * aShower.x() / aShower.position().mag();
522 pp1[1] = aShower.energy() * aShower.y() / aShower.position().mag();
523 pp1[2] = aShower.energy() * aShower.z() / aShower.position().mag();
524 pp1[3] = aShower.energy();
525 theta1 = ( aShower.position() - posG ).theta();
526 phi1 = ( aShower.position() - posG ).phi();
527 if ( phi1 < 0 ) { phi1 = twopi + phi1; }
528 thtgap1 = aShower.ThetaGap();
529 phigap1 = aShower.PhiGap();
530 RecEmcID seed, nseed;
531 seed = aShower.getShowerId();
532 nseed = aShower.NearestSeed();
533 eseed = aShower.eSeed();
534 e3x3 = aShower.e3x3();
535 e5x5 = aShower.e5x5();
536 etof2x1 = aShower.getETof2x1();
537 etof2x3 = aShower.getETof2x3();
538 if ( aShower.getCluster() )
539 { cluster2ndMoment = aShower.getCluster()->getSecondMoment(); }
540 secondMoment = aShower.secondMoment();
541 latMoment = aShower.latMoment();
542 a20Moment = aShower.a20Moment();
543 a42Moment = aShower.a42Moment();
544 if ( nseed.is_valid() ) { enseed = fHitMap[nseed].getEnergy(); }
545 else { enseed = 0; }
546
547 dphi1 = phi1 - tphi;
548 if ( dphi1 < -pi ) dphi1 = dphi1 + twopi;
549 if ( dphi1 > pi ) dphi1 = dphi1 - twopi;
550 }
551
552 if ( iShowerVec != fShowerVec.end() )
553 {
554 iShowerVec++;
555 if ( iShowerVec != fShowerVec.end() )
556 {
557 aShower = *iShowerVec;
558 pp2[0] = aShower.energy() * aShower.x() / aShower.position().mag();
559 pp2[1] = aShower.energy() * aShower.y() / aShower.position().mag();
560 pp2[2] = aShower.energy() * aShower.z() / aShower.position().mag();
561 pp2[3] = aShower.energy();
562 }
563 }
564
565 // for pi0
566 if ( fShowerVec.size() >= 2 )
567 {
568 RecEmcShowerVec::iterator iShowerVec1, iShowerVec2;
569 iShowerVec1 = fShowerVec.begin();
570 iShowerVec2 = fShowerVec.begin() + 1;
571 double e1 = ( *iShowerVec1 ).energy();
572 double e2 = ( *iShowerVec2 ).energy();
573 double angle = ( *iShowerVec1 ).position().angle( ( *iShowerVec2 ).position() );
574 mpi0 = sqrt( 2 * e1 * e2 * ( 1 - cos( angle ) ) );
575 }
576 m_tuple->write();
577 }
578 }
579#endif
580
581 return StatusCode::SUCCESS;
582}
583
584// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
585StatusCode EmcRec::finalize() {
586 EmcRecParameter::Kill();
587 if ( fCluster2Shower ) delete fCluster2Shower;
588
589 MsgStream log( msgSvc(), name() );
590 log << MSG::INFO << "in finalize()" << endmsg;
591 return StatusCode::SUCCESS;
592}
DECLARE_COMPONENT
DECLARE_COMPONENT(BesBdkRc)
e1
Double_t e1
Definition DataBase/tau_mode.c:7
e2
Double_t e2
Definition DataBase/tau_mode.c:7
pi
double pi
Definition EvtPhokharaDef.hh:17
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
ntupleSvc
INTupleSvc * ntupleSvc()
Definition ServiceAccessor.h:51
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:12
EmcID::crystal_id
static Identifier crystal_id(const unsigned int barrel_ec, const unsigned int theta_module, const unsigned int phi_module)
For a single crystal.
Definition EmcID.cxx:63
EmcID::barrel_ec
static unsigned int barrel_ec(const Identifier &id)
Values of different levels (failure returns 0).
Definition EmcID.cxx:36
EmcID::theta_module
static unsigned int theta_module(const Identifier &id)
Definition EmcID.cxx:41
EmcID::phi_module
static unsigned int phi_module(const Identifier &id)
Definition EmcID.cxx:46
EmcRecParameter::SetPositionMode
void SetPositionMode(std::vector< std::string > &mode)
Definition EmcRecParameter.cxx:661
EmcRecParameter::GetInstance
static EmcRecParameter & GetInstance()
Definition EmcRecParameter.cxx:494
EmcRecTDS::RegisterToTDS
StatusCode RegisterToTDS(RecEmcHitMap &aHitMap, RecEmcClusterMap &aClusterMap, RecEmcShowerMap &aShowerMap)
EmcRecTDS::CheckRegister
StatusCode CheckRegister()
Definition EmcRecTDS.cxx:204
EmcRec::initialize
StatusCode initialize()
Definition EmcRec.cxx:60
EmcRec::finalize
StatusCode finalize()
Definition EmcRec.cxx:585
EmcRec::EmcRec
EmcRec(const std::string &name, ISvcLocator *pSvcLocator)
Definition EmcRec.cxx:36
EmcRec::execute
StatusCode execute()
Definition EmcRec.cxx:153
IEmcCalibConstSvc::getPartID
virtual unsigned int getPartID(int Index) const =0
IEmcCalibConstSvc::getIxtalNumber
virtual int getIxtalNumber(int No) const =0
IEmcCalibConstSvc::getPhiIndex
virtual unsigned int getPhiIndex(int Index) const =0
IEmcCalibConstSvc::getThetaIndex
virtual unsigned int getThetaIndex(int Index) const =0
IEmcCalibConstSvc::getIndex
virtual int getIndex(unsigned int PartId, unsigned int ThetaIndex, unsigned int PhiIndex) const =0
Identifier::is_valid
bool is_valid() const
Check if id is in a valid state.
Definition DetectorDescription/Identifier/include/Identifier/Identifier.h:159
RawDataUtil::EmcTime
static double EmcTime(int timeChannel)
Definition Event/RawEvent/include/RawEvent/RawDataUtil.h:17
RawDataUtil::EmcCharge
static double EmcCharge(int measure, int chargeChannel)
Definition Event/RawEvent/include/RawEvent/RawDataUtil.h:22
RecEmcCluster::getSecondMoment
double getSecondMoment() const
Definition RecEmcCluster.cxx:207
RecEmcDigit::Assign
void Assign(const RecEmcID &CellId, const RecEmcADC &ADC, const RecEmcTDC &TDC)
Definition RecEmcDigit.cxx:70
RecEmcShower::ThetaGap
int ThetaGap() const
Definition RecEmcShower.cxx:240
RecEmcShower::NearestSeed
RecEmcID NearestSeed() const
Definition RecEmcShower.cxx:236
RecEmcShower::PhiGap
int PhiGap() const
Definition RecEmcShower.cxx:244