DotsConnection Class Reference

#include <DotsConnection.h>

Inheritance diagram for DotsConnection:

Classes
struct	cell
struct	edge

Public Member Functions
	DotsConnection (const std::string &name, ISvcLocator *pSvcLocator)
StatusCode	initialize ()
StatusCode	execute ()
StatusCode	finalize ()
bool	dotComp (dot *dot1, dot *dot2)
vector< const MdcDigi * >	getMdcDigiVec (cell &aCell)
int	numSharedDots (cell &a, cell &b)
bool	sameCell (cell &a, cell &b)
bool	forward (cell &a, cell &b)
bool	isNeighbour (cell &a, cell &b)
void	print (cell &aCell)
bool	makeOrder (Line &l)
bool	makeEdge (cell &a, cell &b, edge &e)
vector< const MdcDigi * >	getMdcDigiVec (edge &aEdge)
vector< vector< int > >	CellAutomaton (vector< int > &v_cellIdx, vector< cell > &v_cell, vector< edge > &v_edge)
void	stateAutomaton (cell &aCell, vector< cell > &v_cell, vector< edge > &v_edge)
void	nextCellFinding (cell &aCell, vector< cell > &v_cell, vector< edge > &v_edge, map< int, pair< double, double > > &setWireHitPos, map< int, int > &setWireEntry)
void	makeDotCellEdge (vector< int > &vecWireIdx)
void	longestPathFromCA (vector< int > &vecWireIdx, vector< vector< int > > &v_path)
void	longestPathFromCA (set< int > &setWireIdx, vector< vector< int > > &v_path)

Public Attributes
int	myNCycle
vector< int >	myCellMaxState
map< int, dot >	myMapDots
vector< cell >	myVecCells
vector< int >	myVecGoodCellIdx
vector< edge >	myVecEdges

Detailed Description

Definition at line 208 of file DotsConnection.h.

Constructor & Destructor Documentation

DotsConnection()

DotsConnection::DotsConnection	(	const std::string &	name,
		ISvcLocator *	pSvcLocator )

Definition at line 27 of file DotsConnection.cxx.

                                                                              :
    Algorithm(name, pSvcLocator)
{
    myMdcCalibFunSvc=NULL;
 
    // Part 1: Declare the properties
    //declareProperty("MyInt", m_myInt);
    declareProperty("Ntuple", myNtProd=0);
    declareProperty("driftTimeUpLimit",  myDriftTimeUpLimit = 400);
    declareProperty("MdcHitChi2Cut", myMdcHitChi2Cut = 100);
    declareProperty("ChiCut_circle", myChiCut_circle=5);
    declareProperty("NmaxDeact_circle", myNmaxDeact_circle=1);
    declareProperty("ChiCut_helix", myChiCut_helix=5);
    declareProperty("NmaxDeact_helix", myNmaxDeact_helix=1);
    declareProperty("Debug", myDebug=0);
    declareProperty("DebugNb", myDebugNb=0);
    declareProperty("Chi2CutDiverge", myChi2CutDiverge=99999999);
    //declareProperty("MCparID", myMCparID=1);//1: e, 2: mu, 3: pi, 4: K, 5: p
    declareProperty("IniHelix", myIniHelix=1);
    declareProperty("nBinTanl",     myNBinTanl= 100);
    declareProperty("nBinDz",       myNBinDz  = 100);
    declareProperty("tanlRange",    myTanlRange = 3.0);
    declareProperty("dzRange",      myDzRange = 50);
    declareProperty("useWireCrossPoint",       myUseWireCrossPoint = false);
    declareProperty("wireCrossPointAvgMode",   myWireCrossPointAvgMode = false);
    declareProperty("WireCrossPointPersistent",myWireCrossPointPersistent = false);
    declareProperty("IdealTracking",myRunIdealTracking= false);
    declareProperty("UseInTrkerHits",myUseInTrkerHits = true);
    declareProperty("TrkFollow",myRunTrkFollow= true);
}

Member Function Documentation

CellAutomaton()

vector< vector< int > > DotsConnection::CellAutomaton	(	vector< int > &	v_cellIdx,
		vector< cell > &	v_cell,
		vector< edge > &	v_edge )

Definition at line 7036 of file DotsConnection.cxx.

{
    bool debug=false; debug=true;
 
    vector<int> v_rootCellIdx;
 
    vector<vector<int> > v_path;
    if(v_cellIdx.size()==0) return v_path;
 
    // --- initialize
    for(vector<int>::iterator it_cellIdx=v_cellIdx.begin(); it_cellIdx!=v_cellIdx.end(); it_cellIdx++) {
        cell& aCell=v_cell.at(*it_cellIdx);
        aCell.state=-1;
        aCell.color=0;
        if(aCell.set_inEdge.empty()) v_rootCellIdx.push_back(*it_cellIdx);
    }
    myNCycle=0;
    myCellMaxState.clear();
 
    // --- state calculation with depth-first-visit
    for(vector<int>::iterator it_cellIdx=v_cellIdx.begin(); it_cellIdx!=v_cellIdx.end(); it_cellIdx++) {
        cell& aCell=v_cell.at(*it_cellIdx);
        //if(aCell.color==2) continue;
        if(aCell.color==0) {
            if(debug) cout<<setw(6)<<""<<"== start state calculation with cell "<<aCell.id<<endl;
            stateAutomaton(aCell, v_cell, v_edge);
        }
    }
    if(debug) cout<<setw(6)<<""<<"The graph has a cycle? "<<myNCycle<<endl;
    int maxState = 0;
    if(myCellMaxState.size()) {
        maxState = v_cell.at(myCellMaxState.at(0)).state;
        if(debug) cout<<setw(6)<<""<<"maximum state: "<<v_cell.at(myCellMaxState.at(0)).state<<", "<<myCellMaxState.size()<<" cells"<<endl;
    }
 
    // --- form track segments (fit to circles)
    for(vector<int>::iterator it_cell_maxStat=myCellMaxState.begin(); it_cell_maxStat!=myCellMaxState.end(); it_cell_maxStat++) 
        //for(int iRootCell=0; iRootCell<=v_rootCellIdx.size(); iRootCell++)
    {
        if(debug) cout<<setw(6)<<""<<"~~~~> to find a path: "<<endl;
        cell& aRootCell = v_cell.at(*it_cell_maxStat);
        //cell& aRootCell = v_cell.at(v_rootCellIdx.at(iRootCell));
 
        map<int, pair<double,double> > setWireHitPos;
        map<int, int> setWireEntry;
        if(aRootCell.state==1) {
            for(vector<Line>::iterator it_line=aRootCell.vecFittedLine.begin(); it_line!=aRootCell.vecFittedLine.end(); it_line++) {
                double rho=(*it_line).rho;
                double theta=(*it_line).theta;
                vector<double>& vS=(*it_line).S;
                for(int i=0; i<NUMBER_OF_DOTS_PER_CELL; i++) {
                    int wireIdx = aRootCell.dots[i]->id;
                    double s=vS.at(i);
                    double x, y;
                    getXY(rho,theta,s,x,y);
                    map<int, int>::iterator it_entry; // setWireEntry;
                    it_entry = setWireEntry.find(wireIdx);
                    if(it_entry!=setWireEntry.end()) {
                        it_entry->second++;
                        map<int, pair<double,double> >::iterator it_pos=setWireHitPos.find(wireIdx);
                        it_pos->second.first+=x;
                        it_pos->second.second+=y;
                    }
                    else {
                        setWireEntry[wireIdx]=1;
                        pair<double,double> xy(x,y);
                        setWireHitPos[wireIdx]=xy;
                    }
                }
            }
        } 
        else nextCellFinding(aRootCell, v_cell, v_edge, setWireHitPos, setWireEntry);
        cout<<setw(6)<<""<<"a path found with "<<setWireEntry.size()<<" dots: "; //<<endl;
        vector<int> aTrkPath;
        //vector< pair<double, double> > pos_hits;
        //vector<double> phi_LR;
        for(map<int, int>::iterator it_entry=setWireEntry.begin(); it_entry!=setWireEntry.end(); it_entry++) {
            int wireIdx = it_entry->first;
            aTrkPath.push_back(wireIdx);
            int layer = myMdcGeomSvc->Wire(wireIdx)->Layer();
            int wire  = myMdcGeomSvc->Wire(wireIdx)->Cell();
            cout<<"("<<layer<<","<<wire<<") ";
            //map<int, pair<double,double> >::iterator it_pos=setWireHitPos.find(wireIdx);
            //it_pos->second.first/=it_entry->second;
            //it_pos->second.second/=it_entry->second;
            ////cout<<"hit ("<<layer<<","<<wire<<"), filled "<<it_entry->second<<" times, mean x,y="<<it_pos->second.first<<","<<it_pos->second.second<<endl;
            ////pos_hits.push_back(it_pos->second);
            //double phiLR=atan2(it_pos->second.second-0.1*myRLayer[layer]*sin(myWirePhi[layer][wire]), 
            //      it_pos->second.first -0.1*myRLayer[layer]*cos(myWirePhi[layer][wire]));
            //if(myMapMdcDigiDd[wireIdx]<0.15) phiLR=-9999.;
            //phi_LR.push_back(phiLR);
        }
        cout<<endl;
        v_path.push_back(aTrkPath);
        //vector<double> par_Taubin = myDotsHelixFitter.CircleFitByTaubin(pos_hits);
        //vector<double> par_circle = myDotsHelixFitter.circleParConversion(par_Taubin);
        //cout<<"fitted circle: dr="<<par_circle.at(0)<<", phi0="<<par_circle.at(1)<<", kappa="<<par_circle.at(2)<<endl;
        //struct trkCandi oneTrkCandi=getTrkCandi(aTrkPath);
        //int nDropHits;
        //cout<<"the path fitted to circle with initial LR : "<<endl;
        // ---------------------------- trkCandi,   chiCut, useIniHelix,hitSel,setHitFitFlag, layerMin, layerMax, use_additional_points, usePhiAmbi, vPhiAmb
        //nDropHits =   circleFitTaubin(oneTrkCandi, 5,     false,      2,     true,          -1,       50,       false,                 true,       phi_LR);
        //nDropHits =   circleFitTaubin(oneTrkCandi, 5, false, 2, true);
    }
 
    return v_path;
}// CellAutomaton

