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