BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
MucRecTrkExt.cxx
Go to the documentation of this file.
1#include "EmcRawEvent/EmcDigi.h"
2#include "EmcRecEventModel/RecEmcEventModel.h" //2006.11.11
3#include "EventModel/EventHeader.h"
4#include "ExtEvent/RecExtTrack.h"
5#include "GaudiKernel/IDataManagerSvc.h"
6#include "GaudiKernel/IDataProviderSvc.h"
7#include "GaudiKernel/INTupleSvc.h"
8#include "GaudiKernel/IPartPropSvc.h"
9#include "GaudiKernel/ISvcLocator.h"
10#include "GaudiKernel/MsgStream.h"
11#include "GaudiKernel/PropertyMgr.h"
12#include "GaudiKernel/SmartDataPtr.h"
13#include "Identifier/Identifier.h"
14#include "Identifier/MucID.h"
15#include "McTruth/EmcMcHit.h"
16#include "McTruth/McKine.h"
17#include "McTruth/McParticle.h"
18#include "McTruth/McRelTableDefs.h"
19#include "McTruth/MdcMcHit.h"
20#include "McTruth/MucMcHit.h"
21#include "McTruth/TofMcHit.h"
22#include "MdcRawEvent/MdcDigi.h"
23#include "MdcRecEvent/RecMdcTrack.h"
24#include "MucRawEvent/MucDigi.h"
25#include "ReconEvent/ReconEvent.h"
26#include "TofRawEvent/TofDigi.h"
27
28#include "MucGeomSvc/IMucGeomSvc.h"
29#include "MucGeomSvc/MucConstant.h"
30#include "MucGeomSvc/MucGeoGap.h"
31#include "MucGeomSvc/MucGeoGeneral.h"
32#include "MucGeomSvc/MucGeoStrip.h"
33#include "MucRecEvent/MucRecHit.h"
34#include "MucRecEvent/MucRecHitContainer.h"
35#include "MucRecEvent/RecMucTrack.h"
36#include "MucRecTrkExt.h"
37#include "MucRecTrkExtParameter.h"
38#include <fstream>
39#include <iomanip>
40#include <iostream>
41#include <vector>
42
43using namespace std;
44using namespace Event;
45
46/////////////////////////////////////////////////////////////////////////////
47DECLARE_COMPONENT( MucRecTrkExt )
48MucRecTrkExt::MucRecTrkExt( const std::string& name, ISvcLocator* pSvcLocator )
49 : Algorithm( name, pSvcLocator ) {
50 // Declare the properties
51 declareProperty( "ExtTrackSeedMode",
52 m_ExtTrackSeedMode = 1 ); // 1: My ext from Mc track, 2: from Ext track, 3:
53 // from MUC for cosmic ray
54 declareProperty( "CompareWithMcHit", m_CompareWithMcHit = 0 );
55 declareProperty( "FittingMethod", m_fittingMethod = 1 );
56 declareProperty( "EmcShowerSeedMode",
57 m_EmcShowerSeedMode = 0 ); // 0: Not use EmcShower as seed, 1: use...
58 declareProperty( "MucHitSeedMode",
59 m_MucHitSeedMode = 0 ); // 0: Not use MucHit as seed,1 use...
60 declareProperty( "ConfigFile", m_configFile = "MucConfig.xml" );
61 declareProperty( "Blind", m_Blind = false );
62 declareProperty( "NtOutput", m_NtOutput = 0 ); // NTuple save to root or not
63 declareProperty( "MsOutput", m_MsOutput = false ); // Debug message cout or not
64 declareProperty( "FilterFile", m_filter_filename );
65}
66
67// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
68StatusCode MucRecTrkExt::initialize() {
69 MsgStream log( msgSvc(), name() );
70 log << MSG::INFO << "in initialize()" << endmsg;
71
72 // Get the Particle Properties Service
73 IPartPropSvc* p_PartPropSvc;
74 static const bool CREATEIFNOTTHERE( true );
75 StatusCode PartPropStatus = service( "PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE );
76 if ( !PartPropStatus.isSuccess() || 0 == p_PartPropSvc )
77 {
78 log << MSG::WARNING << " Could not initialize Particle Properties Service" << endmsg;
79 return PartPropStatus;
80 }
81 m_particleTable = p_PartPropSvc->PDT();
82
83 m_totalEvent = 0;
84
85 m_NDigisTotal = 0;
86 m_NHitsTotal = 0;
87 m_NHitsFoundTotal = 0;
88 m_NHitsLostTotal = 0;
89 m_NHitsMisFoundTotal = 0;
90 m_NHitsLostByMdcTotal = 0;
91 m_NHitsLostByExtTotal = 0;
92
93 m_NTracksTotal = 0;
94 m_NTracksRecoTotal = 0;
95 m_NTracksLostByMdcTotal = 0;
96 m_NTracksLostByExtTotal = 0;
97 m_NTracksMdcGhostTotal = 0;
98
99 for ( int i = 0; i < 20; i++ ) m_NHitsLost.push_back( 0 );
100 for ( int i = 0; i < 10; i++ ) m_NHitsLostInGap.push_back( 0 );
101
102 IMucGeomSvc* mucGeomSvc;
103 StatusCode sc = service( "MucGeomSvc", mucGeomSvc );
104 if ( sc == StatusCode::SUCCESS )
105 {
106 // cout <<"dump"<<endl;
107 mucGeomSvc->Dump();
108 // cout<<"Hi, event routine is running"<<endl;
109 }
110 else { return StatusCode::FAILURE; }
111 m_MucRecHitContainer = 0;
112
113 //--------------book ntuple------------------
114 // NTuplePtr nt1(ntupleSvc(), "MyTuples/1");
115 if ( m_NtOutput >= 1 )
116 {
117 NTuplePtr nt1( ntupleSvc(), "FILE401/T" );
118
119 if ( nt1 ) { m_tuple = nt1; }
120 else
121 {
122 // m_tuple = ntupleSvc()->book ("MyTuples/1", CLID_RowWiseTuple, "MdcTrkRecon
123 // N-Tuple");
124 m_tuple = ntupleSvc()->book( "FILE401/T", CLID_RowWiseTuple, "MucTrkRecon N-Tuple" );
125 if ( m_tuple )
126 {
127 sc = m_tuple->addItem( "posx", m_posx );
128 sc = m_tuple->addItem( "posy", m_posy );
129 sc = m_tuple->addItem( "posz", m_posz );
130 sc = m_tuple->addItem( "posx_ext", m_posx_ext );
131 sc = m_tuple->addItem( "posy_ext", m_posy_ext );
132 sc = m_tuple->addItem( "posz_ext", m_posz_ext );
133 sc = m_tuple->addItem( "posxsigma", m_posx_sigma );
134 sc = m_tuple->addItem( "posysigma", m_posy_sigma );
135 sc = m_tuple->addItem( "poszsigma", m_posz_sigma );
136 sc = m_tuple->addItem( "Depth", m_depth );
137 sc = m_tuple->addItem( "Distance", m_Distance );
138 sc = m_tuple->addItem( "MaxHits", m_MaxHits );
139 sc = m_tuple->addItem( "Chi2", m_Chi2 );
140 sc = m_tuple->addItem( "Dist_x", m_Dist_x );
141 sc = m_tuple->addItem( "Dist_y", m_Dist_y );
142 sc = m_tuple->addItem( "Dist_z", m_Dist_z );
143 sc = m_tuple->addItem( "px_mc", m_px_mc );
144 sc = m_tuple->addItem( "py_mc", m_py_mc );
145 sc = m_tuple->addItem( "pz_mc", m_pz_mc );
146 sc = m_tuple->addItem( "emctrack", m_emctrack );
147 sc = m_tuple->addItem( "muchasdigi", m_muc_digi );
148
149 sc = m_tuple->addItem( "part", m_part );
150 sc = m_tuple->addItem( "seg", m_seg );
151 sc = m_tuple->addItem( "gap", m_gap );
152 sc = m_tuple->addItem( "strip", m_strip );
153 sc = m_tuple->addItem( "diff", m_diff );
154 sc = m_tuple->addItem( "distance", m_distance );
155 sc = m_tuple->addItem( "run", m_run );
156 sc = m_tuple->addItem( "event", m_event );
157 }
158 else
159 { // did not manage to book the N tuple....
160 log << MSG::ERROR << " Cannot book N-tuple:" << long( m_tuple ) << endmsg;
161 // return StatusCode::FAILURE;
162 }
163 }
164 }
165
166 // load the filter event
167 if ( m_filter_filename == "" )
168 {
169 m_filter_filename = getenv( "MUCRECALGROOT" );
170 m_filter_filename += "/share/filter.txt";
171 }
172
173 if ( m_filter_filename.size() )
174 {
175 std::ifstream infile( m_filter_filename.c_str() );
176
177 while ( !infile.eof() )
178 {
179 FilterEvent filterevt;
180 infile >> filterevt.bossver >> filterevt.runid >> filterevt.eventid;
181 if ( ( !infile.good() ) || ( infile.eof() ) ) { break; }
182
183 m_filter_event.push_back( filterevt );
184 // cout << filterevt.bossver << " "
185 // << filterevt.runid << " "
186 // << filterevt.eventid << std::endl;
187 }
188 }
189
190 //--------------end of book ntuple------------------
191
192 return StatusCode::SUCCESS;
193}
194
195// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
196StatusCode MucRecTrkExt::execute() {
197
198 MsgStream log( msgSvc(), name() );
199 log << MSG::INFO << "in execute()" << endmsg;
200
201 StatusCode sc = StatusCode::SUCCESS;
202
203 // Part 1: Get the event header, print out event and run number
204 int event, run;
205
206 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
207 if ( !eventHeader )
208 {
209 log << MSG::FATAL << "Could not find Event Header" << endmsg;
210 // return( StatusCode::FAILURE);
211 }
212 m_totalEvent++;
213 log << MSG::INFO << "Event: " << m_totalEvent
214 << "\tcurrent run: " << eventHeader->runNumber()
215 << "\tcurrent event: " << eventHeader->eventNumber() << endmsg;
216
217 event = eventHeader->eventNumber();
218 run = eventHeader->runNumber();
219
220 string release = getenv( "BES_RELEASE" );
221 // if(release=="6.6.5"&&run==35896&&event==7448){return ( StatusCode::SUCCESS);}
222 // if(release=="6.6.5"&&run==37241&&event==1690731){return ( StatusCode::SUCCESS);}
223 // filter the event
224 for ( std::vector<FilterEvent>::iterator it = m_filter_event.begin();
225 it != m_filter_event.end(); ++it )
226 {
227 const FilterEvent& fe = ( *it );
228 if ( release == fe.bossver && run == fe.runid && event == fe.eventid )
229 {
230 cout << "SKIP: " << fe.bossver << " " << fe.runid << " " << fe.eventid << std::endl;
231 return StatusCode::SUCCESS;
232 }
233 }
234
235 // Part 2: Retrieve MC truth
236 if ( m_CompareWithMcHit == 1 )
237 {
238 SmartDataPtr<McParticleCol> mcParticleCol( eventSvc(), "/Event/MC/McParticleCol" );
239 if ( !mcParticleCol )
240 {
241 log << MSG::FATAL << "Could not find McParticle" << endmsg;
242 // return( StatusCode::FAILURE);
243 }
244
245 McParticleCol::iterator iter_mc = mcParticleCol->begin();
246 // do not loop; just get first particle info
247
248 // if(!(*iter_mc)->primaryParticle()) continue;
249 int pid = ( *iter_mc )->particleProperty();
250 int charge = 0;
251 double mass = -1.0;
252
253 if ( pid > 0 )
254 {
255 if ( m_particleTable->particle( pid ) )
256 {
257 charge = (int)m_particleTable->particle( pid )->charge();
258 mass = m_particleTable->particle( pid )->mass();
259 }
260 else
261 log << MSG::ERROR << "no this particle id in particle table, please check data"
262 << endmsg;
263 }
264 else if ( pid < 0 )
265 {
266 if ( m_particleTable->particle( -pid ) )
267 {
268 charge = (int)m_particleTable->particle( -pid )->charge();
269 charge *= -1;
270 mass = m_particleTable->particle( -pid )->mass();
271 }
272 else
273 log << MSG::ERROR << "no this particle id in particle table, please check data"
274 << endmsg;
275 }
276 else { log << MSG::WARNING << "wrong particle id, please check data" << endmsg; }
277
278 // if(!charge) {
279 // log << MSG::WARNING
280 // << "neutral particle charge = 0!!! ...just skip it ! " << endmsg;
281 // continue;
282 // }
283
284 HepLorentzVector initialMomentum = ( *iter_mc )->initialFourMomentum();
285 HepLorentzVector initialPos = ( *iter_mc )->initialPosition();
286 if ( m_NtOutput >= 1 )
287 {
288 m_px_mc = initialMomentum.px();
289 m_py_mc = initialMomentum.py();
290 m_pz_mc = initialMomentum.pz();
291 }
292 // cout<<"mc mom: "<<m_px_mc<<endl;
293
294 log << MSG::INFO << " particleId = " << ( *iter_mc )->particleProperty() << endmsg;
295 }
296
297 // Part 6: Retrieve MUC digi
298 int digiId = 0;
299 SmartDataPtr<MucDigiCol> mucDigiCol( eventSvc(), "/Event/Digi/MucDigiCol" );
300 if ( !mucDigiCol )
301 {
302 log << MSG::FATAL << "Could not find MUC digi" << endmsg;
303 return ( StatusCode::FAILURE );
304 }
305 MucDigiCol::iterator iter4 = mucDigiCol->begin();
306 for ( ; iter4 != mucDigiCol->end(); iter4++, digiId++ ) {}
307 log << MSG::INFO << "Total MUC digis:\t" << digiId << endmsg;
308 if ( digiId == 0 ) return ( StatusCode::SUCCESS );
309
310 // Retrieve MdcMcHit
311 /*
312 SmartDataPtr<Event::MdcMcHitCol> aMdcMcHitCol(eventSvc(),"/Event/MC/MdcMcHitCol");
313 if (!aMdcMcHitCol) {
314 log << MSG::WARNING << "Could not find MdcMcHitCol" << endmsg;
315 //return( StatusCode::FAILURE);
316 }
317 //log << MSG::DEBUG << "MdcMcHitCol contains " << aMdcMcHitCol->size() << " Hits " << endmsg;
318
319 // Retrieve TofMcHit
320 SmartDataPtr<Event::TofMcHitCol> aTofMcHitCol(eventSvc(),"/Event/MC/TofMcHitCol");
321 if (!aTofMcHitCol) {
322 log << MSG::WARNING << "Could not find TofMcHitCol" << endmsg;
323 //return( StatusCode::FAILURE);
324 }
325 //log << MSG::DEBUG << "TofMcHitCol contains " << aTofMcHitCol->size() << " Hits " << endmsg;
326
327 // Retrieve EmcMcHit
328 SmartDataPtr<Event::EmcMcHitCol> aEmcMcHitCol(eventSvc(),"/Event/MC/EmcMcHitCol");
329 if (!aEmcMcHitCol) {
330 log << MSG::WARNING << "Could not find EmcMcHitCol" << endmsg;
331 //return( StatusCode::FAILURE);
332 }
333 //log << MSG::DEBUG << "EmcMcHitCol contains " << aEmcMcHitCol->size() << " Hits " << endmsg;
334 */
335
336 Identifier mucID;
337
338 // McPartToMucHitTab assocMcMuc;
339 // assocMcMuc.init();
340
341 if ( m_CompareWithMcHit )
342 {
343 McPartToMucHitTab assocMcMuc;
344 assocMcMuc.init();
345
346 SmartDataPtr<McParticleCol> mcParticleCol( eventSvc(), "/Event/MC/McParticleCol" );
347 if ( !mcParticleCol )
348 {
349 log << MSG::FATAL << "Could not find McParticle" << endmsg;
350 // return( StatusCode::FAILURE);
351 }
352 McParticleCol::iterator iter_mc = mcParticleCol->begin();
353
354 // Retrieve MucMcHit
355 SmartDataPtr<Event::MucMcHitCol> aMucMcHitCol( eventSvc(), "/Event/MC/MucMcHitCol" );
356 if ( !aMucMcHitCol )
357 {
358 log << MSG::WARNING << "Could not find MucMcHitCol" << endmsg;
359 // return( StatusCode::FAILURE);
360 }
361
362 log << MSG::DEBUG << "MucMcHitCol contains " << aMucMcHitCol->size() << " Hits " << endmsg;
363 MucMcHitCol::iterator iter_MucMcHit = aMucMcHitCol->begin();
364 for ( ; iter_MucMcHit != aMucMcHitCol->end(); iter_MucMcHit++ )
365 {
366 mucID = ( *iter_MucMcHit )->identify();
367
368 log << MSG::DEBUG << " MucMcHit "
369 << " : "
370 << " part " << MucID::barrel_ec( mucID ) << " seg " << MucID::segment( mucID )
371 << " gap " << MucID::layer( mucID ) << " strip " << MucID::channel( mucID )
372 << " Track Id " << ( *iter_MucMcHit )->getTrackIndex() << " pos x "
373 << ( *iter_MucMcHit )->getPositionX() << " pos y "
374 << ( *iter_MucMcHit )->getPositionY() << " pos z "
375 << ( *iter_MucMcHit )->getPositionZ() << endmsg;
376
377 McParticle* assocMcPart = 0;
378 for ( iter_mc = mcParticleCol->begin(); iter_mc != mcParticleCol->end(); iter_mc++ )
379 {
380 if ( ( *iter_mc )->trackIndex() == ( *iter_MucMcHit )->getTrackIndex() )
381 {
382 assocMcPart = *iter_mc;
383 break;
384 }
385 }
386 if ( assocMcPart == 0 )
387 {
388 log << MSG::WARNING << "No Corresponding Mc Particle found for this MucMcHit"
389 << endmsg;
390 }
391
392 MucMcHit* assocMucMcHit = *iter_MucMcHit;
393 McPartToMucHitRel* relMcMuc = 0;
394 relMcMuc = new McPartToMucHitRel( assocMcPart, assocMucMcHit );
395 if ( relMcMuc == 0 ) log << MSG::DEBUG << "relMcMuc not created " << endmsg;
396 else
397 {
398 bool addstat = assocMcMuc.addRelation( relMcMuc );
399 if ( !addstat ) delete relMcMuc;
400 }
401 }
402
403 log << MSG::DEBUG << " Fill McPartToMucHitTab, size " << assocMcMuc.size() << endmsg;
404 iter_mc = mcParticleCol->begin();
405 for ( ; iter_mc != mcParticleCol->end(); iter_mc++ )
406 {
407 log << MSG::DEBUG << " Track index " << ( *iter_mc )->trackIndex()
408 << " particleId = " << ( *iter_mc )->particleProperty() << endmsg;
409
410 vector<McPartToMucHitRel*> vecMucMcHit = assocMcMuc.getRelByFirst( *iter_mc );
411 vector<McPartToMucHitRel*>::iterator iter_MucMcHit = vecMucMcHit.begin();
412 for ( ; iter_MucMcHit != vecMucMcHit.end(); iter_MucMcHit++ )
413 {
414 mucID = ( *iter_MucMcHit )->getSecond()->identify();
415
416 log << MSG::DEBUG
417 //<< " McPartToMucHitTab " << iter_assocMcMuc
418 << " MC Particle index " << ( *iter_MucMcHit )->getFirst()->trackIndex()
419 << " contains "
420 << " MucMcHit "
421 << " part " << MucID::barrel_ec( mucID ) << " seg " << MucID::segment( mucID )
422 << " gap " << MucID::layer( mucID ) << " strip "
423 << MucID::channel( mucID )
424 //<< " posX " << (*iter_MucMcHit)->getSecond()->getPositionX()
425 //<< " posY " << (*iter_MucMcHit)->getSecond()->getPositionY()
426 //<< " posZ " << (*iter_MucMcHit)->getSecond()->getPositionZ()
427 << endmsg;
428 }
429 }
430
431 // assocMcPart = new McParticle(**iter_mc);
432 // MucMcHit *assocMucMcHit = new MucMcHit(mucID, (*iter_MucMcHit)->getTrackIndex(),
433 // (*iter_MucMcHit)->getPositionX(),
434 // (*iter_MucMcHit)->getPositionY(),
435 //(*iter_MucMcHit)->getPositionZ(),
436 // (*iter_MucMcHit)->getPx(), (*iter_MucMcHit)->getPy(),
437 //(*iter_MucMcHit)->getPz() );
438
439 // Retrieve McPartToMucHitTab
440 // SmartDataPtr<McPartToMucHitTab>
441 // aMcPartToMucHitTab(eventSvc(),"/Event/MC/McPartToMucHitTab"); if (!aMcPartToMucHitTab) {
442 // log << MSG::ERROR << "Could not find McPartToMucHitTab" << endmsg;
443 // return( StatusCode::FAILURE);
444 //}
445 }
446
447 //
448 // Read in Muc Digi;
449 if ( !m_MucRecHitContainer ) { m_MucRecHitContainer = new MucRecHitContainer(); }
450 m_MucRecHitContainer->Clear();
451 MucRecHitCol* aMucRecHitCol = new MucRecHitCol();
452 m_MucRecHitContainer->SetMucRecHitCol( aMucRecHitCol );
453
454 SmartIF<IDataManagerSvc> dataManSvc( eventSvc() );
455 DataObject* mucRecHitCol;
456 eventSvc()->findObject( "/Event/Recon/MucRecHitCol", mucRecHitCol );
457 if ( mucRecHitCol != NULL )
458 {
459 dataManSvc->clearSubTree( "/Event/Recon/MucRecHitCol" );
460 eventSvc()->unregisterObject( "/Event/Recon/MucRecHitCol" );
461 }
462
463 sc = eventSvc()->registerObject( "/Event/Recon/MucRecHitCol", aMucRecHitCol );
464 if ( !sc )
465 {
466 log << MSG::ERROR << "/Event/Recon/MucRecHitCol not registerd!" << endmsg;
467 return ( StatusCode::FAILURE );
468 }
469
470 log << MSG::INFO << "Add digis" << endmsg;
471 MucDigiCol::iterator iter_Muc = mucDigiCol->begin();
472 int mucDigiId = 0;
473 for ( ; iter_Muc != mucDigiCol->end(); iter_Muc++, mucDigiId++ )
474 {
475 mucID = ( *iter_Muc )->identify();
476 int part = MucID::barrel_ec( mucID );
477 int seg = MucID::segment( mucID );
478 int gap = MucID::layer( mucID );
479 int strip = MucID::channel( mucID );
480 // m_MucRecHitContainer->AddHit(mucID);
481 m_MucRecHitContainer->AddHit( part, seg, gap, strip );
482
483 log << MSG::DEBUG << " digit" << mucDigiId << " : "
484 << " part " << part << " seg " << seg << " gap " << gap << " strip " << strip
485 << endmsg;
486 }
487
488 //
489 // Create track Container;
490 RecMucTrackCol* aRecMucTrackCol = new RecMucTrackCol();
491
492 // Register RecMucTrackCol to TDS
493 DataObject* aReconEvent;
494 eventSvc()->findObject( "/Event/Recon", aReconEvent );
495 if ( aReconEvent == NULL )
496 {
497 // then register Recon
498 aReconEvent = new ReconEvent();
499 StatusCode sc = eventSvc()->registerObject( "/Event/Recon", aReconEvent );
500 if ( sc != StatusCode::SUCCESS )
501 {
502 log << MSG::FATAL << "Could not register ReconEvent" << endmsg;
503 return ( StatusCode::FAILURE );
504 }
505 }
506 StatusCode fr = eventSvc()->findObject( "/Event/Recon", aReconEvent );
507 if ( fr != StatusCode::SUCCESS )
508 {
509 log << MSG::WARNING << "Could not find register ReconEvent, will create it" << endmsg;
510 StatusCode sc = eventSvc()->registerObject( "/Event/Recon", aReconEvent );
511 if ( sc != StatusCode::SUCCESS )
512 {
513 log << MSG::FATAL << "Could not register ReconEvent" << endmsg;
514 return ( StatusCode::FAILURE );
515 }
516 }
517
518 DataObject* mucTrackCol;
519 eventSvc()->findObject( "/Event/Recon/RecMucTrackCol", mucTrackCol );
520 if ( mucTrackCol != NULL )
521 {
522 dataManSvc->clearSubTree( "/Event/Recon/RecMucTrackCol" );
523 eventSvc()->unregisterObject( "/Event/Recon/RecMucTrackCol" );
524 }
525
526 sc = eventSvc()->registerObject( "/Event/Recon/RecMucTrackCol", aRecMucTrackCol );
527 if ( sc != StatusCode::SUCCESS )
528 {
529 log << MSG::FATAL << "Could not register MUC track collection" << endmsg;
530 return ( StatusCode::FAILURE );
531 }
532
533 // check RecMucTrackCol registered
534 SmartDataPtr<RecMucTrackCol> findRecMucTrackCol( eventSvc(), "/Event/Recon/RecMucTrackCol" );
535 if ( !findRecMucTrackCol )
536 {
537 log << MSG::FATAL << "Could not find RecMucTrackCol" << endmsg;
538 return ( StatusCode::FAILURE );
539 }
540 aRecMucTrackCol->clear();
541
542 // m_ExtTrackSeedMode 1: ext from MC track, 2: use Ext track
543
544 // Retrieve Ext track Col from TDS
545 // Check ExtTrackCol in TDS.
546 SmartDataPtr<RecExtTrackCol> aExtTrackCol( eventSvc(), "/Event/Recon/RecExtTrackCol" );
547 if ( !aExtTrackCol )
548 {
549 log << MSG::WARNING << "Can't find ExtTrackCol in TDS!" << endmsg;
550 // return( StatusCode::FAILURE);
551 }
552
553 SmartDataPtr<RecMdcTrackCol> aMdcTrackCol( eventSvc(), "/Event/Recon/RecMdcTrackCol" );
554 if ( !aMdcTrackCol )
555 {
556 log << MSG::WARNING << "Can't find MdcTrackCol in TDS!" << endmsg;
557 // return( StatusCode::FAILURE);
558 }
559
560 // Retrieve Emc track col from TDS 2006.11.11 liangyt
561 // Check EmcRecShowerCol in TDS.
562 SmartDataPtr<RecEmcShowerCol> aShowerCol( eventSvc(), "/Event/Recon/RecEmcShowerCol" );
563 if ( !aShowerCol )
564 {
565 log << MSG::WARNING << "Could not find RecEmcShowerCol" << endmsg;
566 // return( StatusCode::FAILURE);
567 }
568
569 // EmcRecShowerCol::iterator iShowerCol;
570 // for(iShowerCol=aShowerCol->begin();
571 // iShowerCol!=aShowerCol->end();
572 // iShowerCol++){
573 // cout<<"check EmcRecShowerCol registered:"<<endl;
574 // <<"shower position: "<<(*iShowerCol)->Position()<<endl;
575 // }
576
577 int muctrackId = 0;
578
579 log << MSG::INFO << "Add tracks info, ExtTrackSeedMode = " << m_ExtTrackSeedMode
580 << ", EmcShowerSeedMode = " << m_EmcShowerSeedMode << endmsg;
581 // if (m_ExtTrackSeedMode == 1 || !aExtTrackCol) {
582 if ( m_ExtTrackSeedMode == 1 )
583 {
584 SmartDataPtr<McParticleCol> mcParticleCol( eventSvc(), "/Event/MC/McParticleCol" );
585 if ( !mcParticleCol )
586 {
587 log << MSG::FATAL << "Could not find McParticle" << endmsg;
588 // return( StatusCode::FAILURE);
589 return ( StatusCode::SUCCESS );
590 }
591 McParticleCol::iterator iter_mc = mcParticleCol->begin();
592
593 int trackIndex = -99;
594 for ( int iTrack = 0; iter_mc != mcParticleCol->end(); iter_mc++, iTrack++ )
595 {
596 if ( !( *iter_mc )->primaryParticle() ) continue;
597 int pid = ( *iter_mc )->particleProperty();
598 int charge = 0;
599 double mass = -1.0;
600 if ( pid > 0 )
601 {
602 if ( m_particleTable->particle( pid ) )
603 {
604 charge = (int)m_particleTable->particle( pid )->charge();
605 mass = m_particleTable->particle( pid )->mass();
606 }
607 else
608 log << MSG::ERROR << "no this particle id in particle table, please check data"
609 << endmsg;
610 }
611 else if ( pid < 0 )
612 {
613 if ( m_particleTable->particle( -pid ) )
614 {
615 charge = (int)m_particleTable->particle( -pid )->charge();
616 charge *= -1;
617 mass = m_particleTable->particle( -pid )->mass();
618 }
619 else
620 log << MSG::ERROR << "no this particle id in particle table, please check data"
621 << endmsg;
622 }
623 else { log << MSG::WARNING << "wrong particle id, please check data" << endmsg; }
624
625 if ( !charge )
626 {
627 log << MSG::WARNING << " neutral particle charge = 0!!! ...just skip it !" << endmsg;
628 continue;
629 }
630
631 trackIndex = ( *iter_mc )->trackIndex();
632 log << MSG::DEBUG << "iTrack " << iTrack << " index " << trackIndex
633 << " particleId = " << ( *iter_mc )->particleProperty() << endmsg;
634
635 RecMucTrack* aTrack = new RecMucTrack();
636 aTrack->setTrackId( trackIndex );
637 aTrack->setId( muctrackId );
638
639 HepLorentzVector initialMomentum = ( *iter_mc )->initialFourMomentum();
640 HepLorentzVector initialPos = ( *iter_mc )->initialPosition();
641 float theta0 = initialMomentum.theta();
642 float phi0 = initialMomentum.phi();
643 // float dirX = sin(theta0) * cos(phi0);
644 // float dirY = sin(theta0) * sin(phi0);
645 // float dirZ = cos(theta0);
646 float x0 = initialPos.x();
647 float y0 = initialPos.y();
648 float z0 = initialPos.z();
649
650 Hep3Vector startPos( x0, y0, z0 );
651 Hep3Vector startP( initialMomentum.px(), initialMomentum.py(), initialMomentum.pz() );
652 log << MSG::DEBUG << "startP " << startP << " startPos " << startPos << endmsg;
653 Hep3Vector endPos( 0, 0, 0 ), endP;
654
655 aTrack->GetMdcExtTrack( startPos, startP, charge, endPos, endP );
656 aTrack->SetMdcPos( x0, y0, z0 );
657 aTrack->SetMdcMomentum( startP.x(), startP.y(), startP.z() );
658 aTrack->SetExtTrackID( trackIndex );
659 aTrack->SetExtMucPos( endPos.x(), endPos.y(), endPos.z() );
660 aTrack->SetExtMucMomentum( endP.x(), endP.y(), endP.z() );
661 aTrack->SetMucPos( endPos.x(), endPos.y(), endPos.z() );
662 aTrack->SetMucMomentum( endP.x(), endP.y(), endP.z() );
663 aTrack->SetCurrentPos( endPos.x(), endPos.y(), endPos.z() );
664 aTrack->SetCurrentDir( endP.x(), endP.y(), endP.z() );
665 // aTrack->SetMucVertexPos( x0, y0, z0);
666 // aTrack->SetMucVertexDir( dirX, dirY, dirZ );
667 // aTrack->SetCurrentPos( x0, y0, z0);
668 // aTrack->SetCurrentDir( dirX, dirY, dirZ );
669
670 // log << MSG::DEBUG
671 //<< " pdg " << (*iter_mc)->particleProperty()
672 //<< " mucPos " << aTrack->GetMucVertexPos()
673 //<< " mucDir " << aTrack->GetMucVertexDir()
674 //<< endmsg;
675 log << MSG::INFO << "Add track, trkIndex: " << trackIndex << "\tmucTrkId: " << muctrackId
676 << endmsg;
677 aRecMucTrackCol->add( aTrack );
678 muctrackId++;
679 }
680 }
681 else if ( m_ExtTrackSeedMode == 2 )
682 {
683 if ( !aExtTrackCol || !aMdcTrackCol )
684 {
685 log << MSG::WARNING << "Can't find ExtTrackCol or MdcTrackCol in TDS!" << endmsg;
686 return StatusCode::SUCCESS;
687 }
688
689 int trackIndex = -99;
690 double kdep = -99;
691 double krechi = 0.;
692 int kdof = 0;
693 int kbrLay = -1;
694 int kecLay = -1;
695 // log << MSG::INFO << "MdcTrackCol:\t " << aMdcTrackCol->size() << "\tExtTrackCol:\t " <<
696 // aExtTrackCol->size() << endmsg;
697 RecExtTrackCol::iterator iter_ExtTrack = aExtTrackCol->begin();
698 RecMdcTrackCol::iterator iter_MdcTrack = aMdcTrackCol->begin();
699 int iExtTrack = 0;
700 for ( ; iter_ExtTrack != aExtTrackCol->end() && iter_MdcTrack != aMdcTrackCol->end();
701 iter_ExtTrack++, iter_MdcTrack++, iExtTrack++ )
702 {
703 trackIndex = ( *iter_ExtTrack )->GetTrackId();
704 log << MSG::DEBUG << "iExtTrack " << iExtTrack << " index " << trackIndex << " MucPos "
705 << iExtTrack << ( *iter_ExtTrack )->mucPosition().x() << " "
706 << ( *iter_ExtTrack )->mucPosition().y() << " "
707 << ( *iter_ExtTrack )->mucPosition().z() << " r "
708 << ( *iter_ExtTrack )->mucPosition().r() << endmsg;
709
710 if ( ( *iter_ExtTrack )->mucPosition().x() == -99 &&
711 ( *iter_ExtTrack )->mucPosition().y() == -99 &&
712 ( *iter_ExtTrack )->mucPosition().z() == -99 )
713 {
714 log << MSG::INFO << "Bad ExtTrack, trackIndex: " << trackIndex << ", skip" << endmsg;
715 continue;
716 }
717
718 // added by LI Chunhua 2013/02/01
719 krechi = ( *iter_ExtTrack )->MucKalchi2();
720 kdof = ( *iter_ExtTrack )->MucKaldof();
721 kdep = ( *iter_ExtTrack )->MucKaldepth();
722 kbrLay = ( *iter_ExtTrack )->MucKalbrLastLayer();
723 kecLay = ( *iter_ExtTrack )->MucKalecLastLayer();
724 if ( kdof <= 0 ) krechi = 0.;
725 else krechi = krechi / kdof;
726 kdep = kdep / 10.;
727 //*********************************
728 RecMucTrack* aTrack = new RecMucTrack();
729 aTrack->setTrackId( trackIndex );
730 aTrack->setId( muctrackId );
731
732 aTrack->setType( 0 ); // 0 for use seed from mdc ext;
733
734 aTrack->SetExtTrack( *iter_ExtTrack );
735
736 // added by LI Chunhua 2013/02/01
737 aTrack->setkalRechi2( krechi );
738 aTrack->setkalDof( kdof );
739 aTrack->setkalDepth( kdep );
740 aTrack->setkalbrLastLayer( kbrLay );
741 aTrack->setkalecLastLayer( kecLay );
742 //******************************
743 // Hep3Vector mdcPos = (*iter_ExtTrack)->GetMdcPosition();
744 Hep3Vector mdcMomentum = ( *iter_MdcTrack )->p3();
745 Hep3Vector mucPos = ( *iter_ExtTrack )->mucPosition();
746 Hep3Vector mucMomentum = ( *iter_ExtTrack )->mucMomentum();
747
748 // aTrack->SetMdcPos( mdcPos.x(), mdcPos.y(), mdcPos.z() );
749 aTrack->SetMdcMomentum( mdcMomentum.x(), mdcMomentum.y(), mdcMomentum.z() );
750 aTrack->SetExtMucPos( mucPos.x(), mucPos.y(), mucPos.z() );
751 // cout << "ExtTrack MucPos " << aTrack->GetExtMucPos() << endl;
752 aTrack->SetExtMucMomentum( mucMomentum.x(), mucMomentum.y(), mucMomentum.z() );
753 aTrack->SetMucPos( mucPos.x(), mucPos.y(), mucPos.z() );
754 aTrack->SetMucMomentum( mucMomentum.x(), mucMomentum.y(), mucMomentum.z() );
755 aTrack->SetCurrentPos( mucPos.x(), mucPos.y(), mucPos.z() );
756 aTrack->SetCurrentDir( mucMomentum.x(), mucMomentum.y(), mucMomentum.z() );
757 aTrack->SetRecMode( 0 ); // mdc ext mode
758
759 log << MSG::INFO << "Add track, trkIndex: " << trackIndex << "\tmucTrkId: " << muctrackId
760 << endmsg;
761 aRecMucTrackCol->add( aTrack );
762 muctrackId++;
763 }
764 }
765 else if ( m_ExtTrackSeedMode == 3 )
766 { // cosmic ray
767
768 if ( !aMdcTrackCol )
769 {
770 log << MSG::WARNING << "Can't find MdcTrackCol in TDS!" << endmsg;
771 return StatusCode::SUCCESS;
772 // return( StatusCode::FAILURE);
773 }
774 log << MSG::INFO << "Mdc track size: " << aMdcTrackCol->size() << endmsg;
775
776 int trackIndex = -99;
777 for ( RecMdcTrackCol::iterator iter_mdc1 = aMdcTrackCol->begin();
778 iter_mdc1 != aMdcTrackCol->end(); iter_mdc1++ )
779 {
780 // if((*iter_mdc1)->getFi0() > 0.5*pi && (*iter_mdc1)->getFi0() < 1.5*pi){
781 // cout<<"this is down track"<<endl;
782 // }
783
784 trackIndex = ( *iter_mdc1 )->trackId();
785 HepVector helix = ( *iter_mdc1 )->helix();
786 // cout<<"helix: "<<helix<<endl;
787
788 float x0 = helix[0] * cos( helix[1] );
789 float y0 = helix[0] * sin( helix[1] );
790 float z0 = helix[3];
791
792 // float dx = 1/(sqrt(1+helix[4]*helix[4])) * (-1* sin(helix[1]));
793 // float dy = 1/(sqrt(1+helix[4]*helix[4])) * cos(helix[1]);
794 // float dz = 1/(sqrt(1+helix[4]*helix[4])) * helix[4];
795
796 float dx = -1 * sin( helix[1] );
797 float dy = cos( helix[1] );
798 float dz = helix[4];
799
800 // cout<<"pos: "<<x0<<" "<<y0<<" "<<z0<<" dir: "<<dx<<" "<<dy<<" "<<dz<<endl;
801
802 RecMucTrack* aTrack = new RecMucTrack();
803 aTrack->setTrackId( trackIndex );
804 aTrack->setId( muctrackId );
805 muctrackId++;
806
807 aTrack->setType( 3 ); // 3 for use seed from mdc without magnet field;
808
809 Hep3Vector mucPos( x0, y0, z0 );
810 Hep3Vector mucMomentum( dx, dy, dz );
811
812 aTrack->SetExtMucPos( x0, y0, z0 );
813 aTrack->SetExtMucMomentum( dx, dy, dz );
814
815 aTrack->SetMucPos( x0, y0, z0 );
816 aTrack->SetMucMomentum( dx, dy, dz );
817 aTrack->SetCurrentPos( x0, y0, z0 );
818 aTrack->SetCurrentDir( dx, dy, dz );
819 aTrack->SetRecMode( 0 ); // mdc ext mode
820
821 log << MSG::INFO << "Add track, trkIndex: " << trackIndex << "\tmucTrkId: " << muctrackId
822 << endmsg;
823 aRecMucTrackCol->add( aTrack );
824 muctrackId++;
825 }
826 }
827 else { log << MSG::INFO << "ExtTrackSeedMode error!" << endmsg; }
828
829 // Rec from Emc extend track
830 if ( m_EmcShowerSeedMode == 1 )
831 {
832 int trackIndex = 999;
833 RecEmcShowerCol::iterator iShowerCol;
834 for ( iShowerCol = aShowerCol->begin(); iShowerCol != aShowerCol->end(); iShowerCol++ )
835 {
836 // cout<<"shower id: "<<(*iShowerCol)->ShowerId()<<" "
837 // <<"shower energy: "<<(*iShowerCol)->Energy()<<" "
838 // <<"shower position: "<<(*iShowerCol)->Position()<<endl;
839
840 RecMucTrack* aTrack = new RecMucTrack();
841 aTrack->setTrackId( trackIndex );
842 aTrack->setId( muctrackId );
843 aTrack->setType( 1 ); // 1 for use seed from emc hits;
844
845 Hep3Vector mucPos = ( *iShowerCol )->position();
846 Hep3Vector mucMomentum = ( *iShowerCol )->position();
847
848 aTrack->SetExtMucPos( mucPos.x(), mucPos.y(), mucPos.z() );
849 // cout << "EmcTrack MucPos " << aTrack->GetExtMucPos() << endl;
850 aTrack->SetExtMucMomentum( mucMomentum.x(), mucMomentum.y(), mucMomentum.z() );
851 aTrack->SetMucPos( mucPos.x(), mucPos.y(), mucPos.z() );
852 aTrack->SetMucMomentum( mucMomentum.x(), mucMomentum.y(), mucMomentum.z() );
853 aTrack->SetCurrentPos( mucPos.x(), mucPos.y(), mucPos.z() );
854 aTrack->SetCurrentDir( mucMomentum.x(), mucMomentum.y(), mucMomentum.z() );
855 aTrack->SetRecMode( 1 ); // emc ext mode
856 // cout<<" MucRecTrkExt 111 trackIndex :"<<aTrack->GetIndex()<<endl;
857
858 log << MSG::INFO << "Add track, trkIndex: " << trackIndex << "\tmucTrkId: " << muctrackId
859 << endmsg;
860 aRecMucTrackCol->add( aTrack );
861 muctrackId++;
862
863 m_emcrec = 1;
864 }
865 }
866 log << MSG::DEBUG << " track container filled " << endmsg;
867
868 // All input filled, begin track finding;
869 log << MSG::INFO << "Start tracking..." << endmsg;
870 int runVerbose = 1;
871 RecMucTrack* aTrack = 0;
872 for ( int iTrack = 0; iTrack < (int)aRecMucTrackCol->size(); iTrack++ )
873 {
874 log << MSG::DEBUG << "iTrack " << iTrack << endmsg;
875 aTrack = ( *aRecMucTrackCol )[iTrack];
876 // cout<<"in MucRecTrkExt trackIndex :"<<aTrack->GetIndex()<<endl;
877
878 Hep3Vector currentPos = aTrack->GetCurrentPos();
879 Hep3Vector currentDir = aTrack->GetCurrentDir();
880 if ( currentPos.mag() < kMinor )
881 {
882 log << MSG::WARNING << "No MUC intersection in track " << iTrack << endmsg;
883 continue;
884 }
885 // log << MSG::INFO << " pos " << currentPos << " dir " << currentDir << endmsg;
886
887 int firstHitFound[2] = { 0, 0 }; // Has the fist position in this orient determined? if so,
888 // could CorrectDirection()
889 int firstHitGap[2] = { -1, -1 }; // When the fist position in this orient determined, the
890 // gap it is on
891 for ( int partSeq = 0; partSeq < (int)MucID::getPartNum(); partSeq++ )
892 {
893 int iPart = kPartSeq[partSeq];
894 for ( int iGap = 0; iGap < (int)MucID::getGapNum( iPart ); iGap++ )
895 {
896 // log << MSG::INFO << iPart << iGap << endmsg;
897 int seg = -1;
898 int orient = MucGeoGeneral::Instance()->GetGap( iPart, 0, iGap )->Orient();
899 ;
900
901 float xInsct, yInsct, zInsct;
902 aTrack->Project( iPart, iGap, xInsct, yInsct, zInsct, seg );
903 if ( m_MsOutput )
904 cout << "part " << iPart << " gap " << iGap << " x " << xInsct << " y " << yInsct
905 << " z " << zInsct << " seg " << seg << endl;
906
907 if ( seg == -1 ) continue;
908
909 aTrack->SetCurrentInsct( xInsct, yInsct, zInsct );
910
911 for ( int iDeltaSeg = 0; iDeltaSeg < kNSegSearch; iDeltaSeg++ )
912 {
913 int iSeg = seg + kDeltaSeg[iDeltaSeg]; // also find on neighbor seg;
914 if ( iSeg < 0 ) iSeg += MucID::getSegNum( iPart );
915 if ( iSeg > (int)MucID::getSegNum( iPart ) - 1 ) iSeg -= MucID::getSegNum( iPart );
916
917 // log << MSG::DEBUG << iPart << iSeg << iGap
918 // << "NHits " << m_MucRecHitContainer->GetGapHitCount(iPart, seg, iGap) << endmsg;
919
920 //----------now change window size(depond on mom)
921 // float window = kWindowSize[iPart][iGap];
922 Hep3Vector mom_mdc = aTrack->getMdcMomentum();
923 float window = getWindowSize( mom_mdc, iPart, iSeg, iGap );
924
925 if ( firstHitFound[orient] != 1 )
926 window *= kFirstHitWindowRatio; // if no hits have been found on this orient,
927 // expand the window.
928 for ( int iHit = 0; iHit < m_MucRecHitContainer->GetGapHitCount( iPart, iSeg, iGap );
929 iHit++ )
930 {
931 log << MSG::DEBUG << "iSeg " << iSeg << " iHit " << iHit << endmsg;
932 MucRecHit* pHit = m_MucRecHitContainer->GetHit( iPart, iSeg, iGap, iHit );
933 // cout<< "strip " << pHit->Strip() << " center " << pHit->GetCenterPos() << endl;
934
935 if ( !pHit )
936 {
937 log << MSG::WARNING << "MucRecTrkExt: null pointer to pHit" << endmsg;
938 continue;
939 }
940 else
941 {
942 // cout<<"in MucRecTrkExt: pHit exist : "<<iPart<<" "<<iSeg<<" "<<iGap<<"
943 // "<<iHit<<" mode:"<<pHit->GetHitMode()<<endl;
944 // Get the displacement of hit pHit to intersection
945 // float dX = aTrack->GetHitDistance(pHit);
946 float dX = aTrack->GetHitDistance2( pHit ); // not abs value now!
947 log << MSG::DEBUG << "distance = " << setw( 8 ) << dX << " size " << setw( 4 )
948 << window << endmsg;
949
950 if ( m_NtOutput >= 2 )
951 { // too many info
952 m_distance = dX;
953 m_part = iPart;
954 m_seg = iSeg;
955 m_gap = iGap;
956 m_strip = pHit->Strip();
957 m_diff = -99;
958 m_tuple->write();
959 }
960
961 if ( fabs( dX ) < window )
962 {
963 // Attach the hit if it exists
964 // cout << "meet window " << pHit->GetSoftID() << endl;
965 //****************if this if emc track, we abondon used hit in
966 // mdc*******************
967 // if(m_emcrec == 1 )
968 if ( aTrack->GetRecMode() == 0 )
969 {
970 pHit->SetHitMode( 1 ); // mdc ext
971 aTrack->AttachHit( pHit );
972 // cout<<"in MucRecTrkExt: trackmode==0 so mdc ext "<<iPart<<" "<<iSeg<<"
973 // "<<iGap<<" "<<iHit<<endl;
974 }
975 if ( aTrack->GetRecMode() == 1 )
976 {
977 // cout<<"in MucRecTrkExt: HitMode = "<<pHit->GetHitMode()<<" part:
978 // "<<iPart<<" "<<iSeg<<" "<<iGap<<"
979 // "<<iHit<<endl;
980 if ( pHit->GetHitMode() != 1 )
981 {
982 aTrack->AttachHit( pHit ); // this hit has not been used by mdc ext
983 pHit->SetHitMode( 2 ); // emc ext
984 }
985 }
986
987 // push back distance between ext and hits
988 aTrack->pushExtDistHits( dX ); // 2009-03-12
989
990 if ( firstHitGap[orient] == -1 ) firstHitGap[orient] = iGap;
991 firstHitFound[orient] = 1;
992 // cout << "You could correct directon in orient " << orient << endl;
993
994 log << MSG::DEBUG << " part " << iPart << " seg " << iSeg << " gap " << iGap
995 << " strip " << setw( 2 ) << pHit->Strip() << " attatched" << endmsg;
996 log << MSG::DEBUG << "current total hits " << aTrack->GetTotalHits() << endmsg;
997 }
998 else
999 {
1000 m_NHitsLostInGap[iGap]++;
1001 // log << MSG::DEBUG cout << " part " << iPart << " seg " << iSeg << " gap " <<
1002 // iGap
1003 // << " strip " << setw(2) << pHit->GetSoftID().GetStrip()
1004 // << " not attached !" << endmsg;
1005 }
1006 } // end pHit else
1007 } // end for iHit
1008 aTrack->CalculateInsct( iPart, iSeg, iGap );
1009 } // end for iDeltaSeg
1010
1011 // When correct dir in the orient is permitted and this gap is not the gap first hit
1012 // locates.
1013 if ( m_ExtTrackSeedMode != 3 && !m_Blind )
1014 { // for cosmic ray, we do not correct dir and pos...
1015 if ( firstHitFound[orient] && firstHitGap[orient] != iGap ) aTrack->CorrectDir();
1016 aTrack->CorrectPos();
1017 }
1018 // cout << "Current Intersection " << aTrack->GetCurrentInsct() << endl;
1019 // cout << "Current Direction " << aTrack->GetCurrentDir() << endl;
1020 aTrack->AttachInsct( aTrack->GetCurrentInsct() );
1021 } // end for iGap
1022 } // end for iSeg
1023 aTrack->LineFit( m_fittingMethod );
1024 aTrack->ComputeTrackInfo( m_fittingMethod );
1025 log << MSG::INFO << "Fit track done! trackIndex: " << aTrack->trackId()
1026 << ", mucId: " << aTrack->id() << ", RecMode: " << aTrack->GetRecMode() << endmsg;
1027
1028 if ( m_NtOutput >= 1 )
1029 {
1030 m_depth = aTrack->depth();
1031 m_Distance = aTrack->distance();
1032 m_MaxHits = aTrack->maxHitsInLayer();
1033 m_Chi2 = aTrack->chi2();
1034 m_Dist_x = aTrack->GetExtMucDistance().x();
1035 m_Dist_y = aTrack->GetExtMucDistance().y();
1036 m_Dist_z = aTrack->GetExtMucDistance().z();
1037 m_posx_ext = aTrack->GetMucStripPos().x();
1038 m_posy_ext = aTrack->GetMucStripPos().y();
1039 m_posz_ext = aTrack->GetMucStripPos().z();
1040
1041 m_emctrack = m_emcrec;
1042 }
1043
1044 // test distance of line/quad fitting
1045 if ( m_NtOutput >= 2 )
1046 {
1047 vector<MucRecHit*> attachedHits = aTrack->GetHits();
1048 vector<MucRecHit*> expectedHits = aTrack->GetExpectedHits();
1049 vector<float> distanceHits = aTrack->getDistHits();
1050 vector<float> distanceHits_quad = aTrack->getQuadDistHits();
1051 vector<float> distanceHits_ext = aTrack->getExtDistHits();
1052
1053 for ( int i = 0; i < expectedHits.size(); i++ )
1054 {
1055 MucRecHit* ihit = expectedHits[i];
1056 // cout<<"expected Hits: "<<ihit->Part()<<" "<<ihit->Seg()<<" "<<ihit->Gap()<<"
1057 // "<<ihit->Strip()<<endl;
1058 }
1059
1060 for ( int j = 0; j < attachedHits.size(); j++ )
1061 {
1062 MucRecHit* jhit = attachedHits[j];
1063 // if(attachedHits.size()!=distanceHits.size())cout<<"attached hits size
1064 // no same as disthits!!!"<<endl;
1065 if ( attachedHits.size() == distanceHits.size() )
1066 { // same gap, cale distance
1067 m_part = jhit->Part();
1068 m_seg = jhit->Seg();
1069 m_gap = jhit->Gap();
1070 m_strip = jhit->Strip();
1071 m_distance = distanceHits[j];
1072 m_Distance = distanceHits_quad[j];
1073 m_diff = -9999;
1074 // cout<<"distance = "<<m_dist<<endl;
1075 // if(distanceHits.size() == distanceHits_ext.size())
1076 // cout<<"ext dist: "<<distanceHits_ext[j]<<endl;
1077 m_tuple->write();
1078 }
1079 }
1080 } // end if output
1081
1082 int nHitsAttached = aTrack->GetTotalHits();
1083 // m_NHitsAttachedTotal += nHitsAttached;
1084 // if(mucDigiCol->size() - recHitNum != 0)
1085 log << MSG::DEBUG << "track Index " << aTrack->trackId() << setw( 2 )
1086 << mucDigiCol->size() - nHitsAttached << " of " << setw( 2 ) << mucDigiCol->size()
1087 << " lost " << endmsg;
1088 // m_NHitsLost[int(mucDigiCol->size() - nHitsAttached)]++;
1089 } // end for iTrack
1090
1091 // if (aRecMucTrackCol->size() == 0) m_NHitsLost[0]++;
1092 // m_NHitsTotal += mucDigiCol->size();
1093
1094 //****************** look up unrec hits to do recon with MUC info ***************??
1095
1096 if ( m_MucHitSeedMode == 1 )
1097 {
1098 MucRecHit *pHit = 0, *pHit0 = 0, *pHit1 = 0;
1099 int count, orient;
1100
1101 for ( int iPart = 0; iPart < (int)MucID::getPartNum(); iPart++ )
1102 {
1103 for ( int iSeg = 0; iSeg < (int)MucID::getSegNum( iPart ); iSeg++ )
1104 {
1105 bool hit0 = false, hit1 = false;
1106 int firstgap0 = -1, firstgap1 = -1;
1107 int nStrip0 = 0, nStrip1 = 0;
1108 Hep3Vector posHit0, posHit1;
1109 for ( int iGap = 0; iGap < (int)MucID::getGapNum( iPart ); iGap++ )
1110 {
1111 count = m_MucRecHitContainer->GetGapHitCount( iPart, iSeg, iGap );
1112 for ( int iHit0 = 0; iHit0 < count; iHit0++ )
1113 {
1114 // cout << "iHit0 = " << iHit0 << endl;
1115 pHit = m_MucRecHitContainer->GetHit( iPart, iSeg, iGap, iHit0 );
1116 if ( !pHit )
1117 {
1118 log << MSG::FATAL << "MucRecRoadFinder-E10 "
1119 << " part " << iPart << " seg " << iSeg << " gap " << iGap << " null pointer"
1120 << endl;
1121 }
1122 else
1123 {
1124 // cout<< "pHit Mode is : " << pHit->GetHitMode() << endl;
1125 if ( pHit->GetHitMode() == -1 )
1126 {
1127 orient = MucGeoGeneral::Instance()->GetGap( iPart, 0, iGap )->Orient();
1128 if ( orient == 1 && hit0 == false )
1129 {
1130 hit0 = true;
1131 firstgap0 = iGap;
1132 }
1133 if ( iGap == firstgap0 )
1134 {
1135 posHit0 += pHit->GetCenterPos();
1136 nStrip0++;
1137 }
1138
1139 if ( orient == 0 && hit1 == false )
1140 {
1141 hit1 = true;
1142 firstgap1 = iGap;
1143 }
1144 if ( iGap == firstgap1 )
1145 {
1146 posHit1 += pHit->GetCenterPos();
1147 nStrip1++;
1148 }
1149 // if(orient == 1 && hit0 == false){ //orient = 1
1150 // hit0 = true;
1151 // pHit0 = pHit; //pHit0 is the first hit of first gap in this
1152 // segment with orient 1
1153 // }
1154 // if(orient == 0 && hit1 == false){
1155 // hit1 = true;
1156 // pHit1 = pHit; //pHit0 is the first hit of first gap in this
1157 // segment with orient 0
1158 // }
1159 }
1160 } // pHit0 exist;
1161 } // iHit0
1162 } // iGap
1163
1164 // if hit0 hit1 exist, make a seed
1165 if ( hit0 && hit1 )
1166 {
1167 posHit0 /= nStrip0;
1168 posHit1 /= nStrip1;
1169 // cout<< "pHit0 position is " << posHit0 <<" "<<nStrip0<<endl;
1170 // cout<< "pHit1 position is " << posHit1 <<" "<<nStrip1<<endl;
1171
1172 int trackIndex = 999;
1173 RecMucTrack* aTrack = new RecMucTrack();
1174 aTrack->setTrackId( trackIndex );
1175 aTrack->setId( muctrackId );
1176 muctrackId++;
1177
1178 aTrack->setType( 2 ); // 2 for use seed from Muc Information
1179
1180 Hep3Vector mucPos, mucMomentum;
1181 if ( iPart == 1 )
1182 {
1183 mucMomentum.set( posHit0.x(), posHit0.y(), posHit1.z() );
1184 mucPos = mucMomentum * 0.9; // shorten mucPos, otherwise maybe no intersection
1185 // point found for the track and this gap
1186 }
1187 else
1188 {
1189 mucMomentum.set( posHit0.x(), posHit1.y(), posHit0.z() * 0.5 + posHit1.z() * 0.5 );
1190 mucPos = mucMomentum * 0.9; // shorten mucPos, otherwise maybe no intersection
1191 // point found for the track and this gap
1192 }
1193
1194 aTrack->SetExtMucPos( mucPos.x(), mucPos.y(), mucPos.z() );
1195 // cout << "EmcTrack MucPos " << aTrack->GetExtMucPos() << endl;
1196 aTrack->SetExtMucMomentum( mucMomentum.x(), mucMomentum.y(), mucMomentum.z() );
1197 aTrack->SetMucPos( mucPos.x(), mucPos.y(), mucPos.z() );
1198 aTrack->SetMucMomentum( mucMomentum.x(), mucMomentum.y(), mucMomentum.z() );
1199 aTrack->SetCurrentPos( mucPos.x(), mucPos.y(), mucPos.z() );
1200 aTrack->SetCurrentDir( mucMomentum.x(), mucMomentum.y(), mucMomentum.z() );
1201 // cout<<" MucRecTrkExt 111 trackIndex :"<<aTrack->GetIndex()<<endl;
1202 aTrack->SetRecMode( 2 );
1203 aRecMucTrackCol->add( aTrack );
1204
1205 TrackFinding( aTrack );
1206 }
1207
1208 if ( m_NtOutput >= 1 )
1209 {
1210 m_depth = aTrack->depth();
1211 m_Distance = aTrack->distance();
1212 m_MaxHits = aTrack->maxHitsInLayer();
1213 m_Chi2 = aTrack->chi2();
1214 m_Dist_x = aTrack->GetExtMucDistance().x();
1215 m_Dist_y = aTrack->GetExtMucDistance().y();
1216 m_Dist_z = aTrack->GetExtMucDistance().z();
1217 m_posx_ext = aTrack->GetMucStripPos().x();
1218 m_posy_ext = aTrack->GetMucStripPos().y();
1219 m_posz_ext = aTrack->GetMucStripPos().z();
1220
1221 m_emctrack = 2;
1222 }
1223
1224 } // iSeg
1225 } // iPart
1226 } // end if SeedMode == 1
1227
1228 //************** end of look up unrec hits to do recon with MUC info ************??
1229 int nTracksTotal = 0; // mcParticleCol->size();
1230 int nTracksFound = 0;
1231 int nTracksLost = 0;
1232 int nTracksLostByExt = 0;
1233 int nTracksMisFound = 0;
1234
1235 int nDigisTotal = 0; // mucDigiCol->size();
1236 int nHitsTotal = 0; // assocMcMuc.size();
1237 int nHitsFound = 0;
1238 int nHitsLost = 0;
1239 int nHitsMisFound = 0;
1240
1241 /*
1242 if (m_CompareWithMcHit) {
1243 //
1244 // Compare RecMucTracks with McPartToMucHitTab;
1245
1246 log << MSG::DEBUG << " *******************************" << endmsg;
1247 log << MSG::DEBUG << " Compare Mc Info with rec tracks" << endmsg;
1248 log << MSG::DEBUG << " McParticleCol size " << mcParticleCol->size() << endmsg;
1249 log << MSG::DEBUG << " McPartToMcuHitTab size " << assocMcMuc.size() << endmsg;
1250
1251 iter_mc = mcParticleCol->begin();
1252 for (;iter_mc != mcParticleCol->end(); iter_mc++) {
1253 log << MSG::DEBUG << " McParticle index " << (*iter_mc)->getTrackIndex()
1254 << " particleId = " << (*iter_mc)->particleProperty()
1255 << endmsg;
1256
1257 bool foundRecoTrack = false;
1258 log << MSG::DEBUG << " Reconed RecMucTrackCol size " << aRecMucTrackCol->size() <<
1259 endmsg; aTrack = 0; for(int iTrack = 0; iTrack < (int)aRecMucTrackCol->size(); iTrack++) {
1260 log << MSG::DEBUG << "iTrack " << iTrack << endmsg;
1261 aTrack = (*aRecMucTrackCol)[iTrack];
1262 if (aTrack->GetExtTrackID() == (*iter_mc)->getTrackIndex()) {
1263 log << MSG::DEBUG << " found iTrack " << iTrack
1264 << " index " << aTrack->GetExtTrackID()
1265 << " equals to " << " McParticle index " << (*iter_mc)->getTrackIndex()
1266 << endmsg;
1267 foundRecoTrack = true;
1268 break;
1269 }
1270 }
1271 if (foundRecoTrack == false) {
1272 nTracksLost++;
1273 m_TrackLostByMdc.push_back(eventHeader->eventNumber());
1274 log << MSG::DEBUG << " this McParticle is not found in RecMucTrackCol,"
1275 << "->not intialized from ExtTrack-> not constructed in Mdc" << endmsg;
1276 }
1277 else {
1278 nTracksFound++;
1279 }
1280
1281 int nHitsFoundInTrack = 0;
1282 int nHitsLostInTrack = 0;
1283
1284 vector< McPartToMucHitRel* > vecMucMcHit = assocMcMuc.getRelByFirst(*iter_mc);
1285 log << MSG::DEBUG << " McParticle " << (*iter_mc)->getTrackIndex()
1286 << " contains " << vecMucMcHit.size() << " MucMcHit " << endmsg;
1287 if (!aTrack) {
1288 nHitsLostInTrack = vecMucMcHit.size();
1289 m_NHitsLostByMdcTotal += nHitsLostInTrack;
1290 log << MSG::DEBUG << " This Mc particle has no corresponding Reco Track, "
1291 << " all of its MucMcHits lost" << endmsg;
1292 }
1293 else if ((aTrack->GetExtMucPos()).mag() < 0.1) {
1294 nHitsLostInTrack = vecMucMcHit.size();
1295 m_NHitsLostByExtTotal += nHitsLostInTrack;
1296 nTracksLostByExt++;
1297 m_TrackLostByExt.push_back(eventHeader->eventNumber());
1298 log << MSG::DEBUG << " This Mc particle 's Ext track's intersection with Muc is
1299 zero,"
1300 << " all of its MucMcHits lost" << endmsg;
1301 }
1302 else {
1303 vector< McPartToMucHitRel* >::iterator iter_MucMcHit = vecMucMcHit.begin();
1304 for (; iter_MucMcHit != vecMucMcHit.end(); iter_MucMcHit++) {
1305 mucID = (*iter_MucMcHit)->getSecond()->identify();
1306 int part = MucID::barrel_ec(mucID);
1307 int seg = MucID::segment(mucID);
1308 int gap = MucID::layer(mucID);
1309 int strip = MucID::channel(mucID);
1310
1311 log << MSG::DEBUG
1312 //<< " McPartToMucHitTab " << iter_assocMcMuc
1313 << " McParticle index "
1314 << (*iter_MucMcHit)->getFirst()->getTrackIndex()
1315 << " contains " << " MucMcHit "
1316 << " part " << part
1317 << " seg " << seg
1318 << " gap " << gap
1319 << " strip " << strip
1320 //<< " posX " << (*iter_MucMcHit)->getSecond()->getPositionX()
1321 // << " posY " << (*iter_MucMcHit)->getSecond()->getPositionY()
1322 //<< " posZ " << (*iter_MucMcHit)->getSecond()->getPositionZ()
1323 << endmsg;
1324
1325 if (aTrack->HasHit(part, seg, gap, strip)) {
1326 nHitsFoundInTrack++;
1327 log << MSG::DEBUG << " This MucMchit found on this Reco track" << endmsg;
1328 }
1329 else {
1330 nHitsLostInTrack++;
1331 log << MSG::DEBUG << " This MucMchit not found on this Reco track, it is lost "
1332 << endmsg;
1333 }
1334 }
1335 }
1336
1337 nHitsFound += nHitsFoundInTrack;
1338 nHitsLost += nHitsLostInTrack;
1339
1340 m_NHitsLost[nHitsLost]++;
1341 if (nHitsLost > 0) {
1342 m_TrackLostHit.push_back(eventHeader->eventNumber());
1343 m_TrackLostHitNumber.push_back(nHitsLost);
1344 }
1345
1346 log << MSG::DEBUG << " In " << vecMucMcHit.size() << " MucMcHit : "
1347 << nHitsFoundInTrack << " found in Reco track, "
1348 << nHitsLostInTrack << " lost " << endmsg;
1349 }
1350
1351 // Reversely, Compare McPartToMucHitTab with RecMucTracks;
1352 log << MSG::DEBUG << " *******************************" << endmsg;
1353 log << MSG::DEBUG << " Compare rec tracks with Mc Info" << endmsg;
1354 log << MSG::DEBUG << " Reconed RecMucTrackCol size " << aRecMucTrackCol->size() <<
1355 endmsg; for(int iTrack = 0; iTrack < (int)aRecMucTrackCol->size(); iTrack++) { log <<
1356 MSG::DEBUG << "iTrack " << iTrack << endmsg; aTrack = (*aRecMucTrackCol)[iTrack];
1357 //cout << "MucPos " << aTrack->GetMucPos() << " MucMomentum " <<
1358 aTrack->GetMucMomentum() << endl;
1359
1360 bool foundMcParticle = false;
1361 iter_mc = mcParticleCol->begin();
1362 for (;iter_mc != mcParticleCol->end(); iter_mc++) {
1363 if ((*iter_mc)->getTrackIndex() == aTrack->GetExtTrackID()) {
1364 log << MSG::DEBUG << " found McParticle index " << (*iter_mc)->getTrackIndex()
1365 << " equals to " << " Reconed track " << iTrack
1366 << " index " << aTrack->GetExtTrackID()
1367 << endmsg;
1368 foundMcParticle = true;
1369 break;
1370 }
1371 }
1372 if (foundMcParticle == false) {
1373 nTracksMisFound++;
1374 log << MSG::DEBUG << " This Reconed track has no corresponding McParticle, where is
1375 it from? " << endmsg;
1376 }
1377
1378 int nHitsFoundInMcParticle = 0;
1379 int nHitsMisFoundInMcParticle = 0;
1380 vector< MucRecHit* > vecMucRecHit = aTrack->GetHits();
1381 log << MSG::DEBUG << " Reconed Track " << iTrack
1382 << " index " << aTrack->GetExtTrackID()
1383 << " contains " << aTrack->GetTotalHits() << " MucRecMcHit " << endmsg;
1384 vector< MucRecHit* >::iterator iter_MucRecHit = vecMucRecHit.begin();
1385 for (; iter_MucRecHit != vecMucRecHit.end(); iter_MucRecHit++) {
1386 int part = (*iter_MucRecHit)->Part();
1387 int seg = (*iter_MucRecHit)->Seg();
1388 int gap = (*iter_MucRecHit)->Gap();
1389 int strip = (*iter_MucRecHit)->Strip();
1390
1391 log << MSG::DEBUG
1392 << " Reconed track index "
1393 << aTrack->GetExtTrackID()
1394 << " contains " << " MucRecHit "
1395 << " part " << part
1396 << " seg " << seg
1397 << " gap " << gap
1398 << " strip " << strip;
1399
1400 bool foundMcHit = false;
1401 vector< McPartToMucHitRel* > vecMucMcHit = assocMcMuc.getRelByFirst(*iter_mc);
1402 vector< McPartToMucHitRel* >::iterator iter_MucMcHit = vecMucMcHit.begin();
1403 for (; iter_MucMcHit != vecMucMcHit.end(); iter_MucMcHit++) {
1404 mucID = (*iter_MucMcHit)->getSecond()->identify();
1405 if (part == MucID::barrel_ec(mucID) &&
1406 seg == MucID::segment(mucID) &&
1407 gap == MucID::layer(mucID) &&
1408 strip == MucID::channel(mucID)) {
1409 foundMcHit = true;
1410 break;
1411 }
1412 }
1413
1414 if (foundMcHit == true) {
1415 nHitsFoundInMcParticle++;
1416 log << MSG::DEBUG << " This MucRecHit belongs to this Mc Particle " << endmsg;
1417 }
1418 else {
1419 nHitsMisFoundInMcParticle++;
1420 log << MSG::DEBUG << " This MucRecHit does not belong to this Reco track, it should
1421 not be attached " << endmsg;
1422 }
1423 }
1424
1425 //nHitsFound += nHitsFoundInMcParticle;
1426 nHitsMisFound += nHitsMisFoundInMcParticle;
1427
1428 log << MSG::DEBUG << " In " << vecMucRecHit.size() << " MucRecHit : "
1429 << nHitsFoundInMcParticle << " found in Mc Particle, "
1430 << nHitsMisFoundInMcParticle << " not found in Mc Particle, mis attached "
1431 << endmsg;
1432 }
1433 }
1434 */
1435 /*
1436 log << MSG::INFO << " This event contains " << nTracksTotal << " Mc Particle, "
1437 << nTracksFound << " tracks found, "
1438 << nTracksLost << " tracks lost, "
1439 << nTracksMisFound << " tracks mis found "
1440 << endmsg;
1441 log << MSG::INFO << " This event contains " << nDigisTotal << " Muc Digis " << endmsg;
1442 log << MSG::INFO << " This event contains " << nHitsTotal << " MucMcHits, "
1443 << nHitsFound << " mc hits found, "
1444 << nHitsLost << " mc hits lost, "
1445 << nHitsMisFound << " mc hits mis found "
1446 << endmsg;
1447 */
1448 m_NDigisTotal += nDigisTotal;
1449 m_NHitsTotal += nHitsTotal;
1450 m_NHitsFoundTotal += nHitsFound;
1451 m_NHitsLostTotal += nHitsLost;
1452 m_NHitsMisFoundTotal += nHitsMisFound;
1453 // m_NHitsLost[nHitsLost]++;
1454
1455 m_NTracksTotal += nTracksTotal;
1456 m_NTracksRecoTotal += nTracksFound;
1457 m_NTracksLostByMdcTotal += nTracksLost;
1458 m_NTracksLostByExtTotal += nTracksLostByExt;
1459 m_NTracksMdcGhostTotal += nTracksMisFound;
1460 if ( aRecMucTrackCol->size() > 0 )
1461 {
1462 RecMucTrack* aRecMucTrack = ( *aRecMucTrackCol )[0];
1463 // m_posx = 0.0; m_posy = 0.0; m_posz = 0.0;
1464 if ( m_NtOutput >= 1 )
1465 {
1466 m_posx = aRecMucTrack->getMucPos().x();
1467 m_posy = aRecMucTrack->getMucPos().y();
1468 m_posz = aRecMucTrack->getMucPos().z();
1469
1470 m_posx_sigma = aRecMucTrack->getMucPosSigma().x();
1471 m_posy_sigma = aRecMucTrack->getMucPosSigma().y();
1472 m_posz_sigma = aRecMucTrack->getMucPosSigma().z();
1473 }
1474 // cout << m_posx << " " << m_posy << " " << m_posz << endl;
1475 // m_tuple->write();
1476 }
1477
1478 if ( m_NtOutput >= 1 ) m_tuple->write();
1479 return StatusCode::SUCCESS;
1480}
1481
1482// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
1483StatusCode MucRecTrkExt::finalize() {
1484
1485 MsgStream log( msgSvc(), name() );
1486 log << MSG::INFO << "in finalize()" << endmsg;
1487
1488 log << MSG::INFO << " In " << m_NDigisTotal << " total MucDigis " << endmsg;
1489 log << MSG::INFO << " In " << m_NHitsTotal << " total MucMcHits " << endmsg;
1490 log << MSG::INFO << m_NHitsFoundTotal << " hits found in Reco tracks, " << endmsg;
1491 log << MSG::INFO << m_NHitsLostTotal << " hits lost (not found), of which "
1492 << m_NHitsLostByMdcTotal << " lost because Mdc track not constructed "
1493 << m_NHitsLostByExtTotal << " lost because Ext track not intersect with muc " << endmsg;
1494 log << MSG::INFO << m_NHitsMisFoundTotal
1495 << " hits mis found (not belong to this track, but mis attached)" << endmsg;
1496
1497 log << MSG::INFO << " In " << m_NTracksTotal << " total Mc tracks " << endmsg;
1498 log << MSG::INFO << m_NTracksRecoTotal << " tracks found " << endmsg;
1499 log << MSG::INFO << m_NTracksLostByMdcTotal
1500 << " tracks lost because Mdc track not constructed " << endmsg;
1501 log << MSG::INFO << m_NTracksLostByExtTotal
1502 << " tracks lost because Ext track not intersect with muc " << endmsg;
1503 log << MSG::INFO << m_NTracksMdcGhostTotal << " tracks are Mdc ghost tracks " << endmsg;
1504
1505 /*
1506 for(int i = 0; i < 20; i++) log << MSG::DEBUG << "lost " << i << " hits track " <<
1507 m_NHitsLost[i] << endmsg; for(int i = 0; i < 9; i++) log << MSG::DEBUG << "lost on gap " <<
1508 i << " is " << m_NHitsLostInGap[i] << endmsg; cout << m_TrackLostHit.size() << " track has
1509 hit lost, their event id : " << endl; for (int i = 0; i < m_TrackLostHit.size(); i++) { cout
1510 << m_TrackLostHit[i] << " : " << m_TrackLostHitNumber[i] << endl;
1511 }
1512 cout << m_TrackLostByMdc.size() << " tracks lost by Mdc, their event id : " << endl;
1513 for (int i = 0; i < m_TrackLostByMdc.size(); i++) {
1514 cout << m_TrackLostByMdc[i] << endl;
1515 }
1516 cout << m_TrackLostByExt.size() << " tracks lost by Ext no intersection with muc, their event
1517 id : " << endl; for (int i = 0; i < m_TrackLostByExt.size(); i++) { cout <<
1518 m_TrackLostByExt[i] << endl;
1519 }
1520 */
1521
1522 return StatusCode::SUCCESS;
1523}
1524
1525void MucRecTrkExt::TrackFinding( RecMucTrack* aTrack ) {
1526 MsgStream log( msgSvc(), name() );
1527
1528 Hep3Vector currentPos = aTrack->GetCurrentPos();
1529 Hep3Vector currentDir = aTrack->GetCurrentDir();
1530 // if(currentPos.mag() < kMinor) {
1531 // log << MSG::WARNING << "No MUC intersection in track " << endmsg;
1532 // continue;
1533 // }
1534
1535 int firstHitFound[2] = { 0, 0 }; // Has the fist position in this orient determined? if so,
1536 // could CorrectDirection()
1537 int firstHitGap[2] = { -1, -1 }; // When the fist position in this orient determined, the gap
1538 // it is on
1539 for ( int partSeq = 0; partSeq < (int)MucID::getPartNum(); partSeq++ )
1540 {
1541 int iPart = kPartSeq[partSeq];
1542 for ( int iGap = 0; iGap < (int)MucID::getGapNum( iPart ); iGap++ )
1543 {
1544 int seg = -1;
1545 int orient = MucGeoGeneral::Instance()->GetGap( iPart, 0, iGap )->Orient();
1546 ;
1547 float xInsct, yInsct, zInsct;
1548 aTrack->Project( iPart, iGap, xInsct, yInsct, zInsct, seg );
1549 log << MSG::INFO << "part " << iPart << " gap " << iGap << " x " << xInsct << " y "
1550 << yInsct << " z " << zInsct << " seg " << seg << endmsg;
1551 // cout<<"project: gap="<<iGap<<" seg="<<seg<<" "<<xInsct<<" "<<yInsct<<"
1552 // "<<zInsct<<endl;
1553
1554 if ( seg == -1 )
1555 {
1556 // log << MSG::DEBUG << "no intersection with part " << iPart
1557 // << " gap " << iGap << " in this track !" << endl;
1558 continue;
1559 }
1560 aTrack->SetCurrentInsct( xInsct, yInsct, zInsct );
1561
1562 for ( int iDeltaSeg = 0; iDeltaSeg < kNSegSearch; iDeltaSeg++ )
1563 {
1564 int iSeg = seg + kDeltaSeg[iDeltaSeg]; // also find on neighbor seg;
1565 if ( iSeg < 0 ) iSeg += MucID::getSegNum( iPart );
1566 if ( iSeg > (int)MucID::getSegNum( iPart ) - 1 ) iSeg -= MucID::getSegNum( iPart );
1567
1568 // float window = kWindowSize[iPart][iGap];
1569 Hep3Vector mom_mdc = aTrack->getMdcMomentum();
1570 float window = getWindowSize( mom_mdc, iPart, iSeg, iGap );
1571
1572 if ( firstHitFound[orient] != 1 )
1573 window *= kFirstHitWindowRatio; // if no hits have been found on this orient, expand
1574 // the window.
1575
1576 for ( int iHit = 0; iHit < m_MucRecHitContainer->GetGapHitCount( iPart, iSeg, iGap );
1577 iHit++ )
1578 {
1579 log << MSG::DEBUG << "iSeg " << iSeg << " iHit " << iHit << endmsg;
1580 MucRecHit* pHit = m_MucRecHitContainer->GetHit( iPart, iSeg, iGap, iHit );
1581 // cout<< "strip " << pHit->Strip() << " center " << pHit->GetCenterPos() << endl;
1582
1583 if ( !pHit )
1584 {
1585 log << MSG::WARNING << "MucRecTrkExt: null pointer to pHit" << endmsg;
1586 continue;
1587 }
1588 else
1589 {
1590 // Get the displacement of hit pHit to intersection
1591 float dX = aTrack->GetHitDistance2( pHit );
1592 log << MSG::DEBUG << "distance = " << setw( 8 ) << dX << " size " << setw( 4 )
1593 << window << endmsg;
1594
1595 // cout<<"dX= "<<dX<<" window="<<window<<endl;
1596 if ( fabs( dX ) < window )
1597 {
1598 // Attach the hit if it exists
1599 // cout << "meet window " << pHit->GetSoftID() << endl;
1600 //****************if this if emc track, we abondon used hit in
1601 // mdc*******************
1602 // if(m_emcrec == 1 )
1603 if ( aTrack->GetRecMode() == 0 )
1604 {
1605 pHit->SetHitMode( 1 ); // mdc ext
1606 aTrack->AttachHit( pHit );
1607 // cout<<"in MucRecTrkExt: trackmode==0 so mdc ext "<<iPart<<" "<<iSeg<<"
1608 // "<<iGap<<" "<<iHit<<endl;
1609 }
1610 if ( aTrack->GetRecMode() == 1 )
1611 {
1612 // cout<<"in MucRecTrkExt: HitMode = "<<pHit->GetHitMode()<<" part: "<<iPart<<"
1613 // "<<iSeg<<" "<<iGap<<"
1614 // "<<iHit<<endl;
1615 if ( pHit->GetHitMode() != 1 )
1616 {
1617 aTrack->AttachHit( pHit ); // this hit has not been used by mdc ext
1618 pHit->SetHitMode( 2 ); // emc ext
1619 }
1620 }
1621 if ( aTrack->GetRecMode() == 2 )
1622 {
1623 // cout<<"in MucRecTrkExt: HitMode = "<<pHit->GetHitMode()<<" part: "<<iPart<<"
1624 // "<<iSeg<<" "<<iGap<<"
1625 // "<<iHit<<endl;
1626 if ( pHit->GetHitMode() == -1 )
1627 {
1628 aTrack->AttachHit( pHit ); // this hit has not been used by mdc ext
1629 pHit->SetHitMode( 2 ); // emc ext
1630 }
1631 }
1632
1633 if ( firstHitGap[orient] == -1 ) firstHitGap[orient] = iGap;
1634 firstHitFound[orient] = 1;
1635 // cout << "You could correct directon in orient " << orient << endl;
1636 // cout<< " part " << iPart << " seg " << iSeg << " gap " << iGap
1637 // << " strip " << setw(2) << pHit->Strip() << " attatched" << endl;
1638 log << MSG::DEBUG << " part " << iPart << " seg " << iSeg << " gap " << iGap
1639 << " strip " << setw( 2 ) << pHit->Strip() << " attatched" << endmsg;
1640 log << MSG::DEBUG << "current total hits " << aTrack->GetTotalHits() << endmsg;
1641 }
1642 else
1643 {
1644 m_NHitsLostInGap[iGap]++;
1645 // log << MSG::DEBUG cout << " part " << iPart << " seg " << iSeg << " gap " <<
1646 // iGap
1647 // << " strip " << setw(2) << pHit->GetSoftID().GetStrip()
1648 // << " not attached !" << endmsg;
1649 }
1650 } // end pHit else
1651 } // end iHit
1652
1653 aTrack->CalculateInsct( iPart, iSeg, iGap );
1654 } // end DeltaSeg
1655
1656 // When correct dir in the orient is permitted and this gap is not the gap first hit
1657 // locates.
1658 if ( firstHitFound[orient] && firstHitGap[orient] != iGap ) aTrack->CorrectDir();
1659 aTrack->CorrectPos();
1660 // cout << "Current Intersection " << aTrack->GetCurrentInsct() << endl;
1661 // cout << "Current Direction " << aTrack->GetCurrentDir() << endl;
1662 aTrack->AttachInsct( aTrack->GetCurrentInsct() );
1663 } // end iGap
1664 } // end iPart
1665
1666 aTrack->LineFit( m_fittingMethod );
1667 aTrack->ComputeTrackInfo( m_fittingMethod );
1668
1669 // cout<<" in TrackFinding: depth= "<<aTrack->GetDepth()<<endl;
1670} // End TrackFinding()
1671
1672// the function to get window size with (mom, part, seg, gap); getMdcMomentum
1673float MucRecTrkExt::getWindowSize( Hep3Vector mom, int part, int seg, int gap ) {
1674 int mom_id;
1675 if ( mom.mag() < 0.6 ) mom_id = 0;
1676 else if ( mom.mag() < 0.8 ) mom_id = 1;
1677 else if ( mom.mag() < 1 ) mom_id = 2;
1678 else mom_id = 3;
1679
1680 return kWindowSize_Mom_Gap[mom_id][gap];
1681}
1682
1683// END
DECLARE_COMPONENT
DECLARE_COMPONENT(BesBdkRc)
release
struct sembuf release
Definition BesClient.cxx:112
mass
double mass
Definition CosmicGenerator.cxx:128
count
DOUBLE_PRECISION count[3]
Definition EvtPhokharaDef.hh:63
sin
double sin(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:185
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
kMinor
const double kMinor
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MucGeomSvc/MucConstant.h:8
MucRecHitCol
ObjectVector< MucRecHit > MucRecHitCol
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MucRecEvent/MucRecHit.h:125
RecMucTrackCol
ObjectVector< RecMucTrack > RecMucTrackCol
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MucRecEvent/RecMucTrack.h:357
kNSegSearch
const int kNSegSearch
Definition MucRecRoadFinderParameter.h:68
kPartSeq
const int kPartSeq[3]
Definition MucRecRoadFinderParameter.h:63
kDeltaSeg
const int kDeltaSeg[kNSegSearch]
Definition MucRecRoadFinderParameter.h:69
kWindowSize_Mom_Gap
const float kWindowSize_Mom_Gap[4][9]
Definition MucRecTrkExtParameter.h:13
kFirstHitWindowRatio
const float kFirstHitWindowRatio
Definition MucRecTrkExtParameter.h:21
ntupleSvc
INTupleSvc * ntupleSvc()
Definition ServiceAccessor.h:51
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:12
Event::RelTable::init
void init()
Initialize the internal pointer to an ObjectList of relations.
Definition Event/RelTable/include/RelTable/RelTable.h:33
Event::RelTable::size
unsigned long size() const
This method returns the number of relations in the table.
Definition Event/RelTable/include/RelTable/RelTable.h:226
Event::RelTable::addRelation
bool addRelation(Relation< T1, T2 > *rel)
Definition Event/RelTable/include/RelTable/RelTable.h:127
Event::RelTable::getRelByFirst
std::vector< Relation< T1, T2 > * > getRelByFirst(const T1 *pobj) const
Definition Event/RelTable/include/RelTable/RelTable.h:146
IMucGeomSvc::Dump
virtual void Dump()=0
MucGeoGeneral::GetGap
MucGeoGap * GetGap(const int part, const int seg, const int gap) const
Get a pointer to the gap identified by (part,seg,gap).
Definition MucGeoGeneral.cxx:435
MucGeoGeneral::Instance
static MucGeoGeneral * Instance()
Get a pointer to the single instance of MucGeoGeneral.
Definition MucGeoGeneral.cxx:428
MucID::barrel_ec
static int barrel_ec(const Identifier &id)
Values of different levels.
Definition MucID.cxx:38
MucID::layer
static int layer(const Identifier &id)
Definition MucID.cxx:58
MucID::channel
static int channel(const Identifier &id)
Definition MucID.cxx:68
MucID::segment
static int segment(const Identifier &id)
Definition MucID.cxx:48
MucID::getPartNum
static value_type getPartNum()
Definition MucID.cxx:131
MucID::getSegNum
static value_type getSegNum(int part)
Definition MucID.cxx:134
MucID::getGapNum
static value_type getGapNum(int part)
Definition MucID.cxx:141
MucRecHit::GetCenterPos
Hep3Vector GetCenterPos() const
Get Center position of the strip (global coords).
Definition MucRecHit.cxx:93
MucRecTrkExt::TrackFinding
void TrackFinding(RecMucTrack *aTrack)
Definition MucRecTrkExt.cxx:1525
MucRecTrkExt::execute
StatusCode execute()
Definition MucRecTrkExt.cxx:196
MucRecTrkExt::getWindowSize
float getWindowSize(Hep3Vector mom, int part, int seg, int gap)
Definition MucRecTrkExt.cxx:1673
MucRecTrkExt::MucRecTrkExt
MucRecTrkExt(const std::string &name, ISvcLocator *pSvcLocator)
Definition MucRecTrkExt.cxx:48
MucRecTrkExt::finalize
StatusCode finalize()
Definition MucRecTrkExt.cxx:1483
MucRecTrkExt::initialize
StatusCode initialize()
Definition MucRecTrkExt.cxx:68
RecMucTrack::GetTotalHits
int GetTotalHits() const
How many hits in all does this track contain?
Definition RecMucTrack.cxx:1462
RecMucTrack::LineFit
void LineFit(int fittingMethod)
Line fit with hits on a seg with max hits.
Definition RecMucTrack.cxx:1007
RecMucTrack::CorrectPos
void CorrectPos()
Correct current position of this trajectory.
Definition RecMucTrack.cxx:545
RecMucTrack::CorrectDir
void CorrectDir()
Correct direction of this trajectory.
Definition RecMucTrack.cxx:538
RecMucTrack::ComputeTrackInfo
void ComputeTrackInfo(int fittingmethod)
Get corresponding monte carlo track pointer.
Definition RecMucTrack.cxx:1679
RecMucTrack::SetMdcMomentum
void SetMdcMomentum(const float px, const float py, const float pz)
set start moment of the track in Mdc.
Definition RecMucTrack.cxx:209
RecMucTrack::getMucPos
Hep3Vector getMucPos() const
start position of this track in Muc.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MucRecEvent/RecMucTrack.h:163
RecMucTrack::AttachInsct
void AttachInsct(Hep3Vector insct)
Attach the intersection to this trajectory.
Definition RecMucTrack.cxx:534
RecMucTrack::SetExtMucPos
void SetExtMucPos(const float x, const float y, const float z)
Definition RecMucTrack.cxx:234
RecMucTrack::SetExtTrackID
void SetExtTrackID(int id)
set Ext track id. for compute from mdc myself
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MucRecEvent/RecMucTrack.h:103
RecMucTrack::SetMdcPos
void SetMdcPos(const float x, const float y, const float z)
set start position of the track in Mdc.
Definition RecMucTrack.cxx:205
RecMucTrack::SetCurrentDir
void SetCurrentDir(const float x, const float y, const float z)
set current direction of the trajectory.
Definition RecMucTrack.cxx:246
RecMucTrack::SetExtTrack
void SetExtTrack(RecExtTrack *extTrack)
set Ext track point.
Definition RecMucTrack.cxx:261
RecMucTrack::setTrackId
void setTrackId(const int trackId)
set the index for this track.
Definition RecMucTrack.cxx:199
RecMucTrack::GetExtMucDistance
Hep3Vector GetExtMucDistance() const
Distance match of the ext track with muc track in first layer.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MucRecEvent/RecMucTrack.h:223
RecMucTrack::GetMdcExtTrack
void GetMdcExtTrack(Hep3Vector mdcStartPos, Hep3Vector mdcStartMomentum, int charge, Hep3Vector &mucStartPos, Hep3Vector &mucStartMomentum)
compute ext track myself from mdc.
Definition RecMucTrack.cxx:278
RecMucTrack::getMdcMomentum
Hep3Vector getMdcMomentum() const
momentum of this track in Mdc
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MucRecEvent/RecMucTrack.h:171
RecMucTrack::CalculateInsct
Hep3Vector CalculateInsct(const int part, const int seg, const int gap)
Calculate intersection from all hits attached on this gap.
Definition RecMucTrack.cxx:469
RecMucTrack::Project
void Project(const int &part, const int &gap, float &x, float &y, float &z, int &seg)
Where does the trajectory of this track intersect a specific gap?
Definition RecMucTrack.cxx:354
RecMucTrack::SetMucMomentum
void SetMucMomentum(const float px, const float py, const float pz)
set start moment of the track in Muc.
Definition RecMucTrack.cxx:227
RecMucTrack::GetCurrentInsct
Hep3Vector GetCurrentInsct() const
Current intersection.
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MucRecEvent/RecMucTrack.h:180
RecMucTrack::SetExtMucMomentum
void SetExtMucMomentum(const float px, const float py, const float pz)
set start moment of ext track in Muc.
Definition RecMucTrack.cxx:238
RecMucTrack::GetHits
vector< MucRecHit * > GetHits() const
Get all hits on this track.
Definition RecMucTrack.cxx:1654
RecMucTrack::AttachHit
void AttachHit(MucRecHit *hit)
set corresponding monte carlo track pointer.
Definition RecMucTrack.cxx:325
RecMucTrack::SetCurrentInsct
void SetCurrentInsct(const float x, const float y, const float z)
set current intersection of the trajectory with strip plane.
Definition RecMucTrack.cxx:250
RecMucTrack::SetCurrentPos
void SetCurrentPos(const float x, const float y, const float z)
set current position of the trajectory.
Definition RecMucTrack.cxx:242
RecMucTrack::GetExpectedHits
vector< MucRecHit * > GetExpectedHits() const
Definition RecMucTrack.cxx:1656
RecMucTrack::SetMucPos
void SetMucPos(const float x, const float y, const float z)
set start position of the track in Muc. (after line fit and correction)
Definition RecMucTrack.cxx:213
RecMucTrack::GetHitDistance2
float GetHitDistance2(const MucRecHit *hit)
no abs value
Definition RecMucTrack.cxx:437
Event::McPartToMucHitRel
Event::Relation< Event::McParticle, Event::MucMcHit > McPartToMucHitRel
Definition Event/McTruth/include/McTruth/McRelTableDefs.h:48
Event::McPartToMucHitTab
Event::RelTable< Event::McParticle, Event::MucMcHit > McPartToMucHitTab
Definition Event/McTruth/include/McTruth/McRelTableDefs.h:47