4#include "GaudiKernel/MsgStream.h"
6#include "GaudiKernel/SmartDataPtr.h"
7#include "GaudiKernel/IPartPropSvc.h"
8#include "GaudiKernel/INTupleSvc.h"
11#include "HepPDT/ParticleDataTable.hh"
28 Algorithm(name, pSvcLocator)
30 myMdcCalibFunSvc=NULL;
34 declareProperty(
"Ntuple", myNtProd=0);
35 declareProperty(
"driftTimeUpLimit", myDriftTimeUpLimit = 400);
36 declareProperty(
"MdcHitChi2Cut", myMdcHitChi2Cut = 100);
37 declareProperty(
"ChiCut_circle", myChiCut_circle=5);
38 declareProperty(
"NmaxDeact_circle", myNmaxDeact_circle=1);
39 declareProperty(
"ChiCut_helix", myChiCut_helix=5);
40 declareProperty(
"NmaxDeact_helix", myNmaxDeact_helix=1);
41 declareProperty(
"Debug", myDebug=0);
42 declareProperty(
"DebugNb", myDebugNb=0);
43 declareProperty(
"Chi2CutDiverge", myChi2CutDiverge=99999999);
45 declareProperty(
"IniHelix", myIniHelix=1);
46 declareProperty(
"nBinTanl", myNBinTanl= 100);
47 declareProperty(
"nBinDz", myNBinDz = 100);
48 declareProperty(
"tanlRange", myTanlRange = 3.0);
49 declareProperty(
"dzRange", myDzRange = 50);
50 declareProperty(
"useWireCrossPoint", myUseWireCrossPoint =
false);
51 declareProperty(
"wireCrossPointAvgMode", myWireCrossPointAvgMode =
false);
52 declareProperty(
"WireCrossPointPersistent",myWireCrossPointPersistent =
false);
53 declareProperty(
"IdealTracking",myRunIdealTracking=
false);
54 declareProperty(
"UseInTrkerHits",myUseInTrkerHits =
true);
55 declareProperty(
"TrkFollow",myRunTrkFollow=
true);
65 MsgStream log(
msgSvc(), name());
66 log << MSG::INFO <<
" DotsConnection initialize()" << endreq;
74 StatusCode sc = service (
"MdcGeomSvc", imdcGeomSvc);
75 myMdcGeomSvc =
dynamic_cast<MdcGeomSvc*
> (imdcGeomSvc);
77 return( StatusCode::FAILURE);
91 for(
int i=0; i<nLayer; i++)
94 myNWire[i]=aLayer->
NCell();
95 myRLayer[i]=aLayer->
Radius();
98 double nomShift = aLayer->
nomShift();
99 if(nomShift>0) myWireFlag[i]=1;
100 else if(nomShift<0) myWireFlag[i]=-1;
101 else myWireFlag[i]=0;
102 if(i>0&&myWireFlag[i]!=lastWireFlag) {
105 myAreaLayer[i]=iArea;
106 lastWireFlag=myWireFlag[i];
110 nCellTot+=myNWire[i];
111 for(
int j=0; j<myNWire[i]; j++)
119 int wireidx = aWire->
Id();
120 myMdcDigiGroup[wireidx].
SetWire(aWire);
122 myWirePhi[i][j]=aWire->
Forward().phi();
127 int iInnerLayer = i-1;
129 for(
int k=0; k<myNWire[iInnerLayer]; k++)
132 if(k_last<0) k_last=myNWire[iInnerLayer]-1;
133 double dphi_last = dPhi(myWirePhi[iInnerLayer][k_last], myWirePhi[i][j]);
134 double dphi = dPhi(myWirePhi[iInnerLayer][k], myWirePhi[i][j]);
135 if(dphi_last<0&&dphi>0) {
136 myInnerWire[i][j][0]=k_last;
137 myInnerWire[i][j][1]=k;
150 for(
int i=0; i<nLayer; i++)
152 for(
int j=0; j<myNWire[i]; j++)
154 int iOuterLayer = i+1;
155 if(iOuterLayer<nLayer) {
156 for(
int k=0; k<myNWire[iOuterLayer]; k++)
159 if(k_last<0) k_last=myNWire[iOuterLayer]-1;
160 double dphi_last = dPhi(myWirePhi[iOuterLayer][k_last], myWirePhi[i][j]);
161 double dphi = dPhi(myWirePhi[iOuterLayer][k], myWirePhi[i][j]);
162 if(dphi_last<0&&dphi>0) {
163 myOuterWire[i][j][0]=k_last;
164 myOuterWire[i][j][1]=k;
173 cout<<
"Total "<<nCellTot<<
" cells"<<endl;
180 int n_wire = innerWires_tem.size();
181 cout<<
" Layer "<<layer<<
", wire "<<wire<<
", inner neighbours: "<<endl;
182 for(
int i=0; i<n_wire; i++) {
183 int wireIdx = innerWires_tem[i];
184 layer=myMdcGeomSvc->
Wire(wireIdx)->
Layer();
185 wire=myMdcGeomSvc->
Wire(wireIdx)->
Cell();
186 cout<<
"(L"<<layer<<
", W"<<wire<<
") ";
192 StatusCode Cgem_sc = service (
"CgemGeomSvc", ISvc);
193 if(Cgem_sc.isSuccess()) myCgemGeomSvc=
dynamic_cast<CgemGeomSvc *
>(ISvc);
194 else cout<<
"DotsConnection::initialize(): can not get CgemGeomSvc"<<endl;
198 sc = service (
"RawDataProviderSvc", irawDataProviderSvc);
200 if ( sc.isFailure() ){
201 log << MSG::FATAL << name()<<
" Could not load RawDataProviderSvc!" << endreq;
202 return StatusCode::FAILURE;
207 sc = service(
"MdcCalibFunSvc", imdcCalibSvc);
208 if ( sc.isSuccess() ){
212 cout<<
"DotsConnection::initialize(): can not get MdcCalibFunSvc"<<endl;
220 NTuplePtr nt_ptr(
ntupleSvc(),
"TestDotsHelixFitter/trkPar");
221 if( nt_ptr ) myNtHelixFitter = nt_ptr;
224 myNtHelixFitter =
ntupleSvc()->book(
"TestDotsHelixFitter/trkPar",CLID_ColumnWiseTuple,
"trkPar");
225 if( myNtHelixFitter ) {
226 myNtHelixFitter->addItem (
"run", myRUN);
227 myNtHelixFitter->addItem (
"evt", myEVT);
228 myNtHelixFitter->addItem (
"pid", myPID);
229 myNtHelixFitter->addItem (
"Npar", myNPar, 0,5);
230 myNtHelixFitter->addIndexedItem(
"HelixMC", myNPar, myArrayHelixMC);
231 myNtHelixFitter->addIndexedItem(
"HelixFitted", myNPar, myArrayHelixFitted);
232 myNtHelixFitter->addItem (
"NhitCircle", myNHitsCircle, 0,100);
233 myNtHelixFitter->addIndexedItem(
"LayerHitsCircle", myNHitsCircle, myLayerHitsCircle);
234 myNtHelixFitter->addIndexedItem(
"ChiHitsCircle", myNHitsCircle, myChiHitsCircle);
235 myNtHelixFitter->addItem (
"Nhit", myNHits, 0,100);
236 myNtHelixFitter->addIndexedItem(
"LayerHits", myNHits, myLayerHits);
237 myNtHelixFitter->addIndexedItem(
"ChiHits", myNHits, myChiHits);
238 myNtHelixFitter->addItem (
"NXhit", myNXHits, 0,100);
239 myNtHelixFitter->addItem (
"NVhit", myNVHits, 0,100);
240 myNtHelixFitter->addItem (
"NXCluster", myNXCluster, 0,100);
241 myNtHelixFitter->addItem (
"NVCluster", myNVCluster, 0,100);
242 myNtHelixFitter->addItem (
"TrkIdBest", myTrkIdBest);
243 myNtHelixFitter->addItem (
"NHitsBestTrk", myNHitsBestTrk);
244 myNtHelixFitter->addItem (
"NSameHitsBestTrk", myNSameHitsBestTrk);
250 if(myNtProd&2 || myNtProd&4) {
251 NTuplePtr nt_ptr(
ntupleSvc(),
"TestDotsHelixFitter/lineTrkSegHough");
252 if( nt_ptr ) myNtTrkSeg = nt_ptr;
255 myNtTrkSeg =
ntupleSvc()->book(
"TestDotsHelixFitter/lineTrkSegHough",CLID_ColumnWiseTuple,
"hough");
257 myNtTrkSeg->addItem (
"run", myRUN);
258 myNtTrkSeg->addItem (
"evt", myEVT);
259 myNtTrkSeg->addItem (
"nhits", myNt_nhits, 0,100);
260 myNtTrkSeg->addIndexedItem(
"layer", myNt_nhits, myNt_layer);
261 myNtTrkSeg->addIndexedItem(
"wire", myNt_nhits, myNt_wire);
262 myNtTrkSeg->addIndexedItem(
"xhit", myNt_nhits, myNt_Xhit);
263 myNtTrkSeg->addIndexedItem(
"yhit", myNt_nhits, myNt_Yhit);
264 myNtTrkSeg->addIndexedItem(
"DDhit", myNt_nhits, myNt_DDhit);
265 myNtTrkSeg->addItem (
"rho", myNt_rho);
266 myNtTrkSeg->addItem (
"theta", myNt_theta);
280 cout<<
"DotsConnection::initialize() myIniHelix = "<<myIniHelix<<endl;
283 myRoughTanlDzMap = TH2D(
"roughTanlDzMap",
"roughTanlDzMap",myNBinTanl,-1.*myTanlRange,myTanlRange,myNBinDz,-1.*myDzRange,myDzRange);
286 myThetaRange=acos(-1);
287 myRoughRhoThetaMap=TH2D(
"roughRhoThetaMap",
"roughRhoThetaMap",100,0,myThetaRange,100,-1.*myRhoRange,myRhoRange);
289 return StatusCode::SUCCESS;
297 MsgStream log(
msgSvc(), name());
298 log << MSG::INFO <<
"DotsConnection execute()" << endreq;
308 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),
"/Event/EventHeader");
310 log << MSG::WARNING <<
"Could not find Event Header" << endreq;
311 return StatusCode::FAILURE;
313 myRun = eventHeader->runNumber();
314 myEvt = eventHeader->eventNumber();
316 cout<<endl<<
"------------------------------ "<<endl;
317 cout<<
"run:"<<myRun<<
" , event: "<<myEvt<<endl;
320 cout<<endl<<
"------------------------------ "<<endl;
321 cout<<
"DotsConnection: run:"<<myRun<<
" , event: "<<myEvt<<endl;
339 bool bookTrkCol = registerRecMdcTrack();
342 cout<<
"DotsConnection::execute(): failed to register RecMdcTrackCol!"<<endl;
343 return StatusCode::FAILURE;
350 myVecTrkCandidates.clear();
355 if(vecDigi.size()<2000) {
357 fillMdcDigiMap(vecDigi);
358 buildMdcDigiNeighbors(myMapMdcDigi, 3, 3);
368 else log << MSG::WARNING <<
"Too Many MDC digis: "<<vecDigi.size() << endreq;
372 if(myRunIdealTracking) {
373 getMcFinalChargedStates();
375 associateDigisToMcParticles();
383 return StatusCode::SUCCESS;
391 MsgStream log(
msgSvc(), name());
392 log << MSG::INFO <<
"DotsConnection finalize()" << endreq;
395 cout<<
"DotsConnection finalize(): " << endl;
401 return StatusCode::SUCCESS;
405double DotsConnection::dPhi(
double phi1,
double phi2)
410 if(dphi<-M_PI||dphi>
M_PI) dphi=atan2(
sin(dphi),
cos(dphi));
420 double pos_x, pos_y, pos_z;
425 IPartPropSvc* p_PartPropSvc;
426 static const bool CREATEIFNOTTHERE(
true);
427 StatusCode PartPropStatus = service(
"PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
428 if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) {
429 cout<<
" Could not initialize Particle Properties Service" << endl;
431 HepPDT::ParticleDataTable* p_particleTable = p_PartPropSvc->PDT();
434 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),
"/Event/MC/McParticleCol");
436 McParticleCol::iterator iter_mc = mcParticleCol->begin();
437 for (;iter_mc != mcParticleCol->end(); iter_mc++ ) {
438 if(!(*iter_mc)->primaryParticle())
continue;
439 int pid = (*iter_mc)->particleProperty();
440 int pid_abs =
abs(pid);
441 if( p_particleTable->particle(pid) ){
442 name = p_particleTable->particle(pid)->name();
443 charge= p_particleTable->particle(pid)->charge();
444 }
else if( p_particleTable->particle(-pid) ){
445 name =
"anti " + p_particleTable->particle(-pid)->name();
446 charge = (-1.)*p_particleTable->particle(-pid)->charge();
448 pos_x = (*iter_mc)->initialPosition().x();
449 pos_y = (*iter_mc)->initialPosition().y();
450 pos_z = (*iter_mc)->initialPosition().z();
451 px = (*iter_mc)->initialFourMomentum().px();
452 py = (*iter_mc)->initialFourMomentum().py();
453 pz = (*iter_mc)->initialFourMomentum().pz();
454 if(pid_abs==11||pid_abs==13||pid_abs==211||pid_abs==321||pid_abs==2212)
break;
460 Hep3Vector p3Truth(px, py, pz);
463 helix_truth.pivot(
origin);
465 cout<<
"DotsConnection::getMCHelix() finds a valid MC particle "<<name<<
" at "<<posTruth<<
" with p3="<<p3Truth<<endl;
466 cout<<
"DotsConnection::getMCHelix() Helix = "<<helix_truth.a()<<endl;
472void DotsConnection::getMcFinalChargedStates()
479 double pos_x, pos_y, pos_z;
485 IPartPropSvc* p_PartPropSvc;
486 static const bool CREATEIFNOTTHERE(
true);
487 StatusCode PartPropStatus = service(
"PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
488 if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) {
489 cout<<
" Could not initialize Particle Properties Service" << endl;
491 HepPDT::ParticleDataTable* p_particleTable = p_PartPropSvc->PDT();
494 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),
"/Event/MC/McParticleCol");
496 McParticleCol::iterator iter_mc = mcParticleCol->begin();
497 if(myDebug) cout<<
"MC charged particles: "<<endl;
498 for (;iter_mc != mcParticleCol->end(); iter_mc++ )
500 int pdgid = (*iter_mc)->particleProperty();
501 if(myDebug) cout<<
"idx, pdg, from "<<(*iter_mc)->trackIndex()<<
", "<<pdgid<<
", "<<((*iter_mc)->mother()).trackIndex()<<endl;
502 if( !( (*iter_mc)->primaryParticle() || (*iter_mc)->decayFromGenerator() || (*iter_mc)->decayInFlight() ) )
continue;
503 int pid = (*iter_mc)->particleProperty();
504 int mpid = ((*iter_mc)->mother()).particleProperty();
505 int pid_abs =
abs(pid);
507 if(pid_abs==11||pid_abs==13||pid_abs==211||pid_abs==321||pid_abs==2212)
509 if( p_particleTable->particle(pid) )
511 name = p_particleTable->particle(pid)->name();
512 charge= p_particleTable->particle(pid)->charge();
514 else if( p_particleTable->particle(-pid) )
516 name =
"anti " + p_particleTable->particle(-pid)->name();
517 charge = (-1.)*p_particleTable->particle(-pid)->charge();
519 pos_x = (*iter_mc)->initialPosition().x();
520 pos_y = (*iter_mc)->initialPosition().y();
521 pos_z = (*iter_mc)->initialPosition().z();
522 px = (*iter_mc)->initialFourMomentum().px();
523 py = (*iter_mc)->initialFourMomentum().py();
524 pz = (*iter_mc)->initialFourMomentum().pz();
527 Hep3Vector p3Truth(px, py, pz);
529 helix_truth.pivot(
origin);
531 double dphi=helix_truth.IntersectCylinder(myRLayer[42]/10.);
532 if(dphi==0) dphi=
M_PI*charge;
534 double lenOuter = helix_truth.flightLength(posOuter);
535 double lenInner = helix_truth.flightLength(posTruth);
536 double lenInMdc = lenOuter-lenInner;
537 myVecTrkLenFirstHalf.push_back(lenInMdc);
539 int trkIdx = (*iter_mc)->trackIndex();
540 int mtrkIdx = ((*iter_mc)->mother()).trackIndex();
541 myVecMCTrkId.push_back(trkIdx);
542 myVecPDG.push_back(pid);
543 vector<double> helix;
544 helix.push_back(helix_truth.dr());
545 helix.push_back(helix_truth.phi0());
546 helix.push_back(helix_truth.kappa());
547 helix.push_back(helix_truth.dz());
548 helix.push_back(helix_truth.tanl());
549 myVecHelix.push_back(helix);
552 cout<<Form(
"trk idx %d, pdg code %d, mother %d", trkIdx, pid, mpid)<<endl;
553 cout<<Form(
"p3 = (%.6f, %.6f, %.6f)", px, py, pz)<<endl;
554 cout<<Form(
"pos = (%.6f, %.6f, %.6f)", pos_x, pos_y, pos_z)<<endl;
555 cout<<Form(
"dphi = %.6f", dphi)<<endl;
556 cout<<Form(
"lenOuter = %.6f, lenInner = %.6f", lenOuter,lenInner)<<endl;
557 cout<<
"******************************************************"<<endl;
572void DotsConnection::resetFCSVec()
574 myVecMCTrkId.clear();
576 myVecTrkLenFirstHalf.clear();
578 myVecCgemMcHit.clear();
579 myVecCgemXcluster.clear();
580 myVecCgemVcluster.clear();
583void DotsConnection::associateDigisToMcParticles()
588 double T0 = getEventStartTime();
591 SmartDataPtr<Event::CgemMcHitCol> cgemMcHitCol(eventSvc(),
"/Event/MC/CgemMcHitCol");
592 if(!cgemMcHitCol) cout<<
"DotsConnection::associateDigisToMcParticles() does not find CgemMcHitCol!"<<endl;
596 RecCgemClusterCol::iterator iter_cluster_begin;
597 RecCgemClusterCol::iterator iter_cluster;
598 SmartDataPtr<RecCgemClusterCol> aCgemClusterCol(eventSvc(),
"/Event/Recon/RecCgemClusterCol");
599 if(!aCgemClusterCol) cout<<
"DotsConnection::associateDigisToMcParticles() does not find RecCgemClusterCol!"<<endl;
600 else iter_cluster_begin=aCgemClusterCol->begin();
604 SmartDataPtr<Event::MdcMcHitCol> mdcMcHitCol(eventSvc(),
"/Event/MC/MdcMcHitCol");
605 if(!mdcMcHitCol) cout<<
"DotsConnection::associateDigisToMcParticles() does not find MdcMcHitCol!"<<endl;
608 cout<<
"MDC MC hits obtained, n="<<mdcMcHitCol->size()<<endl;
609 Event::MdcMcHitCol::iterator iter_mdcMcHit = mdcMcHitCol->begin();
610 for(; iter_mdcMcHit!=mdcMcHitCol->end(); iter_mdcMcHit++ )
612 string creatorProcess = (*iter_mdcMcHit)->getCreatorProcess();
613 int trkId = (*iter_mdcMcHit)->getTrackIndex();
614 int isSec = (*iter_mdcMcHit)->getIsSecondary();
615 double trkLen = (*iter_mdcMcHit)->getFlightLength();
616 double px=(*iter_mdcMcHit)->getMomentumX();
617 double py=(*iter_mdcMcHit)->getMomentumY();
618 double pz=(*iter_mdcMcHit)->getMomentumZ();
619 double posx=(*iter_mdcMcHit)->getPositionX();
620 double posy=(*iter_mdcMcHit)->getPositionY();
621 double posz=(*iter_mdcMcHit)->getPositionZ();
622 int pdgCode = (*iter_mdcMcHit)->getCurrentTrackPID();
623 int digiIdx = (*iter_mdcMcHit)->getDigiIdx();
627 Hep3Vector pos(posx,posy,0);
628 Hep3Vector p3(px,py,0);
629 double dotProd = pos*p3;
632 cout<<
" MDC MC hit from "<<creatorProcess
635 <<
", len="<<trkLen/10
636 <<
", p3=("<<px<<
", "<<py<<
", "<<pz<<
")"
637 <<
", pos=("<<posx<<
", "<<posy<<
", "<<posz<<
")"
638 <<
", layer "<<layer<<
" wire "<<wire
639 <<
", dotProd="<<dotProd
640 <<
", pdg code = "<<pdgCode
641 <<
", digi index = "<<digiIdx
647 uint32_t getDigiFlag = 0;
650 cout<<
"DotsConnection::associateDigisToMcParticles() get "<<mdcDigiVec.size()<<
" MDC digis from RawDataProviderSvc"<<endl;
652 clearMdcDigiPointer();
653 vector<MdcDigi*>::iterator iter_mdcDigi = mdcDigiVec.begin();
654 for(; iter_mdcDigi!=mdcDigiVec.end(); iter_mdcDigi++)
663 double tprop = myMdcCalibFunSvc->
getTprop(layer, 0);
664 double charge = (*iter_mdcDigi)->getChargeChannel();
665 double T0Walk = myMdcCalibFunSvc->
getT0(layer,wire) + myMdcCalibFunSvc->
getTimeWalk(layer, charge);
667 double driftT = rawTime - T0 - TOF - T0Walk - tprop;
668 if(driftT>myDriftTimeUpLimit)
continue;
670 myMdcDigiPointer[layer][wire]=*iter_mdcDigi;
675 int nMcTrk = myVecMCTrkId.size();
676 for(
int i=0; i<nMcTrk; i++)
679 int trkIdx = myVecMCTrkId[i];
680 double trkLenInMdc = myVecTrkLenFirstHalf[i];
682 cout<<endl<<
"* MC particle "<<trkIdx<<
", pdg code "<<myVecPDG[i]<<endl;
685 double trkLenMcHit_min=9999.;
686 double pos_innermost[3],p3_innermost[3];
689 if(myUseInTrkerHits&&aCgemClusterCol) {
690 int CgemCluster[3][2]={0,0,0,0,0,0};
691 myVecCgemXcluster.clear();
692 myVecCgemVcluster.clear();
693 myVecCgem1DCluster.clear();
695 Event::CgemMcHitCol::iterator iter_CgemMcHit = cgemMcHitCol->begin();
696 for(; iter_CgemMcHit!=cgemMcHitCol->end(); iter_CgemMcHit++ )
698 string creatorProcess = (*iter_CgemMcHit)->GetCreatorProcess();
700 cout<<
" CGEM MC hit from process "<<creatorProcess;
701 if(creatorProcess==
"Generator"||creatorProcess==
"Decay")
704 if((*iter_CgemMcHit)->GetTrackID()!=trkIdx) {
705 if(myDebug) cout <<
" trkIdx diff "<<endl;
709 if((*iter_CgemMcHit)->GetIsSecondary()) {
710 if(myDebug) cout <<
" is secondary "<<endl;
713 double px=(*iter_CgemMcHit)->GetMomentumXOfPrePoint();
714 double py=(*iter_CgemMcHit)->GetMomentumYOfPrePoint();
715 double pz=(*iter_CgemMcHit)->GetMomentumZOfPrePoint();
716 double posx=(*iter_CgemMcHit)->GetPositionXOfPrePoint();
717 double posy=(*iter_CgemMcHit)->GetPositionYOfPrePoint();
718 double posz=(*iter_CgemMcHit)->GetPositionZOfPrePoint();
719 Hep3Vector pos(posx,posy,0);
720 Hep3Vector p3(px,py,0);
721 double dotProd = pos*p3;
725 if(myDebug) cout <<
" coming back "<<endl;
728 double trkLenMcHit = ((*iter_CgemMcHit)->GetFlightLengthPrePoint())/10.;
729 if(myDebug) cout <<
" trkLenMcHit, trkLenInMdc ="<<trkLenMcHit<<
", "<<trkLenInMdc<<endl;
730 if(trkLenMcHit>1.5*trkLenInMdc) {
731 if(myDebug) cout <<
" trkLenMcHit>1.5*trkLenInMdc "<<endl;
735 if((pos.perp()/10.)<trkLenMcHit_min) {
736 trkLenMcHit_min=pos.perp()/10.;
737 pos_innermost[0]=posx;
738 pos_innermost[1]=posy;
739 pos_innermost[2]=posz;
745 HepPoint3D aPivot(posx/10., posy/10., posz/10.);
746 Hep3Vector aP3(px/1000., py/1000., pz/1000.);
747 double charge=myVecHelix[i][2]/fabs(myVecHelix[i][2]);
751 cout<<
"R="<<pos.perp()/10.<<
", z="<<posz/10.<<
", helix = ("<<aHelix.dr()<<
", "<<aHelix.phi0()<<
", "<<aHelix.kappa()<<
", "<<aHelix.dz()<<
", "<<aHelix.tanl()<<
")"<<endl;
754 const vector<int> & vec_xCluster = (*iter_CgemMcHit)->GetXclusterIdxVec();
755 int nXCluster=vec_xCluster.size();
756 for(
int j=0; j<nXCluster; j++)
758 iter_cluster=iter_cluster_begin+vec_xCluster[j];
759 int clusterId = (*iter_cluster)->getclusterid();
760 int layer=(*iter_cluster)->getlayerid();
761 int sheet=(*iter_cluster)->getsheetid();
763 cout<<
" find one X-cluster on layer "<<layer<<
" (phi="<<(*iter_cluster)->getrecphi()<<
") ";
764 if(CgemCluster[layer][0]==0)
766 myVecCgemXcluster.push_back(*iter_cluster);
767 myVecCgem1DCluster.push_back(*iter_cluster);
768 myVecCgemXCluIdx[layer][sheet].push_back(clusterId);
770 CgemCluster[layer][0]++;
772 cout<<
", associated. ";
776 const vector<int> & vec_vCluster = (*iter_CgemMcHit)->GetVclusterIdxVec();
777 int nVCluster=vec_vCluster.size();
778 for(
int j=0; j<nVCluster; j++)
780 iter_cluster=iter_cluster_begin+vec_vCluster[j];
781 int clusterId = (*iter_cluster)->getclusterid();
782 int layer=(*iter_cluster)->getlayerid();
783 int sheet=(*iter_cluster)->getsheetid();
785 cout<<
" find one V-cluster on layer "<<layer;
786 if(CgemCluster[layer][1]==0)
788 myVecCgemVcluster.push_back(*iter_cluster);
789 myVecCgem1DCluster.push_back(*iter_cluster);
790 myVecCgemVCluIdx[layer][sheet].push_back(clusterId);
791 CgemCluster[layer][1]++;
793 cout<<
", associated. ";
797 if(myDebug) cout<<endl;
801 myVecMdcDigi.clear();
802 int nMdcXHits(0), nMdcVHits(0);
822 Event::MdcMcHitCol::iterator iter_mdcMcHit = mdcMcHitCol->begin();
824 for(; iter_mdcMcHit!=mdcMcHitCol->end(); iter_mdcMcHit++ )
826 int trkId = (*iter_mdcMcHit)->getTrackIndex();
827 if(trkId!=trkIdx)
continue;
828 if((*iter_mdcMcHit)->getDigiIdx()==-9999)
continue;
831 if(!myUseInTrkerHits&&layer<8)
continue;
833 int isSec = (*iter_mdcMcHit)->getIsSecondary();
840 double trkLenMcHit = ((*iter_mdcMcHit)->getFlightLength())/10.;
841 if(trkLenMcHit>1.5*trkLenInMdc) {
845 double px=(*iter_mdcMcHit)->getMomentumX();
846 double py=(*iter_mdcMcHit)->getMomentumY();
847 double pz=(*iter_mdcMcHit)->getMomentumZ();
848 double posx=(*iter_mdcMcHit)->getPositionX();
849 double posy=(*iter_mdcMcHit)->getPositionY();
850 double posz=(*iter_mdcMcHit)->getPositionZ();
851 Hep3Vector pos(posx,posy,0);
852 Hep3Vector p3(px,py,0);
853 double dotProd = pos*p3;
860 if(trkLenMcHit<trkLenMcHit_min) {
861 trkLenMcHit_min=trkLenMcHit;
862 pos_innermost[0]=posx;
863 pos_innermost[1]=posy;
864 pos_innermost[2]=posz;
870 HepPoint3D aPivot(posx/10., posy/10., posz/10.);
871 Hep3Vector aP3(px/1000., py/1000., pz/1000.);
872 double charge=myVecHelix[i][2]/fabs(myVecHelix[i][2]);
876 cout<<
"R="<<pos.perp()/10.<<
", z="<<posz/10.<<
", helix = ("<<aHelix.dr()<<
", "<<aHelix.phi0()<<
", "<<aHelix.kappa()<<
", "<<aHelix.dz()<<
", "<<aHelix.tanl()<<
") layer "<<layer<<endl;
881 const MdcDigi* aMdcDigiPt = myMdcDigiPointer[layer][wire];
884 myVecMdcDigi.push_back(aMdcDigiPt);
885 myMdcDigiPointer[layer][wire]=NULL;
887 cout<<
" MDC digi on layer "<<layer<<
" wire "<<wire<<
" associated"<<endl;
888 if(layer<8||(layer>19&&layer<36)) nMdcVHits++;
892 if(myDebug) cout<<
" digi association finished"<<endl;
899 if((myVecCgemVcluster.size()+nMdcVHits)>=2&&(myVecCgemXcluster.size()+nMdcXHits+myVecCgemVcluster.size()+nMdcVHits)>5)
903 myDotsHelixFitter.
setDChits(myVecMdcDigi,T0);
906 HepVector a_helix(5,0);
909 a_helix(1)=myVecHelix[i][0];
910 a_helix(2)=myVecHelix[i][1];
911 a_helix(3)=myVecHelix[i][2];
912 a_helix(4)=myVecHelix[i][3];
913 a_helix(5)=myVecHelix[i][4];
915 else if(myIniHelix==2)
917 HepPoint3D posTruth_mdc(pos_innermost[0]/10., pos_innermost[1]/10., pos_innermost[2]/10.);
918 Hep3Vector p3Truth(p3_innermost[0]/1000., p3_innermost[1]/1000., p3_innermost[2]/1000.);
919 double charge=myVecHelix[i][2]/fabs(myVecHelix[i][2]);
922 a_helix(1)=ini_helix.dr();
923 a_helix(2)=ini_helix.phi0();
924 a_helix(3)=ini_helix.kappa();
925 a_helix(4)=ini_helix.dz();
926 a_helix(5)=ini_helix.tanl();
928 cout<<
"helix at IP: "<<myVecHelix[i][0]<<
", "<<myVecHelix[i][1]<<
", "<<myVecHelix[i][2]<<
", "<<myVecHelix[i][3]<<
", "<<myVecHelix[i][4]<<endl;
929 cout<<
"helix at innermost hit: "<<a_helix(1)<<
", "<<a_helix(2)<<
", "<<a_helix(3)<<
", "<<a_helix(4)<<
", "<<a_helix(5)<<endl;
930 cout<<
"innermost hit position: r="<<posTruth_mdc.perp()<<
", phi="<<posTruth_mdc.phi()<<endl;
933 else if(myIniHelix==3)
939 HepPoint3D posTruth_mdc(pos_innermost[0]/10., pos_innermost[1]/10., pos_innermost[2]/10.);
940 Hep3Vector p3Truth(p3_innermost[0]/1000., p3_innermost[1]/1000., p3_innermost[2]/1000.);
941 double charge=myVecHelix[i][2]/fabs(myVecHelix[i][2]);
944 if(myVecCgemXcluster.size()+nMdcXHits>=3) {
946 a_helix(1)=circleParVec[0];
947 a_helix(2)=circleParVec[1];
948 a_helix(3)=circleParVec[2];
951 a_helix(1)=ini_helix.dr();
952 a_helix(2)=ini_helix.phi0();
953 a_helix(3)=ini_helix.kappa();
955 a_helix(4)=ini_helix.dz();
956 a_helix(5)=ini_helix.tanl();
977 if(myNtProd&1 && nIter==1) {
981 vector<const RecCgemCluster*>::iterator it_cluster=myVecCgem1DCluster.begin();
983 for(; it_cluster!=myVecCgem1DCluster.end(); it_cluster++, i_cluster++)
985 int flag=(*it_cluster)->getflag();
986 if(flag!=0)
continue;
987 int layer=(*it_cluster)->getlayerid();
989 if(myNHitsCircle<100)
991 myLayerHitsCircle[myNHitsCircle]=layer;
992 myChiHitsCircle[myNHitsCircle] =vecChiCgem[i_cluster];
998 vector<const MdcDigi*>::iterator it_digi=myVecMdcDigi.begin();
1002 for(; it_digi!=myVecMdcDigi.end(); it_digi++, i_digi++)
1007 if(myWireFlag[layer]!=0)
continue;
1008 if(myNHitsCircle<100)
1010 myLayerHitsCircle[myNHitsCircle]=layer;
1011 myChiHitsCircle[myNHitsCircle] =vecChiMdc[i_digi];
1018 if(myDotsHelixFitter.
deactiveHits(myChiCut_circle,myNmaxDeact_circle)==0)
break;
1019 if(nIter>100)
break;
1022 cout<<
"initial helix "<<ini_helix.a()<<endl;
1023 cout<<
"fitFlag="<<fitFlag<<endl;
1024 cout<<
"nIter="<<nIter<<endl;
1025 cout<<
"fitted circle "<<myDotsHelixFitter.
getHelix()<<endl;
1026 cout<<
"circle chi2="<<myDotsHelixFitter.
getChi2()<<endl;
1027 cout<<
"circle error="<<myDotsHelixFitter.
getEa()<<endl<<endl;
1037 if(fitFlag!=0)
break;
1038 if(myDotsHelixFitter.
deactiveHits(myChiCut_helix,myNmaxDeact_helix)==0)
break;
1043 if(fitFlag!=0)
break;
1044 if(myDotsHelixFitter.
activeHits(myChiCut_helix)==0)
break;
1050 cout<<
"fitFlag="<<fitFlag<<endl;
1051 cout<<
"fitted helix "<<myDotsHelixFitter.
getHelix()<<endl;
1052 cout<<
"helix chi2="<<myDotsHelixFitter.
getChi2()<<endl;
1053 cout<<
"helix error="<<myDotsHelixFitter.
getEa()<<endl;
1059 cout<<
" Didn't found enough hits!!!!"<<endl;
1076 myNXCluster=myVecCgemXcluster.size();
1077 myNVCluster=myVecCgemVcluster.size();
1079 myArrayHelixMC[0]=myVecHelix[i][0];
1081 myArrayHelixMC[1]=myVecHelix[i][1];
1082 myArrayHelixMC[2]=myVecHelix[i][2];
1083 myArrayHelixMC[3]=myVecHelix[i][3];
1084 myArrayHelixMC[4]=myVecHelix[i][4];
1087 myNSameHitsBestTrk=0;
1092 saveARecMdcTrack(i);
1097 vector<const RecCgemCluster*>::iterator it_cluster=myVecCgem1DCluster.begin();
1099 for(; it_cluster!=myVecCgem1DCluster.end(); it_cluster++, i_cluster++)
1101 int flag=(*it_cluster)->getflag();
1102 int layer=(*it_cluster)->getlayerid();
1103 if(flag==1) layer=-1*(layer+1);
1106 myLayerHits[myNHits]=layer;
1107 myChiHits[myNHits] =vecChiCgem[i_cluster];
1113 vector<const MdcDigi*>::iterator it_digi=myVecMdcDigi.begin();
1117 for(; it_digi!=myVecMdcDigi.end(); it_digi++, i_digi++)
1119 if(vecIsActMdc[i_digi]==0)
continue;
1124 myLayerHits[myNHits]=layer;
1125 myChiHits[myNHits] =vecChiMdc[i_digi];
1134 if(myNtProd&1) myNtHelixFitter->write();
1148 uint32_t getDigiFlag = 0;
1157 cout<<
"DotsConnection::getMdcDigiVec() get "<<mdcDigiVec.size()<<
" MDC digis from RawDataProviderSvc"<<endl;
1160 vector<const MdcDigi*> aMdcDigiVec;
1161 vector<MdcDigi*>::iterator iter_mdcDigi = mdcDigiVec.begin();
1162 for(; iter_mdcDigi!=mdcDigiVec.end(); iter_mdcDigi++) {
1170 aMdcDigiVec.push_back(*iter_mdcDigi);
1184double DotsConnection::getEventStartTime()
1189 SmartDataPtr<RecEsTimeCol> evTimeCol(eventSvc(),
"/Event/Recon/RecEsTimeCol");
1191 RecEsTimeCol::iterator iter_evt = evTimeCol->begin();
1192 if(iter_evt != evTimeCol->end()){
1194 T0 = (*iter_evt)->getTest();
1197 else cout<<
"error: DotsConnection::getEventStartTime() failed to access event start time"<<endl;
1203double DotsConnection::getRoughDD(
const MdcDigi* aMdcDigi)
1211 double tprop = myMdcCalibFunSvc->
getTprop(layer, 0);
1213 double T0Walk = myMdcCalibFunSvc->
getT0(layer,wire) + myMdcCalibFunSvc->
getTimeWalk(layer, charge);
1215 double driftT = rawTime - myEvtT0 - TOF - T0Walk - tprop;
1218 double entranceAngle = 0;
1219 double driftDist = 0.1 * myMdcCalibFunSvc->
driftTimeToDist(driftT,layer,wire,lrCalib,entranceAngle);
1225void DotsConnection::testDotsHelixFitterAllHits()
1235 double T0 = getEventStartTime();
1240 myDotsHelixFitter.
setDChits(vecDigi,T0);
1243 vector<const MdcDigi*>::iterator iter_mdcDigi = vecDigi.begin();
1244 map<double, const MdcDigi*> aMapDcDigi;
1245 for(; iter_mdcDigi!=vecDigi.end(); iter_mdcDigi++)
1247 aMapDcDigi[myDotsHelixFitter.
getFlightLength(*iter_mdcDigi)] = *iter_mdcDigi;
1250 vector<const RecCgemCluster*> vecCgemCluster = getCgemClusterVec();
1258 HepVector a = ini_helix.
a();
1259 HepVector a_new = myDotsHelixFitter.
getHelix();
1260 for(
int i=0; i<5; i++)
1262 myArrayHelixMC[i]=a[i];
1263 myArrayHelixFitted[i]=a_new[i];
1265 myNtHelixFitter->write();
1268bool DotsConnection::getCgemClusterIntersection(
int idx_Xclu,
int idx_Vclu,
double& z)
1270 RecCgemClusterCol::iterator iter_Xclu = myIterCgemClusterBegin+idx_Xclu;
1271 int flag = (*iter_Xclu)->getflag();
1272 if(flag!=0)
return false;
1273 int layer = (*iter_Xclu)->getlayerid();
1274 int sheet = (*iter_Xclu)->getsheetid();
1276 RecCgemClusterCol::iterator iter_Vclu = myIterCgemClusterBegin+idx_Vclu;
1277 int flagV = (*iter_Vclu)->getflag();
1278 if(flagV!=1||(*iter_Vclu)->getlayerid()!=layer||(*iter_Vclu)->getsheetid()!=sheet)
return false;
1280 double phi = (*iter_Xclu)->getrecphi();
1281 double V_loc = (*iter_Vclu)->getrecv();
1292vector<const RecCgemCluster*> DotsConnection::getCgemClusterVec(
int view)
1294 vector<const RecCgemCluster*> aVecCgemCluster;
1295 SmartDataPtr<RecCgemClusterCol> aCgemClusterCol(eventSvc(),
"/Event/Recon/RecCgemClusterCol");
1298 RecCgemClusterCol::iterator iter_cluster=aCgemClusterCol->begin();
1299 int nCluster = aCgemClusterCol->size();
1312 for(; iter_cluster!=aCgemClusterCol->end(); iter_cluster++)
1324 int flag = (*iter_cluster)->getflag();
1325 if(view==0||view==1)
1327 if(flag==view) aVecCgemCluster.push_back(*iter_cluster);
1330 if(flag==0||flag==1) aVecCgemCluster.push_back(*iter_cluster);
1335 else cout<<
"DotsConnection::getCgemClusterVec() does not find RecCgemClusterCol!"<<endl;
1336 return aVecCgemCluster;
1346bool DotsConnection::getCgemClusters()
1348 for(
int i=0; i<3; i++)
1350 myVecCgemXClusterIdx[i].clear();
1351 myVecCgemVClusterIdx[i].clear();
1352 myVecCgemXClusterTrkIdx[i].clear();
1353 myVecCgemVClusterTrkIdx[i].clear();
1356 SmartDataPtr<RecCgemClusterCol> aCgemClusterCol(eventSvc(),
"/Event/Recon/RecCgemClusterCol");
1360 myIterCgemClusterBegin=aCgemClusterCol->begin();
1361 myIterCgemClusterEnd=aCgemClusterCol->end();
1362 myIterClu=myIterCgemClusterBegin;
1363 int nCluster = aCgemClusterCol->size();
1364 for(
int i=0; i<nCluster; i++, myIterClu++)
1366 int flag = (*myIterClu)->getflag();
1367 int layer= (*myIterClu)->getlayerid();
1368 if(flag==0) myVecCgemXClusterIdx[layer].push_back(i);
1369 else if(flag==1) myVecCgemVClusterIdx[layer].push_back(i);
1371 for(
int i=0; i<3; i++) {
1372 int size = myVecCgemXClusterIdx[i].size();
1373 myVecCgemXClusterTrkIdx[i].resize(size, -1);
1374 size = myVecCgemVClusterIdx[i].size();
1375 myVecCgemVClusterTrkIdx[i].resize(size, -1);
1382void DotsConnection::testDotsHelixFitterPartHits()
1392 double T0 = getEventStartTime();
1394 myDotsHelixFitter.
setT0(T0);
1400 vector<const MdcDigi*>::iterator iter_mdcDigi = vecDigi.begin();
1401 map<double, const MdcDigi*> aMapDcDigi;
1402 for(; iter_mdcDigi!=vecDigi.end(); iter_mdcDigi++)
1404 aMapDcDigi[myDotsHelixFitter.
getFlightLength(*iter_mdcDigi)] = *iter_mdcDigi;
1409 vector<const MdcDigi*> smallVecDigi;
1410 map<double, const MdcDigi*>::iterator iter_digi = aMapDcDigi.begin();
1412 const MdcDigi* digiUdStudy=NULL;
1413 for(; iter_digi!=aMapDcDigi.end(); iter_digi++)
1415 const MdcDigi* aDigi = iter_digi->second;
1416 int layer = myDotsHelixFitter.
getLayer(aDigi);
1417 if(layer==layerUdStudy)
1422 if(layer>layerUdStudy&&nLayerCount>=0&&nLayerCount<layerUsed)
1424 smallVecDigi.push_back(aDigi);
1427 if(nLayerCount==layerUsed)
break;
1429 if( digiUdStudy!=NULL && nLayerCount==layerUsed )
1431 myDotsHelixFitter.
setDChits(smallVecDigi,T0);
1454void DotsConnection::clearMdcDigiPointer()
1456 for(
int i=0; i<43; i++)
1457 for(
int j=0; j<288; j++)
1458 myMdcDigiPointer[i][j]=NULL;
1467bool DotsConnection::registerRecMdcTrack()
1469 MsgStream log(
msgSvc(), name());
1471 IDataProviderSvc* eventSvc = NULL;
1472 service(
"EventDataSvc", eventSvc);
1474 log << MSG::INFO <<
"makeTds:event Svc has been found" << endreq;
1476 log << MSG::FATAL <<
"makeTds:Could not find eventSvc" << endreq;
1479 IDataManagerSvc *dataManSvc =
dynamic_cast<IDataManagerSvc*
>(eventSvc);;
1482 myRecMdcTrackCol = NULL;
1483 DataObject *aRecMdcTrackCol;
1484 eventSvc->findObject(
"/Event/Recon/RecMdcTrackCol",aRecMdcTrackCol);
1485 if(aRecMdcTrackCol != NULL)
1487 myRecMdcTrackCol =
dynamic_cast<RecMdcTrackCol*
> (aRecMdcTrackCol);
1493 sc = eventSvc->registerObject(
"/Event/Recon/RecMdcTrackCol", myRecMdcTrackCol);
1494 if(sc.isFailure()) {
1495 log << MSG::FATAL <<
" Could not register RecMdcTrack collection" <<endreq;
1498 log << MSG::INFO <<
"RecMdcTrackCol registered successfully!" <<endreq;
1501 DataObject *aRecMdcHitCol;
1502 eventSvc->findObject(
"/Event/Recon/RecMdcHitCol",aRecMdcHitCol);
1503 if(aRecMdcHitCol != NULL) {
1506 myRecMdcHitCol=
dynamic_cast<RecMdcHitCol*
> (aRecMdcHitCol);
1511 sc = eventSvc->registerObject(
"/Event/Recon/RecMdcHitCol", myRecMdcHitCol);
1512 if(sc.isFailure()) {
1513 log << MSG::FATAL <<
" Could not register RecMdcHit collection" <<endreq;
1516 log << MSG::INFO <<
"RecMdcHitCol registered successfully!" <<endreq;
1526bool DotsConnection::saveARecMdcTrack(
int iTrk)
1532 int trackId = myRecMdcTrackCol->size();
1536 HepVector aHelixVec = myDotsHelixFitter.
getHelix();
1537 helixPar[0]=aHelixVec[0];
1538 helixPar[1]=aHelixVec[1];
1539 helixPar[2]=aHelixVec[2];
1540 helixPar[3]=aHelixVec[3];
1541 helixPar[4]=aHelixVec[4];
1546 myArrayHelixFitted[0]=aHelixVec[0];
1547 myArrayHelixFitted[1]=aHelixVec[1];
1548 myArrayHelixFitted[2]=aHelixVec[2];
1549 myArrayHelixFitted[3]=aHelixVec[3];
1550 myArrayHelixFitted[4]=aHelixVec[4];
1553 int q = helixPar[2]>0? 1:-1;
1555 double pxy = aHelix.
pt();
1556 double px = aHelix.
momentum(0).x();
1557 double py = aHelix.
momentum(0).y();
1558 double pz = aHelix.
momentum(0).z();
1559 double p = aHelix.
momentum(0).mag();
1560 double theta = aHelix.
direction(0).theta();
1564 double r = poca.perp();
1565 HepSymMatrix Ea = aHelix.
Ea();
1567 double errorMat[15];
1569 for (
int ie = 0 ; ie < 5 ; ie ++){
1570 for (
int je = ie ; je < 5 ; je ++){
1571 errorMat[k] = Ea[ie][je];
1575 double chisq = myDotsHelixFitter.
getChi2();
1577 recMdcTrack->
setPxy(pxy);
1578 recMdcTrack->
setPx(px);
1579 recMdcTrack->
setPy(py);
1580 recMdcTrack->
setPz(pz);
1581 recMdcTrack->
setP(p);
1583 recMdcTrack->
setPhi(phi);
1585 recMdcTrack->
setX(poca.x());
1586 recMdcTrack->
setY(poca.y());
1587 recMdcTrack->
setZ(poca.z());
1588 recMdcTrack->
setR(sqrt(poca.x()*poca.x() + poca.y()*poca.y()));
1598 recMdcTrack->
setStat(myVecPDG[iTrk]);
1600 int maxLayerId = -1;
1601 int minLayerId = 43;
1603 double fltLen = -0.00001;
1604 int layerMaxFltLen=-1;
1609 vector<const RecCgemCluster*>::iterator it_cluster=myVecCgem1DCluster.begin();
1611 for(; it_cluster!=myVecCgem1DCluster.end(); it_cluster++, i_cluster++)
1613 int flag=(*it_cluster)->getflag();
1614 int layer=(*it_cluster)->getlayerid();
1615 int sheet=(*it_cluster)->getsheetid();
1616 if(fabs(vecChiCgem[i_cluster]+9999)<1e-10)
1617 cout<<
"vecChiCgem="<<vecChiCgem[i_cluster]<<endl;
1618 if(flag==0) myVecCgemXCluChi[layer][sheet].push_back(vecChiCgem[i_cluster]);
1619 else if(flag==1) myVecCgemVCluChi[layer][sheet].push_back(vecChiCgem[i_cluster]);
1622 map<int,int> clusterFitStat;
1623 SmartDataPtr<RecCgemClusterCol> aCgemClusterCol(eventSvc(),
"/Event/Recon/RecCgemClusterCol");
1624 if(!aCgemClusterCol) cout<<
"DotsConnection::saveARecMdcTrack() does not find RecCgemClusterCol!"<<endl;
1626 RecCgemClusterCol::iterator iter_cluster_begin=aCgemClusterCol->begin();
1627 RecCgemClusterCol::iterator iter_cluster=iter_cluster_begin;
1628 for(; iter_cluster!=aCgemClusterCol->end(); iter_cluster++)
1630 int flag = (*iter_cluster)->getflag();
1631 if(flag!=2)
continue;
1632 int layer=(*iter_cluster)->getlayerid();
1633 int sheet=(*iter_cluster)->getsheetid();
1634 int idXClu = (*iter_cluster)->getclusterflagb();
1636 vector<int>::iterator
iter=myVecCgemXCluIdx[layer][sheet].begin();
1638 for(;
iter!=myVecCgemXCluIdx[layer][sheet].end();
iter++, i_cluster++)
1639 if((*
iter)==idXClu && fabs(myVecCgemXCluChi[layer][sheet][i_cluster]+9999)>1e-10)
1641 if(!matchX)
continue;
1642 int idVClu = (*iter_cluster)->getclusterflage();
1644 iter=myVecCgemVCluIdx[layer][sheet].begin();
1645 for(;
iter!=myVecCgemVCluIdx[layer][sheet].end();
iter++)
1646 if((*
iter)==idVClu) matchV=
true;
1651 clusterRefVec.push_back(recCgemCluster);
1652 clusterFitStat[clusterid] = 1;
1653 if(maxLayerId<layer)
1667 int nMdcHits=aRecMdcHitVec.size();
1669 vector<RecMdcHit>::iterator iter_recMdcHit = aRecMdcHitVec.begin();
1670 for(; iter_recMdcHit!=aRecMdcHitVec.end(); iter_recMdcHit++)
1672 if(iter_recMdcHit->getChisqAdd()>myMdcHitChi2Cut)
1676 recMdcHit->
setId(hitId);
1679 myRecMdcHitCol->push_back(recMdcHit);
1680 SmartRef<RecMdcHit> refHit(recMdcHit);
1681 hitRefVec.push_back(refHit);
1687 setMdcIdt.insert(layer*1000+wire);
1688 if(layer>maxLayerId)
1692 if(layer<minLayerId)
1696 if(fltLen<recMdcHit->getFltLen()) {
1698 layerMaxFltLen=layer;
1703 cout<<
"track "<<trackId<<
", "<<nMdcHitsKept<<
"/"<<nMdcHits<<
" MDC hits kept, "<<clusterRefVec.size()<<
" 2D-clusters kept"<<endl;
1706 if(maxLayerId>=0&&maxLayerId<3) {
1710 if(fltLen>0) fiTerm=-fltLen*
sin(theta)/aHelix.
radius();
1712 cout<<
"fltLen<0!: maxLayerId="<<maxLayerId<<
", n_cluster="<<clusterRefVec.size()<<
", nMdcHitsKept="<<nMdcHitsKept<<endl;
1713 cout<<
"myVecCgem1DCluster.size="<<myVecCgem1DCluster.size()<<endl;
1720 cout<<
"fiTerm, layer, fltLen, pos= "<<fiTerm<<
", "<<layerMaxFltLen<<
", "<<fltLen<<
", "<<posOut<<endl;
1724 std::sort(clusterRefVec.begin(),clusterRefVec.end(),
sortCluster);
1728 int n_CgemCluster(0), n_MdcHits(0), n_hits_tot(0);
1729 int n_sameCgemCluster(0), n_sameMdcHits(0), n_sameHit_tot(0);
1731 RecMdcTrackCol::iterator iter_bestTrk;
1732 RecMdcTrackCol::iterator iter_trk = myRecMdcTrackCol->begin();
1733 for(; iter_trk!=myRecMdcTrackCol->end(); iter_trk++)
1735 int stat =(*iter_trk)->stat();
1738 n_sameCgemCluster=0;
1740 n_CgemCluster=aClusterRefVec.size();
1741 ClusterRefVec::iterator itCluster=aClusterRefVec.begin();
1742 for(; itCluster!=aClusterRefVec.end(); itCluster++) {
1743 int clusterid = (*itCluster)->getclusterid();
1744 if(clusterFitStat.count(clusterid)) n_sameCgemCluster++;
1748 HitRefVec aHitRefVec = (*iter_trk)->getVecHits();
1749 n_MdcHits=aHitRefVec.size();
1750 HitRefVec::iterator it_hit = aHitRefVec.begin();
1751 for(; it_hit!=aHitRefVec.end(); it_hit++) {
1755 if(setMdcIdt.find(layid*1000+localwid)!=setMdcIdt.end()) n_sameMdcHits++;
1757 if((n_sameMdcHits+2*n_sameCgemCluster)>n_sameHit_tot) {
1758 n_sameHit_tot=n_sameMdcHits+2*n_sameCgemCluster;
1759 n_hits_tot=n_MdcHits+2*n_CgemCluster;
1760 iter_bestTrk=iter_trk;
1761 bestTrkId=(*iter_trk)->trackId();
1765 myTrkIdBest=bestTrkId;
1766 myNHitsBestTrk=n_hits_tot;
1767 myNSameHitsBestTrk=n_sameHit_tot;
1770 int stat = (*iter_bestTrk)->stat();
1771 (*iter_bestTrk)->setStat(stat*100000*myVecPDG[iTrk]/
abs(myVecPDG[iTrk])+myVecPDG[iTrk]);
1775 myRecMdcTrackCol->push_back(recMdcTrack);
1780void DotsConnection::clearMdcNeighbours()
1788 for(
int i=0; i<6796; i++)
if(myMdcDigiGroup[i].GetNoNeiborHits()||myMdcDigiGroup[i].HasHit()) myMdcDigiGroup[i].
Reset();
1798void DotsConnection::fillMdcDigiMap(vector<const MdcDigi*>& vecMdcDigiPnt)
1800 myMapMdcDigi.clear();
1801 vector<const MdcDigi*>::iterator iter_mdcDigi = vecMdcDigiPnt.begin();
1802 for(; iter_mdcDigi!=vecMdcDigiPnt.end(); iter_mdcDigi++)
1808 int wireid = myMdcGeomSvc->
Wire(layer_id,wire_id)->
Id();
1811 myMapMdcDigi[wireid]=*iter_mdcDigi;
1812 double dd = getRoughDD(*iter_mdcDigi);
1813 myMapMdcDigiDd[wireid]=dd;
1824void DotsConnection::buildMdcDigiNeighbors(map<int, const MdcDigi*>& aMapMdcDigi,
int SameLRange,
int DiffLRange)
1826 clearMdcNeighbours();
1828 map<int, const MdcDigi*>::iterator iter_mapDigi = aMapMdcDigi.begin();
1829 for(; iter_mapDigi!=aMapMdcDigi.end(); iter_mapDigi++)
1831 int wireid=iter_mapDigi->first;
1832 myMdcDigiGroup[wireid].
AddHit(iter_mapDigi->second);
1840 for(
int range=0;range<SameLRange;range++)
1842 int prevWire = wirePrev_tem;
1843 if(aMapMdcDigi.find(prevWire)!=aMapMdcDigi.end()) myMdcDigiGroup[prevWire].
AddNext(wireid, range);
1848 int nextWire = wireNext_tem;
1849 if(aMapMdcDigi.find(nextWire)!=aMapMdcDigi.end()) myMdcDigiGroup[nextWire].
AddPrev(wireid, range);
1859 for(
int range=0;range<DiffLRange;range++)
1862 vector<int> innerWires = innerWires_tem;
1864 for(
unsigned int i=0;i<innerWires.size();i++){
1866 if(aMapMdcDigi.find(innerWires[i])!=aMapMdcDigi.end()) myMdcDigiGroup[innerWires[i]].
AddOuter(wireid, range);
1869 vector<int> tem2 = innerWires_tem;
1875 for(
unsigned int k=0;k<tem1.size();k++){
1876 for(
unsigned int j=0;j<tem2.size();j++){
1877 if(tem1.at(k)==tem2.at(j)){
1883 innerWires_tem.push_back(tem1.at(k));
1891 vector<int> outerWires = outerWires_tem;
1894 for(
unsigned int i=0;i<outerWires.size();i++){
1896 if(aMapMdcDigi.find(outerWires[i])!=aMapMdcDigi.end()) myMdcDigiGroup[outerWires[i]].
AddInner(wireid, range);
1900 vector<int> tem2 = outerWires_tem;
1906 for(
unsigned int k=0;k<tem1.size();k++){
1907 for(
unsigned int j=0;j<tem2.size();j++){
1908 if(tem1.at(k)==tem2.at(j)){
1914 outerWires_tem.push_back(tem1.at(k));
1926void DotsConnection::fillMdcDigiBuildNeighbors(vector<const MdcDigi*> vecMdcDigiPnt,
int SameLRange,
int DiffLRange)
1928 clearMdcNeighbours();
1933 vector<const MdcDigi*>::iterator iter_mdcDigi = vecMdcDigiPnt.begin();
1934 for(; iter_mdcDigi!=vecMdcDigiPnt.end(); iter_mdcDigi++)
1940 int wireid = myMdcGeomSvc->
Wire(layer_id,wire_id)->
Id();
1942 myMdcDigiGroup[wireid].
AddHit((*iter_mdcDigi));
1943 myMapMdcDigi[wireid]=*iter_mdcDigi;
1952 for(
int range=0;range<SameLRange;range++)
1954 int prevWire = wirePrev_tem;
1955 myMdcDigiGroup[prevWire].
AddNext(wireid);
1960 int nextWire = wireNext_tem;
1961 myMdcDigiGroup[nextWire].
AddPrev(wireid);
1971 for(
int range=0;range<DiffLRange;range++)
1974 vector<int> innerWires = innerWires_tem;
1976 for(
unsigned int i=0;i<innerWires.size();i++){
1978 myMdcDigiGroup[innerWires[i]].
AddOuter(wireid);
1981 vector<int> tem2 = innerWires_tem;
1987 for(
unsigned int k=0;k<tem1.size();k++){
1988 for(
unsigned int j=0;j<tem2.size();j++){
1989 if(tem1.at(k)==tem2.at(j)){
1995 innerWires_tem.push_back(tem1.at(k));
2003 vector<int> outerWires = outerWires_tem;
2006 for(
unsigned int i=0;i<outerWires.size();i++){
2008 myMdcDigiGroup[outerWires[i]].
AddInner(wireid);
2012 vector<int> tem2 = outerWires_tem;
2018 for(
unsigned int k=0;k<tem1.size();k++){
2019 for(
unsigned int j=0;j<tem2.size();j++){
2020 if(tem1.at(k)==tem2.at(j)){
2026 outerWires_tem.push_back(tem1.at(k));
2046 map<int, const MdcDigi*>::iterator iter_mapDigi = myMapMdcDigi.begin();
2047 for(; iter_mapDigi!=myMapMdcDigi.end(); iter_mapDigi++)
2049 int wireIdx=iter_mapDigi->first;
2063 if(myMdcDigiGroup[wireIdx].IsOuterEnd()) myOuterEnds.push_back(wireIdx);
2074void DotsConnection::findOuterEnds()
2076 myOuterEnds.clear();
2077 myNeighbourSeeds.clear();
2079 vector<int> wireIdxVec_noNbLayer;
2095 map<int, const MdcDigi*>::iterator iter_mapDigi = myMapMdcDigi.begin();
2096 for(; iter_mapDigi!=myMapMdcDigi.end(); iter_mapDigi++)
2098 int wireIdx=iter_mapDigi->first;
2112 if(myMdcDigiGroup[wireIdx].IsOuterEnd()) myOuterEnds.push_back(wireIdx);
2114 if((myMdcDigiGroup[wireIdx].GetNTypesNeighborHits()==0&&myMdcDigiGroup[wireIdx].GetNhitsNearby()!=0)||myMdcDigiGroup[wireIdx].isCircuEnd()) wireIdxVec_noNbLayer.push_back(wireIdx);
2118 myNeighbourSeeds=myOuterEnds;
2120 myNeighbourSeeds.insert(myNeighbourSeeds.begin(), wireIdxVec_noNbLayer.begin(), wireIdxVec_noNbLayer.end());
2129void DotsConnection::getMdcHitCluster()
2133 int layer, nAxial, nStereo, nStereo1, nStereo2;
2134 vector<vector<const MdcDigi*> > MdcHitClusters;
2135 vector<int> vecWireId_tmp;
2138 while(!myNeighbourSeeds.empty())
2146 vector<const MdcDigi*> aMdcHitCluster;
2147 vector<int> vecMdcXDigiIdx;
2148 vector<int> vecMdcV1DigiIdx;
2149 vector<int> vecMdcV2DigiIdx;
2150 vector<int> vecMdcVDigiIdx;
2154 int wireIdx = myNeighbourSeeds.back();
2155 myNeighbourSeeds.pop_back();
2165 vecWireId_tmp.push_back(wireIdx);
2166 while(!vecWireId_tmp.empty())
2168 wireIdx=vecWireId_tmp.back();
2169 vecWireId_tmp.pop_back();
2170 if(!myMdcDigiGroup[wireIdx].Used())
2172 layer = myMdcGeomSvc->
Wire(wireIdx)->
Layer();
2173 if(layer<layer_min) layer_min=layer;
2174 if(layer>layer_max) layer_max=layer;
2175 if(myWireFlag[layer]==0) {
2176 aTrkCandi.mdcHitIdx[layer].push_back(nAxial);
2178 vecMdcXDigiIdx.push_back(wireIdx);
2179 if(layer<xlayer_min) xlayer_min=layer;
2180 if(layer>xlayer_max) xlayer_max=layer;
2183 aTrkCandi.mdcHitIdx[layer].push_back(nStereo);
2184 vecMdcVDigiIdx.push_back(wireIdx);
2186 if(myWireFlag[layer]==1) {
2188 vecMdcV1DigiIdx.push_back(wireIdx);
2190 else if(myWireFlag[layer]==-1) {
2192 vecMdcV2DigiIdx.push_back(wireIdx);
2195 aMdcHitCluster.push_back(myMapMdcDigi[wireIdx]);
2199 vector<int> neighbours = myMdcDigiGroup[wireIdx].
GetNeiborHits();
2200 for(
int j=0; j<neighbours.size(); j++) {
2202 if(!myMdcDigiGroup[neighbours[j]].Used()) vecWireId_tmp.push_back(neighbours[j]);
2207 else if(method==1) {
2208 layer = myMdcGeomSvc->
Wire(wireIdx)->
Layer();
2210 vector<int> nbSameLayer;
2211 nbSameLayer.push_back(wireIdx);
2212 int nbWireId=myMdcDigiGroup[wireIdx].
GetPrev();
2213 while(nbWireId>=0) {
2214 if(!myMdcDigiGroup[nbWireId].Used()) nbSameLayer.push_back(nbWireId);
2215 nbWireId=myMdcDigiGroup[nbWireId].
GetPrev();
2217 nbWireId=myMdcDigiGroup[wireIdx].
GetNext();
2218 while(nbWireId>=0) {
2219 if(!myMdcDigiGroup[nbWireId].Used()) nbSameLayer.push_back(nbWireId);
2220 nbWireId=myMdcDigiGroup[nbWireId].
GetNext();
2224 int nSame=nbSameLayer.size();
2226 for(
int i=0; i<nSame; i++)
2229 wireIdx=nbSameLayer[i];
2230 if(!myMdcDigiGroup[wireIdx].Used())
2233 if(layer<layer_min) layer_min=layer;
2234 if(layer>layer_max) layer_max=layer;
2235 if(myWireFlag[layer]==0) {
2236 aTrkCandi.mdcHitIdx[layer].push_back(nAxial);
2238 vecMdcXDigiIdx.push_back(wireIdx);
2239 if(layer<xlayer_min) xlayer_min=layer;
2240 if(layer>xlayer_max) xlayer_max=layer;
2243 aTrkCandi.mdcHitIdx[layer].push_back(nStereo);
2244 vecMdcVDigiIdx.push_back(wireIdx);
2246 if(myWireFlag[layer]==1) {
2248 vecMdcV1DigiIdx.push_back(wireIdx);
2250 else if(myWireFlag[layer]==-1) {
2252 vecMdcV2DigiIdx.push_back(wireIdx);
2255 aMdcHitCluster.push_back(myMapMdcDigi[wireIdx]);
2262 if(layer>=8&&layer<=42) {
2264 set<int> idxSetInnerLayer;
2266 vector<int> seedsInnerLayer;
2267 for(
int j=0; j<nbSameLayer.size(); j++) {
2268 const vector<int>& InnerNbs=myMdcDigiGroup[nbSameLayer[j]].
GetInner();
2269 for(
int k=0; k<InnerNbs.size(); k++) {
2270 wireIdx=InnerNbs[k];
2271 if(myMdcDigiGroup[wireIdx].Used())
continue;
2272 idxSetInnerLayer.insert(wireIdx);
2275 int nbWireId=myMdcDigiGroup[wireIdx].
GetPrev();
2276 while(nbWireId>=0) {
2277 if(!myMdcDigiGroup[nbWireId].Used()) idxSetInnerLayer.insert(nbWireId);
2278 nbWireId=myMdcDigiGroup[nbWireId].
GetPrev();
2280 nbWireId=myMdcDigiGroup[wireIdx].
GetNext();
2281 while(nbWireId>=0) {
2282 if(!myMdcDigiGroup[nbWireId].Used()) idxSetInnerLayer.insert(nbWireId);
2283 nbWireId=myMdcDigiGroup[nbWireId].
GetNext();
2287 nHits=idxSetInnerLayer.size();
2288 set<int>::iterator it=idxSetInnerLayer.begin();
2289 int first(-100), last(-100),
prev(-100);
2290 for(
int iHit=0; iHit<nHits; iHit++, it++) {
2291 if(iHit==0)
first=*it;
2292 if(iHit==nHits-1) last =*it;
2296 seedsInnerLayer.push_back(prev);
2302 if(first==0&&last+1==myNWire[layer]) nGap--;
2303 else seedsInnerLayer.push_back(last);
2305 nbSameLayer.clear();
2306 nbSameLayer.insert(nbSameLayer.end(), idxSetInnerLayer.begin(), idxSetInnerLayer.end());
2309 myNeighbourSeeds.insert(myNeighbourSeeds.end(), seedsInnerLayer.begin(), seedsInnerLayer.end());
2312 cout<<
" layer "<<layer<<
", "<<nHits<<
" hits, "<<nGap<<
" gaps, "<<endl;
2320 if(!aMdcHitCluster.empty()) {
2321 MdcHitClusters.push_back(aMdcHitCluster);
2322 aTrkCandi.mdcXDigiWireIdx=vecMdcXDigiIdx;
2323 int nXMdcDigi = vecMdcXDigiIdx.size();
2324 aTrkCandi.mdcXDigiChi.clear();
2325 aTrkCandi.mdcXDigiChi.resize(nXMdcDigi, 9999.0);
2326 aTrkCandi.mdcXDigiFitted.clear();
2327 aTrkCandi.mdcXDigiFitted.resize(nXMdcDigi, 1);
2328 aTrkCandi.mdcV1DigiWireIdx=vecMdcV1DigiIdx;
2329 aTrkCandi.mdcV2DigiWireIdx=vecMdcV2DigiIdx;
2330 aTrkCandi.mdcVDigiWireIdx=vecMdcVDigiIdx;
2331 int nVMdcDigi = vecMdcVDigiIdx.size();
2332 aTrkCandi.mdcVDigiChi.clear();
2333 aTrkCandi.mdcVDigiChi.resize(nVMdcDigi, 9999.0);
2334 aTrkCandi.mdcVDigiFitted.clear();
2335 aTrkCandi.mdcVDigiFitted.resize(nVMdcDigi, 1);
2336 aTrkCandi.xLayerMin=xlayer_min;
2337 aTrkCandi.xLayerMax=xlayer_max;
2338 aTrkCandi.layerMin=layer_min;
2339 aTrkCandi.layerMax=layer_max;
2340 if((nAxial+nStereo1+nStereo2)<3) aTrkCandi.isGood=0;
2341 else aTrkCandi.isGood=1;
2342 int trkIdx = myVecTrkCandidates.size();
2343 aTrkCandi.trkIdx = trkIdx;
2344 myVecTrkCandidates.push_back(aTrkCandi);
2361 vector<vector<const MdcDigi*> >::iterator
iter = MdcHitClusters.begin();
2363 for(;
iter!=MdcHitClusters.end();
iter++, ipath++) {
2364 cout<<endl<<
" --------> hit-cluster "<<ipath<<
":" <<endl;
2365 vector<const MdcDigi*> vecDigi = *
iter;
2366 vector<const MdcDigi*>::iterator iter_digi = vecDigi.begin();
2368 for(; iter_digi!=vecDigi.end(); iter_digi++, iHit++) {
2373 cout<<
" (layer "<<layer<<
", wire "<<wire<<
", view "<<myWireFlag[layer]<<
") ";
2374 if(iHit%5==0) cout<<endl;
2385void DotsConnection::mdcHitClustering(
const vector<int>& vecWireIdx, vector<vector<int> >& vecHitClu)
2387 int nHit=vecWireIdx.size();
2388 vector<int> vecUseFlag;
2389 for(
int i=0; i<nHit; i++) {
2390 int wireIdx = vecWireIdx.at(i);
2391 vecUseFlag.push_back(myMdcDigiGroup[wireIdx].getUseFlag());
2395 vector<int> vecWireId_tmp;
2396 for(
int i=0; i<nHit; i++)
2398 int wireIdx = vecWireIdx.at(i);
2399 if(myMdcDigiGroup[wireIdx].getUseFlag()==2)
2401 vector<int> aHitClu;
2402 vecWireId_tmp.push_back(wireIdx);
2403 while(!vecWireId_tmp.empty())
2405 wireIdx=vecWireId_tmp.back();
2406 vecWireId_tmp.pop_back();
2407 if(myMdcDigiGroup[wireIdx].getUseFlag()==2)
2408 aHitClu.push_back(wireIdx);
2410 vector<int> neighbours = myMdcDigiGroup[wireIdx].
GetNeiborHits();
2411 for(
int j=0; j<neighbours.size(); j++)
2412 if(myMdcDigiGroup[neighbours[j]].getUseFlag()==2) vecWireId_tmp.push_back(neighbours[j]);
2414 vecHitClu.push_back(aHitClu);
2418 for(
int i=0; i<nHit; i++) {
2419 int wireIdx = vecWireIdx.at(i);
2420 myMdcDigiGroup[wireIdx].
SetUsedFlag(vecUseFlag.at(i));
2428 int nHit=vecWireIdx.size();
2429 int layer, wireIdx, nAxial(0), nStereo(0), nStereo1(0), nStereo2(0);
2434 vector<int> vecMdcXDigiIdx;
2435 vector<int> vecMdcXDigiFitFlag;
2436 vector<int> vecMdcV1DigiIdx;
2437 vector<int> vecMdcV2DigiIdx;
2438 vector<int> vecMdcVDigiIdx;
2439 for(
int i=0; i<nHit; i++) {
2440 wireIdx=vecWireIdx[i];
2441 layer = myMdcGeomSvc->Wire(wireIdx)->Layer();
2442 if(layer<layer_min) layer_min=layer;
2443 if(layer>layer_max) layer_max=layer;
2444 if(myWireFlag[layer]==0) {
2445 aTrkCandi.mdcHitIdx[layer].push_back(nAxial);
2447 vecMdcXDigiIdx.push_back(wireIdx);
2449 if(myMdcDigiGroup[wireIdx].isOnNode()) fitFlag=0;
2450 vecMdcXDigiFitFlag.push_back(fitFlag);
2451 if(layer<xlayer_min) xlayer_min=layer;
2452 if(layer>xlayer_max) xlayer_max=layer;
2455 aTrkCandi.mdcHitIdx[layer].push_back(nStereo);
2456 vecMdcVDigiIdx.push_back(wireIdx);
2458 if(myWireFlag[layer]==1) {
2460 vecMdcV1DigiIdx.push_back(wireIdx);
2462 else if(myWireFlag[layer]==-1) {
2464 vecMdcV2DigiIdx.push_back(wireIdx);
2468 aTrkCandi.mdcXDigiWireIdx=vecMdcXDigiIdx;
2469 int nXMdcDigi = vecMdcXDigiIdx.size();
2470 aTrkCandi.mdcXDigiChi.clear();
2471 aTrkCandi.mdcXDigiChi.resize(nXMdcDigi, 9999.0);
2472 aTrkCandi.mdcXDigiFitted.clear();
2473 aTrkCandi.mdcXDigiFitted.resize(nXMdcDigi, 1);
2475 aTrkCandi.mdcV1DigiWireIdx=vecMdcV1DigiIdx;
2476 aTrkCandi.mdcV2DigiWireIdx=vecMdcV2DigiIdx;
2477 aTrkCandi.mdcVDigiWireIdx=vecMdcVDigiIdx;
2478 int nVMdcDigi = vecMdcVDigiIdx.size();
2479 aTrkCandi.mdcVDigiChi.clear();
2480 aTrkCandi.mdcVDigiChi.resize(nVMdcDigi, 9999.0);
2481 aTrkCandi.mdcVDigiFitted.clear();
2482 aTrkCandi.mdcVDigiFitted.resize(nVMdcDigi, 1);
2483 aTrkCandi.xLayerMin=xlayer_min;
2484 aTrkCandi.xLayerMax=xlayer_max;
2485 aTrkCandi.layerMin=layer_min;
2486 aTrkCandi.layerMax=layer_max;
2487 if((nAxial+nStereo1+nStereo2)<3) aTrkCandi.isGood=0;
2488 else aTrkCandi.isGood=1;
2493void DotsConnection::findCgemTrkSeg()
2496 vector<int> vecCgemXClusterIdx[3];
2497 vector<int> vecCgemVClusterIdx[3];
2498 vector<vector<pair<int, double> > > vecCgemVCluIdxZ[3];
2499 for(
int i=0; i<3; i++) {
2502 size = myVecCgemVClusterIdx[i].size();
2503 for(
int j=0; j<size; j++) {
2504 if(myVecCgemVClusterTrkIdx[i][j]==-1) {
2505 vecCgemVClusterIdx[i].push_back(myVecCgemVClusterIdx[i][j]);
2509 size = myVecCgemXClusterIdx[i].size();
2510 for(
int j=0; j<size; j++) {
2511 if(myVecCgemXClusterTrkIdx[i][j]==-1) {
2512 int idx_Xclu = myVecCgemXClusterIdx[i][j];
2513 vecCgemXClusterIdx[i].push_back(idx_Xclu);
2514 vector<pair<int, double> > vecVidxZ;
2515 for(
int k=0; k<vecCgemVClusterIdx[i].size(); k++) {
2516 int idx_Vclu = vecCgemVClusterIdx[i][k];
2518 if(getCgemClusterIntersection(idx_Xclu,idx_Vclu,z_XV)) {
2519 pair<int, double> idxZ(idx_Vclu, z_XV);
2520 vecVidxZ.push_back(idxZ);
2523 vecCgemVCluIdxZ[i].push_back(vecVidxZ);
2531 for(
int i2=0; i2<vecCgemXClusterIdx[1].size(); i2++) {
2532 int idx2=vecCgemXClusterIdx[1][i2];
2536 for(
int i1=0; i1<vecCgemXClusterIdx[0].size(); i1++) {
2540 int idx1=vecCgemXClusterIdx[0][i1];
2542 for(
int i3=0; i3<vecCgemXClusterIdx[2].size(); i3++) {
2546 int idx3=vecCgemXClusterIdx[2][i3];
2550 aTrkCandi.CgemXClusterID.push_back(idx2);
2551 aTrkCandi.CgemXClusterID.push_back(idx1);
2553 cout<<
"--- Try a triple-CGEM_Xcluster circle fit: "<<idx3<<
", "<<idx2<<
", "<<idx1<<endl;
2554 int nDropHits = circleFitTaubin(aTrkCandi, 5);
2556 if(isGoodCgemCircleCandi(aTrkCandi)) {
2558 cout<<
"a good CGEM cluster circle candidate found"<<endl;
2561 for(
int i=0; i<3; i++) {
2564 s[i] = fabs(dphi*circle_fittedHelix.
radius());
2567 double dz_min=9999.;
2568 double tanl_dzmin=9999.;
2569 int v_dzmin[3]={-1,-1,-1};
2571 for(
int v1=0; v1<vecCgemVCluIdxZ[0][i1].size(); v1++) {
2572 double z1=vecCgemVCluIdxZ[0][i1][v1].second;
2573 for(
int v2=0; v2<vecCgemVCluIdxZ[1][i2].size(); v2++) {
2574 double z2=vecCgemVCluIdxZ[1][i2][v2].second;
2575 for(
int v3=0; v3<vecCgemVCluIdxZ[2][i3].size(); v3++) {
2576 double z3=vecCgemVCluIdxZ[2][i3][v3].second;
2577 double r21=(z2-z1)/(
s[1]-
s[0]);
2578 double r32=(z3-z2)/(
s[2]-
s[1]);
2579 double s_sum(0), z_sum(0);
2580 double s2_sum(0), sz_sum(0);
2581 s_sum=
s[0]+
s[1]+
s[2];
2583 s2_sum=
s[0]*
s[0]+
s[1]*
s[1]+
s[2]*
s[2];
2584 sz_sum=
s[0]*z1+
s[1]*z2+
s[2]*z3;
2585 double tanl_lfit = (n_sz*sz_sum-s_sum*z_sum)/(n_sz*s2_sum-s_sum*s_sum);
2586 double dz_lfit = (s2_sum*z_sum-s_sum*sz_sum)/(n_sz*s2_sum-s_sum*s_sum);
2588 cout<<
"r21, r32, del_r, dz, tanl = "<<r21<<
", "<<r32<<
", "<<fabs(r21-r32)<<
", "<<dz_lfit<<
", "<<tanl_lfit<<endl;
2591 if(fabs(r21-r32)>0.08) {
2592 if(myDebugNb==2) cout<<
" dropped "<<endl;
2595 if(fabs(dz_lfit)<fabs(dz_min)) {
2597 tanl_dzmin=tanl_lfit;
2605 if(v_dzmin[0]!=-1) {
2606 aTrkCandi.CgemVClusterID.push_back(vecCgemVCluIdxZ[2][i3][v_dzmin[2]].first);
2607 aTrkCandi.CgemVClusterID.push_back(vecCgemVCluIdxZ[1][i2][v_dzmin[1]].first);
2608 aTrkCandi.CgemVClusterID.push_back(vecCgemVCluIdxZ[0][i1][v_dzmin[0]].first);
2609 int fitFlag = helixFit(aTrkCandi, 5);
2612 updateCgemFitItem(aTrkCandi, 1);
2614 int trkIdx = myVecTrkCandidates.size();
2615 aTrkCandi.trkIdx = trkIdx;
2616 myVecTrkCandidates.push_back(aTrkCandi);
2629void DotsConnection::printTrkCandi(
struct trkCandi& aTrkCandi)
2631 cout<<
"Trk candidate id "<<aTrkCandi.
trkIdx<<
", good? "<<aTrkCandi.
isGood<<endl;
2633 cout<<
" nXHits="<<nXHits<<
":";
2634 for(
int i=0; i<nXHits; i++) {
2636 if((i+1)%20==0) cout<<endl<<
" ";
2641int DotsConnection::nXHitsActive(
struct trkCandi& aTrkCandi)
2645 for(
int i=0; i<nXHits; i++)
if(aTrkCandi.
mdcXDigiFitted[i]>0) nXActive++;
2653 for(
int i=0; i<nXHits; i++) {
2655 int layer = myMdcGeomSvc->
Wire(wireIdx)->
Layer();
2656 if(layer>=layerMin && layer<=layerMax) {
2664int DotsConnection::setXHitsInBigSetFitFlagUncertain(
struct trkCandi& aTrkCandi,
bool set)
2668 for(
int i=0; i<nXHits; i++) {
2671 if(myMdcDigiGroup[wireIdx].isInBigSet()) {
2708vector<const MdcDigi*> DotsConnection::getXDigiVec(
struct trkCandi& aTrkCandi,
int sel, vector<int>* vecFitFlag)
2712 vector<const MdcDigi*> aXDigiVec;
2719 if(sel==1&&(*
iter)<0)
continue;
2722 aXDigiVec.push_back(myMapMdcDigi[mdcDigiIdx]);
2723 int layer = myMdcGeomSvc->
Wire(mdcDigiIdx)->
Layer();
2724 if(vecFitFlag!=NULL) {
2726 if(fitFlag<0) fitFlag=0;
2727 vecFitFlag->push_back(fitFlag);
2731 myVecDigiIdx.push_back(i);
2743vector<const MdcDigi*> DotsConnection::getVDigiVec(
struct trkCandi& aTrkCandi,
int sel)
2745 vector<const MdcDigi*> aVDigiVec;
2751 if(sel==1&&(*
iter)<0)
continue;
2754 aVDigiVec.push_back(myMapMdcDigi[mdcDigiIdx]);
2755 myVecDigiIdx.push_back(-(i+1));
2760vector<const MdcDigi*> DotsConnection::getDigiVec(
struct trkCandi& aTrkCandi,
int sel)
2762 myVecDigiIdx.clear();
2763 vector<const MdcDigi*> aDigiVec = getXDigiVec(aTrkCandi, sel);
2764 vector<const MdcDigi*> vDigiVec = getVDigiVec(aTrkCandi, sel);
2765 aDigiVec.insert(aDigiVec.end(), vDigiVec.begin(), vDigiVec.end());
2769vector<const RecCgemCluster*> DotsConnection::getXCluVec(
struct trkCandi& aTrkCandi,
int sel)
2771 vector<const RecCgemCluster*> aVecCgemCluster;
2778 if(sel==2)
continue;
2780 myIterClu=myIterCgemClusterBegin+cluIdx;
2781 aVecCgemCluster.push_back(*myIterClu);
2782 myVecCluIdx.push_back(i);
2784 return aVecCgemCluster;
2787vector<const RecCgemCluster*> DotsConnection::getVCluVec(
struct trkCandi& aTrkCandi,
int sel)
2789 vector<const RecCgemCluster*> aVecCgemCluster;
2796 if(sel==2)
continue;
2798 myIterClu=myIterCgemClusterBegin+cluIdx;
2799 aVecCgemCluster.push_back(*myIterClu);
2800 myVecCluIdx.push_back(-(i+1));
2802 return aVecCgemCluster;
2805vector<const RecCgemCluster*> DotsConnection::getCluVec(
struct trkCandi& aTrkCandi,
int sel)
2807 myVecCluIdx.clear();
2808 vector<const RecCgemCluster*> aVecCgemCluster = getXCluVec(aTrkCandi, sel);
2809 vector<const RecCgemCluster*> vVecCgemCluster = getVCluVec(aTrkCandi, sel);
2810 aVecCgemCluster.insert(aVecCgemCluster.end(), vVecCgemCluster.begin(), vVecCgemCluster.end());
2811 return aVecCgemCluster;
2814void DotsConnection::tagRedudantHits(
struct trkCandi& aTrkCandi)
2817 cout<<
"DotsConnection::tagRedudantHits: "<<endl;
2818 vector<vector<int> > hitSets[43];
2819 vector<vector<int> > outerHitSets[43];
2820 vector<vector<int> > innerHitSets[43];
2821 vector<int> isNode[43];
2822 vector<int> vecWireIdx;
2823 vector<vector<int> > trackPaths;
2825 int iLayerHasHits=0;
2826 for(
int layer=8; layer<43; layer++)
2828 int nHits = aTrkCandi.
mdcHitIdx[layer].size();
2831 map<int, int> idxSet;
2832 for(
int iHit=0; iHit<nHits; iHit++) {
2833 int idx = aTrkCandi.
mdcHitIdx[layer][iHit];
2837 int wire = myMdcGeomSvc->
Wire(wireIdx)->
Cell();
2839 idxSet[wire]=wireIdx;
2841 map<int, int>::iterator it=idxSet.begin();
2842 int first(-100), last(-100),
prev(-100);
2845 cout<<
" layer "<<layer<<
", "<<nHits<<
" hits: ";
2846 for(
int iHit=0; iHit<nHits; iHit++, it++) {
2847 if(iHit==0)
first=it->first;
2853 if((it->first-prev)>1) {
2856 hitSets[layer].push_back(vecWireIdx);
2861 cout<<
" "<<it->first;
2863 vecWireIdx.push_back(it->second);
2865 if(first==0&&last+1==myNWire[layer]) {
2867 if(hitSets[layer].size()>0) {
2868 vector<int> firstVec = hitSets[layer][0];
2869 hitSets[layer].erase(hitSets[layer].begin());
2870 vecWireIdx.insert(vecWireIdx.end(), firstVec.begin(), firstVec.end());
2874 hitSets[layer].push_back(vecWireIdx);
2877 cout<<
" ===> "<<nGap<<
" gaps"<<endl;
2878 int nSet=hitSets[layer].size();
2879 isNode[layer].clear();
2880 isNode[layer].resize(nSet,0);
2881 if(iLayerHasHits==1) {
2882 for(
int iSet=0; iSet<nSet; iSet++) {
2883 trackPaths.push_back(hitSets[layer][iSet]);
2886 }
else if(iLayerHasHits>1) {
2887 int nSetIn = hitSets[layer-1].size();
2888 vector<vector<int> > newTrkPaths;
2889 vector<vector<int> > innerNb;
2890 for(
int i2=0; i2<nSet; i2++) {
2892 innerNb.push_back(aVec);
2894 for(
int i1=0; i1<nSetIn; i1++) {
2895 vector<int> outerNb;
2897 cout<<
" inner-hit-set "<<i1<<
" has outer-Nb-hit-set: ";
2898 for(
int iSet=0; iSet<nSet; iSet++) {
2899 if(isInOutNeighbours(hitSets[layer-1][i1], hitSets[layer][iSet])) {
2900 outerNb.push_back(iSet);
2903 innerNb[iSet].push_back(i1);
2907 cout<<
" ("<<outerNb.size()<<
" in total)"<<endl;
2908 if(outerNb.size()>1) {
2909 isNode[layer-1][i1]=1;
2911 cout<<
" innner set "<<i1<<
" is a node"<<endl;
2913 outerHitSets[layer-1].push_back(outerNb);
2916 int aWireIdx=hitSets[layer-1][i1][0];
2917 vector<vector<int> >::iterator it_path=trackPaths.begin();
2918 for(; it_path!=trackPaths.end(); it_path++) {
2919 vector<int> aPath = *it_path;
2920 vector<int>::iterator it_hit = find(aPath.begin(), aPath.end(), aWireIdx);
2921 if(it_hit!=aPath.end()) {
2923 for(
int iSet=0; iSet<outerNb.size(); iSet++) {
2924 vector<int> aNewPath=aPath;
2925 aNewPath.insert(aNewPath.end(), hitSets[layer][outerNb[iSet]].begin(), hitSets[layer][outerNb[iSet]].end());
2926 if(iSet==0) *it_path=aNewPath;
2927 else newTrkPaths.push_back(aNewPath);
2933 cout<<
" "<<nTrk<<
" trk path(s) with this inner-hit-set"<<endl;
2936 for(
int i2=0; i2<nSet; i2++) {
2937 if(innerNb[i2].size()>1) {
2938 isNode[layer][i2]=1;
2940 cout<<
" set "<<i2<<
" is a node"<<endl;
2944 trackPaths.insert(trackPaths.end(), newTrkPaths.begin(), newTrkPaths.end());
2946 cout<<
" "<<trackPaths.size()<<
" possible track path(s) "<<endl;
2952 vector<int>* mdcDigiFitFlag;
2953 if(myWireFlag[layer]==0) mdcDigiFitFlag=&(aTrkCandi.
mdcXDigiFitted);
2955 for(
int iHit=0; iHit<nHits; iHit++) {
2956 int idx = aTrkCandi.
mdcHitIdx[layer][iHit];
2957 mdcDigiFitFlag->at(idx)=0;
2963 cout<<trackPaths.size()<<
" possible track path(s) found "<<endl;
2965 for(
int l=8; l<43; l++) {
2966 int nSet = hitSets[l].size();
2967 for(
int i=0; i<nSet; i++) {
2968 if(isNode[l][i]==1) {
2969 int nHit = hitSets[l][i].size();
2970 for(
int j=0; j<nHit; j++) {
2971 int wireIdx = hitSets[l][i][j];
2977 if(trackPaths.size()>1) {
2979 int i_bestPath = -1;
2980 int nMax_goodHits=0;
2981 vector<struct trkCandi> vecTrkCandi;
2982 for(vector<vector<int> >::iterator it_path=trackPaths.begin(); it_path!=trackPaths.end(); it_path++, i_path++) {
2984 cout<<
"circlr fit of path "<<i_path<<
": "<<endl;
2985 struct trkCandi aTrkCandi=getTrkCandi(*it_path);
2986 vecTrkCandi.push_back(aTrkCandi);
2987 int nDropHits = circleFitTaubin(vecTrkCandi.back(), 5,
false, 2,
true);
2989 if(nMax_goodHits<nFittedXhits) {
2990 nMax_goodHits=nFittedXhits;
2995 cout<<
"the best path is "<<i_bestPath<<endl;
2999bool DotsConnection::split(
struct trkCandi& aTrkCandi,
int idxCandi)
3001 const int nHitsBig=3;
3002 const int nLayerHitsBig=3;
3003 vector<vector<int> > hitSets[43];
3004 vector<int> gapSize[43];
3005 vector<int> gapSize_small[43];
3006 vector<int> bigSets[43];
3007 vector<int> inFatCluster[43];
3008 vector<vector<int> > bigInnerSets[43];
3009 vector<vector<int> > bigOuterSets[43];
3010 vector<vector<int> > outerHitSets[43];
3011 vector<vector<int> > innerHitSets[43];
3012 vector<int> isNode[43];
3013 vector<int> vecWireIdx;
3014 vector<vector<int> > trackPaths;
3015 vector<vector<int> > fatPaths;
3016 vector<vector<int> > thinPaths;
3035 int iLayerHasHits=0;
3036 for(
int layer=8; layer<43; layer++)
3038 int nHits = aTrkCandi.
mdcHitIdx[layer].size();
3047 map<int, int> idxSet;
3048 for(
int iHit=0; iHit<nHits; iHit++) {
3049 int idx = aTrkCandi.
mdcHitIdx[layer][iHit];
3053 int wire = myMdcGeomSvc->
Wire(wireIdx)->
Cell();
3055 idxSet[wire]=wireIdx;
3057 map<int, int>::iterator it=idxSet.begin();
3058 int first(-100), last(-100),
prev(-100);
3061 cout<<
" layer "<<layer<<
", "<<nHits<<
" hits: ";
3062 for(
int iHit=0; iHit<nHits; iHit++, it++) {
3063 if(iHit==0)
first=it->first;
3068 gap_size=it->first-
prev-1;
3078 hitSets[layer].push_back(vecWireIdx);
3079 gapSize[layer].push_back(gap_size);
3084 cout<<
" "<<it->first;
3086 vecWireIdx.push_back(it->second);
3088 gap_size=myNWire[layer]-last-1+
first;
3089 if(first==0&&last+1==myNWire[layer]) {
3091 if(hitSets[layer].size()>0) {
3092 vector<int> firstVec = hitSets[layer][0];
3098 hitSets[layer].erase(hitSets[layer].begin());
3099 vecWireIdx.insert(vecWireIdx.end(), firstVec.begin(), firstVec.end());
3100 gap_size=gapSize[layer][0];
3101 gapSize[layer].erase(gapSize[layer].begin());
3102 if(myDebugNb==2) cout<<
" (comb. w. 1st set) ";
3108 hitSets[layer].push_back(vecWireIdx);
3109 gapSize[layer].push_back(gap_size);
3111 int nSet=hitSets[layer].size();
3113 cout<<
" ===> "<<nGap<<
" gaps "<<endl;
3114 cout<<
" N_gapSize = "<<gapSize[layer].size()<<endl;
3117 for(
int iGap=0; iGap<nSet; iGap++)
3119 int gapSizeNext=gapSize[layer][iGap];
3120 int iPrevGap=iGap-1;
3121 if(iPrevGap<0) iPrevGap+=nSet;
3122 int gapSizePrev=gapSize[layer][iPrevGap];
3123 if(myDebugNb==2) cout<<
" gap "<<iGap<<
" size : "<<gapSizeNext;
3124 if(
abs(gapSizePrev)<
abs(gapSizeNext)) {
3126 gapSize_small[layer].push_back(-1*
abs(gapSizePrev));
3129 else gapSize_small[layer].push_back(
abs(gapSizeNext));
3130 if(myDebugNb==2) cout<<
" -> "<<gapSize_small[layer].back()<<endl;
3136 isNode[layer].clear(); isNode[layer].resize(nSet,0);
3137 inFatCluster[layer].clear(); inFatCluster[layer].resize(nSet,0);
3138 if(iLayerHasHits==1) {
3139 for(
int iSet=0; iSet<nSet; iSet++) {
3141 if(hitSets[layer][iSet].size()>nHitsBig) {
3142 bigSets[layer].push_back(iSet);
3143 vector<int> aFatPath;
3144 aFatPath.push_back(layer);
3145 aFatPath.push_back(iSet);
3146 fatPaths.push_back(aFatPath);
3150 }
else if(iLayerHasHits>1) {
3151 int nSetIn = hitSets[layer-1].size();
3153 vector<vector<int> > newFatPaths;
3154 vector<vector<int> > innerNb;
3155 vector<vector<int> > innerBigNb;
3156 for(
int iSet=0; iSet<nSet; iSet++) {
3157 if(hitSets[layer][iSet].size()>nHitsBig) bigSets[layer].push_back(iSet);
3159 innerNb.push_back(aVec);
3160 innerBigNb.push_back(aVec);
3164 for(
int i1=0; i1<nSetIn; i1++) {
3165 vector<int> outerNbSetIdx;
3166 vector<int> bigOuterNbSetIdx;
3168 cout<<
" inner-hit-set "<<i1<<
" has outer-Nb-hit-set: ";
3170 for(
int iSet=0; iSet<nSet; iSet++) {
3171 if(isInOutNeighbours(hitSets[layer-1][i1], hitSets[layer][iSet])) {
3174 outerNbSetIdx.push_back(iSet);
3175 innerNb[iSet].push_back(i1);
3176 if(hitSets[layer][iSet].size()>nHitsBig) {
3179 bigOuterNbSetIdx.push_back(iSet);
3181 if(hitSets[layer-1][i1].size()>nHitsBig) innerBigNb[iSet].push_back(i1);
3186 cout<<
" ("<<outerNbSetIdx.size()<<
" in total)"<<endl;
3187 if(outerNbSetIdx.size()>1) {
3188 isNode[layer-1][i1]=1;
3190 cout<<
" innner set "<<i1<<
" is a node"<<endl;
3192 outerHitSets[layer-1].push_back(outerNbSetIdx);
3193 bigOuterSets[layer-1].push_back(bigOuterNbSetIdx);
3215 if(hitSets[layer-1][i1].size()>nHitsBig) {
3216 int nBigPath = fatPaths.size();
3217 for(
int iBigPath=0; iBigPath<nBigPath; iBigPath++)
3219 int last_layer=fatPaths[iBigPath][0]+fatPaths[iBigPath].size()-2;
3220 if(last_layer==(layer-1)&&fatPaths[iBigPath].back()==i1) {
3221 vector<int> aBigPath=fatPaths[iBigPath];
3222 for(
int iSet=0; iSet<bigOuterNbSetIdx.size(); iSet++) {
3223 int setIdx = bigOuterNbSetIdx[iSet];
3224 vector<int> aNewBigPath = aBigPath;
3225 aNewBigPath.push_back(setIdx);
3226 if(iSet==0) fatPaths[iBigPath]=aNewBigPath;
3227 else newFatPaths.push_back(aNewBigPath);
3233 bigInnerSets[layer]=innerBigNb;
3234 innerHitSets[layer]=innerNb;
3235 for(
int iSet=0; iSet<bigSets[layer].size(); iSet++) {
3236 int setIdx = bigSets[layer][iSet];
3237 if(innerBigNb[setIdx].size()==0) {
3238 vector<int> aFatPath;
3239 aFatPath.push_back(layer);
3240 aFatPath.push_back(setIdx);
3241 newFatPaths.push_back(aFatPath);
3244 fatPaths.insert(fatPaths.end(), newFatPaths.begin(), newFatPaths.end());
3245 for(
int i2=0; i2<nSet; i2++)
3247 if(innerNb[i2].size()>1) {
3248 isNode[layer][i2]=1;
3250 cout<<
" set "<<i2<<
" is a node"<<endl;
3285 if(myDebugNb==2) cout<<
" In total, "<<fatPaths.size()<<
" fat path(s) found "<<endl;
3288 int nFatPath = fatPaths.size();
3289 for(
int iFat=0; iFat<nFatPath; iFat++) {
3290 int nLayer = fatPaths[iFat].size()-1;
3291 if(nLayer>nLayerHitsBig) {
3293 cout<<
" a fat path with nLayer="<<nLayer<<
" :";
3294 for(
int iL=0; iL<nLayer; iL++) {
3295 int l=fatPaths[iFat][0]+iL;
3296 int setIdx = fatPaths[iFat][iL+1];
3298 cout<<
" L"<<l<<
",Set"<<setIdx<<
" ";
3299 inFatCluster[l][setIdx]=1;
3308 for(
int l=8; l<43; l++) {
3309 int nBigThisLayer = 0;
3310 int nSet= inFatCluster[l].size();
3312 for(
int iSet=0; iSet<nSet; iSet++)
3313 if(inFatCluster[l][iSet]!=0) {
3317 if(myDebugNb==2) cout<<
" in big cluster:"<<endl;
3318 if(nBigThisLayer==1)
if(myDebugNb==2) cout<<
" layer "<<l<<
" set";
3322 if(nBigThisLayer>0)
if(myDebugNb==2) cout<<endl;
3330 cout<<
" === form track paths === "<<endl;
3332 for(
int l=8; l<43; l++) {
3333 vector<vector<int> > newTrkPaths;
3335 int nSet = hitSets[l].size();
3336 if(nSet>0) iLayerHasHits++;
3339 cout<<
" layer "<<l<<
", "<<nSet<<
" hit-sets";
3405 if(iLayerHasHits>1) {
3406 int nPath=trackPaths.size();
3407 int nSetIn=hitSets[l-1].size();
3408 for(
int i1=0; i1<nSetIn; i1++)
3410 if(inFatCluster[l-1][i1]!=0)
continue;
3411 int aWireIdx=hitSets[l-1][i1][0];
3412 for(
int ipath=0; ipath<nPath; ipath++) {
3413 vector<int> aPath = trackPaths[ipath];
3414 int pathSize = aPath.size();
3415 vector<int>::iterator it_hit = find(aPath.begin(), aPath.end(), aWireIdx);
3416 if(it_hit!=aPath.end()) {
3418 vector<int>& outerNbSetIdx=outerHitSets[l-1][i1];
3433 for(
int i2=0; i2<outerNbSetIdx.size(); i2++) {
3434 int setIdx=outerNbSetIdx[i2];
3435 if(inFatCluster[l][setIdx]!=0)
continue;
3437 vector<int> aNewPath=aPath;
3438 aNewPath.insert(aNewPath.end(), hitSets[l][setIdx].begin(), hitSets[l][setIdx].end());
3454 if(nBr==1) trackPaths[ipath]=aNewPath;
3455 else newTrkPaths.push_back(aNewPath);
3458 cout<<
" inner-hit-set "<<i1<<
" matches trk path "<<ipath<<
", "<<nBr<<
" branches found"<<endl;
3464 for(
int i=0; i<nSet; i++) {
3465 int nHit = hitSets[l][i].size();
3467 if(isNode[l][i]==1) {
3468 for(
int j=0; j<nHit; j++) {
3469 int wireIdx = hitSets[l][i][j];
3475 for(
int j=0; j<nHit; j++) {
3476 int wireIdx = hitSets[l][i][j];
3482 if(inFatCluster[l][i]==0) {
3483 int nNewPath = newTrkPaths.size();
3484 if(iLayerHasHits==1) newTrkPaths.push_back(hitSets[l][i]);
3486 int nInNb=innerHitSets[l][i].size();
3489 for(
int i1=0; i1<innerHitSets[l][i].size(); i1++) {
3490 if(inFatCluster[l-1][innerHitSets[l][i][i1]]==0) nInSmallNb++;
3493 newTrkPaths.push_back(hitSets[l][i]);
3504 if(myDebugNb==2&&nNewPath<newTrkPaths.size()) cout<<
" hit-set "<<i<<
" taken as a start of a brand new trk path (root)"<<endl;
3508 trackPaths.insert(trackPaths.end(), newTrkPaths.begin(), newTrkPaths.end());
3510 cout<<
" ---> "<<trackPaths.size()<<
" possible track path(s) "<<endl;
3514 cout<<
" In total, "<<trackPaths.size()<<
" possible track path(s) found "<<endl;
3518 if(trackPaths.size()>1) {
3520 int i_bestPath = -1;
3521 int nMax_goodHits=0;
3524 myVecNhitPathId.clear();
3525 myVecTrkPaths.clear();
3526 for(vector<vector<int> >::iterator it_path=trackPaths.begin(); it_path!=trackPaths.end(); it_path++, i_path++) {
3528 cout<<
" ~~~ circlr fit of path "<<i_path<<
": "<<endl;
3529 struct trkCandi oneTrkCandi=getTrkCandi(*it_path);
3531 myVecTrkPaths.push_back(oneTrkCandi);
3532 int nDropHits = circleFitTaubin(myVecTrkPaths.back(), 5,
false, 2,
true);
3534 if(nFittedXhits<4)
continue;
3536 pair<int, int> nHitPathId(nFittedXhits, myVecTrkPaths.size()-1);
3537 myVecNhitPathId.push_back(nHitPathId);
3538 if(nMax_goodHits<nFittedXhits) {
3539 nMax_goodHits=nFittedXhits;
3544 int nPath = myVecNhitPathId.size();
3549 myVecNhitPathId.erase(myVecNhitPathId.begin()+i_bestPath);
3552 int trkIdx = myVecTrkCandidates.size();
3556 aTrkCandi.
trkIdx = trkIdx;
3559 myVecTrkPaths[i_bestPath].trkIdx = idxCandi;
3562 myVecTrkCandidates.insert(myVecTrkCandidates.begin()+idxCandi, myVecTrkPaths[i_bestPath]);
3563 for(
size_t idxx = idxCandi; idxx < myVecTrkCandidates.size();idxx ++) myVecTrkCandidates[trkIdx].trkIdx = idxx;
3567 cout<<
" ----- the best path is "<<i_bestPath<<endl;
3568 cout<<
" save the original trk candidate as the last one "<<endl;
3569 cout<<
" save the best path as the trk candidate in processing"<<endl;
3578void DotsConnection::tagClusterHits(
struct trkCandi& aTrkCandi)
3581 int nContiLayTooManyHits=0;
3582 for(
int layer=8; layer<43; layer++) {
3583 int nHits = aTrkCandi.
mdcHitIdx[layer].size();
3586 nContiLayTooManyHits++;
3588 if(nHits<3||layer==42) {
3589 if(nContiLayTooManyHits>=3)
3592 cout<<
" nContiLayTooManyHits>=3 to set fit flag=0"<<endl;
3593 for(
int iL=layer; iL>layer-nContiLayTooManyHits; iL--)
3595 vector<int>& hitId = aTrkCandi.
mdcHitIdx[iL];
3596 if(myWireFlag[iL]==0) {
3597 for(
int i=0; i<hitId.size(); i++) {
3603 for(
int i=0; i<hitId.size(); i++) {
3610 nContiLayTooManyHits=0;
3613 cout<<
"tagClusterHits: layer "<<layer<<
", nHits="<<nHits<<
", nContiLayTooManyHits="<<nContiLayTooManyHits<<endl;
3625bool DotsConnection::isInOutNeighbours(
const vector<int>& innerSet,
const vector<int>& outerSet)
3628 int nInner = innerSet.size();
3629 vector<int>::const_iterator
iter;
3631 for(
int i=0; i<nInner; i++) {
3632 int wireIdx=innerSet[i];
3635 const vector<int>& vectmp=myMdcDigiGroup[wireIdx].
GetOuter();
3636 vector<int>::const_iterator it=vectmp.begin();
3637 for(; it!=vectmp.end(); it++) {
3639 iter = find(outerSet.begin(), outerSet.end(), *it);
3641 if(
iter!=outerSet.end())
3653bool DotsConnection::isInOutNeighbours(
struct trkCandi& innerTrkCandi,
struct trkCandi& outerTrkCandi,
int gap)
3657 vector<int>* innerHitIdxVec;
3658 int type = myWireFlag[innerTrkCandi.
layerMax];
3662 vector<int> vecNeighbours1;
3663 for(
int i=0; i<innerHitIdxVec->size(); i++)
3665 int wireIdx = innerHitIdxVec->at(i);
3669 vecNeighbours1.insert(vecNeighbours1.end(), vectmp.begin(), vectmp.end());
3673 vector<int>* outerHitIdxVec;
3674 type = myWireFlag[outerTrkCandi.
layerMin];
3678 vector<int> vecNeighbours2;
3679 for(
int i=0; i<outerHitIdxVec->size(); i++)
3681 int wireIdx = outerHitIdxVec->at(i);
3694 vector<int> vectmp2;
3695 int nwire=vectmp.size();
3696 for(
int j=0; j<nwire; j++) {
3698 vectmp2.insert(vectmp2.end(), vectmp3.begin(), vectmp3.end());
3703 vecNeighbours2.insert(vecNeighbours2.end(), vectmp.begin(), vectmp.end());
3707 vector<int>::iterator
iter;
3708 for(
int i=0; i<vecNeighbours1.size(); i++)
3710 iter = find(vecNeighbours2.begin(), vecNeighbours2.end(), vecNeighbours1[i]);
3711 if(
iter!=vecNeighbours2.end()) {
3720void DotsConnection::combineTrkCandi(
struct trkCandi& trkCandi_1,
struct trkCandi& trkCandi_2)
3740 layer = myMdcGeomSvc->
Wire(wireIdx)->
Layer();
3741 trkCandi_1.
mdcHitIdx[layer].push_back(nXhits_1);
3761 layer = myMdcGeomSvc->
Wire(wireIdx)->
Layer();
3778 layer = myMdcGeomSvc->
Wire(wireIdx)->
Layer();
3795 layer = myMdcGeomSvc->
Wire(wireIdx)->
Layer();
3796 trkCandi_1.
mdcHitIdx[layer].push_back(nVhits_1);
3804void DotsConnection::rmPathInTrkCandi(
struct trkCandi& aTrkCandi,
struct trkCandi& path)
3806 for(
int i=0; i<43; i++) aTrkCandi.
mdcHitIdx[i].clear();
3814 cout<<
"rmPathInTrkCandi: "<<vecMdcXDigiIdx.size()<<
" X-hits to ";
3816 vector<int> vecMdcXDigiIdx_2_fitted;
3817 int n2=vecMdcXDigiIdx_2.size();
3819 for(
int i=0; i<
n2; i++) {
3820 wireIdx=vecMdcXDigiIdx_2[i];
3821 if(path.
mdcXDigiFitted[i]>0||myMdcDigiGroup[wireIdx].isOnNode()==0) vecMdcXDigiIdx_2_fitted.push_back(wireIdx);
3823 vector<int> vecMdcXDigiIdx_new = getUncommonVec(vecMdcXDigiIdx, vecMdcXDigiIdx_2_fitted);
3825 int nXMdcDigi = vecMdcXDigiIdx_new.size();
3826 for(
int i=0; i<nXMdcDigi; i++) {
3827 wireIdx = vecMdcXDigiIdx_new[i];
3828 layer = myMdcGeomSvc->
Wire(wireIdx)->
Layer();
3829 aTrkCandi.
mdcHitIdx[layer].push_back(i);
3838 cout<<nXMdcDigi<<endl;
3839 double nXRm=vecMdcXDigiIdx.size()-nXMdcDigi;
3847 vector<int> vecMdcVDigiIdx_2_fitted;
3848 n2=vecMdcVDigiIdx_2.size();
3849 for(
int i=0; i<
n2; i++) {
3850 wireIdx=vecMdcVDigiIdx_2[i];
3851 if(path.
mdcVDigiFitted[i]>0||myMdcDigiGroup[wireIdx].isOnNode()==0) vecMdcVDigiIdx_2_fitted.push_back(vecMdcVDigiIdx_2[i]);
3853 vector<int> vecMdcVDigiIdx_new = getUncommonVec(vecMdcVDigiIdx, vecMdcVDigiIdx_2_fitted);
3854 int nVMdcDigi = vecMdcVDigiIdx_new.size();
3856 cout<<
" "<<vecMdcVDigiIdx.size()<<
" V-hits to "<<nVMdcDigi<<endl;
3857 double nVRm=vecMdcVDigiIdx.size()-nVMdcDigi;
3858 for(
int i=0; i<nVMdcDigi; i++) {
3859 wireIdx = vecMdcVDigiIdx_new[i];
3860 layer = myMdcGeomSvc->
Wire(wireIdx)->
Layer();
3861 aTrkCandi.
mdcHitIdx[layer].push_back(i);
3871 if((nXMdcDigi+nVMdcDigi)<3 || nXRm==0)
3884vector<int> DotsConnection::getUncommonVec(vector<int> v1, vector<int>& v2)
3887 vector<int>::iterator it;
3888 for(
int i=0; i<
n; i++) {
3889 it=find(v1.begin(), v1.end(), v2[i]);
3890 if(it!=v1.end()) v1.erase(it);
3895void DotsConnection::setMinusVec(set<int>&
s, vector<int>&
v)
3898 for(
int i=0; i<
n; i++) {
3899 set<int>::iterator it_s =
s.find(
v.at(i));
3900 if(it_s!=
s.end())
s.erase(it_s);
3904void DotsConnection::resetXDigiFitItem(
struct trkCandi& aTrkCandi)
3913void DotsConnection::resetVDigiFitItem(
struct trkCandi& aTrkCandi)
3922void DotsConnection::updateDigiFitItem(
struct trkCandi& aTrkCandi)
3928 for(
int i=0; i<vecActive.size(); i++)
3930 int idx = myVecDigiIdx[i];
3943void DotsConnection::updateCgemFitItem(
struct trkCandi& aTrkCandi,
int assignCgemClu)
3948 for(
int i=0; i<vecActive.size(); i++)
3950 int idx = myVecCluIdx[i];
3955 if(vecActive[i]==-1) {
3961 else if(assignCgemClu>0&&cluIdx>=0){
3962 myIterClu=myIterCgemClusterBegin+cluIdx;
3963 int layer= (*myIterClu)->getlayerid();
3964 int size=myVecCgemXClusterIdx[layer].size();
3965 for(
int j=0; j<size; j++) {
3966 if(myVecCgemXClusterIdx[layer][j]==cluIdx) {
3967 myVecCgemXClusterTrkIdx[layer][j]=aTrkCandi.
trkIdx;
3978 if(vecActive[i]==-1) {
3984 else if(assignCgemClu>0&&cluIdx>=0){
3985 myIterClu=myIterCgemClusterBegin+cluIdx;
3986 int layer= (*myIterClu)->getlayerid();
3987 int size=myVecCgemVClusterIdx[layer].size();
3988 for(
int j=0; j<size; j++) {
3989 if(myVecCgemVClusterIdx[layer][j]==cluIdx) {
3990 myVecCgemVClusterTrkIdx[layer][j]=aTrkCandi.
trkIdx;
4000void DotsConnection::dropRedundantMdcXDigi(
struct trkCandi& aTrkCandi)
4002 vector<const MdcDigi*> aGoodXDigiVec;
4003 for(
int i=8; i<43; i++)
4008 aGoodXDigiVec.push_back(myMapMdcDigi[wireIdx]);
4011 vector<int> idxVec = aTrkCandi.
mdcHitIdx[i];
4017void DotsConnection::calcuCiclePar3CgemClu(
struct trkCandi& aTrkCandi)
4020 double x[3]={0,0,0};
4021 double y[3]={0,0,0};
4025 myIterClu=myIterCgemClusterBegin+cluIdx;
4026 double phi_cluster = (*myIterClu)->getrecphi();
4027 int layer= (*myIterClu)->getlayerid();
4029 x[i]=r*
cos(phi_cluster);
4030 y[i]=r*
sin(phi_cluster);
4033 x0 = ((pow(x[2],2)+pow(y[2],2))*(y[1]-y[0])+(pow(x[1],2)+pow(y[1],2))*(y[0]-y[2])+(pow(x[0],2)+pow(y[0],2))*(y[2]-y[1]))/(2.*((x[2] - x[0])*(y[1] - y[0]) - (x[1] -x[0])*(y[2] - y[0])));
4034 y0 = ((pow(x[2],2)+pow(y[2],2))*(
x[1]-
x[0])+(pow(x[1],2)+pow(y[1],2))*(
x[0]-
x[2])+(pow(x[0],2)+pow(y[0],2))*(
x[2]-
x[1]))/(2.*((y[2] - y[0])*(x[1] - x[0]) - (y[1] -y[0])*(x[2] - x[0])));
4035 r0 = sqrt(pow(x[1]-x0, 2) + pow(y[1]-y0,2));
4036 double xmean=(
x[0]+
x[1]+
x[2])/3.0;
4037 double ymean=(y[0]+y[1]+y[2])/3.0;
4038 double dr_fit = sqrt(x0*x0+y0*y0)-r0;
4039 double phi0_fit = atan2(y0,x0);
4041 double myAlpha=aHelix.
alpha();
4042 double kappa_fit= myAlpha/r0;
4043 Hep3Vector pos_mean(xmean, ymean);
4044 Hep3Vector pos_center(x0, y0);
4045 Hep3Vector
vec_z=pos_mean.cross(pos_center);
4051 kappa_fit = -kappa_fit;
4053 while(phi0_fit< 0 ) phi0_fit+=2*
M_PI;
4054 while(phi0_fit> 2*
M_PI) phi0_fit-=2*
M_PI;
4056 aTrkCandi.
a_helix[1]=phi0_fit;
4057 aTrkCandi.
a_helix[2]=kappa_fit;
4059 cout<<__FUNCTION__<<
": circle par = "<<dr_fit<<
", "<<phi0_fit<<
", "<<kappa_fit<<endl;
4069std::vector<std::pair<double, double> > DotsConnection::findWireCrossPoint(
struct trkCandi& aTrkCandi, vector<double>* ref_track_param){
4070 std::vector<const MdcDigi *> usedDigis=getVDigiVec(aTrkCandi,1);
4097 const int prevLayers[]={23,27,31};
4098 const int nextLayers[]={24,28,32};
4099 const int allNotedLayers[]={23,24,27,28,31,32};
4101 std::vector<const MdcGeoWire *> crossWires[6];
4102 std::cout<<
"---findWireCrossPoint: getting digi layer and wires"<<
'\n';
4103 for (
size_t ii=0; ii<usedDigis.size(); ++ii){
4104 const MdcDigi * digi=usedDigis[ii];
4110 for (
size_t i=0;i<
sizeof(allNotedLayers)/
sizeof(int) ;++i){
4111 if (layer==allNotedLayers[i]) {
4116 if (layerIdx == -1)
continue;
4118 crossWires[layerIdx].push_back(wir);
4121 for (
size_t ii=0; ii<
sizeof(allNotedLayers)/
sizeof(int); ii++) {
4122 int layerGroupId=ii/2;
4123 if (crossWires[ii].empty()){
4124 crossWires[(ii xor (
size_t)1)].clear();
4127 std::vector<std::pair<double, double> > res;
4128 for (
size_t ii=0; ii<
sizeof(allNotedLayers)/
sizeof(int)/2; ii++){
4129 std::vector<const MdcGeoWire *> & wirePrev=crossWires[ii*2];
4130 std::vector<const MdcGeoWire *> & wirePost=crossWires[ii*2+1];
4131 typedef std::vector<const MdcGeoWire *>::iterator WireIt;
4132 std::vector<std::pair<double, double> > conjs;
4133 for (WireIt it1=wirePrev.begin(); it1!=wirePrev.end(); ++it1){
4134 for(WireIt it2=wirePost.begin(); it2!=wirePost.end(); ++it2){
4145 cout<<
"--- findWireCrossPoint found point: ";
4147 cout<<
"("<<p.x()<<
','<<p.y()<<
','<<p.z()<<
": layer.wire "<<(*it1)->Layer()<<
"."<<(*it1)->Id() <<
" : "<<(*it2)->Layer()<<
"."<<(*it2)->Id()<<
") ";
4148 conjs.push_back(std::make_pair(p.x(),p.y()));
4151 if (conjs.empty())
continue;
4153 std::pair<double, double> res2=std::make_pair(0.,0.);
4154 typedef std::vector<std::pair<double, double> >::iterator
PairIt;
4156 for (PairIt it1=conjs.begin(); it1!=conjs.end(); it1++){
4157 if (!isfinite(it1->first) || !isfinite(it1->second))
continue;
4158 if (!myWireCrossPointAvgMode)res.push_back(*it1);
4159 res2.first+= it1->first;
4160 res2.second+= it1->second;
4163 if (nCount==0)
continue;
4164 res2.first /=nCount;
4165 res2.second/=nCount;
4166 if (myWireCrossPointAvgMode)res.push_back(res2);
4187int DotsConnection::circleFitTaubin(
struct trkCandi& aTrkCandi,
double chiCut,
bool useIniHelix,
int hitSel,
bool setHitFitFlag,
int layerMin,
int layerMax,
bool use_additional_points,
bool usePhiAmbi,
const vector<double>& vPhiAmb)
4189 myVecDigiIdx.clear();
4190 vector<int> aXDigiFitFlagVec;
4191 vector<const MdcDigi*> aXDigiVec = getXDigiVec(aTrkCandi, hitSel, &aXDigiFitFlagVec);
4203 myDotsHelixFitter.
setDChits(aXDigiVec, myEvtT0);
4205 if(setHitFitFlag&&myDotsHelixFitter.
setDChitsFitFlag(aXDigiFitFlagVec))
if(myDebugNb==2) cout<<
"set fit flag successfully "<<endl;
4208 if(myDotsHelixFitter.
setDChitsPhiAmbi(vPhiAmb))
if(myDebugNb==2) cout<<
"usePhiAmbi"<<endl;
4211 myVecCluIdx.clear();
4212 vector<const RecCgemCluster*> aVecCgemCluster = getXCluVec(aTrkCandi,hitSel);
4217 HepVector a_helix(5,0);
4218 a_helix(1)=aTrkCandi.
a_helix[0];
4219 a_helix(2)=aTrkCandi.
a_helix[1];
4220 a_helix(3)=aTrkCandi.
a_helix[2];
4225 vector<double> circleParVec;
4226 if(use_additional_points) {
4231 std::vector<std::pair<double, double> > extraPoints=findWireCrossPoint(aTrkCandi);
4234 cout<<
" --- extraPoint for Taubin fit: ";
4235 for (
int i=0; i < extraPoints.size(); i++){
4236 cout<<
"("<<extraPoints[i].first<<
","<<extraPoints[i].second<<
") ";
4243 use_additional_points=
false;
4245 cout<<
" --- Circle par from Taubin fit: "<<circleParVec[0]
4246 <<
", "<<circleParVec[1]
4247 <<
", "<<circleParVec[2]
4248 <<
", chi2 = "<<myDotsHelixFitter.
getChi2()
4251 if(fabs(chiCut-9999.0)<0.001)
break;
4252 if(myDotsHelixFitter.
deactiveHits(chiCut, 1)==0)
break;
4256 cout<<
"DotsConnection::"<<__FUNCTION__<<
": "<<nDeactived<<
" hits deactived "<<endl;
4260 int n_reactive = myDotsHelixFitter.
activeHits(20, 1);
4261 if(n_reactive==0)
break;
4263 n_active+=n_reactive;
4266 cout<<
" +++ Circle par from Taubin fit: "<<circleParVec[0]
4267 <<
", "<<circleParVec[1]
4268 <<
", "<<circleParVec[2]
4269 <<
", chi2 = "<<myDotsHelixFitter.
getChi2()
4275 cout<<
"DotsConnection::"<<__FUNCTION__<<
": "<<n_active<<
" uncertain hits actived "<<endl;
4279 int n_reactive = myDotsHelixFitter.
activeHits(chiCut, 0);
4280 if(n_reactive==0)
break;
4281 n_active_2+=n_reactive;
4285 cout<<
" +++ Circle par from Taubin fit: "<<circleParVec[0]
4286 <<
", "<<circleParVec[1]
4287 <<
", "<<circleParVec[2]
4288 <<
", chi2 = "<<myDotsHelixFitter.
getChi2()
4294 cout<<
"DotsConnection::"<<__FUNCTION__<<
": "<<n_active_2<<
" hits reactived "<<endl;
4298 if(myDotsHelixFitter.
deactiveHits(chiCut, 1)==0)
break;
4299 if(fabs(chiCut-9999.0)<0.001)
break;
4300 else nDeactived_2++;
4303 cout<<
" --- Circle par from Taubin fit: "<<circleParVec[0]
4304 <<
", "<<circleParVec[1]
4305 <<
", "<<circleParVec[2]
4306 <<
", chi2 = "<<myDotsHelixFitter.
getChi2()
4311 cout<<
"DotsConnection::"<<__FUNCTION__<<
": "<<nDeactived_2<<
" hits deactived "<<endl;
4312 nDeactived=nDeactived-n_active-n_active_2+nDeactived_2;
4314 aTrkCandi.
a_helix[0]=circleParVec[0];
4315 aTrkCandi.
a_helix[1]=circleParVec[1];
4316 aTrkCandi.
a_helix[2]=circleParVec[2];
4322void DotsConnection::circleFit(
struct trkCandi& aTrkCandi)
4324 myVecDigiIdx.clear();
4325 vector<const MdcDigi*> aXDigiVec = getXDigiVec(aTrkCandi);
4328 myDotsHelixFitter.
setDChits(aXDigiVec, myEvtT0);
4331 HepVector a_helix(5,0);
4332 a_helix(1)=aTrkCandi.
a_helix[0];
4333 a_helix(2)=aTrkCandi.
a_helix[1];
4334 a_helix(3)=aTrkCandi.
a_helix[2];
4340 HepVector new_a = myDotsHelixFitter.
getHelix();
4342 cout<<
"Circle par from LS fit: "<<new_a[0]
4345 <<
", chi2 = "<<myDotsHelixFitter.
getChi2()
4353 vector<const MdcDigi*> aDigiVec = getDigiVec(aTrkCandi, 3);
4354 myDotsHelixFitter.
setDChits(aDigiVec, myEvtT0);
4355 vector<const RecCgemCluster*> aVecCgemCluster = getCluVec(aTrkCandi, 2);
4359 HepVector a_helix(5,0);
4360 a_helix(1)=aTrkCandi.
a_helix[0];
4361 a_helix(2)=aTrkCandi.
a_helix[1];
4362 a_helix(3)=aTrkCandi.
a_helix[2];
4363 a_helix(4)=aTrkCandi.
a_helix[3];
4364 a_helix(5)=aTrkCandi.
a_helix[4];
4376 if(fitFlag!=0)
break;
4387 if(fitFlag!=0)
break;
4389 int nReActi = myDotsHelixFitter.
activeHits(chiCut);
4390 if(nReActi==0)
break;
4391 nReActiTot+=nReActi;
4396 HepVector aHelix = myDotsHelixFitter.
getHelix();
4397 double chi2 = myDotsHelixFitter.
getChi2();
4399 for(
int i=0; i<5; i++) aTrkCandi.
a_helix[i]=aHelix[i];
4400 aTrkCandi.
chi2=chi2;
4407 cout<<
"fitFlag="<<fitFlag<<endl;
4408 cout<<nDropTot<<
" hits dropped, "<<nReActiTot<<
" hits reactived "<<endl;
4409 cout<<
"helix chi2="<<myDotsHelixFitter.
getChi2()<<endl;
4410 cout<<
"fitted helix: ("<<aHelix(1)<<
", "<<aHelix(2)<<
", "<<aHelix(3)<<
", "<<aHelix(4)<<
", "<<aHelix(5)<<
")"<<endl;
4417bool DotsConnection::matchMdcDigi(
struct trkCandi& aTrkCandi,
int idx)
4420 int nMaxXHits = nXHitsActive(aTrkCandi);
4421 for(; direct<=1; direct+=2)
4423 if(aTrkCandi.
isGood!=1)
continue;
4425 if(direct<0) layer = aTrkCandi.
layerMin;
4427 int layerToFind = layer+direct;
4429 cout<<
" try to match "<<direct<<
" MDC hits: "<<endl;
4430 while(layerToFind>=8&&layerToFind<=42)
4432 int gap =
abs(layer-layerToFind)-1;
4436 cout<<
" layerToFind = "<<layerToFind<<
", gap="<<gap<<endl;
4437 vector<struct trkCandi>::iterator iter_other = myVecTrkCandidates.begin();
4442 vector<struct trkCandi>::iterator iter_matchBest;
4443 for(; iter_other!=myVecTrkCandidates.end(); iter_other++, idxTrk++)
4446 if(idxTrk==idx)
continue;
4448 if(direct<0) layerFound = (*iter_other).layerMax;
4449 else layerFound = (*iter_other).layerMin;
4450 if(layerFound!=layerToFind)
continue;
4451 if((*iter_other).isGood==-1)
continue;
4453 cout<<
" layer "<<layerToFind<<
" is found"<<endl;
4455 struct trkCandi& innerTrkCandi = *iter_other;
4456 struct trkCandi& outerTrkCandi = aTrkCandi;
4457 if( isInOutNeighbours(innerTrkCandi,outerTrkCandi,gap)==
false )
continue;
4460 struct trkCandi& outerTrkCandi = *iter_other;
4461 struct trkCandi& innerTrkCandi = aTrkCandi;
4462 if( isInOutNeighbours(innerTrkCandi,outerTrkCandi,gap)==
false )
continue;
4465 cout<<
" "<<gap<<
"-layer-gap neighboured ";
4466 if(direct<0) cout<<
"inner ";
4467 else cout<<
"outer ";
4468 cout<<
"segment found: "<<endl;}
4471 vector<int>::iterator iter_digi = (*iter_other).mdcXDigiWireIdx.begin();
4472 int n_matchedXhits=0;
4474 for(; iter_digi!=(*iter_other).mdcXDigiWireIdx.end(); iter_digi++)
4476 int wireIdx = *iter_digi;
4477 int layer = myMdcGeomSvc->
Wire(wireIdx)->
Layer();
4478 if(myWireFlag[layer]!=0)
continue;
4479 int wire = myMdcGeomSvc->
Wire(wireIdx)->
Cell();
4482 cout<<
" (layer "<<layer<<
", wire "<<wire<<
", view "<<myWireFlag[layer]<<
", chi "<<myDotsHelixFitter.
getDigiChi()<<
") "<<endl;
4484 if(fabs(myDotsHelixFitter.
getDigiChi())<500) n_matchedXhits++;
4486 if(n_matchedXhits>nmax_match) {
4487 nmax_match=n_matchedXhits;
4488 iter_matchBest = iter_other;
4494 cout<<
" "<<n_neighbSeg<<
" neighboured trk segments found, and the best one is to be combined"<<endl;
4496 combineTrkCandi(aTrkCandi, *iter_matchBest);
4499 cout<<
"nmax_match="<<nmax_match<<
", n_Xhits="<<n_Xhits<<endl;
4501 bool useIniHelix=
true;
4502 if(n_Xhits>nMaxXHits) {
4505 cout<<
" circle fit without ini Helix, with all hits"<<endl;
4507 else if(myDebugNb==2) cout<<
" circle fit with ini Helix and all hits"<<endl;
4508 circleFitTaubin(aTrkCandi, 5, useIniHelix, 1,
false);
4509 updateDigiFitItem(aTrkCandi);
4511 else if(n_Xhits>0) {
4514 cout<<
" circle fit without ini Helix, rejected hits"<<endl;
4515 circleFitTaubin(aTrkCandi, 5,
false, 2,
true);
4518 if(nXHitAct<nMaxXHits) {
4520 cout<<
" "<<nXHitAct<<
" hits kept, < "<<nMaxXHits<<
" "<<endl;
4521 cout<<
" circle fit with inner layer hits: "<<endl;}
4522 setFitFlagUncertain(aTrkCandi, 20, 60);
4523 circleFitTaubin(aTrkCandi, 5,
false, 2,
true);
4526 cout<<
" "<<nXHitAct<<
" hits kept"<<endl;
4528 updateDigiFitItem(aTrkCandi);
4530 double dr=aTrkCandi.
a_helix[0];
4534 cout<<
"a bad circleFitTaubin! "<<endl;
4536 resetXDigiFitItem(aTrkCandi);
4540 if(direct<0) layer = aTrkCandi.
layerMin;
4542 layerToFind = layer+direct;
4546 else layerToFind +=direct;
4554bool DotsConnection::matchCgemXCluster(
struct trkCandi& aTrkCandi)
4558 cout<<
"as layerMin = "<<aTrkCandi.
layerMin<<
" > 10, no CGEM cluster matching tried "<<endl;
4561 int nCluMatchedLastTime=0;
4565 cout<<
" ~~ the "<<nIter+1<<
" time to match Cgem X-Cluster: "<<endl;
4568 double minResidX[3]={9999., 9999., 9999.};
4569 int idxClusterMinResid[3]={-1, -1, -1};
4571 int nCluExpt(0), nLayerWithCluNearby(0);
4573 for(
int i=2; i>=0; i--)
4577 phi_trk = pos.phi();
4578 int nClu = myVecCgemXClusterIdx[i].size();
4579 for(
int j=0; j<nClu; j++)
4581 myIterClu=myIterCgemClusterBegin+myVecCgemXClusterIdx[i][j];
4582 double phi_CC = (*myIterClu)->getrecphi();
4583 double del_phi = phi_CC-phi_trk;
4584 if(del_phi<-M_PI||del_phi>
M_PI) del_phi=atan2(
sin(del_phi),
cos(del_phi));
4586 if(fabs(del_X)<10.0)
4589 cout<<
"CGEM layer "<<i<<
", a X-cluster with dX="<<del_X*10<<
" mm found"<<endl;
4590 if(fabs(del_X)<fabs(minResidX[i])) {
4592 idxClusterMinResid[i]=myVecCgemXClusterIdx[i][j];
4596 if(fabs(minResidX[i])<1.0)
4601 if(fabs(minResidX[i])<10.) nLayerWithCluNearby++;
4605 int nDropHits = circleFitTaubin(aTrkCandi, 5,
true, 2,
true);
4607 cout<<
"After matching CGEM and circle fit, nDropHits="<<nDropHits<<endl;
4608 if(nDropHits>nCluMatched) {
4609 double circlePar[3];
4610 circlePar[0]=aTrkCandi.
a_helix[0];
4611 circlePar[1]=aTrkCandi.
a_helix[1];
4612 circlePar[2]=aTrkCandi.
a_helix[2];
4614 cout<<
"nDropHits>nAdd => refit without initial helix: "<<endl;
4615 int nDropHits_2 = circleFitTaubin(aTrkCandi, 5,
false, 2,
true);
4617 cout<<
"Without initial helix, nDropHits="<<nDropHits_2<<endl;
4618 if(nDropHits_2<nDropHits) nDropHits=nDropHits_2;
4620 aTrkCandi.
a_helix[0]=circlePar[0];
4621 aTrkCandi.
a_helix[1]=circlePar[1];
4622 aTrkCandi.
a_helix[2]=circlePar[2];
4623 nDropHits = circleFitTaubin(aTrkCandi, 5,
true, 2,
true);
4626 updateCgemFitItem(aTrkCandi);
4627 int nXCluDrop = NXcluToDrop(aTrkCandi);
4628 while(nDropHits>nCluMatched||nXCluDrop>0) {
4630 cout<<
" nDropHits="<<nDropHits<<
", nCluMatched="<<nCluMatched<<
", nXCluDrop="<<nXCluDrop<<
"=> drop one CGEM X-cluster"<<endl;
4632 nDropHits = circleFitTaubin(aTrkCandi, 5,
true, 2,
true);
4633 updateCgemFitItem(aTrkCandi);
4634 nXCluDrop = NXcluToDrop(aTrkCandi);
4636 if(nCluMatched<=0)
break;
4639 if(nLayerWithCluNearby==nCluMatched)
break;
4641 if(nCluMatched<=nCluMatchedLastTime)
break;
4643 nCluMatchedLastTime = nCluMatched;
4665void DotsConnection::matchCgemClusterMdcDigi()
4668 vector<struct trkCandi>::iterator iter_trkCandi = myVecTrkCandidates.begin();
4670 for(; iter_trkCandi!=myVecTrkCandidates.end(); iter_trkCandi++, i_candi++)
4673 cout<<endl<<
" ----> trk candi "<<i_candi<<
": "<<endl;
4674 if((*iter_trkCandi).isGood!=1)
continue;
4675 if((*iter_trkCandi).mdcXDigiWireIdx.size()>=3)
4678 tagRedudantHits(*iter_trkCandi);
4679 int nDropHits = circleFitTaubin(*iter_trkCandi, 5,
false, 2,
true);
4680 updateDigiFitItem(*iter_trkCandi);
4687 double circlePar[3];
4688 circlePar[0]=(*iter_trkCandi).a_helix[0];
4689 circlePar[1]=(*iter_trkCandi).a_helix[1];
4690 circlePar[2]=(*iter_trkCandi).a_helix[2];
4691 if(myDotsHelixFitter.
getNActiveHits()<4||circlePar[0]!=circlePar[0])
4694 cout<<
"a bad circleFitTaubin! "<<endl;
4695 (*iter_trkCandi).isGood=0;
4696 resetXDigiFitItem(*iter_trkCandi);
4829 matchMdcDigi(*iter_trkCandi, i_candi);
4831 if((*iter_trkCandi).isGood==0)
continue;
4836 matchCgemXCluster(*iter_trkCandi);
4922 cout<<
"CGEM V-cluster candidates (s,z): "<<endl;
4924 vector<int>::iterator iter_XClu = iter_trkCandi->CgemXClusterID.begin();
4926 myRoughTanlDzMap.Reset();
4927 map<int,pair<double,double> > map_Vclu_sz;
4928 for(; iter_XClu!=iter_trkCandi->CgemXClusterID.end(); iter_XClu++)
4930 int idx_Xclu = *iter_XClu;
4931 if(idx_Xclu<0)
continue;
4932 myIterClu=myIterCgemClusterBegin+idx_Xclu;
4933 int layer = (*myIterClu)->getlayerid();
4935 cout<<
"layer "<<layer;
4938 if(fabs(dphi)<1e-10) {
4940 cout<<
"DotsConnection::"<<__FUNCTION__<<
": error, no intersection between track and CGEM layer "<<layer<<endl;
4943 double s = fabs(dphi*circle_fittedHelix.
radius());
4945 int nV = myVecCgemVClusterIdx[layer].size();
4946 for(
int iV=0; iV<nV; iV++) {
4947 int idx_Vclu = myVecCgemVClusterIdx[layer][iV];
4948 if(getCgemClusterIntersection(idx_Xclu,idx_Vclu,z_XV)) {
4949 pair<double, double> sz(
s,z_XV);
4950 fillTanlDzMap(
s,z_XV);
4952 map_Vclu_sz[idx_Vclu]=sz;
4954 cout<<
" ("<<
s<<
", "<<z_XV<<
")";
4966 vector<vector<pair<double, double> > > vec_vecsz;
4967 vector<int>& mdcVDigiWireIdx = iter_trkCandi->mdcVDigiWireIdx;
4968 int nVdigi = mdcVDigiWireIdx.size();
4969 vector<int>::iterator iter_digi = mdcVDigiWireIdx.begin();
4971 cout<<
"MDC V-hits candidates (s,z): "<<endl;
4973 for( ; i_digi<nVdigi; i_digi++,iter_digi++)
4975 int wireIdx = *iter_digi;
4976 const MdcDigi* aMdcDigi = myMapMdcDigi[wireIdx];
4977 vector<pair<double, double> > vec_sz = myDotsHelixFitter.
getSZ(aMdcDigi);
4978 vec_vecsz.push_back(vec_sz);
4979 vector<pair<double, double> >::iterator iter_sz = vec_sz.begin();
4980 for(; iter_sz!=vec_sz.end(); iter_sz++) {
4981 fillTanlDzMap(iter_sz->first,iter_sz->second);
4983 if(vec_sz.size()>0) nsz_mdc++;
4985 if((i_digi+1)%3==0)
if(myDebugNb==2) cout<<endl;
4987 if((i_digi+1)%3!=1)
if(myDebugNb==2) cout<<endl;
4990 if((nsz_cgem+nsz_mdc)<3) {
4992 cout<<
" no sufficient sz pair! stop tanl, dz, V-hits association. "<<endl;
4993 (*iter_trkCandi).isGood=0;
4999 double x_peak(0.), y_peak(0.);
5000 double x_weight(0.), y_weight(0.);
5001 getWeightedPeak(myRoughTanlDzMap, x_peak, y_peak, x_weight, y_weight);
5004 cout<<
"from tanL-dz Hough map: dz="<<y_weight<<
", tanL = "<<x_weight<<endl;
5007 myRoughTanlDzMap.Reset();
5008 double s_sum(0), z_sum(0);
5009 double s2_sum(0), sz_sum(0);
5013 int cgemVCluIdx[3]={-1,-1,-1};
5014 double cgemVClu_dz_min[3]={9999.,9999.,9999.};
5015 double cgemVClu_s_dzmin[3]={9999.,9999.,9999.};
5016 double cgemVClu_z_dzmin[3]={9999.,9999.,9999.};
5018 cout<<
"CGEM V-cluster candidates: "<<endl;
5019 for(map<
int,pair<double,double> >::iterator i_map=map_Vclu_sz.begin(); i_map!=map_Vclu_sz.end(); i_map++)
5021 double s = (*i_map).second.first;
5022 double z = (*i_map).second.second;
5023 double z_hough = x_weight*
s+y_weight;
5024 double delta_z = z-z_hough;
5026 cout<<
" delta_z = "<<delta_z<<
" ";
5027 double deltaZcut=10.;
5028 if(fabs(delta_z)<deltaZcut) {
5029 int idx_Vclu = (*i_map).first;
5030 myIterClu=myIterCgemClusterBegin+idx_Vclu;
5031 int layer = (*myIterClu)->getlayerid();
5032 if(fabs(delta_z)<fabs(cgemVClu_dz_min[layer])) {
5033 cgemVClu_dz_min[layer]=delta_z;
5034 cgemVClu_s_dzmin[layer]=
s;
5035 cgemVClu_z_dzmin[layer]=z;
5036 cgemVCluIdx[layer]=idx_Vclu;
5042 for(
int layer=2; layer>=0; layer--) {
5043 if(cgemVCluIdx[layer]!=-1) {
5045 fillTanlDzMap(cgemVClu_s_dzmin[layer],cgemVClu_z_dzmin[layer]);
5046 s_sum+=cgemVClu_s_dzmin[layer];
5047 s2_sum+=cgemVClu_s_dzmin[layer]*cgemVClu_s_dzmin[layer];
5048 z_sum+=cgemVClu_z_dzmin[layer];
5049 sz_sum+=cgemVClu_s_dzmin[layer]*cgemVClu_z_dzmin[layer];
5056 cout<<
"MDC V-hits candidates: "<<endl;
5057 resetVDigiFitItem(*iter_trkCandi);
5058 vector<vector<pair<double, double> > >::iterator iter_vvsz = vec_vecsz.begin();
5060 for(; iter_vvsz!=vec_vecsz.end(); iter_vvsz++, i_digi++)
5062 vector<pair<double, double> >& vec_sz = *iter_vvsz;
5063 vector<pair<double, double> >::iterator iter_sz = vec_sz.begin();
5064 double delta_z_min = 9999.0;
5065 double z_dzmin(0), s_dzmin(0);
5066 for(; iter_sz!=vec_sz.end(); iter_sz++) {
5067 double s = iter_sz->first;
5068 double z = iter_sz->second;
5069 double z_hough = x_weight*
s+y_weight;
5070 double delta_z = z-z_hough;
5071 if(fabs(delta_z)<fabs(delta_z_min)) {
5072 delta_z_min=delta_z;
5078 cout<<
"delta_z_min="<<delta_z_min;
5079 double deltaZcut_mdc = 10.0;
5080 if(fabs(delta_z_min)>deltaZcut_mdc) {
5083 iter_trkCandi->mdcVDigiFitted[i_digi]=-1;
5086 iter_trkCandi->mdcVDigiFitted[i_digi]=1;
5087 fillTanlDzMap(s_dzmin, z_dzmin);
5089 cout<<
"("<<s_dzmin<<
","<<z_dzmin<<
")";
5092 s2_sum+=s_dzmin*s_dzmin;
5093 sz_sum+=s_dzmin*z_dzmin;
5098 if((i_digi+1)%3==0) cout<<endl;}
5101 if((i_digi+1)%3!=1) cout<<endl;
5104 getWeightedPeak(myRoughTanlDzMap, x_peak, y_peak, x_weight, y_weight);
5106 cout<<
"After rough selection of sz of V-hits: "<<endl;
5108 cout<<
"from tanL-dz Hough map: dz="<<y_weight<<
", tanL = "<<x_weight<<endl;}
5111 double tanl_lfit = (n_sz*sz_sum-s_sum*z_sum)/(n_sz*s2_sum-s_sum*s_sum);
5112 double dz_lfit = (s2_sum*z_sum-s_sum*sz_sum)/(n_sz*s2_sum-s_sum*s_sum);
5114 cout<<
"from s-z line fit: dz="<<dz_lfit<<
", tanL = "<<tanl_lfit<<endl;
5116 iter_trkCandi->a_helix[3]=dz_lfit;
5117 iter_trkCandi->a_helix[4]=tanl_lfit;
5121 for(
int i=0; i<3; i++) {
5123 cgemVClu_dz_min[i]=9999.;
5124 cgemVClu_s_dzmin[i]=9999.;
5125 cgemVClu_z_dzmin[i]=9999.;
5128 cout<<
"CGEM V-cluster candidates: "<<endl;
5129 for(map<
int,pair<double,double> >::iterator i_map=map_Vclu_sz.begin(); i_map!=map_Vclu_sz.end(); i_map++)
5131 double s = (*i_map).second.first;
5132 double z = (*i_map).second.second;
5133 double z_hough = tanl_lfit*
s+dz_lfit;
5134 double delta_z = z-z_hough;
5136 cout<<
" delta_z = "<<delta_z<<
" ";
5137 double deltaZcut=10.;
5138 if(fabs(delta_z)<deltaZcut) {
5139 int idx_Vclu = (*i_map).first;
5140 myIterClu=myIterCgemClusterBegin+idx_Vclu;
5141 int layer = (*myIterClu)->getlayerid();
5142 if(fabs(delta_z)<fabs(cgemVClu_dz_min[layer])) {
5143 cgemVClu_dz_min[layer]=delta_z;
5144 cgemVClu_s_dzmin[layer]=
s;
5145 cgemVClu_z_dzmin[layer]=z;
5146 cgemVCluIdx[layer]=idx_Vclu;
5152 for(
int layer=2; layer>=0; layer--) {
5153 if(cgemVCluIdx[layer]!=-1) {
5154 iter_trkCandi->CgemVClusterID.push_back(cgemVCluIdx[layer]);
5161 cout<<
"MDC V-hits candidates: "<<endl;
5162 resetVDigiFitItem(*iter_trkCandi);
5164 iter_vvsz = vec_vecsz.begin();
5166 for(; iter_vvsz!=vec_vecsz.end(); iter_vvsz++, i_digi++)
5168 vector<pair<double, double> >& vec_sz = *iter_vvsz;
5169 vector<pair<double, double> >::iterator iter_sz = vec_sz.begin();
5170 double delta_z_min = 9999.0;
5171 double z_dzmin(0), s_dzmin(0);
5172 for(; iter_sz!=vec_sz.end(); iter_sz++) {
5173 double s = iter_sz->first;
5174 double z = iter_sz->second;
5175 double z_hough = tanl_lfit*
s+dz_lfit;
5176 double delta_z = z-z_hough;
5177 if(fabs(delta_z)<fabs(delta_z_min)) {
5178 delta_z_min=delta_z;
5184 cout<<
"delta_z_min="<<delta_z_min;
5185 double deltaZcut_mdc = 10.0;
5186 if(fabs(delta_z_min)>deltaZcut_mdc) {
5189 iter_trkCandi->mdcVDigiFitted[i_digi]=-1;
5192 iter_trkCandi->mdcVDigiFitted[i_digi]=1;
5194 cout<<
"("<<s_dzmin<<
","<<z_dzmin<<
")";
5199 if((i_digi+1)%3==0) cout<<endl;}
5202 if((i_digi+1)%3!=1) cout<<endl;
5206 int fitFlag = helixFit(*iter_trkCandi, 5);
5209 updateDigiFitItem(*iter_trkCandi);
5210 updateCgemFitItem(*iter_trkCandi, 1);
5225void DotsConnection::processTrkCandi()
5230 bool multiPath, isGoodTrk;
5232 vector<struct trkCandi>::iterator iter_trkCandi = myVecTrkCandidates.begin();
5234 for(; iter_trkCandi!=myVecTrkCandidates.end(); i_candi++, iter_trkCandi=myVecTrkCandidates.begin()+i_candi)
5237 cout<<endl<<
" ----> process trk candi "<<i_candi<<
": "<<endl;
5238 struct trkCandi& aTrkCandi=myVecTrkCandidates[i_candi];
5241 if(aTrkCandi.
isGood!=1)
continue;
5244 multiPath = split(aTrkCandi, i_candi);
5252 struct trkCandi& aTrkCandi_new=myVecTrkCandidates[i_candi];
5253 if(aTrkCandi_new.
isGood!=1)
break;
5255 isGoodTrk = completeTrkCandi(aTrkCandi_new, i_candi);
5259 int nPath=myVecNhitPathId.size();
5261 int nMax_goodHits=0;
5264 for(
int i=0; i<nPath; i++) {
5265 int nFittedXhits = myVecNhitPathId[i].first;
5266 if(nMax_goodHits<nFittedXhits) {
5267 nMax_goodHits=nFittedXhits;
5268 i_bestPath=myVecNhitPathId[i].second;
5272 myVecTrkPaths[i_bestPath].trkIdx = i_candi;
5273 myVecTrkCandidates[i_candi] = myVecTrkPaths[i_bestPath];
5277 myVecNhitPathId.erase(myVecNhitPathId.begin()+j_bestPath);
5279 if(myDebugNb==2) cout<<
" ~~~~~ try another path "<<endl;
5284 myVecTrkCandidates[i_candi+1].isGood=0;
5285 if(myDebugNb==2) cout<<
" ~~~~~ no more path "<<endl;
5292 rmPathInTrkCandi(myVecTrkCandidates[i_candi+1], myVecTrkCandidates[i_candi]);
5310void DotsConnection::processMdcHitCluter()
5312 bool debug=
false; debug=
true;
5313 int nMdcClu = myVecTrkCandidates.size();
5314 for(
int iMdcClu=0; iMdcClu<nMdcClu; iMdcClu++)
5316 if(debug) cout<<
"processMdcHitCluter "<<iMdcClu<<endl;
5317 struct trkCandi& aTrkCandi = myVecTrkCandidates.at(iMdcClu);
5326 vector<int> vecHit[8];
5327 for(
int layer=8; layer<43; layer++)
5329 int iArea = myAreaLayer[layer];
5330 int nHits = aTrkCandi.
mdcHitIdx[layer].size();
5331 for(
int iHit=0; iHit<nHits; iHit++)
5333 int idx = aTrkCandi.
mdcHitIdx[layer][iHit];
5337 vecHit[iArea].push_back(wireIdx);
5342 vector<struct trkCandi> vTrkSeg[8];
5343 for(
int iArea=2; iArea<=7; iArea++)
5345 if(vecHit[iArea].size()==0)
continue;
5346 vector<vector<int> > vecHitClu;
5347 mdcHitClustering(vecHit[iArea], vecHitClu);
5348 int nHitClu = vecHitClu.size();
5349 if(debug) cout<<endl<<
"Area "<<iArea<<
" has "<<vecHit[iArea].size()<<
" hits, "<<nHitClu<<
" hit-cluster(s)"<<endl;
5350 for(
int iHitClu=0; iHitClu<nHitClu; iHitClu++)
5352 if(debug) cout<<setw(4)<<
""<<
"hit-cluster "<<iHitClu<<
":"<<endl;
5353 vector<int>& aHitClu = vecHitClu.at(iHitClu);
5354 int nHit=aHitClu.size();
5355 set<int> aSetHits(aHitClu.begin(),aHitClu.end());
5361 vector<int> aPath(aSetHits.begin(), aSetHits.end());
5362 if(aPath.size()>3 && (iArea==2||iArea==7) ) {
5363 struct trkCandi oneTrkCandi=getTrkCandi(aPath);
5364 int nDropHits = circleFitTaubin(oneTrkCandi, 5,
false, 2,
true);
5366 if(nFittedXhits>3&&nDropHits<=3) {
5368 vTrkSeg[iArea].push_back(oneTrkCandi);
5369 setMinusVec(aSetHits,aPath);
5375 vector<vector<int> > VPathCA;
5377 if(VPathCA.size()>0) {
5378 for(
int iPath=0; iPath<VPathCA.size(); iPath++) {
5379 aPath=VPathCA.at(iPath);
5380 struct trkCandi oneTrkCandi=getTrkCandi(aPath);
5381 if(iArea==2||iArea==7) {
5382 int nDropHits = circleFitTaubin(oneTrkCandi, 5,
false, 2,
true);
5384 if(nFittedXhits>=3&&nDropHits<=3) {
5386 vTrkSeg[iArea].push_back(oneTrkCandi);
5387 setMinusVec(aSetHits,aPath);
5394 if(debug) cout<<
"No good circle fit from CA path, fill all hits as one (bad) trk candidate"<<endl;
5395 vector<int> aPath(aSetHits.begin(), aSetHits.end());
5396 vTrkSeg[iArea].push_back(getTrkCandi(aPath));
5401 if(aSetHits.size()==0)
break;
5413bool DotsConnection::completeTrkCandi(
struct trkCandi& aTrkCandi,
int idxCandi)
5416 if(myDebugNb==2) cout<<
" ~~~~~ completeTrkCandi "<<idxCandi<<endl;
5422 int iTryCircleFit=0;
5423 int nDropMin=9999;
int iCircleFitBest=-1;
5424 bool needUseWireCrossPoint=myUseWireCrossPoint;
5425 for(; iTryCircleFit<2; iTryCircleFit++) {
5426 int nDropHits = circleFitTaubin(aTrkCandi, 5,
false, 2,
true,-1,50,needUseWireCrossPoint);
5427 needUseWireCrossPoint=
false;
5428 if(nDropMin>nDropHits) {
5430 iCircleFitBest=iTryCircleFit;
5436 cout<<
"nDropHits="<<nDropHits<<
" < NActiveHits=myDotsHelixFitter.getNActiveHits(), nDropInnerXHits="<<nDropInnerXHits<<
", nActOuterXHits="<<nActOuterXHits<<endl;
5437 if(iTryCircleFit==0&&nActOuterXHits>0&&nDropInnerXHits>nActOuterXHits) {
5439 cout<<
"N_drop_innerXhits="<<nDropInnerXHits<<
" >= N_keep_outerXhits="<<nActOuterXHits<<
", try circle fit without outer X-hits!"<<endl;
5440 setFitFlagUncertain(aTrkCandi, 20, 60);
5445 else if(iTryCircleFit==0) {
5447 int nXInBig=setXHitsInBigSetFitFlagUncertain(aTrkCandi,
false);
5448 if(nXInBig>0&&nXInBig<nDropMin) {
5450 cout<<
"nDropHits>=NActiveHits, try circle fit without hits in big cluster"<<endl;
5451 nXInBig=setXHitsInBigSetFitFlagUncertain(aTrkCandi,
true);
5462 updateDigiFitItem(aTrkCandi);
5463 double circlePar[3];
5464 circlePar[0]=aTrkCandi.
a_helix[0];
5465 circlePar[1]=aTrkCandi.
a_helix[1];
5466 circlePar[2]=aTrkCandi.
a_helix[2];
5467 if(myDotsHelixFitter.
getNActiveHits()<4||circlePar[0]!=circlePar[0])
5470 cout<<
"a bad circleFitTaubin! "<<endl;
5472 resetXDigiFitItem(aTrkCandi);
5477 matchMdcDigi(aTrkCandi, idxCandi);
5478 if(aTrkCandi.
isGood==0)
return false;
5481 matchCgemXCluster(aTrkCandi);
5484 if(selectVHits(aTrkCandi)) {
5485 int fitFlag = helixFit(aTrkCandi, 5);
5488 if(myDebugNb==2) cout<<
"a good Helix fit"<<endl;
5489 updateDigiFitItem(aTrkCandi);
5490 updateCgemFitItem(aTrkCandi, 1);
5498bool DotsConnection::selectVHits(
struct trkCandi& aTrkCandi)
5503 cout<<
"CGEM V-cluster candidates (s,z): "<<endl;
5505 vector<int>::iterator iter_XClu = aTrkCandi.
CgemXClusterID.begin();
5507 myRoughTanlDzMap.Reset();
5508 map<int,pair<double,double> > map_Vclu_sz;
5511 int idx_Xclu = *iter_XClu;
5512 if(idx_Xclu<0)
continue;
5513 myIterClu=myIterCgemClusterBegin+idx_Xclu;
5514 int layer = (*myIterClu)->getlayerid();
5516 cout<<
"layer "<<layer;
5519 if(fabs(dphi)<1e-10) {
5521 cout<<
"DotsConnection::"<<__FUNCTION__<<
": error, no intersection between track and CGEM layer "<<layer<<endl;
5524 double s = fabs(dphi*circle_fittedHelix.
radius());
5526 int nV = myVecCgemVClusterIdx[layer].size();
5527 for(
int iV=0; iV<nV; iV++) {
5528 int idx_Vclu = myVecCgemVClusterIdx[layer][iV];
5529 if(getCgemClusterIntersection(idx_Xclu,idx_Vclu,z_XV)) {
5530 pair<double, double> sz(
s,z_XV);
5531 fillTanlDzMap(
s,z_XV);
5533 map_Vclu_sz[idx_Vclu]=sz;
5535 cout<<
" ("<<
s<<
", "<<z_XV<<
")";
5547 vector<vector<pair<double, double> > > vec_vecsz;
5549 int nVdigi = mdcVDigiWireIdx.size();
5550 vector<int>::iterator iter_digi = mdcVDigiWireIdx.begin();
5552 cout<<
"MDC V-hits candidates (s,z): "<<endl;
5554 for( ; i_digi<nVdigi; i_digi++,iter_digi++)
5556 int wireIdx = *iter_digi;
5557 const MdcDigi* aMdcDigi = myMapMdcDigi[wireIdx];
5558 vector<pair<double, double> > vec_sz = myDotsHelixFitter.
getSZ(aMdcDigi);
5559 vec_vecsz.push_back(vec_sz);
5560 vector<pair<double, double> >::iterator iter_sz = vec_sz.begin();
5561 for(; iter_sz!=vec_sz.end(); iter_sz++) {
5562 fillTanlDzMap(iter_sz->first,iter_sz->second);
5564 if(vec_sz.size()>0) nsz_mdc++;
5566 if((i_digi+1)%3==0) cout<<endl;
5568 if(myDebugNb==2)
if((i_digi+1)%3!=1) cout<<endl;
5571 if((nsz_cgem+nsz_mdc)<3) {
5573 cout<<
" no sufficient sz pair! stop tanl, dz, V-hits association. "<<endl;
5580 double x_peak(0.), y_peak(0.);
5581 double x_weight(0.), y_weight(0.);
5582 getWeightedPeak(myRoughTanlDzMap, x_peak, y_peak, x_weight, y_weight);
5585 cout<<
"from tanL-dz Hough map: dz="<<y_weight<<
", tanL = "<<x_weight<<endl;
5588 myRoughTanlDzMap.Reset();
5589 double s_sum(0), z_sum(0);
5590 double s2_sum(0), sz_sum(0);
5594 int cgemVCluIdx[3]={-1,-1,-1};
5595 double cgemVClu_dz_min[3]={9999.,9999.,9999.};
5596 double cgemVClu_s_dzmin[3]={9999.,9999.,9999.};
5597 double cgemVClu_z_dzmin[3]={9999.,9999.,9999.};
5599 cout<<
"CGEM V-cluster candidates: "<<endl;
5600 for(map<
int,pair<double,double> >::iterator i_map=map_Vclu_sz.begin(); i_map!=map_Vclu_sz.end(); i_map++)
5602 double s = (*i_map).second.first;
5603 double z = (*i_map).second.second;
5604 double z_hough = x_weight*
s+y_weight;
5605 double delta_z = z-z_hough;
5607 cout<<
" delta_z = "<<delta_z<<
" ";
5608 double deltaZcut=10.;
5609 if(fabs(delta_z)<deltaZcut) {
5610 int idx_Vclu = (*i_map).first;
5611 myIterClu=myIterCgemClusterBegin+idx_Vclu;
5612 int layer = (*myIterClu)->getlayerid();
5613 if(fabs(delta_z)<fabs(cgemVClu_dz_min[layer])) {
5614 cgemVClu_dz_min[layer]=delta_z;
5615 cgemVClu_s_dzmin[layer]=
s;
5616 cgemVClu_z_dzmin[layer]=z;
5617 cgemVCluIdx[layer]=idx_Vclu;
5623 for(
int layer=2; layer>=0; layer--) {
5624 if(cgemVCluIdx[layer]!=-1) {
5626 fillTanlDzMap(cgemVClu_s_dzmin[layer],cgemVClu_z_dzmin[layer]);
5627 s_sum+=cgemVClu_s_dzmin[layer];
5628 s2_sum+=cgemVClu_s_dzmin[layer]*cgemVClu_s_dzmin[layer];
5629 z_sum+=cgemVClu_z_dzmin[layer];
5630 sz_sum+=cgemVClu_s_dzmin[layer]*cgemVClu_z_dzmin[layer];
5637 cout<<
"MDC V-hits candidates: "<<endl;
5638 resetVDigiFitItem(aTrkCandi);
5639 vector<vector<pair<double, double> > >::iterator iter_vvsz = vec_vecsz.begin();
5641 for(; iter_vvsz!=vec_vecsz.end(); iter_vvsz++, i_digi++)
5643 vector<pair<double, double> >& vec_sz = *iter_vvsz;
5644 vector<pair<double, double> >::iterator iter_sz = vec_sz.begin();
5645 double delta_z_min = 9999.0;
5646 double z_dzmin(0), s_dzmin(0);
5647 for(; iter_sz!=vec_sz.end(); iter_sz++) {
5648 double s = iter_sz->first;
5649 double z = iter_sz->second;
5650 double z_hough = x_weight*
s+y_weight;
5651 double delta_z = z-z_hough;
5652 if(fabs(delta_z)<fabs(delta_z_min)) {
5653 delta_z_min=delta_z;
5659 cout<<
"delta_z_min="<<delta_z_min;
5660 double deltaZcut_mdc = 10.0;
5661 if(fabs(delta_z_min)>deltaZcut_mdc) {
5668 fillTanlDzMap(s_dzmin, z_dzmin);
5670 cout<<
"("<<s_dzmin<<
","<<z_dzmin<<
")";
5673 s2_sum+=s_dzmin*s_dzmin;
5674 sz_sum+=s_dzmin*z_dzmin;
5679 if((i_digi+1)%3==0) cout<<endl;}
5682 if((i_digi+1)%3!=1) cout<<endl;
5685 getWeightedPeak(myRoughTanlDzMap, x_peak, y_peak, x_weight, y_weight);
5687 cout<<
"After rough selection of sz of V-hits: "<<endl;
5689 cout<<
"from tanL-dz Hough map: dz="<<y_weight<<
", tanL = "<<x_weight<<endl;
5693 double tanl_lfit = (n_sz*sz_sum-s_sum*z_sum)/(n_sz*s2_sum-s_sum*s_sum);
5694 double dz_lfit = (s2_sum*z_sum-s_sum*sz_sum)/(n_sz*s2_sum-s_sum*s_sum);
5696 cout<<
"from s-z line fit: dz="<<dz_lfit<<
", tanL = "<<tanl_lfit<<endl;
5699 aTrkCandi.
a_helix[4]=tanl_lfit;
5703 for(
int i=0; i<3; i++) {
5705 cgemVClu_dz_min[i]=9999.;
5706 cgemVClu_s_dzmin[i]=9999.;
5707 cgemVClu_z_dzmin[i]=9999.;
5710 cout<<
"CGEM V-cluster candidates: "<<endl;
5711 for(map<
int,pair<double,double> >::iterator i_map=map_Vclu_sz.begin(); i_map!=map_Vclu_sz.end(); i_map++)
5713 double s = (*i_map).second.first;
5714 double z = (*i_map).second.second;
5715 double z_hough = tanl_lfit*
s+dz_lfit;
5716 double delta_z = z-z_hough;
5718 cout<<
" delta_z = "<<delta_z<<
" ";
5719 double deltaZcut=10.;
5720 if(fabs(delta_z)<deltaZcut) {
5721 int idx_Vclu = (*i_map).first;
5722 myIterClu=myIterCgemClusterBegin+idx_Vclu;
5723 int layer = (*myIterClu)->getlayerid();
5724 if(fabs(delta_z)<fabs(cgemVClu_dz_min[layer])) {
5725 cgemVClu_dz_min[layer]=delta_z;
5726 cgemVClu_s_dzmin[layer]=
s;
5727 cgemVClu_z_dzmin[layer]=z;
5728 cgemVCluIdx[layer]=idx_Vclu;
5734 for(
int layer=2; layer>=0; layer--) {
5735 if(cgemVCluIdx[layer]!=-1) {
5743 cout<<
"MDC V-hits candidates: "<<endl;
5744 resetVDigiFitItem(aTrkCandi);
5746 iter_vvsz = vec_vecsz.begin();
5748 for(; iter_vvsz!=vec_vecsz.end(); iter_vvsz++, i_digi++)
5750 vector<pair<double, double> >& vec_sz = *iter_vvsz;
5751 vector<pair<double, double> >::iterator iter_sz = vec_sz.begin();
5752 double delta_z_min = 9999.0;
5753 double z_dzmin(0), s_dzmin(0);
5754 for(; iter_sz!=vec_sz.end(); iter_sz++) {
5755 double s = iter_sz->first;
5756 double z = iter_sz->second;
5757 double z_hough = tanl_lfit*
s+dz_lfit;
5758 double delta_z = z-z_hough;
5759 if(fabs(delta_z)<fabs(delta_z_min)) {
5760 delta_z_min=delta_z;
5766 cout<<
"delta_z_min="<<delta_z_min;
5767 double deltaZcut_mdc = 10.0;
5768 if(fabs(delta_z_min)>deltaZcut_mdc) {
5776 cout<<
"("<<s_dzmin<<
","<<z_dzmin<<
")";
5781 if((i_digi+1)%3==0) cout<<endl;}
5784 if((i_digi+1)%3!=1) cout<<endl;
5786 if(nVhits>2)
return true;
5791int DotsConnection::fillTanlDzMap(
double s,
double z,
int vote)
5795 double dx = (2.0*myTanlRange)/(myNBinTanl*vote);
5796 dx = dx>1e-6?dx:1e-6;
5797 double x = -1.0*myTanlRange + dx/2;
5798 while(x<myTanlRange)
5800 double y = -
s*
x + z;
5802 cut = fabs(y)<myDzRange;
5805 myRoughTanlDzMap.Fill(x,y);
5813void DotsConnection::getWeightedPeak(TH2D& h,
double& x_peak,
double& y_peak,
double& x_weight,
double& y_weight,
int x_ext,
int y_ext)
5815 int ix_max, iy_max, iz_max;
5816 h.GetMaximumBin(ix_max, iy_max, iz_max);
5817 int nx=h.GetXaxis()->GetNbins();
5818 int ny=h.GetYaxis()->GetNbins();
5821 if(ix_max==0&&iy_max==0) {
5823 for(
int i=0; i<nx; i++)
5826 for(
int j=0; j<ny; j++)
5828 double n_tmp = h.GetBinContent(i,j);
5834 x_peak=h.GetXaxis()->GetBinCenter(ix_max);
5835 y_peak=h.GetYaxis()->GetBinCenter(iy_max);
5836 x_weight=0.; y_weight=0.;
5838 if(x_ext>=0&&y_ext>=0) {
5839 for(
int i1=ix_max-x_ext; i1<=ix_max+x_ext; i1++)
5841 double x_tmp = h.GetXaxis()->GetBinCenter(i1);
5846 for(
int i2=iy_max-y_ext; i2<=iy_max+y_ext; i2++)
5848 double n_tmp = h.GetBinContent(i1_tmp,i2);
5849 if(i2<1||i2>ny) n_tmp=0.;
5850 if(i1<1||i1>nx) n_tmp=0.;
5852 double y_tmp = h.GetYaxis()->GetBinCenter(i2);
5853 x_weight+=x_tmp*n_tmp;
5854 y_weight+=y_tmp*n_tmp;
5858 x_weight=x_weight/
weight;
5859 y_weight=y_weight/
weight;
5867void DotsConnection::getRhoTheta_bisection(vector<const MdcDigi*>& aVecMdcDigi)
5869 bool debug=
false; debug=
true;
5871 int nHits=aVecMdcDigi.size();
5873 vector<double> vx, vy, vr;
5874 vector<const MdcDigi*>::iterator i_digi = aVecMdcDigi.begin();
5875 for(
int ihit=0; i_digi!=aVecMdcDigi.end(); i_digi++,ihit++)
5880 int wireid = myMdcGeomSvc->
Wire(layer,wire)->
Id();
5881 vx.push_back(0.1*myRLayer[layer]*
cos(myWirePhi[layer][wire]));
5882 vy.push_back(0.1*myRLayer[layer]*
sin(myWirePhi[layer][wire]));
5883 vr.push_back(myMapMdcDigiDd[wireid]);
5884 if( (myNtProd&2) && ihit<100) {
5885 myNt_Xhit[ihit]=vx.back();
5886 myNt_Yhit[ihit]=vy.back();
5887 myNt_DDhit[ihit]=vr.back();
5892 double rho_center = 0.0;
5893 double theta_center = myThetaRange/2.0;
5894 double bin_rho=myRhoRange;
5895 double bin_theta=myThetaRange/2.0;
5896 double theta, rho, theta_min, rho_min, rho_max;
5900 cout<<
" ------ hough start : rho_center, theta_center = "<<rho_center<<
", "<<theta_center<<endl;
5902 theta_min=theta_center-bin_theta;
5903 rho_min=rho_center-bin_rho;
5906 int Npass[2][2]={0,0,0,0};
5907 for(
int iHit=0; iHit<nHits; iHit++) {
5908 int pass[2][2]={0,0,0,0};
5909 for(
int i=0; i<3; i++) {
5910 theta=theta_center+(i-1)*bin_theta;
5911 for(
int sign=-1; sign<2; sign+=2) {
5912 rho=vx[iHit]*
cos(theta)+vy[iHit]*
sin(theta)+sign*vr[iHit];
5913 int i_rho=(int)floor((rho-rho_min)/bin_rho);
5914 if(i_rho>=0&&i_rho<2) {
5915 if(i<2) pass[0][i_rho]=1;
5916 if(i>0) pass[1][i_rho]=1;
5920 for(
int i=0; i<3; i++) {
5921 rho=rho_center+(i-1)*bin_rho;
5922 for(
int sign=-1; sign<2; sign+=2) {
5923 double a=vx[iHit]*vx[iHit]+vy[iHit]*vy[iHit];
5924 double b=-2.0*(rho+sign*vr[iHit])*vx[iHit];
5925 double c=pow(rho+sign*vr[iHit],2)-vy[iHit]*vy[iHit];
5929 for(
int sign2=-1; sign2<2; sign2+=2) {
5930 double costh=(-b+sign2*d)/(2.0*a);
5931 if(fabs(costh)<=1) {
5933 int i_theta=(int)floor((theta-theta_min)/bin_theta);
5934 if(i_theta>=0&&i_theta<2) {
5935 if(i<2) pass[i_theta][0]=1;
5936 if(i>2) pass[i_theta][1]=1;
5943 Npass[0][0]+=pass[0][0];
5944 Npass[0][1]+=pass[0][1];
5945 Npass[1][0]+=pass[1][0];
5946 Npass[1][1]+=pass[1][1];
5953 for(
int i=0; i<4; i++) {
5956 cout<<
"["<<ii<<
","<<jj<<
"]="<<Npass[ii][jj]<<
" ";
5957 if(nAccu_peak<Npass[ii][jj]) {
5958 nAccu_peak=Npass[ii][jj];
5977 theta_center=(quarter_peak/2-1)*0.5*bin_theta;
5978 rho_center=(quarter_peak%2-1)*0.5*bin_rho;
5984 cout<<
"hough iter "<<i_iter<<
": rho_center, theta_center = "<<rho_center<<
", "<<theta_center<<endl;
5988 if(i_iter>10)
break;
5993 double x_sum(0), y_sum(0), x2_sum(0), y2_sum(0), xy_sum(0);
5997 x_sum=0.; y_sum=0.; x2_sum=0.; y2_sum=0.; xy_sum=0.;
5998 for(
int iHit=0; iHit<nHits; iHit++) {
6001 x2_sum+=pow(vx[iHit],2);
6002 y2_sum+=pow(vy[iHit],2);
6003 xy_sum+=vx[iHit]*vy[iHit];
6005 double Y=x_sum*x_sum-nHits*x2_sum;
6006 double X=nHits*xy_sum-x_sum*y_sum;
6008 rho=(x_sum*xy_sum-x2_sum*y_sum)/sqrt(X*X+Y*Y);
6020 myNt_rho=rho_center;
6021 myNt_theta=theta_center;
6022 myNtTrkSeg->write();
6026bool DotsConnection::LineFit_TLS(vector<const MdcDigi*>& aVecMdcDigi, vector<double>& vRho, vector<double>& vTheta, vector<vector<int> >& vvLR, vector<vector<double> >& vvS,
bool ini,
double rho_ini,
double theta_ini,
const vector<int>& LR_ini)
6028 bool debug=
false; debug=
true;
6030 bool goodLine=
false;
6032 const int nHitMax=5;
6034 int nHits=aVecMdcDigi.size();
6035 if(nHits<3)
return false;
6036 if(nHits>nHitMax) nHits=nHitMax;
6039 vector<double> vx, vy, vr;
6041 vector<const MdcDigi*>::iterator i_digi = aVecMdcDigi.begin();
6042 vector<int> vLRFloat;
6044 cout<<setw(6)<<
""<<
"---> Try fitting "<<nHits<<
" hits (L,W,dd): ";
6046 for(
int ihit=0; i_digi!=aVecMdcDigi.end(); i_digi++,ihit++)
6051 int wireid = myMdcGeomSvc->
Wire(layer,wire)->
Id();
6052 vx.push_back(0.1*myRLayer[layer]*
cos(myWirePhi[layer][wire]));
6053 vy.push_back(0.1*myRLayer[layer]*
sin(myWirePhi[layer][wire]));
6054 vr.push_back(myMapMdcDigiDd[wireid]);
6055 if(debug) cout<<
"("<<setw(2)<<layer <<
", "<<setw(3)<<wire<<
", "<<setw(8)<<vr.back()<<
") ";
6056 if( (myNtProd&2) && ihit<100) {
6057 myNt_layer[ihit]=layer;
6058 myNt_wire[ihit]=wire;
6059 myNt_Xhit[ihit]=vx.back();
6060 myNt_Yhit[ihit]=vy.back();
6061 myNt_DDhit[ihit]=vr.back();
6064 int LR=LR_ini.at(ihit);
6066 if(debug) cout<<
"LR="<<setw(2)<<LR<<
"; ";
6067 if(LR==0) vLRFloat.push_back(ihit);
6069 else vLR.push_back(0);
6071 if(debug) cout<<endl;
6087 double theta_TLS, rho_TLS, theta0_TLS(0.), rho0_TLS(0.), theta_last, rho_last;
6088 double x_sum, y_sum, x_mean, y_mean, dxdy_sum, dx2mdy2_sum, rho_hit[nHitMax];
6093 int nLRFloat = nHits;
6096 nLRFloat=vLRFloat.size();
6097 theta0_TLS=theta_ini;
6100 if(debug) cout<<setw(6)<<
""<<
"ini rho, theta = "<<rho0_TLS<<
", "<<theta0_TLS<<endl;
6102 int nLR=pow(2,nLRFloat);
6103 for(
int iLR=iLR_0; iLR<=nLR; iLR++) {
6106 for(
int j=0; j<nLRFloat; j++) {
6107 int ibHit = pow(2,j);
6109 if(ini) i_LR=vLRFloat.at(j);
6110 if( (iLR-1) & ibHit ) vLR[i_LR]= 1;
6116 theta_TLS=theta0_TLS;
6119 double drho2_tls=99999.;
6120 double drho2_min=999999.;
6121 double theta_TLS_drmin, rho_TLS_drmin;
6122 int i_lastMin(0), iPeriod(0), period(0);
6123 double drho2_last=0.;
6124 set<double> set_drho2;
6128 for(
int j=0; j<nHits; j++)
6132 double x=vx[j]+vr[j]*
cos(theta_TLS)*vLR[j];
6133 double y=vy[j]+vr[j]*
sin(theta_TLS)*vLR[j];
6136 rho_hit[j]=
x*
cos(theta_TLS)+y*
sin(theta_TLS);
6137 drho2_tls+=pow(rho_hit[j]-rho_TLS, 2);
6141 dxdy_sum=0.; dx2mdy2_sum=0.;
6142 for(
int j=0; j<nHits; j++)
6146 double x=vx[j]+vr[j]*
cos(theta_TLS)*vLR[j];
6147 double y=vy[j]+vr[j]*
sin(theta_TLS)*vLR[j];
6151 dx2mdy2_sum+=dy*dy-dx*dx;
6153 double Y_TLS=-2*dxdy_sum;
6154 double X_TLS=dx2mdy2_sum;
6155 theta_TLS=0.5*atan2(Y_TLS,X_TLS);
6156 rho_TLS=x_mean*
cos(theta_TLS)+y_mean*
sin(theta_TLS);
6159 theta0_TLS=theta_TLS;
6162 theta0_TLS+=acos(-1.);
6170 while(fabs(theta_TLS-theta0_TLS+acos(-1.))<fabs(theta_TLS-theta0_TLS))
6174 theta_TLS+=acos(-1.);
6178 while(fabs(theta_TLS-theta0_TLS-acos(-1.))<fabs(theta_TLS-theta0_TLS))
6182 theta_TLS-=acos(-1.);
6187 if(fabs(drho2_tls-drho2_last)<0.0001)
6189 if(debug) cout<<setw(6)<<
""<<
"Line TLS fit converges at iter "<<setw(2)<<
iter<<
" ";
6192 if(drho2_min>drho2_tls) {
6193 drho2_min=drho2_tls;
6194 theta_TLS_drmin=theta_TLS;
6195 rho_TLS_drmin=rho_TLS;
6204 bool findDrho2=
false;
6205 set<double>::iterator it_set=set_drho2.begin();
6206 for(; it_set!=set_drho2.end(); it_set++) {
6207 if(fabs(*it_set-drho2_tls)<0.0001) {
6214 if(debug) cout<<setw(6)<<
""<<setw(38)<<
"Repeat a result, break!";
6220 else set_drho2.insert(drho2_tls);
6225 if(debug) cout<<setw(6)<<
""<<setw(38)<<
"Too many iterations, break!";
6226 drho2_tls=drho2_min;
6227 rho_TLS=rho_TLS_drmin;
6228 theta_TLS=theta_TLS_drmin;
6231 theta_last=theta_TLS;
6232 drho2_last=drho2_tls;
6237 if(isGood_Drho2_line(drho2_tls)) {
6239 vRho.push_back(rho_TLS);
6240 vTheta.push_back(theta_TLS);
6241 vvLR.push_back(vLR);
6243 for(
int j=0; j<nHits; j++)
6245 double s=getS(rho_TLS,theta_TLS,vx[j],vy[j]);
6252 while(theta_TLS<-0.5*acos(-1.)) {
6253 theta_TLS+=acos(-1.);
6256 while(theta_TLS>0.5*acos(-1.)) {
6257 theta_TLS-=acos(-1.);
6260 cout<<setw(6)<<
""<<
"rho="<<setw(8)<<rho_TLS<<
", theta="<<setw(8)<<theta_TLS;
6262 for(
int j=0; j<nHits; j++) cout<<setw(3)<<std::right<<vLR[j];
6263 cout<<
"}, drho2="<<setw(8)<<drho2_tls;
6264 if(isGood_Drho2_line(drho2_tls)) {
6267 for(
int j=0; j<nHits; j++) cout<<setw(4)<<vS[j]<<
" ";
6282 myNtTrkSeg->write();
6289vector<vector<double> > DotsConnection::getTangentCircles(vector<const MdcDigi*>& aVecMdcDigi)
6291 vector<vector<double> > circles;
6292 int nHits=aVecMdcDigi.size();
6293 if(nHits<3)
return circles;
6295 vector<double> xhit, yhit, DDhit;
6296 vector<const MdcDigi*>::iterator i_digi = aVecMdcDigi.begin();
6297 for(
int ihit=0; ihit<3&&i_digi!=aVecMdcDigi.end(); i_digi++,ihit++)
6302 int wireid = myMdcGeomSvc->
Wire(layer,wire)->
Id();
6303 xhit.push_back(0.1*myRLayer[layer]*
cos(myWirePhi[layer][wire]));
6304 yhit.push_back(0.1*myRLayer[layer]*
sin(myWirePhi[layer][wire]));
6305 DDhit.push_back(myMapMdcDigiDd[wireid]);
6319 double xmean=(xhit[0]+xhit[1]+xhit[2])/3.0;
6320 double ymean=(yhit[0]+yhit[1]+yhit[2])/3.0;
6322 for(
int ic=0; ic<8; ic++) {
6323 double v11 = 2 * xhit[1] - 2 * xhit[0];
6324 double v12 = 2 * yhit[1] - 2 * yhit[0];
6325 double v13 = xhit[0] * xhit[0] - xhit[1] * xhit[1] + yhit[0] * yhit[0] - yhit[1] * yhit[1] - DDhit[0] * DDhit[0] + DDhit[1] * DDhit[1];
6326 double v14 = 2 * LR[ic][1] * DDhit[1] - 2 * LR[ic][0] * DDhit[0];
6328 double v21 = 2 * xhit[2] - 2 * xhit[1];
6329 double v22 = 2 * yhit[2] - 2 * yhit[1];
6330 double v23 = xhit[1] * xhit[1] - xhit[2] * xhit[2] + yhit[1] * yhit[1] - yhit[2] * yhit[2] - DDhit[1] * DDhit[1] + DDhit[2] * DDhit[2];
6331 double v24 = 2 * LR[ic][2] * DDhit[2] - 2 * LR[ic][1] * DDhit[1];
6333 double w12 = v12 / v11;
6334 double w13 = v13 / v11;
6335 double w14 = v14 / v11;
6337 double w22 = v22 / v21 - w12;
6338 double w23 = v23 / v21 - w13;
6339 double w24 = v24 / v21 - w14;
6341 double P = -w23 / w22;
6342 double Q = w24 / w22;
6343 double M = -w12 *
P - w13;
6344 double n = w14 - w12 * Q;
6346 double aa =
n *
n + Q * Q - 1;
6347 double bb = 2 * M *
n - 2 *
n * xhit[0] + 2 *
P * Q - 2 * Q * yhit[0] + 2 * LR[ic][0] * DDhit[0];
6348 double cc = xhit[0] * xhit[0] + M * M - 2 * M * xhit[0] +
P *
P + yhit[0] * yhit[0] - 2 *
P * yhit[0] - DDhit[0] * DDhit[0];
6350 double D = bb * bb - 4 * aa * cc;
6351 double rs = (-bb - sqrt(D)) / (2 * aa);
6353 double xs = M +
n * rs;
6354 double ys =
P + Q * rs;
6360 vector<double> circle;
6361 circle.push_back(xs);
6362 circle.push_back(ys);
6363 circle.push_back(rs);
6364 circle.push_back(xmean);
6365 circle.push_back(ymean);
6366 circles.push_back(circle);
6372vector<double> DotsConnection::getBestTangentCircle(vector<const MdcDigi*>& aVecMdcDigi)
6374 vector<double> circle;
6376 int nHits=aVecMdcDigi.size();
6377 if(nHits<=3)
return circle;
6411 int nFill=min(nHits,100);
6412 for(
int ihit=0; ihit<nFill; ihit++)
6414 Identifier id = aVecMdcDigi.at(ihit)->identify();
6417 int wireid = myMdcGeomSvc->
Wire(layer,wire)->
Id();
6418 double x=0.1*myRLayer[layer]*
cos(myWirePhi[layer][wire]);
6419 double y=0.1*myRLayer[layer]*
sin(myWirePhi[layer][wire]);
6420 double r=myMapMdcDigiDd[wireid];
6421 myNt_layer[ihit]=layer;
6422 myNt_wire[ihit]=wire;
6430 myNtTrkSeg->write();
6438void DotsConnection::fillVecCgemClu(
struct trkCandi& aTrkCandi)
6442 vector<int>::iterator iter_XClu = aTrkCandi.
CgemXClusterID.begin();
6445 int idx_Xclu = *iter_XClu;
6446 if(idx_Xclu<0)
continue;
6447 myIterClu=myIterCgemClusterBegin+idx_Xclu;
6448 int layer = (*myIterClu)->getlayerid();
6449 int sheet = (*myIterClu)->getsheetid();
6450 myVecCgemXCluIdx[layer][sheet].push_back(idx_Xclu);
6454 vector<int>::iterator iter_VClu = aTrkCandi.
CgemVClusterID.begin();
6457 int idx_clu = *iter_VClu;
6458 if(idx_clu<0)
continue;
6459 myIterClu=myIterCgemClusterBegin+idx_clu;
6460 int layer = (*myIterClu)->getlayerid();
6461 int sheet = (*myIterClu)->getsheetid();
6462 myVecCgemVCluIdx[layer][sheet].push_back(idx_clu);
6468int DotsConnection::NXcluToDrop(
struct trkCandi& aTrkCandi)
6471 vector<int>::iterator iter_XClu = aTrkCandi.
CgemXClusterID.begin();
6474 int idx_Xclu = *iter_XClu;
6475 if(idx_Xclu<0) nToDrop++;
6481bool DotsConnection::isGoodCgemCircleCandi(
struct trkCandi& aTrkCandi)
6484 HepVector a_helix(5,0);
6485 a_helix(1)=aTrkCandi.
a_helix[0];
6486 a_helix(2)=aTrkCandi.
a_helix[1];
6487 a_helix(3)=aTrkCandi.
a_helix[2];
6490 int n_sign[2]={0,0};
6493 Hep3Vector pos_center(center.x(), center.y());
6496 for(
int i=0; i<size; i++,
iter++) {
6503 myIterClu=myIterCgemClusterBegin+cluIdx;
6504 int layer = (*myIterClu)->getlayerid();
6506 double phi_cluster = (*myIterClu)->getrecphi();
6507 double x_clu = rCGEM*
cos(phi_cluster);
6508 double y_clu = rCGEM*
sin(phi_cluster);
6509 Hep3Vector pos_clu(x_clu, y_clu);
6510 Hep3Vector
vec_z=pos_clu.cross(pos_center);
6511 if(
vec_z.z()<0) n_sign[0]++;
6512 else if(
vec_z.z()>0) n_sign[1]++;
6515 if(n_sign[0]+n_sign[1]>0 && n_sign[0]*n_sign[1]==0)
return true;
6519bool DotsConnection::saveARecMdcTrack(
struct trkCandi& aTrkCandi)
6521 if(aTrkCandi.
isGood!=1)
return false;
6526 int trackId = myRecMdcTrackCol->size();
6530 HepVector a_helix(5,0);
6531 for(
int i=0; i<5; i++) {
6532 helixPar[i]=aTrkCandi.
a_helix[i];
6533 a_helix(i+1)=helixPar[i];
6540 int q = helixPar[2]>0? 1:-1;
6541 double pxy = aHelix.pt();
6542 double px = aHelix.momentum(0).x();
6543 double py = aHelix.momentum(0).y();
6544 double pz = aHelix.momentum(0).z();
6545 double p = aHelix.momentum(0).mag();
6546 double theta = aHelix.direction(0).theta();
6547 double phi = aHelix.direction(0).phi();
6550 double r = poca.perp();
6551 HepSymMatrix Ea = aHelix.Ea();
6553 double errorMat[15];
6555 for (
int ie = 0 ; ie < 5 ; ie ++){
6556 for (
int je = ie ; je < 5 ; je ++){
6557 errorMat[k] = Ea[ie][je];
6561 double chisq = aTrkCandi.
chi2;
6563 recMdcTrack->
setPxy(pxy);
6564 recMdcTrack->
setPx(px);
6565 recMdcTrack->
setPy(py);
6566 recMdcTrack->
setPz(pz);
6567 recMdcTrack->
setP(p);
6569 recMdcTrack->
setPhi(phi);
6571 recMdcTrack->
setX(poca.x());
6572 recMdcTrack->
setY(poca.y());
6573 recMdcTrack->
setZ(poca.z());
6574 recMdcTrack->
setR(sqrt(poca.x()*poca.x() + poca.y()*poca.y()));
6584 int maxLayerId = -1;
6585 int minLayerId = 43;
6587 double fltLen = -0.00001;
6588 int layerMaxFltLen=-1;
6590 fillVecCgemClu(aTrkCandi);
6592 map<int,int> clusterFitStat;
6593 myIterClu=myIterCgemClusterBegin;
6594 for(; myIterClu!=myIterCgemClusterEnd; myIterClu++)
6596 int flag = (*myIterClu)->getflag();
6597 if(flag!=2)
continue;
6598 int layer=(*myIterClu)->getlayerid();
6599 int sheet=(*myIterClu)->getsheetid();
6600 int idXClu = (*myIterClu)->getclusterflagb();
6602 vector<int>::iterator
iter=myVecCgemXCluIdx[layer][sheet].begin();
6604 for(;
iter!=myVecCgemXCluIdx[layer][sheet].end();
iter++, i_cluster++)
6608 if(!matchX)
continue;
6609 int idVClu = (*myIterClu)->getclusterflage();
6611 iter=myVecCgemVCluIdx[layer][sheet].begin();
6612 for(;
iter!=myVecCgemVCluIdx[layer][sheet].end();
iter++)
6613 if((*
iter)==idVClu) matchV=
true;
6618 clusterRefVec.push_back(recCgemCluster);
6619 clusterFitStat[clusterid] = 1;
6620 if(maxLayerId<layer)
6630 vector<RecMdcHit> aRecMdcHitVec = aTrkCandi.
RecMdcHitVec;
6631 vector<RecMdcHit>::iterator iter_recMdcHit = aRecMdcHitVec.begin();
6632 for(; iter_recMdcHit!=aRecMdcHitVec.end(); iter_recMdcHit++)
6634 if(iter_recMdcHit->getChisqAdd()>myMdcHitChi2Cut)
6638 recMdcHit->
setId(hitId);
6641 myRecMdcHitCol->push_back(recMdcHit);
6642 SmartRef<RecMdcHit> refHit(recMdcHit);
6643 hitRefVec.push_back(refHit);
6648 if(layer>maxLayerId)
6652 if(layer<minLayerId)
6656 if(fltLen<recMdcHit->getFltLen()) {
6658 layerMaxFltLen=layer;
6664 if(maxLayerId>=0&&maxLayerId<3) {
6666 fltLen=aHelix.flightLength(rmax);
6668 if(fltLen>0) fiTerm=-fltLen*
sin(theta)/aHelix.radius();
6670 cout<<
"fltLen<0!: maxLayerId="<<maxLayerId<<
", n_cluster="<<clusterRefVec.size()<<endl;
6671 cout<<
"myVecCgem1DCluster.size="<<myVecCgem1DCluster.size()<<endl;
6678 std::sort(clusterRefVec.begin(),clusterRefVec.end(),
sortCluster);
6682 myRecMdcTrackCol->push_back(recMdcTrack);
6687int DotsConnection::saveGoodTrkCandi()
6691 vector<struct trkCandi>::iterator iter_trkCandi = myVecTrkCandidates.begin();
6692 for(; iter_trkCandi!=myVecTrkCandidates.end(); iter_trkCandi++, i_candi++)
6694 if((*iter_trkCandi).isGood!=1)
continue;
6696 if(saveARecMdcTrack(*iter_trkCandi)) nSaved++;
6703 vector<const MdcDigi*> aMdcDigiVec;
6706 for(
int i=0; i<NUMBER_OF_DOTS_PER_CELL; i++)
6709 aMdcDigiVec.push_back(myMapMdcDigi[aCell.
dots[i]->
id]);
6733 for(
int i=0; i<NUMBER_OF_DOTS_PER_CELL; i++) {
6735 for(
int j=0; j<NUMBER_OF_DOTS_PER_CELL; j++) {
6757 for(
int i=0; i<NUMBER_OF_DOTS_PER_CELL; i++)
6783 bool twoDotsShare=
false;
6797 return twoDotsShare;
6802 cout<<
" cell "<<setw(5)<<aCell.
id;
6804 cout<<
": (L, W, id) = ";
6805 for(
int i=0; i<NUMBER_OF_DOTS_PER_CELL; i++) {
6809 cout<<
"("<<setw(3)<<layer<<
","<<setw(3)<<Cell<<
", "<<setw(5)<<aCell.
dots[i]->
id<<
") ";
6819 int a_shared[NUMBER_OF_DOTS_PER_CELL]={0,0,0};
6820 int b_shared[NUMBER_OF_DOTS_PER_CELL]={0,0,0};
6823 vector<int> viDotShared_a;
6824 vector<int> viDotShared_b;
6829 for(
int i=0; i<NUMBER_OF_DOTS_PER_CELL; i++) {
6831 for(
int j=0; j<NUMBER_OF_DOTS_PER_CELL; j++) {
6834 comDots.insert(id_a);
6835 viDotShared_a.push_back(i);
6836 viDotShared_b.push_back(j);
6846 cout<<setw(6)<<
""<<
"cannot make edge due to n_shared="<<n_shared<<endl;
6860 double diff_theta_min=9999.;
6861 double diff_rho_min =9999.;
6864 int iL_a_match(-1), iL_b_match(-1);
6865 int nReverse_b_match=0;
6867 for(
int iL_a=0; itL_a!=a.
vecFittedLine.end(); itL_a++, iL_a++) {
6868 double rho_a = (*itL_a).rho;
6869 double theta_a = (*itL_a).theta;
6871 for(
int iL_b=0; itL_b!=b.
vecFittedLine.end(); itL_b++, iL_b++) {
6872 double rho_b = (*itL_b).rho;
6873 double theta_b = (*itL_b).theta;
6875 while( fabs(theta_b-acos(-1.)-theta_a) < fabs(theta_b-theta_a) ) {
6880 while( fabs(theta_b+acos(-1.)-theta_a) < fabs(theta_b-theta_a) ) {
6885 if(diff_theta_min>fabs(theta_b-theta_a)) {
6886 diff_theta_min=fabs(theta_b-theta_a);
6887 diff_rho_min = fabs(rho_b-rho_a);
6889 iL_a_match=iL_a; iL_b_match=iL_b;
6890 nReverse_b_match=nReverse;
6894 if(diff_theta_min>0.1)
6897 cout<<setw(6)<<
""<<
"theta unmatched!"<<endl;
6900 else cout<<setw(6)<<
""<<
"theta matched!"<<endl;
6905 int nShared=viDotShared_a.size();
6906 cout<<setw(6)<<
""<<
"shared LR in a: (";
6907 for(
int i=0; i<nShared; i++) cout<<setw(3)<<LR_a.at(viDotShared_a.at(i));
6909 int reverse=1;
if(nReverse_b_match%2!=0) {reverse=-1; cout<<setw(6)<<
""<<
"LR reversed for b"<<endl;}
6910 cout<<setw(6)<<
""<<
"shared LR in b: (";
6911 for(
int i=0; i<nShared; i++) cout<<setw(3)<<reverse*LR_b.at(viDotShared_b.at(i));
6914 const double shortDD=0.16;
6915 for(
int i=0; i<nShared; i++) {
6916 if( reverse*LR_b.at(viDotShared_b.at(i)) != LR_a.at(viDotShared_a.at(i)) ) {
6917 if(myMapMdcDigiDd[b.
dots[viDotShared_b.at(i)]->
id]<shortDD)
continue;
6922 if(LR_match) cout<<setw(6)<<
""<<
"LR matched"<<endl;
6924 cout<<setw(6)<<
""<<
"LR unmatched"<<endl;
6929 map<int, int> mapWidxLR;
6930 vector<int>::iterator it =LR_a.begin();
6931 int idot_a_uncom=-1;
6932 for(
int i=0; i<NUMBER_OF_DOTS_PER_CELL; i++, it++) {
6933 mapWidxLR[a.
dots[i]->
id]=*it;
6941 int idot_b_uncom=-1;
6942 for(
int i=0; i<NUMBER_OF_DOTS_PER_CELL; i++, it++) {
6943 mapWidxLR[b.
dots[i]->
id]=(*it)*reverse;
6964 vector<dot*>::iterator it_dot=e.
dots.begin();
6965 for(; it_dot!=e.
dots.end(); it_dot++) {
6966 int idx = (*it_dot)->id;
6968 if(myMapMdcDigiDd[idx]>=shortDD) LR_i=mapWidxLR[idx];
6969 LR_ini.push_back(LR_i);
6977 vector<double> vRho; vector<double> vTheta; vector<vector<int> > vvLR; vector<vector<double> > vvS;
6979 bool hasGdLine = LineFit_TLS(aVecMdcDigi, vRho,vTheta,vvLR,vvS,
true,rho_ini,theta_ini,LR_ini);
6981 cout<<setw(6)<<
""<<
" no good line fit for the edge"<<endl;
6990 int n_rho = vRho.size();
6991 cout<<setw(6)<<
""<<
"an good edge with "<<n_rho<<
" fitted lines"<<endl;
6992 for(
int i_line=0; i_line<n_rho; i_line++) {
6994 aLine.
rho=vRho[i_line];
6995 aLine.
theta=vTheta[i_line];
6996 aLine.
LR=vvLR[i_line];
6997 aLine.
S=vvS[i_line];
7001 if(aLine.
order.front()>1&&aLine.
order.back()>1)
7004 cout<<setw(6)<<
""<<
"good S order : ";
7005 for(
int i=0; i<aLine.
order.size(); i++) cout<<setw(2)<<aLine.
order.at(i);
7007 if(aLine.
order.front()==3) {
7013 else if(aLine.
order.front()!=2) {
7014 cout<<setw(6)<<
""<<
" something wrong about the dots' order"<<endl;
7026 vector<const MdcDigi*> aMdcDigiVec;
7028 vector<dot*>::iterator it_dot=aEdge.
dots.begin();
7029 for(; it_dot!=aEdge.
dots.end(); it_dot++)
7031 aMdcDigiVec.push_back(myMapMdcDigi[(*it_dot)->id]);
7038 bool debug=
false; debug=
true;
7040 vector<int> v_rootCellIdx;
7042 vector<vector<int> > v_path;
7043 if(v_cellIdx.size()==0)
return v_path;
7046 for(vector<int>::iterator it_cellIdx=v_cellIdx.begin(); it_cellIdx!=v_cellIdx.end(); it_cellIdx++) {
7047 cell& aCell=v_cell.at(*it_cellIdx);
7050 if(aCell.
set_inEdge.empty()) v_rootCellIdx.push_back(*it_cellIdx);
7056 for(vector<int>::iterator it_cellIdx=v_cellIdx.begin(); it_cellIdx!=v_cellIdx.end(); it_cellIdx++) {
7057 cell& aCell=v_cell.at(*it_cellIdx);
7059 if(aCell.
color==0) {
7060 if(debug) cout<<setw(6)<<
""<<
"== start state calculation with cell "<<aCell.
id<<endl;
7064 if(debug) cout<<setw(6)<<
""<<
"The graph has a cycle? "<<
myNCycle<<endl;
7075 if(debug) cout<<setw(6)<<
""<<
"~~~~> to find a path: "<<endl;
7076 cell& aRootCell = v_cell.at(*it_cell_maxStat);
7079 map<int, pair<double,double> > setWireHitPos;
7080 map<int, int> setWireEntry;
7081 if(aRootCell.
state==1) {
7083 double rho=(*it_line).rho;
7084 double theta=(*it_line).theta;
7085 vector<double>& vS=(*it_line).S;
7086 for(
int i=0; i<NUMBER_OF_DOTS_PER_CELL; i++) {
7087 int wireIdx = aRootCell.
dots[i]->
id;
7090 getXY(rho,theta,
s,
x,y);
7091 map<int, int>::iterator it_entry;
7092 it_entry = setWireEntry.find(wireIdx);
7093 if(it_entry!=setWireEntry.end()) {
7095 map<int, pair<double,double> >::iterator it_pos=setWireHitPos.find(wireIdx);
7096 it_pos->second.first+=
x;
7097 it_pos->second.second+=y;
7100 setWireEntry[wireIdx]=1;
7101 pair<double,double> xy(
x,y);
7102 setWireHitPos[wireIdx]=xy;
7107 else nextCellFinding(aRootCell, v_cell, v_edge, setWireHitPos, setWireEntry);
7108 cout<<setw(6)<<
""<<
"a path found with "<<setWireEntry.size()<<
" dots: ";
7109 vector<int> aTrkPath;
7112 for(map<int, int>::iterator it_entry=setWireEntry.begin(); it_entry!=setWireEntry.end(); it_entry++) {
7113 int wireIdx = it_entry->first;
7114 aTrkPath.push_back(wireIdx);
7115 int layer = myMdcGeomSvc->
Wire(wireIdx)->
Layer();
7116 int wire = myMdcGeomSvc->
Wire(wireIdx)->
Cell();
7117 cout<<
"("<<layer<<
","<<wire<<
") ";
7129 v_path.push_back(aTrkPath);
7146 bool debug=
false; debug=
true;
7149 cout<<
"[ discover";
print(aCell);
7150 cout<<setw(5)<<
""<<
"loop next cells:"<<endl;
7154 for(set<int>::iterator it_edgeIdx=outEdges.begin(); it_edgeIdx!=outEdges.end(); it_edgeIdx++) {
7155 cell& aNextCell = v_cell.at(v_edge.at(*it_edgeIdx).twoCells.second);
7156 if(debug) cout<<setw(5)<<
""<<
"cell "<<aNextCell.
id<<endl;
7159 if(aNextCell.
state>state_max) state_max=aNextCell.
state;
7162 aCell.
state+=state_max;
7175 cout<<
"finish cell "<<aCell.
id<<
" with state "<<aCell.
state;
7184 bool debug=
false; debug=
true;
7188 cout<<setw(5)<<
""<<
"cell "<<aCell.
id<<
", state="<<aCell.
state;
7189 if(outEdges.size()) cout<<
", to fine the next cell";
7190 else cout<<
", no next cell"<<endl;
7192 for(set<int>::iterator it_edgeIdx=outEdges.begin(); it_edgeIdx!=outEdges.end(); it_edgeIdx++) {
7193 edge& aEdge = v_edge.at(*it_edgeIdx);
7195 double rho=(*it_line).rho;
7196 double theta=(*it_line).theta;
7197 vector<double>& vS=(*it_line).S;
7198 vector<dot*>::iterator it_dot=aEdge.
dots.begin();
7199 for(
int i=0; it_dot!=aEdge.
dots.end(); it_dot++, i++) {
7200 int wireIdx = (*it_dot)->id;
7203 getXY(rho,theta,
s,
x,y);
7204 map<int, int>::iterator it_entry;
7205 it_entry = setWireEntry.find(wireIdx);
7206 if(it_entry!=setWireEntry.end()) {
7208 map<int, pair<double,double> >::iterator it_pos=setWireHitPos.find(wireIdx);
7209 it_pos->second.first+=
x;
7210 it_pos->second.second+=y;
7213 setWireEntry[wireIdx]=1;
7214 pair<double,double> xy(
x,y);
7215 setWireHitPos[wireIdx]=xy;
7220 cell& aNextCell = v_cell.at(v_edge.at(*it_edgeIdx).twoCells.second);
7223 if(debug) cout<<
", cell "<<aNextCell.
id<<
" found, state = "<<aNextCell.
state<<endl;
7224 nextCellFinding(aNextCell, v_cell, v_edge, setWireHitPos, setWireEntry);
7240 int nHit=vecWireIdx.size();
7241 for(
int i=0; i<nHit; i++) {
7242 int wireIdx=vecWireIdx.at(i);
7252 map<int, dot>::iterator it_mapDots =
myMapDots.begin();
7253 for(; it_mapDots!=
myMapDots.end(); it_mapDots++)
7255 int wireIdx=it_mapDots->first;
7256 int layer= myMdcGeomSvc->
Wire(wireIdx)->
Layer();
7257 int Cell = myMdcGeomSvc->
Wire(wireIdx)->
Cell();
7260 vector<int> vec_neib = myMdcDigiGroup[wireIdx].
GetNeiborHits(1);
7261 for(vector<int>::iterator it_nb1=vec_neib.begin(); it_nb1!=vec_neib.end(); it_nb1++)
7264 int wireIdx1=*it_nb1;
7266 int layer1= myMdcGeomSvc->
Wire(wireIdx1)->
Layer();
7267 if(myWireFlag[layer1]!=myWireFlag[layer])
continue;
7268 int Cell1 = myMdcGeomSvc->
Wire(wireIdx1)->
Cell();
7271 vector<int>::iterator it_nb2=it_nb1; it_nb2++;
7273 if(it_nb2==vec_neib.end())
break;
7274 int wireIdx2=*it_nb2;
7276 int layer2= myMdcGeomSvc->
Wire(wireIdx2)->
Layer();
7277 if(myWireFlag[layer2]!=myWireFlag[layer]) {it_nb2++;
continue;}
7278 int Cell2 = myMdcGeomSvc->
Wire(wireIdx2)->
Cell();
7279 if(debug) cout<<
" === try cell (id, L, W): "
7280 <<
"("<<setw(5)<<wireIdx1<<
","<<setw(3)<<layer1<<
","<<setw(3)<<Cell1<<
") "
7281 <<
"("<<setw(5)<<wireIdx<<
","<<setw(3)<<layer<<
","<<setw(3)<<Cell<<
") "
7282 <<
"("<<setw(5)<<wireIdx2<<
","<<setw(3)<<layer2<<
","<<setw(3)<<Cell2<<
") "
7290 int wireMid = aCell.
dots[1]->
id;
7291 int nCell_wireMid =
myMapDots[wireMid].vecCellIdx.size();
7292 bool cellExist=
false;
7293 if(debug) cout<<
" wireMid="<<wireMid<<
" with "<<nCell_wireMid<<
" cells"<<endl;
7294 for(
int iCell=0; iCell<nCell_wireMid; iCell++) {
7295 if(debug) cout<<
" i-cell "<<iCell<<endl;
7299 if(debug) cout<<
" the cell exists"<<endl;
7302 else if(debug) cout<<
" the cell is new"<<endl;
7310 vector<double> vRho; vector<double> vTheta; vector<vector<int> > vvLR; vector<vector<double> > vvS;
7311 bool hasGdLine = LineFit_TLS(aVecMdcDigi, vRho,vTheta,vvLR,vvS);
7313 int n_rho = vRho.size();
7314 if(debug) cout<<
" the cell is good with N_lines="<<setw(2)<<n_rho<<endl;
7315 for(
int i_line=0; i_line<n_rho; i_line++) {
7317 aLine.
rho=vRho[i_line];
7318 aLine.
theta=vTheta[i_line];
7319 aLine.
LR=vvLR[i_line];
7320 aLine.
S=vvS[i_line];
7321 if(debug) cout<<
" S={"
7322 <<setw(5)<<aLine.
S[0]<<
" "
7323 <<setw(5)<<aLine.
S[1]<<
" "
7324 <<setw(5)<<aLine.
S[2]<<
"}"
7327 if(aLine.
S[0]<aLine.
S[1]&&aLine.
S[2]<aLine.
S[1]) {
7328 if(debug) cout<<
" line "<<i_line<<
" has a min S_middle"<<endl;
7331 if(aLine.
S[0]>aLine.
S[1]&&aLine.
S[2]>aLine.
S[1]) {
7332 if(debug) cout<<
" line "<<i_line<<
" has a max S_middle"<<endl;
7339 if(fabs(aLine.
S[0])<fabs(aLine.
S[2])) {
7340 aLine.
order.push_back(0);
7341 aLine.
order.push_back(1);
7342 aLine.
order.push_back(2);
7345 aLine.
order.push_back(2);
7346 aLine.
order.push_back(1);
7347 aLine.
order.push_back(0);
7349 if(debug) cout<<
" line "<<i_line<<
" hits' order "<<aLine.
order[0]<<
" "<<aLine.
order[1]<<
" "<<aLine.
order[2]<<endl;
7352 myVecCells.back().vecFittedLine.push_back(aLine);
7354 if(
myVecCells.back().vecFittedLine.size()>0) {
7356 cout<<setw(6)<<
""<<
"# Good Cell "<<
myVecCells.back().id<<endl;
7370 vector<int>::iterator i2_goodCell;
7375 for(
int idot=0; idot<3; idot++) {
7377 vector<int>& vCellIdx=
myMapDots[cell_a.
dots[idot]->
id].vecGoodCellIdx[ivc];
7378 i2_goodCell=vCellIdx.begin();
7379 for(; i2_goodCell!=vCellIdx.end(); i2_goodCell++) {
7380 int iGdCell2=*i2_goodCell;
7381 if(iGdCell2<iGdCell)
continue;
7384 cout<<
" ===> compare two cells: "<<endl;
7385 cout<<setw(6)<<
""<<
"Good cell "<<iGdCell<<
", ";
print(cell_a);
7386 cout<<setw(6)<<
""<<
"Good cell "<<*i2_goodCell<<
", ";
print(cell_b);
7391 if(debug) cout<<setw(6)<<
""<<
"neighbour with two shared dots in the middle"<<endl;
7399 idx_cell_fwd = *i1_goodCell;
7400 if(debug) cout<<setw(6)<<
""<<
"two cells' order reversed"<<endl;
7409 if(debug) cout<<setw(6)<<
""<<
"cell "<<cell_fwd->
id<<
" tagged as a new forward cell of cell "<<cell_bkd->
id<<endl;
7411 if(
makeEdge(cell_a,cell_b,aNewEdge)) {
7418 if(debug) cout<<setw(6)<<
""<<
"a good out edge added to cell "<<aNewEdge.
twoCells.first<<
" and "<<aNewEdge.
twoCells.second<<endl;
7434 int n_path = v_path.size();
7441 cout<<setw(6)<<
""<<n_path<<
" longest paths found"<<endl;
7443 else cout<<setw(6)<<
""<<
"no path found"<<endl;
7448 vector<int> vecWireIdx(setWireIdx.begin(),setWireIdx.end());
double sin(const BesAngle a)
double cos(const BesAngle a)
double P(RecMdcKalTrack *trk)
bool sortCluster(const RecCgemCluster *clusterA, const RecCgemCluster *clusterB)
bool global_show_this_event_detail
bool sortCluster(const RecCgemCluster *clusterA, const RecCgemCluster *clusterB)
double abs(const EvtComplex &c)
****INTEGER imax DOUBLE PRECISION m_pi *DOUBLE PRECISION m_amfin DOUBLE PRECISION m_Chfin DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_sinw2 DOUBLE PRECISION m_GFermi DOUBLE PRECISION m_MfinMin DOUBLE PRECISION m_ta2 INTEGER m_out INTEGER m_KeyFSR INTEGER m_KeyQCD *COMMON c_Semalib $ !copy of input $ !CMS energy $ !beam mass $ !final mass $ !beam charge $ !final charge $ !smallest final mass $ !Z mass $ !Z width $ !EW mixing angle $ !Gmu Fermi $ alphaQED at q
*********Class see also m_nmax DOUBLE PRECISION m_MasPhot DOUBLE PRECISION m_phsu DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_r2 DOUBLE PRECISION m_WtMass INTEGER m_nmax INTEGER m_Nevgen INTEGER m_IsFSR INTEGER m_MarTot *COMMON c_KarFin $ !Output file $ !Event serial number $ !alpha QED at Thomson limit $ !minimum energy at CMS for remooval $ !infrared cut
*********Class see also m_nmax DOUBLE PRECISION m_MasPhot DOUBLE PRECISION m_phsu DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_r2 DOUBLE PRECISION m_WtMass INTEGER m_nmax INTEGER m_Nevgen INTEGER m_IsFSR INTEGER m_MarTot *COMMON c_KarFin $ !Output file $ !Event serial number $ !alpha QED at Thomson limit $ !minimum energy at CMS for remooval $ !infrared dimensionless $ !dummy photon IR regulator $ !crude photon multiplicity enhancement factor *EVENT $ !MC crude volume of PhhSpace *Sfactors $ !YFS formfactor IR part only $ !YFS formfactor non IR finite part $ !mass weight
**********Class see also m_nmax DOUBLE PRECISION m_amel DOUBLE PRECISION m_x2 DOUBLE PRECISION m_alfinv DOUBLE PRECISION m_Xenph INTEGER m_KeyWtm INTEGER m_idyfs DOUBLE PRECISION m_zini DOUBLE PRECISION m_q2 DOUBLE PRECISION m_Wt_KF DOUBLE PRECISION m_WtCut INTEGER m_KFfin *COMMON c_KarLud $ !Input CMS energy[GeV] $ !CMS energy after beam spread beam strahlung[GeV] $ !Beam energy spread[GeV] $ !z boost due to beam spread $ !electron beam mass *ff pair spectrum $ !minimum v
std::vector< MdcDigi * > MdcDigiVec
ObjectVector< RecMdcHit > RecMdcHitCol
SmartRefVector< RecCgemCluster > ClusterRefVec
ObjectVector< RecMdcTrack > RecMdcTrackCol
SmartRefVector< RecMdcHit > HitRefVec
double getZFromPhiV(double phi, double V, int checkXRange=1) const
bool OnThePlane(double phi, double z) const
CgemGeoReadoutPlane * getReadoutPlane(int iLayer, int iSheet) const
void stateAutomaton(cell &aCell, vector< cell > &v_cell, vector< edge > &v_edge)
void longestPathFromCA(vector< int > &vecWireIdx, vector< vector< int > > &v_path)
DotsConnection(const std::string &name, ISvcLocator *pSvcLocator)
bool forward(cell &a, cell &b)
bool dotComp(dot *dot1, dot *dot2)
int numSharedDots(cell &a, cell &b)
void nextCellFinding(cell &aCell, vector< cell > &v_cell, vector< edge > &v_edge, map< int, pair< double, double > > &setWireHitPos, map< int, int > &setWireEntry)
bool makeEdge(cell &a, cell &b, edge &e)
bool sameCell(cell &a, cell &b)
vector< int > myVecGoodCellIdx
vector< cell > myVecCells
vector< int > myCellMaxState
vector< edge > myVecEdges
map< int, dot > myMapDots
void makeDotCellEdge(vector< int > &vecWireIdx)
vector< const MdcDigi * > getMdcDigiVec(cell &aCell)
bool isNeighbour(cell &a, cell &b)
vector< vector< int > > CellAutomaton(vector< int > &v_cellIdx, vector< cell > &v_cell, vector< edge > &v_edge)
void setInitialHelix(KalmanFit::Helix aHelix)
void setChi2Diverge(double cut=1000000)
int activeHits(double chi_cut=10, int sel=0)
activate myMdcDigiIsActive's hit based on myChiVecMdcDigi
int getLayer(const MdcDigi *aDcDigi)
void setDChits(vector< const MdcDigi * > aVecMdcDigi, double T0)
void setCgemClusters(vector< const RecCgemCluster * > aVecCgemCluster)
void useAxialHitsOnly(bool x=true)
bool getPosAtCgem(int layer, HepPoint3D &pos)
int deactiveHits(double chi_cut=10, int nMax=1, int layerMin=0, int layerMax=50)
bool setDChitsPhiAmbi(const vector< double > &vec)
vector< double > calculateCirclePar_Taubin(bool useIniHelix=false, int layerMin=-1, int layerMax=50, bool useExtraPoints=false)
double getRmidGapCgem(int i)
vector< double > getVecChiCgemCluster()
int getNActiveOuterXHits()
const HepSymMatrix & getEa()
KalmanFit::Helix getClassHelix()
void setModeWireCrossPointPersistent(bool mode)
vector< pair< double, double > > getSZ(const MdcDigi *aDcDigi)
return S, the track length in XY; and Z of a hit on a stereo wire
bool setDChitsFitFlag(vector< int > &vec)
set myMdcDigiIsActive from vec
void fitCircleOnly(bool x=true)
void setExtraPoints(vector< pair< double, double > > extraPoints)
void usePhiAmbi(bool x=true)
vector< RecMdcHit > makeRecMdcHitVec(int sel=1)
vector< double > getVecChiMdcDigi()
vector< int > getVecMdcDigiIsAct()
double IntersectCylinder(double r)
void updateDcDigiInfo(const MdcDigi *aDcDigi)
based on input aDcDigi, modifies myFlightLength, myLeftRight, myPosOnWire, myDocaFromDigi,...
vector< int > getVecIsActiveCgemCluster()
void setPxy(const double pxy)
void setTrackId(const int trackId)
void setPy(const double py)
void setZ(const double z)
void setX(const double x)
void setError(double err[15])
void setTheta(const double theta)
void setStat(const int stat)
void setP(const double p)
void setHelix(double helix[5])
void setPoca(double poca[3])
void setR(const double r)
void setCharge(const int charge)
void setY(const double y)
void setChi2(const double chi)
void setPhi(const double phi)
void setPz(const double pz)
void setPx(const double px)
double flightLength(HepPoint3D &hit) const
Hep3Vector momentum(double dPhi=0.) const
returns momentum vector after rotating angle dPhi in phi direction.
const HepSymMatrix & Ea(void) const
returns error matrix.
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
double radius(void) const
returns radious of helix.
const HepVector & a(void) const
returns helix parameters.
const HepPoint3D & pivot(void) const
returns pivot position.
Hep3Vector direction(double dPhi=0.) const
returns direction vector after rotating angle dPhi in phi direction.
double getT0(int layid, int cellid) const
double driftTimeToDist(double drifttime, int layid, int cellid, int lr, double entrance=0.0) const
double getTprop(int lay, double z) const
double getTimeWalk(int layid, double Q) const
vector< int > GetInner(int i=0)
void AddOuter(int wireID, int gap=0)
void SetWire(const MdcGeoWire *aWire)
vector< int > GetOuter(int i=0)
void AddInner(int wireID, int gap=0)
void AddPrev(int wireID, int gap=0)
void SetUsedFlag(int flag=1)
void AddHit(const MdcDigi *aHit)
void AddNext(int wireID, int gap=0)
vector< int > GetNeiborHits()
double Radius(void) const
double nomShift(void) const
int GetNeighborIDType2(void) const
vector< int > GetNeighborIDType4(void) const
HepPoint3D Forward(void) const
int GetNeighborIDType1(void) const
vector< int > GetNeighborIDType3(void) const
const MdcGeoWire *const Wire(unsigned id)
HepPoint3D getAdjacentIntersectionPointSlower(const MdcGeoWire *wireA, const MdcGeoWire *wireB)
const int getSuperLayerSize()
const MdcGeoLayer *const Layer(unsigned id)
static int layer(const Identifier &id)
Values of different levels (failure returns 0)
static int wire(const Identifier &id)
MdcDigiVec & getMdcDigiVec(uint32_t control=0)
static double MdcTime(int timeChannel)
virtual Identifier identify() const
unsigned int getChargeChannel() const
unsigned int getTimeChannel() const
int getclusterid(void) const
int getlayerid(void) const
const Identifier getMdcId(void) const
const double getFltLen(void) const
void setVecClusters(ClusterRefVec vecclusters)
void setPivot(const HepPoint3D &pivot)
void setVecHits(HitRefVec vechits)
void setFiTerm(double fiterm)
CalibPairCol::const_iterator PairIt
vector< Line > vecFittedLine
dot * dots[NUMBER_OF_DOTS_PER_CELL]
vector< Line > vecFittedLine
pair< int, int > twoCells
pair< int, int > twoLines
vector< RecMdcHit > RecMdcHitVec
vector< int > mdcV2DigiWireIdx
vector< int > mdcXDigiFitted
vector< double > mdcVDigiChi
vector< int > mdcXDigiWireIdx
vector< int > CgemXClusterID
vector< double > mdcXDigiChi
vector< int > mdcVDigiWireIdx
vector< int > mdcVDigiFitted
vector< int > mdcHitIdx[43]
vector< int > mdcV1DigiWireIdx
vector< int > CgemVClusterID