Referenced by longestPathFromCA().

dotComp()

bool DotsConnection::dotComp ( dot * dot1, dot * dot2 )

inline

Definition at line 215 of file DotsConnection.h.

                                           {
            bool less=false;
            int layer1=myMdcGeomSvc->Wire(dot1->id)->Layer();
            int layer2=myMdcGeomSvc->Wire(dot2->id)->Layer();
            if(layer1<layer2) less=true;
            else if(layer1==layer2) {
                int wCell1=myMdcGeomSvc->Wire(dot1->id)->Cell();
                int wCell2=myMdcGeomSvc->Wire(dot2->id)->Cell();
                if(wCell1<wCell2) {
                    if(fabs(wCell1-wCell2)<0.5*myNWire[layer1]) less=true;
                }
                else if(fabs(wCell1-wCell2)>0.5*myNWire[layer1]) less=true;
            }
            return less;
        };

Referenced by forward().

execute()

StatusCode DotsConnection::execute

(

)

Definition at line 293 of file DotsConnection.cxx.

                                   {
    cout.precision(4);
 
    // Part 1: Get the messaging service, print where you are
    MsgStream log(msgSvc(), name());
    log << MSG::INFO << "DotsConnection execute()" << endreq;
 
    // Part 2: Print out the different levels of messages
    /*log << MSG::DEBUG << "A DEBUG message" << endreq;
      log << MSG::INFO << "An INFO message" << endreq;
      log << MSG::WARNING << "A WARNING message" << endreq;
      log << MSG::ERROR << "An ERROR message" << endreq;
      log << MSG::FATAL << "A FATAL error message" << endreq;
      */
 
    SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
    if (!eventHeader) {
        log << MSG::WARNING << "Could not find Event Header" << endreq;
        return StatusCode::FAILURE;
    }
    myRun = eventHeader->runNumber();
    myEvt = eventHeader->eventNumber();
    if(myDebug) {
        cout<<endl<<"------------------------------ "<<endl;
        cout<<"run:"<<myRun<<"  , event: "<<myEvt<<endl;
    }
    if(myDebugNb==2) {
        cout<<endl<<"------------------------------ "<<endl;
        cout<<"DotsConnection: run:"<<myRun<<"  , event: "<<myEvt<<endl;
    }
 
 
    //if(!(myEvt==240|| myEvt==1374 || myEvt==5084 || myEvt==9015) ) 
    //if(myEvt!=6746) 
    //  return StatusCode::SUCCESS;
    //if(myRun==763) 
    //if(myEvt!=3946) return StatusCode::SUCCESS;// run for specific events
    //if(myEvt!=3946) return StatusCode::SUCCESS;// run for specific events
 
    //testDotsHelixFitterAllHits();
    //testDotsHelixFitterPartHits();
 
    // --- event start time
    getEventStartTime();
 
    // --- register RecMdcTrack/HitCol if not there
    bool bookTrkCol = registerRecMdcTrack();
    if(!bookTrkCol) 
    {
        cout<<"DotsConnection::execute(): failed to register RecMdcTrackCol!"<<endl;
        return StatusCode::FAILURE;
    }
 
    //if(myEvt!=2376) return StatusCode::SUCCESS;// run for specific events
 
    // --- tracking starting with grouping neighboured digis
    if(myRunTrkFollow) {
        myVecTrkCandidates.clear();
        //cout<<"test 3"<<endl;
        //myMdcDigiGroup[503].GetWire()->GetNeighborIDType1();
        //cout<<"test 3 done "<<endl;
        vector<const MdcDigi*> vecDigi = getMdcDigiVec();
        if(vecDigi.size()<2000) {
            // fillMdcDigiBuildNeighbors(vecDigi,1,1); // expired
            fillMdcDigiMap(vecDigi);
            buildMdcDigiNeighbors(myMapMdcDigi, 3, 3);
            findOuterEnds();
            getMdcHitCluster();
            getCgemClusters();
            //matchCgemClusterMdcDigi();
            processTrkCandi();
            //processMdcHitCluter();
            findCgemTrkSeg();
            saveGoodTrkCandi();
        }
        else log << MSG::WARNING << "Too Many MDC digis: "<<vecDigi.size() << endreq;
    }
 
    // --- ideal tracking
    if(myRunIdealTracking) {
        getMcFinalChargedStates();
        //myWatch_tot.Start(myWatch_reset);
        associateDigisToMcParticles();
        //myWatch_reset=kFALSE;
        //myWatch_tot.Stop();
        //cout<<"associateDigisToMcParticles used CPU time: "<<myWatch_tot.CpuTime()<<" s"<<endl;
        //cout<<"      finding "<<myWatch_finding.CpuTime()<<" s"<<endl;
        //cout<<"      fitting "<<myWatch_fitting.CpuTime()<<" s"<<endl;
    }
 
    return StatusCode::SUCCESS;
}

finalize()

StatusCode DotsConnection::finalize

(

)

Definition at line 388 of file DotsConnection.cxx.

                                    {
 
    // Part 1: Get the messaging service, print where you are
    MsgStream log(msgSvc(), name());
    log << MSG::INFO << "DotsConnection finalize()" << endreq;
 
    if(myDebugNb==2) {
        cout<<"DotsConnection finalize(): " << endl;
        //cout<<"    associateDigisToMcParticles used CPU time: "<<myWatch_tot.CpuTime()<<" s"<<endl;
        //cout<<"        finding "<<myWatch_finding.CpuTime()<<" s"<<endl;
        //cout<<"        fitting "<<myWatch_fitting.CpuTime()<<" s"<<endl;
    }
 
    return StatusCode::SUCCESS;
}

forward()

bool DotsConnection::forward	(	cell &	a,
		cell &	b )

Definition at line 6752 of file DotsConnection.cxx.

