00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014 #include <math.h>
00015 #include "CLHEP/String/Strings.h"
00016
00017 #include "MdcGeomSvc/MdcGeomSvc.h"
00018 #include "MdcTables/MdcTables.h"
00019
00020 #include "TrkReco/TrkReco.h"
00021 #include "TrkReco/TMDC.h"
00022 #include "TrkReco/TMDCWire.h"
00023 #include "TrkReco/TMDCWireHit.h"
00024 #include "TrkReco/TPerfectFinder.h"
00025 #include "TrkReco/TConformalFinder0.h"
00026 #include "TrkReco/TConformalFinder.h"
00027 #include "TrkReco/TCurlFinder.h"
00028 #include "TrkReco/TFastFinder.h"
00029 #include "TrkReco/TTrack.h"
00030 #include "TrkReco/TTrackMC.h"
00031 #include "TrkReco/TTrackHEP.h"
00032 #include "TrkReco/TMLink.h"
00033 #include "TrkReco/TTrackBase.h"
00034
00035 #include "EvTimeEvent/RecEsTime.h"
00036 #include "MdcCalibFunSvc/IMdcCalibFunSvc.h"
00037 #include "MdcCalibFunSvc/MdcCalibFunSvc.h"
00038
00039 #include "GaudiKernel/MsgStream.h"
00040 #include "GaudiKernel/AlgFactory.h"
00041 #include "GaudiKernel/ISvcLocator.h"
00042 #include "GaudiKernel/IDataManagerSvc.h"
00043 #include "GaudiKernel/SmartDataPtr.h"
00044 #include "GaudiKernel/IDataProviderSvc.h"
00045 #include "GaudiKernel/PropertyMgr.h"
00046 #include "GaudiKernel/IMessageSvc.h"
00047 #include "GaudiKernel/IDataProviderSvc.h"
00048 #include "GaudiKernel/IJobOptionsSvc.h"
00049 #include "GaudiKernel/Bootstrap.h"
00050 #include "GaudiKernel/StatusCode.h"
00051
00052 #include "Identifier/Identifier.h"
00053 #include "Identifier/MdcID.h"
00054 #include "ReconEvent/ReconEvent.h"
00055 #include "EventModel/EventHeader.h"
00056 #include "MdcRawEvent/MdcDigi.h"
00057 #include "McTruth/McParticle.h"
00058 #include "McTruth/MdcMcHit.h"
00059 #include "RawEvent/RawDataUtil.h"
00060 #include "RawDataProviderSvc/IRawDataProviderSvc.h"
00061
00062 #include "RawDataProviderSvc/MdcRawDataProvider.h"
00063
00064
00065 #include <vector>
00066 #include <iostream>
00067 #include <fstream>
00068
00069 #include "BesTimerSvc/IBesTimerSvc.h"
00070 #include "BesTimerSvc/BesTimerSvc.h"
00071 #include "BesTimerSvc/BesTimer.h"
00072
00073
00074 using namespace std;
00075 using namespace Event;
00076
00077 int TrkRecoTest = 0;
00078 float TrkRecoHelixFitterChisqMax = 0;
00079
00080 TrkReco::TrkReco(const std::string& name, ISvcLocator* pSvcLocator) :
00081 Algorithm(name, pSvcLocator),
00082 _cdc(0),
00083 _perfectFinder(0),
00084 _rkfitter("range fitter",
00085 false,
00086 0,
00087 true),
00088 _confFinder(0),
00089 _curlFinder(0),
00090 _nEvents(0) {
00091
00093 initPara();
00094
00095 TrkRecoHelixFitterChisqMax = b_helixFitterChisqMax;
00096 }
00097
00098
00099 StatusCode TrkReco::initialize(){
00100 MsgStream log(msgSvc(), name());
00101 log << MSG::INFO << "in initialize()" << endreq;
00102
00103 StatusCode sc;
00104
00106 IRawDataProviderSvc* irawDataProviderSvc;
00107 sc = service ("RawDataProviderSvc", irawDataProviderSvc);
00108 _rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*> (irawDataProviderSvc);
00109 if ( sc.isFailure() ){
00110 log << MSG::FATAL << "Could not load RawDataProviderSvc!" << endreq;
00111 return StatusCode::FAILURE;
00112 }
00113
00115
00116
00117
00118
00119
00120
00121
00123 IMdcCalibFunSvc* imdcCalibSvc;
00124 sc = service ("MdcCalibFunSvc",imdcCalibSvc);
00125 _mdcCalibFunSvc = dynamic_cast<MdcCalibFunSvc*> (imdcCalibSvc);
00126 if ( sc.isFailure() ){
00127 log << MSG::FATAL << "Could not load MdcCalibFunSvc!" << endreq;
00128 }
00129
00131 if (b_tuple) InitTuple();
00132
00133 #ifdef TRKRECO_DEBUG
00134 b_debugLevel = 2;
00135 #endif
00136
00137
00138
00139
00140
00141 TrkRecoHelixFitterChisqMax = b_helixFitterChisqMax;
00142
00143
00144
00145
00146
00147
00148
00149 if ((! _confFinder) && (b_conformalFinder == 0)) {
00150 _confFinder = new TConformalFinder0(b_conformalMaxSigma,
00151 b_conformalFraction,
00152 b_conformalStereoZ3,
00153 b_conformalStereoZ4,
00154 b_conformalStereoChisq3,
00155 b_conformalStereoChisq4,
00156 b_conformalStereoMaxSigma,
00157 b_conformalFittingCorrections,
00158 b_salvageLevel,
00159 b_doConformalFinderCosmic);
00160 }
00161 else if (! _confFinder) {
00162 _confFinder = new TConformalFinder(b_doConformalFastFinder,
00163 b_doConformalSlowFinder,
00164 b_conformalPerfectSegmentFinding,
00165 b_conformalMaxSigma,
00166 b_conformalStereoMaxSigma,
00167 b_salvageLevel,
00168 b_conformalMinNLinksForSegment,
00169 b_conformalMinNCores,
00170 b_conformalMinNSegments,
00171 b_conformalSalvageLoadWidth,
00172 b_conformalStereoMode,
00173 b_conformalStereoLoadWidth,
00174 b_conformalStereoSzSegmentDistance,
00175 b_conformalStereoSzLinkDistance,
00176 b_conformalFittingFlag);
00177 }
00178 _confFinder->debugLevel(b_debugLevel);
00179 _confFinder->doStereo(b_doConformalFinderStereo);
00180
00181
00182 if (b_doSalvage == 2) _confFinder->doSalvage(true);
00183
00184
00185 if (! _curlFinder)
00186 _curlFinder = new TCurlFinder((unsigned)min_segment,
00187 (unsigned)min_salvage,
00188 bad_distance_for_salvage,
00189 good_distance_for_salvage,
00190 (unsigned)min_sequence,
00191 (unsigned)min_fullwire,
00192 range_for_axial_search,
00193 range_for_stereo_search,
00194 (unsigned)superlayer_for_stereo_search,
00195 range_for_axial_last2d_search,
00196 range_for_stereo_last2d_search,
00197 trace2d_distance,
00198 trace2d_first_distance,
00199 trace3d_distance,
00200 (unsigned)determine_one_track,
00201 selector_max_impact,
00202 selector_max_sigma,
00203 selector_strange_pz,
00204 selector_replace_dz,
00205 (unsigned)stereo_2dfind,
00206 (unsigned)merge_exe,
00207 merge_ratio,
00208 merge_z_diff,
00209 mask_distance,
00210 ratio_used_wire,
00211 range_for_stereo1,
00212 range_for_stereo2,
00213 range_for_stereo3,
00214 range_for_stereo4,
00215 range_for_stereo5,
00216 range_for_stereo6,
00217 z_cut,
00218 z_diff_for_last_attend,
00219 (unsigned)svd_reconstruction,
00220 min_svd_electrons,
00221 (unsigned)on_correction,
00222 (unsigned)output_2dtracks,
00223 (unsigned)curl_version,
00224
00225 minimum_seedLength,
00226 minimum_2DTrackLength,
00227 minimum_3DTrackLength,
00228 minimum_closeHitsLength,
00229 MIN_RADIUS_OF_STRANGE_TRACK,
00230 ULTIMATE_MIN_RADIUS_OF_STRANGE_TRACK);
00231 _curlFinder->debugLevel(b_debugLevel);
00232
00233
00234 _trackManager.maxMomentum(b_momentumCut);
00235 _trackManager.debugLevel(b_debugLevel);
00236 _trackManager.fittingFlag(b_fittingFlag);
00237
00238
00239
00240
00241
00242
00243
00244 TrkRecoTest = b_test;
00245
00246 if (b_timeTest && b_tuple) {
00247 StatusCode sc = service( "BesTimerSvc", m_timersvc);
00248 if( sc.isFailure() ) {
00249 log << MSG::WARNING << " Unable to locate BesTimerSvc" << endreq;
00250 return StatusCode::FAILURE;
00251 }
00252 m_timer[1] = m_timersvc->addItem("Execution");
00253 m_timer[1]->propName("nExecution");
00254 }
00255
00256 return StatusCode::SUCCESS;
00257 }
00258
00259 TMDC *
00260 TrkReco::cdcInit(void) {
00261 _cdcVersion = dtostring(b_cdcVersion);
00262 if (! _cdc) _cdc = TMDC::getTMDC(_cdcVersion);
00263 _cdc->debugLevel(b_debugLevel);
00264
00265
00266 float fudgeFactor = b_fudgeFactor;
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277 fudgeFactor = 1.0;
00278
00279 _cdc->fudgeFactor(fudgeFactor);
00280 return _cdc;
00281 }
00282
00283
00284 StatusCode TrkReco::finalize() {
00285
00286 if (b_timeTest && b_tuple) m_timersvc->print();
00287
00288
00289 if (_cdc) _cdc->fastClear();
00290 if (b_doConformalFinder && _confFinder) _confFinder->clear();
00291 if (b_doCurlFinder && _curlFinder) _curlFinder->clear();
00292
00293
00294
00295
00296
00297 if (_confFinder) delete _confFinder;
00298 if (_curlFinder) delete _curlFinder;
00299
00300 return StatusCode::SUCCESS;
00301 }
00302 StatusCode TrkReco::beginRun(){
00304 IMdcGeomSvc* imdcGeomSvc;
00305 StatusCode sc1 =Gaudi::svcLocator()->service("MdcGeomSvc", imdcGeomSvc);
00306 _mdcGeomSvc = dynamic_cast<MdcGeomSvc*> (imdcGeomSvc);
00307 if (sc1.isFailure()) {
00308 return( StatusCode::FAILURE);
00309 }
00310 _cdc = cdcInit();
00311
00312 return StatusCode::SUCCESS;
00313 }
00314 void
00315 TrkReco::disp_stat(const char * m) {
00316 std::string msg = m;
00317
00318 }
00319
00320
00321 StatusCode TrkReco::execute() {
00322 MsgStream log(msgSvc(), name());
00323 log << MSG::INFO << "in execute()" << endreq;
00324 StatusCode sc;
00325 if (b_timeTest && b_tuple) m_timer[1]->start();
00326
00327
00328 if (b_goodTrk && b_tuple)
00329 for (int ii=0;ii<43;ii++){
00330 for (int jj=0;jj<288;jj++){
00331 havedigi[ii][jj]= -99;
00332 }
00333 }
00334
00335
00336
00337
00338
00339
00340 if(!b_RungeKuttaCorrection){
00341 IDataManagerSvc *dataManSvc;
00342 DataObject *aTrackCol;
00343 DataObject *aRecHitCol;
00344 if(!m_combineTracking){
00345 dataManSvc = dynamic_cast<IDataManagerSvc*> (eventSvc().get());
00346 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
00347 if(aTrackCol != NULL) {
00348 dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
00349 eventSvc()->unregisterObject("/Event/Recon/RecMdcTrackCol");
00350 }
00351 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
00352 if(aRecHitCol != NULL) {
00353 dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
00354 eventSvc()->unregisterObject("/Event/Recon/RecMdcHitCol");
00355 }
00356 }
00357 DataObject *aReconEvent;
00358 eventSvc()->findObject("/Event/Recon",aReconEvent);
00359 if(!aReconEvent) {
00360 ReconEvent* recevt = new ReconEvent;
00361 StatusCode sc = eventSvc()->registerObject("/Event/Recon",recevt );
00362 if(sc!=StatusCode::SUCCESS) {
00363 log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
00364 return( StatusCode::FAILURE);
00365 }
00366 }
00367 RecMdcTrackCol* trkcol;
00368 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
00369 if (aTrackCol) {
00370 trkcol = dynamic_cast<RecMdcTrackCol*> (aTrackCol);
00371 }else{
00372 trkcol= new RecMdcTrackCol;
00373 sc = eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trkcol);
00374 if(!sc.isSuccess()) {
00375 log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
00376 return StatusCode::FAILURE;
00377 }
00378 }
00379 RecMdcHitCol* hitcol;
00380 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
00381 if (aRecHitCol) {
00382 hitcol = dynamic_cast<RecMdcHitCol*> (aRecHitCol);
00383 }else{
00384 hitcol= new RecMdcHitCol;
00385 sc = eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitcol);
00386 if(!sc.isSuccess()) {
00387 log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
00388 return StatusCode::FAILURE;
00389 }
00390 }
00391
00392 t0_bes = 0.;
00393 t0Sta = -1;
00394 if (useESTime) {
00395 SmartDataPtr<RecEsTimeCol> aevtimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
00396
00397 if (aevtimeCol && aevtimeCol->size() ) {
00398 RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
00399 t0_bes = (*iter_evt)->getTest();
00400 t0Sta = (*iter_evt)->getStat();
00401 }else{
00402 log << MSG::WARNING << "Could not find EsTimeCol" << endreq;
00403 return StatusCode::SUCCESS;
00404 }
00405 }
00406
00407
00408 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00409 if (!eventHeader) {
00410 log << MSG::FATAL << "Could not find Event Header" << endreq;
00411 return( StatusCode::FAILURE);
00412 }
00413 _nEvents = eventHeader->eventNumber();
00414 if (b_tuple) std::cout << "TrkReco ... processing ev# " << _nEvents << std::endl;
00415
00416
00417 int digiId;
00418 uint32_t getDigiFlag = 0;
00419 if(m_dropHot || m_combineTracking)getDigiFlag |= MdcRawDataProvider::b_dropHot;
00420 if(m_keepBadTdc) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
00421 if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
00422 MdcDigiVec mdcDigiVec = _rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
00423 if (0 == mdcDigiVec.size()){
00424 log << MSG::WARNING << " No hits in MdcDigiVec" << endreq;
00425 return StatusCode::SUCCESS;
00426 }
00427
00428
00429
00430
00431
00432 MdcDigiVec::iterator iter1 = mdcDigiVec.begin();
00433 MC_DIGI_SIZE = mdcDigiVec.size();
00434
00435 digiId = 0;
00436 Identifier mdcId;
00437 int layerId;
00438 int wireId;
00439
00440
00441 unsigned nt = MdcRecWirhitCol::getMdcRecWirhitCol()->size();
00442
00443 MdcRecWirhitCol::getMdcRecWirhitCol()->clear();
00444
00445 for (;iter1 != mdcDigiVec.end(); iter1++, digiId++) {
00446
00447 mdcId = (*iter1)->identify();
00448 layerId = MdcID::layer(mdcId);
00449 wireId = MdcID::wire(mdcId);
00450
00451
00452
00453
00454
00455
00456
00457
00458 if (b_goodTrk && b_tuple) havedigi[layerId][wireId] = (*iter1)->getTrackIndex();
00459
00460 MdcRec_wirhit mhit;
00461
00462 mhit.id = digiId;
00463 mhit.geo = _mdcGeomSvc->Wire(layerId,wireId);
00464
00465 double tof;
00466 tof = 1.28 * mhit.geo->Lyr()->Radius()/299.8;
00467 mhit.tdc = RawDataUtil::MdcTime((*iter1)->getTimeChannel());
00468
00469
00470 mhit.tdc -= t0_bes;
00471
00472
00473 _trackManager.sett0bes(t0_bes);
00474
00475
00476 mhit.adc = (*iter1)->getChargeChannel();
00477
00478
00479
00480 double timewalk=_mdcCalibFunSvc->getTimeWalk(layerId,mhit.adc);
00481 double T0 = _mdcCalibFunSvc->getT0(layerId,wireId);
00482 double drifttime = mhit.tdc-tof-timewalk-T0;
00483 double dist2 = _mdcCalibFunSvc->driftTimeToDist(drifttime,layerId, wireId, 2, 0);
00484 mhit.erddl = _mdcCalibFunSvc->getSigma(layerId, 2 , dist2, 0.0,0.0,0.0,mhit.adc);
00485
00486 mhit.ddl = dist2/10.;
00487 mhit.ddr = mhit.ddl;
00488 mhit.erddl = mhit.erddl/10.;
00489 mhit.erddr = mhit.erddl;
00490
00491 mhit.lr = 2;
00492 mhit.stat = 0;
00493 mhit.stat = mhit.stat |= 1048576;
00494 mhit.stat = mhit.stat |= 2097152;
00495 mhit.stat = mhit.stat |= 4194304;
00496 mhit.stat = mhit.stat |= 1073741824;
00497
00498
00499
00500
00501 MdcRecWirhitCol::getMdcRecWirhitCol()->push_back(mhit);
00502
00503 if (b_tuple) {
00504 t10_tdc = mhit.tdc;
00505 t10_adc = mhit.adc;
00506 t10_drift = mhit.ddl;
00507 t10_dDrift = mhit.erddl;
00508 t10_lyrId = layerId;
00509 t10_localId = wireId;
00510 m_tuple10->write();
00511 }
00512 }
00513
00514 unsigned nT0Reset = 0;
00515
00516
00517 TrkReco_start:
00518
00519
00520 clear();
00521
00522
00523 if (mdcDigiVec.size() > 2000){
00524 log << MSG::WARNING << " Too many hits in MdcDigiVec" << endreq;
00525 return StatusCode::SUCCESS;
00526 }
00527
00528
00529 _cdc->update(b_doMCAnalysis);
00530
00531
00532 unsigned mask = 0;
00533 if (! b_useAllHits) mask = WireHitFindingValid;
00534 const AList<TMDCWireHit> & axialHits = _cdc->axialHits(mask);
00535 const AList<TMDCWireHit> & stereoHits = _cdc->stereoHits(mask);
00536 const AList<TMDCWireHit> & allHits = _cdc->hits(mask);
00537
00538
00539
00540 AList<TTrack> tracks;
00541 AList<TTrack> tracks2D;
00542
00543
00544 if (b_doPerfectFinder) {
00545 _perfectFinder->doit(axialHits, stereoHits, tracks, tracks2D);
00546 _trackManager.append(tracks);
00547 }
00548
00549 else {
00550
00551 if (b_doConformalFinder) {
00552
00553
00554 if (b_doT0Reset) {
00555 if (b_nT0ResetMax > nT0Reset)
00556 ((TConformalFinder *) _confFinder)->doT0Reset(true);
00557 else
00558 ((TConformalFinder *) _confFinder)->doT0Reset(false);
00559 }
00560
00561 _confFinder->doit(axialHits, stereoHits, tracks, tracks2D);
00562
00563
00564 if (b_doT0Reset) {
00565 ++nT0Reset;
00566 if (((TConformalFinder *) _confFinder)->T0ResetDone())
00567 goto TrkReco_start;
00568 }
00569
00570
00571
00572
00573 _trackManager.append(tracks);
00574 _trackManager.append2D(tracks2D);
00575 if (b_conformalFinder == 0) {
00576 if (b_doSalvage == 1) _trackManager.salvage(allHits);
00577 if (b_doSalvage) _trackManager.mask();
00578 }
00579 }
00580
00581
00582 if (b_doCurlFinder) {
00583 if ((! b_doSalvage) && (b_conformalFinder == 0))
00584 _trackManager.maskCurlHits(axialHits,
00585 stereoHits,
00586 _trackManager.tracks());
00587 AList<TTrack> confTracks = _trackManager.tracks();
00588 tracks.append(confTracks);
00589 _curlFinder->doit(axialHits, stereoHits, tracks, tracks2D);
00590 tracks.remove(confTracks);
00591
00592 }
00593
00594
00595
00596
00597
00598 if (b_doCurlFinder) {
00599 _trackManager.append(tracks);
00600 _trackManager.append2D(tracks2D);
00601 }
00602
00603
00604
00605 _trackManager.merge();
00606
00607
00608
00609
00610
00611 }
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627 if (b_tuple) FillTuple();
00628
00629
00630 int t_tkStat = 2;
00631 if (!b_doConformalFinder && b_doCurlFinder) t_tkStat = 3;
00632 StatusCode scTds = _trackManager.makeTds(trkcol, hitcol, t_tkStat,b_RungeKuttaCorrection,m_CalibFlag);
00633 if (scTds != StatusCode::SUCCESS) return( StatusCode::FAILURE);
00634
00635
00636 if (b_doT0Determination) {
00637 _trackManager.determineT0(b_doT0Determination, b_nTracksForT0);
00638 if (b_tuple) t3_t0Rec = _trackManager.paraT0();
00639 }
00640 }
00641 else if(b_RungeKuttaCorrection==1){
00642 if (useESTime) {
00643 SmartDataPtr<RecEsTimeCol> aevtimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
00644 if (aevtimeCol && aevtimeCol->size() ) {
00645 RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
00646 t0_bes = (*iter_evt)->getTest();
00647 t0Sta = (*iter_evt)->getStat();
00648 }else{
00649 log << MSG::WARNING << "Could not find EsTimeCol" << endreq;
00650 return StatusCode::SUCCESS;
00651 }
00652 }
00653
00654 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00655 if (!eventHeader) {
00656 log << MSG::FATAL << "Could not find Event Header" << endreq;
00657 return( StatusCode::FAILURE);
00658 }
00659 _nEvents = eventHeader->eventNumber();
00660 if (b_tuple) std::cout << "TrkReco ... processing ev# " << _nEvents << std::endl;
00661 AList<TTrack> rktracks;
00662 int digiId;
00663 uint32_t getDigiFlag = 0;
00664 if(m_dropHot || m_combineTracking)getDigiFlag |= MdcRawDataProvider::b_dropHot;
00665 if(m_keepBadTdc) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
00666 if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
00667 MdcDigiVec mdcDigiVec = _rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
00668 if (0 == mdcDigiVec.size()){
00669 log << MSG::WARNING << " No hits in MdcDigiVec" << endreq;
00670 return StatusCode::SUCCESS;
00671 }
00672 MdcDigiVec::iterator iter1 = mdcDigiVec.begin();
00673 MC_DIGI_SIZE = mdcDigiVec.size();
00674
00675 digiId = 0;
00676 Identifier mdcId;
00677 int layerId;
00678 int wireId;
00679
00680 unsigned nt = MdcRecWirhitCol::getMdcRecWirhitCol()->size();
00681 MdcRecWirhitCol::getMdcRecWirhitCol()->clear();
00682
00683 for (;iter1 != mdcDigiVec.end(); iter1++, digiId++) {
00684 mdcId = (*iter1)->identify();
00685 layerId = MdcID::layer(mdcId);
00686 wireId = MdcID::wire(mdcId);
00687 if (b_goodTrk && b_tuple) havedigi[layerId][wireId] = (*iter1)->getTrackIndex();
00688
00689 MdcRec_wirhit mhit;
00690
00691 mhit.id = digiId;
00692 mhit.geo = _mdcGeomSvc->Wire(layerId,wireId);
00693 double tof;
00694 tof = 1.28 * mhit.geo->Lyr()->Radius()/299.8;
00695 mhit.tdc = RawDataUtil::MdcTime((*iter1)->getTimeChannel());
00696 mhit.timechannel=(*iter1)->getTimeChannel();
00697 mhit.tdc -= t0_bes;
00698
00699 _trackManager.sett0bes(t0_bes);
00700 mhit.adc = (*iter1)->getChargeChannel();
00701 double timewalk=_mdcCalibFunSvc->getTimeWalk(layerId,mhit.adc);
00702 double T0 = _mdcCalibFunSvc->getT0(layerId,wireId);
00703 double drifttime = mhit.tdc-tof-timewalk-T0;
00704 double dist2 = _mdcCalibFunSvc->driftTimeToDist(drifttime,layerId, wireId, 2, 0);
00705 mhit.erddl = _mdcCalibFunSvc->getSigma(layerId, 2 , dist2, 0.0,0.0,0.0,mhit.adc);
00706 mhit.ddl = dist2/10.;
00707 mhit.ddr = mhit.ddl;
00708 mhit.erddl = mhit.erddl/10.;
00709 mhit.erddr = mhit.erddl;
00710
00711 mhit.lr = 2;
00712 mhit.stat = 0;
00713 mhit.stat = mhit.stat |= 1048576;
00714 mhit.stat = mhit.stat |= 2097152;
00715 mhit.stat = mhit.stat |= 4194304;
00716 mhit.stat = mhit.stat |= 1073741824;
00717 MdcRecWirhitCol::getMdcRecWirhitCol()->push_back(mhit);
00718
00719 if (b_tuple) {
00720 t10_tdc = mhit.tdc;
00721 t10_adc = mhit.adc;
00722 t10_drift = mhit.ddl;
00723 t10_dDrift = mhit.erddl;
00724 t10_lyrId = layerId;
00725 t10_localId = wireId;
00726 m_tuple10->write();
00727 }
00728 }
00729
00730 unsigned nT0Reset = 0;
00731
00732 clear();
00733 if (mdcDigiVec.size() > 2000){
00734 log << MSG::WARNING << " Too many hits in MdcDigiVec" << endreq;
00735 return StatusCode::SUCCESS;
00736 }
00737 int counter=0;
00738 _cdc->update(b_doMCAnalysis);
00739
00740 unsigned mask = 0;
00741 if (! b_useAllHits) mask = WireHitFindingValid;
00742 const AList<TMDCWireHit> & axialHits = _cdc->axialHits(mask);
00743 const AList<TMDCWireHit> & stereoHits = _cdc->stereoHits(mask);
00744 const AList<TMDCWireHit> & allHits = _cdc->hits(mask);
00745 AList<TMLink> _allHits[3];
00746 TConformalFinder0::conformalTransformationDriftCircle(ORIGIN, axialHits, _allHits[0]);
00747 TConformalFinder0::conformalTransformationDriftCircle(ORIGIN, stereoHits, _allHits[1]);
00748 _allHits[2].append(_allHits[0]);
00749 _allHits[2].append(_allHits[1]);
00750
00751
00752 SmartDataPtr<RecMdcTrackCol> mdcTracks(eventSvc(),EventModel::Recon::RecMdcTrackCol);
00753 if( ! mdcTracks )
00754 {
00755 log << MSG::ERROR << "Unable to retrieve RecMdcTrackCol" << endreq;
00756 return StatusCode::FAILURE;
00757 } else {
00758 log << MSG::DEBUG << "RecMdcTrackCol retrieved of size "<< mdcTracks->size() << endreq;
00759 int ntrk=0;
00760 for(RecMdcTrackCol::iterator it=mdcTracks->begin(); it!=mdcTracks->end(); it++)
00761 {ntrk++;
00762 TRunge *r = new TRunge(**it);
00763 TTrack *t1 = new TTrack(*r);
00764 HitRefVec gothits = (*it)->getVecHits();
00765 HitRefVec::iterator it_gothit = gothits.begin();
00766 unsigned stat=(*it)->stat();
00767 int nhit=0;
00768 for( ; it_gothit != gothits.end(); it_gothit++){
00769 for(int i=0;i<_allHits[2].length();i++){
00770 int lyrraw=_allHits[2][i]->wire()->layerId();
00771 int wireraw=_allHits[2][i]->wire()->localId();
00772 int g_layer = MdcID::layer((*it_gothit)->getMdcId());
00773 int g_wire = MdcID::wire((*it_gothit)->getMdcId());
00774 if(lyrraw==g_layer&&g_wire==wireraw){
00775 nhit++;
00776 t1->append(*_allHits[2][i]);
00777 }
00778 }
00779 }
00780 TRunge *rr = new TRunge(*t1);
00781 TRunge *tt = new TRunge(*t1);
00782 int err=_rkfitter.fit(*rr);
00783 int nhits=rr->links().length();
00784 int ndrop=0;
00785 int nmax=0;
00786 if(nhits<=10)nmax=0;
00787 if(nhits>10)nmax=(int)nhit*0.3;
00788 for(int ii=0;ii<nmax;ii++){
00789
00790 ndrop=rr->DropWorst();
00791 if(ndrop)err=_rkfitter.fit(*rr);
00792 if(err)break;
00793 }
00794 if(err==-2) counter++;
00795 if(m_CalibFlag==1) {
00796 tt->removeLinks();
00797 tt->append(rr->links());
00798 for(int i=0;i<43;i++){
00799 err= _rkfitter.fit(*tt,0,i);
00800 }
00801 rr->removeLinks();
00802 rr->append(tt->links());
00803 TTrack *t =new TTrack(*rr);
00804 if(err==0){
00805 rktracks.append(*t);
00806 t->setFinderType(100+stat);
00807 }
00808 }
00809 if(m_CalibFlag==0){
00810 TTrack *t =new TTrack(*rr);
00811 if(err==0){
00812
00813 t->setFinderType(100+stat);
00814 rktracks.append(*t);
00815 }
00816 }
00817 delete r ;
00818 delete t1;
00819 delete rr;
00820 delete tt;
00821 }
00822 }
00823
00824 _trackManager.append(rktracks);
00825 IDataManagerSvc *dataManSvc;
00826 DataObject *aTrackCol;
00827 DataObject *aRecHitCol;
00828 if(!m_combineTracking){
00829 dataManSvc = dynamic_cast<IDataManagerSvc*> (eventSvc().get());
00830 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
00831 if(aTrackCol != NULL) {
00832 dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
00833 eventSvc()->unregisterObject("/Event/Recon/RecMdcTrackCol");
00834 }
00835 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
00836 if(aRecHitCol != NULL) {
00837 dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
00838 eventSvc()->unregisterObject("/Event/Recon/RecMdcHitCol");
00839 }
00840 }
00841 DataObject *aReconEvent;
00842 eventSvc()->findObject("/Event/Recon",aReconEvent);
00843 if(!aReconEvent) {
00844 ReconEvent* recevt = new ReconEvent;
00845 StatusCode sc = eventSvc()->registerObject("/Event/Recon",recevt );
00846 if(sc!=StatusCode::SUCCESS) {
00847 log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
00848 return( StatusCode::FAILURE);
00849 }
00850 }
00851 RecMdcTrackCol* trkcol;
00852 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
00853 if (aTrackCol) {
00854 trkcol = dynamic_cast<RecMdcTrackCol*> (aTrackCol);
00855 }else{
00856 trkcol= new RecMdcTrackCol;
00857 sc = eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trkcol);
00858 if(!sc.isSuccess()) {
00859 log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
00860 return StatusCode::FAILURE;
00861 }
00862 }
00863 RecMdcHitCol* hitcol;
00864 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
00865 if (aRecHitCol) {
00866 hitcol = dynamic_cast<RecMdcHitCol*> (aRecHitCol);
00867 }else{
00868 hitcol= new RecMdcHitCol;
00869 sc = eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitcol);
00870 if(!sc.isSuccess()) {
00871 log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
00872 return StatusCode::FAILURE;
00873 }
00874 }
00875 if (b_tuple) FillTuple();
00876
00877 int t_tkStat = 2;
00878 if (!b_doConformalFinder && b_doCurlFinder) t_tkStat = 3;
00879 StatusCode scTds = _trackManager.makeTds(trkcol, hitcol, t_tkStat,b_RungeKuttaCorrection,m_CalibFlag);
00880 if (b_doT0Determination) {
00881 _trackManager.determineT0(b_doT0Determination, b_nTracksForT0);
00882 if (b_tuple) t3_t0Rec = _trackManager.paraT0();
00883 }
00884 for (unsigned i = 0; i < 3; i++) {
00885 if (i == 2)
00886 HepAListDeleteAll(_allHits[i]);
00887 else
00888 _allHits[i].removeAll();
00889 }
00890
00891 rktracks.removeAll();
00892 }
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905 return StatusCode::SUCCESS;
00906 }
00907
00908 void
00909 TrkReco::mcInformation(void) {
00910
00911
00912
00913 const AList<TTrack> & allTracks = _trackManager.tracksFinal();
00914
00915 unsigned nHep = TTrackHEP::list().length();
00916 unsigned nTrk = allTracks.length();
00917 unsigned * N = (unsigned *) malloc(nHep * sizeof(unsigned));
00918 for (unsigned i = 0; i < nHep; i++) N[i] = 0;
00919
00920
00921 for (unsigned i = 0; i < nTrk; i++) {
00922 TTrackMC * mc = allTracks[i]->_mc;
00923 if (! mc) {
00924 mc = new TTrackMC(* allTracks[i]);
00925 _mcTracks.append(mc);
00926 allTracks[i]->_mc = mc;
00927 }
00928
00929 mc->update();
00930 if (mc->hepId() != -1)
00931 if (mc->charge())
00932 ++N[mc->hepId()];
00933 }
00934
00935
00936 for (unsigned i = 0; i < nHep; i++) {
00937 if (N[i] < 2) {
00938 for (unsigned j = 0; j < nTrk; j++)
00939 if (allTracks[j]->_mc->hepId() == i)
00940 allTracks[j]->_mc->_quality += TTrackUnique;
00941 }
00942 }
00943
00944
00945 for (unsigned i = 0; i < nTrk; i++) {
00946 unsigned & quality = allTracks[i]->_mc->_quality;
00947 if ((quality & TTrackGhost) && (quality & TTrackUnique))
00948 quality += TTrackGood;
00949 }
00950
00951
00952 free(N);
00953 }
00954
00955 void
00956 TrkReco::clear(void) {
00957
00958
00959 HepAListDeleteAll(_mcTracks);
00960
00961
00962 if (b_doPerfectFinder) _perfectFinder->clear();
00963 if (b_doConformalFinder) _confFinder->clear();
00964 if (b_doCurlFinder) _curlFinder->clear();
00965 _trackManager.clear();
00966 _trackManager.clearTables();
00967 }
00968
00969 void
00970 TrkReco::fastClear(void) {
00971 clear();
00972 }
00973
00974 bool
00975 TrkReco::mcEvent(void) const {
00976
00977
00978
00979
00980
00981
00982
00983 return true;
00984 }
00985
00986
00987
00988
00989
00990
00991 void
00992 TrkReco::initPara(void){
00993 _rkfitter.setBesFromGdml();
00994 declareProperty("mcPar", b_mcPar = 0);
00995 declareProperty("mcHit", b_mcHit = 0);
00996 declareProperty("tuple", b_tuple = 0);
00997 declareProperty("goodTrk", b_goodTrk = 0);
00998 declareProperty("timeTest", b_timeTest = 0);
00999
01000 declareProperty("cdcVersion", b_cdcVersion = 1.0);
01001 declareProperty("fudgeFactor", b_fudgeFactor = 1.0);
01002
01003 declareProperty("useESTime", useESTime = 1.0);
01004
01005 declareProperty("debugLevel", b_debugLevel = 0);
01006 declareProperty("useAllHits", b_useAllHits = 0);
01007 declareProperty("doT0Reset", b_doT0Reset = 0);
01008 declareProperty("nT0ResetMax", b_nT0ResetMax = 1);
01009 declareProperty("doMCAnalysis", b_doMCAnalysis = 1);
01010 declareProperty("helixFitterChisqMax", b_helixFitterChisqMax = 0.);
01011
01012 declareProperty("RungeKuttaCorrection", b_RungeKuttaCorrection = 0);
01013 declareProperty("doPerfectFinder", b_doPerfectFinder = 0);
01014 declareProperty("perfectFitting", b_perfectFitting = 0);
01015
01016 declareProperty("conformalFinder", b_conformalFinder = 1);
01017 declareProperty("doConformalFinder", b_doConformalFinder = 0);
01018 declareProperty("doConformalFastFinder", b_doConformalFastFinder = 1);
01019 declareProperty("doConformalSlowFinder", b_doConformalSlowFinder = 1);
01020 declareProperty("conformalPerfectSegmentFinding", b_conformalPerfectSegmentFinding = 0);
01021 declareProperty("conformalFittingFlag", b_conformalFittingFlag = 4);
01022 declareProperty("conformalMaxSigma", b_conformalMaxSigma = 30.);
01023 declareProperty("conformalMinNLinksForSegment", b_conformalMinNLinksForSegment = 2);
01024 declareProperty("conformalMinNCores", b_conformalMinNCores = 2);
01025 declareProperty("conformalMinNSegments", b_conformalMinNSegments = 2);
01026 declareProperty("salvageLevel", b_salvageLevel = 30.);
01027 declareProperty("conformalSalvageLoadWidth", b_conformalSalvageLoadWidth = 1);
01028 declareProperty("conformalStereoMode", b_conformalStereoMode = 1);
01029 declareProperty("conformalStereoLoadWidth", b_conformalStereoLoadWidth = 5);
01030 declareProperty("conformalStereoMaxSigma", b_conformalStereoMaxSigma = 30.);
01031 declareProperty("conformalStereoSzSegmentDistance", b_conformalStereoSzSegmentDistance = 10.);
01032 declareProperty("conformalStereoSzLinkDistance", b_conformalStereoSzLinkDistance = 15.);
01033
01034 declareProperty("doConformalFinderStereo", b_doConformalFinderStereo = 1);
01035 declareProperty("doConformalFinderCosmic", b_doConformalFinderCosmic = 0);
01036 declareProperty("conformalFraction", b_conformalFraction = 0.7);
01037 declareProperty("conformalStereoZ3", b_conformalStereoZ3 = 20.);
01038 declareProperty("conformalStereoZ4", b_conformalStereoZ4 = 20.);
01039 declareProperty("conformalStereoChisq3", b_conformalStereoChisq3 = 30.);
01040 declareProperty("conformalStereoChisq4", b_conformalStereoChisq4 = 30.);
01041 declareProperty("conformalFittingCorrections", b_conformalFittingCorrections = 0);
01042
01043 declareProperty("doCurlFinder", b_doCurlFinder = 1);
01044 declareProperty("doSalvage", b_doSalvage = 2);
01045 declareProperty("doMerge", b_doMerge = 0);
01046 declareProperty("momentumCut", b_momentumCut = 0.0);
01047 declareProperty("doT0Determination", b_doT0Determination = 7);
01048 declareProperty("nTracksForT0", b_nTracksForT0 = 2);
01049 declareProperty("sortMode", b_sortMode = 0);
01050 declareProperty("test", b_test = 0);
01051 declareProperty("fittingFlag", b_fittingFlag = 0);
01052 declareProperty("doAssociation", b_doAssociation = 1);
01053 declareProperty("associateSigma", b_associateSigma = 60);
01054
01055 declareProperty("dropHot", m_dropHot= 0);
01056 declareProperty("combineTracking",m_combineTracking=0);
01057 declareProperty("keepBadTdc",m_keepBadTdc=0);
01058 declareProperty("keepUnmatch",m_keepUnmatch=0);
01059
01060 declareProperty("CalibFlag",m_CalibFlag=0);
01061
01062 declareProperty("min_segment", min_segment = 5);
01063 declareProperty("min_salvage", min_salvage = 10);
01064 declareProperty("bad_distance_for_salvage", bad_distance_for_salvage = 1.0);
01065 declareProperty("good_distance_for_salvage", good_distance_for_salvage = 0.2);
01066 declareProperty("min_sequence", min_sequence = 6);
01067 declareProperty("min_fullwire", min_fullwire = 5);
01068 declareProperty("range_for_axial_search", range_for_axial_search = 1.5);
01069 declareProperty("range_for_stereo_search", range_for_stereo_search = 2.5);
01070 declareProperty("superlayer_for_stereo_search", superlayer_for_stereo_search = 3);
01071 declareProperty("range_for_axial_last2d_search", range_for_axial_last2d_search = 1.5);
01072 declareProperty("range_for_stereo_last2d_search", range_for_stereo_last2d_search = 2.0);
01073 declareProperty("trace2d_distance", trace2d_distance = 35.0);
01074 declareProperty("trace2d_first_distance", trace2d_first_distance = 0.5);
01075 declareProperty("trace3d_distance", trace3d_distance = 30.0);
01076 declareProperty("determine_one_track", determine_one_track = 0);
01077 declareProperty("selector_max_impact", selector_max_impact = 4.0);
01078 declareProperty("selector_max_sigma", selector_max_sigma = 36.0);
01079 declareProperty("selector_strange_pz", selector_strange_pz = 10.0);
01080 declareProperty("selector_replace_dz", selector_replace_dz = 15.0);
01081 declareProperty("stereo_2dfind", stereo_2dfind = 0);
01082 declareProperty("merge_exe", merge_exe = 1);
01083 declareProperty("merge_ratio", merge_ratio = 0.1);
01084 declareProperty("merge_z_diff", merge_z_diff = 10.0);
01085 declareProperty("mask_distance", mask_distance = 0.5);
01086 declareProperty("ratio_used_wire", ratio_used_wire = 0.3);
01087 declareProperty("range_for_stereo1", range_for_stereo1 = 2.4);
01088 declareProperty("range_for_stereo2", range_for_stereo2 = 2.7);
01089 declareProperty("range_for_stereo3", range_for_stereo3 = 2.9);
01090 declareProperty("range_for_stereo4", range_for_stereo4 = 3.4);
01091 declareProperty("range_for_stereo5", range_for_stereo5 = 4.1);
01092 declareProperty("range_for_stereo6", range_for_stereo6 = 5.0);
01093 declareProperty("z_cut", z_cut = 50.0);
01094 declareProperty("z_diff_for_last_attend", z_diff_for_last_attend = 1.5);
01095 declareProperty("on_correction", on_correction = 4);
01096 declareProperty("output_2dtracks", output_2dtracks = 1);
01097 declareProperty("curl_version", curl_version = 0);
01098
01099 declareProperty("minimum_seedLength", minimum_seedLength = 5);
01100 declareProperty("minimum_2DTrackLength", minimum_2DTrackLength = 7);
01101 declareProperty("minimum_3DTrackLength", minimum_3DTrackLength = 10);
01102 declareProperty("minimum_closeHitsLength", minimum_closeHitsLength = 5);
01103 declareProperty("MIN_RADIUS_OF_STRANGE_TRACK",MIN_RADIUS_OF_STRANGE_TRACK = 50);
01104 declareProperty("ULTIMATE_MIN_RADIUS_OF_STRANGE_TRACK",ULTIMATE_MIN_RADIUS_OF_STRANGE_TRACK = 5);
01105
01106 }
01107
01108 void
01109 TrkReco::InitTuple(void){
01110 m_tuple=ntupleSvc()->book("FILE101/rec3D",CLID_ColumnWiseTuple,"MdcRecEvent");
01111 m_tuple->addItem ("mc_phi", t_mcphi);
01112 m_tuple->addItem ("mc_theta", t_mctheta);
01113 m_tuple->addItem ("mc_ptot", t_mcptot);
01114 m_tuple->addItem ("mc_pt", t_mcpt);
01115 m_tuple->addItem ("mc_pz", t_mcpz);
01116 m_tuple->addItem ("mc_t0", t_mct0);
01117 m_tuple->addItem ("nDigi", t_nDigi);
01118
01119 m_tuple->addItem ("dr", t_dr);
01120 m_tuple->addItem ("dz", t_dz);
01121 m_tuple->addItem ("pt", t_pt);
01122 m_tuple->addItem ("ptot", t_ptot);
01123 m_tuple->addItem ("tanlmd",t_tanlmd);
01124 m_tuple->addItem ("phi", t_phi);
01125 m_tuple->addItem ("radius",t_radius);
01126 m_tuple->addItem ("chi2", t_chi2);
01127 m_tuple->addItem ("nHits", t_nHits);
01128 m_tuple->addItem ("nCores",t_nCores);
01129 m_tuple->addItem ("nSegs", t_nSegs);
01130 m_tuple->addItem ("ndf", t_ndf);
01131
01132 m_tuple->addItem ("dpt", t_dpt);
01133 m_tuple->addItem ("dptot", t_dptot);
01134 m_tuple->addItem ("dlmd", t_dlmd);
01135 m_tuple->addItem ("dphi", t_dphi);
01136
01137 m_tuple->addItem ("t0", t_t0);
01138 m_tuple->addItem ("t0_sta",t0_sta);
01139
01140 m_tuple->addItem ("evtNo", t_evtNo);
01141 m_tuple->addItem ("length", t_length);
01142 m_tuple->addItem ("length2", t_length2);
01143
01144 m_tuple->addItem ("gd_theta", t_good_theta);
01145 m_tuple->addItem ("gd_nLayers", t_gdNLayers);
01146 m_tuple->addItem ("mc_nLayers", t_mcNLayers);
01147 m_tuple->addItem ("best_nLayers",t_bestNLayers);
01148 m_tuple->addItem ("best_mc_nLayers",t_bestMcNLayers);
01149
01150 m_tuple3=ntupleSvc()->book("FILE101/raw",CLID_ColumnWiseTuple,"Raw Data");
01151 m_tuple3->addItem ("mc_t0", t3_mct0);
01152 m_tuple3->addItem ("mc_theta", t3_mctheta);
01153 m_tuple3->addItem ("mc_phi", t3_mcphi);
01154 m_tuple3->addItem ("mc_ptot", t3_mcptot);
01155 m_tuple3->addItem ("mc_pt", t3_mcpt);
01156 m_tuple3->addItem ("mc_pid", t3_mcpid);
01157 m_tuple3->addItem ("evtNo", t3_evtNo);
01158
01159 m_tuple31=ntupleSvc()->book("FILE101/rawEvt",CLID_ColumnWiseTuple,"Raw Data");
01160 m_tuple31->addItem ("nDigi", t3_nDigi);
01161 m_tuple31->addItem ("goodLength", t3_goodLength);
01162 m_tuple31->addItem ("length", t3_length);
01163 m_tuple31->addItem ("t0_rec", t3_t0Rec);
01164 m_tuple31->addItem ("t0", t3_t0);
01165 m_tuple31->addItem ("t0_sta", t3_t0Sta);
01166 m_tuple31->addItem ("finalLength", t3_finalLength);
01167
01168 m_tuple31->addItem ("mc_theta", t3_mctheta);
01169 m_tuple31->addItem ("mc_ptot", t3_mcptot);
01170 m_tuple31->addItem ("mc_pt", t3_mcpt);
01171 m_tuple31->addItem ("evtNo", t3_evtNo);
01172
01173 m_tuple2=ntupleSvc()->book("FILE101/rec2D",CLID_ColumnWiseTuple,"MdcRecEvent 2D");
01174 m_tuple2->addItem ("mc_theta", t2_mctheta);
01175 m_tuple2->addItem ("mc_pt", t2_mcpt);
01176 m_tuple2->addItem ("nDigi", t2_nDigi);
01177 m_tuple2->addItem ("nHits", t2_nHits);
01178 m_tuple2->addItem ("nSegs", t2_nSegs);
01179 m_tuple2->addItem ("chi2", t2_chi2);
01180 m_tuple2->addItem ("evtNo", t2_evtNo);
01181 m_tuple2->addItem ("ndf", t2_ndf);
01182 m_tuple2->addItem ("length", t2_length);
01183 m_tuple2->addItem ("length2", t2_length2);
01184 m_tuple2->addItem ("radius", t2_radius);
01185
01186 m_tuple4=ntupleSvc()->book("FILE101/hit",CLID_ColumnWiseTuple,"MdcRecEvent Hits");
01187 m_tuple4->addItem ("d_cal", t4_Dist);
01188 m_tuple4->addItem ("d_meas", t4_drift);
01189 m_tuple4->addItem ("e_meas", t4_dDrift);
01190 m_tuple4->addItem ("mc_drift", t4_mcDrift);
01191 m_tuple4->addItem ("mcLR", t4_mcLR);
01192 m_tuple4->addItem ("pull", t4_pull);
01193 m_tuple4->addItem ("lyrId", t4_lyrId);
01194 m_tuple4->addItem ("localId", t4_localId);
01195 m_tuple4->addItem ("LR", t4_LR);
01196 m_tuple4->addItem ("tdc", t4_tdc);
01197 m_tuple4->addItem ("z", t4_z);
01198 m_tuple4->addItem ("bz", t4_bz);
01199 m_tuple4->addItem ("fz", t4_fz);
01200 m_tuple4->addItem ("fy", t4_fy);
01201 m_tuple4->addItem ("phi", t4_phi);
01202 m_tuple4->addItem ("nHits", t4_nHits);
01203
01204 m_tuple5=ntupleSvc()->book("FILE101/char",CLID_ColumnWiseTuple,"MdcRecEvent Charge");
01205 m_tuple5->addItem ("ptotPos", t5_ptotPos);
01206 m_tuple5->addItem ("ptotNeg", t5_ptotNeg);
01207 m_tuple5->addItem ("drPos", t5_drPos);
01208 m_tuple5->addItem ("drNeg", t5_drNeg);
01209 m_tuple5->addItem ("dzPos", t5_dzPos);
01210 m_tuple5->addItem ("dzNeg", t5_dzNeg);
01211
01212 m_tuple6=ntupleSvc()->book("FILE101/urec",CLID_ColumnWiseTuple,"MdcRecEvent urec");
01213 m_tuple6->addItem ("length2", u_length2);
01214 m_tuple6->addItem ("mc_ptot", u_mcptot);
01215 m_tuple6->addItem ("mc_pt", u_mcpt);
01216 m_tuple6->addItem ("mc_theta", u_mctheta);
01217 m_tuple6->addItem ("nDigi", u_nDigi);
01218 m_tuple6->addItem ("evtNo", u_evtNo);
01219 m_tuple6->addItem ("mc_t0", u_mct0);
01220 m_tuple6->addItem ("t0", ut_t0);
01221 m_tuple6->addItem ("t0_sta", ut0_sta);
01222
01223 if (b_timeTest && b_tuple) {
01224 m_tuple7=ntupleSvc()->book("FILE101/time",CLID_ColumnWiseTuple,"MdcRecEvent time");
01225 m_tuple7->addItem ("time", ti_eventTime);
01226 m_tuple7->addItem ("recNum", ti_recTrkNum);
01227 m_tuple7->addItem ("evtNo", ti_evtNo);
01228 m_tuple7->addItem ("nHits", ti_nHits);
01229 m_tuple7->addItem ("nDigi", ti_nDigi);
01230 }
01231
01232 m_tuple9=ntupleSvc()->book("FILE101/seg",CLID_ColumnWiseTuple,"MdcRecEvent segments");
01233 m_tuple9->addItem ("times", t9_times);
01234 m_tuple9->addItem ("nLinks", t9_nLinks);
01235 m_tuple9->addItem ("nUsed", t9_nUsed);
01236 m_tuple9->addItem ("nSL", t9_nSL);
01237 m_tuple9->addItem ("mctheta", t9_mctheta);
01238
01239 m_tuple10=ntupleSvc()->book("FILE101/rawHit",CLID_ColumnWiseTuple,"MdcRecEvent mchitCOL");
01240 m_tuple10->addItem ("tdc", t10_tdc);
01241 m_tuple10->addItem ("adc", t10_adc);
01242 m_tuple10->addItem ("Drift", t10_drift);
01243 m_tuple10->addItem ("dDrift", t10_dDrift);
01244 m_tuple10->addItem ("lyrId", t10_lyrId);
01245 m_tuple10->addItem ("localId", t10_localId);
01246 }
01247
01248 void
01249 TrkReco::FillTuple(void) {
01250 MsgStream log(msgSvc(), name());
01251 log << MSG::INFO << "fill Tuple()" << endreq;
01252
01254 AList<TTrack> tracks;
01255 AList<TTrack> tracks2D;
01256 tracks = _trackManager.tracks();
01257 tracks2D = _trackManager.tracks2D();
01258
01259
01260 if(tracks.length()==0) cout<<"zslength: 3D length=0, and the 2D length is"<<tracks2D.length()<<endl;
01261
01263 if (b_timeTest && b_tuple) {
01264 m_timer[1]->stop();
01265 ti_eventTime = m_timer[1]->elapsed();
01266 ti_nDigi = MC_DIGI_SIZE;
01267 ti_recTrkNum = tracks.length();
01268 ti_evtNo = _nEvents;
01269 for (unsigned i = 0; i < ti_recTrkNum; ++i)
01270 ti_nHits += tracks[i]->nLinks();
01271 m_tuple7->write();
01272 }
01273
01275 CLHEP::Hep3Vector MC_TRACK_VEC;
01276 CLHEP::HepLorentzVector MC_TRACK_LRZVEC;
01277 float MC_TRACK_NUM;
01278 double MC_EVENT_TIME;
01279
01280 int digiId;
01281 int num_par = 0;
01282 MC_EVENT_TIME = -1;
01283 if (b_mcPar) {
01284 SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
01285 if (mcParticleCol) {
01286 Event::McParticleCol::iterator iter_mc = mcParticleCol->begin();
01287 digiId = 0;
01288 for (;iter_mc != mcParticleCol->end(); iter_mc++, digiId++) {
01289 if((*iter_mc)->statusFlags()==8320||(*iter_mc)->statusFlags()==128) {
01290 MC_EVENT_TIME = (*iter_mc)->initialFourMomentum().t();
01291 break;
01292 }
01293 }
01294
01295 iter_mc = mcParticleCol->begin();
01296 MC_TRACK_LRZVEC = (*iter_mc)->initialFourMomentum();
01297 MC_TRACK_VEC = (*iter_mc)->initialFourMomentum().vect();
01298 MC_TRACK_NUM = mcParticleCol->size();
01299
01300 #ifdef TRKRECO_DEBUG
01301 iter_mc = mcParticleCol->begin();
01302 digiId = 0;
01303 for (;iter_mc != mcParticleCol->end(); iter_mc++, digiId++) {
01304 log << MSG::INFO << "MDC digit No: " << digiId << endreq;
01305 log << MSG::INFO
01306 << " particleId = " << (*iter_mc)->particleProperty()
01307 << " theta = " << (*iter_mc)->initialFourMomentum().theta()
01308 <<" phi= "<< (*iter_mc)->initialFourMomentum().phi()
01309 <<" px= "<< (*iter_mc)->initialFourMomentum().px()
01310 <<" py= "<< (*iter_mc)->initialFourMomentum().py()
01311 <<" pz= "<< (*iter_mc)->initialFourMomentum().pz()
01312 << endreq;
01313 }
01314 #endif
01315 digiId = 0;
01316 iter_mc = mcParticleCol->begin();
01317 for (;iter_mc != mcParticleCol->end(); iter_mc++, digiId++) {
01318 CLHEP::HepLorentzVector LRZVEC = (*iter_mc)->initialFourMomentum();
01319 CLHEP::Hep3Vector VEC = (*iter_mc)->initialFourMomentum().vect();
01320
01321
01322 t3_mctheta = LRZVEC.theta();
01323 t3_mcphi = VEC.phi();
01324 t3_mcptot = VEC.mag()/1000.;
01325 t3_mcpt = VEC.perp()/1000.;
01326 t3_mct0 = (*iter_mc)->initialFourMomentum().t();
01327 t3_mcpid = (*iter_mc)->particleProperty();
01328 t3_evtNo = _nEvents;
01329 if((t3_mcpid==211||t3_mcpid==-211||t3_mcpid==321)&&fabs(cos(t3_mctheta))<0.93) ++num_par;
01330 m_tuple3->write();
01331 }
01332 }
01333 }
01334
01336 if (tracks.length()==0) {
01337 u_length2 = tracks2D.length();
01338 u_mct0 = MC_EVENT_TIME;
01339 u_mcptot = MC_TRACK_VEC.mag()/1000.;
01340 u_mcpt = MC_TRACK_VEC.perp()/1000.;
01341 u_mctheta = MC_TRACK_LRZVEC.theta();
01342 u_nDigi = MC_DIGI_SIZE;
01343 u_evtNo = _nEvents;
01344 ut0_sta = -1;
01345 ut_t0 = t0_bes;
01346 ut0_sta = t0Sta;
01347 m_tuple6->write();
01348 }
01349
01350
01351 float dDot = 2;
01352 int flagi = -1, flagj = -1;
01353 HepVector3D ppp1, ppp2;
01354 float pppdDot = 2;
01355 for (unsigned i = 0; i < tracks.length(); i++){
01356 if(tracks.length()<2) break;
01357 if(tracks[i]->charge()>0) {
01358 ppp1 = tracks[i]->p().unit();
01359 for (unsigned j = 0; j < tracks.length(); j++){
01360 if (tracks[j]->charge()<0) {
01361 ppp2 = tracks[j]->p().unit();
01362 pppdDot = ppp1.dot(ppp2);
01363 if (pppdDot<dDot) {
01364 flagi = i;
01365 flagj = j;
01366 dDot = pppdDot;
01367 }
01368 }
01369 }
01370 }
01371 }
01372 if (flagi != -1 && flagj != -1 && dDot < 0) {
01373 t5_ptotPos = tracks[flagi]->ptot();
01374 t5_ptotNeg = tracks[flagj]->ptot();
01375 t5_drPos = tracks[flagi]->helix().dr();
01376 t5_drNeg = tracks[flagj]->helix().dr();
01377 t5_dzPos = tracks[flagi]->helix().dz();
01378 t5_dzNeg = tracks[flagj]->helix().dz();
01379 m_tuple5->write();
01380 }
01381
01382 unsigned nGood = 0;
01383 for (unsigned i = 0; i < tracks.length(); i++){
01384 for(unsigned k = 0; k < tracks[i]->links().length(); k++){
01385 t4_Dist = tracks[i]->links()[k]->distance();
01386 t4_drift = tracks[i]->links()[k]->drift();
01387 t4_dDrift= tracks[i]->links()[k]->dDrift();
01388 t4_pull = tracks[i]->links()[k]->pull();
01389 t4_LR = 2;
01390 t4_LR = tracks[i]->links()[k]->leftRight();
01391 if (t4_LR == 0) t4_LR = -1;
01392 t4_lyrId = tracks[i]->links()[k]->wire()->layerId();
01393 t4_localId = tracks[i]->links()[k]->wire()->localId();
01394 t4_tdc = tracks[i]->links()[k]->hit()->reccdc()->tdc;
01395 t4_z = tracks[i]->links()[k]->positionOnTrack().z();
01396 t4_bz = tracks[i]->links()[k]->wire()->backwardPosition().z();
01397 t4_fz = tracks[i]->links()[k]->wire()->forwardPosition().z();
01398 t4_fy = tracks[i]->links()[k]->wire()->forwardPosition().y();
01399 t4_nHits = tracks[i]->links().length();
01400 unsigned lyrID = tracks[i]->links()[k]->wire()->layerId();
01401 float phi0 = _cdc->layer(lyrID)->offset();
01402 int nWir = (int) _cdc->layer(lyrID)->nWires();
01403 t4_phi = phi0 + t4_localId*2*pi/nWir;
01404 if(t4_phi<0) t4_phi+=2*pi;
01405
01406 if (b_mcHit) {
01407 SmartDataPtr<Event::MdcMcHitCol> mcMdcMcHitCol(eventSvc(),"/Event/MC/MdcMcHitCol");
01408 if (mcMdcMcHitCol) {
01409 Event::MdcMcHitCol::iterator iter_mchit1 = mcMdcMcHitCol->begin();
01410 digiId=0;
01411 for (;iter_mchit1 != mcMdcMcHitCol->end(); iter_mchit1++, digiId++) {
01412 const Identifier ident = (*iter_mchit1)->identify();
01413 if(MdcID::layer(ident) != t4_lyrId) continue;
01414 if(MdcID::wire(ident)==t4_localId) {
01415 t4_mcDrift = (*iter_mchit1)->getDriftDistance();
01416 t4_mcLR = (*iter_mchit1)->getPositionFlag();
01417 if (t4_mcLR == 0) t4_mcLR = -1;
01418 break;
01419 }
01420 }
01421 }
01422 }
01423
01424
01425 m_tuple4->write();
01426 }
01427
01428 unsigned segSize = tracks[i]->segments().length();
01429 for(unsigned kk = 0; kk < segSize; ++kk) {
01430 unsigned segLength = tracks[i]->segments()[kk]->links().length();
01431 AList<TMLink> ll = tracks[i]->segments()[kk]->links();
01432 int nSmall = 0;
01433 for(unsigned nn = 0; nn < segLength; ++nn) {
01434 double tmpX = ll[nn]->positionD().x();
01435 double tmpY = ll[nn]->positionD().y();
01436 double tmpA = tracks[i]->segments()[kk]->lineTsf().x();
01437 double tmpB = tracks[i]->segments()[kk]->lineTsf().y();
01438 double tmpC = tracks[i]->segments()[kk]->lineTsf().z();
01439 double dis = fabs(tmpA * tmpX + tmpB * tmpY + tmpC) / sqrt(tmpA * tmpA + tmpB * tmpB);
01440 double idealDis = maxdDistance(ll[nn]);
01441 if(fabs(dis-ll[nn]->cDrift())/idealDis<0.001 && nSmall<2) {
01442 ++nSmall;
01443 continue;
01444 }
01445 t9_times += fabs(dis-ll[nn]->cDrift())/idealDis;
01446 }
01447 t9_nSL = ll[0]->wire()->superLayerId();
01448 t9_nLinks = segLength;
01449 t9_mctheta = MC_TRACK_LRZVEC.theta();
01450 t9_times = t9_times / (t9_nLinks - nSmall);
01451 m_tuple9->write();
01452 }
01453
01454 t_mcphi = MC_TRACK_VEC.phi();
01455 t_mctheta = MC_TRACK_LRZVEC.theta();
01456 t_mcptot = MC_TRACK_VEC.mag()/1000.;
01457 t_mcpt = MC_TRACK_VEC.perp()/1000.;
01458 t_mcpz = MC_TRACK_LRZVEC.pz()/1000.;
01459 t_mct0 = MC_EVENT_TIME;
01460 t_nDigi = MC_DIGI_SIZE;
01461
01462 t_dr = tracks[i]->helix().dr();
01463 t_dz = tracks[i]->helix().dz();
01464 t_pt = tracks[i]->pt();
01465 t_ptot = tracks[i]->ptot();
01466 t_tanlmd = tracks[i]->helix().tanl();
01467 t_phi = tracks[i]->helix().phi0();
01468 t_chi2 = tracks[i]->chi2();
01469 t_nHits = tracks[i]->nLinks();
01470 t_nCores = tracks[i]->nCores();
01471 t_nSegs = tracks[i]->segments().length();
01472 t_ndf = tracks[i]->ndf();
01473 t_radius = tracks[i]->radius();
01474 t_evtNo = _nEvents;
01475 t_length = tracks.length();
01476 t_length2 = tracks2D.length();
01477
01478 t_dpt = tracks[i]->pt() - t_mcpt;
01479 t_dptot = tracks[i]->ptot() - t_mcptot;
01480 t_dlmd = atan(tracks[i]->helix().tanl()) - (pi/2 - t_mctheta);
01481 t_dphi = tracks[i]->helix().phi0() - t_mcphi;
01482
01483 t_t0 = t0_bes;
01484 t0_sta = t0Sta;
01485
01486 if (b_mcPar) {
01487 if (b_goodTrk && b_tuple) {
01488 unsigned mcTrackSize = MC_TRACK_NUM;
01489 unsigned recTrackSize = tracks.length();
01490 const unsigned matchSize = 20;
01491
01492
01493 int mLayers[matchSize];
01494 for (int ii=0; ii<matchSize; ii++)
01495 mLayers[ii] = 0;
01496 for(int jj = 0; jj < 43; ++jj) {
01497 int tmp[matchSize];
01498 for(unsigned kk = 0; kk < matchSize; ++kk)
01499 tmp[kk] = 0;
01500 for(int kk = 0; kk < 288; ++kk) {
01501 if (havedigi[jj][kk]<0) continue;
01502 tmp[havedigi[jj][kk]] = 1;
01503 }
01504 for(int kk = 0; kk < matchSize; ++kk)
01505 if(tmp[kk] == 1) ++mLayers[kk];
01506 }
01507
01508 unsigned trackSize = tracks[i]->nLinks();
01509 int trkIndex[43];
01510 int nLayers[matchSize];
01511 for (int j = 0; j < 43; ++j)
01512 trkIndex[j] = -99;
01513 for (int j = 0; j < matchSize; ++j)
01514 nLayers[j] = 0;
01515
01516 for (int j = 0; j < trackSize; ++j) {
01517 unsigned lId = tracks[i]->links()[j]->wire()->layerId();
01518 unsigned loId = tracks[i]->links()[j]->wire()->localId();
01519 if (havedigi[lId][loId] >= 0) trkIndex[lId] = havedigi[lId][loId];
01520 }
01521 for (int j = 0; j < 43; ++j)
01522 if (trkIndex[j] >= 0) ++nLayers[trkIndex[j]];
01523
01524 for (int j = 0; j < matchSize; ++j) {
01525
01526 if (nLayers[j]==0) continue;
01527 if ((float)nLayers[j]/(float)mLayers[j] > 0.51) {
01528 t_gdNLayers = nLayers[j];
01529 t_mcNLayers = mLayers[j];
01530 t_good_theta = MC_TRACK_LRZVEC.theta();
01531 ++nGood;
01532 break;
01533 }
01534 }
01535 if (t_good_theta == 0.) {
01536 int tmpLayers = 0;
01537 int tmpi = -1;
01538 for (int j = 0; j < matchSize; ++j) {
01539 if (nLayers[j]==0) continue;
01540 if (nLayers[j] > tmpLayers) {
01541 tmpLayers = nLayers[j];
01542 tmpi = j;
01543 }
01544 }
01545 if (tmpi != -1) {
01546 t_bestNLayers = nLayers[tmpi];
01547 t_bestMcNLayers = mLayers[tmpi];
01548 }
01549 else {
01550 t_bestNLayers = -1;
01551 t_bestMcNLayers = -1;
01552 }
01553 }
01554 else {
01555 t_bestNLayers = -2;
01556 t_bestMcNLayers = -2;
01557 }
01558 }
01559 }
01560 m_tuple->write();
01561 }
01562
01563 t3_mctheta = MC_TRACK_LRZVEC.theta();
01564 t3_mcptot = MC_TRACK_VEC.mag()/1000.;
01565 t3_mcpt = MC_TRACK_VEC.perp()/1000.;
01566
01567 t3_nDigi = MC_DIGI_SIZE;
01568 t3_t0 = t0_bes;
01569 t3_t0Sta = t0Sta;
01570 t3_goodLength = nGood;
01571 t3_length = tracks.length();
01572 t3_finalLength = num_par;
01573 m_tuple31->write();
01574
01575 for (unsigned i = 0; i < tracks2D.length(); i++) {
01576 t2_mctheta = MC_TRACK_LRZVEC.theta();
01577 t2_mcpt = MC_TRACK_VEC.perp()/1000.;
01578
01579 t2_nDigi = MC_DIGI_SIZE;
01580 t2_nHits = tracks2D[i]->nLinks();
01581 t2_nSegs = tracks2D[i]->segments().length();
01582 t2_chi2 = tracks2D[i]->chi2();
01583 t2_ndf = tracks2D[i]->ndf();
01584 t2_radius = tracks2D[i]->radius();
01585 t2_evtNo = _nEvents;
01586 t2_length = tracks.length();
01587 t2_length2 = tracks2D.length();
01588 m_tuple2->write();
01589 }
01590 }
01591
01592 double
01593 TrkReco::maxdDistance(TMLink *l) const{
01594 unsigned sl = l->wire()->superLayerId();
01595 unsigned lId = l->wire()->layerId();
01596 double _averR[11] = {9.7, 14.5, 22.14, 28.62, 35.1, 42.39, 48.87, 55.35, 61.83, 69.12, 74.79};
01597 double _averR2[43] = { 7.885, 9.07, 10.29, 11.525,
01598 12.72, 13.875, 15.01, 16.16,
01599 19.7, 21.3, 22.955, 24.555,
01600 26.215, 27.82, 29.465, 31.06,
01601 32.69, 34.265, 35.875, 37.455,
01602 39.975, 41.52, 43.12, 44.76,
01603 46.415, 48.02, 49.685, 51.37,
01604 53.035, 54.595, 56.215, 57.855,
01605 59.475, 60.995, 62.565, 64.165,
01606 66.68, 68.285, 69.915, 71.57,
01607 73.21, 74.76, 76.345};
01608 double radius = _averR2[lId];
01609 const double singleSigma = l->dDrift();
01610 return 2 * singleSigma / (radius * radius);
01611 }
01612
01613
01614
01615
01616
01617
01618
01619
01620
01621
01622
01623
01624
01625
01626
01627
01628
01629
01630
01631
01632
01633
01634
01635
01636
01637
01638
01639
01640
01641
01642
01643
01644
01645
01646
01647
01648
01649
01650
01651
01652
01653
01654
01655
01656
01657
01658
01659
01660
01661
01662
01663
01664
01665
01666
01667
01668
01669
01670
01671
01672
01673
01674
01675
01676
01677
01678
01679
01680
01681
01682
01683
01684
01685
01686
01687
01688
01689
01690
01691
01692
01693
01694
01695
01696
01697
01698
01699
01700
01701
01702
01703
01704
01705
01706
01707
01708
01709
01710
01711
01712
01713
01714
01715
01716
01717
01718
01719
01720
01721
01722
01723
01724
01725
01726
01727
01728
01729
01730
01731
01732
01733
01734
01735
01736
01737
01738
01739
01740
01741
01742
01743
01744
01745
01746
01747
01748
01749
01750
01751
01752
01753
01754
01755
01756
01757
01758
01759
01760
01761
01762
01763
01764
01765
01766
01767
01768
01769
01770
01771
01772
01773
01774
01775
01776
01777
01778
01779
01780
01781
01782
01783
01784
01785
01786
01787
01788
01789
01790
01791
01792
01793
01794
01795
01796