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