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