BOSS 6.6.4.p01
BESIII Offline Software System
Loading...
Searching...
No Matches
FTFinder.cxx
Go to the documentation of this file.
1// -*- C++ -*-
2//
3// Package: MdcFastTrkAlg
4// Module: FTFinder
5//
6// Description: fast track finder class for MdcFastTrkAlg
7//
8//
9
10#include <cmath>
11#include <iostream>
12#include <vector>
13
14#include "GaudiKernel/MsgStream.h"
15#include "GaudiKernel/ThreadGaudi.h"
16#include "GaudiKernel/AlgFactory.h"
17#include "GaudiKernel/ISvcLocator.h"
18#include "GaudiKernel/IMessageSvc.h"
19#include "GaudiKernel/IDataProviderSvc.h"
20#include "GaudiKernel/SmartDataPtr.h"
21#include "GaudiKernel/SmartRef.h"
22#include "GaudiKernel/SmartRefVector.h"
23#include "GaudiKernel/IDataProviderSvc.h"
24#include "GaudiKernel/PropertyMgr.h"
25#include "GaudiKernel/IJobOptionsSvc.h"
26#include "GaudiKernel/IMessageSvc.h"
27#include "GaudiKernel/Bootstrap.h"
28#include "GaudiKernel/StatusCode.h"
29#include "GaudiKernel/IDataManagerSvc.h"
30#include "GaudiKernel/IPartPropSvc.h"
31
32#include "MdcGeomSvc/IMdcGeomSvc.h"
33#include "MdcGeomSvc/MdcGeoWire.h"
34#include "MdcGeomSvc/MdcGeoLayer.h"
35#include "MdcGeomSvc/MdcGeoSuper.h"
36
37#include "Identifier/Identifier.h"
38#include "Identifier/MdcID.h"
39
40#include "EventModel/EventHeader.h"
41#include "MdcRawEvent/MdcDigi.h"
42#include "RawEvent/RawDataUtil.h"
43#include "EvTimeEvent/RecEsTime.h"
44#include "RawDataProviderSvc/IRawDataProviderSvc.h"
45#include "TrigEvent/TrigData.h"
46#include "TrigEvent/TrigGTD.h"
47#include "TrigEvent/TrigMdc.h"
48#include "TrigEvent/TrigGTDProvider.h"
49
50
51#include "MdcRecEvent/RecMdcTrack.h"
52#include "MdcRecEvent/RecMdcHit.h"
53#include "ReconEvent/ReconEvent.h"
54
55#include "MdcFastTrkAlg/FTFinder.h"
56#include "MdcFastTrkAlg/FTWire.h"
57#include "MdcFastTrkAlg/FTLayer.h"
58#include "MdcFastTrkAlg/FTSuperLayer.h"
59#include "MdcFastTrkAlg/FTSegment.h"
60#include "MdcFastTrkAlg/FTTrack.h"
61#include "MdcFastTrkAlg/MdcParameter.h"
62
63#ifndef OnlineMode
64#include "McTruth/McParticle.h"
65
66#include "AIDA/IAxis.h"
67#include "AIDA/IHistogram1D.h"
68#include "AIDA/IHistogram2D.h"
69#include "AIDA/IHistogram3D.h"
70#include "AIDA/IHistogramFactory.h"
71#endif
72
73
74using namespace Event;
75//...Globals...
76
77//MCTruth
78
79#ifndef OnlineMode
80
81extern NTuple::Item<long> g_ntrkMC, g_eventNo;
82extern NTuple::Array<float> g_theta0MC, g_phi0MC;
83extern NTuple::Array<float> g_pxMC, g_pyMC, g_pzMC, g_ptMC;
84
85//recon
86
87extern NTuple::Item<long> g_ntrk;
88extern NTuple::Array<float> g_px, g_py, g_pz, g_pt, g_p;
89extern NTuple::Array<float> g_phi, g_theta;
90extern NTuple::Array<float> g_vx, g_vy, g_vz;
91extern NTuple::Array<float> g_dr, g_phi0, g_kappa, g_dz, g_tanl;
92extern NTuple::Item<float> g_estime;
93
94extern IHistogram1D* g_ncellMC;
95extern IHistogram1D* g_ncell;
96extern IHistogram1D* g_naxialhit;
97extern IHistogram1D* g_nstereohit;
98extern IHistogram1D* g_nhit;
99extern IHistogram2D* g_hitmap[20];
100
101extern int num_2Dtrk, num_3Dtrk, num_finaltrk;
102
103#endif
104
105
106//float FTTrack::xtCoEff = 0.;
107//int FTTrack::evtTiming = 0;
108//float FTTrack::Testime = 0;
109
110//FTFinder * FTFinder::_tFinder = NULL;
111
112MdcParameter* FTFinder::param = MdcParameter::instance();
113
114//FTFinder *
115//FTFinder::GetPointer(void)
116//{
117// if (!_tFinder) _tFinder = new FTFinder;
118// return _tFinder;
119//}
120
121FTFinder::FTFinder()
122 :// findEventVertex(1),
123 // evtTimeCorr(1),
124
125 // minPt(0.07),
126 // minDr(7.5),
127 // tOffSet(0.),
128 // xtCoEff(0.0344),
129 // doIt(1),
130#ifndef OnlineMode
131 // mkMdst(true),
132#endif
133 // mkTds(true),
134 tOffSet(0.),
135 t0Estime(-999.),
136 tEstFlag(0),
137
138 _wire(NULL),
139 _layer(NULL),
140 _superLayer(NULL),
141 _tracks(*(new FTList<FTTrack *>(20))),
142 _linkedSegments(new FTList<FTSegment *>(6)),
143 _axialSegCollect(10),
144 _vx(0.),
145 _vy(0.),
146 _vz(0.),
147 _ExpNo(0),
148 _RunNo(0),
149 m_total_trk(0),
150 pivot(0,0,0)
151{
152 StatusCode scmgn = Gaudi::svcLocator()->service("MagneticFieldSvc",m_pmgnIMF);
153 if(scmgn!=StatusCode::SUCCESS) {
154 std::cout<< "Unable to open Magnetic field service"<<std::endl;
155 }
156}
157
158void
159FTFinder::init()
160{
161 if (!param->_doIt) return;
162// param->updateAlpha();
163
164 // Get the Particle Properties Service
165 IPartPropSvc* p_PartPropSvc;
166 static const bool CREATEIFNOTTHERE(true);
167 StatusCode PartPropStatus = Gaudi::svcLocator()->service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
168 if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) {
169 // log << MSG::ERROR << " Could not initialize Particle Properties Service" << endreq;
170 std::cerr << "Could not initialize Particle Properties Service" << std::endl;
171 // return PartPropStatus;
172 return;
173 }
174 m_particleTable = p_PartPropSvc->PDT();
175
176 // IRawDataProviderSvc* m_rawDataProviderSvc;
177 StatusCode RawData = Gaudi::svcLocator()-> service ("RawDataProviderSvc", m_rawDataProviderSvc, CREATEIFNOTTHERE);
178 if ( !RawData.isSuccess() ){
179 std::cerr << "Could not load RawDataProviderSvc!" << m_rawDataProviderSvc << endreq;
180 return;
181 }
182
183// IMdcGeomSvc* mdcGeomSvc;
184// StatusCode sc = Gaudi::svcLocator()->service("MdcGeomSvc", mdcGeomSvc);
185//
186// if (sc == StatusCode::SUCCESS) {
187// //mdcGeomSvc->Dump();
188//
189// } else {
190// return ;
191// }
192//
193// const int Nwire = mdcGeomSvc->getWireSize();
194// const int Nlyr = mdcGeomSvc->getLayerSize();
195// const int Nsup = mdcGeomSvc->getSuperLayerSize();
196//
197// if(!Nwire || !Nlyr){
198// std::cerr << "FTFINDER::GEOCDC not found (please put cdctable before l4)" << std::endl;
199// std::cerr << "JOB will stop" << std::endl;
200// exit(-1);
201// }
202//
203// if (!_wire) _wire = (FTWire *) malloc((Nwire+1) * sizeof(FTWire));
204// if (!_layer) _layer = (FTLayer *) malloc(Nlyr * sizeof(FTLayer));
205// if (!_superLayer) _superLayer =
206// (FTSuperLayer *) malloc(Nsup * sizeof(FTSuperLayer));
207//
208// if (!_wire || !_layer || !_superLayer){
209// std::cerr << "FTFINDER::Cannot allocate geometries" << std::endl;
210// std::cerr << "JOB will stop" << std::endl;
211// exit(-1);
212// }
213//
214// int superLayerID = 0;
215// int layerID = 0;
216// int localLayerID = 0;
217// int localWireID = 0;
218// int localID = 0;
219// int wireID;
220// MdcGeoLayer * layer_back = NULL;
221// MdcGeoSuper * superLayer_back = NULL;
222// int k = 0;
223// int Nlayer[12];
224// int Nlocal[12];
225// int NlocalWireID[44];
226//
227// for (wireID = 0;wireID <= Nwire; wireID++){
228// MdcGeoLayer * layer = (wireID==Nwire) ? NULL : mdcGeomSvc->Wire(wireID)->Lyr();
229// if (layer != layer_back){
230// layer_back = layer;
231// MdcGeoSuper * superLayer = (wireID==Nwire) ? NULL : mdcGeomSvc->Layer(layerID)->Sup();
232// if (superLayer != superLayer_back){
233// superLayer_back = superLayer;
234// Nlayer[k] = localLayerID;
235// Nlocal[k] = localID;
236// localLayerID = 0;
237// k++;
238// }
239// NlocalWireID[layerID] = localWireID;
240// localID = 0;
241// localWireID = 0;
242// layerID++;
243// localLayerID++;
244// }
245// localID++;
246// localWireID++;
247// }
248//
249//
250// superLayerID = -1;
251// layerID = -1;
252// localLayerID = 0;
253// localID = 0;
254// layer_back = NULL;
255// superLayer_back = NULL;
256// for (wireID = 0;wireID < Nwire; wireID++){
257// MdcGeoLayer * layer = mdcGeomSvc->Wire(wireID)->Lyr();
258// if (layer != layer_back){
259// layer_back = layer;
260// MdcGeoSuper * superLayer = mdcGeomSvc->Layer(layerID+1)->Sup();
261// if (superLayer != superLayer_back){
262// superLayer_back = superLayer;
263// // initialize super-layer
264// superLayerID++;
265// new(_superLayer+superLayerID) FTSuperLayer(wireID,
266// Nlocal[superLayerID+1],
267// layerID+1,
268// Nlayer[superLayerID+1],
269// superLayerID);
270// localLayerID=0;
271// }
272// // initialize layer
273// layerID++;
274// double slantWithSymbol = (mdcGeomSvc->Layer(layerID)->Slant())*(mdcGeomSvc->Layer(layerID)->Sup()->Type());
275// // slant in new MdcGeomSvc include symbol
276// new(_layer+layerID) FTLayer(0.1*layer->Radius(), layer->Slant(),
277// 0.1*(+layer->Length()/2), 0.1*(-layer->Length()/2), 0.1*layer->Offset(),
278// layerID, localLayerID++, NlocalWireID[layerID+1],
279// _superLayer[superLayerID]);
280// localID = 0;
281// }
282// // initialize wire
283// const MdcGeoWire * wire = mdcGeomSvc->Wire(wireID);
284// if (superLayerID == 2 || superLayerID == 3 ||
285// superLayerID == 4 || superLayerID == 9 ||
286// superLayerID == 10){ // axial wire
287// new(_wire+wireID) FTWire(0.1*(float)0.5*(wire->Backward().x()+wire->Forward().x()),
288// 0.1*(float)0.5*(wire->Backward().y()+wire->Forward().y()),
289// 0.,0.,
290// _layer[layerID], localID++, _wire+Nwire);
291//
292// }else{ // stereo wire
293// new(_wire+wireID) FTWire(0.1*(float)wire->Backward().x(),
294// 0.1*(float)wire->Backward().y(),
295// 0.1*(float)wire->Forward().x() - 0.1*(float)wire->Backward().x(),
296// 0.1*(float)wire->Forward().y() - 0.1*(float)wire->Backward().y(),
297// _layer[layerID], localID++, _wire+Nwire);
298// }
299// }
300//
301// // make virtual wire object for the pointer of boundary's neighbor
302// new(_wire+Nwire) FTWire();
303//
304// for (int i = 0; i < Nwire; i++){
305// (*(_wire+i)).initNeighbor();
306// }
307//
308// _widbase[0] = 0;
309// for (int k = 1; k < 43; k++){
310// _widbase[k] = _widbase[k-1] + (*(_layer+k-1)).NWire();
311// }
312
313 //minPt = (float)param->_minPt;
314 //minDr = (float)param->_minDr;
315 //xtCoEff = param->_xtCoEff;
316}
317
318void
319FTFinder::term()
320{
321 if (param->_doIt) clear();
322 delete &_tracks;
323 delete _linkedSegments;
324 if (!param->_doIt) return;
325 free(_wire);
326 free(_layer);
327 free(_superLayer);
328}
329
330void
331FTFinder::begin_run(){
332
333 eventNo = 0;
334 runNo=0;
335 expflag=0;
336 if (!param->_doIt) return;
337 _ExpNo = 0;
338 _RunNo = 0;
339
340 IMdcGeomSvc* mdcGeomSvc;
341 StatusCode sc = Gaudi::svcLocator()->service("MdcGeomSvc", mdcGeomSvc);
342
343 if (sc == StatusCode::SUCCESS) {
344 //mdcGeomSvc->Dump();
345
346 } else {
347 return ;
348 }
349
350 const int Nwire = mdcGeomSvc->getWireSize();
351 const int Nlyr = mdcGeomSvc->getLayerSize();
352 const int Nsup = mdcGeomSvc->getSuperLayerSize();
353
354 if(!Nwire || !Nlyr){
355 std::cerr << "FTFINDER::GEOCDC not found (please put cdctable before l4)" << std::endl;
356 std::cerr << "JOB will stop" << std::endl;
357 exit(-1);
358 }
359
360 if (!_wire) _wire = (FTWire *) malloc((Nwire+1) * sizeof(FTWire));
361 if (!_layer) _layer = (FTLayer *) malloc(Nlyr * sizeof(FTLayer));
362 if (!_superLayer) _superLayer =
363 (FTSuperLayer *) malloc(Nsup * sizeof(FTSuperLayer));
364
365 if (!_wire || !_layer || !_superLayer){
366 std::cerr << "FTFINDER::Cannot allocate geometries" << std::endl;
367 std::cerr << "JOB will stop" << std::endl;
368 exit(-1);
369 }
370
371 int superLayerID = 0;
372 int layerID = 0;
373 int localLayerID = 0;
374 int localWireID = 0;
375 int localID = 0;
376 int wireID;
377 MdcGeoLayer * layer_back = NULL;
378 MdcGeoSuper * superLayer_back = NULL;
379 int k = 0;
380 int Nlayer[12];
381 int Nlocal[12];
382 int NlocalWireID[44];
383
384 for (wireID = 0;wireID <= Nwire; wireID++){
385 MdcGeoLayer * layer = (wireID==Nwire) ? NULL : mdcGeomSvc->Wire(wireID)->Lyr();
386 if (layer != layer_back){
387 layer_back = layer;
388 MdcGeoSuper * superLayer = (wireID==Nwire) ? NULL : mdcGeomSvc->Layer(layerID)->Sup();
389 if (superLayer != superLayer_back){
390 superLayer_back = superLayer;
391 Nlayer[k] = localLayerID;
392 Nlocal[k] = localID;
393 localLayerID = 0;
394 k++;
395 }
396 NlocalWireID[layerID] = localWireID;
397 localID = 0;
398 localWireID = 0;
399 layerID++;
400 localLayerID++;
401 }
402 localID++;
403 localWireID++;
404 }
405
406
407 superLayerID = -1;
408 layerID = -1;
409 localLayerID = 0;
410 localID = 0;
411 layer_back = NULL;
412 superLayer_back = NULL;
413 for (wireID = 0;wireID < Nwire; wireID++){
414 MdcGeoLayer * layer = mdcGeomSvc->Wire(wireID)->Lyr();
415 if (layer != layer_back){
416 layer_back = layer;
417 MdcGeoSuper * superLayer = mdcGeomSvc->Layer(layerID+1)->Sup();
418 if (superLayer != superLayer_back){
419 superLayer_back = superLayer;
420 // initialize super-layer
421 superLayerID++;
422 new(_superLayer+superLayerID) FTSuperLayer(wireID,
423 Nlocal[superLayerID+1],
424 layerID+1,
425 Nlayer[superLayerID+1],
426 superLayerID);
427 localLayerID=0;
428 }
429 // initialize layer
430 layerID++;
431 double slantWithSymbol = (mdcGeomSvc->Layer(layerID)->Slant())*(mdcGeomSvc->Layer(layerID)->Sup()->Type());
432 // slant in new MdcGeomSvc include symbol
433 new(_layer+layerID) FTLayer(0.1*layer->Radius(), layer->Slant(),
434 0.1*(+layer->Length()/2), 0.1*(-layer->Length()/2), 0.1*layer->Offset(),
435 layerID, localLayerID++, NlocalWireID[layerID+1],
436 _superLayer[superLayerID]);
437 localID = 0;
438 }
439 // initialize wire
440 const MdcGeoWire * wire = mdcGeomSvc->Wire(wireID);
441 if (superLayerID == 2 || superLayerID == 3 ||
442 superLayerID == 4 || superLayerID == 9 ||
443 superLayerID == 10){ // axial wire
444 new(_wire+wireID) FTWire(0.1*(float)0.5*(wire->Backward().x()+wire->Forward().x()),
445 0.1*(float)0.5*(wire->Backward().y()+wire->Forward().y()),
446 0.,0.,
447 _layer[layerID], localID++, _wire+Nwire);
448
449 }else{ // stereo wire
450 new(_wire+wireID) FTWire(0.1*(float)wire->Backward().x(),
451 0.1*(float)wire->Backward().y(),
452 0.1*(float)wire->Forward().x() - 0.1*(float)wire->Backward().x(),
453 0.1*(float)wire->Forward().y() - 0.1*(float)wire->Backward().y(),
454 _layer[layerID], localID++, _wire+Nwire);
455 }
456 }
457
458 // make virtual wire object for the pointer of boundary's neighbor
459 new(_wire+Nwire) FTWire();
460
461 for (int i = 0; i < Nwire; i++){
462 (*(_wire+i)).initNeighbor();
463 }
464
465 _widbase[0] = 0;
466 for (int k = 1; k < 43; k++){
467 _widbase[k] = _widbase[k-1] + (*(_layer+k-1)).NWire();
468 }
469}
470
471void
472FTFinder::event() {
473 IMessageSvc *msgSvc;
474 Gaudi::svcLocator()->service("MessageSvc", msgSvc);
475
476 MsgStream log(msgSvc, "FTFinder");
477
478 IDataProviderSvc* eventSvc = NULL;
479 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
480
481 if (!param->_doIt) return;
482 //--
483 // clear old information
484 //--
485 clear();
486
487 //check whether Recon already registered
488 DataObject *aReconEvent;
489 eventSvc->findObject("/Event/Recon",aReconEvent);
490 if(!aReconEvent) {
491 // register ReconEvent Data Object to TDS;
492 ReconEvent* recevt = new ReconEvent;
493 StatusCode sc = eventSvc->registerObject("/Event/Recon",recevt );
494 if(sc!=StatusCode::SUCCESS) {
495 log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
496 return;
497 }
498 }
499 //register Event start time
500 IDataManagerSvc *dataManSvc = dynamic_cast<IDataManagerSvc*> (eventSvc);
501 DataObject *aEsTimeEvent;
502 eventSvc->findObject("/Event/Recon/RecEsTimeCol",aEsTimeEvent);
503 if(aEsTimeEvent != NULL) {
504 dataManSvc->clearSubTree("/Event/Recon/RecEsTimeCol");
505 eventSvc->unregisterObject("/Event/Recon/RecEsTimeCol");
506 }
507
508 RecEsTimeCol *aRecEsTimeCol = new RecEsTimeCol;
509 StatusCode est;
510 est = eventSvc->registerObject("/Event/Recon/RecEsTimeCol",aRecEsTimeCol);
511 if( est.isFailure() ) {
512 log << MSG::FATAL << "Could not register RecEsTimeCol" << endreq;
513 return;
514 }
515 log << MSG::DEBUG << "RecEsTimeCol registered successfully!" <<endreq;
516 //--
517 // update wirehit information
518 //--
519 updateMdc();
520
521 //--
522 // segment finding
523 //--
524 for (int i = 0; i^11; i++){
525 (*(_superLayer+i)).mkSegmentList();
526 }
527
528 //--
529 // reduce noise and get start time from segment fit
530 for(int j=0;j<10;j++){
531 (*(_superLayer+j)).reduce_noise(tEstime);
532 }
533
534 getTestime();
535
536 //--
537 // register t0Estime to TDS
538 //--
539 //if(t0Estime!=-999&&t0Estime<24) registT0();
540
541 //--
542 // axial segment linking
543 //--
544 mkTrackList();
545
546 //--
547 // 2D track fitting
548 //--
549 //(*_superLayer).reAppendSalvage();
550
551 int n = _tracks.length();
552 log << MSG::DEBUG << "number of 2D tracks: " << n <<endreq;
553#ifndef OnlineMode
554 num_2Dtrk+=n;
555#endif
556 for (int i = 0; i^n; i++){
557 if (!_tracks[i]->r_phiFit()){
558 delete _tracks[i];
559 log<<MSG::DEBUG<<"===========>deleted 2D track : "<< i+1 <<endreq;
560 n = _tracks.remove(i);
561 }
562 }
563
564 if(t0Estime!=-999){
565 //begin to make Event start time
566 RecEsTime *arecestime = new RecEsTime;
567 arecestime -> setTest(t0Estime);
568 arecestime -> setStat(tEstFlag);
569 aRecEsTimeCol->push_back(arecestime);
570 }
571 if (!n){
572 makeTds();
573 return;
574 }
575
576 ///--
577 /// vertex(r-phi) fit and event timing correction
578 //--
579 int vtx_flag = VertexFit(0);
580 evtTiming = (param->_evtTimeCorr) ? CorrectEvtTiming() : 0;
581
582 ///--
583 /// 2D track re-fitting
584 //--
585 if(param->_hitscut==1){
586 for (int i = 0; i^n; i++){
587 _tracks[i]->r_phiReFit(_vx, _vy, vtx_flag);
588 }
589 }
590
591 //--
592 //cut bad hits from 2D track re-fitting
593 if(param->_hitscut==2){
594 for (int i = 0; i^n; i++){
595 for(int j = 0; j<2; j++){
596 i_rPhiFit= _tracks[i]->r_phi2Fit(_vx, _vy, vtx_flag);
597 if(i_rPhiFit!=-99) _tracks[i]->r_phi3Fit(i_rPhiFit,_vx, _vy, vtx_flag);
598 }
599 _tracks[i]->r_phi4Fit(_vx, _vy, vtx_flag);
600 }
601 }
602
603 //--
604 // stereo segment linking
605 //--
606 mkTrack3D();
607
608 //--
609 // final track fittng
610 //--
611 n = _tracks.length();
612 log<< MSG::DEBUG <<"number of 3D tracks: " << n << endreq;
613
614#ifndef OnlineMode
615 num_3Dtrk+=n;
616#endif
617
618 for (int i = 0; i^n; i++) {
619
620#ifndef OnlineMode
621 log<<MSG::DEBUG<<"=======>3D track: "<<i<<endreq;
622 _tracks[i]->printout();
623#endif
624
625 if (_tracks[i]->get_nhits() < 18) {
626 delete _tracks[i];
627 log<< MSG::DEBUG <<"================>deleted 3D track : "<< i+1 <<endreq;
628 n = _tracks.remove(i);
629 }
630 }
631
632 if (!n){
633 makeTds();
634 return;
635 }
636
637 for (int i = 0; i^n; i++){
638 _tracks[i]->s_zFit();
639 }
640
641 //--
642 // find primary event vertex
643 //--
644 if (param->_findEventVertex){
645 VertexFit(1);
646 }
647
648 n = _tracks.length();
649 log<<MSG::DEBUG<<"final number of tracks: " << n << endreq;
650
651#ifndef OnlineMode
652 num_finaltrk+=n;
653 if (param->_mkMdst) makeMdst();
654#endif
655
656 //--
657 // output tracking result into TDS
658 // by wangdy
659 //--
660 if (param->_mkTds) makeTds();
661
662}
663
664int
665FTFinder::updateMdc(void){
666 IMessageSvc *msgSvc;
667 Gaudi::svcLocator()->service("MessageSvc", msgSvc);
668
669 MsgStream log(msgSvc, "FTFinder");
670
671 IDataProviderSvc* eventSvc = NULL;
672 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
673
674
675 if (eventSvc) {
676 //log << MSG::INFO << "Could find event Svc" << endreq;
677 } else {
678 log << MSG::FATAL << "Could not find Event Header" << endreq;
679 return 0;
680 }
681
682
683 // Part 1: Get the event header, print out event and run number
684
685 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc,"/Event/EventHeader");
686 if (!eventHeader) {
687 log << MSG::FATAL << "Could not find Event Header" << endreq;
688 return( StatusCode::FAILURE);
689 }
690 log << MSG::DEBUG << "MdcFastTrkAlg: retrieved event: " << eventHeader->eventNumber() << " run: " << eventHeader->runNumber() << endreq;
691 eventNo = eventHeader->eventNumber();
692 runNo = eventHeader->runNumber();
693 if(runNo>0) expflag=1;
694 else expflag=0;
695 int digiId;
696
697#ifndef OnlineMode
698 g_eventNo = eventHeader->eventNumber();
699
700 // Part 2: Retrieve MC truth
701 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc,"/Event/MC/McParticleCol");
702 if (!mcParticleCol) {
703 log << MSG::WARNING << "Could not find McParticle" << endreq;
704 // return( StatusCode::FAILURE);
705 }
706 else{
707 McParticleCol ::iterator iter_mc = mcParticleCol->begin();
708 digiId = 0;
709 int ntrkMC = 0;
710 int charge = 0;
711 for (;iter_mc != mcParticleCol->end(); iter_mc++, digiId++) {
712 log << MSG::DEBUG << "MDC digit No: " << digiId << endreq;
713 log << MSG::DEBUG
714 << " particleId = " << (*iter_mc)->particleProperty ()
715 << endreq;
716 int statusFlags = (*iter_mc)->statusFlags();
717 int pid = (*iter_mc)->particleProperty();
718 if(pid >0) {
719 if(m_particleTable->particle( pid )){
720 charge = (int)m_particleTable->particle( pid )->charge();
721 }
722 } else if ( pid <0 ) {
723 if(m_particleTable->particle( -pid )){
724 charge = (int)m_particleTable->particle( -pid )->charge();
725 charge *= -1;
726 }
727 } else {
728 log << MSG::WARNING << "wrong particle id, please check data" <<endreq;
729 }
730
731 if(charge==0 || abs(cos((*iter_mc)->initialFourMomentum().theta()))>0.93) continue;
732
733 g_theta0MC[ntrkMC] = (*iter_mc)->initialFourMomentum().theta();
734 g_phi0MC[ntrkMC] = (*iter_mc)->initialFourMomentum().phi();
735 g_pxMC[ntrkMC] = (*iter_mc)->initialFourMomentum().px();
736 g_pyMC[ntrkMC] = (*iter_mc)->initialFourMomentum().py();
737 g_pzMC[ntrkMC] = (*iter_mc)->initialFourMomentum().pz();
738 g_ptMC[ntrkMC] = 0.001 * sqrt(g_pxMC[ntrkMC]*g_pxMC[ntrkMC] + g_pyMC[ntrkMC]*g_pyMC[ntrkMC]);
739 ntrkMC++;
740 }
741 g_ntrkMC = ntrkMC;
742 }
743#endif
744
745#ifdef MultiThread
746 //Part 3: Retrieve MDC digi
747 SmartDataPtr<MdcDigiCol> mdcDigiVec(eventSvc,"/Event/Digi/MdcDigiCol");
748 if (!mdcDigiVec) {
749 log << MSG::FATAL << "Could not find MDC digi" << endreq;
750 return( StatusCode::FAILURE);
751 }
752 MdcDigiCol::iterator iter1 = mdcDigiVec->begin();
753 digiId = 0;
754 Identifier mdcId;
755 int layerId;
756 int wireId;
757 for (;iter1 != mdcDigiVec->end(); iter1++, digiId++) {
758#else
759 //--use RawDataProvider to get MdcDigi
760 // bool m_dropRecHits=true;//or false
761 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec();
762 MdcDigiVec::iterator iter1 = mdcDigiVec.begin();
763 digiId = 0;
764 Identifier mdcId;
765 int layerId;
766 int wireId;
767 for (;iter1 != mdcDigiVec.end(); iter1++, digiId++) {
768#endif
769 mdcId = (*iter1)->identify();
770 layerId = MdcID::layer(mdcId);
771 wireId = MdcID::wire(mdcId);
772#ifndef OnlineMode
773 log << MSG::DEBUG << "MDC digit No: " << digiId << endreq;
774 log << MSG::DEBUG
775 << " time_channel = " << RawDataUtil::MdcTime((*iter1)->getTimeChannel())
776 << " charge_channel = " << RawDataUtil::MdcCharge((*iter1)->getChargeChannel())
777 << " layerId = " << layerId
778 << " wireId = " << wireId
779 << endreq;
780#endif
781
782 const int localwid = wireId;
783 const int wid = localwid + _widbase[layerId];
784
785#ifndef OnlineMode
786 g_ncellMC->fill((float)wid, 1.0);
787#endif
788 //... skip invalid wireID's
789 if (wid < 0 || wid > 6795) continue;
790 FTWire & w = *(_wire + wid);
791 float tdc = RawDataUtil::MdcTime((*iter1)->getTimeChannel());
792 float adc = RawDataUtil::MdcCharge((*iter1)->getChargeChannel());
793
794#ifndef OnlineMode
795 if (eventNo == 466) {
796 g_hitmap[0]->fill(w.x(), w.y());
797 }
798
799 if (eventNo == 538) {
800 g_hitmap[1]->fill(w.x(), w.y());
801 }
802
803 if (eventNo == 408) {
804 g_hitmap[2]->fill(w.x(), w.y());
805 }
806
807 if (eventNo == 499) {
808 g_hitmap[3]->fill(w.x(), w.y());
809 }
810
811 if (eventNo == 603) {
812 g_hitmap[4]->fill(w.x(), w.y());
813 }
814
815 if (eventNo == 938) {
816 g_hitmap[5]->fill(w.x(), w.y());
817 }
818#endif
819 // tdc/adc calibration
820 //new tdc include drift time and flight time
821 float time = tdc-sqrt(w.x()*w.x()+w.y()*w.y())/30;
822 w.time(time);
823 w.wireId(wid);
824 w.setAdc(adc);
825
826 //... check adc validation
827 const FTLayer & lyr = w.layer();
828
829 //... save to the array
830 //w.distance(time*40*0.0001); //original
831 w.distance(0.0); //by X.-R. Lu
832 w.setChi2(999);
833 w.state(FTWireHit);
834 FTSuperLayer & sup = (FTSuperLayer &) lyr.superLayer();
835 sup.appendHit(&w);
836 }
837
838 return 1;
839}
840
841
842void
843FTFinder::clear(void){
844 evtTiming = 0;
845
846 if ( _superLayer ) {
847 for (int i = 0; i^11; i++) (*(_superLayer+i)).clear();
848 }
849 _tracks.deleteAll();
850 _axialSegCollect.deleteAll();
851 _vx = -99999.;
852 _vy = -99999.;
853 _vz = -99999.;
854 for (int i = 0; i < 10; i++) {
855 tEstime[i].clear();
856 }
857}
858
859int
860FTFinder::getTestime(void)
861{
862 IMessageSvc *msgSvc;
863 Gaudi::svcLocator()->service("MessageSvc", msgSvc);
864
865 MsgStream log(msgSvc, "FTFinder");
866
867 IDataProviderSvc* eventSvc = NULL;
868 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
869
870 float sumT=0,estime=0;
871 int n=0;
872 t0Estime=-999;
873 // for(int i=0;i<4;i++){
874 for(int i=0;i<10;i++){
875 for(int j=0;j<tEstime[i].length();j++){
876 if(tEstime[i][j]!=0){
877 sumT+=tEstime[i][j];
878 n++;
879 }
880 }
881 }
882 if(n!=0){
883 estime=sumT/n;
884 estime+=_t0cal;
885 if(estime>-1){ // && estime<50){
886 //if(m_nbunch==3){
887 if(expflag==0){
888 int nbunch=((int)(estime-tOffSet))/_bunchtime;
889 if(((int)(estime-tOffSet))%(int)_bunchtime>_bunchtime/2) nbunch=nbunch+1;
890 t0Estime=nbunch*_bunchtime+tOffSet;
891 //cout<<"_bunchtime: "<< _bunchtime<<" , t0Estime : "<< t0Estime<<endl;
892 }
893 else{
894 t0Estime=estime;
895 }
896 //if(m_nbunch==6){
897 /* int nbunch=((int)(estime-tOffSet))/4;
898 if(((int)(estime-tOffSet))%4>2) nbunch=nbunch+1;
899 t0Estime=nbunch*4+tOffSet;
900 */
901 // if(t0Estime>-1 && t0Estime<24) FTTrack::Testime=t0Estime;
902 int trigtimming=0;
903 // Retrieve trigger timming information
904 // timing system: TOF:1, MDC:2, EMC:3, NONE:0
905 SmartDataPtr<TrigData> trigData(eventSvc,"/Event/Trig/TrigData");
906 if (trigData) {
907 trigtimming=trigData->getTimingType();
908 log << MSG::INFO <<"Timing type: "<< trigData->getTimingType() << endreq;
909 }
910 if(trigtimming==1) tEstFlag=117;
911 if(trigtimming==2) tEstFlag=127;
912 if(trigtimming==3) tEstFlag=137;
913 if(trigtimming==0) tEstFlag=107;
914#ifndef OnlineMode
915 g_estime=estime;
916#endif
917 }
918 }
919}
920
921void
922FTFinder::mkTrackList(void)
923{
924
925 IMessageSvc *msgSvc;
926 Gaudi::svcLocator()->service("MessageSvc", msgSvc);
927
928 MsgStream log(msgSvc, "FTFinder");
929
930 FTTrack ** tracks = (FTTrack **)malloc( 6 * sizeof(FTTrack *));
931 int * Nsame = (int *)malloc( 6 * sizeof(int));
932 int ntrack_cand = 0;
933 int i,j,k,l;
934
935 //for (int i =0; i < tracks().lenght(); i++) {
936 // FTTrack* fttrk = tracks()[i];
937 // fttrk->setFTFinder(this);
938 //}
939
940 FTList<FTSegment *> inner_segments(10);
941 FTList<FTSegment *> outer_segments(10);
942
943 linkAxialSuperLayer234(inner_segments);
944 linkAxialSuperLayer910(outer_segments);
945
946 _axialSegCollect.append(inner_segments);
947 _axialSegCollect.append(outer_segments);
948
949 int innerN = inner_segments.length();
950 int outerN = outer_segments.length();
951
952 log << MSG::DEBUG << innerN <<" segments in inner axial layers!"<<endreq;
953 log << MSG::DEBUG << outerN <<" segments in outer axial layers!"<<endreq;
954
955 for (l = 0; l^innerN; l++){
956 if (inner_segments[l]->track()) continue;
957 FTSegment *inner = inner_segments[l];
958 FTTrack *track_cand = NULL;
959
960#ifndef OnlineMode
961 inner->printout();
962#endif
963
964 if(outerN==0){
965 track_cand = linkAxialSegments(&inner, NULL);
966 if(!track_cand) continue;
967 _linkedSegments = new FTList<FTSegment *>(5);
968 // track_cand->printout();
969 // append this to _tracks
970 FTList<FTSegment *>& segments = track_cand->axial_segments();
971 int n=segments.length();
972 for (i = 0; i < n; i++){
973 segments[i]->track(track_cand);
974 }
975 _tracks.append(track_cand);
976 track_cand->setFTFinder(this);
977 ntrack_cand++;
978 continue;
979 }
980
981 for(k = 0; k^outerN; k++) {
982 if (outer_segments[k]->track()) { //already used
983 if(k==outerN-1) { //if the last outer segment
984 if (inner_segments[l]->track()) continue; //inner segment never used
985 track_cand = linkAxialSegments(&inner, NULL);
986 if(!track_cand) continue;
987 _linkedSegments = new FTList<FTSegment *>(5);
988 // append this to _tracks
989 FTList<FTSegment *>& segments = track_cand->axial_segments();
990 int n=segments.length();
991 for (i = 0; i < n; i++){
992 segments[i]->track(track_cand);
993 }
994 _tracks.append(track_cand);
995 track_cand->setFTFinder(this);
996 ntrack_cand++;
997 }
998 continue;
999 }
1000
1001 FTSegment *outer = outer_segments[k];
1002 track_cand = linkAxialSegments(&inner, &outer);
1003
1004 if (!track_cand ){ // link failed
1005 if( k == outerN-1 ) { //if the last outer segment
1006 if (inner_segments[l]->track()) continue; //inner segment never used
1007 track_cand = linkAxialSegments(&inner, NULL);
1008 }
1009 if(!track_cand) continue;
1010 }
1011 // else {
1012 // innerN = inner_segments.remove(l);
1013 // outerN = outer_segments.remove(k);
1014 // }
1015 //track_cand->printout();
1016
1017 ntrack_cand++;
1018 _linkedSegments = new FTList<FTSegment *>(5);
1019
1020 //compair this to appended track candidate
1021
1022 FTList<FTSegment *>& segments = track_cand->axial_segments();
1023 int n = segments.length();
1024
1025 //compare this to the appended track candidate
1026 if(inner_segments[l]->track()) { //already used
1027 FTTrack * trk_tmp = inner_segments[l]->track();
1028 int nTrks = _tracks.length();
1029 int cache_i=0;
1030 for (i = 0; i^nTrks; i++){
1031 if (_tracks[i] == trk_tmp){
1032 cache_i = i;
1033 break;
1034 }
1035 }
1036 FTList<FTSegment *> & segments2 = trk_tmp->axial_segments();
1037 int n2 = segments2.length();
1038 if(n >= n2 && track_cand->chi2_kappa_tmp() < trk_tmp->chi2_kappa_tmp()){
1039 // swap this and appended one
1040 for (i = 0; i < n2; i++){
1041 segments2[i]->track(NULL);
1042 }
1043 for (i = 0; i < n; i++){
1044 segments[i]->track(track_cand);
1045 }
1046 _tracks.replace(cache_i,track_cand);
1047 track_cand->setFTFinder(this);
1048 delete trk_tmp;
1049 }else{
1050 delete track_cand;
1051 }
1052 }else{
1053 // append this
1054 for (i = 0; i < n; i++){
1055 segments[i]->track(track_cand);
1056 }
1057 _tracks.append(track_cand);
1058 track_cand->setFTFinder(this);
1059 }
1060 } // end of loop of outerN
1061 } //end of loop innerN
1062
1063 free(tracks);
1064 free(Nsame);
1065
1066#ifndef OnlineMode
1067 int n=_tracks.length();
1068 for(int i=0; i<n; i++){
1069 FTList<FTSegment *> & segments = _tracks[i]->axial_segments();
1070 log << MSG::DEBUG <<" loop connected axial track: " << i <<endreq;
1071 int l=segments.length();
1072 for(int j=0; j<l; j++){
1073 segments[j]->printout();
1074 }
1075 }
1076#endif
1077
1078}
1079
1080FTTrack *
1081FTFinder::linkAxialSegments(FTSegment **inner, FTSegment ** outer)
1082{
1083 float chi2_kappa = param->_chi2_kappa;
1084 _linkedSegments->clear();
1085 if(!outer) {
1086 // allow only 2, 3, 4, axial segments in one track
1087 int n = (*inner)->wireHits().length();
1088 _linkedSegments->append(*inner);
1089 if(n >= 7) return (new FTTrack(*_linkedSegments, (*inner)->kappa(), chi2_kappa));
1090 else return NULL;
1091 }
1092 int n = _linkedSegments->append(*outer);
1093 float SigmaK = (*outer)->kappa();
1094 float SigmaRR = (*outer)->r();
1095 SigmaRR *= SigmaRR;
1096 float SigmaKRR = SigmaK*SigmaRR;
1097 float SigmaKKRR = SigmaK*SigmaKRR;
1098 FTSegment & s = * (*_linkedSegments)[n-1];
1099 FTSegment * innerSegment = NULL;
1100 float SigmaK_cache, SigmaRR_cache, SigmaKRR_cache, SigmaKKRR_cache;
1101
1102 float Min_chi2 = param->_Min_chi2;
1103 float inX = s.incomingX();
1104 float inY = s.incomingY();
1105 //const FTWire & innerBoundHit = * s.innerBoundHits().first();
1106 //float in_r = innerBoundHit.layer().r();
1107 //float incomingPhi = innerBoundHit.phi();
1108 float in_r = s.innerBoundHits().first()->layer().r();
1109 float incomingPhi = s.incomingPhi();
1110
1111 FTSegment* next = *inner;
1112 const FTWire & NextOuterBoundHit = * next->outerBoundHits().first();
1113
1114 //float deltaPhi =fabs(incomingPhi - next->outgoingPhi());
1115 //if (deltaPhi > param->_deltaPhi && deltaPhi < (2*M_PI-param->_deltaPhi)) return NULL;
1116 float outX = next->outgoingX();
1117 float outY = next->outgoingY();
1118
1119 //////////////////////////////////////////////////////
1120 float _trk_d = -2*(-1. / 2.99792458 /m_pmgnIMF->getReferField())/s.kappa();
1121 float _angle1 = asin(NextOuterBoundHit.layer().r()/_trk_d);
1122 float _angle2 = asin(s.outerBoundHits().first()->layer().r()/_trk_d);
1123 float _ang_diff = _angle2 - _angle1;
1124 float _diff = s.outgoingPhi() - next->outgoingPhi();
1125 _diff = _diff - (int(_diff/M_PI))*2*M_PI;
1126 float _require = _ang_diff - _diff;
1127 //cut of connecting inner and outer segments
1128 if (_require < -0.10 || _require > 0.11) return NULL;
1129
1130 float SegK = next->kappa();
1131 float SegRR = next->r();
1132 SegRR *= SegRR;
1133 const float out_r = NextOuterBoundHit.layer().r();
1134 // kappa = -2. * alpha * ((Vout X Vin)_z / |Vin|*|Vout|) / |Vin-Vout|
1135 float GapK = 2.*(-1. / 2.99792458 /m_pmgnIMF->getReferField())*(inX*outY-inY*outX) /
1136 (in_r*out_r*sqrt((inX-outX)*(inX-outX)+(inY-outY)*(inY-outY)));
1137 //float GapRR = (currentLayer==j+1||currentLayer==j+2) ? 0.5*(in_r+out_r) : in_r+out_r;
1138 float GapRR = 0.5*(in_r+out_r);
1139 GapRR*=GapRR;
1140 float SigmaK_tmp = (SigmaK + SegK + GapK);
1141 float SigmaRR_tmp = SigmaRR + SegRR + GapRR;
1142 float SigmaKRR_tmp = SigmaKRR + SegK*SegRR + GapK*GapRR;
1143 float SigmaKKRR_tmp = SigmaKKRR + SegK*SegK*SegRR + GapK*GapK*GapRR;
1144 float MuK_tmp = SigmaK_tmp/(2*n+1);
1145 float chi2 = (MuK_tmp*MuK_tmp*SigmaRR_tmp
1146 - 2.*MuK_tmp*SigmaKRR_tmp + SigmaKKRR_tmp)/(2*n+1);
1147 if ((chi2-chi2_kappa) < Min_chi2){
1148 Min_chi2 = chi2;
1149 innerSegment = next;
1150 SigmaK_cache = SigmaK_tmp;
1151 SigmaRR_cache = SigmaRR_tmp;
1152 SigmaKRR_cache = SigmaKRR_tmp;
1153 SigmaKKRR_cache = SigmaKKRR_tmp;
1154 }
1155 if (innerSegment){
1156 n = _linkedSegments->append(*inner);
1157 SigmaK = SigmaK_cache;
1158 SigmaRR = SigmaRR_cache;
1159 SigmaKRR = SigmaKRR_cache;
1160 SigmaKKRR = SigmaKKRR_cache;
1161 chi2_kappa = Min_chi2;
1162
1163 if (n > 1) return (new FTTrack(*_linkedSegments,SigmaK/(2*n-1),chi2_kappa));
1164 }else{
1165 // allow only 3, 4, 5 axial segments in one track
1166 /*n = (*inner)->wireHits().length();
1167 _linkedSegments->clear();
1168 _linkedSegments->append(*inner);
1169 if(n > 10) return (new FTTrack(*_linkedSegments, SegK, Min_chi2));*/
1170 return NULL;
1171 }
1172 //if (fabs(SigmaK/(2*n-1)) > 1.2/minPt) return NULL;
1173 return NULL;
1174
1175}
1176
1177void
1178FTFinder::linkAxialSuperLayer234(FTList<FTSegment *> & inner_segments)
1179{
1180 FTList<FTSegment *> _segments34(10);
1181 FTList<FTSegment *>& SuperLayer3Segments = _superLayer[3].segments();
1182 FTList<FTSegment *>& SuperLayer4Segments = _superLayer[4].segments();
1183 linkAxialSegments_step(SuperLayer3Segments, SuperLayer4Segments,
1184 _segments34, param->_D_phi2, param->_chi2_2);
1185 _segments34.append(SuperLayer3Segments);
1186 SuperLayer3Segments.clear();
1187
1188 FTList<FTSegment *>& SuperLayer2Segments = _superLayer[2].segments();
1189 linkAxialSegments_step(SuperLayer2Segments, _segments34,
1190 inner_segments, param->_D_phi1, param->_chi2_1);
1191 inner_segments.append(_segments34);
1192
1193 //zoujh: connect the 2&4 superlayers leap over 3
1194 int n = SuperLayer2Segments.length();
1195 int m = SuperLayer4Segments.length();
1196 for (int i = 0; i^n; i++) {
1197 FTSegment * inner = SuperLayer2Segments[i];
1198 //inner->printout();
1199 float in_outerPhi = inner->outgoingPhi();
1200 for (int j = 0; j^m; j++) {
1201 FTSegment * outer = SuperLayer4Segments[j];
1202 // outer->printout();
1203 float out_innerPhi = outer->incomingPhi();
1204 float D_phi = fabs(in_outerPhi - out_innerPhi);
1205 if (D_phi > M_PI) D_phi = 2*M_PI - D_phi;
1206 /////////////////////////////////////////////////////
1207 if (D_phi > M_PI/12.5) continue;
1208 float inX = inner->outgoingX();
1209 float inY = inner->outgoingY();
1210 float outX = outer->incomingX();
1211 float outY = outer->incomingY();
1212 float in_r = inner->outerBoundHits().first()->layer().r();
1213 float out_r = outer->innerBoundHits().first()->layer().r();
1214 float GapK = 2.*(-1. / 2.99792458 /m_pmgnIMF->getReferField())*(inY*outX-inX*outY) /
1215 (in_r*out_r*sqrt((inX-outX)*(inX-outX)+(inY-outY)*(inY-outY)));
1216 //float cache_in = ((inner->kappa()+outer->kappa()+GapK)/3.0 - inner->kappa())/in_r;
1217 float cache_in = ((outer->kappa()+GapK)/3.0 - inner->kappa()*2.0/3.0)/in_r;
1218 float cache_out = ((inner->kappa()+GapK)/3.0 - outer->kappa()*2.0/3.0)/out_r;
1219 float cache_gap = ((inner->kappa()+outer->kappa())/3.0-GapK*2.0/3.0)*2.0/(in_r+out_r);
1220 float chi2_z = cache_in*cache_in + cache_out*cache_out + cache_gap*cache_gap;
1221 if (chi2_z > param->_D_phi1) continue;
1222 /////////////////////////////////////////////////////
1223 inner->connect_outer(outer->wireHits(),outer->outerBoundHits());
1224 inner->update();
1225 inner_segments.append(inner);
1226 n = SuperLayer2Segments.remove(i);
1227 m = SuperLayer4Segments.remove(j);
1228 delete outer;
1229 break;
1230 }
1231 }
1232 inner_segments.append(SuperLayer2Segments);
1233 SuperLayer2Segments.clear();
1234 inner_segments.append(SuperLayer4Segments);
1235 SuperLayer4Segments.clear();
1236}
1237
1238void
1239FTFinder::linkAxialSuperLayer910(FTList<FTSegment *> & outer_segments)
1240{
1241 FTList<FTSegment *>& SuperLayer9Segments = _superLayer[9].segments();
1242 FTList<FTSegment *>& SuperLayer10Segments = _superLayer[10].segments();
1243 linkAxialSegments_step(SuperLayer9Segments, SuperLayer10Segments,
1244 outer_segments, param->_D_phi3, param->_chi2_3);
1245 outer_segments.append(SuperLayer9Segments);
1246 SuperLayer9Segments.clear();
1247 outer_segments.append(SuperLayer10Segments);
1248 SuperLayer10Segments.clear();
1249}
1250
1251void
1252FTFinder::linkAxialSegments_step(FTList<FTSegment *>& InnerSegments,
1253 FTList<FTSegment *>& OuterSegments,
1254 FTList<FTSegment *>& Connected,
1255 float maxDphi, float maxChi2)
1256{
1257 IMessageSvc *msgSvc;
1258 Gaudi::svcLocator()->service("MessageSvc", msgSvc);
1259 MsgStream log(msgSvc, "FTFinder");
1260
1261 int n = InnerSegments.length();
1262 int m = OuterSegments.length();
1263 for (int i = 0; i^n; i++){
1264 FTSegment * inner = InnerSegments[i];
1265#ifndef OnlineMode
1266 // log<<MSG::DEBUG<<"linking: "<<endreq;
1267 // inner->printout();
1268#endif
1269 const FTLayer & in_outerBound = inner->outerBoundHits().first()->layer();
1270 float in_outerPhi = inner->outgoingPhi();
1271 float min_Dphi = M_PI/2;
1272 int min_Dphi_index = -1;
1273 for (int j = 0; j^m; j++) {
1274 FTSegment * outer = OuterSegments[j];
1275#ifndef OnlineMode
1276 // log<<MSG::DEBUG<<"segment: "<<j<<endreq;
1277 // outer->printout();
1278#endif
1279 float D_phi = fabs(in_outerPhi - outer->incomingPhi());
1280 if (D_phi > M_PI) D_phi = 2*M_PI - D_phi;
1281 if (D_phi > min_Dphi) continue;
1282 float inX = inner->incomingX();
1283 float inY = inner->incomingY();
1284 float outX = outer->outgoingX();
1285 float outY = outer->outgoingY();
1286 float in_r = inner->innerBoundHits().first()->layer().r();
1287 float out_r = outer->outerBoundHits().first()->layer().r();
1288 float allK = 2.*(-1. / 2.99792458 /m_pmgnIMF->getReferField())*(inY*outX-inX*outY) /
1289 (in_r*out_r*sqrt((inX-outX)*(inX-outX)+(inY-outY)*(inY-outY)));
1290 //float cache_in = ((inner->kappa() + outer->kappa() + allK)/3.0 - inner->kappa())/in_r;
1291 float cache_in = ((outer->kappa() + allK)/3.0 - inner->kappa()*2/3.0)/in_r;
1292 float cache_out = ((inner->kappa() + allK)/3.0 - outer->kappa()*2/3.0)/out_r;
1293 float cache_all = ((inner->kappa()+outer->kappa())/3.0-allK*2/3.0)*2.0/(in_r+out_r);
1294 float chi2_z = cache_in*cache_in + cache_out*cache_out + cache_all*cache_all;
1295 log<<MSG::DEBUG<<"D_phi: "<< D_phi <<" chi2_z: "<< chi2_z <<" maxChi2: " <<maxChi2 <<endreq;
1296 if (chi2_z > maxChi2) continue;
1297 min_Dphi = D_phi;
1298 min_Dphi_index = j;
1299 }
1300 if (min_Dphi_index < 0) continue;
1301 log<<MSG::DEBUG<<"min_Dphi: "<< min_Dphi <<endreq;
1302 FTSegment * outer = OuterSegments[min_Dphi_index];
1303 const FTLayer & out_innerBound = outer->innerBoundHits().first()->layer();
1304 switch (out_innerBound.layerId() - in_outerBound.layerId()) {
1305 case 1:
1306 if (min_Dphi > maxDphi) continue;
1307 break;
1308 case 2:
1309 if (min_Dphi > maxDphi*1.5) continue;
1310 break;
1311 case 3:
1312 if (min_Dphi > maxDphi*2.25) continue;
1313 break;
1314 default:
1315 continue;
1316 }
1317 inner->connect_outer(outer->wireHits(), outer->outerBoundHits());
1318 inner->update();
1319 Connected.append(inner);
1320 n = InnerSegments.remove(i);
1321 delete outer;
1322 m = OuterSegments.remove(min_Dphi_index);
1323 log<<MSG::DEBUG<<"DONE!!"<<endreq;
1324 }
1325 return;
1326}
1327
1328void
1329FTFinder::mkTrack3D(void){
1330 IMessageSvc *msgSvc;
1331 Gaudi::svcLocator()->service("MessageSvc", msgSvc);
1332 MsgStream log(msgSvc, "FTFinder");
1333
1334 int n = _tracks.length();
1335 // select segment candidate and cache sz
1336 FTList<FTSegment *> multi_trk_cand(20);
1337 for (int i=8; i>=0 ; i-- ){
1338 //if (i==2 || i==3 || i==4) continue;
1339 if (i == 4) i -= 3;
1340 FTList<FTSegment *> & segments = _superLayer[i].segments();
1341 int m = segments.length();
1342 for (int j = 0; j^m; j++){
1343 FTSegment * s = segments[j];
1344#ifndef OnlineMode
1345 s->printout();
1346#endif
1347 int nTrack = 0;
1348 FTTrack * t;
1349 for (int k = 0; k^n; k++){
1350#ifndef OnlineMode
1351 log<<MSG::DEBUG<<"coupling to track: " << k <<endreq;
1352 //_tracks[k]->printout();
1353#endif
1354 if (s->update3D(_tracks[k])){ // calculate s and z
1355 t = _tracks[k];
1356 nTrack++;
1357 }
1358 }
1359 switch(nTrack){
1360 case 0:
1361 continue;
1362 case 1:
1363 t->append_stereo_cache(s);
1364 break;
1365 default:
1366 multi_trk_cand.append(s);
1367 break;
1368 }
1369 }
1370 // whose relation between this is unique
1371 for (int j = 0; j^n; j++) _tracks[j]->updateSZ();
1372 }
1373 // link segments by tanLambda & dz
1374 for (int i = 0; i^n; i++) _tracks[i]->linkStereoSegments();
1375 n = multi_trk_cand.length();
1376 for (int i = 0; i^n; i++) multi_trk_cand[i]->linkStereoSegments();
1377}
1378
1379int
1380FTFinder::VertexFit2D()
1381{
1382 int nTrks = _tracks.length();
1383 if (nTrks < 2){
1384 _vx = -99999.;
1385 _vy = -99999.;
1386 _vz = -99999.;
1387 return 0;
1388 }
1389 FTList<float> px(10);
1390 FTList<float> py(10);
1391 FTList<float> dx(10);
1392 FTList<float> dy(10);
1393 FTList<double> weight(10);
1394 for (int i = 0; i < nTrks; i++){
1395 const Lpav & la = _tracks[i]->lpav();
1396 HepVector a = la.Hpar(pivot);
1397
1398 const float dr_i = a(1);
1399 if (fabs(a(1)) > 1.5) continue;
1400 const float px_i = - sin(a(2));
1401 const float py_i = cos(a(2));
1402
1403 //const float dx_i = dr_i*py_i;
1404 //const float dy_i = -dr_i*px_i;
1405
1406 float weight_i = la.chisq()/(la.nc()*0.02);
1407
1408 px.append(px_i);
1409 py.append(py_i);
1410 dx.append(dr_i*py_i);
1411 dy.append(-dr_i*px_i);
1412 weight.append(exp(-weight_i*weight_i));
1413 }
1414 if (dx.length() < 2){
1415 _vx = -99999.;
1416 _vy = -99999.;
1417 _vz = -99999.;
1418 return 0;
1419 }
1420 double S_weight = 0.;
1421 for (int i = dx.length() - 1; i; i--){
1422 const float px_i = px[i];
1423 const float py_i = py[i];
1424 const float dx_i = dx[i];
1425 const float dy_i = dy[i];
1426 const float weight_i = weight[i];
1427 for (int j = 0; j < i; j++){
1428 const float px_j = px[j];
1429 const float py_j = py[j];
1430 //const float weight_j = weight[j];
1431
1432 const float ddx = dx[j] - dx_i;
1433 const float ddy = dx[j] - dy_i;
1434
1435 const float tmp_par = px_i*py_j - px_j*py_i;
1436 const float par = (py_j*ddx-px_j*ddy)/tmp_par;
1437
1438 double weight_ij = weight_i*weight[j];
1439 S_weight += weight_i*weight[j];
1440 if (tmp_par==0 ||
1441 par < -0.5 || (py_i*ddx-px_i*ddy)/tmp_par < -0.5 ||
1442 fabs((px_i*px_j + py_i*py_j)/
1443 sqrt((px_i*px_i+py_i*py_i)*(px_j*px_j+py_j*py_j)))>0.86){
1444 _vx += (dx_i+0.5*ddx)*weight_ij;
1445 _vy += (dy_i+0.5*ddy)*weight_ij;
1446 }else{
1447 _vx += (dx_i+par*px_i)*weight_ij;
1448 _vy += (dy_i+par*py_i)*weight_ij;
1449 }
1450 }
1451 }
1452
1453 int rtn_flag = 0;
1454 if (S_weight == 0.){
1455 _vx = -99999.;
1456 _vy = -99999.;
1457 _vz = -99999.;
1458 }else{
1459 _vx /= S_weight;
1460 _vy /= S_weight;
1461 _vz = -99999.;
1462 rtn_flag = 1;
1463 }
1464 return rtn_flag;
1465
1466}
1467
1468int
1469FTFinder::VertexFit(int z_flag)
1470{
1471 _vx = 0.;
1472 _vy = 0.;
1473 _vz = 0.;
1474 if (!z_flag){
1475 return VertexFit2D();
1476 }
1477 int n = _tracks.length();
1478 if (n < 2){
1479 _vx = -99999.;
1480 _vy = -99999.;
1481 _vz = -99999.;
1482 return 0;
1483 }
1484 FTList<float> px(10);
1485 FTList<float> py(10);
1486 FTList<float> pz(10);
1487 FTList<float> dx(10);
1488 FTList<float> dy(10);
1489 FTList<float> dz(10);
1490 FTList<float> pmag2(10);
1491 FTList<float> weight(10);
1492 FTList<float> weight_z(10);
1493 FTList<float> sigma2_r(10);
1494 FTList<float> sigma2_z(10);
1495
1496 for (int i = 0; i < n; i++){
1497 const Lpav & la = _tracks[i]->lpav();
1498 if (la.nc() <= 3) continue;
1499 const zav & za = _tracks[i]->Zav();
1500 if(za.nc() > 2 && za.b() > -900){
1501 pmag2.append(1+za.b()*za.b());
1502 pz.append(za.a());
1503 dz.append(za.b());
1504 sigma2_z.append(za.chisq()/(za.nc()-2));
1505 weight_z.append(exp(-(za.chisq()/(za.nc()-2))));
1506 }else{
1507 continue;
1508 }
1509
1510 HepVector a = la.Hpar(pivot);
1511
1512 const float dr_i = a(1);
1513 const float px_i = - std::sin(a(2));
1514 const float py_i = std::cos(a(2));
1515
1516 px.append(px_i);
1517 py.append(py_i);
1518 dx.append(dr_i*py_i);
1519 dy.append(-dr_i*px_i);
1520 sigma2_r.append(std::sqrt(la.chisq())/(la.nc()-3));
1521 weight.append(exp(-(std::sqrt(la.chisq())/(la.nc()-3))));
1522 }
1523
1524 n = dx.length();
1525 if (n < 2){
1526 _vx = -9999.;
1527 _vy = -9999.;
1528 _vz = -9999.;
1529 return 0;
1530 }
1531
1532 FTList<float> vtx_x(20);
1533 FTList<float> vtx_y(20);
1534 FTList<float> vtx_z(20);
1535 FTList<float> weight2(20);
1536 FTList<float> weight2_z(20);
1537 FTList<float> vtx_chi2(20);
1538 unsigned n_vtx(0);
1539 for (int i = n - 1; i; i--){
1540 for (int j = 0; j < i; j++){
1541 // min. chi2 fit w/ line approximation
1542 const float pij_x = py[i]*pz[j]-pz[i]*py[j];
1543 const float pij_y = pz[i]*px[j]-px[i]*pz[j];
1544 const float pij_z = px[i]*py[j]-py[i]*px[j];
1545
1546 if (pij_x == 0.0f && pij_y == 0.0f && pij_z == 0.0f ) continue;
1547
1548 const float sr = sigma2_r[i]+sigma2_r[j];
1549 const float sz = sigma2_z[i]+sigma2_z[j];
1550
1551 const float ddx = dx[i]-dx[j];
1552 const float ddy = dy[i]-dy[j];
1553 const float ddz = dz[i]-dz[j];
1554
1555 const float pij_mag2 = pij_x*pij_x/(sr*sz)+pij_y*pij_y/(sr*sz)+pij_z*pij_z/(sr*sr);
1556
1557 const float d_i = (pij_x*(py[j]*ddz-pz[j]*ddy)/(sr*sz)+
1558 pij_y*(pz[j]*ddx-px[j]*ddz)/(sr*sz)+
1559 pij_z*(px[j]*ddy-py[j]*ddx)/(sr*sr))/pij_mag2;
1560
1561 const float d_j = (pij_x*(py[i]*ddz-pz[i]*ddy)/(sr*sz)+
1562 pij_y*(pz[i]*ddx-px[i]*ddz)/(sr*sz)+
1563 pij_z*(px[i]*ddy-py[i]*ddx)/(sr*sr))/pij_mag2;
1564
1565 const float vtx_x_i = dx[i] + px[i]*d_i;
1566 const float vtx_y_i = dy[i] + py[i]*d_i;
1567 const float vtx_z_i = dz[i] + pz[i]*d_i;
1568
1569 const float vtx_x_j = dx[j] + px[j]*d_j;
1570 const float vtx_y_j = dy[j] + py[j]*d_j;
1571 const float vtx_z_j = dz[j] + pz[j]*d_j;
1572
1573 const float weight_ij = weight[i]+weight[j];
1574 vtx_x.append((weight[j]*vtx_x_i + weight[i]*vtx_x_j)/weight_ij);
1575 vtx_y.append((weight[j]*vtx_y_i + weight[i]*vtx_y_j)/weight_ij);
1576 vtx_z.append((weight_z[j]*vtx_z_i + weight_z[i]*vtx_z_j)/(weight_z[i]+weight_z[j]));
1577
1578 weight2.append(exp(-sr));
1579 weight2_z.append(exp(-sz));
1580 vtx_chi2.append(((vtx_x_i-vtx_x_j)*(vtx_x_i-vtx_x_j)+(vtx_y_i-vtx_y_j)*(vtx_y_i-vtx_y_j))/sr +
1581 (vtx_z_i-vtx_z_j)*(vtx_z_i-vtx_z_j)/sz);
1582 n_vtx++;
1583 }
1584 }
1585 n = vtx_chi2.length();
1586 float S_weight(0.0f);
1587 float S_weight_z(0.0f);
1588 for (int i = 0; i != n; i++){
1589 if (vtx_chi2[i] > 10.) continue;
1590 float w(std::exp(-vtx_chi2[i]));
1591 _vx += vtx_x[i]*weight2[i]*w;
1592 _vy += vtx_y[i]*weight2[i]*w;
1593 _vz += vtx_z[i]*weight2_z[i]*w;
1594 S_weight += weight2[i]*w;
1595 S_weight_z += weight2_z[i]*w;
1596 }
1597 int rtn_flag = 0;
1598 if (S_weight <= 0.){
1599 _vx = -9999.;
1600 _vy = -9999.;
1601 }else{
1602 _vx /= S_weight;
1603 _vy /= S_weight;
1604 rtn_flag = 1;
1605 }
1606 if (!z_flag) return rtn_flag;
1607 if (S_weight_z <= 0.){
1608 _vz = -9999.;
1609 }else{
1610 _vz /= S_weight_z;
1611 rtn_flag++;
1612 }
1613 return rtn_flag;
1614}
1615
1616int
1617FTFinder::findBestVertex(void)
1618{
1619 int nTrks = _tracks.length();
1620 if (nTrks < 2){
1621 _vx = -9999.;
1622 _vy = -9999.;
1623 _vz = -9999.;
1624 return 0;
1625 }
1626 float min_dr = 9999.;
1627 float phi0 = 9999.;
1628 for (int i = 0; i < nTrks; i++){
1629 HepVector a = _tracks[i]->lpav().Hpar(pivot);
1630 if (fabs(a(1)) < fabs(min_dr)){
1631 min_dr = a(1);
1632 phi0 = a(2);
1633 }
1634 }
1635 _vx = min_dr*cos(phi0);
1636 _vy = min_dr*sin(phi0);
1637 return 1;
1638}
1639
1640Hep3Vector
1641FTFinder::vertex(void) const{
1642 return Hep3Vector(_vx,_vy,_vz);
1643}
1644
1645int
1646FTFinder::CorrectEvtTiming(void)
1647{
1648 int nTrks = _tracks.length();
1649 float weight_sum = 0.;
1650 float dt_sum2 = 0.;
1651 for (int i = 0; i^nTrks; i++){
1652 float dt_sum = 0.;
1653 float dtt_sum = 0.;
1654 int nHits = 0;
1655 const Lpav & la = _tracks[i]->lpav();
1656 FTList<FTSegment *>& axial_sgmnts = _tracks[i]->axial_segments();
1657 int m = axial_sgmnts.length();
1658 for (int j = 0; j^m; j++){
1659 FTList<FTWire *>& hits = axial_sgmnts[j]->wireHits();
1660 int l = hits.length();
1661 for (int k = 0; k^l; k++){
1662 FTWire & h = * hits[k];
1663 const float x = h.x();
1664 const float y = h.y();
1665 double d0 = fabs(la.d((double)x,(double)y));
1666 if (d0 >= 0.47f*h.layer().csize()) continue;
1667 nHits++;
1668 float dt = x2t(h.layer(), d0) - h.time();
1669 // float dt = d0/(40*0.0001) - h.time();
1670 dt_sum += dt;
1671 dtt_sum += (dt*dt);
1672 }
1673 }
1674 if (!nHits) continue;
1675 float weight_t = exp(-(dtt_sum - (dt_sum*dt_sum/nHits))/(nHits*1600));
1676 weight_sum += (nHits*weight_t);
1677 dt_sum2 += (dt_sum*weight_t);
1678 }
1679 // cout<<"event correction time: "<< int(dt_sum2/weight_sum)<<endl;
1680 return int(dt_sum2/weight_sum);
1681
1682}
1683
1684#ifndef OnlineMode
1685void
1686FTFinder::makeMdst(void)
1687{
1688 // static const double pi_mass(0.13956995);
1689 const int nTrks = _tracks.length();
1690 int Ntable(0);
1691
1692 for (int i = 0; i<nTrks; i++){
1693
1694 const FTTrack & trk = * _tracks[i];
1695 FTList<FTSegment *> &axialSegments = trk.axial_segments();
1696 FTList<FTSegment *> &stereoSegments = trk.stereo_segments();
1697 int axialSegmentsN = axialSegments.length();
1698 int stereoSegmentsN = stereoSegments.length();
1699 for (int i = 0; i < axialSegmentsN ; i++) {
1700 //FTSuperLayer& superLayer= axialSegments[i]->superLayer();
1701 FTList<FTWire *> &wires = axialSegments[i]->wireHits();
1702 int wiresN = wires.length();
1703 for (int j=0; j < wiresN; j++) {
1704 //int id = wires[j]->localId();
1705 g_ncell->fill(wires[j]->getWireId());
1706 }
1707 }
1708
1709 for (int i = 0; i < stereoSegmentsN ; i++) {
1710 //FTSuperLayer& superLayer = stereoSegments[i]->superLayer();
1711 FTList<FTWire *> &wires = stereoSegments[i]->wireHits();
1712 int wiresN = wires.length();
1713 for (int j=0; j < wiresN; j++) {
1714 //int id = wires[j]->localId();
1715 g_ncell->fill(wires[j]->getWireId());
1716 }
1717 }
1718
1719 const Helix& h = * trk.helix();
1720 if (h.tanl() < -9000.) continue;
1721
1722 Ntable++;
1723 Hep3Vector p(h.momentum());
1724 HepPoint3D v(h.x());
1725 //float m_charge = (h.kappa() > 0) ? 1. : -1.;
1726 float m_px = p.x();
1727 g_px[Ntable-1] = p.x();
1728 float m_py = p.y();
1729 g_py[Ntable-1] = p.y();
1730 float m_pz = p.z();
1731 g_pz[Ntable-1] = p.z();
1732 g_p[Ntable-1] = p.mag();
1733 g_phi[Ntable-1] = atan2(m_py, m_px);
1734
1735 g_dr[Ntable-1] = h.dr();
1736 g_phi0[Ntable-1] = h.phi0();
1737 g_kappa[Ntable-1] = h.kappa();
1738 g_pt[Ntable-1] = 1/fabs(h.kappa());
1739 g_dz[Ntable-1] = h.dz();
1740 g_tanl[Ntable-1] = h.tanl();
1741 g_theta[Ntable-1] = atan2(1/fabs(h.kappa()),(double)(m_pz));
1742 g_vx[Ntable-1] = v(0);
1743 g_vy[Ntable-1] = v(1);
1744 g_vz[Ntable-1] = v(2);
1745 }
1746 g_ntrk = Ntable;
1747
1748}
1749#endif
1750
1751
1752
1753//-------------------------------
1754// output tracking results to TDS by wangdy
1755//-------------------------------
1756
1757StatusCode
1758FTFinder::makeTds(void)
1759{
1760 IMessageSvc *msgSvc;
1761 Gaudi::svcLocator()->service("MessageSvc", msgSvc);
1762 MsgStream log(msgSvc, "FTFinder");
1763
1764 IDataProviderSvc* eventSvc = NULL;
1765 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
1766
1767 if (eventSvc) {
1768 //log << MSG::DEBUG << "makeTds:event Svc has been found" << endreq;
1769 } else {
1770 log << MSG::FATAL << "makeTds:Could not find eventSvc" << endreq;
1771 return StatusCode::FAILURE ;
1772 }
1773
1774 RecMdcTrackCol* trkcol = new RecMdcTrackCol;
1775 RecMdcHitCol* hitcol = new RecMdcHitCol;
1776
1777 //make RecMdcTrackCol
1778#ifndef OnlineMode
1779 log << MSG::DEBUG << "beginning to make RecMdcTrackCol" <<endreq;
1780#endif
1781 int ntrk = tracks().length();
1782 int trkid=0;
1783 for (int i =0; i < ntrk; i++) {
1784 RecMdcTrack* trk = new RecMdcTrack;
1785 FTTrack* fttrk = tracks()[i];
1786
1787#ifndef OnlineMode
1788 // fttrk->printout();
1789#endif
1790
1791 if (fttrk->helix()->tanl() < -9000.){
1792 log << MSG::DEBUG << "deleted trackId = " << i+1 << " due to tanl = " << fttrk->helix()->tanl() << endreq;
1793 delete trk;
1794 continue;
1795 }
1796 // int trackindex = i;
1797 HepPoint3D pos = fttrk->helix()->pivot();
1798 int charge = -1;
1799 HepVector m_a(5,0);
1800 m_a = fttrk->helix()->a();
1801 m_a[2]=-m_a[2];
1802 HepSymMatrix m_ea = fttrk->helix()->Ea();
1803 float fiterm = fttrk->lpav().phi(77.0);
1804 float chi2lpav = fttrk->lpav().chisq();
1805 float chi2zav = fttrk->Zav().chisq();
1806
1807// note: this alg itself can not provide error matrix,so m_ea above is empty;
1808// the following error matrix values are based on BESII results
1809// the values of online Bhabha events are multiplied by factor 10
1810// wangdy 20050418
1811 m_ea[0][0] = 0.0085;
1812 m_ea[1][1] = 0.000011;
1813 m_ea[2][2] = 0.0018;
1814 m_ea[3][3] = 1.2;
1815 m_ea[4][4] = 0.00026;
1816 m_ea[1][0] = m_ea[0][1] = -0.00029;
1817 m_ea[2][0] = m_ea[0][2] = charge*0.004;
1818 m_ea[2][1] = m_ea[1][2] = charge*0.00012;
1819 m_ea[3][0] = m_ea[0][3] = -0.017;
1820 m_ea[3][1] = m_ea[1][3] = 0.0;
1821 m_ea[3][2] = m_ea[2][3] = 0.0;
1822 m_ea[4][0] = m_ea[0][4] = 0.0;
1823 m_ea[4][1] = m_ea[1][4] = 0.0;
1824 m_ea[4][2] = m_ea[2][4] = 0.0;
1825 m_ea[4][3] = m_ea[3][4] = -0.018;
1826
1827 trk->setTrackId(trkid);
1828 trk->setHelix(m_a);
1829 trk->setPivot(pos);
1830 trk->setError(m_ea);
1831 trk->setFiTerm(fiterm);
1832 trk->setChi2(chi2lpav+chi2zav);
1833#ifndef OnlineMode
1834 log <<MSG::DEBUG << " trackId = " << i << endreq;
1835 log <<MSG::DEBUG <<"fast-tracking kappa "<<m_a[2]
1836 <<" fast-tracking tanl "<<m_a[4]
1837 <<endreq;
1838 log <<MSG::DEBUG <<"push_backed kappa "<<trk->helix(2)
1839 <<" push_backed tanl "<< trk->helix(4)
1840 <<endreq;
1841
1842 log << MSG::DEBUG << "beginning to make hitlist and RecMdcHitCol " <<endreq;
1843#endif
1844
1845 HitRefVec hitrefvec;
1846
1847 int hitindex = 0;
1848
1849 //axial segments part
1850 FTList<FTSegment *>& seglist_ax = fttrk->axial_segments();
1851 int nseg_ax = seglist_ax.length();
1852 int ntrackhits = 0;
1853 for (int iseg_ax = 0; iseg_ax < nseg_ax; iseg_ax++) {
1854 FTList<FTWire *>& hitlist_ax = seglist_ax[iseg_ax]->wireHits();
1855 int nhitlist_ax = hitlist_ax.length();
1856 ntrackhits += nhitlist_ax;
1857 for( int ihit_ax = 0; ihit_ax < nhitlist_ax; ihit_ax++,hitindex++) {
1858 FTWire wire_ax = *(hitlist_ax[ihit_ax]);
1859 double dd = wire_ax.distance();
1860 double adc = RawDataUtil::MdcChargeChannel(wire_ax.getAdc());
1861 double tdc = wire_ax.time();
1862 double chi2 = wire_ax.getChi2();
1863 const Identifier mdcid = MdcID::wire_id(wire_ax.layer().layerId(),wire_ax.localId());
1864
1865 RecMdcHit* hit = new RecMdcHit;
1866 hit->setId(hitindex);
1867 hit->setDriftDistLeft( dd );
1868 hit->setDriftDistRight( dd );
1869 hit->setAdc( adc );
1870 hit->setTdc( tdc );
1871 hit->setMdcId( mdcid );
1872 hit->setChisqAdd( chi2 );
1873#ifndef OnlineMode
1874 log << MSG::DEBUG << "Hit DriftDistLeft " << hit->getDriftDistLeft() << endreq;
1875 log << MSG::DEBUG << "MDC Hit ADC " << hit->getAdc() << endreq;
1876 log << MSG::DEBUG << "Hit MDC Identifier " << hit->getMdcId() << endreq;
1877 log << MSG::DEBUG << "Hit Chisq " << hit->getChisqAdd() << endreq;
1878#endif
1879 hitcol->push_back(hit);
1880 SmartRef<RecMdcHit> refhit(hit);
1881 hitrefvec.push_back(refhit);
1882 }
1883 }
1884
1885 //stereo segments part
1886 FTList<FTSegment *>& seglist_st = fttrk->stereo_segments();
1887 int nseg_st = seglist_st.length();
1888 int ntrackster = 0;
1889
1890 for (int iseg_st = 0; iseg_st < nseg_st; iseg_st++) {
1891 FTList<FTWire *>& hitlist_st = seglist_st[iseg_st]->wireHits();
1892 int nhitlist_st = hitlist_st.length();
1893 ntrackhits += nhitlist_st;
1894 ntrackster += nhitlist_st;
1895 for( int ihit_st = 0; ihit_st < nhitlist_st; ihit_st++,hitindex++) {
1896 FTWire wire_st = *(hitlist_st[ihit_st]);
1897 double dd = wire_st.distance();
1898 //double adc = wire_st.getAdc();
1899 double adc = RawDataUtil::MdcChargeChannel(wire_st.getAdc());
1900 double tdc = wire_st.time();
1901 double chi2 = wire_st.getChi2();
1902 const Identifier mdcid = MdcID::wire_id(wire_st.layer().layerId(),wire_st.localId());
1903
1904 RecMdcHit* hit = new RecMdcHit;
1905 hit->setId(hitindex);
1906 hit->setDriftDistLeft( dd );
1907 hit->setDriftDistRight( dd );
1908 hit->setAdc( adc );
1909 hit->setTdc( tdc);
1910 hit->setMdcId( mdcid );
1911 hit->setChisqAdd( chi2 );
1912#ifndef OnlineMode
1913 log << MSG::DEBUG << "Z Hit DriftDistLeft " << hit->getDriftDistLeft() << endreq;
1914 log << MSG::DEBUG << "Z MDC Hit ADC " << hit->getAdc() << endreq;
1915 log << MSG::DEBUG << "Z Hit MDC Identifier " << hit->getMdcId() << endreq;
1916 log << MSG::DEBUG << "Z Hit Chisq " << hit->getChisqAdd() << endreq;
1917#endif
1918 hitcol->push_back(hit);
1919 SmartRef<RecMdcHit> refhit(hit);
1920 hitrefvec.push_back(refhit);
1921 }
1922 }
1923 trk->setNhits(ntrackhits);
1924 trk->setNster(ntrackster);
1925 trk->setVecHits(hitrefvec);
1926 trkcol->push_back(trk);
1927 trkid++;
1928#ifndef OnlineMode
1929 g_naxialhit->fill((float)(ntrackhits-ntrackster), 1.0);
1930 g_nstereohit->fill((float)ntrackster, 1.0);
1931 g_nhit->fill((float)ntrackhits, 1.0);
1932#endif
1933 }
1934
1935 IDataManagerSvc *dataManSvc = dynamic_cast<IDataManagerSvc*> (eventSvc);
1936 DataObject *aTrackCol;
1937 eventSvc->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
1938 if(aTrackCol != NULL) {
1939 dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
1940 eventSvc->unregisterObject("/Event/Recon/RecMdcTrackCol");
1941 }
1942 DataObject *aRecHitCol;
1943 eventSvc->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
1944 if(aRecHitCol != NULL) {
1945 dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
1946 eventSvc->unregisterObject("/Event/Recon/RecMdcHitCol");
1947 }
1948
1949 //register RecMdcHitCol
1950 StatusCode hitsc;
1951 hitsc = eventSvc->registerObject("/Event/Recon/RecMdcHitCol", hitcol);
1952 if( hitsc.isFailure() ) {
1953 log << MSG::FATAL << "Could not register RecMdcHitCol" << endreq;
1954 return hitsc;
1955 }
1956#ifndef OnlineMode
1957 log << MSG::DEBUG << "RecMdcHitCol registered successfully!" <<endreq;
1958
1959 RecMdcTrackCol::iterator bef = trkcol->begin();
1960 for( ; bef != trkcol->end(); bef++) {
1961 log <<MSG::DEBUG <<" registered kappa"<<(*bef)->helix(2)
1962 <<" registered tanl"<< (*bef)->helix(4)
1963 <<endreq;
1964 }
1965#endif
1966
1967 //register RecMdcTrackCol
1968 StatusCode trksc;
1969 trksc = eventSvc->registerObject("/Event/Recon/RecMdcTrackCol", trkcol);
1970 if( trksc.isFailure() ) {
1971 log << MSG::FATAL << "Could not register RecMdcHit" << endreq;
1972 return trksc;
1973 }
1974#ifndef OnlineMode
1975 log << MSG::DEBUG << "RecMdcTrackCol registered successfully!" <<endreq;
1976
1977 //check the result:RecMdcHitCol
1978 SmartDataPtr<RecMdcHitCol> newhitCol(eventSvc,"/Event/Recon/RecMdcHitCol");
1979 if (!newhitCol) {
1980 log << MSG::FATAL << "Could not find RecMdcHit" << endreq;
1981 return( StatusCode::FAILURE);
1982 }
1983 log << MSG::DEBUG << "Begin to check RecMdcHitCol"<<endreq;
1984 RecMdcHitCol::iterator iter_hit = newhitCol->begin();
1985 for( ; iter_hit != newhitCol->end(); iter_hit++){
1986 log << MSG::DEBUG << "retrieved MDC Hit:"
1987 << " DDL " <<(*iter_hit)->getDriftDistLeft()
1988 << " DDR " <<(*iter_hit)->getDriftDistRight()
1989 << " ADC " <<(*iter_hit)->getAdc() << endreq;
1990 }
1991
1992 //check the result:RecMdcTrackCol
1993 SmartDataPtr<RecMdcTrackCol> newtrkCol(eventSvc,"/Event/Recon/RecMdcTrackCol");
1994 if (!newtrkCol) {
1995 log << MSG::FATAL << "Could not find RecMdcTrackCol" << endreq;
1996 return( StatusCode::FAILURE);
1997 }
1998 log << MSG::DEBUG << "Begin to check RecMdcTrackCol"<<endreq;
1999 RecMdcTrackCol::iterator iter_trk = newtrkCol->begin();
2000 for( ; iter_trk != newtrkCol->end(); iter_trk++){
2001 log << MSG::DEBUG << "retrieved MDC track:"
2002 << "Track Id: " << (*iter_trk)->trackId()
2003 << " Pivot: " << (*iter_trk)->poca()[0] << " "
2004 << (*iter_trk)->poca()[1] << " " << (*iter_trk)->poca()[2]
2005 << endreq ;
2006
2007 /*<< "Phi0: " << (*iter_trk)->getFi0() << " Error of Phi0 "
2008 << (*iter_trk)->getError()(2,2) << endreq
2009 /<< "kappa: " << (*iter_trk)->getCpa() << " Error of kappa "
2010 << (*iter_trk)->getError()(3,3) << endreq
2011 << "Tanl: " << (*iter_trk)->getTanl() << " Error of Tanl "
2012 << (*iter_trk)->getError()(5,5) << endreq
2013 << "Chisq of fit: "<< (*iter_trk)->getChisq()
2014 << " Phi terminal: "<< (*iter_trk)->getFiTerm()
2015 << endreq
2016 << "Number of hits: "<< (*iter_trk)->getNhits()
2017 << " Number of stereo hits " << (*iter_trk)->getNster()
2018 << endreq;*/
2019
2020 log << MSG::DEBUG << "hitList of this track:" << endreq;
2021 HitRefVec gothits = (*iter_trk)->getVecHits();
2022 HitRefVec::iterator it_gothit = gothits.begin();
2023 for( ; it_gothit != gothits.end(); it_gothit++){
2024 log << MSG::DEBUG << "hits Id: "<<(*it_gothit)->getId()
2025 << " hits DDL&DDR " <<(*it_gothit)->getDriftDistLeft()
2026 << " hits MDC IDentifier " <<(*it_gothit)->getMdcId()
2027 << endreq
2028 << " hits TDC " <<(*it_gothit)->getTdc()
2029 << " hits ADC " <<(*it_gothit)->getAdc() << endreq;
2030 }
2031
2032 }
2033#endif
2034
2035 return StatusCode::SUCCESS;
2036}
2037
dt
TGraphErrors * dt
Definition: AbsCor.cxx:72
sin
double sin(const BesAngle a)
Definition: BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213
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
exp
EvtComplex exp(const EvtComplex &c)
Definition: EvtComplex.hh:252
abs
double abs(const EvtComplex &c)
Definition: EvtComplex.hh:212
g_nstereohit
IHistogram1D * g_nstereohit
Definition: ntupleItem.h:43
g_pz
NTuple::Array< float > g_pz
Definition: FTFinder.cxx:88
g_theta
NTuple::Array< float > g_theta
Definition: FTFinder.cxx:89
num_3Dtrk
int num_3Dtrk
Definition: FTFinder.cxx:101
g_pxMC
NTuple::Array< float > g_pxMC
Definition: ntupleItem.h:24
g_pzMC
NTuple::Array< float > g_pzMC
Definition: FTFinder.cxx:83
g_ncellMC
IHistogram1D * g_ncellMC
Definition: ntupleItem.h:40
num_2Dtrk
int num_2Dtrk
Definition: MdcFastTrkAlg.cxx:34
g_pyMC
NTuple::Array< float > g_pyMC
Definition: FTFinder.cxx:83
g_dr
NTuple::Array< float > g_dr
Definition: ntupleItem.h:33
g_hitmap
IHistogram2D * g_hitmap[20]
Definition: ntupleItem.h:45
g_kappa
NTuple::Array< float > g_kappa
Definition: FTFinder.cxx:91
g_eventNo
NTuple::Item< long > g_eventNo
Definition: FTFinder.cxx:81
g_ptMC
NTuple::Array< float > g_ptMC
Definition: FTFinder.cxx:83
g_estime
NTuple::Item< float > g_estime
Definition: ntupleItem.h:34
g_dz
NTuple::Array< float > g_dz
Definition: FTFinder.cxx:91
g_tanl
NTuple::Array< float > g_tanl
Definition: FTFinder.cxx:91
g_ntrk
NTuple::Item< long > g_ntrk
Definition: ntupleItem.h:28
g_ntrkMC
NTuple::Item< long > g_ntrkMC
Definition: ntupleItem.h:22
g_phi0MC
NTuple::Array< float > g_phi0MC
Definition: FTFinder.cxx:82
num_finaltrk
int num_finaltrk
Definition: FTFinder.cxx:101
g_py
NTuple::Array< float > g_py
Definition: FTFinder.cxx:88
g_px
NTuple::Array< float > g_px
Definition: ntupleItem.h:30
g_ncell
IHistogram1D * g_ncell
Definition: ntupleItem.h:41
g_pt
NTuple::Array< float > g_pt
Definition: FTFinder.cxx:88
g_nhit
IHistogram1D * g_nhit
Definition: ntupleItem.h:44
g_p
NTuple::Array< float > g_p
Definition: FTFinder.cxx:88
g_phi0
NTuple::Array< float > g_phi0
Definition: FTFinder.cxx:91
g_vx
NTuple::Array< float > g_vx
Definition: ntupleItem.h:32
g_naxialhit
IHistogram1D * g_naxialhit
Definition: ntupleItem.h:42
g_vy
NTuple::Array< float > g_vy
Definition: FTFinder.cxx:90
g_theta0MC
NTuple::Array< float > g_theta0MC
Definition: ntupleItem.h:23
g_vz
NTuple::Array< float > g_vz
Definition: FTFinder.cxx:90
g_phi
NTuple::Array< float > g_phi
Definition: ntupleItem.h:31
FTWireHit
#define FTWireHit
Definition: FTWire.h:12
s
XmlRpcServer s
Definition: HelloServer.cpp:11
weight
*********Class see also m_nmax DOUBLE PRECISION m_MasPhot DOUBLE PRECISION m_phsu DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_r2 DOUBLE PRECISION m_WtMass INTEGER m_nmax INTEGER m_Nevgen INTEGER m_IsFSR INTEGER m_MarTot *COMMON c_KarFin $ !Output file $ !Event serial number $ !alpha QED at Thomson limit $ !minimum energy at CMS for remooval $ !infrared dimensionless $ !dummy photon IR regulator $ !crude photon multiplicity enhancement factor *EVENT $ !MC crude volume of PhhSpace *Sfactors $ !YFS formfactor IR part only $ !YFS formfactor non IR finite part $ !mass weight
Definition: KarFin.h:34
v
**********Class see also m_nmax DOUBLE PRECISION m_amel DOUBLE PRECISION m_x2 DOUBLE PRECISION m_alfinv DOUBLE PRECISION m_Xenph INTEGER m_KeyWtm INTEGER m_idyfs DOUBLE PRECISION m_zini DOUBLE PRECISION m_q2 DOUBLE PRECISION m_Wt_KF DOUBLE PRECISION m_WtCut INTEGER m_KFfin *COMMON c_KarLud $ !Input CMS energy[GeV] $ !CMS energy after beam spread beam strahlung[GeV] $ !Beam energy spread[GeV] $ !z boost due to beam spread $ !electron beam mass *ff pair spectrum $ !minimum v
Definition: KarLud.h:35
m_pz
NTuple::Item< double > m_pz
Definition: MdcHistItem.h:78
MdcDigiVec
std::vector< MdcDigi * > MdcDigiVec
Definition: MdcRawDataProvider.h:11
RecEsTimeCol
ObjectVector< RecEsTime > RecEsTimeCol
Definition: RecEsTime.h:53
RecMdcHitCol
ObjectVector< RecMdcHit > RecMdcHitCol
Definition: RecMdcHit.h:99
RecMdcTrackCol
ObjectVector< RecMdcTrack > RecMdcTrackCol
Definition: RecMdcTrack.h:79
HitRefVec
SmartRefVector< RecMdcHit > HitRefVec
Definition: RecMdcTrack.h:22
n2
int n2
Definition: SD0Tag.cxx:55
msgSvc
IMessageSvc * msgSvc()
Definition: ServiceAccessor.h:13
M_PI
#define M_PI
Definition: TConstant.h:4
DstMdcTrack::setTrackId
void setTrackId(const int trackId)
Definition: DstMdcTrack.h:82
DstMdcTrack::setNster
void setNster(const int ns)
Definition: DstMdcTrack.h:98
DstMdcTrack::setError
void setError(double err[15])
Definition: DstMdcTrack.cxx:112
DstMdcTrack::setHelix
void setHelix(double helix[5])
Definition: DstMdcTrack.cxx:98
DstMdcTrack::helix
const HepVector helix() const
......
Definition: DstMdcTrack.cxx:151
DstMdcTrack::setChi2
void setChi2(const double chi)
Definition: DstMdcTrack.h:96
FTFinder::t0Estime
float t0Estime
Definition: FTFinder.h:161
FTFinder::vertex
CLHEP::Hep3Vector vertex(void) const
returns event primary vertex
Definition: FTFinder.cxx:1641
FTFinder::superLayer
FTSuperLayer * superLayer(int id) const
returns superlayer
Definition: FTFinder.h:196
FTFinder::i_rPhiFit
int i_rPhiFit
Definition: FTFinder.h:166
FTFinder::event
void event()
track finder core
Definition: FTFinder.cxx:472
FTFinder::tEstFlag
int tEstFlag
Definition: FTFinder.h:167
FTFinder::init
void init()
initializer(creates geometry)
Definition: FTFinder.cxx:159
FTFinder::tOffSet
float tOffSet
Definition: FTFinder.h:162
FTFinder::evtTiming
float evtTiming
Definition: FTFinder.h:163
FTFinder::tEstime
FTList< float > tEstime[10]
Definition: FTFinder.h:168
FTFinder::eventNo
int eventNo
Definition: FTFinder.h:158
FTFinder::begin_run
void begin_run()
begin run function(reads constants)
Definition: FTFinder.cxx:331
FTFinder::tracks
FTList< FTTrack * > & tracks(void) const
returns track list
Definition: FTFinder.h:202
FTFinder::term
void term()
terminator
Definition: FTFinder.cxx:319
FTFinder::runNo
int runNo
Definition: FTFinder.h:159
FTFinder::expflag
int expflag
Definition: FTFinder.h:160
FTFinder::getWireId
int getWireId(FTWire *) const
returns wire ID for given FTWire object
Definition: FTFinder.h:208
FTFinder::FTFinder
FTFinder()
Constructors and destructor.
Definition: FTFinder.cxx:121
FTFinder::pivot
const HepPoint3D pivot
Definition: FTFinder.h:164
FTFinder::x2t
float x2t(const FTLayer &l, const float x) const
convert x to t
Definition: FTFinder.h:214
FTLayer::r
const float r(void) const
returns r form origin
Definition: FTLayer.h:141
FTLayer::layerId
const int layerId(void) const
returns layer ID
Definition: FTLayer.h:112
FTLayer::superLayer
const FTSuperLayer & superLayer(void) const
returns super-layer
Definition: FTLayer.h:183
FTLayer::csize
double csize(void) const
returns cell size
Definition: FTLayer.h:190
FTList
Definition: FTList.h:16
FTList::first
T & first(void) const
returns the first object in the list
Definition: FTList.h:224
FTList::deleteAll
void deleteAll(void)
delete all object and clear(allocated memory remains same)
Definition: FTList.h:266
FTList::replace
void replace(int i, T src)
replace index-th object by the object src
Definition: FTList.h:157
FTList::clear
void clear(void)
clear lists but the allocated memory remains same
Definition: FTList.h:182
FTList::remove
int remove(int &)
remove objects by index and returns decremented index and length
Definition: FTList.h:140
FTList::length
int length(void) const
returns the length of the list
Definition: FTList.h:245
FTList::append
int append(const T &x)
append an object into the end of the list
Definition: FTList.h:127
FTSegment::incomingPhi
float incomingPhi(void) const
returns phi of incoming position
Definition: FTSegment.h:264
FTSegment::printout
void printout(void)
Definition: FTSegment.cxx:193
FTSegment::outgoingY
float outgoingY(void) const
returns y of outgoing position
Definition: FTSegment.h:243
FTSegment::outgoingPhi
float outgoingPhi(void) const
returns phi of outgoing position
Definition: FTSegment.h:280
FTSegment::innerBoundHits
FTList< FTWire * > & innerBoundHits(void) const
returns innerBoundHits
Definition: FTSegment.h:222
FTSegment::incomingY
float incomingY(void) const
returns y of incoming position
Definition: FTSegment.h:257
FTSegment::update
void update(void)
update information for axial segment
Definition: FTSegment.cxx:83
FTSegment::wireHits
FTList< FTWire * > & wireHits(void) const
returns wire-hit list
Definition: FTSegment.h:215
FTSegment::outgoingX
float outgoingX(void) const
returns x of outgoing position
Definition: FTSegment.h:236
FTSegment::kappa
float kappa(void) const
returns kappa(axial)
Definition: FTSegment.h:313
FTSegment::connect_outer
void connect_outer(const FTList< FTWire * > &, const FTList< FTWire * > &)
connect short segments
Definition: FTSegment.h:170
FTSegment::outerBoundHits
FTList< FTWire * > & outerBoundHits(void) const
returns outerBoundHits
Definition: FTSegment.h:229
FTSegment::incomingX
float incomingX(void) const
returns x of incoming position
Definition: FTSegment.h:250
FTSuperLayer::segments
FTList< FTSegment * > & segments(void) const
returns segement list
Definition: FTSuperLayer.h:215
FTSuperLayer::appendHit
void appendHit(FTWire *)
append wireHit to the list of hits
Definition: FTSuperLayer.h:233
FTTrack::helix
Helix * helix(void) const
returns helix parameters
Definition: FTTrack.h:251
FTTrack::axial_segments
FTList< FTSegment * > & axial_segments(void) const
returns axial segments
Definition: FTTrack.h:191
FTTrack::setFTFinder
void setFTFinder(FTFinder *)
Definition: FTTrack.h:323
FTTrack::Zav
const zav & Zav(void) const
returns zav
Definition: FTTrack.h:244
FTTrack::stereo_segments
FTList< FTSegment * > & stereo_segments(void) const
returns stereo_segments
Definition: FTTrack.h:198
FTTrack::lpav
const Lpav & lpav(void) const
returns lpav
Definition: FTTrack.h:237
FTTrack::chi2_kappa_tmp
float chi2_kappa_tmp(void) const
returns sigmaKappa at linking
Definition: FTTrack.h:230
FTWire
Definition: FTWire.h:43
FTWire::getAdc
float getAdc(void) const
Definition: FTWire.h:591
FTWire::y
const float y(void) const
returns position y
Definition: FTWire.h:373
FTWire::distance
float distance(void) const
returns drift distance
Definition: FTWire.h:401
FTWire::getChi2
float getChi2(void) const
Definition: FTWire.h:600
FTWire::state
unsigned state(void) const
returns state
Definition: FTWire.h:480
FTWire::initNeighbor
void initNeighbor(void)
initNeighbor
Definition: FTWire.h:255
FTWire::layer
const FTLayer & layer(void) const
returns layer
Definition: FTWire.h:359
FTWire::setAdc
void setAdc(float adc)
wangdy add:set Adc value
Definition: FTWire.h:595
FTWire::localId
const int localId(void) const
returns local ID
Definition: FTWire.h:380
FTWire::x
const float x(void) const
returns position x
Definition: FTWire.h:366
FTWire::wireId
void wireId(int wireID)
set wireId
Definition: FTWire.h:580
FTWire::time
float time(void) const
rerurns TDC time(after t0 subtraction)
Definition: FTWire.h:457
FTWire::setChi2
void setChi2(float chi2)
set residual fit chi2
Definition: FTWire.h:604
Helix::momentum
Hep3Vector momentum(double dPhi=0.) const
returns momentum vector after rotating angle dPhi in phi direction.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-08/src/Helix.cxx:184
Helix::Ea
const HepSymMatrix & Ea(void) const
returns error matrix.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-08/TrackUtil/Helix.h:269
Helix::x
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-08/src/Helix.cxx:129
Helix::dr
double dr(void) const
returns an element of parameters.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-08/TrackUtil/Helix.h:227
Helix::a
const HepVector & a(void) const
returns helix parameters.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-08/TrackUtil/Helix.h:263
Helix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-08/TrackUtil/Helix.h:209
IMagneticFieldSvc::getReferField
virtual double getReferField()=0
IMdcGeomSvc::getLayerSize
virtual const int getLayerSize()=0
IMdcGeomSvc::Layer
virtual const MdcGeoLayer *const Layer(unsigned id)=0
IMdcGeomSvc::Wire
virtual const MdcGeoWire *const Wire(unsigned id)=0
IMdcGeomSvc::getWireSize
virtual const int getWireSize()=0
IMdcGeomSvc::getSuperLayerSize
virtual const int getSuperLayerSize()=0
IRawDataProviderSvc::getMdcDigiVec
virtual MdcDigiVec & getMdcDigiVec(uint32_t control=0)=0
Lpar::d
double d(double x, double y) const
Definition: KalFitAlg/KalFitAlg-00-07-55-p03/KalFitAlg/lpav/Lpar.h:208
Lpar::Hpar
HepVector Hpar(const HepPoint3D &pivot) const
Definition: KalFitAlg/KalFitAlg-00-07-55-p03/KalFitAlg/lpav/Lpar.h:269
Lpar::phi
double phi(double r, int dir=0) const
Definition: KalFitAlg/KalFitAlg-00-07-55-p03/src/lpav/Lpar.cxx:187
MdcGeoLayer::Radius
double Radius(void) const
Definition: MdcGeoLayer.h:160
MdcGeoLayer::Slant
double Slant(void) const
Definition: MdcGeoLayer.h:158
MdcGeoLayer::Length
double Length(void) const
Definition: MdcGeoLayer.h:161
MdcGeoLayer::Sup
MdcGeoSuper * Sup(void) const
Definition: MdcGeoLayer.h:174
MdcGeoLayer::Offset
double Offset(void) const
Definition: MdcGeoLayer.h:166
MdcGeoSuper::Type
int Type(void) const
Definition: MdcGeoSuper.h:35
MdcGeoWire::Lyr
MdcGeoLayer * Lyr(void) const
Definition: MdcGeoWire.h:140
MdcGeoWire::Forward
HepPoint3D Forward(void) const
Definition: MdcGeoWire.h:129
MdcGeoWire::Backward
HepPoint3D Backward(void) const
Definition: MdcGeoWire.h:128
MdcID::wire_id
static Identifier wire_id(int wireType, int layer, int wire)
For a single wire.
Definition: MdcID.cxx:77
MdcID::layer
static int layer(const Identifier &id)
Values of different levels (failure returns 0)
Definition: MdcID.cxx:49
MdcID::wire
static int wire(const Identifier &id)
Definition: MdcID.cxx:54
MdcParameter::_findEventVertex
const int _findEventVertex
Definition: MdcParameter.h:35
MdcParameter::_mkMdst
const bool _mkMdst
Definition: MdcParameter.h:49
MdcParameter::_D_phi1
const float _D_phi1
Definition: MdcParameter.h:67
MdcParameter::_D_phi2
const float _D_phi2
Definition: MdcParameter.h:67
MdcParameter::_chi2_2
const float _chi2_2
Definition: MdcParameter.h:77
MdcParameter::instance
static MdcParameter * instance()
Definition: MdcParameter.cxx:6
MdcParameter::_chi2_kappa
const float _chi2_kappa
Definition: MdcParameter.h:59
MdcParameter::_chi2_1
const float _chi2_1
Definition: MdcParameter.h:77
MdcParameter::_doIt
const int _doIt
Definition: MdcParameter.h:47
MdcParameter::_hitscut
const int _hitscut
Definition: MdcParameter.h:81
MdcParameter::_Min_chi2
const float _Min_chi2
Definition: MdcParameter.h:61
MdcParameter::_evtTimeCorr
const int _evtTimeCorr
Definition: MdcParameter.h:37
MdcParameter::_chi2_3
const float _chi2_3
Definition: MdcParameter.h:77
MdcParameter::_D_phi3
const float _D_phi3
Definition: MdcParameter.h:67
MdcParameter::_mkTds
const bool _mkTds
Definition: MdcParameter.h:51
RawDataUtil::MdcChargeChannel
static int MdcChargeChannel(double charge)
Definition: RawDataUtil.h:12
RawDataUtil::MdcTime
static double MdcTime(int timeChannel)
Definition: RawDataUtil.h:8
RawDataUtil::MdcCharge
static double MdcCharge(int chargeChannel)
Definition: RawDataUtil.h:11
RecMdcHit::setMdcId
void setMdcId(Identifier mdcid)
Definition: RecMdcHit.h:67
RecMdcHit::getChisqAdd
const double getChisqAdd(void) const
Definition: RecMdcHit.h:46
RecMdcHit::setDriftDistLeft
void setDriftDistLeft(double ddl)
Definition: RecMdcHit.h:60
RecMdcHit::getAdc
const double getAdc(void) const
Definition: RecMdcHit.h:51
RecMdcHit::getMdcId
const Identifier getMdcId(void) const
Definition: RecMdcHit.h:49
RecMdcHit::setTdc
void setTdc(double tdc)
Definition: RecMdcHit.h:68
RecMdcHit::getDriftDistLeft
const double getDriftDistLeft(void) const
Definition: RecMdcHit.h:42
RecMdcHit::setAdc
void setAdc(double adc)
Definition: RecMdcHit.h:69
RecMdcHit::setChisqAdd
void setChisqAdd(double pChisq)
Definition: RecMdcHit.h:64
RecMdcHit::setDriftDistRight
void setDriftDistRight(double ddr)
Definition: RecMdcHit.h:61
RecMdcHit::setId
void setId(int id)
Definition: RecMdcHit.h:58
RecMdcTrack::setPivot
void setPivot(const HepPoint3D &pivot)
Definition: RecMdcTrack.h:68
RecMdcTrack::setVecHits
void setVecHits(HitRefVec vechits)
Definition: RecMdcTrack.h:69
RecMdcTrack::setNhits
void setNhits(int nhits)
Definition: RecMdcTrack.h:67
RecMdcTrack::setFiTerm
void setFiTerm(double fiterm)
Definition: RecMdcTrack.h:66
Event
Definition: Event.h:21
EvtCyclic3::next
Index next(Index i)
Definition: EvtCyclic3.cc:107
g_nstereohit
IHistogram1D * g_nstereohit
Definition: ntupleItem.h:43
g_pz
NTuple::Array< float > g_pz
Definition: ntupleItem.h:30
g_theta
NTuple::Array< float > g_theta
Definition: ntupleItem.h:31
g_pxMC
NTuple::Array< float > g_pxMC
Definition: ntupleItem.h:24
g_pzMC
NTuple::Array< float > g_pzMC
Definition: ntupleItem.h:24
g_ncellMC
IHistogram1D * g_ncellMC
Definition: ntupleItem.h:40
g_pyMC
NTuple::Array< float > g_pyMC
Definition: ntupleItem.h:24
g_dr
NTuple::Array< float > g_dr
Definition: ntupleItem.h:33
g_hitmap
IHistogram2D * g_hitmap[20]
Definition: ntupleItem.h:45
g_kappa
NTuple::Array< float > g_kappa
Definition: ntupleItem.h:33
g_eventNo
NTuple::Item< long > g_eventNo
Definition: ntupleItem.h:22
g_ptMC
NTuple::Array< float > g_ptMC
Definition: ntupleItem.h:24
g_estime
NTuple::Item< float > g_estime
Definition: ntupleItem.h:34
g_dz
NTuple::Array< float > g_dz
Definition: ntupleItem.h:33
g_tanl
NTuple::Array< float > g_tanl
Definition: ntupleItem.h:33
g_ntrk
NTuple::Item< long > g_ntrk
Definition: ntupleItem.h:28
g_ntrkMC
NTuple::Item< long > g_ntrkMC
Definition: ntupleItem.h:22
g_phi0MC
NTuple::Array< float > g_phi0MC
Definition: ntupleItem.h:23
g_py
NTuple::Array< float > g_py
Definition: ntupleItem.h:30
g_px
NTuple::Array< float > g_px
Definition: ntupleItem.h:30
g_ncell
IHistogram1D * g_ncell
Definition: ntupleItem.h:41
g_pt
NTuple::Array< float > g_pt
Definition: ntupleItem.h:30
g_nhit
IHistogram1D * g_nhit
Definition: ntupleItem.h:44
g_p
NTuple::Array< float > g_p
Definition: ntupleItem.h:30
g_phi0
NTuple::Array< float > g_phi0
Definition: ntupleItem.h:33
g_vx
NTuple::Array< float > g_vx
Definition: ntupleItem.h:32
g_naxialhit
IHistogram1D * g_naxialhit
Definition: ntupleItem.h:42
g_vy
NTuple::Array< float > g_vy
Definition: ntupleItem.h:32
g_theta0MC
NTuple::Array< float > g_theta0MC
Definition: ntupleItem.h:23
g_vz
NTuple::Array< float > g_vz
Definition: ntupleItem.h:32
g_phi
NTuple::Array< float > g_phi
Definition: ntupleItem.h:31
t
int t()
Definition: t.c:1