BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EventAssemblyAlg.cxx
Go to the documentation of this file.
1#include "EventAssemblyAlg.h"
2
3#include "GaudiKernel/IDataManagerSvc.h"
4#include "GaudiKernel/IDataProviderSvc.h"
5#include "GaudiKernel/ISvcLocator.h"
6#include "GaudiKernel/MsgStream.h"
7#include "GaudiKernel/PropertyMgr.h"
8#include "GaudiKernel/SmartDataPtr.h"
9// #include "GaudiKernel/INTupleSvc.h"
10// #include "GaudiKernel/IHistogramSvc.h"
11
12#include "EvTimeEvent/RecEsTime.h"
13#include "EventModel/EventHeader.h"
14#include "EvtRecEvent/EvtRecEvent.h"
15#include "EvtRecEvent/EvtRecObject.h"
16#include "EvtRecEvent/EvtRecTrack.h"
17#include "Identifier/Identifier.h"
18#include "Identifier/TofID.h"
19
20DECLARE_COMPONENT( EventAssemblyAlg )
21/////////////////////////////////////////////////////////////////////////////
22EventAssemblyAlg::EventAssemblyAlg( const std::string& name, ISvcLocator* pSvcLocator )
23 : Algorithm( name, pSvcLocator ) {
24 // Declare the properties
25 declareProperty( "MsgFlag", myMsgFlag = false );
26 declareProperty( "ActiveCutFlag", myActiveCutFlag = true );
27 declareProperty( "EmcThetaMatchCut", myEmcThetaCut = -1 );
28 declareProperty( "EmcPhiMatchCut", myEmcPhiCut = -1 );
29 declareProperty( "EmcThetaNSigmaCut", myEmcThetaNSigmaCut = 5 );
30 declareProperty( "EmcPhiNSigmaCut", myEmcPhiNSigmaCut = 5 );
31 declareProperty( "ParticleId", myParticleId );
32 declareProperty( "Output", m_Output = 0 );
33}
34
35//////////////////////////////////////////
36StatusCode EventAssemblyAlg::initialize() {
37 MsgStream log( msgSvc(), name() );
38 log << MSG::INFO << "in initialize()" << endmsg;
39
40 if ( myEmcThetaCut < 0 ) myEmcThetaCut = 0.1745;
41 if ( myEmcPhiCut < 0 ) myEmcPhiCut = 0.2618;
42
43 log << MSG::INFO << "EmcThetaCut,EmcPhiCut: " << myEmcThetaCut << ", " << myEmcPhiCut
44 << endmsg;
45 log << MSG::INFO << "emcThetaNSigmaCut, emcPhiNSigmaCut: " << myEmcThetaNSigmaCut << ", "
46 << myEmcPhiNSigmaCut << endmsg;
47
48 if ( m_Output == 1 )
49 {
50 NTuplePtr nt1( ntupleSvc(), "FILE998/match" );
51 if ( nt1 ) myNTuple1 = nt1;
52 else
53 {
54 myNTuple1 = ntupleSvc()->book( "FILE998/match", CLID_ColumnWiseTuple, "match" );
55 if ( myNTuple1 )
56 {
57 myNTuple1->addItem( "exted", myExted );
58 myNTuple1->addItem( "matched", myMatched );
59 myNTuple1->addItem( "Ematched", myEnergyMatched );
60 }
61 else
62 {
63 log << MSG::ERROR << " Cannot book N-tuple:" << long( myNTuple1 ) << endmsg;
64 // return StatusCode::FAILURE;
65 }
66 }
67 }
68 return StatusCode::SUCCESS;
69}
70
71///////////////////////////////////////
72StatusCode EventAssemblyAlg::execute() {
73 MsgStream log( msgSvc(), name() );
74 log << MSG::INFO << "in execute()" << endmsg;
75
76 int numOfChargedTrks = 0;
77 int numOfKalTrks = 0;
78 int numOfDedx = 0;
79 int numOfExt = 0;
80 int numOfTof = 0;
81 int numOfShower = 0;
82 int numOfMatchedShower = 0;
83 int numOfMucTrks = 0;
84
85 // Part 1: Get the event header, print out event and run number
86 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
87 if ( !eventHeader )
88 {
89 log << MSG::FATAL << "Could not find Event Header" << endmsg;
90 return ( StatusCode::FAILURE );
91 }
92 log << MSG::INFO << ": retrieved event: " << eventHeader->eventNumber()
93 << " run: " << eventHeader->runNumber() << endmsg;
94
95 // Part 2: Retrieve Rec data
96 SmartDataPtr<RecMdcTrackCol> recMdcTrackCol( eventSvc(), EventModel::Recon::RecMdcTrackCol );
97 RecMdcTrackCol::iterator beginOfMdcTrkCol = RecMdcTrackCol::iterator( NULL );
98 if ( !recMdcTrackCol ) { log << MSG::INFO << "Could not find RecMdcTrackCol!" << endmsg; }
99 else
100 {
101 beginOfMdcTrkCol = recMdcTrackCol->begin();
102 numOfChargedTrks = recMdcTrackCol->size();
103 }
104
105 SmartDataPtr<RecMdcKalTrackCol> mdcKalTrackCol( eventSvc(),
106 EventModel::Recon::RecMdcKalTrackCol );
107 RecMdcKalTrackCol::iterator beginOfMdcKalTrackCol = RecMdcKalTrackCol::iterator( NULL );
108 if ( !mdcKalTrackCol ) { log << MSG::INFO << "Could not find MdcKalTrackCol!" << endmsg; }
109 else
110 {
111 beginOfMdcKalTrackCol = mdcKalTrackCol->begin();
112 numOfKalTrks = mdcKalTrackCol->size();
113 }
114
115 SmartDataPtr<RecMdcDedxCol> mdcDedxCol( eventSvc(), EventModel::Recon::RecMdcDedxCol );
116 RecMdcDedxCol::iterator beginOfDedxCol = RecMdcDedxCol::iterator( NULL );
117 if ( !mdcDedxCol ) { log << MSG::INFO << "Could not find MdcDedxCol!" << endmsg; }
118 else
119 {
120 beginOfDedxCol = mdcDedxCol->begin();
121 numOfDedx = mdcDedxCol->size();
122 }
123
124 SmartDataPtr<RecExtTrackCol> extTrackCol( eventSvc(), EventModel::Recon::RecExtTrackCol );
125 RecExtTrackCol::iterator beginOfExtTrackCol = RecExtTrackCol::iterator( NULL );
126 if ( !extTrackCol ) { log << MSG::INFO << "Could not find RecExtTrackCol!" << endmsg; }
127 else
128 {
129 beginOfExtTrackCol = extTrackCol->begin();
130 numOfExt = extTrackCol->size();
131 }
132
133 SmartDataPtr<RecTofTrackCol> TofTrackCol( eventSvc(), EventModel::Recon::RecTofTrackCol );
134 RecTofTrackCol::iterator beginOfTofTrackCol = RecTofTrackCol::iterator( NULL );
135 if ( !TofTrackCol ) { log << MSG::INFO << "Could not find RecTofTrackCol!" << endmsg; }
136 else
137 {
138 beginOfTofTrackCol = TofTrackCol->begin();
139 numOfTof = TofTrackCol->size();
140 }
141
142 SmartDataPtr<RecEmcShowerCol> emcRecShowerCol( eventSvc(),
143 EventModel::Recon::RecEmcShowerCol );
144 RecEmcShowerCol::iterator beginOfEmcTrackCol = RecEmcShowerCol::iterator( NULL );
145 if ( !emcRecShowerCol ) { log << MSG::INFO << "Could not find RecEmcShowerCol!" << endmsg; }
146 else
147 {
148 beginOfEmcTrackCol = emcRecShowerCol->begin();
149 numOfShower = emcRecShowerCol->size();
150 }
151
152 SmartDataPtr<RecMucTrackCol> mucTrackCol( eventSvc(), EventModel::Recon::RecMucTrackCol );
153 RecMucTrackCol::iterator beginOfMucTrackCol = RecMucTrackCol::iterator( NULL );
154 if ( !mucTrackCol ) { log << MSG::INFO << "Could not find RecMucTrackCol!" << endmsg; }
155 else
156 {
157 beginOfMucTrackCol = mucTrackCol->begin();
158 numOfMucTrks = mucTrackCol->size();
159 }
160
161 // Some match flag
162
163 int MaxHit = numOfChargedTrks + numOfShower + numOfTof;
164
165 vector<bool> dEdxMatched;
166 vector<bool> kalTrkMatched;
167 vector<bool> extMatched;
168 vector<bool> TofMatched;
169 vector<bool> emcMatched;
170 vector<bool> mucMatched;
171
172 for ( int i = 0; i < MaxHit; i++ )
173 {
174 dEdxMatched.push_back( false );
175 kalTrkMatched.push_back( false );
176 extMatched.push_back( false );
177 TofMatched.push_back( false );
178 emcMatched.push_back( false );
179 mucMatched.push_back( false );
180 }
181
182 // Part 3: Match tracks
183 // RecTrackListCol* aNewRecTrkListCol = new RecTrackListCol;
184 EvtRecTrackCol* aNewEvtRecTrackCol = new EvtRecTrackCol;
185 if ( numOfChargedTrks + numOfShower ) // if mdc trackCol and emc trackCol neither
186 // exists,there is no match job.
187 {
188 for ( int i = 0; i < numOfChargedTrks; i++ ) // match charged tracks
189 {
190 EvtRecTrack* aNewEvtRecTrack = new EvtRecTrack;
191 int mdcTrkID = ( *( beginOfMdcTrkCol + i ) )->trackId();
192 aNewEvtRecTrack->setTrackId( mdcTrkID );
193 aNewEvtRecTrack->setMdcTrack( *( beginOfMdcTrkCol + i ) );
194 // aNewEvtRecTrack->setMdcTrkIdx(i);
195 for ( int j = 0; j < numOfKalTrks; j++ ) // match MdcKalTrack
196 if ( !kalTrkMatched[j] && mdcTrkID == ( *( beginOfMdcKalTrackCol + j ) )->trackId() )
197 {
198 aNewEvtRecTrack->setMdcKalTrack( *( beginOfMdcKalTrackCol + j ) );
199 // aNewRecTrkList->setMdcKalTrkIdx(j);
200 kalTrkMatched[j] = true;
201 }
202
203 for ( int j = 0; j < numOfDedx; j++ ) // match dEdx
204 if ( !dEdxMatched[j] && mdcTrkID == ( *( beginOfDedxCol + j ) )->trackId() )
205 {
206 aNewEvtRecTrack->setMdcDedx( *( beginOfDedxCol + j ) );
207 // aNewRecTrkList->setDedxTrkIdx(j);
208 dEdxMatched[j] = true;
209 }
210
211 for ( int j = 0; j < numOfExt; j++ ) // match ExtTrack
212 if ( !extMatched[j] && mdcTrkID == ( *( beginOfExtTrackCol + j ) )->trackId() )
213 {
214 aNewEvtRecTrack->setExtTrack( *( beginOfExtTrackCol + j ) );
215 // aNewRecTrkList->setExtTrkIdx(j);
216 extMatched[j] = true;
217 }
218 for ( int j = 0; j < numOfTof; j++ ) // match Tof
219 if ( !TofMatched[j] && mdcTrkID == ( *( beginOfTofTrackCol + j ) )->trackID() )
220 {
221 aNewEvtRecTrack->addTofTrack( *( beginOfTofTrackCol + j ) );
222 // aNewRecTrkList->setTofTrkIdx(j);
223 TofMatched[j] = true;
224 }
225 for ( int j = 0; j < numOfMucTrks; j++ ) // match MucTrack
226 if ( !mucMatched[j] && mdcTrkID == ( *( beginOfMucTrackCol + j ) )->trackId() )
227 {
228 aNewEvtRecTrack->setMucTrack( *( beginOfMucTrackCol + j ) );
229 // aNewRecTrkList->setMucTrkIdx(j);
230 mucMatched[j] = true;
231 }
232
233 if ( m_Output == 1 )
234 {
235 myExted = -1.;
236 myMatched = -1.;
237 myEnergyMatched = -1.;
238 }
239 if ( aNewEvtRecTrack->isExtTrackValid() &&
240 ( aNewEvtRecTrack->extTrack() )->emcVolumeNumber() != -1 ) // start to match Emc
241 // showers
242 {
243 if ( m_Output == 1 ) myExted = 1.;
244 // resieve ext information @ Emc
245 Hep3Vector extPos = ( aNewEvtRecTrack->extTrack() )->emcPosition();
246 double extTheta = extPos.theta();
247 double extPhi = extPos.phi();
248
249 double dAngleMin = 9999.;
250 int indexOfEmcNearest = -1;
251
252 for ( int j = 0; j < numOfShower; j++ ) // start to find nearest emc shower
253 {
254 if ( !emcMatched[j] )
255 {
256 HepPoint3D emcPot = ( *( beginOfEmcTrackCol + j ) )->position();
257 Hep3Vector emcPos( emcPot.x(), emcPot.y(), emcPot.z() );
258 double dAngle = extPos.angle( emcPos );
259 if ( dAngle < dAngleMin )
260 {
261 dAngleMin = dAngle;
262 indexOfEmcNearest = j;
263 }
264 }
265 }
266 if ( indexOfEmcNearest >= 0 ) // if nearest Emc shower is found
267 {
268 HepPoint3D emcPot = ( *( beginOfEmcTrackCol + indexOfEmcNearest ) )->position();
269 double theta = emcPot.theta();
270 double phi = emcPot.phi();
271 double deltaTheta = ( extTheta - theta );
272 double deltaPhi = fmod( extPhi - phi + 5 * M_PI, 2 * M_PI ) - M_PI;
273 if ( myActiveCutFlag )
274 {
275 double tanLamda = ( *( beginOfMdcTrkCol + i ) )->helix( 4 );
276 double kappa = ( *( beginOfMdcTrkCol + i ) )->helix( 2 );
277 double p = sqrt( 1 + tanLamda * tanLamda ) / fabs( kappa );
278 double aThetaCut = thetaCut( p, myEmcThetaNSigmaCut, myParticleId );
279 double aPhiCUt = phiCut( p, myEmcPhiNSigmaCut, myParticleId );
280 if ( fabs( deltaTheta ) < aThetaCut && fabs( deltaPhi ) < aPhiCUt )
281 {
282 aNewEvtRecTrack->setEmcShower( *( beginOfEmcTrackCol + indexOfEmcNearest ) );
283 // aNewRecTrkList->setEmcTrkIdx(indexOfEmcNearest);
284 emcMatched[indexOfEmcNearest] = true;
285 numOfMatchedShower++;
286 if ( m_Output == 1 ) myMatched = 1.;
287 }
288 }
289 else if ( fabs( deltaTheta ) < myEmcThetaCut &&
290 fabs( deltaPhi ) < myEmcPhiCut ) // if the nearest shower is
291 // in the fixed match
292 // window,it will be matched.
293 {
294 aNewEvtRecTrack->setEmcShower( *( beginOfEmcTrackCol + indexOfEmcNearest ) );
295 // aNewRecTrkList->setEmcTrkIdx(indexOfEmcNearest);
296 emcMatched[indexOfEmcNearest] = true;
297 numOfMatchedShower++;
298 if ( m_Output == 1 ) myMatched = 1.;
299 }
300 }
301 }
302 aNewEvtRecTrackCol->push_back( aNewEvtRecTrack );
303 if ( m_Output == 1 )
304 {
305 myEnergyMatched = 0.;
306 myNTuple1->write();
307 }
308 } // end charged track match
309
310 int id = 0;
311 for ( int i = 0; i < numOfShower;
312 i++ ) // If Emc shower isn't macthed, it will be added to the end of RecTrackListCol
313 {
314 if ( !emcMatched[i] )
315 {
316 EvtRecTrack* aNewEvtRecTrack = new EvtRecTrack;
317 aNewEvtRecTrack->setTrackId( id + numOfChargedTrks );
318 aNewEvtRecTrack->setEmcShower( *( beginOfEmcTrackCol + i ) );
319 // aNewRecTrkList->setEmcTrkIdx(i);
320
321 // EMC and TOF match
322 HepPoint3D emcPos( ( *( beginOfEmcTrackCol + i ) )->position() );
323 double emcPhi = emcPos.phi();
324 emcPhi = emcPhi < 0 ? emcPhi + CLHEP::twopi : emcPhi;
325
326 int module = 9999, nphi1 = 9999, nphi2 = 9999, nphi1_mrpc = 9999;
327 module = ( *( beginOfEmcTrackCol + i ) )->module();
328 if ( module == 1 )
329 { // barrel
330 nphi1 = (int)( emcPhi * 88. / CLHEP::twopi ); // tof layer1 number
331 nphi2 = (int)( emcPhi * 88. / CLHEP::twopi + 0.5 ); // tof layer2 number
332 nphi2 += 88;
333 }
334 else
335 { // endcap
336 nphi1 = (int)( emcPhi * 48. / CLHEP::twopi + 0.5 );
337
338 // new mrpc detector
339 double fi_endtof_mrpc = atan2( emcPos.y(), emcPos.x() ) - 3.141592654 / 2.;
340 if ( fi_endtof_mrpc < 0 ) fi_endtof_mrpc += 2 * 3.141592654;
341 nphi1_mrpc = ( (int)( fi_endtof_mrpc / ( 3.141593 / 18 ) ) ) +
342 1; // We have 36 modules in each endcap (divided again into two layers)!
343 // The numbering starts with 1.
344
345 if ( nphi1_mrpc % 2 == 0 )
346 {
347 nphi1_mrpc = nphi1_mrpc / 2;
348 } // Calculate the correct number as in each endcap we have two layers of modules
349 // furnished with numbers 1 to 18.
350 else { nphi1_mrpc = ( nphi1_mrpc + 1 ) / 2; }
351
352 if ( emcPos.z() > 0 )
353 {
354 ( nphi1_mrpc -= 19 );
355 nphi1_mrpc = nphi1_mrpc * ( -1 );
356 } // Consider that for the east endcap the module at x=0;y=max has the number 18 (the
357 // numbering is != to west endcap)
358 // nphi1_mrpc is now the number of the nearest ToF module!
359 }
360
361 int dPhiMin = 9999;
362 int partid_dPhiMin = 7;
363 bool mrpcmatch = false;
364
365 int indexOfTofNearest = -1;
366 for ( int j = 0; j < numOfTof; j++ )
367 { // start to find nearest Tof shower
368 if ( !TofMatched[j] )
369 {
370
371 Identifier id( ( *( beginOfTofTrackCol + j ) )->tofID() );
372 bool is_mrpc = TofID::is_mrpc( id );
373
374 if ( !is_mrpc )
375 {
376 int barrel_ec = TofID::barrel_ec( id );
377 int layer = TofID::layer( id );
378 int im = TofID::phi_module( id );
379 im += layer * 88;
380
381 if ( module != barrel_ec ) continue;
382
383 int dPhi;
384 if ( layer == 0 ) { dPhi = abs( im - nphi1 ); }
385 else { dPhi = abs( im - nphi2 ); }
386 if ( module == 1 ) { dPhi = dPhi >= 44 ? 88 - dPhi : dPhi; }
387 else { dPhi = dPhi >= 24 ? 48 - dPhi : dPhi; }
388
389 if ( dPhi < dPhiMin )
390 {
391 dPhiMin = dPhi;
392 indexOfTofNearest = j;
393 }
394 } // close if(!is_mrpc)
395 else // is_mrpc
396 {
397 int barrel_ec = TofID::barrel_ec( id );
398 int module_mrpc = TofID::module( id );
399
400 int dphi;
401 dphi = abs( module_mrpc - nphi1_mrpc );
402
403 if ( dphi < dPhiMin )
404 {
405 dPhiMin = dphi;
406 partid_dPhiMin = barrel_ec;
407 indexOfTofNearest = j;
408 mrpcmatch = true;
409 }
410 } // close else is_mrpc
411
412 } // close if( !TofMatched[j] )
413 } // close for(int j=0;j<numOfTof;j++)
414
415 /*
416 if(0)
417 {
418 cout << "Event Assembly Algorithm" <<endl;
419 cout << "EMC Shower number = " << i << endl;
420 cout << "mrpcmatch ? " << mrpcmatch << endl;
421 cout << "indexOfTofNearest " << indexOfTofNearest <<endl;
422 cout << "dphiMin " << dPhiMin << endl;
423 cout << "partid_mrpc " << partid_dPhiMin <<endl;
424 cout << "partid_emc " << module <<endl <<endl;
425 }
426 */
427
428 // match in z-direction
429 if ( indexOfTofNearest >= 0 )
430 {
431 HepPoint3D tofPos( 0, 0.867,
432 ( *( beginOfTofTrackCol + indexOfTofNearest ) )->zrhit() );
433 double matWindow = 0;
434 double eShower = ( *( beginOfEmcTrackCol + i ) )->e5x5(); // energy of a 5x5 cluster
435 if ( eShower > 0 ) { matWindow = 0.01277 + 0.004268 / sqrt( eShower ); }
436 }
437
438 if ( indexOfTofNearest >= 0 && dPhiMin <= 1 )
439 { // if nearest Tof shower is found
440 // calculate matching window for barrel
441 double eShower = ( *( beginOfEmcTrackCol + i ) )->e5x5();
442 double matWindow = 0;
443 if ( eShower > 0 ) { matWindow = 0.01277 + 0.004268 / sqrt( eShower ); }
444 matWindow = 0.2; // temp code
445
446 HepPoint3D tofPos( 0, 0.867,
447 ( *( beginOfTofTrackCol + indexOfTofNearest ) )->zrhit() );
448 if ( fabs( tofPos.theta() - emcPos.theta() ) < matWindow || module != 1 )
449 {
450 aNewEvtRecTrack->addTofTrack( *( beginOfTofTrackCol + indexOfTofNearest ) );
451 TofMatched[indexOfTofNearest] = true;
452 // cout << "EventAssembly : New neutral Toftrack added " << endl;
453 }
454 } // close if(indexOfTofNearest>=0 && dPhiMin<=1)
455
456 aNewEvtRecTrackCol->push_back( aNewEvtRecTrack );
457 id++;
458 } // close if(!)
459 } // close for(int i=0;i<numOfShower;i++) <---This are the EMC showers
460 } // end match
461
462 // match statistics
463 if ( myMsgFlag )
464 cout << "Charged tracks: " << numOfChargedTrks << " Ext track: " << numOfExt
465 << " Shower:" << numOfShower << " Matched shower:" << numOfMatchedShower << endl;
466
467 // Part 4: register RecTrackListCol and EventList in TDS
468 DataObject* aDataObject1;
469 eventSvc()->findObject( EventModel::EvtRec::Event, aDataObject1 );
470 if ( aDataObject1 == NULL )
471 {
472 EvtRecObject* aRecObject = new EvtRecObject;
473 StatusCode sc = eventSvc()->registerObject( EventModel::EvtRec::Event, aRecObject );
474 if ( sc.isFailure() )
475 {
476 log << MSG::FATAL << "Could not register EvtRecObject in TDS!" << endmsg;
477 return StatusCode::FAILURE;
478 }
479 }
480
481 DataObject* aDataObject2;
482 eventSvc()->findObject( EventModel::EvtRec::EvtRecEvent, aDataObject2 );
483 if ( aDataObject2 == NULL )
484 {
485 EvtRecEvent* aRecEvent = new EvtRecEvent;
486 aRecEvent->setTotalTracks( numOfChargedTrks + numOfShower - numOfMatchedShower );
487 aRecEvent->setTotalCharged( numOfChargedTrks );
488 aRecEvent->setTotalNeutral( numOfShower - numOfMatchedShower );
489 StatusCode sc = eventSvc()->registerObject( EventModel::EvtRec::EvtRecEvent, aRecEvent );
490 if ( sc.isFailure() )
491 {
492 log << MSG::FATAL << "Could not register EvtRecEvent in TDS!" << endmsg;
493 return ( StatusCode::FAILURE );
494 }
495 }
496 else
497 {
498 EvtRecEvent* aRecEvent = dynamic_cast<EvtRecEvent*>( aDataObject2 );
499 aRecEvent->setTotalTracks( numOfChargedTrks + numOfShower - numOfMatchedShower );
500 aRecEvent->setTotalCharged( numOfChargedTrks );
501 aRecEvent->setTotalNeutral( numOfShower - numOfMatchedShower );
502 }
503
504 DataObject* aDataObject3;
505 eventSvc()->findObject( EventModel::EvtRec::EvtRecTrackCol, aDataObject3 );
506 if ( aDataObject3 != NULL )
507 {
508 // IDataManagerSvc* dataMgrSvc = dynamic_cast<IDataManagerSvc*> (eventSvc());
509 SmartIF<IDataManagerSvc> dataMgrSvc( eventSvc() );
510 dataMgrSvc->clearSubTree( EventModel::EvtRec::EvtRecTrackCol );
511 eventSvc()->unregisterObject( EventModel::EvtRec::EvtRecTrackCol );
512 }
513
514 StatusCode sc =
515 eventSvc()->registerObject( EventModel::EvtRec::EvtRecTrackCol, aNewEvtRecTrackCol );
516 if ( sc.isFailure() )
517 {
518 log << MSG::FATAL << "Could not register EvtRecTrackCol in TDS!" << endmsg;
519 return ( StatusCode::FAILURE );
520 }
521
522 // Part 5: check RecTrackListCol in TDS
523 SmartDataPtr<EvtRecTrackCol> evtRecTrackCol( eventSvc(),
524 EventModel::EvtRec::EvtRecTrackCol );
525 if ( !evtRecTrackCol )
526 { log << MSG::FATAL << "Could not find RecTrackListCol in TDS!" << endmsg; }
527
528 EvtRecTrackCol::iterator itOfRecTrkListCol = evtRecTrackCol->begin();
529 if ( myMsgFlag )
530 for ( int i = 0; itOfRecTrkListCol != evtRecTrackCol->end(); itOfRecTrkListCol++, i++ )
531 {
532 cout << "******************** Track " << i << ": *********************" << endl;
533 cout << "Track ID: " << ( *itOfRecTrkListCol )->trackId() << endl;
534 if ( ( *itOfRecTrkListCol )->isMdcTrackValid() )
535 {
536 RecMdcTrack* mdcTrk = ( *itOfRecTrkListCol )->mdcTrack();
537 double kappa = mdcTrk->helix( 2 );
538 double tanl = mdcTrk->helix( 4 );
539 cout << "Mdc kappa, tanl: " << kappa << ", " << tanl << endl;
540 }
541 if ( ( *itOfRecTrkListCol )->isExtTrackValid() )
542 {
543 RecExtTrack* extTrack = ( *itOfRecTrkListCol )->extTrack();
544 Hep3Vector emcPos = extTrack->emcPosition();
545 cout << "Ext Emc pos:" << emcPos.x() << "," << emcPos.y() << "," << emcPos.z() << endl;
546 }
547 if ( ( *itOfRecTrkListCol )->isEmcShowerValid() )
548 {
549 RecEmcShower* emcTrack = ( *itOfRecTrkListCol )->emcShower();
550 HepPoint3D emcPos = emcTrack->position();
551 double x = emcPos.x();
552 double y = emcPos.y();
553 double z = emcPos.z();
554 cout << "Emc rec pos:" << x << "," << y << "," << z << endl;
555 }
556 }
557
558 return StatusCode::SUCCESS;
559}
560
561// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
562StatusCode EventAssemblyAlg::finalize() {
563 MsgStream log( msgSvc(), name() );
564 log << MSG::INFO << "in finalize()" << endmsg;
565
566 return StatusCode::SUCCESS;
567}
568
569//////////////////////////////
570double EventAssemblyAlg::thetaCut( double p, double nSigma, int parId ) {
571 if ( p < 0.3 ) p = 0.3;
572 if ( parId < 1 || parId > 3 ) parId = 3;
573 if ( nSigma <= 0 ) nSigma = 5.;
574 double cut = myEmcThetaCut;
575 switch ( parId )
576 {
577 case 1: cut = fabs( -0.00159 ) + ( exp( -3.566 * p - 4.14 ) + 0.006804 ) * nSigma; break;
578 case 2: cut = fabs( -0.00145 ) + ( exp( -4.268 * p - 3.305 ) + 0.01129 ) * nSigma; break;
579 case 3: cut = fabs( -0.00161 ) + ( exp( -5.955 * p - 3.052 ) + 0.01751 ) * nSigma; break;
580 }
581
582 return cut;
583}
584
585///////////////////////////////
586double EventAssemblyAlg::phiCut( double p, double nSigma, int parId ) {
587 if ( p < 0.3 ) p = 0.3;
588 if ( parId < 1 || parId > 3 ) parId = 3;
589 if ( nSigma <= 0 ) nSigma = 5.;
590 double cut = myEmcPhiCut;
591 switch ( parId )
592 {
593 case 1:
594 cut = fabs( exp( -2.187 * p - 2.945 ) + 0.03236 ) +
595 ( exp( -2.977 * p - 3.558 ) + 0.005566 ) * nSigma;
596 break;
597 case 2:
598 cut = fabs( exp( -2.216 * p - 2.13 ) + 0.03314 ) +
599 ( exp( -2.436 * p - 3.392 ) + 0.005049 ) * nSigma;
600 break;
601 case 3:
602 cut = fabs( exp( -1.63 * p - 2.931 ) + 0.03223 ) +
603 ( exp( -3.192 * p - 2.756 ) + 0.01533 ) * nSigma;
604 break;
605 }
606
607 return cut;
608}
DECLARE_COMPONENT
DECLARE_COMPONENT(BesBdkRc)
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition CosmicGenerator.h:36
EvtRecTrackCol
ObjectVector< EvtRecTrack > EvtRecTrackCol
Definition Event/EvtRecEvent/include/EvtRecEvent/EvtRecTrack.h:104
exp
EvtComplex exp(const EvtComplex &c)
Definition EvtComplex.hh:209
cut
*********Class see also m_nmax DOUBLE PRECISION m_MasPhot DOUBLE PRECISION m_phsu DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_r2 DOUBLE PRECISION m_WtMass INTEGER m_nmax INTEGER m_Nevgen INTEGER m_IsFSR INTEGER m_MarTot *COMMON c_KarFin $ !Output file $ !Event serial number $ !alpha QED at Thomson limit $ !minimum energy at CMS for remooval $ !infrared cut
Definition KarFin.h:27
ntupleSvc
INTupleSvc * ntupleSvc()
Definition ServiceAccessor.h:51
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:12
M_PI
#define M_PI
Definition TConstant.h:4
DstMdcTrack::helix
const HepVector helix() const
......
Definition DstMdcTrack.cxx:147
EventAssemblyAlg::EventAssemblyAlg
EventAssemblyAlg(const std::string &name, ISvcLocator *pSvcLocator)
Definition EventAssemblyAlg.cxx:22
EventAssemblyAlg::initialize
StatusCode initialize()
Definition EventAssemblyAlg.cxx:36
EvtRecTrack::setMdcKalTrack
void setMdcKalTrack(const RecMdcKalTrack *trk)
Definition Event/EvtRecEvent/include/EvtRecEvent/EvtRecTrack.h:76
EvtRecTrack::setEmcShower
void setEmcShower(const RecEmcShower *shower)
Definition Event/EvtRecEvent/include/EvtRecEvent/EvtRecTrack.h:81
EvtRecTrack::addTofTrack
void addTofTrack(const SmartRef< RecTofTrack > trk)
Definition Event/EvtRecEvent/include/EvtRecEvent/EvtRecTrack.h:80
TofID::is_mrpc
static bool is_mrpc(const Identifier &id)
Definition TofID.cxx:98
TofID::phi_module
static int phi_module(const Identifier &id)
Definition TofID.cxx:65
TofID::barrel_ec
static int barrel_ec(const Identifier &id)
Values of different levels (failure returns 0).
Definition TofID.cxx:54
TofID::layer
static int layer(const Identifier &id)
Definition TofID.cxx:59
TofID::module
static int module(const Identifier &id)
Definition TofID.cxx:114