BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
TMDC.cxx
Go to the documentation of this file.
1//-----------------------------------------------------------------------------
2// $Id: TMDC.cxx,v 1.19 2022/01/28 22:14:13 maqm Exp $
3//-----------------------------------------------------------------------------
4// Filename : TMDC.cc
5// Section : Tracking MDC
6// Owner : Yoshi Iwasaki
7// Email : yoshihito.iwasaki@kek.jp
8//-----------------------------------------------------------------------------
9// Description : A class to represent MDC.
10// See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
11//-----------------------------------------------------------------------------
12
13#ifdef TRKRECO_DEBUG_DETAIL
14# ifndef TRKRECO_DEBUG
15# define TRKRECO_DEBUG
16# endif
17#endif
18#include <iostream>
19// #include "CLHEP/String/Strings.h"
20// #include "belle.h"
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"
30// #include MDC_H
31// #include BELLETDF_H
32// #include "MdcRecGeo/MdcRecGeo.h"
33#include "MdcGeomSvc/IMdcGeomSvc.h"
34#include "MdcTables/MdcTables.h"
35// #include "tables/bestdf.h"
36
37// #include "GaudiKernel/Algorithm.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"
47// #include "MDCRawEvent/MdcDigi.h"
48// #include "TOFRawEvent/TofDigi.h"
49// #include "EMCRawEvent/CalDigi.h"
50// #include "McTruth/McKine.h"
51
52#include <iostream>
53#include <vector>
54
55using namespace std;
56using namespace Event;
57
59/*zsl
60extern "C" {
61 void calcdc_driftdist_(int *,
62 int *,
63 int *,
64 float[3],
65 float[3],
66 float *,
67 float *,
68 float *);
69 void calcdc_tof2_(int *, float *, float *, float *);
70};
71*/
72
73std::string TMDC::name( void ) const { return "TMDC"; }
74
75std::string TMDC::version( void ) const { return "2.18"; }
76
77TMDC* TMDC::_cdc = 0;
78
79TMDC* TMDC::getTMDC( const std::string& version ) {
80 if ( !_cdc ) _cdc = new TMDC( version );
81 return _cdc;
82}
83
84TMDC* TMDC::getTMDC( void ) {
85 if ( !_cdc ) _cdc = new TMDC( "1.0" );
86 return _cdc;
87}
88
89TMDC::TMDC( const std::string& version )
90 : _debugLevel( 0 )
91 , _fudgeFactor( 1. )
92 , _cdcVersion( version )
93 ,
94 // _cdcVersion(-2000),
95 // _nWires(BsCouTab(GEOMDC_WIRE)),
96 // _nSuperLayers(BsCouTab(GEOMDC_SUPER)),
97 // _nLayers(BsCouTab(GEOMDC_LAYER)),
98 // _geo(MdcRecGeo::getMdcRecGeo()),
99 // _nWires(_geo->getWireSize()),
100 // _nSuperLayers(_geo->getSuperLayerSize()),
101 // _nLayers(_geo->getLayerSize()),
102 _nWires( 6796 )
103 , _nSuperLayers( 11 )
104 , _nLayers( 43 )
105 , _newCdc( version == "-2000" ) {
106
107 GMDC = this;
108
109 //...Get MDC information...
110 std::cout << "TMDC ... MDC version = " << _cdcVersion << std::endl;
111
112 //...Check GEOMDC_WIRE...
113 /* if (BsCouTab(GEOMDC_LAYER) == 0) {
114 std::cout << "TMDC !!! Bank GEOMDC_LAYER not found" << std::endl;
115 std::cout << " !!! mission impossible" << std::endl;
116 std::cout << " !!! TrkReco will crash!" << std::endl;
117 std::cout << " !!! Please stop execution" << std::endl;
118 return;
119 }
120 */
121
122 IMdcGeomSvc* mdcGeomSvc;
123 StatusCode sc = Gaudi::svcLocator()->service( "MdcGeomSvc", mdcGeomSvc );
124 if ( sc == StatusCode::SUCCESS )
125 {
126 //...Create super-layers...
127 for ( unsigned i = 0; i < _nSuperLayers; i++ ) _superLayers.append( new AList<TMDCLayer> );
128
129 //...Loop over all wires layers...
130 unsigned nWires = 0;
131 for ( unsigned i = 0; i < _nLayers; i++ )
132 {
133 // struct geocdc_layer * l =
134 // (struct geocdc_layer *)
135 // BsGetEnt(GEOMDC_LAYER, i + 1, BBS_No_Index);
136
137 const MdcGeoLayer* l = mdcGeomSvc->Layer( i );
138 TMDCLayer* layer = new TMDCLayer( l );
139 _layers.append( layer );
140 _superLayers[layer->superLayerId()]->append( layer );
141
142 //...Loop over all wires in a layer...
143 for ( unsigned j = 0; j < l->NCell(); j++ )
144 {
145 const MdcGeoWire* w = mdcGeomSvc->Wire( l->Wirst() + j );
146 TMDCWire* tw = new TMDCWire( w, layer );
147 _wires.append( tw );
148 layer->append( tw );
149
150 ++nWires;
151 }
152 }
153 }
154 else { return; }
155}
156
157void TMDC::dump( const std::string& msg ) const {
158 // if ( msg.index("name") != -1
159 // || msg.index("version") != -1
160 // || msg.index("detail") != -1
161 // || msg == "") {
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;
168
169 std::string tab( " " );
170
171 if ( msg == "" || msg.find( "geometry" ) != std::string::npos )
172 {
173
174 //...Get information..."
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;
182
183 std::cout << " layer information" << std::endl;
184 for ( unsigned i = 0; i < _nLayers; i++ ) _layers[i]->dump( "", tab );
185
186 std::cout << " wire information" << std::endl;
187 for ( unsigned i = 0; i < _nWires; i++ ) ( _wires[i] )->dump( "neighbor", tab );
188
189 return;
190 }
191 if ( msg.find( "hits" ) != std::string::npos )
192 {
193 std::cout << " hits : " << _hits.length() << std::endl;
194 for ( unsigned i = 0; i < _hits.length(); i++ ) { _hits[i]->dump( "state mc", tab ); }
195 }
196 if ( msg.find( "axialHits" ) != std::string::npos )
197 {
198 std::cout << " hits : " << _axialHits.length() << std::endl;
199 for ( unsigned i = 0; i < _axialHits.length(); i++ )
200 { _axialHits[i]->dump( "state mc", tab ); }
201 }
202 if ( msg.find( "stereoHits" ) != std::string::npos )
203 {
204 std::cout << " hits : " << _stereoHits.length() << std::endl;
205 for ( unsigned i = 0; i < _stereoHits.length(); i++ )
206 { _stereoHits[i]->dump( "state mc", tab ); }
207 }
208}
209
210const TMDCWire* const TMDC::wire( unsigned layerId, int localId ) const {
211 if ( layerId < _nLayers ) return _layers[layerId]->wire( localId );
212
213 return NULL;
214}
215
216const TMDCWire* const TMDC::wire( const HepPoint3D& p ) const {
217 float r = p.mag();
218 float phi = p.phi();
219 return wire( r, phi );
220}
221
222const TMDCWire* const TMDC::wire( float r, float p ) const {
223
224 std::cout << "r,phi = " << r << "," << p << std::endl;
225
226 unsigned id = 25;
227 bool ok = false;
228 const TMDCLayer* l;
229 while ( !ok )
230 {
231 l = layer( id );
232 if ( !l ) return NULL;
233
234 const MdcGeoLayer* geo = l->geocdc();
235 if ( geo->Radius() / 10 + geo->RCSiz2() / 10 < r ) ++id;
236 else if ( geo->Radius() / 10 - geo->RCSiz1() / 10 > r ) --id;
237 else ok = true;
238 }
239 // float dPhi = 2. * M_PI / float(l->nWires());
240 // if (l->geocdc()->m_offset > 0.) p -= dPhi / 2.;
241 float dPhi = 2. * M_PI / float( l->nWires() );
242 p -= l->geocdc()->Offset() / dPhi;
243 unsigned localId = unsigned( phi( p ) / dPhi );
244 return l->wire( localId );
245}
246
247void TMDC::clear( void ) {
248 unsigned i = 0;
249 while ( TMDCWire* w = _wires[i++] ) w->clear();
250
251 _hitWires.removeAll();
252 _axialHits.removeAll();
253 _stereoHits.removeAll();
254 HepAListDeleteAll( _hits );
255 HepAListDeleteAll( _hitsMC );
256 HepAListDeleteAll( _badHits );
257
259}
260
261void TMDC::fastClear( void ) {
262 unsigned i = 0;
263 while ( TMDCWire* w = _hitWires[i++] ) w->clear();
264 i = 0;
265 while ( TMDCWireHit* h = _badHits[i++] ) ( (TMDCWire*)h->wire() )->clear();
266
267 _hitWires.removeAll();
268 _axialHits.removeAll();
269 _stereoHits.removeAll();
270 HepAListDeleteAll( _hits );
271 HepAListDeleteAll( _hitsMC );
272 HepAListDeleteAll( _badHits );
273
275}
276
277void TMDC::update( bool mcAnalysis ) {
278 //...Already updated?...
279 // zsl if (updated()) return;
280
281 //...Clear old information...
282 fastClear();
283
284 //...Loop over RECMDC_WIRHIT bank...
285 // unsigned nReccdc = BsCouTab(RECMDC_WIRHIT);
286 unsigned nReccdc = MdcRecWirhitCol::getMdcRecWirhitCol()->size();
287#ifdef TRKRECO_DEBUG
288 cout << "size of MdcRecWirhit:" << nReccdc << endl;
289#endif
290 for ( unsigned i = 0; i < nReccdc; i++ )
291 {
292 // struct reccdc_wirhit * h =
293 // (struct reccdc_wirhit *)
294 // BsGetEnt(RECMDC_WIRHIT, i + 1, BBS_No_Index);
296 //...Check validity...
297 // zsl if (! (h->stat & WireHitFindingValid)) continue;
298
299 //...Obtain a pointer to GEOMDC...
300 // struct geocdc_wire * g =
301 // (struct geocdc_wire *)
302 // BsGetEnt(GEOMDC_WIRE, h->m_geo, BBS_No_Index);
303 const MdcGeoWire* g = h->geo;
304 // if (_newCdc)
305 // if (g->m_ID < 256)
306 // continue;
307
308 //...Get a pointer to a TMDCWire...
309 // TMDCWire * w = _wires[g->m_ID - 1];
310 TMDCWire* w = _wires[g->Id()];
311 _hitWires.append( w );
312
313 //...Create TMDCWireHit...
314 TMDCWireHit* hit = new TMDCWireHit( w, h, _fudgeFactor );
315 _hits.append( hit );
316
317 //... update the wirehit pointer in the TMDCWire (by zang shilei)
318 w->hit( hit );
319 // if (w->hit()) cout<<"test .......... right~~~~~"<<endl;
320
321 //...Axial or stereo...
322 if ( w->axial() ) _axialHits.append( hit );
323 else _stereoHits.append( hit );
324 }
325
326 //...Hit classfication...
328
329 //...MC information...
330 if ( mcAnalysis ) updateMC();
331
332 //...Update information...
334}
335
336void TMDC::updateMC( void ) {
337
338 //...Create TTrackHEP...
339 TTrackHEP::update();
340
341 //...Loop over DATMDC_MCWIRHIT bank...
342 unsigned n = 0;
343 // for (unsigned i = 0; i < BsCouTab(DATMDC_MCWIRHIT); i++) {
344 for ( unsigned i = 0; i < MdcDatMcwirhitCol::getMdcDatMcwirhitCol()->size(); i++ )
345 {
346 // struct datcdc_mcwirhit * h =
347 // (struct datcdc_mcwirhit *)
348 // BsGetEnt(DATMDC_MCWIRHIT, i + 1, BBS_No_Index);
350 //...Get a pointer to RECMDC_WIRHIT...
351 // reccdc_wirhit * whp =
352 // (reccdc_wirhit *) BsGetEnt(RECMDC_WIRHIT, h->m_dat, BBS_No_Index);
353 MdcRec_wirhit* whp = h->dat->rec;
354 //...Get TrkReco objects...
355 TMDCWireHit* wh = 0;
356 TMDCWire* w = 0;
357 if ( whp )
358 {
359 if ( whp->stat & WireHitFindingValid )
360 {
361 unsigned n = _hits.length();
362 unsigned j = ( whp->id < n ) ? whp->id : n;
363 while ( j )
364 {
365 --j;
366 if ( _hits[j]->reccdc() == whp )
367 {
368 wh = _hits[j];
369 w = _wires[wh->wire()->id()];
370 break;
371 }
372 }
373 }
374 }
375 if ( !w )
376 {
377 // geocdc_wire * g =
378 // (geocdc_wire *) BsGetEnt(GEOMDC_WIRE, h->m_geo, BBS_No_Index);
379 // w = _wires[g->m_ID - 1];
380 const MdcGeoWire* g = h->geo;
381 }
382
383 //...Create TMDCWireHitMC...
384 TMDCWireHitMC* hit = new TMDCWireHitMC( w, wh, h );
385 _hitsMC.append( hit );
386 w->hit( hit );
387 if ( wh ) wh->mc( hit );
388
389 //...TTrackHEP...
390 // TTrackHEP * hep = TTrackHEP::list()[h->m_hep - 1];
391 TTrackHEP* hep = TTrackHEP::list()[h->hep->id];
392 hit->_hep = hep;
393 if ( hep ) hep->_hits.append( hit );
394 else
395 {
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
400 // cout << " h->m_hep, h->m_hep -1 = " << h->m_hep;
401 // std::cout << ", " << h->m_hep - 1 << std::endl;
402 std::cout << " TTrackHEP list length = ";
403 std::cout << TTrackHEP::list().length() << std::endl;
404// BsShwDat(GEN_HEPEVT);
405// BsShwDat(DATMDC_MCWIRHIT);
406#endif
407 }
408 }
409}
410
412 unsigned n = _hits.length();
413
414 for ( unsigned i = 0; i < n; i++ )
415 {
416 TMDCWireHit* h = _hits[i];
417 const TMDCWire* w = h->wire();
418 unsigned state = h->state();
419
420 //...Cache pointers to a neighbor...
421 const TMDCWire* neighbor[6];
422 for ( unsigned j = 0; j < 6; j++ ) neighbor[j] = w->neighbor( j );
423
424 // output to test ...
425 // cout<<w->localId()<<endl;
426 // if(neighbor[0]) cout<<"inner local:"<<neighbor[0]->localId()<<",
427 //"<<neighbor[1]->localId()<<endl; if(neighbor[5]) cout<<"outer
428 // local:"<<neighbor[4]->localId()<<", "<<neighbor[5]->localId()<<endl;
429
430 //...Decide hit pattern...
431 unsigned pattern = 0;
432 for ( unsigned j = 0; j < 6; j++ )
433 {
434 if ( neighbor[j] )
435 if ( neighbor[j]->hit() ) pattern += ( 1 << j );
436 }
437 state |= ( pattern << WireHitNeighborHit );
438
439 //...Check isolation...
440 const TMDCWireHit* hr1 = neighbor[2]->hit();
441 const TMDCWireHit* hl1 = neighbor[3]->hit();
442 if ( ( hr1 == 0 ) && ( hl1 == 0 ) ) { state |= WireHitIsolated; }
443 else
444 {
445 const TMDCWireHit* hr2 = neighbor[2]->neighbor( 2 )->hit();
446 const TMDCWireHit* hl2 = neighbor[3]->neighbor( 3 )->hit();
447 if ( ( hr2 == 0 ) && ( hr1 != 0 ) && ( hl1 == 0 ) ||
448 ( hl2 == 0 ) && ( hl1 != 0 ) && ( hr1 == 0 ) )
449 state |= WireHitIsolated;
450 }
451
452 //...Check continuation...
453 unsigned superLayer = w->superLayerId();
454 bool previous = false;
455 bool next = false;
456 if ( neighbor[0] == 0 ) previous = true;
457 else
458 {
459 if ( ( neighbor[0]->hit() ) || neighbor[1]->hit() ) previous = true;
460 }
461 if ( neighbor[5] == 0 ) next = true;
462 else
463 {
464 if ( ( neighbor[4]->hit() ) || neighbor[5]->hit() ) next = true;
465 }
466 // if (previous && next) state |= WireHitContinuous;
467 if ( previous || next ) state |= WireHitContinuous;
468
469 //...Solve LR locally...
470 if ( ( pattern == 34 ) || ( pattern == 42 ) || ( pattern == 40 ) || ( pattern == 10 ) ||
471 ( pattern == 35 ) || ( pattern == 50 ) )
472 state |= WireHitPatternRight;
473 else if ( ( pattern == 17 ) || ( pattern == 21 ) || ( pattern == 20 ) ||
474 ( pattern == 5 ) || ( pattern == 19 ) || ( pattern == 49 ) )
475 state |= WireHitPatternLeft;
476
477 //...Store it...
478 h->state( state );
479 }
480}
481
482const AList<TMDCWireHit>& TMDC::axialHits( unsigned mask ) const {
483 if ( !mask ) return _axialHits;
484 else if ( mask == WireHitFindingValid ) return _axialHits;
485 std::cout << "TMDC::axialHits !!! unsupported mask given" << std::endl;
486 return _axialHits;
487}
488
489const AList<TMDCWireHit>& TMDC::stereoHits( unsigned mask ) const {
490 if ( !mask ) return _stereoHits;
491 else if ( mask == WireHitFindingValid ) return _stereoHits;
492 std::cout << "TMDC::stereoHits !!! unsupported mask given" << std::endl;
493 return _stereoHits;
494}
495
496const AList<TMDCWireHit>& TMDC::hits( unsigned mask ) const {
497 if ( !mask ) return _hits;
498 else if ( mask == WireHitFindingValid ) return _hits;
499 std::cout << "TMDC::hits !!! unsupported mask given" << std::endl;
500 return _hits;
501}
502
504 if ( !updated() ) update();
505 if ( _badHits.length() ) return _badHits;
506
507 //...Loop over RECMDC_WIRHIT bank...
508 // unsigned nReccdc = BsCouTab(RECMDC_WIRHIT);
509 unsigned nReccdc = MdcRecWirhitCol::getMdcRecWirhitCol()->size();
510 for ( unsigned i = 0; i < nReccdc; i++ )
511 {
512 // struct reccdc_wirhit * h =
513 // (struct reccdc_wirhit *)
514 // BsGetEnt(RECMDC_WIRHIT, i + 1, BBS_No_Index);
516
517 //...Check validity...
518 if ( h->stat & WireHitFindingValid ) continue;
519
520 //...Obtain a pointer to GEOMDC...
521 // geocdc_wire * g =
522 // (geocdc_wire *) BsGetEnt(GEOMDC_WIRE, h->m_geo, BBS_No_Index);
523 const MdcGeoWire* g = h->geo;
524
525 //...Get a pointer to a TMDCWire...
526 // TMDCWire * w = _wires[g->m_ID - 1];
527 TMDCWire* w = _wires[g->Id()];
528
529 //...Create TMDCWireHit...
530 _badHits.append( new TMDCWireHit( w, h, _fudgeFactor ) );
531 }
532
533 return _badHits;
534}
535
536std::string TMDC::wireName( unsigned wireId ) {
537 if ( superLayerId( wireId ) % 2 )
538 return std::to_string( layerId( wireId ) ) + "=" + std::to_string( localId( wireId ) );
539 return std::to_string( layerId( wireId ) ) + "-" + std::to_string( localId( wireId ) );
540}
541
542std::string TMDC::wireName( const MdcGeoWire* const w ) {
543 if ( !w ) return std::string( "no such a wire" );
544 // unsigned id = w->m_ID - 1;
545 unsigned id = w->Id();
546 return wireName( id );
547}
548
549// std::string
550// TMDC::wireName(const Geocdc_wire * const w) {
551// if (! w) return std::string("no such a wire");
552// unsigned id = w->get_ID() - 1;
553// return wireName(id);
554// }
555
556unsigned TMDC::layerId( unsigned id ) {
557 // int length = _wires.length();
558 // if(id < length) return _wires[id]->geocdc()->Layer();
559 // return 9999;
560
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;
565
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;
570
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;
575
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;
580
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;
585
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;
590
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;
595
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;
600
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;
605
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;
610
611 else if ( id < 6284 ) return 40;
612 else if ( id < 6572 ) return 41;
613 else if ( id < 6860 ) return 42;
614
615 return 9999;
616}
617
618unsigned TMDC::layerId( const MdcGeoWire* const w ) {
619 if ( !w ) return 9999;
620 // zsl unsigned id = w->Id();
621 // return layerId(id);
622 return w->Layer();
623}
624
625// unsigned
626// TMDC::layerId(const Geocdc_wire * const w) {
627// if (! w) return 9999;
628// unsigned id = w->get_ID() - 1;
629// return layerId(id);
630// }
631
632unsigned TMDC::localId( unsigned id ) {
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;
637
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;
642
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;
647
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;
652
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;
657
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;
662
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;
667
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;
672
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;
677
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;
682
683 else if ( id < 6284 ) return id - 5996;
684 else if ( id < 6572 ) return id - 6284;
685 else if ( id < 6860 ) return id - 6572;
686
687 return 9999;
688
689 /* if (id < 384) return id % 64;
690 else if (id < 624) return (id - 384) % 80;
691 else if (id < 1200) return (id - 624) % 96;
692 else if (id < 1584) return (id - 1200) % 128;
693 else if (id < 2304) return (id - 1584) % 144;
694 else if (id < 2944) return (id - 2304) % 160;
695 else if (id < 3904) return (id - 2944) % 176;
696 else if (id < 4736) return (id - 3904) % 208;
697 else if (id < 5936) return (id - 4736) % 240;
698 else if (id < 6960) return (id - 5936) % 256;
699 else if (id < 8400) return (id - 6960) % 288;
700
701 return 9999;
702 */
703}
704
705unsigned TMDC::localId( const MdcGeoWire* const w ) {
706 if ( !w ) return 9999;
707 // unsigned id = w->Id();
708 // return localId(id);
709 return w->Cell();
710}
711
712// unsigned
713// TMDC::localId(const Geocdc_wire * const w) {
714// if (! w) return 9999;
715// unsigned id = w->get_ID() - 1;
716// return localId(id);
717// }
718
719unsigned TMDC::superLayerId( unsigned id ) {
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;
731
732 return 9999;
733 /* if (id < 384) return 0;
734 else if (id < 624) return 1;
735 else if (id < 1200) return 2;
736 else if (id < 1584) return 3;
737 else if (id < 2304) return 4;
738 else if (id < 2944) return 5;
739 else if (id < 3904) return 6;
740 else if (id < 4736) return 7;
741 else if (id < 5936) return 8;
742 else if (id < 6960) return 9;
743 else if (id < 8400) return 10;
744
745 return 9999;
746 */
747}
748
749unsigned TMDC::superLayerId( const MdcGeoWire* const w ) {
750 if ( !w ) return 9999;
751 unsigned id = w->Id();
752 return superLayerId( id );
753}
754
755// unsigned
756// TMDC::superLayerId(const Geocdc_wire * const w) {
757// if (! w) return 9999;
758// unsigned id = w->get_ID() - 1;
759// return superLayerId(id);
760// }
761
762unsigned TMDC::superLayerId( const MdcGeoLayer* const l ) {
763 if ( !l ) return 9999;
764 unsigned id = l->Id();
765
766 if ( id < 44 ) return id / 4;
767 return 9999;
768 /* if (id < 6) return 0;
769 else if (id < 9) return 1;
770 else if (id < 15) return 2;
771 else if (id < 18) return 3;
772 else if (id < 23) return 4;
773 else if (id < 27) return 5;
774 else if (id < 32) return 6;
775 else if (id < 36) return 7;
776 else if (id < 41) return 8;
777 else if (id < 45) return 9;
778 else if (id < 50) return 10;
779
780 return 9999;
781 */
782}
783
784unsigned TMDC::localLayerId( unsigned id ) {
785 if ( id < 6860 ) return layerId( id ) % 4;
786 return 9999;
787
788 /* if (id < 64) return 0;
789 else if (id < 128) return 1;
790 else if (id < 192) return 2;
791 else if (id < 256) return 3;
792 else if (id < 320) return 4;
793 else if (id < 384) return 5;
794
795 else if (id < 464) return 0;
796 else if (id < 544) return 1;
797 else if (id < 624) return 2;
798
799 else if (id < 720) return 0;
800 else if (id < 816) return 1;
801 else if (id < 912) return 2;
802 else if (id < 1008) return 3;
803 else if (id < 1104) return 4;
804 else if (id < 1200) return 5;
805
806 else if (id < 1328) return 0;
807 else if (id < 1456) return 1;
808 else if (id < 1584) return 2;
809
810 else if (id < 1728) return 0;
811 else if (id < 1872) return 1;
812 else if (id < 2016) return 2;
813 else if (id < 2160) return 3;
814 else if (id < 2304) return 4;
815
816 else if (id < 2464) return 0;
817 else if (id < 2624) return 1;
818 else if (id < 2784) return 2;
819 else if (id < 2944) return 3;
820
821 else if (id < 3136) return 0;
822 else if (id < 3328) return 1;
823 else if (id < 3520) return 2;
824 else if (id < 3712) return 3;
825 else if (id < 3904) return 4;
826
827 else if (id < 4112) return 0;
828 else if (id < 4320) return 1;
829 else if (id < 4528) return 2;
830 else if (id < 4736) return 3;
831
832 else if (id < 4976) return 0;
833 else if (id < 5216) return 1;
834 else if (id < 5456) return 2;
835 else if (id < 5696) return 3;
836 else if (id < 5936) return 4;
837
838 else if (id < 6192) return 0;
839 else if (id < 6448) return 1;
840 else if (id < 6704) return 2;
841 else if (id < 6960) return 3;
842
843 else if (id < 7248) return 0;
844 else if (id < 7536) return 1;
845 else if (id < 7824) return 2;
846 else if (id < 8112) return 3;
847 else if (id < 8400) return 4;
848
849 return 9999;
850 */
851}
852
853unsigned TMDC::localLayerId( const MdcGeoWire* const w ) {
854 if ( !w ) return 9999;
855 return w->Layer() % 4;
856
857 // unsigned id = w->Id();
858 // return localLayerId(id);
859}
860
861// unsigned
862// TMDC::localLayerId(const Geocdc_wire * const w) {
863// if (! w) return 9999;
864// unsigned id = w->get_ID() - 1;
865// return localLayerId(id);
866// }
867
868unsigned TMDC::localLayerId( const MdcGeoLayer* const l ) {
869 if ( !l ) return 9999;
870 unsigned id = l->Id();
871 if ( id < 44 ) return id % 4;
872 return 9999;
873
874 /* if (id < 1) return 0;
875 else if (id < 2) return 1;
876 else if (id < 3) return 2;
877 else if (id < 4) return 3;
878 else if (id < 5) return 4;
879 else if (id < 6) return 5;
880
881 else if (id < 7) return 0;
882 else if (id < 8) return 1;
883 else if (id < 9) return 2;
884
885 else if (id < 10) return 0;
886 else if (id < 11) return 1;
887 else if (id < 12) return 2;
888 else if (id < 13) return 3;
889 else if (id < 14) return 4;
890 else if (id < 15) return 5;
891
892 else if (id < 16) return 0;
893 else if (id < 17) return 1;
894 else if (id < 18) return 2;
895
896 else if (id < 19) return 0;
897 else if (id < 20) return 1;
898 else if (id < 21) return 2;
899 else if (id < 22) return 3;
900 else if (id < 23) return 4;
901
902 else if (id < 24) return 0;
903 else if (id < 25) return 1;
904 else if (id < 26) return 2;
905 else if (id < 27) return 3;
906
907 else if (id < 28) return 0;
908 else if (id < 29) return 1;
909 else if (id < 30) return 2;
910 else if (id < 31) return 3;
911 else if (id < 32) return 4;
912
913 else if (id < 33) return 0;
914 else if (id < 34) return 1;
915 else if (id < 35) return 2;
916 else if (id < 36) return 3;
917
918 else if (id < 37) return 0;
919 else if (id < 38) return 1;
920 else if (id < 39) return 2;
921 else if (id < 40) return 3;
922 else if (id < 41) return 4;
923
924 else if (id < 42) return 0;
925 else if (id < 43) return 1;
926 else if (id < 44) return 2;
927 else if (id < 45) return 3;
928
929 else if (id < 46) return 0;
930 else if (id < 47) return 1;
931 else if (id < 48) return 2;
932 else if (id < 49) return 3;
933 else if (id < 50) return 4;
934
935 return 9999;
936 */
937}
938
939unsigned TMDC::axialStereoLayerId( const MdcGeoLayer* const l ) {
940 if ( !l ) return 9999;
941 unsigned id = l->Id();
942
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;
947
948 return 9999;
949}
950
951unsigned TMDC::layerId( unsigned as, unsigned id ) {
952
953 //...Axial case...
954 if ( as == 0 )
955 {
956 if ( id < 12 ) return id + 8;
957 else if ( id < 19 ) return id + 24;
958 return 9999;
959 }
960
961 //...Stereo case...
962 if ( id < 8 ) return id;
963 else if ( id < 24 ) return id + 12;
964 return 9999;
965}
966
967std::string TMDC::wireName( const MdcRec_wirhit& h ) {
968 // geocdc_wire * g = (geocdc_wire *) BsGetEnt(GEOMDC_WIRE,
969 // h.m_geo,
970 // BBS_No_Index);
971 const MdcGeoWire* g = h.geo;
972 return wireName( g );
973}
974
975void TMDC::driftDistance( TMLink& l, const TTrack& t, unsigned flag, float t0Offset ) {
976
977 //...No correction...
978 if ( flag == 0 )
979 {
980 if ( l.hit() )
981 {
982 l.drift( 0, l.hit()->drift( 0 ) );
983 l.drift( 1, l.hit()->drift( 1 ) );
984 l.dDrift( 0, l.hit()->dDrift( 0 ) );
985 l.dDrift( 1, l.hit()->dDrift( 1 ) );
986 }
987 else
988 {
989 l.drift( 0, 0. );
990 l.drift( 1, 0. );
991 l.dDrift( 0, 0. );
992 l.dDrift( 1, 0. );
993 }
994
995 return;
996 }
997
998 //...TOF correction...
999 float tof = 0.;
1000 if ( flag && 1 )
1001 {
1002 int imass = 3;
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] );
1007 // zsl calcdc_tof2_(& imass, & p, & s, & tof);
1008 }
1009
1010 //...T0 correction....
1011 if ( !( flag && 2 ) ) t0Offset = 0.;
1012
1013 //...Propagation corrections...
1014 const TMDCWireHit& h = *l.hit();
1015 int wire = h.wire()->id();
1016 HepVector3D tp = t.helix().momentum( l.dPhi() );
1017 // float p[3] = {tp.x(), tp.y(), tp.z()};
1018 double p[3] = { tp.x(), tp.y(), tp.z() };
1019 const HepPoint3D& onWire = l.positionOnWire();
1020 // maqm float x[3] = {onWire.x(), onWire.y(), onWire.z()};
1021 double x[3] = { onWire.x(), onWire.y(), onWire.z() };
1022 float time = h.reccdc()->tdc + t0Offset - tof;
1023 float dist = 1.e6;// add initialize 25-05-15
1024 float edist = 1.e9;// add initialize 25-05-15
1025 int prop = ( flag & 4 );
1026
1027 //...Calculation with left side...
1028 int side = -1;
1029 if ( side == 0 ) side = -1;
1030 // zsl calcdc_driftdist_(& prop,
1031 // & wire,
1032 // & side,
1033 // p,
1034 // x,
1035 // & time,
1036 // & dist,
1037 // & edist);
1038 l.drift( 0, dist );
1039 l.dDrift( 0, edist );
1040
1041 //...Calculation with left side...
1042 side = 1;
1043 // zsl calcdc_driftdist_(& prop,
1044 // & wire,
1045 // & side,
1046 // p,
1047 // x,
1048 // & time,
1049 // & dist,
1050 // & edist);
1051 l.drift( 1, dist );
1052 l.dDrift( 1, edist );
1053
1054 //...tan(lambda) correction...
1055 if ( flag && 8 )
1056 {
1057 float tanl = abs( p[2] ) / tp.perp();
1058 float c;
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;
1062 else c = 0.84;
1063
1064 l.dDrift( 0, l.dDrift( 0 ) * c );
1065 l.dDrift( 1, l.dDrift( 1 ) * c );
1066 }
1067}
HepGeom::Vector3D< double > HepVector3D
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
const Int_t n
Double_t time
double w
XmlRpcServer s
HepGeom::Point3D< double > HepPoint3D
#define M_PI
Definition TConstant.h:4
TMDC * GMDC
Definition TMDC.cxx:58
virtual const MdcGeoLayer *const Layer(unsigned id)=0
virtual const MdcGeoWire *const Wire(unsigned id)=0
static vector< MdcDat_mcwirhit > * getMdcDatMcwirhitCol(void)
static vector< MdcRec_wirhit > * getMdcRecWirhitCol(void)
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()'.
Definition TMDCLayer.cxx:56
const MdcGeoLayer * geocdc(void) const
returns a pointer to GEOMDC_WIR.
unsigned nWires(void) const
returns # of wires.
struct MdcRec_wirhit * reccdc(void) const
returns a pointer to RECMDC_WIRHIT.
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.
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.
Definition TMDCWire.cxx:94
void fastClear(void)
clears TMDC information.
Definition TMDC.cxx:261
static unsigned localId(unsigned wireId)
Definition TMDC.cxx:632
static unsigned axialStereoLayerId(const MdcGeoLayer *const)
Definition TMDC.cxx:939
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.
Definition TMDC.cxx:336
static void driftDistance(TMLink &link, const TTrack &track, unsigned correctionFlag=0, float T0Offset=0.)
Definition TMDC.cxx:975
std::string cdcVersion(void) const
returns MDC version.
static unsigned superLayerId(unsigned wireId)
Definition TMDC.cxx:719
std::string name(void) const
returns name.
Definition TMDC.cxx:73
static unsigned layerId(unsigned wireId)
Definition TMDC.cxx:556
const AList< TMDCWireHit > & hits(unsigned mask=0) const
returns a list of TMDCWireHit. 'update()' must be called before calling this function.
Definition TMDC.cxx:496
static TMDC * getTMDC(void)
Definition TMDC.cxx:84
const AList< TMDCWireHit > & axialHits(unsigned mask=0) const
returns a list of axial hits. 'update()' must be called before calling this function.
Definition TMDC.cxx:482
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)
Definition TMDC.cxx:536
std::string version(void) const
returns version.
Definition TMDC.cxx:75
const AList< TMDCWireHit > & badHits(void)
returns bad hits(finding invalid hits).
Definition TMDC.cxx:503
const AList< TMDCWireHit > & stereoHits(unsigned mask=0) const
returns a list of stereo hits. 'update()' must be called before calling this function.
Definition TMDC.cxx:489
static unsigned localLayerId(unsigned wireId)
Definition TMDC.cxx:784
void dump(const std::string &message) const
dumps debug information.
Definition TMDC.cxx:157
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.
Definition TMDC.cxx:247
void classification(void)
classify hits.
Definition TMDC.cxx:411
A class to represent a GEN_HEPEVT particle in tracking.
static const AList< TTrackHEP > & list(void)
returns a list of TTrackHEP's.
Definition TTrackHEP.cxx:70
A class to represent a track in tracking.
virtual void update(void)
updates an object.
Definition TUpdater.cxx:33
virtual bool updated(void) const
returns true if an object is updated.
Definition TUpdater.cxx:52
virtual void clear(void)
clears an object.
int t()
Definition t.c:1