CGEM BOSS 6.6.5.h
BESIII Offline Software System
Loading...
Searching...
No Matches
Hough.cxx
Go to the documentation of this file.
1#include "GaudiKernel/SmartDataPtr.h"
2#include "GaudiKernel/Bootstrap.h"
3#include "GaudiKernel/IDataManagerSvc.h"
4#include "GaudiKernel/IDataProviderSvc.h"
5#include "GaudiKernel/IPartPropSvc.h"
6
7#include "HoughTransAlg/Hough.h"
8#include "EventModel/EventHeader.h"
9#include "EvTimeEvent/RecEsTime.h"
10#include "McTruth/McParticle.h"
11#include "HepPDT/ParticleDataTable.hh"
12#include "MdcTrkRecon/MdcTrackParams.h"
13#include "MdcRecoUtil/Pdt.h"
14#include "Identifier/Identifier.h"
15#include "Identifier/CgemID.h"
16#include "Identifier/MdcID.h"
17#include "MdcGeom/MdcDetector.h"
18#include "MdcData/MdcRecoHitOnTrack.h"
19#include "CgemData/CgemHitOnTrack.h"
20#include "TrkBase/TrkRep.h"
21#include "TrkBase/TrkDifTraj.h"
22#include "TrkBase/TrkExchangePar.h"
23
24#include <math.h>
25#include <iostream>
26#include <fstream>
27
28class TrkRep;
29using namespace std;
30
31HoughFinder::HoughFinder(const std::string& name, ISvcLocator* pSvcLocator):Algorithm(name, pSvcLocator)
32{
33 declareProperty("debug", m_debug = 0);
34 declareProperty("cgem", m_cgem = 1);
35 declareProperty("mcTruth", m_mcTruth = 0);
36 declareProperty("fillNTuple", m_fillNTuple= 0);
37 declareProperty("trackCharge", m_trackCharge = 0);
38 declareProperty("findPeakMethod",m_findPeakMethod=0);
39
40 //declareProperty("nStep_xy", m_nStep_xy = 1);
41 //declareProperty("xBin_xy", m_xBin_xy);
42 //declareProperty("yBin_xy", m_yBin_xy);
43 //declareProperty("xExtend_xy", m_xExtend_xy);
44 //declareProperty("yExtend_xy", m_yExtend_xy);
45 //declareProperty("nVote_xy", m_nVote_xy);
46 //declareProperty("nRMS_xy", m_nRMS_xy);
47 //declareProperty("xInclude_xy", m_xInclude_xy = 2);
48 //declareProperty("yInclude_xy", m_yInclude_xy = 2);
49
50 //declareProperty("nStep_sz", m_nStep_sz = 1);
51 //declareProperty("xBin_sz", m_xBin_sz);
52 //declareProperty("yBin_sz", m_yBin_sz);
53 //declareProperty("xExtend_sz", m_xExtend_sz);
54 //declareProperty("yExtend_sz", m_yExtend_sz);
55 //declareProperty("nVote_sz", m_nVote_sz);
56 //declareProperty("nRMS_sz", m_nRMS_sz);
57 //declareProperty("xInclude_sz", m_xInclude_sz = 2);
58 //declareProperty("yInclude_sz", m_yInclude_sz = 2);
59
60 //read raw data setup
61 declareProperty("keepBadTdc", m_keepBadTdc = 0);
62 declareProperty("dropHot", m_dropHot = 0);
63 declareProperty("keepUnmatch", m_keepUnmatch = 0);
64 declareProperty("driftTimeUpLimit", m_driftTimeUpLimit = 1500);
65 //declareProperty("",m_ = 0);
66 declareProperty("pdtFile", m_pdtFile = "pdt.table");
67
68 // tests by llwang
69 declareProperty("nBinTheta", m_nBinTheta=100 );
70 declareProperty("nBinRho", m_nBinRho=50 );
71 declareProperty("rhoRange", m_rhoRange = 0.1);
72 declareProperty("ExtPeak_theta",m_ExtPeak_theta= 4);
73 declareProperty("ExtPeak_rho", m_ExtPeak_rho= 4);
74 declareProperty("shareHitRate", m_shareHitRate = 0.7);
75
76 declareProperty("nBinTanl", m_nBinTanl = 100);
77 declareProperty("nBinDz", m_nBinDz = 100);
78 declareProperty("ExtPeak_tanl", m_ExtPeak_tanl = 1);
79 declareProperty("ExtPeak_dz", m_ExtPeak_dz = 1);
80 declareProperty("filter", m_filter= 0);
81 declareProperty("eventFile", m_evtFile= "EventList.txt");
82 declareProperty("fitFlag", m_fitFlag= 3);
83 declareProperty("storeFlag", storeFlag=1);
84 declareProperty("checkHits", m_checkHits=1);
85 declareProperty("Layer", Layer=20);
86 declareProperty("clearTrack", m_clearTrack=1);
87 declareProperty("useCgemInGlobalFit",m_useCgemInGlobalFit=3);
88 declareProperty("printFlag", printFlag=0);
89 declareProperty("removeNOuterHits", m_removeNOuterHits = 0);
90
91 // --- cuts to combine tracks
92 declareProperty("drCut_combine", m_drCut=2.0);
93 declareProperty("phi0Cut_combine", m_phi0Cut=0.2);
94 declareProperty("kappaCut_combine",m_kappaCut=1.5);
95 declareProperty("dzCut_combine", m_dzCut=5.0);
96 declareProperty("tanlCut_combine", m_tanlCut=0.5);
97
98 // --- cuts for hits saved in RecMdcTrack
99 declareProperty("chi2CutHits", m_chi2CutHits=10);
100
101 m_maxFireLayer = -1;
102 m_totEvtProcessed=0;
103}
104
105StatusCode HoughFinder::initialize()
106{
107 MsgStream log(msgSvc(), name());
108 log << MSG::INFO << "HoughFinder::initialize()" << endreq;
109 StatusCode sc;
110
111 // --- RawData
112 IRawDataProviderSvc* irawDataProviderSvc;
113 sc = service ("RawDataProviderSvc", irawDataProviderSvc);
114 m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*> (irawDataProviderSvc);
115 if ( sc.isFailure() ){
116 log << MSG::ERROR << name()<<" Could not load RawDataProviderSvc!" << endreq;
117 return StatusCode::FAILURE;
118 }
119
120
121 // --- MDC Geometry
122 IMdcGeomSvc* imdcGeomSvc;
123 sc = service ("MdcGeomSvc", imdcGeomSvc);
124 m_mdcGeomSvc = dynamic_cast<MdcGeomSvc*> (imdcGeomSvc);
125 if ( sc.isFailure() ){
126 log << MSG::ERROR << "Could not load MdcGeoSvc!" << endreq;
127 return StatusCode::FAILURE;
128 }
129
130 // --- MdcCalibFunSvc
131 IMdcCalibFunSvc* imdcCalibSvc;
132 sc = service ("MdcCalibFunSvc", imdcCalibSvc);
133 m_mdcCalibFunSvc = dynamic_cast<MdcCalibFunSvc*>(imdcCalibSvc);
134 if ( sc.isFailure() ){
135 log << MSG::ERROR << "Could not load MdcCalibFunSvc!" << endreq;
136 return StatusCode::FAILURE;
137 }
138
139 if(m_cgem){
140 // --- CGEM geometry
141 ICgemGeomSvc* icgemGeomSvc;
142 sc = service ("CgemGeomSvc", icgemGeomSvc);
143 m_cgemGeomSvc = dynamic_cast<CgemGeomSvc*> (icgemGeomSvc);
144 if ( sc.isFailure() ){
145 log << MSG::ERROR << "Could not load CgemGeomSvc!" << endreq;
146 return StatusCode::FAILURE;
147 }
148
149 // --- CgemCalibFunSvc
150 ICgemCalibFunSvc* icgemCalibSvc;
151 sc = service ("CgemCalibFunSvc", icgemCalibSvc);
152 m_cgemCalibFunSvc = dynamic_cast<CgemCalibFunSvc*>(icgemCalibSvc);
153 if ( sc.isFailure() ){
154 log << MSG::ERROR << "Could not load CgemCalibFunSvc!" << endreq;
155 return StatusCode::FAILURE;
156 }
157
158 }
159
160 // --- MagneticFieldSvc
161 sc = service ("MagneticFieldSvc",m_pIMF);
162 if( sc.isFailure() ) {
163 log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
164 return StatusCode::FAILURE;
165 }
166 m_bfield = new BField(m_pIMF);
167 //log << MSG::INFO << "field z = "<<m_bfield->bFieldNominal()<< endreq;
168 m_trkContextEv = new TrkContextEv(m_bfield);
169
170 m_mdcDetector = new MdcDetector();
171 //read pdt
172 Pdt::readMCppTable(m_pdtFile);
173
174 // --- Mdc Print Service
175 //IMdcPrintSvc* imdcPrintSvc;
176 //sc = service ("MdcPrintSvc", imdcPrintSvc);
177 //m_mdcPrintSvc = dynamic_cast<MdcPrintSvc*> (imdcPrintSvc);
178 //if ( sc.isFailure() ){
179 //log << MSG::FATAL << "Could not load MdcPrintSvc!" << endreq;
180 //return StatusCode::FAILURE;
181 //}
182
183 // --- Bes Timer Service
184 //sc = service( "BesTimerSvc", m_timersvc);
185 //if( sc.isFailure() ) {
186 //log << MSG::WARNING << " Unable to locate BesTimerSvc" << endreq;
187 //return StatusCode::FAILURE;
188 //}
189 //m_timer_all = m_timersvc->addItem("Execution");
190 //m_timer_all->propName("nExecution");
191
192 HoughHit::setMdcGeomSvc(m_mdcGeomSvc);
193 HoughHit::setMdcCalibFunSvc(m_mdcCalibFunSvc);
194 HoughHit::setCgemGeomSvc(m_cgemGeomSvc);
195 HoughHit::setCgemCalibFunSvc(m_cgemCalibFunSvc);
196 HoughHit::setMdcDetector(m_mdcDetector);
197 HoughTrack::m_clearTrack = m_clearTrack;
198 HoughTrack::m_useCgemInGlobalFit = m_useCgemInGlobalFit;
199
200 m_roughRhoThetaMap=TH2D("roughRhoThetaMap","roughRhoThetaMap",m_nBinTheta,0.,M_PI, m_nBinRho, -1.*m_rhoRange, m_rhoRange);
201 //m_tanlRange = 0.93/sqrt(1-0.93*0.93);//FIXME
202 m_tanlRange = 3.0;//FIXME
203 m_dzRange = 50;
204 m_roughTanlDzMap = TH2D("roughTanlDzMap","roughTanlDzMap",m_nBinTanl,-1.*m_tanlRange,m_tanlRange,m_nBinDz,-1.*m_dzRange,m_dzRange);
205
206 const int nPt=12;
207 double rho[nPt] = {0.07, 0.056,0.045,0.038,0.0335,0.03,0.0198,0.0148,0.0074,0.00376,0.00303,0.00157};
208 double cut1_Cgem_99[12]={-2.14,-1.36, -0.6 , -0.46, -0.35, -0.59, -0.32, -0.26, -0.22, -0.21, -0.21, -0.211};
209 double cut2_Cgem_99[12]={ 0.5 , 1.77, 1.99, 1.94, 1.99, 1.9 , 0.31, 0.27, 0.24, 0.23, 0.23, 0.222};
210 double cut1_ODC1_99[12]={-1.28,-0.71, -0.58, -0.62, -0.64, -0.68, -0.18, -0.14, -0.11, -0.1 , -0.1 , -0.11 };
211 double cut2_ODC1_99[12]={ 0.9 , 0.74, 0.42, 0.37, 0.32, 0.37, 0.28, 0.25, 0.24, 0.24, 0.24, 0.23};
212 double cut1_ODC2_99[12]={ -2.14, -1.25, -1.28, -0.86, -0.47, -0.78, -0.36, -0.44, -0.61, -0.67, -0.64, -0.82 };
213 double cut2_ODC2_99[12]={ -1.35, 0.55 , 0.53 , 0.83 , 0.85 , 0.83 , 0.38 , 0.55 , 0.49 , 0.46 , 0.49 , 0.4};
214 m_cut1_cgem = new TGraph(nPt, rho, cut1_Cgem_99);
215 m_cut2_cgem = new TGraph(nPt, rho, cut2_Cgem_99);
216 m_cut1_ODC1 = new TGraph(nPt, rho, cut1_ODC1_99);
217 m_cut2_ODC1 = new TGraph(nPt, rho, cut2_ODC1_99);
218 m_cut1_ODC2 = new TGraph(nPt, rho, cut1_ODC2_99);
219 m_cut2_ODC2 = new TGraph(nPt, rho, cut2_ODC2_99);
220
221
222 int s(0),l(0),npt(0);
223 double momenta[100] = {0};
224 double cuts[2][43][100] = {0};
225 string cutFilePath = getenv("HOUGHTRANSALGROOT");
226 cutFilePath += "/share/cut.txt";
227 ifstream fin(cutFilePath.c_str(), ios::in);
228 //cout<<"cutfiledir:"<<cutFilePath.c_str()<<endl;
229 //ifstream fin("/workfs/bes/huangzhen/workarea/CgemBoss/6.6.5.c/Reconstruction/HoughFinder/HoughFinder-00-00-13/share/cut.txt",ios::in);
230 if (!fin.good())
231 cout << "Error in " << __FILE__<<" when opening "<<cutFilePath.c_str()<<endl;
232 bool firstline(true);
233 string strcom;
234 while(std::getline(fin, strcom)){
235 if(firstline){
236 fin >> npt
237 >> momenta[0] >> momenta[1] >> momenta[2] >> momenta[3] >> momenta[4] >> momenta[5] >> momenta[6] >> momenta[7] >> momenta[8] >> momenta[9]
238 >> momenta[10] >> momenta[11] >> momenta[12] >> momenta[13] >> momenta[14] >> momenta[15] >> momenta[16] >> momenta[17] >> momenta[18] >> momenta[19]
239 >> momenta[20] >> momenta[21] >> momenta[22] >> momenta[23] >> momenta[24] >> momenta[25] >> momenta[26] >> momenta[27] >> momenta[28] >> momenta[29]
240 >> momenta[30] >> momenta[31] >> momenta[32] >> momenta[33] >> momenta[34] >> momenta[35] >> momenta[36] >> momenta[37] >> momenta[38] >> momenta[39]
241 >> momenta[40] >> momenta[41] >> momenta[42] >> momenta[43] >> momenta[44] >> momenta[45] >> momenta[46] >> momenta[47] >> momenta[48] >> momenta[49]
242 >> momenta[50] >> momenta[51] >> momenta[52] >> momenta[53] >> momenta[54] >> momenta[55] >> momenta[56] >> momenta[57] >> momenta[58] >> momenta[59]
243 >> momenta[60] >> momenta[61] >> momenta[62] >> momenta[63] >> momenta[64] >> momenta[65] >> momenta[66] >> momenta[67] >> momenta[68] >> momenta[69]
244 >> momenta[70] >> momenta[71] >> momenta[72] >> momenta[73] >> momenta[74] >> momenta[75] >> momenta[76] >> momenta[77] >> momenta[78] >> momenta[79]
245 >> momenta[80] >> momenta[81] >> momenta[82] >> momenta[83] >> momenta[84] >> momenta[85] >> momenta[86] >> momenta[87] >> momenta[88] >> momenta[89]
246 >> momenta[90] >> momenta[91] >> momenta[92] >> momenta[93] >> momenta[94] >> momenta[95] >> momenta[96] >> momenta[97] >> momenta[98] >> momenta[99]
247 ;
248 firstline = false;
249 continue;
250 }
251 fin >> s >> l;
252 fin >> cuts[s][l][0] >> cuts[s][l][1] >> cuts[s][l][2] >> cuts[s][l][3] >> cuts[s][l][4] >> cuts[s][l][5] >> cuts[s][l][6] >> cuts[s][l][7] >> cuts[s][l][8] >> cuts[s][l][9]
253 >> cuts[s][l][10] >> cuts[s][l][11] >> cuts[s][l][12] >> cuts[s][l][13] >> cuts[s][l][14] >> cuts[s][l][15] >> cuts[s][l][16] >> cuts[s][l][17] >> cuts[s][l][18] >> cuts[s][l][19]
254 >> cuts[s][l][20] >> cuts[s][l][21] >> cuts[s][l][22] >> cuts[s][l][23] >> cuts[s][l][24] >> cuts[s][l][25] >> cuts[s][l][26] >> cuts[s][l][27] >> cuts[s][l][28] >> cuts[s][l][29]
255 >> cuts[s][l][30] >> cuts[s][l][31] >> cuts[s][l][32] >> cuts[s][l][33] >> cuts[s][l][34] >> cuts[s][l][35] >> cuts[s][l][36] >> cuts[s][l][37] >> cuts[s][l][38] >> cuts[s][l][39]
256 >> cuts[s][l][40] >> cuts[s][l][41] >> cuts[s][l][42] >> cuts[s][l][43] >> cuts[s][l][44] >> cuts[s][l][45] >> cuts[s][l][46] >> cuts[s][l][47] >> cuts[s][l][48] >> cuts[s][l][49]
257 >> cuts[s][l][50] >> cuts[s][l][51] >> cuts[s][l][52] >> cuts[s][l][53] >> cuts[s][l][54] >> cuts[s][l][55] >> cuts[s][l][56] >> cuts[s][l][57] >> cuts[s][l][58] >> cuts[s][l][59]
258 >> cuts[s][l][60] >> cuts[s][l][61] >> cuts[s][l][62] >> cuts[s][l][63] >> cuts[s][l][64] >> cuts[s][l][65] >> cuts[s][l][66] >> cuts[s][l][67] >> cuts[s][l][68] >> cuts[s][l][69]
259 >> cuts[s][l][70] >> cuts[s][l][71] >> cuts[s][l][72] >> cuts[s][l][73] >> cuts[s][l][74] >> cuts[s][l][75] >> cuts[s][l][76] >> cuts[s][l][77] >> cuts[s][l][78] >> cuts[s][l][79]
260 >> cuts[s][l][80] >> cuts[s][l][81] >> cuts[s][l][82] >> cuts[s][l][83] >> cuts[s][l][84] >> cuts[s][l][85] >> cuts[s][l][86] >> cuts[s][l][87] >> cuts[s][l][88] >> cuts[s][l][89]
261 >> cuts[s][l][90] >> cuts[s][l][91] >> cuts[s][l][92] >> cuts[s][l][93] >> cuts[s][l][94] >> cuts[s][l][95] >> cuts[s][l][96] >> cuts[s][l][97] >> cuts[s][l][98] >> cuts[s][l][99]
262 ;
263 if(fin.peek()<0)break;
264 }
265 for(int is=0;is<=1;is++){
266 for(int il=0;il<43;il++){
267 TGraph* gr = new TGraph();
268 int point(0);
269 //cout<<il<<" ";
270 for(int ipt=0;ipt<npt;ipt++){
271 //cout<<cuts[is][il][ipt]<<" ";
272 if(fabs(cuts[is][il][ipt])<1e-6)continue;
273 //cut[is][il]->SetPoint(point,momenta[ipt]/1000.,cuts[is][il][ipt]);
274 gr->SetPoint(point,momenta[ipt]/1000.,cuts[is][il][ipt]);
275 point++;
276 }
277 if(point>0)HoughTrack::m_cut[is][il] = gr;
278 //cout<<endl;
279 }
280 }
281
282
283 if(m_fillNTuple) sc = bookTuple();
284 if ( sc.isFailure() ){
285 log << MSG::FATAL << "Could not book Tuple !" << endreq;
286 return StatusCode::FAILURE;
287 }
288
289
290 myDotsHelixFitter.initialize();
291
292 return StatusCode::SUCCESS;
293}
294
295StatusCode HoughFinder::execute()
296{
297 MsgStream log(msgSvc(), name());
298 log << MSG::INFO << "HoughFinder::execute()" << endreq;
299 cout.precision(6);
300
301
302 //cout<<"line "<<__LINE__<<endl;
303 // --- Get event header
304 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
305 if(!eventHeader){
306 log << MSG::ERROR << "Could not find Event Header" << endreq;
307 return StatusCode::FAILURE;
308 }
309 m_run = eventHeader->runNumber();
310 m_event = eventHeader->eventNumber();
311 m_totEvtProcessed++; //if(m_totEvtProcessed%100==0) cout<<"HoughFinder: "<<m_totEvtProcessed<<" events processed."<<endl;
312
313 //cout<<"line "<<__LINE__<<endl;
314 // --- event start time
315 SmartDataPtr<RecEsTimeCol> evTimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
316 if(!evTimeCol){
317 log << MSG::ERROR << "Could not retrieve RecEsTimeCol" << endreq;
318 return StatusCode::FAILURE;
319 }
320 m_bunchT0 = 0;
321 RecEsTimeCol::iterator iter_evt = evTimeCol->begin();
322 if (iter_evt != evTimeCol->end()){
323 m_bunchT0 = (*iter_evt)->getTest()*1.e-9;//m_bunchT0-s, getTest-ns
324 //cout<<"this event T0 = "<<m_bunchT0<<endl;
325 }
326
327 // --- clear, register and store reconstructed track
328 registerTrack(m_trackList_tds,m_hitList_tds);
329
330
331 //cout<<"line "<<__LINE__<<endl;
332 if(m_filter){
333 ifstream lowPt_Evt;
334 lowPt_Evt.open(m_evtFile.c_str());
335 vector<int> evtlist;
336 int evtNum;
337 while( lowPt_Evt >> evtNum) {
338 evtlist.push_back(evtNum);
339 }
340 vector<int>::iterator iter_evt = find(evtlist.begin(),evtlist.end(),m_event);
341 if(iter_evt == evtlist.end()){
342 setFilterPassed(false);
343 return StatusCode::SUCCESS;
344 }
345 else{
346 setFilterPassed(true);
347 }
348 }
349 //if(m_event<2000||m_event>3000){
350 //setFilterPassed(false);
351 //return StatusCode::SUCCESS;
352 //}
353
354 if(m_debug)
355 {
356 //cout<<endl;
357 cout<<"===================================================================================================="<<endl;
358 cout<<"run:"<<m_run<<" , event: "<<m_event<<endl;
359 //cout<<m_event<<" , ";
360 //cout<<endl;
361 }
362
363 //cout<<"line "<<__LINE__<<endl;
364 // --- prepare hits for reconstruction
365 int nHit = makeHoughHitList();
366
367 //cout<<"line "<<__LINE__<<endl;
368 int nMcHit;
369 if(m_run<0&&m_mcTruth){
370 int nMcHit = getMcHitCol();
371 //cout<<"nMcHit "<<nMcHit<<endl;
372 if(printFlag)cout<<"@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ McParticle @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"<<endl;
373 getMcParticleCol();
374 if(printFlag)cout<<"@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ McParticle @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"<<endl;
375 if(printFlag)cout<<endl;
376 if(printFlag)cout<<endl;
377 if(printFlag)cout<<endl;
378 //halfCircle(m_mcTrackCol,m_mcHitCol);
379 }
380
381 //cout<<"line "<<__LINE__<<endl;
382 searchCircle();
383 //cout<<"ntracks:"<<m_houghTrackList.size()<<endl;
384 //cout<<"line "<<__LINE__<<endl;
385 checkHot(m_houghTrackList); // remove tracks with too many shared hits (if >70%)
386 //cout<<"line "<<__LINE__<<endl;
387 solveSharedHits(m_XHoughHitList); // solve hits sharing by smaller distance
388 //checkTrack(m_houghTrackList);
389 //cout<<"line "<<__LINE__<<endl;
390 associateVHits();
391 //cout<<"line "<<__LINE__<<endl;
392 if(printFlag)cout<<"################################################## HoughTrack ##################################################"<<endl;
393 if(storeFlag==1)storeTracks(m_trackList_tds, m_hitList_tds, m_houghTrackList);
394 else
395 if(storeFlag==2)storeTrack(m_houghTrackList,m_trackList_tds,m_hitList_tds);
396 else
397 if(storeFlag==3)storeRecTracks(m_trackList_tds, m_hitList_tds, m_houghTrackList);
398 if(printFlag)cout<<"################################################## HoughTrack ##################################################"<<endl;
399 if(printFlag)cout<<endl;
400 if(printFlag)cout<<endl;
401 if(printFlag)cout<<endl;
402 //cout<<"line "<<__LINE__<<endl;
403
404 if(m_fillNTuple)dumpHit();
405 if(m_fillNTuple==1)dumpHoughTrack();
406 if(m_fillNTuple==1)dumpHoughEvent();
407 if(m_fillNTuple==2)dumpTdsTrack();
408 if(m_fillNTuple==2)dumpTdsEvent();
409 //cout<<"line "<<__LINE__<<endl;
410
411 if(printFlag)cout<<"$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ TdsTrack $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$"<<endl;
412 if(printFlag)printTdsTrack();
413 if(printFlag)cout<<"$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ TdsTrack $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$"<<endl;
414 if(printFlag)cout<<endl;
415 if(printFlag)cout<<endl;
416 if(printFlag)cout<<endl;
417
418 if(printFlag)cout<<"&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& HoughHit &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&"<<endl;
419 if(printFlag)printHitList();
420 if(printFlag)cout<<"&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& HoughHit &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&"<<endl;
421 if(printFlag)cout<<endl;
422 if(printFlag)cout<<endl;
423 if(printFlag)cout<<endl;
424
425 clearTrack(m_houghTrackList);
426 //cout<<"line "<<__LINE__<<endl;
427
428 clearMemory();
429 //cout<<"line "<<__LINE__<<endl;
430 if(m_debug)cout<<endl;
431
432 return StatusCode::SUCCESS;
433}// <-- execute()
434
435StatusCode HoughFinder::finalize()
436{
437 MsgStream log(msgSvc(), name());
438 log << MSG::INFO << "HoughFinder::finalize()" << endreq;
439 return StatusCode::SUCCESS;
440}
441
442bool largerPt(HoughTrack trk1, HoughTrack trk2)
443{
444 return fabs(trk1.pt())>fabs(trk2.pt());
445}
446
447int HoughFinder::makeHoughHitList()
448{
449 MsgStream log(msgSvc(), name());
450 int nHit(0);
451 if(!(m_mdcHitCol.empty())) m_mdcHitCol.clear();
452 if(!(m_houghHitList.empty())) m_houghHitList.clear();
453 if(!(m_XHoughHitList.empty())) m_XHoughHitList.clear();
454 if(!(m_VHoughHitList.empty())) m_VHoughHitList.clear();
455
456 if(m_cgem){
457 // --- RecCgemCluster
458 SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
459 if(!recCgemClusterCol){
460 log << MSG::ERROR << "Could not retrieve RecCgemClusterCol" << endreq;
461 }else{
462 RecCgemClusterCol::iterator cgemClusterIter = recCgemClusterCol->begin();
463 m_recCgemClusterColBegin=cgemClusterIter;
464 for (;cgemClusterIter != recCgemClusterCol->end(); cgemClusterIter++,nHit++){
465 const RecCgemCluster* recCgemCluster = (*cgemClusterIter);
466 //m_recCgemCluster = (*cgemClusterIter);
467 HoughHit hit(recCgemCluster, m_bunchT0, nHit);
468 m_houghHitList.push_back(hit);
469 if(hit.getLayer()>m_maxFireLayer)m_maxFireLayer = hit.getLayer();
470 }
471 }
472 }
473
474 uint32_t getDigiFlag = 0;
475 if(m_dropHot)getDigiFlag |= MdcRawDataProvider::b_dropHot;
476 if(m_keepBadTdc) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
477 if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
478 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
479
480 MdcDigiVec::iterator mdcDigiIter = mdcDigiVec.begin();
481 for (;mdcDigiIter != mdcDigiVec.end(); mdcDigiIter++,nHit++){
482 const MdcDigi* mdcDigi = (*mdcDigiIter);
483 HoughHit hit(mdcDigi, m_bunchT0, nHit);
484 if(fabs(hit.driftTime())>m_driftTimeUpLimit) continue;
485 MdcHit *mdcHit= new MdcHit(mdcDigi,m_mdcDetector);
486 hit.setMdcHit(mdcHit);
487
488 m_mdcHitCol.push_back(mdcHit);
489 m_houghHitList.push_back(hit);
490 if(hit.getLayer()>m_maxFireLayer)m_maxFireLayer = hit.getLayer();
491 }
492
493 for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end();hitIt++){
494 int flag=hitIt->getFlag();
495 if(hitIt->getHitType()==HoughHit::CGEMHIT){
496 if(flag==0)m_XHoughHitList.push_back(&(*hitIt));
497 if(flag==2)m_VHoughHitList.push_back(&(*hitIt));
498 }
499 if(hitIt->getHitType()==HoughHit::MDCHIT){
500 if(flag==0)m_XHoughHitList.push_back(&(*hitIt));
501 else m_VHoughHitList.push_back(&(*hitIt));
502 }
503 }
504 //for(vector<HoughHit*>::iterator hitIt = m_XHoughHitList.begin(); hitIt != m_XHoughHitList.end(); hitIt++){
505 //(*hitIt)->print();
506 //}
507 //for(vector<HoughHit*>::iterator hitIt = m_VHoughHitList.begin(); hitIt != m_VHoughHitList.end(); hitIt++){
508 //(*hitIt)->print();
509 //}
510
511 return nHit;
512}
513
514int HoughFinder::getMcHitCol()
515{
516 MsgStream log(msgSvc(), name());
517 int nMcHit(0);
518 if(!(m_mcHitCol.empty()))m_mcHitCol.clear();
519
520 // --- get cgemMcHit
521 if(m_cgem){
522 SmartDataPtr<CgemMcHitCol> cgemMcHitCol(eventSvc(), "/Event/MC/CgemMcHitCol");
523 if (!cgemMcHitCol){
524 log << MSG::ERROR << "Could not retrieve CgemMcHitCol" << endreq;
525 }else{
526 CgemMcHitCol::iterator cgemMcHitIt = cgemMcHitCol->begin();
527 for(;cgemMcHitIt != cgemMcHitCol->end();cgemMcHitIt++,nMcHit++){
528 const CgemMcHit* cgemMcHit = (*cgemMcHitIt);
529 HoughHit hit(cgemMcHit,m_bunchT0,nMcHit);
530 m_mcHitCol.push_back(hit);
531 HoughHit* phit = &(*m_mcHitCol.rbegin());
532
533 // --- match fired strip ID between CemgMcHit and RecCgemCluster
534 vector<int> xFireStripID;
535 vector<int> vFireStripID;
536 int sheet;
537 TArrayI identifier = phit->getCgemMcHit()->GetIdentifier();
538 for(int ii = 0; ii < identifier.GetSize(); ii++){
539 //const Identifier ident = Identifier(identifier.GetAt(ii));
540 const Identifier ident = Identifier(identifier.At(ii));
541 if(CgemID::is_xstrip(ident))xFireStripID.push_back(CgemID::strip(ident));
542 else vFireStripID.push_back(CgemID::strip(ident));
543 sheet = CgemID::sheet(ident);
544 //cout
545 //<< " hitID:" << phit->getHitID()
546 //<< " layerID=" << CgemID::layer(ident)
547 //<< " sheetID=" << CgemID::sheet(ident)
548 //<< " isXstrip=" << CgemID::is_xstrip(ident)
549 //<< " stripID=" << CgemID::strip(ident)
550 //<<endl;
551 }
552 sort(xFireStripID.begin(),xFireStripID.end());
553 sort(vFireStripID.begin(),vFireStripID.end());
554 //vector<int>::iterator uniqueX = unique(xFireStripID.begin(),xFireStrip ID.end());
555 //vector<int>::iterator uniqueV = unique(vFireStripID.begin(),vFireStrip ID.end());
556 //xFireStripID.erase(uniqueX,xFireStripID.end());
557 //vFireStripID.erase(uniqueV,vFireStripID.end());
558 //for(vector<int>::iterator i = vFireStripID.begin();i!=vFireStripID.end();i++)cout<<*i<<endl;
559 //cout<<endl;
560
561 bool findX(false), findV(false);
562 const RecCgemCluster* xCluster = NULL;
563 const RecCgemCluster* vCluster = NULL;
564 for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end();hitIt++){
565 if(hitIt->getHitType()!=HoughHit::CGEMHIT)continue;
566 const RecCgemCluster* recCgemCluster = hitIt->getCgemCluster();
567 if(recCgemCluster->getlayerid() != phit->getLayer())continue;
568 if(!(recCgemCluster->getflag()==0||recCgemCluster->getflag()==1))continue;
569 //if(recCgemCluster->getsheetid() != sheet)continue;
570 //if(recCgemCluster->getTrkId() != *(phit->getTrkID().begin())continue;
571 if(!findX){
572 findX = recCgemCluster->getflag() == 0
573 && xFireStripID.size() > 0
574 && recCgemCluster->getclusterflagb() == *(xFireStripID.begin())
575 && recCgemCluster->getclusterflage() == *(xFireStripID.rbegin());
576 if(findX){
577 xCluster = recCgemCluster;
578 hitIt->setPairHit(phit);
579 }
580 }
581 if(!findV){
582 findV = recCgemCluster->getflag() == 1
583 && vFireStripID.size() > 0
584 && recCgemCluster->getclusterflagb() == *(vFireStripID.begin())
585 && recCgemCluster->getclusterflage() == *(vFireStripID.rbegin());
586 if(findV){
587 vCluster = recCgemCluster;
588 hitIt->setPairHit(phit);
589 }
590 }
591
592 }
593
594 /// --- put match RecCgemCluster into CemgMcHit
595 if(findX&&findV){
596 for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end();hitIt++){
597 if(hitIt->getHitType()!=HoughHit::CGEMHIT)continue;
598 const RecCgemCluster* recCgemCluster = hitIt->getCgemCluster();
599 if(recCgemCluster->getlayerid() != phit->getLayer())continue;
600 if(recCgemCluster->getflag()==0||recCgemCluster->getflag()==1)continue;
601 //if(recCgemCluster->getsheetid() != sheet)continue;
602 //if(recCgemCluster->getTrkId() != *(phit->getTrkID().begin())continue;
603 if(recCgemCluster->getclusterflagb() == xCluster->getclusterid()
604 && recCgemCluster->getclusterflage() == vCluster->getclusterid()){
605 phit->setCgemCluster(recCgemCluster);
606 phit->setWire(recCgemCluster->getclusterid());
607 phit->setFlag(recCgemCluster->getflag());
608 phit->setPairHit(&(*hitIt));
609 hitIt->setPairHit(phit);
610 }
611 }
612 }
613 else if(findX&&!findV){
614 phit->setCgemCluster(xCluster);
615 phit->setWire(xCluster->getclusterid());
616 phit->setFlag(xCluster->getflag());
617 }
618 else if(!findX&&findV){
619 phit->setCgemCluster(vCluster);
620 phit->setWire(vCluster->getclusterid());
621 phit->setFlag(vCluster->getflag());
622 }
623 }
624 }
625 }
626
627 // --- get mdcMcHit
628 SmartDataPtr<MdcMcHitCol> mdcMcHitCol(eventSvc(),"/Event/MC/MdcMcHitCol");
629 if(!mdcMcHitCol){
630 log << MSG::ERROR << "Could not retrieve MdcMcHitCol" << endreq;
631 }else{
632 MdcMcHitCol::iterator mdcMcHitIt = mdcMcHitCol->begin();
633 for(;mdcMcHitIt != mdcMcHitCol->end();mdcMcHitIt++,nMcHit++){
634 const MdcMcHit* mdcMcHit = (*mdcMcHitIt);
635 HoughHit hit(mdcMcHit,m_bunchT0,nMcHit);
636 m_mcHitCol.push_back(hit);
637 HoughHit* phit = &(*m_mcHitCol.rbegin());
638 for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end();hitIt++){
639 if(hitIt->getHitType()!=HoughHit::MDCHIT)continue;
640 const MdcDigi* mdcDigi = hitIt->getDigi();
641 if(mdcDigi->identify() == phit->getMdcMcHit()->identify()){
642 phit->setDigi(mdcDigi);
643 phit->setPairHit(&(*hitIt));
644 hitIt->setPairHit(phit);
645 }
646 }
647 }
648 }
649
650 for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end();hitIt++){
651 for(HitVector_Iterator mcHitIt = m_mcHitCol.begin(); mcHitIt != m_mcHitCol.end();mcHitIt++){
652 if(hitIt->getLayer()!=mcHitIt->getLayer())continue;
653 if(hitIt->getHitType() == HoughHit::CGEMHIT && mcHitIt->getHitType() == HoughHit::CGEMMCHIT){
654 if(mcHitIt->getCgemCluster()==NULL)continue;
655 int recClusterID(-1);
656 int mcClusterID(-2);
657 int recFlag = hitIt->getCgemCluster()->getflag();
658 int mcFlag = mcHitIt->getCgemCluster()->getflag();;
659 if(recFlag==mcFlag){
660 recClusterID = hitIt->getCgemCluster()->getclusterid();
661 mcClusterID = mcHitIt->getCgemCluster()->getclusterid();
662 }
663
664 if(recFlag==1&&mcFlag==2){
665 recClusterID = hitIt->getCgemCluster()->getclusterid();
666 mcClusterID = mcHitIt->getCgemCluster()->getclusterflage();
667 }
668
669 if(recFlag==0&&mcFlag==2){
670 recClusterID = hitIt->getCgemCluster()->getclusterid();
671 mcClusterID = mcHitIt->getCgemCluster()->getclusterflagb();
672 }
673
674 if(recFlag==2&&mcFlag==1){
675 recClusterID = hitIt->getCgemCluster()->getclusterflage();
676 mcClusterID = mcHitIt->getCgemCluster()->getclusterid();
677 }
678
679 if(recFlag==2&&mcFlag==0){
680 recClusterID = hitIt->getCgemCluster()->getclusterflagb();
681 mcClusterID = mcHitIt->getCgemCluster()->getclusterid();
682 }
683
684 if(recClusterID==mcClusterID){
685 hitIt->setPairHit(&(*mcHitIt));
686 hitIt->setCgemMcHit(mcHitIt->getCgemMcHit());
687 //mcHitIt->addPairHit(&(*hitIt));
688
689 //cout<<hitIt->getLayer()<<" "<<mcHitIt->getLayer()<<" "<<hitIt->getWire()<<" "<<mcHitIt->getWire()<<endl;
690 //cout<<hitIt->getDigi()->identify()<<" "<<mcHitIt->getMdcMcHit()->identify()<<endl;
691 //cout<<hitIt->getLayer()<<" "<<mcHitIt->getLayer()<<" ";
692 //cout<<hitIt->getWire()<<" "<<mcHitIt->getWire()<<endl;
693 //cout<<recFlag<<" "<<mcFlag<<" ";
694 //cout<<recClusterID<<" "<<mcClusterID<<endl;
695
696 //cout<<hitIt->getHitPosition();
697 //cout<<endl;
698 //cout<<mcHitIt->getHitPosition();
699 //cout<<endl;
700 //cout<<endl;
701 }
702 }
703 else if(hitIt->getHitType() == HoughHit::MDCHIT && mcHitIt->getHitType() == HoughHit::MDCMCHIT){
704 if(mcHitIt->getDigi()==NULL)continue;
705 //if(hitIt->getLayer() == mcHitIt->getLayer() && hitIt->getWire() == mcHitIt->getWire())
706 if(hitIt->getDigi()->identify() == mcHitIt->getDigi()->identify())
707 {
708 hitIt->setPairHit(&(*mcHitIt));
709 hitIt->setMdcMcHit(mcHitIt->getMdcMcHit());
710 //mcHitIt->addPairHit(&(*hitIt));
711 //cout<<hitIt->getLayer()<<" "<<hitIt->getWire()<<endl;
712 //cout<<mcHitIt->getLayer()<<" "<<mcHitIt->getWire()<<endl;
713 //cout<<endl;
714 //cout<<hitIt->getDigi()->identify()<<" "<<mcHitIt->getMdcMcHit()->identify()<<endl;
715 }
716 }
717 }
718 }
719 //for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end();hitIt++){
720 //if(hitIt->getHitType()!=HoughHit::MDCHIT)continue;
721 //cout<<hitIt->getLayer()<<" "<<hitIt->getWire()<<" ";
722 //if(hitIt->getPairHit()!=NULL){
723 //cout<<hitIt->getPairHit()->getHitID()<<" ";
724 //cout<<hitIt->getPairHit()->getLayer()<<" "<<hitIt->getPairHit()->getWire()<<" ";
725 //}
726 //cout<<endl;
727 //}
728
729 return nMcHit;
730}
731
732int HoughFinder::getMcParticleCol(){
733 const int maxCell[43] = {40,44,48,56, 64,72,80,80,
734 76,76,88,88, 100,100,112,112, 128,128,140,140,
735 160,160,160,160, 176,176,176,176, 208,208,208,208, 240,240,240,240,
736 256,256,256,256, 288,288,288 };
737 m_mcTrackCol.clear();
738 MsgStream log(msgSvc(), name());
739 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
740 if (!mcParticleCol) {
741 log << MSG::ERROR << "Could not find McParticle" << endreq;
742 }else{
743 bool psipDecay(false);
744 McParticleCol::iterator iter_mc = mcParticleCol->begin();
745 for (;iter_mc != mcParticleCol->end(); iter_mc++){
746 int trackIndex = (*iter_mc)->trackIndex();
747 int pid = (*iter_mc)->particleProperty();
748 bool primaryParticle = (*iter_mc)->primaryParticle();
749 bool leafParticle = (*iter_mc)->leafParticle();
750 bool decayFromGenerator = (*iter_mc)->decayFromGenerator();
751 bool decayInFlight = (*iter_mc)->decayInFlight();
752 HepLorentzVector initialPosition = (*iter_mc)->initialPosition();
753 HepLorentzVector initialFourMomentum = (*iter_mc)->initialFourMomentum();
754 //HepLorentzVector finalPosition = (*iter_mc)->finalPosition();
755 //HepLorentzVector finalFourMomentum = (*iter_mc)->finalFourMomentum();
756 unsigned int statusFlags = (*iter_mc)->statusFlags();
757 //string process = (*iter_mc)->getProcess();
758 //if(primaryParticle)continue;
759 //if(!decayFromGenerator)continue;
760 //if(pid==100443)psipDecay = true;
761 //if(!psipDecay) continue;
762 //if(!(fabs(pid)==11||fabs(pid)==211))continue;
763
764 string name;
765 int charge(0);
766 if( pid == 0 ) {
767 //cout << "Wrong particle id" <<endl;
768 continue;
769 }else{
770 IPartPropSvc* p_PartPropSvc;
771 static const bool CREATEIFNOTTHERE(true);
772 StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
773 if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) {
774 log << MSG::ERROR << "Could not initialize Particle Properties Service" << endreq;
775 }
776 HepPDT::ParticleDataTable* p_particleTable = p_PartPropSvc->PDT();
777 if( p_particleTable->particle(pid) ){
778 name = p_particleTable->particle(pid)->name();
779 charge = p_particleTable->particle(pid)->charge();
780 }else if( p_particleTable->particle(-pid) ){
781 name = "anti " + p_particleTable->particle(-pid)->name();
782 charge = (-1)*p_particleTable->particle(-pid)->charge();
783 }
784 }
785
786 //cout
787 //<<setw(8)<<trackIndex
788 //<<setw(12)<<name
789 //<<setw(8)<<pid
790 //<<setw(8)<<primaryParticle
791 //<<setw(8)<<leafParticle
792 //<<setw(8)<<decayFromGenerator
793 //<<setw(8)<<decayInFlight
794 //<<setw(8)<<statusFlags
795 //<<endl;
796
797 HoughTrack track(charge,initialPosition.v(),initialFourMomentum.v(),trackIndex);
798 //HoughTrack* track = new HoughTrack(charge,initialPosition.v(),initialFourMomentum.v(),trackIndex);
799 HepPoint3D pivot(0,0,0);
800 track.pivot(pivot);
801
802 int nHot(0);
803 vector<HoughHit*> hotList;
804 vector<HoughHit>& hitList = m_mcHitCol;
805 for(HitVector_Iterator hitIt = hitList.begin(); hitIt != hitList.end();hitIt++){
806 vector<int> trkID = hitIt->getTrkID();
807 if(trkID[0]==track.getTrkID()){
808 if(hitIt->getHitType()==HoughHit::CGEMMCHIT){
809 track.addHitPnt(&(*hitIt));
810 track.addVecStereoHitPnt(&(*hitIt));
811 hotList.push_back(&(*hitIt));
812 }
813 if(hitIt->getHitType()==HoughHit::MDCMCHIT){
814 if(hitIt->getFlag()==0)track.addHitPnt(&(*hitIt));
815 else track.addVecStereoHitPnt(&(*hitIt));
816 hotList.push_back(&(*hitIt));
817 }
818 nHot++;
819 }
820 }
821 if(nHot>0)m_mcTrackCol.push_back(track);
822
823 ///*
824 bool classifyHotByHotLayer(false);
825 bool classifyHotByTrackCenter(true);
826 if(classifyHotByHotLayer){
827 int lastHitLayer(0);
828 int deltaLayer(0);
829 int halfCircle(0);
830 HoughHit* lastLayerHit;
831 vector<HoughHit*> hitOnLayer;
832 int nMdcHot(0);
833 for(vector<HoughHit*>::iterator hitIt = hotList.begin(); hitIt != hotList.end(); hitIt++){
834 if((*hitIt)->getHitType()!=HoughHit::MDCHIT)continue;
835 if(nMdcHot==0){//the first hit on track
836 hitOnLayer.clear();
837 hitOnLayer.push_back(*hitIt);
838 lastLayerHit = *hitIt;
839 lastHitLayer = (*hitIt)->getLayer();
840 deltaLayer = 0;
841 halfCircle = 1;
842 }else{
843 if((*hitIt)->getLayer()==lastHitLayer){//find hits on the same layer
844 hitOnLayer.push_back(*hitIt);
845 lastHitLayer = (*hitIt)->getLayer();
846 }else{//another layer hit
847 if(deltaLayer*((*hitIt)->getLayer()-lastHitLayer)>=0){//not the turn layer
848 for(vector<HoughHit*>::iterator iter = hitOnLayer.begin(); iter != hitOnLayer.end(); iter++){
849 (*iter)->setHalfCircle(halfCircle);
850 lastLayerHit = *iter;
851 }
852 }else{//the turn layer
853 halfCircle++;
854 int maxWireGap(0);
855 int turnPoint(0);
856 int deltaWire(0);
857 for(vector<HoughHit*>::iterator iter = hitOnLayer.begin(); iter != hitOnLayer.end(); iter++){
858 deltaWire = fabs((*iter)->getWire()-lastLayerHit->getWire());
859 if(deltaWire>(maxCell[(*iter)->getLayer()])/2) deltaWire = maxCell[(*iter)->getLayer()] - deltaWire;
860 if(deltaWire>maxWireGap){
861 maxWireGap = deltaWire;
862 turnPoint = iter - hitOnLayer.begin();
863 }
864 lastLayerHit = *iter;
865 }
866 deltaWire = fabs((*hitIt)->getWire()-lastLayerHit->getWire());
867 if(deltaWire>(maxCell[(*hitIt)->getLayer()])/2) deltaWire = maxCell[(*hitIt)->getLayer()] - deltaWire;
868 if(deltaWire>maxWireGap){
869 maxWireGap = deltaWire;
870 turnPoint = hitOnLayer.size();
871 }
872
873 if(maxWireGap>2){
874 for(vector<HoughHit*>::iterator iter = hitOnLayer.begin(); iter != hitOnLayer.end(); iter++){
875 int shift = iter-hitOnLayer.begin();
876 if(shift<turnPoint)(*iter)->setHalfCircle(halfCircle-1);
877 else (*iter)->setHalfCircle(halfCircle);
878 }
879 }else{
880 for(vector<HoughHit*>::iterator iter = hitOnLayer.begin(); iter != hitOnLayer.end(); iter++){
881 int shift = iter-hitOnLayer.begin();
882 if(shift<hitOnLayer.size()/2)(*iter)->setHalfCircle(halfCircle-1);
883 else (*iter)->setHalfCircle(halfCircle);
884 }
885 }
886 }
887 hitOnLayer.clear();
888 hitOnLayer.push_back(*hitIt);
889 lastHitLayer = (*hitIt)->getLayer();
890 deltaLayer = (*hitIt)->getLayer()-lastHitLayer;
891 }
892
893 }
894 nMdcHot++;
895 }
896 }
897 //*/
898 ///*
899 if(classifyHotByTrackCenter){
900 int cgemHalf(0), mdcHalf(0);
901 int cgemSign(0), mdcSign(0);
902 int half(0), sign(0);
903 for(vector<HoughHit*>::iterator hitIt = hotList.begin(); hitIt != hotList.end(); hitIt++){
904 //(*hitIt)->print();
905 if((*hitIt)->getHitType() == HoughHit::CGEMMCHIT){
906 if(track.judgeHalf(*hitIt) != cgemSign){
907 cgemHalf++;
908 cgemSign = track.judgeHalf(*hitIt);
909 half = cgemHalf;
910 sign = cgemSign;
911 }
912 }
913 else if((*hitIt)->getHitType() == HoughHit::MDCMCHIT){
914 if(track.judgeHalf(*hitIt) != mdcSign){
915 mdcHalf++;
916 mdcSign = track.judgeHalf(*hitIt);
917 half = mdcHalf;
918 sign = mdcSign;
919 }
920 }
921 (*hitIt)->setHalfCircle(half);
922 }
923 }
924 if(hotList.size()>0){
925 if(printFlag)track.print();
926 //track.printHot();
927 //cout<<endl;
928 }
929 for(vector<HoughHit*>::iterator hotIt = hotList.begin(); hotIt != hotList.end();hotIt++){
930 HoughHit* hitIt = *hotIt;
931 if(printFlag)hitIt->print();
932 //cout<<hitIt->getHitID();
933 //cout<<" ("<<hitIt->getLayer()<<", "<<hitIt->getWire()<<") "<<hitIt->getFlag();
934 //cout<<endl;
935 }
936 if(printFlag)if(hotList.size())cout<<endl;
937 }
938 sort(m_mcTrackCol.begin(),m_mcTrackCol.end(),largerPt);
939 //cout<<"N MCtrack: "<<m_mcTrackCol.size()<<endl;
940 for(vector<HoughTrack>::iterator trkIter=m_mcTrackCol.begin();trkIter!=m_mcTrackCol.end();trkIter++){
941 //trkIter->print();
942 }
943 }
944 return m_mcTrackCol.size();
945}
946
947int HoughFinder::fillHistogram(HoughHit* hit, TH2D* hitMap, int charge, int vote)
948{
949 int nBin = hitMap->GetXaxis()->GetNbins();
950 double xMin = hitMap->GetXaxis()->GetXmin();
951 double xMax = hitMap->GetXaxis()->GetXmax();
952 double yMin = hitMap->GetYaxis()->GetXmin();
953 double yMax = hitMap->GetYaxis()->GetXmax();
954 double dx = (xMax-xMin)/(nBin*vote)>1e-6?(xMax-xMin)/(nBin*vote):1e-6;
955 double x = xMin + dx/2;
956
957 int flag = hit->getFlag();
958 int nPoint(0);
959 if(flag==0){
960 double xHit = hit->getHitPosition().x();
961 double yHit = hit->getHitPosition().y();
962 double rHit = hit->getDriftDist();//drift distance
963 double denominator = xHit*xHit+yHit*yHit-rHit*rHit;
964
965 //double X1(0),Y1(0),X2(0),Y2(0),R(0);
966 double X1 = 2*xHit/denominator;
967 double Y1 = 2*yHit/denominator;
968 double X2 = 2*xHit/denominator;
969 double Y2 = 2*yHit/denominator;
970 double R = 2*rHit/denominator;
971
972 while(x<xMax){
973 double y1 = X1*cos(x) + Y1*sin(x) + R;
974 double y2 = X2*cos(x) + Y2*sin(x) - R;
975 double slope1 = -X1*sin(x) + Y1*cos(x);
976 double slope2 = -X2*sin(x) + Y2*cos(x);
977 double hitOnCircle1 = charge*slope1*y1;
978 double hitOnCircle2 = charge*slope2*y2;
979
980 bool cut(true);
981 cut = yMin<y1 && y1<yMax;
982 cut = cut && hitOnCircle1 <=0;
983 if(cut){
984 hitMap->Fill(x,y1);
985 nPoint++;
986 }
987
988 cut = yMin<y2 && y2<yMax;
989 cut = cut && hitOnCircle2 <= 0;
990 if(cut){
991 //if(hit.getHitType() == HoughHit::CGEMHIT)continue;// comment by llwang
992 //int bin = hitMap->FindBin(x,y2);
993 //if(hitMap->GetBinContent(bin)>0)continue;
994 hitMap->Fill(x,y2);
995 nPoint++;
996 }
997 x += dx;
998 }
999 }else{
1000 vector<HoughHit::S_Z> sz = hit->getSZ();
1001 vector<HoughHit::S_Z>::iterator iter = sz.begin();
1002 for(;iter!=sz.end();iter++){
1003 double S = iter->first;
1004 double Z = iter->second;
1005 double dy = -S*dx;
1006 double y = -dy + Z;
1007 while(x<xMax){
1008 x += dx;
1009 y += dy;
1010 //y = S*cos(x) + Z*sin(x);
1011 bool cut(true);
1012 cut = yMin<y && y<yMax;
1013 if(cut){
1014 hitMap->Fill(x,y);
1015 nPoint++;
1016 }
1017 }
1018 }
1019 }
1020 return nPoint;
1021}
1022
1023int HoughFinder::fillHistogram(HoughHit* hit, TH2D* hitMap, HoughTrack* trkCandi, int vote)
1024{
1025 int charge = trkCandi->getCharge();
1026 int nBin = hitMap->GetXaxis()->GetNbins();
1027 double xMin = hitMap->GetXaxis()->GetXmin();
1028 double xMax = hitMap->GetXaxis()->GetXmax();
1029 double yMin = hitMap->GetYaxis()->GetXmin();
1030 double yMax = hitMap->GetYaxis()->GetXmax();
1031 double dx = (xMax-xMin)/(nBin*vote)>1e-6?(xMax-xMin)/(nBin*vote):1e-6;// step size
1032 double x = xMin + dx/2;// first step
1033
1034 int flag = hit->getFlag();
1035 int nPoint(0);
1036 if(flag==0) // X-view hits
1037 {
1038 bool firstHalf = trkCandi->judgeHalf(hit)==1;
1039 if(!firstHalf) return 0;
1040 double xHit = hit->getHitPosition().x();
1041 double yHit = hit->getHitPosition().y();
1042 double rHit = hit->getDriftDist();//drift distance
1043 double denominator = xHit*xHit+yHit*yHit-rHit*rHit;
1044
1045 //double X1(0),Y1(0),X2(0),Y2(0),R(0);
1046 double X1 = 2*xHit/denominator;
1047 double Y1 = 2*yHit/denominator;
1048 double X2 = 2*xHit/denominator;
1049 double Y2 = 2*yHit/denominator;
1050 double R = 2*rHit/denominator;
1051
1052 while(x<xMax){
1053 double y1 = X1*cos(x) + Y1*sin(x) + R;
1054 double y2 = X2*cos(x) + Y2*sin(x) - R;
1055 double slope1 = -X1*sin(x) + Y1*cos(x);
1056 double slope2 = -X2*sin(x) + Y2*cos(x);
1057 double hitOnCircle1 = charge*slope1*y1;
1058 double hitOnCircle2 = charge*slope2*y2;
1059
1060 bool cut(true);
1061 cut = yMin<y1 && y1<yMax;
1062 //cut = cut && hitOnCircle1 <=0;
1063 if(cut){
1064 hitMap->Fill(x,y1);
1065 nPoint++;
1066 }
1067
1068 cut = yMin<y2 && y2<yMax;
1069 //cut = cut && hitOnCircle2 <= 0;
1070 if(cut){
1071 //if(hit.getHitType() == HoughHit::CGEMHIT)continue;// comment by llwang
1072 //int bin = hitMap->FindBin(x,y2);
1073 //if(hitMap->GetBinContent(bin)>0)continue;
1074 hitMap->Fill(x,y2);
1075 nPoint++;
1076 }
1077 x += dx;
1078 }
1079 }else{ // V-view hits
1080 //if(hit.getPairHit()==NULL)return nPoint;// FIXME
1081 //if(hit.getPairHit()->getHalfCircle()!=1)return nPoint;// FIXME
1082 int nTangency = trkCandi->calculateZ_S(hit);
1083 //cout<<nTangency<<endl;
1084 if(nTangency==0){
1085 //hit.setUse(-1);// why? FIXME
1086 }
1087 else
1088 {
1089 vector<HoughHit::S_Z> sz = hit->getSZ();
1090 vector<HoughHit::S_Z>::iterator iter = sz.begin();
1091 for(;iter!=sz.end();iter++)
1092 {
1093 double S = iter->first;
1094 double Z = iter->second;
1095 while(x<xMax)
1096 {
1097 double y = -S*x + Z;
1098 //y = S*cos(x) + Z*sin(x);
1099 bool cut(true);
1100 cut = yMin<y && y<yMax;
1101 if(cut)
1102 {
1103 hitMap->Fill(x,y);
1104 nPoint++;
1105 }
1106 x += dx;
1107 }
1108 }
1109 }
1110 }
1111 return nPoint;
1112}
1113
1114int HoughFinder::fillHistogram(vector<HoughHit*> &hitList, TH2D* hitMap, int charge, int vote, HoughTrack* trkCandi )
1115{
1116 int nHitsFilled = 0;
1117 m_VHoughHitListOnSZmap.clear();
1118 std::vector<HoughHit*>::iterator iter = hitList.begin();
1119 for(; iter!= hitList.end(); iter++)
1120 {
1121 int used = (*iter)->getUse();
1122 //(*iter)->print();
1123 if(used==0) {
1124 int nPoint = 0;
1125 if(trkCandi!=NULL){
1126 //if((*iter)->getPairHit()==NULL) continue;// FIXME
1127 //if((*iter)->getPairHit()->getHalfCircle()!=1)continue;// FIXME
1128 nPoint=fillHistogram((*iter),hitMap,trkCandi,vote);
1129 }
1130 else{
1131 //if((*iter)->getPairHit()==NULL) continue;// FIXME
1132 //if((*iter)->getPairHit()->getHalfCircle()!=1)continue;// FIXME
1133 nPoint=fillHistogram(*iter,hitMap,charge,vote);
1134 }
1135 if(nPoint>0) {
1136 nHitsFilled++;
1137 m_VHoughHitListOnSZmap.push_back(*iter);
1138 }
1139 }
1140 }
1141 return nHitsFilled;
1142}
1143
1144//////////////////////////////////////////////////---Hough Transform---//////////////////////////////////////////////////
1145/*
1146//Depth First Search Peak
1147int HoughFinder::findTrack(vector<HoughHit> &hitList, vector<HoughTrack*>& trackVector, double x, double y, double dx, double dy, int charge, int mapType, int step)
1148{
1149 step++;
1150 if(mapType==0){
1151 m_nStep = m_nStep_xy;
1152 m_xBin = m_xBin_xy;
1153 m_yBin = m_yBin_xy;
1154 m_xExtend = m_xExtend_xy;
1155 m_yExtend = m_yExtend_xy;
1156 m_nVote = m_nVote_xy;
1157 m_nRMS = m_nRMS_xy;
1158 m_xInclude = m_xInclude_xy;
1159 m_yInclude = m_yInclude_xy;
1160 }else{
1161 m_nStep = m_nStep_sz;
1162 m_xBin = m_xBin_sz;
1163 m_yBin = m_yBin_sz;
1164 m_xExtend = m_xExtend_sz;
1165 m_yExtend = m_yExtend_sz;
1166 m_nVote = m_nVote_sz;
1167 m_nRMS = m_nRMS_sz;
1168 m_xInclude = m_xInclude_sz;
1169 m_yInclude = m_yInclude_sz;
1170 }
1171
1172 // --- build maps
1173 map<int,TH2D*> hitMapList;
1174 TH2D* houghMap = buildMap(hitList,hitMapList,x,y,dx,dy,charge,mapType,step);
1175 //cout<<"hitMapList.size()+++"<<hitMapList.size()<<" "<<charge<<endl;
1176 // --- find peaks
1177 findPeak(hitList,hitMapList,houghMap,trackVector,charge,mapType,step);
1178 // --- clear maps
1179 clearMap(houghMap,hitMapList);
1180 //cout<<"hitMapList.size()---"<<hitMapList.size()<<endl;
1181 //cout<<endl<<endl;
1182 return trackVector.size();
1183}
1184
1185
1186TH2D* HoughFinder::buildMap(vector<HoughHit> &hitList, map<int,TH2D*>& hitMapList, double x, double y, double dx, double dy, int charge, int mapType, int step)
1187{
1188 int vote = m_nVote[step-1];
1189 int xBin = m_xBin[step-1];
1190 int yBin = m_yBin[step-1];
1191 double xExtend = 0.5+m_xExtend[step-1];
1192 double yExtend = 0.5+m_yExtend[step-1];
1193 if(m_findPeakMethod==1){
1194 vote = m_nVote[step-m_nStep-1];
1195 xBin = m_xBin[step-m_nStep-1];
1196 yBin = m_yBin[step-m_nStep-1];
1197 xExtend = 0.5+m_xExtend[step-m_nStep-1];
1198 yExtend = 0.5+m_yExtend[step-m_nStep-1];
1199 }
1200 double xMin = x-dx*xExtend;
1201 double xMax = x+dx*xExtend;
1202 double yMin = y-dy*yExtend;
1203 double yMax = y+dy*yExtend;
1204 stringstream ssname;
1205 ssname<<"step"<<step;
1206 string sname = ssname.str();
1207 const char * name = sname.c_str();
1208 TH2D* houghMap = new TH2D(name, name, xBin, xMin, xMax, yBin, yMin, yMax);
1209 HitVector_Iterator hitIt = hitList.begin();
1210 for(; hitIt != hitList.end();hitIt++){
1211 //hitIt->print();
1212 if(hitIt->getFlag()==1 && hitIt->getHitType() == HoughHit::CGEMHIT)continue;
1213 if(mapType==0 && hitIt->getFlag() != 0)continue;
1214 if(mapType!=0 && hitIt->getFlag() == 0)continue;
1215 int hitID = hitIt->getHitID();
1216 stringstream ssname;
1217 ssname<<"step"<<step<<",hitID"<<hitID;
1218 string sname = ssname.str();
1219 const char * name = sname.c_str();
1220 TH2D* hitMap = new TH2D(name, name, xBin, xMin, xMax, yBin, yMin, yMax);
1221 int nPoint = fillHistogram((*hitIt),hitMap,charge,vote);
1222 if(nPoint>0){
1223 houghMap->Add(hitMap);
1224 hitMapList[hitID] = hitMap;
1225 }else delete hitMap;
1226 }
1227 return houghMap;
1228}
1229
1230int HoughFinder::findPeak(vector<HoughHit> &hitList, map<int,TH2D*>& hitMapList, TH2D* houghMap, vector<HoughTrack*>& trackVector, int charge, int mapType, int step)
1231{
1232 int trkID = trackVector.size();
1233 int binx(0),biny(0),binz(0);
1234 houghMap->GetMaximumBin(binx,biny,binz);
1235 int peakHeight = houghMap->GetBinContent(binx,biny);
1236 int nRMS = m_nRMS[step-1];
1237 if(m_findPeakMethod==1)nRMS = m_nRMS[step-m_nStep-1];
1238 double peakCut = mapDev(houghMap,nRMS);
1239 int nPeaks(0);
1240 int nRemove(0);
1241 while(peakHeight>peakCut){
1242 houghMap->GetMaximumBin(binx,biny,binz);
1243 double x = houghMap->GetXaxis()->GetBinCenter(binx);
1244 double y = houghMap->GetYaxis()->GetBinCenter(biny);
1245 double dx = houghMap->GetXaxis()->GetBinWidth(binx);
1246 double dy = houghMap->GetYaxis()->GetBinWidth(biny);
1247 //cout<<x<<" "<<y<<endl;
1248 centroid(houghMap,m_xInclude,m_yInclude,x,y);
1249 //cout<<x<<" "<<y<<endl;
1250 if(step<m_nStep){
1251 findTrack(hitList,trackVector,x,y,dx,dy,charge,mapType,step);
1252 nRemove = removeHitMap(hitMapList,houghMap,trackVector);
1253 if(nRemove==0)nRemove = removeHitMap(hitMapList,houghMap,binx,biny);
1254 }else{
1255 if(mapType == 0){
1256 HoughTrack* track = new HoughTrack(charge,x,y,dx,dy,trkID);
1257 //track->print();
1258 int nHot = track->findHot(hitList,mapType);
1259 if(nHot>=3){
1260 trackVector.push_back(track);
1261 nRemove = removeHitMap(hitMapList,houghMap,track);
1262 if(nRemove==0)nRemove = removeHitMap(hitMapList,houghMap,binx,biny);
1263 trkID++;
1264 }else{
1265 nRemove = removeHitMap(hitMapList,houghMap,binx,biny);
1266 }
1267
1268 }else if(trackVector.size()>0){
1269 HoughTrack* track = *(trackVector.begin());
1270 double tanl = track->tanl();
1271 double dz = track->dz();
1272 double dTanl = track->getDTanl();
1273 double dDz = track->getDDz();
1274 vector<HoughHit> hot;
1275 int nHot1 = track->findHot(hitList,hot,mapType);
1276 track->setTanl(x);
1277 track->setDz(y);
1278 track->setDTanl(dx);
1279 track->setDDz(dy);
1280 track->updateHelix();
1281 int nHot2 = track->findHot(hitList,hot,mapType);
1282 if(nHot1>nHot2){
1283 track->setTanl(tanl);
1284 track->setDz(dz);
1285 track->setDTanl(dTanl);
1286 track->setDDz(dDz);
1287 track->updateHelix();
1288 }
1289 nRemove = removeHitMap(hitMapList,houghMap,binx,biny);
1290 }
1291 }
1292 //double xKurtosis = houghMap->GetKurtosis(1);
1293 //double yKurtosis = houghMap->GetKurtosis(2);
1294 houghMap->GetMaximumBin(binx,biny,binz);
1295 peakHeight = houghMap->GetBinContent(binx,biny);
1296 peakCut = mapDev(houghMap,nRMS);
1297 nPeaks++;
1298 //if(mapType!=0 && nPeaks>0)break;
1299 if(nRemove==0)break;
1300 }
1301 return nPeaks;
1302}
1303
1304//remove map of hits that on the peak cell
1305int HoughFinder::removeHitMap(map<int,TH2D*>& hitMapList, TH2D* houghMap,int binx, int biny)
1306{
1307 bool debug(0);
1308 int nRemove(0);
1309 //if(debug){
1310 //cout<<endl;
1311 //for(map<int,TH2D*>::iterator iter = hitMapList.begin();iter != hitMapList.end();iter++){
1312 //cout<<"hit->"<<iter->first<<" ";
1313 //if(iter->second==NULL)cout<<"erased"<<endl;
1314 //else cout<<iter->second->GetBinContent(binx,biny)<<endl;
1315 //}
1316 //cout<<endl;
1317 //}
1318
1319 //int binx(0),biny(0),binz(0);
1320 //houghMap->GetMaximumBin(binx,biny,binz);
1321 map<int,TH2D*>::iterator iter = hitMapList.begin();
1322 for(;iter != hitMapList.end();iter++){
1323 if(debug)cout<<"hit->"<<iter->first<<" ";
1324 if(iter->second!=NULL){
1325 if(iter->second->GetBinContent(binx,biny)>0){
1326 houghMap->Add(iter->second,-1);
1327 delete iter->second;
1328 iter->second = NULL;
1329 //hitMapList.erase(iter);
1330 //hitMapList.erase(iter->first);
1331 nRemove++;
1332 if(debug)cout<<" erasing";
1333 }
1334 }
1335 else{
1336 if(debug)cout<<" erased";
1337 }
1338 if(debug)cout<<endl;
1339 }
1340
1341 //if(debug){
1342 //cout<<endl;
1343 //for(map<int,TH2D*>::iterator iter = hitMapList.begin();iter != hitMapList.end();iter++){
1344 //cout<<"hit->"<<iter->first<<" ";
1345 //if(iter->second==NULL)cout<<"erased"<<endl;
1346 //else cout<<iter->second->GetBinContent(binx,biny)<<endl;
1347 //}
1348 //cout<<endl;
1349 //}
1350
1351 return nRemove;
1352}
1353
1354//remove map of hits that on the track
1355int HoughFinder::removeHitMap(map<int,TH2D*>& hitMapList, TH2D* houghMap, HoughTrack* track)
1356{
1357 bool debug(1);
1358 int nRemove(0);
1359 //if(debug){
1360 //cout<<endl;
1361 //for(map<int,TH2D*>::iterator iter = hitMapList.begin();iter != hitMapList.end();iter++){
1362 //cout<<"hit->"<<iter->first<<" ";
1363 //if(iter->second==NULL)cout<<"erased"<<endl;
1364 //else{
1365 //int bin = iter->second->FindBin(track->getAngle(),track->getRho());
1366 //cout<<iter->second->GetBinContent(bin)<<endl;
1367 //}
1368 //}
1369 //cout<<endl;
1370 //}
1371
1372 vector<HoughHit>* hot = track->getHot();
1373 HitVector_Iterator hotIt = hot->begin();
1374 for(;hotIt!=hot->end();hotIt++){
1375 int hitID = hotIt->getHitID();
1376 if(debug)cout<<"hit->"<<hitID<<" ";
1377 map<int,TH2D*>::iterator iter = hitMapList.find(hitID);
1378 if(iter->second!=NULL){
1379 if(iter != hitMapList.end()){
1380 houghMap->Add(iter->second,-1);
1381 delete iter->second;
1382 iter->second = NULL;
1383 //hitMapList.erase(iter);
1384 //hitMapList.erase(iter->first);
1385 nRemove++;
1386 if(debug)cout<<"erasing"<<endl;
1387 }
1388 }else{
1389 if(debug)cout<<"erased"<<endl;
1390 }
1391 }
1392
1393 //if(debug){
1394 //cout<<endl;
1395 //for(map<int,TH2D*>::iterator iter = hitMapList.begin();iter != hitMapList.end();iter++){
1396 //cout<<"hit->"<<iter->first<<" ";
1397 //if(iter->second==NULL)cout<<"erased"<<endl;
1398 //else{
1399 //int bin = iter->second->FindBin(track->getAngle(),track->getRho());
1400 //cout<<iter->second->GetBinContent(bin)<<endl;
1401 //}
1402 //}
1403 //cout<<endl;
1404 //}
1405 return nRemove;
1406}
1407
1408//remove map of hits that on all tracks
1409int HoughFinder::removeHitMap(map<int,TH2D*>& hitMapList, TH2D* houghMap, vector<HoughTrack*>& trackVector)
1410{
1411 int nRemove(0);
1412 TrackVector_Iterator trkIter = trackVector.begin();
1413 for(;trkIter!=trackVector.end();trkIter++){
1414 if((*trkIter)->getFlag() == 1)continue;
1415 vector<HoughHit>* hot = (*trkIter)->getHot();
1416 HitVector_Iterator iter = hot->begin();
1417 for(;iter!=hot->end();iter++){
1418 int hitID = iter->getHitID();
1419 map<int,TH2D*>::iterator iter = hitMapList.find(hitID);
1420 if(iter != hitMapList.end()){
1421 houghMap->Add(iter->second,-1);
1422 delete iter->second;
1423 //iter->second = NULL;
1424 hitMapList.erase(iter);
1425 nRemove++;
1426 }
1427 }
1428 }
1429 //for(map<int,TH2D*>::iterator iter = hitMapList.begin();iter != hitMapList.end();iter++){
1430 //if(iter->second!=NULL)cout<<iter->first<<" "<<iter->second->GetBinContent(binx,biny)<<endl;
1431 //else cout<<iter->first<<" "<<" erase"<<endl;
1432 //}
1433 //cout<<endl;
1434 return nRemove;
1435}
1436
1437int HoughFinder::clearMap(TH2D* houghMap, map<int,TH2D*>& hitMapList)
1438{
1439 int nRemove(0);
1440 map<int,TH2D*>::iterator iter = hitMapList.begin();
1441 for(;iter != hitMapList.end();iter++){
1442 if(iter->second != NULL){
1443 delete iter->second;
1444 iter->second = NULL;
1445 //hitMapList.erase(iter);
1446 nRemove++;
1447 }
1448 }
1449 hitMapList.clear();
1450
1451 delete houghMap;
1452 houghMap = NULL;
1453
1454 return nRemove;
1455}
1456
1457
1458int HoughFinder::centroid(TH2D* houghMap, int xInclude, int yInclude, double &xc, double &yc)
1459{
1460 double X(0),Y(0),W(0);
1461 int binx(0),biny(0),binz(0);
1462 houghMap->GetMaximumBin(binx,biny,binz);
1463 int xMin = binx-xInclude>0?binx-xInclude:0;
1464 int xMax = binx+xInclude<houghMap->GetNbinsX()?binx+xInclude:houghMap->GetNbinsX();
1465 int yMin = biny-yInclude>0?biny-yInclude:0;
1466 int yMax = biny+yInclude<houghMap->GetNbinsY()?biny+yInclude:houghMap->GetNbinsY();
1467 //cout<<houghMap->GetNbinsX()<<" "<<houghMap->GetNbinsY()<<endl;
1468 for(int i=binx-xInclude;i<=binx+xInclude;i++){
1469 for(int j=biny-yInclude;j<=biny+yInclude;j++){
1470 double x = houghMap->GetXaxis()->GetBinCenter(i);
1471 double y = houghMap->GetYaxis()->GetBinCenter(j);
1472 int w = houghMap->GetBinContent(i,j);
1473 X += w*x;
1474 Y += w*y;
1475 W += w;
1476 }
1477 }
1478 if(W>0){
1479 xc = X/W;
1480 yc = Y/W;
1481 }
1482 return W;
1483}
1484
1485double HoughFinder::mapDev(TH2D* houghMap,double nRMS)
1486{
1487 int xBins = houghMap->GetXaxis()->GetNbins();
1488 int yBins = houghMap->GetYaxis()->GetNbins();
1489 int count_use=0;
1490 double sum=0;
1491 double sum2=0;
1492 for(int i=0;i<xBins;i++){
1493 for(int j=0;j<yBins;j++){
1494 int count=houghMap->GetBinContent(i+1,j+1);
1495 if(count!=0){
1496 count_use++;
1497 sum+=count;
1498 sum2+=count*count;
1499 }
1500 }
1501 }
1502 double aver(0),rms(0);
1503 if(count_use != 0){
1504 aver=sum/count_use;
1505 if(count_use==1) rms = sqrt(sum2/count_use-aver*aver);
1506 else rms = sqrt((sum2-count_use*aver*aver)/(count_use-1));
1507 }
1508 double cut = aver + nRMS*rms;
1509 return cut;
1510}
1511//*/
1512/*
1513//Breadth First Search Peak
1514int HoughFinder::findTrack(vector<HoughTrack>& trackVector, int charge, int mapType, int step)
1515{
1516 step++;
1517 vector<HoughTrack>& trkVector;
1518 TrackVector_Iterator trkIter = trackVector.begin();
1519 for(;trkIter!=trackVector.end();trkIter++){
1520 int xBin = 10;
1521 int yBin = 10;
1522 double xPeak = trkIter->getAlpha();
1523 double yPeak = trkIter->getRho();
1524 double dx(0),dy(0);
1525 if(mapType==0){
1526 dx = trkIter->getDAngle();
1527 dy = trkIter->getDRho();
1528 }else{
1529 dx = trkIter->getDTanl();
1530 dy = trkIter->getDDz();
1531 }
1532 findTrack(trackVector,xPeak-dx,xPeak+dx,yPeak-dy,yPeak+dy,charge,mapType,m_nStep+step);
1533 }
1534 trackVector.swap(trkVector);
1535 if(step<m_nStep)findTrack(trackVector);
1536}
1537
1538//Breadth First Search Peak
1539int HoughFinder::findTrack(vector<HoughTrack>& trackVector, int charge, int mapType, int step)
1540{
1541 while(step<m_nStep){
1542 step++;
1543 vector<HoughTrack>& trkVector;
1544 trackVector.swap(trkVector);
1545 trackVector.clear();
1546 TrackVector_Iterator trkIter = trkVector.begin();
1547 for(;trkIter!=trkVector.end();trkIter++){
1548 int xBin = 10;
1549 int yBin = 10;
1550 double xPeak = trkIter->getAlpha();
1551 double yPeak = trkIter->getRho();
1552 double dx(0),dy(0);
1553 if(mapType==0){
1554 dx = trkIter->getDAngle();
1555 dy = trkIter->getDRho();
1556 }else{
1557 dx = trkIter->getDTanl();
1558 dy = trkIter->getDDz();
1559 }
1560 findTrack(trackVector,xPeak-dx,xPeak+dx,yPeak-dy,yPeak+dy,charge,mapType,m_nStep+step);
1561 }
1562 }
1563 //trackVector.swap(trkVector);
1564 if(step<m_nStep)findTrack(trackVector);
1565}
1566//*/
1567
1568bool moreHot(HoughTrack trk1, HoughTrack trk2)
1569{
1570 int nHot1 = trk1.getHotList().size();
1571 int nHot2 = trk2.getHotList().size();
1572 return nHot1>nHot2;
1573 //return trk1.getHot().size() > trk2.getHot().size();
1574}
1575
1576//check track by hots, remove the shorter track with too many shared hits
1577// --- v1 ---
1578/*
1579int HoughFinder::checkHot(vector<HoughTrack>& trackVector)
1580{
1581 int nDelete(0);
1582 std::sort(trackVector.begin(),trackVector.end(),moreHot);// sort track candidates from more hits to less
1583 for(vector<HoughTrack>::iterator trkIT1 = trackVector.begin(); trkIT1!=trackVector.end(); trkIT1++){
1584 if((trkIT1)->getFlag() == 1)continue;
1585 vector<HoughHit>* hot1 = (trkIT1)->getHot();
1586 trkIT1->dropRedundentCgemXHits();
1587 //cout<<"hit size "<<hot1->size()<<endl;
1588 if(hot1->size()<3){
1589 (trkIT1)->setFlag(1);
1590 trkIT1->releaseSelHits();
1591 continue;
1592 }
1593
1594 for(vector<HoughTrack>::iterator trkIT2 = trkIT1+1; trkIT2!=trackVector.end();trkIT2++){
1595 if((trkIT2)->getFlag() == 1)continue;
1596 vector<HoughHit>* hot2 = (trkIT2)->getHot();
1597 if(hot2->size()<3){
1598 (trkIT2)->setFlag(1);
1599 trkIT2->releaseSelHits();
1600 continue;
1601 }
1602
1603 int shareHit(0);
1604 for(HitVector_Iterator hitIt1 = hot1->begin(); hitIt1 != hot1->end(); hitIt1++){
1605 for(HitVector_Iterator hitIt2 = hot2->begin(); hitIt2 != hot2->end(); hitIt2++){
1606 //cout<<"id2, id1 = "<<hitIt2->getHitID()<<", "<<hitIt1->getHitID()<<endl;
1607 if(hitIt2->getHitID()==hitIt1->getHitID())shareHit++;
1608 }
1609 }
1610 //cout<<"shareHit = "<<shareHit<<endl;
1611 if(shareHit>m_shareHitRate*hot2->size()){
1612 //delete *trkIT2;
1613 //trackVector.erase(trkIT2);
1614 //trkIT2--;
1615 (trkIT2)->setFlag(1);
1616 trkIT2->releaseSelHits();
1617 nDelete++;
1618 cout<<"too many shared hits ("<<shareHit<<"/"<<hot2->size()
1619 <<") in track "<<(trkIT2)->getTrkID()
1620 <<" with track "<<trkIT1->getTrkID()<<endl;
1621 }
1622 }
1623 }
1624 return nDelete;
1625}
1626//*/
1627
1628// --- v2 ---
1629int HoughFinder::checkHot(vector<HoughTrack>& trackVector)
1630{
1631 //m_debug=true;
1632 int nDelete(0);
1633 if(trackVector.size()==0)return 0;
1634 std::sort(trackVector.begin(),trackVector.end(),moreHot);// sort track candidates from more hits to less
1635
1636 vector<HoughTrack>::iterator trkIT1 = trackVector.end();
1637 //vector<HoughTrack>::iterator trkIT1 = trackVector.rbegin();
1638 bool loop = true;
1639
1640 while(loop)
1641 {
1642 trkIT1--;
1643 if(trkIT1==trackVector.begin()) loop=false;
1644 if((trkIT1)->getFlag() == 1) continue;
1645 trkIT1->dropRedundentCgemXHits();
1646 //cout<<"hit size "<<hot1->size()<<endl;
1647 int nHits = trkIT1->getVecHitPnt().size();
1648 if(nHits<3){
1649 (trkIT1)->setFlag(1);
1650 trkIT1->releaseSelHits();
1651 if(m_debug) {
1652 cout<<"too less hits ("<<nHits<<") in track "<<(trkIT1)->getTrkID()<<endl;
1653 }
1654 continue;
1655 }
1656 int nShared = trkIT1->getNHitsShared();
1657
1658 //cout<<"shareHit = "<<shareHit<<endl;
1659 if(nShared>m_shareHitRate*nHits){
1660 trkIT1->setFlag(1);
1661 trkIT1->releaseSelHits();
1662 nDelete++;
1663 if(m_debug) {
1664 cout<<"too many shared hits ("<<nShared<<"/"<<nHits
1665 <<") in track "<<(trkIT1)->getTrkID()
1666 <<endl;
1667 }
1668 }
1669 }
1670 //m_debug=false;
1671 return nDelete;
1672}
1673
1674int HoughFinder::checkTrack(vector<HoughTrack>& trackVector)
1675{
1676 int nDelete(0);
1677 double dDr(999);
1678 double dPhi0(999);
1679 double dKappa(999);
1680 double dDz(999);
1681 double dTanl(999);
1682
1683 //m_debug=true;
1684
1685 std::sort(trackVector.begin(),trackVector.end(),moreHot);// sort track candidates from more hits to less
1686 for(vector<HoughTrack>::iterator trkIT1 = trackVector.begin(); trkIT1!=trackVector.end(); trkIT1++){
1687 if((trkIT1)->getFlag() == 1)continue;
1688 int charge = (trkIT1)->getCharge();
1689 double dr = (trkIT1)->dr();
1690 double phi0 = (trkIT1)->phi0();
1691 double kappa = (trkIT1)->kappa();
1692 double dz = (trkIT1)->dz();
1693 double tanl = (trkIT1)->tanl();
1694 for(vector<HoughTrack>::iterator trkIT2 = trkIT1+1; trkIT2!=trackVector.end();trkIT2++){
1695 if((trkIT2)->getFlag() == 1)continue;
1696 bool sameTrack(false);
1697 dDr = fabs(dr - (trkIT2)->dr());
1698 dPhi0 = fabs(phi0 - (trkIT2)->phi0());
1699 if((trkIT2)->getCharge() == charge){
1700 dKappa = fabs(kappa - (trkIT2)->kappa());
1701 dTanl = fabs(tanl - (trkIT2)->tanl());
1702 }else{
1703 dKappa = fabs(fabs(kappa) - fabs((trkIT2)->kappa()));
1704 dTanl = fabs(fabs(tanl) - fabs((trkIT2)->tanl()));
1705 }
1706 double r = fabs(dz)<fabs((trkIT2)->dz())?(trkIT1)->radius():(trkIT2)->radius();
1707 double k = dz<fabs((trkIT2)->dz())?tanl:(trkIT2)->tanl();
1708 dDz = fabs(dz - (trkIT2)->dz());
1709 int nCircle = fabs(dDz)/(2*M_PI*r*k);
1710 dDz = fabs(dDz - nCircle*(2*M_PI*r*k));
1711 if(dDz>M_PI*r*k)dDz = fabs(dDz - 2*M_PI*r*k);
1712 if(dPhi0<m_phi0Cut&&dKappa<m_kappaCut&&dTanl<m_tanlCut) sameTrack=true;
1713 //sameTrack = sameTrack&&(dDr<m_drCut);//FIXME
1714 //sameTrack = sameTrack&&(dPhi0<m_phi0Cut);//FIXME
1715 //sameTrack = sameTrack&&(dKappa<m_kappaCut);//FIXME
1716 //sameTrack = sameTrack&&(dDz<m_dzCut);//FIXME
1717 //sameTrack = sameTrack&&(dTanl<m_tanlCut);//FIXME
1718 if(sameTrack){
1719 //delete *trkIT2;
1720 //trackVector.erase(trkIT2);
1721 //trkIT2--;
1722 (trkIT2)->setFlag(1);
1723 if(m_debug) {
1724 cout<<"dDr, dPhi0, dKappa, dDz, dTanl "
1725 <<setw(10)<<dDr
1726 <<setw(10)<<dPhi0
1727 <<setw(10)<<dKappa
1728 <<setw(10)<<dDz
1729 <<setw(10)<<dTanl
1730 <<endl;
1731 cout<<"similar track "<<trkIT2->getTrkID()<<" is removed from hough track list"<<endl;
1732 }
1733
1734 // --- associate the hits of track 2 to track 1 if within pi/2 // to-do
1735 nDelete++;
1736 }
1737 }
1738 }
1739 if(m_debug) cout<<"checkTrack(): "<<nDelete<<" similar tracks are removed from hough track list"<<endl;
1740 //m_debug=false;
1741 return nDelete;
1742}
1743
1744/*
1745//fit tracks on xy plane, calculate S and Z of hits as well as add them to tracks, then calculate tanl and dz by Hough Transform and fit tracks again
1746int HoughFinder::completeTrack(vector<HoughTrack*>& trackVector)
1747{
1748 int nTrk(0);
1749 for(TrackVector_Iterator trackIT = trackVector.begin();trackIT!=trackVector.end();trackIT++){
1750 if((*trackIT)->getFlag() == 1)continue;
1751
1752 //cout<<"begin circlr fit "<<endl;
1753 TrkErrCode trkErrCode = (*trackIT)->fitCircle(m_mdcDetector,m_trkContextEv,m_bunchT0);
1754 //cout<<"finish circlr fit "<<endl;
1755 if(trkErrCode.failure()){
1756 (*trackIT)->setFlag(1);
1757 continue;
1758 }
1759 //cout<<"finish circlr fit "<<endl;
1760
1761 vector<HoughHit> stereoHot;
1762 for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end();hitIt++){
1763 HoughHit hit(*hitIt);
1764 if((*trackIT)->calculateZ_S(hit)>0)stereoHot.push_back(hit);
1765 }
1766 if(stereoHot.size()<2){//FIXME
1767 (*trackIT)->setFlag(1);
1768 continue;
1769 }
1770 //cout<<"finish stereo hot finding "<<endl;
1771
1772 double tanl = 0;
1773 double dz = 0;
1774 //double dTanl = 0.93/sqrt(1-0.93*0.93);
1775 double dTanl = 3.0;
1776 double dDz = 50.0;
1777 int mapType = 1;
1778 int step = 0;
1779 int charge = 0;
1780 vector<HoughTrack*> trkVector;
1781 trkVector.push_back(*trackIT);
1782 //cout<<"before trk finding "<<endl;
1783 findTrack(stereoHot,trkVector,tanl,dz,dTanl,dDz,charge,mapType,step);
1784 //cout<<"end trk finding "<<endl;
1785 TrackVector_Iterator trkIT = trkVector.begin();
1786 HoughTrack* trk = (*trkIT);
1787 int nHot = trk->findHot(stereoHot,mapType);
1788 if(nHot<2){//FIXME
1789 (*trackIT)->setFlag(1);
1790 continue;
1791 }
1792 trkErrCode = trk->fitHelix(m_mdcDetector,m_trkContextEv,m_bunchT0,m_mdcHitCol);
1793 if(trkErrCode.failure()){
1794 (*trackIT)->setFlag(1);
1795 continue;
1796 }
1797 nTrk++;
1798 }
1799 return nTrk;
1800}
1801//*/
1802StatusCode HoughFinder::registerTrack(RecMdcTrackCol*& recMdcTrackCol ,RecMdcHitCol*& recMdcHitCol)
1803{
1804 MsgStream log(msgSvc(), name());
1805 StatusCode sc;
1806 IDataProviderSvc* eventSvc = NULL;
1807 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
1808 if (eventSvc) {
1809 log << MSG::INFO << "makeTds:event Svc has been found" << endreq;
1810 } else {
1811 log << MSG::FATAL << "makeTds:Could not find eventSvc" << endreq;
1812 return StatusCode::FAILURE;
1813 //return StatusCode::SUCCESS;
1814 }
1815 IDataManagerSvc *dataManSvc = dynamic_cast<IDataManagerSvc*>(eventSvc);;
1816
1817
1818 // --- register RecMdcTrack
1819 recMdcTrackCol = new RecMdcTrackCol;
1820 DataObject *aRecMdcTrackCol;
1821 eventSvc->findObject("/Event/Recon/RecMdcTrackCol",aRecMdcTrackCol);
1822 if(aRecMdcTrackCol != NULL) {
1823 dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
1824 eventSvc->unregisterObject("/Event/Recon/RecMdcTrackCol");
1825 }
1826
1827 sc = eventSvc->registerObject("/Event/Recon/RecMdcTrackCol", recMdcTrackCol);
1828 if(sc.isFailure()) {
1829 log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
1830 return StatusCode::FAILURE;
1831 }
1832 log << MSG::INFO << "RecMdcTrackCol registered successfully!" <<endreq;
1833
1834
1835
1836 recMdcHitCol = new RecMdcHitCol;
1837 DataObject *aRecMdcHitCol;
1838 eventSvc->findObject("/Event/Recon/RecMdcHitCol",aRecMdcHitCol);
1839 if(aRecMdcHitCol != NULL) {
1840 dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
1841 eventSvc->unregisterObject("/Event/Recon/RecMdcHitCol");
1842 }
1843
1844 sc = eventSvc->registerObject("/Event/Recon/RecMdcHitCol", recMdcHitCol);
1845 if(sc.isFailure()) {
1846 log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
1847 return StatusCode::FAILURE;
1848 }
1849 log << MSG::INFO << "RecMdcHitCol registered successfully!" <<endreq;
1850
1851
1852 return StatusCode::SUCCESS;
1853}//end of stroeTracks
1854
1855int HoughFinder::storeTrack(vector<HoughTrack>& trackVector, RecMdcTrackCol*& recMdcTrackCol ,RecMdcHitCol*& recMdcHitCol)
1856{
1857 MdcTrackParams mdcTrackParams;
1858 //mdcTrackParams.lPrint=8;
1859 mdcTrackParams.lRemoveInActive=1;
1860 mdcTrackParams.maxGapLength=99;
1861 mdcTrackParams.nHitDeleted=99;
1862 MdcTrackList mdcTrackList(mdcTrackParams);
1863 mdcTrackList.setNoInner(true);
1864 vector<MdcTrack*> mdcTrackVector;
1865 for(vector<HoughTrack>::iterator trkIT = trackVector.begin();trkIT!=trackVector.end();trkIT++){
1866 if((trkIT)->getFlag() == 1)continue;
1867 if((trkIT)->getFlag() == -1)continue;
1868 if((trkIT)->getFlag() == -2)continue;
1869 //(trkIT)->sortHot();
1870 MdcTrack* mdcTrack = new MdcTrack((trkIT)->getTrkRecoTrk());
1871 mdcTrackList.append(mdcTrack);
1872 mdcTrackVector.push_back(mdcTrack);
1873 }
1874 int nDeleted(0);
1875 mdcTrackList.setNoInner(true);
1876 if(m_checkHits)nDeleted = mdcTrackList.arbitrateHits();
1877 if(nDeleted>0)cout<<"nDeleted "<<nDeleted<<endl;
1878
1879 int trkID=0;
1880 for(int l=0;l<mdcTrackList.length();l++){
1881 mdcTrackList[l]->storeTrack(trkID,recMdcTrackCol,recMdcHitCol,4);
1882 trkID++;
1883 }
1884
1885 vector<MdcTrack*>::iterator iter = mdcTrackVector.begin();
1886 for(;iter!=mdcTrackVector.end();iter++){
1887 if(m_clearTrack)delete *iter;
1888 }
1889 return trkID;
1890}
1891
1892bool sortCluster(const RecCgemCluster* clusterA , const RecCgemCluster* clusterB){
1893 return clusterA->getlayerid()<clusterB->getlayerid();
1894}
1895
1896int HoughFinder::storeRecTracks(RecMdcTrackCol* trackList, RecMdcHitCol* hitList, vector<HoughTrack>& trackVector)
1897{
1898 bool debug=false; //debug=true;
1899 int trackId(0);
1900 int tkStat =4;
1901 if(debug) cout<<"N Rectrack: "<<trackVector.size()<<endl;
1902 sort(trackVector.begin(),trackVector.end(),largerPt);
1903
1904 for(vector<HoughTrack>::iterator trkIT = trackVector.begin(); trkIT!=trackVector.end(); trkIT++)
1905 {
1906
1907 if(debug) cout<<"trk flag: "<<(trkIT)->getFlag()<<endl;
1908
1909 //if((trkIT)->getFlag() == 1)continue;
1910 //if((trkIT)->getFlag() == -1)continue;
1911 //if((trkIT)->getFlag() == -2)continue;
1912 if((trkIT)->getFlag()!=0) continue;
1913
1914 if(debug)
1915 {
1916 cout<<"the following track will be saved: "<<endl;
1917 trkIT->print();
1918 }
1919
1920 //new a Rec. Track MdcTrack
1921 RecMdcTrack* recMdcTrack = new RecMdcTrack();
1922
1923 recMdcTrack->setTrackId(trackId);
1924
1925 double helixPar[5];
1926 helixPar[0]=trkIT->getDr();
1927 helixPar[1]=trkIT->getPhi0();
1928 helixPar[2]=trkIT->getKappa();
1929 helixPar[3]=trkIT->getDz();
1930 helixPar[4]=trkIT->getTanl();//helixPar[4]+=0.1;
1931 recMdcTrack->setHelix(helixPar);
1932
1933 int q = trkIT->getCharge();
1934 double pxy = trkIT->pt();
1935 double px = trkIT->momentum(0).x();
1936 double py = trkIT->momentum(0).y();
1937 double pz = trkIT->momentum(0).z();
1938 double p = trkIT->momentum(0).mag();
1939 double theta = trkIT->direction(0).theta();
1940 double phi = trkIT->direction(0).phi();
1941 HepPoint3D poca = trkIT->x(0);
1942 HepPoint3D pivot = trkIT->pivot();
1943 double r = poca.perp();
1944 HepSymMatrix Ea = trkIT->Ea();
1945 //cout<<"Ea="<<Ea<<endl;
1946 double errorMat[15];
1947 int k = 0;
1948 for (int ie = 0 ; ie < 5 ; ie ++){
1949 for (int je = ie ; je < 5 ; je ++){
1950 errorMat[k] = Ea[ie][je];
1951 k++;
1952 }
1953 }
1954 // ---> a Helix for DotsHelixFitter
1955 HepVector a(5,0);
1956 a(1) =trkIT->getDr();
1957 a(2) =trkIT->getPhi0();
1958 a(3) =trkIT->getKappa();
1959 a(4) =trkIT->getDz();
1960 a(5) =trkIT->getTanl();
1961 KalmanFit::Helix aHelix(pivot,a,Ea); // <-- a Helix for DotsHelixFitter ---
1962 //myDotsHelixFitter.setInitialHelix(aHelix);// set a Helix for DotsHelixFitter
1963 double chisq = trkIT->getChi2();
1964 recMdcTrack->setCharge(q);
1965 recMdcTrack->setPxy(pxy);
1966 recMdcTrack->setPx(px);
1967 recMdcTrack->setPy(py);
1968 recMdcTrack->setPz(pz);
1969 recMdcTrack->setP(p);
1970 recMdcTrack->setTheta(theta);
1971 recMdcTrack->setPhi(phi);
1972 recMdcTrack->setPoca(poca);
1973 recMdcTrack->setX(poca.x());//poca
1974 recMdcTrack->setY(poca.y());
1975 recMdcTrack->setZ(poca.z());
1976 recMdcTrack->setR(sqrt(poca.x()*poca.x() + poca.y()*poca.y()));
1977 HepPoint3D apivot(0.,0.,0.);
1978 recMdcTrack->setPivot(apivot);
1979 recMdcTrack->setVX0(0.);//pivot
1980 recMdcTrack->setVY0(0.);
1981 recMdcTrack->setVZ0(0.);
1982 recMdcTrack->setError(errorMat);
1983 recMdcTrack->setError(Ea);
1984 recMdcTrack->setChi2(chisq);
1985 recMdcTrack->setStat(tkStat);
1986
1987 //-----------hit list-------------
1988 int hitId = 0;
1989 HitRefVec hitRefVec;
1990 ClusterRefVec clusterRefVec;
1991 //vector<int> vecHits;
1992 map<int,int> clusterFitStat;
1993 //int maxLayer = 0;
1994 int maxLayerId = -1;
1995 int minLayerId = 43;
1996 double fiTerm = 999.;
1997 HoughHit* fiTermHot = NULL;
1998
1999 //vector<HoughHit*> hotList = trkIT->getVecHitPnt();
2000 vector<HoughHit*> hotList = trkIT->getHotList(2);
2001 for(vector<HoughHit*>::iterator hotIt = hotList.begin(); hotIt != hotList.end();hotIt++)
2002 {
2003 //cout<<"handle hit "<<hitId<<endl;
2004 HoughHit* hot(*hotIt);
2005 //int layer = (**hotIt).getLayer();
2006 //int wire = (**hotIt).getWire();
2007 //cout<<"("<<setw(2)<<layer<<","<<setw(3)<<wire<<") ";
2008
2009 /*
2010 if(hot->getHitType()==HoughHit::MDCHIT)
2011 {
2012 //if(hot->getUse()!=1) continue;
2013 //cout<<"create a new RecMdcHit"<<endl;
2014
2015 // --- make a RecMdcHit from digi, aHelix using DotsHelixFitter
2016 RecMdcHit* recMdcHit = new RecMdcHit(myDotsHelixFitter.makeRecMdcHit(hot->getDigi(), aHelix));
2017 recMdcHit->setId(hitId);
2018 recMdcHit->setTrkId(trackId);
2019
2020 //
2021 // --> make a RecMdcHit from a HoughHit (parameters from Hough)
2022 RecMdcHit* recMdcHit = new RecMdcHit;
2023 recMdcHit->setId(hitId);
2024 recMdcHit->setTrkId(trackId);
2025 HepPoint3D hitPoint = hot->getHitPosition();
2026 double dPhi = trkIT->dPhi(hitPoint);
2027 HepPoint3D position = trkIT->x(dPhi);
2028 Hep3Vector direction = trkIT->direction(dPhi);
2029 //cout<<"get 1"<<endl;
2030
2031 double ddl = -1*hot->getDriftDist();
2032 double ddr = hot->getDriftDist();
2033 //double erddl = hot->hitSigma();
2034 //double erddr = hot->hitSigma();
2035 double pChisq = hot->getResidual()*hot->getResidual();//FIXME
2036 int lr = 0;
2037 int stat = hot->getUse();
2038 stat = 1; // llwang
2039 //cout<<"hit stat = "<<stat<<endl;
2040 //cout<<"get 2"<<endl;
2041 Identifier mdcId = hot->getDigi()->identify();
2042 //cout<<"get 2.1"<<endl;
2043 double tdc = hot->getDigi()->getTimeChannel();
2044 //cout<<"get 2.2"<<endl;
2045 double adc = hot->getDigi()->getChargeChannel();
2046 //cout<<"get 2.3"<<endl;
2047 double driftT = hot->driftTime();
2048 //cout<<"get 2.4"<<endl;
2049 double doca = hot->getDriftDist() +hot->getResidual();
2050 //cout<<"get 2.5"<<endl;
2051 double entra = direction.phi()-position.phi();//FIXME
2052 //cout<<"get 2.6"<<endl;
2053 while(entra<-0.5*M_PI)entra+= M_PI;
2054 while(entra> 0.5*M_PI)entra -= M_PI;
2055 //cout<<"get 3"<<endl;
2056 double zhit = hitPoint.z();
2057 double fltLen = trkIT->flightLength(hitPoint);
2058 //cout<<"after get "<<hitId<<endl;
2059
2060 recMdcHit->setDriftDistLeft(ddl);
2061 recMdcHit->setDriftDistRight(ddr);
2062 //recMdcHit->setErrDriftDistLeft(erddl);
2063 //recMdcHit->setErrDriftDistRight(erddr);
2064 recMdcHit->setChisqAdd(pChisq);
2065 recMdcHit->setFlagLR(lr);
2066 recMdcHit->setStat(stat);
2067 recMdcHit->setMdcId(mdcId);
2068 recMdcHit->setTdc(tdc);
2069 recMdcHit->setAdc(adc);
2070 recMdcHit->setDriftT(driftT);
2071 recMdcHit->setDoca(doca);
2072 recMdcHit->setEntra(entra);
2073 recMdcHit->setZhit(zhit);
2074 recMdcHit->setFltLen(fltLen);
2075 //cout<<"after set "<<hitId<<endl;
2076 // <--- make a RecMdcHit from a HoughHit
2077
2078 hitList->push_back(recMdcHit);
2079 SmartRef<RecMdcHit> refHit(recMdcHit);
2080 hitRefVec.push_back(refHit);
2081
2082 if(hot->getLayer()>maxLayerId){
2083 maxLayerId = hot->getLayer();
2084 fiTermHot = hot;
2085 }
2086 if(hot->getLayer()<minLayerId){
2087 minLayerId = hot->getLayer();
2088 }
2089 //if(hot->getUse()==1){
2090 recMdcHit->setStat(1);
2091 //}else{
2092 // recMdcHit->setStat(0);
2093 //}
2094 //cout<<"finish a DC hit "<<hitId<<endl;
2095 hitId++;
2096 }*/
2097 if(hot->getHitType()==HoughHit::CGEMHIT)
2098 {
2099 //hot->print();
2100 //cout<<endl;
2101 const RecCgemCluster* recCgemCluster = hot->getCgemCluster();
2102 int clusterid = hot->getCgemCluster()->getclusterid();
2103 //int stat = hot->getUse();
2104 int stat = 1;
2105 clusterFitStat[clusterid] = stat;
2106 clusterRefVec.push_back(recCgemCluster);
2107 if(debug) cout<<" cluster "<<clusterid<<" kept! "<<endl;
2108 if(hot->getLayer()>maxLayerId)
2109 {
2110 maxLayerId = hot->getLayer();
2111 fiTermHot = hot;
2112 }
2113 if(hot->getLayer()<minLayerId)
2114 {
2115 minLayerId = hot->getLayer();
2116 }
2117 }
2118 }// loop hotList (vector<HoughHit*>) of a hough track
2119
2120 // --- loop RecMdcHitVec in Hough track
2121 double fltLen = -1;
2122 vector<RecMdcHit>* RecMdcHitVecPnt = trkIT->getRecMdcHitVec();
2123 int nMdcHits=RecMdcHitVecPnt->size();
2124 int nMdcHitsKept=0;
2125 vector<RecMdcHit>::iterator iter_recMdcHit = RecMdcHitVecPnt->begin();
2126 for(; iter_recMdcHit!=RecMdcHitVecPnt->end(); iter_recMdcHit++)
2127 {
2128 RecMdcHit* recMdcHit = new RecMdcHit(*iter_recMdcHit);
2129 if(recMdcHit->getChisqAdd()>m_chi2CutHits) // skip hit with chi2>m_chi2CutHits
2130 {
2131 //cout<<"skip hit with chi2="<<recMdcHit->getChisqAdd()<<endl;
2132 delete recMdcHit;
2133 continue;
2134 }
2135 recMdcHit->setId(hitId);
2136 recMdcHit->setTrkId(trackId);
2137 recMdcHit->setStat(1);
2138 hitList->push_back(recMdcHit);
2139 SmartRef<RecMdcHit> refHit(recMdcHit);
2140 hitRefVec.push_back(refHit);
2141 nMdcHitsKept++;
2142
2143 Identifier mdcid = recMdcHit->getMdcId();
2144 int layer = MdcID::layer(mdcid);
2145 int wire = MdcID::wire(mdcid);
2146 if(debug) cout<<" hit at layer "<<layer<<" wire "<<wire<<" kept! "<<endl;
2147 if(layer>maxLayerId)
2148 {
2149 maxLayerId = layer;
2150 fltLen=recMdcHit->getFltLen();
2151 //fiTermHot = hot;
2152 }
2153 if(layer<minLayerId)
2154 {
2155 minLayerId = layer;
2156 }
2157 hitId++;
2158 }
2159 if(debug) cout<<"track "<<trackId<<", "<<nMdcHitsKept<<"/"<<nMdcHits<<" hits kept"<<endl;
2160
2161 //fi terminal angle
2162 if (fiTermHot!=NULL){
2163 //fiTerm=(1./sqrt(1.+tanDip*tanDip))*fiTermHot->fltLen()*helix.omega();
2164 HepPoint3D point = fiTermHot->getHitPosition();
2165 fiTerm = trkIT->flightArc(point)/trkIT->radius();
2166 }
2167 if(fltLen>0) fiTerm=-fltLen*sin(theta)/trkIT->radius();
2168 recMdcTrack->setFiTerm(fiTerm);
2169
2170 // --- setN* functions called in setVec* functions
2171 recMdcTrack->setVecHits(hitRefVec);
2172 std::sort(clusterRefVec.begin(),clusterRefVec.end(),sortCluster);
2173 recMdcTrack->setVecClusters(clusterRefVec,clusterFitStat);
2174 trackList->push_back(recMdcTrack);
2175 trackId++;
2176 }// loop HoughTrack list
2177 return trackList->size();
2178
2179}
2180
2181int HoughFinder::storeTracks(RecMdcTrackCol* trackList, RecMdcHitCol* hitList, vector<HoughTrack>& trackVector){
2182 bool debug=false; //debug=true;
2183 int trackId(0);
2184 int tkStat =4;
2185 if(debug) cout<<"N Rectrack: "<<trackVector.size()<<endl;
2186 sort(trackVector.begin(),trackVector.end(),largerPt);
2187
2188 for(vector<HoughTrack>::iterator trkIT = trackVector.begin();trkIT!=trackVector.end();trkIT++){
2189 if(debug) cout<<"trk flag: "<<(trkIT)->getFlag()<<endl;
2190 if((trkIT)->getFlag() == 1)continue;
2191 if((trkIT)->getFlag() == -1)continue;
2192 if((trkIT)->getFlag() == -2)continue;
2193 if(trkIT->getTrkRecoTrk()==NULL)
2194 {
2195 if(debug) cout<<"trk getTrkRecoTrk NULL!"<<endl;
2196 continue;
2197 }
2198 //get the results of the fit to this track
2199 const TrkFit* fitresult = trkIT->getTrkRecoTrk()->fitResult();
2200 //check the fit worked
2201 if (fitresult == 0)
2202 {
2203 if(debug) cout<<"no fit result!"<<endl;
2204 continue;
2205 }
2206 if(debug)
2207 {
2208 cout<<"the following track will be saved: "<<endl;
2209 trkIT->print();
2210 }
2211 if(printFlag)trkIT->print();
2212 if(printFlag)trkIT->printHot();
2213 //cout<<endl;
2214
2215 //new a Rec. Track MdcTrack
2216 RecMdcTrack* recMdcTrack = new RecMdcTrack();
2217 const TrkHitList* aList = trkIT->getTrkRecoTrk()->hits();
2218 //some track info
2219 const BField& theField = trkIT->getTrkRecoTrk()->bField();
2220 double Bz = theField.bFieldZ();
2221 //Fit info
2222 int hitId = 0;
2223 int nHits=0;
2224 int nSt=0;
2225 //int nAct=0; int nSeg=0;
2226 //int maxlay = 0; int minlay = 43; int seg[11];/* 11 slayer */ int segme;
2227 //****************************
2228 //* active flag open this
2229 //****************************
2230 //* if (allHit>0){
2231 //* nHits = aList->nHit();//add inActive hit also
2232 //* }else{
2233 //* nHits = fitresult->nMdc();// = nActive()
2234 //* }
2235 //****************************
2236 //for 5.1.0 use all hits
2237 nHits = aList->nHit();
2238 //****************************
2239 // use active only
2240 // nHits = fitresult->nMdc();// = nActive()
2241 //****************************
2242
2243 int q = fitresult->charge();
2244 double chisq = fitresult->chisq();
2245 int nDof = fitresult->nDof();
2246
2247 //if(nHits-nDof!=5) cout<<"HoughFinder::storeTracks: nHits-nDof!=5, nHits="<<nHits<<" , nDof="<<nDof<<", chi2="<<chisq<<" , evt="<<m_event<<endl;
2248
2249 //track parameters at closest approach to beamline
2250 double fltLenPoca = 0.0;
2251 TrkExchangePar helix = fitresult->helix(fltLenPoca);
2252 //std::cout<< __FILE__ << " phi0 " << helix.phi0() << " omega " <<helix.omega()<<std::endl;
2253 double phi0 = helix.phi0();
2254 double tanDip = helix.tanDip();
2255 //double z0 = helix.z0();
2256 double d0 = helix.d0();
2257 //momenta and position
2258 //CLHEP::Hep3Vector mom = fitresult->momentum(fltLenPoca);
2259 HepPoint3D poca = fitresult->position(fltLenPoca);
2260
2261 double helixPar[5];
2262 //dro =-d0
2263 helixPar[0] = -d0;
2264 //phi0 = phi0 - pi/2 [0,2pi)
2265 double tphi0 = phi0 - Constants::pi/2.;
2266 if (tphi0 < 0. ) tphi0 += Constants::pi * 2.;
2267 helixPar[1] = tphi0;
2268 //kappa = q/pxy
2269 double pxy = fitresult->pt();
2270 if (pxy == 0.) helixPar[2] = 9999.;
2271 else helixPar[2] = q/fabs(pxy);
2272 if(pxy>9999.) helixPar[2] = 0.00001;
2273 //dz = z0
2274 helixPar[3] = helix.z0();
2275 //tanl
2276 helixPar[4] = tanDip;
2277 //error
2278 //V(Y)=S * V(X) * ST , mS = S , mVy = V(Y) , helix.covariance() = V(X)
2279 HepSymMatrix mS(helix.params().num_row(),0);
2280 mS[0][0]=-1.;//dr0=-d0
2281 mS[1][1]=1.;
2282 mS[2][2]=-333.567/Bz;//pxy -> cpa
2283 mS[3][3]=1.;
2284 mS[4][4]=1.;
2285 HepSymMatrix mVy = helix.covariance().similarity(mS);
2286 double errorMat[15];
2287 int k = 0;
2288 for (int ie = 0 ; ie < 5 ; ie ++){
2289 for (int je = ie ; je < 5 ; je ++) {
2290 errorMat[k] = mVy[ie][je];
2291 k++;
2292 }
2293 }
2294 double p,px,py,pz;
2295 px = pxy * (-sin(helixPar[1]));
2296 py = pxy * cos(helixPar[1]);
2297 pz = pxy * helixPar[4];
2298 p = sqrt(pxy*pxy + pz*pz);
2299 //theta, phi
2300 double theta = acos(pz/p);
2301 double phi = atan2(py,px);
2302 recMdcTrack->setTrackId(trackId);
2303 recMdcTrack->setHelix(helixPar);
2304 recMdcTrack->setCharge(q);
2305 recMdcTrack->setPxy(pxy);
2306 recMdcTrack->setPx(px);
2307 recMdcTrack->setPy(py);
2308 recMdcTrack->setPz(pz);
2309 recMdcTrack->setP(p);
2310 recMdcTrack->setTheta(theta);
2311 recMdcTrack->setPhi(phi);
2312 recMdcTrack->setPoca(poca);
2313 recMdcTrack->setX(poca.x());//poca
2314 recMdcTrack->setY(poca.y());
2315 recMdcTrack->setZ(poca.z());
2316 recMdcTrack->setR(sqrt(poca.x()*poca.x() + poca.y()*poca.y()));
2317 HepPoint3D pivot(0.,0.,0.);
2318 recMdcTrack->setPivot(pivot);
2319 recMdcTrack->setVX0(0.);//pivot
2320 recMdcTrack->setVY0(0.);
2321 recMdcTrack->setVZ0(0.);
2322 recMdcTrack->setError(errorMat);
2323 recMdcTrack->setError(mVy);
2324 recMdcTrack->setChi2(chisq);
2325 recMdcTrack->setStat(tkStat);
2326 //recMdcTrack->setNhits(nHits);
2327 //-----------hit list-------------
2328 HitRefVec hitRefVec;
2329 ClusterRefVec clusterRefVec;
2330 vector<int> vecHits;
2331 map<int,int> clusterFitStat;
2332 //for fiterm
2333 int maxLayer = 0;
2334 //cout<<__FILE__ <<" "<<__LINE__ <<endl;
2335 for(TrkHitList::hot_iterator hot = aList->begin();hot!=aList->end();hot++){
2336 if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack) /*|| typeid(*hot)==typeid(TrkHitOnTrk)*/ ){
2337 const MdcRecoHitOnTrack* recoHot;
2338 recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
2339 int layerId = recoHot->mdcHit()->layernumber();
2340 int wireId = recoHot->mdcHit()->wirenumber();
2341 if(maxLayer < layerId)maxLayer = layerId;
2342 }
2343 }
2344 //cout<<"maxLayer:"<<maxLayer<<endl;
2345 int maxLayerId = -1;
2346 int minLayerId = 43;
2347 double fiTerm = 999.;
2348 const MdcRecoHitOnTrack* fiTermHot = NULL;
2349 TrkHitList::hot_iterator hot = aList->begin();
2350 for (;hot!=aList->end();hot++){
2351 //cout<<__FILE__ <<" "<<__LINE__ <<" "<< typeid(*hot).name()<<endl;
2352 if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack) /*|| typeid(*hot)==typeid(TrkHitOnTrk)*/ ){
2353 const MdcRecoHitOnTrack* recoHot;
2354 recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
2355 if(!(recoHot->isActive())) continue;// added 2020-05-26 wangll
2356 int layerId = recoHot->mdcHit()->layernumber();
2357 if(layerId>(maxLayer - m_removeNOuterHits))continue;
2358 //if (!recoHot->mdcHit()) return;
2359 RecMdcHit* recMdcHit = new RecMdcHit;
2360 //id
2361 recMdcHit->setId(hitId);
2362 //---------BESIII----------
2363 //phiWire - phiHit <0 flagLR=0 left driftleft<0 ;
2364 //phiWire - phiHit >0 flagLR=1 right driftright>0;
2365 //flag = 2 ambig;
2366 //-----Babar wireAmbig----
2367 //-1 -> right, 0 -> don't know, +1 -> left
2368 //+1 phiWire-phiHit<0 Left,
2369 //-1 phiWire-phiHit>0 Right,
2370 //0 don't know
2371 //ambig() w.r.t track direction
2372 //wireAmbig() w.r.t. wire location
2373 double hotWireAmbig = recoHot->wireAmbig();
2374 double driftDist = fabs(recoHot->drift());
2375 double sigma = recoHot->hitRms();
2376 double doca = fabs(recoHot->dcaToWire());
2377 int flagLR=2;
2378 if ( hotWireAmbig == 1){
2379 flagLR = 0; //left driftDist <0
2380 doca *= -1.;//2012-07-19
2381 }else if( hotWireAmbig == -1){
2382 flagLR = 1; //right driftDist >0
2383 }else if( hotWireAmbig == 0){
2384 flagLR = 2; //don't know
2385 }
2386 recMdcHit->setFlagLR(flagLR);
2387 recMdcHit->setDriftDistLeft((-1)*driftDist);
2388 recMdcHit->setErrDriftDistLeft(sigma);
2389 recMdcHit->setDriftDistRight(driftDist);
2390 recMdcHit->setErrDriftDistRight(sigma);
2391 //Mdc Id
2392 Identifier mdcId = recoHot->mdcHit()->digi()->identify();
2393 recMdcHit->setMdcId(mdcId);
2394 //corrected ADC
2395
2396 //contribution to chisquare
2397 //fill chisq
2398 double res=999.,rese=999.;
2399 if (recoHot->resid(res,rese,false)){
2400 }else{}
2401 double deltaChi=0;
2402 recoHot->getFitStuff(deltaChi);//yzhang open 2010-09-20
2403 recMdcHit->setChisqAdd( deltaChi * deltaChi );
2404 //set tracks containing this hit
2405 recMdcHit->setTrkId(trackId);
2406 //doca
2407 recMdcHit->setDoca(doca);//doca sign left<0
2408 //entrance angle
2409
2410 recMdcHit->setEntra(recoHot->entranceAngle());
2411 //z of hit
2412 HepPoint3D pos; Hep3Vector dir;
2413 double fltLen = recoHot->fltLen();
2414 recMdcHit->setFltLen(fltLen);
2415 recoHot->getParentRep()->traj().getInfo(fltLen,pos,dir);
2416 recMdcHit->setZhit(pos.z());
2417 double tof = recoHot->getParentRep()->arrivalTime(fltLen);
2418 recMdcHit->setTdc(recoHot->mdcHit()->tdcIndex());
2419 recMdcHit->setAdc(recoHot->mdcHit()->adcIndex());
2420 //drift time
2421 recMdcHit->setDriftT(recoHot->mdcHit()->driftTime(tof,pos.z()));
2422 //for fiterm
2423 int wireId = recoHot->mdcHit()->wirenumber();
2424 //std::cout<<" MdcTrack::store() ("<<layerId<<","<<wireId<<") fltLen "<<fltLen<<" wireAmbig "<<hotWireAmbig<<" y "<<pos.y()<<std::endl;
2425 // <<recoHot->entranceAngle()<<std::endl;
2426 if (layerId >= maxLayerId){
2427 maxLayerId = layerId;
2428 fiTermHot = recoHot;
2429 }
2430 if (layerId < minLayerId){
2431 minLayerId = layerId;
2432 }
2433 // status flag
2434 if (recoHot->isActive()) {
2435 recMdcHit->setStat(1);
2436 }else{
2437 recMdcHit->setStat(0);
2438 }
2439 // for 5.1.0 use all hits
2440 if (recoHot->layer()->view()) {
2441 ++nSt;
2442 }
2443 hitList->push_back(recMdcHit);
2444 SmartRef<RecMdcHit> refHit(recMdcHit);
2445 hitRefVec.push_back(refHit);
2446 vecHits.push_back(mdcId.get_value());
2447 ++hitId;
2448 }else if(typeid(*hot)==typeid(CgemHitOnTrack)){
2449 const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hot));
2450 const RecCgemCluster* recCgemCluster = recoHot->baseHit();
2451 int clusterid = recCgemCluster->getclusterid();
2452 int stat = recoHot->isActive();
2453 if(stat==0) continue;// added 2020-05-26 wangll
2454 //cout<<"stat:"<<stat<<endl;
2455 clusterFitStat[clusterid] = stat;
2456 //cout<<clusterid<<" "<<stat<<endl;
2457 //cout<<"clusterid:"<< recCgemCluster->getclusterid()<<" layer:"<<recCgemCluster->getlayerid()<<endl;
2458 clusterRefVec.push_back(recCgemCluster);
2459 }
2460 }
2461 std::sort(clusterRefVec.begin(),clusterRefVec.end(),sortCluster);
2462 //fi terminal angle
2463 if (fiTermHot!=NULL){
2464 fiTerm=(1./sqrt(1.+tanDip*tanDip))*fiTermHot->fltLen()*helix.omega();
2465 }
2466 recMdcTrack->setFiTerm(fiTerm);
2467 // number of stereo hits contained
2468 //recMdcTrack->setNhits(hitId);
2469 //recMdcTrack->setNster(nSt);
2470 //recMdcTrack->setNdof(nDof);
2471 //recMdcTrack->setFirstLayer(maxLayerId);
2472 //recMdcTrack->setLastLayer(minLayerId);
2473 //recMdcTrack->setFirstLayer(minLayerId);
2474 //recMdcTrack->setLastLayer(maxLayerId);
2475 //recMdcTrack->setClusterFitStat(clusterFitStat);
2476 //recMdcTrack->setVecClusters(clusterRefVec);
2477 //recMdcTrack->setNcluster(clusterRefVec.size());
2478
2479 // --- setN* functions called in setVec* functions
2480 recMdcTrack->setVecHits(hitRefVec);
2481 recMdcTrack->setVecClusters(clusterRefVec,clusterFitStat);
2482 trackList->push_back(recMdcTrack);
2483 trackId++;
2484 }
2485 //cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
2486 /*
2487 int trackId(0);
2488 for(vector<HoughTrack>::iterator trkIT = trackVector.begin();trkIT!=trackVector.end();trkIT++){
2489//TrkErrCode trkErrCode = trkIT->fitCircle(m_mdcDetector,m_trkContextEv,m_bunchT0);
2490HoughTrack* houghTrack = &(*trkIT);
2491if(houghTrack->getFlag() == 1)continue;
2492if(houghTrack->getFlag() < 0)continue;
2493//houghTrack->print();
2494//houghTrack->printHot();
2495trackId=trackList->size();
2496//cout<<"start store trk "<<trackId;
2497RecMdcTrack* recMdcTrack = new RecMdcTrack();
2498//vector<HoughHit>* hotList = houghTrack->getHot();
2499//vector<HoughHit*> hotList = houghTrack->getVecHitPnt();
2500vector<HoughHit*> hotList = houghTrack->getHotList();
2501//cout<<" with "<<hotList.size()<<" hits";//<<endl;
2502double helixPar[5];
2503helixPar[0] = houghTrack->dr();
2504helixPar[1] = houghTrack->phi0();
2505helixPar[2] = houghTrack->kappa();
2506helixPar[3] = houghTrack->dz();
2507helixPar[4] = houghTrack->tanl();
2508int q = houghTrack->getCharge();
2509double pxy = houghTrack->pt();
2510double px = houghTrack->momentum(0).x();
2511double py = houghTrack->momentum(0).y();
2512double pz = houghTrack->momentum(0).z();
2513double p = houghTrack->momentum(0).mag();
2514double theta = houghTrack->direction(0).theta();
2515double phi = houghTrack->direction(0).phi();
2516HepPoint3D poca = houghTrack->x(0);
2517HepPoint3D pivot = houghTrack->pivot();
2518double r = poca.perp();
2519HepSymMatrix Ea = houghTrack->Ea();
2520//cout<<" phi, pxy, pz = "<<phi<<", "<<pxy<<", "<<pz<<endl;
2521double errorMat[15];
2522int k = 0;
2523for (int ie = 0 ; ie < 5 ; ie ++){
2524for (int je = ie ; je < 5 ; je ++){
2525errorMat[k] = Ea[ie][je];
2526k++;
2527}
2528}
2529double chisq = houghTrack->getChi2();
2530int tkStat = 4+houghTrack->getCircleFitStat();
2531//if(tkStat==4) cout<<" track "<<trackId<<" has no good circle fit"<<endl;
2532int nHits = hotList.size();
2533
2534recMdcTrack->setTrackId(trackId);
2535recMdcTrack->setHelix(helixPar);
2536recMdcTrack->setCharge(q);
2537recMdcTrack->setPxy(pxy);
2538recMdcTrack->setPx(px);
2539recMdcTrack->setPy(py);
2540recMdcTrack->setPz(pz);
2541recMdcTrack->setP(p);
2542recMdcTrack->setTheta(theta);
2543recMdcTrack->setPhi(phi);
2544recMdcTrack->setPoca(poca);
2545recMdcTrack->setX(poca.x());//poca
2546recMdcTrack->setY(poca.y());
2547recMdcTrack->setZ(poca.z());
2548recMdcTrack->setR(r);
2549recMdcTrack->setPivot(pivot);
2550recMdcTrack->setVX0(0.);//pivot
2551recMdcTrack->setVY0(0.);
2552recMdcTrack->setVZ0(0.);
2553recMdcTrack->setError(errorMat);
2554recMdcTrack->setError(Ea);
2555recMdcTrack->setChi2(chisq);
2556recMdcTrack->setStat(tkStat);
2557recMdcTrack->setNhits(nHits);
2558//-----------hit list-------------
2559int hitId = 0;
2560int nSt(0);
2561int maxLayerId = -1;
2562int minLayerId = 43;
2563double fiTerm = 999.;
2564HoughHit fiTermHot;
2565
2566HitRefVec hitRefVec;
2567ClusterRefVec clusterRefVec;
2568map<int,int> clusterFitStat;
2569for(vector<HoughHit*>::iterator hotIt = hotList.begin(); hotIt != hotList.end();hotIt++)
2570{
2571 //cout<<"handle hit "<<hitId<<endl;
2572 HoughHit* hot(*hotIt);
2573 //int layer = (**hotIt).getLayer();
2574 //int wire = (**hotIt).getWire();
2575 //cout<<"("<<setw(2)<<layer<<","<<setw(3)<<wire<<") ";
2576 if(hot->getHitType()==HoughHit::MDCHIT){
2577 //cout<<"create a new RecMdcHit"<<endl;
2578 RecMdcHit* recMdcHit = new RecMdcHit;
2579
2580 recMdcHit->setId(hitId);
2581 recMdcHit->setTrkId(trackId);
2582 HepPoint3D hitPoint = hot->getHitPosition();
2583 double dPhi = houghTrack->dPhi(hitPoint);
2584 HepPoint3D position = houghTrack->x(dPhi);
2585 Hep3Vector direction = houghTrack->direction(dPhi);
2586 //cout<<"get 1"<<endl;
2587
2588 double ddl = -1*hot->getDriftDist();
2589 double ddr = hot->getDriftDist();
2590 //double erddl = hot->hitSigma();
2591 //double erddr = hot->hitSigma();
2592 double pChisq = hot->getResidual()*hot->getResidual();//FIXME
2593 int lr = 0;
2594 int stat = hot->getUse();
2595 stat = 1; // llwang
2596 //cout<<"hit stat = "<<stat<<endl;
2597 //cout<<"get 2"<<endl;
2598 Identifier mdcId = hot->getDigi()->identify();
2599 //cout<<"get 2.1"<<endl;
2600 double tdc = hot->getDigi()->getTimeChannel();
2601 //cout<<"get 2.2"<<endl;
2602 double adc = hot->getDigi()->getChargeChannel();
2603 //cout<<"get 2.3"<<endl;
2604 double driftT = hot->driftTime();
2605 //cout<<"get 2.4"<<endl;
2606 double doca = hot->getDriftDist() +hot->getResidual();
2607 //cout<<"get 2.5"<<endl;
2608 double entra = direction.phi()-position.phi();//FIXME
2609 //cout<<"get 2.6"<<endl;
2610 while(entra<-0.5*M_PI)entra+= M_PI;
2611 while(entra> 0.5*M_PI)entra -= M_PI;
2612 //cout<<"get 3"<<endl;
2613 double zhit = hitPoint.z();
2614 double fltLen = houghTrack->flightLength(hitPoint);
2615 //cout<<"after get "<<hitId<<endl;
2616
2617 recMdcHit->setDriftDistLeft(ddl);
2618 recMdcHit->setDriftDistRight(ddr);
2619 //recMdcHit->setErrDriftDistLeft(erddl);
2620 //recMdcHit->setErrDriftDistRight(erddr);
2621 recMdcHit->setChisqAdd(pChisq);
2622 recMdcHit->setFlagLR(lr);
2623 recMdcHit->setStat(stat);
2624 recMdcHit->setMdcId(mdcId);
2625 recMdcHit->setTdc(tdc);
2626 recMdcHit->setAdc(adc);
2627 recMdcHit->setDriftT(driftT);
2628 recMdcHit->setDoca(doca);
2629 recMdcHit->setEntra(entra);
2630 recMdcHit->setZhit(zhit);
2631 recMdcHit->setFltLen(fltLen);
2632 //cout<<"after set "<<hitId<<endl;
2633
2634 hitList->push_back(recMdcHit);
2635 SmartRef<RecMdcHit> refHit(recMdcHit);
2636 hitRefVec.push_back(refHit);
2637
2638 if(hot->getLayer()>maxLayerId){
2639 maxLayerId = hot->getLayer();
2640 fiTermHot = *hot;
2641 }
2642 if(hot->getLayer()<minLayerId){
2643 minLayerId = hot->getLayer();
2644 }
2645 if(hot->getUse()==1){
2646 recMdcHit->setStat(1);
2647 }else{
2648 recMdcHit->setStat(0);
2649 }
2650 if(hot->getFlag()!=0)++nSt;
2651 //cout<<"finish a DC hit "<<hitId<<endl;
2652 }else if(hot->getHitType()==HoughHit::CGEMHIT){
2653 //hot->print();
2654 //cout<<endl;
2655 const RecCgemCluster* recCgemCluster = hot->getCgemCluster();
2656 int clusterid = hot->getCgemCluster()->getclusterid();
2657 //int stat = hot->getUse();
2658 int stat = 1;
2659 clusterFitStat[clusterid] = stat;
2660 clusterRefVec.push_back(recCgemCluster);
2661 }
2662 hitId++;
2663}
2664
2665//HepPoint3D point = fiTermHot.getHitPosition();
2666//fiTerm = houghTrack->flightArc(point)/houghTrack->radius();
2667//recMdcTrack->setFiTerm(fiTerm);
2668recMdcTrack->setNdof(nHits-5);
2669recMdcTrack->setNster(nSt);
2670recMdcTrack->setFirstLayer(minLayerId);
2671recMdcTrack->setLastLayer(maxLayerId);
2672//cout<<"push back a RecMdcTrack with "<<recMdcTrack->getNhits()<<" hits, "<<recMdcTrack->getNcluster()<<" clusters"<<endl;
2673recMdcTrack->setVecHits(hitRefVec);
2674//cout<<"before setVecClusters"<<endl;
2675//recMdcTrack->setVecClusters(clusterRefVec,clusterFitStat);
2676//recMdcTrack->setVecClusters(clusterRefVec);
2677recMdcTrack->setVecClusters(clusterRefVec,clusterFitStat);
2678//cout<<"clusterRefVec.size = "<<clusterRefVec.size()<<endl;
2679//cout<<"end setVecClusters"<<endl;
2680recMdcTrack->setNcluster(clusterRefVec.size());
2681//cout<<"push back a RecMdcTrack with "<<recMdcTrack->getNhits()<<" hits, "<<recMdcTrack->getNcluster()<<" clusters"<<endl;
2682trackList->push_back(recMdcTrack);
2683//trackId++;
2684//cout<<helixPar[2]<<endl;
2685}
2686return trackList->size();
2687//*/
2688return trackList->size();
2689}
2690
2691void HoughFinder::printHitList()
2692{
2693 //cout<<"nHit ==========> "<<m_houghHitList.size()<<endl;
2694 cout<<"========== truth =========="<<endl;
2695 for(HitVector_Iterator hitIt=m_mcHitCol.begin();hitIt!=m_mcHitCol.end();hitIt++){
2696 hitIt->print();
2697 }
2698 cout<<endl;
2699 cout<<"========== Digi =========="<<endl;
2700 for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end();hitIt++)
2701 {
2702 hitIt->print();
2703 }
2704 cout<<endl;
2705 /*
2706 {
2707 vector<HoughHit>& hitList = m_mcHitCol;
2708 for(HitVector_Iterator hitIt = hitList.begin(); hitIt != hitList.end();hitIt++){
2709 cout<<setw(4)<<hitIt->getHitID();
2710 cout<<"("<<setw(2)<<hitIt->getLayer()<<", "<<setw(3)<<hitIt->getWire()<<") ";
2711 cout<<"("<<setw(10)<<hitIt->getHitPosition().x()<<", "<<setw(10)<<hitIt->getHitPosition().y()<<", "<<setw(10)<<hitIt->getHitPosition().z()<<") ";
2712 vector<int> trkID = hitIt->getTrkID();
2713 cout<<"trkID:"<<setw(4)<<trkID[0];
2714 cout<<"half:"<<setw(4)<<hitIt->getHalfCircle();
2715 if(hitIt->getHitType()==HoughHit::CGEMMCHIT){
2716 if(hitIt->getCgemCluster()!=NULL){
2717 //cout<<hitIt->getFlag()<<" ";
2718 //if(hitIt->getFlag()==2)
2719 cout<<setw(5)<<hitIt->getCgemCluster()->getclusterid();
2720 cout<<"["<<setw(3)<<hitIt->getCgemCluster()->getclusterflagb()<<", "<<setw(3)<<hitIt->getCgemCluster()->getclusterflage()<<"] ";
2721 }
2722 //else cout<<"NULL ";
2723 }
2724 else{
2725 if(hitIt->getDigi()!=NULL){
2726 cout<<setw(5)<<hitIt->getDigi()->getTrackIndex();
2727 cout<<hitIt->getDigi()->identify();
2728 }
2729 //else cout<<"NULL ";
2730
2731 }
2732 cout<<endl;
2733 }
2734 }
2735 cout<<endl;
2736 //*/
2737 /*
2738 {
2739 vector<HoughHit>& hitList = m_houghHitList;
2740 cout<<hitList.size()<<endl;
2741 for(HitVector_Iterator hitIt = hitList.begin(); hitIt != hitList.end();hitIt++){
2742 cout<<setw(4)<<hitIt->getHitID();
2743 cout<<"("<<setw(2)<<hitIt->getLayer()<<", "<<setw(3)<<hitIt->getWire()<<") ";
2744 cout<<"("<<setw(10)<<hitIt->getHitPosition().x()<<", "<<setw(10)<<hitIt->getHitPosition().y()<<", "<<setw(10)<<hitIt->getHitPosition().z()<<") ";
2745 cout<<setw(4)<<hitIt->getFlag();
2746 if(hitIt->getHitType()==HoughHit::CGEMHIT){
2747 cout<<setw(4)<<hitIt->getCgemCluster()->getclusterid();
2748 cout<<"["<<setw(4)<<hitIt->getCgemCluster()->getclusterflagb()<<", "<<setw(4)<<hitIt->getCgemCluster()->getclusterflage()<<"] ";
2749 }
2750 else{
2751 cout<<setw(5)<<hitIt->getDigi()->getTrackIndex();
2752 cout<<hitIt->getDigi()->identify()<<" ";
2753 //cout<<"("<<hitIt->getPairHit()->getLayer()<<", "<<hitIt->getPairHit()->getWire()<<") ";
2754 //if(hitIt->getPairHit()!=NULL)cout<<hitIt->getDriftDist()<<" "<<hitIt->getPairHit()->getDriftDist()<<" "<<hitIt->getDriftDist()-hitIt->getPairHit()->getDriftDist()<<endl;
2755 }
2756 if(hitIt->getPairHit()!=NULL)cout<<hitIt->getPairHit()->getHitID();
2757 //else cout<<"NULL";
2758 cout<<endl;
2759 }
2760 }
2761 cout<<endl;
2762 //*/
2763}
2764
2765
2766int HoughFinder::searchCircle()
2767{
2768 //m_debug=false; // m_debug=true;
2769 if(m_debug) {
2770 cout<<"***************************************"<<endl;
2771 cout<<"start circle search : "<<endl;
2772 }
2773 bool tryCgemLeftOnly=false;
2774 while(true)
2775 {
2776 int charge=0;
2777 // --- fill rough map
2778 int nXhitsLeft=0;
2779 m_roughRhoThetaMap.Reset();
2780 //int charge=0;
2781 int vote=1;
2782
2783 nXhitsLeft = fillHistogram(m_XHoughHitList,&m_roughRhoThetaMap,charge,vote);
2784 if(m_debug) {
2785 cout<<"------------------------------------------------------------------------------------ "<<endl;
2786 cout<<"nXhitsLeft = "<<nXhitsLeft<<endl;
2787 }
2788
2789 if(nXhitsLeft<3) {
2790 if(!tryCgemLeftOnly) {
2791 int nCgemHitsLeft = activeUnusedCgemHitsOnly(m_XHoughHitList);
2792 tryCgemLeftOnly=true;
2793 if(nCgemHitsLeft>=3) {
2794 if(m_debug) cout<<"Now try CGEM left only"<<endl;
2795 continue;
2796 }
2797 }
2798 break;
2799 }
2800
2801 // --- find the peak position on rough Hough map
2802 double x_peak(0.), y_peak(0.);
2803 double x_weight(0.), y_weight(0.);
2804 getWeightedPeak(m_roughRhoThetaMap, x_peak, y_peak, x_weight, y_weight);
2805 if(m_debug) cout<<"rough x_weight, y_weight = "<<x_weight<<", "<<y_weight<<endl;
2806 //cout<<endl;
2807
2808 // --- range for fine map
2809 x_peak = x_weight;
2810 y_peak = y_weight;
2811 double x_min=x_peak-0.5*m_ExtPeak_theta*M_PI/m_nBinTheta;
2812 double x_max=x_peak+0.5*m_ExtPeak_theta*M_PI/m_nBinTheta;
2813 double y_min=y_peak-0.5*m_ExtPeak_rho*0.2/m_nBinRho;
2814 double y_max=y_peak+0.5*m_ExtPeak_rho*0.2/m_nBinRho;
2815
2816 // --- binning for fine map
2817 int nFineBin_Theta = nFineBinTheta(y_peak);
2818 int nFineBin_Rho = nFineBinRho(y_peak);
2819 if(m_debug) cout<<"nThetaFine, nRhoFine = "<<nFineBin_Theta<<", "<<nFineBin_Rho<<endl;
2820
2821 int trkFoundThisTry = 0;
2822 int NHitTried = 0;
2823 // --- loop -/+ hypotheses
2824 for(charge=-1; charge<2; charge+=2)
2825 {
2826
2827 if(m_debug) cout<<"/******** Try iCharge = "<<charge<<" **********/"<<endl;
2828
2829 // create a rough track candidate
2830 int trkID = m_houghTrackList.size();
2831 double binTheta(0.), binRho(0.);
2832 HoughTrack track(charge, x_peak, y_peak, binTheta, binRho, trkID);
2833 if(m_debug) cout<<"center = "<<track.center()<<endl;
2834 // --- fill the fine map
2835 m_fineRhoThetaMap.Reset();
2836 m_fineRhoThetaMap.SetBins(nFineBin_Theta,x_min,x_max,nFineBin_Rho,y_min,y_max);
2837 vote=1;
2838 int nHitFilled(0);
2839 //track.print();
2840 if(fabs(y_weight)<0.02) // > 300 MeV
2841 {
2842 //cout<<"fill high pt histogram"<<endl;
2843 nHitFilled = fillHistogram(m_XHoughHitList,&m_fineRhoThetaMap,0,vote);
2844 }
2845 else
2846 {
2847 //cout<<"fill low pt histogram"<<endl;
2848 nHitFilled = fillHistogram(m_XHoughHitList,&m_fineRhoThetaMap,charge,vote);
2849 }
2850 if(nHitFilled<3) {
2851 if(m_debug) cout<<"nHitFilled<3 -> continue"<<endl;
2852 //delete track;
2853 continue;
2854 }
2855 double theta_peak, rho_peak, theta_weight, rho_weight;
2856 getWeightedPeak(m_fineRhoThetaMap, theta_peak, rho_peak, theta_weight, rho_weight);
2857
2858 // update the track candidate
2859 track.update(theta_weight,rho_weight);
2860 //cout<<" ============== before fitting ================="<<endl;
2861 //cout<<"fine x_weight, y_weight = "<<theta_weight<<", "<<rho_weight<<endl;
2862 //track.print();
2863
2864 // --- select X-view hits
2865 //vector<HoughHit*> XHitList_sel;
2866 int flag=0;// only X-view hits
2867 track.resetNhitHalf();
2868 int nXsel(0);
2869 if(fabs(y_weight)<0.02) // > 300 MeV
2870 {
2871 track.findXHot(m_XHoughHitList,0);
2872 //if(track->getNhitFirstHalf()<track->getNhitSecondHalf())
2873 if(track.getNhitUnusedFirstHalf()<track.getNhitUnusedSecondHalf())
2874 {
2875 if(m_debug) cout<<"--- the second half has more hits, try next charge: "<<endl;
2876 //delete track;
2877 continue;
2878 }
2879 else
2880 {
2881 track.clearHits();
2882 //XHitList_sel.clear();
2883 if(m_debug) cout<<" --- the first half has more hits, go ahead: "<<endl;
2884 }
2885 }
2886 nXsel = track.findXHot(m_XHoughHitList, charge);
2887 //track.printHot();
2888 NHitTried = track.getNTried();
2889 if(m_debug) cout<<"sel "<<nXsel<<" hits, "<<NHitTried<<" new hits tried"<<endl;
2890
2891 //if(track.isAGoodCircle())
2892 if(nXsel>=3)
2893 {
2894 //track.printHot();
2895 //track.print();
2896 //TrkContextEv* trkContextEv = new TrkContextEv(m_bfield);
2897 //if(m_fitFlag>0)TrkErrCode trkErrCode = track.fitCircle(m_mdcDetector,trkContextEv,m_bunchT0);
2898 int fitFlag_circle = 0;
2899 if(m_fitFlag>0)
2900 {
2901 //TrkErrCode trkErrCode = track.fitCircle(m_mdcDetector,m_trkContextEv,m_bunchT0);
2902 //if(m_debug)
2903 //cout<<"circle helix from TrkFitter"<<track.a()<<", chi2="<<track.getChi2()<<endl;
2904 //cout<<"m_bunchT0 = "<<m_bunchT0<<endl;
2905 //track.update(theta_weight,rho_weight);
2906 //cout<<"to start circle fit"<<endl;
2907 fitFlag_circle = track.fitCircle(&myDotsHelixFitter,m_bunchT0*1e9);
2908 //cout<<"circle fit ends"<<endl;
2909 if(fitFlag_circle==0)
2910 {
2911 HepVector a_circle_fit = myDotsHelixFitter.getHelix();
2912 track.updateCirclePar(a_circle_fit[0], a_circle_fit[1], a_circle_fit[2]);
2913 track.clearHits();
2914 int nXsel_new = track.findXHot(m_XHoughHitList, charge);
2915 if(m_debug) {
2916 cout<<" sel "<<nXsel_new<<" hits after fitCircle"<<endl;
2917 track.printHot();
2918 }
2919 NHitTried = track.getNTried();
2920 //cout<<" ============== after fitting ================="<<endl;
2921 //track.printHot();
2922 //track.print();
2923 }
2924 //else if(m_debug) cout<<"fitCircle failed "<<endl;
2925 }
2926 //track.print();
2927
2928 //*/
2929 //track.markUsedHot(XHitList_sel);
2930 if(track.isAGoodCircle())
2931 {
2932 track.markUsedHot(1);
2933 //track.setMdcHit( &m_mdcHitCol);
2934 m_houghTrackList.push_back(track);
2935 if(m_debug) cout<<"save this circle trk"<<endl;
2936 //track.print();
2937 //track.printHot();
2938 trkFoundThisTry++;
2939 break;
2940 }
2941 //else
2942 //{
2943 // if(m_debug) cout<<"drop this circle trk"<<endl;
2944 // track.markUsedHot(-1);
2945 //}
2946 }
2947 //else
2948 //{
2949 // --- drop the circle track
2950 if(m_debug) cout<<"drop this circle trk"<<endl;
2951 //track.markUsedHot(XHitList_sel, -1);
2952 track.markUsedHot(-1);
2953 //delete track;
2954 //}
2955 }// loop charge
2956 if(NHitTried==0)
2957 {
2958 if(m_debug) cout<<"no new hits tried in this try"<<endl;
2959 if(!tryCgemLeftOnly) {
2960 int nCgemHitsLeft = activeUnusedCgemHitsOnly(m_XHoughHitList);
2961 tryCgemLeftOnly=true;
2962 if(m_debug) cout<<"check nCgemHitsLeft="<<nCgemHitsLeft<<endl;
2963 if(nCgemHitsLeft>=3) {
2964 //cout<<"Now try CGEM left only"<<endl;
2965 continue;
2966 }
2967 }
2968 if(m_debug) cout<<"stop circle trk finding "<<endl;
2969 break;
2970 }
2971 }// while
2972
2973 //m_debug=false;
2974 return 1;
2975}
2976
2977void HoughFinder::getWeightedPeak(TH2D& h, double& x_peak, double& y_peak, double& x_weight, double& y_weight, int x_ext, int y_ext)
2978{
2979 int ix_max, iy_max, iz_max;
2980 h.GetMaximumBin(ix_max, iy_max, iz_max);
2981 int nx=h.GetXaxis()->GetNbins();
2982 int ny=h.GetYaxis()->GetNbins();
2983 //cout<<"peak at "<<ix_max<<", "<<iy_max<<" in a map "<<nx<<" * "<<ny<<endl;
2984
2985 if(ix_max==0&&iy_max==0) {
2986 //cout<<"nx, ny = "<<nx<<", "<<ny<<endl;
2987 for(int i=0; i<nx; i++)
2988 {
2989 //cout<<endl;
2990 for(int j=0; j<ny; j++)
2991 {
2992 double n_tmp = h.GetBinContent(i,j);
2993 //cout<<n_tmp<<" ";
2994 }
2995 }
2996 }
2997
2998 x_peak=h.GetXaxis()->GetBinCenter(ix_max);
2999 y_peak=h.GetYaxis()->GetBinCenter(iy_max);
3000 x_weight=0.; y_weight=0.;
3001 double weight(0.);
3002 if(x_ext>=0&&y_ext>=0) {
3003 for(int i1=ix_max-x_ext; i1<=ix_max+x_ext; i1++)
3004 {
3005 double x_tmp = h.GetXaxis()->GetBinCenter(i1);
3006 int i1_tmp = i1;
3007 //if(i1<1) i1_tmp=m_nBinTheta+i1;
3008 //if(i1>1) i1_tmp=m_nBinTheta+i1;
3009 //if(i1==m_nBinTheta+1) i1_tmp=1;
3010 for(int i2=iy_max-y_ext; i2<=iy_max+y_ext; i2++)
3011 {
3012 double n_tmp = h.GetBinContent(i1_tmp,i2);
3013 if(i2<1||i2>ny) n_tmp=0.;
3014 if(i1<1||i1>nx) n_tmp=0.;
3015 weight+=n_tmp;
3016 double y_tmp = h.GetYaxis()->GetBinCenter(i2);
3017 x_weight+=x_tmp*n_tmp;
3018 y_weight+=y_tmp*n_tmp;
3019 //cout<<"i1,i2="<<i1<<","<<i2<<endl;
3020 }
3021 }
3022 x_weight=x_weight/weight;
3023 y_weight=y_weight/weight;
3024 }
3025 else {
3026 x_weight=x_peak;
3027 y_weight=y_peak;
3028 }
3029}
3030
3031
3032int HoughFinder::nFineBinTheta(double rho)
3033{
3034 double rho_rough = fabs(rho);
3035 double k_theta = (600.-1200.)/0.0335;
3036 int nThetaFine = int(k_theta*rho_rough+1200);
3037 if(nThetaFine<500) nThetaFine=500;
3038 //cout<<"rho_rough, nThetaFine = "<<rho_rough<<", "<<nThetaFine<<endl;
3039 //cout<<"nThetaFine = "<<nThetaFine<<endl;
3040 nThetaFine=ceil(1.0*nThetaFine/m_nBinTheta*m_ExtPeak_theta);
3041 return nThetaFine;
3042}
3043
3044int HoughFinder::nFineBinRho(double rho)
3045{
3046 double rho_rough = fabs(rho);
3047 int nRhoFine = 50;
3048 if(rho_rough<0.0198)
3049 {
3050 double k1_rho = (850.-1200.)/0.0198;
3051 nRhoFine=int(k1_rho*rho_rough+1200);
3052 }
3053 else if(rho_rough<0.03)
3054 {
3055 double k2_rho = (300.-850.)/(0.03-0.0198);
3056 nRhoFine=int(k2_rho*(rho_rough-0.0198)+850.);
3057 }
3058 else
3059 {
3060 double k3_rho = (100.-300.)/(0.056-0.03);
3061 nRhoFine=int(k3_rho*(rho_rough-0.03)+300.);
3062 if(nRhoFine<50) nRhoFine=50;
3063 }
3064 //cout<<"nRhoFine = "<<nRhoFine<<endl;
3065 nRhoFine=ceil(1.0*nRhoFine/m_nBinRho*m_ExtPeak_rho);
3066
3067 return nRhoFine;
3068}
3069
3070void HoughFinder::XhitCutWindow(double rho, int ilayer, double charge, double& cut1, double &cut2)
3071{
3072 rho=fabs(rho);
3073 if(rho>0.07) rho=0.07;
3074 TGraph* g_cut1, *g_cut2;
3075 if(ilayer<=3) {
3076 g_cut1=m_cut1_cgem;
3077 g_cut2=m_cut2_cgem;
3078 }
3079 else if(ilayer<=19) {
3080 g_cut1=m_cut1_ODC1;
3081 g_cut2=m_cut2_ODC1;
3082 }
3083 else if(ilayer<=42) {
3084 g_cut1=m_cut1_ODC2;
3085 g_cut2=m_cut2_ODC2;
3086 }
3087 // charge < 0
3088 cut1=g_cut1->Eval(rho);
3089 cut2=g_cut2->Eval(rho);
3090 if(charge>0) {// charge > 0
3091 double cut=cut1;
3092 cut1=-1*cut2;
3093 cut2=-1*cut;
3094 }
3095 else {// charge = 0
3096 double cut=max(fabs(cut1),fabs(cut2));
3097 cut1=-1*cut;
3098 cut2=cut;
3099 }
3100 cut1=1.1*cut1;
3101 cut2=1.1*cut2;
3102}
3103
3104void HoughFinder::clearTrack(vector<HoughTrack>& trackVector)
3105{
3106 /*
3107 vector<HoughTrack*>::iterator it = trackVector.begin();
3108 for(; it!=trackVector.end(); it++)
3109 {
3110 HoughTrack* trk = (*it);
3111 delete trk;
3112 }
3113 */
3114 trackVector.clear();
3115}
3116
3117void HoughFinder::clearMemory()
3118{
3119 //cout<<"m_mdcHitCol.size()"<<m_mdcHitCol.size()<<endl;
3120 vector<MdcHit*>::iterator imdcHit = m_mdcHitCol.begin();
3121 for(;imdcHit != m_mdcHitCol.end();imdcHit++){
3122 delete (*imdcHit);
3123 }
3124
3125 for(vector<HoughTrack>::iterator trkIT = m_houghTrackList.begin();trkIT!=m_houghTrackList.end();trkIT++){
3126 trkIT->clearMemory();
3127 //TrkRecoTrk* trkRecoTrk = trkIT->getTrkRecoTrk();
3128 //delete trkRecoTrk;
3129 }
3130}
3131
3132
3133int HoughFinder::activeUnusedCgemHitsOnly(vector<HoughHit*>& hitPntList)
3134{
3135 int nCgemLeft(0);
3136 vector<HoughHit*>::iterator iter = hitPntList.begin();
3137 for(; iter!=hitPntList.end(); iter++)
3138 {
3139 //if((*iter)->getPairHit()==NULL) continue;// FIXME
3140 //if((*iter)->getPairHit()->getHalfCircle()!=1)continue;// FIXME
3141 int useOld = (*iter)->getUse();
3142 int iLayer = (*iter)->getLayer();
3143 if(iLayer<3) {
3144 if(useOld==-1||useOld==0) {
3145 (*iter)->setUse(0);
3146 nCgemLeft++;
3147 }
3148 }
3149 else
3150 if(useOld==0) (*iter)->setUse(-1);
3151 }
3152
3153 return nCgemLeft;
3154}
3155
3156// re-associate shard hits, llwang 2018-12-26
3157void HoughFinder::solveSharedHits(vector<HoughHit*>& hitList)
3158{
3159 vector<HoughHit*>::iterator itHit = hitList.begin();
3160 for(; itHit!=hitList.end(); itHit++)
3161 {
3162 //vector<HoughTrack*> trkPntList = it->getTrkPntVec();
3163 //int nTrkSharing = trkPntList.size();
3164 vector<int> trkIdList = (*itHit)->getTrkID();
3165 int nTrkSharing = trkIdList.size();
3166 //cout<<"nTrkSharing = "<<nTrkSharing<<endl;
3167 if(nTrkSharing>1) {
3168 //cout<<"found one hit shared by "<<nTrkSharing<<" trks: layer "<<(*itHit)->getLayer()
3169 //<<", pos("<<(*itHit)->getHitPosition().x()<<","<<(*itHit)->getHitPosition().y()<<") with ";
3170 int trkId_minDelD=-1;
3171 double minResid = 99.0;
3172 //vector<int> trkIdList_toDel;
3173 vector<double> vecRes = (*itHit)->getVecResid();
3174 int nTrk = vecRes.size();
3175 //cout<<nTrk<<" residuals: ";
3176 for(int i=0; i<nTrk; i++)// find out the minimum residual
3177 {
3178 //cout<<" trk "<<trkIdList[i]<<":"<<vecRes[i];
3179 if(minResid>fabs(vecRes[i])) {
3180 minResid=fabs(vecRes[i]);
3181 trkId_minDelD=trkIdList[i];
3182 }
3183 }
3184 //cout<<endl;
3185 //cout<<" trkId_minDelD, minResid = "<<trkId_minDelD<<", "<<minResid<<endl;
3186 for(int i=0; i<nTrk; i++)// remove other associations
3187 {
3188 if(trkIdList[i]!=trkId_minDelD) {
3189 (*itHit)->dropTrkID(trkIdList[i]);// remove the trk from the hit
3190 vector<HoughTrack>::iterator itTrk = getHoughTrkIt(m_houghTrackList,trkIdList[i]);
3191 if(itTrk!=m_houghTrackList.end()) itTrk->dropHitPnt(&(*(*itHit)));// remove the hit from the trk
3192 }
3193 }
3194 trkIdList = (*itHit)->getTrkID();
3195 nTrkSharing = trkIdList.size();
3196 //cout<<"new nTrkSharing = "<<nTrkSharing<<endl;
3197 }// if(nTrkSharing>1)
3198 }// loop hits
3199}
3200
3201
3202vector<HoughTrack>::iterator HoughFinder::getHoughTrkIt(vector<HoughTrack>& houghTrkList, int trkId)
3203{
3204 vector<HoughTrack>::iterator it = houghTrkList.begin();
3205 for(; it!=houghTrkList.end(); it++)
3206 {
3207 if(it->getTrkID()==trkId) break;
3208 }
3209 return it;
3210}
3211
3212
3213int HoughFinder::associateVHits()
3214{
3215 //m_debug=true;
3216 //int trkSize = m_houghTrackList.size();
3217 //cout<<"found "<<trkSize<<" circles on x-y plane"<<endl;
3218 if(m_debug)
3219 {
3220 cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
3221 cout<<" start V hits association: "<<endl;
3222 }
3223 //cout<<__FILE__<<" "<<__LINE__<<endl;
3224 for(vector<HoughTrack>::iterator trkIT = m_houghTrackList.begin(); trkIT!=m_houghTrackList.end(); trkIT++)// loop circles
3225 {
3226 //trkIT->print();
3227 //trkIT->printHot();
3228 //HoughTrack* trkCandi = m_houghTrackList[i];
3229 if(trkIT->getFlag()!=0) continue; // skip bad circle candidates
3230 int nXVhit = 0;
3231 m_roughTanlDzMap.Reset();
3232 int vote = 3;
3233 //nXVhit = fillHistogram(m_VHoughHitList,&m_roughTanlDzMap,0,vote,trkCandi);
3234 nXVhit = fillHistogram(m_VHoughHitList,&m_roughTanlDzMap,0,vote,&(*trkIT));// fill tanLambda dz map
3235 if(m_debug)
3236 cout<<nXVhit<<" nXVhits filled in the sz Hough map"<<endl;
3237 if(nXVhit<2)
3238 {
3239 trkIT->setFlag(3); // not enough V hits
3240 continue;// skip circle with Vhits candidates <2
3241 }
3242
3243 double x_peak(0.), y_peak(0.);
3244 double x_weight(0.), y_weight(0.);
3245 getWeightedPeak(m_roughTanlDzMap, x_peak, y_peak, x_weight, y_weight);
3246 //cout<<"rough x_peak, y_peak = "<<x_peak<<", "<<y_peak<<endl;
3247 //cout<<"rough x_weight, y_weight = "<<x_weight<<", "<<y_weight<<endl;
3248 //int binx(0),biny(0),binz(0);
3249 //m_roughTanlDzMap.GetMaximumBin(binx,biny,binz);
3250 //int peakHeight = m_roughTanlDzMap.GetBinContent(binx,biny);
3251 //cout<<"peakHeight:"<<peakHeight<<endl;
3252 //cout<<endl;
3253 //cout<<endl;
3254 //x_peak = x_weight;
3255 //y_peak = y_weight;
3256
3257 //-----------------------------------------------------------------
3258 //binning study
3259 trkIT->setTanl(x_weight);
3260 trkIT->setDz(y_weight);
3261 trkIT->updateHelix();// update pivot and a
3262 int flag = 1;
3263 int charge = trkIT->getCharge();
3264
3265 //int nHot = trkIT->findHot(m_VHoughHitList, flag, charge,3,m_event);
3266 //trkIT->dropRedundentCgemVHits();
3267 //trkIT->setMdcHit( &m_mdcHitCol);
3268 //trkIT->print();
3269 //int l(20);
3270 //int l = Layer; l = 35;
3271 //while(l<36){
3272 // //while(1)
3273 // //cout<<"----------------------------------------------------------------"<<endl;
3274 // //cout<<"----"<<endl;
3275 // double dr = trkIT->dr();
3276 // double phi0 = trkIT->phi0();
3277 // double kappa = trkIT->kappa();
3278 // double dz = trkIT->dz();
3279 // double tanl = trkIT->tanl();
3280 // double cutFactor=1.0;
3281 // if(nXVhit<4) cutFactor=3;
3282 // //int nHot = trkIT->findVHot(m_VHoughHitList,charge,l,cutFactor);
3283 // int nHot = trkIT->findVHot(m_VHoughHitListOnSZmap,charge,l,cutFactor);
3284 // trkIT->dropRedundentCgemVHits();
3285 // if(m_debug)
3286 // {
3287 // cout<<"track id "<<trkIT->getTrkID()<<", maxLayer="<<l<<", nVhits="<<nHot<<endl;
3288 // trkIT->print();
3289 // }
3290 // //cout<<nHot<<endl;
3291 // //trkIT->print();
3292 // vector<HoughHit*> hot = trkIT->getHotList();
3293 // //cout<<"nHot: "<<hot.size()<<endl;
3294 // //cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
3295 // //trkIT->print();
3296 // //trkIT->printHot();
3297 //
3298 // TrkErrCode trkErrCode;
3299 // if(m_fitFlag>1) {
3300 // trkErrCode = trkIT->fitHelix(m_mdcDetector,m_bfield,m_bunchT0,hot,l);
3301 // //trkIT->fitHelix(&myDotsHelixFitter,m_bunchT0*1e9,m_recCgemClusterColBegin);
3302 // }
3303 // //trkIT->print();
3304 // //cout<<"----------------------------------------------------------------"<<endl;
3305 // bool isFitDivergency(false);
3306 // isFitDivergency = isFitDivergency||fabs(trkIT->dr() -dr )>999.;
3307 // isFitDivergency = isFitDivergency||fabs(trkIT->phi0() -phi0 )>999.;
3308 // isFitDivergency = isFitDivergency||fabs(trkIT->kappa() -kappa)>999.;
3309 // isFitDivergency = isFitDivergency||fabs(trkIT->dz() -dz )>999.;
3310 // isFitDivergency = isFitDivergency||fabs(trkIT->tanl() -tanl )>999.;
3311 // if(isFitDivergency){
3312 // trkIT->setFlag(-2);
3313 // if(m_debug)cout<<"helix fit diverged!"<<endl;
3314 // break;
3315 // }
3316 // if(trkErrCode.failure()){
3317 // //if(m_debug)cout<<__FILE__<<" line"<<__LINE__<<":helix fit failure!"<<endl;
3318 // trkIT->setFlag(-1);
3319 // if(m_debug){
3320 // trkErrCode.print(std::cout);cout<<endl;//cout<<endl;
3321 // cout<<"helix fit failed!"<<endl;
3322 // }
3323 // }
3324 // else{
3325 // trkIT->setFlag(0);
3326 // if(m_mcTruth)findMcTrack(&(*trkIT));
3327 // if(m_debug)cout<<"helix fit good!"<<endl;
3328 // }
3329 // if(m_debug)
3330 // {
3331 // cout<<"track info after fit: "<<endl;
3332 // trkIT->print();
3333 // }
3334 //
3335 // l++;
3336 // if(l>m_maxFireLayer)break;
3337 // //break;
3338 // //if(nHot==0)break;
3339 //}// while
3340 //
3341 //if(trkIT->getFlag()!=0)
3342 //{
3343 // //cout<<"trk flag !=0 !"<<endl;"
3344 // continue;
3345 //}
3346 ////cout<<"nHot: "<<nHot<<endl;
3347 ////trkIT->print();
3348 ////trkIT->printHot();
3349 //double dr = trkIT->dr();
3350 //double phi0 = trkIT->phi0();
3351 //double kappa = trkIT->kappa();
3352 //double dz = trkIT->dz();
3353 //double tanl = trkIT->tanl();
3354 //
3355 //double cutFactor=1.0;
3356 //if(nXVhit<4) cutFactor=3;
3357 ////int nHot = trkIT->findVHot(m_VHoughHitList,charge,l,cutFactor);
3358 //int nHot = trkIT->findVHot(m_VHoughHitListOnSZmap,charge,l,cutFactor);
3359 //trkIT->dropRedundentCgemVHits();
3360 //if(m_debug) cout<<"after helix fit, nVhits="<<nHot<<endl;
3361 //vector<HoughHit*> hot = trkIT->getHotList();
3362 //TrkErrCode trkErrCode = trkIT->fitHelix(m_mdcDetector,m_trkContextEv,m_bunchT0,m_mdcHitCol,hot);
3363 //
3364 //bool isFitDivergency(false);
3365 //isFitDivergency = isFitDivergency||fabs(trkIT->dr() -dr )>999.;
3366 //isFitDivergency = isFitDivergency||fabs(trkIT->phi0() -phi0 )>999.;
3367 //isFitDivergency = isFitDivergency||fabs(trkIT->kappa() -kappa)>999.;
3368 //isFitDivergency = isFitDivergency||fabs(trkIT->dz() -dz )>999.;
3369 //isFitDivergency = isFitDivergency||fabs(trkIT->tanl() -tanl )>999.;
3370 //if(isFitDivergency){
3371 // trkIT->setFlag(-2);
3372 // if(m_debug)cout<<"helix fit diverged!"<<endl;
3373 // break;
3374 //}
3375 //if(trkErrCode.failure()){
3376 // //if(m_debug)cout<<__FILE__<<" line"<<__LINE__<<":helix fit failure!"<<endl;
3377 // trkIT->setFlag(-1);
3378 // if(m_debug)cout<<"helix fit failed!"<<endl;
3379 //}
3380 //else{
3381 // trkIT->setFlag(0);
3382 // if(m_mcTruth)findMcTrack(&(*trkIT));
3383 // if(m_debug)cout<<"helix fit good!"<<endl;
3384 //}
3385 //if(trkIT->getFlag()!=0)
3386 //{
3387 // if(m_debug) cout<<"helix flag !=0 !"<<endl;
3388 //}
3389
3390
3391 int helixFitFlag(0);
3392 int nHot_last=0;
3393 int l=35;
3394 //while(1)
3395 for(int j=0; j<2; j++)
3396 {
3397 double cutFactor=1.0;
3398 if(nXVhit<4) cutFactor=3;
3399 int nHot = trkIT->findVHot(m_VHoughHitListOnSZmap,charge,l,cutFactor);
3400 trkIT->dropRedundentCgemVHits();
3401 helixFitFlag = trkIT->fitHelix(&myDotsHelixFitter,m_bunchT0*1e9,m_recCgemClusterColBegin);// trkIT's m_dr etc. updated already
3402 if(helixFitFlag==0)
3403 {
3404 trkIT->updateHelix();// update pivot and a
3405 //if(nHot_last!=nHot)
3406 //{
3407 //nHot_last=nHot;
3408 //continue;
3409 //}
3410 //else break;
3411 }
3412 else break;
3413 }
3414 if(helixFitFlag!=0) // if fit fails
3415 {
3416 trkIT->setFlag(helixFitFlag);
3417 continue; // skip bad helix fit
3418 }
3419 else // if fit successes
3420 {
3421 vector<RecMdcHit> aRecMdcHitVec=myDotsHelixFitter.makeRecMdcHitVec(1);
3422 trkIT->setRecMdcHitVec(aRecMdcHitVec);
3423 //trkIT->setRecMdcHitVec(myDotsHelixFitter.makeRecMdcHitVec(1));
3424 //cout<<"Helix fit success and aRecMdcHitVec filled "<<endl;
3425 //trkIT->printRecMdcHitVec();
3426 //trkIT->printHot();
3427 }
3428
3429
3430
3431 //cout<<"getId() "<<m_trkContextEv->getId()<<endl;
3432 //trkIT->print();
3433 //cout<<trkIT->dr() -dr <<endl;
3434 //cout<<trkIT->phi0() -phi0 <<endl;
3435 //cout<<trkIT->kappa() -kappa<<endl;
3436 //cout<<trkIT->dz() -dz <<endl;
3437 //cout<<trkIT->tanl() -tanl <<endl;
3438 //trkIT->printHot();
3439 }// loop circles
3440 //cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
3441 //m_debug=false;
3442 return 0;
3443}
3444
3445void HoughFinder::findMcTrack(HoughTrack* track)
3446{
3447 TH1I trkID("track index","",100,0,100);
3448 vector<HoughHit*> hotList = track->getHotList();
3449 for(vector<HoughHit*>::iterator hotIt = hotList.begin(); hotIt != hotList.end(); hotIt++){
3450 if((*hotIt)->getPairHit()!=NULL)trkID.Fill(((*hotIt)->getPairHit()->getTrkID())[0]);
3451 }
3452 int trkid = trkID.GetMaximumBin()-1;
3453 for(vector<HoughTrack>::iterator trkIter = m_mcTrackCol.begin();trkIter!=m_mcTrackCol.end();trkIter++){
3454 if(trkIter->getTrkID()==trkid){
3455 if(m_run<0&&m_mcTruth)track->setMcTrack(&(*trkIter));
3456 break;
3457 }
3458 }
3459}
3460
3461int HoughFinder::nFineBinTanl(double tanl)
3462{
3463}
3464
3465int HoughFinder::nFineBinDz(double tanl)
3466{
3467}
3468
3469void HoughFinder::XVhitCutWindow(double tanl, int ilayer, double charge, double& cut1, double &cut2)
3470{
3471}
3472
3473//if(m_hist)ntuple_hit->write();
3474StatusCode HoughFinder::bookTuple(){
3475 MsgStream log(msgSvc(), name());
3476
3477 NTuplePtr nt1(ntupleSvc(), "mdcHoughFinder/hit");
3478 if ( nt1 ){
3479 ntuple_hit= nt1;
3480 } else {
3481 ntuple_hit= ntupleSvc()->book("mdcHoughFinder/hit", CLID_ColumnWiseTuple, "hit");
3482 if(ntuple_hit){
3483 ntuple_hit->addItem("hit_run", m_hit_run);
3484 ntuple_hit->addItem("hit_event", m_hit_event);
3485
3486 ntuple_hit->addItem("hit_nhit", m_hit_nhit, 0, 10000);
3487 ntuple_hit->addItem("hit_hitID", m_hit_nhit, m_hit_hitID);
3488 ntuple_hit->addItem("hit_hitType", m_hit_nhit, m_hit_hitType);
3489 ntuple_hit->addItem("hit_layer", m_hit_nhit, m_hit_layer);
3490 ntuple_hit->addItem("hit_wire", m_hit_nhit, m_hit_wire);
3491 ntuple_hit->addItem("hit_flag", m_hit_nhit, m_hit_flag);
3492 ntuple_hit->addItem("hit_halfCircle", m_hit_nhit, m_hit_halfCircle);
3493 ntuple_hit->addItem("hit_x", m_hit_nhit, m_hit_x);
3494 ntuple_hit->addItem("hit_y", m_hit_nhit, m_hit_y);
3495 ntuple_hit->addItem("hit_z", m_hit_nhit, m_hit_z);
3496 ntuple_hit->addItem("hit_drift", m_hit_nhit, m_hit_drift);
3497
3498 ntuple_hit->addItem("mcHit_hitID", m_hit_nhit, m_mcHit_hitID);
3499 ntuple_hit->addItem("mcHit_hitType", m_hit_nhit, m_mcHit_hitType);
3500 ntuple_hit->addItem("mcHit_layer", m_hit_nhit, m_mcHit_layer);
3501 ntuple_hit->addItem("mcHit_wire", m_hit_nhit, m_mcHit_wire);
3502 ntuple_hit->addItem("mcHit_flag", m_hit_nhit, m_mcHit_flag);
3503 ntuple_hit->addItem("mcHit_halfCircle", m_hit_nhit, m_mcHit_halfCircle);
3504 ntuple_hit->addItem("mcHit_x", m_hit_nhit, m_mcHit_x);
3505 ntuple_hit->addItem("mcHit_y", m_hit_nhit, m_mcHit_y);
3506 ntuple_hit->addItem("mcHit_z", m_hit_nhit, m_mcHit_z);
3507 ntuple_hit->addItem("mcHit_drift", m_hit_nhit, m_mcHit_drift);
3508 } else { log << MSG::ERROR << "Cannot book tuple mdcHoughFinder/hit" <<endmsg;
3509 return StatusCode::FAILURE;
3510 }
3511 }
3512
3513 NTuplePtr nt2(ntupleSvc(), "mdcHoughFinder/track");
3514 if ( nt2 ){
3515 ntuple_track = nt2;
3516 } else {
3517 ntuple_track = ntupleSvc()->book("mdcHoughFinder/track", CLID_ColumnWiseTuple, "track");
3518 if(ntuple_track){
3519 ntuple_track->addItem("trk_run", m_trk_run);
3520 ntuple_track->addItem("trk_event", m_trk_event);
3521 ntuple_track->addItem("trk_nTrack", m_trk_nTrack);
3522 ntuple_track->addItem("trk_trackID", m_trk_trackID);
3523 ntuple_track->addItem("trk_charge", m_trk_charge);
3524 ntuple_track->addItem("trk_flag", m_trk_flag);
3525 ntuple_track->addItem("trk_angle", m_trk_angle);
3526 ntuple_track->addItem("trk_rho", m_trk_rho);
3527 ntuple_track->addItem("trk_dAngle", m_trk_dAngle);
3528 ntuple_track->addItem("trk_dRho", m_trk_dRho);
3529 ntuple_track->addItem("trk_dTanl", m_trk_dTanl);
3530 ntuple_track->addItem("trk_dDz", m_trk_dDz);
3531 ntuple_track->addItem("trk_Xc", m_trk_Xc);
3532 ntuple_track->addItem("trk_Yc", m_trk_Yc);
3533 ntuple_track->addItem("trk_R", m_trk_R);
3534 ntuple_track->addItem("trk_dr", m_trk_dr);
3535 ntuple_track->addItem("trk_phi0", m_trk_phi0);
3536 ntuple_track->addItem("trk_kappa", m_trk_kappa);
3537 ntuple_track->addItem("trk_dz", m_trk_dz);
3538 ntuple_track->addItem("trk_tanl", m_trk_tanl);
3539 ntuple_track->addItem("trk_pxy", m_trk_pxy);
3540 ntuple_track->addItem("trk_px", m_trk_px);
3541 ntuple_track->addItem("trk_py", m_trk_py);
3542 ntuple_track->addItem("trk_pz", m_trk_pz);
3543 ntuple_track->addItem("trk_p", m_trk_p);
3544 ntuple_track->addItem("trk_phi", m_trk_phi);
3545 ntuple_track->addItem("trk_theta", m_trk_theta);
3546 ntuple_track->addItem("trk_cosTheta", m_trk_cosTheta);
3547 ntuple_track->addItem("trk_vx", m_trk_vx);
3548 ntuple_track->addItem("trk_vy", m_trk_vy);
3549 ntuple_track->addItem("trk_vz", m_trk_vz);
3550 ntuple_track->addItem("trk_vr", m_trk_vr);
3551 ntuple_track->addItem("trk_chi2", m_trk_chi2);
3552 ntuple_track->addItem("trk_fiTerm", m_trk_fiTerm);
3553 ntuple_track->addItem("trk_nhit", m_trk_nhit);
3554 ntuple_track->addItem("trk_ncluster", m_trk_ncluster);
3555 ntuple_track->addItem("trk_stat", m_trk_stat);
3556 ntuple_track->addItem("trk_ndof", m_trk_ndof);
3557 ntuple_track->addItem("trk_nster", m_trk_nster);
3558 ntuple_track->addItem("trk_nlayer", m_trk_nlayer);
3559 ntuple_track->addItem("trk_firstLayer", m_trk_firstLayer);
3560 ntuple_track->addItem("trk_lastLayer", m_trk_lastLayer);
3561 ntuple_track->addItem("trk_nCgemXClusters", m_trk_nCgemXClusters);
3562 ntuple_track->addItem("trk_nCgemVClusters", m_trk_nCgemVClusters);
3563
3564 ntuple_track->addItem("trk_nHot", m_trk_nHot, 0, 10000);
3565 ntuple_track->addItem("hot_hitID", m_trk_nHot, m_hot_hitID);
3566 ntuple_track->addItem("hot_hitType", m_trk_nHot, m_hot_hitType);
3567 ntuple_track->addItem("hot_layer", m_trk_nHot, m_hot_layer);
3568 ntuple_track->addItem("hot_wire", m_trk_nHot, m_hot_wire);
3569 ntuple_track->addItem("hot_flag", m_trk_nHot, m_hot_flag);
3570 ntuple_track->addItem("hot_halfCircle", m_trk_nHot, m_hot_halfCircle);
3571 ntuple_track->addItem("hot_x", m_trk_nHot, m_hot_x);
3572 ntuple_track->addItem("hot_y", m_trk_nHot, m_hot_y);
3573 ntuple_track->addItem("hot_z", m_trk_nHot, m_hot_z);
3574 ntuple_track->addItem("hot_drift", m_trk_nHot, m_hot_drift);
3575 ntuple_track->addItem("hot_residual", m_trk_nHot, m_hot_residual);
3576
3577 ntuple_track->addItem("mcHot_hitID", m_trk_nHot, m_mcHot_hitID);
3578 ntuple_track->addItem("mcHot_hitType", m_trk_nHot, m_mcHot_hitType);
3579 ntuple_track->addItem("mcHot_layer", m_trk_nHot, m_mcHot_layer);
3580 ntuple_track->addItem("mcHot_wire", m_trk_nHot, m_mcHot_wire);
3581 ntuple_track->addItem("mcHot_flag", m_trk_nHot, m_mcHot_flag);
3582 ntuple_track->addItem("mcHot_halfCircle", m_trk_nHot, m_mcHot_halfCircle);
3583 ntuple_track->addItem("mcHot_x", m_trk_nHot, m_mcHot_x);
3584 ntuple_track->addItem("mcHot_y", m_trk_nHot, m_mcHot_y);
3585 ntuple_track->addItem("mcHot_z", m_trk_nHot, m_mcHot_z);
3586 ntuple_track->addItem("mcHot_drift", m_trk_nHot, m_mcHot_drift);
3587
3588 ntuple_track->addItem("mcTrk_trackID", m_mcTrk_trackID);
3589 ntuple_track->addItem("mcTrk_charge", m_mcTrk_charge);
3590 ntuple_track->addItem("mcTrk_flag", m_mcTrk_flag);
3591 ntuple_track->addItem("mcTrk_angle", m_mcTrk_angle);
3592 ntuple_track->addItem("mcTrk_rho", m_mcTrk_rho);
3593 ntuple_track->addItem("mcTrk_dAngle", m_mcTrk_dAngle);
3594 ntuple_track->addItem("mcTrk_dRho", m_mcTrk_dRho);
3595 ntuple_track->addItem("mcTrk_dTanl", m_mcTrk_dTanl);
3596 ntuple_track->addItem("mcTrk_dDz", m_mcTrk_dDz);
3597 ntuple_track->addItem("mcTrk_Xc", m_mcTrk_Xc);
3598 ntuple_track->addItem("mcTrk_Yc", m_mcTrk_Yc);
3599 ntuple_track->addItem("mcTrk_R", m_mcTrk_R);
3600 ntuple_track->addItem("mcTrk_dr", m_mcTrk_dr);
3601 ntuple_track->addItem("mcTrk_phi0", m_mcTrk_phi0);
3602 ntuple_track->addItem("mcTrk_kappa", m_mcTrk_kappa);
3603 ntuple_track->addItem("mcTrk_dz", m_mcTrk_dz);
3604 ntuple_track->addItem("mcTrk_tanl", m_mcTrk_tanl);
3605 ntuple_track->addItem("mcTrk_pxy", m_mcTrk_pxy);
3606 ntuple_track->addItem("mcTrk_px", m_mcTrk_px);
3607 ntuple_track->addItem("mcTrk_py", m_mcTrk_py);
3608 ntuple_track->addItem("mcTrk_pz", m_mcTrk_pz);
3609 ntuple_track->addItem("mcTrk_p", m_mcTrk_p);
3610 ntuple_track->addItem("mcTrk_phi", m_mcTrk_phi);
3611 ntuple_track->addItem("mcTrk_theta", m_mcTrk_theta);
3612 ntuple_track->addItem("mcTrk_cosTheta", m_mcTrk_cosTheta);
3613 ntuple_track->addItem("mcTrk_vx", m_mcTrk_vx);
3614 ntuple_track->addItem("mcTrk_vy", m_mcTrk_vy);
3615 ntuple_track->addItem("mcTrk_vz", m_mcTrk_vz);
3616 ntuple_track->addItem("mcTrk_vr", m_mcTrk_vr);
3617 ntuple_track->addItem("mcTrk_chi2", m_mcTrk_chi2);
3618 ntuple_track->addItem("mcTrk_fiTerm", m_mcTrk_fiTerm);
3619 ntuple_track->addItem("mcTrk_nhit", m_mcTrk_nhit);
3620 ntuple_track->addItem("mcTrk_ncluster", m_mcTrk_ncluster);
3621 ntuple_track->addItem("mcTrk_stat", m_mcTrk_stat);
3622 ntuple_track->addItem("mcTrk_ndof", m_mcTrk_ndof);
3623 ntuple_track->addItem("mcTrk_nster", m_mcTrk_nster);
3624 ntuple_track->addItem("mcTrk_nlayer", m_mcTrk_nlayer);
3625 ntuple_track->addItem("mcTrk_firstLayer", m_mcTrk_firstLayer);
3626 ntuple_track->addItem("mcTrk_lastLayer", m_mcTrk_lastLayer);
3627 ntuple_track->addItem("mcTrk_nCgemXClusters", m_mcTrk_nCgemXClusters);
3628 ntuple_track->addItem("mcTrk_nCgemVClusters", m_mcTrk_nCgemVClusters);
3629
3630 ntuple_track->addItem("mcTrk_nHot", m_mcTrk_nHot, 0, 10000);
3631 ntuple_track->addItem("mcTrkHot_hitID", m_mcTrk_nHot, m_mcTrkHot_hitID);
3632 ntuple_track->addItem("mcTrkHot_hitType", m_mcTrk_nHot, m_mcTrkHot_hitType);
3633 ntuple_track->addItem("mcTrkHot_layer", m_mcTrk_nHot, m_mcTrkHot_layer);
3634 ntuple_track->addItem("mcTrkHot_wire", m_mcTrk_nHot, m_mcTrkHot_wire);
3635 ntuple_track->addItem("mcTrkHot_flag", m_mcTrk_nHot, m_mcTrkHot_flag);
3636 ntuple_track->addItem("mcTrkHot_halfCircle", m_mcTrk_nHot, m_mcTrkHot_halfCircle);
3637 ntuple_track->addItem("mcTrkHot_x", m_mcTrk_nHot, m_mcTrkHot_x);
3638 ntuple_track->addItem("mcTrkHot_y", m_mcTrk_nHot, m_mcTrkHot_y);
3639 ntuple_track->addItem("mcTrkHot_z", m_mcTrk_nHot, m_mcTrkHot_z);
3640 ntuple_track->addItem("mcTrkHot_drift", m_mcTrk_nHot, m_mcTrkHot_drift);
3641 } else {
3642 log << MSG::ERROR << "Cannot book tuple mdcHoughFinder/track" <<endmsg;
3643 return StatusCode::FAILURE;
3644 }
3645 }
3646
3647 NTuplePtr nt3(ntupleSvc(), "mdcHoughFinder/event");
3648 if ( nt3 ){
3649 ntuple_event = nt3;
3650 } else {
3651 ntuple_event = ntupleSvc()->book("mdcHoughFinder/event", CLID_ColumnWiseTuple, "event");
3652 if(ntuple_event){
3653 ntuple_event->addItem("evt_run", m_evt_run);
3654 ntuple_event->addItem("evt_event", m_evt_event);
3655 ntuple_event->addItem("evt_nXCluster", m_evt_nXCluster);
3656 ntuple_event->addItem("evt_nVCluster", m_evt_nVCluster);
3657 ntuple_event->addItem("evt_nXVCluster", m_evt_nXVCluster);
3658 ntuple_event->addItem("evt_nCgemCluster", m_evt_nCGEMCluster);
3659 ntuple_event->addItem("evt_nAxialHit", m_evt_nAxialHit);
3660 ntuple_event->addItem("evt_nStereoHit", m_evt_nStereoHit);
3661 ntuple_event->addItem("evt_nODCHit", m_evt_nODCHit);
3662 ntuple_event->addItem("evt_nHit", m_evt_nHit);
3663
3664 ntuple_event->addItem("evt_nTrack", m_evt_nTrack, 0, 100);
3665 ntuple_event->addItem("evtTrk_trackID", m_evt_nTrack, m_evtTrk_trackID);
3666 ntuple_event->addItem("evtTrk_charge", m_evt_nTrack, m_evtTrk_charge);
3667 ntuple_event->addItem("evtTrk_flag", m_evt_nTrack, m_evtTrk_flag);
3668 ntuple_event->addItem("evtTrk_angle", m_evt_nTrack, m_evtTrk_angle);
3669 ntuple_event->addItem("evtTrk_rho", m_evt_nTrack, m_evtTrk_rho);
3670 ntuple_event->addItem("evtTrk_dAngle", m_evt_nTrack, m_evtTrk_dAngle);
3671 ntuple_event->addItem("evtTrk_dRho", m_evt_nTrack, m_evtTrk_dRho);
3672 ntuple_event->addItem("evtTrk_dTanl", m_evt_nTrack, m_evtTrk_dTanl);
3673 ntuple_event->addItem("evtTrk_dDz", m_evt_nTrack, m_evtTrk_dDz);
3674 ntuple_event->addItem("evtTrk_Xc", m_evt_nTrack, m_evtTrk_Xc);
3675 ntuple_event->addItem("evtTrk_Yc", m_evt_nTrack, m_evtTrk_Yc);
3676 ntuple_event->addItem("evtTrk_R", m_evt_nTrack, m_evtTrk_R);
3677 ntuple_event->addItem("evtTrk_dr", m_evt_nTrack, m_evtTrk_dr);
3678 ntuple_event->addItem("evtTrk_phi0", m_evt_nTrack, m_evtTrk_phi0);
3679 ntuple_event->addItem("evtTrk_kappa", m_evt_nTrack, m_evtTrk_kappa);
3680 ntuple_event->addItem("evtTrk_dz", m_evt_nTrack, m_evtTrk_dz);
3681 ntuple_event->addItem("evtTrk_tanl", m_evt_nTrack, m_evtTrk_tanl);
3682 ntuple_event->addItem("evtTrk_pxy", m_evt_nTrack, m_evtTrk_pxy);
3683 ntuple_event->addItem("evtTrk_px", m_evt_nTrack, m_evtTrk_px);
3684 ntuple_event->addItem("evtTrk_py", m_evt_nTrack, m_evtTrk_py);
3685 ntuple_event->addItem("evtTrk_pz", m_evt_nTrack, m_evtTrk_pz);
3686 ntuple_event->addItem("evtTrk_p", m_evt_nTrack, m_evtTrk_p);
3687 ntuple_event->addItem("evtTrk_phi", m_evt_nTrack, m_evtTrk_phi);
3688 ntuple_event->addItem("evtTrk_theta", m_evt_nTrack, m_evtTrk_theta);
3689 ntuple_event->addItem("evtTrk_cosTheta", m_evt_nTrack, m_evtTrk_cosTheta);
3690 ntuple_event->addItem("evtTrk_vx", m_evt_nTrack, m_evtTrk_vx);
3691 ntuple_event->addItem("evtTrk_vy", m_evt_nTrack, m_evtTrk_vy);
3692 ntuple_event->addItem("evtTrk_vz", m_evt_nTrack, m_evtTrk_vz);
3693 ntuple_event->addItem("evtTrk_vr", m_evt_nTrack, m_evtTrk_vr);
3694 ntuple_event->addItem("evtTrk_chi2", m_evt_nTrack, m_evtTrk_chi2);
3695 ntuple_event->addItem("evtTrk_fiTerm", m_evt_nTrack, m_evtTrk_fiTerm);
3696 ntuple_event->addItem("evtTrk_nhit", m_evt_nTrack, m_evtTrk_nhit);
3697 ntuple_event->addItem("evtTrk_ncluster", m_evt_nTrack, m_evtTrk_ncluster);
3698 ntuple_event->addItem("evtTrk_stat", m_evt_nTrack, m_evtTrk_stat);
3699 ntuple_event->addItem("evtTrk_ndof", m_evt_nTrack, m_evtTrk_ndof);
3700 ntuple_event->addItem("evtTrk_nster", m_evt_nTrack, m_evtTrk_nster);
3701 ntuple_event->addItem("evtTrk_nlayer", m_evt_nTrack, m_evtTrk_nlayer);
3702 ntuple_event->addItem("evtTrk_firstLayer", m_evt_nTrack, m_evtTrk_firstLayer);
3703 ntuple_event->addItem("evtTrk_lastLayer", m_evt_nTrack, m_evtTrk_lastLayer);
3704 ntuple_event->addItem("evtTrk_nCgemXClusters", m_evt_nTrack, m_evtTrk_nCgemXClusters);
3705 ntuple_event->addItem("evtTrk_nCgemVClusters", m_evt_nTrack, m_evtTrk_nCgemVClusters);
3706
3707 ntuple_event->addItem("mcEvtTrk_trackID", m_evt_nTrack, m_mcEvtTrk_trackID);
3708 ntuple_event->addItem("mcEvtTrk_charge", m_evt_nTrack, m_mcEvtTrk_charge);
3709 ntuple_event->addItem("mcEvtTrk_flag", m_evt_nTrack, m_mcEvtTrk_flag);
3710 ntuple_event->addItem("mcEvtTrk_angle", m_evt_nTrack, m_mcEvtTrk_angle);
3711 ntuple_event->addItem("mcEvtTrk_rho", m_evt_nTrack, m_mcEvtTrk_rho);
3712 ntuple_event->addItem("mcEvtTrk_dAngle", m_evt_nTrack, m_mcEvtTrk_dAngle);
3713 ntuple_event->addItem("mcEvtTrk_dRho", m_evt_nTrack, m_mcEvtTrk_dRho);
3714 ntuple_event->addItem("mcEvtTrk_dTanl", m_evt_nTrack, m_mcEvtTrk_dTanl);
3715 ntuple_event->addItem("mcEvtTrk_dDz", m_evt_nTrack, m_mcEvtTrk_dDz);
3716 ntuple_event->addItem("mcEvtTrk_Xc", m_evt_nTrack, m_mcEvtTrk_Xc);
3717 ntuple_event->addItem("mcEvtTrk_Yc", m_evt_nTrack, m_mcEvtTrk_Yc);
3718 ntuple_event->addItem("mcEvtTrk_R", m_evt_nTrack, m_mcEvtTrk_R);
3719 ntuple_event->addItem("mcEvtTrk_dr", m_evt_nTrack, m_mcEvtTrk_dr);
3720 ntuple_event->addItem("mcEvtTrk_phi0", m_evt_nTrack, m_mcEvtTrk_phi0);
3721 ntuple_event->addItem("mcEvtTrk_kappa", m_evt_nTrack, m_mcEvtTrk_kappa);
3722 ntuple_event->addItem("mcEvtTrk_dz", m_evt_nTrack, m_mcEvtTrk_dz);
3723 ntuple_event->addItem("mcEvtTrk_tanl", m_evt_nTrack, m_mcEvtTrk_tanl);
3724 ntuple_event->addItem("mcEvtTrk_pxy", m_evt_nTrack, m_mcEvtTrk_pxy);
3725 ntuple_event->addItem("mcEvtTrk_px", m_evt_nTrack, m_mcEvtTrk_px);
3726 ntuple_event->addItem("mcEvtTrk_py", m_evt_nTrack, m_mcEvtTrk_py);
3727 ntuple_event->addItem("mcEvtTrk_pz", m_evt_nTrack, m_mcEvtTrk_pz);
3728 ntuple_event->addItem("mcEvtTrk_p", m_evt_nTrack, m_mcEvtTrk_p);
3729 ntuple_event->addItem("mcEvtTrk_phi", m_evt_nTrack, m_mcEvtTrk_phi);
3730 ntuple_event->addItem("mcEvtTrk_theta", m_evt_nTrack, m_mcEvtTrk_theta);
3731 ntuple_event->addItem("mcEvtTrk_cosTheta", m_evt_nTrack, m_mcEvtTrk_cosTheta);
3732 ntuple_event->addItem("mcEvtTrk_vx", m_evt_nTrack, m_mcEvtTrk_vx);
3733 ntuple_event->addItem("mcEvtTrk_vy", m_evt_nTrack, m_mcEvtTrk_vy);
3734 ntuple_event->addItem("mcEvtTrk_vz", m_evt_nTrack, m_mcEvtTrk_vz);
3735 ntuple_event->addItem("mcEvtTrk_vr", m_evt_nTrack, m_mcEvtTrk_vr);
3736 ntuple_event->addItem("mcEvtTrk_chi2", m_evt_nTrack, m_mcEvtTrk_chi2);
3737 ntuple_event->addItem("mcEvtTrk_fiTerm", m_evt_nTrack, m_mcEvtTrk_fiTerm);
3738 ntuple_event->addItem("mcEvtTrk_nhit", m_evt_nTrack, m_mcEvtTrk_nhit);
3739 ntuple_event->addItem("mcEvtTrk_ncluster", m_evt_nTrack, m_mcEvtTrk_ncluster);
3740 ntuple_event->addItem("mcEvtTrk_stat", m_evt_nTrack, m_mcEvtTrk_stat);
3741 ntuple_event->addItem("mcEvtTrk_ndof", m_evt_nTrack, m_mcEvtTrk_ndof);
3742 ntuple_event->addItem("mcEvtTrk_nster", m_evt_nTrack, m_mcEvtTrk_nster);
3743 ntuple_event->addItem("mcEvtTrk_nlayer", m_evt_nTrack, m_mcEvtTrk_nlayer);
3744 ntuple_event->addItem("mcEvtTrk_firstLayer", m_evt_nTrack, m_mcEvtTrk_firstLayer);
3745 ntuple_event->addItem("mcEvtTrk_lastLayer", m_evt_nTrack, m_mcEvtTrk_lastLayer);
3746 ntuple_event->addItem("mcEvtTrk_nCgemXClusters", m_evt_nTrack, m_mcEvtTrk_nCgemXClusters);
3747 ntuple_event->addItem("mcEvtTrk_nCgemVClusters", m_evt_nTrack, m_mcEvtTrk_nCgemVClusters);
3748 } else {
3749 log << MSG::ERROR << "Cannot book tuple mdcHoughFinder/event" <<endmsg;
3750 return StatusCode::FAILURE;
3751 }
3752 }
3753 return StatusCode::SUCCESS;
3754}
3755
3756int HoughFinder::dumpHit()
3757{
3758 m_hit_run = m_run;
3759 m_hit_event = m_event;
3760 m_hit_nhit = m_houghHitList.size();
3761 int i(0);
3762 for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end(); hitIt++){
3763 m_hit_hitID[i] = hitIt->getHitID();
3764 m_hit_hitType[i] = hitIt->getHitType();
3765 m_hit_layer[i] = hitIt->getLayer();
3766 m_hit_wire[i] = hitIt->getWire();
3767 m_hit_flag[i] = hitIt->getFlag();
3768 m_hit_halfCircle[i] = hitIt->getHalfCircle();
3769 m_hit_x[i] = hitIt->getHitPosition().x();
3770 m_hit_y[i] = hitIt->getHitPosition().y();
3771 m_hit_z[i] = hitIt->getHitPosition().z();
3772 m_hit_drift[i] = hitIt->getDriftDist();
3773
3774 if(hitIt->getPairHit()!=NULL){
3775 m_mcHit_hitID[i] = hitIt->getPairHit()->getHitID();
3776 m_mcHit_hitType[i] = hitIt->getPairHit()->getHitType();
3777 m_mcHit_layer[i] = hitIt->getPairHit()->getLayer();
3778 m_mcHit_wire[i] = hitIt->getPairHit()->getWire();
3779 m_mcHit_flag[i] = hitIt->getPairHit()->getFlag();
3780 m_mcHit_halfCircle[i] = hitIt->getPairHit()->getHalfCircle();
3781 m_mcHit_x[i] = hitIt->getPairHit()->getHitPosition().x();
3782 m_mcHit_y[i] = hitIt->getPairHit()->getHitPosition().y();
3783 m_mcHit_z[i] = hitIt->getPairHit()->getHitPosition().z();
3784 m_mcHit_drift[i] = hitIt->getPairHit()->getDriftDist();
3785 }
3786 i++;
3787 }
3788 ntuple_hit->write();
3789 return m_houghHitList.size();
3790}
3791
3792int HoughFinder::dumpHoughTrack()
3793{
3794 if(m_houghTrackList.size()==0)return 0;
3795 m_trk_run = m_run;
3796 m_trk_event = m_event;
3797 m_trk_nTrack = m_houghTrackList.size();
3798 HepPoint3D positionOntrack(999,999,999), positionOnDetector(999,999,999);
3799 for(vector<HoughTrack>::iterator trkIt = m_houghTrackList.begin(); trkIt != m_houghTrackList.end(); trkIt++){
3800 m_trk_trackID = trkIt->getTrkID();
3801 m_trk_charge = trkIt->getCharge();
3802 m_trk_flag = trkIt->getFlag();
3803 m_trk_angle = trkIt->getAngle();
3804 m_trk_rho = trkIt->getRho();
3805 m_trk_dAngle = trkIt->getDAngle();
3806 m_trk_dRho = trkIt->getDRho();
3807 m_trk_dTanl = trkIt->getDTanl();
3808 m_trk_dDz = trkIt->getDDz();
3809 m_trk_Xc = trkIt->center().x();
3810 m_trk_Yc = trkIt->center().y();
3811 m_trk_R = trkIt->radius();
3812 m_trk_dr = trkIt->dr();
3813 m_trk_phi0 = trkIt->phi0();
3814 m_trk_kappa = trkIt->kappa();
3815 m_trk_dz = trkIt->dz();
3816 m_trk_tanl = trkIt->tanl();
3817 m_trk_pxy = trkIt->pt();
3818 m_trk_px = trkIt->momentum(0).x();
3819 m_trk_py = trkIt->momentum(0).y();
3820 m_trk_pz = trkIt->momentum(0).z();
3821 m_trk_p = trkIt->momentum(0).mag();
3822 m_trk_phi = trkIt->direction(0).phi();
3823 m_trk_theta = trkIt->direction(0).theta();
3824 m_trk_cosTheta = cos(trkIt->direction(0).theta());
3825 m_trk_vx = trkIt->x(0).x();
3826 m_trk_vy = trkIt->x(0).y();
3827 m_trk_vz = trkIt->x(0).z();
3828 m_trk_vr = trkIt->x(0).perp();
3829 m_trk_chi2 = trkIt->getChi2();
3830 //m_trk_chi2[i] = trkIt->
3831 //m_trk_fiTerm[i] = trkIt->
3832 //m_trk_nhit[i] = trkIt->
3833 //m_trk_ncluster[i] = trkIt->
3834 //m_trk_stat[i] = trkIt->
3835 //m_trk_ndof[i] = trkIt->
3836 //m_trk_nster[i] = trkIt->
3837 //m_trk_nlayer[i] = trkIt->
3838 //m_trk_firstLayer[i] = trkIt->
3839 //m_trk_lastLayer[i] = trkIt->
3840 //m_trk_nCgemXClusters[i] = trkIt->
3841 //m_trk_nCgemVClusters[i] = trkIt->
3842
3843 vector<HoughHit*> hotList = trkIt->getHotList();
3844 m_trk_nHot = hotList.size();
3845 int i(0);
3846 for(vector<HoughHit*>::iterator hotIt = hotList.begin(); hotIt != hotList.end(); hotIt++){
3847 m_hot_hitID[i] = (*hotIt)->getHitID();
3848 m_hot_hitType[i] = (*hotIt)->getHitType();
3849 m_hot_layer[i] = (*hotIt)->getLayer();
3850 m_hot_wire[i] = (*hotIt)->getWire();
3851 m_hot_flag[i] = (*hotIt)->getFlag();
3852 m_hot_halfCircle[i] = (*hotIt)->getHalfCircle();
3853 m_hot_x[i] = (*hotIt)->getHitPosition().x();
3854 m_hot_y[i] = (*hotIt)->getHitPosition().y();
3855 m_hot_z[i] = (*hotIt)->getHitPosition().z();
3856 m_hot_drift[i] = (*hotIt)->getDriftDist();
3857 m_hot_residual[i] = (*hotIt)->residual(&(*trkIt),positionOntrack, positionOnDetector);
3858
3859 if((*hotIt)->getPairHit()!=NULL){
3860 m_mcHot_hitID[i] = (*hotIt)->getPairHit()->getHitID();
3861 m_mcHot_hitType[i] = (*hotIt)->getPairHit()->getHitType();
3862 m_mcHot_layer[i] = (*hotIt)->getPairHit()->getLayer();
3863 m_mcHot_wire[i] = (*hotIt)->getPairHit()->getWire();
3864 m_mcHot_flag[i] = (*hotIt)->getPairHit()->getFlag();
3865 m_mcHot_halfCircle[i] = (*hotIt)->getPairHit()->getHalfCircle();
3866 m_mcHot_x[i] = (*hotIt)->getPairHit()->getHitPosition().x();
3867 m_mcHot_y[i] = (*hotIt)->getPairHit()->getHitPosition().y();
3868 m_mcHot_z[i] = (*hotIt)->getPairHit()->getHitPosition().z();
3869 m_mcHot_drift[i] = (*hotIt)->getPairHit()->getDriftDist();
3870 }
3871 i++;
3872 }
3873
3874 if(trkIt->getMcTrack()!=NULL){
3875 m_mcTrk_trackID = trkIt->getMcTrack()->getTrkID();
3876 m_mcTrk_charge = trkIt->getMcTrack()->getCharge();
3877 m_mcTrk_flag = trkIt->getMcTrack()->getFlag();
3878 m_mcTrk_angle = trkIt->getMcTrack()->getAngle();
3879 m_mcTrk_rho = trkIt->getMcTrack()->getRho();
3880 m_mcTrk_dAngle = trkIt->getMcTrack()->getDAngle();
3881 m_mcTrk_dRho = trkIt->getMcTrack()->getDRho();
3882 m_mcTrk_dTanl = trkIt->getMcTrack()->getDTanl();
3883 m_mcTrk_dDz = trkIt->getMcTrack()->getDDz();
3884 m_mcTrk_Xc = trkIt->getMcTrack()->center().x();
3885 m_mcTrk_Yc = trkIt->getMcTrack()->center().y();
3886 m_mcTrk_R = trkIt->getMcTrack()->radius();
3887 m_mcTrk_dr = trkIt->getMcTrack()->dr();
3888 m_mcTrk_phi0 = trkIt->getMcTrack()->phi0();
3889 m_mcTrk_kappa = trkIt->getMcTrack()->kappa();
3890 m_mcTrk_dz = trkIt->getMcTrack()->dz();
3891 m_mcTrk_tanl = trkIt->getMcTrack()->tanl();
3892 m_mcTrk_pxy = trkIt->getMcTrack()->pt();
3893 m_mcTrk_px = trkIt->getMcTrack()->momentum(0).x();
3894 m_mcTrk_py = trkIt->getMcTrack()->momentum(0).y();
3895 m_mcTrk_pz = trkIt->getMcTrack()->momentum(0).z();
3896 m_mcTrk_p = trkIt->getMcTrack()->momentum(0).mag();
3897 m_mcTrk_phi = trkIt->getMcTrack()->direction(0).phi();
3898 m_mcTrk_theta = trkIt->getMcTrack()->direction(0).theta();
3899 m_mcTrk_cosTheta = cos(trkIt->getMcTrack()->direction(0).theta());
3900 m_mcTrk_vx = trkIt->getMcTrack()->x(0).x();
3901 m_mcTrk_vy = trkIt->getMcTrack()->x(0).y();
3902 m_mcTrk_vz = trkIt->getMcTrack()->x(0).z();
3903 m_mcTrk_vr = trkIt->getMcTrack()->x(0).perp();
3904 //m_mcTrk_chi2 = trkIt->getMcTrack()->
3905 //m_mcTrk_fiTerm = trkIt->getMcTrack()->
3906 //m_mcTrk_nhit = trkIt->getMcTrack()->
3907 //m_mcTrk_ncluster = trkIt->getMcTrack()->
3908 //m_mcTrk_stat = trkIt->getMcTrack()->
3909 //m_mcTrk_ndof = trkIt->getMcTrack()->
3910 //m_mcTrk_nster = trkIt->getMcTrack()->
3911 //m_mcTrk_nlayer = trkIt->getMcTrack()->
3912 //m_mcTrk_firstLayer = trkIt->getMcTrack()->
3913 //m_mcTrk_lastLayer = trkIt->getMcTrack()->
3914 //m_mcTrk_nCgemXClusters = trkIt->getMcTrack()->
3915 //m_mcTrk_nCgemVClusters = trkIt->getMcTrack()->
3916
3917 vector<HoughHit*> mcHotList = trkIt->getMcTrack()->getHotList();
3918 m_mcTrk_nHot = mcHotList.size();
3919 int j(0);
3920 for(vector<HoughHit*>::iterator mcHotIt = mcHotList.begin(); mcHotIt != mcHotList.end(); mcHotIt++){
3921 m_mcTrkHot_hitID[j] = (*mcHotIt)->getHitID();
3922 m_mcTrkHot_hitType[j] = (*mcHotIt)->getHitType();
3923 m_mcTrkHot_layer[j] = (*mcHotIt)->getLayer();
3924 m_mcTrkHot_wire[j] = (*mcHotIt)->getWire();
3925 m_mcTrkHot_flag[j] = (*mcHotIt)->getFlag();
3926 m_mcTrkHot_halfCircle[j] = (*mcHotIt)->getHalfCircle();
3927 m_mcTrkHot_x[j] = (*mcHotIt)->getHitPosition().x();
3928 m_mcTrkHot_y[j] = (*mcHotIt)->getHitPosition().y();
3929 m_mcTrkHot_z[j] = (*mcHotIt)->getHitPosition().z();
3930 m_mcTrkHot_drift[j] = (*mcHotIt)->getDriftDist();
3931 j++;
3932 }
3933 }
3934 ntuple_track->write();
3935 }
3936 return m_houghTrackList.size();
3937}
3938
3939int HoughFinder::dumpHoughEvent()
3940{
3941 if(m_houghTrackList.size()==0)return 0;
3942 m_evt_run = m_run;
3943 m_evt_event = m_event;
3944 m_evt_nTrack = m_houghTrackList.size();
3945 int i(0);
3946 for(vector<HoughTrack>::iterator trkIt = m_houghTrackList.begin(); trkIt != m_houghTrackList.end(); trkIt++){
3947
3948 m_evtTrk_trackID[i] = trkIt->getTrkID();
3949 m_evtTrk_charge[i] = trkIt->getCharge();
3950 m_evtTrk_flag[i] = trkIt->getFlag();
3951 m_evtTrk_angle[i] = trkIt->getAngle();
3952 m_evtTrk_rho[i] = trkIt->getRho();
3953 m_evtTrk_dAngle[i] = trkIt->getDAngle();
3954 m_evtTrk_dRho[i] = trkIt->getDRho();
3955 m_evtTrk_dTanl[i] = trkIt->getDTanl();
3956 m_evtTrk_dDz[i] = trkIt->getDDz();
3957 m_evtTrk_Xc[i] = trkIt->center().x();
3958 m_evtTrk_Yc[i] = trkIt->center().y();
3959 m_evtTrk_R[i] = trkIt->radius();
3960 m_evtTrk_dr[i] = trkIt->dr();
3961 m_evtTrk_phi0[i] = trkIt->phi0();
3962 m_evtTrk_kappa[i] = trkIt->kappa();
3963 m_evtTrk_dz[i] = trkIt->dz();
3964 m_evtTrk_tanl[i] = trkIt->tanl();
3965 m_evtTrk_pxy[i] = trkIt->pt();
3966 m_evtTrk_px[i] = trkIt->momentum(0).x();
3967 m_evtTrk_py[i] = trkIt->momentum(0).y();
3968 m_evtTrk_pz[i] = trkIt->momentum(0).z();
3969 m_evtTrk_p[i] = trkIt->momentum(0).mag();
3970 m_evtTrk_phi[i] = trkIt->direction(0).phi();
3971 m_evtTrk_theta[i] = trkIt->direction(0).theta();
3972 m_evtTrk_cosTheta[i] = cos(trkIt->direction(0).theta());
3973 m_evtTrk_vx[i] = trkIt->x(0).x();
3974 m_evtTrk_vy[i] = trkIt->x(0).y();
3975 m_evtTrk_vz[i] = trkIt->x(0).z();
3976 m_evtTrk_vr[i] = trkIt->x(0).perp();
3977
3978 //m_evtTrk_chi2[i] = trkIt->
3979 //m_evtTrk_fiTerm[i] = trkIt->
3980 //m_evtTrk_nhit[i] = trkIt->
3981 //m_evtTrk_ncluster[i] = trkIt->
3982 //m_evtTrk_stat[i] = trkIt->
3983 //m_evtTrk_ndof[i] = trkIt->
3984 //m_evtTrk_nster[i] = trkIt->
3985 //m_evtTrk_nlayer[i] = trkIt->
3986 //m_evtTrk_firstLayer[i] = trkIt->
3987 //m_evtTrk_lastLayer[i] = trkIt->
3988 //m_evtTrk_nCgemXClusters[i] = trkIt->
3989 //m_evtTrk_nCgemVClusters[i] = trkIt->
3990 if(trkIt->getMcTrack()!=NULL){
3991 m_mcEvtTrk_trackID[i] = trkIt->getMcTrack()->getTrkID();
3992 m_mcEvtTrk_charge[i] = trkIt->getMcTrack()->getCharge();
3993 m_mcEvtTrk_flag[i] = trkIt->getMcTrack()->getFlag();
3994 m_mcEvtTrk_angle[i] = trkIt->getMcTrack()->getAngle();
3995 m_mcEvtTrk_rho[i] = trkIt->getMcTrack()->getRho();
3996 m_mcEvtTrk_dAngle[i] = trkIt->getMcTrack()->getDAngle();
3997 m_mcEvtTrk_dRho[i] = trkIt->getMcTrack()->getDRho();
3998 m_mcEvtTrk_dTanl[i] = trkIt->getMcTrack()->getDTanl();
3999 m_mcEvtTrk_dDz[i] = trkIt->getMcTrack()->getDDz();
4000 m_mcEvtTrk_Xc[i] = trkIt->getMcTrack()->center().x();
4001 m_mcEvtTrk_Yc[i] = trkIt->getMcTrack()->center().y();
4002 m_mcEvtTrk_R[i] = trkIt->getMcTrack()->radius();
4003 m_mcEvtTrk_dr[i] = trkIt->getMcTrack()->dr();
4004 m_mcEvtTrk_phi0[i] = trkIt->getMcTrack()->phi0();
4005 m_mcEvtTrk_kappa[i] = trkIt->getMcTrack()->kappa();
4006 m_mcEvtTrk_dz[i] = trkIt->getMcTrack()->dz();
4007 m_mcEvtTrk_tanl[i] = trkIt->getMcTrack()->tanl();
4008 m_mcEvtTrk_pxy[i] = trkIt->getMcTrack()->pt();
4009 m_mcEvtTrk_px[i] = trkIt->getMcTrack()->momentum(0).x();
4010 m_mcEvtTrk_py[i] = trkIt->getMcTrack()->momentum(0).y();
4011 m_mcEvtTrk_pz[i] = trkIt->getMcTrack()->momentum(0).z();
4012 m_mcEvtTrk_p[i] = trkIt->getMcTrack()->momentum(0).mag();
4013 m_mcEvtTrk_phi[i] = trkIt->getMcTrack()->direction(0).phi();
4014 m_mcEvtTrk_theta[i] = trkIt->getMcTrack()->direction(0).theta();
4015 m_mcEvtTrk_cosTheta[i] = cos(trkIt->getMcTrack()->direction(0).theta());
4016 m_mcEvtTrk_vx[i] = trkIt->getMcTrack()->x(0).x();
4017 m_mcEvtTrk_vy[i] = trkIt->getMcTrack()->x(0).y();
4018 m_mcEvtTrk_vz[i] = trkIt->getMcTrack()->x(0).z();
4019 m_mcEvtTrk_vr[i] = trkIt->getMcTrack()->x(0).perp();
4020 //m_mcEvtTrk_chi2[i] = trkIt->getMcTrack()->
4021 //m_mcEvtTrk_fiTerm[i] = trkIt->getMcTrack()->
4022 //m_mcEvtTrk_nhit[i] = trkIt->getMcTrack()->
4023 //m_mcEvtTrk_ncluster[i] = trkIt->getMcTrack()->
4024 //m_mcEvtTrk_stat[i] = trkIt->getMcTrack()->
4025 //m_mcEvtTrk_ndof[i] = trkIt->getMcTrack()->
4026 //m_mcEvtTrk_nster[i] = trkIt->getMcTrack()->
4027 //m_mcEvtTrk_nlayer[i] = trkIt->getMcTrack()->
4028 //m_mcEvtTrk_firstLayer[i] = trkIt->getMcTrack()->
4029 //m_mcEvtTrk_lastLayer[i] = trkIt->getMcTrack()->
4030 //m_mcEvtTrk_nCgemXClusters[i] = trkIt->getMcTrack()->
4031 //m_mcEvtTrk_nCgemVClusters[i] = trkIt->getMcTrack()->
4032 }
4033 i++;
4034 }
4035 ntuple_event->write();
4036 return m_houghTrackList.size();
4037}
4038
4039int HoughFinder::dumpTdsTrack()
4040{
4041 SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
4042 if(!recMdcTrackCol)return 0;
4043 m_evt_run = m_run;
4044 m_evt_event = m_event;
4045 m_evt_nTrack = recMdcTrackCol->size();
4046 for(RecMdcTrackCol::iterator iter = recMdcTrackCol->begin();iter!=recMdcTrackCol->end();iter++){
4047 m_trk_trackID = (*iter)->trackId() ;
4048 m_trk_charge = (*iter)->charge() ;
4049 //m_trk_flag = (*iter)->
4050 //m_trk_angle = (*iter)->
4051 //m_trk_rho = (*iter)->
4052 //m_trk_dAngle = (*iter)->
4053 //m_trk_dRho = (*iter)->
4054 //m_trk_dTanl = (*iter)->
4055 //m_trk_dDz = (*iter)->
4056 //m_trk_Xc = (*iter)->
4057 //m_trk_Yc = (*iter)->
4058 //m_trk_R = (*iter)->
4059 m_trk_dr = (*iter)->helix(0) ;
4060 m_trk_phi0 = (*iter)->helix(1) ;
4061 m_trk_kappa = (*iter)->helix(2) ;
4062 m_trk_dz = (*iter)->helix(3) ;
4063 m_trk_tanl = (*iter)->helix(4) ;
4064 m_trk_pxy = (*iter)->pxy() ;
4065 m_trk_px = (*iter)->px() ;
4066 m_trk_py = (*iter)->py() ;
4067 m_trk_pz = (*iter)->pz() ;
4068 m_trk_p = (*iter)->p() ;
4069 m_trk_phi = (*iter)->theta() ;
4070 m_trk_theta = (*iter)->phi() ;
4071 m_trk_cosTheta = cos((*iter)->theta()) ;
4072 m_trk_vx = (*iter)->x() ;
4073 m_trk_vy = (*iter)->y() ;
4074 m_trk_vz = (*iter)->z() ;
4075 m_trk_vr = (*iter)->r() ;
4076 m_trk_chi2 = (*iter)->chi2() ;
4077 m_trk_fiTerm = (*iter)->getFiTerm() ;
4078 m_trk_nhit = (*iter)->getNhits() ;
4079 m_trk_ncluster = (*iter)->getNcluster() ;
4080 m_trk_stat = (*iter)->stat() ;
4081 m_trk_ndof = (*iter)->ndof() ;
4082 m_trk_nster = (*iter)->nster() ;
4083 m_trk_nlayer = (*iter)->nlayer() ;
4084 m_trk_firstLayer = (*iter)->firstLayer() ;
4085 m_trk_lastLayer = (*iter)->lastLayer() ;
4086 m_trk_nCgemXClusters = (*iter)->nCgemXCluster();
4087 m_trk_nCgemVClusters = (*iter)->nCgemVCluster();
4088
4089 HitRefVec hitRefVec = (*iter)->getVecHits();
4090 HitRefVec::iterator hotIt = hitRefVec.begin();
4091 ClusterRefVec clusterRefVec = (*iter)->getVecClusters();
4092 ClusterRefVec::iterator clusterIt = clusterRefVec.begin();
4093 m_trk_nHot = clusterRefVec.size() + hitRefVec.size();
4094 int i(0);
4095 for(clusterIt=clusterRefVec.begin();clusterIt!=clusterRefVec.end();clusterIt++){
4096 for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end(); hitIt++){
4097 if(hitIt->getHitType()!=HoughHit::CGEMHIT)continue;
4098 if((*clusterIt)->getclusterid()==(hitIt)->getCgemCluster()->getclusterid()){
4099 m_hot_hitID[i] = (hitIt)->getHitID();
4100 m_hot_hitType[i] = (hitIt)->getHitType();
4101 m_hot_layer[i] = (hitIt)->getLayer();
4102 m_hot_wire[i] = (hitIt)->getWire();
4103 m_hot_flag[i] = (hitIt)->getFlag();
4104 m_hot_halfCircle[i] = (hitIt)->getHalfCircle();
4105 m_hot_x[i] = (hitIt)->getHitPosition().x();
4106 m_hot_y[i] = (hitIt)->getHitPosition().y();
4107 m_hot_z[i] = (hitIt)->getHitPosition().z();
4108 m_hot_drift[i] = (hitIt)->getDriftDist();
4109 m_hot_residual[i] = -999;
4110
4111 if((hitIt)->getPairHit()!=NULL){
4112 m_mcHot_hitID[i] = (hitIt)->getPairHit()->getHitID();
4113 m_mcHot_hitType[i] = (hitIt)->getPairHit()->getHitType();
4114 m_mcHot_layer[i] = (hitIt)->getPairHit()->getLayer();
4115 m_mcHot_wire[i] = (hitIt)->getPairHit()->getWire();
4116 m_mcHot_flag[i] = (hitIt)->getPairHit()->getFlag();
4117 m_mcHot_halfCircle[i] = (hitIt)->getPairHit()->getHalfCircle();
4118 m_mcHot_x[i] = (hitIt)->getPairHit()->getHitPosition().x();
4119 m_mcHot_y[i] = (hitIt)->getPairHit()->getHitPosition().y();
4120 m_mcHot_z[i] = (hitIt)->getPairHit()->getHitPosition().z();
4121 m_mcHot_drift[i] = (hitIt)->getPairHit()->getDriftDist();
4122 }
4123 i++;
4124 }
4125 }
4126 }
4127 for(hotIt=hitRefVec.begin();hotIt!=hitRefVec.end();hotIt++){
4128 int layer = MdcID::layer((*hotIt)->getMdcId());
4129 int wire = MdcID::wire((*hotIt)->getMdcId());
4130 for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end(); hitIt++){
4131 if(hitIt->getHitType()!=HoughHit::MDCHIT)continue;
4132 if(layer == (hitIt)->getLayer() && wire == (hitIt)->getWire()){
4133 m_hot_hitID[i] = (hitIt)->getHitID();
4134 m_hot_hitType[i] = (hitIt)->getHitType();
4135 m_hot_layer[i] = (hitIt)->getLayer();
4136 m_hot_wire[i] = (hitIt)->getWire();
4137 m_hot_flag[i] = (hitIt)->getFlag();
4138 m_hot_halfCircle[i] = (hitIt)->getHalfCircle();
4139 m_hot_x[i] = (hitIt)->getHitPosition().x();
4140 m_hot_y[i] = (hitIt)->getHitPosition().y();
4141 m_hot_z[i] = (hitIt)->getHitPosition().z();
4142 m_hot_drift[i] = (hitIt)->getDriftDist();
4143 m_hot_residual[i] = 999;
4144
4145 if((hitIt)->getPairHit()!=NULL){
4146 m_mcHot_hitID[i] = (hitIt)->getPairHit()->getHitID();
4147 m_mcHot_hitType[i] = (hitIt)->getPairHit()->getHitType();
4148 m_mcHot_layer[i] = (hitIt)->getPairHit()->getLayer();
4149 m_mcHot_wire[i] = (hitIt)->getPairHit()->getWire();
4150 m_mcHot_flag[i] = (hitIt)->getPairHit()->getFlag();
4151 m_mcHot_halfCircle[i] = (hitIt)->getPairHit()->getHalfCircle();
4152 m_mcHot_x[i] = (hitIt)->getPairHit()->getHitPosition().x();
4153 m_mcHot_y[i] = (hitIt)->getPairHit()->getHitPosition().y();
4154 m_mcHot_z[i] = (hitIt)->getPairHit()->getHitPosition().z();
4155 m_mcHot_drift[i] = (hitIt)->getPairHit()->getDriftDist();
4156 }
4157 }
4158 }
4159 i++;
4160 }
4161
4162 int maxSameHitNumber(0);
4163 vector<HoughTrack>::iterator sameTrkIt = m_mcTrackCol.end();
4164 vector<HoughTrack>::iterator trkIter = m_mcTrackCol.begin();
4165 for(;trkIter!=m_mcTrackCol.end();trkIter++){
4166 int sameHitNumber(0);
4167 vector<HoughHit*> hitList = trkIter->getHotList();
4168 for(vector<HoughHit*>::iterator mchotIt = hitList.begin(); mchotIt != hitList.end();mchotIt++){
4169 if((*mchotIt)->getHitType()==HoughHit::CGEMHIT){
4170 for(clusterIt=clusterRefVec.begin();clusterIt!=clusterRefVec.end();clusterIt++){
4171 if((*clusterIt)->getclusterid()==(*mchotIt)->getCgemCluster()->getclusterid()){
4172 sameHitNumber++;
4173 }
4174 }
4175 }else if((*mchotIt)->getHitType()==HoughHit::MDCHIT){
4176 for(hotIt=hitRefVec.begin();hotIt!=hitRefVec.end();hotIt++){
4177 int layer = MdcID::layer((*hotIt)->getMdcId());
4178 int wire = MdcID::wire((*hotIt)->getMdcId());
4179 if(layer == (*mchotIt)->getLayer() && wire == (*mchotIt)->getWire()){
4180 sameHitNumber++;
4181 }
4182 }
4183 }
4184 }
4185 if(sameHitNumber>maxSameHitNumber){
4186 maxSameHitNumber = sameHitNumber;
4187 sameTrkIt = trkIter;
4188 }
4189
4190 }
4191
4192 if(sameTrkIt!=m_mcTrackCol.end()){
4193 m_mcTrk_trackID = sameTrkIt->getTrkID();
4194 m_mcTrk_charge = sameTrkIt->getCharge();
4195 m_mcTrk_flag = sameTrkIt->getFlag();
4196 m_mcTrk_angle = sameTrkIt->getAngle();
4197 m_mcTrk_rho = sameTrkIt->getRho();
4198 m_mcTrk_dAngle = sameTrkIt->getDAngle();
4199 m_mcTrk_dRho = sameTrkIt->getDRho();
4200 m_mcTrk_dTanl = sameTrkIt->getDTanl();
4201 m_mcTrk_dDz = sameTrkIt->getDDz();
4202 m_mcTrk_Xc = sameTrkIt->center().x();
4203 m_mcTrk_Yc = sameTrkIt->center().y();
4204 m_mcTrk_R = sameTrkIt->radius();
4205 m_mcTrk_dr = sameTrkIt->dr();
4206 m_mcTrk_phi0 = sameTrkIt->phi0();
4207 m_mcTrk_kappa = sameTrkIt->kappa();
4208 m_mcTrk_dz = sameTrkIt->dz();
4209 m_mcTrk_tanl = sameTrkIt->tanl();
4210 m_mcTrk_pxy = sameTrkIt->pt();
4211 m_mcTrk_px = sameTrkIt->momentum(0).x();
4212 m_mcTrk_py = sameTrkIt->momentum(0).y();
4213 m_mcTrk_pz = sameTrkIt->momentum(0).z();
4214 m_mcTrk_p = sameTrkIt->momentum(0).mag();
4215 m_mcTrk_phi = sameTrkIt->direction(0).phi();
4216 m_mcTrk_theta = sameTrkIt->direction(0).theta();
4217 m_mcTrk_cosTheta = cos(sameTrkIt->direction(0).theta());
4218 m_mcTrk_vx = sameTrkIt->x(0).x();
4219 m_mcTrk_vy = sameTrkIt->x(0).y();
4220 m_mcTrk_vz = sameTrkIt->x(0).z();
4221 m_mcTrk_vr = sameTrkIt->x(0).perp();
4222 //m_mcTrk_chi2 = sameTrkIt->getChi2();
4223 //m_mcTrk_fiTerm = sameTrkIt->
4224 //m_mcTrk_nhit = (sameTrkIt->getHotList()).size();
4225 //m_mcTrk_ncluster = sameTrkIt->
4226 //m_mcTrk_stat = sameTrkIt->
4227 //m_mcTrk_ndof = sameTrkIt->
4228 //m_mcTrk_nster = sameTrkIt->
4229 //m_mcTrk_nlayer = sameTrkIt->
4230 //m_mcTrk_firstLayer = sameTrkIt->
4231 //m_mcTrk_lastLayer = sameTrkIt->
4232 //m_mcTrk_nC emXClusters = sameTrkIt->
4233 //m_mcTrk_nCgemVClusters = sameTrkIt->
4234 vector<HoughHit*> mcHotList = sameTrkIt->getHotList();
4235 int j(0);
4236 for(vector<HoughHit*>::iterator mcHotIt = mcHotList.begin(); mcHotIt != mcHotList.end(); mcHotIt++){
4237 m_mcTrkHot_hitID[j] = (*mcHotIt)->getHitID();
4238 m_mcTrkHot_hitType[j] = (*mcHotIt)->getHitType();
4239 m_mcTrkHot_layer[j] = (*mcHotIt)->getLayer();
4240 m_mcTrkHot_wire[j] = (*mcHotIt)->getWire();
4241 m_mcTrkHot_flag[j] = (*mcHotIt)->getFlag();
4242 m_mcTrkHot_halfCircle[j] = (*mcHotIt)->getHalfCircle();
4243 m_mcTrkHot_x[j] = (*mcHotIt)->getHitPosition().x();
4244 m_mcTrkHot_y[j] = (*mcHotIt)->getHitPosition().y();
4245 m_mcTrkHot_z[j] = (*mcHotIt)->getHitPosition().z();
4246 m_mcTrkHot_drift[j] = (*mcHotIt)->getDriftDist();
4247 j++;
4248 }
4249 }
4250 ntuple_track->write();
4251 }
4252 return recMdcTrackCol->size();
4253}
4254
4255int HoughFinder::dumpTdsEvent()
4256{
4257 SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
4258 if(!recMdcTrackCol)return 0;
4259 m_evt_run = m_run;
4260 m_evt_event = m_event;
4261 m_evt_nTrack = recMdcTrackCol->size();
4262 int i(0);
4263 for(RecMdcTrackCol::iterator iter = recMdcTrackCol->begin();iter!=recMdcTrackCol->end();iter++){
4264 m_evtTrk_trackID[i] = (*iter)->trackId() ;
4265 m_evtTrk_charge[i] = (*iter)->charge() ;
4266 //m_evtTrk_flag[i] = (*iter)->
4267 //m_evtTrk_angle[i] = (*iter)->
4268 //m_evtTrk_rho[i] = (*iter)->
4269 //m_evtTrk_dAngle[i] = (*iter)->
4270 //m_evtTrk_dRho[i] = (*iter)->
4271 //m_evtTrk_dTanl[i] = (*iter)->
4272 //m_evtTrk_dDz[i] = (*iter)->
4273 //m_evtTrk_Xc[i] = (*iter)->
4274 //m_evtTrk_Yc[i] = (*iter)->
4275 //m_evtTrk_R[i] = (*iter)->
4276 m_evtTrk_dr[i] = (*iter)->helix(0) ;
4277 m_evtTrk_phi0[i] = (*iter)->helix(1) ;
4278 m_evtTrk_kappa[i] = (*iter)->helix(2) ;
4279 m_evtTrk_dz[i] = (*iter)->helix(3) ;
4280 m_evtTrk_tanl[i] = (*iter)->helix(4) ;
4281 m_evtTrk_pxy[i] = (*iter)->pxy() ;
4282 m_evtTrk_px[i] = (*iter)->px() ;
4283 m_evtTrk_py[i] = (*iter)->py() ;
4284 m_evtTrk_pz[i] = (*iter)->pz() ;
4285 m_evtTrk_p[i] = (*iter)->p() ;
4286 m_evtTrk_phi[i] = (*iter)->theta() ;
4287 m_evtTrk_theta[i] = (*iter)->phi() ;
4288 m_evtTrk_cosTheta[i] = cos((*iter)->theta()) ;
4289 m_evtTrk_vx[i] = (*iter)->x() ;
4290 m_evtTrk_vy[i] = (*iter)->y() ;
4291 m_evtTrk_vz[i] = (*iter)->z() ;
4292 m_evtTrk_vr[i] = (*iter)->r() ;
4293 m_evtTrk_chi2[i] = (*iter)->chi2() ;
4294 m_evtTrk_fiTerm[i] = (*iter)->getFiTerm() ;
4295 m_evtTrk_nhit[i] = (*iter)->getNhits() ;
4296 m_evtTrk_ncluster[i] = (*iter)->getNcluster() ;
4297 m_evtTrk_stat[i] = (*iter)->stat() ;
4298 m_evtTrk_ndof[i] = (*iter)->ndof() ;
4299 m_evtTrk_nster[i] = (*iter)->nster() ;
4300 m_evtTrk_nlayer[i] = (*iter)->nlayer() ;
4301 m_evtTrk_firstLayer[i] = (*iter)->firstLayer() ;
4302 m_evtTrk_lastLayer[i] = (*iter)->lastLayer() ;
4303 m_evtTrk_nCgemXClusters[i] = (*iter)->nCgemXCluster();
4304 m_evtTrk_nCgemVClusters[i] = (*iter)->nCgemVCluster();
4305
4306 HitRefVec hitRefVec = (*iter)->getVecHits();
4307 HitRefVec::iterator hitIt = hitRefVec.begin();
4308 ClusterRefVec clusterRefVec = (*iter)->getVecClusters();
4309 ClusterRefVec::iterator clusterIt = clusterRefVec.begin();
4310
4311 int maxSameHitNumber(0);
4312 vector<HoughTrack>::iterator sameTrkIt = m_mcTrackCol.end();
4313 vector<HoughTrack>::iterator trkIter = m_mcTrackCol.begin();
4314 for(;trkIter!=m_mcTrackCol.end();trkIter++){
4315 int sameHitNumber(0);
4316 vector<HoughHit*> hitList = trkIter->getHotList();
4317 for(vector<HoughHit*>::iterator hotIt = hitList.begin(); hotIt != hitList.end();hotIt++){
4318 if((*hotIt)->getHitType()==HoughHit::CGEMHIT){
4319 for(;clusterIt!=clusterRefVec.end();clusterIt++){
4320 if((*clusterIt)->getclusterid()==(*hotIt)->getCgemCluster()->getclusterid()){
4321 sameHitNumber++;
4322 }
4323 }
4324 }else if((*hotIt)->getHitType()==HoughHit::MDCHIT){
4325 for(;hitIt!=hitRefVec.end();hitIt++){
4326 int layer = MdcID::layer((*hitIt)->getMdcId());
4327 int wire = MdcID::wire((*hitIt)->getMdcId());
4328 if(layer == (*hotIt)->getLayer() && wire == (*hotIt)->getWire()){
4329 sameHitNumber++;
4330 }
4331 }
4332 }
4333 }
4334 if(sameHitNumber>maxSameHitNumber){
4335 maxSameHitNumber = sameHitNumber;
4336 sameTrkIt = trkIter;
4337 }
4338
4339 }
4340
4341 if(sameTrkIt!=m_mcTrackCol.end()){
4342 m_mcEvtTrk_trackID[i] = sameTrkIt->getTrkID();
4343 m_mcEvtTrk_charge[i] = sameTrkIt->getCharge();
4344 m_mcEvtTrk_flag[i] = sameTrkIt->getFlag();
4345 m_mcEvtTrk_angle[i] = sameTrkIt->getAngle();
4346 m_mcEvtTrk_rho[i] = sameTrkIt->getRho();
4347 m_mcEvtTrk_dAngle[i] = sameTrkIt->getDAngle();
4348 m_mcEvtTrk_dRho[i] = sameTrkIt->getDRho();
4349 m_mcEvtTrk_dTanl[i] = sameTrkIt->getDTanl();
4350 m_mcEvtTrk_dDz[i] = sameTrkIt->getDDz();
4351 m_mcEvtTrk_Xc[i] = sameTrkIt->center().x();
4352 m_mcEvtTrk_Yc[i] = sameTrkIt->center().y();
4353 m_mcEvtTrk_R[i] = sameTrkIt->radius();
4354 m_mcEvtTrk_dr[i] = sameTrkIt->dr();
4355 m_mcEvtTrk_phi0[i] = sameTrkIt->phi0();
4356 m_mcEvtTrk_kappa[i] = sameTrkIt->kappa();
4357 m_mcEvtTrk_dz[i] = sameTrkIt->dz();
4358 m_mcEvtTrk_tanl[i] = sameTrkIt->tanl();
4359 m_mcEvtTrk_pxy[i] = sameTrkIt->pt();
4360 m_mcEvtTrk_px[i] = sameTrkIt->momentum(0).x();
4361 m_mcEvtTrk_py[i] = sameTrkIt->momentum(0).y();
4362 m_mcEvtTrk_pz[i] = sameTrkIt->momentum(0).z();
4363 m_mcEvtTrk_p[i] = sameTrkIt->momentum(0).mag();
4364 m_mcEvtTrk_phi[i] = sameTrkIt->direction(0).phi();
4365 m_mcEvtTrk_theta[i] = sameTrkIt->direction(0).theta();
4366 m_mcEvtTrk_cosTheta[i] = cos(sameTrkIt->direction(0).theta());
4367 m_mcEvtTrk_vx[i] = sameTrkIt->x(0).x();
4368 m_mcEvtTrk_vy[i] = sameTrkIt->x(0).y();
4369 m_mcEvtTrk_vz[i] = sameTrkIt->x(0).z();
4370 m_mcEvtTrk_vr[i] = sameTrkIt->x(0).perp();
4371 //m_mcEvtTrk_chi2[i] = sameTrkIt->getChi2();
4372 //m_mcEvtTrk_fiTerm[i] = sameTrkIt->
4373 //m_mcEvtTrk_nhit[i] = (sameTrkIt->getHotList()).size();
4374 //m_mcEvtTrk_ncluster[i] = sameTrkIt->
4375 //m_mcEvtTrk_stat[i] = sameTrkIt->
4376 //m_mcEvtTrk_ndof[i] = sameTrkIt->
4377 //m_mcEvtTrk_nster[i] = sameTrkIt->
4378 //m_mcEvtTrk_nlayer[i] = sameTrkIt->
4379 //m_mcEvtTrk_firstLayer[i] = sameTrkIt->
4380 //m_mcEvtTrk_lastLayer[i] = sameTrkIt->
4381 //m_mcEvtTrk_nCgemXClusters[i] = sameTrkIt->
4382 //m_mcEvtTrk_nCgemVClusters[i] = sameTrkIt->
4383 }
4384 i++;
4385 }
4386 ntuple_event->write();
4387 return recMdcTrackCol->size();
4388}
4389
4390int HoughFinder::printTdsTrack()
4391{
4392 SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
4393 if(!recMdcTrackCol)return 0;
4394 for(RecMdcTrackCol::iterator iter = recMdcTrackCol->begin();iter!=recMdcTrackCol->end();iter++){
4395 double dr = (*iter)->helix(0) ;
4396 double phi0 = (*iter)->helix(1) ;
4397 double kappa = (*iter)->helix(2) ;
4398 double dz = (*iter)->helix(3) ;
4399 double tanl = (*iter)->helix(4) ;
4400 cout<<setw(12)<<"dr:" <<setw(15)<<dr
4401 <<setw(12)<<"phi0:" <<setw(15)<<phi0
4402 <<setw(12)<<"kappa:" <<setw(15)<<kappa
4403 <<setw(12)<<"dz:" <<setw(15)<<dz
4404 <<setw(12)<<"tanl:" <<setw(15)<<tanl
4405 <<endl;
4406
4407 ClusterRefVec clusterRefVec = (*iter)->getVecClusters();
4408 ClusterRefVec::iterator clusterIt = clusterRefVec.begin();
4409 for(clusterIt=clusterRefVec.begin();clusterIt!=clusterRefVec.end();clusterIt++){
4410 //cout<<(*clusterIt)->;
4411 }
4412
4413 HitRefVec hitRefVec = (*iter)->getVecHits();
4414 HitRefVec::iterator hotIt = hitRefVec.begin();
4415 //cout<<"("<<"layer"<<","<<"wire"<<") "
4416 cout<<"("<<"l "<<","<<" w "<<") "
4417 //<<"Stat "
4418 //<<"FlagLR "
4419 <<"S "
4420 <<"F "
4421 <<"DriftLeft "
4422 <<" DriftRight "
4423 <<"ErrDrift_L "
4424 <<"ErrDrift_R "
4425 <<"ChisqAdd "
4426 //<<" Tdc "
4427 //<<" Adc "
4428 <<" DriftT "
4429 <<" Doca "
4430 <<" Entra "
4431 <<" Zhit "
4432 <<" FltLen "
4433 <<endl;
4434 for(hotIt=hitRefVec.begin();hotIt!=hitRefVec.end();hotIt++){
4435 int layer = MdcID::layer((*hotIt)->getMdcId());
4436 int wire = MdcID::wire((*hotIt)->getMdcId());
4437 cout<<"("<<setw( 2)<<layer<<","<<setw( 3)<<wire<<") "
4438 <<setw( 3)<<(*hotIt)->getStat()
4439 <<setw( 3)<<(*hotIt)->getFlagLR()
4440 <<setw(12)<<(*hotIt)->getDriftDistLeft()
4441 <<setw(12)<<(*hotIt)->getDriftDistRight()
4442 <<setw(12)<<(*hotIt)->getErrDriftDistLeft()
4443 <<setw(12)<<(*hotIt)->getErrDriftDistRight()
4444 <<setw(12)<<(*hotIt)->getChisqAdd()
4445 //<<setw( 6)<<(*hotIt)->getTdc()
4446 //<<setw( 6)<<(*hotIt)->getAdc()
4447 <<setw(10)<<(*hotIt)->getDriftT()
4448 <<setw(12)<<(*hotIt)->getDoca()
4449 <<setw(12)<<(*hotIt)->getEntra()
4450 <<setw(10)<<(*hotIt)->getZhit()
4451 <<setw(8 )<<(*hotIt)->getFltLen()
4452 <<endl;
4453 }
4454 cout<<endl;
4455 }
4456 return recMdcTrackCol->size();
4457}
4458
4459//cout<<__FILE__<<" "<<__LINE__<<endl;
sin
double sin(const BesAngle a)
Definition BesAngle.h:210
cos
double cos(const BesAngle a)
Definition BesAngle.h:213
sortCluster
bool sortCluster(const RecCgemCluster *clusterA, const RecCgemCluster *clusterB)
Definition CgemMdcFitAlg.cxx:768
x
Double_t x[10]
Definition DataBase/tau_mode.c:57
gr
TGraph * gr
Definition DataBase/tau_mode.c:69
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition EvtStreamInputIterator.hh:101
s
XmlRpcServer s
Definition HelloServer.cpp:11
HitVector_Iterator
vector< HoughHit >::iterator HitVector_Iterator
Definition HoughHit.h:183
moreHot
bool moreHot(HoughTrack trk1, HoughTrack trk2)
Definition Hough.cxx:1568
largerPt
bool largerPt(HoughTrack trk1, HoughTrack trk2)
Definition Hough.cxx:442
sortCluster
bool sortCluster(const RecCgemCluster *clusterA, const RecCgemCluster *clusterB)
Definition Hough.cxx:1892
q
****INTEGER imax DOUBLE PRECISION m_pi *DOUBLE PRECISION m_amfin DOUBLE PRECISION m_Chfin DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_sinw2 DOUBLE PRECISION m_GFermi DOUBLE PRECISION m_MfinMin DOUBLE PRECISION m_ta2 INTEGER m_out INTEGER m_KeyFSR INTEGER m_KeyQCD *COMMON c_Semalib $ !copy of input $ !CMS energy $ !beam mass $ !final mass $ !beam charge $ !final charge $ !smallest final mass $ !Z mass $ !Z width $ !EW mixing angle $ !Gmu Fermi $ alphaQED at q
Definition KKsem.h:33
cut
*********Class see also m_nmax DOUBLE PRECISION m_MasPhot DOUBLE PRECISION m_phsu DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_r2 DOUBLE PRECISION m_WtMass INTEGER m_nmax INTEGER m_Nevgen INTEGER m_IsFSR INTEGER m_MarTot *COMMON c_KarFin $ !Output file $ !Event serial number $ !alpha QED at Thomson limit $ !minimum energy at CMS for remooval $ !infrared cut
Definition KarFin.h:27
weight
*********Class see also m_nmax DOUBLE PRECISION m_MasPhot DOUBLE PRECISION m_phsu DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_r2 DOUBLE PRECISION m_WtMass INTEGER m_nmax INTEGER m_Nevgen INTEGER m_IsFSR INTEGER m_MarTot *COMMON c_KarFin $ !Output file $ !Event serial number $ !alpha QED at Thomson limit $ !minimum energy at CMS for remooval $ !infrared dimensionless $ !dummy photon IR regulator $ !crude photon multiplicity enhancement factor *EVENT $ !MC crude volume of PhhSpace *Sfactors $ !YFS formfactor IR part only $ !YFS formfactor non IR finite part $ !mass weight
Definition KarFin.h:34
MdcDigiVec
std::vector< MdcDigi * > MdcDigiVec
Definition MdcRawDataProvider.h:11
RecMdcHitCol
ObjectVector< RecMdcHit > RecMdcHitCol
Definition RecMdcHit.h:99
ClusterRefVec
SmartRefVector< RecCgemCluster > ClusterRefVec
Definition RecMdcTrack.h:28
RecMdcTrackCol
ObjectVector< RecMdcTrack > RecMdcTrackCol
Definition RecMdcTrack.h:102
HitRefVec
SmartRefVector< RecMdcHit > HitRefVec
Definition RecMdcTrack.h:26
ntupleSvc
INTupleSvc * ntupleSvc()
Definition ServiceAccessor.h:55
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:13
M_PI
#define M_PI
Definition TConstant.h:4
X1
#define X1
X2
#define X2
CgemHitOnTrack::baseHit
const RecCgemCluster * baseHit() const
Definition CgemHitOnTrack.h:119
CgemID::strip
static int strip(const Identifier &id)
Definition CgemID.cxx:83
CgemID::sheet
static int sheet(const Identifier &id)
Definition CgemID.cxx:77
CgemID::is_xstrip
static bool is_xstrip(const Identifier &id)
Definition CgemID.cxx:64
Constants::pi
static const double pi
Definition Constants.h:38
DotsHelixFitter::getHelix
HepVector getHelix()
Definition DotsHelixFitter.h:55
DotsHelixFitter::makeRecMdcHitVec
vector< RecMdcHit > makeRecMdcHitVec(int sel=1)
Definition DotsHelixFitter.cxx:1187
DstMdcTrack::setPxy
void setPxy(const double pxy)
Definition DstMdcTrack.h:86
DstMdcTrack::setTrackId
void setTrackId(const int trackId)
Definition DstMdcTrack.h:84
DstMdcTrack::setPy
void setPy(const double py)
Definition DstMdcTrack.h:88
DstMdcTrack::setZ
void setZ(const double z)
Definition DstMdcTrack.h:95
DstMdcTrack::setX
void setX(const double x)
Definition DstMdcTrack.h:93
DstMdcTrack::setError
void setError(double err[15])
Definition DstMdcTrack.cxx:118
DstMdcTrack::setTheta
void setTheta(const double theta)
Definition DstMdcTrack.h:91
DstMdcTrack::setStat
void setStat(const int stat)
Definition DstMdcTrack.h:97
DstMdcTrack::setP
void setP(const double p)
Definition DstMdcTrack.h:90
DstMdcTrack::setHelix
void setHelix(double helix[5])
Definition DstMdcTrack.cxx:104
DstMdcTrack::setPoca
void setPoca(double poca[3])
Definition DstMdcTrack.cxx:111
DstMdcTrack::setR
void setR(const double r)
Definition DstMdcTrack.h:96
DstMdcTrack::setCharge
void setCharge(const int charge)
Definition DstMdcTrack.h:85
DstMdcTrack::setY
void setY(const double y)
Definition DstMdcTrack.h:94
DstMdcTrack::setChi2
void setChi2(const double chi)
Definition DstMdcTrack.h:98
DstMdcTrack::setPhi
void setPhi(const double phi)
Definition DstMdcTrack.h:92
DstMdcTrack::setPz
void setPz(const double pz)
Definition DstMdcTrack.h:89
DstMdcTrack::setPx
void setPx(const double px)
Definition DstMdcTrack.h:87
Event::CgemMcHit::GetIdentifier
TArrayI GetIdentifier() const
Definition CgemMcHit.h:122
Event::MdcMcHit::identify
Identifier identify() const
Definition MdcMcHit.cxx:24
Helix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition Reconstruction/TrackUtil/TrackUtil-00-00-12/TrackUtil/Helix.h:218
HoughFinder::printTdsTrack
int printTdsTrack()
Definition Hough.cxx:4390
HoughFinder::execute
StatusCode execute()
Definition Hough.cxx:295
HoughFinder::storeFlag
int storeFlag
Definition Hough.h:66
HoughFinder::Layer
int Layer
Definition Hough.h:63
HoughFinder::activeUnusedCgemHitsOnly
int activeUnusedCgemHitsOnly(vector< HoughHit * > &hitPntList)
Definition Hough.cxx:3133
HoughFinder::checkTrack
int checkTrack(vector< HoughTrack > &trackVector)
Definition Hough.cxx:1674
HoughFinder::storeRecTracks
int storeRecTracks(RecMdcTrackCol *trackList, RecMdcHitCol *hitList, vector< HoughTrack > &trackVector)
Definition Hough.cxx:1896
HoughFinder::makeHoughHitList
int makeHoughHitList()
Definition Hough.cxx:447
HoughFinder::finalize
StatusCode finalize()
Definition Hough.cxx:435
HoughFinder::solveSharedHits
void solveSharedHits(vector< HoughHit * > &hitList)
Definition Hough.cxx:3157
HoughFinder::printFlag
int printFlag
Definition Hough.h:68
HoughFinder::storeTrack
int storeTrack(vector< HoughTrack > &trackVector, RecMdcTrackCol *&recMdcTrackCol, RecMdcHitCol *&recMdcHitCol)
Definition Hough.cxx:1855
HoughFinder::storeTracks
int storeTracks(RecMdcTrackCol *trackList, RecMdcHitCol *hitList, vector< HoughTrack > &trackVector)
Definition Hough.cxx:2181
HoughFinder::getHoughTrkIt
vector< HoughTrack >::iterator getHoughTrkIt(vector< HoughTrack > &houghTrkList, int trkId)
Definition Hough.cxx:3202
HoughFinder::bookTuple
StatusCode bookTuple()
Definition Hough.cxx:3474
HoughFinder::getMcHitCol
int getMcHitCol()
Definition Hough.cxx:514
HoughFinder::m_useCgemInGlobalFit
int m_useCgemInGlobalFit
Definition Hough.h:64
HoughFinder::initialize
StatusCode initialize()
Definition Hough.cxx:105
HoughFinder::clearTrack
void clearTrack(vector< HoughTrack > &trackVector)
Definition Hough.cxx:3104
HoughFinder::HoughFinder
HoughFinder(const std::string &name, ISvcLocator *pSvcLocator)
Definition Hough.cxx:31
HoughFinder::getMcParticleCol
int getMcParticleCol()
Definition Hough.cxx:732
HoughFinder::checkHot
int checkHot(vector< HoughTrack > &trackVector)
Definition Hough.cxx:1629
HoughFinder::printHitList
void printHitList()
Definition Hough.cxx:2691
HoughFinder::fillHistogram
int fillHistogram(HoughHit *hit, TH2D *hitMap, int charge, int vote)
Definition Hough.cxx:947
HoughFinder::m_clearTrack
int m_clearTrack
Definition Hough.h:65
HoughFinder::registerTrack
StatusCode registerTrack(RecMdcTrackCol *&trackList_tds, RecMdcHitCol *&hitList_tds)
Definition Hough.cxx:1802
HoughFinder::clearMemory
void clearMemory()
Definition Hough.cxx:3117
HoughHit::setCgemCluster
void setCgemCluster(const RecCgemCluster *cgemCluster)
Definition HoughHit.h:75
HoughHit::getSZ
vector< S_Z > getSZ() const
Definition HoughHit.h:63
HoughHit::setMdcGeomSvc
static void setMdcGeomSvc(MdcGeomSvc *mdcGeomSvc)
Definition HoughHit.h:99
HoughHit::getCgemCluster
const RecCgemCluster * getCgemCluster() const
Definition HoughHit.h:41
HoughHit::getDriftDist
double getDriftDist() const
Definition HoughHit.h:54
HoughHit::setMdcHit
void setMdcHit(MdcHit *mdcHit)
Definition HoughHit.h:135
HoughHit::setDigi
void setDigi(const MdcDigi *mdcDigi)
Definition HoughHit.h:76
HoughHit::getMdcMcHit
const MdcMcHit * getMdcMcHit() const
Definition HoughHit.h:44
HoughHit::getHitPosition
HepPoint3D getHitPosition() const
Definition HoughHit.h:51
HoughHit::getHitType
HitType getHitType() const
Definition HoughHit.h:40
HoughHit::getCgemMcHit
const CgemMcHit * getCgemMcHit() const
Definition HoughHit.h:43
HoughHit::setPairHit
void setPairHit(HoughHit *pairHit)
Definition HoughHit.h:116
HoughHit::setMdcCalibFunSvc
static void setMdcCalibFunSvc(MdcCalibFunSvc *mdcCalibFunSvc)
Definition HoughHit.h:100
HoughHit::getFlag
int getFlag() const
Definition HoughHit.h:47
HoughHit::print
void print()
Definition HoughHit.cxx:327
HoughHit::driftTime
double driftTime()
Definition HoughHit.cxx:296
HoughHit::setCgemCalibFunSvc
static void setCgemCalibFunSvc(CgemCalibFunSvc *cgemCalibFunSvc)
Definition HoughHit.h:102
HoughHit::setWire
void setWire(int wire)
Definition HoughHit.h:80
HoughHit::getLayer
int getLayer() const
Definition HoughHit.h:45
HoughHit::setFlag
void setFlag(int flag)
Definition HoughHit.h:81
HoughHit::setCgemGeomSvc
static void setCgemGeomSvc(CgemGeomSvc *cgemGeomSvc)
Definition HoughHit.h:101
HoughHit::setMdcDetector
static void setMdcDetector(const MdcDetector *mdcDetector)
Definition HoughHit.h:103
HoughHit::getWire
int getWire() const
Definition HoughHit.h:46
HoughHit::CGEMMCHIT
@ CGEMMCHIT
Definition HoughHit.h:27
HoughHit::CGEMHIT
@ CGEMHIT
Definition HoughHit.h:27
HoughHit::MDCMCHIT
@ MDCMCHIT
Definition HoughHit.h:27
HoughTrack::m_useCgemInGlobalFit
static int m_useCgemInGlobalFit
Definition HoughTrack.h:146
HoughTrack::m_clearTrack
static int m_clearTrack
Definition HoughTrack.h:145
HoughTrack::addVecStereoHitPnt
void addVecStereoHitPnt(HoughHit *aHitPnt)
Definition HoughTrack.h:135
HoughTrack::getTrkID
int getTrkID() const
Definition HoughTrack.h:49
HoughTrack::judgeHalf
int judgeHalf(HoughHit *hit)
Definition HoughTrack.cxx:596
HoughTrack::setMcTrack
void setMcTrack(HoughTrack *mcTrack)
Definition HoughTrack.h:142
HoughTrack::getHotList
vector< HoughHit * > getHotList(int type=2)
Definition HoughTrack.cxx:2127
HoughTrack::m_cut
static TGraph * m_cut[2][43]
Definition HoughTrack.h:16
HoughTrack::addHitPnt
void addHitPnt(HoughHit *aHitPnt)
Definition HoughTrack.h:84
Identifier::get_value
value_type get_value() const
Definition Identifier.h:163
MdcHitOnTrack::wireAmbig
int wireAmbig() const
Definition MdcHitOnTrack.cxx:415
MdcHitOnTrack::dcaToWire
double dcaToWire() const
Definition MdcHitOnTrack.cxx:123
MdcHitOnTrack::drift
double drift() const
Definition MdcHitOnTrack.h:75
MdcHitOnTrack::entranceAngle
double entranceAngle() const
Definition MdcHitOnTrack.cxx:172
MdcHitOnTrack::layer
const MdcLayer * layer() const
Definition MdcHitOnTrack.cxx:354
MdcHit::digi
const MdcDigi * digi() const
Definition MdcHit.h:55
MdcHit::layernumber
unsigned layernumber() const
Definition MdcHit.h:61
MdcHit::wirenumber
unsigned wirenumber() const
Definition MdcHit.h:62
MdcHit::adcIndex
unsigned adcIndex() const
Definition MdcHit.h:64
MdcHit::driftTime
double driftTime(double tof, double z) const
Definition MdcHit.cxx:142
MdcHit::tdcIndex
unsigned tdcIndex() const
Definition MdcHit.h:63
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
MdcLayer::view
int view(void) const
Definition MdcLayer.h:28
MdcRecoHitOnTrack::mdcHit
const MdcHit * mdcHit() const
Definition MdcRecoHitOnTrack.cxx:38
MdcTrackListBase::setNoInner
void setNoInner(bool noInnerFlag)
Definition MdcTrackListBase.h:63
Pdt::readMCppTable
static void readMCppTable(std::string filenm)
Definition Pdt.cxx:291
RawDataProviderSvc::getMdcDigiVec
MdcDigiVec & getMdcDigiVec(uint32_t control=0)
Definition RawDataProviderSvc.h:29
RawData::identify
virtual Identifier identify() const
Definition RawData.cxx:15
RecCgemCluster::getclusterid
int getclusterid(void) const
Definition RecCgemCluster.h:54
RecCgemCluster::getlayerid
int getlayerid(void) const
Definition RecCgemCluster.h:56
RecCgemCluster::getflag
int getflag(void) const
Definition RecCgemCluster.h:58
RecCgemCluster::getclusterflagb
int getclusterflagb(void) const
Definition RecCgemCluster.h:73
RecCgemCluster::getclusterflage
int getclusterflage(void) const
Definition RecCgemCluster.h:74
RecMdcHit::setMdcId
void setMdcId(Identifier mdcid)
Definition RecMdcHit.h:67
RecMdcHit::setErrDriftDistRight
void setErrDriftDistRight(double erddr)
Definition RecMdcHit.h:63
RecMdcHit::setFltLen
void setFltLen(double fltLen)
Definition RecMdcHit.h:74
RecMdcHit::getChisqAdd
const double getChisqAdd(void) const
Definition RecMdcHit.h:46
RecMdcHit::setErrDriftDistLeft
void setErrDriftDistLeft(double erddl)
Definition RecMdcHit.h:62
RecMdcHit::setDriftDistLeft
void setDriftDistLeft(double ddl)
Definition RecMdcHit.h:60
RecMdcHit::getMdcId
const Identifier getMdcId(void) const
Definition RecMdcHit.h:49
RecMdcHit::setDoca
void setDoca(double doca)
Definition RecMdcHit.h:71
RecMdcHit::setStat
void setStat(int stat)
Definition RecMdcHit.h:66
RecMdcHit::setTdc
void setTdc(double tdc)
Definition RecMdcHit.h:68
RecMdcHit::setAdc
void setAdc(double adc)
Definition RecMdcHit.h:69
RecMdcHit::setFlagLR
void setFlagLR(int lr)
Definition RecMdcHit.h:65
RecMdcHit::setChisqAdd
void setChisqAdd(double pChisq)
Definition RecMdcHit.h:64
RecMdcHit::getFltLen
const double getFltLen(void) const
Definition RecMdcHit.h:56
RecMdcHit::setZhit
void setZhit(double zhit)
Definition RecMdcHit.h:73
RecMdcHit::setDriftT
void setDriftT(double driftT)
Definition RecMdcHit.h:70
RecMdcHit::setDriftDistRight
void setDriftDistRight(double ddr)
Definition RecMdcHit.h:61
RecMdcHit::setTrkId
void setTrkId(int trkid)
Definition RecMdcHit.h:59
RecMdcHit::setId
void setId(int id)
Definition RecMdcHit.h:58
RecMdcHit::setEntra
void setEntra(double entra)
Definition RecMdcHit.h:72
RecMdcTrack::setVecClusters
void setVecClusters(ClusterRefVec vecclusters)
Definition RecMdcTrack.cxx:86
RecMdcTrack::setPivot
void setPivot(const HepPoint3D &pivot)
Definition RecMdcTrack.h:79
RecMdcTrack::setVZ0
void setVZ0(double z0)
Definition RecMdcTrack.h:76
RecMdcTrack::setVY0
void setVY0(double y0)
Definition RecMdcTrack.h:75
RecMdcTrack::setVecHits
void setVecHits(HitRefVec vechits)
Definition RecMdcTrack.cxx:80
RecMdcTrack::setFiTerm
void setFiTerm(double fiterm)
Definition RecMdcTrack.h:77
RecMdcTrack::setVX0
void setVX0(double x0)
Definition RecMdcTrack.h:74
Trajectory::getInfo
virtual void getInfo(double fltLen, HepPoint3D &pos, Hep3Vector &direction) const =0
TrkAbsFit::print
virtual void print(std::ostream &ostr) const =0
TrkAbsFit::pt
virtual double pt(double fltL=0.) const =0
TrkAbsFit::chisq
virtual double chisq() const =0
TrkAbsFit::charge
virtual int charge() const =0
TrkAbsFit::nDof
virtual int nDof() const =0
TrkAbsFit::traj
virtual const TrkDifTraj & traj() const =0
TrkAbsFit::position
virtual HepPoint3D position(double fltL) const =0
TrkExchangePar::phi0
double phi0() const
Definition TrkExchangePar.h:43
TrkExchangePar::omega
double omega() const
Definition TrkExchangePar.h:44
TrkExchangePar::z0
double z0() const
Definition TrkExchangePar.h:45
TrkExchangePar::params
const HepVector & params() const
Definition TrkExchangePar.h:48
TrkExchangePar::d0
double d0() const
Definition TrkExchangePar.h:42
TrkExchangePar::tanDip
double tanDip() const
Definition TrkExchangePar.h:46
TrkExchangePar::covariance
const HepSymMatrix & covariance() const
Definition TrkExchangePar.h:50
TrkFit::helix
virtual TrkExchangePar helix(double fltL) const =0
TrkHitList::begin
hot_iterator begin() const
Definition TrkHitList.h:45
TrkHitList::nHit
unsigned nHit() const
Definition TrkHitList.h:44
TrkHitList::end
hot_iterator end() const
Definition TrkHitList.h:46
TrkHitOnTrk::fltLen
double fltLen() const
Definition TrkHitOnTrk.h:91
TrkHitOnTrk::resid
double resid(bool exclude=false) const
Definition TrkHitOnTrk.cxx:191
TrkHitOnTrk::hitRms
double hitRms() const
Definition TrkHitOnTrk.h:89
TrkHitOnTrk::getParentRep
const TrkRep * getParentRep() const
Definition TrkHitOnTrk.h:73
TrkHitOnTrk::isActive
bool isActive() const
Definition TrkHitOnTrk.h:200
TrkHitOnTrk::getFitStuff
TrkErrCode getFitStuff(HepVector &derivs, double &deltaChi) const
Definition TrkHitOnTrk.cxx:248
TrkRep::arrivalTime
virtual double arrivalTime(double fltL) const
Definition TrkRep.cxx:192