{
    bool forward=false;
    //set<int>::iterator it_dot_a = a.set_dotIdx.begin();
    //set<int>::iterator it_dot_b = b.set_dotIdx.begin();
    for(int i=0; i<NUMBER_OF_DOTS_PER_CELL; i++) 
    {
        //if(*it_dot_a == *it_dot_b) 
        if(a.dots[i]->id == b.dots[i]->id) 
        {
            //it_dot_a++; it_dot_b++;
            continue;
        }
        //else if(*it_dot_a < *it_dot_b) 
        else if(dotComp(a.dots[i], b.dots[i])) 
        {
            //int layer=myMdcGeomSvc->Wire(*it_dot_a)->Layer();
            //if(layer==myMdcGeomSvc->Wire(*it_dot_b)->Layer()) 
            //{// same layer
            //  int Cell_a=myMdcGeomSvc->Wire(*it_dot_a)->Cell();
            //  if(Cell_a+1!=myMdcGeomSvc->Wire(*it_dot_b)->Cell()) break;// at the boundary
            //}
            forward=true;
            break;
        }
    }
    return forward;
}

Referenced by forward().

getMdcDigiVec() [1/2]

vector< const MdcDigi * > DotsConnection::getMdcDigiVec

(

cell &

aCell

)

Definition at line 6701 of file DotsConnection.cxx.

{
    vector<const MdcDigi*> aMdcDigiVec;
    //set<int>::iterator it = aCell.set_dotIdx.begin();
    //for(; it != aCell.set_dotIdx.end(); it++) 
    for(int i=0; i<NUMBER_OF_DOTS_PER_CELL; i++)
    {
        //aMdcDigiVec.push_back(myMapMdcDigi[*it]);
        aMdcDigiVec.push_back(myMapMdcDigi[aCell.dots[i]->id]);
    }
    return aMdcDigiVec;
}

Referenced by execute(), getMdcDigiVec(), getMdcDigiVec(), makeDotCellEdge(), and makeEdge().

getMdcDigiVec() [2/2]

vector< const MdcDigi * > DotsConnection::getMdcDigiVec

(

edge &

aEdge

)

Definition at line 7024 of file DotsConnection.cxx.

{
    vector<const MdcDigi*> aMdcDigiVec;
    //set<int>::iterator it = aEdge.set_dotIdx.begin(); for(; it != aEdge.set_dotIdx.end(); it++) 
    vector<dot*>::iterator it_dot=aEdge.dots.begin();
    for(; it_dot!=aEdge.dots.end(); it_dot++) 
    {
        aMdcDigiVec.push_back(myMapMdcDigi[(*it_dot)->id]);
    }
    return aMdcDigiVec;
}

initialize()

StatusCode DotsConnection::initialize

(

)

