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