BOSS 7.0.2
BESIII Offline Software System
Loading...
Searching...
No Matches
MdcHoughFinder.cxx
Go to the documentation of this file.
1#include "MdcHoughFinder/MdcHoughFinder.h"
2#include "GaudiKernel/MsgStream.h"
3#include "GaudiKernel/AlgFactory.h"
4#include "GaudiKernel/ISvcLocator.h"
5#include "GaudiKernel/SmartDataPtr.h"
6#include "GaudiKernel/IDataProviderSvc.h"
7#include "GaudiKernel/PropertyMgr.h"
8#include "GaudiKernel/IService.h"
9#include "EventModel/EventHeader.h"
10#include "McTruth/DecayMode.h"
11#include "McTruth/McParticle.h"
12#include "McTruth/MdcMcHit.h"
13#include "RawDataProviderSvc/IRawDataProviderSvc.h"
14#include "RawDataProviderSvc/RawDataProviderSvc.h"
15#include "MdcRawEvent/MdcDigi.h"
16#include "Identifier/Identifier.h"
17#include "Identifier/MdcID.h"
18#include "MdcGeomSvc/IMdcGeomSvc.h"
19#include "MdcGeomSvc/MdcGeoWire.h"
20#include "MdcGeomSvc/MdcGeoLayer.h"
21
22#include <vector>
23#include <iostream>
24#include <math.h>
25#include <map>
26
27#include "EvTimeEvent/RecEsTime.h"
28#include "MdcCalibFunSvc/IMdcCalibFunSvc.h"
29#include "RawDataProviderSvc/MdcRawDataProvider.h"
30#include "TrkBase/TrkFit.h"
31#include "TrkBase/TrkHitList.h"
32#include "TrkBase/TrkExchangePar.h"
33#include "TrkFitter/TrkHelixMaker.h"
34#include "TrkFitter/TrkCircleMaker.h"
35#include "TrkFitter/TrkLineMaker.h"
36#include "TrkFitter/TrkHelixFitter.h"
37#include "MdcxReco/Mdcxprobab.h"
38#include "MdcData/MdcHit.h"
39#include "MdcData/MdcRecoHitOnTrack.h"
40#include "MdcTrkRecon/MdcTrack.h"
41#include "MdcPrintSvc/IMdcPrintSvc.h"
42#include "MdcGeom/EntranceAngle.h"
43#include "TrackUtil/Helix.h"
44#include "GaudiKernel/IPartPropSvc.h"
45#include "MdcRecoUtil/Pdt.h"
46#include "RawEvent/RawDataUtil.h"
47
48#include "TTree.h"
49#include "TH2D.h"
50
51using namespace std;
52using namespace Event;
53
54/////////////////////////////////////////////////////////////////////////////
55
56MdcHoughFinder::MdcHoughFinder(const std::string& name, ISvcLocator* pSvcLocator) :
57 Algorithm(name, pSvcLocator)
58{
59 // Declare the properties
60 declareProperty("trackNumMc", m_trackNum_Mc_set=1);
61 declareProperty("method", m_method=0);
62 declareProperty("debug", m_debug = 0);
63 declareProperty("data", m_data= 0);
64 declareProperty("binx", m_binx= 100);
65 declareProperty("biny", m_biny= 100);
66 declareProperty("peakCellNum", m_peakCellNum);
67 declareProperty("ndev", m_ndev= 3);
68 declareProperty("f_pro", m_fpro= 0.8);
69 declareProperty("f_hit_pro", m_hit_pro= 0.8);
70 //declareProperty("debug", m_debug = 0);
71 declareProperty("pdtFile", m_pdtFile = "pdt.table");
72 declareProperty("pickHits", m_pickHits = true);
73 declareProperty("pid", m_pid = 0);
74 //declareProperty("helixHitsSigma", m_helixHitsSigma);
75 declareProperty("combineTracking",m_combineTracking = true);
76 declareProperty("getDigiFlag", m_getDigiFlag = 0);
77 declareProperty("maxMdcDigi", m_maxMdcDigi = 0);
78 declareProperty("keepBadTdc", m_keepBadTdc = 0);
79 declareProperty("dropHot", m_dropHot= 0);
80 declareProperty("keepUnmatch", m_keepUnmatch = 0);
81 declareProperty("minMdcDigi", m_minMdcDigi = 0);
82
83}
84
85StatusCode MdcHoughFinder::beginRun(){
86
87 //Initailize MdcDetector
88 m_gm = MdcDetector::instance(0);
89 if(NULL == m_gm) return StatusCode::FAILURE;
90
91
92 return StatusCode::SUCCESS;
93}
94
95// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
96StatusCode MdcHoughFinder::initialize(){
97
98 MsgStream log(msgSvc(), name());
99 log << MSG::INFO << "in initialize()" << endreq;
100
101 StatusCode sc;
102
103 //for(int i=0;i<43;i++) TrkHelixFitter::nSigmaCut[i] = m_helixHitsSigma[i];
104
105 //nfailure=0;
106 IPartPropSvc* p_PartPropSvc;
107 static const bool CREATEIFNOTTHERE(true);
108 sc = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
109 if (!sc.isSuccess() || 0 == p_PartPropSvc) {
110 log << MSG::ERROR << " Could not initialize PartPropSvc" << endreq;
111 return sc;
112 }
113 m_particleTable = p_PartPropSvc->PDT();
114
115 // RawData
116 IRawDataProviderSvc* irawDataProviderSvc;
117 sc = service ("RawDataProviderSvc", irawDataProviderSvc);
118 m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*> (irawDataProviderSvc);
119 if ( sc.isFailure() ){
120 log << MSG::FATAL << name()<<" Could not load RawDataProviderSvc!" << endreq;
121 return StatusCode::FAILURE;
122
123 }
124
125 // Geometry
126 IMdcGeomSvc* imdcGeomSvc;
127 sc = service ("MdcGeomSvc", imdcGeomSvc);
128 m_mdcGeomSvc = dynamic_cast<MdcGeomSvc*> (imdcGeomSvc);
129 if ( sc.isFailure() ){
130 log << MSG::FATAL << "Could not load MdcGeoSvc!" << endreq;
131 return StatusCode::FAILURE;
132 }
133
134 sc = service ("MagneticFieldSvc",m_pIMF);
135 if(sc!=StatusCode::SUCCESS) {
136 log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
137 }
138 m_bfield = new BField(m_pIMF);
139 log << MSG::INFO << "field z = "<<m_bfield->bFieldNominal()<< endreq;
140
141 m_context = new TrkContextEv(m_bfield);
142
143 Pdt::readMCppTable(m_pdtFile);
144
145 //Get MdcCalibFunSvc
146 IMdcCalibFunSvc* imdcCalibSvc;
147 sc = service ("MdcCalibFunSvc", imdcCalibSvc);
148 m_mdcCalibFunSvc = dynamic_cast<MdcCalibFunSvc*>(imdcCalibSvc);
149 if ( sc.isFailure() ){
150 log << MSG::FATAL << "Could not load MdcCalibFunSvc!" << endreq;
151 return StatusCode::FAILURE;
152 }
153
154 //initialize ntuplehit
155 NTuplePtr nt(ntupleSvc(), "MdcHoughFinder/hit");
156 if ( nt ){
157 ntuplehit = nt;
158 } else {
159 ntuplehit = ntupleSvc()->book("MdcHoughFinder/hit", CLID_ColumnWiseTuple, "hit");
160 if(ntuplehit){
161 ntuplehit->addItem("eventNum", m_eventNum);
162 ntuplehit->addItem("runNum", m_runNum);
163 ntuplehit->addItem("costaCut", m_cosCut);
164 // test the first layer in hough space
165 // ntuplehit->addItem("Layer1Num", m_nLayerNum,0, 43);
166 // ntuplehit->addItem("layer1X", m_nLayerNum, m_layer1X);
167 // ntuplehit->addItem("layer1Y", m_nLayerNum, m_layer1Y);
168 // ntuplehit->addItem("layer1U", m_nLayerNum, m_layer1U);
169 // ntuplehit->addItem("layer1V", m_nLayerNum, m_layer1V);
170
171
172 ntuplehit->addItem("nHitMc", m_nHit_Mc,0, 6796);
173 ntuplehit->addItem("layerMc", m_nHit_Mc, m_layer_Mc);
174 ntuplehit->addItem("cellMc", m_nHit_Mc, m_cell_Mc);
175
176
177 ntuplehit->addItem("nHit", m_nHit,0, 6796);
178 ntuplehit->addItem("layerNhit", 43, m_layerNhit);
179 // ntuplehit->addItem("hitCol", m_nHit, m_hitCol);
180 ntuplehit->addItem("nCross", m_nCross,0, 125000);
181 ntuplehit->addItem("rho", m_nCross, m_rho);
182 ntuplehit->addItem("theta", m_nCross, m_theta);
183
184 ntuplehit->addItem("hitSignal", m_nHit, m_hitSignal);
185 ntuplehit->addItem("layer", m_nHit, m_layer);
186 ntuplehit->addItem("cell", m_nHit, m_cell);
187 ntuplehit->addItem("x_east", m_nHit, m_x_east);
188 ntuplehit->addItem("y_east", m_nHit, m_y_east);
189 ntuplehit->addItem("z_east", m_nHit, m_z_east);
190 ntuplehit->addItem("x_west", m_nHit, m_x_west);
191 ntuplehit->addItem("y_west", m_nHit, m_y_west);
192 ntuplehit->addItem("z_west", m_nHit, m_z_west);
193 ntuplehit->addItem("x", m_nHit, m_x);
194 ntuplehit->addItem("y", m_nHit, m_y);
195 ntuplehit->addItem("z", m_nHit, m_z);
196 ntuplehit->addItem("u", m_nHit, m_u);
197 ntuplehit->addItem("v", m_nHit, m_v);
198 ntuplehit->addItem("p", m_nHit, m_p);
199 ntuplehit->addItem("slant", m_nHit, m_slant);
200 // ntuplehit->addItem("z_stereo", m_stereohit, m_zStereo);
201 // ntuplehit->addItem("z_num", m_stereohit, m_z_num);
202 ntuplehit->addItem("maxCount", m_maxCount);
203
204 ntuplehit->addItem("peakColumn", m_npeak,0,100);
205 ntuplehit->addItem("peakWidth", m_npeak, m_peakWidth);
206 ntuplehit->addItem("peakHigh", m_npeak, m_peakHigh);
207 ntuplehit->addItem("peakArea", m_npeak, m_peakArea);
208 ntuplehit->addItem("peakAreaLeast", m_areaLeast);
209 ntuplehit->addItem("peakAreaLeastNum", m_areaLeastNum);
210 ntuplehit->addItem("peakHitSelect", m_areaSelectHit);
211 ntuplehit->addItem("peakHitSelectSignal", m_areaSelectHit_signal);
212 // ntuplehit->addItem("mc_x", m_Mcnhit, mc_x);
213 // ntuplehit->addItem("mc_y", m_Mcnhit, mc_y);
214 // ntuplehit->addItem("mc_z", m_Mcnhit, mc_z);
215 ntuplehit->addItem("circleCenterX", m_x_circle);
216 ntuplehit->addItem("circleCenterY", m_y_circle);
217 ntuplehit->addItem("circleR", m_r);
218
219
220 ntuplehit->addItem("zStereoNum", m_zStereoNum,0,1000);
221 ntuplehit->addItem("zStereo", m_zStereoNum, m_zStereo );
222 ntuplehit->addItem("l", m_zStereoNum, m_l);
223
224 ntuplehit->addItem("trackNum_Mc", m_trackNum_Mc, 0,100);
225 ntuplehit->addItem("trackNum", m_trackNum, 0,100);
226 ntuplehit->addItem("d0", m_trackNum, m_d0);
227 ntuplehit->addItem("phi0", m_trackNum, m_phi0);
228 ntuplehit->addItem("omega", m_trackNum, m_omega);
229 ntuplehit->addItem("z0", m_trackNum, m_z0);
230 ntuplehit->addItem("tanl", m_trackNum, m_tanl);
231
232 ntuplehit->addItem("pt_rec", m_trackNum, m_pt);
233 ntuplehit->addItem("pt2_rec", m_trackNum, m_pt2);
234 ntuplehit->addItem("pz_rec", m_trackNum, m_pz);
235 ntuplehit->addItem("p_rec", m_trackNum, m_pxyz);
236
237 ntuplehit->addItem("nHitSignal", m_nHitSignal);
238 ntuplehit->addItem("nHitAxialSignal", m_nHitAxialSignal);
239 ntuplehit->addItem("nHitSignal_select", m_trackNum, m_nHitSignal_select);
240 ntuplehit->addItem("nHitAxialSignal_selcet", m_trackNum, m_nHitAxialSignal_select);
241 ntuplehit->addItem("fitFailure", m_trackNum, m_nFitFailure);
242 ntuplehit->addItem("nFit", m_trackNum, m_2d_nFit);
243 ntuplehit->addItem("3dnFit", m_trackNum, m_3d_nFit);
244 ntuplehit->addItem("nHitAxial", m_nHitAxial);
245 ntuplehit->addItem("nFitSucess", m_nFitSucess);
246 // ntuplehit->addItem("nHitAxialSelect", m_nHitAxialSelect);
247 // ntuplehit->addItem("nHitSelect", m_trackNum, m_nHitSelect);
248
249 ntuplehit->addItem("pidTruth", m_pidTruth);
250 ntuplehit->addItem("costaTruth", m_costaTruth);
251 ntuplehit->addItem("phiTruth", m_phi0Truth);
252 ntuplehit->addItem("drTruth", m_drTruth);
253 ntuplehit->addItem("dzTruth", m_dzTruth);
254 ntuplehit->addItem("ptTruth", m_ptTruth);
255 ntuplehit->addItem("pzTruth", m_pzTruth);
256 ntuplehit->addItem("pTruth", m_pTruth);
257 ntuplehit->addItem("qTruth", m_qTruth);
258
259
260 } else { log << MSG::ERROR << "Cannot book tuple MdcHoughFinder/hit" << endmsg;
261 return StatusCode::FAILURE;
262 }
263 }
264
265
266 if(m_debug>2)TrkHelixFitter::m_debug = true;
267
268 return StatusCode::SUCCESS;
269
270
271}
272
273
274// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
275StatusCode MdcHoughFinder::execute() {
276
277 MsgStream log(msgSvc(), name());
278 log << MSG::INFO << "in execute()" << endreq;
279
280 //event start time
281 SmartDataPtr<RecEsTimeCol> evTimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
282 if (!evTimeCol) {
283 log << MSG::WARNING<< "Could not retrieve RecEsTimeCol , use t0=0" << endreq;
284 m_bunchT0=0.;
285 }else{
286 RecEsTimeCol::iterator iter_evt = evTimeCol->begin();
287 if (iter_evt != evTimeCol->end()){
288 m_bunchT0 = (*iter_evt)->getTest()*1.e-9;//m_bunchT0-s, getTest-ns
289 }
290 }
291
292 //if(m_debug)TrkHelixFitter::m_debug = true;
293 // Get the event header, print out event and run number
294
295 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
296 if (!eventHeader) {
297 log << MSG::FATAL << "Could not find Event Header" << endreq;
298 return( StatusCode::FAILURE);
299 }
300
301
302 DataObject *aTrackCol;
303 DataObject *aRecHitCol;
304 SmartIF<IDataManagerSvc> dataManSvc(eventSvc());
305 if(!m_combineTracking){
306 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
307 if(aTrackCol != NULL) {
308 dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
309 eventSvc()->unregisterObject("/Event/Recon/RecMdcTrackCol");
310 }
311 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
312 if(aRecHitCol != NULL) {
313 dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
314 eventSvc()->unregisterObject("/Event/Recon/RecMdcHitCol");
315 }
316 }
317
318 RecMdcTrackCol* trackList;
319 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
320 if (aTrackCol) {
321 trackList = dynamic_cast<RecMdcTrackCol*> (aTrackCol);
322 }else{
323 trackList = new RecMdcTrackCol;
324 StatusCode sc = eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trackList);
325 if(!sc.isSuccess()) {
326 log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
327 return StatusCode::FAILURE;
328 }
329 }
330 int nTdsTk = (int) trackList->size();
331 int nFitSucess = 0;
332
333 RecMdcHitCol* hitList;
334 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
335 if (aRecHitCol) {
336 hitList = dynamic_cast<RecMdcHitCol*> (aRecHitCol);
337 }else{
338 hitList = new RecMdcHitCol;
339 StatusCode sc = eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitList);
340 if(!sc.isSuccess()) {
341 log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
342 return StatusCode::FAILURE;
343 }
344 }
345
346 //pick the first layer
347
348 // m_nLayerNum=43;
349 // for(int i=0;i<43;i++){
350 // const MdcGeoWire *geowir_first_layer =m_mdcGeomSvc->Wire(i,1);
351 // HepPoint3D eastP = geowir_first_layer->Backward()/10.0;//mm 2 cm
352 // HepPoint3D westP = geowir_first_layer->Forward() /10.0;//mm 2 cm
353 // m_layer1X[i]=(eastP.x()+westP.x())/2;
354 // m_layer1Y[i]=(eastP.y()+westP.y())/2;
355 // m_layer1U[i]=CFtrans((eastP.x()+westP.x())/2.,(eastP.y()+westP.y())/2.);
356 // m_layer1V[i]=CFtrans((eastP.y()+westP.y())/2.,(eastP.x()+westP.x())/2.);
357 //
358 // cout<<"("<<i<<","<<"1"<<")"<<":"<<"("<<m_layer1X[i]<<","<<m_layer1Y[i]<<")"<<endl;
359 // }
360
361 t_eventNum=eventHeader->eventNumber();
362 t_runNum=eventHeader->runNumber();
363 m_eventNum=t_eventNum;
364 m_runNum=t_runNum;
365 log << MSG::INFO << "MdcHoughFinder: retrieved event: " << eventHeader->eventNumber() << " run: " << eventHeader->runNumber() << endreq;
366
367 // get mc particle digi
368 int mc_particle=GetMcInfo();
369
370
371 if(abs(m_costaTruth)<=0.83){m_cosCut=1;}
372 else{m_cosCut=0;}
373
374 //int n_curve=0;
375 //double vec_ptTruth=m_ptTruth;
376 //double n_costaTruth=m_costaTruth;
377 //double z_flight=2*M_PI*(771.4/2)*(vec_ptTruth/tan(n_costaTruth));
378 //if(vec_ptTruth<0.12&&(z_flight<1193.)){
379 ///// n_curve++;
380 //// m_curve=n_curve;
381 //}
382 vector<double> vec_u;
383 vector<double> vec_v;
384 vector<double> vec_p;
385 vector<double> vec_x_east;
386 vector<double> vec_y_east;
387 vector<double> vec_z_east;
388 vector<double> vec_x_west;
389 vector<double> vec_y_west;
390 vector<double> vec_z_west;
391 vector<double> vec_z_Mc;
392 vector<double> vec_x;
393 vector<double> vec_y;
394 vector<double> vec_z;
395 vector<int> vec_layer;
396 vector<int> vec_wire;
397 vector<int> vec_slant;
398 vector<double> vec_zStereo;
399 vector<double> l;
400
401 vector<int> vec_layer_Mc;
402 vector<int> vec_wire_Mc;
403 t_maxP=-999.;
404 t_minP=999.;
405
406 // retrieve Mdc truth hits
407 int t_nHit_Mc;
408 int digiId_Mc=0;
409 int nHitAxial_Mc=0;
410
411 SmartDataPtr<Event::MdcMcHitCol> mcMdcMcHitCol(eventSvc(),"/Event/MC/MdcMcHitCol");
412 if (!mcMdcMcHitCol) {
413 log << MSG::INFO << "Could not find MdcMcHit" << endreq;
414 }else{
415 Event::MdcMcHitCol::iterator iterMcHit = mcMdcMcHitCol->begin();
416
417 for(; iterMcHit != mcMdcMcHitCol->end(); ++iterMcHit,digiId_Mc++) {
418 Identifier mdcid = (*iterMcHit)->identify();
419 int layerId_Mc = MdcID::layer(mdcid);
420 int wireId_Mc = MdcID::wire(mdcid);
421
422 vec_layer_Mc.push_back(layerId_Mc);
423 vec_wire_Mc.push_back(wireId_Mc);
424
425 // cout<<"mcHit: "<<"("<<layerId_Mc<<","<<wireId_Mc<<")"<<endl;
426 //vec_layer.push_back(layerId_Mc);
427 //vec_wire.push_back(wireId_Mc);
428 m_layer_Mc[digiId_Mc]=layerId_Mc;
429 m_cell_Mc[digiId_Mc]=wireId_Mc;
430
431 if(m_data==0){
432 if ((layerId_Mc>=8&&layerId_Mc<=19)||(layerId_Mc>=36)) {
433 nHitAxial_Mc++;}
434 if(m_debug>0) {cout<<"("<<layerId_Mc<<","<<wireId_Mc<<")"<<"nHitAxial_Mc: "<<nHitAxial_Mc<<endl;}
435 m_layer[digiId_Mc]=layerId_Mc;
436 m_cell[digiId_Mc]=wireId_Mc;
437 const MdcGeoWire *geowir =m_mdcGeomSvc->Wire(layerId_Mc,wireId_Mc);
438 HepPoint3D eastP = geowir->Backward()/10.0;//mm 2 cm
439 HepPoint3D westP = geowir->Forward() /10.0;//mm 2 cm
440
441 m_x_east[digiId_Mc]=eastP.x();
442 m_y_east[digiId_Mc]=eastP.y();
443 m_z_east[digiId_Mc]=eastP.z();
444 m_x_west[digiId_Mc]=westP.x();
445 m_y_west[digiId_Mc]=westP.y();
446 m_z_west[digiId_Mc]=westP.z();
447
448 vec_x_east.push_back(eastP.x());
449 vec_y_east.push_back(eastP.y());
450 vec_z_east.push_back(eastP.z());
451 vec_x_west.push_back(westP.x());
452 vec_y_west.push_back(westP.y());
453 vec_z_west.push_back(westP.z());
454
455 double x = (*iterMcHit)->getPositionX()/10.;//mm 2 cm
456 double y = (*iterMcHit)->getPositionY()/10.;//mm 2 cm
457 double z = (*iterMcHit)->getPositionZ()/10.;//mm 2 cm
458
459 double u=CFtrans(x,y);
460 double v=CFtrans(y,x);
461 double p=sqrt(u*u+v*v);
462
463 vec_x.push_back(x);
464 vec_y.push_back(y);
465 vec_z.push_back(z);
466 vec_u.push_back(u);
467 vec_v.push_back(v);
468 vec_p.push_back(p);
469
470 m_x[digiId_Mc]=x;
471 m_y[digiId_Mc]=y;
472 m_z[digiId_Mc]=z;
473 m_u[digiId_Mc]=u;
474 m_v[digiId_Mc]=v;
475 m_p[digiId_Mc]=p;
476
477 int layer= layerId_Mc;
478 vec_slant.push_back(SlantId(layer));
479 m_slant[digiId_Mc]=SlantId(layer);
480 if (m_slant!=0)
481 {cout<<"layer: "<<layerId_Mc<<"wire: "<<wireId_Mc<<"x_east: "<<m_x_east[digiId_Mc]<<"y_east: "<<m_y_east[digiId_Mc]<<endl;}
482 }
483 }
484 }
485
486 t_nHit_Mc=digiId_Mc;
487 m_nHit_Mc=digiId_Mc;
488 if(m_data==0) {m_nHit=digiId_Mc;}
489
490 //m_nHitAxial=nHitAxial_Mc;
491 // m_nFitFailure=0;
492
493 // for(int i =0;i<t_nHit_Mc;i++){
494 // if(t_maxP<vec_p[i]) {t_maxP=vec_p[i];}
495 // }
496 //cout<<"t_maxP:"<<t_maxP<<" "<<endl;
497
498 // for(int i =0;i<t_nHit_Mc;i++){
499 // if(t_minP>vec_p[i]) {t_minP=vec_p[i];}
500 // }
501 //cout<<"t_minP:"<<t_minP<<" "<<endl;
502
503
504 // retrieve Mdc digi vector form RawDataProviderSvc
505 vector<int> vec_hitSignal;
506 vector<int> vec_track_index;
507 int track_index_max=0;
508 uint32_t getDigiFlag = 0;
509 if(m_dropHot || m_combineTracking)getDigiFlag |= MdcRawDataProvider::b_dropHot;
510 if(m_keepBadTdc) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
511 if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
512 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
513 // MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec();
514 cout<<"MdcDigiVec: "<<mdcDigiVec.size()<<endl;
515 MdcDigiVec::iterator iter1 = mdcDigiVec.begin();
516 int t_nHit;
517 int digiId = 0;
518 int nHitAxial = 0;
519 int nHitLayer[43];
520 int nHitSignal=0;
521 int nHitAxialSignal=0;
522
523 for (;iter1 != mdcDigiVec.end(); iter1++, digiId++) {
524 int track_index=(*iter1)->getTrackIndex();
525// cout<<"mc track_index: "<<track_index<<endl;
526 if(track_index>=1000) track_index-=1000;
527 vec_track_index.push_back(track_index);
528 if(track_index>=0){
529 nHitSignal++;
530 m_hitSignal[digiId]=1;
531 }
532 //m_hitCol[digiId]=digiId;
533 Identifier mdcId= (*iter1)->identify();
534 int layerId = MdcID::layer(mdcId);
535 int wireId = MdcID::wire(mdcId);
536
537
538 //cout<<"layer: "<<layerId<<" "<<"wire: "<<wireId<<"hitSignal: "<<vec_hitSignal[digiId]<<endl;
539
540 if ((layerId>=8&&layerId<=19)||(layerId>=36)) {
541 nHitAxial++;}
542 if (((layerId>=8&&layerId<=19)||(layerId>=36))&&track_index>=0) {
543 nHitAxialSignal++;}
544
545 vec_layer.push_back(layerId);
546 vec_wire.push_back(wireId);
547
548 nHitLayer[layerId]++;
549 m_layerNhit[layerId]=nHitLayer[layerId];
550
551 const MdcGeoWire *geowir =m_mdcGeomSvc->Wire(layerId,wireId);
552 HepPoint3D eastP = geowir->Backward()/10.0;//mm 2 cm
553 HepPoint3D westP = geowir->Forward() /10.0;//mm 2 cm
554
555 vec_x_east.push_back(eastP.x());
556 vec_y_east.push_back(eastP.y());
557 vec_z_east.push_back(eastP.z());
558 vec_x_west.push_back(westP.x());
559 vec_y_west.push_back(westP.y());
560 vec_z_west.push_back(westP.z());
561
562 m_x_east[digiId]=eastP.x();
563 m_y_east[digiId]=eastP.y();
564 m_z_east[digiId]=eastP.z();
565 m_x_west[digiId]=westP.x();
566 m_y_west[digiId]=westP.y();
567 m_z_west[digiId]=westP.z();
568 if(m_debug>0) {cout<<"event: "<<t_eventNum<<" "<<"layer: "<<layerId<<" "<<"wireId: "<<wireId<<" "<<"zeast: "<<vec_z_east.at(digiId)<<" "<<"zwest: "<<vec_z_west.at(digiId)<<" "<<endl;}
569 //conformal trans:
570 double x =(eastP.x()+westP.x())/2.;
571 double y =(eastP.y()+westP.y())/2.;
572 vec_x.push_back((vec_x_east[digiId]+vec_x_west[digiId])/2);
573 vec_y.push_back((vec_y_east[digiId]+vec_y_west[digiId])/2);
574 m_x[digiId]=(vec_x_east[digiId]+vec_x_west[digiId])/2;
575 m_y[digiId]=(vec_y_east[digiId]+vec_y_west[digiId])/2;
576 double u=CFtrans(x,y);
577 double v=CFtrans(y,x);
578 // cout<< "u: "<<u<<"v: "<<v<<"digiId: "<<digiId<<endl;
579 //cout<<"layer: "<<layerId<<" "<<"wireId: "<<wireId<<" "<<"x: "<<x<<" "<<"y: "<<y<<" "<<endl;
580 vec_p.push_back(sqrt(u*u+v*v));
581 vec_u.push_back(u);
582 vec_v.push_back(v);
583 m_u[digiId]=u;
584 m_v[digiId]=v;
585 m_p[digiId]=sqrt(u*u+v*v);
586
587 m_layer[digiId]=geowir->Layer();
588 m_cell[digiId]=geowir->Cell();
589 int layer= layerId;
590 vec_slant.push_back(SlantId(layer));
591 m_slant[digiId]=SlantId(layer);
592 //cout<<"layer: "<<layer<<" "<<"slant: "<<m_slant[digiId]<<endl;
593 // if (m_slant!=0)
594 // {cout<<"layer: "<<layerId<<"wire: "<<wireId<<"x_east: "<<m_x_east[digiId]<<"y_east: "<<m_y_east[digiId]<<endl;}
595 }
596 for(int i=0;i<digiId;i++){
597 if(track_index_max<=vec_track_index[i]) track_index_max=vec_track_index[i];
598 }
599 track_index_max++;
600 m_trackNum_Mc=track_index_max;
601// cout<<"mc truth has :"<<track_index_max<<" mc track."<<endl;
602 vector< vector <int> > mc_track_hit(track_index_max,vector<int>() );
603 for(int i=0;i<track_index_max;i++){
604 for(int j=0;j<digiId;j++){
605 if(vec_track_index[j]==i) mc_track_hit[i].push_back(j);
606 }
607 }
608 //composory set mc_track number
609 track_index_max=m_trackNum_Mc_set;
610 //debug mc track hit
611 //for(int i=0;i<track_index_max;i++){
612 // for(int j=0;j<mc_track_hit[i].size();j++){
613 // int hitId=mc_track_hit[i][j];
614 // cout<<"track: "<<i<<" has hit: "<<mc_track_hit[i][j]<<"("<<vec_layer[hitId]<<","<<vec_wire[hitId]<<")"<<endl;
615 // }
616 //}
617
618 //t_nHit=digiId;
619 m_nHit=digiId;
620 m_nHitAxial=nHitAxial;
621 //m_nFitFailure=0;
622
623 m_nHitSignal=nHitSignal;
624 m_nHitAxialSignal=nHitAxialSignal;
625
626 cout<<"hit number: "<<digiId<<endl;
627
628
629 // get the max&min amplitude of the curve
630 for(int i =0;i<m_nHit;i++){
631 if(t_maxP<vec_p[i]) {t_maxP=vec_p[i];}
632 }
633
634 for(int i =0;i<m_nHit;i++){
635 if(t_minP>vec_p[i]) {t_minP=vec_p[i];}
636 }
637 t_maxP=t_maxP+0.01;
638
639
640 //double peakArea[100];
641 //m_npeak=100;
642 //for( int bin_i=0;bin_i<10;bin_i++){
643 // for(int bin_j=0;bin_j<10;bin_j++){
644 // int binx=m_binx+100*bin_i;
645 // int biny=m_biny+100*bin_j;
646 int nhit;
647 if(m_data==0) {
648 nhit=m_nHit_Mc;}
649 else{
650 nhit=m_nHit;}
651
652 binX=m_binx;
653 binY=m_biny;
654 double binwidth=M_PI/binX;
655 double binhigh=2*t_maxP/binY;
656 TH2D *h1=new TH2D("h1","h1",binX,0,M_PI,binY,-t_maxP,t_maxP);
657
658 //method 0
659 //2-point combine
660 vector<double> vec_rho;
661 vector<double> vec_theta;
662 vector< vector<int> > vec_hitNum(2,vector<int>()); //store 2 hits in a cross point
663 int numCross=(int)(nhit*(nhit-1)/2);
664 if(m_method==0){
665 RhoTheta(numCross,nhit,vec_u,vec_v,vec_rho,vec_theta,vec_hitNum); //calculate cross point
666 FillRhoTheta(h1,vec_theta,vec_rho,numCross); //fill histgram method by rhotheta
667 }
668 // method 1
669 vector< vector< vector<int> > > vec_selectNum(binX,vector< vector<int> >(binY,vector<int>() ) );
670 if(m_method==1){
671 FillHist(h1,vec_u,vec_v,nhit,vec_selectNum);
672 }
673
674
675 vector<bool> vec_truthHit(nhit,false);
676 vector< vector <int> > countij(binX,vector <int> (binY,0) );
677 // vector< vector < vector<int> > > vec_selectNum(binX,vector < vector<int> >(binY,vector<int>() ));
678 vector< vector< vector<int> > > vec_selectHit(binX,vector< vector<int> >(binY,vector<int>() ) );
679 if(m_method==2){
680 FillCells(h1,nhit,binX,binY,vec_u,vec_v,vec_layer,vec_wire,countij,vec_selectNum,vec_selectHit);
681 }
682
683 //cout h1
684 //for(int i=0;i<binX;i++){
685 // for(int j=0;j<binY;j++){
686 // int count =h1->GetBinContent(i+1,j+1);
687 //// for(int k=0;k<count;k++){
688 // cout<<"("<<i<<","<<j<<")"<<": "<<count<<endl;
689 //// }
690 // }
691 //}
692
693
694 //int nHitSelect=0;
695 //int nHitSignal_select=0;
696 //int nHitAxialSignal_select=0;
697 //find peak
698 int max_count=0;
699 int max_binx=1;
700 int max_biny=1;
701 for(int i=0;i<binX;i++){
702 for(int j=0;j<binY;j++){
703 int count=h1->GetBinContent(i+1,j+1);
704 if(max_count<count) {
705 max_count=count;
706 max_binx=i+1;
707 max_biny=j+1;
708 }
709 }
710 }
711 // cout<<"("<<max_binx<<","<<max_biny<<","<<max_count<<")"<<"numCross: "<<numCross<<endl;
712 // h1->Draw("lego2");
713
714 //int peak_cellNum=1; // to determine how many cell should be embraced in one orientation
715 //int count_sum[1000];
716 //for(int i=0;i<1000;i++){
717 // count_sum[i]=0;
718 // for(int alphax=max_binx-i;alphax<=max_binx+i;alphax++){
719 // for(int rhox=max_biny-i;rhox<=max_biny+i;rhox++){
720 // int ax;
721 // if (alphax<0) { ax=alphax+binX; }
722 // else if (alphax>=binX) { ax=alphax-binX; }
723 // else { ax=alphax; }
724 // count_sum[i]+=h1->GetBinContent(ax,rhox);
725 // }
726 // }
727 // if(count_sum[i]>m_fpro*numCross){
728 // peak_cellNum=i;
729 // cout<<"i: "<<i<<endl;
730 // cout<<"count_sum: "<<count_sum[i]<<endl;
731 // cout<<"pro: "<<(double)count_sum[i]/(double)numCross<<endl;
732 // break;
733 // }
734 //}
735
736
737 /*
738 //return hit number
739 vector<bool> vec_hitSelect(m_nHit,false);
740 int bin_left=max_binx-peak_cellNum;
741 int bin_right=max_binx+peak_cellNum;
742 int bin_up=max_biny+peak_cellNum;
743 int bin_down=max_biny-peak_cellNum;
744 double area_left=(bin_left-1)*binwidth;
745 double area_right=(bin_right)*binwidth;
746 double area_down=(bin_down-1)*binhigh-t_maxP;
747 double area_up=(bin_up)*binhigh-t_maxP;
748 // cout<<"("<<area_left<<","<<area_right<<","<<area_down<<","<<area_up<<")"<<endl;
749 double peak_width=(2*peak_cellNum+1)*binwidth;
750 double peak_high=(2*peak_cellNum+1)*binhigh;
751 cout<<"the area of the peak is: "<<peak_width*peak_high<<endl;
752 //test if the hit is in this area
753 for(int i=0;i<numCross;i++){
754 if(vec_theta[i]>=area_left&&vec_theta[i]<area_right&&vec_rho[i]>=area_down&&vec_rho[i]<area_up) {
755 int hitNum_1=vec_hitNum[0][i];
756 int hitNum_2=vec_hitNum[1][i];
757 vec_hitSelect[hitNum_1]=true;
758 vec_hitSelect[hitNum_2]=true;
759 }
760 }
761 int selectHitNum=0;
762 for(int i=0;i<m_nHit;i++){
763 if(vec_hitSelect[i]==1) {selectHitNum++;}
764 // cout<<"if hit is select: "<<i<<": "<<vec_hitSelect[i]<<endl;
765 }
766 cout<<"how many Hit have been selected: "<<selectHitNum<<" "<<(double)selectHitNum/(double)m_nHit;
767 m_areaSelectHit=(double)selectHitNum/(double)m_nHit;
768
769
770 // // the smallest cell has been found
771 // //cout<<"width: "<<peak_width<<" "<<"high: "<<peak_high<<endl;
772 // int column=(int)(bin_i*10+bin_j);
773 // peakArea[column]=peak_width*peak_high;
774 // m_peakWidth[column]=peak_width;
775 // m_peakHigh[column]=peak_high;
776 /// m_peakArea[column]=peak_width*peak_high;
777 // int columnx_split=(int)(column/10)*100+100;
778 // int columny_split=(int)(column%10)*100+100;
779 // //cout<<"column: "<<column<<" "<<"("<<columnx_split<<","<<columny_split<<")"<<peakArea[column]<<endl;
780 delete h1;
781 // }
782 // }
783 // double area_least=peakArea[0];
784 // for(int i=0;i<100;i++){
785 // if(area_least>peakArea[i]){
786 // area_least=peakArea[i];
787 // }
788 // }
789 // int num_area_least=0;
790 // for(int i=0;i<100;i++){
791 // if(area_least==peakArea[i]){
792 // num_area_least=i;
793 // }
794 // }
795 // int binx_split=(int)(num_area_least/10)*100+100;
796 // int biny_split=(int)(num_area_least%10)*100+100;
797 // m_areaLeast=area_least;
798 // m_areaLeastNum=num_area_least;
799 //cout<<"the max peakArea: "<<num_area_least<<"("<<binx_split<<","<<biny_split<<")"<<"peakarea is : "<<area_least<<endl;
800 */
801
802
803
804
805 // select hit
806 int count_use=0;
807 double sum=0;
808 double sum2=0;
809 for(int i=0;i<binX;i++){
810 for(int j=0;j<binY;j++){
811 int count=h1->GetBinContent(i+1,j+1);
812 // cout<<"("<<i<<","<<j<<")"<<" "<<"count: "<<count<<" ";
813 if(count!=0){
814 count_use++;
815 sum+=count;
816 sum2+=count*count;
817 }
818 }
819 }
820 double f_n=m_ndev;
821 cout<<"count_use: "<<count_use<<endl;
822 double aver=sum/count_use;
823 double dev=sqrt(sum2/count_use-aver*aver);
824 int min_counts=(int)(aver+f_n*dev);
825 cout<<"aver: "<<aver<<" "<<"dev: "<<dev<<" min: "<<min_counts<<endl;
826
827
828 vector< vector <int> > hough_trans_CS(binX,vector <int> (binY,0) );
829 vector< vector <int> > hough_trans_CS_peak(binX,vector <int> (binY,0) );
830 // int hough_trans_CS[binX][binY];
831 // int hough_trans_CS_peak[binX][binY];
832 for(int i=0;i<binX;i++){
833 for(int j=0;j<binY;j++){
834 hough_trans_CS[i][j]=h1->GetBinContent(i+1,j+1);
835 }
836 }
837
838 int peak_num=0;
839 for (int r=0; r<binY; r++) {
840 for (int a=0; a<binX; a++) {
841 long max_around=0;
842 for (int ar=a-1; ar<=a+1; ar++) {
843 for (int rx=r-1; rx<=r+1; rx++) {
844 int ax;
845 if (ar<0) { ax=ar+binX; }
846 else if (ar>=binX) { ax=ar-binX; }
847 else { ax=ar; }
848 if ( (ax!=a || rx!=r) && rx>=0 && rx<binY) {
849 if (hough_trans_CS[ax][rx]>max_around) { max_around=hough_trans_CS[ax][rx]; }
850 }
851 }
852 }
853 if (hough_trans_CS[a][r]<=min_counts) { hough_trans_CS_peak[a][r]=0; }
854 else if (hough_trans_CS[a][r]<max_around) { hough_trans_CS_peak[a][r]=0; }
855 else if (hough_trans_CS[a][r]==max_around) { hough_trans_CS_peak[a][r]=1; }
856 else if (hough_trans_CS[a][r]>max_around) { hough_trans_CS_peak[a][r]=2; }
857 if(hough_trans_CS_peak[a][r]!=0) {
858 // cout<<" peak: "<<"("<<a<<","<<r<<")"<<":"<<hough_trans_CS_peak[a][r]<<" "<<hough_trans_CS[a][r]<<endl;
859 peak_num++;
860 }
861 }
862 }
863 //cout<<"peak_num: "<<peak_num<<endl;
864
865 // //find double-point-peaks in parameter space
866 int peak_num_2=0;
867 for (int r=0; r<binY; r++) {
868 for (int a=0; a<binX; a++) {
869 if (hough_trans_CS_peak[a][r]==1) {
870 hough_trans_CS_peak[a][r]=2;
871 for (int ar=a-1; ar<=a+1; ar++) {
872 for (int rx=r-1; rx<=r+1; rx++) {
873 int ax;
874 if (ar<0) { ax=ar+binX; }
875 else if (ar>=binX) { ax=ar-binX; }
876 else { ax=ar; }
877 if ( (ax!=a || rx!=r) && rx>=0 && rx<binY) {
878 if (hough_trans_CS_peak[ax][rx]==1) { hough_trans_CS_peak[ax][rx]=0; }
879 }
880 }
881 }
882 }
883 if(hough_trans_CS_peak[a][r]!=0) {
884 peak_num_2++;
885 //cout<<" peak: "<<"("<<a<<","<<r<<")"<<":"<<hough_trans_CS_peak[a][r]<<" "<<hough_trans_CS[a][r]<<endl;
886 // for(int i=0;i<hough_trans_CS[a][r];i++){
887 // int hitNum=selectNum[a][r][i];
888 // cout<<"hit: "<<hitNum<<"("<<layer[hitNum]<<","<<cell[hitNum]<<")"<<endl;
889 // }
890 }
891 }
892 }
893 //cout<<"peak_mum_2: "<<peak_num_2<<endl;
894
895 //fill the histogram again
896 TH2D *h2 = new TH2D("h2","test2",binX,0,3.14,binY,-t_maxP,t_maxP);
897 for(int i=0;i<binX;i++){
898 for(int j=0;j<binY;j++){
899 if(hough_trans_CS_peak[i][j]==2){
900 h2->SetBinContent(i+1,j+1,hough_trans_CS[i][j]);
901 }
902 }
903 }
904
905 vector<int> peakList1[3];
906 for(int n=0;n<peak_num_2;n++){
907 for (int r=0; r<binY; r++) {
908 for (int a=0; a<binX; a++) {
909 if (hough_trans_CS_peak[a][r]==2) {
910 peakList1[0].push_back(a+1);
911 peakList1[1].push_back(r+1);
912 peakList1[2].push_back(hough_trans_CS[a][r]);
913 }
914 }
915 }
916 }
917 cout<<"finish peak?"<<endl;
918 if(m_method==0) cout<<"numCross: "<<numCross<<endl;
919 // //sort peaks by height
920 int npeak1=peak_num_2;
921 int n_max;
922 for (int na=0; na<npeak1-1; na++) {
923 n_max=na;
924 for (int nb=na+1; nb<npeak1; nb++) {
925 if (peakList1[2][n_max]<peakList1[2][nb]) { n_max=nb; }
926 }
927 if (n_max!=na) { // swap
928 float swap[3];
929 for (int i=0; i<3; i++) {
930 swap[i]=peakList1[i][na];
931 peakList1[i][na]=peakList1[i][n_max];
932 peakList1[i][n_max]=swap[i];
933 }
934 }
935 }
936 cout<<"npeak1: "<<npeak1<<endl;
937 for(int n=0;n<npeak1;n++){
938 int bintheta=peakList1[0][n];
939 int binrho=peakList1[1][n];
940 int count=peakList1[2][n];
941 for(int i=0;i<count;i++){
942 //cout<<"("<<bintheta<<","<<binrho<<","<<count<<")"<<endl;
943 //cout<<"n"<<n<<": "<<"("<<bintheta<<","<<binrho<<","<<count<<"):"<<vec_selectNum[bintheta-1][binrho-1][i]<<endl;
944 }
945 //cout<<"peakList1: "<<n<<" "<<"bina: "<<peakList1[0][n]<<" "<<"binr: "<<peakList1[1][n]<<" "<<"houghpeak: "<<peakList1[2][n]<<endl;
946 }
947
948 typedef std::map<int, int> M2;
949 typedef std::map<int, M2> M1;
950 M2 peak_combine_num;
951 M1 hit_combine;
952
953 int peak_track=0;
954
955 if(m_method==2){
956 //combine aroud
957 int m=0;
958 int n=0;
959 int addnum=999;
960 vector<bool> peaktrack(npeak1,true);
961 while(addnum!=0)
962 {
963 addnum=0;
964 double peak_cellNum[43];
965 int peakX[43];
966 int peakY[43];
967 double bin_left[43];
968 double bin_right[43];
969 double bin_up[43];
970 double bin_down[43];
971 double area_left[43];
972 double area_right[43];
973 //double area_mid[43];
974 double area_down[43];
975 double area_up[43];
976 // for(int i=0;i<npeak1;i++) { peaktrack[i]=true; }
977 for(int iLayer=0; iLayer<m_gm->nLayer(); iLayer++){
978 peak_cellNum[iLayer]=m_peakCellNum[iLayer];
979 peakX[iLayer]=peakList1[0][n]; //start from the highest peak
980 peakY[iLayer]=peakList1[1][n];
981 bin_left[iLayer]=peakX[iLayer]-peak_cellNum[iLayer];
982 bin_right[iLayer]=peakX[iLayer]+peak_cellNum[iLayer];
983 bin_up[iLayer]=peakY[iLayer]+peak_cellNum[iLayer];
984 bin_down[iLayer]=peakY[iLayer]-peak_cellNum[iLayer];
985 area_left[iLayer]=(bin_left[iLayer]-1)*binwidth;
986 area_right[iLayer]=(bin_right[iLayer])*binwidth;
987 // area_mid=(area_left+area_right)/2.;
988 area_down[iLayer]=(bin_down[iLayer]-1)*binhigh-t_maxP;
989 area_up[iLayer]=(bin_up[iLayer])*binhigh-t_maxP;
990 }
991 int count_peak=0;
992 for(int k=0;k<nhit;k++){
993 int layer=vec_layer[k];
994 double lineLeft=vec_u[k]*cos(area_left[layer])+vec_v.at(k)*sin(area_right[layer]);
995 double lineRight=vec_u[k]*cos(area_left[layer])+vec_v.at(k)*sin(area_right[layer]);
996 // double lineMid=vec_u[k]*cos(area_mid)+vec_v[k]*sin(area_mid);
997 if((lineLeft<area_up[layer])&&(lineLeft>area_down[layer])||(lineRight<area_up[layer])&&(lineRight>area_down[layer])||((lineLeft-area_up[layer])*(area_down[layer]-lineRight)>0)){
998 // hit_combine[m].push_back(k);
999 // if(m==0){
1000 // nHitSelect++;
1001 // if(vec_track_index[k]>=0){nHitSignal_select++;}
1002 // if(vec_track_index[k]>=0&&vec_slant[k]==0){nHitAxialSignal_select++;}
1003 // }
1004 hit_combine[m][count_peak]=k;
1005 //vec_truthHit[m][k]=1;
1006 count_peak++;
1007 }
1008 //else {vec_truthHit[m][k]=0;}
1009 peak_combine_num[m]=count_peak;
1010 //int sizepeak=hit_combine[0].size();
1011 //cout<<"peak1size: "<<sizepeak<<endl;
1012 }
1013 for(int i=n+1;i<npeak1;i++){
1014 if(peaktrack[i]==false) continue;
1015 int peaktheta=peakList1[0][i];
1016 int peakrho=peakList1[1][i];
1017 int peakhitNum=peakList1[2][i];
1018 int count_same_hit=0;
1019 for(int j=0;j<peakhitNum;j++){
1020 //for(int k=0;k<hit_combine[m].size();k++)
1021 for(int k=0;k<peak_combine_num[m];k++){
1022 if(vec_selectNum[peaktheta-1][peakrho-1][j]==hit_combine[m][k]){
1023 count_same_hit++;
1024 }
1025 }
1026 }
1027 double f_hit=m_hit_pro;
1028 if(count_same_hit>f_hit*peakhitNum){
1029 peaktrack[i]=false;
1030 }
1031 }
1032 for(int i=n+1;i<npeak1;i++){
1033 if(peaktrack[i]==true){
1034 addnum=i;
1035 break;
1036 }
1037 }
1038 if(addnum!=0) m++;
1039 cout<<"peak_m: "<<m+1<<endl;
1040 n=n+addnum;
1041 }
1042
1043
1044 for(int i=0;i<npeak1;i++){
1045 cout<<"track"<<i<<": "<<"("<<peakList1[0][i]<<","<<peakList1[1][i]<<","<<peakList1[2][i]<<")"<<" "<<"truth: "<<peaktrack[i]<<endl;
1046 if(peaktrack[i]==true){
1047 peak_track++;
1048 }
1049 }
1050 for( int i=0;i<peak_track;i++){
1051 for(int j =0;j<peak_combine_num[i];j++){
1052 int hit_number=hit_combine[i][j];
1053 cout<<" peak "<<i<<" has select hits: "<<vec_layer[hit_number]<<" "<<vec_wire[hit_number]<<endl;
1054 }
1055 cout<<"has collect :"<<peak_combine_num[i]<<" hits "<<endl;
1056 }
1057 cout<<"event"<<t_eventNum<<": "<<"has found: "<<peak_track<<" track."<<endl;
1058 m_trackNum=peak_track;
1059 //m_nHitSelect=nHitSelect;
1060 //m_nHitSignal_select=nHitSignal_select;
1061 //m_nHitAxialSignal_select=nHitAxialSignal_select;
1062
1063 //according mc and track
1064 vector< vector<int> > rec_mc_num(peak_track,vector<int>(track_index_max,0) );
1065 vector< vector<int> > rec_mc_num_axial(peak_track,vector<int>(track_index_max,0) );
1066 if(track_index_max!=peak_track) cout<<"not match track number !"<<endl;
1067 for(int mc_track_num=0;mc_track_num<track_index_max;mc_track_num++){
1068 for(int i=0;i<peak_track;i++){
1069 for(int j=0;j<peak_combine_num[i];j++){
1070 for(int k=0;k<mc_track_hit[mc_track_num].size();k++){
1071 if(hit_combine[i][j]==mc_track_hit[mc_track_num][k]){
1072 rec_mc_num[i][mc_track_num]++;
1073 int hit_num=mc_track_hit[mc_track_num][k];
1074 if(vec_slant[hit_num]==0) rec_mc_num_axial[i][mc_track_num]++;
1075 }
1076 }
1077 }
1078 }
1079 }
1080 vector<int> rec_mc(peak_track,999);
1081 for(int i=0;i<peak_track;i++){
1082 for(int mc_track_num=0;mc_track_num<track_index_max;mc_track_num++){
1083 if(rec_mc_num[i][mc_track_num]>0.5*peak_combine_num[i]) rec_mc[i]=mc_track_num;
1084 cout<<"trackrec: "<<i<<" trackmc: "<<mc_track_num<<" "<<rec_mc_num[i][mc_track_num]<<" "<<peak_combine_num[i]<<endl;
1085 }
1086 }
1087 for(int i=0;i<peak_track;i++){
1088 cout<<"rec :"<<i<<"belong to mc track: "<<rec_mc[i]<<endl;
1089 }
1090 for(int i=0;i<peak_track;i++){
1091 int mc_track_num=rec_mc[i];
1092 if(mc_track_num!=999) {
1093 cout<<"debug mc_track_num: "<<mc_track_num<<endl;
1094 m_nHitSignal_select[i]=rec_mc_num[i][mc_track_num]; //????????? mc_track_num=999?
1095 m_nHitAxialSignal_select[i]=rec_mc_num_axial[i][mc_track_num]; //????????? mc_track_num=999?
1096 cout<<"m_nHitSignal_select: "<<m_nHitSignal_select[i]<<endl;
1097 }
1098 }
1099
1100 //cout<<"peak1: "<<endl;
1101 //for(int i=0;i<peak_track;i++){
1102 // for(int j=0;j<peak_combine_num[i];j++){
1103 // int hitNumtemp=hit_combine[i][j];
1104 // vec_truthHit[hitNumtemp]=i;
1105 // cout<<" "<<"("<<vec_layer[hitNumtemp]<<","<<vec_wire[hitNumtemp]<<")"<<": "<<vec_track_index[hitNumtemp]<<endl;
1106 // }
1107 //}
1108
1109 // if(peak_track==2){
1110 // cout<<"peak2: "<<endl;
1111 // for(int i=0;i<hit_combine[1].size();i++){
1112 // int hitNumtemp=hit_combine[1][i];
1113 // vec_truthHit[hitNumtemp]=1;
1114 // cout<<" "<<"("<<vec_layer[hitNumtemp]<<","<<vec_wire[hitNumtemp]<<")"<<": "<<vec_track_index[hitNumtemp]<<endl;
1115 // }
1116 // }
1117
1118
1119 }
1120
1121
1122 //yzhang 2015-03-03 delete
1123// //return hit number
1124// if(m_method==0){
1125// vector<int> numCross_select;
1126// vector< vector<bool> > vec_hitSelect(npeak1,vector<bool> (nhit,false) );
1127// vector<int> vec_selectHitNum;
1128// vector<int> vec_selectHitNum_signal;
1129// for(int peakNum=0;peakNum<npeak1;peakNum++){
1130// int peak_cellNum=m_peakCellNum;
1131// // vector<bool> vec_hitSelect(nhit,false);
1132// int peakX=peakList1[0][peakNum];
1133// int peakY=peakList1[1][peakNum];
1134// int bin_left=peakX-peak_cellNum;
1135// int bin_right=peakX+peak_cellNum;
1136// int bin_up=peakY+peak_cellNum;
1137// int bin_down=peakY-peak_cellNum;
1138// double area_left=(bin_left-1)*binwidth;
1139// double area_right=(bin_right)*binwidth;
1140// double area_down=(bin_down-1)*binhigh-t_maxP;
1141// double area_up=(bin_up)*binhigh-t_maxP;
1142// // cout<<"("<<area_left<<","<<area_right<<","<<area_down<<","<<area_up<<")"<<endl;
1143// double peak_width=(2*peak_cellNum+1)*binwidth;
1144// double peak_high=(2*peak_cellNum+1)*binhigh;
1145// //cout<<"the area of the peak is: "<<peak_width*peak_high<<endl;
1146// //test if the hit is in this area
1147// int numCross_select_temp=0;
1148// for(int i=0;i<numCross;i++){
1149// if(vec_theta[i]>=area_left&&vec_theta[i]<area_right&&vec_rho[i]>=area_down&&vec_rho[i]<area_up) {
1150// numCross_select_temp++;
1151// int hitNum_1=vec_hitNum[0][i];
1152// int hitNum_2=vec_hitNum[1][i];
1153// vec_hitSelect[peakNum][hitNum_1]=true;
1154// vec_hitSelect[peakNum][hitNum_2]=true;
1155// }
1156// }
1157// numCross_select.push_back(numCross_select_temp);
1158// int selectHitNum=0;
1159// int selectHitNum_signal=0;
1160// for(int i=0;i<nhit;i++){
1161// if(vec_hitSelect[peakNum][i]==true) {
1162// selectHitNum++;
1163// if(vec_hitSignal[i]==1){
1164// selectHitNum_signal++;
1165// }
1166// // cout<<"if hit is select: "<<i<<": "<<vec_hitSelect[i]<<endl;
1167// }
1168// }
1169// vec_selectHitNum.push_back(selectHitNum);
1170// vec_selectHitNum_signal.push_back(selectHitNum_signal);
1171// // cout<<"hit: "<<selectHitNum<<" "<<(double)selectHitNum/(double)nhit<<endl;
1172// //cout<<"true hit : "<<selectHitNum_signal<<" "<<(double)selectHitNum_signal/(double)nhit<<endl;
1173// m_areaSelectHit=(double)selectHitNum/(double)nhit;
1174// m_areaSelectHit_signal=(double)selectHitNum_signal/(double)nhit;
1175// }
1176//
1177// for(int i=0;i<npeak1;i++){
1178// int selectHitNum=0;
1179// for(int j=0;j<nhit;j++){
1180// if(vec_hitSelect[i][j]==true){
1181// selectHitNum++;
1182// // cout<<"event"<<t_eventNum<<": "<<"peak"<<i<<": "<<"hit: "<<"("<<vec_layer[j]<<","<<vec_wire[j]<<")"<<endl; //cout picked hit
1183// }
1184// }
1185// cout<<"peak"<<i<<": "<<selectHitNum<<" hits"<<endl;
1186// }
1187// // for(int i=0;i<3;i++){
1188// // int bina=peakList1[0][i];
1189// // int binr=peakList1[1][i];
1190// // int hit_peak=peakList1[2][i];
1191// // for(int j=0;j<hit_peak;j++){
1192// // int hitList=vec_selectHit[bina][binr][j];
1193// // cout<<"peaknum: "<<i<<"("<<bina<<","<<binr<<")"<<" "<<peakList1[2][i]<<" hit : "<<hitList<<endl;
1194// // // cout<<"peaknum: "<<i<<"("<<bina<<","<<binr<<")"<<" "<<peakList1[2][i]<<" hit : "<<vec_selectHit[bina][binr][j]<<endl;
1195// // }
1196// // }
1197// //cout<<npeak1<<endl;
1198//
1199// for(int n=0;n<npeak1;n++){
1200// cout<<"peakList1: "<<n<<" "<<"bina: "<<peakList1[0][n]<<" "<<"binr: "<<peakList1[1][n]<<" "<<"houghpeak: "<<peakList1[2][n]<<" "<<"numCross_select:"<<numCross_select[n]<<" "<<" cross: "<<(double)numCross_select[n]/(double)numCross<<" "<<"hit_selectNum: "<<vec_selectHitNum_signal[n]<<" "<<100*(double)vec_selectHitNum_signal[n]/(double)m_nHitSignal<<"% "<<"noise: "<<(double)(vec_selectHitNum[n]-vec_selectHitNum_signal[n])<<endl;
1201// }
1202// }
1203
1204 /*
1205 //combine the neighbor cell
1206 //vector<int> peak_hitList[npeak1];
1207 vector< vector <int> > peak_hitList(npeak1,vector <int> () );
1208 for(int i=0;i<npeak1;i++){
1209 int bina=peakList1[0][i];
1210 int binr=peakList1[1][i];
1211 int hit_peak=peakList1[2][i];
1212 for(int j=0;j<hit_peak;j++){
1213 peak_hitList[i].push_back(vec_selectHit[bina][binr][j]);
1214// cout<<"j: "<<j<<" "<<vec_selectHit[bina][binr][j]<<endl;
1215}
1216
1217for(int a=bina-1;a<=bina+1;a++){
1218for(int rx=binr-1;rx<=binr+1;rx++){
1219int ax;
1220if (a<0) { ax=a+binX; }
1221else if (a>=binX) { ax=a-binX; }
1222else { ax=a; }
1223if ( (ax!=bina || rx!=binr) && rx>=0 && rx<binY) {
1224int hit_around=count[ax][rx];
1225for(int k=0;k<hit_around;k++){
1226bool combine=1;
1227for(int m=0;m<hit_peak;m++){
1228if(vec_selectHit[ax][rx][k]==peak_hitList[i][m]){
1229// cout<<i<<" "<<" "<<ax<<" "<<rx<<" "<<k<<" "<<vec_selectHit[ax][rx][k]<<endl;
1230combine=0;
1231continue;
1232}
1233}
1234if(combine==1) {
1235peak_hitList[i].push_back(vec_selectHit[ax][rx][k]);
1236peakList1[2][i]++;
1237hit_peak++;
1238// cout<<" there is a combine"<<" "<<vec_selectHit[ax][rx][k]<<endl;
1239}
1240}
1241}
1242}
1243}
1244}
1245
1246
1247//for(int n=0;n<npeak1;n++){
1248// cout<<"peakList2: "<<n<<" "<<"a: "<<peakList1[0][n]<<" "<<"r: "<<peakList1[1][n]<<" "<<"houghpeak: "<<peakList1[2][n]<<endl;
1249//}
1250
1251
1252for(int i=0;i<1;i++){
1253int bina=peakList1[0][i];
1254int binr=peakList1[1][i];
1255int hit_peak=peakList1[2][i];
1256for(int j=0;j<hit_peak;j++){
1257int hitList=peak_hitList[i][j];
1258// cout<<"peaknum: "<<i<<"("<<bina<<","<<binr<<")"<<" "<<peakList1[2][i]<<" hit : "<<hitList<<endl;
1259// cout<<"peaknum: "<<i<<"("<<bina<<","<<binr<<")"<<" "<<peakList1[2][i]<<" hit : "<<vec_selectHit[bina][binr][j]<<endl;
1260}
1261}
1262
1263//sort another time
1264int n_max_2;
1265for (int na=0; na<npeak1-1; na++) {
1266n_max_2=na;
1267for (int nb=na+1; nb<npeak1; nb++) {
1268if (peakList1[2][n_max_2]<peakList1[2][nb]) { n_max_2=nb; }
1269}
1270if (n_max_2!=na) { // swap
1271float swap[3];
1272for (int i=0; i<3; i++) {
1273swap[i]=peakList1[i][na];
1274peakList1[i][na]=peakList1[i][n_max_2];
1275peakList1[i][n_max_2]=swap[i];
1276}
1277}
1278}
1279*/
1280
1281// for(int n=0;n<npeak1;n++){
1282// cout<<"peakList3: "<<n<<" "<<"a: "<<peakList1[0][n]<<" "<<"r: "<<peakList1[1][n]<<" "<<"houghpeak: "<<peakList1[2][n]<<endl;
1283// }
1284
1285// int max_count=-999;
1286// for(int i=0;i<binX;i++){
1287// for(int j=0;j<binY;j++){
1288// if(max_count<count[i][j]) {max_count=count[i][j];}
1289// }
1290// }
1291
1292// m_maxCount=max_count;
1293//if(m_debug>0) cout<<"maxcount: "<<max_count<<endl;
1294
1295// for(int i=0;i<binX;i++){
1296// for(int j=0;j<binY;j++){
1297// // cout<<"("<<i<<","<<j<<")"<<" "<<"count: "<<count[i][j]<<" "<<endl;
1298// for(int l=0,m=vec_selectNum[i][j].size();l<m;++l){
1299// if(count[i][j]==max_count){
1300// //int layerId=vec_layer.at(vec_selectNum[i][j].at(l));
1301// //if ((layerId>=8&&layerId<=19)||(layerId>=36)) {
1302// // nHitAxialSelect++;
1303// //}
1304// vec_truthHit[vec_selectNum[i][j][l]]=true;
1305//
1306// }
1307// }
1308// }
1309//}
1310// for(int i=0;i<n;i++){
1311// cout<<"hit:("<<vec_layer[i]<<","<<vec_wire[i]<<")"<<"is"<<vec_truthHit[i]<<endl;
1312// }
1313// cout<<"selected true hits is: "<<endl;
1314
1315//for(int i=0;i<n;i++){
1316// if(vec_truthHit[i]==true){
1317// cout<<"hit:("<<vec_layer[i]<<","<<vec_wire[i]<<")"<<"is"<<vec_truthHit[i]<<" "<<"hitSignal:"<<vec_hitSignal[i]<<endl;
1318// }
1319//}
1320
1321
1322//int nHitAxialSelect=0;
1323//int nHitSelect=0;
1324//int nHitSignal_select=0;
1325//int nHitAxialSignal_select=0;
1326//for(int i=0;i<n;i++){
1327// if(vec_truthHit[i]==false) continue;
1328// nHitSelect++;
1329// if(vec_hitSignal[i]==1) {nHitSignal_select++;}
1330// int layerId=vec_layer[i];
1331// if ((layerId>=8&&layerId<=19)||(layerId>=36)) {
1332// nHitAxialSelect++;
1333// if(vec_hitSignal[i]==1) {nHitAxialSignal_select++;}
1334// }
1335//}
1336//
1337//m_nHitAxialSelect=nHitAxialSelect;
1338//m_nHitSelect=nHitSelect;
1339//
1340//m_nHitSignal_select=nHitSignal_select;
1341//m_nHitAxialSignal_select=nHitAxialSignal_select;
1342
1343// backtransform all track-
1344
1345double d0=-999.;
1346double phi0=-999.;
1347double omega=-999.;
1348double z0=0;
1349double tanl=0;
1350
1351for(track_fit=0;track_fit<peak_track;track_fit++){
1352 for(int i=0;i<nhit;i++){
1353 vec_truthHit[i]=false;
1354 }
1355 cout<<"track: "<<track_fit<<" has select: "<<peak_combine_num[track_fit]<<" hit ."<<endl;
1356 for(int i=0;i<peak_combine_num[track_fit];i++){
1357 int hitNum=hit_combine[track_fit][i];
1358 vec_truthHit[hitNum]=true;
1359 }
1360
1361 int nHitAxialSelect_temp=10;
1362 int leastSquare=LeastSquare(nhit,nHitAxialSelect_temp,vec_x,vec_y,vec_slant,vec_layer,vec_truthHit,d0,phi0,omega);
1363
1364 if(leastSquare==0){
1365
1366 TrkExchangePar tt(d0,phi0,omega,z0,tanl);
1367 TrkCircleMaker circlefactory;
1368 float chisum =0.;
1369 TrkRecoTrk* newTrack = circlefactory.makeTrack(tt, chisum, *m_context, m_bunchT0);
1370 bool permitFlips = true;
1371 bool lPickHits = m_pickHits;
1372 circlefactory.setFlipAndDrop(*newTrack, permitFlips, lPickHits);
1373 // combine hit list
1374 int nDigi = digiToHots(newTrack,vec_truthHit);
1375 int nRecTk = 0;
1376 //if(m_debug>0) newTrack->printAll(std::cout);
1377 //fit
1378 TrkHitList* qhits = newTrack->hits();
1379 TrkErrCode err=qhits->fit();
1380 //cout<<"++++++++++++++++++"<<err.failure()<<endl;
1381 const TrkFit* newFit = newTrack->fitResult();
1382 bool fitSuccess = false;
1383 //test fit result
1384 if (!newFit || (newFit->nActive()<3)) {
1385 if(m_debug>0){
1386 cout << "evt "<<t_eventNum<<" global fit failed ";
1387 if(newFit) cout <<" nAct "<<newFit->nActive();
1388 cout<<endl;
1389 }
1390 //return StatusCode::SUCCESS;
1391 }else{
1392 nRecTk = 1;
1393 fitSuccess = true;
1394 m_nEvtSuccess++;
1395 // if(m_debug>0) cout <<"evt "<<t_eventNum<< " global fit success"<<endl;
1396 if(m_debug>0) newTrack->print(std::cout);
1397
1398 }
1399
1400 vector<int> nfit2d(peak_track,0);
1401 if(m_debug>0) cout<<" hot list:"<<endl;
1402 TrkHotList::hot_iterator hotIter= qhits->hotList().begin();
1403 while (hotIter!=qhits->hotList().end()) {
1404 nfit2d[track_fit]++;
1405 cout <<"(" <<((MdcHit*)(hotIter->hit()))->layernumber()
1406 <<","<<((MdcHit*)(hotIter->hit()))->wirenumber()
1407 <<":"<<hotIter->isActive()<<") ";
1408 hotIter++;
1409 }
1410
1411 m_2d_nFit[track_fit]=nfit2d[track_fit];
1412
1413 //-------------try to out put parameters---
1414 double x_cirtemp;
1415 double y_cirtemp;
1416 double r_temp;
1417 if(err.failure()==0){
1418 //cout<<"++++++++++++++++++"<<err.failure()<<endl;
1419 TrkExchangePar par = newFit->helix(0.);
1420 d0=par.d0();
1421 phi0=par.phi0();
1422 omega=par.omega();
1423 // cout<<"d0: "<<par.d0()<<" "<<"d0temp: "<<-d0_temp<<endl;
1424 // cout<<"phi0: "<<par.phi0()<<" "<<"phi0temp: "<<phi0_temp<<endl;
1425 // cout<<"w: "<<par.omega()<<" "<<"omegatemp: "<<omega_temp<<endl;
1426 // vector<double> b; // vector<double> x_stereo;
1427 r_temp=1./-par.omega();
1428 x_cirtemp=(r_temp-par.d0())*cos(par.phi0()-M_PI/2.);
1429 y_cirtemp=(r_temp-par.d0())*sin(par.phi0()-M_PI/2.);
1430
1431 m_pt[track_fit]=r_temp/333.567;
1432 cout<<"pt_rec: "<<m_pt[track_fit]<<endl;
1433 }
1434 else{
1435 m_nFitFailure[track_fit]=2;
1436 }
1437 //caculate z position
1438
1439 int z_stereoNum= Zposition(nhit,vec_slant,x_cirtemp,y_cirtemp,r_temp,vec_x_east,vec_x_west,vec_y_east, vec_y_west,vec_z_east,vec_z_west, vec_layer, vec_wire,vec_z,vec_zStereo,l,vec_truthHit);
1440 //cout<<"z_stereoNum: "<<z_stereoNum<<endl;
1441 //cout<<"if zposition is finish: "<<endl;
1442
1443 m_zStereoNum=z_stereoNum;
1444 for(int i=0;i<z_stereoNum;i++){
1445
1446 m_zStereoNum=z_stereoNum;
1447 m_zStereo[i]=vec_zStereo[i];
1448 // cout<<"eventNum: "<<t_eventNum<<" "<<"1111111111111111"<<endl;
1449 m_l[i]=l[i];
1450 if(m_debug>0) cout<<" l: "<<m_l[i]<<" "<<"z: "<<vec_zStereo[i]<<endl;
1451 }
1452 // cout<<"333333333333333333"<<endl;
1453
1454 Linefit(z_stereoNum,vec_zStereo,l,tanl,z0);
1455
1456 cout<<"tanl: "<<tanl<<endl;
1457 cout<<"z0: "<<z0<<endl;
1458 z0=dz_mc;
1459 tanl=tanl_mc;
1460 cout<<"tanltruth: "<<tanl<<endl;
1461 cout<<"z0truth: "<<z0<<endl;
1462 //-------------------------------------------5 parameter fit-----------------------
1463
1464 TrkExchangePar tt2(d0,phi0,omega,z0,tanl);
1465 TrkHelixMaker helixfactory;
1466 chisum =0.;
1467 TrkRecoTrk* newTrack2 = helixfactory.makeTrack(tt2, chisum, *m_context, m_bunchT0);
1468 permitFlips = true;
1469 lPickHits = m_pickHits;
1470 helixfactory.setFlipAndDrop(*newTrack2, permitFlips, lPickHits);
1471 // combine hit list
1472 int nDigi2 = digiToHots2(newTrack2,vec_truthHit);
1473 // int nRecTk = 0;
1474 //if(m_debug>0) newTrack->printAll(std::cout);
1475 //fit
1476 TrkHitList* qhits2 = newTrack2->hits();
1477 TrkErrCode err2=qhits2->fit();
1478 //cout<<"++++++++++++++++++"<<err.failure()<<endl;
1479 const TrkFit* newFit2 = newTrack2->fitResult();
1480 // bool fitSuccess = false;
1481
1482 //test fit result
1483 if (!newFit2 || (newFit2->nActive()<5)) {
1484 //fit failed
1485 if(m_debug>0){
1486 cout << "evt "<<t_eventNum<<" global fit failed ";
1487 if(newFit2) cout <<" nAct "<<newFit2->nActive();
1488 cout<<endl;
1489 }
1490 //return StatusCode::SUCCESS;
1491 }else{
1492 //fit success
1493 nFitSucess++;
1494 MdcTrack* mdcTrack = new MdcTrack(newTrack2);
1495 int tkStat = 4;//track find by Hough set stat=4
1496 int tkId = nTdsTk + track_fit;
1497 mdcTrack->storeTrack(tkId, trackList, hitList, tkStat);
1498 nRecTk = 1;
1499 fitSuccess = true;
1500 m_nEvtSuccess++;
1501 // if(m_debug>0) cout <<"evt "<<t_eventNum<< " global fit success"<<endl;
1502 if(m_debug>0) newTrack2->print(std::cout);
1503
1504 }
1505
1506 vector<int> nfit3d(peak_track,0);
1507 if(m_debug>0) cout<<" hot list:"<<endl;
1508 TrkHotList::hot_iterator hotIter2= qhits2->hotList().begin();
1509 while (hotIter2!=qhits2->hotList().end()) {
1510 nfit3d[track_fit]++;
1511 cout <<"(" <<((MdcHit*)(hotIter2->hit()))->layernumber()
1512 <<","<<((MdcHit*)(hotIter2->hit()))->wirenumber()
1513 <<":"<<hotIter2->isActive()<<") ";
1514 hotIter2++;
1515 }
1516
1517 m_3d_nFit[track_fit]=nfit3d[track_fit];
1518
1519 // -------------try to out put parameters---
1520 if(err2.failure()==0){
1521 //cout<<"++++++++++++++++++"<<err2.failure()<<endl;
1522 TrkExchangePar par2 = newFit2->helix(0.);
1523
1524 //cout<<"d0: "<<par2.d0()<<endl;
1525 //cout<<"phi0: "<<par2.phi0()<<endl;
1526 //cout<<"w: "<<par2.omega()<<endl;
1527 //cout<<"z: "<<par2.z0()<<endl;
1528 //cout<<"tanl: "<<par2.tanDip()<<endl;
1529
1530 r_temp=1./-par2.omega();
1531 x_cirtemp=(r_temp-par2.d0())*cos(par2.phi0()-M_PI/2.);
1532 y_cirtemp=(r_temp-par2.d0())*sin(par2.phi0()-M_PI/2.);
1533
1534 m_d0[track_fit]=par2.d0();
1535 m_phi0[track_fit]=par2.phi0();
1536 m_omega[track_fit]=par2.omega();
1537 m_z0[track_fit]=par2.z0();
1538 m_tanl[track_fit]=par2.tanDip();
1539
1540 double pxy=r_temp/333.567;
1541 double pz=pxy*par2.tanDip();
1542 double pxyz=sqrt(pxy*pxy+pz*pz);
1543 m_pt2[track_fit]=pxy;
1544 m_pz[track_fit]=pxy*par2.tanDip();
1545 m_pxyz[track_fit]=pxyz;
1546 cout<<"track"<<track_fit<<": "<<"pt_rec2: "<<m_pt2[track_fit]<<endl;
1547 cout<<"eventNum: "<<t_eventNum<<" "<<"track"<<track_fit<<": "<<"pz_rec: "<<m_pz[track_fit]<<endl;
1548 }
1549 else{
1550 m_nFitFailure[track_fit]=3;
1551 }
1552 // vector<double> b; // vector<double> x_stereo;
1553 // r_temp=1./-par.omega();
1554 // x_cirtemp=(r_temp-par.d0())*cos(par.phi0()-M_PI/2.);
1555 // y_cirtemp=(r_temp-par.d0())*sin(par.phi0()-M_PI/2.);
1556 // if(m_data==0) {
1557 // m_mc_pt=r_temp/333.567;
1558 // cout<<"pt_rec: "<<m_mc_pt<<endl;
1559 // }
1560 // if(m_data==1) {m_pt=r_temp/333.567;}
1561 //
1562 // else{
1563 // if(m_data==0) {m_mc_nfitfailure=2;}
1564 // if(m_data==1){m_nFitFailure=2;}
1565 // }
1566
1567
1568 //------------------------------------clean-------------------------
1569 }
1570 //vector<double>().swap(vec_u);
1571 //vector<double>().swap(vec_v);
1572 //vector<double>().swap(vec_p);
1573 //vector<double>().swap(vec_x);
1574 //vector<double>().swap(vec_y);
1575 //vector<double>().swap(vec_z);
1576 //vector<double>().swap(vec_x_east);
1577 //vector<double>().swap(vec_x_west);
1578 //vector<double>().swap(vec_y_east);
1579 //vector<double>().swap(vec_y_west);
1580 //vector<double>().swap(vec_z_east);
1581 //vector<double>().swap(vec_z_west);
1582 //vector<double>().swap(vec_zStereo);
1583 //vector<double>().swap(l);
1584 //vector<int>().swap(vec_layer);
1585 //vector<int>().swap(vec_wire);
1586}
1587m_nFitSucess=nFitSucess;
1588cout<<" Hough finder find "<<nFitSucess<<" tot "<<nTdsTk<< endl;
1589
1590delete h1;
1591delete h2;
1592cout<<"break: ????"<<endl;
1593ntuplehit->write();
1594cout<<"finish event "<<t_eventNum<<endl;
1595cout<<endl;
1596cout<<endl;
1597return StatusCode::SUCCESS;
1598}
1599
1600// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
1601StatusCode MdcHoughFinder::finalize() {
1602 MsgStream log(msgSvc(), name());
1603 log << MSG::INFO << "in finalize()" << endreq;
1604 return StatusCode::SUCCESS;
1605}
1606int MdcHoughFinder::digiToHots(TrkRecoTrk* newTrack,vector<bool> vec_truthHit){
1607 TrkHitList* trkHitList = newTrack->hits();
1608 //HepPoint3D poca; Hep3Vector tempDir;
1609 //getInfo(helix,0,poca,tempDir);
1610
1611 //if(m_debug>3)std::cout<<"track poca "<<poca<<std::endl;
1612 // new MdcRecoHitOnTracks
1613 //SmartDataPtr<MdcDigiCol> mdcDigiCol(eventSvc(),"/Event/Digi/MdcDigiCol");
1614 //if (!mdcDigiCol) {
1615 // cout << "Could not find MdcDigiCol" << endl;
1616 // return 0;
1617 //}
1618
1619 //if(mdcDigiCol->size()<5) return mdcDigiCol->size();
1620
1621 //if(m_debug>0) std::cout<<" nDigi = "<<mdcDigiCol->size()<<std::endl;
1622 int digiId=0;
1623 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec();
1624 MdcDigiVec::iterator iter1 = mdcDigiVec.begin();
1625 for (;iter1 != mdcDigiVec.end(); iter1++, digiId++) {
1626 // MdcDigiVec::iterator iter = mdcDigiCol->begin();
1627 //for (int i=0;iter != mdcDigiCol->end(); iter++,i++ ) {
1628 if(vec_truthHit[digiId]==false) continue;
1629 const MdcDigi* aDigi = *iter1;
1630 MdcHit *hit = new MdcHit(aDigi, m_gm);
1631 Identifier id = aDigi->identify();
1632 int layer = MdcID::layer(id);
1633 int wire = MdcID::wire(id);
1634
1635 if ((layer>=0&&layer<=7)||(layer>=20&&layer<=35)) {continue;}
1636 // if ((layer>=0&&layer<=7)||(layer>=20)) {continue;}
1637 //if (layer<=35) {continue;}
1638
1639 if(m_debug>0) std::cout<<" ("<<layer<<","<<wire<<") "<<std::endl;
1640 hit->setCalibSvc(m_mdcCalibFunSvc);
1641 hit->setCountPropTime(true); //hit->setCountPropTime(m_countPropTime);
1642 // if(m_debug>0) cout<<"hitx: "<<hit->x()<<" "<<"hity: "<<hit->y()<<" "<<endl;
1643 //hit->setZTruth(t_mcZ[layer][wire]);
1644
1645 // new fltLen and ambig
1646 int newAmbig = 0;
1647 // calc. position of this point
1648 //m_hitCol->push_back(hit);
1649 MdcRecoHitOnTrack temp(*hit, newAmbig, m_bunchT0*1.e9);//m_bunchT0 nano second here
1650 MdcHitOnTrack* newhot = &temp;
1651 // cout<<"ambig: "<<newhot->ambig()<<" xxxxxxx"<<endl;
1652 double fltLen = m_gm->Layer(layer)->rMid();
1653 newhot->setFltLen(fltLen);
1654 //newhot->setFltLen(0);
1655 //newhot->setHitRms(0.02);
1656 //if(m_debug>1) std::cout<<name()<<" ("<<layer<<","<<wire<<") fltLen "<<newhot->fltLen()<<std::endl;
1657 if(m_debug>2) std::cout<<" ("<<layer<<","<<wire<<",rt "<<hit->rawTime()<<",dt "<<hit->driftTime(m_bunchT0,0)<<") T0 " << m_mdcCalibFunSvc->getT0(layer,wire) <<" timewalk "<< m_mdcCalibFunSvc->getTimeWalk(layer, aDigi->getChargeChannel())<<endl;
1658 //if(hit->driftTime(m_bunchT0,0)>800) continue;
1659 if(hit->driftDist(m_bunchT0,0)>1.2) continue;
1660 trkHitList->appendHot(newhot);
1661 //}//end of loop hits
1662 }
1663 if(m_debug>1) std::cout<<std::endl;
1664 return trkHitList->nHit();
1665}
1666
1667double MdcHoughFinder::CFtrans(double x,double y){
1668 return x/(x*x+y*y);
1669}
1670
1671
1672// MsgStream log(msgSvc(), name());
1673// log << MSG::INFO << "in finalize()" << endreq;
1674
1675// return StatusCode::SUCCESS;
1676
1677
1678int MdcHoughFinder::SlantId(int layer){
1679 if ((layer>=0&&layer<=7)||(layer>=20&&layer<=35)) {
1680 if ((layer>=0&&layer<=3)||(layer>=20&&layer<=23)||(layer>=28&&layer<=31)){
1681 // m_slant[digiId]=-1;
1682 return -1;
1683 }else{
1684 return 1;
1685 }
1686 } else{
1687 return 0;
1688 }
1689}
1690
1691
1692int MdcHoughFinder::GetMcInfo(){
1693 //t_t0Truth=-1;
1694 MsgStream log(msgSvc(), name());
1695 StatusCode sc;
1696 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
1697 if (!mcParticleCol) {
1698 log << MSG::ERROR << "Could not find McParticle" << endreq;
1699 return -999;
1700 }
1701 int itk = 0;
1702 int t_qTruth = 0;
1703 int t_pidTruth = -999;
1704 int t_McTkId = -999;
1705 t_nTrkMC=0;
1706 Helix* mchelix;
1707 McParticleCol::iterator iter_mc = mcParticleCol->begin();
1708 for (;iter_mc != mcParticleCol->end(); iter_mc++ ) {
1709 if(!(*iter_mc)->primaryParticle()) continue;
1710 t_pidTruth = (*iter_mc)->particleProperty();
1711
1712 if(m_debug>2) cout<< "tk "<<itk<< " pid="<< t_pidTruth<< endl;
1713 if((m_pid!=0) && (t_pidTruth != m_pid)){ continue; }
1714
1715 //t_t0Truth=(*iter_mc)->initialPosition().t();
1716
1717 if(itk>=10) break;
1718
1719 int pid = t_pidTruth;
1720 if( pid == 0 ) {
1721 log << MSG::WARNING << "Wrong particle id" <<endreq;
1722 continue;
1723 }else{
1724 IPartPropSvc* p_PartPropSvc;
1725 static const bool CREATEIFNOTTHERE(true);
1726 StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
1727 if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) {
1728 std::cout<< " Could not initialize Particle Properties Service" << std::endl;
1729 }
1730 HepPDT::ParticleDataTable* p_particleTable = p_PartPropSvc->PDT();
1731 std::string name;
1732 if( p_particleTable->particle(pid) ){
1733 name = p_particleTable->particle(pid)->name();
1734 t_qTruth = p_particleTable->particle(pid)->charge();
1735 }else if( p_particleTable->particle(-pid) ){
1736 name = "anti " + p_particleTable->particle(-pid)->name();
1737 t_qTruth = (-1)*p_particleTable->particle(-pid)->charge();
1738 }
1739 if(m_debug>2) std::cout << " particle "<<name<<" charge= "<<t_qTruth<<std::endl;
1740 }
1741
1742 t_McTkId = (*iter_mc)->trackIndex();
1743
1744 double px = (*iter_mc)->initialFourMomentum().px();//GeV
1745 double py = (*iter_mc)->initialFourMomentum().py();//GeV
1746 double pz = (*iter_mc)->initialFourMomentum().pz();//GeV
1747
1748 m_pzTruth=pz;
1749
1750 mchelix = new Helix((*iter_mc)->initialPosition().v(), (*iter_mc)->initialFourMomentum().v(), t_qTruth);//cm
1751 mchelix->pivot( HepPoint3D(0.,0.,0.) );
1752
1753 if(m_debug>0){
1754 std::cout<<"Truth tk "<<itk<<" pid "<<pid<<" charge "<<t_qTruth
1755 << " helix "<< mchelix->a()<<" p "<<mchelix->momentum(0.)<<std::endl;
1756 }
1757 t_nTrkMC++;
1758 itk++;
1759 }
1760 if(t_nTrkMC!=1) {
1761 std::cout<<"WARNING run "<<t_runNum<<" evt "<<t_eventNum<<" not single event. nTrkMc="<<t_nTrkMC<<std::endl;
1762 return -999;
1763 }else{
1764 if(m_debug>2) std::cout<<"nTrkMc=1"<<std::endl;
1765 }
1766 //--fill truth in evt
1767
1768 m_pidTruth = t_pidTruth;
1769 m_costaTruth = mchelix->cosTheta();
1770 m_phi0Truth = mchelix->phi0();
1771 m_drTruth = mchelix->dr();
1772 m_dzTruth = mchelix->dz();
1773 m_ptTruth = mchelix->pt();
1774 m_pTruth = mchelix->momentum(0.).mag();
1775 m_qTruth = t_qTruth;
1776 dz_mc =mchelix->dz();
1777 tanl_mc =mchelix->tanl();
1778 return 1;
1779}
1780
1781
1782int MdcHoughFinder::LeastSquare(int n,int nHitAxialSelect,vector<double> vec_x,vector<double> vec_y, vector<int> vec_slant,vector<int> vec_layer,vector<bool> vec_truthHit,double &d0,double &phi0, double &omega){
1783 //if(nHitAxialSelect<3) {
1784 //m_nFitFailure=1;
1785 // return 1;
1786 // }
1787 //else{
1788 double x_sum=0;
1789 double y_sum=0;
1790 double x2_sum=0;
1791 double y2_sum=0;
1792 double x3_sum=0;
1793 double y3_sum=0;
1794 double x2y_sum=0;
1795 double xy2_sum=0;
1796 double xy_sum=0;
1797 double a1=0;
1798 double a2=0;
1799 double b1=0;
1800 double b2=0;
1801 double c1=0;
1802 double c2=0;
1803 double x_aver,y_aver,r2;
1804 for(int i=0;i<n;i++){
1805 if(vec_truthHit[i]==false) continue;
1806 if(vec_slant[i]!=0) continue;
1807 //cout<<"hitNum: "<<i<<" "<<"layer:"<<" "<<vec_layer[i]<<" "<<endl;
1808 x_sum=x_sum+vec_x[i];
1809 y_sum=y_sum+vec_y[i];
1810 x2_sum=x2_sum+vec_x[i]*vec_x[i];
1811 y2_sum=y2_sum+vec_y[i]*vec_y[i];
1812 x3_sum=x3_sum+vec_x[i]*vec_x[i]*vec_x[i];
1813 y3_sum=y3_sum+vec_y[i]*vec_y[i]*vec_y[i];
1814 x2y_sum=x2y_sum+vec_x[i]*vec_x[i]*vec_y[i];
1815 xy2_sum=xy2_sum+vec_x[i]*vec_y[i]*vec_y[i];
1816 xy_sum=xy_sum+vec_x[i]*vec_y[i];
1817 }
1818 a1=x2_sum-x_sum*x_sum/n;
1819 a2=xy_sum-x_sum*y_sum/n;
1820 b1=a2;
1821 b2=y2_sum-y_sum*y_sum/n;
1822 c1=(x3_sum+xy2_sum-x_sum*(x2_sum+y2_sum)/n)/2.;
1823 c2=(y3_sum+x2y_sum-y_sum*(x2_sum+y2_sum)/n)/2.;
1824
1825 x_aver=x_sum/n;
1826 y_aver=y_sum/n;
1827
1828 m_x_circle =(b1*c2-b2*c1)/(b1*b1-a1*b2);
1829 m_y_circle =(b1*c1-a1*c2)/(b1*b1-a1*b2);
1830
1831 double x_cirtemp =(b1*c2-b2*c1)/(b1*b1-a1*b2);
1832 double y_cirtemp =(b1*c1-a1*c2)/(b1*b1-a1*b2);
1833
1834 r2=(x2_sum+y2_sum-2*x_cirtemp*x_sum-2*y_cirtemp*y_sum)/n+x_cirtemp*x_cirtemp+y_cirtemp*y_cirtemp;
1835
1836 m_r=sqrt(r2);
1837
1838 double r_temp=sqrt(r2);
1839 double d0_temp=sqrt(x_cirtemp*x_cirtemp+y_cirtemp*y_cirtemp)-r_temp;
1840 //cout<<"rtemp: "<<r_temp<<" "<<endl;
1841 //cout<<"xtemp: "<<x_cirtemp<<" "<<endl;
1842 //cout<<"ytemp: "<<y_cirtemp<<" "<<endl;
1843 //cout<<"d0temp: "<<-d0_temp<<" "<<endl;
1844 double pt_temp=r_temp/333.567;
1845 cout<<"pt_temp: "<<pt_temp<<endl;
1846 //--------------------------------------------------add drift and fit---------------------------
1847
1848 d0=-d0_temp;
1849 // double phi0=atan2(y_cirtemp,x_cirtemp);
1850 // double omega=-1./r_temp;
1851 phi0=atan2(y_cirtemp,x_cirtemp)+M_PI/2.;
1852 omega=-1./r_temp;
1853 double z0=0.;
1854 double tanl=0.;
1855
1856 // double phi0_temp=phi0;
1857 // double omega_temp=omega;
1858 //d0=0;
1859 //phi0=-3.08301;
1860 //omega=-0.0117116;
1861 //double z0=0.0;
1862 //double tanl=0.0;
1863 if(m_debug<0)std::cout<<" before global fit Helix PatPar :"<<d0<<","<<phi0<<","<<omega<<","<<z0<<","<<tanl<< std::endl;
1864 return 0;
1865 //}
1866}
1867
1868int MdcHoughFinder::Zposition(int n,vector<int> vec_slant,double x_cirtemp,double y_cirtemp,double r_temp,vector<double> vec_x_east,vector<double> vec_x_west,vector<double> vec_y_east,vector<double> vec_y_west,vector<double> vec_z_east,vector<double> vec_z_west,vector<int> vec_layer,vector<int> vec_wire,vector<double> vec_z,vector<double>& vec_zStereo,vector<double>& l,vector<bool> vec_truthHit){
1869 double x_stereo;
1870 double y_stereo;
1871 // double theta;
1872 int stereoId=0;
1873 //vector<double> z_stereo;
1874 // vector<double> vec_zStereo;
1875 int nDropDelta=0;
1876 for(int i=0;i<n;i++){
1877 if(vec_truthHit[i]==false) continue;
1878 double k;
1879 double b;
1880 if(vec_slant[i]!=0){
1881 if(vec_x_east[i]!=vec_x_west[i]){
1882 k=(vec_y_east[i]-vec_y_west[i])/(vec_x_east[i]-vec_x_west[i]);
1883 b=-k*vec_x_east[i]+vec_y_east[i];
1884 // cout<<"eventNum: "<<t_eventNum<<" "<<"layerid: "<<" "<<m_layer[i]<<"cellid: "<<m_cell[i]<<" "<<"k: "<<k<<" "<<"b: "<<b<<" "<<m_y_east[i]<<" "<<m_y_west[i]<<" "<<m_x_east[i]<<" "<<m_x_west[i]<<" "<<endl;
1885 double delta=4*(k*b-k*y_cirtemp-x_cirtemp)*(k*b-k*y_cirtemp-x_cirtemp)-4*(k*k+1)*(x_cirtemp*x_cirtemp+b*b+y_cirtemp*y_cirtemp-2*b*y_cirtemp-r_temp*r_temp);
1886 //cout<<"layer: "<<vec_layer[i]<<"wire: "<<vec_wire[i]<<endl;
1887 //cout<<"x_east: "<<vec_x_east[i]<<" "<<"x_west: "<<vec_x_west[i]<<endl;
1888 //cout<<"k: "<<k<<" "<<"b: "<<b<<endl;
1889 //cout<<"x_cirtemp: "<<x_cirtemp<<" "<<"y_cirtemp: "<<y_cirtemp<<endl;
1890 if(delta<0) {
1891 nDropDelta++;
1892 continue;
1893 }
1894 //cout<<"delta: "<<delta<<endl;
1895 double x1=(2*x_cirtemp+2*k*y_cirtemp-2*k*b+sqrt(delta))/(2*(k*k+1));
1896 double x2=(2*x_cirtemp+2*k*y_cirtemp-2*k*b-sqrt(delta))/(2*(k*k+1));
1897 double x_middle=(vec_x_east[i]+vec_x_west[i])/2.;
1898 // if(abs(x1-vec_x[i])>abs(x2-vec_x[i])) {x_stereo=x2;}
1899 if(abs(x1-x_middle)>abs(x2-x_middle)) {x_stereo=x2;}
1900 else {x_stereo=x1;}
1901 y_stereo=k*x_stereo+b;
1902 }
1903 else{
1904 //cout<<"layer: "<<vec_layer[i]<<"wire: "<<vec_wire[i]<<endl;
1905 //cout<<"x_east: "<<vec_x_east[i]<<" "<<"x_west: "<<vec_x_west[i]<<endl;
1906 x_stereo=vec_x_east[i];
1907 double y1=sqrt(r_temp*r_temp-(x_cirtemp-x_stereo)*(x_cirtemp-x_stereo))+y_cirtemp;
1908 double y2=y_cirtemp-sqrt(r_temp*r_temp-(x_cirtemp-x_stereo)*(x_cirtemp-x_stereo));
1909 double y_middle=(vec_y_east[i]+vec_y_west[i])/2.;
1910 if(abs(y1-y_middle)>abs(y2-y_middle)) {y_stereo=y2;}
1911 else{y_stereo=y1;}
1912 }
1913 // cout<<"layer: "<<vec_layer[i]<<"wire: "<<vec_wire[i]<<"x_stereo: "<<x_stereo<<" "<<"y_stereo: "<<y_stereo<<endl;
1914
1915 double theta_temp;
1916 if(x_cirtemp==0||x_stereo-x_cirtemp==0){
1917 theta_temp=0;
1918 }
1919 else{
1920 theta_temp=M_PI-atan2(y_stereo-y_cirtemp,x_stereo-x_cirtemp)+atan2(y_cirtemp,x_cirtemp);
1921 if(theta_temp>2*M_PI){
1922 theta_temp=theta_temp-2*M_PI;
1923 }
1924 if(theta_temp<0){
1925 theta_temp=theta_temp+2*M_PI;
1926 }
1927 }
1928 //cout<<"thetatemp: "<<theta_temp<<endl;
1929 //if(theta_temp>2*M_PI) {
1930 // theta_temp=theta_temp-2*M_PI;
1931 //}
1932 //if(theta_temp<0.) {
1933 // theta_temp=theta_temp+2*M_PI;
1934 //}
1935 //cout<<"debug: if thetatemp is cout over"<<endl;
1936 l.push_back(r_temp*theta_temp);
1937 //cout<<"debug: if l is pushback"<<endl;
1938 //l.push_back(M_PI-atan2(y_stereo,x_stereo)-atan2(y_cirtemp,x_cirtemp));
1939 //l.push_back(r_temp*(M_PI-atan2(y_stereo,x_stereo)-atan2(y_cirtemp,x_cirtemp)));
1940
1941 double d1=sqrt((x_stereo-vec_x_west[i])*(x_stereo-vec_x_west[i])+(y_stereo-vec_y_west[i])*(y_stereo-vec_y_west[i]));
1942 double d2=sqrt((vec_x_east[i]-vec_x_west[i])*(vec_x_east[i]-vec_x_west[i])+(vec_y_east[i]-vec_y_west[i])*(vec_y_east[i]-vec_y_west[i]));
1943 vec_zStereo.push_back(vec_z_west[i]-(vec_z_west[i]-vec_z_east[i])*d1/d2);
1944 //cout<<"debug: if vec_zStereo is pushback"<<endl;
1945 //if(stereoId==0) {vec_zStereo.push_back(vec_z_west[i]-(vec_z_west[i]-vec_z_east[i])*d1/d2);}
1946 //if(stereoId>0){
1947 // if(vec_layer[i]==vec_layer.at(i-1)) {
1948 // vec_zStereo[(stereoId-1)]=((z_stereo.at(stereoId)+z_stereo.at(stereoId-1))/2.);
1949 // vec_zStereo.push_back((z_stereo.at(stereoId)+z_stereo.at(stereoId-1))/2.);
1950 // }
1951 // else{
1952 // vec_zStereo.push_back(vec_z_west[i]-(vec_z_west[i]-vec_z_east[i])*d1/d2);
1953 // }
1954 //}
1955
1956 // double delphi=atan2(y_stereo,x_stereo)-par.phi0();
1957 //cout<<"eventNum: "<<t_eventNum<<" "<<"layerid: "<<" "<<m_layer[i]<<"cellid: "<<" "<<m_cell[i]<<" "<<"z_stereo: "<<z_stereo<<" "<<d1<<" "<<d2<<" "<<endl;
1958 //cout<<"eventNum: "<<t_eventNum<<" "<<"layerid: "<<" "<<m_layer[i]<<"cellid: "<<m_cell[i]<<" "<<"k: "<<k<<" "<<"b: "<<b<<" "<<m_y_east[i]<<" "<<m_y_west[i]<<" "<<m_x_east[i]<<" "<<m_x_west[i]<<" "<<"x_stereo: "<<" "<<x_stereo<<" "<<"x: "<<vec_x[i]<<"y_stereo: "<<" "<<y_stereo<<" "<<"y: "<<vec_y[i]<<" "<<"zstereo: "<<z_stereo<<endl;
1959 // cout<<"layerid: "<<" "<<vec_layer[i]<<" "<<"cellid: "<<vec_wire[i]<<" "<<"x_stereo: "<<" "<<x_stereo<<" "<<"y_stereo: "<<" "<<y_stereo<<" "<<"zstereo: "<<z_stereo.at(stereoId)<<" "<<vec_zStereo.at(stereoId)<<" "<<"ztruth: "<<vec_z[i]<<endl;
1960
1961 // cout<<"layer: "<<vec_layer[i]<<"wire: "<<vec_wire[i]<<"z: "<<vec_zStereo[stereoId]<<" "<<"l: "<<l[stereoId]<<" "<<vec_truthHit[i]<<endl;
1962 stereoId=stereoId+1;
1963 //cout<<"theta: "<<theta<<endl;
1964 }
1965 else {k=b=0;}
1966 }
1967 cout<<"nDropDelta: "<<nDropDelta<<endl;
1968 // cout<<"if for is finished : "<<endl;
1969 //int digiId=0;
1970 //for(int i=0;i<n;i++){
1971 // if(vec_slant[i]!=0){
1972 // cout<<"layerid: "<<" "<<vec_layer[i]<<" "<<"cellid: "<<vec_wire[i]<<" "<<"zstereo: "<<z_stereo.at(digiId)<<" "<<vec_zStereo.at(digiId)<<" "<<"ztruth: "<<vec_z[i]<<endl;
1973 // digiId++;
1974 // }
1975 //}
1976 //double sum_zstereo=0;
1977 //for(int i=0;i<stereoId;i++){
1978 // sum_zstereo=sum_zstereo+z_stereo[i];
1979 //}
1980 //cout<<"sum_zstereo: "<<sum_zstereo<<endl;
1981 return stereoId;
1982}
1983
1984
1985void MdcHoughFinder::Linefit(int z_stereoNum,vector<double> vec_zStereo,vector<double> l,double& tanl,double& z0){
1986
1987 double l_sum=0;
1988 double z_sum=0;
1989 double l2_sum=0;
1990 double z2_sum=0;
1991 double lz_sum=0;
1992 for(int i=0;i<z_stereoNum;i++){
1993 // cout<<"event: "<<t_eventNum<<" "<<"z: "<<vec_zStereo[i]<<" "<<"l: "<<l[i]<<endl;
1994 l_sum=l_sum+l[i];
1995 z_sum=z_sum+vec_zStereo[i];
1996 l2_sum=l2_sum+l[i]*l[i];
1997 // cout<<l2_sum<<endl;
1998 z2_sum=z2_sum+vec_zStereo[i]*vec_zStereo[i];
1999 lz_sum=lz_sum+l[i]*vec_zStereo[i];
2000 }
2001 double b=(l2_sum*z_sum-lz_sum*l_sum)/(z_stereoNum*l2_sum-l_sum);
2002 double k=(z_stereoNum*lz_sum-l_sum*z_sum)/(z_stereoNum*l2_sum-l_sum*l_sum);
2003 z0=b;
2004 tanl=k;
2005 cout<<k<<" "<<b<<endl;
2006}
2007
2008
2009
2010int MdcHoughFinder::digiToHots2(TrkRecoTrk* newTrack,vector<bool> vec_truthHit){
2011 TrkHitList* trkHitList = newTrack->hits();
2012 //HepPoint3D poca; Hep3Vector tempDir;
2013 //getInfo(helix,0,poca,tempDir);
2014
2015 //if(m_debug>3)std::cout<<"track poca "<<poca<<std::endl;
2016 // new MdcRecoHitOnTracks
2017 //SmartDataPtr<MdcDigiCol> mdcDigiCol(eventSvc(),"/Event/Digi/MdcDigiCol");
2018 //if (!mdcDigiCol) {
2019 // cout << "Could not find MdcDigiCol" << endl;
2020 // return 0;
2021 //}
2022
2023 //if(mdcDigiCol->size()<5) return mdcDigiCol->size();
2024
2025 //if(m_debug>0) std::cout<<" nDigi = "<<mdcDigiCol->size()<<std::endl;
2026 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec();
2027 MdcDigiVec::iterator iter1 = mdcDigiVec.begin();
2028 int digiId=0;
2029 for (;iter1 != mdcDigiVec.end(); iter1++, digiId++) {
2030 //MdcDigiVec::iterator iter = mdcDigiCol->begin();
2031 //for (int i=0;iter != mdcDigiCol->end(); iter++,i++ ) {
2032 if(vec_truthHit[digiId]==false) continue;
2033 const MdcDigi* aDigi = *iter1;
2034 MdcHit *hit = new MdcHit(aDigi, m_gm);
2035 Identifier id = aDigi->identify();
2036 int layer = MdcID::layer(id);
2037 int wire = MdcID::wire(id);
2038
2039 // if ((layer>=0&&layer<=7)||(layer>=20&&layer<=35)) {continue;}
2040 // if ((layer>=0&&layer<=7)||(layer>=20)) {continue;}
2041 //if (layer<=35) {continue;}
2042
2043 if(m_debug>0) std::cout<<" ("<<layer<<","<<wire<<") "<<std::endl;
2044 hit->setCalibSvc(m_mdcCalibFunSvc);
2045 hit->setCountPropTime(true); //hit->setCountPropTime(m_countPropTime);
2046 // if(m_debug>0) cout<<"hitx: "<<hit->x()<<" "<<"hity: "<<hit->y()<<" "<<endl;
2047 //hit->setZTruth(t_mcZ[layer][wire]);
2048
2049 // new fltLen and ambig
2050 int newAmbig = 0;
2051 // calc. position of this point
2052 //m_hitCol->push_back(hit);
2053 MdcRecoHitOnTrack temp(*hit, newAmbig, m_bunchT0*1.e9);//m_bunchT0 nano second here
2054 MdcHitOnTrack* newhot = &temp;
2055 // cout<<"ambig: "<<newhot->ambig()<<" xxxxxxx"<<endl;
2056 double fltLen = m_gm->Layer(layer)->rMid();
2057 newhot->setFltLen(fltLen);
2058 //newhot->setFltLen(0);
2059 //newhot->setHitRms(0.02);
2060 //if(m_debug>1) std::cout<<name()<<" ("<<layer<<","<<wire<<") fltLen "<<newhot->fltLen()<<std::endl;
2061 if(m_debug>2) std::cout<<" ("<<layer<<","<<wire<<",rt "<<hit->rawTime()<<",dt "<<hit->driftTime(m_bunchT0,0)<<") T0 " << m_mdcCalibFunSvc->getT0(layer,wire) <<" timewalk "<< m_mdcCalibFunSvc->getTimeWalk(layer, aDigi->getChargeChannel())<<endl;
2062 //if(hit->driftTime(m_bunchT0,0)>800) continue;
2063 if(hit->driftDist(m_bunchT0,0)>1.2) continue;
2064 trkHitList->appendHot(newhot);
2065 //}//end of loop hits
2066 }
2067 if(m_debug>1) std::cout<<std::endl;
2068 return trkHitList->nHit();
2069}
2070void MdcHoughFinder:: Multi_array(int binX,int binY){
2071 int **count=new int*[binX];
2072 for(int i=0;i<binX;i++){
2073 count[i]=new int [binY];
2074 }
2075 //vector<int> vec_selectNum[binX][binY];
2076 int ***selectNum=new int**[binX];
2077 for(int i=0;i<binX;i++){
2078 selectNum[i]=new int* [binY];
2079 for(int j=0;j<binY;j++){
2080 //selectNum[i][j]=new int[];
2081 selectNum[i][j]=new int [count[i][j]];
2082 }
2083 }
2084}
2085void MdcHoughFinder::FillCells(TH2D *h1,int n,int binX,int binY,vector<double> vec_u,vector<double> vec_v,vector<int> vec_layer,vector<int> vec_wire,vector< vector <int> >& countij,vector< vector < vector<int> > >& vec_selectNum,vector< vector < vector<int> > >& vec_selectHit){
2086 double binWidth=M_PI/binX;
2087 double binHigh=2*t_maxP/binY;
2088
2089 // vector<int>** selectNum=new vector<int>*[binX];
2090 // for(int i=0;i<binX;i++){
2091 // selectNum[i]=new vector<int>[binY];
2092 // for(int j=0;j<binY;j++){
2093 // selectNum[i][j]=new vector<int>;
2094 // }
2095 // }
2096 // TH2D *h = new TH2D("hough","hough",binX,0,3.14,binY,-max_p,max_p);
2097
2098 for(int i=0;i<binX;i++){
2099 for(int j=0;j<binY;j++){
2100 int count=0;
2101 double binLeft=i*binWidth;
2102 double binRight=(i+1)*binWidth;
2103 double binDown=-t_maxP+j*binHigh;
2104 double binUp=-t_maxP+(j+1)*binHigh;
2105 double binMid=(binLeft+binRight)/2;
2106 for(int k=0;k<n;k++){
2107 double lineLeft=vec_u[k]*cos(binLeft)+vec_v.at(k)*sin(binRight);
2108 double lineRight=vec_u[k]*cos(binLeft)+vec_v.at(k)*sin(binRight);
2109 double lineMid=vec_u[k]*cos(binMid)+vec_v[k]*sin(binMid);
2110 if(((lineLeft<binUp)&&(lineLeft>binDown)) || ((lineRight<binUp)&&(lineRight>binDown)) || (((lineLeft-binUp)*(binDown-lineRight)>0))){
2111 vec_selectNum[i][j].push_back(k);
2112 int layerId=vec_layer[k];
2113 int cellId=vec_wire[k];
2114 int hit_temp=(int)layerId*1000+cellId;
2115 vec_selectHit[i][j].push_back(hit_temp);
2116 count++;
2117 }
2118 }
2119 countij[i][j]=count;
2120
2121 //for(int l=0,m=vec_selectNum[i][j].size();l<m;++l){
2122 // {// cout << "section " << "(" << i << ", " << j << ") statistics: " << count[i][j] << endl;
2123 // cout << "("<<i<<","<<j<<")"<<": "<<vec_selectNum[i][j][l] << ", "<<endl;
2124 // // trueTrk[vec_selectNum[i][j].at(l)]=true;
2125 // }
2126 //}
2127 if(m_debug>0) cout<<"countij: "<<countij[i][j]<<endl;
2128 }
2129 }
2130 // for(int i=0;i<binX;i++){
2131 // for(int j=0;j<binY;j++){
2132 // for(int k=0;k<countij[i][j];k++){
2133 // cout<<"("<<i<<","<<j<<")"<<vec_selectNum[i][j][k]<<endl;
2134 // }
2135 // }
2136 // }
2137
2138 for(int i=0;i<binX;i++){
2139 for(int j=0;j<binY;j++){
2140 h1->SetBinContent(i+1,j+1,countij[i][j]);
2141 }
2142 }
2143}
2144void MdcHoughFinder::RhoTheta(int numCross,int nhit,vector<double> vec_u,vector<double> vec_v,vector<double>& vec_rho,vector<double>& vec_theta,vector< vector<int> >& vec_hitNum){
2145 for(int i=0;i<nhit;i++){
2146 for(int j=i+1;j<nhit;j++){
2147 double k,b,x_cross,y_cross;
2148 double rho_temp,theta_temp;
2149 k=(vec_v[i]-vec_v[j])/(vec_u[i]-vec_u[j]);
2150 b=vec_v[i]-k*vec_u[i];
2151 x_cross=-b/(k+1/k);
2152 y_cross=b/(k*k+1);
2153 rho_temp=sqrt(x_cross*x_cross+y_cross*y_cross);
2154 theta_temp=atan2(y_cross,x_cross);
2155 if(theta_temp<0) {
2156 theta_temp=theta_temp+M_PI;
2157 rho_temp=-rho_temp;
2158 }
2159 vec_rho.push_back(rho_temp);
2160 vec_theta.push_back(theta_temp);
2161 vec_hitNum[0].push_back(i);
2162 vec_hitNum[1].push_back(j);
2163 //cout<<"u:"<<vec_u[i]<<" "<<"v:"<<vec_v[i]<<" "<<"("<<vec_theta[i]<<","<<vec_rho[i]<<")"<<endl;
2164 }
2165 }
2166 //cout<<"numCross: "<<numCross<<endl;
2167 m_nCross=numCross;
2168
2169 for(int i=0;i<numCross;i++){
2170 m_rho[i]=vec_rho[i];
2171 m_theta[i]=vec_theta[i];
2172 //cout<<" "<<vec_theta[i]<<" "<<vec_rho[i]<<endl;
2173 }
2174}
2175void MdcHoughFinder::FillHist(TH2D *h1,vector<double> vec_u,vector<double> vec_v,int nhit,vector< vector < vector< int > > > vec_selectNum){
2176 for(int n=0;n<nhit;n++){
2177 for(int i=0;i<binX;i++){
2178 double binwidth=M_PI/binX;
2179 double binhigh=2*t_maxP/binY;
2180 double binMid=(i+0.5)*binwidth;
2181 double rho=vec_u[n]*cos(binMid)+vec_v[n]*sin(binMid);
2182 int j=(int)((rho+t_maxP)/binhigh);
2183 int count=0;
2184 count=h1->GetBinContent(i+1,j+1);
2185 count++;
2186 h1->SetBinContent(i+1,j+1,count);
2187 }
2188 }
2189 //for(int i=0;i<binX;i++){
2190 // double binwidth=M_PI/binX;
2191 // double binhigh=2*t_maxP/binY;
2192 // double binMid=(i+0.5)*binwidth;
2193 // for(int n=0;n<nhit;n++){
2194 // double rho=vec_u[n]*cos(binMid)+vec_v[n]*sin(binMid);
2195 // int j=(int)((rho+t_maxP)/binhigh);
2196 // int count=h1->GetBinContent(i+1,j+1);
2197 // count++;
2198 // h1->SetBinContent(i+1,j+1,count);
2199 // // cout<<"i: "<<i<<" "<<"j: "<<j<<" "<<"n: "<<n<<endl;
2200 // vec_selectNum[i][j].push_back(n);
2201 // }
2202 //}
2203}
2204void MdcHoughFinder::FillRhoTheta(TH2D *h1 ,vector<double> vec_theta,vector<double> vec_rho,int numCross){
2205 for(int i=0;i<numCross;i++){
2206 double binwidth=M_PI/binX;
2207 double binhigh=2*t_maxP/binY;
2208 int binxNum=vec_theta[i]/binwidth+1;
2209 int binyNum=(vec_rho[i]+t_maxP)/binhigh+1;
2210 int count =h1->GetBinContent(binxNum,binyNum);
2211 count++;
2212 h1->SetBinContent(binxNum,binyNum,count);
2213 }
2214}
n
const Int_t n
Definition: DataBase/tau_mode.c:65
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
count
DOUBLE_PRECISION count[2]
Definition: EvtPhokharaDef.hh:68
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition: Gam4pikp.cxx:37
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
RecMdcTrackCol
ObjectVector< RecMdcTrack > RecMdcTrackCol
Definition: InstallArea/include/MdcRecEvent/MdcRecEvent/RecMdcTrack.h:79
v
**********Class see also m_nmax DOUBLE PRECISION m_amel DOUBLE PRECISION m_x2 DOUBLE PRECISION m_alfinv DOUBLE PRECISION m_Xenph INTEGER m_KeyWtm INTEGER m_idyfs DOUBLE PRECISION m_zini DOUBLE PRECISION m_q2 DOUBLE PRECISION m_Wt_KF DOUBLE PRECISION m_WtCut INTEGER m_KFfin *COMMON c_KarLud $ !Input CMS energy[GeV] $ !CMS energy after beam spread beam strahlung[GeV] $ !Beam energy spread[GeV] $ !z boost due to beam spread $ !electron beam mass *ff pair spectrum $ !minimum v
Definition: KarLud.h:35
M_PI
#define M_PI
Definition: TConstant.h:4
Helix::momentum
Hep3Vector momentum(double dPhi=0.) const
returns momentum vector after rotating angle dPhi in phi direction.
Definition: Reconstruction/TrackUtil/TrackUtil-00-00-08/src/Helix.cxx:184
Helix::dr
double dr(void) const
returns an element of parameters.
Definition: InstallArea/include/TrackUtil/TrackUtil/Helix.h:227
Helix::a
const HepVector & a(void) const
returns helix parameters.
Definition: InstallArea/include/TrackUtil/TrackUtil/Helix.h:263
Helix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition: InstallArea/include/TrackUtil/TrackUtil/Helix.h:209
MdcCalibFunSvc::getT0
double getT0(int layid, int cellid) const
Definition: MdcCalibFunSvc.cxx:507
MdcCalibFunSvc::getTimeWalk
double getTimeWalk(int layid, double Q) const
Definition: MdcCalibFunSvc.cxx:514
MdcDetector::Layer
const MdcLayer * Layer(unsigned id) const
Definition: InstallArea/include/MdcGeom/MdcGeom/MdcDetector.h:33
MdcDetector::instance
static MdcDetector * instance()
Definition: MdcDetector.cxx:21
MdcGeomSvc::Wire
const MdcGeoWire *const Wire(unsigned id)
Definition: MdcGeomSvc.cxx:768
MdcHit::setCalibSvc
void setCalibSvc(const MdcCalibFunSvc *calibSvc)
Definition: MdcHit.cxx:136
MdcHit::driftTime
double driftTime(double tof, double z) const
Definition: MdcHit.cxx:142
MdcHit::setCountPropTime
void setCountPropTime(const bool count)
Definition: InstallArea/include/MdcData/MdcData/MdcHit.h:86
MdcHit::driftDist
double driftDist(double, int, double, double, double) const
Definition: MdcHit.cxx:156
MdcHoughFinder::execute
StatusCode execute()
Definition: MdcHoughFinder.cxx:275
MdcHoughFinder::initialize
StatusCode initialize()
Definition: MdcHoughFinder.cxx:96
MdcHoughFinder::beginRun
StatusCode beginRun()
Definition: MdcHoughFinder.cxx:85
MdcHoughFinder::MdcHoughFinder
MdcHoughFinder(const std::string &name, ISvcLocator *pSvcLocator)
Definition: MdcHoughFinder.cxx:56
MdcHoughFinder::finalize
StatusCode finalize()
Definition: MdcHoughFinder.cxx:1601
MdcID::layer
static int layer(const Identifier &id)
Values of different levels (failure returns 0)
Definition: MdcID.cxx:49
MdcID::wire
static int wire(const Identifier &id)
Definition: MdcID.cxx:54
MdcTrack::storeTrack
void storeTrack(int trackId, RecMdcTrackCol *trackList, RecMdcHitCol *hitList, int tkStat)
Definition: MdcTrack.cxx:143
Pdt::readMCppTable
static void readMCppTable(std::string filenm)
RawData::identify
virtual Identifier identify() const
Definition: RawData.cxx:15
RawData::getChargeChannel
unsigned int getChargeChannel() const
Definition: RawData.cxx:45
TrkFit::nActive
virtual int nActive() const =0
TrkFit::helix
virtual TrkExchangePar helix(double fltL) const =0
TrkHitList::fit
TrkErrCode fit()
Definition: TrkHitList.cxx:59
TrkHitList::appendHot
TrkHitOnTrk * appendHot(const TrkHitOnTrk *theHot)
Definition: TrkHitList.cxx:91
TrkHitList::hotList
const TrkHotList & hotList() const
Definition: TrkHitList.cxx:51
TrkRecoTrk::print
virtual void print(std::ostream &) const
TrkRecoTrk::fitResult
const TrkFit * fitResult() const
Definition: TrkRecoTrk.cxx:387
TrkSimpleMaker::setFlipAndDrop
bool setFlipAndDrop(TrkRecoTrk &, bool allowFlips, bool allowDrops) const
TrkSimpleMaker::makeTrack
TrkRecoTrk * makeTrack(const TrkExchangePar &helix, const double chi2, const TrkContext &, double trackT0) const