Definition at line 60 of file DotsConnection.cxx.

                                     {
 
    //bool debug=false; //debug=true;
 
    // Part 1: Get the messaging service, print where you are
    MsgStream log(msgSvc(), name());
    log << MSG::INFO << " DotsConnection initialize()" << endreq;
 
    // Part 2: Print out the property values
    // log << MSG::INFO << "  MyInt =    " << m_myInt << endreq;
    
    // --- Get MdcGeomSvc ---
    IMdcGeomSvc* imdcGeomSvc;
    //StatusCode sc1 =Gaudi::svcLocator()->service("MdcGeomSvc", imdcGeomSvc);
    StatusCode sc = service ("MdcGeomSvc", imdcGeomSvc);
    myMdcGeomSvc =  dynamic_cast<MdcGeomSvc*> (imdcGeomSvc);
    if (sc.isFailure()) {
        return( StatusCode::FAILURE);
    }
    int nWire = myMdcGeomSvc->getWireSize();
    int nLayer= myMdcGeomSvc->getLayerSize();
    int nSuperLayer= myMdcGeomSvc->getSuperLayerSize();
    int nSeg= myMdcGeomSvc->getSegmentNo();
    /*cout<<"nWire = "<<nWire
        <<"  nLayer="<<nLayer
        <<"  nSuperLayer="<<nSuperLayer
        <<"  nSeg="<<nSeg
        <<endl;*/
    int nCellTot=0;
    int lastWireFlag=-1;
    int iArea=0;
    for(int i=0; i<nLayer; i++)
    {
        const MdcGeoLayer* aLayer = myMdcGeomSvc->Layer(i);
        myNWire[i]=aLayer->NCell();
        myRLayer[i]=aLayer->Radius();
        //cout<<"layer "<<i<<", "<<aLayer->NCell()<<" cells, slant "<<aLayer->Slant()<<", R="<<aLayer->Radius()<<endl;
 
        double nomShift = aLayer->nomShift();
        if(nomShift>0)      myWireFlag[i]=1;
        else if(nomShift<0) myWireFlag[i]=-1;
        else                myWireFlag[i]=0;
        if(i>0&&myWireFlag[i]!=lastWireFlag) {
            iArea++;
        }
        myAreaLayer[i]=iArea;
        lastWireFlag=myWireFlag[i];
        //cout<<"layer "<<i<<" in area "<<iArea<<endl;
 
        //nCellTot+=aLayer->NCell();
        nCellTot+=myNWire[i];
        for(int j=0; j<myNWire[i]; j++)
        {
            const MdcGeoWire* aWire = myMdcGeomSvc->Wire(i,j);
            //cout<<"        wire "<<j<<", BWpos ="<<aWire->Backward()
            //  <<", "<<aWire->BWirePos()
            //  <<", FWpos ="<<aWire->Forward()
            //  <<", "<<aWire->FWirePos()
            //  <<endl;
            int wireidx = aWire->Id();
            myMdcDigiGroup[wireidx].SetWire(aWire);// set wire pointer
            //myMdcDigiGroupPnt[i][j]=&(myMdcDigiGroup[wireidx]);
            myWirePhi[i][j]=aWire->Forward().phi();
            //cout<<"Mdc layer "<<i<<", wire "<<j<<", phi="<<myWirePhi[i][j]<<endl;
            //int layer = aWire->Layer();
            //int cell  = aWire->Cell();
            //cout<<"wire idx "<<wireidx<<",  WireLayer "<<layer<<", cell "<<cell<<", previous idx:"<<aWire->GetNeighborIDType1()<<endl;
            int iInnerLayer = i-1;
            if(iInnerLayer>=0) {
                for(int k=0; k<myNWire[iInnerLayer]; k++)
                {
                    int k_last = k-1;
                    if(k_last<0) k_last=myNWire[iInnerLayer]-1;
                    double dphi_last = dPhi(myWirePhi[iInnerLayer][k_last], myWirePhi[i][j]);
                    double dphi = dPhi(myWirePhi[iInnerLayer][k], myWirePhi[i][j]);
                    if(dphi_last<0&&dphi>0) {
                        myInnerWire[i][j][0]=k_last;
                        myInnerWire[i][j][1]=k;
                        //cout<<"k_last, k ="<<k_last<<", "<<k<<endl;
                        break;
                    }
                }
            }
        }
    }
 
    //cout<<"test "<<endl;
    //myMdcDigiGroup[503].GetWire()->GetNeighborIDType1();
    //cout<<"test done "<<endl;
    
    for(int i=0; i<nLayer; i++)
    {
        for(int j=0; j<myNWire[i]; j++)
        {
            int iOuterLayer = i+1;
            if(iOuterLayer<nLayer) {
                for(int k=0; k<myNWire[iOuterLayer]; k++)
                {
                    int k_last = k-1;// index of previous wire
                    if(k_last<0) k_last=myNWire[iOuterLayer]-1;
                    double dphi_last = dPhi(myWirePhi[iOuterLayer][k_last], myWirePhi[i][j]);
                    double dphi = dPhi(myWirePhi[iOuterLayer][k], myWirePhi[i][j]);
                    if(dphi_last<0&&dphi>0) {
                        myOuterWire[i][j][0]=k_last;
                        myOuterWire[i][j][1]=k;
                        //cout<<"k_last, k ="<<k_last<<", "<<k<<endl;
                        //cout<<"phi="<<myWirePhi[i][j]<<", outer phi = "<<myWirePhi[iOuterLayer][k_last]<<", "<<myWirePhi[iOuterLayer][k]<<endl;
                        break;
                    }
                }
            }
        }
    }
    cout<<"Total "<<nCellTot<<" cells"<<endl;
    // ----------------------
    
    // --- some specific tests
    int layer=14;
    int wire=50;
    vector<int> innerWires_tem = myMdcGeomSvc->Wire(layer,wire)->GetNeighborIDType3();
    int n_wire = innerWires_tem.size();
    cout<<" Layer "<<layer<<", wire "<<wire<<", inner neighbours: "<<endl;
    for(int i=0; i<n_wire; i++) {
        int wireIdx = innerWires_tem[i];
        layer=myMdcGeomSvc->Wire(wireIdx)->Layer();
        wire=myMdcGeomSvc->Wire(wireIdx)->Cell();
        cout<<"(L"<<layer<<", W"<<wire<<") ";
    }
    cout<<endl;
 
 
    ICgemGeomSvc* ISvc;
    StatusCode Cgem_sc = service ("CgemGeomSvc", ISvc);
    if(Cgem_sc.isSuccess()) myCgemGeomSvc=dynamic_cast<CgemGeomSvc *>(ISvc);
    else cout<<"DotsConnection::initialize(): can not get CgemGeomSvc"<<endl;
 
    // --- Initialize RawDataProviderSvc ---
    IRawDataProviderSvc* irawDataProviderSvc;
    sc = service ("RawDataProviderSvc", irawDataProviderSvc);
    myRawDataProviderSvc = dynamic_cast<RawDataProviderSvc*> (irawDataProviderSvc);
    if ( sc.isFailure() ){
        log << MSG::FATAL << name()<<" Could not load RawDataProviderSvc!" << endreq;
        return StatusCode::FAILURE;
    } 
 
    // --- MdcCalibFunSvc ---
    IMdcCalibFunSvc* imdcCalibSvc;                                                                                                      
    sc = service("MdcCalibFunSvc", imdcCalibSvc);
    if ( sc.isSuccess() ){ 
        myMdcCalibFunSvc = dynamic_cast<MdcCalibFunSvc*>(imdcCalibSvc);
    }
    else {
        cout<<"DotsConnection::initialize(): can not get MdcCalibFunSvc"<<endl;
    }
    // --- initialize DotsHelixFitter ---
    myDotsHelixFitter.initialize();
    myDotsHelixFitter.setChi2Diverge(myChi2CutDiverge);
    myDotsHelixFitter.setModeWireCrossPointPersistent(myWireCrossPointPersistent);
    if(myNtProd&1) 
    {
        NTuplePtr nt_ptr(ntupleSvc(),"TestDotsHelixFitter/trkPar");
        if( nt_ptr ) myNtHelixFitter = nt_ptr;
        else
        {
            myNtHelixFitter = ntupleSvc()->book("TestDotsHelixFitter/trkPar",CLID_ColumnWiseTuple,"trkPar");
            if( myNtHelixFitter ) {
                myNtHelixFitter->addItem       ("run",         myRUN);
                myNtHelixFitter->addItem       ("evt",         myEVT);
                myNtHelixFitter->addItem       ("pid",         myPID);
                myNtHelixFitter->addItem       ("Npar",        myNPar,  0,5);
                myNtHelixFitter->addIndexedItem("HelixMC",     myNPar,  myArrayHelixMC);
                myNtHelixFitter->addIndexedItem("HelixFitted", myNPar,  myArrayHelixFitted);
                myNtHelixFitter->addItem       ("NhitCircle",        myNHitsCircle, 0,100);
                myNtHelixFitter->addIndexedItem("LayerHitsCircle",   myNHitsCircle, myLayerHitsCircle);
                myNtHelixFitter->addIndexedItem("ChiHitsCircle",     myNHitsCircle, myChiHitsCircle);
                myNtHelixFitter->addItem       ("Nhit",        myNHits, 0,100);
                myNtHelixFitter->addIndexedItem("LayerHits",   myNHits, myLayerHits);
                myNtHelixFitter->addIndexedItem("ChiHits",     myNHits, myChiHits);
                myNtHelixFitter->addItem       ("NXhit",       myNXHits, 0,100);
                myNtHelixFitter->addItem       ("NVhit",       myNVHits, 0,100);
                myNtHelixFitter->addItem       ("NXCluster",   myNXCluster, 0,100);
                myNtHelixFitter->addItem       ("NVCluster",   myNVCluster, 0,100);
                myNtHelixFitter->addItem       ("TrkIdBest",        myTrkIdBest);
                myNtHelixFitter->addItem       ("NHitsBestTrk",     myNHitsBestTrk);
                myNtHelixFitter->addItem       ("NSameHitsBestTrk", myNSameHitsBestTrk);
                //cout<<"myNtHelixFitter added !"<<endl;
            }
        }
    }
 
    if(myNtProd&2 || myNtProd&4) {
        NTuplePtr nt_ptr(ntupleSvc(),"TestDotsHelixFitter/lineTrkSegHough");
        if( nt_ptr ) myNtTrkSeg = nt_ptr;
        else
        {
            myNtTrkSeg = ntupleSvc()->book("TestDotsHelixFitter/lineTrkSegHough",CLID_ColumnWiseTuple,"hough");
            if( myNtTrkSeg ) {
                myNtTrkSeg->addItem       ("run",   myRUN);
                myNtTrkSeg->addItem       ("evt",   myEVT);
                myNtTrkSeg->addItem       ("nhits", myNt_nhits,  0,100);
                myNtTrkSeg->addIndexedItem("layer", myNt_nhits,  myNt_layer);
                myNtTrkSeg->addIndexedItem("wire",  myNt_nhits,  myNt_wire);
                myNtTrkSeg->addIndexedItem("xhit",  myNt_nhits,  myNt_Xhit);
                myNtTrkSeg->addIndexedItem("yhit",  myNt_nhits,  myNt_Yhit);
                myNtTrkSeg->addIndexedItem("DDhit", myNt_nhits,  myNt_DDhit);
                myNtTrkSeg->addItem       ("rho",   myNt_rho);
                myNtTrkSeg->addItem       ("theta", myNt_theta);
            }
        }
    }
 
    // --- get MdcUtilitySvc
    /*
       sc = service("MdcUtilitySvc", myMdcUtilitySvc);
       if( sc != StatusCode::SUCCESS ){
       log << MSG::FATAL << "can not use MdcUtilitySvc" << endreq;
       }*/
 
    //myWatch_reset=kTRUE;
 
    cout<<"DotsConnection::initialize() myIniHelix = "<<myIniHelix<<endl;
 
    // --- Hough maps
    myRoughTanlDzMap = TH2D("roughTanlDzMap","roughTanlDzMap",myNBinTanl,-1.*myTanlRange,myTanlRange,myNBinDz,-1.*myDzRange,myDzRange);
 
    myRhoRange=78;// in cm
    myThetaRange=acos(-1); // pi
    myRoughRhoThetaMap=TH2D("roughRhoThetaMap","roughRhoThetaMap",100,0,myThetaRange,100,-1.*myRhoRange,myRhoRange);
 
    return StatusCode::SUCCESS;
}

isNeighbour()

bool DotsConnection::isNeighbour	(	cell &	a,
		cell &	b )

Definition at line 6781 of file DotsConnection.cxx.

{
    bool twoDotsShare=false;
 
    //if(a.dots[0]->id==b.dots[1]->id) {
    //  if(a.dots[1]->id==b.dots[0]->id && a.dots[2]->id!=b.dots[2]->id) twoDotsShare=true;
    //  else if(a.dots[2]->id==b.dots[0]->id && a.dots[1]->id!=b.dots[2]->id) twoDotsShare=true;
    //}
    //else if(a.dots[0]->id==b.dots[2]->id) {
    //  if(a.dots[1]->id==b.dots[0]->id && a.dots[2]->id!=b.dots[1]->id) twoDotsShare=true;
    //  else if(a.dots[2]->id==b.dots[0]->id && a.dots[1]->id!=b.dots[1]->id) twoDotsShare=true;
    //}
 
    if(a.dots[0]->id!=b.dots[2]->id && ( (a.dots[1]->id==b.dots[0]->id&&a.dots[2]->id==b.dots[1]->id )||( a.dots[1]->id==b.dots[1]->id&&a.dots[2]->id==b.dots[0]->id)) ) twoDotsShare=true;
    else if(a.dots[2]->id!=b.dots[0]->id && ( (a.dots[0]->id==b.dots[1]->id&&a.dots[1]->id==b.dots[2]->id )||( a.dots[0]->id==b.dots[2]->id&&a.dots[1]->id==b.dots[1]->id)) ) twoDotsShare=true;
 
    return twoDotsShare;
}

