13#ifdef TRKRECO_DEBUG_DETAIL
21#include "TrkReco/TMDC.h"
22#include "TrkReco/TMDCLayer.h"
23#include "TrkReco/TMDCUtil.h"
24#include "TrkReco/TMDCWire.h"
25#include "TrkReco/TMDCWireHit.h"
26#include "TrkReco/TMDCWireHitMC.h"
27#include "TrkReco/TMLink.h"
28#include "TrkReco/TTrack.h"
29#include "TrkReco/TTrackHEP.h"
33#include "MdcGeomSvc/IMdcGeomSvc.h"
34#include "MdcTables/MdcTables.h"
38#include "EventModel/Event.h"
39#include "GaudiKernel/Bootstrap.h"
40#include "GaudiKernel/IDataProviderSvc.h"
41#include "GaudiKernel/IMessageSvc.h"
42#include "GaudiKernel/ISvcLocator.h"
43#include "GaudiKernel/MsgStream.h"
44#include "GaudiKernel/PropertyMgr.h"
45#include "GaudiKernel/SmartDataPtr.h"
46#include "GaudiKernel/StatusCode.h"
80 if ( !_cdc ) _cdc =
new TMDC(
version );
85 if ( !_cdc ) _cdc =
new TMDC(
"1.0" );
89TMDC::TMDC(
const std::string& version )
92 , _cdcVersion( version )
103 , _nSuperLayers( 11 )
105 , _newCdc( version ==
"-2000" ) {
110 std::cout <<
"TMDC ... MDC version = " << _cdcVersion << std::endl;
123 StatusCode sc = Gaudi::svcLocator()->service(
"MdcGeomSvc", mdcGeomSvc );
124 if ( sc == StatusCode::SUCCESS )
127 for (
unsigned i = 0; i < _nSuperLayers; i++ ) _superLayers.append(
new AList<TMDCLayer> );
131 for (
unsigned i = 0; i < _nLayers; i++ )
139 _layers.append(
layer );
143 for (
unsigned j = 0; j < l->
NCell(); j++ )
162 if ( msg.find(
"name" ) != std::string::npos || msg.find(
"version" ) != std::string::npos ||
163 msg.find(
"detail" ) != std::string::npos || msg ==
"" )
164 { std::cout <<
name() <<
"(" <<
version() <<
") "; }
165 if ( msg.find(
"detail" ) != std::string::npos || msg.find(
"state" ) != std::string::npos )
166 { std::cout <<
"Debug Level=" << _debugLevel; }
167 std::cout << std::endl;
169 std::string tab(
" " );
171 if ( msg ==
"" || msg.find(
"geometry" ) != std::string::npos )
175 unsigned nLayer = _layers.length();
176 std::cout <<
" version : " <<
version() << std::endl;
177 std::cout <<
" cdc version: " <<
cdcVersion() << std::endl;
178 std::cout <<
" # of wires : " << _wires.length() << std::endl;
179 std::cout <<
" # of layers: " << nLayer << std::endl;
180 std::cout <<
" super layer information" << std::endl;
181 std::cout <<
" # of super layers = " << _superLayers.length() << std::endl;
183 std::cout <<
" layer information" << std::endl;
184 for (
unsigned i = 0; i < _nLayers; i++ ) _layers[i]->
dump(
"", tab );
186 std::cout <<
" wire information" << std::endl;
187 for (
unsigned i = 0; i < _nWires; i++ ) ( _wires[i] )->dump(
"neighbor", tab );
191 if ( msg.find(
"hits" ) != std::string::npos )
193 std::cout <<
" hits : " << _hits.length() << std::endl;
194 for (
unsigned i = 0; i < _hits.length(); i++ ) { _hits[i]->dump(
"state mc", tab ); }
196 if ( msg.find(
"axialHits" ) != std::string::npos )
198 std::cout <<
" hits : " << _axialHits.length() << std::endl;
199 for (
unsigned i = 0; i < _axialHits.length(); i++ )
200 { _axialHits[i]->dump(
"state mc", tab ); }
202 if ( msg.find(
"stereoHits" ) != std::string::npos )
204 std::cout <<
" hits : " << _stereoHits.length() << std::endl;
205 for (
unsigned i = 0; i < _stereoHits.length(); i++ )
206 { _stereoHits[i]->dump(
"state mc", tab ); }
224 std::cout <<
"r,phi = " << r <<
"," << p << std::endl;
232 if ( !l )
return NULL;
235 if ( geo->
Radius() / 10 + geo->
RCSiz2() / 10 < r ) ++id;
236 else if ( geo->
Radius() / 10 - geo->
RCSiz1() / 10 > r ) --id;
241 float dPhi = 2. *
M_PI / float( l->
nWires() );
243 unsigned localId = unsigned(
phi( p ) / dPhi );
249 while (
TMDCWire*
w = _wires[i++] )
w->clear();
251 _hitWires.removeAll();
252 _axialHits.removeAll();
253 _stereoHits.removeAll();
254 HepAListDeleteAll( _hits );
255 HepAListDeleteAll( _hitsMC );
256 HepAListDeleteAll( _badHits );
263 while (
TMDCWire*
w = _hitWires[i++] )
w->clear();
267 _hitWires.removeAll();
268 _axialHits.removeAll();
269 _stereoHits.removeAll();
270 HepAListDeleteAll( _hits );
271 HepAListDeleteAll( _hitsMC );
272 HepAListDeleteAll( _badHits );
288 cout <<
"size of MdcRecWirhit:" << nReccdc << endl;
290 for (
unsigned i = 0; i < nReccdc; i++ )
311 _hitWires.append(
w );
322 if (
w->axial() ) _axialHits.append( hit );
323 else _stereoHits.append( hit );
361 unsigned n = _hits.length();
362 unsigned j = ( whp->
id <
n ) ? whp->
id :
n;
366 if ( _hits[j]->reccdc() == whp )
385 _hitsMC.append( hit );
387 if ( wh ) wh->
mc( hit );
393 if ( hep ) hep->_hits.append( hit );
396 std::cout <<
"TMDC::updateMC !!! mission impossible" << std::endl;
397 std::cout <<
" This error will cause TrkReco crush";
398 std::cout << std::endl;
399#ifdef TRKRECO_DEBUG_DETAIL
402 std::cout <<
" TTrackHEP list length = ";
412 unsigned n = _hits.length();
414 for (
unsigned i = 0; i <
n; i++ )
418 unsigned state = h->
state();
422 for (
unsigned j = 0; j < 6; j++ ) neighbor[j] =
w->neighbor( j );
431 unsigned pattern = 0;
432 for (
unsigned j = 0; j < 6; j++ )
435 if ( neighbor[j]->hit() ) pattern += ( 1 << j );
447 if ( ( hr2 == 0 ) && ( hr1 != 0 ) && ( hl1 == 0 ) ||
448 ( hl2 == 0 ) && ( hl1 != 0 ) && ( hr1 == 0 ) )
454 bool previous =
false;
456 if ( neighbor[0] == 0 ) previous =
true;
459 if ( ( neighbor[0]->hit() ) || neighbor[1]->hit() ) previous =
true;
461 if ( neighbor[5] == 0 ) next =
true;
464 if ( ( neighbor[4]->hit() ) || neighbor[5]->hit() ) next =
true;
470 if ( ( pattern == 34 ) || ( pattern == 42 ) || ( pattern == 40 ) || ( pattern == 10 ) ||
471 ( pattern == 35 ) || ( pattern == 50 ) )
473 else if ( ( pattern == 17 ) || ( pattern == 21 ) || ( pattern == 20 ) ||
474 ( pattern == 5 ) || ( pattern == 19 ) || ( pattern == 49 ) )
483 if ( !mask )
return _axialHits;
485 std::cout <<
"TMDC::axialHits !!! unsupported mask given" << std::endl;
490 if ( !mask )
return _stereoHits;
492 std::cout <<
"TMDC::stereoHits !!! unsupported mask given" << std::endl;
497 if ( !mask )
return _hits;
499 std::cout <<
"TMDC::hits !!! unsupported mask given" << std::endl;
505 if ( _badHits.length() )
return _badHits;
510 for (
unsigned i = 0; i < nReccdc; i++ )
530 _badHits.append(
new TMDCWireHit(
w, h, _fudgeFactor ) );
538 return std::to_string(
layerId( wireId ) ) +
"=" + std::to_string(
localId( wireId ) );
539 return std::to_string(
layerId( wireId ) ) +
"-" + std::to_string(
localId( wireId ) );
543 if ( !
w )
return std::string(
"no such a wire" );
545 unsigned id =
w->Id();
561 if (
id < 40 )
return 0;
562 else if (
id < 84 )
return 1;
563 else if (
id < 132 )
return 2;
564 else if (
id < 188 )
return 3;
566 else if (
id < 252 )
return 4;
567 else if (
id < 324 )
return 5;
568 else if (
id < 404 )
return 6;
569 else if (
id < 484 )
return 7;
571 else if (
id < 560 )
return 8;
572 else if (
id < 636 )
return 9;
573 else if (
id < 724 )
return 10;
574 else if (
id < 812 )
return 11;
576 else if (
id < 912 )
return 12;
577 else if (
id < 1012 )
return 13;
578 else if (
id < 1124 )
return 14;
579 else if (
id < 1236 )
return 15;
581 else if (
id < 1364 )
return 16;
582 else if (
id < 1492 )
return 17;
583 else if (
id < 1632 )
return 18;
584 else if (
id < 1772 )
return 19;
586 else if (
id < 1932 )
return 20;
587 else if (
id < 2092 )
return 21;
588 else if (
id < 2252 )
return 22;
589 else if (
id < 2412 )
return 23;
591 else if (
id < 2604 )
return 24;
592 else if (
id < 2796 )
return 25;
593 else if (
id < 2988 )
return 26;
594 else if (
id < 3180 )
return 27;
596 else if (
id < 3388 )
return 28;
597 else if (
id < 3596 )
return 29;
598 else if (
id < 3804 )
return 30;
599 else if (
id < 4012 )
return 31;
601 else if (
id < 4252 )
return 32;
602 else if (
id < 4492 )
return 33;
603 else if (
id < 4732 )
return 34;
604 else if (
id < 4972 )
return 35;
606 else if (
id < 5228 )
return 36;
607 else if (
id < 5484 )
return 37;
608 else if (
id < 5740 )
return 38;
609 else if (
id < 5996 )
return 39;
611 else if (
id < 6284 )
return 40;
612 else if (
id < 6572 )
return 41;
613 else if (
id < 6860 )
return 42;
619 if ( !
w )
return 9999;
633 if (
id < 40 )
return id;
634 else if (
id < 84 )
return id - 40;
635 else if (
id < 132 )
return id - 84;
636 else if (
id < 188 )
return id - 132;
638 else if (
id < 252 )
return id - 188;
639 else if (
id < 324 )
return id - 252;
640 else if (
id < 404 )
return id - 324;
641 else if (
id < 484 )
return id - 404;
643 else if (
id < 560 )
return id - 484;
644 else if (
id < 636 )
return id - 560;
645 else if (
id < 724 )
return id - 636;
646 else if (
id < 812 )
return id - 724;
648 else if (
id < 912 )
return id - 812;
649 else if (
id < 1012 )
return id - 912;
650 else if (
id < 1124 )
return id - 1012;
651 else if (
id < 1236 )
return id - 1124;
653 else if (
id < 1364 )
return id - 1236;
654 else if (
id < 1492 )
return id - 1364;
655 else if (
id < 1632 )
return id - 1492;
656 else if (
id < 1772 )
return id - 1632;
658 else if (
id < 1932 )
return id - 1772;
659 else if (
id < 2092 )
return id - 1932;
660 else if (
id < 2252 )
return id - 2092;
661 else if (
id < 2412 )
return id - 2252;
663 else if (
id < 2604 )
return id - 2412;
664 else if (
id < 2796 )
return id - 2604;
665 else if (
id < 2988 )
return id - 2796;
666 else if (
id < 3180 )
return id - 2988;
668 else if (
id < 3388 )
return id - 3180;
669 else if (
id < 3596 )
return id - 3388;
670 else if (
id < 3804 )
return id - 3596;
671 else if (
id < 4012 )
return id - 3804;
673 else if (
id < 4252 )
return id - 4012;
674 else if (
id < 4492 )
return id - 4252;
675 else if (
id < 4732 )
return id - 4492;
676 else if (
id < 4972 )
return id - 4732;
678 else if (
id < 5228 )
return id - 4972;
679 else if (
id < 5484 )
return id - 5228;
680 else if (
id < 5740 )
return id - 5484;
681 else if (
id < 5996 )
return id - 5740;
683 else if (
id < 6284 )
return id - 5996;
684 else if (
id < 6572 )
return id - 6284;
685 else if (
id < 6860 )
return id - 6572;
706 if ( !
w )
return 9999;
720 if (
id < 188 )
return 0;
721 else if (
id < 484 )
return 1;
722 else if (
id < 812 )
return 2;
723 else if (
id < 1236 )
return 3;
724 else if (
id < 1772 )
return 4;
725 else if (
id < 2412 )
return 5;
726 else if (
id < 3180 )
return 6;
727 else if (
id < 4012 )
return 7;
728 else if (
id < 4972 )
return 8;
729 else if (
id < 5996 )
return 9;
730 else if (
id < 6860 )
return 10;
750 if ( !
w )
return 9999;
751 unsigned id =
w->Id();
763 if ( !l )
return 9999;
764 unsigned id = l->
Id();
766 if (
id < 44 )
return id / 4;
785 if (
id < 6860 )
return layerId(
id ) % 4;
854 if ( !
w )
return 9999;
855 return w->Layer() % 4;
869 if ( !l )
return 9999;
870 unsigned id = l->
Id();
871 if (
id < 44 )
return id % 4;
940 if ( !l )
return 9999;
941 unsigned id = l->
Id();
943 if (
id < 8 )
return id;
944 else if (
id < 20 )
return id - 8;
945 else if (
id < 36 )
return id - 12;
946 else if (
id < 43 )
return id - 24;
956 if (
id < 12 )
return id + 8;
957 else if (
id < 19 )
return id + 24;
962 if (
id < 8 )
return id;
963 else if (
id < 24 )
return id + 12;
1003 float tl =
t.helix().a()[4];
1004 float f = sqrt( 1. + tl * tl );
1005 float s = fabs(
t.helix().curv() ) * fabs( l.
dPhi() ) *
f;
1006 float p =
f / fabs(
t.helix().a()[2] );
1011 if ( !(
flag && 2 ) ) t0Offset = 0.;
1018 double p[3] = { tp.x(), tp.y(), tp.z() };
1021 double x[3] = { onWire.x(), onWire.y(), onWire.z() };
1025 int prop = (
flag & 4 );
1029 if ( side == 0 ) side = -1;
1057 float tanl =
abs( p[2] ) / tp.perp();
1059 if ( ( tanl >= 0.0 ) && ( tanl < 0.5 ) ) c = -0.48 * tanl + 1.3;
1060 else if ( ( tanl >= 0.5 ) && ( tanl < 1.0 ) ) c = -0.28 * tanl + 1.2;
1061 else if ( ( tanl >= 1.0 ) && ( tanl < 1.5 ) ) c = -0.16 * tanl + 1.08;
HepGeom::Vector3D< double > HepVector3D
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
#define WireHitFindingValid
#define WireHitContinuous
#define WireHitPatternLeft
#define WireHitPatternRight
#define WireHitNeighborHit
HepGeom::Point3D< double > HepPoint3D
virtual const MdcGeoLayer *const Layer(unsigned id)=0
virtual const MdcGeoWire *const Wire(unsigned id)=0
static vector< MdcDat_mcwirhit > * getMdcDatMcwirhitCol(void)
double Radius(void) const
double RCSiz2(void) const
double RCSiz1(void) const
double Offset(void) const
static vector< MdcRec_wirhit > * getMdcRecWirhitCol(void)
A class to represent a wire layer.
unsigned superLayerId(void) const
returns super layer id.
const TMDCWire *const wire(int id) const
returns a pointer to a wire. 'id' can be negative or 'id' can be greater than 'nWires()'.
const MdcGeoLayer * geocdc(void) const
returns a pointer to GEOMDC_WIR.
unsigned nWires(void) const
returns # of wires.
A class to represent a MC wire hit in MDC.
struct MdcRec_wirhit * reccdc(void) const
returns a pointer to RECMDC_WIRHIT.
unsigned state(void) const
returns state.
float dDrift(unsigned) const
returns drift distance error.
float drift(unsigned) const
returns drift distance.
const TMDCWire *const wire(void) const
returns a pointer to a TMDCWire.
const TMDCWireHitMC *const mc(void) const
returns a pointer to TMDCWireHitMC.
A class to represent a wire in MDC.
unsigned id(void) const
returns id.
const TMDCWireHit *const hit(void) const
returns a pointer to a TMDCWireHit.
const TMDCWire *const neighbor(unsigned) const
returns a pointer to a neighbor wire.
void fastClear(void)
clears TMDC information.
static unsigned localId(unsigned wireId)
static unsigned axialStereoLayerId(const MdcGeoLayer *const)
const TMDCWire *const wire(unsigned wireId) const
returns a pointer to a wire. 0 will be returned if 'wireId' is invalid.
void updateMC(void)
updates TMDC information for MC.
static void driftDistance(TMLink &link, const TTrack &track, unsigned correctionFlag=0, float T0Offset=0.)
std::string cdcVersion(void) const
returns MDC version.
static unsigned superLayerId(unsigned wireId)
std::string name(void) const
returns name.
static unsigned layerId(unsigned wireId)
const AList< TMDCWireHit > & hits(unsigned mask=0) const
returns a list of TMDCWireHit. 'update()' must be called before calling this function.
static TMDC * getTMDC(void)
const AList< TMDCWireHit > & axialHits(unsigned mask=0) const
returns a list of axial hits. 'update()' must be called before calling this function.
unsigned nWires(void) const
const TMDCLayer *const layer(unsigned id) const
returns a pointer to a layer. 0 will be returned if 'id' is invalid.
static std::string wireName(unsigned wireId)
std::string version(void) const
returns version.
static float phi(float phi)
const AList< TMDCWireHit > & badHits(void)
returns bad hits(finding invalid hits).
const AList< TMDCWireHit > & stereoHits(unsigned mask=0) const
returns a list of stereo hits. 'update()' must be called before calling this function.
static unsigned localLayerId(unsigned wireId)
void dump(const std::string &message) const
dumps debug information.
const AList< TMDCLayer > *const superLayer(unsigned id) const
returns a pointer to a super-layer. 0 will be returned if 'id' is invalid.
void clear(void)
clears all TMDC information.
void classification(void)
classify hits.
A class to relate TMDCWireHit and TTrack objects.
const HepPoint3D & positionOnWire(void) const
returns the closest point on wire to a track.
const TMDCWireHit * hit(void) const
returns a pointer to a hit.
float dDrift(void) const
returns/sets drift distance error.
double dPhi(void) const
returns dPhi to the closest point.
float drift(void) const
returns/sets drift distance.
A class to represent a GEN_HEPEVT particle in tracking.
static const AList< TTrackHEP > & list(void)
returns a list of TTrackHEP's.
A class to represent a track in tracking.
virtual void update(void)
updates an object.
virtual bool updated(void) const
returns true if an object is updated.
virtual void clear(void)
clears an object.