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