Referenced by makeDotCellEdge().

longestPathFromCA() [1/2]

void DotsConnection::longestPathFromCA	(	set< int > &	setWireIdx,
		vector< vector< int > > &	v_path )

Definition at line 7446 of file DotsConnection.cxx.

{
    vector<int> vecWireIdx(setWireIdx.begin(),setWireIdx.end());
    longestPathFromCA(vecWireIdx, v_path);
}

longestPathFromCA() [2/2]

void DotsConnection::longestPathFromCA	(	vector< int > &	vecWireIdx,
		vector< vector< int > > &	v_path )

Definition at line 7428 of file DotsConnection.cxx.

{
    makeDotCellEdge(vecWireIdx);
 
    v_path = CellAutomaton(myVecGoodCellIdx,myVecCells,myVecEdges);
 
    int n_path = v_path.size();
    //for(int i_path=0; i_path<n_path; i_path++) {
    //  vector<int>& aPath = v_path.at(i_path);
    //}
    if(n_path>0) {
        //path = v_path.at(0);
        //if(n_path>1) 
        cout<<setw(6)<<""<<n_path<<" longest paths found"<<endl;
    }
    else cout<<setw(6)<<""<<"no path found"<<endl;
}

Referenced by longestPathFromCA().

makeDotCellEdge()

void DotsConnection::makeDotCellEdge

(

vector< int > &

vecWireIdx

)

Definition at line 7231 of file DotsConnection.cxx.

{
    bool debug=false; //debug=true;
    myMapDots.clear();
    myVecCells.clear();
    myVecGoodCellIdx.clear(); 
    myVecEdges.clear();
 
    // --- make dots and put into myMapDots
    int nHit=vecWireIdx.size();
    for(int i=0; i<nHit; i++) {
        int wireIdx=vecWireIdx.at(i);
        dot aDot;
        aDot.detector = dot::MDC;
        aDot.digi.mdcDigi=myMapMdcDigi[wireIdx];
        aDot.id=wireIdx;
        aDot.cluId=-1;
        myMapDots[wireIdx]=aDot;
    }
 
    // --- make Cells, fill myVecCells, myVecGoodCellIdx, associate cells to wireMids
    map<int, dot>::iterator it_mapDots = myMapDots.begin();
    for(; it_mapDots!=myMapDots.end(); it_mapDots++)
    {
        int wireIdx=it_mapDots->first;
        int layer= myMdcGeomSvc->Wire(wireIdx)->Layer(); //myWireFlag[layer]
        int Cell = myMdcGeomSvc->Wire(wireIdx)->Cell();
 
        //vector<int> vec_neib = myMdcDigiGroup[wireIdx].GetNeiborHits();// get all neighbours, but causes cell duplication
        vector<int> vec_neib = myMdcDigiGroup[wireIdx].GetNeiborHits(1);// get all neighbours and nbs nb, but causes cell duplication
        for(vector<int>::iterator it_nb1=vec_neib.begin(); it_nb1!=vec_neib.end(); it_nb1++)
        {
            // --- get a neighbour
            int wireIdx1=*it_nb1;
            if(myMapDots.find(wireIdx1)==myMapDots.end()) continue;// good dot
            int layer1= myMdcGeomSvc->Wire(wireIdx1)->Layer();
            if(myWireFlag[layer1]!=myWireFlag[layer]) continue;// same stereo angle
            int Cell1 = myMdcGeomSvc->Wire(wireIdx1)->Cell();
 
            // --- get another neighbour 
            vector<int>::iterator it_nb2=it_nb1; it_nb2++; 
            while(true) {// loops
                if(it_nb2==vec_neib.end()) break; // end of loops
                int wireIdx2=*it_nb2;
                if(myMapDots.find(wireIdx2)==myMapDots.end()) {it_nb2++; continue;}// good dot
                int layer2= myMdcGeomSvc->Wire(wireIdx2)->Layer();
                if(myWireFlag[layer2]!=myWireFlag[layer]) {it_nb2++; continue;}// same stereo angle
                int Cell2 = myMdcGeomSvc->Wire(wireIdx2)->Cell();
                if(debug) cout<<" === try cell (id, L, W): "
                    <<"("<<setw(5)<<wireIdx1<<","<<setw(3)<<layer1<<","<<setw(3)<<Cell1<<") "
                        <<"("<<setw(5)<<wireIdx<<","<<setw(3)<<layer<<","<<setw(3)<<Cell<<") "
                        <<"("<<setw(5)<<wireIdx2<<","<<setw(3)<<layer2<<","<<setw(3)<<Cell2<<") "
                        <<endl;
                // --- create a cell
                cell aCell;
                aCell.dots[0]=&(myMapDots[wireIdx1]);
                aCell.dots[1]=&(myMapDots[wireIdx]);
                aCell.dots[2]=&(myMapDots[wireIdx2]);
                // --- check if the cell has already exist
                int wireMid = aCell.dots[1]->id;// wire index of the middle hit
                int nCell_wireMid = myMapDots[wireMid].vecCellIdx.size();
                bool cellExist=false;
                if(debug) cout<<" wireMid="<<wireMid<<" with "<<nCell_wireMid<<" cells"<<endl;
                for(int iCell=0; iCell<nCell_wireMid; iCell++) {
                    if(debug) cout<<" i-cell "<<iCell<<endl;
                    if(sameCell(aCell, myVecCells.at(myMapDots[wireMid].vecCellIdx[iCell]))) { cellExist=true; break; }
                }
                if(cellExist) { 
                    if(debug) cout<<" the cell exists"<<endl;
                    it_nb2++; continue;
                }// drop a duplicated cell
                else if(debug) cout<<" the cell is new"<<endl;
                // --- a new cell found
                aCell.id=myVecCells.size();
                myMapDots.find(wireMid)->second.vecCellIdx.push_back(myVecCells.size());// associate the cell to wireMid
                myVecCells.push_back(aCell);
                myVecCells.back().idInGood=-1;
                // --- fit the cell
                vector<const MdcDigi*> aVecMdcDigi=getMdcDigiVec(aCell);
                vector<double> vRho; vector<double> vTheta; vector<vector<int> > vvLR; vector<vector<double> > vvS;
                bool hasGdLine = LineFit_TLS(aVecMdcDigi, vRho,vTheta,vvLR,vvS);
                if(hasGdLine) {
                    int n_rho = vRho.size();
                    if(debug) cout<<" the cell is good with N_lines="<<setw(2)<<n_rho<<endl;
                    for(int i_line=0; i_line<n_rho; i_line++) {
                        Line aLine;
                        aLine.rho=vRho[i_line];
                        aLine.theta=vTheta[i_line];
                        aLine.LR=vvLR[i_line];
                        aLine.S=vvS[i_line];
                        if(debug) cout<<" S={"
                            <<setw(5)<<aLine.S[0]<<" "
                                <<setw(5)<<aLine.S[1]<<" "
                                <<setw(5)<<aLine.S[2]<<"}"
                                <<endl;
                        // --- only keep the line fit if the middle hit has a middle S
                        if(aLine.S[0]<aLine.S[1]&&aLine.S[2]<aLine.S[1]) {
                            if(debug) cout<<" line "<<i_line<<" has a min S_middle"<<endl;
                            continue;
                        }
                        if(aLine.S[0]>aLine.S[1]&&aLine.S[2]>aLine.S[1]) {
                            if(debug) cout<<" line "<<i_line<<" has a max S_middle"<<endl;
                            continue;
                        }
                        //if(aLine.S[0]*aLine.S[2]<0) {
                        //  if(debug) cout<<" line "<<i_line<<" has different sign for S"<<endl;
                        //  continue;
                        //}
                        if(fabs(aLine.S[0])<fabs(aLine.S[2])) {
                            aLine.order.push_back(0);
                            aLine.order.push_back(1);
                            aLine.order.push_back(2);
                        }
                        else {
                            aLine.order.push_back(2);
                            aLine.order.push_back(1);
                            aLine.order.push_back(0);
                        }
                        if(debug) cout<<" line "<<i_line<<" hits' order "<<aLine.order[0]<<" "<<aLine.order[1]<<" "<<aLine.order[2]<<endl;
 
                        // --- save the good fitted line to the cell
                        myVecCells.back().vecFittedLine.push_back(aLine);
                    }
                    if(myVecCells.back().vecFittedLine.size()>0) {
                        if(debug) cout<<"Good Cell "<<myVecGoodCellIdx.size()<<endl;
                        cout<<setw(6)<<""<<"# Good Cell "<<myVecCells.back().id<<endl;
                        myVecCells.back().idInGood=myVecGoodCellIdx.size();
                        myMapDots[wireMid].vecGoodCellIdx[0].push_back(myVecGoodCellIdx.size());// associate the good cell to wireMid
                        myVecGoodCellIdx.push_back(myVecCells.size()-1);// good cells
                    }
                }
                it_nb2++;
            }// loop nb2
        }// loop nb1
    }// loop the map of dots to make cells
 
    // --- make edges
    debug=true;
    vector<int>::iterator i1_goodCell = myVecGoodCellIdx.begin();
    vector<int>::iterator i2_goodCell;
    int iGdCell=0;
    for(; i1_goodCell!=myVecGoodCellIdx.end(); i1_goodCell++, iGdCell++) {// ---> loop good cells
        // --- set Cells' id, state, weight
        cell& cell_a = myVecCells.at(*i1_goodCell);
        for(int idot=0; idot<3; idot++) {// loop the two end-dots of the good cell (non-middle)
            int ivc=0;
            vector<int>& vCellIdx=myMapDots[cell_a.dots[idot]->id].vecGoodCellIdx[ivc];
            i2_goodCell=vCellIdx.begin();
            for(; i2_goodCell!=vCellIdx.end(); i2_goodCell++) {// loop another good cells
                int iGdCell2=*i2_goodCell;
                if(iGdCell2<iGdCell) continue;//reduce number of comparisons
                cell& cell_b = myVecCells.at(myVecGoodCellIdx.at(*i2_goodCell));
                if(debug) {
                    cout<<" ===> compare two cells: "<<endl;
                    cout<<setw(6)<<""<<"Good cell "<<iGdCell<<", ";      print(cell_a); 
                    cout<<setw(6)<<""<<"Good cell "<<*i2_goodCell<<", "; print(cell_b); 
                }
                //if( (cell_b.dots[0]->id==cell_a.dots[1]->id&&cell_b.dots[2]->id!=cell_a.dots[2-idot]->id) || (cell_b.dots[2]->id==cell_a.dots[1]->id&&cell_b.dots[0]->id!=cell_a.dots[2-idot]->id) )
                if(isNeighbour(cell_a,cell_b))
                {
                    if(debug) cout<<setw(6)<<""<<"neighbour with two shared dots in the middle"<<endl;
                    cell* cell_bkd;
                    cell* cell_fwd; 
                    int idx_cell_fwd;
                    if(cell_a.dots[1]->id>cell_b.dots[1]->id) 
                    {
                        cell_bkd = &cell_b;
                        cell_fwd = &cell_a;
                        idx_cell_fwd = *i1_goodCell;
                        if(debug) cout<<setw(6)<<""<<"two cells' order reversed"<<endl;
                    }
                    else {
                        cell_bkd = &cell_a;
                        cell_fwd = &cell_b;
                        idx_cell_fwd = myVecGoodCellIdx.at(*i2_goodCell);
                    }
                    if(cell_bkd->set_fwdCell.find(idx_cell_fwd)==cell_bkd->set_fwdCell.end()) {// a new forward cell
                        cell_bkd->set_fwdCell.insert(idx_cell_fwd);
                        if(debug) cout<<setw(6)<<""<<"cell "<<cell_fwd->id<<" tagged as a new forward cell of cell "<<cell_bkd->id<<endl;
                        edge aNewEdge; 
                        if(makeEdge(cell_a,cell_b,aNewEdge)) {
                            cell& cell_first = myVecCells.at(aNewEdge.twoCells.first);
                            cell_first.set_outEdge.insert(myVecEdges.size());
                            cell& cell_second = myVecCells.at(aNewEdge.twoCells.second);
                            cell_second.set_inEdge.insert(myVecEdges.size());
                            aNewEdge.id=myVecEdges.size();
                            myVecEdges.push_back(aNewEdge);
                            if(debug) cout<<setw(6)<<""<<"a good out edge added to cell "<<aNewEdge.twoCells.first<<" and "<<aNewEdge.twoCells.second<<endl;
                        }
                    }
                }
            }
        }// loop the two end-dots of the good cell (non-middle)
    }// <--- loop good cells
 
}// makeDotCellEdge(vector<int>& vecWireIdx)

