00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include <cmath>
00011 #include <iostream>
00012 #include <vector>
00013
00014 #include "GaudiKernel/MsgStream.h"
00015 #include "GaudiKernel/ThreadGaudi.h"
00016 #include "GaudiKernel/AlgFactory.h"
00017 #include "GaudiKernel/ISvcLocator.h"
00018 #include "GaudiKernel/IMessageSvc.h"
00019 #include "GaudiKernel/IDataProviderSvc.h"
00020 #include "GaudiKernel/SmartDataPtr.h"
00021 #include "GaudiKernel/SmartRef.h"
00022 #include "GaudiKernel/SmartRefVector.h"
00023 #include "GaudiKernel/IDataProviderSvc.h"
00024 #include "GaudiKernel/PropertyMgr.h"
00025 #include "GaudiKernel/IJobOptionsSvc.h"
00026 #include "GaudiKernel/IMessageSvc.h"
00027 #include "GaudiKernel/Bootstrap.h"
00028 #include "GaudiKernel/StatusCode.h"
00029 #include "GaudiKernel/IDataManagerSvc.h"
00030 #include "GaudiKernel/IPartPropSvc.h"
00031
00032 #include "MdcGeomSvc/IMdcGeomSvc.h"
00033 #include "MdcGeomSvc/MdcGeoWire.h"
00034 #include "MdcGeomSvc/MdcGeoLayer.h"
00035 #include "MdcGeomSvc/MdcGeoSuper.h"
00036
00037 #include "Identifier/Identifier.h"
00038 #include "Identifier/MdcID.h"
00039
00040 #include "EventModel/EventHeader.h"
00041 #include "MdcRawEvent/MdcDigi.h"
00042 #include "RawEvent/RawDataUtil.h"
00043 #include "EvTimeEvent/RecEsTime.h"
00044 #include "RawDataProviderSvc/IRawDataProviderSvc.h"
00045 #include "TrigEvent/TrigData.h"
00046 #include "TrigEvent/TrigGTD.h"
00047 #include "TrigEvent/TrigMdc.h"
00048 #include "TrigEvent/TrigGTDProvider.h"
00049
00050
00051 #include "MdcRecEvent/RecMdcTrack.h"
00052 #include "MdcRecEvent/RecMdcHit.h"
00053 #include "ReconEvent/ReconEvent.h"
00054
00055 #include "MdcFastTrkAlg/FTFinder.h"
00056 #include "MdcFastTrkAlg/FTWire.h"
00057 #include "MdcFastTrkAlg/FTLayer.h"
00058 #include "MdcFastTrkAlg/FTSuperLayer.h"
00059 #include "MdcFastTrkAlg/FTSegment.h"
00060 #include "MdcFastTrkAlg/FTTrack.h"
00061 #include "MdcFastTrkAlg/MdcParameter.h"
00062
00063 #ifndef OnlineMode
00064 #include "McTruth/McParticle.h"
00065
00066 #include "AIDA/IAxis.h"
00067 #include "AIDA/IHistogram1D.h"
00068 #include "AIDA/IHistogram2D.h"
00069 #include "AIDA/IHistogram3D.h"
00070 #include "AIDA/IHistogramFactory.h"
00071 #endif
00072
00073
00074 using namespace Event;
00075
00076
00077
00078
00079 #ifndef OnlineMode
00080
00081 extern NTuple::Item<long> g_ntrkMC, g_eventNo;
00082 extern NTuple::Array<float> g_theta0MC, g_phi0MC;
00083 extern NTuple::Array<float> g_pxMC, g_pyMC, g_pzMC, g_ptMC;
00084
00085
00086
00087 extern NTuple::Item<long> g_ntrk;
00088 extern NTuple::Array<float> g_px, g_py, g_pz, g_pt, g_p;
00089 extern NTuple::Array<float> g_phi, g_theta;
00090 extern NTuple::Array<float> g_vx, g_vy, g_vz;
00091 extern NTuple::Array<float> g_dr, g_phi0, g_kappa, g_dz, g_tanl;
00092 extern NTuple::Item<float> g_estime;
00093
00094 extern IHistogram1D* g_ncellMC;
00095 extern IHistogram1D* g_ncell;
00096 extern IHistogram1D* g_naxialhit;
00097 extern IHistogram1D* g_nstereohit;
00098 extern IHistogram1D* g_nhit;
00099 extern IHistogram2D* g_hitmap[20];
00100
00101 extern int num_2Dtrk, num_3Dtrk, num_finaltrk;
00102
00103 #endif
00104
00105
00106
00107
00108
00109
00110
00111
00112 MdcParameter* FTFinder::param = MdcParameter::instance();
00113
00114
00115
00116
00117
00118
00119
00120
00121 FTFinder::FTFinder()
00122 :
00123
00124
00125
00126
00127
00128
00129
00130 #ifndef OnlineMode
00131
00132 #endif
00133
00134 tOffSet(0.),
00135 t0Estime(-999.),
00136 tEstFlag(0),
00137
00138 _wire(NULL),
00139 _layer(NULL),
00140 _superLayer(NULL),
00141 _tracks(*(new FTList<FTTrack *>(20))),
00142 _linkedSegments(new FTList<FTSegment *>(6)),
00143 _axialSegCollect(10),
00144 _vx(0.),
00145 _vy(0.),
00146 _vz(0.),
00147 _ExpNo(0),
00148 _RunNo(0),
00149 m_total_trk(0),
00150 pivot(0,0,0)
00151 {
00152 StatusCode scmgn = Gaudi::svcLocator()->service("MagneticFieldSvc",m_pmgnIMF);
00153 if(scmgn!=StatusCode::SUCCESS) {
00154 std::cout<< "Unable to open Magnetic field service"<<std::endl;
00155 }
00156 }
00157
00158 void
00159 FTFinder::init()
00160 {
00161 if (!param->_doIt) return;
00162
00163
00164
00165 IPartPropSvc* p_PartPropSvc;
00166 static const bool CREATEIFNOTTHERE(true);
00167 StatusCode PartPropStatus = Gaudi::svcLocator()->service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
00168 if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) {
00169
00170 std::cerr << "Could not initialize Particle Properties Service" << std::endl;
00171
00172 return;
00173 }
00174 m_particleTable = p_PartPropSvc->PDT();
00175
00176
00177 StatusCode RawData = Gaudi::svcLocator()-> service ("RawDataProviderSvc", m_rawDataProviderSvc, CREATEIFNOTTHERE);
00178 if ( !RawData.isSuccess() ){
00179 std::cerr << "Could not load RawDataProviderSvc!" << m_rawDataProviderSvc << endreq;
00180 return;
00181 }
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316 }
00317
00318 void
00319 FTFinder::term()
00320 {
00321 if (param->_doIt) clear();
00322 delete &_tracks;
00323 delete _linkedSegments;
00324 if (!param->_doIt) return;
00325 free(_wire);
00326 free(_layer);
00327 free(_superLayer);
00328 }
00329
00330 void
00331 FTFinder::begin_run(){
00332
00333 eventNo = 0;
00334 runNo=0;
00335 expflag=0;
00336 if (!param->_doIt) return;
00337 _ExpNo = 0;
00338 _RunNo = 0;
00339
00340 IMdcGeomSvc* mdcGeomSvc;
00341 StatusCode sc = Gaudi::svcLocator()->service("MdcGeomSvc", mdcGeomSvc);
00342
00343 if (sc == StatusCode::SUCCESS) {
00344
00345
00346 } else {
00347 return ;
00348 }
00349
00350 const int Nwire = mdcGeomSvc->getWireSize();
00351 const int Nlyr = mdcGeomSvc->getLayerSize();
00352 const int Nsup = mdcGeomSvc->getSuperLayerSize();
00353
00354 if(!Nwire || !Nlyr){
00355 std::cerr << "FTFINDER::GEOCDC not found (please put cdctable before l4)" << std::endl;
00356 std::cerr << "JOB will stop" << std::endl;
00357 exit(-1);
00358 }
00359
00360 if (!_wire) _wire = (FTWire *) malloc((Nwire+1) * sizeof(FTWire));
00361 if (!_layer) _layer = (FTLayer *) malloc(Nlyr * sizeof(FTLayer));
00362 if (!_superLayer) _superLayer =
00363 (FTSuperLayer *) malloc(Nsup * sizeof(FTSuperLayer));
00364
00365 if (!_wire || !_layer || !_superLayer){
00366 std::cerr << "FTFINDER::Cannot allocate geometries" << std::endl;
00367 std::cerr << "JOB will stop" << std::endl;
00368 exit(-1);
00369 }
00370
00371 int superLayerID = 0;
00372 int layerID = 0;
00373 int localLayerID = 0;
00374 int localWireID = 0;
00375 int localID = 0;
00376 int wireID;
00377 MdcGeoLayer * layer_back = NULL;
00378 MdcGeoSuper * superLayer_back = NULL;
00379 int k = 0;
00380 int Nlayer[12];
00381 int Nlocal[12];
00382 int NlocalWireID[44];
00383
00384 for (wireID = 0;wireID <= Nwire; wireID++){
00385 MdcGeoLayer * layer = (wireID==Nwire) ? NULL : mdcGeomSvc->Wire(wireID)->Lyr();
00386 if (layer != layer_back){
00387 layer_back = layer;
00388 MdcGeoSuper * superLayer = (wireID==Nwire) ? NULL : mdcGeomSvc->Layer(layerID)->Sup();
00389 if (superLayer != superLayer_back){
00390 superLayer_back = superLayer;
00391 Nlayer[k] = localLayerID;
00392 Nlocal[k] = localID;
00393 localLayerID = 0;
00394 k++;
00395 }
00396 NlocalWireID[layerID] = localWireID;
00397 localID = 0;
00398 localWireID = 0;
00399 layerID++;
00400 localLayerID++;
00401 }
00402 localID++;
00403 localWireID++;
00404 }
00405
00406
00407 superLayerID = -1;
00408 layerID = -1;
00409 localLayerID = 0;
00410 localID = 0;
00411 layer_back = NULL;
00412 superLayer_back = NULL;
00413 for (wireID = 0;wireID < Nwire; wireID++){
00414 MdcGeoLayer * layer = mdcGeomSvc->Wire(wireID)->Lyr();
00415 if (layer != layer_back){
00416 layer_back = layer;
00417 MdcGeoSuper * superLayer = mdcGeomSvc->Layer(layerID+1)->Sup();
00418 if (superLayer != superLayer_back){
00419 superLayer_back = superLayer;
00420
00421 superLayerID++;
00422 new(_superLayer+superLayerID) FTSuperLayer(wireID,
00423 Nlocal[superLayerID+1],
00424 layerID+1,
00425 Nlayer[superLayerID+1],
00426 superLayerID);
00427 localLayerID=0;
00428 }
00429
00430 layerID++;
00431 double slantWithSymbol = (mdcGeomSvc->Layer(layerID)->Slant())*(mdcGeomSvc->Layer(layerID)->Sup()->Type());
00432
00433 new(_layer+layerID) FTLayer(0.1*layer->Radius(), layer->Slant(),
00434 0.1*(+layer->Length()/2), 0.1*(-layer->Length()/2), 0.1*layer->Offset(),
00435 layerID, localLayerID++, NlocalWireID[layerID+1],
00436 _superLayer[superLayerID]);
00437 localID = 0;
00438 }
00439
00440 const MdcGeoWire * wire = mdcGeomSvc->Wire(wireID);
00441 if (superLayerID == 2 || superLayerID == 3 ||
00442 superLayerID == 4 || superLayerID == 9 ||
00443 superLayerID == 10){
00444 new(_wire+wireID) FTWire(0.1*(float)0.5*(wire->Backward().x()+wire->Forward().x()),
00445 0.1*(float)0.5*(wire->Backward().y()+wire->Forward().y()),
00446 0.,0.,
00447 _layer[layerID], localID++, _wire+Nwire);
00448
00449 }else{
00450 new(_wire+wireID) FTWire(0.1*(float)wire->Backward().x(),
00451 0.1*(float)wire->Backward().y(),
00452 0.1*(float)wire->Forward().x() - 0.1*(float)wire->Backward().x(),
00453 0.1*(float)wire->Forward().y() - 0.1*(float)wire->Backward().y(),
00454 _layer[layerID], localID++, _wire+Nwire);
00455 }
00456 }
00457
00458
00459 new(_wire+Nwire) FTWire();
00460
00461 for (int i = 0; i < Nwire; i++){
00462 (*(_wire+i)).initNeighbor();
00463 }
00464
00465 _widbase[0] = 0;
00466 for (int k = 1; k < 43; k++){
00467 _widbase[k] = _widbase[k-1] + (*(_layer+k-1)).NWire();
00468 }
00469 }
00470
00471 void
00472 FTFinder::event() {
00473 IMessageSvc *msgSvc;
00474 Gaudi::svcLocator()->service("MessageSvc", msgSvc);
00475
00476 MsgStream log(msgSvc, "FTFinder");
00477
00478 IDataProviderSvc* eventSvc = NULL;
00479 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
00480
00481 if (!param->_doIt) return;
00482
00483
00484
00485 clear();
00486
00487
00488 DataObject *aReconEvent;
00489 eventSvc->findObject("/Event/Recon",aReconEvent);
00490 if(!aReconEvent) {
00491
00492 ReconEvent* recevt = new ReconEvent;
00493 StatusCode sc = eventSvc->registerObject("/Event/Recon",recevt );
00494 if(sc!=StatusCode::SUCCESS) {
00495 log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
00496 return;
00497 }
00498 }
00499
00500 IDataManagerSvc *dataManSvc = dynamic_cast<IDataManagerSvc*> (eventSvc);
00501 DataObject *aEsTimeEvent;
00502 eventSvc->findObject("/Event/Recon/RecEsTimeCol",aEsTimeEvent);
00503 if(aEsTimeEvent != NULL) {
00504 dataManSvc->clearSubTree("/Event/Recon/RecEsTimeCol");
00505 eventSvc->unregisterObject("/Event/Recon/RecEsTimeCol");
00506 }
00507
00508 RecEsTimeCol *aRecEsTimeCol = new RecEsTimeCol;
00509 StatusCode est;
00510 est = eventSvc->registerObject("/Event/Recon/RecEsTimeCol",aRecEsTimeCol);
00511 if( est.isFailure() ) {
00512 log << MSG::FATAL << "Could not register RecEsTimeCol" << endreq;
00513 return;
00514 }
00515 log << MSG::DEBUG << "RecEsTimeCol registered successfully!" <<endreq;
00516
00517
00518
00519 updateMdc();
00520
00521
00522
00523
00524 for (int i = 0; i^11; i++){
00525 (*(_superLayer+i)).mkSegmentList();
00526 }
00527
00528
00529
00530 for(int j=0;j<10;j++){
00531 (*(_superLayer+j)).reduce_noise(tEstime);
00532 }
00533
00534 getTestime();
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544 mkTrackList();
00545
00546
00547
00548
00549
00550
00551 int n = _tracks.length();
00552 log << MSG::DEBUG << "number of 2D tracks: " << n <<endreq;
00553 #ifndef OnlineMode
00554 num_2Dtrk+=n;
00555 #endif
00556 for (int i = 0; i^n; i++){
00557 if (!_tracks[i]->r_phiFit()){
00558 delete _tracks[i];
00559 log<<MSG::DEBUG<<"===========>deleted 2D track : "<< i+1 <<endreq;
00560 n = _tracks.remove(i);
00561 }
00562 }
00563
00564 if(t0Estime!=-999){
00565
00566 RecEsTime *arecestime = new RecEsTime;
00567 arecestime -> setTest(t0Estime);
00568 arecestime -> setStat(tEstFlag);
00569 aRecEsTimeCol->push_back(arecestime);
00570 }
00571 if (!n){
00572 makeTds();
00573 return;
00574 }
00575
00578
00579 int vtx_flag = VertexFit(0);
00580 evtTiming = (param->_evtTimeCorr) ? CorrectEvtTiming() : 0;
00581
00584
00585 if(param->_hitscut==1){
00586 for (int i = 0; i^n; i++){
00587 _tracks[i]->r_phiReFit(_vx, _vy, vtx_flag);
00588 }
00589 }
00590
00591
00592
00593 if(param->_hitscut==2){
00594 for (int i = 0; i^n; i++){
00595 for(int j = 0; j<2; j++){
00596 i_rPhiFit= _tracks[i]->r_phi2Fit(_vx, _vy, vtx_flag);
00597 if(i_rPhiFit!=-99) _tracks[i]->r_phi3Fit(i_rPhiFit,_vx, _vy, vtx_flag);
00598 }
00599 _tracks[i]->r_phi4Fit(_vx, _vy, vtx_flag);
00600 }
00601 }
00602
00603
00604
00605
00606 mkTrack3D();
00607
00608
00609
00610
00611 n = _tracks.length();
00612 log<< MSG::DEBUG <<"number of 3D tracks: " << n << endreq;
00613
00614 #ifndef OnlineMode
00615 num_3Dtrk+=n;
00616 #endif
00617
00618 for (int i = 0; i^n; i++) {
00619
00620 #ifndef OnlineMode
00621 log<<MSG::DEBUG<<"=======>3D track: "<<i<<endreq;
00622 _tracks[i]->printout();
00623 #endif
00624
00625 if (_tracks[i]->get_nhits() < 18) {
00626 delete _tracks[i];
00627 log<< MSG::DEBUG <<"================>deleted 3D track : "<< i+1 <<endreq;
00628 n = _tracks.remove(i);
00629 }
00630 }
00631
00632 if (!n){
00633 makeTds();
00634 return;
00635 }
00636
00637 for (int i = 0; i^n; i++){
00638 _tracks[i]->s_zFit();
00639 }
00640
00641
00642
00643
00644 if (param->_findEventVertex){
00645 VertexFit(1);
00646 }
00647
00648 n = _tracks.length();
00649 log<<MSG::DEBUG<<"final number of tracks: " << n << endreq;
00650
00651 #ifndef OnlineMode
00652 num_finaltrk+=n;
00653 if (param->_mkMdst) makeMdst();
00654 #endif
00655
00656
00657
00658
00659
00660 if (param->_mkTds) makeTds();
00661
00662 }
00663
00664 int
00665 FTFinder::updateMdc(void){
00666 IMessageSvc *msgSvc;
00667 Gaudi::svcLocator()->service("MessageSvc", msgSvc);
00668
00669 MsgStream log(msgSvc, "FTFinder");
00670
00671 IDataProviderSvc* eventSvc = NULL;
00672 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
00673
00674
00675 if (eventSvc) {
00676
00677 } else {
00678 log << MSG::FATAL << "Could not find Event Header" << endreq;
00679 return 0;
00680 }
00681
00682
00683
00684
00685 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc,"/Event/EventHeader");
00686 if (!eventHeader) {
00687 log << MSG::FATAL << "Could not find Event Header" << endreq;
00688 return( StatusCode::FAILURE);
00689 }
00690 log << MSG::DEBUG << "MdcFastTrkAlg: retrieved event: " << eventHeader->eventNumber() << " run: " << eventHeader->runNumber() << endreq;
00691 eventNo = eventHeader->eventNumber();
00692 runNo = eventHeader->runNumber();
00693 if(runNo>0) expflag=1;
00694 else expflag=0;
00695 int digiId;
00696
00697 #ifndef OnlineMode
00698 g_eventNo = eventHeader->eventNumber();
00699
00700
00701 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc,"/Event/MC/McParticleCol");
00702 if (!mcParticleCol) {
00703 log << MSG::WARNING << "Could not find McParticle" << endreq;
00704
00705 }
00706 else{
00707 McParticleCol ::iterator iter_mc = mcParticleCol->begin();
00708 digiId = 0;
00709 int ntrkMC = 0;
00710 int charge = 0;
00711 for (;iter_mc != mcParticleCol->end(); iter_mc++, digiId++) {
00712 log << MSG::DEBUG << "MDC digit No: " << digiId << endreq;
00713 log << MSG::DEBUG
00714 << " particleId = " << (*iter_mc)->particleProperty ()
00715 << endreq;
00716 int statusFlags = (*iter_mc)->statusFlags();
00717 int pid = (*iter_mc)->particleProperty();
00718 if(pid >0) {
00719 if(m_particleTable->particle( pid )){
00720 charge = (int)m_particleTable->particle( pid )->charge();
00721 }
00722 } else if ( pid <0 ) {
00723 if(m_particleTable->particle( -pid )){
00724 charge = (int)m_particleTable->particle( -pid )->charge();
00725 charge *= -1;
00726 }
00727 } else {
00728 log << MSG::WARNING << "wrong particle id, please check data" <<endreq;
00729 }
00730
00731 if(charge==0 || abs(cos((*iter_mc)->initialFourMomentum().theta()))>0.93) continue;
00732
00733 g_theta0MC[ntrkMC] = (*iter_mc)->initialFourMomentum().theta();
00734 g_phi0MC[ntrkMC] = (*iter_mc)->initialFourMomentum().phi();
00735 g_pxMC[ntrkMC] = (*iter_mc)->initialFourMomentum().px();
00736 g_pyMC[ntrkMC] = (*iter_mc)->initialFourMomentum().py();
00737 g_pzMC[ntrkMC] = (*iter_mc)->initialFourMomentum().pz();
00738 g_ptMC[ntrkMC] = 0.001 * sqrt(g_pxMC[ntrkMC]*g_pxMC[ntrkMC] + g_pyMC[ntrkMC]*g_pyMC[ntrkMC]);
00739 ntrkMC++;
00740 }
00741 g_ntrkMC = ntrkMC;
00742 }
00743 #endif
00744
00745 #ifdef MultiThread
00746
00747 SmartDataPtr<MdcDigiCol> mdcDigiVec(eventSvc,"/Event/Digi/MdcDigiCol");
00748 if (!mdcDigiVec) {
00749 log << MSG::FATAL << "Could not find MDC digi" << endreq;
00750 return( StatusCode::FAILURE);
00751 }
00752 MdcDigiCol::iterator iter1 = mdcDigiVec->begin();
00753 digiId = 0;
00754 Identifier mdcId;
00755 int layerId;
00756 int wireId;
00757 for (;iter1 != mdcDigiVec->end(); iter1++, digiId++) {
00758 #else
00759
00760
00761 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec();
00762 MdcDigiVec::iterator iter1 = mdcDigiVec.begin();
00763 digiId = 0;
00764 Identifier mdcId;
00765 int layerId;
00766 int wireId;
00767 for (;iter1 != mdcDigiVec.end(); iter1++, digiId++) {
00768 #endif
00769 mdcId = (*iter1)->identify();
00770 layerId = MdcID::layer(mdcId);
00771 wireId = MdcID::wire(mdcId);
00772 #ifndef OnlineMode
00773 log << MSG::DEBUG << "MDC digit No: " << digiId << endreq;
00774 log << MSG::DEBUG
00775 << " time_channel = " << RawDataUtil::MdcTime((*iter1)->getTimeChannel())
00776 << " charge_channel = " << RawDataUtil::MdcCharge((*iter1)->getChargeChannel())
00777 << " layerId = " << layerId
00778 << " wireId = " << wireId
00779 << endreq;
00780 #endif
00781
00782 const int localwid = wireId;
00783 const int wid = localwid + _widbase[layerId];
00784
00785 #ifndef OnlineMode
00786 g_ncellMC->fill((float)wid, 1.0);
00787 #endif
00788
00789 if (wid < 0 || wid > 6795) continue;
00790 FTWire & w = *(_wire + wid);
00791 float tdc = RawDataUtil::MdcTime((*iter1)->getTimeChannel());
00792 float adc = RawDataUtil::MdcCharge((*iter1)->getChargeChannel());
00793
00794 #ifndef OnlineMode
00795 if (eventNo == 466) {
00796 g_hitmap[0]->fill(w.x(), w.y());
00797 }
00798
00799 if (eventNo == 538) {
00800 g_hitmap[1]->fill(w.x(), w.y());
00801 }
00802
00803 if (eventNo == 408) {
00804 g_hitmap[2]->fill(w.x(), w.y());
00805 }
00806
00807 if (eventNo == 499) {
00808 g_hitmap[3]->fill(w.x(), w.y());
00809 }
00810
00811 if (eventNo == 603) {
00812 g_hitmap[4]->fill(w.x(), w.y());
00813 }
00814
00815 if (eventNo == 938) {
00816 g_hitmap[5]->fill(w.x(), w.y());
00817 }
00818 #endif
00819
00820
00821 float time = tdc-sqrt(w.x()*w.x()+w.y()*w.y())/30;
00822 w.time(time);
00823 w.wireId(wid);
00824 w.setAdc(adc);
00825
00826
00827 const FTLayer & lyr = w.layer();
00828
00829
00830
00831 w.distance(0.0);
00832 w.setChi2(999);
00833 w.state(FTWireHit);
00834 FTSuperLayer & sup = (FTSuperLayer &) lyr.superLayer();
00835 sup.appendHit(&w);
00836 }
00837
00838 return 1;
00839 }
00840
00841
00842 void
00843 FTFinder::clear(void){
00844 evtTiming = 0;
00845
00846 if ( _superLayer ) {
00847 for (int i = 0; i^11; i++) (*(_superLayer+i)).clear();
00848 }
00849 _tracks.deleteAll();
00850 _axialSegCollect.deleteAll();
00851 _vx = -99999.;
00852 _vy = -99999.;
00853 _vz = -99999.;
00854 for (int i = 0; i < 10; i++) {
00855 tEstime[i].clear();
00856 }
00857 }
00858
00859 int
00860 FTFinder::getTestime(void)
00861 {
00862 IMessageSvc *msgSvc;
00863 Gaudi::svcLocator()->service("MessageSvc", msgSvc);
00864
00865 MsgStream log(msgSvc, "FTFinder");
00866
00867 IDataProviderSvc* eventSvc = NULL;
00868 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
00869
00870 float sumT=0,estime=0;
00871 int n=0;
00872 t0Estime=-999;
00873
00874 for(int i=0;i<10;i++){
00875 for(int j=0;j<tEstime[i].length();j++){
00876 if(tEstime[i][j]!=0){
00877 sumT+=tEstime[i][j];
00878 n++;
00879 }
00880 }
00881 }
00882 if(n!=0){
00883 estime=sumT/n;
00884 estime+=_t0cal;
00885 if(estime>-1){
00886
00887 if(expflag==0){
00888 int nbunch=((int)(estime-tOffSet))/_bunchtime;
00889 if(((int)(estime-tOffSet))%(int)_bunchtime>_bunchtime/2) nbunch=nbunch+1;
00890 t0Estime=nbunch*_bunchtime+tOffSet;
00891
00892 }
00893 else{
00894 t0Estime=estime;
00895 }
00896
00897
00898
00899
00900
00901
00902 int trigtimming=0;
00903
00904
00905 SmartDataPtr<TrigData> trigData(eventSvc,"/Event/Trig/TrigData");
00906 if (trigData) {
00907 trigtimming=trigData->getTimingType();
00908 log << MSG::INFO <<"Timing type: "<< trigData->getTimingType() << endreq;
00909 }
00910 if(trigtimming==1) tEstFlag=117;
00911 if(trigtimming==2) tEstFlag=127;
00912 if(trigtimming==3) tEstFlag=137;
00913 if(trigtimming==0) tEstFlag=107;
00914 #ifndef OnlineMode
00915 g_estime=estime;
00916 #endif
00917 }
00918 }
00919 }
00920
00921 void
00922 FTFinder::mkTrackList(void)
00923 {
00924
00925 IMessageSvc *msgSvc;
00926 Gaudi::svcLocator()->service("MessageSvc", msgSvc);
00927
00928 MsgStream log(msgSvc, "FTFinder");
00929
00930 FTTrack ** tracks = (FTTrack **)malloc( 6 * sizeof(FTTrack *));
00931 int * Nsame = (int *)malloc( 6 * sizeof(int));
00932 int ntrack_cand = 0;
00933 int i,j,k,l;
00934
00935
00936
00937
00938
00939
00940 FTList<FTSegment *> inner_segments(10);
00941 FTList<FTSegment *> outer_segments(10);
00942
00943 linkAxialSuperLayer234(inner_segments);
00944 linkAxialSuperLayer910(outer_segments);
00945
00946 _axialSegCollect.append(inner_segments);
00947 _axialSegCollect.append(outer_segments);
00948
00949 int innerN = inner_segments.length();
00950 int outerN = outer_segments.length();
00951
00952 log << MSG::DEBUG << innerN <<" segments in inner axial layers!"<<endreq;
00953 log << MSG::DEBUG << outerN <<" segments in outer axial layers!"<<endreq;
00954
00955 for (l = 0; l^innerN; l++){
00956 if (inner_segments[l]->track()) continue;
00957 FTSegment *inner = inner_segments[l];
00958 FTTrack *track_cand = NULL;
00959
00960 #ifndef OnlineMode
00961 inner->printout();
00962 #endif
00963
00964 if(outerN==0){
00965 track_cand = linkAxialSegments(&inner, NULL);
00966 if(!track_cand) continue;
00967 _linkedSegments = new FTList<FTSegment *>(5);
00968
00969
00970 FTList<FTSegment *>& segments = track_cand->axial_segments();
00971 int n=segments.length();
00972 for (i = 0; i < n; i++){
00973 segments[i]->track(track_cand);
00974 }
00975 _tracks.append(track_cand);
00976 track_cand->setFTFinder(this);
00977 ntrack_cand++;
00978 continue;
00979 }
00980
00981 for(k = 0; k^outerN; k++) {
00982 if (outer_segments[k]->track()) {
00983 if(k==outerN-1) {
00984 if (inner_segments[l]->track()) continue;
00985 track_cand = linkAxialSegments(&inner, NULL);
00986 if(!track_cand) continue;
00987 _linkedSegments = new FTList<FTSegment *>(5);
00988
00989 FTList<FTSegment *>& segments = track_cand->axial_segments();
00990 int n=segments.length();
00991 for (i = 0; i < n; i++){
00992 segments[i]->track(track_cand);
00993 }
00994 _tracks.append(track_cand);
00995 track_cand->setFTFinder(this);
00996 ntrack_cand++;
00997 }
00998 continue;
00999 }
01000
01001 FTSegment *outer = outer_segments[k];
01002 track_cand = linkAxialSegments(&inner, &outer);
01003
01004 if (!track_cand ){
01005 if( k == outerN-1 ) {
01006 if (inner_segments[l]->track()) continue;
01007 track_cand = linkAxialSegments(&inner, NULL);
01008 }
01009 if(!track_cand) continue;
01010 }
01011
01012
01013
01014
01015
01016
01017 ntrack_cand++;
01018 _linkedSegments = new FTList<FTSegment *>(5);
01019
01020
01021
01022 FTList<FTSegment *>& segments = track_cand->axial_segments();
01023 int n = segments.length();
01024
01025
01026 if(inner_segments[l]->track()) {
01027 FTTrack * trk_tmp = inner_segments[l]->track();
01028 int nTrks = _tracks.length();
01029 int cache_i=0;
01030 for (i = 0; i^nTrks; i++){
01031 if (_tracks[i] == trk_tmp){
01032 cache_i = i;
01033 break;
01034 }
01035 }
01036 FTList<FTSegment *> & segments2 = trk_tmp->axial_segments();
01037 int n2 = segments2.length();
01038 if(n >= n2 && track_cand->chi2_kappa_tmp() < trk_tmp->chi2_kappa_tmp()){
01039
01040 for (i = 0; i < n2; i++){
01041 segments2[i]->track(NULL);
01042 }
01043 for (i = 0; i < n; i++){
01044 segments[i]->track(track_cand);
01045 }
01046 _tracks.replace(cache_i,track_cand);
01047 track_cand->setFTFinder(this);
01048 delete trk_tmp;
01049 }else{
01050 delete track_cand;
01051 }
01052 }else{
01053
01054 for (i = 0; i < n; i++){
01055 segments[i]->track(track_cand);
01056 }
01057 _tracks.append(track_cand);
01058 track_cand->setFTFinder(this);
01059 }
01060 }
01061 }
01062
01063 free(tracks);
01064 free(Nsame);
01065
01066 #ifndef OnlineMode
01067 int n=_tracks.length();
01068 for(int i=0; i<n; i++){
01069 FTList<FTSegment *> & segments = _tracks[i]->axial_segments();
01070 log << MSG::DEBUG <<" loop connected axial track: " << i <<endreq;
01071 int l=segments.length();
01072 for(int j=0; j<l; j++){
01073 segments[j]->printout();
01074 }
01075 }
01076 #endif
01077
01078 }
01079
01080 FTTrack *
01081 FTFinder::linkAxialSegments(FTSegment **inner, FTSegment ** outer)
01082 {
01083 float chi2_kappa = param->_chi2_kappa;
01084 _linkedSegments->clear();
01085 if(!outer) {
01086
01087 int n = (*inner)->wireHits().length();
01088 _linkedSegments->append(*inner);
01089 if(n >= 7) return (new FTTrack(*_linkedSegments, (*inner)->kappa(), chi2_kappa));
01090 else return NULL;
01091 }
01092 int n = _linkedSegments->append(*outer);
01093 float SigmaK = (*outer)->kappa();
01094 float SigmaRR = (*outer)->r();
01095 SigmaRR *= SigmaRR;
01096 float SigmaKRR = SigmaK*SigmaRR;
01097 float SigmaKKRR = SigmaK*SigmaKRR;
01098 FTSegment & s = * (*_linkedSegments)[n-1];
01099 FTSegment * innerSegment = NULL;
01100 float SigmaK_cache, SigmaRR_cache, SigmaKRR_cache, SigmaKKRR_cache;
01101
01102 float Min_chi2 = param->_Min_chi2;
01103 float inX = s.incomingX();
01104 float inY = s.incomingY();
01105
01106
01107
01108 float in_r = s.innerBoundHits().first()->layer().r();
01109 float incomingPhi = s.incomingPhi();
01110
01111 FTSegment* next = *inner;
01112 const FTWire & NextOuterBoundHit = * next->outerBoundHits().first();
01113
01114
01115
01116 float outX = next->outgoingX();
01117 float outY = next->outgoingY();
01118
01120 float _trk_d = -2*(-1. / 2.99792458 /m_pmgnIMF->getReferField())/s.kappa();
01121 float _angle1 = asin(NextOuterBoundHit.layer().r()/_trk_d);
01122 float _angle2 = asin(s.outerBoundHits().first()->layer().r()/_trk_d);
01123 float _ang_diff = _angle2 - _angle1;
01124 float _diff = s.outgoingPhi() - next->outgoingPhi();
01125 _diff = _diff - (int(_diff/M_PI))*2*M_PI;
01126 float _require = _ang_diff - _diff;
01127
01128 if (_require < -0.10 || _require > 0.11) return NULL;
01129
01130 float SegK = next->kappa();
01131 float SegRR = next->r();
01132 SegRR *= SegRR;
01133 const float out_r = NextOuterBoundHit.layer().r();
01134
01135 float GapK = 2.*(-1. / 2.99792458 /m_pmgnIMF->getReferField())*(inX*outY-inY*outX) /
01136 (in_r*out_r*sqrt((inX-outX)*(inX-outX)+(inY-outY)*(inY-outY)));
01137
01138 float GapRR = 0.5*(in_r+out_r);
01139 GapRR*=GapRR;
01140 float SigmaK_tmp = (SigmaK + SegK + GapK);
01141 float SigmaRR_tmp = SigmaRR + SegRR + GapRR;
01142 float SigmaKRR_tmp = SigmaKRR + SegK*SegRR + GapK*GapRR;
01143 float SigmaKKRR_tmp = SigmaKKRR + SegK*SegK*SegRR + GapK*GapK*GapRR;
01144 float MuK_tmp = SigmaK_tmp/(2*n+1);
01145 float chi2 = (MuK_tmp*MuK_tmp*SigmaRR_tmp
01146 - 2.*MuK_tmp*SigmaKRR_tmp + SigmaKKRR_tmp)/(2*n+1);
01147 if ((chi2-chi2_kappa) < Min_chi2){
01148 Min_chi2 = chi2;
01149 innerSegment = next;
01150 SigmaK_cache = SigmaK_tmp;
01151 SigmaRR_cache = SigmaRR_tmp;
01152 SigmaKRR_cache = SigmaKRR_tmp;
01153 SigmaKKRR_cache = SigmaKKRR_tmp;
01154 }
01155 if (innerSegment){
01156 n = _linkedSegments->append(*inner);
01157 SigmaK = SigmaK_cache;
01158 SigmaRR = SigmaRR_cache;
01159 SigmaKRR = SigmaKRR_cache;
01160 SigmaKKRR = SigmaKKRR_cache;
01161 chi2_kappa = Min_chi2;
01162
01163 if (n > 1) return (new FTTrack(*_linkedSegments,SigmaK/(2*n-1),chi2_kappa));
01164 }else{
01165
01166
01167
01168
01169
01170 return NULL;
01171 }
01172
01173 return NULL;
01174
01175 }
01176
01177 void
01178 FTFinder::linkAxialSuperLayer234(FTList<FTSegment *> & inner_segments)
01179 {
01180 FTList<FTSegment *> _segments34(10);
01181 FTList<FTSegment *>& SuperLayer3Segments = _superLayer[3].segments();
01182 FTList<FTSegment *>& SuperLayer4Segments = _superLayer[4].segments();
01183 linkAxialSegments_step(SuperLayer3Segments, SuperLayer4Segments,
01184 _segments34, param->_D_phi2, param->_chi2_2);
01185 _segments34.append(SuperLayer3Segments);
01186 SuperLayer3Segments.clear();
01187
01188 FTList<FTSegment *>& SuperLayer2Segments = _superLayer[2].segments();
01189 linkAxialSegments_step(SuperLayer2Segments, _segments34,
01190 inner_segments, param->_D_phi1, param->_chi2_1);
01191 inner_segments.append(_segments34);
01192
01193
01194 int n = SuperLayer2Segments.length();
01195 int m = SuperLayer4Segments.length();
01196 for (int i = 0; i^n; i++) {
01197 FTSegment * inner = SuperLayer2Segments[i];
01198
01199 float in_outerPhi = inner->outgoingPhi();
01200 for (int j = 0; j^m; j++) {
01201 FTSegment * outer = SuperLayer4Segments[j];
01202
01203 float out_innerPhi = outer->incomingPhi();
01204 float D_phi = fabs(in_outerPhi - out_innerPhi);
01205 if (D_phi > M_PI) D_phi = 2*M_PI - D_phi;
01207 if (D_phi > M_PI/12.5) continue;
01208 float inX = inner->outgoingX();
01209 float inY = inner->outgoingY();
01210 float outX = outer->incomingX();
01211 float outY = outer->incomingY();
01212 float in_r = inner->outerBoundHits().first()->layer().r();
01213 float out_r = outer->innerBoundHits().first()->layer().r();
01214 float GapK = 2.*(-1. / 2.99792458 /m_pmgnIMF->getReferField())*(inY*outX-inX*outY) /
01215 (in_r*out_r*sqrt((inX-outX)*(inX-outX)+(inY-outY)*(inY-outY)));
01216
01217 float cache_in = ((outer->kappa()+GapK)/3.0 - inner->kappa()*2.0/3.0)/in_r;
01218 float cache_out = ((inner->kappa()+GapK)/3.0 - outer->kappa()*2.0/3.0)/out_r;
01219 float cache_gap = ((inner->kappa()+outer->kappa())/3.0-GapK*2.0/3.0)*2.0/(in_r+out_r);
01220 float chi2_z = cache_in*cache_in + cache_out*cache_out + cache_gap*cache_gap;
01221 if (chi2_z > param->_D_phi1) continue;
01223 inner->connect_outer(outer->wireHits(),outer->outerBoundHits());
01224 inner->update();
01225 inner_segments.append(inner);
01226 n = SuperLayer2Segments.remove(i);
01227 m = SuperLayer4Segments.remove(j);
01228 delete outer;
01229 break;
01230 }
01231 }
01232 inner_segments.append(SuperLayer2Segments);
01233 SuperLayer2Segments.clear();
01234 inner_segments.append(SuperLayer4Segments);
01235 SuperLayer4Segments.clear();
01236 }
01237
01238 void
01239 FTFinder::linkAxialSuperLayer910(FTList<FTSegment *> & outer_segments)
01240 {
01241 FTList<FTSegment *>& SuperLayer9Segments = _superLayer[9].segments();
01242 FTList<FTSegment *>& SuperLayer10Segments = _superLayer[10].segments();
01243 linkAxialSegments_step(SuperLayer9Segments, SuperLayer10Segments,
01244 outer_segments, param->_D_phi3, param->_chi2_3);
01245 outer_segments.append(SuperLayer9Segments);
01246 SuperLayer9Segments.clear();
01247 outer_segments.append(SuperLayer10Segments);
01248 SuperLayer10Segments.clear();
01249 }
01250
01251 void
01252 FTFinder::linkAxialSegments_step(FTList<FTSegment *>& InnerSegments,
01253 FTList<FTSegment *>& OuterSegments,
01254 FTList<FTSegment *>& Connected,
01255 float maxDphi, float maxChi2)
01256 {
01257 IMessageSvc *msgSvc;
01258 Gaudi::svcLocator()->service("MessageSvc", msgSvc);
01259 MsgStream log(msgSvc, "FTFinder");
01260
01261 int n = InnerSegments.length();
01262 int m = OuterSegments.length();
01263 for (int i = 0; i^n; i++){
01264 FTSegment * inner = InnerSegments[i];
01265 #ifndef OnlineMode
01266
01267
01268 #endif
01269 const FTLayer & in_outerBound = inner->outerBoundHits().first()->layer();
01270 float in_outerPhi = inner->outgoingPhi();
01271 float min_Dphi = M_PI/2;
01272 int min_Dphi_index = -1;
01273 for (int j = 0; j^m; j++) {
01274 FTSegment * outer = OuterSegments[j];
01275 #ifndef OnlineMode
01276
01277
01278 #endif
01279 float D_phi = fabs(in_outerPhi - outer->incomingPhi());
01280 if (D_phi > M_PI) D_phi = 2*M_PI - D_phi;
01281 if (D_phi > min_Dphi) continue;
01282 float inX = inner->incomingX();
01283 float inY = inner->incomingY();
01284 float outX = outer->outgoingX();
01285 float outY = outer->outgoingY();
01286 float in_r = inner->innerBoundHits().first()->layer().r();
01287 float out_r = outer->outerBoundHits().first()->layer().r();
01288 float allK = 2.*(-1. / 2.99792458 /m_pmgnIMF->getReferField())*(inY*outX-inX*outY) /
01289 (in_r*out_r*sqrt((inX-outX)*(inX-outX)+(inY-outY)*(inY-outY)));
01290
01291 float cache_in = ((outer->kappa() + allK)/3.0 - inner->kappa()*2/3.0)/in_r;
01292 float cache_out = ((inner->kappa() + allK)/3.0 - outer->kappa()*2/3.0)/out_r;
01293 float cache_all = ((inner->kappa()+outer->kappa())/3.0-allK*2/3.0)*2.0/(in_r+out_r);
01294 float chi2_z = cache_in*cache_in + cache_out*cache_out + cache_all*cache_all;
01295 log<<MSG::DEBUG<<"D_phi: "<< D_phi <<" chi2_z: "<< chi2_z <<" maxChi2: " <<maxChi2 <<endreq;
01296 if (chi2_z > maxChi2) continue;
01297 min_Dphi = D_phi;
01298 min_Dphi_index = j;
01299 }
01300 if (min_Dphi_index < 0) continue;
01301 log<<MSG::DEBUG<<"min_Dphi: "<< min_Dphi <<endreq;
01302 FTSegment * outer = OuterSegments[min_Dphi_index];
01303 const FTLayer & out_innerBound = outer->innerBoundHits().first()->layer();
01304 switch (out_innerBound.layerId() - in_outerBound.layerId()) {
01305 case 1:
01306 if (min_Dphi > maxDphi) continue;
01307 break;
01308 case 2:
01309 if (min_Dphi > maxDphi*1.5) continue;
01310 break;
01311 case 3:
01312 if (min_Dphi > maxDphi*2.25) continue;
01313 break;
01314 default:
01315 continue;
01316 }
01317 inner->connect_outer(outer->wireHits(), outer->outerBoundHits());
01318 inner->update();
01319 Connected.append(inner);
01320 n = InnerSegments.remove(i);
01321 delete outer;
01322 m = OuterSegments.remove(min_Dphi_index);
01323 log<<MSG::DEBUG<<"DONE!!"<<endreq;
01324 }
01325 return;
01326 }
01327
01328 void
01329 FTFinder::mkTrack3D(void){
01330 IMessageSvc *msgSvc;
01331 Gaudi::svcLocator()->service("MessageSvc", msgSvc);
01332 MsgStream log(msgSvc, "FTFinder");
01333
01334 int n = _tracks.length();
01335
01336 FTList<FTSegment *> multi_trk_cand(20);
01337 for (int i=8; i>=0 ; i-- ){
01338
01339 if (i == 4) i -= 3;
01340 FTList<FTSegment *> & segments = _superLayer[i].segments();
01341 int m = segments.length();
01342 for (int j = 0; j^m; j++){
01343 FTSegment * s = segments[j];
01344 #ifndef OnlineMode
01345 s->printout();
01346 #endif
01347 int nTrack = 0;
01348 FTTrack * t;
01349 for (int k = 0; k^n; k++){
01350 #ifndef OnlineMode
01351 log<<MSG::DEBUG<<"coupling to track: " << k <<endreq;
01352
01353 #endif
01354 if (s->update3D(_tracks[k])){
01355 t = _tracks[k];
01356 nTrack++;
01357 }
01358 }
01359 switch(nTrack){
01360 case 0:
01361 continue;
01362 case 1:
01363 t->append_stereo_cache(s);
01364 break;
01365 default:
01366 multi_trk_cand.append(s);
01367 break;
01368 }
01369 }
01370
01371 for (int j = 0; j^n; j++) _tracks[j]->updateSZ();
01372 }
01373
01374 for (int i = 0; i^n; i++) _tracks[i]->linkStereoSegments();
01375 n = multi_trk_cand.length();
01376 for (int i = 0; i^n; i++) multi_trk_cand[i]->linkStereoSegments();
01377 }
01378
01379 int
01380 FTFinder::VertexFit2D()
01381 {
01382 int nTrks = _tracks.length();
01383 if (nTrks < 2){
01384 _vx = -99999.;
01385 _vy = -99999.;
01386 _vz = -99999.;
01387 return 0;
01388 }
01389 FTList<float> px(10);
01390 FTList<float> py(10);
01391 FTList<float> dx(10);
01392 FTList<float> dy(10);
01393 FTList<double> weight(10);
01394 for (int i = 0; i < nTrks; i++){
01395 const Lpav & la = _tracks[i]->lpav();
01396 HepVector a = la.Hpar(pivot);
01397
01398 const float dr_i = a(1);
01399 if (fabs(a(1)) > 1.5) continue;
01400 const float px_i = - sin(a(2));
01401 const float py_i = cos(a(2));
01402
01403
01404
01405
01406 float weight_i = la.chisq()/(la.nc()*0.02);
01407
01408 px.append(px_i);
01409 py.append(py_i);
01410 dx.append(dr_i*py_i);
01411 dy.append(-dr_i*px_i);
01412 weight.append(exp(-weight_i*weight_i));
01413 }
01414 if (dx.length() < 2){
01415 _vx = -99999.;
01416 _vy = -99999.;
01417 _vz = -99999.;
01418 return 0;
01419 }
01420 double S_weight = 0.;
01421 for (int i = dx.length() - 1; i; i--){
01422 const float px_i = px[i];
01423 const float py_i = py[i];
01424 const float dx_i = dx[i];
01425 const float dy_i = dy[i];
01426 const float weight_i = weight[i];
01427 for (int j = 0; j < i; j++){
01428 const float px_j = px[j];
01429 const float py_j = py[j];
01430
01431
01432 const float ddx = dx[j] - dx_i;
01433 const float ddy = dx[j] - dy_i;
01434
01435 const float tmp_par = px_i*py_j - px_j*py_i;
01436 const float par = (py_j*ddx-px_j*ddy)/tmp_par;
01437
01438 double weight_ij = weight_i*weight[j];
01439 S_weight += weight_i*weight[j];
01440 if (tmp_par==0 ||
01441 par < -0.5 || (py_i*ddx-px_i*ddy)/tmp_par < -0.5 ||
01442 fabs((px_i*px_j + py_i*py_j)/
01443 sqrt((px_i*px_i+py_i*py_i)*(px_j*px_j+py_j*py_j)))>0.86){
01444 _vx += (dx_i+0.5*ddx)*weight_ij;
01445 _vy += (dy_i+0.5*ddy)*weight_ij;
01446 }else{
01447 _vx += (dx_i+par*px_i)*weight_ij;
01448 _vy += (dy_i+par*py_i)*weight_ij;
01449 }
01450 }
01451 }
01452
01453 int rtn_flag = 0;
01454 if (S_weight == 0.){
01455 _vx = -99999.;
01456 _vy = -99999.;
01457 _vz = -99999.;
01458 }else{
01459 _vx /= S_weight;
01460 _vy /= S_weight;
01461 _vz = -99999.;
01462 rtn_flag = 1;
01463 }
01464 return rtn_flag;
01465
01466 }
01467
01468 int
01469 FTFinder::VertexFit(int z_flag)
01470 {
01471 _vx = 0.;
01472 _vy = 0.;
01473 _vz = 0.;
01474 if (!z_flag){
01475 return VertexFit2D();
01476 }
01477 int n = _tracks.length();
01478 if (n < 2){
01479 _vx = -99999.;
01480 _vy = -99999.;
01481 _vz = -99999.;
01482 return 0;
01483 }
01484 FTList<float> px(10);
01485 FTList<float> py(10);
01486 FTList<float> pz(10);
01487 FTList<float> dx(10);
01488 FTList<float> dy(10);
01489 FTList<float> dz(10);
01490 FTList<float> pmag2(10);
01491 FTList<float> weight(10);
01492 FTList<float> weight_z(10);
01493 FTList<float> sigma2_r(10);
01494 FTList<float> sigma2_z(10);
01495
01496 for (int i = 0; i < n; i++){
01497 const Lpav & la = _tracks[i]->lpav();
01498 if (la.nc() <= 3) continue;
01499 const zav & za = _tracks[i]->Zav();
01500 if(za.nc() > 2 && za.b() > -900){
01501 pmag2.append(1+za.b()*za.b());
01502 pz.append(za.a());
01503 dz.append(za.b());
01504 sigma2_z.append(za.chisq()/(za.nc()-2));
01505 weight_z.append(exp(-(za.chisq()/(za.nc()-2))));
01506 }else{
01507 continue;
01508 }
01509
01510 HepVector a = la.Hpar(pivot);
01511
01512 const float dr_i = a(1);
01513 const float px_i = - std::sin(a(2));
01514 const float py_i = std::cos(a(2));
01515
01516 px.append(px_i);
01517 py.append(py_i);
01518 dx.append(dr_i*py_i);
01519 dy.append(-dr_i*px_i);
01520 sigma2_r.append(std::sqrt(la.chisq())/(la.nc()-3));
01521 weight.append(exp(-(std::sqrt(la.chisq())/(la.nc()-3))));
01522 }
01523
01524 n = dx.length();
01525 if (n < 2){
01526 _vx = -9999.;
01527 _vy = -9999.;
01528 _vz = -9999.;
01529 return 0;
01530 }
01531
01532 FTList<float> vtx_x(20);
01533 FTList<float> vtx_y(20);
01534 FTList<float> vtx_z(20);
01535 FTList<float> weight2(20);
01536 FTList<float> weight2_z(20);
01537 FTList<float> vtx_chi2(20);
01538 unsigned n_vtx(0);
01539 for (int i = n - 1; i; i--){
01540 for (int j = 0; j < i; j++){
01541
01542 const float pij_x = py[i]*pz[j]-pz[i]*py[j];
01543 const float pij_y = pz[i]*px[j]-px[i]*pz[j];
01544 const float pij_z = px[i]*py[j]-py[i]*px[j];
01545
01546 if (pij_x == 0.0f && pij_y == 0.0f && pij_z == 0.0f ) continue;
01547
01548 const float sr = sigma2_r[i]+sigma2_r[j];
01549 const float sz = sigma2_z[i]+sigma2_z[j];
01550
01551 const float ddx = dx[i]-dx[j];
01552 const float ddy = dy[i]-dy[j];
01553 const float ddz = dz[i]-dz[j];
01554
01555 const float pij_mag2 = pij_x*pij_x/(sr*sz)+pij_y*pij_y/(sr*sz)+pij_z*pij_z/(sr*sr);
01556
01557 const float d_i = (pij_x*(py[j]*ddz-pz[j]*ddy)/(sr*sz)+
01558 pij_y*(pz[j]*ddx-px[j]*ddz)/(sr*sz)+
01559 pij_z*(px[j]*ddy-py[j]*ddx)/(sr*sr))/pij_mag2;
01560
01561 const float d_j = (pij_x*(py[i]*ddz-pz[i]*ddy)/(sr*sz)+
01562 pij_y*(pz[i]*ddx-px[i]*ddz)/(sr*sz)+
01563 pij_z*(px[i]*ddy-py[i]*ddx)/(sr*sr))/pij_mag2;
01564
01565 const float vtx_x_i = dx[i] + px[i]*d_i;
01566 const float vtx_y_i = dy[i] + py[i]*d_i;
01567 const float vtx_z_i = dz[i] + pz[i]*d_i;
01568
01569 const float vtx_x_j = dx[j] + px[j]*d_j;
01570 const float vtx_y_j = dy[j] + py[j]*d_j;
01571 const float vtx_z_j = dz[j] + pz[j]*d_j;
01572
01573 const float weight_ij = weight[i]+weight[j];
01574 vtx_x.append((weight[j]*vtx_x_i + weight[i]*vtx_x_j)/weight_ij);
01575 vtx_y.append((weight[j]*vtx_y_i + weight[i]*vtx_y_j)/weight_ij);
01576 vtx_z.append((weight_z[j]*vtx_z_i + weight_z[i]*vtx_z_j)/(weight_z[i]+weight_z[j]));
01577
01578 weight2.append(exp(-sr));
01579 weight2_z.append(exp(-sz));
01580 vtx_chi2.append(((vtx_x_i-vtx_x_j)*(vtx_x_i-vtx_x_j)+(vtx_y_i-vtx_y_j)*(vtx_y_i-vtx_y_j))/sr +
01581 (vtx_z_i-vtx_z_j)*(vtx_z_i-vtx_z_j)/sz);
01582 n_vtx++;
01583 }
01584 }
01585 n = vtx_chi2.length();
01586 float S_weight(0.0f);
01587 float S_weight_z(0.0f);
01588 for (int i = 0; i != n; i++){
01589 if (vtx_chi2[i] > 10.) continue;
01590 float w(std::exp(-vtx_chi2[i]));
01591 _vx += vtx_x[i]*weight2[i]*w;
01592 _vy += vtx_y[i]*weight2[i]*w;
01593 _vz += vtx_z[i]*weight2_z[i]*w;
01594 S_weight += weight2[i]*w;
01595 S_weight_z += weight2_z[i]*w;
01596 }
01597 int rtn_flag = 0;
01598 if (S_weight <= 0.){
01599 _vx = -9999.;
01600 _vy = -9999.;
01601 }else{
01602 _vx /= S_weight;
01603 _vy /= S_weight;
01604 rtn_flag = 1;
01605 }
01606 if (!z_flag) return rtn_flag;
01607 if (S_weight_z <= 0.){
01608 _vz = -9999.;
01609 }else{
01610 _vz /= S_weight_z;
01611 rtn_flag++;
01612 }
01613 return rtn_flag;
01614 }
01615
01616 int
01617 FTFinder::findBestVertex(void)
01618 {
01619 int nTrks = _tracks.length();
01620 if (nTrks < 2){
01621 _vx = -9999.;
01622 _vy = -9999.;
01623 _vz = -9999.;
01624 return 0;
01625 }
01626 float min_dr = 9999.;
01627 float phi0 = 9999.;
01628 for (int i = 0; i < nTrks; i++){
01629 HepVector a = _tracks[i]->lpav().Hpar(pivot);
01630 if (fabs(a(1)) < fabs(min_dr)){
01631 min_dr = a(1);
01632 phi0 = a(2);
01633 }
01634 }
01635 _vx = min_dr*cos(phi0);
01636 _vy = min_dr*sin(phi0);
01637 return 1;
01638 }
01639
01640 Hep3Vector
01641 FTFinder::vertex(void) const{
01642 return Hep3Vector(_vx,_vy,_vz);
01643 }
01644
01645 int
01646 FTFinder::CorrectEvtTiming(void)
01647 {
01648 int nTrks = _tracks.length();
01649 float weight_sum = 0.;
01650 float dt_sum2 = 0.;
01651 for (int i = 0; i^nTrks; i++){
01652 float dt_sum = 0.;
01653 float dtt_sum = 0.;
01654 int nHits = 0;
01655 const Lpav & la = _tracks[i]->lpav();
01656 FTList<FTSegment *>& axial_sgmnts = _tracks[i]->axial_segments();
01657 int m = axial_sgmnts.length();
01658 for (int j = 0; j^m; j++){
01659 FTList<FTWire *>& hits = axial_sgmnts[j]->wireHits();
01660 int l = hits.length();
01661 for (int k = 0; k^l; k++){
01662 FTWire & h = * hits[k];
01663 const float x = h.x();
01664 const float y = h.y();
01665 double d0 = fabs(la.d((double)x,(double)y));
01666 if (d0 >= 0.47f*h.layer().csize()) continue;
01667 nHits++;
01668 float dt = x2t(h.layer(), d0) - h.time();
01669
01670 dt_sum += dt;
01671 dtt_sum += (dt*dt);
01672 }
01673 }
01674 if (!nHits) continue;
01675 float weight_t = exp(-(dtt_sum - (dt_sum*dt_sum/nHits))/(nHits*1600));
01676 weight_sum += (nHits*weight_t);
01677 dt_sum2 += (dt_sum*weight_t);
01678 }
01679
01680 return int(dt_sum2/weight_sum);
01681
01682 }
01683
01684 #ifndef OnlineMode
01685 void
01686 FTFinder::makeMdst(void)
01687 {
01688
01689 const int nTrks = _tracks.length();
01690 int Ntable(0);
01691
01692 for (int i = 0; i<nTrks; i++){
01693
01694 const FTTrack & trk = * _tracks[i];
01695 FTList<FTSegment *> &axialSegments = trk.axial_segments();
01696 FTList<FTSegment *> &stereoSegments = trk.stereo_segments();
01697 int axialSegmentsN = axialSegments.length();
01698 int stereoSegmentsN = stereoSegments.length();
01699 for (int i = 0; i < axialSegmentsN ; i++) {
01700
01701 FTList<FTWire *> &wires = axialSegments[i]->wireHits();
01702 int wiresN = wires.length();
01703 for (int j=0; j < wiresN; j++) {
01704
01705 g_ncell->fill(wires[j]->getWireId());
01706 }
01707 }
01708
01709 for (int i = 0; i < stereoSegmentsN ; i++) {
01710
01711 FTList<FTWire *> &wires = stereoSegments[i]->wireHits();
01712 int wiresN = wires.length();
01713 for (int j=0; j < wiresN; j++) {
01714
01715 g_ncell->fill(wires[j]->getWireId());
01716 }
01717 }
01718
01719 const Helix& h = * trk.helix();
01720 if (h.tanl() < -9000.) continue;
01721
01722 Ntable++;
01723 Hep3Vector p(h.momentum());
01724 HepPoint3D v(h.x());
01725
01726 float m_px = p.x();
01727 g_px[Ntable-1] = p.x();
01728 float m_py = p.y();
01729 g_py[Ntable-1] = p.y();
01730 float m_pz = p.z();
01731 g_pz[Ntable-1] = p.z();
01732 g_p[Ntable-1] = p.mag();
01733 g_phi[Ntable-1] = atan2(m_py, m_px);
01734
01735 g_dr[Ntable-1] = h.dr();
01736 g_phi0[Ntable-1] = h.phi0();
01737 g_kappa[Ntable-1] = h.kappa();
01738 g_pt[Ntable-1] = 1/fabs(h.kappa());
01739 g_dz[Ntable-1] = h.dz();
01740 g_tanl[Ntable-1] = h.tanl();
01741 g_theta[Ntable-1] = atan2(1/fabs(h.kappa()),(double)(m_pz));
01742 g_vx[Ntable-1] = v(0);
01743 g_vy[Ntable-1] = v(1);
01744 g_vz[Ntable-1] = v(2);
01745 }
01746 g_ntrk = Ntable;
01747
01748 }
01749 #endif
01750
01751
01752
01753
01754
01755
01756
01757 StatusCode
01758 FTFinder::makeTds(void)
01759 {
01760 IMessageSvc *msgSvc;
01761 Gaudi::svcLocator()->service("MessageSvc", msgSvc);
01762 MsgStream log(msgSvc, "FTFinder");
01763
01764 IDataProviderSvc* eventSvc = NULL;
01765 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
01766
01767 if (eventSvc) {
01768
01769 } else {
01770 log << MSG::FATAL << "makeTds:Could not find eventSvc" << endreq;
01771 return StatusCode::FAILURE ;
01772 }
01773
01774 RecMdcTrackCol* trkcol = new RecMdcTrackCol;
01775 RecMdcHitCol* hitcol = new RecMdcHitCol;
01776
01777
01778 #ifndef OnlineMode
01779 log << MSG::DEBUG << "beginning to make RecMdcTrackCol" <<endreq;
01780 #endif
01781 int ntrk = tracks().length();
01782 int trkid=0;
01783 for (int i =0; i < ntrk; i++) {
01784 RecMdcTrack* trk = new RecMdcTrack;
01785 FTTrack* fttrk = tracks()[i];
01786
01787 #ifndef OnlineMode
01788
01789 #endif
01790
01791 if (fttrk->helix()->tanl() < -9000.){
01792 log << MSG::DEBUG << "deleted trackId = " << i+1 << " due to tanl = " << fttrk->helix()->tanl() << endreq;
01793 delete trk;
01794 continue;
01795 }
01796
01797 HepPoint3D pos = fttrk->helix()->pivot();
01798 int charge = -1;
01799 HepVector m_a(5,0);
01800 m_a = fttrk->helix()->a();
01801 m_a[2]=-m_a[2];
01802 HepSymMatrix m_ea = fttrk->helix()->Ea();
01803 float fiterm = fttrk->lpav().phi(77.0);
01804 float chi2lpav = fttrk->lpav().chisq();
01805 float chi2zav = fttrk->Zav().chisq();
01806
01807
01808
01809
01810
01811 m_ea[0][0] = 0.0085;
01812 m_ea[1][1] = 0.000011;
01813 m_ea[2][2] = 0.0018;
01814 m_ea[3][3] = 1.2;
01815 m_ea[4][4] = 0.00026;
01816 m_ea[1][0] = m_ea[0][1] = -0.00029;
01817 m_ea[2][0] = m_ea[0][2] = charge*0.004;
01818 m_ea[2][1] = m_ea[1][2] = charge*0.00012;
01819 m_ea[3][0] = m_ea[0][3] = -0.017;
01820 m_ea[3][1] = m_ea[1][3] = 0.0;
01821 m_ea[3][2] = m_ea[2][3] = 0.0;
01822 m_ea[4][0] = m_ea[0][4] = 0.0;
01823 m_ea[4][1] = m_ea[1][4] = 0.0;
01824 m_ea[4][2] = m_ea[2][4] = 0.0;
01825 m_ea[4][3] = m_ea[3][4] = -0.018;
01826
01827 trk->setTrackId(trkid);
01828 trk->setHelix(m_a);
01829 trk->setPivot(pos);
01830 trk->setError(m_ea);
01831 trk->setFiTerm(fiterm);
01832 trk->setChi2(chi2lpav+chi2zav);
01833 #ifndef OnlineMode
01834 log <<MSG::DEBUG << " trackId = " << i << endreq;
01835 log <<MSG::DEBUG <<"fast-tracking kappa "<<m_a[2]
01836 <<" fast-tracking tanl "<<m_a[4]
01837 <<endreq;
01838 log <<MSG::DEBUG <<"push_backed kappa "<<trk->helix(2)
01839 <<" push_backed tanl "<< trk->helix(4)
01840 <<endreq;
01841
01842 log << MSG::DEBUG << "beginning to make hitlist and RecMdcHitCol " <<endreq;
01843 #endif
01844
01845 HitRefVec hitrefvec;
01846
01847 int hitindex = 0;
01848
01849
01850 FTList<FTSegment *>& seglist_ax = fttrk->axial_segments();
01851 int nseg_ax = seglist_ax.length();
01852 int ntrackhits = 0;
01853 for (int iseg_ax = 0; iseg_ax < nseg_ax; iseg_ax++) {
01854 FTList<FTWire *>& hitlist_ax = seglist_ax[iseg_ax]->wireHits();
01855 int nhitlist_ax = hitlist_ax.length();
01856 ntrackhits += nhitlist_ax;
01857 for( int ihit_ax = 0; ihit_ax < nhitlist_ax; ihit_ax++,hitindex++) {
01858 FTWire wire_ax = *(hitlist_ax[ihit_ax]);
01859 double dd = wire_ax.distance();
01860 double adc = RawDataUtil::MdcChargeChannel(wire_ax.getAdc());
01861 double tdc = wire_ax.time();
01862 double chi2 = wire_ax.getChi2();
01863 const Identifier mdcid = MdcID::wire_id(wire_ax.layer().layerId(),wire_ax.localId());
01864
01865 RecMdcHit* hit = new RecMdcHit;
01866 hit->setId(hitindex);
01867 hit->setDriftDistLeft( dd );
01868 hit->setDriftDistRight( dd );
01869 hit->setAdc( adc );
01870 hit->setTdc( tdc );
01871 hit->setMdcId( mdcid );
01872 hit->setChisqAdd( chi2 );
01873 #ifndef OnlineMode
01874 log << MSG::DEBUG << "Hit DriftDistLeft " << hit->getDriftDistLeft() << endreq;
01875 log << MSG::DEBUG << "MDC Hit ADC " << hit->getAdc() << endreq;
01876 log << MSG::DEBUG << "Hit MDC Identifier " << hit->getMdcId() << endreq;
01877 log << MSG::DEBUG << "Hit Chisq " << hit->getChisqAdd() << endreq;
01878 #endif
01879 hitcol->push_back(hit);
01880 SmartRef<RecMdcHit> refhit(hit);
01881 hitrefvec.push_back(refhit);
01882 }
01883 }
01884
01885
01886 FTList<FTSegment *>& seglist_st = fttrk->stereo_segments();
01887 int nseg_st = seglist_st.length();
01888 int ntrackster = 0;
01889
01890 for (int iseg_st = 0; iseg_st < nseg_st; iseg_st++) {
01891 FTList<FTWire *>& hitlist_st = seglist_st[iseg_st]->wireHits();
01892 int nhitlist_st = hitlist_st.length();
01893 ntrackhits += nhitlist_st;
01894 ntrackster += nhitlist_st;
01895 for( int ihit_st = 0; ihit_st < nhitlist_st; ihit_st++,hitindex++) {
01896 FTWire wire_st = *(hitlist_st[ihit_st]);
01897 double dd = wire_st.distance();
01898
01899 double adc = RawDataUtil::MdcChargeChannel(wire_st.getAdc());
01900 double tdc = wire_st.time();
01901 double chi2 = wire_st.getChi2();
01902 const Identifier mdcid = MdcID::wire_id(wire_st.layer().layerId(),wire_st.localId());
01903
01904 RecMdcHit* hit = new RecMdcHit;
01905 hit->setId(hitindex);
01906 hit->setDriftDistLeft( dd );
01907 hit->setDriftDistRight( dd );
01908 hit->setAdc( adc );
01909 hit->setTdc( tdc);
01910 hit->setMdcId( mdcid );
01911 hit->setChisqAdd( chi2 );
01912 #ifndef OnlineMode
01913 log << MSG::DEBUG << "Z Hit DriftDistLeft " << hit->getDriftDistLeft() << endreq;
01914 log << MSG::DEBUG << "Z MDC Hit ADC " << hit->getAdc() << endreq;
01915 log << MSG::DEBUG << "Z Hit MDC Identifier " << hit->getMdcId() << endreq;
01916 log << MSG::DEBUG << "Z Hit Chisq " << hit->getChisqAdd() << endreq;
01917 #endif
01918 hitcol->push_back(hit);
01919 SmartRef<RecMdcHit> refhit(hit);
01920 hitrefvec.push_back(refhit);
01921 }
01922 }
01923 trk->setNhits(ntrackhits);
01924 trk->setNster(ntrackster);
01925 trk->setVecHits(hitrefvec);
01926 trkcol->push_back(trk);
01927 trkid++;
01928 #ifndef OnlineMode
01929 g_naxialhit->fill((float)(ntrackhits-ntrackster), 1.0);
01930 g_nstereohit->fill((float)ntrackster, 1.0);
01931 g_nhit->fill((float)ntrackhits, 1.0);
01932 #endif
01933 }
01934
01935 IDataManagerSvc *dataManSvc = dynamic_cast<IDataManagerSvc*> (eventSvc);
01936 DataObject *aTrackCol;
01937 eventSvc->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
01938 if(aTrackCol != NULL) {
01939 dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
01940 eventSvc->unregisterObject("/Event/Recon/RecMdcTrackCol");
01941 }
01942 DataObject *aRecHitCol;
01943 eventSvc->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
01944 if(aRecHitCol != NULL) {
01945 dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
01946 eventSvc->unregisterObject("/Event/Recon/RecMdcHitCol");
01947 }
01948
01949
01950 StatusCode hitsc;
01951 hitsc = eventSvc->registerObject("/Event/Recon/RecMdcHitCol", hitcol);
01952 if( hitsc.isFailure() ) {
01953 log << MSG::FATAL << "Could not register RecMdcHitCol" << endreq;
01954 return hitsc;
01955 }
01956 #ifndef OnlineMode
01957 log << MSG::DEBUG << "RecMdcHitCol registered successfully!" <<endreq;
01958
01959 RecMdcTrackCol::iterator bef = trkcol->begin();
01960 for( ; bef != trkcol->end(); bef++) {
01961 log <<MSG::DEBUG <<" registered kappa"<<(*bef)->helix(2)
01962 <<" registered tanl"<< (*bef)->helix(4)
01963 <<endreq;
01964 }
01965 #endif
01966
01967
01968 StatusCode trksc;
01969 trksc = eventSvc->registerObject("/Event/Recon/RecMdcTrackCol", trkcol);
01970 if( trksc.isFailure() ) {
01971 log << MSG::FATAL << "Could not register RecMdcHit" << endreq;
01972 return trksc;
01973 }
01974 #ifndef OnlineMode
01975 log << MSG::DEBUG << "RecMdcTrackCol registered successfully!" <<endreq;
01976
01977
01978 SmartDataPtr<RecMdcHitCol> newhitCol(eventSvc,"/Event/Recon/RecMdcHitCol");
01979 if (!newhitCol) {
01980 log << MSG::FATAL << "Could not find RecMdcHit" << endreq;
01981 return( StatusCode::FAILURE);
01982 }
01983 log << MSG::DEBUG << "Begin to check RecMdcHitCol"<<endreq;
01984 RecMdcHitCol::iterator iter_hit = newhitCol->begin();
01985 for( ; iter_hit != newhitCol->end(); iter_hit++){
01986 log << MSG::DEBUG << "retrieved MDC Hit:"
01987 << " DDL " <<(*iter_hit)->getDriftDistLeft()
01988 << " DDR " <<(*iter_hit)->getDriftDistRight()
01989 << " ADC " <<(*iter_hit)->getAdc() << endreq;
01990 }
01991
01992
01993 SmartDataPtr<RecMdcTrackCol> newtrkCol(eventSvc,"/Event/Recon/RecMdcTrackCol");
01994 if (!newtrkCol) {
01995 log << MSG::FATAL << "Could not find RecMdcTrackCol" << endreq;
01996 return( StatusCode::FAILURE);
01997 }
01998 log << MSG::DEBUG << "Begin to check RecMdcTrackCol"<<endreq;
01999 RecMdcTrackCol::iterator iter_trk = newtrkCol->begin();
02000 for( ; iter_trk != newtrkCol->end(); iter_trk++){
02001 log << MSG::DEBUG << "retrieved MDC track:"
02002 << "Track Id: " << (*iter_trk)->trackId()
02003 << " Pivot: " << (*iter_trk)->poca()[0] << " "
02004 << (*iter_trk)->poca()[1] << " " << (*iter_trk)->poca()[2]
02005 << endreq ;
02006
02007
02008
02009
02010
02011
02012
02013
02014
02015
02016
02017
02018
02019
02020 log << MSG::DEBUG << "hitList of this track:" << endreq;
02021 HitRefVec gothits = (*iter_trk)->getVecHits();
02022 HitRefVec::iterator it_gothit = gothits.begin();
02023 for( ; it_gothit != gothits.end(); it_gothit++){
02024 log << MSG::DEBUG << "hits Id: "<<(*it_gothit)->getId()
02025 << " hits DDL&DDR " <<(*it_gothit)->getDriftDistLeft()
02026 << " hits MDC IDentifier " <<(*it_gothit)->getMdcId()
02027 << endreq
02028 << " hits TDC " <<(*it_gothit)->getTdc()
02029 << " hits ADC " <<(*it_gothit)->getAdc() << endreq;
02030 }
02031
02032 }
02033 #endif
02034
02035 return StatusCode::SUCCESS;
02036 }
02037