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