Referenced by longestPathFromCA().

makeEdge()

bool DotsConnection::makeEdge	(	cell &	a,
		cell &	b,
		edge &	e )

Definition at line 6816 of file DotsConnection.cxx.

{
    // --- get the set of shared dots
    int a_shared[NUMBER_OF_DOTS_PER_CELL]={0,0,0};
    int b_shared[NUMBER_OF_DOTS_PER_CELL]={0,0,0};
    int n_shared=0;
    set<int> comDots;
    vector<int> viDotShared_a;
    vector<int> viDotShared_b;
    //set<int>& dots_a = a.set_dotIdx;
    //set<int>& dots_b = b.set_dotIdx;
    //std::set_intersection(dots_a.begin(), dots_a.end(), dots_b.begin(), dots_b.end(), 
    //      std::inserter(comDots, comDots.begin()));
    for(int i=0; i<NUMBER_OF_DOTS_PER_CELL; i++) {
        int id_a=a.dots[i]->id;
        for(int j=0; j<NUMBER_OF_DOTS_PER_CELL; j++) {
            int id_b=b.dots[j]->id;
            if(id_a==id_b) {
                comDots.insert(id_a);
                viDotShared_a.push_back(i);
                viDotShared_b.push_back(j);
                //e.set_dotIdx.insert(id_a);
                e.dots.push_back(a.dots[i]);// fill the common hits as the 1st, 2ed hits
                a_shared[i]=1;
                b_shared[j]=1;
                n_shared++;
            }
        }
    }
    if(n_shared!=2) {
        cout<<setw(6)<<""<<"cannot make edge due to n_shared="<<n_shared<<endl;
        return false;
    }
 
    // --- get the indices of shared dots
    //set<int>::iterator it_dots_a=dots_a.begin();
    //for(int i=0; it_dots_a!=dots_a.end(); it_dots_a++, i++) if(comDots.find(*it_dots_a)!=comDots.end()) viDotShared_a.push_back(i);
    //for(int i=0; i<NUMBER_OF_DOTS_PER_CELL; i++) if(comDots.find(a.dots[i]->id)!=comDots.end()) viDotShared_a.push_back(i);
 
    //for(int i=0; i<NUMBER_OF_DOTS_PER_CELL; i++) if(comDots.find(b.dots[i]->id)!=comDots.end()) viDotShared_b.push_back(i);
    //set<int>::iterator it_dots_b=dots_b.begin();
    //for(int i=0; it_dots_b!=dots_b.end(); it_dots_b++, i++) if(comDots.find(*it_dots_b)!=comDots.end()) viDotShared_b.push_back(i);
 
    // --- compare rho, theta
    double diff_theta_min=9999.;
    double diff_rho_min  =9999.;
    //vector<Line>::iterator itL_a_match = a.vecFittedLine.end();
    //vector<Line>::iterator itL_b_match = b.vecFittedLine.end();
    int iL_a_match(-1), iL_b_match(-1);
    int nReverse_b_match=0;
    vector<Line>::iterator itL_a=a.vecFittedLine.begin();
    for(int iL_a=0; itL_a!=a.vecFittedLine.end(); itL_a++, iL_a++) {
        double rho_a   = (*itL_a).rho;
        double theta_a = (*itL_a).theta;
        vector<Line>::iterator itL_b=b.vecFittedLine.begin();
        for(int iL_b=0; itL_b!=b.vecFittedLine.end(); itL_b++, iL_b++) {
            double rho_b   = (*itL_b).rho;
            double theta_b = (*itL_b).theta;
            int nReverse=0; 
            while( fabs(theta_b-acos(-1.)-theta_a) < fabs(theta_b-theta_a) ) {
                theta_b-=acos(-1.);
                rho_b*=-1;
                nReverse++;
            }
            while( fabs(theta_b+acos(-1.)-theta_a) < fabs(theta_b-theta_a) ) {
                theta_b+=acos(-1.);
                rho_b*=-1;
                nReverse++;
            }
            if(diff_theta_min>fabs(theta_b-theta_a)) {
                diff_theta_min=fabs(theta_b-theta_a);
                diff_rho_min = fabs(rho_b-rho_a);
                //itL_a_match=itL_a; itL_b_match=itL_b;
                iL_a_match=iL_a; iL_b_match=iL_b;
                nReverse_b_match=nReverse;
            }
        }// loop lines of cell b
    }// loop lines of cell a
    if(diff_theta_min>0.1) // || diff_rho_min>1.0) 
    {
        //cout<<"theta or rho unmatched!"<<endl;
        cout<<setw(6)<<""<<"theta unmatched!"<<endl;
        return false;
    }
    else cout<<setw(6)<<""<<"theta matched!"<<endl;
 
    // --- compare LR of shared dots
    vector<int>& LR_a=a.vecFittedLine.at(iL_a_match).LR;
    vector<int>& LR_b=b.vecFittedLine.at(iL_b_match).LR;
    int nShared=viDotShared_a.size();
    cout<<setw(6)<<""<<"shared LR in a: (";
    for(int i=0; i<nShared; i++) cout<<setw(3)<<LR_a.at(viDotShared_a.at(i));
    cout<<")"<<endl;
    int reverse=1; if(nReverse_b_match%2!=0) {reverse=-1; cout<<setw(6)<<""<<"LR reversed for b"<<endl;}
    cout<<setw(6)<<""<<"shared LR in b: (";
    for(int i=0; i<nShared; i++) cout<<setw(3)<<reverse*LR_b.at(viDotShared_b.at(i));
    cout<<")"<<endl;
    bool LR_match=true;
    const double shortDD=0.16;// to be optimized
    for(int i=0; i<nShared; i++) {
        if( reverse*LR_b.at(viDotShared_b.at(i)) != LR_a.at(viDotShared_a.at(i)) ) {
            if(myMapMdcDigiDd[b.dots[viDotShared_b.at(i)]->id]<shortDD) continue;// to be optimized
            LR_match=false;
            break;
        }
    }
    if(LR_match) cout<<setw(6)<<""<<"LR matched"<<endl;
    else {
        cout<<setw(6)<<""<<"LR unmatched"<<endl;
        return false;
    }
 
    // --- fill mapWidxLR and uncommon dots of the edge
    map<int, int> mapWidxLR;
    vector<int>::iterator it =LR_a.begin();
    int idot_a_uncom=-1;
    for(int i=0; i<NUMBER_OF_DOTS_PER_CELL; i++, it++) {
        mapWidxLR[a.dots[i]->id]=*it;
        if(a_shared[i]==0) 
        {
            e.dots.push_back(a.dots[i]);// fill the uncommon hit of cell a as the 3rd hit
            idot_a_uncom=i;
        }
    }
    it =LR_b.begin();
    int idot_b_uncom=-1;
    for(int i=0; i<NUMBER_OF_DOTS_PER_CELL; i++, it++) {
        mapWidxLR[b.dots[i]->id]=(*it)*reverse;
        if(b_shared[i]==0) 
        {
            e.dots.push_back(b.dots[i]);// fill the uncommon hit of cell b as the 4th hit
            idot_b_uncom=i;
        }
    }
 
    // --- get all dots
    //std::set_union(dots_a.begin(), dots_a.end(), dots_b.begin(), dots_b.end(),
    //      std::inserter(e.set_dotIdx, e.set_dotIdx.begin()) );
    //set<int>::iterator it_dot = e.set_dotIdx.begin();
    //for(; it_dot!=e.set_dotIdx.end(); it_dot++) {
    //  int idx = *it_dot;
    //  if(comDots.find(idx)!=comDots.end())    e.vDotStatus.push_back(edge::Common);
    //  else {
    //      if(dots_a.find(idx)!=dots_a.end()) e.vDotStatus.push_back(edge::Prev);
    //      else if(dots_b.find(idx)!=dots_b.end()) e.vDotStatus.push_back(edge::Next);
    //  }
    //}
    vector<int> LR_ini;
    vector<dot*>::iterator it_dot=e.dots.begin();
    for(; it_dot!=e.dots.end(); it_dot++) {
        int idx = (*it_dot)->id;
        int LR_i=0;
        if(myMapMdcDigiDd[idx]>=shortDD) LR_i=mapWidxLR[idx];
        LR_ini.push_back(LR_i);
    }
 
    // --- fit
    //double rho_ini  =0.5*(a.vecFittedLine.at(iL_a_match).rho  +b.vecFittedLine.at(iL_b_match).rho);
    //double theta_ini=0.5*(a.vecFittedLine.at(iL_a_match).theta+b.vecFittedLine.at(iL_b_match).theta);
    double rho_ini   = a.vecFittedLine.at(iL_a_match).rho  ;
    double theta_ini = a.vecFittedLine.at(iL_a_match).theta;
    vector<double> vRho; vector<double> vTheta; vector<vector<int> > vvLR; vector<vector<double> > vvS;
    vector<const MdcDigi*> aVecMdcDigi=getMdcDigiVec(e);
    bool hasGdLine = LineFit_TLS(aVecMdcDigi, vRho,vTheta,vvLR,vvS,true,rho_ini,theta_ini,LR_ini);
    if(!hasGdLine) {
        cout<<setw(6)<<""<<" no good line fit for the edge"<<endl;
        return false;
    }
 
    // --- make an edge
    e.twoCells.first=a.id;
    e.twoCells.second=b.id;
    e.twoLines.first=iL_a_match;
    e.twoLines.second=iL_b_match;
    int n_rho = vRho.size();
    cout<<setw(6)<<""<<"an good edge with "<<n_rho<<" fitted lines"<<endl;
    for(int i_line=0; i_line<n_rho; i_line++) {
        Line aLine;
        aLine.rho=vRho[i_line];
        aLine.theta=vTheta[i_line];
        aLine.LR=vvLR[i_line];
        aLine.S=vvS[i_line];
        // --- save the good fitted line to the edge
        if(makeOrder(aLine)) 
        {
            if(aLine.order.front()>1&&aLine.order.back()>1)
            {
                e.vecFittedLine.push_back(aLine);
                cout<<setw(6)<<""<<"good S order : ";
                for(int i=0; i<aLine.order.size(); i++) cout<<setw(2)<<aLine.order.at(i);
                cout<<endl;
                if(aLine.order.front()==3) {     // uncommon hit of cell b
                    e.twoCells.first=b.id;
                    e.twoCells.second=a.id;
                    e.twoLines.first=iL_b_match;
                    e.twoLines.second=iL_a_match;
                }
                else if(aLine.order.front()!=2) {// uncommon hit of cell a
                    cout<<setw(6)<<""<<" something wrong about the dots' order"<<endl;
                }
            }
        }
    }// loop fitted lines
 
    if(e.vecFittedLine.size()>0) return true;
    else return false;
}

