BOSS 7.0.3
BESIII Offline Software System
Loading...
Searching...
No Matches
Hough.cxx
Go to the documentation of this file.
1#include "MdcHoughFinder/Hough.h"
2#include "GaudiKernel/MsgStream.h"
3#include "GaudiKernel/AlgFactory.h"
4#include "GaudiKernel/ISvcLocator.h"
5#include "GaudiKernel/SmartDataPtr.h"
6#include "GaudiKernel/IDataProviderSvc.h"
7#include "GaudiKernel/IPartPropSvc.h"
8#include "GaudiKernel/PropertyMgr.h"
9#include "GaudiKernel/IService.h"
10#include "GaudiKernel/NTuple.h"
11#include "GaudiKernel/INTupleSvc.h"
12#include "EventModel/EventHeader.h"
13#include "MdcRawEvent/MdcDigi.h"
14#include "Identifier/Identifier.h"
15#include "Identifier/MdcID.h"
16
17#include <iostream>
18#include <math.h>
19
20#include "EvTimeEvent/RecEsTime.h"
21#include "MdcGeom/EntranceAngle.h"
22#include "RawEvent/RawDataUtil.h"
23#include "MdcData/MdcHit.h"
24
25#include "McTruth/DecayMode.h"
26#include "McTruth/McParticle.h"
27#include "TrackUtil/Helix.h"
28#include "MdcRecoUtil/Pdt.h"
29
30#include "TrkBase/TrkFit.h"
31#include "TrkBase/TrkHitList.h"
32#include "TrkBase/TrkExchangePar.h"
33#include "TrkFitter/TrkHelixMaker.h"
34#include "TrkFitter/TrkCircleMaker.h"
35#include "TrkFitter/TrkLineMaker.h"
36#include "TrkFitter/TrkHelixFitter.h"
37#include "MdcxReco/Mdcxprobab.h"
38#include "MdcData/MdcRecoHitOnTrack.h"
39#include "MdcPrintSvc/IMdcPrintSvc.h"
40
41#include "MdcHoughFinder/HoughTrackList.h"
42#include "MdcHoughFinder/HoughTrack.h"
43#include <math.h>
44
45#include "EvTimeEvent/RecEsTime.h"
46//#include "MdcCalibFunSvc/IMdcCalibFunSvc.h"
47#include "McTruth/MdcMcHit.h"
48#include "TrkBase/TrkFit.h"
49#include "TrkBase/TrkHitList.h"
50#include "TrkBase/TrkExchangePar.h"
51#include "TrkFitter/TrkHelixMaker.h"
52#include "TrkFitter/TrkCircleMaker.h"
53#include "TrkFitter/TrkLineMaker.h"
54#include "TrkFitter/TrkHelixFitter.h"
55#include "MdcPrintSvc/IMdcPrintSvc.h"
56#include "MdcGeom/EntranceAngle.h"
57#include "TrackUtil/Helix.h"
58#include "GaudiKernel/IPartPropSvc.h"
59#include "MdcRecoUtil/Pdt.h"
60#include "RawEvent/RawDataUtil.h"
61#include "MdcData/MdcHit.h"
62#include "MdcTrkRecon/MdcMap.h"
63
64#include "TTree.h"
65#include "TH2D.h"
66
67#include "BesTimerSvc/IBesTimerSvc.h"
68#include "BesTimerSvc/BesTimerSvc.h"
69#include "BesTimerSvc/BesTimer.h"
70
71
72
73using namespace std;
74using namespace Event;
75
76MdcHoughFinder::MdcHoughFinder(const std::string& name, ISvcLocator* pSvcLocator) :
77 Algorithm(name, pSvcLocator)
78{
79 // Declare the properties
80 declareProperty("debug", m_debug = 0);
81 declareProperty("debugMap", m_debugMap = 0);
82 declareProperty("debug2D", m_debug2D = 0);
83 declareProperty("debugTrack", m_debugTrack = 0);
84 declareProperty("debugStereo", m_debugStereo= 0);
85 declareProperty("debugZs", m_debugZs= 0);
86 declareProperty("debug3D", m_debug3D= 0);
87 declareProperty("debugArbHit", m_debugArbHit= 0);
88 declareProperty("hist", m_hist = 1);
89 declareProperty("filter", m_filter= 0);
90 //read raw data setup
91 declareProperty("keepBadTdc", m_keepBadTdc = 0);
92 declareProperty("dropHot", m_dropHot= 0);
93 declareProperty("keepUnmatch", m_keepUnmatch = 0);
94 // combine pattsf
95 declareProperty("combineTracking",m_combineTracking = false);
96 declareProperty("removeBadTrack", m_removeBadTrack = 0);
97 declareProperty("dropTrkDrCut", m_dropTrkDrCut= 1.);
98 declareProperty("dropTrkDzCut", m_dropTrkDzCut= 10.);
99 declareProperty("dropTrkPtCut", m_dropTrkPtCut= 0.12);
100 declareProperty("dropTrkChi2Cut", m_dropTrkChi2Cut = 10000.);
101 //input setup
102 declareProperty("inputType", m_inputType = 0);
103 //set MdcHoughFinder map
104 declareProperty("mapCharge", m_mapCharge= -1); //0 use whole ; 1 only half
105 declareProperty("useHalf", m_useHalf= 0); //0 use whole ; 1 only half
106 declareProperty("mapHitStyle", m_mapHit= 0); //0 : all ; 1 :axial
107 declareProperty("nbinTheta", m_nBinTheta = 100);
108 declareProperty("nbinRho", m_nBinRho = 100);
109 declareProperty("rhoRange", m_rhoRange = 0.1);
110 declareProperty("peakWidth", m_peakWidth= 3);
111 declareProperty("peakHigh", m_peakHigh= 1);
112 declareProperty("hitPro", m_hitPro= 0.4);
113
114 declareProperty("recpos", m_recpos= 1);
115 declareProperty("recneg", m_recneg= 1);
116 declareProperty("combineSelf", m_combine= 1);
117 declareProperty("z0CutCompareHough", m_z0Cut_compareHough= 10);
118
119 //split drift circle
120 declareProperty("n1", m_npart= 100);
121 declareProperty("n2", m_n2= 16);
122
123 declareProperty("d1", m_d1= 0.2);
124 declareProperty("d2", m_d2= 0.2);
125
126 declareProperty("pdtFile", m_pdtFile = "pdt.table");
127 declareProperty("eventFile", m_evtFile= "EventList");
128
129}
130//
131StatusCode MdcHoughFinder::beginRun(){
132 //Initailize MdcDetector
133 Global::m_gm = MdcDetector::instance(0);
134 if(NULL == Global::m_gm) return StatusCode::FAILURE;
135
136 return StatusCode::SUCCESS;
137}
138
139// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
140StatusCode MdcHoughFinder::initialize(){
141
142 MsgStream log(msgSvc(), name());
143 log << MSG::INFO << "in initialize()" << endreq;
144
145 StatusCode sc;
146
147 //particle
148 IPartPropSvc* p_PartPropSvc;
149 static const bool CREATEIFNOTTHERE(true);
150 sc = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
151 if (!sc.isSuccess() || 0 == p_PartPropSvc) {
152 log << MSG::ERROR << " Could not initialize PartPropSvc" << endreq;
153 return sc;
154 }
155 m_particleTable = p_PartPropSvc->PDT();
156
157 // RawData
158 IRawDataProviderSvc* irawDataProviderSvc;
159 sc = service ("RawDataProviderSvc", irawDataProviderSvc);
160 m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*> (irawDataProviderSvc);
161 if ( sc.isFailure() ){
162 log << MSG::FATAL << name()<<" Could not load RawDataProviderSvc!" << endreq;
163 return StatusCode::FAILURE;
164 }
165
166 // Geometry
167 IMdcGeomSvc* imdcGeomSvc;
168 sc = service ("MdcGeomSvc", imdcGeomSvc);
169 m_mdcGeomSvc = dynamic_cast<MdcGeomSvc*> (imdcGeomSvc);
170 if ( sc.isFailure() ){
171 log << MSG::FATAL << "Could not load MdcGeoSvc!" << endreq;
172 return StatusCode::FAILURE;
173 }
174
175 // magneticfield
176 sc = service ("MagneticFieldSvc",m_pIMF);
177 if(sc!=StatusCode::SUCCESS) {
178 log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
179 }
180 m_bfield = new BField(m_pIMF);
181 log << MSG::INFO << "field z = "<<m_bfield->bFieldNominal()<< endreq;
182 m_context = new TrkContextEv(m_bfield);
183
184 //read pdt
185 Pdt::readMCppTable(m_pdtFile);
186
187 //Get MdcCalibFunSvc
188 IMdcCalibFunSvc* imdcCalibSvc;
189 sc = service ("MdcCalibFunSvc", imdcCalibSvc);
190 m_mdcCalibFunSvc = dynamic_cast<MdcCalibFunSvc*>(imdcCalibSvc);
191 if ( sc.isFailure() ){
192 log << MSG::FATAL << "Could not load MdcCalibFunSvc!" << endreq;
193 return StatusCode::FAILURE;
194 }
195
196 //initialize MdcPrintSvc
197 IMdcPrintSvc* imdcPrintSvc;
198 sc = service ("MdcPrintSvc", imdcPrintSvc);
199 m_mdcPrintSvc = dynamic_cast<MdcPrintSvc*> (imdcPrintSvc);
200 if ( sc.isFailure() ){
201 log << MSG::FATAL << "Could not load MdcPrintSvc!" << endreq;
202 return StatusCode::FAILURE;
203 }
204
205 //time
206 sc = service( "BesTimerSvc", m_timersvc);
207 if( sc.isFailure() ) {
208 log << MSG::WARNING << " Unable to locate BesTimerSvc" << endreq;
209 return StatusCode::FAILURE;
210 }
211 m_timer_all = m_timersvc->addItem("Execution");
212 m_timer_all->propName("nExecution");
213
214 //initialize static
215 HoughHit::setMdcCalibFunSvc(m_mdcCalibFunSvc);
216 HoughHit::setMdcGeomSvc(m_mdcGeomSvc);
217 HoughHit::setBunchTime(m_bunchT0);
218 Hough2D::setContext(m_context);
219 Hough2D::setCalib(m_mdcCalibFunSvc);
220 Hough3D::setContext(m_context);
221 Hough3D::setCalib(m_mdcCalibFunSvc);
222
223 HoughHit::_npart=m_npart;
224 HoughMap::m_useHalfCir=m_useHalf;
225 HoughMap::m_N1=m_npart;
226 HoughMap::m_N2=m_n2;
227
228 HoughPeak::m_dcut1=m_d1;
229 HoughPeak::m_dcut2=m_d2;
230
231 if(m_debugMap> 0) HoughMap::m_debug = 1;
232 if(m_debug2D > 0) Hough2D::m_debug = 1;
233 if(m_debug3D > 0) Hough3D::m_debug = 1;
234 if(m_debugTrack> 0) HoughTrack::m_debug = 1;
235 if(m_debugStereo> 0) HoughStereo::m_debug = 1;
236 if(m_debugZs> 0) HoughZsFit::m_debug = 1;
237 if(m_debug3D > 4) TrkHelixFitter::m_debug = 1;
238
239 if(m_hist) sc = bookTuple();
240 if ( sc.isFailure() ){
241 log << MSG::FATAL << "Could not book Tuple !" << endreq;
242 return StatusCode::FAILURE;
243 }
244 return StatusCode::SUCCESS;
245}
246
247// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
248StatusCode MdcHoughFinder::execute() {
249
250 MsgStream log(msgSvc(), name());
251 log << MSG::INFO << "in execute()" << endreq;
252 cout.precision(6);
253
254
255 m_timer_all->start();
256
257 //event start time
258 SmartDataPtr<RecEsTimeCol> evTimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
259 if (!evTimeCol) {
260 log << MSG::WARNING << "Could not find EvTimeCol" << endreq;
261 return StatusCode::SUCCESS;
262 }else{
263 RecEsTimeCol::iterator iter_evt = evTimeCol->begin();
264 if (iter_evt != evTimeCol->end()){
265 m_bunchT0 = (*iter_evt)->getTest()*1.e-9;//m_bunchT0-s, getTest-ns
266 }
267 }
268 HoughHit::setBunchTime(m_bunchT0);
269
270 //-- Get the event header, print out event and run number
271 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
272 if (!eventHeader) {
273 log << MSG::FATAL << "Could not find Event Header" << endreq;
274 return StatusCode::FAILURE;
275 }
276
277 t_eventNum=eventHeader->eventNumber();
278 t_runNum=eventHeader->runNumber();
279 if(m_hist){
280 m_eventNum=eventHeader->eventNumber();
281 m_runNum=eventHeader->runNumber();
282 }
283 if(m_debug>0) cout<<"begin evt "<<eventHeader->eventNumber()<<endl;
284
285 //prepare tds
286 RecMdcTrackCol* trackList_tds;
287 RecMdcHitCol* hitList_tds;
288 vector<RecMdcTrack*> vec_trackPatTds;
289 int nTdsTk = storeTracks(trackList_tds,hitList_tds,vec_trackPatTds);
290 //print track in pattsf with bad vertex
291 if(m_debug>0){
292 RecMdcTrackCol::iterator iteritrk_pattsf = vec_trackPatTds.begin();
293 for(;iteritrk_pattsf!=vec_trackPatTds.end();iteritrk_pattsf++){
294 cout<<"in PATTSF LOST: "<<(*iteritrk_pattsf)->helix(0)<<" "<<(*iteritrk_pattsf)->helix(1)<<" "<<(*iteritrk_pattsf)->helix(2)<<" "<<(*iteritrk_pattsf)->helix(3)<<" "<<(*iteritrk_pattsf)->helix(4)<<" chi2 "<<(*iteritrk_pattsf)->chi2()<<endl;
295 }
296 }
297
298 //for arbi hits
299 MdcTrackParams m_par;
300 if(m_debugArbHit>0 ) m_par.lPrint=8;
301 m_par.lRemoveInActive=1;
302 //m_par.lUseQualCuts=0;
303 //m_par.maxGapLength=99;
304 //m_par.maxChisq=99;
305 //m_par.minHits=99;
306
307 // if filter read eventNum in file
308 if(m_filter){
309 ifstream lowPt_Evt;
310 lowPt_Evt.open(m_evtFile.c_str());
311 vector<int> evtlist;
312 int evtNum;
313 while( lowPt_Evt >> evtNum) {
314 evtlist.push_back(evtNum);
315 }
316 vector<int>::iterator iter_evt = find(evtlist.begin(),evtlist.end(),t_eventNum);
317 if( iter_evt == evtlist.end() ) { setFilterPassed(false); return StatusCode::SUCCESS; }
318 else setFilterPassed(true);
319 }
320
321 if(m_inputType == -1) GetMcInfo();
322
323 //-- Get MDC digi vector
324 if(m_debug>0) cout<<"step1 : prepare digi "<<endl;
325 MdcDigiVec mdcDigiVec = prepareDigi();
326
327 //-- Create MdcHoughFinder hit list
328 bool debugTruth = false;
329 if(m_inputType == -1) debugTruth= true;
330 if(m_debug>0) cout<<"step2 : hits-> hough hit list "<<endl;
331 HoughHitList houghHitList(mdcDigiVec);
332 // if( mdcDigiVec.size() < 10 ) return StatusCode::SUCCESS;
333 if( houghHitList.nHit() < 10 || houghHitList.nHit()>500 ) return StatusCode::SUCCESS;
334 if(m_debug>0) houghHitList.printAll();
335 if(debugTruth) houghHitList.addTruthInfo(g_tkParTruth);
336
337 vector<MdcHit*> mdcHitCol_neg;
338 vector<MdcHit*> mdcHitCol_pos;
339 MdcDigiVec::iterator iter = mdcDigiVec.begin();
340 for (;iter != mdcDigiVec.end(); iter++) {
341 const MdcDigi* digi = (*iter);
342 if( HoughHit(digi).driftTime()>1000 || HoughHit(digi).driftTime()<=0 ) continue;
343 MdcHit *mdcHit= new MdcHit(digi, Global::m_gm);
344 MdcHit *mdcHit_pos= new MdcHit(digi, Global::m_gm);
345 mdcHitCol_neg.push_back(mdcHit);
346 mdcHitCol_pos.push_back(mdcHit_pos);
347 }
348
349 HoughMap *m_map = new HoughMap(-1,houghHitList,m_mapHit,m_nBinTheta,m_nBinRho,-1.*m_rhoRange,m_rhoRange,m_peakWidth,m_peakHigh,m_hitPro);//, 2dOr3d);
350
351 // if(m_hist) dumpHoughHitList(houghHitList);
352 if(m_hist) m_nHit = houghHitList.nHit();
353 if(m_hist) dumpHoughMap(*m_map);
354
355 //track
356 if(m_debug>0) cout<<"step3 : neg track list "<<endl;
357 vector<HoughTrack*> vec_track_neg;
358 vector<HoughTrack*> vec_track2D_neg;
359 MdcTrackList mdcTrackList_neg(m_par) ;
360 if(m_recneg){
361 HoughTrackList trackList_neg(*m_map);
362 int trackList_size = trackList_neg.getTrackNum();
363 vector< vector<int> > stat_2d(2,vector<int>(trackList_size,0) );
364 vector< vector<int> > stat_3d(2,vector<int>(trackList_size,0) );
365 int ifit=0;
366 int ifit3D=0;
367 for(int itrack=0;itrack<trackList_size;itrack++){
368 if(m_debug>0) cout<<"begin track: "<<itrack<<endl;
369 //suppose charge -
370 if(m_debug>0) cout<<" suppose charge -1 "<<endl;
371 HoughTrack &track_neg = trackList_neg.getTrack(itrack);
372 track_neg.setMdcHit( &mdcHitCol_neg);
373 track_neg.setHoughHitList(houghHitList.getList());
374 track_neg.setCharge(-1);
375 stat_2d[0][itrack] = track_neg.fit2D(m_bunchT0);
376 int track_charge_2d = track_neg.trackCharge2D();
377 if(m_debug>0) cout<<" charge -1 stat2d "<<stat_2d[0][itrack]<<" "<<track_charge_2d<<endl;
378 if( stat_2d[0][itrack] == 0 || track_charge_2d ==0 ) continue;
379 //vec_track2D_neg.push_back( &track_neg );
380 //fill 2D inf
381 ifit++;
382 //3D fit
383 int nHit3d = track_neg.find_stereo_hit();
384 int npair= track_neg.find_pair_hit();
385 //cout<<"npair "<<npair<<endl;
386
387 int track_charge_3d = track_neg.trackCharge3D();
388 if(m_debug>0) cout<<" nhitstereo -1 "<<nHit3d<<" "<<track_charge_3d<<endl;
389 if( nHit3d <3 || track_charge_3d==0 ) continue;
390
391 //choose fit method
392 if( npair==0 ) stat_3d[0][itrack] = track_neg.fit3D();
393 else stat_3d[0][itrack] = track_neg.fit3D_inner();
394 //else stat_3d[0][itrack] = track_neg.fit3D();
395
396 if(m_debug>0) cout<<" charge -1 stat3d "<<stat_3d[0][itrack]<<" "<<endl;
397 if( stat_3d[0][itrack]==0 ) continue;
398 vec_track_neg.push_back( &track_neg );
399 //fill 3D inf
400 /*
401 if(m_hist){
402 m_tanl_neg[ifit3D] = track_neg.getTanl_zs();
403 m_tanl3D_neg[ifit3D] = track_neg.getTanl();
404 m_z0_neg[ifit3D] = track_neg.getZ0_zs();
405 m_z0_3D_neg[ifit3D] = track_neg.getZ0();
406 m_pro_neg[ifit3D] = track_neg.getPro();
407 m_pt3D_neg[ifit3D] = track_neg.getPt3D();
408 m_nHit3D_neg[ifit3D] = track_neg.getNfit3D();
409 m_chi2_3D_neg[ifit3D] = track_neg.getChi2_3D();
410 }
411 ifit3D++;
412 int axial_use=0;
413 int stereo_use=0;
414 int cir0=0;
415 int cirmore=0;
416 int ncir_track;
417 for(int ister =0;ister<track_neg.getHoughHitList().size();ister++){
418 HoughRecHit *rechit = &(track_neg.getHoughHitList().at(ister));
419 if( rechit->getSlayerType()==0 ) {
420 double deltaD = rechit->getDeltaD();
421 double flt= rechit->getFltLen();
422 if(m_hist){
423 m_axial_layer[axial_use]=rechit->getLayerId();
424 m_axial_wire[axial_use]=rechit->getWireId();
425 m_axial_deltaD[axial_use] = deltaD;
426 m_axial_flt[axial_use] = flt;
427 }
428 //judge which circle
429 int cirlist = rechit->getCirList();
430 if(cirlist == 0) cir0++;
431 else cirmore++;
432 axial_use++;
433 }
434 if( rechit->getSlayerType()!=0 ) {
435 //stereo hit info
436 double z0= rechit->getzAmb(0);
437 double z1= rechit->getzAmb(1);
438 double s0= rechit->getsAmb(0);
439 double s1= rechit->getsAmb(1);
440 int ambig = rechit->getAmbig();
441 int ambigTruth = rechit->getLrTruth();
442 int style = rechit->getStyle();
443 //cout<<"ambig ambigTruth "<<ambig<<" "<<ambigTruth<<endl;
444 double z = rechit->getzPos();
445 double s = rechit->getsPos();
446 double zTruth;
447 double sTruth;
448 if(ambigTruth == 1 ) {zTruth = z0;sTruth=s0;}
449 else if(ambigTruth == -1 ) {zTruth = z1;sTruth=s1;}
450 //calcu dist of hitz to zsline
451 double deltaZ = zTruth - ( sTruth * track_neg.getTanl_zs()+track_neg.getZ0_zs());
452 //else cout<<" not get ambig truth info "<<endl;
453 int nsol = rechit->getnsol();
454 double disToCir= rechit->getDisToCir();
455 int multicir = rechit->getCirList();
456 if(m_hist){
457 m_stereo_layer[stereo_use]=rechit->getLayerId();
458 m_stereo_wire[stereo_use]=rechit->getWireId();
459 m_stereo_style[stereo_use]=style;
460 m_stereo_z0[stereo_use]=z0;
461 m_stereo_z1[stereo_use]=z1;
462 m_stereo_s0[stereo_use]=s0;
463 m_stereo_s1[stereo_use]=s1;
464 m_stereo_z[stereo_use]=z;
465 m_stereo_s[stereo_use]=s;
466 m_stereo_zTruth[stereo_use]=zTruth;
467 m_stereo_sTruth[stereo_use]=sTruth;
468 m_stereo_deltaZ[stereo_use]=deltaZ;
469 m_stereo_nsol[stereo_use]=nsol;
470 m_stereo_ambig[stereo_use]=ambig;
471 m_stereo_ambig_truth[stereo_use]=ambigTruth;
472 m_stereo_disToCir[stereo_use]=disToCir;
473 m_stereo_cirlist[stereo_use]=multicir;
474 }
475 stereo_use++;
476 }
477 }
478 if( cir0>cirmore ) ncir_track = 0;
479 else ncir_track = 1;
480 if(m_hist){
481 m_evt_stereo= t_eventNum;
482 m_run_stereo= t_runNum;
483 m_cos_stereo= t_cos;
484 m_tanlTruth_stereo= t_tanl;
485 m_ncir_stereo= ncir_track;
486 m_npair_stereo= npair;
487 m_tanl_stereo= track_neg.getTanl_zs();
488 m_z0_stereo= track_neg.getZ0_zs();
489 m_tanl3D_stereo= track_neg.getTanl();
490 m_z03D_stereo= track_neg.getZ0();
491 m_nHit_axial= axial_use;
492 m_nHit_stereo = stereo_use;
493 ntuplehit->write();
494 }
495 */
496 }
497
498 std::sort ( vec_track_neg.begin(),vec_track_neg.end(),more_pt);
499 //track for clear
500 vector<MdcTrack*> vec_MdcTrack_neg;
501 for( unsigned int i =0;i<vec_track_neg.size();i++){
502 MdcTrack *mdcTrack = new MdcTrack(vec_track_neg[i]->getTrk());
503 vec_MdcTrack_neg.push_back(mdcTrack);
504 if(m_debug>0) cout<<"trackneg: "<<i<<" pt "<<vec_track_neg[i]->getPt3D()<<endl;
505 if(m_debug>0) vec_track_neg[i]->print();
506 }
507 if(m_debug>0) cout<<"step4 : judge neg track list "<<endl;
508 judgeHit(mdcTrackList_neg,vec_MdcTrack_neg);
509 if(m_debug>0) cout<<"finish - charge "<<endl;
510 }
511
512 //do rec pos charge
513
514 HoughMap *m_map2 = new HoughMap(1,houghHitList,m_mapHit,m_nBinTheta,m_nBinRho,-1.*m_rhoRange,m_rhoRange,m_peakWidth,m_peakHigh,m_hitPro);//, 2dOr3d);
515 if(m_debug>0) cout<<"step5 : pos track list "<<endl;
516 vector<HoughTrack*> vec_track_pos;
517 MdcTrackList mdcTrackList_pos(m_par) ;
518 if(m_recpos){
519 HoughTrackList tracklist_pos(*m_map2);
520 int tracklist_size2 = tracklist_pos.getTrackNum();
521 vector< vector<int> > stat_2d2(2,vector<int>(tracklist_size2,0) );
522 vector< vector<int> > stat_3d2(2,vector<int>(tracklist_size2,0) );
523 int ifit=0;
524 int ifit3D=0;
525 for(int itrack=0;itrack<tracklist_size2;itrack++){
526 //suppose charge +
527 if(m_debug>0) cout<<" suppose charge +1 "<<endl;
528 HoughTrack &track_pos = tracklist_pos.getTrack(itrack);
529 track_pos.setMdcHit( &mdcHitCol_pos);
530 track_pos.setHoughHitList(houghHitList.getList());
531 track_pos.setCharge(1);
532 stat_2d2[0][itrack] = track_pos.fit2D(m_bunchT0);
533 int track_charge_2d = track_pos.trackCharge2D();
534 if(m_debug>0) cout<<" charge +1 stat2d "<<stat_2d2[1][itrack]<<" "<<track_charge_2d<<endl;
535 if( stat_2d2[0][itrack] == 0 || track_charge_2d==0 ) continue;
536 //fill 2d inf
537 ifit++;
538 //3d fit
539 int nHit3d = track_pos.find_stereo_hit();
540 int npair= track_pos.find_pair_hit();
541 //cout<<"npair "<<npair<<endl;
542
543 int track_charge_3d = track_pos.trackCharge3D();
544 if(m_debug>0) cout<<" nhitstereo +1 "<<nHit3d<<" "<<track_charge_3d<<endl;
545 if( nHit3d <3 || track_pos.trackCharge3D()==0 ) continue;
546 //choose fit method
547 if( npair==0 ) stat_3d2[0][itrack] = track_pos.fit3D();
548 else stat_3d2[0][itrack] = track_pos.fit3D_inner();
549 //else stat_3d2[0][itrack] = track_neg.fit3D();
550
551 if(m_debug>0) cout<<" charge +1 stat3d "<<stat_3d2[1][itrack]<<" "<<endl;
552 if( stat_3d2[0][itrack]==0 ) continue;
553 vec_track_pos.push_back( &track_pos );
554 //fill 3d inf
555 ifit3D++;
556 }
557
558 //sort pos tracklist
559 std::sort ( vec_track_pos.begin(),vec_track_pos.end(),more_pt);
560
561 // clear pos track
562 vector<MdcTrack*> vec_MdcTrack_pos;
563 for( unsigned int i =0;i<vec_track_pos.size();i++){
564 MdcTrack *mdcTrack = new MdcTrack(vec_track_pos[i]->getTrk());
565 vec_MdcTrack_pos.push_back(mdcTrack);
566 if(m_debug>0) cout<<"trackpos : "<<i<<" pt "<<vec_track_pos[i]->getPt3D()<<endl;
567 if(m_debug>0) vec_track_pos[i]->print();
568 }
569 if(m_debug>0) cout<<"step6 : judge pos track list "<<endl;
570 judgeHit(mdcTrackList_pos,vec_MdcTrack_pos);
571 }
572
573 //combine hough itself
574 // if(m_debug>0) cout<<"step7 : combine neg&pos track list "<<endl;
575 if(m_combine){
576 compareHough(mdcTrackList_neg);
577 compareHough(mdcTrackList_pos);
578 }
579 if( mdcTrackList_neg.length()!=0 && mdcTrackList_pos.length()!=0 ) judgeChargeTrack(mdcTrackList_neg,mdcTrackList_pos);
580 //add tracklist
581 mdcTrackList_neg+=mdcTrackList_pos; //neg -> all charge
582
583 //compare pattsf&hough self
584 if(m_combineTracking) nTdsTk = comparePatTsf(mdcTrackList_neg,trackList_tds);
585
586 // store tds
587 if(m_debug>0) cout<<"step8 : store tds "<<endl;
588 if(m_debug>0) cout<<"store tds "<<endl;
589 int tkId = nTdsTk ; //trkid
590 for( unsigned int i =0;i<mdcTrackList_neg.length();i++){
591 if(m_debug>0) cout<<"- charge size i : "<<i<<" "<<mdcTrackList_neg.length()<<endl;
592 int tkStat = 4;//track find by houghspace set stat=4
593 mdcTrackList_neg[i]->storeTrack(tkId, trackList_tds, hitList_tds, tkStat);
594 tkId++;
595 delete mdcTrackList_neg[i];
596 }
597 if(m_debug>0) m_mdcPrintSvc->printRecMdcTrackCol();
598
599 m_timer_all->stop();
600 double t_teventTime = m_timer_all->elapsed();
601 if(m_hist) m_time = t_teventTime;
602
603 if( m_hist ) ntuple_evt->write();
604
605 delete m_map;
606 delete m_map2;
607 if(m_debug>0) cout<<"after delete map "<<endl;
608 for(int ihit=0;ihit<mdcHitCol_neg.size();ihit++){
609 delete mdcHitCol_neg[ihit];
610 delete mdcHitCol_pos[ihit];
611 }
612
613 if(m_debug>0) cout<<"after delete hit"<<endl;
614 //clearMem(mdcTrackList_neg,mdcTrackList_pos);
615 if(m_debug>0) cout<<"end event "<<endl;
616
617 return StatusCode::SUCCESS;
618}
619
620// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
621StatusCode MdcHoughFinder::finalize() {
622 MsgStream log(msgSvc(), name());
623 delete m_bfield ;
624 log << MSG::INFO<< "in finalize()" << endreq;
625
626 return StatusCode::SUCCESS;
627}
628
629int MdcHoughFinder::GetMcInfo(){
630 MsgStream log(msgSvc(), name());
631 StatusCode sc;
632 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
633 if (!mcParticleCol) {
634 log << MSG::ERROR << "Could not find McParticle" << endreq;
635 return StatusCode::FAILURE;
636 }
637
638
639 g_tkParTruth.clear();//clear track param.
640
641 //t_t0Truth=-1;
642 int t_qTruth = 0;
643 int t_pidTruth = -999;
644 int nMcTk=0;
645 McParticleCol::iterator iter_mc = mcParticleCol->begin();
646 for (;iter_mc != mcParticleCol->end(); iter_mc++ ) {
647 t_pidTruth = (*iter_mc)->particleProperty();
648 // cout<<" pid "<<t_pidTruth<<endl;
649 if(!(*iter_mc)->primaryParticle()) continue;
650 //if((m_pid!=0) && (t_pidTruth != m_pid)){ continue; }
651 //t_t0Truth=(*iter_mc)->initialPosition().t();
652 std::string name;
653 int pid = t_pidTruth;
654 if( pid == 0 ) {
655 log << MSG::WARNING << "Wrong particle id" <<endreq;
656 continue;
657 }else{
658 IPartPropSvc* p_PartPropSvc;
659 static const bool CREATEIFNOTTHERE(true);
660 StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
661 if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) {
662 std::cout<< " Could not initialize Particle Properties Service" << std::endl;
663 }
664 HepPDT::ParticleDataTable* p_particleTable = p_PartPropSvc->PDT();
665 if( p_particleTable->particle(pid) ){
666 name = p_particleTable->particle(pid)->name();
667 t_qTruth = p_particleTable->particle(pid)->charge();
668 }else if( p_particleTable->particle(-pid) ){
669 name = "anti " + p_particleTable->particle(-pid)->name();
670 t_qTruth = (-1)*p_particleTable->particle(-pid)->charge();
671 }
672 }
673
674 //helix
675 if(m_debug>1) std::cout<<__FILE__<<" "<<__LINE__<<" new helix with pos "<<(*iter_mc)->initialPosition().v()<<" mom "<<(*iter_mc)->initialFourMomentum().v()<<" q "<<t_qTruth<<std::endl;
676 Helix mchelix = Helix((*iter_mc)->initialPosition().v(), (*iter_mc)->initialFourMomentum().v(), t_qTruth);//cm
677 mchelix.pivot( HepPoint3D(0.,0.,0.) );
678
679 int mcTkId = (*iter_mc)->trackIndex();
680
681 pair<int, const HepVector> p(mcTkId,mchelix.a());
682 g_tkParTruth.insert(p);
683 // exchange to rho theta test negtive charge
684 double rho , theta;
685 if( mchelix.phi0()>M_PI) { rho=-1./fabs(mchelix.radius()); theta = mchelix.phi0()-M_PI;}
686 else {rho =1./fabs(mchelix.radius()) ; theta = mchelix.phi0();}
687 t_d0=mchelix.dr();
688 t_z0=mchelix.dz();
689 t_pt=mchelix.pt();
690 t_p=mchelix.momentum(0.).mag();
691 t_tanl=mchelix.tanl();
692 t_cos= mchelix.cosTheta();
693
694 //m_pzTruth[nMcTk]=mchelix.momentum(0.).z();
695 //m_pidTruth[nMcTk] = t_pidTruth;
696 //m_costaTruth[nMcTk] = mchelix.cosTheta();
697 //m_ptTruth[nMcTk] = mchelix.pt();
698 //m_pTruth[nMcTk] = mchelix.momentum(0.).mag();
699 //m_qTruth[nMcTk] = t_qTruth;
700 //m_drTruth[nMcTk] = mchelix.dr();
701 //m_phi0Truth[nMcTk] = mchelix.phi0();
702 //m_omegaTruth[nMcTk] =1./(mchelix.radius()); //Q
703 //m_dzTruth[nMcTk] = mchelix.dz();
704 //m_tanl_mc[nMcTk] =mchelix.tanl();
705 //m_rho_mc[nMcTk] = rho;
706 //m_theta_mc[nMcTk] = theta;
707
708 if(m_hist){
709 m_pzTruth=mchelix.momentum(0.).z();
710 m_pidTruth = t_pidTruth;
711 m_costaTruth = mchelix.cosTheta();
712 m_ptTruth = mchelix.pt();
713 m_pTruth = mchelix.momentum(0.).mag();
714 m_qTruth = t_qTruth;
715 m_drTruth = mchelix.dr();
716 m_phi0Truth = mchelix.phi0();
717 m_omegaTruth =1./(mchelix.radius()); //Q
718 m_dzTruth = mchelix.dz();
719 m_tanl_mc =mchelix.tanl();
720 m_rho_mc = rho;
721 m_theta_mc = theta;
722 }
723
724 // if(m_debug>0){
725 // std::cout<<"Truth tk "<<nMcTk<<" " <<name<<" pid "<<pid<<" charge "<<t_qTruth << " helix "<< mchelix.a()<<" p "<<mchelix.momentum(0.)<<" p "<<mchelix.momentum(0.).mag()<<" pt "<<mchelix.momentum(0.).perp()<<" pz "<<mchelix.momentum(0.).z()<<std::endl;
726 //
727 // cout<<"ptTruth "<<mchelix.pt()<<endl;
728 // cout<<"tanlTruth"<<mchelix.tanl()<<endl;
729 // cout<<"d0Truth"<<mchelix.dr()<<endl;
730 // }
731 // nMcTk++;
732 }
733 // m_nTrkMC = nMcTk;
734 //if(m_debug>2) m_mdcPrintSvc->printMdcMcHitCol();
735
736 g_firstHit.set(0,0,0);
737 SmartDataPtr<MdcMcHitCol> mdcMcHitCol(eventSvc(),"/Event/MC/MdcMcHitCol");
738 if(mdcMcHitCol){
739 MdcMcHitCol::iterator iter_mchit = mdcMcHitCol->begin();
740 for (;iter_mchit != mdcMcHitCol->end(); iter_mchit++) {
741 if((*iter_mchit)->getTrackIndex()==0) {
742 g_firstHit.setX((*iter_mchit)->getPositionX()/10.);//mm to cm
743 g_firstHit.setY((*iter_mchit)->getPositionY()/10.);
744 g_firstHit.setZ((*iter_mchit)->getPositionZ()/10.);
745 break;
746 }
747 }
748 }else{
749 std::cout<<__FILE__<<" Could not get MdcMcHitCol "<<std::endl;
750 return StatusCode::FAILURE;
751 }
752
753 //add truth info. to HoughHit in list
754 return StatusCode::SUCCESS;
755}
756
757MdcDigiVec MdcHoughFinder::prepareDigi(){
758 // retrieve Mdc digi vector form RawDataProviderSvc
759 uint32_t getDigiFlag = 0;
760 if(m_dropHot || m_combineTracking )getDigiFlag |= MdcRawDataProvider::b_dropHot;
761 if(m_keepBadTdc) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
762 if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
763
764 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
765 return mdcDigiVec;
766}
767
768void MdcHoughFinder::dumpHoughHitList(const HoughHitList& houghhitlist){
769
770 int num_first_half=0;
771 vector<double> vtemp,utemp;
772 double k,b,theta,rho,x_cross,y_cross;
773 std::vector<HoughHit>::const_iterator iter = houghhitlist.getList().begin();
774 int nsig_axial=0;
775 int exit_multiturn=0;
776 for(int iHit=0;iter!= houghhitlist.getList().end(); iter++,iHit++){
777 const HoughHit h = (*iter);
778 if( m_debug>0 ) h.printTruth();
779 //if( h.getCirList()!=0 ) continue;
780 m_layer[iHit] = h.getLayerId();
781 m_cell[iHit] = h.getWireId();
782 m_x[iHit] = h.getMidX();
783 m_y[iHit] = h.getMidY();
784 m_u[iHit] = h.getU();
785 m_v[iHit] = h.getV();
786 m_drift[iHit] = h.getDriftDist();
787 m_r[iHit] = h.getR();
788 m_slant[iHit] = h.getSlayerType();
789 m_uTruth[iHit] = h.getUTruth();
790 m_vTruth[iHit] = h.getVTruth();
791 m_xTruth[iHit] = h.getPointTruth().x();
792 m_yTruth[iHit] = h.getPointTruth().y();
793 m_driftTruth[iHit] = h.getDriftDistTruth();
794 m_rTruth[iHit] = h.getRTruth();
795 int type =0; // type : 0 1 2
796 type = h.digi()->getTrackIndex(); // signal noise : 0 1
797 if ( type <0 ) type=1;
798 if ( type ==0 && h.getCirList()>0 ) type = 2; // signal && 2nd half : 2
799 m_type[iHit] = type;
800 if ( h.digi()->getTrackIndex()>0 && h.getCirList()>1 ) exit_multiturn = 1; // signal && 2nd half : 2
801
802 num_first_half++;
803 if( h.getCirList()==0 && h.digi()->getTrackIndex()>=0 && h.getStyle()!=-999 && h.getSlayerType()==0 ) nsig_axial++;
804
805 if ( h.getSlayerType()==0 && h.getCirList()==0 && h.digi()->getTrackIndex()>=0 && h.getStyle()!=-999 && utemp.size()<10 ) //axial,1st,signal,size<10
806 {
807 utemp.push_back(h.getUTruth());
808 vtemp.push_back(h.getVTruth());
809 }
810 m_nSig_Axial=nsig_axial;
811 m_exit_multiturn = exit_multiturn;
812
813
814 Leastfit(utemp,vtemp,k,b);
815 //calcu truth
816 //k,b from truth
817 x_cross = -b/(k+1/k);
818 y_cross = b/(1+k*k);
819 rho=sqrt(x_cross*x_cross+y_cross*y_cross);
820 theta=atan2(y_cross,x_cross);
821 if(theta<0) {
822 theta=theta+M_PI;
823 rho=-rho;
824 }
825
826 m_rho_line =rho;
827 m_theta_line=theta;
828
829
830 }
831 // m_nHit = houghhitlist.nHit();
832 m_nHit = num_first_half;
833}
834
835void MdcHoughFinder::dumpHoughMap(const HoughMap& houghmap){
836 //m_MapMax = houghmap.MAX;
837
838 // //dump max bin
839 // int BINX=m_nBinTheta;
840 // int BINY=m_nBinRho;
841 // m_xybinNum=BINX*BINY;
842 // double **s = houghmap.getS();
843 // double max = 0;
844 // for(int i=0;i<BINX;i++){
845 // for(int j=0;j<BINY;j++){
846 // if( s[i][j]>max ) max=s[i][j];
847 // m_xybinS[BINY*i+j] = s[i][j];
848 // }
849 // }
850 // m_xybinMax = max;
851
852 //method -- narrow
853 // int BINX=m_n2;
854 // int BINY=m_n2;
855 // m_xybinNum=BINX*BINY;
856 // double **s = houghmap.getS2();
857
858 // m_theta_left = houghmap.Theta_left;
859 // m_theta_right= houghmap.Theta_right;
860 // m_rho_down= houghmap.Rho_down;
861 // m_rho_up= houghmap.Rho_up;
862 // m_theta= houghmap.Theta;
863 // m_rho= houghmap.Rho;
864 // m_height = houghmap.Height;
865
866 // m_aver= houghmap.Aver;
867 // m_sigma= houghmap.Sigma;
868
869 // dump tracklist
870 m_nMapPk=houghmap.getPeakNumber();
871 //for(int iPk=0;iPk< houghmap.getPeakNumber(); iPk++){
872 // m_PeakOrder[iPk] = houghmap.getPeak(iPk).getPeakNum();
873 // m_PeakRho[iPk] = houghmap.getPeak(iPk).getRho();
874 //cout<<"rho: "<<houghmap.getPeak(iPk).getRho()<<endl;
875 // m_PeakTheta[iPk]= houghmap.getPeak(iPk).getTheta();
876 //cout<<"theta : "<<houghmap.getPeak(iPk).getTheta()<<endl;
877 //m_PeakHeight[iPk]= houghmap.getPeak(iPk).peakHeight();
878 // m_PeakHit[iPk]= houghmap.getPeak(iPk).getHoughHitList().size();
879 // m_PeakHitA[iPk] = houghmap.getPeak(iPk).getHitNumA(6);
880 //}
881 // m_nMapTrk=houghmap.getTrackNumber();
882 //cout<<"npk ntrk "<<m_nMapPk<<" "<<m_nMapTrk<<endl;
883 //for(int iTk=0;iTk< houghmap.getTrackNumber(); iTk++){
884 // m_TrackRho[iTk] = houghmap.getTrack(iTk).getRho();
885 // m_TrackTheta[iTk]= houghmap.getTrack(iTk).getTheta();
886 // m_TrackHitA[iTk] = houghmap.getTrack(iTk).getHitNumA(6);
887 //}
888 //maxlayer_slant
889 // m_nMaxLayerSlant = houghmap.get_maxlayer_slant().size();
890 // //cout<<"m_nMaxLayerSlant "<<m_nMaxLayerSlant<<endl;
891 // for(int inmax=0;inmax< houghmap.get_maxlayer_slant().size(); inmax++){
892 // m_MaxLayerSlant[inmax]=houghmap.get_maxlayer_slant().at(inmax);
893 // m_MaxLayer[inmax]=houghmap.get_maxlayer();
894 // }
895 // //nomaxlayer_slant
896 // m_nNoMaxLayerSlant = houghmap.get_nomaxlayer_slant().size();
897 // //cout<<"m_nNoMaxLayerSlant "<<m_nNoMaxLayerSlant<<endl;
898 // for(int inmax=0;inmax< houghmap.get_nomaxlayer_slant().size(); inmax++){
899 // m_NoMaxLayerSlant[inmax]=houghmap.get_nomaxlayer_slant().at(inmax);
900 // m_NoMaxLayerid[inmax]=houghmap.get_nomaxlayerid().at(inmax);
901 // }
902 //
903 // std::vector<HoughHit>::const_iterator iter = houghmap.getHitList().getList().begin();
904 // for(int iHit=0;iter!= houghmap.getHitList().getList().end(); iter++,iHit++){
905 // const HoughHit h = (*iter);
906 // //if( h.getSlayerType()!=0 || h.getCirList()!=0 ) continue;
907 // if( h.getSlayerType()!=0 ) continue;
908 // m_deltaD[iHit] = h.getDeltaD();
909 // m_flt[iHit] = h.getFltLen();
910 // }
911}
912void MdcHoughFinder::diffMap(const HoughMap& houghmap,const HoughMap& houghmap2){
913 //map 1
914 m_nMap1Pk=houghmap.getPeakNumber();
915 for(int iPk=0;iPk< houghmap.getPeakNumber(); iPk++){
916 m_PkRho1[iPk] = houghmap.getPeak(iPk).getRho();
917 m_PkTheta1[iPk]= houghmap.getPeak(iPk).getTheta();
918 }
919 m_nMap1Tk=houghmap.getTrackNumber();
920 for(int iTk=0;iTk< houghmap.getTrackNumber(); iTk++){
921 m_TkRho1[iTk] = houghmap.getTrack(iTk).getRho();
922 m_TkTheta1[iTk]= houghmap.getTrack(iTk).getTheta();
923 }
924 //map 2
925 m_nMap2Pk=houghmap2.getPeakNumber();
926 for(int iPk=0;iPk< houghmap2.getPeakNumber(); iPk++){
927 m_PkRho2[iPk] = houghmap2.getPeak(iPk).getRho();
928 m_PkTheta2[iPk]= houghmap2.getPeak(iPk).getTheta();
929 }
930 m_nMap2Tk=houghmap2.getTrackNumber();
931 for(int iTk=0;iTk< houghmap2.getTrackNumber(); iTk++){
932 m_TkRho2[iTk] = houghmap2.getTrack(iTk).getRho();
933 m_TkTheta2[iTk]= houghmap2.getTrack(iTk).getTheta();
934 }
935}
936
937void MdcHoughFinder::Leastfit(vector<double> u,vector<double> v,double& k ,double& b){
938 int N = u.size();
939 if( N <= 2 ) return;
940 double *x=new double[N];
941 double *y=new double[N];
942 for(int i=0;i<N;i++){
943 x[i]=u[i];
944 y[i]=v[i];
945 }
946
947 TF1 *fpol1=new TF1("fpol1","pol1",-50,50);
948 TGraph *tg=new TGraph(N,x,y);
949 tg->Fit("fpol1","QN");
950 double ktemp =fpol1->GetParameter(1);
951 double btemp =fpol1->GetParameter(0);
952 k = ktemp;
953 b = btemp;
954 delete []x;
955 delete []y;
956 delete fpol1;
957 delete tg;
958}
959
960bool more_pt(const HoughTrack* tracka,const HoughTrack* trackb){
961 return tracka->getPt3D() > trackb->getPt3D();
962}
963
964int MdcHoughFinder::storeTracks(RecMdcTrackCol*& trackList_tds ,RecMdcHitCol*& hitList_tds,vector<RecMdcTrack*>& vec_trackPatTds){
965 MsgStream log(msgSvc(), name());
966 StatusCode sc;
967 // tds
968 DataObject *aTrackCol;
969 DataObject *aRecHitCol;
970 SmartIF<IDataManagerSvc> dataManSvc(eventSvc());
971 if(!m_combineTracking){
972 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
973 if(aTrackCol != NULL) {
974 dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
975 eventSvc()->unregisterObject("/Event/Recon/RecMdcTrackCol");
976 }
977 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
978 if(aRecHitCol != NULL) {
979 dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
980 eventSvc()->unregisterObject("/Event/Recon/RecMdcHitCol");
981 }
982 }
983
984 //RecMdcTrackCol* trackList_tds;
985 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
986 if (aTrackCol) {
987 trackList_tds = dynamic_cast<RecMdcTrackCol*> (aTrackCol);
988 }else{
989 trackList_tds = new RecMdcTrackCol;
990 StatusCode sc = eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trackList_tds);
991 if(!sc.isSuccess()) {
992 log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
993 return StatusCode::FAILURE;
994 }
995 }
996 //RecMdcHitCol* hitList_tds;
997 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
998 if (aRecHitCol) {
999 hitList_tds = dynamic_cast<RecMdcHitCol*> (aRecHitCol);
1000 }else{
1001 hitList_tds = new RecMdcHitCol;
1002 StatusCode sc = eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitList_tds);
1003 if(!sc.isSuccess()) {
1004 log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
1005 return StatusCode::FAILURE;
1006 }
1007 }
1008 //remove bad quality or low pt track(s)
1009 if(m_removeBadTrack && m_combineTracking) {
1010 if( m_debug>0 ) cout<<"PATTSF collect "<<trackList_tds->size()<<" track. "<<endl;
1011 if(trackList_tds->size()!=0){
1012 RecMdcTrackCol::iterator iter_pat = trackList_tds->begin();
1013 for(;iter_pat!=trackList_tds->end();iter_pat++){
1014 double d0=(*iter_pat)->helix(0);
1015 double kap = (*iter_pat)->helix(2);
1016 double pt = 0.00001;
1017 if(fabs(kap)>0.00001) pt = fabs(1./kap);
1018 double dz=(*iter_pat)->helix(3);
1019 double chi2=(*iter_pat)->chi2();
1020 if( m_debug>0) cout<<"d0: "<<d0<<" z0: "<<dz<<" pt "<<pt<<" chi2 "<<chi2<<endl;
1021 //if( (fabs(d0)>m_dropTrkDrCut || fabs(dz)>m_dropTrkDzCut || chi2>m_dropTrkChi2Cut) && pt<=m_dropTrkPtCut)
1022 if( pt<=m_dropTrkPtCut ){
1023 vec_trackPatTds.push_back(*iter_pat);
1024 //remove hits on track
1025 HitRefVec rechits = (*iter_pat)->getVecHits();
1026 HitRefVec::iterator hotIter = rechits.begin();
1027 while (hotIter!=rechits.end()) {
1028 Identifier id = (*hotIter)->getMdcId();
1029 int layer = MdcID::layer(id);
1030 int wire = MdcID::wire(id);
1031 if(m_debug>0) cout <<"remove hit " << (*hotIter)->getId()<<"("<<layer<<","<<wire<<")"<<endl;
1032 hitList_tds->remove(*hotIter);
1033 hotIter++;
1034 }
1035 //if( (fabs(d0)>3*m_dropTrkDrCut || fabs(dz)>3*m_dropTrkDzCut) ){ // drop pattsf fate track
1036 // trackList_tds->remove(*iter_pat);
1037 // iter_pat--;
1038 //}
1039 }
1040 if( m_debug>0 ) cout<<"after delete trackcol size : "<<trackList_tds->size()<<endl;
1041 }
1042 } else {
1043 if(m_debug>0) cout<<" PATTSF find 0 track. "<<endl;
1044 }
1045 }//end of remove bad quality high pt track
1046
1047 int nTdsTk=trackList_tds->size();
1048 return nTdsTk;
1049}//end of stroeTracks
1050
1051void MdcHoughFinder::clearMem(MdcTrackList& list1,MdcTrackList& list2){
1052 cout<<"length in clearMem "<<list1.length()<<" "<<list2.length()<<endl;
1053 cout<<"in list "<<vectrk_for_clean.size()<<endl;
1054
1055 for(unsigned int i=0;i<vectrk_for_clean.size();i++){
1056 int sametrk=0;
1057 for(unsigned int j=0;j<list1.length();j++){
1058 cout<<"-i j "<<i<<" "<<j<<" "<<vectrk_for_clean[i]<<" "<<&(list1[j]->track())<<endl;
1059 if(vectrk_for_clean[i] == &(list1[j]->track()) ) {sametrk=1;cout<<"not delete new trk "<<endl;}
1060 //delete list1[j];
1061 }
1062 for(unsigned int k=0;k<list2.length();k++){
1063 cout<<"+i k "<<i<<" "<<k<<" "<<vectrk_for_clean[i]<<" "<<&(list2[k]->track())<<endl;
1064 if(vectrk_for_clean[i] == &(list2[k]->track()) ) {sametrk=1;cout<<"not delete new trk "<<endl;}
1065 //delete list2[k];
1066 }
1067 if( sametrk ==0 ) {
1068 cout<<"delete i "<<i<<endl;
1069 delete vectrk_for_clean[i];
1070 vectrk_for_clean[i] = NULL;
1071 }
1072 }
1073 vectrk_for_clean.clear();
1074}
1075
1076int MdcHoughFinder::judgeChargeTrack(MdcTrackList& list1 ,MdcTrackList& list2){
1077 if(m_debug>0) cout<<"in judgeChargeTrack"<<endl;
1078 if(m_debug>0) cout<<"ntrack length neg : "<<list1.length()<<endl;
1079 if(m_debug>0) cout<<"ntrack length pos : "<<list2.length()<<endl;
1080 for(int itrack1=0; itrack1<list1.nTrack(); itrack1++){
1081 MdcTrack* track_1 = list1[itrack1];
1082 TrkExchangePar par1 = track_1->track().fitResult()->helix(0.);
1083 double d0_1 = par1.d0();
1084 double phi0_1 = par1.phi0();
1085 double omega_1 = par1.omega();
1086 double cr=fabs(1./par1.omega());
1087 double cx= sin(par1.phi0()) *(par1.d0()+1/par1.omega()) * -1.; //z axis verse,x_babar = -x_belle
1088 double cy= -1.*cos(par1.phi0()) *(par1.d0()+1/par1.omega()) * -1;//???
1089 double z0_1 = par1.z0();
1090 for(int itrack2=0; itrack2<list2.nTrack(); itrack2++){
1091 MdcTrack* track_2 = list2[itrack2];
1092 TrkExchangePar par2 = track_2->track().fitResult()->helix(0.);
1093 double d0_2 = par2.d0();
1094 double phi0_2 = par2.phi0();
1095 double omega_2 = par2.omega();
1096 double cR=fabs(1./par2.omega());
1097 double cX= sin(par2.phi0()) *(par2.d0()+1/par2.omega()) * -1.; //z axis verse,x_babar = -x_belle
1098 double cY= -1.*cos(par2.phi0()) *(par2.d0()+1/par2.omega()) * -1;//???
1099 double z0_2= par2.z0();
1100
1101 double bestDiff = 1.0e+20;
1102 if(fabs(cr)*(1.-0.25) <= fabs(cR) && fabs(cR) <= fabs(cr)*(1.+0.25)){
1103 if(fabs(cx-cX) <= fabs(cr)*(0.25)&& fabs(cy-cY) <= fabs(cr)*(0.25) ){
1104 double diff = fabs((fabs(cr)-fabs(cR))*(cx-cX)*(cy-cY));
1105 if(diff < bestDiff){
1106 bestDiff = diff;
1107 }
1108 }
1109 }
1110
1111 if(bestDiff != 1.0e20) {
1112 if(m_debug>0) {
1113 cout<<"There is ambiguous +- charge in the track list "<<endl;
1114 cout<<"z0_1 : "<<z0_1<<" z0_2 : "<<z0_2<<endl;
1115 }
1116 //iter_hough2 = vec_track.erase(iter_hough2); iter_hough2--;
1117 if( fabs(z0_1) >= fabs(z0_2) ) { list1.remove( track_1 ); itrack1--; break;}
1118 if( fabs(z0_1) < fabs(z0_2) ) { list2.remove( track_2 ); itrack2--; break;}
1119 }
1120 if(bestDiff == 1.0e20) { if(m_debug>0) cout<<" no (+-) track in hough"<<endl; }
1121 }
1122 }
1123}
1124
1125void MdcHoughFinder::compareHough(MdcTrackList& mdcTrackList){
1126 if(m_debug>0) cout<<"in compareHough "<<endl;
1127 if(m_debug>0) cout<<"ntrack : "<<mdcTrackList.length()<<endl;
1128 for(int itrack=0; itrack<mdcTrackList.length(); itrack++){
1129 MdcTrack* track = mdcTrackList[itrack];
1130 TrkExchangePar par = track->track().fitResult()->helix(0.);
1131 int nhit = track->track().hots()->nHit();
1132 double pt = (1./par.omega())/333.567;
1133 if(m_debug>0) cout<<"i Track : "<<itrack<<" nHit: "<<nhit<<" pt: "<<pt<<endl;
1134 }
1135 // vector<HoughTrack>::iterator iter_hough = vec_track.begin();
1136 // vector<HoughTrack>::iterator iter_hough2 = vec_track.begin()+1;
1137 for(int itrack1=0; itrack1<mdcTrackList.length(); itrack1++){
1138 MdcTrack* track_1 = mdcTrackList[itrack1];
1139 TrkExchangePar par1 = track_1->track().fitResult()->helix(0.);
1140 int nhit1 = track_1->track().hots()->nHit();
1141 double d0_1 = par1.d0();
1142 double phi0_1 = par1.phi0();
1143 double omega_1 = par1.omega();
1144 double z0_1= par1.z0();
1145 double cr=fabs(1./par1.omega());
1146 double cx= sin(par1.phi0()) *(par1.d0()+1/par1.omega()) * -1.; //z axis verse,x_babar = -x_belle
1147 double cy= -1.*cos(par1.phi0()) *(par1.d0()+1/par1.omega()) * -1;//???
1148 if(m_debug>0) cout<<"circle 1 : "<<cr<<","<<cx<<","<<cy<<endl;
1149
1150 for(int itrack2=itrack1+1; itrack2<mdcTrackList.length(); itrack2++){
1151 MdcTrack* track_2 = mdcTrackList[itrack2];
1152 TrkExchangePar par2 = track_2->track().fitResult()->helix(0.);
1153 int nhit2 = track_2->track().hots()->nHit();
1154 double d0_2 = par2.d0();
1155 double phi0_2 = par2.phi0();
1156 double omega_2 = par2.omega();
1157 double z0_2= par2.z0();
1158 double cR=fabs(1./par2.omega());
1159 double cX= sin(par2.phi0()) *(par2.d0()+1/par2.omega()) * -1.; //z axis verse,x_babar = -x_belle
1160 double cY= -1.*cos(par2.phi0()) *(par2.d0()+1/par2.omega()) * -1;//???
1161 if(m_debug>0) cout<<"circle 2 : "<<cR<<","<<cX<<","<<cY<<endl;
1162
1163 double bestDiff = 1.0e+20;
1164 if(fabs(cr)*(1.-0.25) <= fabs(cR) && fabs(cR) <= fabs(cr)*(1.+0.25)){
1165 if(fabs(cx-cX) <= fabs(cr)*(0.25)&& fabs(cy-cY) <= fabs(cr)*(0.25) ){
1166 double diff = fabs((fabs(cr)-fabs(cR))*(cx-cX)*(cy-cY));
1167 if(diff < bestDiff){
1168 bestDiff = diff;
1169 }
1170 }
1171 }
1172
1173 double zdiff = z0_1-z0_2;
1174 if(bestDiff != 1.0e20 && fabs(zdiff)<25.){
1175 if(m_debug>0) cout<<"z0_1 : "<<z0_1<<" z0_2 : "<<z0_2<<endl;
1176 //if( fabs(z0_1)>fabs(z0_2) && fabs(z0_1)>m_z0Cut_compareHough)
1177 if( nhit1<nhit2 && (fabs(z0_1)>fabs(z0_2) && fabs(z0_1)>m_z0Cut_compareHough) ){ //FIXME
1178 if(m_debug>0) cout<<"remove track1 "<<1./par1.omega()/333.567<<endl;
1179 //remove 1
1180 mdcTrackList.remove( track_1 );
1181 itrack1--;
1182 break;
1183 }
1184 else{
1185 //remove 2
1186 if(m_debug>0) cout<<"remove track2 "<<1./par2.omega()/333.567<<endl;
1187 //remove 1
1188 mdcTrackList.remove( track_2 );
1189 itrack2--;
1190 }
1191 }
1192 if(bestDiff == 1.0e20) { if(m_debug>0) cout<<" no same track in hough"<<endl; }
1193 }
1194 }
1195}
1196int MdcHoughFinder::comparePatTsf(MdcTrackList& mdcTrackList, RecMdcTrackCol* trackList_tds){
1197 // vector<HoughTrack>::iterator iter_hough = vec_track.begin();
1198 // int itrk_hough=0;
1199 // for(;iter_hough != vec_track.end();iter_hough++)
1200 // if(m_debug>0) cout<<"being hough trk "<<itrk_hough<<endl;
1201 // itrk_hough++;
1202 // double cr= (*iter_hough).getCirR();
1203 // double cx= (*iter_hough).getCirX();
1204 // double cy= (*iter_hough).getCirY();
1205 for(int itrack1=0; itrack1<mdcTrackList.length(); itrack1++){
1206 MdcTrack* track_1 = mdcTrackList[itrack1];
1207 TrkExchangePar par1 = track_1->track().fitResult()->helix(0.);
1208 int nhit1 = track_1->track().hots()->nHit();
1209 double d0_1 = par1.d0();
1210 double phi0_1 = par1.phi0();
1211 double omega_1 = par1.omega();
1212 double z0_1= par1.z0();
1213 double cr=fabs(1./par1.omega());
1214 double cx= sin(par1.phi0()) *(par1.d0()+1/par1.omega()) * -1.; //z axis verse,x_babar = -x_belle
1215 double cy= -1.*cos(par1.phi0()) *(par1.d0()+1/par1.omega()) * -1;//???
1216
1217 //track2
1218 double bestDiff = 1.0e+20;
1219 double ptDiff = 0.0;
1220 double cR, cX, cY;
1221 vector<double> a0,a1,a2,a3,a4;
1222 RecMdcTrackCol::iterator iterBest;
1223 RecMdcTrackCol::iterator iteritrk = trackList_tds->begin();
1224 int itrk=0;
1225 for(;iteritrk!=trackList_tds->end();iteritrk++){
1226 if(m_debug>0) cout<<"being pattsf trk "<<itrk<<endl;
1227 double pt=(*iteritrk)->pxy();
1228 a0.push_back( (*iteritrk)->helix(0) );
1229 a1.push_back( (*iteritrk)->helix(1) );
1230 a2.push_back( (*iteritrk)->helix(2) );
1231 a3.push_back( (*iteritrk)->helix(3) );
1232 a4.push_back( (*iteritrk)->helix(4) );
1233 cR=(-333.567)/a2[itrk];
1234 cX=(cR+a0[itrk])*cos(a1[itrk]);
1235 cY=(cR+a0[itrk])*sin(a1[itrk]);
1236 if(m_debug>0)
1237 {
1238 cout<<" compare PATTSF and HOUGH "<<endl;
1239 cout<<" fabs(cX) "<< fabs(cX)<< "fabs(cx) "<<fabs(cx)<<endl;
1240 cout<<" fabs(cY) "<< fabs(cY)<< "fabs(cy) "<<fabs(cy)<<endl;
1241 cout<<" fabs(cR) "<< fabs(cR)<< "fabs(cr) "<<fabs(cr)<<endl;
1242 cout<<" fabs(cx-cX) "<< fabs(cx-cX)<< "fabs(cy-cY) "<<fabs(cy-cY)<<" fabs(cr-cR) "<< fabs(cr-cR)<<endl;
1243 }
1244
1245 if(fabs(cr)*(1.-0.25) <= fabs(cR) && fabs(cR) <= fabs(cr)*(1.+0.25)){
1246 if(fabs(cx-cX) <= fabs(cr)*(0.25)&& fabs(cy-cY) <= fabs(cr)*(0.25) ){
1247 double diff = fabs((fabs(cr)-fabs(cR))*(cx-cX)*(cy-cY));
1248 if(m_debug>0) cout<<" diff "<<diff<<" pt "<<pt<<endl;
1249 if( fabs(pt) > ptDiff){
1250 ptDiff= pt;
1251 iterBest=iteritrk;
1252 bestDiff = diff;
1253 }
1254 //if(diff < bestDiff){
1255 // bestDiff = diff;
1256 // // bestIndex = itrk;
1257 // iterBest=iteritrk;
1258 //}
1259 }
1260 }
1261 itrk++;
1262 }
1263 if(bestDiff != 1.0e20) {
1264 if(m_debug>0) cout<<" same track pattsf & hough , then compare nhit. "<<endl;
1265 double d0_pattsf =(*iterBest)->helix(0);
1266 double z0_pattsf =(*iterBest)->helix(3);
1267 double d0_hough = d0_1;
1268 double z0_hough = z0_1;
1269 int use_pattsf(1),use_hough(1);
1270 if( fabs(d0_pattsf)>m_dropTrkDrCut || fabs(z0_pattsf)>m_dropTrkDzCut ) use_pattsf=0;
1271 if( fabs(d0_hough)>m_dropTrkDrCut || fabs(z0_hough)>m_dropTrkDzCut ) use_hough=0;
1272 if( use_pattsf==0 && use_hough==1 ) {
1273 trackList_tds->remove(*iterBest);
1274 if(m_debug>0) cout<<" use houghTrack, vertex pattsf wrong"<<endl;
1275 }
1276 if( use_pattsf==1 && use_hough==0 ) {
1277 mdcTrackList.remove( track_1 );
1278 itrack1--;
1279 if(m_debug>0) cout<<" use pattsfTrack, vertex hough wrong"<<endl;
1280 }
1281 if( (use_pattsf==0 && use_hough==0) || (use_pattsf==1 && use_hough==1) ){
1282 int nhit_pattsf=( (*iterBest)->ndof()+5 );
1283 int nhit_hough = nhit1;
1284 if(m_debug>0) cout<<" nhit "<<nhit_pattsf<<" "<<nhit_hough<<endl;
1285 if( nhit_pattsf>nhit_hough ) {
1286 mdcTrackList.remove( track_1 );
1287 itrack1--;
1288 if(m_debug>0) cout<<" use pattsfTrack "<<endl;
1289 }
1290 else{
1291 trackList_tds->remove(*iterBest);
1292 if(m_debug>0) cout<<" use houghTrack "<<endl;
1293 }
1294 }
1295 }
1296 // if(m_debug>0) cout<<" d0"<<d0_pattsf<<" "<<d0_hough<<endl;
1297 // if(m_debug>0) cout<<" z0"<<z0_pattsf<<" "<<z0_hough<<endl;
1298 // if( fabs(d0_pattsf) >= fabs(d0_hough) && fabs(z0_pattsf) >= fabs(z0_hough) ) {
1299 // trackList_tds->remove(*iterBest);
1300 // if(m_debug>0) cout<<" use houghTrack "<<endl;
1301 // }
1302 // else{
1303 // iter_hough = vec_track.erase(iter_hough); iter_hough--;
1304 // if(m_debug>0) cout<<" use pattsfTrack "<<endl;
1305 // }
1306
1307 if(bestDiff == 1.0e20) { if(m_debug>0) cout<<" no same track in pattsf"<<endl; }
1308 //if(bestDiff != 1.0e20) { iter_hough = vec_track.erase(iter_hough); iter_hough--;}
1309 //if(bestDiff == 1.0e20) { if(m_debug>0) cout<<" no merge "<<endl; }
1310 }
1311 //cout<<"size after combine "<<trackList_tds->size()<<endl;
1312 if(trackList_tds->size()!=0){
1313 int digi=0;
1314 RecMdcTrackCol::iterator iter_pat = trackList_tds->begin();
1315 for(;iter_pat!=trackList_tds->end();iter_pat++,digi++){
1316 (*iter_pat)->setTrackId(digi);
1317 // cout<<"digi "<<(*iter_pat)->trackId()<<endl;
1318 }
1319 }
1320 return trackList_tds->size();
1321}
1322//int MdcHoughFinder::judgeHit(MdcTrackList& list,const vector<HoughTrack>& trackCol)
1323int MdcHoughFinder::judgeHit(MdcTrackList& list,vector<MdcTrack*>& trackCol){
1324 if(m_debug>0) cout<<"in judgHit: "<<endl;
1325 for( unsigned int i =0;i<trackCol.size();i++){
1326 //MdcTrack *mdcTrack = new MdcTrack(trackCol[i].getTrk());
1327 //list.append(mdcTrack);
1328 list.append(trackCol[i]);
1329 }
1330 int nDeleted = list.arbitrateHits();
1331 return nDeleted;
1332}
1333StatusCode MdcHoughFinder::bookTuple(){
1334 MsgStream log(msgSvc(), name());
1335 NTuplePtr nt1(ntupleSvc(), "MdcHoughFinder/evt");
1336 if ( nt1 ){
1337 ntuple_evt = nt1;
1338 } else {
1339 ntuple_evt = ntupleSvc()->book("MdcHoughFinder/evt", CLID_ColumnWiseTuple, "evt");
1340 if(ntuple_evt){
1341 ntuple_evt->addItem("eventNum", m_eventNum);
1342 ntuple_evt->addItem("runNum", m_runNum);
1343
1344 ntuple_evt->addItem("nHit", m_nHit,0, 6796);
1345
1346 //ntuple_evt->addItem("layer", m_nHit, m_layer);
1347 //ntuple_evt->addItem("cell", m_nHit, m_cell);
1348 //ntuple_evt->addItem("x", m_nHit, m_x);
1349 //ntuple_evt->addItem("y", m_nHit, m_y);
1350 //ntuple_evt->addItem("u", m_nHit, m_u);
1351 //ntuple_evt->addItem("v", m_nHit, m_v);
1352 //ntuple_evt->addItem("drift", m_nHit, m_drift);
1353 //ntuple_evt->addItem("r", m_nHit, m_r);
1354
1355 //ntuple_evt->addItem("uTruth", m_nHit, m_uTruth);
1356 //ntuple_evt->addItem("vTruth", m_nHit, m_vTruth);
1357 //ntuple_evt->addItem("xTruth", m_nHit, m_xTruth);
1358 //ntuple_evt->addItem("yTruth", m_nHit, m_yTruth);
1359 //ntuple_evt->addItem("driftTruth", m_nHit, m_driftTruth);
1360 //ntuple_evt->addItem("rTruth", m_nHit, m_rTruth);
1361 //ntuple_evt->addItem("type", m_nHit, m_type);
1362
1363 //ntuple_evt->addItem("deltaD", m_nHit, m_deltaD);
1364 //ntuple_evt->addItem("flt", m_nHit, m_flt);
1365 //ntuple_evt->addItem("slant", m_nHit, m_slant);
1366
1367 //ntuple_evt->addItem("nSigAxial", m_nSig_Axial, 0, 6796);
1368 //ntuple_evt->addItem("xybinNum", m_xybinNum, 0,10000000);
1369 //ntuple_evt->addItem("xybinMax", m_xybinMax);
1370 //ntuple_evt->addItem("xybinNL", m_xybinNum, m_xybinNL);
1371 //ntuple_evt->addItem("xybinRL", m_xybinNum, m_xybinRL);
1372 //ntuple_evt->addItem("xybinS", m_xybinNum, m_xybinS);
1373 //ntuple_evt->addItem("theta_left", m_theta_left);
1374 //ntuple_evt->addItem("theta_right", m_theta_right);
1375 //ntuple_evt->addItem("rho_down", m_rho_down);
1376 //ntuple_evt->addItem("rho_up", m_rho_up);
1377 ////by houghmap maxbin value 0~M_PI
1378 //ntuple_evt->addItem("theta", m_theta);
1379 //ntuple_evt->addItem("rho", m_rho);
1380 //ntuple_evt->addItem("height", m_height);
1381 //ntuple_evt->addItem("aver", m_aver);
1382 //ntuple_evt->addItem("sigma", m_sigma);
1383
1384 ////multiturn
1385 //ntuple_evt->addItem("multiturn", m_exit_multiturn);
1386
1387 ////diff two charge
1388 //ntuple_evt->addItem("nMap1pk", m_nMap1Pk, 0,100);
1389 //ntuple_evt->addItem("nMap1tk", m_nMap1Tk, 0,100);
1390 //ntuple_evt->addItem("nMap2pk", m_nMap2Pk, 0,100);
1391 //ntuple_evt->addItem("nMap2tk", m_nMap2Tk, 0,100);
1392 //ntuple_evt->addItem("map1_pkrho", m_nMap1Pk, m_PkRho1);
1393 //ntuple_evt->addItem("map1_pktheta", m_nMap1Pk, m_PkTheta1);
1394 //ntuple_evt->addItem("map1_tkrho", m_nMap1Tk, m_TkRho1);
1395 //ntuple_evt->addItem("map1_tktheta", m_nMap1Tk, m_TkTheta1);
1396 //ntuple_evt->addItem("map2_pkrho", m_nMap2Pk, m_PkRho2);
1397 //ntuple_evt->addItem("map2_pktheta", m_nMap2Pk, m_PkTheta2);
1398 //ntuple_evt->addItem("map2_tkrho", m_nMap2Tk, m_TkRho2);
1399 //ntuple_evt->addItem("map2_tktheta", m_nMap2Tk, m_TkTheta2);
1400
1401 ////MC truth
1402 //ntuple_evt->addItem("thetaline", m_theta_line);
1403 //ntuple_evt->addItem("rholine", m_rho_line);
1404
1405 // ntuple_evt->addItem("nmax", m_nMaxLayerSlant,0,200);
1406 // ntuple_evt->addItem("maxslant", m_nMaxLayerSlant, m_MaxLayerSlant);
1407 // ntuple_evt->addItem("maxslantlayer", m_nMaxLayerSlant, m_MaxLayer);
1408 // ntuple_evt->addItem("nnomax", m_nNoMaxLayerSlant,0,1000);
1409 // ntuple_evt->addItem("nomaxslant", m_nNoMaxLayerSlant, m_NoMaxLayerSlant);
1410 // ntuple_evt->addItem("nomaxlayerid", m_nNoMaxLayerSlant, m_NoMaxLayerid);
1411
1412 // ntuple_evt->addItem("MapMax", m_MapMax );
1413 ntuple_evt->addItem("nMapPk", m_nMapPk, 0,100);
1414 // ntuple_evt->addItem("nPeakNum", m_nMapPk, m_PeakOrder);
1415 // ntuple_evt->addItem("nPeakRho", m_nMapPk, m_PeakRho);
1416 // ntuple_evt->addItem("nPeakTheta", m_nMapPk, m_PeakTheta);
1417 // ntuple_evt->addItem("nPeakHeight", m_nMapPk, m_PeakHeight);
1418 // ntuple_evt->addItem("nPeakHit", m_nMapPk, m_PeakHit);
1419 // ntuple_evt->addItem("nPeakHitA", m_nMapPk, m_PeakHitA);
1420 // ntuple_evt->addItem("nMapTrk", m_nMapTrk, 0,1000);
1421 // ntuple_evt->addItem("nTrkRho", m_nMapTrk, m_TrackRho);
1422 // ntuple_evt->addItem("nTrkTheta", m_nMapTrk, m_TrackTheta);
1423 // ntuple_evt->addItem("nTrkHitA", m_nMapTrk, m_TrackHitA);
1424
1425 //rec - charge
1426 //global 2d fit
1427 // ntuple_evt->addItem("nTrk2D_neg", m_nTrk2D_neg, 0,100);
1428 // ntuple_evt->addItem("pt2D_neg", m_nTrk2D_neg, m_pt2D_neg);
1429 // ntuple_evt->addItem("nHit2D_neg", m_nTrk2D_neg, m_nHit2D_neg);
1430 // ntuple_evt->addItem("chi2_2D_neg", m_nTrk2D_neg, m_chi2_2D_neg);
1431 // //global 3d fit
1432 // ntuple_evt->addItem("nTrk3D_neg", m_nTrk3D_neg, 0,100);
1433 // ntuple_evt->addItem("tanl_neg", m_nTrk3D_neg, m_tanl_neg);
1434 // ntuple_evt->addItem("tanl3D_neg", m_nTrk3D_neg, m_tanl3D_neg);
1435 // ntuple_evt->addItem("z0_neg", m_nTrk3D_neg, m_z0_neg);
1436 // ntuple_evt->addItem("z0_3D_neg", m_nTrk3D_neg, m_z0_3D_neg);
1437 // ntuple_evt->addItem("pro_neg", m_nTrk3D_neg, m_pro_neg);
1438 // ntuple_evt->addItem("pt3D_neg", m_nTrk3D_neg, m_pt3D_neg);
1439 // ntuple_evt->addItem("nHit3D_neg", m_nTrk3D_neg, m_nHit3D_neg);
1440 // ntuple_evt->addItem("chi2_3D_neg", m_nTrk3D_neg, m_chi2_3D_neg);
1441 //
1442 // //truth inf
1443 // // ntuple_evt->addItem("nTrkMC", m_nTrkMC,0,10);
1444 // ntuple_evt->addItem("pidTruth", m_pidTruth);
1445 // ntuple_evt->addItem("costaTruth", m_costaTruth);
1446 // ntuple_evt->addItem("ptTruth", m_ptTruth);
1447 // ntuple_evt->addItem("pzTruth", m_pzTruth);
1448 // ntuple_evt->addItem("pTruth", m_pTruth);
1449 // ntuple_evt->addItem("qTruth", m_qTruth);
1450 // ntuple_evt->addItem("drTruth", m_drTruth);
1451 // ntuple_evt->addItem("phiTruth", m_phi0Truth);
1452 // ntuple_evt->addItem("omegaTruth", m_omegaTruth);
1453 // ntuple_evt->addItem("dzTruth", m_dzTruth);
1454 // ntuple_evt->addItem("tanlTruth", m_tanl_mc);
1455 // ntuple_evt->addItem("rhoTruth", m_rho_mc);
1456 // ntuple_evt->addItem("thetaTruth", m_theta_mc);
1457
1458
1459 ntuple_evt->addItem("time", m_time);
1460
1461 } else { log << MSG::ERROR << "Cannot book tuple MdcHoughFinder/hough" <<endmsg;
1462 return StatusCode::FAILURE;
1463 }
1464 }
1465 /*
1466 NTuplePtr nt2(ntupleSvc(), "MdcHoughFinder/hit");
1467 if ( nt2 ){
1468 ntuplehit= nt2;
1469 } else {
1470 ntuplehit= ntupleSvc()->book("MdcHoughFinder/hit", CLID_ColumnWiseTuple, "hit");
1471 if(ntuplehit){
1472 ntuplehit->addItem("evt_stereo", m_evt_stereo);
1473 ntuplehit->addItem("run_stereo", m_run_stereo);
1474 ntuplehit->addItem("cos_stereo", m_cos_stereo);
1475 ntuplehit->addItem("tanlTruth_stereo", m_tanlTruth_stereo);
1476 ntuplehit->addItem("ncir_stereo", m_ncir_stereo);
1477 ntuplehit->addItem("npair_stereo", m_npair_stereo);
1478 ntuplehit->addItem("tanl_stereo", m_tanl_stereo);
1479 ntuplehit->addItem("tanl3D_stereo", m_tanl3D_stereo);
1480 ntuplehit->addItem("z0_stereo", m_z0_stereo);
1481 ntuplehit->addItem("z03D_stereo", m_z03D_stereo);
1482 ntuplehit->addItem("nHit_axial", m_nHit_axial, 0, 6796);
1483 ntuplehit->addItem("axial_layer", m_nHit_axial, m_axial_layer);
1484 ntuplehit->addItem("axial_wire", m_nHit_axial, m_axial_wire);
1485 ntuplehit->addItem("axial_deltaD", m_nHit_axial, m_axial_deltaD);
1486 ntuplehit->addItem("axial_flt", m_nHit_axial, m_axial_flt);
1487 ntuplehit->addItem("nHit_stereo", m_nHit_stereo, 0, 6796);
1488 ntuplehit->addItem("stereo_layer", m_nHit_stereo, m_stereo_layer);
1489 ntuplehit->addItem("stereo_wire", m_nHit_stereo, m_stereo_wire);
1490 ntuplehit->addItem("stereo_style", m_nHit_stereo, m_stereo_style);
1491 ntuplehit->addItem("stereo_z0", m_nHit_stereo, m_stereo_z0);
1492 ntuplehit->addItem("stereo_z1", m_nHit_stereo, m_stereo_z1);
1493 ntuplehit->addItem("stereo_s0", m_nHit_stereo, m_stereo_s0);
1494 ntuplehit->addItem("stereo_s1", m_nHit_stereo, m_stereo_s1);
1495 ntuplehit->addItem("stereo_z", m_nHit_stereo, m_stereo_z);
1496 ntuplehit->addItem("stereo_s", m_nHit_stereo, m_stereo_s);
1497 ntuplehit->addItem("stereo_zTruth", m_nHit_stereo, m_stereo_zTruth);
1498 ntuplehit->addItem("stereo_sTruth", m_nHit_stereo, m_stereo_sTruth);
1499 ntuplehit->addItem("stereo_deltaZ", m_nHit_stereo, m_stereo_deltaZ);
1500 ntuplehit->addItem("stereo_nsol", m_nHit_stereo, m_stereo_nsol);
1501 ntuplehit->addItem("stereo_ambig", m_nHit_stereo, m_stereo_ambig);
1502 ntuplehit->addItem("stereo_ambig_truth", m_nHit_stereo, m_stereo_ambig_truth);
1503 ntuplehit->addItem("stereo_disToCir", m_nHit_stereo, m_stereo_disToCir);
1504 ntuplehit->addItem("stereo_cirlist", m_nHit_stereo, m_stereo_cirlist);
1505 }
1506 }
1507 */
1508
1509 return StatusCode::SUCCESS;
1510}
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition: EvtStreamInputIterator.hh:100
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition: Gam4pikp.cxx:37
more_pt
bool more_pt(const HoughTrack *tracka, const HoughTrack *trackb)
Definition: Hough.cxx:960
sin
double sin(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:213
vectrk_for_clean
vector< TrkRecoTrk * > vectrk_for_clean
Definition: Hough3D.cxx:4
RecMdcHitCol
ObjectVector< RecMdcHit > RecMdcHitCol
Definition: InstallArea/include/MdcRecEvent/MdcRecEvent/RecMdcHit.h:99
RecMdcTrackCol
ObjectVector< RecMdcTrack > RecMdcTrackCol
Definition: InstallArea/include/MdcRecEvent/MdcRecEvent/RecMdcTrack.h:79
HitRefVec
SmartRefVector< RecMdcHit > HitRefVec
Definition: InstallArea/include/MdcRecEvent/MdcRecEvent/RecMdcTrack.h:22
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
more_pt
bool more_pt(const HoughTrack *tracka, const HoughTrack *trackb)
Definition: Hough.cxx:960
M_PI
#define M_PI
Definition: TConstant.h:4
BesTimer::start
void start(void)
Definition: BesTimer.cxx:27
BesTimer::stop
void stop(void)
Definition: BesTimer.cxx:39
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::radius
double radius(void) const
returns radious of helix.
Definition: InstallArea/include/TrackUtil/TrackUtil/Helix.h:215
Helix::dr
double dr(void) const
returns an element of parameters.
Definition: InstallArea/include/TrackUtil/TrackUtil/Helix.h:227
Helix::a
const HepVector & a(void) const
returns helix parameters.
Definition: InstallArea/include/TrackUtil/TrackUtil/Helix.h:263
Helix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition: InstallArea/include/TrackUtil/TrackUtil/Helix.h:209
Hough2D::setContext
static void setContext(TrkContextEv *context)
Definition: InstallArea/include/MdcHoughFinder/MdcHoughFinder/Hough2D.h:54
Hough2D::setCalib
static void setCalib(const MdcCalibFunSvc *mdcCalibFunSvc)
Definition: InstallArea/include/MdcHoughFinder/MdcHoughFinder/Hough2D.h:55
Hough3D::setContext
static void setContext(TrkContextEv *context)
Definition: InstallArea/include/MdcHoughFinder/MdcHoughFinder/Hough3D.h:55
Hough3D::setCalib
static void setCalib(const MdcCalibFunSvc *mdcCalibFunSvc)
Definition: InstallArea/include/MdcHoughFinder/MdcHoughFinder/Hough3D.h:56
HoughHitList::addTruthInfo
int addTruthInfo(std::map< int, const HepVector > mcTkPars)
Definition: HoughHitList.cxx:207
HoughHitList::getList
const std::vector< HoughHit > & getList() const
Definition: HoughHitList.cxx:344
HoughHit::setMdcGeomSvc
static void setMdcGeomSvc(MdcGeomSvc *geomSvc)
Definition: InstallArea/include/MdcHoughFinder/MdcHoughFinder/HoughHit.h:36
HoughHit::printTruth
void printTruth() const
Definition: HoughHit.cxx:194
HoughHit::setMdcCalibFunSvc
static void setMdcCalibFunSvc(const MdcCalibFunSvc *calibSvc)
Definition: InstallArea/include/MdcHoughFinder/MdcHoughFinder/HoughHit.h:35
HoughTrack::trackCharge2D
int trackCharge2D()
Definition: HoughTrack.cxx:820
HoughTrack::find_stereo_hit
int find_stereo_hit()
Definition: HoughTrack.cxx:697
HoughTrack::fit2D
int fit2D(double bunchtime)
Definition: HoughTrack.cxx:190
HoughTrack::trackCharge3D
int trackCharge3D()
Definition: HoughTrack.cxx:838
HoughTrack::setMdcHit
void setMdcHit(const vector< MdcHit * > *mdchit)
Definition: InstallArea/include/MdcHoughFinder/MdcHoughFinder/HoughTrack.h:91
HoughTrack::fit3D_inner
int fit3D_inner()
Definition: HoughTrack.cxx:312
HoughTrack::find_pair_hit
int find_pair_hit()
Definition: HoughTrack.cxx:946
IBesTimerSvc::addItem
virtual BesTimer * addItem(const std::string &name)=0
MdcDetector::instance
static MdcDetector * instance()
Definition: MdcDetector.cxx:21
MdcHoughFinder::execute
StatusCode execute()
Definition: Hough.cxx:248
MdcHoughFinder::initialize
StatusCode initialize()
Definition: Hough.cxx:140
MdcHoughFinder::beginRun
StatusCode beginRun()
Definition: Hough.cxx:131
MdcHoughFinder::MdcHoughFinder
MdcHoughFinder(const std::string &name, ISvcLocator *pSvcLocator)
Definition: Hough.cxx:76
MdcHoughFinder::finalize
StatusCode finalize()
Definition: Hough.cxx:621
MdcHoughFinder::bookTuple
StatusCode bookTuple()
Definition: Hough.cxx:1333
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
MdcPrintSvc::printRecMdcTrackCol
void printRecMdcTrackCol() const
Definition: MdcPrintSvc.cxx:109
MdcTrackListBase::remove
void remove(MdcTrack *atrack)
Definition: MdcTrackListBase.cxx:608
Pdt::readMCppTable
static void readMCppTable(std::string filenm)
RawData::getTrackIndex
int getTrackIndex() const
Definition: RawData.cxx:50
TrkFit::helix
virtual TrkExchangePar helix(double fltL) const =0
TrkHotList::nHit
virtual int nHit(TrkEnums::TrkViewInfo view=TrkEnums::bothView) const =0
TrkRecoTrk::fitResult
const TrkFit * fitResult() const
Definition: TrkRecoTrk.cxx:387