00001 #include "MdcHoughFinder/MdcHoughFinder.h"
00002 #include "GaudiKernel/MsgStream.h"
00003 #include "GaudiKernel/AlgFactory.h"
00004 #include "GaudiKernel/ISvcLocator.h"
00005 #include "GaudiKernel/SmartDataPtr.h"
00006 #include "GaudiKernel/IDataProviderSvc.h"
00007 #include "GaudiKernel/PropertyMgr.h"
00008 #include "GaudiKernel/IService.h"
00009 #include "EventModel/EventHeader.h"
00010 #include "McTruth/DecayMode.h"
00011 #include "McTruth/McParticle.h"
00012 #include "McTruth/MdcMcHit.h"
00013 #include "RawDataProviderSvc/IRawDataProviderSvc.h"
00014 #include "RawDataProviderSvc/RawDataProviderSvc.h"
00015 #include "MdcRawEvent/MdcDigi.h"
00016 #include "Identifier/Identifier.h"
00017 #include "Identifier/MdcID.h"
00018 #include "MdcGeomSvc/IMdcGeomSvc.h"
00019 #include "MdcGeomSvc/MdcGeoWire.h"
00020 #include "MdcGeomSvc/MdcGeoLayer.h"
00021
00022 #include <vector>
00023 #include <iostream>
00024 #include <math.h>
00025 #include <map>
00026
00027 #include "EvTimeEvent/RecEsTime.h"
00028 #include "MdcCalibFunSvc/IMdcCalibFunSvc.h"
00029 #include "RawDataProviderSvc/MdcRawDataProvider.h"
00030 #include "TrkBase/TrkFit.h"
00031 #include "TrkBase/TrkHitList.h"
00032 #include "TrkBase/TrkExchangePar.h"
00033 #include "TrkFitter/TrkHelixMaker.h"
00034 #include "TrkFitter/TrkCircleMaker.h"
00035 #include "TrkFitter/TrkLineMaker.h"
00036 #include "TrkFitter/TrkHelixFitter.h"
00037 #include "MdcxReco/Mdcxprobab.h"
00038 #include "MdcData/MdcHit.h"
00039 #include "MdcData/MdcRecoHitOnTrack.h"
00040 #include "MdcTrkRecon/MdcTrack.h"
00041 #include "MdcPrintSvc/IMdcPrintSvc.h"
00042 #include "MdcGeom/EntranceAngle.h"
00043 #include "TrackUtil/Helix.h"
00044 #include "GaudiKernel/IPartPropSvc.h"
00045 #include "MdcRecoUtil/Pdt.h"
00046 #include "RawEvent/RawDataUtil.h"
00047
00048 #include "TTree.h"
00049 #include "TH2D.h"
00050
00051 using namespace std;
00052 using namespace Event;
00053
00055
00056 MdcHoughFinder::MdcHoughFinder(const std::string& name, ISvcLocator* pSvcLocator) :
00057 Algorithm(name, pSvcLocator)
00058 {
00059
00060 declareProperty("trackNumMc", m_trackNum_Mc_set=1);
00061 declareProperty("method", m_method=0);
00062 declareProperty("debug", m_debug = 0);
00063 declareProperty("data", m_data= 0);
00064 declareProperty("binx", m_binx= 100);
00065 declareProperty("biny", m_biny= 100);
00066 declareProperty("peakCellNum", m_peakCellNum);
00067 declareProperty("ndev", m_ndev= 3);
00068 declareProperty("f_pro", m_fpro= 0.8);
00069 declareProperty("f_hit_pro", m_hit_pro= 0.8);
00070
00071 declareProperty("pdtFile", m_pdtFile = "pdt.table");
00072 declareProperty("pickHits", m_pickHits = true);
00073 declareProperty("pid", m_pid = 0);
00074
00075 declareProperty("combineTracking",m_combineTracking = true);
00076 declareProperty("getDigiFlag", m_getDigiFlag = 0);
00077 declareProperty("maxMdcDigi", m_maxMdcDigi = 0);
00078 declareProperty("keepBadTdc", m_keepBadTdc = 0);
00079 declareProperty("dropHot", m_dropHot= 0);
00080 declareProperty("keepUnmatch", m_keepUnmatch = 0);
00081 declareProperty("minMdcDigi", m_minMdcDigi = 0);
00082
00083 }
00084
00085 StatusCode MdcHoughFinder::beginRun(){
00086
00087
00088 m_gm = MdcDetector::instance(0);
00089 if(NULL == m_gm) return StatusCode::FAILURE;
00090
00091
00092 return StatusCode::SUCCESS;
00093 }
00094
00095
00096 StatusCode MdcHoughFinder::initialize(){
00097
00098 MsgStream log(msgSvc(), name());
00099 log << MSG::INFO << "in initialize()" << endreq;
00100
00101 StatusCode sc;
00102
00103
00104
00105
00106 IPartPropSvc* p_PartPropSvc;
00107 static const bool CREATEIFNOTTHERE(true);
00108 sc = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
00109 if (!sc.isSuccess() || 0 == p_PartPropSvc) {
00110 log << MSG::ERROR << " Could not initialize PartPropSvc" << endreq;
00111 return sc;
00112 }
00113 m_particleTable = p_PartPropSvc->PDT();
00114
00115
00116 IRawDataProviderSvc* irawDataProviderSvc;
00117 sc = service ("RawDataProviderSvc", irawDataProviderSvc);
00118 m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*> (irawDataProviderSvc);
00119 if ( sc.isFailure() ){
00120 log << MSG::FATAL << name()<<" Could not load RawDataProviderSvc!" << endreq;
00121 return StatusCode::FAILURE;
00122
00123 }
00124
00125
00126 IMdcGeomSvc* imdcGeomSvc;
00127 sc = service ("MdcGeomSvc", imdcGeomSvc);
00128 m_mdcGeomSvc = dynamic_cast<MdcGeomSvc*> (imdcGeomSvc);
00129 if ( sc.isFailure() ){
00130 log << MSG::FATAL << "Could not load MdcGeoSvc!" << endreq;
00131 return StatusCode::FAILURE;
00132 }
00133
00134 sc = service ("MagneticFieldSvc",m_pIMF);
00135 if(sc!=StatusCode::SUCCESS) {
00136 log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
00137 }
00138 m_bfield = new BField(m_pIMF);
00139 log << MSG::INFO << "field z = "<<m_bfield->bFieldNominal()<< endreq;
00140
00141 m_context = new TrkContextEv(m_bfield);
00142
00143 Pdt::readMCppTable(m_pdtFile);
00144
00145
00146 IMdcCalibFunSvc* imdcCalibSvc;
00147 sc = service ("MdcCalibFunSvc", imdcCalibSvc);
00148 m_mdcCalibFunSvc = dynamic_cast<MdcCalibFunSvc*>(imdcCalibSvc);
00149 if ( sc.isFailure() ){
00150 log << MSG::FATAL << "Could not load MdcCalibFunSvc!" << endreq;
00151 return StatusCode::FAILURE;
00152 }
00153
00154
00155 NTuplePtr nt(ntupleSvc(), "MdcHoughFinder/hit");
00156 if ( nt ){
00157 ntuplehit = nt;
00158 } else {
00159 ntuplehit = ntupleSvc()->book("MdcHoughFinder/hit", CLID_ColumnWiseTuple, "hit");
00160 if(ntuplehit){
00161 ntuplehit->addItem("eventNum", m_eventNum);
00162 ntuplehit->addItem("runNum", m_runNum);
00163 ntuplehit->addItem("costaCut", m_cosCut);
00164
00165
00166
00167
00168
00169
00170
00171
00172 ntuplehit->addItem("nHitMc", m_nHit_Mc,0, 6796);
00173 ntuplehit->addItem("layerMc", m_nHit_Mc, m_layer_Mc);
00174 ntuplehit->addItem("cellMc", m_nHit_Mc, m_cell_Mc);
00175
00176
00177 ntuplehit->addItem("nHit", m_nHit,0, 6796);
00178 ntuplehit->addItem("layerNhit", 43, m_layerNhit);
00179
00180 ntuplehit->addItem("nCross", m_nCross,0, 125000);
00181 ntuplehit->addItem("rho", m_nCross, m_rho);
00182 ntuplehit->addItem("theta", m_nCross, m_theta);
00183
00184 ntuplehit->addItem("hitSignal", m_nHit, m_hitSignal);
00185 ntuplehit->addItem("layer", m_nHit, m_layer);
00186 ntuplehit->addItem("cell", m_nHit, m_cell);
00187 ntuplehit->addItem("x_east", m_nHit, m_x_east);
00188 ntuplehit->addItem("y_east", m_nHit, m_y_east);
00189 ntuplehit->addItem("z_east", m_nHit, m_z_east);
00190 ntuplehit->addItem("x_west", m_nHit, m_x_west);
00191 ntuplehit->addItem("y_west", m_nHit, m_y_west);
00192 ntuplehit->addItem("z_west", m_nHit, m_z_west);
00193 ntuplehit->addItem("x", m_nHit, m_x);
00194 ntuplehit->addItem("y", m_nHit, m_y);
00195 ntuplehit->addItem("z", m_nHit, m_z);
00196 ntuplehit->addItem("u", m_nHit, m_u);
00197 ntuplehit->addItem("v", m_nHit, m_v);
00198 ntuplehit->addItem("p", m_nHit, m_p);
00199 ntuplehit->addItem("slant", m_nHit, m_slant);
00200
00201
00202 ntuplehit->addItem("maxCount", m_maxCount);
00203
00204 ntuplehit->addItem("peakColumn", m_npeak,0,100);
00205 ntuplehit->addItem("peakWidth", m_npeak, m_peakWidth);
00206 ntuplehit->addItem("peakHigh", m_npeak, m_peakHigh);
00207 ntuplehit->addItem("peakArea", m_npeak, m_peakArea);
00208 ntuplehit->addItem("peakAreaLeast", m_areaLeast);
00209 ntuplehit->addItem("peakAreaLeastNum", m_areaLeastNum);
00210 ntuplehit->addItem("peakHitSelect", m_areaSelectHit);
00211 ntuplehit->addItem("peakHitSelectSignal", m_areaSelectHit_signal);
00212
00213
00214
00215 ntuplehit->addItem("circleCenterX", m_x_circle);
00216 ntuplehit->addItem("circleCenterY", m_y_circle);
00217 ntuplehit->addItem("circleR", m_r);
00218
00219
00220 ntuplehit->addItem("zStereoNum", m_zStereoNum,0,1000);
00221 ntuplehit->addItem("zStereo", m_zStereoNum, m_zStereo );
00222 ntuplehit->addItem("l", m_zStereoNum, m_l);
00223
00224 ntuplehit->addItem("trackNum_Mc", m_trackNum_Mc, 0,100);
00225 ntuplehit->addItem("trackNum", m_trackNum, 0,100);
00226 ntuplehit->addItem("d0", m_trackNum, m_d0);
00227 ntuplehit->addItem("phi0", m_trackNum, m_phi0);
00228 ntuplehit->addItem("omega", m_trackNum, m_omega);
00229 ntuplehit->addItem("z0", m_trackNum, m_z0);
00230 ntuplehit->addItem("tanl", m_trackNum, m_tanl);
00231
00232 ntuplehit->addItem("pt_rec", m_trackNum, m_pt);
00233 ntuplehit->addItem("pt2_rec", m_trackNum, m_pt2);
00234 ntuplehit->addItem("pz_rec", m_trackNum, m_pz);
00235 ntuplehit->addItem("p_rec", m_trackNum, m_pxyz);
00236
00237 ntuplehit->addItem("nHitSignal", m_nHitSignal);
00238 ntuplehit->addItem("nHitAxialSignal", m_nHitAxialSignal);
00239 ntuplehit->addItem("nHitSignal_select", m_trackNum, m_nHitSignal_select);
00240 ntuplehit->addItem("nHitAxialSignal_selcet", m_trackNum, m_nHitAxialSignal_select);
00241 ntuplehit->addItem("fitFailure", m_trackNum, m_nFitFailure);
00242 ntuplehit->addItem("nFit", m_trackNum, m_2d_nFit);
00243 ntuplehit->addItem("3dnFit", m_trackNum, m_3d_nFit);
00244 ntuplehit->addItem("nHitAxial", m_nHitAxial);
00245 ntuplehit->addItem("nFitSucess", m_nFitSucess);
00246
00247
00248
00249 ntuplehit->addItem("pidTruth", m_pidTruth);
00250 ntuplehit->addItem("costaTruth", m_costaTruth);
00251 ntuplehit->addItem("phiTruth", m_phi0Truth);
00252 ntuplehit->addItem("drTruth", m_drTruth);
00253 ntuplehit->addItem("dzTruth", m_dzTruth);
00254 ntuplehit->addItem("ptTruth", m_ptTruth);
00255 ntuplehit->addItem("pzTruth", m_pzTruth);
00256 ntuplehit->addItem("pTruth", m_pTruth);
00257 ntuplehit->addItem("qTruth", m_qTruth);
00258
00259
00260 } else { log << MSG::ERROR << "Cannot book tuple MdcHoughFinder/hit" << endmsg;
00261 return StatusCode::FAILURE;
00262 }
00263 }
00264
00265
00266 if(m_debug>2)TrkHelixFitter::m_debug = true;
00267
00268 return StatusCode::SUCCESS;
00269
00270
00271 }
00272
00273
00274
00275 StatusCode MdcHoughFinder::execute() {
00276
00277 MsgStream log(msgSvc(), name());
00278 log << MSG::INFO << "in execute()" << endreq;
00279
00280
00281 SmartDataPtr<RecEsTimeCol> evTimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
00282 if (!evTimeCol) {
00283 log << MSG::WARNING<< "Could not retrieve RecEsTimeCol , use t0=0" << endreq;
00284 m_bunchT0=0.;
00285 }else{
00286 RecEsTimeCol::iterator iter_evt = evTimeCol->begin();
00287 if (iter_evt != evTimeCol->end()){
00288 m_bunchT0 = (*iter_evt)->getTest()*1.e-9;
00289 }
00290 }
00291
00292
00293
00294
00295 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00296 if (!eventHeader) {
00297 log << MSG::FATAL << "Could not find Event Header" << endreq;
00298 return( StatusCode::FAILURE);
00299 }
00300
00301
00302 DataObject *aTrackCol;
00303 DataObject *aRecHitCol;
00304 SmartIF<IDataManagerSvc> dataManSvc(eventSvc());
00305 if(!m_combineTracking){
00306 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
00307 if(aTrackCol != NULL) {
00308 dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
00309 eventSvc()->unregisterObject("/Event/Recon/RecMdcTrackCol");
00310 }
00311 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
00312 if(aRecHitCol != NULL) {
00313 dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
00314 eventSvc()->unregisterObject("/Event/Recon/RecMdcHitCol");
00315 }
00316 }
00317
00318 RecMdcTrackCol* trackList;
00319 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
00320 if (aTrackCol) {
00321 trackList = dynamic_cast<RecMdcTrackCol*> (aTrackCol);
00322 }else{
00323 trackList = new RecMdcTrackCol;
00324 StatusCode sc = eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trackList);
00325 if(!sc.isSuccess()) {
00326 log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
00327 return StatusCode::FAILURE;
00328 }
00329 }
00330 int nTdsTk = (int) trackList->size();
00331 int nFitSucess = 0;
00332
00333 RecMdcHitCol* hitList;
00334 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
00335 if (aRecHitCol) {
00336 hitList = dynamic_cast<RecMdcHitCol*> (aRecHitCol);
00337 }else{
00338 hitList = new RecMdcHitCol;
00339 StatusCode sc = eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitList);
00340 if(!sc.isSuccess()) {
00341 log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
00342 return StatusCode::FAILURE;
00343 }
00344 }
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361 t_eventNum=eventHeader->eventNumber();
00362 t_runNum=eventHeader->runNumber();
00363 m_eventNum=t_eventNum;
00364 m_runNum=t_runNum;
00365 log << MSG::INFO << "MdcHoughFinder: retrieved event: " << eventHeader->eventNumber() << " run: " << eventHeader->runNumber() << endreq;
00366
00367
00368 int mc_particle=GetMcInfo();
00369
00370
00371 if(abs(m_costaTruth)<=0.83){m_cosCut=1;}
00372 else{m_cosCut=0;}
00373
00374
00375
00376
00377
00378
00381
00382 vector<double> vec_u;
00383 vector<double> vec_v;
00384 vector<double> vec_p;
00385 vector<double> vec_x_east;
00386 vector<double> vec_y_east;
00387 vector<double> vec_z_east;
00388 vector<double> vec_x_west;
00389 vector<double> vec_y_west;
00390 vector<double> vec_z_west;
00391 vector<double> vec_z_Mc;
00392 vector<double> vec_x;
00393 vector<double> vec_y;
00394 vector<double> vec_z;
00395 vector<int> vec_layer;
00396 vector<int> vec_wire;
00397 vector<int> vec_slant;
00398 vector<double> vec_zStereo;
00399 vector<double> l;
00400
00401 vector<int> vec_layer_Mc;
00402 vector<int> vec_wire_Mc;
00403 t_maxP=-999.;
00404 t_minP=999.;
00405
00406
00407 int t_nHit_Mc;
00408 int digiId_Mc=0;
00409 int nHitAxial_Mc=0;
00410
00411 SmartDataPtr<Event::MdcMcHitCol> mcMdcMcHitCol(eventSvc(),"/Event/MC/MdcMcHitCol");
00412 if (!mcMdcMcHitCol) {
00413 log << MSG::INFO << "Could not find MdcMcHit" << endreq;
00414 }else{
00415 Event::MdcMcHitCol::iterator iterMcHit = mcMdcMcHitCol->begin();
00416
00417 for(; iterMcHit != mcMdcMcHitCol->end(); ++iterMcHit,digiId_Mc++) {
00418 Identifier mdcid = (*iterMcHit)->identify();
00419 int layerId_Mc = MdcID::layer(mdcid);
00420 int wireId_Mc = MdcID::wire(mdcid);
00421
00422 vec_layer_Mc.push_back(layerId_Mc);
00423 vec_wire_Mc.push_back(wireId_Mc);
00424
00425
00426
00427
00428 m_layer_Mc[digiId_Mc]=layerId_Mc;
00429 m_cell_Mc[digiId_Mc]=wireId_Mc;
00430
00431 if(m_data==0){
00432 if ((layerId_Mc>=8&&layerId_Mc<=19)||(layerId_Mc>=36)) {
00433 nHitAxial_Mc++;}
00434 if(m_debug>0) {cout<<"("<<layerId_Mc<<","<<wireId_Mc<<")"<<"nHitAxial_Mc: "<<nHitAxial_Mc<<endl;}
00435 m_layer[digiId_Mc]=layerId_Mc;
00436 m_cell[digiId_Mc]=wireId_Mc;
00437 const MdcGeoWire *geowir =m_mdcGeomSvc->Wire(layerId_Mc,wireId_Mc);
00438 HepPoint3D eastP = geowir->Backward()/10.0;
00439 HepPoint3D westP = geowir->Forward() /10.0;
00440
00441 m_x_east[digiId_Mc]=eastP.x();
00442 m_y_east[digiId_Mc]=eastP.y();
00443 m_z_east[digiId_Mc]=eastP.z();
00444 m_x_west[digiId_Mc]=westP.x();
00445 m_y_west[digiId_Mc]=westP.y();
00446 m_z_west[digiId_Mc]=westP.z();
00447
00448 vec_x_east.push_back(eastP.x());
00449 vec_y_east.push_back(eastP.y());
00450 vec_z_east.push_back(eastP.z());
00451 vec_x_west.push_back(westP.x());
00452 vec_y_west.push_back(westP.y());
00453 vec_z_west.push_back(westP.z());
00454
00455 double x = (*iterMcHit)->getPositionX()/10.;
00456 double y = (*iterMcHit)->getPositionY()/10.;
00457 double z = (*iterMcHit)->getPositionZ()/10.;
00458
00459 double u=CFtrans(x,y);
00460 double v=CFtrans(y,x);
00461 double p=sqrt(u*u+v*v);
00462
00463 vec_x.push_back(x);
00464 vec_y.push_back(y);
00465 vec_z.push_back(z);
00466 vec_u.push_back(u);
00467 vec_v.push_back(v);
00468 vec_p.push_back(p);
00469
00470 m_x[digiId_Mc]=x;
00471 m_y[digiId_Mc]=y;
00472 m_z[digiId_Mc]=z;
00473 m_u[digiId_Mc]=u;
00474 m_v[digiId_Mc]=v;
00475 m_p[digiId_Mc]=p;
00476
00477 int layer= layerId_Mc;
00478 vec_slant.push_back(SlantId(layer));
00479 m_slant[digiId_Mc]=SlantId(layer);
00480 if (m_slant!=0)
00481 {cout<<"layer: "<<layerId_Mc<<"wire: "<<wireId_Mc<<"x_east: "<<m_x_east[digiId_Mc]<<"y_east: "<<m_y_east[digiId_Mc]<<endl;}
00482 }
00483 }
00484 }
00485
00486 t_nHit_Mc=digiId_Mc;
00487 m_nHit_Mc=digiId_Mc;
00488 if(m_data==0) {m_nHit=digiId_Mc;}
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505 vector<int> vec_hitSignal;
00506 vector<int> vec_track_index;
00507 int track_index_max=0;
00508 uint32_t getDigiFlag = 0;
00509 if(m_dropHot || m_combineTracking)getDigiFlag |= MdcRawDataProvider::b_dropHot;
00510 if(m_keepBadTdc) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
00511 if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
00512 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
00513
00514 cout<<"MdcDigiVec: "<<mdcDigiVec.size()<<endl;
00515 MdcDigiVec::iterator iter1 = mdcDigiVec.begin();
00516 int t_nHit;
00517 int digiId = 0;
00518 int nHitAxial = 0;
00519 int nHitLayer[43];
00520 int nHitSignal=0;
00521 int nHitAxialSignal=0;
00522
00523 for (;iter1 != mdcDigiVec.end(); iter1++, digiId++) {
00524 int track_index=(*iter1)->getTrackIndex();
00525
00526 if(track_index>=1000) track_index-=1000;
00527 vec_track_index.push_back(track_index);
00528 if(track_index>=0){
00529 nHitSignal++;
00530 m_hitSignal[digiId]=1;
00531 }
00532
00533 Identifier mdcId= (*iter1)->identify();
00534 int layerId = MdcID::layer(mdcId);
00535 int wireId = MdcID::wire(mdcId);
00536
00537
00538
00539
00540 if ((layerId>=8&&layerId<=19)||(layerId>=36)) {
00541 nHitAxial++;}
00542 if (((layerId>=8&&layerId<=19)||(layerId>=36))&&track_index>=0) {
00543 nHitAxialSignal++;}
00544
00545 vec_layer.push_back(layerId);
00546 vec_wire.push_back(wireId);
00547
00548 nHitLayer[layerId]++;
00549 m_layerNhit[layerId]=nHitLayer[layerId];
00550
00551 const MdcGeoWire *geowir =m_mdcGeomSvc->Wire(layerId,wireId);
00552 HepPoint3D eastP = geowir->Backward()/10.0;
00553 HepPoint3D westP = geowir->Forward() /10.0;
00554
00555 vec_x_east.push_back(eastP.x());
00556 vec_y_east.push_back(eastP.y());
00557 vec_z_east.push_back(eastP.z());
00558 vec_x_west.push_back(westP.x());
00559 vec_y_west.push_back(westP.y());
00560 vec_z_west.push_back(westP.z());
00561
00562 m_x_east[digiId]=eastP.x();
00563 m_y_east[digiId]=eastP.y();
00564 m_z_east[digiId]=eastP.z();
00565 m_x_west[digiId]=westP.x();
00566 m_y_west[digiId]=westP.y();
00567 m_z_west[digiId]=westP.z();
00568 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;}
00569
00570 double x =(eastP.x()+westP.x())/2.;
00571 double y =(eastP.y()+westP.y())/2.;
00572 vec_x.push_back((vec_x_east[digiId]+vec_x_west[digiId])/2);
00573 vec_y.push_back((vec_y_east[digiId]+vec_y_west[digiId])/2);
00574 m_x[digiId]=(vec_x_east[digiId]+vec_x_west[digiId])/2;
00575 m_y[digiId]=(vec_y_east[digiId]+vec_y_west[digiId])/2;
00576 double u=CFtrans(x,y);
00577 double v=CFtrans(y,x);
00578
00579
00580 vec_p.push_back(sqrt(u*u+v*v));
00581 vec_u.push_back(u);
00582 vec_v.push_back(v);
00583 m_u[digiId]=u;
00584 m_v[digiId]=v;
00585 m_p[digiId]=sqrt(u*u+v*v);
00586
00587 m_layer[digiId]=geowir->Layer();
00588 m_cell[digiId]=geowir->Cell();
00589 int layer= layerId;
00590 vec_slant.push_back(SlantId(layer));
00591 m_slant[digiId]=SlantId(layer);
00592
00593
00594
00595 }
00596 for(int i=0;i<digiId;i++){
00597 if(track_index_max<=vec_track_index[i]) track_index_max=vec_track_index[i];
00598 }
00599 track_index_max++;
00600 m_trackNum_Mc=track_index_max;
00601
00602 vector< vector <int> > mc_track_hit(track_index_max,vector<int>() );
00603 for(int i=0;i<track_index_max;i++){
00604 for(int j=0;j<digiId;j++){
00605 if(vec_track_index[j]==i) mc_track_hit[i].push_back(j);
00606 }
00607 }
00608
00609 track_index_max=m_trackNum_Mc_set;
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619 m_nHit=digiId;
00620 m_nHitAxial=nHitAxial;
00621
00622
00623 m_nHitSignal=nHitSignal;
00624 m_nHitAxialSignal=nHitAxialSignal;
00625
00626 cout<<"hit number: "<<digiId<<endl;
00627
00628
00629
00630 for(int i =0;i<m_nHit;i++){
00631 if(t_maxP<vec_p[i]) {t_maxP=vec_p[i];}
00632 }
00633
00634 for(int i =0;i<m_nHit;i++){
00635 if(t_minP>vec_p[i]) {t_minP=vec_p[i];}
00636 }
00637 t_maxP=t_maxP+0.01;
00638
00639
00640
00641
00642
00643
00644
00645
00646 int nhit;
00647 if(m_data==0) {
00648 nhit=m_nHit_Mc;}
00649 else{
00650 nhit=m_nHit;}
00651
00652 binX=m_binx;
00653 binY=m_biny;
00654 double binwidth=M_PI/binX;
00655 double binhigh=2*t_maxP/binY;
00656 TH2D *h1=new TH2D("h1","h1",binX,0,M_PI,binY,-t_maxP,t_maxP);
00657
00658
00659
00660 vector<double> vec_rho;
00661 vector<double> vec_theta;
00662 vector< vector<int> > vec_hitNum(2,vector<int>());
00663 int numCross=(int)(nhit*(nhit-1)/2);
00664 if(m_method==0){
00665 RhoTheta(numCross,nhit,vec_u,vec_v,vec_rho,vec_theta,vec_hitNum);
00666 FillRhoTheta(h1,vec_theta,vec_rho,numCross);
00667 }
00668
00669 vector< vector< vector<int> > > vec_selectNum(binX,vector< vector<int> >(binY,vector<int>() ) );
00670 if(m_method==1){
00671 FillHist(h1,vec_u,vec_v,nhit,vec_selectNum);
00672 }
00673
00674
00675 vector<bool> vec_truthHit(nhit,false);
00676 vector< vector <int> > countij(binX,vector <int> (binY,0) );
00677
00678 vector< vector< vector<int> > > vec_selectHit(binX,vector< vector<int> >(binY,vector<int>() ) );
00679 if(m_method==2){
00680 FillCells(h1,nhit,binX,binY,vec_u,vec_v,vec_layer,vec_wire,countij,vec_selectNum,vec_selectHit);
00681 }
00682
00683
00684
00685
00686
00688
00690
00691
00692
00693
00694
00695
00696
00697
00698 int max_count=0;
00699 int max_binx=1;
00700 int max_biny=1;
00701 for(int i=0;i<binX;i++){
00702 for(int j=0;j<binY;j++){
00703 int count=h1->GetBinContent(i+1,j+1);
00704 if(max_count<count) {
00705 max_count=count;
00706 max_binx=i+1;
00707 max_biny=j+1;
00708 }
00709 }
00710 }
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
00758
00759
00760
00761
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806 int count_use=0;
00807 double sum=0;
00808 double sum2=0;
00809 for(int i=0;i<binX;i++){
00810 for(int j=0;j<binY;j++){
00811 int count=h1->GetBinContent(i+1,j+1);
00812
00813 if(count!=0){
00814 count_use++;
00815 sum+=count;
00816 sum2+=count*count;
00817 }
00818 }
00819 }
00820 double f_n=m_ndev;
00821 cout<<"count_use: "<<count_use<<endl;
00822 double aver=sum/count_use;
00823 double dev=sqrt(sum2/count_use-aver*aver);
00824 int min_counts=(int)(aver+f_n*dev);
00825 cout<<"aver: "<<aver<<" "<<"dev: "<<dev<<" min: "<<min_counts<<endl;
00826
00827
00828 vector< vector <int> > hough_trans_CS(binX,vector <int> (binY,0) );
00829 vector< vector <int> > hough_trans_CS_peak(binX,vector <int> (binY,0) );
00830
00831
00832 for(int i=0;i<binX;i++){
00833 for(int j=0;j<binY;j++){
00834 hough_trans_CS[i][j]=h1->GetBinContent(i+1,j+1);
00835 }
00836 }
00837
00838 int peak_num=0;
00839 for (int r=0; r<binY; r++) {
00840 for (int a=0; a<binX; a++) {
00841 long max_around=0;
00842 for (int ar=a-1; ar<=a+1; ar++) {
00843 for (int rx=r-1; rx<=r+1; rx++) {
00844 int ax;
00845 if (ar<0) { ax=ar+binX; }
00846 else if (ar>=binX) { ax=ar-binX; }
00847 else { ax=ar; }
00848 if ( (ax!=a || rx!=r) && rx>=0 && rx<binY) {
00849 if (hough_trans_CS[ax][rx]>max_around) { max_around=hough_trans_CS[ax][rx]; }
00850 }
00851 }
00852 }
00853 if (hough_trans_CS[a][r]<=min_counts) { hough_trans_CS_peak[a][r]=0; }
00854 else if (hough_trans_CS[a][r]<max_around) { hough_trans_CS_peak[a][r]=0; }
00855 else if (hough_trans_CS[a][r]==max_around) { hough_trans_CS_peak[a][r]=1; }
00856 else if (hough_trans_CS[a][r]>max_around) { hough_trans_CS_peak[a][r]=2; }
00857 if(hough_trans_CS_peak[a][r]!=0) {
00858
00859 peak_num++;
00860 }
00861 }
00862 }
00863
00864
00865
00866 int peak_num_2=0;
00867 for (int r=0; r<binY; r++) {
00868 for (int a=0; a<binX; a++) {
00869 if (hough_trans_CS_peak[a][r]==1) {
00870 hough_trans_CS_peak[a][r]=2;
00871 for (int ar=a-1; ar<=a+1; ar++) {
00872 for (int rx=r-1; rx<=r+1; rx++) {
00873 int ax;
00874 if (ar<0) { ax=ar+binX; }
00875 else if (ar>=binX) { ax=ar-binX; }
00876 else { ax=ar; }
00877 if ( (ax!=a || rx!=r) && rx>=0 && rx<binY) {
00878 if (hough_trans_CS_peak[ax][rx]==1) { hough_trans_CS_peak[ax][rx]=0; }
00879 }
00880 }
00881 }
00882 }
00883 if(hough_trans_CS_peak[a][r]!=0) {
00884 peak_num_2++;
00885
00886
00887
00888
00889
00890 }
00891 }
00892 }
00893
00894
00895
00896 TH2D *h2 = new TH2D("h2","test2",binX,0,3.14,binY,-t_maxP,t_maxP);
00897 for(int i=0;i<binX;i++){
00898 for(int j=0;j<binY;j++){
00899 if(hough_trans_CS_peak[i][j]==2){
00900 h2->SetBinContent(i+1,j+1,hough_trans_CS[i][j]);
00901 }
00902 }
00903 }
00904
00905 vector<int> peakList1[3];
00906 for(int n=0;n<peak_num_2;n++){
00907 for (int r=0; r<binY; r++) {
00908 for (int a=0; a<binX; a++) {
00909 if (hough_trans_CS_peak[a][r]==2) {
00910 peakList1[0].push_back(a+1);
00911 peakList1[1].push_back(r+1);
00912 peakList1[2].push_back(hough_trans_CS[a][r]);
00913 }
00914 }
00915 }
00916 }
00917 cout<<"finish peak?"<<endl;
00918 if(m_method==0) cout<<"numCross: "<<numCross<<endl;
00919
00920 int npeak1=peak_num_2;
00921 int n_max;
00922 for (int na=0; na<npeak1-1; na++) {
00923 n_max=na;
00924 for (int nb=na+1; nb<npeak1; nb++) {
00925 if (peakList1[2][n_max]<peakList1[2][nb]) { n_max=nb; }
00926 }
00927 if (n_max!=na) {
00928 float swap[3];
00929 for (int i=0; i<3; i++) {
00930 swap[i]=peakList1[i][na];
00931 peakList1[i][na]=peakList1[i][n_max];
00932 peakList1[i][n_max]=swap[i];
00933 }
00934 }
00935 }
00936 cout<<"npeak1: "<<npeak1<<endl;
00937 for(int n=0;n<npeak1;n++){
00938 int bintheta=peakList1[0][n];
00939 int binrho=peakList1[1][n];
00940 int count=peakList1[2][n];
00941 for(int i=0;i<count;i++){
00942
00943
00944 }
00945
00946 }
00947
00948 typedef std::map<int, int> M2;
00949 typedef std::map<int, M2> M1;
00950 M2 peak_combine_num;
00951 M1 hit_combine;
00952
00953 int peak_track=0;
00954
00955 if(m_method==2){
00956
00957 int m=0;
00958 int n=0;
00959 int addnum=999;
00960 vector<bool> peaktrack(npeak1,true);
00961 while(addnum!=0)
00962 {
00963 addnum=0;
00964 double peak_cellNum[43];
00965 int peakX[43];
00966 int peakY[43];
00967 double bin_left[43];
00968 double bin_right[43];
00969 double bin_up[43];
00970 double bin_down[43];
00971 double area_left[43];
00972 double area_right[43];
00973
00974 double area_down[43];
00975 double area_up[43];
00976
00977 for(int iLayer=0; iLayer<m_gm->nLayer(); iLayer++){
00978 peak_cellNum[iLayer]=m_peakCellNum[iLayer];
00979 peakX[iLayer]=peakList1[0][n];
00980 peakY[iLayer]=peakList1[1][n];
00981 bin_left[iLayer]=peakX[iLayer]-peak_cellNum[iLayer];
00982 bin_right[iLayer]=peakX[iLayer]+peak_cellNum[iLayer];
00983 bin_up[iLayer]=peakY[iLayer]+peak_cellNum[iLayer];
00984 bin_down[iLayer]=peakY[iLayer]-peak_cellNum[iLayer];
00985 area_left[iLayer]=(bin_left[iLayer]-1)*binwidth;
00986 area_right[iLayer]=(bin_right[iLayer])*binwidth;
00987
00988 area_down[iLayer]=(bin_down[iLayer]-1)*binhigh-t_maxP;
00989 area_up[iLayer]=(bin_up[iLayer])*binhigh-t_maxP;
00990 }
00991 int count_peak=0;
00992 for(int k=0;k<nhit;k++){
00993 int layer=vec_layer[k];
00994 double lineLeft=vec_u[k]*cos(area_left[layer])+vec_v.at(k)*sin(area_right[layer]);
00995 double lineRight=vec_u[k]*cos(area_left[layer])+vec_v.at(k)*sin(area_right[layer]);
00996
00997 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)){
00998
00999
01000
01001
01002
01003
01004 hit_combine[m][count_peak]=k;
01005
01006 count_peak++;
01007 }
01008
01009 peak_combine_num[m]=count_peak;
01010
01011
01012 }
01013 for(int i=n+1;i<npeak1;i++){
01014 if(peaktrack[i]==false) continue;
01015 int peaktheta=peakList1[0][i];
01016 int peakrho=peakList1[1][i];
01017 int peakhitNum=peakList1[2][i];
01018 int count_same_hit=0;
01019 for(int j=0;j<peakhitNum;j++){
01020
01021 for(int k=0;k<peak_combine_num[m];k++){
01022 if(vec_selectNum[peaktheta-1][peakrho-1][j]==hit_combine[m][k]){
01023 count_same_hit++;
01024 }
01025 }
01026 }
01027 double f_hit=m_hit_pro;
01028 if(count_same_hit>f_hit*peakhitNum){
01029 peaktrack[i]=false;
01030 }
01031 }
01032 for(int i=n+1;i<npeak1;i++){
01033 if(peaktrack[i]==true){
01034 addnum=i;
01035 break;
01036 }
01037 }
01038 if(addnum!=0) m++;
01039 cout<<"peak_m: "<<m+1<<endl;
01040 n=n+addnum;
01041 }
01042
01043
01044 for(int i=0;i<npeak1;i++){
01045 cout<<"track"<<i<<": "<<"("<<peakList1[0][i]<<","<<peakList1[1][i]<<","<<peakList1[2][i]<<")"<<" "<<"truth: "<<peaktrack[i]<<endl;
01046 if(peaktrack[i]==true){
01047 peak_track++;
01048 }
01049 }
01050 for( int i=0;i<peak_track;i++){
01051 for(int j =0;j<peak_combine_num[i];j++){
01052 int hit_number=hit_combine[i][j];
01053 cout<<" peak "<<i<<" has select hits: "<<vec_layer[hit_number]<<" "<<vec_wire[hit_number]<<endl;
01054 }
01055 cout<<"has collect :"<<peak_combine_num[i]<<" hits "<<endl;
01056 }
01057 cout<<"event"<<t_eventNum<<": "<<"has found: "<<peak_track<<" track."<<endl;
01058 m_trackNum=peak_track;
01059
01060
01061
01062
01063
01064 vector< vector<int> > rec_mc_num(peak_track,vector<int>(track_index_max,0) );
01065 vector< vector<int> > rec_mc_num_axial(peak_track,vector<int>(track_index_max,0) );
01066 if(track_index_max!=peak_track) cout<<"not match track number !"<<endl;
01067 for(int mc_track_num=0;mc_track_num<track_index_max;mc_track_num++){
01068 for(int i=0;i<peak_track;i++){
01069 for(int j=0;j<peak_combine_num[i];j++){
01070 for(int k=0;k<mc_track_hit[mc_track_num].size();k++){
01071 if(hit_combine[i][j]==mc_track_hit[mc_track_num][k]){
01072 rec_mc_num[i][mc_track_num]++;
01073 int hit_num=mc_track_hit[mc_track_num][k];
01074 if(vec_slant[hit_num]==0) rec_mc_num_axial[i][mc_track_num]++;
01075 }
01076 }
01077 }
01078 }
01079 }
01080 vector<int> rec_mc(peak_track,999);
01081 for(int i=0;i<peak_track;i++){
01082 for(int mc_track_num=0;mc_track_num<track_index_max;mc_track_num++){
01083 if(rec_mc_num[i][mc_track_num]>0.5*peak_combine_num[i]) rec_mc[i]=mc_track_num;
01084 cout<<"trackrec: "<<i<<" trackmc: "<<mc_track_num<<" "<<rec_mc_num[i][mc_track_num]<<" "<<peak_combine_num[i]<<endl;
01085 }
01086 }
01087 for(int i=0;i<peak_track;i++){
01088 cout<<"rec :"<<i<<"belong to mc track: "<<rec_mc[i]<<endl;
01089 }
01090 for(int i=0;i<peak_track;i++){
01091 int mc_track_num=rec_mc[i];
01092 if(mc_track_num!=999) {
01093 cout<<"debug mc_track_num: "<<mc_track_num<<endl;
01094 m_nHitSignal_select[i]=rec_mc_num[i][mc_track_num];
01095 m_nHitAxialSignal_select[i]=rec_mc_num_axial[i][mc_track_num];
01096 cout<<"m_nHitSignal_select: "<<m_nHitSignal_select[i]<<endl;
01097 }
01098 }
01099
01100
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117
01118
01119 }
01120
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139
01140
01141
01142
01143
01144
01145
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161
01162
01163
01164
01165
01166
01167
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197
01198
01199
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212
01213
01214
01215
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226
01227
01228
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288
01289
01290
01291
01292
01293
01294
01295
01296
01297
01298
01299
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311
01312
01313
01314
01315
01316
01317
01318
01319
01320
01321
01322
01323
01324
01325
01326
01327
01328
01329
01330
01331
01332
01333
01334
01335
01336
01337
01338
01339
01340
01341
01342
01343
01344
01345 double d0=-999.;
01346 double phi0=-999.;
01347 double omega=-999.;
01348 double z0=0;
01349 double tanl=0;
01350
01351 for(track_fit=0;track_fit<peak_track;track_fit++){
01352 for(int i=0;i<nhit;i++){
01353 vec_truthHit[i]=false;
01354 }
01355 cout<<"track: "<<track_fit<<" has select: "<<peak_combine_num[track_fit]<<" hit ."<<endl;
01356 for(int i=0;i<peak_combine_num[track_fit];i++){
01357 int hitNum=hit_combine[track_fit][i];
01358 vec_truthHit[hitNum]=true;
01359 }
01360
01361 int nHitAxialSelect_temp=10;
01362 int leastSquare=LeastSquare(nhit,nHitAxialSelect_temp,vec_x,vec_y,vec_slant,vec_layer,vec_truthHit,d0,phi0,omega);
01363
01364 if(leastSquare==0){
01365
01366 TrkExchangePar tt(d0,phi0,omega,z0,tanl);
01367 TrkCircleMaker circlefactory;
01368 float chisum =0.;
01369 TrkRecoTrk* newTrack = circlefactory.makeTrack(tt, chisum, *m_context, m_bunchT0);
01370 bool permitFlips = true;
01371 bool lPickHits = m_pickHits;
01372 circlefactory.setFlipAndDrop(*newTrack, permitFlips, lPickHits);
01373
01374 int nDigi = digiToHots(newTrack,vec_truthHit);
01375 int nRecTk = 0;
01376
01377
01378 TrkHitList* qhits = newTrack->hits();
01379 TrkErrCode err=qhits->fit();
01380
01381 const TrkFit* newFit = newTrack->fitResult();
01382 bool fitSuccess = false;
01383
01384 if (!newFit || (newFit->nActive()<3)) {
01385 if(m_debug>0){
01386 cout << "evt "<<t_eventNum<<" global fit failed ";
01387 if(newFit) cout <<" nAct "<<newFit->nActive();
01388 cout<<endl;
01389 }
01390
01391 }else{
01392 nRecTk = 1;
01393 fitSuccess = true;
01394 m_nEvtSuccess++;
01395
01396 if(m_debug>0) newTrack->print(std::cout);
01397
01398 }
01399
01400 vector<int> nfit2d(peak_track,0);
01401 if(m_debug>0) cout<<" hot list:"<<endl;
01402 TrkHotList::hot_iterator hotIter= qhits->hotList().begin();
01403 while (hotIter!=qhits->hotList().end()) {
01404 nfit2d[track_fit]++;
01405 cout <<"(" <<((MdcHit*)(hotIter->hit()))->layernumber()
01406 <<","<<((MdcHit*)(hotIter->hit()))->wirenumber()
01407 <<":"<<hotIter->isActive()<<") ";
01408 hotIter++;
01409 }
01410
01411 m_2d_nFit[track_fit]=nfit2d[track_fit];
01412
01413
01414 double x_cirtemp;
01415 double y_cirtemp;
01416 double r_temp;
01417 if(err.failure()==0){
01418
01419 TrkExchangePar par = newFit->helix(0.);
01420 d0=par.d0();
01421 phi0=par.phi0();
01422 omega=par.omega();
01423
01424
01425
01426
01427 r_temp=1./-par.omega();
01428 x_cirtemp=(r_temp-par.d0())*cos(par.phi0()-M_PI/2.);
01429 y_cirtemp=(r_temp-par.d0())*sin(par.phi0()-M_PI/2.);
01430
01431 m_pt[track_fit]=r_temp/333.567;
01432 cout<<"pt_rec: "<<m_pt[track_fit]<<endl;
01433 }
01434 else{
01435 m_nFitFailure[track_fit]=2;
01436 }
01437
01438
01439 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);
01440
01441
01442
01443 m_zStereoNum=z_stereoNum;
01444 for(int i=0;i<z_stereoNum;i++){
01445
01446 m_zStereoNum=z_stereoNum;
01447 m_zStereo[i]=vec_zStereo[i];
01448
01449 m_l[i]=l[i];
01450 if(m_debug>0) cout<<" l: "<<m_l[i]<<" "<<"z: "<<vec_zStereo[i]<<endl;
01451 }
01452
01453
01454 Linefit(z_stereoNum,vec_zStereo,l,tanl,z0);
01455
01456 cout<<"tanl: "<<tanl<<endl;
01457 cout<<"z0: "<<z0<<endl;
01458 z0=dz_mc;
01459 tanl=tanl_mc;
01460 cout<<"tanltruth: "<<tanl<<endl;
01461 cout<<"z0truth: "<<z0<<endl;
01462
01463
01464 TrkExchangePar tt2(d0,phi0,omega,z0,tanl);
01465 TrkHelixMaker helixfactory;
01466 chisum =0.;
01467 TrkRecoTrk* newTrack2 = helixfactory.makeTrack(tt2, chisum, *m_context, m_bunchT0);
01468 permitFlips = true;
01469 lPickHits = m_pickHits;
01470 helixfactory.setFlipAndDrop(*newTrack2, permitFlips, lPickHits);
01471
01472 int nDigi2 = digiToHots2(newTrack2,vec_truthHit);
01473
01474
01475
01476 TrkHitList* qhits2 = newTrack2->hits();
01477 TrkErrCode err2=qhits2->fit();
01478
01479 const TrkFit* newFit2 = newTrack2->fitResult();
01480
01481
01482
01483 if (!newFit2 || (newFit2->nActive()<5)) {
01484
01485 if(m_debug>0){
01486 cout << "evt "<<t_eventNum<<" global fit failed ";
01487 if(newFit2) cout <<" nAct "<<newFit2->nActive();
01488 cout<<endl;
01489 }
01490
01491 }else{
01492
01493 nFitSucess++;
01494 MdcTrack* mdcTrack = new MdcTrack(newTrack2);
01495 int tkStat = 4;
01496 int tkId = nTdsTk + track_fit;
01497 mdcTrack->storeTrack(tkId, trackList, hitList, tkStat);
01498 nRecTk = 1;
01499 fitSuccess = true;
01500 m_nEvtSuccess++;
01501
01502 if(m_debug>0) newTrack2->print(std::cout);
01503
01504 }
01505
01506 vector<int> nfit3d(peak_track,0);
01507 if(m_debug>0) cout<<" hot list:"<<endl;
01508 TrkHotList::hot_iterator hotIter2= qhits2->hotList().begin();
01509 while (hotIter2!=qhits2->hotList().end()) {
01510 nfit3d[track_fit]++;
01511 cout <<"(" <<((MdcHit*)(hotIter2->hit()))->layernumber()
01512 <<","<<((MdcHit*)(hotIter2->hit()))->wirenumber()
01513 <<":"<<hotIter2->isActive()<<") ";
01514 hotIter2++;
01515 }
01516
01517 m_3d_nFit[track_fit]=nfit3d[track_fit];
01518
01519
01520 if(err2.failure()==0){
01521
01522 TrkExchangePar par2 = newFit2->helix(0.);
01523
01524
01525
01526
01527
01528
01529
01530 r_temp=1./-par2.omega();
01531 x_cirtemp=(r_temp-par2.d0())*cos(par2.phi0()-M_PI/2.);
01532 y_cirtemp=(r_temp-par2.d0())*sin(par2.phi0()-M_PI/2.);
01533
01534 m_d0[track_fit]=par2.d0();
01535 m_phi0[track_fit]=par2.phi0();
01536 m_omega[track_fit]=par2.omega();
01537 m_z0[track_fit]=par2.z0();
01538 m_tanl[track_fit]=par2.tanDip();
01539
01540 double pxy=r_temp/333.567;
01541 double pz=pxy*par2.tanDip();
01542 double pxyz=sqrt(pxy*pxy+pz*pz);
01543 m_pt2[track_fit]=pxy;
01544 m_pz[track_fit]=pxy*par2.tanDip();
01545 m_pxyz[track_fit]=pxyz;
01546 cout<<"track"<<track_fit<<": "<<"pt_rec2: "<<m_pt2[track_fit]<<endl;
01547 cout<<"eventNum: "<<t_eventNum<<" "<<"track"<<track_fit<<": "<<"pz_rec: "<<m_pz[track_fit]<<endl;
01548 }
01549 else{
01550 m_nFitFailure[track_fit]=3;
01551 }
01552
01553
01554
01555
01556
01557
01558
01559
01560
01561
01562
01563
01564
01565
01566
01567
01568
01569 }
01570
01571
01572
01573
01574
01575
01576
01577
01578
01579
01580
01581
01582
01583
01584
01585
01586 }
01587 m_nFitSucess=nFitSucess;
01588 cout<<" Hough finder find "<<nFitSucess<<" tot "<<nTdsTk<< endl;
01589
01590 delete h1;
01591 delete h2;
01592 cout<<"break: ????"<<endl;
01593 ntuplehit->write();
01594 cout<<"finish event "<<t_eventNum<<endl;
01595 cout<<endl;
01596 cout<<endl;
01597 return StatusCode::SUCCESS;
01598 }
01599
01600
01601 StatusCode MdcHoughFinder::finalize() {
01602 MsgStream log(msgSvc(), name());
01603 log << MSG::INFO << "in finalize()" << endreq;
01604 return StatusCode::SUCCESS;
01605 }
01606 int MdcHoughFinder::digiToHots(TrkRecoTrk* newTrack,vector<bool> vec_truthHit){
01607 TrkHitList* trkHitList = newTrack->hits();
01608
01609
01610
01611
01612
01613
01614
01615
01616
01617
01618
01619
01620
01621
01622 int digiId=0;
01623 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec();
01624 MdcDigiVec::iterator iter1 = mdcDigiVec.begin();
01625 for (;iter1 != mdcDigiVec.end(); iter1++, digiId++) {
01626
01627
01628 if(vec_truthHit[digiId]==false) continue;
01629 const MdcDigi* aDigi = *iter1;
01630 MdcHit *hit = new MdcHit(aDigi, m_gm);
01631 Identifier id = aDigi->identify();
01632 int layer = MdcID::layer(id);
01633 int wire = MdcID::wire(id);
01634
01635 if ((layer>=0&&layer<=7)||(layer>=20&&layer<=35)) {continue;}
01636
01637
01638
01639 if(m_debug>0) std::cout<<" ("<<layer<<","<<wire<<") "<<std::endl;
01640 hit->setCalibSvc(m_mdcCalibFunSvc);
01641 hit->setCountPropTime(true);
01642
01643
01644
01645
01646 int newAmbig = 0;
01647
01648
01649 MdcRecoHitOnTrack temp(*hit, newAmbig, m_bunchT0*1.e9);
01650 MdcHitOnTrack* newhot = &temp;
01651
01652 double fltLen = m_gm->Layer(layer)->rMid();
01653 newhot->setFltLen(fltLen);
01654
01655
01656
01657 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;
01658
01659 if(hit->driftDist(m_bunchT0,0)>1.2) continue;
01660 trkHitList->appendHot(newhot);
01661
01662 }
01663 if(m_debug>1) std::cout<<std::endl;
01664 return trkHitList->nHit();
01665 }
01666
01667 double MdcHoughFinder::CFtrans(double x,double y){
01668 return x/(x*x+y*y);
01669 }
01670
01671
01672
01673
01674
01675
01676
01677
01678 int MdcHoughFinder::SlantId(int layer){
01679 if ((layer>=0&&layer<=7)||(layer>=20&&layer<=35)) {
01680 if ((layer>=0&&layer<=3)||(layer>=20&&layer<=23)||(layer>=28&&layer<=31)){
01681
01682 return -1;
01683 }else{
01684 return 1;
01685 }
01686 } else{
01687 return 0;
01688 }
01689 }
01690
01691
01692 int MdcHoughFinder::GetMcInfo(){
01693
01694 MsgStream log(msgSvc(), name());
01695 StatusCode sc;
01696 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
01697 if (!mcParticleCol) {
01698 log << MSG::ERROR << "Could not find McParticle" << endreq;
01699 return -999;
01700 }
01701 int itk = 0;
01702 int t_qTruth = 0;
01703 int t_pidTruth = -999;
01704 int t_McTkId = -999;
01705 t_nTrkMC=0;
01706 Helix* mchelix;
01707 McParticleCol::iterator iter_mc = mcParticleCol->begin();
01708 for (;iter_mc != mcParticleCol->end(); iter_mc++ ) {
01709 if(!(*iter_mc)->primaryParticle()) continue;
01710 t_pidTruth = (*iter_mc)->particleProperty();
01711
01712 if(m_debug>2) cout<< "tk "<<itk<< " pid="<< t_pidTruth<< endl;
01713 if((m_pid!=0) && (t_pidTruth != m_pid)){ continue; }
01714
01715
01716
01717 if(itk>=10) break;
01718
01719 int pid = t_pidTruth;
01720 if( pid == 0 ) {
01721 log << MSG::WARNING << "Wrong particle id" <<endreq;
01722 continue;
01723 }else{
01724 IPartPropSvc* p_PartPropSvc;
01725 static const bool CREATEIFNOTTHERE(true);
01726 StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
01727 if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) {
01728 std::cout<< " Could not initialize Particle Properties Service" << std::endl;
01729 }
01730 HepPDT::ParticleDataTable* p_particleTable = p_PartPropSvc->PDT();
01731 std::string name;
01732 if( p_particleTable->particle(pid) ){
01733 name = p_particleTable->particle(pid)->name();
01734 t_qTruth = p_particleTable->particle(pid)->charge();
01735 }else if( p_particleTable->particle(-pid) ){
01736 name = "anti " + p_particleTable->particle(-pid)->name();
01737 t_qTruth = (-1)*p_particleTable->particle(-pid)->charge();
01738 }
01739 if(m_debug>2) std::cout << " particle "<<name<<" charge= "<<t_qTruth<<std::endl;
01740 }
01741
01742 t_McTkId = (*iter_mc)->trackIndex();
01743
01744 double px = (*iter_mc)->initialFourMomentum().px();
01745 double py = (*iter_mc)->initialFourMomentum().py();
01746 double pz = (*iter_mc)->initialFourMomentum().pz();
01747
01748 m_pzTruth=pz;
01749
01750 mchelix = new Helix((*iter_mc)->initialPosition().v(), (*iter_mc)->initialFourMomentum().v(), t_qTruth);
01751 mchelix->pivot( HepPoint3D(0.,0.,0.) );
01752
01753 if(m_debug>0){
01754 std::cout<<"Truth tk "<<itk<<" pid "<<pid<<" charge "<<t_qTruth
01755 << " helix "<< mchelix->a()<<" p "<<mchelix->momentum(0.)<<std::endl;
01756 }
01757 t_nTrkMC++;
01758 itk++;
01759 }
01760 if(t_nTrkMC!=1) {
01761 std::cout<<"WARNING run "<<t_runNum<<" evt "<<t_eventNum<<" not single event. nTrkMc="<<t_nTrkMC<<std::endl;
01762 return -999;
01763 }else{
01764 if(m_debug>2) std::cout<<"nTrkMc=1"<<std::endl;
01765 }
01766
01767
01768 m_pidTruth = t_pidTruth;
01769 m_costaTruth = mchelix->cosTheta();
01770 m_phi0Truth = mchelix->phi0();
01771 m_drTruth = mchelix->dr();
01772 m_dzTruth = mchelix->dz();
01773 m_ptTruth = mchelix->pt();
01774 m_pTruth = mchelix->momentum(0.).mag();
01775 m_qTruth = t_qTruth;
01776 dz_mc =mchelix->dz();
01777 tanl_mc =mchelix->tanl();
01778 return 1;
01779 }
01780
01781
01782 int 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){
01783
01784
01785
01786
01787
01788 double x_sum=0;
01789 double y_sum=0;
01790 double x2_sum=0;
01791 double y2_sum=0;
01792 double x3_sum=0;
01793 double y3_sum=0;
01794 double x2y_sum=0;
01795 double xy2_sum=0;
01796 double xy_sum=0;
01797 double a1=0;
01798 double a2=0;
01799 double b1=0;
01800 double b2=0;
01801 double c1=0;
01802 double c2=0;
01803 double x_aver,y_aver,r2;
01804 for(int i=0;i<n;i++){
01805 if(vec_truthHit[i]==false) continue;
01806 if(vec_slant[i]!=0) continue;
01807
01808 x_sum=x_sum+vec_x[i];
01809 y_sum=y_sum+vec_y[i];
01810 x2_sum=x2_sum+vec_x[i]*vec_x[i];
01811 y2_sum=y2_sum+vec_y[i]*vec_y[i];
01812 x3_sum=x3_sum+vec_x[i]*vec_x[i]*vec_x[i];
01813 y3_sum=y3_sum+vec_y[i]*vec_y[i]*vec_y[i];
01814 x2y_sum=x2y_sum+vec_x[i]*vec_x[i]*vec_y[i];
01815 xy2_sum=xy2_sum+vec_x[i]*vec_y[i]*vec_y[i];
01816 xy_sum=xy_sum+vec_x[i]*vec_y[i];
01817 }
01818 a1=x2_sum-x_sum*x_sum/n;
01819 a2=xy_sum-x_sum*y_sum/n;
01820 b1=a2;
01821 b2=y2_sum-y_sum*y_sum/n;
01822 c1=(x3_sum+xy2_sum-x_sum*(x2_sum+y2_sum)/n)/2.;
01823 c2=(y3_sum+x2y_sum-y_sum*(x2_sum+y2_sum)/n)/2.;
01824
01825 x_aver=x_sum/n;
01826 y_aver=y_sum/n;
01827
01828 m_x_circle =(b1*c2-b2*c1)/(b1*b1-a1*b2);
01829 m_y_circle =(b1*c1-a1*c2)/(b1*b1-a1*b2);
01830
01831 double x_cirtemp =(b1*c2-b2*c1)/(b1*b1-a1*b2);
01832 double y_cirtemp =(b1*c1-a1*c2)/(b1*b1-a1*b2);
01833
01834 r2=(x2_sum+y2_sum-2*x_cirtemp*x_sum-2*y_cirtemp*y_sum)/n+x_cirtemp*x_cirtemp+y_cirtemp*y_cirtemp;
01835
01836 m_r=sqrt(r2);
01837
01838 double r_temp=sqrt(r2);
01839 double d0_temp=sqrt(x_cirtemp*x_cirtemp+y_cirtemp*y_cirtemp)-r_temp;
01840
01841
01842
01843
01844 double pt_temp=r_temp/333.567;
01845 cout<<"pt_temp: "<<pt_temp<<endl;
01846
01847
01848 d0=-d0_temp;
01849
01850
01851 phi0=atan2(y_cirtemp,x_cirtemp)+M_PI/2.;
01852 omega=-1./r_temp;
01853 double z0=0.;
01854 double tanl=0.;
01855
01856
01857
01858
01859
01860
01861
01862
01863 if(m_debug<0)std::cout<<" before global fit Helix PatPar :"<<d0<<","<<phi0<<","<<omega<<","<<z0<<","<<tanl<< std::endl;
01864 return 0;
01865
01866 }
01867
01868 int 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){
01869 double x_stereo;
01870 double y_stereo;
01871
01872 int stereoId=0;
01873
01874
01875 int nDropDelta=0;
01876 for(int i=0;i<n;i++){
01877 if(vec_truthHit[i]==false) continue;
01878 double k;
01879 double b;
01880 if(vec_slant[i]!=0){
01881 if(vec_x_east[i]!=vec_x_west[i]){
01882 k=(vec_y_east[i]-vec_y_west[i])/(vec_x_east[i]-vec_x_west[i]);
01883 b=-k*vec_x_east[i]+vec_y_east[i];
01884
01885 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);
01886
01887
01888
01889
01890 if(delta<0) {
01891 nDropDelta++;
01892 continue;
01893 }
01894
01895 double x1=(2*x_cirtemp+2*k*y_cirtemp-2*k*b+sqrt(delta))/(2*(k*k+1));
01896 double x2=(2*x_cirtemp+2*k*y_cirtemp-2*k*b-sqrt(delta))/(2*(k*k+1));
01897 double x_middle=(vec_x_east[i]+vec_x_west[i])/2.;
01898
01899 if(abs(x1-x_middle)>abs(x2-x_middle)) {x_stereo=x2;}
01900 else {x_stereo=x1;}
01901 y_stereo=k*x_stereo+b;
01902 }
01903 else{
01904
01905
01906 x_stereo=vec_x_east[i];
01907 double y1=sqrt(r_temp*r_temp-(x_cirtemp-x_stereo)*(x_cirtemp-x_stereo))+y_cirtemp;
01908 double y2=y_cirtemp-sqrt(r_temp*r_temp-(x_cirtemp-x_stereo)*(x_cirtemp-x_stereo));
01909 double y_middle=(vec_y_east[i]+vec_y_west[i])/2.;
01910 if(abs(y1-y_middle)>abs(y2-y_middle)) {y_stereo=y2;}
01911 else{y_stereo=y1;}
01912 }
01913
01914
01915 double theta_temp;
01916 if(x_cirtemp==0||x_stereo-x_cirtemp==0){
01917 theta_temp=0;
01918 }
01919 else{
01920 theta_temp=M_PI-atan2(y_stereo-y_cirtemp,x_stereo-x_cirtemp)+atan2(y_cirtemp,x_cirtemp);
01921 if(theta_temp>2*M_PI){
01922 theta_temp=theta_temp-2*M_PI;
01923 }
01924 if(theta_temp<0){
01925 theta_temp=theta_temp+2*M_PI;
01926 }
01927 }
01928
01929
01930
01931
01932
01933
01934
01935
01936 l.push_back(r_temp*theta_temp);
01937
01938
01939
01940
01941 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]));
01942 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]));
01943 vec_zStereo.push_back(vec_z_west[i]-(vec_z_west[i]-vec_z_east[i])*d1/d2);
01944
01945
01946
01947
01948
01949
01950
01951
01952
01953
01954
01955
01956
01957
01958
01959
01960
01961
01962 stereoId=stereoId+1;
01963
01964 }
01965 else {k=b=0;}
01966 }
01967 cout<<"nDropDelta: "<<nDropDelta<<endl;
01968
01969
01970
01971
01972
01973
01974
01975
01976
01977
01978
01979
01980
01981 return stereoId;
01982 }
01983
01984
01985 void MdcHoughFinder::Linefit(int z_stereoNum,vector<double> vec_zStereo,vector<double> l,double& tanl,double& z0){
01986
01987 double l_sum=0;
01988 double z_sum=0;
01989 double l2_sum=0;
01990 double z2_sum=0;
01991 double lz_sum=0;
01992 for(int i=0;i<z_stereoNum;i++){
01993
01994 l_sum=l_sum+l[i];
01995 z_sum=z_sum+vec_zStereo[i];
01996 l2_sum=l2_sum+l[i]*l[i];
01997
01998 z2_sum=z2_sum+vec_zStereo[i]*vec_zStereo[i];
01999 lz_sum=lz_sum+l[i]*vec_zStereo[i];
02000 }
02001 double b=(l2_sum*z_sum-lz_sum*l_sum)/(z_stereoNum*l2_sum-l_sum);
02002 double k=(z_stereoNum*lz_sum-l_sum*z_sum)/(z_stereoNum*l2_sum-l_sum*l_sum);
02003 z0=b;
02004 tanl=k;
02005 cout<<k<<" "<<b<<endl;
02006 }
02007
02008
02009
02010 int MdcHoughFinder::digiToHots2(TrkRecoTrk* newTrack,vector<bool> vec_truthHit){
02011 TrkHitList* trkHitList = newTrack->hits();
02012
02013
02014
02015
02016
02017
02018
02019
02020
02021
02022
02023
02024
02025
02026 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec();
02027 MdcDigiVec::iterator iter1 = mdcDigiVec.begin();
02028 int digiId=0;
02029 for (;iter1 != mdcDigiVec.end(); iter1++, digiId++) {
02030
02031
02032 if(vec_truthHit[digiId]==false) continue;
02033 const MdcDigi* aDigi = *iter1;
02034 MdcHit *hit = new MdcHit(aDigi, m_gm);
02035 Identifier id = aDigi->identify();
02036 int layer = MdcID::layer(id);
02037 int wire = MdcID::wire(id);
02038
02039
02040
02041
02042
02043 if(m_debug>0) std::cout<<" ("<<layer<<","<<wire<<") "<<std::endl;
02044 hit->setCalibSvc(m_mdcCalibFunSvc);
02045 hit->setCountPropTime(true);
02046
02047
02048
02049
02050 int newAmbig = 0;
02051
02052
02053 MdcRecoHitOnTrack temp(*hit, newAmbig, m_bunchT0*1.e9);
02054 MdcHitOnTrack* newhot = &temp;
02055
02056 double fltLen = m_gm->Layer(layer)->rMid();
02057 newhot->setFltLen(fltLen);
02058
02059
02060
02061 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;
02062
02063 if(hit->driftDist(m_bunchT0,0)>1.2) continue;
02064 trkHitList->appendHot(newhot);
02065
02066 }
02067 if(m_debug>1) std::cout<<std::endl;
02068 return trkHitList->nHit();
02069 }
02070 void MdcHoughFinder:: Multi_array(int binX,int binY){
02071 int **count=new int*[binX];
02072 for(int i=0;i<binX;i++){
02073 count[i]=new int [binY];
02074 }
02075
02076 int ***selectNum=new int**[binX];
02077 for(int i=0;i<binX;i++){
02078 selectNum[i]=new int* [binY];
02079 for(int j=0;j<binY;j++){
02080
02081 selectNum[i][j]=new int [count[i][j]];
02082 }
02083 }
02084 }
02085 void 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){
02086 double binWidth=M_PI/binX;
02087 double binHigh=2*t_maxP/binY;
02088
02089
02090
02091
02092
02093
02094
02095
02096
02097
02098 for(int i=0;i<binX;i++){
02099 for(int j=0;j<binY;j++){
02100 int count=0;
02101 double binLeft=i*binWidth;
02102 double binRight=(i+1)*binWidth;
02103 double binDown=-t_maxP+j*binHigh;
02104 double binUp=-t_maxP+(j+1)*binHigh;
02105 double binMid=(binLeft+binRight)/2;
02106 for(int k=0;k<n;k++){
02107 double lineLeft=vec_u[k]*cos(binLeft)+vec_v.at(k)*sin(binRight);
02108 double lineRight=vec_u[k]*cos(binLeft)+vec_v.at(k)*sin(binRight);
02109 double lineMid=vec_u[k]*cos(binMid)+vec_v[k]*sin(binMid);
02110 if(((lineLeft<binUp)&&(lineLeft>binDown)) || ((lineRight<binUp)&&(lineRight>binDown)) || (((lineLeft-binUp)*(binDown-lineRight)>0))){
02111 vec_selectNum[i][j].push_back(k);
02112 int layerId=vec_layer[k];
02113 int cellId=vec_wire[k];
02114 int hit_temp=(int)layerId*1000+cellId;
02115 vec_selectHit[i][j].push_back(hit_temp);
02116 count++;
02117 }
02118 }
02119 countij[i][j]=count;
02120
02121
02122
02123
02124
02125
02126
02127 if(m_debug>0) cout<<"countij: "<<countij[i][j]<<endl;
02128 }
02129 }
02130
02131
02132
02133
02134
02135
02136
02137
02138 for(int i=0;i<binX;i++){
02139 for(int j=0;j<binY;j++){
02140 h1->SetBinContent(i+1,j+1,countij[i][j]);
02141 }
02142 }
02143 }
02144 void 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){
02145 for(int i=0;i<nhit;i++){
02146 for(int j=i+1;j<nhit;j++){
02147 double k,b,x_cross,y_cross;
02148 double rho_temp,theta_temp;
02149 k=(vec_v[i]-vec_v[j])/(vec_u[i]-vec_u[j]);
02150 b=vec_v[i]-k*vec_u[i];
02151 x_cross=-b/(k+1/k);
02152 y_cross=b/(k*k+1);
02153 rho_temp=sqrt(x_cross*x_cross+y_cross*y_cross);
02154 theta_temp=atan2(y_cross,x_cross);
02155 if(theta_temp<0) {
02156 theta_temp=theta_temp+M_PI;
02157 rho_temp=-rho_temp;
02158 }
02159 vec_rho.push_back(rho_temp);
02160 vec_theta.push_back(theta_temp);
02161 vec_hitNum[0].push_back(i);
02162 vec_hitNum[1].push_back(j);
02163
02164 }
02165 }
02166
02167 m_nCross=numCross;
02168
02169 for(int i=0;i<numCross;i++){
02170 m_rho[i]=vec_rho[i];
02171 m_theta[i]=vec_theta[i];
02172
02173 }
02174 }
02175 void MdcHoughFinder::FillHist(TH2D *h1,vector<double> vec_u,vector<double> vec_v,int nhit,vector< vector < vector< int > > > vec_selectNum){
02176 for(int n=0;n<nhit;n++){
02177 for(int i=0;i<binX;i++){
02178 double binwidth=M_PI/binX;
02179 double binhigh=2*t_maxP/binY;
02180 double binMid=(i+0.5)*binwidth;
02181 double rho=vec_u[n]*cos(binMid)+vec_v[n]*sin(binMid);
02182 int j=(int)((rho+t_maxP)/binhigh);
02183 int count=0;
02184 count=h1->GetBinContent(i+1,j+1);
02185 count++;
02186 h1->SetBinContent(i+1,j+1,count);
02187 }
02188 }
02189
02190
02191
02192
02193
02194
02195
02196
02197
02198
02199
02200
02201
02202
02203 }
02204 void MdcHoughFinder::FillRhoTheta(TH2D *h1 ,vector<double> vec_theta,vector<double> vec_rho,int numCross){
02205 for(int i=0;i<numCross;i++){
02206 double binwidth=M_PI/binX;
02207 double binhigh=2*t_maxP/binY;
02208 int binxNum=vec_theta[i]/binwidth+1;
02209 int binyNum=(vec_rho[i]+t_maxP)/binhigh+1;
02210 int count =h1->GetBinContent(binxNum,binyNum);
02211 count++;
02212 h1->SetBinContent(binxNum,binyNum,count);
02213 }
02214 }