Referenced by makeDotCellEdge().

makeOrder()

bool DotsConnection::makeOrder ( Line & l )

inline

Definition at line 270 of file DotsConnection.h.

        {
            if(l.S.size()==0) return false;
            map<double, int> map_S_i;
            int n_minus(0), n_plus(0);
            vector<double>::iterator it_S=l.S.begin();
            //cout<<"  S= ";
            for(int i=0; it_S!=l.S.end(); it_S++,i++)
            {
                //cout<<setw(5)<<*it_S<<"  ";
                if((*it_S)<0) n_minus++;
                else n_plus++;
                map_S_i[fabs(*it_S)]=i;
            }
            //cout<<endl;
            if(n_minus==l.S.size()||n_plus==l.S.size()) 
            {
                map<double, int>::iterator it=map_S_i.begin();
                for(; it!=map_S_i.end(); it++) l.order.push_back((*it).second);
                return true;
            }
            else return false;
        };

Referenced by makeEdge().

nextCellFinding()

void DotsConnection::nextCellFinding	(	cell &	aCell,
		vector< cell > &	v_cell,
		vector< edge > &	v_edge,
		map< int, pair< double, double > > &	setWireHitPos,
		map< int, int > &	setWireEntry )

Definition at line 7182 of file DotsConnection.cxx.

