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