00001 #include "MdcHoughFinder/HoughValidUpdate.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 "MdcRawEvent/MdcDigi.h"
00013 #include "Identifier/Identifier.h"
00014 #include "Identifier/MdcID.h"
00015
00016
00017 #include "EvTimeEvent/RecEsTime.h"
00018 #include "RawDataProviderSvc/MdcRawDataProvider.h"
00019 #include "McTruth/MdcMcHit.h"
00020 #include "MdcRecEvent/RecMdcTrack.h"
00021 #include "MdcRecEvent/RecMdcHit.h"
00022 #include "TrkBase/TrkFit.h"
00023 #include "TrkBase/TrkHitList.h"
00024 #include "TrkBase/TrkExchangePar.h"
00025 #include "TrkFitter/TrkHelixMaker.h"
00026 #include "TrkFitter/TrkCircleMaker.h"
00027 #include "TrkFitter/TrkLineMaker.h"
00028 #include "TrkFitter/TrkHelixFitter.h"
00029 #include "MdcxReco/Mdcxprobab.h"
00030 #include "MdcData/MdcHit.h"
00031 #include "MdcData/MdcRecoHitOnTrack.h"
00032 #include "MdcTrkRecon/MdcTrack.h"
00033 #include "MdcGeom/EntranceAngle.h"
00034 #include "TrackUtil/Helix.h"
00035 #include "GaudiKernel/IPartPropSvc.h"
00036 #include "MdcRecoUtil/Pdt.h"
00037 #include "RawEvent/RawDataUtil.h"
00038 #include "TrkBase/TrkErrCode.h"
00039 #include "MdcPrintSvc/MdcPrintSvc.h"
00040
00041 #include "TTree.h"
00042 #include "TF1.h"
00043 #include "TGraph.h"
00044 #include "iomanip"
00045
00046 using namespace std;
00047 using namespace Event;
00048
00050
00051 HoughValidUpdate::HoughValidUpdate(const std::string& name, ISvcLocator* pSvcLocator) :
00052 Algorithm(name, pSvcLocator)
00053 {
00054
00055 declareProperty("drawTuple", m_drawTuple=0);
00056 declareProperty("useMcInfo", m_useMcInfo=0);
00057 declareProperty("method", m_method=0);
00058 declareProperty("debug", m_debug = 0);
00059 declareProperty("data", m_data= 0);
00060 declareProperty("binx", m_binx= 100);
00061 declareProperty("biny", m_biny= 200);
00062 declareProperty("maxrho", m_maxrho= 0.3);
00063 declareProperty("peakCellNumX", m_peakCellNumX);
00064 declareProperty("peakCellNumY", m_peakCellNumY);
00065 declareProperty("ndev1", m_ndev1= 1);
00066 declareProperty("ndev2", m_ndev2= 5);
00067 declareProperty("f_hit_pro", m_hit_pro= 0.4);
00068 declareProperty("pdtFile", m_pdtFile = "pdt.table");
00069 declareProperty("pickHits", m_pickHits = true);
00070 declareProperty("pid", m_pid = 0);
00071 declareProperty("combineTracking",m_combineTracking = false);
00072 declareProperty("getDigiFlag", m_getDigiFlag = 0);
00073 declareProperty("maxMdcDigi", m_maxMdcDigi = 0);
00074 declareProperty("keepBadTdc", m_keepBadTdc = 0);
00075 declareProperty("dropHot", m_dropHot= 0);
00076 declareProperty("keepUnmatch", m_keepUnmatch = 0);
00077 declareProperty("minMdcDigi", m_minMdcDigi = 0);
00078 declareProperty("setnpeak_fit", m_setpeak_fit= -1);
00079 declareProperty("Q", m_q= 0);
00080 declareProperty("eventcut", m_eventcut= -1);
00081 declareProperty("rhocut", m_rhocut= -1.);
00082 declareProperty("mcpar", m_mcpar= 0);
00083 declareProperty("fillTruth", m_fillTruth= 0);
00084 declareProperty("removeBadTrack", m_removeBadTrack= false);
00085 declareProperty("dropTrkDrCut", m_dropTrkDrCut= 1.);
00086 declareProperty("dropTrkDzCut", m_dropTrkDzCut= 10.);
00087 declareProperty("dropTrkPtCut", m_dropTrkPtCut= 99999.);
00088 declareProperty("dropTrkChi2Cut", m_dropTrkChi2Cut = 99999.);
00089 declareProperty("dropTrkChi2NdfCut", m_dropTrkChi2NdfCut = 99999.);
00090 declareProperty("qualityFactor", m_qualityFactor= 3.);
00091 declareProperty("dropTrkNhitCut", m_dropTrkNhitCut= 9);
00092 }
00093
00094 StatusCode HoughValidUpdate::beginRun(){
00095
00096
00097 m_gm = MdcDetector::instance(0);
00098 if(NULL == m_gm) return StatusCode::FAILURE;
00099
00100 return StatusCode::SUCCESS;
00101 }
00102
00103
00104 StatusCode HoughValidUpdate::initialize(){
00105
00106 MsgStream log(msgSvc(), name());
00107 log << MSG::INFO << "in initialize()" << endreq;
00108
00109 StatusCode sc;
00110
00111
00112 IPartPropSvc* p_PartPropSvc;
00113 static const bool CREATEIFNOTTHERE(true);
00114 sc = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
00115 if (!sc.isSuccess() || 0 == p_PartPropSvc) {
00116 log << MSG::ERROR << " Could not initialize PartPropSvc" << endreq;
00117 return sc;
00118 }
00119 m_particleTable = p_PartPropSvc->PDT();
00120
00121
00122 IRawDataProviderSvc* irawDataProviderSvc;
00123 sc = service ("RawDataProviderSvc", irawDataProviderSvc);
00124 m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*> (irawDataProviderSvc);
00125 if ( sc.isFailure() ){
00126 log << MSG::FATAL << name()<<" Could not load RawDataProviderSvc!" << endreq;
00127 return StatusCode::FAILURE;
00128
00129 }
00130
00131
00132 IMdcGeomSvc* imdcGeomSvc;
00133 sc = service ("MdcGeomSvc", imdcGeomSvc);
00134 m_mdcGeomSvc = dynamic_cast<MdcGeomSvc*> (imdcGeomSvc);
00135 if ( sc.isFailure() ){
00136 log << MSG::FATAL << "Could not load MdcGeoSvc!" << endreq;
00137 return StatusCode::FAILURE;
00138 }
00139
00140
00141 IMdcPrintSvc* iMdcPrintSvc;
00142 sc = service ("MdcPrintSvc", iMdcPrintSvc);
00143 m_mdcPrintSvc= dynamic_cast<MdcPrintSvc*> (iMdcPrintSvc);
00144 if ( sc.isFailure() ){
00145 log << MSG::FATAL << "Could not load MdcPrintSvc!" << endreq;
00146 return StatusCode::FAILURE;
00147 }
00148
00149 sc = service ("MagneticFieldSvc",m_pIMF);
00150 if(sc!=StatusCode::SUCCESS) {
00151 log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
00152 }
00153 m_bfield = new BField(m_pIMF);
00154 log << MSG::INFO << "field z = "<<m_bfield->bFieldNominal()<< endreq;
00155
00156 m_context = new TrkContextEv(m_bfield);
00157
00158 Pdt::readMCppTable(m_pdtFile);
00159
00160
00161 IMdcCalibFunSvc* imdcCalibSvc;
00162 sc = service ("MdcCalibFunSvc", imdcCalibSvc);
00163 m_mdcCalibFunSvc = dynamic_cast<MdcCalibFunSvc*>(imdcCalibSvc);
00164 if ( sc.isFailure() ){
00165 log << MSG::FATAL << "Could not load MdcCalibFunSvc!" << endreq;
00166 return StatusCode::FAILURE;
00167 }
00168
00169
00170 if(m_drawTuple){
00171 NTuplePtr nt1(ntupleSvc(), "HoughValidUpdate/hit");
00172 if ( nt1 ){
00173 ntuplehit = nt1;
00174 } else {
00175 ntuplehit = ntupleSvc()->book("HoughValidUpdate/hit", CLID_ColumnWiseTuple, "hit");
00176 if(ntuplehit){
00177 ntuplehit->addItem("eventNum", m_eventNum);
00178 ntuplehit->addItem("runNum", m_runNum);
00179 ntuplehit->addItem("nHitMc", m_nHit_Mc,0, 6796);
00180 ntuplehit->addItem("nHitAxialMc", m_nHitAxial_Mc);
00181 ntuplehit->addItem("layerMc", m_nHit_Mc, m_layer_Mc);
00182 ntuplehit->addItem("cellMc", m_nHit_Mc, m_cell_Mc);
00183
00184 ntuplehit->addItem("nTrkMC", m_nTrkMC,0,10);
00185 ntuplehit->addItem("pidTruth", m_nTrkMC,m_pidTruth);
00186 ntuplehit->addItem("costaTruth", m_nTrkMC,m_costaTruth);
00187 ntuplehit->addItem("ptTruth", m_nTrkMC,m_ptTruth);
00188 ntuplehit->addItem("pzTruth", m_nTrkMC,m_pzTruth);
00189 ntuplehit->addItem("pTruth", m_nTrkMC,m_pTruth);
00190 ntuplehit->addItem("qTruth", m_nTrkMC,m_qTruth);
00191 ntuplehit->addItem("drTruth", m_nTrkMC,m_drTruth);
00192 ntuplehit->addItem("phiTruth", m_nTrkMC,m_phi0Truth);
00193 ntuplehit->addItem("omegaTruth", m_nTrkMC,m_omegaTruth);
00194 ntuplehit->addItem("dzTruth", m_nTrkMC,m_dzTruth);
00195 ntuplehit->addItem("tanlTruth", m_nTrkMC,m_tanl_mc);
00196
00197 ntuplehit->addItem("nHit", m_nHit,0, 6796);
00198 ntuplehit->addItem("nHitDigi", m_nHitDigi);
00199 ntuplehit->addItem("nHitAxial", m_nHitAxial);
00200 ntuplehit->addItem("nHitStereo", m_nHitStereo);
00201 ntuplehit->addItem("layer", m_nHit, m_layer);
00202 ntuplehit->addItem("cell", m_nHit, m_cell);
00203 ntuplehit->addItem("x_east", m_nHit, m_x_east);
00204 ntuplehit->addItem("y_east", m_nHit, m_y_east);
00205 ntuplehit->addItem("z_east", m_nHit, m_z_east);
00206 ntuplehit->addItem("x_west", m_nHit, m_x_west);
00207 ntuplehit->addItem("y_west", m_nHit, m_y_west);
00208 ntuplehit->addItem("z_west", m_nHit, m_z_west);
00209 ntuplehit->addItem("x", m_nHit, m_x);
00210 ntuplehit->addItem("y", m_nHit, m_y);
00211 ntuplehit->addItem("z", m_nHit, m_z);
00212 ntuplehit->addItem("u", m_nHit, m_u);
00213 ntuplehit->addItem("v", m_nHit, m_v);
00214 ntuplehit->addItem("p", m_nHit, m_p);
00215 ntuplehit->addItem("slant", m_nHit, m_slant);
00216 ntuplehit->addItem("layerNhit", 43, m_layerNhit);
00217 ntuplehit->addItem("layerMax", m_layer_max);
00218 ntuplehit->addItem("l_truth", m_nHit, m_ltruth);
00219 ntuplehit->addItem("z_truth", m_nHit, m_ztruth);
00220 ntuplehit->addItem("z_postruth", m_nHit, m_postruth);
00221 ntuplehit->addItem("digi_truth", m_nHit, m_digi_truth);
00222 ntuplehit->addItem("e_delta", m_e_delta);
00223
00224 ntuplehit->addItem("nCross", m_nCross, 0, 125000);
00225 ntuplehit->addItem("rho", m_nCross, m_rho);
00226 ntuplehit->addItem("theta", m_nCross, m_theta);
00227 ntuplehit->addItem("xybinNum", m_xybinNum, 0,100000);
00228 ntuplehit->addItem("xybin", m_xybinNum, m_xybin);
00229 ntuplehit->addItem("xsigma", m_xsigma);
00230 ntuplehit->addItem("ysigma", m_ysigma);
00231
00232 ntuplehit->addItem("npeak_1", m_npeak_1);
00233 ntuplehit->addItem("npeak_2", m_npeak_2);
00234 ntuplehit->addItem("trackFailure", m_trackFailure);
00235 ntuplehit->addItem("npeak_after_tracking", m_npeak_after_tracking, 0,1000);
00236 ntuplehit->addItem("nHitSelect", m_npeak_after_tracking, m_nHitSelect);
00237 ntuplehit->addItem("nHitAxialSelect", m_npeak_after_tracking, m_nHitAxialSelect);
00238 ntuplehit->addItem("nHitStereoSelect", m_npeak_after_tracking, m_nHitStereoSelect);
00239 ntuplehit->addItem("naxiallose", m_npeak_after_tracking, m_axiallose);
00240 ntuplehit->addItem("nstereolose", m_npeak_after_tracking, m_stereolose);
00241
00242 ntuplehit->addItem("trackNum_Mc", m_trackNum_Mc, 0,100);
00243 ntuplehit->addItem("trackNum", m_trackNum, 0,100);
00244 ntuplehit->addItem("trackNum_fit", m_trackNum_fit, 0,100);
00245 ntuplehit->addItem("trackNum_tds", m_trackNum_tds, 0,100);
00246 ntuplehit->addItem("circleCenterX", m_trackNum, m_x_circle);
00247 ntuplehit->addItem("circleCenterY", m_trackNum, m_y_circle);
00248 ntuplehit->addItem("circleR", m_trackNum, m_r);
00249 ntuplehit->addItem("chis_pt", m_trackNum, m_chis_pt);
00250 ntuplehit->addItem("pt_rec", m_trackNum, m_pt);
00251
00252 ntuplehit->addItem("nHitStereo_use",m_nHitStereo_use,0,1000);
00253 ntuplehit->addItem("l_Calcu_Left", m_nHitStereo_use, m_lCalcuLeft);
00254 ntuplehit->addItem("l_Calcu_Right", m_nHitStereo_use, m_lCalcuRight);
00255 ntuplehit->addItem("z_Calcu_Left", m_nHitStereo_use, m_zCalcuLeft);
00256 ntuplehit->addItem("z_Calcu_Right", m_nHitStereo_use, m_zCalcuRight);
00257 ntuplehit->addItem("z_Calcu", m_nHitStereo_use, m_zCalcu);
00258 ntuplehit->addItem("l_Calcu", m_nHitStereo_use, m_lCalcu);
00259 ntuplehit->addItem("z_delta", m_nHitStereo_use, m_delta_z);
00260 ntuplehit->addItem("l_delta", m_nHitStereo_use, m_delta_l);
00261 ntuplehit->addItem("amb", m_nHitStereo_use, m_amb);
00262 ntuplehit->addItem("amb_no_match", m_ambig_no_match);
00263 ntuplehit->addItem("amb_no_match_pro", m_ambig_no_match_propotion);
00264 ntuplehit->addItem("tanl_Cal", m_tanl_Cal);
00265 ntuplehit->addItem("z0_Cal", m_z0_Cal);
00266
00267 ntuplehit->addItem("pt2_rec_final", m_trackNum_fit, m_pt2_rec_final);
00268 ntuplehit->addItem("p_rec_final", m_trackNum_fit, m_p_rec_final);
00269 ntuplehit->addItem("d0_final", m_trackNum_fit, m_d0_final);
00270 ntuplehit->addItem("phi0_final", m_trackNum_fit, m_phi0_final);
00271 ntuplehit->addItem("omega_final", m_trackNum_fit, m_omega_final);
00272 ntuplehit->addItem("z0_final", m_trackNum_fit, m_z0_final);
00273 ntuplehit->addItem("tanl_final", m_trackNum_fit, m_tanl_final);
00274 ntuplehit->addItem("chis_3d_final", m_trackNum_fit, m_chis_3d_final);
00275 ntuplehit->addItem("Q_final", m_trackNum_fit, m_Q);
00276
00277 ntuplehit->addItem("pt2_rec_tds", m_trackNum_tds, m_pt2_rec_tds);
00278 ntuplehit->addItem("p_rec_tds", m_trackNum_tds, m_p_rec_tds);
00279 ntuplehit->addItem("d0_tds", m_trackNum_tds, m_d0_tds);
00280 ntuplehit->addItem("phi0_tds", m_trackNum_tds, m_phi0_tds);
00281 ntuplehit->addItem("omega_tds", m_trackNum_tds, m_omega_tds);
00282 ntuplehit->addItem("z0_tds", m_trackNum_tds, m_z0_tds);
00283 ntuplehit->addItem("tanl_tds", m_trackNum_tds, m_tanl_tds);
00284 ntuplehit->addItem("Q_tds", m_trackNum_tds, m_Q_tds);
00285
00286
00287 ntuplehit->addItem("nhitdis0", m_nhitdis0, 0,2000);
00288 ntuplehit->addItem("nhitdis1", m_nhitdis1, 0,2000);
00289 ntuplehit->addItem("nhitdis2", m_nhitdis2, 0,2000);
00290 ntuplehit->addItem("hitdis0", m_nhitdis0, m_hit_dis0);
00291 ntuplehit->addItem("hitdis1", m_nhitdis1, m_hit_dis1);
00292 ntuplehit->addItem("hitdis2", m_nhitdis2, m_hit_dis2);
00293 ntuplehit->addItem("hitdis1_3d", m_nhitdis1, m_hit_dis1_3d);
00294 ntuplehit->addItem("hitdis2_3d", m_nhitdis2, m_hit_dis2_3d);
00295 ntuplehit->addItem("track_use_intrk1", m_track_use_intrk1);
00296 ntuplehit->addItem("noise_intrk1", m_noise_intrk1);
00297 ntuplehit->addItem("decay", m_decay);
00298 ntuplehit->addItem("outputTrk", m_outputTrk, 0,100);
00299 ntuplehit->addItem("hitOnTrk", m_outputTrk, m_hitOnTrk);
00300
00301
00302 } else { log << MSG::ERROR << "Cannot book tuple HoughValidUpdate/hough" << endmsg;
00303 return StatusCode::FAILURE;
00304 }
00305 }
00306 }
00307
00308 if(m_debug>2) TrkHelixFitter::m_debug = true;
00309
00310 return StatusCode::SUCCESS;
00311
00312
00313 }
00314
00315
00316
00317 StatusCode HoughValidUpdate::execute() {
00318
00319 cout.precision(6);
00320 MsgStream log(msgSvc(), name());
00321 log << MSG::INFO << "in execute()" << endreq;
00322
00323
00324 SmartDataPtr<RecEsTimeCol> evTimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
00325 if (!evTimeCol) {
00326 log << MSG::WARNING<< "Could not retrieve RecEsTimeCol , use t0=0" << endreq;
00327 m_bunchT0=0.;
00328 }else{
00329 RecEsTimeCol::iterator iter_evt = evTimeCol->begin();
00330 if (iter_evt != evTimeCol->end()){
00331 m_bunchT0 = (*iter_evt)->getTest()*1.e-9;
00332 }
00333 }
00334
00335
00336 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00337 if (!eventHeader) {
00338 log << MSG::FATAL << "Could not find Event Header" << endreq;
00339 return( StatusCode::FAILURE);
00340 }
00341
00342 DataObject *aTrackCol;
00343 DataObject *aRecHitCol;
00344 SmartIF<IDataManagerSvc> dataManSvc(eventSvc());
00345 if(!m_combineTracking){
00346 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
00347 if(aTrackCol != NULL) {
00348 dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
00349 eventSvc()->unregisterObject("/Event/Recon/RecMdcTrackCol");
00350 }
00351 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
00352 if(aRecHitCol != NULL) {
00353 dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
00354 eventSvc()->unregisterObject("/Event/Recon/RecMdcHitCol");
00355 }
00356 }
00357
00358 RecMdcTrackCol* trackList;
00359 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
00360 if (aTrackCol) {
00361 trackList = dynamic_cast<RecMdcTrackCol*> (aTrackCol);
00362 }else{
00363 trackList = new RecMdcTrackCol;
00364 StatusCode sc = eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trackList);
00365 if(!sc.isSuccess()) {
00366 log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
00367 return StatusCode::FAILURE;
00368 }
00369 }
00370
00371 RecMdcHitCol* hitList;
00372 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
00373 if (aRecHitCol) {
00374 hitList = dynamic_cast<RecMdcHitCol*> (aRecHitCol);
00375 }else{
00376 hitList = new RecMdcHitCol;
00377 StatusCode sc = eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitList);
00378 if(!sc.isSuccess()) {
00379 log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
00380 return StatusCode::FAILURE;
00381 }
00382 }
00383
00384
00385 if(m_removeBadTrack) {
00386 vector<RecMdcTrack*> vec_trackList;
00387 if( m_debug>0 ) cout<<"PATTSF collect "<<trackList->size()<<" track. "<<endl;
00388 if(trackList->size()!=0){
00389 RecMdcTrackCol::iterator iter_pat = trackList->begin();
00390 for(;iter_pat!=trackList->end();iter_pat++){
00391 double d0=(*iter_pat)->helix(0);
00392 double kap = (*iter_pat)->helix(2);
00393 double pt = 0.00001;
00394 if(fabs(kap)>0.00001) pt = fabs(1./kap);
00395 double dz=(*iter_pat)->helix(4);
00396 double chi2=(*iter_pat)->chi2();
00397 if( m_debug>0) cout<<"d0: "<<d0<<" z0: "<<dz<<" pt "<<pt<<" chi2 "<<chi2<<endl;
00398 if(fabs(d0)>m_dropTrkDrCut || fabs(dz)>m_dropTrkDzCut || chi2>m_dropTrkChi2Cut || pt>m_dropTrkPtCut) {
00399 vec_trackList.push_back(*iter_pat);
00400
00401 HitRefVec rechits = (*iter_pat)->getVecHits();
00402 HitRefVec::iterator hotIter = rechits.begin();
00403 while (hotIter!=rechits.end()) {
00404 if(m_debug>0) cout <<"remove hit " << (*hotIter)->getId() <<endl;
00405 hitList->remove(*hotIter);
00406 hotIter++;
00407 }
00408 trackList->remove(*iter_pat);
00409 iter_pat--;
00410 }
00411
00412 }
00413 } else {
00414 if(m_debug>0) cout<<" PATTSF find 0 track. "<<endl;
00415 }
00416 }
00417
00418 int nTdsTk=trackList->size();
00419
00420
00421 t_eventNum=eventHeader->eventNumber();
00422 t_runNum=eventHeader->runNumber();
00423 if( m_drawTuple ){
00424 m_eventNum=t_eventNum;
00425 m_runNum=t_runNum;
00426 }
00427 if(m_debug>0) cout<<" begin event :"<<t_eventNum<<endl;
00428 if(t_eventNum<=m_eventcut) {
00429 cout<<" eventNum cut "<<t_eventNum<<endl;
00430 return StatusCode::SUCCESS;
00431 }
00432 log << MSG::INFO << "HoughValidUpdate: retrieved event: " << eventHeader->eventNumber() << " run: " << eventHeader->runNumber() << endreq;
00433
00434 vector<double> vec_u;
00435 vector<double> vec_v;
00436 vector<double> vec_p;
00437 vector<double> vec_x_east;
00438 vector<double> vec_y_east;
00439 vector<double> vec_z_east;
00440 vector<double> vec_x_west;
00441 vector<double> vec_y_west;
00442 vector<double> vec_z_west;
00443 vector<double> vec_z_Mc;
00444 vector<double> vec_x;
00445 vector<double> vec_y;
00446 vector<double> vec_z;
00447 vector<int> vec_layer;
00448 vector<int> vec_wire;
00449 vector<int> vec_slant;
00450 vector<double> vec_zStereo;
00451 vector<double> l;
00452
00453 vector<int> vec_layer_Mc;
00454 vector<int> vec_wire_Mc;
00455 vector<int> vec_hit_Mc;
00456 vector<double> vec_ztruth_Mc;
00457 vector<double> vec_flt_truth_mc;
00458 vector<double> vec_pos_flag;
00459 vector<double> vec_phit_MC;
00460
00461
00462 int mc_particle=GetMcInfo();
00463 if(m_debug>2) cout<<"MC particle : "<<mc_particle<<endl;
00464
00465
00466 if(m_fillTruth ==1 && m_useMcInfo) {
00467 int digiId_Mc=0;
00468 int nHitAxial_Mc=0;
00469 SmartDataPtr<Event::MdcMcHitCol> mcMdcMcHitCol(eventSvc(),"/Event/MC/MdcMcHitCol");
00470 if (!mcMdcMcHitCol) {
00471 log << MSG::INFO << "Could not find MdcMcHit" << endreq;
00472 }else{
00473 Event::MdcMcHitCol::iterator iterMcHit = mcMdcMcHitCol->begin();
00474
00475 for(; iterMcHit != mcMdcMcHitCol->end(); ++iterMcHit,digiId_Mc++) {
00476 Identifier mdcid = (*iterMcHit)->identify();
00477 int layerId_Mc = MdcID::layer(mdcid);
00478 int wireId_Mc = MdcID::wire(mdcid);
00479 double hittemp=layerId_Mc*1000+wireId_Mc;
00480 double mc_hit_distance =(*iterMcHit)->getDriftDistance();
00481 double z_mc = (*iterMcHit)->getPositionZ()/10.;
00482 double flt_truth=(z_mc-dz_mc)/tanl_mc;
00483 double pos_flag=(*iterMcHit)->getPositionFlag();
00484
00485
00486 if( (*iterMcHit)->getPositionFlag()==-999 ){
00487
00488 double px= (*iterMcHit)->getPositionX();
00489 double py= (*iterMcHit)->getPositionY();
00490 double pz= (*iterMcHit)->getPositionZ();
00491 double pxyz = sqrt(px*px+py*py+pz*pz);
00492 vec_phit_MC.push_back(pxyz);
00493
00494
00495
00496
00497
00498
00499 }
00500
00501
00502 vec_layer_Mc.push_back(layerId_Mc);
00503 vec_wire_Mc.push_back(wireId_Mc);
00504 vec_hit_Mc.push_back(hittemp);
00505 vec_ztruth_Mc.push_back(z_mc);
00506 vec_flt_truth_mc.push_back(flt_truth);
00507 vec_pos_flag.push_back(pos_flag);
00508
00509 if(m_debug>5) {
00510 cout<<"(" <<layerId_Mc<<","<<wireId_Mc<<"):"<<endl;
00511 cout<<" hit_distance : "<<mc_hit_distance<<endl;
00512 cout<<" position flag : "<<pos_flag<<endl;
00513 cout<<" hit_z_mc : "<<z_mc<<endl;
00514 cout<<" flt_truth : "<<flt_truth<<endl;
00515 }
00516
00517 if(m_data==0){
00518 if ((layerId_Mc>=8&&layerId_Mc<=19)||(layerId_Mc>=36)) {
00519 nHitAxial_Mc++;}
00520 const MdcGeoWire *geowir =m_mdcGeomSvc->Wire(layerId_Mc,wireId_Mc);
00521 HepPoint3D eastP = geowir->Backward()/10.0;
00522 HepPoint3D westP = geowir->Forward() /10.0;
00523
00524 double x = (*iterMcHit)->getPositionX()/10.;
00525 double y = (*iterMcHit)->getPositionY()/10.;
00526 double z = (*iterMcHit)->getPositionZ()/10.;
00527 double u=CFtrans(x,y);
00528 double v=CFtrans(y,x);
00529 double p=sqrt(u*u+v*v);
00530
00531 vec_x_east.push_back(eastP.x());
00532 vec_y_east.push_back(eastP.y());
00533 vec_z_east.push_back(eastP.z());
00534 vec_x_west.push_back(westP.x());
00535 vec_y_west.push_back(westP.y());
00536 vec_z_west.push_back(westP.z());
00537 vec_x.push_back(x);
00538 vec_y.push_back(y);
00539 vec_z.push_back(z);
00540 vec_u.push_back(u);
00541 vec_v.push_back(v);
00542 vec_p.push_back(p);
00543 vec_slant.push_back( SlantId(layerId_Mc) );
00544
00545 m_x_east[digiId_Mc]=eastP.x();
00546 m_y_east[digiId_Mc]=eastP.y();
00547 m_z_east[digiId_Mc]=eastP.z();
00548 m_x_west[digiId_Mc]=westP.x();
00549 m_y_west[digiId_Mc]=westP.y();
00550 m_z_west[digiId_Mc]=westP.z();
00551 m_layer[digiId_Mc]=layerId_Mc;
00552 m_cell[digiId_Mc]=wireId_Mc;
00553 m_x[digiId_Mc]=x;
00554 m_y[digiId_Mc]=y;
00555 m_z[digiId_Mc]=z;
00556 m_u[digiId_Mc]=u;
00557 m_v[digiId_Mc]=v;
00558 m_p[digiId_Mc]=p;
00559 m_slant[digiId_Mc]=SlantId(layerId_Mc);
00560 }
00561 }
00562 }
00563 m_nHit_Mc=digiId_Mc;
00564 m_nHitAxial_Mc=nHitAxial_Mc;
00565 if(0==m_data) {m_nHit=digiId_Mc; m_nHitAxial=nHitAxial_Mc;}
00566 }
00567
00568
00569 vector<int> vec_hit;
00570 vector<int> vec_track_index;
00571 set<int> hit_noise;
00572 uint32_t getDigiFlag = 0;
00573 if(m_dropHot || m_combineTracking)getDigiFlag |= MdcRawDataProvider::b_dropHot;
00574 if(m_keepBadTdc) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
00575 if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
00576
00577 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
00578 if(m_debug>0) cout<<"MdcDigiVec: "<<mdcDigiVec.size()<<endl;
00579 MdcDigiVec::iterator iter1 = mdcDigiVec.begin();
00580 int digiId = 0;
00581 int nHitAxial = 0;
00582 for (;iter1 != mdcDigiVec.end(); iter1++, digiId++) {
00583 Identifier mdcId= (*iter1)->identify();
00584 int layerId = MdcID::layer(mdcId);
00585 int wireId = MdcID::wire(mdcId);
00586 double hittemp=layerId*1000+wireId;
00587 const MdcGeoWire *geowir =m_mdcGeomSvc->Wire(layerId,wireId);
00588 HepPoint3D eastP = geowir->Backward()/10.0;
00589 HepPoint3D westP = geowir->Forward() /10.0;
00590 if ((layerId>=8&&layerId<=19)||(layerId>=36)) nHitAxial++;
00591
00592 int track_index=(*iter1)->getTrackIndex();
00593 if( track_index>=1000 ) track_index-=1000;
00594 if( track_index<0 ) hit_noise.insert(hittemp);
00595 int tchannal=(*iter1)->getTimeChannel();
00596 int qchannal=(*iter1)->getChargeChannel();
00597 if( tchannal!=qchannal && m_debug>0 ) cout<<"("<<layerId<<","<<wireId<<")"<< 0<<endl;
00598 if( m_debug>3 ) cout<<"track_index in digi: "<<"("<<layerId<<","<<wireId<<" "<< track_index<<endl;
00599 if( m_debug>3 ) cout<<"event: "<<t_eventNum<<" "<<"layer: "<<layerId<<" "<<"wireId: "<<wireId<<" "<<"zeast: "<<eastP.z()<<" "<<"zwest: "<<westP.z()<<" "<<endl;
00600
00601
00602 double x =(eastP.x()+westP.x())/2.;
00603 double y =(eastP.y()+westP.y())/2.;
00604 double u=CFtrans(x,y);
00605 double v=CFtrans(y,x);
00606 if(m_debug>3) cout<<"digiId: "<<digiId<<" layer: "<<layerId<<" "<<"wireId: "<<wireId<<" "<<"x: "<<x<<" "<<"y: "<<y<<" "<< "u: "<<u<<"v: "<<v<<endl;
00607
00608 vec_track_index.push_back(track_index);
00609 vec_layer.push_back(layerId);
00610 vec_wire.push_back(wireId);
00611 vec_hit.push_back(hittemp);
00612 vec_x_east.push_back(eastP.x());
00613 vec_y_east.push_back(eastP.y());
00614 vec_z_east.push_back(eastP.z());
00615 vec_x_west.push_back(westP.x());
00616 vec_y_west.push_back(westP.y());
00617 vec_z_west.push_back(westP.z());
00618 vec_x.push_back(x);
00619 vec_y.push_back(y);
00620 vec_p.push_back(sqrt(u*u+v*v));
00621 vec_u.push_back(u);
00622 vec_v.push_back(v);
00623 vec_slant.push_back( SlantId(layerId) );
00624
00625 if(m_drawTuple){
00626 m_x_east[digiId]=eastP.x();
00627 m_y_east[digiId]=eastP.y();
00628 m_z_east[digiId]=eastP.z();
00629 m_x_west[digiId]=westP.x();
00630 m_y_west[digiId]=westP.y();
00631 m_z_west[digiId]=westP.z();
00632 m_layer[digiId]=layerId;
00633 m_cell[digiId]=wireId;
00634 m_x[digiId]=(vec_x_east[digiId]+vec_x_west[digiId])/2;
00635 m_y[digiId]=(vec_y_east[digiId]+vec_y_west[digiId])/2;
00636 m_u[digiId]=u;
00637 m_v[digiId]=v;
00638 m_p[digiId]=sqrt(u*u+v*v);
00639 m_slant[digiId]=SlantId( layerId );
00640 m_layerNhit[layerId]++;
00641 }
00642 }
00643 if(m_drawTuple){
00644 m_nHit=digiId;
00645 m_nHitAxial=nHitAxial;
00646 m_nHitStereo=m_nHit-m_nHitAxial;
00647 }
00648 if(m_debug>0) cout<<"Hit digi number: "<<digiId<<endl;
00649 if(digiId<5||nHitAxial<3) {
00650 if(m_drawTuple) m_trackFailure=1;
00651 if(m_drawTuple) ntuplehit->write();
00652 if(m_debug>0) cout<<"not enough hit "<<endl;
00653 return StatusCode::SUCCESS;
00654 }
00655
00656
00657 vector<int> vec_digitruth(digiId,0);
00658 vector<double> vec_flt_truth(digiId,999.);
00659 vector<double> vec_ztruth(digiId,999.);
00660 vector<double> vec_posflag_truth(digiId,999.);
00661 if(m_useMcInfo){
00662 int if_exit_delta_e=0;
00663 int nhit_digi=0;
00664 for(int iHit=0;iHit<digiId;iHit++){
00665 vector<int>::iterator iter_ihit = find( vec_hit_Mc.begin(),vec_hit_Mc.end(),vec_hit[iHit] );
00666 if( iter_ihit==vec_hit_Mc.end() ) {
00667 vec_digitruth[iHit]=0;
00668 if(m_drawTuple) m_digi_truth[iHit]=0;
00669 if_exit_delta_e=1;
00670 }
00671 else {
00672 vec_digitruth[iHit]=1;
00673 m_digi_truth[iHit]=1;
00674 nhit_digi++;
00675 }
00676 for(int iHit_Mc=0;iHit_Mc<vec_flt_truth_mc.size();iHit_Mc++){
00677 if(vec_hit[iHit]==vec_hit_Mc[iHit_Mc]) {
00678 vec_flt_truth[iHit]=vec_flt_truth_mc[iHit_Mc];
00679 vec_ztruth[iHit]=vec_ztruth_Mc[iHit_Mc];
00680 vec_posflag_truth[iHit]=vec_pos_flag[iHit_Mc];
00681 if(m_drawTuple){
00682 m_ztruth[iHit]=vec_ztruth_Mc[iHit_Mc];
00683 m_ltruth[iHit]=vec_flt_truth_mc[iHit_Mc];
00684 m_postruth[iHit]=vec_pos_flag[iHit_Mc];
00685 }
00686 }
00687 }
00688 if(m_debug>3) cout<<" hitPos: "<<"("<<vec_layer[iHit]<<","<<vec_wire[iHit]<<")"<<" flt_truth: "<<vec_flt_truth[iHit]<<" z_truth: "<<vec_ztruth[iHit]<<" pos_flag: "<<vec_posflag_truth[iHit]<<" digi_truth: "<<vec_digitruth[iHit]<<endl;
00689 }
00690 if( m_drawTuple ){
00691 m_e_delta=if_exit_delta_e;
00692 m_nHitDigi=nhit_digi;
00693 }
00694 }
00695
00696
00697 vector<int>::iterator laymax=max_element( vec_layer.begin(),vec_layer.end() );
00698 if(m_drawTuple) m_layer_max=*laymax;
00699
00700
00701 int nhit;
00702 int nhitaxial;
00703 if(m_data==0 && m_useMcInfo && m_fillTruth) {
00704 nhit=m_nHit_Mc;
00705 nhitaxial=m_nHitAxial_Mc;
00706 }
00707 else{
00708 nhit=digiId;
00709 nhitaxial=nHitAxial;
00710 }
00711
00712 binX=m_binx;
00713 binY=m_biny;
00714 TH2D *h1 = new TH2D("h1","h1",binX,0,M_PI,binY,-m_maxrho,m_maxrho);
00715 TH2D *h2 = new TH2D("h2","h2",binX,0,M_PI,binY,-m_maxrho,m_maxrho);
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
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786 vector< vector< vector<int> > > vec_selectNum(binX,vector< vector<int> >(binY,vector<int>() ) );
00787 M2 peak_hit_num;
00788 M1 peak_hit_list;
00789 int peak_track;
00790 if(m_method==1){
00791 FillHist(h1,vec_u,vec_v,nhit,vec_selectNum);
00792 int max_count=0;
00793 int max_binx=0;
00794 int max_biny=0;
00795 for(int i=0;i<binX;i++){
00796 for(int j=0;j<binY;j++){
00797 int count=h1->GetBinContent(i+1,j+1);
00798 if(max_count<count) {
00799 max_count=count;
00800 max_binx=i+1;
00801 max_biny=j+1;
00802 }
00803 }
00804 }
00805
00806 if(vec_selectNum[max_binx-1][max_biny-1].size() != max_count ) cout<< "ERROR IN VEC!!! "<<endl;
00807 if(m_debug>4) {
00808 cout<<"maxBin: ["<<max_binx-1<<","<<max_biny-1<<"]: "<<vec_selectNum[max_binx-1][max_biny-1].size()<<endl;
00809 for(int i =0;i<vec_selectNum[max_binx-1][max_biny-1].size();i++) cout<<" maxBin select hit: "<<vec_selectNum[max_binx-1][max_biny-1][i]<<endl;
00810 }
00811
00812 if(m_debug>4) {
00813 for(int ibinx=0;ibinx<binX;ibinx++){
00814 for(int ibiny=0;ibiny<binY;ibiny++){
00815 int bincount=h1->GetBinContent(ibinx+1,ibiny+1);
00816 if(vec_selectNum[ibinx][ibiny].size() != bincount ) cout<< "ERROR IN VECTORSELECT filling "<<endl;
00817 if(vec_selectNum[ibinx][ibiny].size() == 0 ) continue;
00818 cout<<"bin:"<<"["<<ibinx<<","<<ibiny<<"]"<<" select:"<<vec_selectNum[ibinx][ibiny].size()<<endl;
00819 for(int i=0;i<vec_selectNum[ibinx][ibiny].size();i++){
00820 int ihit=vec_selectNum[ibinx][ibiny][i];
00821 cout<<" hit: "<<ihit<<" ("<<vec_layer[ihit]<<","<<vec_wire[ihit]<<")"<<endl;
00822 }
00823 }
00824 }
00825 }
00826
00827
00828
00829 int count_use=0;
00830 double sum=0;
00831 double sum2=0;
00832 for(int i=0;i<binX;i++){
00833 for(int j=0;j<binY;j++){
00834 int count=h1->GetBinContent(i+1,j+1);
00835 if(count!=0){
00836 count_use++;
00837 sum+=count;
00838 sum2+=count*count;
00839 }
00840 }
00841 }
00842 double f_n=m_ndev1;
00843 double aver=sum/count_use;
00844 double dev=sqrt(sum2/count_use-aver*aver);
00845 int min_counts=(int)(aver+f_n*dev);
00846 if (min_counts<m_ndev2) min_counts=m_ndev2;
00847 if(m_debug>2) cout<<"aver: "<<aver<<" "<<"dev: "<<dev<<" min: "<<min_counts<<endl;
00848
00849 vector< vector <int> > hough_trans_CS(binX,vector <int> (binY,0) );
00850 vector< vector <int> > hough_trans_CS_peak(binX,vector <int> (binY,0) );
00851 for(int i=0;i<binX;i++){
00852 for(int j=0;j<binY;j++){
00853 hough_trans_CS[i][j]=h1->GetBinContent(i+1,j+1);
00854 }
00855 }
00856
00857 int npeak_1=0;
00858 for (int r=0; r<binY; r++) {
00859 double binHigh=2*m_maxrho/binY;
00860 double rho_peak=r*binHigh-m_maxrho;
00861 for (int a=0; a<binX; a++) {
00862 long max_around=0;
00863 for (int ar=a-1; ar<=a+1; ar++) {
00864 for (int rx=r-1; rx<=r+1; rx++) {
00865 int ax;
00866 if (ar<0) { ax=ar+binX; }
00867 else if (ar>=binX) { ax=ar-binX; }
00868 else { ax=ar; }
00869 if ( (ax!=a || rx!=r) && rx>=0 && rx<binY) {
00870 if (hough_trans_CS[ax][rx]>max_around) { max_around=hough_trans_CS[ax][rx]; }
00871 }
00872 }
00873 }
00874
00875 if (hough_trans_CS[a][r]<=min_counts||fabs(rho_peak)<m_rhocut) { hough_trans_CS_peak[a][r]=0; }
00876 else if (hough_trans_CS[a][r]<max_around) { hough_trans_CS_peak[a][r]=0; }
00877 else if (hough_trans_CS[a][r]==max_around) { hough_trans_CS_peak[a][r]=1; }
00878 else if (hough_trans_CS[a][r]>max_around) { hough_trans_CS_peak[a][r]=2; }
00879 if(hough_trans_CS_peak[a][r]!=0) {
00880 if(m_debug>2) cout<<" possible peak1:"<<"["<<a<<","<<r<<"]"<<": "<<hough_trans_CS_peak[a][r]<<" "<<hough_trans_CS[a][r]<<endl;
00881 npeak_1++;
00882 }
00883 }
00884 }
00885
00886
00887
00888 int npeak_2=0;
00889 for (int r=0; r<binY; r++) {
00890 for (int a=0; a<binX; a++) {
00891 if (hough_trans_CS_peak[a][r]==1) {
00892 hough_trans_CS_peak[a][r]=2;
00893 for (int ar=a-1; ar<=a+1; ar++) {
00894 for (int rx=r-1; rx<=r+1; rx++) {
00895 int ax;
00896 if (ar<0) { ax=ar+binX; }
00897 else if (ar>=binX) { ax=ar-binX; }
00898 else { ax=ar; }
00899 if ( (ax!=a || rx!=r) && rx>=0 && rx<binY) {
00900 if (hough_trans_CS_peak[ax][rx]==1) { hough_trans_CS_peak[ax][rx]=0; }
00901 }
00902 }
00903 }
00904 }
00905 if(hough_trans_CS_peak[a][r]==2) {
00906 double binHigh=2*m_maxrho/binY;
00907 double rho_peak=r*binHigh-m_maxrho;
00908 npeak_2++;
00909 if(m_debug>2){
00910 cout<<" possible peak2: "<<"["<<a<<","<<r<<"]"<<": "<<hough_trans_CS_peak[a][r]<<" "<<hough_trans_CS[a][r]<<" rhopeak: "<<rho_peak<<endl;
00911 for(int i=0;i<hough_trans_CS[a][r];i++){
00912 int hitNum=vec_selectNum[a][r][i];
00913 if(m_debug>2) cout<<" select hit: "<<hitNum<<"("<<vec_layer[hitNum]<<","<<vec_wire[hitNum]<<")"<<endl;
00914 }
00915 }
00916 }
00917 }
00918 }
00919
00920
00921 for(int i=0;i<binX;i++){
00922 for(int j=0;j<binY;j++){
00923 if(hough_trans_CS_peak[i][j]==2){
00924 h2->SetBinContent(i+1,j+1,hough_trans_CS[i][j]);
00925 }
00926 }
00927 }
00928
00929 vector<int> peakList[3];
00930 for(int n=0;n<npeak_2;n++){
00931 for (int r=0; r<binY; r++) {
00932 for (int a=0; a<binX; a++) {
00933 if (hough_trans_CS_peak[a][r]==2) {
00934 peakList[0].push_back(a);
00935 peakList[1].push_back(r);
00936 peakList[2].push_back(hough_trans_CS[a][r]);
00937 }
00938 }
00939 }
00940 }
00941
00942 if(m_drawTuple){
00943 m_npeak_1=npeak_1;
00944 m_npeak_2=npeak_2;
00945 }
00946 if(m_debug>0) {
00947 cout<<"npeak_1: "<<npeak_1<<endl;
00948 cout<<"npeak_2: "<<npeak_2<<endl;
00949 }
00950
00951
00952 int n_max;
00953 for (int na=0; na<npeak_2-1; na++) {
00954 n_max=na;
00955 for (int nb=na+1; nb<npeak_2; nb++) {
00956 if (peakList[2][n_max]<peakList[2][nb]) { n_max=nb; }
00957 }
00958 if (n_max!=na) {
00959 float swap[3];
00960 for (int i=0; i<3; i++) {
00961 swap[i]=peakList[i][na];
00962 peakList[i][na]=peakList[i][n_max];
00963 peakList[i][n_max]=swap[i];
00964 }
00965 }
00966 }
00967
00968 if(npeak_2==0||npeak_2>100){
00969 if(m_debug>0) cout<<" FAILURE IN NPEAK2 !!!!!! "<<endl;
00970 delete h1;
00971 delete h2;
00972 if(m_drawTuple){
00973 m_trackFailure=2;
00974 m_npeak_2=-999;
00975 }
00976 if(m_drawTuple) ntuplehit->write();
00977 return StatusCode::SUCCESS;
00978 }
00979
00980 if(m_debug>2){
00981 for(int n=0;n<npeak_2;n++){
00982 int bintheta=peakList[0][n];
00983 int binrho=peakList[1][n];
00984 int count=peakList[2][n];
00985 cout<<"possible peakList after SORT: "<<n<<": "<<"["<<bintheta<<","<<binrho<<"]: "<<count<<endl;
00986 for(int i=0;i<count;i++){
00987 int hitNum=vec_selectNum[bintheta][binrho][i];
00988 if(m_debug>2) cout<<" "<<" select hit:"<<":"<<hitNum<<" ------ "<<"("<<vec_layer[hitNum]<<","<<vec_wire[hitNum]<<")"<<endl;
00989 }
00990 }
00991 }
00992
00993
00994 peak_track=npeak_2;
00995 Combine(h1,m_peakCellNumX,m_peakCellNumY,vec_selectNum,peak_hit_list,peak_hit_num,peak_track,peakList);
00996 if(m_drawTuple) m_npeak_after_tracking=peak_track;
00997
00998 if(m_debug>0) cout<<"event"<<t_eventNum<<": "<<"has found: "<<peak_track<<" track."<<endl;
00999 for( int i=0;i<peak_track;i++){
01000 if(m_debug>0) cout<<"peak["<<i<<"] collect "<<peak_hit_num[i]<<" hits "<<endl;
01001 int nhitaxialselect=0;
01002 for(int j =0;j<peak_hit_num[i];j++){
01003 int hit_number=peak_hit_list[i][j];
01004 if(vec_slant[hit_number]==0) nhitaxialselect++ ;
01005 if(m_debug>0) cout<<" "<<" collect hits: ("<<vec_layer[hit_number]<<","<<vec_wire[hit_number]<<")"<<endl;
01006 }
01007 if(m_drawTuple){
01008 m_nHitSelect[i]=peak_hit_num[i];
01009 m_nHitAxialSelect[i]=nhitaxialselect;
01010 m_nHitStereoSelect[i]=peak_hit_num[i]-nhitaxialselect;
01011 m_axiallose[i]=m_nHitAxial-m_nHitAxialSelect[i];
01012 m_stereolose[i]=m_nHit-m_nHitSelect[i]-m_axiallose[i];
01013 }
01014 }
01015 }
01016
01017 vector<int> vec_hitbeforefit;
01018 vector<bool> vec_truthHit(nhit,false);
01019 int peak_fit=peak_track;
01020 if(m_setpeak_fit!=-1) {
01021 peak_fit=m_setpeak_fit;
01022 }
01023 if(m_drawTuple) m_trackNum=peak_fit;
01024
01025 vector<int> vec_hit_use;
01026
01027
01028
01029 vector<TrkRecoTrk*> vec_newTrack;
01030 vector<TrkRecoTrk*> vec_track_in_tds;
01031 vector<int> vec_hitOnTrack;
01032 vector<MdcHit*> vec_for_clean;
01033 vector<TrkRecoTrk*> vec_trk_for_clean;
01034 vector<Mytrack> vec_myTrack;
01035 int track_fit_success=0;
01036
01037 for(track_fit=0;track_fit<peak_fit;track_fit++){
01038 double d0_before_least=-999.;
01039 double phi0_before_least=-999.;
01040 double omega_before_least=-999.;
01041 map<int,double> hit_to_circle1;
01042 map<int,double> hit_to_circle2;
01043 map<int,double> hit_to_circle3;
01044 if(m_debug>0) cout<<"Do fitting track: "<<track_fit<<endl;
01045
01046 for(int i=0;i<nhit;i++) vec_truthHit[i]=false;
01047
01048
01049 if(m_debug>1) cout<<"candi track["<<track_fit<<"] collect "<<peak_hit_num[track_fit]<<" hits "<<endl;
01050
01051 for(int i=0;i<peak_hit_num[track_fit];i++){
01052 int hitNum=peak_hit_list[track_fit][i];
01053 int hittemp=vec_layer[hitNum]*1000+vec_wire[hitNum];
01054 vector<int>::iterator iter_ihit = find( vec_hit_use.begin(),vec_hit_use.end(),hittemp );
01055 if( iter_ihit==vec_hit_use.end() ) vec_truthHit[hitNum]=true;
01056 if( m_debug >1 && vec_truthHit[hitNum]==true) cout<<" "<<"collect hit :"<<"("<<vec_layer[hitNum]<<","<<vec_wire[hitNum]<<")"<<" "<<endl;;
01057 }
01058 if( m_debug >1){
01059 for(int i=0;i<nhit;i++){
01060 if(vec_truthHit[i]==false){
01061 cout<<"candi track "<<"uncollect hit: "<<endl;
01062 cout<<" "<<"uncollect hit :"<<"("<<vec_layer[i]<<","<<vec_wire[i]<<")"<<" " <<endl;
01063 }
01064 }
01065 }
01066
01067 int t_nHitAxialSelect=0;
01068 for(int i=0;i<nhit;i++){
01069 if(vec_truthHit[i]==true){
01070 if(vec_slant[i]==0) t_nHitAxialSelect++;
01071 }
01072 }
01073 if(m_debug>1) cout<<"track "<<track_fit<<" with axial select: "<<t_nHitAxialSelect<<endl;
01074 if(t_nHitAxialSelect<3 && m_drawTuple){
01075 m_x_circle[track_fit]=-999.;
01076 m_y_circle[track_fit]=-999.;
01077 m_r[track_fit]=-999.;
01078 m_chis_pt[track_fit] = -999.;
01079 m_pt[track_fit]=-999.;
01080 continue;
01081 }
01082 double circle_x=0;
01083 double circle_y=0;
01084 double circle_r=0;
01085 int leastSquare=LeastSquare(nhit,vec_x,vec_y,vec_slant,vec_layer,vec_truthHit,d0_before_least,phi0_before_least,omega_before_least,circle_x,circle_y,circle_r);
01086
01087 for(int i=0;i<nhit;i++){
01088 int hittemp=vec_layer[i]*1000+vec_wire[i];
01089 double dist_temp=sqrt(pow( (vec_y[i]-circle_y),2)+pow( (vec_x[i]-circle_x),2))-circle_r;
01090 if(dist_temp>10.) vec_truthHit[i]=false;
01091 if(m_debug>1&&vec_truthHit[i]==true) cout<<"IN LEAST1: "<<"("<<vec_layer[i]<<","<<vec_wire[i]<<"):"<<dist_temp<<endl;
01092 }
01093
01094 t_nHitAxialSelect=0;
01095 for(int i=0;i<nhit;i++){
01096 if(vec_truthHit[i]==true){
01097 if(vec_slant[i]==0) t_nHitAxialSelect++;
01098 }
01099 }
01100 if(m_debug>1) cout<<"track "<<track_fit<<"with axial2 select: "<<t_nHitAxialSelect<<endl;
01101 if(t_nHitAxialSelect<3) continue;
01102
01103 int leastSquare2=LeastSquare(nhit,vec_x,vec_y,vec_slant,vec_layer,vec_truthHit,d0_before_least,phi0_before_least,omega_before_least,circle_x,circle_y,circle_r);
01104 for(int i=0;i<nhit;i++){
01105 int hittemp=vec_layer[i]*1000+vec_wire[i];
01106 double dist_temp=sqrt(pow( (vec_y[i]-circle_y),2)+pow( (vec_x[i]-circle_x),2))-circle_r;
01107 hit_to_circle1[hittemp]=dist_temp;
01108 if(m_debug>1&&vec_truthHit[i]==true) cout<<" IN LEAST2: "<<"("<<vec_layer[i]<<","<<vec_wire[i]<<"):"<<dist_temp<<endl;
01109 }
01110
01111 if(m_drawTuple){
01112 m_x_circle[track_fit]=circle_x;
01113 m_y_circle[track_fit]=circle_y;
01114 m_r[track_fit]=circle_r;
01115 m_chis_pt[track_fit]=circle_r/333.567;
01116 }
01117
01118
01119 TrkHitList* qhits[2];
01120 const TrkFit* newFit[2];
01121 TrkHitList* qhits2[2];
01122 const TrkFit* newFit2[2];
01123 TrkRecoTrk* newTrack2[2];
01124
01125 int Q_iq[]={-1,1};
01126 double Q_Chis_2d[]={9999,9999};
01127 double Q_Chis_3d[]={9999,9999};
01128 double Q_z[]={999,999};
01129 double Q_pt2[]={999,999};
01130 double Q_tanl[]={999,999};
01131 int Q_ok[2][2]={0};
01132
01133 int q_start=0;
01134 int q_end=1;
01135 if(m_q==-1) {q_start=0;q_end=0;}
01136 if(m_q==1) {q_start=1;q_end=1;}
01137 for(int i_q=q_start;i_q<=q_end;i_q++){
01138 double d0=d0_before_least;
01139 double phi0=phi0_before_least;
01140 double omega=omega_before_least;
01141 double z0=0;
01142 double tanl=0;
01143 if(m_debug>0 ){
01144 cout<<"BY LEASTSQUARE: "<<endl;
01145 cout<<" d0: "<<d0<<" d0_mc "<<d0_mc<<endl;
01146 cout<<" phi0: "<<phi0<<" phi0_mc "<<phi0_mc<<endl;
01147 cout<<" omega0: "<<omega<<" omega_mc "<<omega_mc<<endl;
01148 }
01149
01150 vector<double> vec_flt_least;
01151 for(int i=0;i<nhit;i++){
01152 double theta_temp;
01153 if(circle_x==0||vec_x[i]-circle_x==0){
01154 theta_temp=0;
01155 }
01156 else{
01157 theta_temp=M_PI-atan2(vec_y[i]-circle_y,vec_x[i]-circle_x)+atan2(circle_y,circle_x);
01158 if(theta_temp>2*M_PI){
01159 theta_temp=theta_temp-2*M_PI;
01160 }
01161 if(theta_temp<0){
01162 theta_temp=theta_temp+2*M_PI;
01163 }
01164 }
01165 if(Q_iq[i_q]==-1) {
01166 double l_temp=circle_r*theta_temp;
01167 vec_flt_least.push_back(l_temp);
01168 }
01169 if(Q_iq[i_q]==1) {
01170 theta_temp=2*M_PI-theta_temp;
01171 double l_temp=circle_r*(theta_temp);
01172 vec_flt_least.push_back(l_temp);
01173 }
01174 }
01175 if(m_debug>2) {
01176 for(int i=0;i<nhit;i++){
01177 cout<<"By least 2d: "<< "("<<vec_layer[i]<<","<<vec_wire[i]<<"):"<<vec_flt_least[i]<<endl;
01178 }
01179 }
01180
01181
01182 if(Q_iq[i_q]==-1){
01183 d0=-d0;
01184 omega=-1./circle_r;
01185 }
01186 if(Q_iq[i_q]==1){
01187 d0=d0;
01188 omega=1./circle_r;
01189 phi0=phi0-M_PI;
01190 }
01191
01192
01193 double x_cirtemp;
01194 double y_cirtemp;
01195 double r_temp;
01196 TrkExchangePar tt(d0,phi0,omega,z0,tanl);
01197 TrkCircleMaker circlefactory;
01198 float chisum =0.;
01199 TrkRecoTrk* newTrack = circlefactory.makeTrack(tt, chisum, *m_context, m_bunchT0);
01200
01201 bool permitFlips = true;
01202 bool lPickHits = m_pickHits;
01203 circlefactory.setFlipAndDrop(*newTrack, permitFlips, lPickHits);
01204 int nDigi = digiToHots(mdcDigiVec,newTrack,vec_truthHit,getDigiFlag,vec_flt_least,hit_to_circle1,vec_for_clean);
01205 if(m_debug>2){
01206 newTrack->printAll(std::cout);
01207 }
01208
01209
01210 qhits[i_q] = newTrack->hits();
01211 TrkErrCode err=qhits[i_q]->fit();
01212 newFit[i_q] = newTrack->fitResult();
01213
01214
01215 if (!newFit[i_q] || (newFit[i_q]->nActive()<3) || err.failure()!=0) {
01216 if(m_debug>0){
01217 cout << "evt "<<t_eventNum<<" global 2d fit failed ";
01218 if(newFit[i_q]) cout <<" nAct "<<newFit[i_q]->nActive();
01219 cout<<"ERR1 failure ="<<err.failure()<<endl;
01220 }
01221
01222 Q_ok[i_q][0]=0;
01223 Q_ok[i_q][1]=0;
01224 Q_Chis_2d[i_q]=-999.;
01225 delete newTrack;
01226 continue;
01227 }else{
01228 Q_ok[i_q][0]=1;
01229 Q_ok[i_q][1]=0;
01230 if(m_debug>0) cout <<"evt "<<t_eventNum<< " global 2d fit success"<<endl;
01231 if(m_debug>2) {
01232 newTrack->print(std::cout);
01233 }
01234 Q_Chis_2d[i_q]=newFit[i_q]->chisq();
01235
01236 TrkExchangePar par = newFit[i_q]->helix(0.);
01237 d0=par.d0();
01238 phi0=par.phi0();
01239 omega=par.omega();
01240 r_temp=Q_iq[i_q]/par.omega();
01241 x_cirtemp = sin(par.phi0()) *(par.d0()+1/par.omega()) * -1.;
01242 y_cirtemp = -1.*cos(par.phi0()) *(par.d0()+1/par.omega()) * -1;
01243
01244 if(m_debug>1) cout<<" circle after globle 2d: "<<"("<<x_cirtemp<<","<<y_cirtemp<<","<<r_temp<<")"<<endl;
01245 if(m_debug>0) cout<<"pt_rec: "<<-1*Q_iq[i_q]/omega/333.567<<endl;
01246
01247 if(m_drawTuple){
01248 m_x_circle[track_fit]=x_cirtemp;
01249 m_y_circle[track_fit]=y_cirtemp;
01250 m_r[track_fit]=r_temp;
01251 m_chis_pt[track_fit] = newFit[i_q]->chisq();
01252 m_pt[track_fit]=r_temp/333.567;
01253 }
01254 }
01255
01256 int nfit2d=0;
01257 if(m_debug>1) cout<<" hot list:"<<endl;
01258 TrkHotList::hot_iterator hotIter= qhits[i_q]->hotList().begin();
01259 int lay=((MdcHit*)(hotIter->hit()))->layernumber();
01260 int wir=((MdcHit*)(hotIter->hit()))->wirenumber();
01261 int hittemp=lay*1000+wir;
01262 while (hotIter!=qhits[i_q]->hotList().end()) {
01263 if(m_debug>1){ cout <<"(" <<((MdcHit*)(hotIter->hit()))->layernumber()
01264 <<","<<((MdcHit*)(hotIter->hit()))->wirenumber()
01265 <<":"<<hotIter->isActive()<<") ";
01266 }
01267 if (hotIter->isActive()==1) nfit2d++;
01268 hotIter++;
01269 }
01270 if(m_debug>0) cout<<"charge "<<i_q<<" In 2Dfit Num Of Active Hits: "<<nfit2d<<endl;
01271
01272
01273 if(m_debug>0 && m_drawTuple){
01274 cout<<" By global 2d charge "<<i_q<<endl;
01275 cout<<" d0: " <<d0<<" "<<m_drTruth[0]<<endl;
01276 cout<<" phi0: " <<phi0<<" "<<m_phi0Truth[0]<<endl;
01277 cout<<" omega: "<<omega<<" "<<m_omegaTruth[0]<<endl;
01278 cout<<" z0: " <<z0<<" "<<m_dzTruth[0]<<endl;
01279 cout<<" tanl: "<<tanl<<" "<<m_tanl_mc[0]<<endl;
01280 }
01281
01282 if(m_debug>2) {
01283 cout<<"least:( "<<circle_x<<","<<circle_y<<","<<circle_r<<endl;
01284 cout<<"2d:( "<<x_cirtemp<<","<<y_cirtemp<<","<<r_temp<<endl;
01285 }
01286 delete newTrack;
01287
01288
01289 vector<double> vec_flt;
01290 vector<double> vec_theta_ontrack;
01291 for(int i=0;i<nhit;i++){
01292 double theta_temp;
01293 if(x_cirtemp==0||vec_x[i]-x_cirtemp==0){
01294 theta_temp=0;
01295 }
01296 else{
01297 theta_temp=M_PI-atan2(vec_y[i]-y_cirtemp,vec_x[i]-x_cirtemp)+atan2(y_cirtemp,x_cirtemp);
01298 if(theta_temp>2*M_PI){
01299 theta_temp=theta_temp-2*M_PI;
01300 }
01301 if(theta_temp<0){
01302 theta_temp=theta_temp+2*M_PI;
01303 }
01304 if(Q_iq[i_q]==-1) {
01305 double l_temp=r_temp*theta_temp;
01306 vec_flt.push_back(l_temp);
01307 vec_theta_ontrack.push_back(theta_temp);
01308
01309 }
01310 if(Q_iq[i_q]==1) {
01311 theta_temp=2*M_PI-theta_temp;
01312 double l_temp=r_temp*(theta_temp);
01313 vec_flt.push_back(l_temp);
01314 vec_theta_ontrack.push_back(theta_temp);
01315
01316 }
01317 }
01318 }
01319
01320 if(m_debug>3){
01321 for(int i=0;i<nhit;i++){
01322 cout<<"i: "<<"theta: "<<vec_theta_ontrack[i]<<" l: "<<vec_flt[i]<<endl;
01323 }
01324 }
01325
01326 for(int i=0;i<nhit;i++){
01327 int hittemp=vec_layer[i]*1000+vec_wire[i];
01328 double dist_temp=sqrt(pow( (vec_y[i]-y_cirtemp),2)+pow( (vec_x[i]-x_cirtemp),2))-r_temp;
01329 hit_to_circle2[hittemp]=dist_temp;
01330 }
01331
01332 vector<double> vec_l_Calcu_Left;
01333 vector<double> vec_l_Calcu_Right;
01334 vector<double> vec_z_Calcu_Left;
01335 vector<double> vec_z_Calcu_Right;
01336 vector<double> vec_z_Calcu;
01337 vector<double> vec_l_Calcu;
01338 vector<double> vec_delta_z;
01339 vector<double> vec_delta_l;
01340
01341 double l_Calcu_Left=-999.;
01342 double l_Calcu_Right=-999.;
01343 double z_Calcu_Left=-999.;
01344 double z_Calcu_Right=-999.;
01345 double z_Calcu=-999.;
01346 double l_Calcu=-999.;
01347 double delta_z=-999.;
01348 double delta_l=-999.;
01349
01350 MdcDigiVec::iterator iter2 = mdcDigiVec.begin();
01351 int zPos;
01352 int digiId = 0;
01353 for (;iter2 != mdcDigiVec.end(); iter2++, digiId++) {
01354 if(vec_truthHit[digiId]==false) continue;
01355 if(vec_slant[digiId]==0) continue;
01356 if(vec_theta_ontrack[digiId]>M_PI) continue;
01357 const MdcDigi* aDigi = *iter2;
01358 MdcHit *hit = new MdcHit(aDigi, m_gm);
01359 hit->setCalibSvc(m_mdcCalibFunSvc);
01360 hit->setCountPropTime(true);
01361 zPos=zPosition(*hit,m_bunchT0,x_cirtemp,y_cirtemp,r_temp,digiId,delta_z,delta_l,l_Calcu_Left,l_Calcu_Right,z_Calcu_Left,z_Calcu_Right,z_Calcu,l_Calcu,Q_iq[i_q]);
01362 if(m_debug>2) cout<<"in ZS calcu hitPos: "<<"("<<vec_layer[digiId]<<","<<vec_wire[digiId]<<")"<<" l_truth: "<<vec_flt_truth[digiId]<<" z_truth: "<<vec_ztruth[digiId]<<" l_Cal: "<<l_Calcu<<" z_Cal: "<<z_Calcu<<endl;
01363 delete hit;
01364 if(zPos==-1||zPos==-2) continue;
01365 vec_l_Calcu_Left.push_back(l_Calcu_Left);
01366 vec_l_Calcu_Right.push_back(l_Calcu_Right);
01367 vec_z_Calcu_Left.push_back(z_Calcu_Left);
01368 vec_z_Calcu_Right.push_back(z_Calcu_Right);
01369 vec_z_Calcu.push_back(z_Calcu);
01370 vec_l_Calcu.push_back(l_Calcu_Right);
01371 vec_delta_z.push_back(delta_z);
01372 vec_delta_l.push_back(delta_l);
01373 }
01374 int nHitUseInZs=vec_delta_z.size();
01375
01376 if(m_drawTuple){
01377 for(int i=0;i<nHitUseInZs;i++){
01378 m_nHitStereo_use=vec_delta_z.size();
01379 m_lCalcuLeft[i]=vec_l_Calcu_Left[i];
01380 m_lCalcuRight[i]=vec_l_Calcu_Right[i];
01381 m_zCalcuLeft[i]=vec_z_Calcu_Left[i];
01382 m_zCalcuRight[i]=vec_z_Calcu_Right[i];
01383 m_lCalcu[i]=vec_l_Calcu[i];
01384 m_zCalcu[i]=vec_z_Calcu[i];
01385 m_delta_z[i]=vec_delta_z[i];
01386 m_delta_l[i]=vec_delta_l[i];
01387 }
01388 }
01389
01390
01391 vector< vector< vector <double> > > point(2, vector< vector<double> >(2,vector<double>() ));
01392 for(int iHit=0;iHit<nHitUseInZs;iHit++){
01393 point[0][0].push_back(vec_l_Calcu_Left[iHit]);
01394 point[0][1].push_back(vec_z_Calcu_Left[iHit]);
01395 point[1][0].push_back(vec_l_Calcu_Right[iHit]);
01396 point[1][1].push_back(vec_z_Calcu_Right[iHit]);
01397 }
01398 vector<int> ambigList;
01399 if(nHitUseInZs!=0){
01400 AmbigChoose(point,nHitUseInZs,ambigList,tanl,z0);
01401 if(m_drawTuple){
01402 m_tanl_Cal=tanl;
01403 m_z0_Cal=z0;
01404 }
01405
01406 if(m_useMcInfo && m_drawTuple){
01407 int ambig_no_match_num=0;
01408 for(int iHit=0;iHit<nHitUseInZs;iHit++){
01409 if(ambigList[iHit]==0) ambigList[iHit]=1;
01410 else ambigList[iHit]=0;
01411 if(ambigList[iHit]!=m_postruth[iHit]) { m_ambig_no_match=1; ambig_no_match_num++; }
01412 }
01413 m_ambig_no_match_propotion=(double)(ambig_no_match_num)/m_nHitStereo_use;
01414 for (int iHit=0;iHit<ambigList.size();iHit++){
01415 m_amb[iHit]=ambigList[iHit];
01416 }
01417 }
01418 }
01419
01420 if(m_debug>0 && m_drawTuple){
01421 cout<<"By 3d track zs fit:"<<endl;
01422 cout<<"d0: "<<d0<<" mc : "<<m_drTruth[0]<<endl;
01423 cout<<"phi0: "<<phi0<<" mc : "<<m_phi0Truth[0]<<endl;
01424 cout<<"omega: "<<omega<<" mc : "<<m_omegaTruth[0]<<endl;
01425 cout<<"z0: "<<z0<<" mc : "<<m_dzTruth[0]<<endl;
01426 cout<<"tanl: "<<tanl<<" mc : "<<m_tanl_mc[0]<<endl;
01427 }
01428
01429 if(m_mcpar==1){
01430 d0=m_drTruth[0];
01431 phi0=m_phi0Truth[0];
01432 if(Q_iq[i_q]==-1) {phi0=m_phi0Truth[0]-3./2.*M_PI;omega=m_omegaTruth[0];}
01433 else {phi0=m_phi0Truth[0]-1./2.*M_PI;omega=-1*m_omegaTruth[0];}
01434 z0=m_dzTruth[0];
01435 tanl=m_tanl_mc[0];
01436 }
01437
01438
01439
01440 if(m_debug>0) cout<< "become 3d fit "<<endl;
01441 TrkExchangePar tt2(d0,phi0,omega,z0,tanl);
01442 TrkHelixMaker helixfactory;
01443 chisum =0.;
01444 newTrack2[i_q] = helixfactory.makeTrack(tt2, chisum, *m_context, m_bunchT0);
01445 vec_trk_for_clean.push_back(newTrack2[i_q]);
01446 permitFlips = true;
01447 lPickHits = m_pickHits;
01448 helixfactory.setFlipAndDrop(*newTrack2[i_q], permitFlips, lPickHits);
01449 int nDigi2 = digiToHots2(mdcDigiVec,newTrack2[i_q],vec_truthHit,vec_flt_truth,getDigiFlag,vec_slant,vec_l_Calcu,vec_flt,vec_theta_ontrack,vec_hitbeforefit,hit_to_circle2,vec_for_clean,tanl);
01450 if(m_debug>2){
01451 cout<<__FILE__<<__LINE__<<"AFTER digiToHots 3d ---------BEFORE 3d fit"<<endl;
01452 newTrack2[i_q]->printAll(std::cout);
01453 }
01454
01455 qhits2[i_q] = newTrack2[i_q]->hits();
01456 TrkErrCode err2=qhits2[i_q]->fit();
01457 newFit2[i_q] = newTrack2[i_q]->fitResult();
01458 int nActiveHit = newTrack2[i_q]->hots()->nActive();
01459 int fitSuccess_3d = 0;
01460 if (!newFit2[i_q] || (newFit2[i_q]->nActive()<5) || err2.failure()!=0) {
01461 fitSuccess_3d=0;
01462 Q_ok[i_q][1]=0;
01463 Q_Chis_3d[i_q]=-999.;
01464 if(m_debug>0){
01465 cout << "evt "<<t_eventNum<<" global 3d fit failed ";
01466 if(!newFit2[i_q]) cout<<" newfit2 point is NULL!!!" <<endl;
01467 if(newFit2[i_q]) cout <<" nAct "<<newFit2[i_q]->nActive();
01468 cout<<"ERR2 failure ="<<err2.failure()<<endl;
01469 }
01470 }else{
01471
01472 TrkExchangePar par2 = newFit2[i_q]->helix(0.);
01473 if( abs( 1/(par2.omega()) )>0.001){
01474 fitSuccess_3d=1;
01475 Q_Chis_3d[i_q]=newFit2[i_q]->chisq();
01476 Q_ok[i_q][1]=1;
01477 if(m_debug>0) cout <<"evt "<<t_eventNum<< " global 3d fit success "<<err2.failure()<<endl;
01478 if(m_debug>2){
01479 cout<<__FILE__<<__LINE__<<"AFTER fit 3d "<<endl;
01480 newTrack2[i_q]->print(std::cout);
01481 }
01482 } else {
01483 fitSuccess_3d=0;
01484 Q_ok[i_q][1]=0;
01485 Q_Chis_3d[i_q]=-999.;
01486 if(m_debug>2) cout<<"3dfit failure of omega "<<endl;
01487 }
01488 bool badQuality = false;
01489
01490
01491 if(fabs(Q_Chis_3d[i_q])>m_qualityFactor*m_dropTrkChi2Cut ){
01492 if(m_debug>0){
01493 std::cout<<__FILE__<<" "<<__LINE__<<" drop track by chi2 "<<Q_Chis_3d[i_q]<<" > "<<m_qualityFactor<<" * "<<m_dropTrkChi2Cut <<std::endl;
01494 }
01495 badQuality=1;
01496 }
01497 if( fabs(par2.d0())>m_qualityFactor*m_dropTrkDrCut) {
01498 if(m_debug>0){
01499 std::cout<<__FILE__<<" "<<__LINE__<<" drop track by d0 "<<par2.d0()<<" > "<<m_qualityFactor<<" * "<<m_dropTrkDrCut <<std::endl;
01500 }
01501 badQuality=1;
01502 }
01503 if( fabs(par2.z0())>m_qualityFactor*m_dropTrkDzCut) {
01504 if(m_debug>0){
01505 std::cout<<__FILE__<<" "<<__LINE__<<" drop track by z0 "<<par2.z0()<<" > "<<m_qualityFactor<<" * "<<m_dropTrkDzCut <<std::endl;
01506 }
01507 badQuality=1;
01508 }
01509 if( (fabs(Q_Chis_3d[i_q])/qhits2[i_q]->nHit()) > m_qualityFactor*m_dropTrkChi2NdfCut) {
01510 if(m_debug>0){
01511 std::cout<<__FILE__<<" "<<__LINE__<<" drop track by chi2/ndf"<<(fabs(Q_Chis_3d[i_q])/qhits2[i_q]->nHit()) <<" > "<<m_qualityFactor<<" * "<<m_dropTrkChi2NdfCut<<std::endl;
01512 }
01513 badQuality=1;
01514 }
01515 if( nActiveHit <= m_dropTrkNhitCut) {
01516 if(m_debug>0){
01517 std::cout<<__FILE__<<" "<<__LINE__<<" drop track by nhit"<<nActiveHit <<" <= "<<m_qualityFactor<<" * "<<m_dropTrkNhitCut<<std::endl;
01518 }
01519 badQuality=1;
01520 }
01521 if(badQuality) {
01522 fitSuccess_3d=0;
01523 Q_ok[i_q][1]=0;
01524 Q_Chis_3d[i_q]=-999.;
01525 if(m_debug>2) cout<<"3dfit failure of bad quality"<<endl;
01526 }
01527
01528 }
01529
01530
01531 if(fitSuccess_3d==1){
01532 TrkExchangePar par2 = newFit2[i_q]->helix(0.);
01533 if(m_debug>0){
01534 cout<<"BY 3d global fit charge "<<i_q<<endl;
01535 cout<<"d0: "<<par2.d0()<<endl;
01536 cout<<"phi0: "<<par2.phi0()<<endl;
01537 cout<<"w: "<<par2.omega()<<endl;
01538 cout<<"z: "<<par2.z0()<<endl;
01539 cout<<"tanl: "<<par2.tanDip()<<endl;
01540 }
01541 r_temp=Q_iq[i_q]/par2.omega();
01542 x_cirtemp = sin(par2.phi0()) *(par2.d0()+1/par2.omega()) * -1.;
01543 y_cirtemp = -1.*cos(par2.phi0()) *(par2.d0()+1/par2.omega()) * -1;
01544 if(m_debug>0) cout<<" circle after globle 3d: "<<"("<<x_cirtemp<<","<<y_cirtemp<<","<<r_temp<<")"<<endl;
01545 double pxy=-1*r_temp/333.567;
01546 double pz=pxy*par2.tanDip();
01547 double pxyz=sqrt(pxy*pxy+pz*pz);
01548
01549 Q_pt2[i_q]=fabs(pxy);
01550 Q_z[i_q]=par2.z0();
01551 Q_tanl[i_q]=par2.tanDip();
01552 } else{
01553 continue;
01554 }
01555
01556
01557 for(int i=0;i<nhit;i++){
01558 int hittemp=vec_layer[i]*1000+vec_wire[i];
01559 double dist_temp=sqrt(pow( (vec_y[i]-y_cirtemp),2)+pow( (vec_x[i]-x_cirtemp),2))-r_temp;
01560 hit_to_circle3[hittemp]=dist_temp;
01561 }
01562 }
01563
01564 if(m_debug>3){
01565 cout<<"chis 2d: "<<Q_Chis_2d[0]<<" "<<Q_Chis_2d[1]<<endl;
01566 cout<<"chis 3d: "<<Q_Chis_3d[0]<<" "<<Q_Chis_3d[1]<<endl;
01567 cout<<"Z: "<<Q_z[0]<<" "<<Q_z[1]<<endl;
01568 cout<<"chis ok-: "<<Q_ok[0][0]<<" "<<Q_ok[0][1]<<endl;
01569 cout<<"chis ok+: "<<Q_ok[1][0]<<" "<<Q_ok[1][1]<<endl;
01570 }
01571
01572 int Q;
01573 int Q_judge[2]={0};
01574 if(Q_ok[0][1]==1) Q_judge[0]=1;
01575 if(Q_ok[1][1]==1) Q_judge[1]=1;
01576 if(Q_judge[0]==1&&Q_judge[1]==0) Q=-1;
01577 if(Q_judge[0]==0&&Q_judge[1]==1) Q=1;
01578 if(Q_judge[0]==1&&Q_judge[1]==1) {
01579 if(fabs(Q_z[0])>fabs(Q_z[1])) Q=1;
01580 else Q=-1;
01581 }
01582 if(Q_judge[0]==0&&Q_judge[1]==0) {Q=0; continue;}
01583 int q_choose=0;
01584 if(Q==-1) q_choose=0;
01585 if(Q== 1) q_choose=1;
01586 TrkExchangePar par_final = newFit2[q_choose]->helix(0.);
01587 if(m_debug>0){
01588 cout<<"after q choose By 3d global fit: "<<Q<<endl;
01589 cout<<"d0: "<<par_final.d0()<<endl;
01590 cout<<"phi0: "<<par_final.phi0()<<endl;
01591 cout<<"w: "<<par_final.omega()<<endl;
01592 cout<<"z: "<<par_final.z0()<<endl;
01593 cout<<"tanl: "<<par_final.tanDip()<<endl;
01594 }
01595 double pxy=-1*Q_iq[q_choose]/par_final.omega()/333.567;
01596 double pz=pxy*par_final.tanDip();
01597 double pxyz=sqrt(pxy*pxy+pz*pz);
01598
01599 if(m_debug>0) cout<<"pt2_rec: "<<pxy<<endl;
01600 if(m_drawTuple){
01601 m_pt2_rec_final[track_fit_success]=pxy;
01602 m_p_rec_final[track_fit_success]=pxyz;
01603 m_d0_final[track_fit_success]=par_final.d0();
01604 m_phi0_final[track_fit_success]=par_final.phi0();
01605 m_omega_final[track_fit_success]=par_final.omega();
01606 m_z0_final[track_fit_success]=par_final.z0();
01607 m_tanl_final[track_fit_success]=par_final.tanDip();
01608 m_Q[track_fit_success]=Q;
01609 }
01610
01611 int nfit3d=0;
01612 if(m_debug>1) cout<<" hot list:"<<endl;
01613 TrkHotList::hot_iterator hotIter2= qhits2[q_choose]->hotList().begin();
01614 while (hotIter2!=qhits2[q_choose]->hotList().end()) {
01615 int lay=((MdcHit*)(hotIter2->hit()))->layernumber();
01616 int wir=((MdcHit*)(hotIter2->hit()))->wirenumber();
01617 int hittemp=lay*1000+wir;
01618 if(m_debug>1){ cout <<"(" <<((MdcHit*)(hotIter2->hit()))->layernumber()
01619 <<","<<((MdcHit*)(hotIter2->hit()))->wirenumber()
01620 <<":"<<hotIter2->isActive()<<") ";
01621 }
01622 if( hotIter2->isActive()==1) {
01623 nfit3d++;
01624 if(track_fit_success==0) vec_hitOnTrack.push_back(hittemp);
01625 }
01626 hotIter2++;
01627 }
01628 if(m_debug>0) cout<<"In 3D fit Num Of Active Hits: "<<nfit3d<<endl;
01629
01630 track_fit_success++;
01631 int choise_temp,choise_no;
01632 if(Q==-1) choise_temp=0,choise_no=1;
01633 if(Q==1) choise_temp=1,choise_no=0;
01634 vec_newTrack.push_back( newTrack2[choise_temp] );
01635 }
01636
01637 for(int iTrack=0;iTrack<vec_newTrack.size();iTrack++){
01638 const TrkFit* newFit_combine = vec_newTrack[iTrack]->fitResult();
01639 TrkExchangePar par_combine= newFit_combine->helix(0.);
01640 double cx_combine=(-333.567/par_combine.omega()+par_combine.d0())*cos(par_combine.phi0());
01641 double cy_combine=(-333.567/par_combine.omega()+par_combine.d0())*sin(par_combine.phi0());
01642 double phi_combine=par_combine.phi0()-M_PI/2;
01643 double pxy_combine=1./par_combine.omega()/333.567;
01644 if(pxy_combine<0) phi_combine+=M_PI;
01645 if(phi_combine<0) phi_combine+=2*M_PI;
01646 if(phi_combine>2*M_PI) phi_combine-=2*M_PI;
01647
01648 Mytrack mytrack;
01649 mytrack.pt=pxy_combine;
01650 mytrack.charge=(int)(fabs(pxy_combine)/pxy_combine);
01651 mytrack.phi=phi_combine;
01652 mytrack.r=-333.567/par_combine.omega();
01653 mytrack.x_cir=cx_combine;
01654 mytrack.y_cir=cy_combine;
01655 mytrack.newTrack=vec_newTrack[iTrack];
01656 mytrack.use_in_tds=true;
01657 mytrack.vec_hit_on_trk=vec_hitOnTrack;
01658
01659 vec_myTrack.push_back(mytrack);
01660
01661 }
01662
01663 std::sort(vec_myTrack.begin(),vec_myTrack.end(),compare);
01664
01665 for(int i=0;i<vec_newTrack.size();i++){
01666 Mytrack *mytrack_i=&(vec_myTrack[i]);
01667 if(mytrack_i->use_in_tds==false) continue;
01668 for(int j=i+1;j<vec_newTrack.size();j++){
01669 Mytrack *mytrack_j=&(vec_myTrack[j]);
01670 if(mytrack_j->use_in_tds==false) continue;
01671 if( fabs(mytrack_j->phi-mytrack_i->phi)<0.1 ) mytrack_j->use_in_tds=false;
01672 if(fabs(mytrack_j->r)*(1.-0.25) <= fabs(mytrack_i->r) && fabs(mytrack_i->r) <= fabs(mytrack_j->r)*(1.+0.25)){
01673 if(fabs(mytrack_j->x_cir-mytrack_i->x_cir) <= fabs(mytrack_j->r)*(0.25)&& fabs(mytrack_j->y_cir-mytrack_i->y_cir) <= fabs(mytrack_j->r)*(0.25) ){
01674 if(mytrack_j->charge==mytrack_i->charge)
01675 mytrack_j->use_in_tds=false;
01676 }
01677 }
01678 }
01679 }
01680
01681 int nTrack_tds=0;
01682 vector<Mytrack> vec_mytrack_in_tds;
01683 for(int i=0;i<track_fit_success;i++){
01684 Mytrack mytrack=vec_myTrack[i];
01685 if(mytrack.use_in_tds==false) continue;
01686 vec_mytrack_in_tds.push_back(mytrack);
01687 nTrack_tds++;
01688 }
01689
01690 if(m_drawTuple){
01691 m_trackNum_tds=nTrack_tds;
01692 m_trackNum_fit=track_fit_success;
01693 }
01694
01695
01696
01697
01698
01699
01700
01701
01702
01703
01704
01705
01706
01707
01708 int outputTrk=0;
01709 int itrack=0;
01710 for(int i=0;i<nTrack_tds;i++){
01711 Mytrack mytrack=vec_mytrack_in_tds[i];
01712 const TrkFit* newFit_tds= mytrack.newTrack->fitResult();
01713 TrkExchangePar par_tds= newFit_tds->helix(0.);
01714
01715
01716
01717 double cr= 1./ par_tds.omega();
01718 double cx= sin(par_tds.phi0()) *(par_tds.d0()+1/par_tds.omega()) * -1.;
01719 double cy= -1.*cos(par_tds.phi0()) *(par_tds.d0()+1/par_tds.omega()) * -1;
01720
01721 unsigned bestIndex = 0;
01722 double bestDiff = 1.0e+20;
01723 double cR, cX, cY;
01724 vector<double> a0,a1,a2,a3,a4;
01725 vector<double> zdelta;
01726 RecMdcTrackCol::iterator iteritrk = trackList->begin();
01727 int itrk=0;
01728 for(;iteritrk!=trackList->end();iteritrk++){
01729 double pt=(*iteritrk)->pxy();
01730 a0.push_back( (*iteritrk)->helix(0) );
01731 a1.push_back( (*iteritrk)->helix(1) );
01732 a2.push_back( (*iteritrk)->helix(2) );
01733 a3.push_back( (*iteritrk)->helix(3) );
01734 a4.push_back( (*iteritrk)->helix(4) );
01735 zdelta.push_back( par_tds.z0()-(*iteritrk)->helix(3) );
01736
01737
01738
01739 cR=(-333.567)/a2[itrk];
01740 cX=(cR+a0[itrk])*cos(a1[itrk]);
01741 cY=(cR+a0[itrk])*sin(a1[itrk]);
01742
01743 if(fabs(cr)*(1.-0.25) <= fabs(cR) && fabs(cR) <= fabs(cr)*(1.+0.25)){
01744 if(fabs(cx-cX) <= fabs(cr)*(0.25)&& fabs(cy-cY) <= fabs(cr)*(0.25) ){
01745 double diff = fabs((fabs(cr)-fabs(cR))*(cx-cX)*(cy-cY));
01746 if(diff < bestDiff){
01747 bestDiff = diff;
01748 bestIndex = itrk;
01749 }
01750 }
01751 }
01752 itrk++;
01753 }
01754
01755 if(bestDiff == 1.0e20) { if(m_debug>0) cout<<" no merge "<<endl; }
01756 else continue;
01757
01758
01759
01760
01761
01762
01763
01764
01765
01766
01767
01768
01769
01770 MdcTrack* mdcTrack;
01771 mdcTrack= new MdcTrack(mytrack.newTrack);
01772 vec_track_in_tds.push_back(mytrack.newTrack);
01773 int tkStat = 4;
01774 int tkId = nTdsTk + itrack;
01775 mdcTrack->storeTrack(tkId, trackList, hitList, tkStat);
01776 if(m_drawTuple) m_hitOnTrk[outputTrk]=mdcTrack->track().hots()->nActive();
01777 outputTrk++;
01778 delete mdcTrack;
01779
01780 double pxy=1/par_tds.omega()/333.567;
01781 double pz=pxy*par_tds.tanDip();
01782 double pxyz=sqrt(pxy*pxy+pz*pz);
01783 if(m_drawTuple){
01784 m_pt2_rec_tds[itrack]=pxy;
01785 m_p_rec_tds[itrack]=pxyz;
01786 m_d0_tds[itrack]=par_tds.d0();
01787 m_phi0_tds[itrack]=par_tds.phi0();
01788 m_omega_tds[itrack]=par_tds.omega();
01789 m_z0_tds[itrack]=par_tds.z0();
01790 m_tanl_tds[itrack]=par_tds.tanDip();
01791 m_Q_tds[itrack]=pxy>0?1:-1;
01792 }
01793 itrack++;
01794 }
01795 if(m_drawTuple) m_outputTrk=outputTrk;
01796 int nTdsTk_temp=trackList->size();
01797 if(m_debug>0) m_mdcPrintSvc->printRecMdcTrackCol();
01798
01799
01800
01801
01802
01803
01804
01805
01806
01807
01808
01809
01810
01811
01812
01813
01814
01815
01816
01817
01818
01819
01820
01821
01822
01823 for(int i=0;i<vec_for_clean.size();i++){
01824 delete vec_for_clean[i];
01825 }
01826 for(int i=0;i<vec_trk_for_clean.size();i++){
01827
01828
01829 vector<TrkRecoTrk*>::iterator iterTrk =find(vec_track_in_tds.begin(),vec_track_in_tds.end(),vec_trk_for_clean[i]);
01830 if(iterTrk ==vec_track_in_tds.end() ) delete vec_trk_for_clean[i];
01831
01832 }
01833 delete h1;
01834 delete h2;
01835
01836 if(m_drawTuple) ntuplehit->write();
01837 if(m_debug>0) cout<<"end event" <<endl;
01838 return StatusCode::SUCCESS;
01839 }
01840
01841
01842 StatusCode HoughValidUpdate::finalize() {
01843 MsgStream log(msgSvc(), name());
01844 delete m_bfield ;
01845 delete m_context ;
01846 log << MSG::INFO << "in finalize()" << endreq;
01847 return StatusCode::SUCCESS;
01848 }
01849 int HoughValidUpdate::digiToHots(MdcDigiVec mdcDigiVec,TrkRecoTrk* newTrack,vector<bool> vec_truthHit,uint32_t getDigiFlag,vector<double> vec_flt_least,map<int,double> hit_to_circle1,vector<MdcHit*>& vec_for_clean){
01850 TrkHitList* trkHitList = newTrack->hits();
01851 int digiId=0;
01852 if(m_debug>0) cout<<"MdcDigiVec(in 2d fit): "<<mdcDigiVec.size()<<endl;
01853 MdcDigiVec::iterator iter = mdcDigiVec.begin();
01854 for (;iter != mdcDigiVec.end(); iter++, digiId++) {
01855 if(vec_truthHit[digiId]==false) continue;
01856 const MdcDigi* aDigi = *iter;
01857 MdcHit *hit = new MdcHit(aDigi, m_gm);
01858 vec_for_clean.push_back(hit);
01859
01860 Identifier id = aDigi->identify();
01861 int layer = MdcID::layer(id);
01862 int wire = MdcID::wire(id);
01863 int hittemp=layer*1000+wire;
01864 if ((layer>=0&&layer<=7)||(layer>=20&&layer<=35)) {continue;}
01865 hit->setCalibSvc(m_mdcCalibFunSvc);
01866 hit->setCountPropTime(true);
01867
01868 int newAmbig = 0;
01869
01870 MdcRecoHitOnTrack temp(*hit, newAmbig, m_bunchT0*1.e9);
01871 MdcHitOnTrack* newhot = &temp;
01872 double fltLen = m_gm->Layer(layer)->rMid();
01873
01874 newhot->setFltLen(vec_flt_least[digiId]);
01875
01876
01877
01878 double ddCut;
01879
01880
01881 if(layer<8) ddCut=1.0;
01882 else ddCut=1.0;
01883 int use_in2d=1;
01884 if(hit->driftDist(m_bunchT0,0)>ddCut||vec_truthHit[digiId]==false||fabs(hit_to_circle1[hittemp])>10.) { use_in2d=0; }
01885 if(m_debug>1) 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())<<" m_bunchT0 "<<m_bunchT0<<" hit to circle: "<<hit_to_circle1[hittemp]<<" dist: "<<hit->driftDist(m_bunchT0,0)<<" use?: "<<use_in2d<<endl;
01886
01887 if(use_in2d==0) continue;
01888 trkHitList->appendHot(newhot);
01889
01890 }
01891 if(m_debug>0) std::cout<<std::endl;
01892 return trkHitList->nHit();
01893 }
01894
01895 double HoughValidUpdate::CFtrans(double x,double y){
01896 return 2*x/(x*x+y*y);
01897 }
01898
01899 int HoughValidUpdate::SlantId(int layer){
01900 if ((layer>=0&&layer<=7)||(layer>=20&&layer<=35)) {
01901 if ((layer>=0&&layer<=3)||(layer>=20&&layer<=23)||(layer>=28&&layer<=31)){
01902
01903 return -1;
01904 }else{
01905 return 1;
01906 }
01907 } else{
01908 return 0;
01909 }
01910 }
01911
01912 int HoughValidUpdate::GetMcInfo(){
01913 MsgStream log(msgSvc(), name());
01914 StatusCode sc;
01915 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
01916 if (!mcParticleCol) {
01917 log << MSG::ERROR << "Could not find McParticle" << endreq;
01918 return -999;
01919 }
01920
01921 int itk = 0;
01922 int t_qTruth = 0;
01923 int t_pidTruth = -999;
01924 int t_McTkId = -999;
01925 Helix* mchelix;
01926 vector<int> vec_decay;
01927 McParticleCol::iterator iter_mc = mcParticleCol->begin();
01928 for (;iter_mc != mcParticleCol->end(); iter_mc++ ) {
01929
01930 t_pidTruth = (*iter_mc)->particleProperty();
01931 bool t_decay= (*iter_mc)->decayInFlight();
01932 if(m_debug>0) cout<<"Decay : "<<t_decay<<endl;
01933 vec_decay.push_back(t_decay);
01934
01935 if(m_debug>2) cout<< "tk "<<itk<< " pid="<< t_pidTruth<< endl;
01936 if((m_pid!=0) && (t_pidTruth != m_pid)){ continue; }
01937 if(itk>=10) break;
01938
01939 int pid = t_pidTruth;
01940 if( pid == 0 ) {
01941 log << MSG::WARNING << "Wrong particle id" <<endreq;
01942 }else{
01943 IPartPropSvc* p_PartPropSvc;
01944 static const bool CREATEIFNOTTHERE(true);
01945 StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
01946 if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) {
01947 std::cout<< " Could not initialize Particle Properties Service" << std::endl;
01948 }
01949 HepPDT::ParticleDataTable* p_particleTable = p_PartPropSvc->PDT();
01950 std::string name;
01951 if( p_particleTable->particle(pid) ){
01952 name = p_particleTable->particle(pid)->name();
01953 t_qTruth = p_particleTable->particle(pid)->charge();
01954 }else if( p_particleTable->particle(-pid) ){
01955 name = "anti " + p_particleTable->particle(-pid)->name();
01956 t_qTruth = (-1)*p_particleTable->particle(-pid)->charge();
01957 }
01958 if(m_debug>2) std::cout << " particle "<<name<<" charge= "<<t_qTruth<<std::endl;
01959 }
01960
01961 t_McTkId = (*iter_mc)->trackIndex();
01962 double px = (*iter_mc)->initialFourMomentum().px();
01963 double py = (*iter_mc)->initialFourMomentum().py();
01964 double pz = (*iter_mc)->initialFourMomentum().pz();
01965
01966
01967 mchelix = new Helix((*iter_mc)->initialPosition().v(), (*iter_mc)->initialFourMomentum().v(), t_qTruth);
01968 mchelix->pivot( HepPoint3D(0.,0.,0.) );
01969
01970
01971 if(m_debug>2){
01972 std::cout<<"Truth tk "<<itk<<" pid "<<pid<<" charge "<<t_qTruth
01973 << " helix "<< mchelix->a()<<" p "<<mchelix->momentum(0.)<<std::endl;
01974 }
01975 if(m_drawTuple){
01976 m_pzTruth[itk]=pz;
01977 m_pidTruth[itk] = t_pidTruth;
01978 m_costaTruth[itk] = mchelix->cosTheta();
01979 m_ptTruth[itk] = mchelix->pt();
01980 m_pTruth[itk] = mchelix->momentum(0.).mag();
01981 m_qTruth[itk] = t_qTruth;
01982 m_drTruth[itk] = mchelix->dr();
01983 m_phi0Truth[itk] = mchelix->phi0();
01984 m_omegaTruth[itk] =1./(mchelix->radius());
01985 m_dzTruth[itk] = mchelix->dz();
01986 m_tanl_mc[itk] =mchelix->tanl();
01987 }
01988 itk++;
01989 if(m_debug>2){
01990 cout<<" d0: " <<" "<<mchelix->dr()<<endl;
01991 cout<<" phi0: " <<" "<<mchelix->phi0()<<endl;
01992 cout<<" omega: "<<" "<<1./(mchelix->radius())<<endl;
01993 cout<<" z0: " <<" "<<mchelix->dz()<<endl;
01994 cout<<" tanl: "<<" "<<mchelix->tanl()<<endl;
01995 cout<<" pt: " <<" "<<mchelix->pt()<<endl;
01996 cout<<" costaTruth: " <<" "<<mchelix->cosTheta()<<endl;
01997 }
01998 delete mchelix;
01999 }
02000 if(m_drawTuple){
02001 vector<int>::iterator iter_idecay=find(vec_decay.begin(),vec_decay.end(),1 );
02002 if (iter_idecay==vec_decay.end() ) m_decay=0;
02003 else m_decay=1;
02004 m_nTrkMC = itk;
02005 }
02006 if(itk!=1) {
02007 if(m_debug>2) std::cout<<"WARNING run "<<t_runNum<<" evt "<<t_eventNum<<" not single event. nTrkMc="<<itk<<std::endl;
02008 return -1;
02009 }else{
02010 if(m_debug>2) std::cout<<"nTrkMc=1"<<std::endl;
02011 return 1;
02012 }
02013
02014 }
02015
02016
02017 int HoughValidUpdate::LeastSquare(int n,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,double& circleX,double& circleY,double& circleR){
02018
02019 double x_sum=0;
02020 double y_sum=0;
02021 double x2_sum=0;
02022 double y2_sum=0;
02023 double x3_sum=0;
02024 double y3_sum=0;
02025 double x2y_sum=0;
02026 double xy2_sum=0;
02027 double xy_sum=0;
02028 double a1=0;
02029 double a2=0;
02030 double b1=0;
02031 double b2=0;
02032 double c1=0;
02033 double c2=0;
02034 double x_aver,y_aver,r2;
02035 int use_in_least=0;
02036 for(int i=0;i<n;i++){
02037 if(vec_truthHit[i]==false) continue;
02038 if(vec_slant[i]!=0) continue;
02039 x_sum=x_sum+vec_x[i];
02040 y_sum=y_sum+vec_y[i];
02041 x2_sum=x2_sum+vec_x[i]*vec_x[i];
02042 y2_sum=y2_sum+vec_y[i]*vec_y[i];
02043 x3_sum=x3_sum+vec_x[i]*vec_x[i]*vec_x[i];
02044 y3_sum=y3_sum+vec_y[i]*vec_y[i]*vec_y[i];
02045 x2y_sum=x2y_sum+vec_x[i]*vec_x[i]*vec_y[i];
02046 xy2_sum=xy2_sum+vec_x[i]*vec_y[i]*vec_y[i];
02047 xy_sum=xy_sum+vec_x[i]*vec_y[i];
02048 use_in_least++;
02049 }
02050 a1=x2_sum-x_sum*x_sum/n;
02051 a2=xy_sum-x_sum*y_sum/n;
02052 b1=a2;
02053 b2=y2_sum-y_sum*y_sum/n;
02054 c1=(x3_sum+xy2_sum-x_sum*(x2_sum+y2_sum)/n)/2.;
02055 c2=(y3_sum+x2y_sum-y_sum*(x2_sum+y2_sum)/n)/2.;
02056
02057 x_aver=x_sum/n;
02058 y_aver=y_sum/n;
02059
02060 double x_cirtemp =(b1*c2-b2*c1)/(b1*b1-a1*b2);
02061 double y_cirtemp =(b1*c1-a1*c2)/(b1*b1-a1*b2);
02062 r2=(x2_sum+y2_sum-2*x_cirtemp*x_sum-2*y_cirtemp*y_sum)/n+x_cirtemp*x_cirtemp+y_cirtemp*y_cirtemp;
02063 double r_temp=sqrt(r2);
02064
02065 circleX = x_cirtemp;
02066 circleY = y_cirtemp;
02067 circleR = r_temp;
02068
02069 double d0_temp=sqrt(x_cirtemp*x_cirtemp+y_cirtemp*y_cirtemp)-r_temp;
02070 double pt_temp=r_temp/333.567;
02071 if(m_debug>0){
02072 cout<<" in least fit : "<<endl;
02073 cout<<"rtemp: "<<r_temp<<" "<<endl;
02074 cout<<"xtemp: "<<x_cirtemp<<" "<<endl;
02075 cout<<"ytemp: "<<y_cirtemp<<" "<<endl;
02076 cout<<"d0temp: "<<d0_temp<<" "<<endl;
02077 cout<<"pt_temp: "<<pt_temp<<endl;
02078 }
02079
02080 d0=d0_temp;
02081
02082 phi0=atan2(y_cirtemp,x_cirtemp)+M_PI/2.;
02083 omega=1/r_temp;
02084 double z0=0.;
02085 double tanl=0.;
02086 if(m_debug>0) std::cout<<" before global fit Helix PatPar :"<<d0<<","<<phi0<<","<<omega<<","<<z0<<","<<tanl<< std::endl;
02087 return 0;
02088
02089 }
02090
02091 int HoughValidUpdate::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){
02092 double x_stereo;
02093 double y_stereo;
02094
02095 int stereoId=0;
02096
02097
02098 int nDropDelta=0;
02099 for(int i=0;i<n;i++){
02100 if(vec_truthHit[i]==false) continue;
02101 double k;
02102 double b;
02103 if(vec_slant[i]!=0){
02104 if(vec_x_east[i]!=vec_x_west[i]){
02105 k=(vec_y_east[i]-vec_y_west[i])/(vec_x_east[i]-vec_x_west[i]);
02106 b=-k*vec_x_east[i]+vec_y_east[i];
02107
02108 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);
02109
02110
02111
02112
02113 if(delta<0) {
02114 nDropDelta++;
02115 continue;
02116 }
02117
02118 double x1=(2*x_cirtemp+2*k*y_cirtemp-2*k*b+sqrt(delta))/(2*(k*k+1));
02119 double x2=(2*x_cirtemp+2*k*y_cirtemp-2*k*b-sqrt(delta))/(2*(k*k+1));
02120 double x_middle=(vec_x_east[i]+vec_x_west[i])/2.;
02121
02122 if(abs(x1-x_middle)>abs(x2-x_middle)) {x_stereo=x2;}
02123 else {x_stereo=x1;}
02124 y_stereo=k*x_stereo+b;
02125 }
02126 else{
02127
02128
02129 x_stereo=vec_x_east[i];
02130 double y1=sqrt(r_temp*r_temp-(x_cirtemp-x_stereo)*(x_cirtemp-x_stereo))+y_cirtemp;
02131 double y2=y_cirtemp-sqrt(r_temp*r_temp-(x_cirtemp-x_stereo)*(x_cirtemp-x_stereo));
02132 double y_middle=(vec_y_east[i]+vec_y_west[i])/2.;
02133 if(abs(y1-y_middle)>abs(y2-y_middle)) {y_stereo=y2;}
02134 else{y_stereo=y1;}
02135 }
02136
02137
02138 double theta_temp;
02139 if(x_cirtemp==0||x_stereo-x_cirtemp==0){
02140 theta_temp=0;
02141 }
02142 else{
02143 theta_temp=M_PI-atan2(y_stereo-y_cirtemp,x_stereo-x_cirtemp)+atan2(y_cirtemp,x_cirtemp);
02144 if(theta_temp>2*M_PI){
02145 theta_temp=theta_temp-2*M_PI;
02146 }
02147 if(theta_temp<0){
02148 theta_temp=theta_temp+2*M_PI;
02149 }
02150 }
02151
02152
02153
02154
02155
02156
02157
02158
02159 l.push_back(r_temp*theta_temp);
02160
02161
02162
02163
02164 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]));
02165 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]));
02166 vec_zStereo.push_back(vec_z_west[i]-(vec_z_west[i]-vec_z_east[i])*d1/d2);
02167
02168
02169
02170
02171
02172
02173
02174
02175
02176
02177
02178
02179
02180
02181
02182
02183
02184
02185 stereoId=stereoId+1;
02186
02187 }
02188 else {k=b=0;}
02189 }
02190 cout<<"nDropDelta: "<<nDropDelta<<endl;
02191
02192
02193
02194
02195
02196
02197
02198
02199
02200
02201
02202
02203
02204 return stereoId;
02205 }
02206
02207
02208 void HoughValidUpdate::Linefit(int z_stereoNum,vector<double> vec_zStereo,vector<double> l,double& tanl,double& z0){
02209
02210 double l_sum=0;
02211 double z_sum=0;
02212 double l2_sum=0;
02213 double z2_sum=0;
02214 double lz_sum=0;
02215 for(int i=0;i<z_stereoNum;i++){
02216
02217 l_sum=l_sum+l[i];
02218 z_sum=z_sum+vec_zStereo[i];
02219 l2_sum=l2_sum+l[i]*l[i];
02220
02221 z2_sum=z2_sum+vec_zStereo[i]*vec_zStereo[i];
02222 lz_sum=lz_sum+l[i]*vec_zStereo[i];
02223 }
02224 double b=(l2_sum*z_sum-lz_sum*l_sum)/(z_stereoNum*l2_sum-l_sum);
02225 double k=(z_stereoNum*lz_sum-l_sum*z_sum)/(z_stereoNum*l2_sum-l_sum*l_sum);
02226 z0=b;
02227 tanl=k;
02228 cout<<k<<" "<<b<<endl;
02229 }
02230
02231
02232
02233 int HoughValidUpdate::digiToHots2(MdcDigiVec mdcDigiVec,TrkRecoTrk* newTrack,vector<bool> vec_truthHit,vector<double> vec_flt_truth,uint32_t getDigiFlag, vector<int> vec_slant, vector<double> vec_l_Calcu,vector<double> vec_flt,vector<double> vec_theta_ontrack,vector<int>& vec_hitbeforefit,map<int,double> hit_to_circle2,vector<MdcHit*>& vec_for_clean,double tanl){
02234 TrkHitList* trkHitList = newTrack->hits();
02235
02236 if(m_debug>0) cout<<"MdcDigiVec(in 3d fit): "<<mdcDigiVec.size()<<endl;
02237 MdcDigiVec::iterator iter1 = mdcDigiVec.begin();
02238 int digiId=0;
02239 int nhit_beforefit_temp=0;
02240 for (;iter1 != mdcDigiVec.end(); iter1++, digiId++) {
02241
02242 if(vec_theta_ontrack[digiId]>M_PI) continue;
02243
02244 const MdcDigi* aDigi = *iter1;
02245 MdcHit *hit = new MdcHit(aDigi, m_gm);
02246 vec_for_clean.push_back(hit);
02247
02248 Identifier id = aDigi->identify();
02249 int layer = MdcID::layer(id);
02250 int wire = MdcID::wire(id);
02251 int hittemp=layer*1000+wire;
02252 hit->setCalibSvc(m_mdcCalibFunSvc);
02253 hit->setCountPropTime(true);
02254
02255
02256 int newAmbig = 0;
02257
02258 MdcRecoHitOnTrack temp(*hit, newAmbig, m_bunchT0*1.e9);
02259 MdcHitOnTrack* newhot = &temp;
02260 double fltLen = m_gm->Layer(layer)->rMid();
02261
02262 newhot->setFltLen(vec_flt[digiId]);
02263
02264
02265
02266
02267
02268
02269 double ddCut;
02270 double distcirCut;
02271
02272
02273 if(layer<8) {ddCut=1.0; distcirCut=5;}
02274 else {ddCut=1.0; distcirCut=2;}
02275 int use_in3d=1;
02276 if(hit->driftDist(m_bunchT0,0)>ddCut||vec_truthHit[digiId]==false||fabs(hit_to_circle2[hittemp])>distcirCut) { use_in3d=0; }
02277 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())<<" vec_truth: "<<vec_truthHit[digiId]<<" theta: "<<vec_theta_ontrack[digiId]<<" ddcut?: "<<hit->driftDist(m_bunchT0,0)<<" distocir: "<<hit_to_circle2[hittemp]<<" use:? "<<use_in3d<<endl;
02278 if(use_in3d==0) continue;
02279
02280 vec_hitbeforefit.push_back(hittemp);
02281 trkHitList->appendHot(newhot);
02282 }
02283 if(m_debug>0) std::cout<<std::endl;
02284 return trkHitList->nHit();
02285 }
02286 void HoughValidUpdate:: Multi_array(int binX,int binY){
02287 int **count=new int*[binX];
02288 for(int i=0;i<binX;i++){
02289 count[i]=new int [binY];
02290 }
02291
02292 int ***selectNum=new int**[binX];
02293 for(int i=0;i<binX;i++){
02294 selectNum[i]=new int* [binY];
02295 for(int j=0;j<binY;j++){
02296
02297 selectNum[i][j]=new int [count[i][j]];
02298 }
02299 }
02300 }
02301
02302 void HoughValidUpdate::RhoTheta(int nhit,vector<double> vec_u,vector<double> vec_v,vector<double>& vec_rho,vector<double>& vec_theta,vector< vector<int> >& vec_hitNum,vector<int> vec_digitruth){
02303 for(int i=0;i<nhit;i++){
02304 if(vec_digitruth[i]==0) continue;
02305 for(int j=i+1;j<nhit;j++){
02306 if(vec_digitruth[j]==0) continue;
02307 double k,b,x_cross,y_cross;
02308 double rho_temp,theta_temp;
02309 k=(vec_v[i]-vec_v[j])/(vec_u[i]-vec_u[j]);
02310 b=vec_v[i]-k*vec_u[i];
02311 x_cross=-b/(k+1/k);
02312 y_cross=b/(k*k+1);
02313 rho_temp=sqrt(x_cross*x_cross+y_cross*y_cross);
02314 theta_temp=atan2(y_cross,x_cross);
02315 if(theta_temp<0) {
02316 theta_temp=theta_temp+M_PI;
02317 rho_temp=-rho_temp;
02318 }
02319 vec_rho.push_back(rho_temp);
02320 vec_theta.push_back(theta_temp);
02321 vec_hitNum[0].push_back(i);
02322 vec_hitNum[1].push_back(j);
02323
02324 }
02325 }
02326 int nCross=vec_rho.size();
02327 m_nCross=vec_rho.size();
02328 if(m_nCross>125000) return ;
02329
02330
02331 for(int i=0;i<nCross;i++){
02332 m_rho[i]=vec_rho[i];
02333 m_theta[i]=vec_theta[i];
02334
02335 }
02336 }
02337 void HoughValidUpdate::RhoThetaAxial(vector<int> vec_slant,int nhit,vector<double> vec_u,vector<double> vec_v,vector<double>& vec_rho,vector<double>& vec_theta,vector< vector<int> >& vec_hitNum,vector<int> vec_digitruth){
02338 for(int i=0;i<nhit;i++){
02339 if(vec_digitruth[i]==0) continue;
02340 if(vec_slant[i]!=0) continue;
02341 for(int j=i+1;j<nhit;j++){
02342 if(vec_slant[j]!=0) continue;
02343 if(vec_digitruth[j]==0) continue;
02344 double k,b,x_cross,y_cross;
02345 double rho_temp,theta_temp;
02346 k=(vec_v[i]-vec_v[j])/(vec_u[i]-vec_u[j]);
02347 b=vec_v[i]-k*vec_u[i];
02348 x_cross=-b/(k+1/k);
02349 y_cross=b/(k*k+1);
02350 rho_temp=sqrt(x_cross*x_cross+y_cross*y_cross);
02351 theta_temp=atan2(y_cross,x_cross);
02352 if(theta_temp<0) {
02353 theta_temp=theta_temp+M_PI;
02354 rho_temp=-rho_temp;
02355 }
02356 vec_rho.push_back(rho_temp);
02357 vec_theta.push_back(theta_temp);
02358 vec_hitNum[0].push_back(i);
02359 vec_hitNum[1].push_back(j);
02360
02361 }
02362 }
02363 int nCross=vec_rho.size();
02364 if( m_drawTuple ) m_nCross=vec_rho.size();
02365
02366 if(m_drawTuple){
02367 for(int i=0;i<nCross;i++){
02368 m_rho[i]=vec_rho[i];
02369 m_theta[i]=vec_theta[i];
02370 }
02371 }
02372
02373 }
02374 void HoughValidUpdate::FillHist(TH2D *h1,vector<double> vec_u,vector<double> vec_v,int nhit,vector< vector < vector< int > > > &vec_selectNum){
02375 for(int n=0;n<nhit;n++){
02376
02377 for(int i=0;i<binX;i++){
02378 double binwidth=M_PI/binX;
02379 double binhigh=2*m_maxrho/binY;
02380 double theta=(i+0.5)*binwidth;
02381 double rho=vec_u[n]*cos(theta)+vec_v[n]*sin(theta);
02382 int j=(int)((rho+m_maxrho)/binhigh);
02383 int count=h1->GetBinContent(i+1,j+1);
02384 count+=1;
02385 h1->SetBinContent(i+1,j+1,count);
02386
02387 vec_selectNum[i][j].push_back(n);
02388 }
02389 }
02390 }
02391 void HoughValidUpdate::FillRhoTheta(TH2D *h1 ,vector<double> vec_theta,vector<double> vec_rho){
02392 for(int i=0;i<vec_theta.size();i++){
02393 double binwidth=M_PI/binX;
02394 double binhigh=2*m_maxrho/binY;
02395 int binxNum=vec_theta[i]/binwidth+1;
02396 int binyNum=(vec_rho[i]+m_maxrho)/binhigh+1;
02397 int count =h1->GetBinContent(binxNum,binyNum);
02398 count++;
02399 h1->SetBinContent(binxNum,binyNum,count);
02400 }
02401 }
02402
02403
02404 int HoughValidUpdate::zPosition(MdcHit& h, double t0, double x_cirtemp,double y_cirtemp, double r_temp,int digiId,double& delta_z,double& delta_l,double& l_Calcu_Left,double& l_Calcu_Right,double& z_Calcu_Left,double& z_Calcu_Right,double& z_Calcu,double& l_Calcu,int& i_q) {
02405
02406
02407
02408
02409 int return_d[2];
02410 int return_ok[2];
02411 return_d[0]=1;
02412 return_d[1]=1;
02413 return_ok[0]=1;
02414 return_ok[1]=1;
02415 for(int ambig=-1;ambig<=1;ambig=ambig+2){
02416
02417
02418 HepPoint3D fp ( h.wire()->getWestPoint()->x(),h.wire()->getWestPoint()->y(),h.wire()->getWestPoint()->z());
02419 HepPoint3D rp ( h.wire()->getEastPoint()->x(),h.wire()->getEastPoint()->y(),h.wire()->getEastPoint()->z());
02420
02421
02422 HepVector3D X = 0.5 * (fp + rp);
02423 if(m_debug>0){
02424 std::cout<<"---------- "<<std::endl;
02425 h.print(std::cout);
02426
02427 std::cout<<"fp rp "<<fp<<" "<<rp<<std::endl;
02428 std::cout<<"Xmid "<<X<<std::endl;
02429 }
02430 HepVector3D xx = HepVector3D(X.x(), X.y(), 0.);
02431
02432
02433
02434
02435
02436
02437
02438
02439
02440
02441
02442
02443
02444
02445
02446
02447
02448
02449 HepPoint3D center (x_cirtemp, y_cirtemp, 0. );
02450
02451 HepVector3D yy = center - xx;
02452 HepVector3D ww = HepVector3D(yy.x(), yy.y(), 0.);
02453 double wwmag2 = ww.mag2();
02454 double wwmag = sqrt(wwmag2);
02455
02456
02457
02458 double driftdist = fabs( h.driftDist(t0,ambig) );
02459 HepVector3D lr(driftdist/wwmag * ww.x(),
02460 driftdist/wwmag * ww.y(), 0.);
02461 if (m_debug >4){
02462 std::cout<<"xc "<<x_cirtemp << " yc "<<y_cirtemp<< " drfitdist "<<driftdist<<std::endl;
02463
02464
02465 std::cout<<"lr "<<lr<<" "<<" ambig "<<ambig
02466 <<" left "<<h.driftDist(0,1)
02467 <<" right "<<h.driftDist(0,-1)<<std::endl;
02468 }
02469
02470
02471 HepPoint3D ORIGIN = HepPoint3D(0., 0., 0.);
02472
02473
02474
02475 if (ambig == 0) lr = ORIGIN;
02476
02477 if (i_q> 0){
02478 if (ambig == -1){
02479 lr = - lr;
02480 }
02481 }else{
02482 if (ambig == 1){
02483 lr = - lr;
02484 }
02485 }
02486 X += lr;
02487
02488
02489
02490 HepPoint3D tmp(-9999., -9999., 0.);
02491 HepVector3D x = HepVector3D(X.x(), X.y(), 0.);
02492 HepVector3D w = x - center;
02493
02494 HepVector3D V = h.wire()->zAxis();
02495 HepVector3D v = HepVector3D(V.x(), V.y(), 0.);
02496 double vmag2 = v.mag2();
02497
02498
02499
02500 double wv = w.dot(v);
02501
02502 double d2 = wv * wv - vmag2 * (w.mag2() - r_temp * r_temp);
02503 if (m_debug >4){
02504 std::cout<<"X_fix "<<X<<" center "<<center<<std::endl;
02505 std::cout<<"V "<<V<<std::endl;
02506 std::cout<<"w "<<w<<" v "<<v<<std::endl;
02507 std::cout<<"wv "<<wv<<endl;
02508 cout<<"d2: "<<d2<<endl;
02509 }
02510
02511
02512 if (d2 < 0.) {
02513
02514 if (m_debug>4){
02515 cout<<"d2: "<<d2<<endl;
02516 std::cout << "in MdcSegInfoSterO !!! stereo: 0. > d2 = " << d2 << " "
02517 << ambig << std::endl;
02518 }
02519
02520 if (ambig==-1) return_d[ambig+1]=0;
02521 else return_d[ambig]=0;
02522 continue;
02523 }
02524 double d = sqrt(d2);
02525
02526
02527
02528 double l[2];
02529 l[0] =-1*( (- wv + d) / vmag2 );
02530 l[1] =-1*( (- wv - d) / vmag2 );
02531
02532 double z[2];
02533
02534 z[0] = X.z() - l[0] * V.z();
02535 z[1] = X.z() - l[1] * V.z();
02536
02537
02538
02539
02540
02541
02542
02543
02544
02545
02546
02547
02548 bool ok[2];
02549 ok[0] = true;
02550 ok[1] = true;
02551
02552
02553
02554 if (m_debug >4){
02555 std::cout << "X.z "<<X.z()<<" V.z "<<V.z()<<std::endl;
02556 std::cout << "l0, l1 = " << l[0] << ", " << l[1] << std::endl;
02557 std::cout << "rpz " << rp.z() << " fpz " << fp.z() << std::endl;
02558 std::cout << "z0 " << z[0] << " z1 " << z[1] << std::endl;
02559 }
02560
02561 if (ambig == 0) {
02562 if(m_debug>4)std::cout<< " ambig = 0 " <<std::endl;
02563 if (z[0] > rp.z()+20. || z[0] < fp.z()-20.) {
02564 ok[0] = false;
02565 }
02566 if (z[1] > rp.z()+20. || z[1] < fp.z()-20.) {
02567 ok[1] = false;
02568 }
02569 } else {
02570 if(m_debug>4)std::cout<< " ambig != 0 " <<std::endl;
02571
02572 if (fabs(z[0]/rp.z())>1.05 ) { ok[0] = false; }
02573
02574 if (fabs(z[1]/rp.z())>1.05 ) { ok[1] = false; }
02575 }
02576 if ((! ok[0]) && (! ok[1])) {
02577 if(m_debug>4&&(ambig!=0)&&fabs(z[1]/rp.z())>1.05)std::cout<< " z[1] bad " << std::endl;
02578 if(m_debug>4&&(ambig!=0)&&fabs(z[0]/rp.z())>1.05)std::cout<< " z[0] bad " << std::endl;
02579
02580 if(m_debug>4) {
02581 std::cout<< " z[1] bad " << "rpz " << rp.z() << " fpz " << fp.z()
02582 << "z0 " << z[0] << " z1 " << z[1] << std::endl;
02583 std::cout<< " !ok[0] && !ok[1] return -2" <<std::endl;
02584 }
02585
02586 if (ambig==-1) return_ok[ambig+1]=0;
02587 else return_ok[ambig]=0;
02588 continue;
02589 }
02590 if (( ok[0]) && ( ok[1])) {
02591
02592 }
02593 unsigned best = 0;
02594 if (ok[1]) best = 1;
02595 HepVector3D p[2];
02596 p[0] = x + l[0] * v;
02597 p[1] = x + l[1] * v;
02598 double cosdPhi = - center.dot((p[best] - center).unit()) / center.mag();
02599 double dPhi;
02600 if(fabs(cosdPhi)<=1.0) {
02601 dPhi = acos(cosdPhi);
02602 } else if (cosdPhi>1.0) {
02603 dPhi = 0.0;
02604 } else {
02605 dPhi = pi;
02606 }
02607 double stemp=r_temp*dPhi;
02608 if( m_debug >1) cout<<" ambig : "<<ambig<<" z: "<<z[best]<<" stemp: "<<stemp<<endl;
02609 double delta_temp=0;
02610 if(m_drawTuple && m_useMcInfo) delta_temp =z[best]-m_ztruth[digiId];
02611 if( m_debug >3) cout<<"deltatemp: "<<delta_temp<<endl;
02612
02613
02614
02615
02616
02617
02618 if(ambig==1){
02619 z_Calcu_Left=z[best];
02620 l_Calcu_Left=stemp;
02621
02622 if( m_debug >1 && m_drawTuple &&m_useMcInfo) cout<<"ambig: "<<ambig<<" zbest : "<<z[best]<<" ztruth: "<<m_ztruth[digiId]<<endl;
02623 if( m_debug >1 && m_drawTuple &&m_useMcInfo) cout<<"ambig: "<<ambig<<" lbest : "<<stemp<<" ltruth: "<<m_ltruth[digiId]<<endl;
02624 if( m_useMcInfo && m_drawTuple && m_postruth[digiId]==1) { z_Calcu=z[best]; l_Calcu=stemp; delta_z= delta_temp; delta_l =l_Calcu-m_ltruth[digiId];}
02625 }
02626 if(ambig==-1){
02627 z_Calcu_Right=z[best];
02628 l_Calcu_Right=stemp;
02629
02630 if( m_debug >1 && m_drawTuple &&m_useMcInfo) cout<<"ambig: "<<ambig<<" zbest : "<<z[best]<<" ztruth: "<<m_ztruth[digiId]<<endl;
02631 if( m_debug >1 && m_drawTuple &&m_useMcInfo) cout<<"ambig: "<<ambig<<" lbest : "<<stemp<<" ltruth: "<<m_ltruth[digiId]<<endl;
02632 if( m_useMcInfo && m_drawTuple && m_postruth[digiId]==0) { z_Calcu=z[best]; l_Calcu=stemp; delta_z= delta_temp; delta_l =l_Calcu-m_ltruth[digiId];}
02633 }
02634
02635
02637
02638
02639
02640
02641
02642
02643
02644
02645
02646
02647
02648
02649
02650
02651
02652
02653
02654
02655
02656
02657
02658
02659
02660
02661
02662
02663
02664
02665
02666
02667
02668
02669
02670
02671
02672
02673
02674
02675
02676
02677
02678
02679
02680
02681
02682
02683 }
02684
02685 if(return_d[0]==0&&return_d[1]==0) return -1;
02686 if(return_ok[0]==0&&return_ok[1]==0) return -2;
02687 return 1;
02688 }
02689 void HoughValidUpdate::AmbigChoose(vector< vector< vector<double> > > point,int n_use,vector<int>& ambigList,double& tanl, double& z0){
02690
02691 int n =n_use;
02692 int ambig_list[4];
02693
02694 double Chis_Least;
02695
02696 double Chis[4];
02697 int ambig_combine=0;
02698
02699 double chi_k[4];
02700 double chi_b[4];
02701 double first_x[3];
02702 double first_y[3];
02703 first_x[0]=0;
02704 first_y[0]=0;
02705 vector< vector <int> > select(4,vector<int>() );
02706 for(int icombine1=0;icombine1<2;icombine1++){
02707 for(int icombine2=0;icombine2<2;icombine2++){
02708 int combine=icombine1*2+icombine2;
02709 first_x[1]=point[icombine1][0][n-1];
02710 first_x[2]=point[icombine2][0][n-2];
02711 first_y[1]=point[icombine1][1][n-1];
02712 first_y[2]=point[icombine2][1][n-2];
02713 double k,b;
02714 double chi2=0;
02715 LeastLine(3,first_x,first_y,k,b,chi2);
02716
02717
02718
02719
02720
02721
02722
02723
02724
02725
02726
02727
02728
02729 for(int iHit=0;iHit<n-2;iHit++){
02730 double dist_to_line[2];
02731 for(int j=0;j<2;j++){
02732 double x=point[j][0][iHit];
02733 double y=point[j][1][iHit];
02734 dist_to_line[j]=( abs(k*x-y+b) ) / sqrt( k*k +1);
02735 }
02736 if(dist_to_line[0]>dist_to_line[1]) select[combine].push_back(1);
02737 else select[combine].push_back(0);
02738
02739 }
02740 select[combine].push_back(icombine2);
02741 select[combine].push_back(icombine1);
02742
02743
02744
02745
02746
02747 double *l_new=new double[n+1];
02748 double *z_new=new double[n+1];
02749 for(int i=0;i<n;i++){
02750 int ambig=select[combine][i];
02751 l_new[i]=point[ambig][0][i];
02752 z_new[i]=point[ambig][1][i];
02753 }
02754 l_new[n]=0;
02755 z_new[n]=0;
02756
02757
02758
02759
02760
02761
02762
02763
02764
02765
02766 double k_least,b_least;
02767 double chi2_least=0;
02768 LeastLine(n+1,l_new,z_new,k_least,b_least,chi2_least);
02769 Chis[combine]=chi2_least;
02770 chi_k[combine]=k_least;
02771 chi_b[combine]=b_least;
02772
02773
02774
02775 delete []l_new;
02776 delete []z_new;
02777
02778
02779
02780
02781
02782
02783
02784
02785
02786
02787
02788
02789
02790
02791 }
02792 }
02793
02794
02795
02796
02797 Chis_Least=Chis[0];
02798 for(int i=0;i<4;i++){
02799 if(Chis_Least>=Chis[i]) {Chis_Least=Chis[i]; ambig_combine=i;}
02800 }
02801 for(int iHit=0;iHit<n;iHit++){
02802 ambigList.push_back(select[ambig_combine][iHit]);
02803
02804
02805
02806 }
02807
02808 tanl=chi_k[ambig_combine];
02809 z0=chi_b[ambig_combine];
02810
02811
02812 }
02813 void HoughValidUpdate::Combine(TH2D *h1,int widthX,int widthY,vector< vector < vector<int> > > vec_selectNum,M1 &peak_hit_list,M2 &peak_hit_num,int &peak_track,vector<int> peakList[]){
02814 int m=0;
02815 int n=0;
02816 int addnum=999;
02817 int npeak1=peak_track;
02818 vector<bool> peaktrack(npeak1,true);
02819 while(addnum!=0)
02820 {
02821
02822 addnum=0;
02823
02824
02825
02826
02827
02828 int peakX=peakList[0][n];
02829 int peakY=peakList[1][n];
02830 for (int px=peakX-widthX; px<=peakX+widthX; px++) {
02831 for (int py=peakY-widthY; py<=peakY+widthY; py++) {
02832 int ax;
02833 if (px<0) { ax=px+binX; }
02834 else if (px>=binX) { ax=px-binX; }
02835 else { ax=px; }
02836 if ( (ax!=peakX || py!=peakY) && py>=0 && py<binY) Combine_two_cell(h1,peakX,peakY,ax,py,vec_selectNum,m,peak_hit_list,peak_hit_num);
02837 }
02838 }
02839
02840
02841
02842
02843 for(int i=n+1;i<npeak1;i++){
02844 if(peaktrack[i]==false) continue;
02845 int peaktheta=peakList[0][i];
02846 int peakrho=peakList[1][i];
02847 int peakhitNum=peakList[2][i];
02848 int count_same_hit=0;
02849 for(int j=0;j<peakhitNum;j++){
02850 for(int k=0;k<peak_hit_num[m];k++){
02851 if(vec_selectNum[peaktheta][peakrho][j]==peak_hit_list[m][k]){
02852 count_same_hit++;
02853 }
02854 }
02855 }
02856 if(m_debug>1) cout<<" pro: "<<"peak: "<<i<<" "<<(double)count_same_hit/peakhitNum<<endl;
02857 double f_hit=m_hit_pro;
02858 if(count_same_hit>f_hit*peakhitNum){
02859 peaktrack[i]=false;
02860
02861 }
02862 }
02863 for(int i=n+1;i<npeak1;i++){
02864 if(peaktrack[i]==true){
02865 addnum=i-n;
02866 break;
02867 }
02868 }
02869 if(addnum!=0) m++;
02870 n=n+addnum;
02871 }
02872 int ntrack=0;
02873 for(int i=0;i<npeak1;i++){
02874
02875 if(peaktrack[i]==true){
02876 ntrack++;
02877 }
02878 }
02879 peak_track=ntrack;
02880 }
02881 void HoughValidUpdate::Combine_two_cell(TH2D *h1,int xPeakCell,int yPeakCell,int xAround,int yAround,vector< vector < vector<int> > > vec_selectNum,int m,M1 &peak_hit_list,M2 &peak_hit_num){
02882 int size_of_around=vec_selectNum[xAround][yAround].size();
02883 if(size_of_around==0) return;
02884
02885
02886
02887
02888 int size_of_peak=vec_selectNum[xPeakCell][yPeakCell].size();
02889
02890
02891
02892
02893 if(peak_hit_list[m].size()==0){
02894 for(int inum=0;inum<size_of_peak;inum++){
02895 peak_hit_list[m][inum]=vec_selectNum[xPeakCell][yPeakCell][inum];
02896 }
02897 }
02898 vector<int> vec_more_select;
02899 for (int i =0;i<size_of_around;i++){
02900 int same_vec_with_hitcombin=0;
02901 for(int j= 0;j<peak_hit_list[m].size();j++){
02902 if(vec_selectNum[xAround][yAround][i]==peak_hit_list[m][j]){
02903 same_vec_with_hitcombin=1;
02904 break;
02905 }
02906 }
02907 if(same_vec_with_hitcombin==0) {
02908
02909 vec_more_select.push_back(vec_selectNum[xAround][yAround][i]);
02910 }
02911 }
02912 int peak_hit_list_size=peak_hit_list[m].size();
02913 for(int i=0;i<vec_more_select.size();i++ ){
02914 peak_hit_list[m][i+peak_hit_list_size]=vec_more_select[i];
02915 }
02916 int peak_hit_list_size_new=peak_hit_list[m].size();
02917 peak_hit_num[m]=peak_hit_list_size_new;
02918
02919 for(int i= 0;i<peak_hit_list_size_new;i++){
02920
02921 }
02922
02923
02924
02925
02926
02927
02928
02929
02930
02931
02932
02933
02934
02935
02936
02937
02938
02939
02940 }
02941
02942
02943
02944 bool compare(const HoughValidUpdate::Mytrack& a,const HoughValidUpdate::Mytrack& b){
02945 return abs(a.pt)>abs(b.pt);
02946 }
02947 int HoughValidUpdate::LeastLine(int nhit,double x[],double y[],double &k,double &b,double& chi2){
02948 double x_sum=0;
02949 double y_sum=0;
02950 double x2_sum=0;
02951 double y2_sum=0;
02952 double xy_sum=0;
02953 for(int i=0;i<nhit;i++){
02954 x_sum=x_sum+x[i];
02955 y_sum=y_sum+y[i];
02956 x2_sum=x2_sum+x[i]*x[i];
02957 y2_sum=y2_sum+y[i]*y[i];
02958 xy_sum=xy_sum+x[i]*y[i];
02959 }
02960 b=(x2_sum*y_sum-xy_sum*x_sum)/(nhit*x2_sum-x_sum*x_sum);
02961 k=(nhit*xy_sum-x_sum*y_sum)/(nhit*x2_sum-x_sum*x_sum);
02962
02963 double yE[3];
02964 for(int i=0;i<nhit;i++){
02965 yE[i]=k*x[i]+b;
02966 double chi2_temp;
02967 if(yE[i]!=0) chi2_temp=(y[i]-yE[i])*(y[i]-yE[i])/yE[i]*yE[i];
02968 else chi2_temp=0;
02969 chi2+=chi2_temp;
02970 }
02971
02972 return 1;
02973 }