{
    bool debug=false; debug=true;
 
    set<int>& outEdges = aCell.set_outEdge;
    if(debug) {
        cout<<setw(5)<<""<<"cell "<<aCell.id<<", state="<<aCell.state;
        if(outEdges.size()) cout<<", to fine the next cell";
        else cout<<", no next cell"<<endl;
    }
    for(set<int>::iterator it_edgeIdx=outEdges.begin(); it_edgeIdx!=outEdges.end(); it_edgeIdx++) {
        edge& aEdge = v_edge.at(*it_edgeIdx);
        for(vector<Line>::iterator it_line=aEdge.vecFittedLine.begin(); it_line!=aEdge.vecFittedLine.end(); it_line++) {
            double rho=(*it_line).rho;
            double theta=(*it_line).theta;
            vector<double>& vS=(*it_line).S;
            vector<dot*>::iterator it_dot=aEdge.dots.begin();
            for(int i=0; it_dot!=aEdge.dots.end(); it_dot++, i++) {
                int wireIdx = (*it_dot)->id;
                double s=vS.at(i);
                double x, y;
                getXY(rho,theta,s,x,y);
                map<int, int>::iterator it_entry; // setWireEntry;
                it_entry = setWireEntry.find(wireIdx);
                if(it_entry!=setWireEntry.end()) {
                    it_entry->second++;
                    map<int, pair<double,double> >::iterator it_pos=setWireHitPos.find(wireIdx);
                    it_pos->second.first+=x;
                    it_pos->second.second+=y;
                }
                else {
                    setWireEntry[wireIdx]=1;
                    pair<double,double> xy(x,y);
                    setWireHitPos[wireIdx]=xy;
                }
            }   
        }
 
        cell& aNextCell = v_cell.at(v_edge.at(*it_edgeIdx).twoCells.second);
        //cout<<setw(5)<<""<<"cell "<<aNextCell.id<<endl;
        if((aNextCell.state+1)==aCell.state) {// find the next cell
            if(debug) cout<<", cell "<<aNextCell.id<<" found, state = "<<aNextCell.state<<endl;
            nextCellFinding(aNextCell, v_cell, v_edge, setWireHitPos, setWireEntry);
            break;
        }
    }// loop out edges
}// nextCellFinding

Referenced by CellAutomaton(), and nextCellFinding().

numSharedDots()

int DotsConnection::numSharedDots	(	cell &	a,
		cell &	b )

Definition at line 6714 of file DotsConnection.cxx.

{
    int nShare=0;
 
    //set<int>::iterator it_dot_a = a.set_dotIdx.begin();
    //set<int>::iterator it_dot_b = b.set_dotIdx.begin();
    //set<int>::iterator it_dot_b_lastMatch = it_dot_b;
    //for(; it_dot_a!=a.set_dotIdx.end(); it_dot_a++) {
    //  it_dot_b = it_dot_b_lastMatch;
    //  for(; it_dot_b!=b.set_dotIdx.end(); it_dot_b++) {
    //      if(*it_dot_a == *it_dot_b) {
    //          nShare++;
    //          it_dot_b_lastMatch=it_dot_b;
    //          it_dot_b_lastMatch++;
    //          break;
    //      }
    //  }
    //}
 
    for(int i=0; i<NUMBER_OF_DOTS_PER_CELL; i++) {
        int id_a=a.dots[i]->id;
        for(int j=0; j<NUMBER_OF_DOTS_PER_CELL; j++) {
            int id_b=b.dots[j]->id;
            if(id_a==id_b) {
                nShare++;
                break;
            }
        }
    }
 
    return nShare;
}

Referenced by sameCell().

print()

void DotsConnection::print

(

cell &

aCell

)

Definition at line 6800 of file DotsConnection.cxx.

{
    cout<<" cell "<<setw(5)<<aCell.id;
    //cout<<": (id, L, W) = ";
    cout<<": (L, W, id) = ";
    for(int i=0; i<NUMBER_OF_DOTS_PER_CELL; i++) {
        if(aCell.dots[i]->detector==dot::MDC) {
            int layer = myMdcGeomSvc->Wire(aCell.dots[i]->id)->Layer();
            int Cell  = myMdcGeomSvc->Wire(aCell.dots[i]->id)->Cell();
            cout<<"("<<setw(3)<<layer<<","<<setw(3)<<Cell<<", "<<setw(5)<<aCell.dots[i]->id<<") ";
        }
    }
    cout<<endl;
}

Referenced by makeDotCellEdge(), and stateAutomaton().

sameCell()

bool DotsConnection::sameCell	(	cell &	a,
		cell &	b )

Definition at line 6747 of file DotsConnection.cxx.

{ 
    return numSharedDots(a, b)==NUMBER_OF_DOTS_PER_CELL ;
}

Referenced by makeDotCellEdge().

stateAutomaton()

void DotsConnection::stateAutomaton	(	cell &	aCell,
		vector< cell > &	v_cell,
		vector< edge > &	v_edge )

Definition at line 7144 of file DotsConnection.cxx.

{
    bool debug=false; debug=true;
    aCell.color=1;
    if(debug) {
        cout<<"[ discover"; print(aCell);
        cout<<setw(5)<<""<<"loop next cells:"<<endl;
    }
    set<int>& outEdges = aCell.set_outEdge;
    int state_max=0;
    for(set<int>::iterator it_edgeIdx=outEdges.begin(); it_edgeIdx!=outEdges.end(); it_edgeIdx++) {
        cell& aNextCell = v_cell.at(v_edge.at(*it_edgeIdx).twoCells.second);
        if(debug) cout<<setw(5)<<""<<"cell "<<aNextCell.id<<endl;
        if(aNextCell.color==0) stateAutomaton(aNextCell, v_cell, v_edge);
        else if(aNextCell.color==1) myNCycle++;
        if(aNextCell.state>state_max) state_max=aNextCell.state;
    }
    aCell.state=1;
    aCell.state+=state_max;
    aCell.color=2;
 
    //if(!(myCellMaxState.size()>0&&aCell.state<v_cell.at(myCellMaxState.at(0)).state)) {
    //  myCellMaxState.push_back(aCell.id);
    //}
    if(myCellMaxState.size()==0) myCellMaxState.push_back(aCell.id);
    else if(v_cell.at(myCellMaxState.at(0)).state==aCell.state) myCellMaxState.push_back(aCell.id);
    else if(v_cell.at(myCellMaxState.at(0)).state<aCell.state) {
        myCellMaxState.clear();
        myCellMaxState.push_back(aCell.id);
    }
    if(debug) {
        cout<<"finish cell "<<aCell.id<<" with state "<<aCell.state;//<<endl;
        if(myCellMaxState.size())
            cout<<"    maximum state: "<<v_cell.at(myCellMaxState.at(0)).state<<", "<<myCellMaxState.size()<<" cells"<<endl;
        cout<<"]"<<endl;
    }
}

Referenced by CellAutomaton(), and stateAutomaton().

Member Data Documentation

myCellMaxState

vector<int> DotsConnection::myCellMaxState

Definition at line 314 of file DotsConnection.h.

Referenced by CellAutomaton(), and stateAutomaton().

myMapDots

map<int, dot> DotsConnection::myMapDots

Definition at line 318 of file DotsConnection.h.

Referenced by makeDotCellEdge().

myNCycle

int DotsConnection::myNCycle

Definition at line 313 of file DotsConnection.h.

Referenced by CellAutomaton(), and stateAutomaton().

myVecCells

vector<cell> DotsConnection::myVecCells

Definition at line 319 of file DotsConnection.h.

Referenced by longestPathFromCA(), and makeDotCellEdge().

myVecEdges

vector<edge> DotsConnection::myVecEdges

Definition at line 321 of file DotsConnection.h.

Referenced by longestPathFromCA(), and makeDotCellEdge().

myVecGoodCellIdx

vector<int> DotsConnection::myVecGoodCellIdx

Definition at line 320 of file DotsConnection.h.

Referenced by longestPathFromCA(), and makeDotCellEdge().

---

The documentation for this class was generated from the following files:

• DotsConnection.h
• DotsConnection.cxx