From 0472c034e768a1dd51e45c5ba4a2bbe470df4229 Mon Sep 17 00:00:00 2001 From: Zhang Yao <zhangyao@ihep.ac.cn> Date: Tue, 12 Jan 2021 11:06:09 +0800 Subject: [PATCH] Update RecGenfitAlg for external lib v98 --- Detector/DetDriftChamber/compact/det.xml | 2 +- .../RecGenfitAlg/src/GenfitTrack.cpp | 95 ++++++++++++------- Reconstruction/RecGenfitAlg/src/GenfitTrack.h | 2 +- .../RecGenfitAlg/src/RecGenfitAlgDC.cpp | 11 ++- .../RecGenfitAlg/src/RecGenfitAlgDC.h | 4 +- .../RecGenfitAlg/src/RecGenfitAlgSDT.cpp | 12 ++- .../RecGenfitAlg/src/RecGenfitAlgSDT.h | 6 +- .../src/TruthTracker/TruthTrackerAlg.cpp | 69 +++++++------- 8 files changed, 116 insertions(+), 85 deletions(-) diff --git a/Detector/DetDriftChamber/compact/det.xml b/Detector/DetDriftChamber/compact/det.xml index 7ee0d9ae..6c2ad0d7 100644 --- a/Detector/DetDriftChamber/compact/det.xml +++ b/Detector/DetDriftChamber/compact/det.xml @@ -16,7 +16,7 @@ </includes> <define> - <constant name="world_size" value="30*m"/> + <constant name="world_size" value="2226*m"/> <constant name="world_x" value="world_size"/> <constant name="world_y" value="world_size"/> <constant name="world_z" value="world_size"/> diff --git a/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp b/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp index b66ecd76..37778c62 100644 --- a/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp +++ b/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp @@ -46,6 +46,8 @@ const int GenfitTrack::s_PDG[2][5] bool sortDCHit(edm4hep::ConstSimTrackerHit hit1,edm4hep::ConstSimTrackerHit hit2) { + std::cout<<"hit1"<<hit1<<std::endl; + std::cout<<"hit2"<<hit2<<std::endl; bool isEarly=hit1.getTime()<hit2.getTime(); return isEarly; } @@ -200,14 +202,29 @@ bool GenfitTrack::addSpacePointTrakerHit(edm4hep::ConstTrackerHit& hit, { edm4hep::Vector3d pos=hit.getPosition(); TVectorD p(3); - p[0]=pos.x; - p[1]=pos.y; - p[2]=pos.z; + p[0]=pos.x*dd4hep::mm; + p[1]=pos.y*dd4hep::mm; + p[2]=pos.z*dd4hep::mm; + GenfitMsg::get()<<MSG::DEBUG<<m_name<<" addSpacePointMeasurement "<<hitID + <<"pos "<<p[0]<<" "<<p[1]<<" "<<p[2]<<" cm"<<endmsg; /// New a SpacepointMeasurement - double data[6]; - for(int i=0;i<6;i++) data[i]=hit.getCovMatrix(i); - TMatrixDSym hitCov(3,data,"F");//?FIXME + double cov[6]; + for(int i=0;i<6;i++) { + cov[i]=hit.getCovMatrix(i); + GenfitMsg::get()<<MSG::DEBUG<<"cov "<<cov[i]<<endmsg; + } + TMatrixDSym hitCov(3);//FIXME? + //in SimpleDigi/src/PlanarDigiAlg.cpp + //cov[0] = u_direction[0]; + //cov[1] = u_direction[1]; + //cov[2] = resU; + //cov[3] = v_direction[0]; + //cov[4] = v_direction[1]; + //cov[5] = resV; + hitCov(0,0)=cov[2]*dd4hep::mm*dd4hep::mm; + hitCov(1,1)=cov[2]*dd4hep::mm*dd4hep::mm; + hitCov(2,2)=cov[2]*dd4hep::mm*dd4hep::mm; genfit::SpacepointMeasurement* sMeas = new genfit::SpacepointMeasurement(p,hitCov,(int) hit.getCellID(),hitID, @@ -215,6 +232,7 @@ bool GenfitTrack::addSpacePointTrakerHit(edm4hep::ConstTrackerHit& hit, genfit::TrackPoint* trackPoint = new genfit::TrackPoint(sMeas,m_track); m_track->insertPoint(trackPoint); + GenfitMsg::get()<<MSG::DEBUG<<"end of addSpacePointMeasurement"<<endmsg; return true; } @@ -237,6 +255,7 @@ bool GenfitTrack::addSpacePointMeasurement(const TVectorD& pos, /// New a SpacepointMeasurement TMatrixDSym hitCov(3); + sigma=sigma*dd4hep::mm; hitCov(0,0)=sigma*sigma; hitCov(1,1)=sigma*sigma; hitCov(2,2)=sigma*sigma; @@ -626,14 +645,14 @@ double GenfitTrack::extrapolateToHit( TVector3& poca, TVector3& pocaDir, ///Add space point measurement from edm4hep::Track to genfit track -int GenfitTrack::addSimTrakerHits(const edm4hep::Track& track, +int GenfitTrack::addSimTrackerHits(const edm4hep::Track& track, const edm4hep::MCRecoTrackerAssociationCollection* assoHits, float sigma){ //A TrakerHit collection std::vector<edm4hep::ConstSimTrackerHit> sortedDCTrackHitCol; - GenfitMsg::get()<<MSG::DEBUG<<" VXD " + GenfitMsg::get()<<MSG::DEBUG<<m_name<<" VXD " <<lcio::ILDDetID::VXD<<" SIT " <<lcio::ILDDetID::SIT<<" SET " <<lcio::ILDDetID::SET<<" FTD " @@ -646,41 +665,48 @@ int GenfitTrack::addSimTrakerHits(const edm4hep::Track& track, UTIL::BitField64 encoder(lcio::ILDCellID0::encoder_string); encoder.setValue(hit.getCellID()); int detID=encoder[lcio::ILDCellID0::subdet]; - GenfitMsg::get()<<MSG::DEBUG<<"detID "<<detID<<std::endl; + //GenfitMsg::get()<<MSG::DEBUG<<m_name<<" "<<iHit<<" hit "<<hit + //<<" detID "<<detID<<endmsg; - if (detID==lcio::ILDDetID::VXD || detID==lcio::ILDDetID::SIT || + if(detID==lcio::ILDDetID::VXD || detID==lcio::ILDDetID::SIT || detID==lcio::ILDDetID::SET || detID==lcio::ILDDetID::FTD){ - if(addSpacePointTrakerHit(hit,hitID)){ - GenfitMsg::get()<<MSG::DEBUG<<"add slicon space point"<<endmsg; - hitID++; - }else{ - GenfitMsg::get()<<MSG::ERROR<<"addSpacePointMeasurementOnTrack " - <<detID<<" "<<hit.getCellID()<<" faieled"<<endmsg; + //if(addSpacePointTrakerHit(hit,hitID)){ + // GenfitMsg::get()<<MSG::DEBUG<<"add slicon space point"<<endmsg; + // hitID++; + //}else{ + // GenfitMsg::get()<<MSG::ERROR<<"addSpacePointTrakerHit" + // <<detID<<" "<<hit.getCellID()<<" faieled"<<endmsg; + //} + }else if(detID==7){ + //if(addSpacePointMeasurement(p,sigma,hit.getCellID(),hitID)){ + // GenfitMsg::get()<<MSG::DEBUG<<"add DC space point"<<endmsg; + // hitID++; + //}else{ + // GenfitMsg::get()<<MSG::ERROR<<"addSpacePointMeasurement" + // <<detID<<" faieled" <<endmsg; + //} + float minTime=FLT_MAX; + edm4hep::ConstSimTrackerHit minTimeSimHit; + //Select the SimTrakerHit with least time + for(int iSimHit=0;iSimHit<(int) assoHits->size();iSimHit++){ + if(assoHits->at(iSimHit).getRec()==hit && + assoHits->at(iSimHit).getSim().getTime()<minTime){ + minTimeSimHit=assoHits->at(iSimHit).getSim(); + } } + //std::cout<<"minTimeSimHit "<<minTimeSimHit<<std::endl; + sortedDCTrackHitCol.push_back(minTimeSimHit); }else{ - GenfitMsg::get()<<MSG::INFO - <<"skip hit with detector ID = "<< detID<<" CellID = " - <<hit.getCellID()<<endmsg; - } - float minTime=FLT_MAX; - edm4hep::ConstSimTrackerHit minTimeSimHit; - //Select the SimTrakerHit with least time - for(int iSimHit=0;iSimHit<assoHits->size();iSimHit++){ - if(assoHits->at(iSimHit).getRec()==hit && - assoHits->at(iSimHit).getSim().getTime()<minTime){ - minTimeSimHit=assoHits->at(iSimHit).getSim(); - } + GenfitMsg::get()<<MSG::DEBUG<<" Skip add this hit!"<<endmsg; } - sortedDCTrackHitCol.push_back(minTimeSimHit); - } - GenfitMsg::get()<<MSG::DEBUG<<"addSpacePointTrakerHit on silicon="<<hitID - <<endmsg; + }//end loop over hit on track + GenfitMsg::get()<<MSG::DEBUG<<" addSimTrakerHits trackerHits_size=" <<track.trackerHits_size()<<endmsg; + ///Add DC hits to track //Sort sim DC hits by time std::sort(sortedDCTrackHitCol.begin(),sortedDCTrackHitCol.end(),sortDCHit); - for(auto dCTrackerHit: sortedDCTrackHitCol){ edm4hep::Vector3d pos=dCTrackerHit.getPosition(); TVectorD p(3); @@ -692,10 +718,11 @@ int GenfitTrack::addSimTrakerHits(const edm4hep::Track& track, GenfitMsg::get()<<MSG::DEBUG<<"add DC space point"<<endmsg; hitID++; }else{ - GenfitMsg::get()<<MSG::ERROR<<"addSpacePointMeasurementOnTrack " + GenfitMsg::get()<<MSG::ERROR<<"addSpacePointMeasurement" <<detID<<" faieled" <<endmsg; } } + GenfitMsg::get()<<MSG::DEBUG<<"GenfitTrack nHitAdded "<<hitID<<endmsg; return hitID; } diff --git a/Reconstruction/RecGenfitAlg/src/GenfitTrack.h b/Reconstruction/RecGenfitAlg/src/GenfitTrack.h index 89a7d87a..6ab8db8e 100644 --- a/Reconstruction/RecGenfitAlg/src/GenfitTrack.h +++ b/Reconstruction/RecGenfitAlg/src/GenfitTrack.h @@ -113,7 +113,7 @@ class GenfitTrack { virtual bool addWireMeasurementOnTrack(edm4hep::Track& track, double sigma); ///Add space point from truth to track - int addSimTrakerHits(const edm4hep::Track& track, + int addSimTrackerHits(const edm4hep::Track& track, const edm4hep::MCRecoTrackerAssociationCollection* assoHits, float sigma);// float nSigmaSelection diff --git a/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.cpp b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.cpp index 9c26a1be..8f8445b8 100644 --- a/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.cpp +++ b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.cpp @@ -194,6 +194,7 @@ StatusCode RecGenfitAlgDC::initialize() StatusCode RecGenfitAlgDC::execute() { + StatusCode sc; m_timer=clock(); info()<<" RecGenfitAlgDC in execute()"<<endmsg; @@ -215,7 +216,7 @@ StatusCode RecGenfitAlgDC::execute() const edm4hep::TrackerHitCollection* didiDCHitsCol=nullptr; didiDCHitsCol=m_digiDCHitsCol.get(); if(nullptr==didiDCHitsCol) { - debug()<<"DigiDCHitsCollection not found"<<endmsg; + debug()<<"DigiDCHitCollection not found"<<endmsg; return StatusCode::SUCCESS; } ///retrieve DC Hit Association @@ -242,13 +243,13 @@ StatusCode RecGenfitAlgDC::execute() double eventStartTime=0; if(!genfitTrack->createGenfitTrackFromEDM4HepTrack(pidType,dcTrack, eventStartTime)){ - debug()<<"createGenfitTrack failed!"<<endmsg; + debug()<<"createGenfitTrackFromEDM4HepTrack failed!"<<endmsg; return StatusCode::SUCCESS; } if(m_useTruthHit){ - if(0==genfitTrack->addSimTrakerHits(dcTrack,assoDCHitsCol, + if(0==genfitTrack->addSimTrackerHits(dcTrack,assoDCHitsCol, m_sigmaHit.value())){ - debug()<<"addSpacePointMeasurementOnTrack failed!"<<endmsg; + debug()<<"addSimTrackerHits failed!"<<endmsg; return StatusCode::FAILURE; } }else{ @@ -293,7 +294,7 @@ StatusCode RecGenfitAlgDC::execute() // //system ("pause"); //} - if(m_tuple) m_tuple->write(); + if(m_tuple) sc=m_tuple->write(); return StatusCode::SUCCESS; } diff --git a/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.h b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.h index 6d3cf11f..4619b468 100644 --- a/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.h +++ b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.h @@ -68,7 +68,7 @@ class RecGenfitAlgDC:public GaudiAlgorithm { "EventHeaderCol", Gaudi::DataHandle::Reader, this}; //Drift chamber rec hit and trac DataHandle<edm4hep::TrackerHitCollection> m_digiDCHitsCol{ - "DigiDCHitsCollection", Gaudi::DataHandle::Reader, this}; + "DigiDCHitCollection", Gaudi::DataHandle::Reader, this}; DataHandle<edm4hep::TrackCollection> m_dcTrackCol{ "DCTrackCollection", Gaudi::DataHandle::Reader, this}; //Mc truth @@ -95,7 +95,7 @@ class RecGenfitAlgDC:public GaudiAlgorithm { Gaudi::Property<std::string> m_readout_name{this, "readout", "DriftChamberHitsCollection"}; Gaudi::Property<int> m_debug{this,"debug",false}; - Gaudi::Property<float> m_sigmaHit{this,"sigmaHit",0.011}; + Gaudi::Property<float> m_sigmaHit{this,"sigmaHit",0.11};//mm Gaudi::Property<float> m_nSigmaHit{this,"nSigmaHit",5}; Gaudi::Property<double> m_initCovResPos{this,"initCovResPos",1}; Gaudi::Property<double> m_initCovResMom{this,"initCovResMom",0.1}; diff --git a/Reconstruction/RecGenfitAlg/src/RecGenfitAlgSDT.cpp b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgSDT.cpp index ae52cb69..f55c7891 100644 --- a/Reconstruction/RecGenfitAlg/src/RecGenfitAlgSDT.cpp +++ b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgSDT.cpp @@ -197,6 +197,7 @@ StatusCode RecGenfitAlgSDT::initialize() StatusCode RecGenfitAlgSDT::execute() { + StatusCode sc=StatusCode::SUCCESS; m_timer=clock(); info()<<" RecGenfitAlgSDT in execute()"<<endmsg; @@ -218,7 +219,7 @@ StatusCode RecGenfitAlgSDT::execute() //const edm4hep::TrackerHitCollection* dcDigiCol=nullptr; //dcDigiCol=m_DCDigiCol.get(); //if(nullptr==dcDigiCol) { - // debug()<<"DigiDCHitsCollection not found"<<endmsg; + // debug()<<"DigiDCHitCollection not found"<<endmsg; // return StatusCode::SUCCESS; //} ///retrieve DC Hit Association @@ -249,9 +250,9 @@ StatusCode RecGenfitAlgSDT::execute() debug()<<"createGenfitTrack failed!"<<endmsg; return StatusCode::SUCCESS; } - if(0==genfitTrack->addSimTrakerHits(sdtTrack,dcHitAssociationCol, + if(0==genfitTrack->addSimTrackerHits(sdtTrack,dcHitAssociationCol, m_sigmaHit.value())){ - debug()<<"addSpacePointMeasurementOnTrack failed!"<<endmsg; + debug()<<"addSimTrackerHits failed!"<<endmsg; return StatusCode::SUCCESS; } if(m_debug) genfitTrack->printSeed(); @@ -280,6 +281,7 @@ StatusCode RecGenfitAlgSDT::execute() ++m_fitSuccess[pidType]; }//end loop over particle type }//end loop over a track + m_nRecTrack++; if(m_tuple) debugEvent(); @@ -289,7 +291,7 @@ StatusCode RecGenfitAlgSDT::execute() // //system ("pause"); //} - if(m_tuple) m_tuple->write(); + if(m_tuple) sc=m_tuple->write(); return StatusCode::SUCCESS; } @@ -334,7 +336,7 @@ void RecGenfitAlgSDT::debugTrack(int pidType,const GenfitTrack* genfitTrack) TMatrixDSym fittedCov; TLorentzVector fittedPos; TVector3 fittedMom; - int fittedState=genfitTrack->getFittedState(fittedPos,fittedMom,fittedCov); + genfitTrack->getFittedState(fittedPos,fittedMom,fittedCov); HelixClass helix;//mm and GeV float pos[3]={float(fittedPos.X()/dd4hep::mm),float(fittedPos.Y()/dd4hep::mm), float(fittedPos.Z()/dd4hep::mm)}; diff --git a/Reconstruction/RecGenfitAlg/src/RecGenfitAlgSDT.h b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgSDT.h index 32da0800..1a1031e1 100644 --- a/Reconstruction/RecGenfitAlg/src/RecGenfitAlgSDT.h +++ b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgSDT.h @@ -67,7 +67,7 @@ class RecGenfitAlgSDT:public GaudiAlgorithm { "EventHeaderCol", Gaudi::DataHandle::Reader, this}; //Drift chamber rec hit and trac DataHandle<edm4hep::TrackerHitCollection> m_DCDigiCol{ - "DigiDCHitsCollection", Gaudi::DataHandle::Reader, this}; + "DigiDCHitCollection", Gaudi::DataHandle::Reader, this}; DataHandle<edm4hep::SimTrackerHitCollection> m_DCSimHitCol{ "DCSimHitsCollection", Gaudi::DataHandle::Reader, this}; DataHandle<edm4hep::TrackCollection> m_DCTrackCol{ @@ -97,8 +97,8 @@ class RecGenfitAlgSDT:public GaudiAlgorithm { dd4hep::DDSegmentation::BitFieldCoder* m_decoder; Gaudi::Property<std::string> m_readout_name{this, "readout", "DriftChamberHitsCollection"}; - Gaudi::Property<int> m_debug{this,"debug",false}; - Gaudi::Property<float> m_sigmaHit{this,"sigmaHit",0.011}; + Gaudi::Property<int> m_debug{this,"debug",0}; + Gaudi::Property<float> m_sigmaHit{this,"sigmaHit",0.11};//mm Gaudi::Property<float> m_nSigmaHit{this,"nSigmaHit",5}; Gaudi::Property<double> m_initCovResPos{this,"initCovResPos",1}; Gaudi::Property<double> m_initCovResMom{this,"initCovResMom",0.1}; diff --git a/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp b/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp index 9e24a601..8053ac44 100644 --- a/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp +++ b/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp @@ -117,39 +117,40 @@ StatusCode TruthTrackerAlg::execute() ///Retrieve silicon Track const edm4hep::TrackCollection* siTrackCol=nullptr; - siTrackCol=m_siSubsetTrackCol.get(); - - ///New SDT track - for(auto siTrack:*siTrackCol){ - edm4hep::Track sdtTrack=sdtTrackCol->create(); - edm4hep::TrackState sdtTrackState; - edm4hep::TrackState siTrackStat=siTrack.getTrackStates(0);//FIXME? - sdtTrackState.location=siTrackStat.location; - sdtTrackState.D0=siTrackStat.D0; - sdtTrackState.phi=siTrackStat.phi; - sdtTrackState.omega=siTrackStat.omega; - sdtTrackState.Z0=siTrackStat.Z0; - sdtTrackState.tanLambda=siTrackStat.tanLambda; - sdtTrackState.referencePoint=siTrackStat.referencePoint; - for(int k=0;k<15;k++){ - sdtTrackState.covMatrix[k]=siTrackStat.covMatrix[k]; - } - sdtTrack.addToTrackStates(sdtTrackState); - sdtTrack.setType(siTrack.getType()); - sdtTrack.setChi2(siTrack.getChi2()); - sdtTrack.setNdf(siTrack.getNdf()); - sdtTrack.setDEdx(siTrack.getDEdx()); - sdtTrack.setDEdxError(siTrack.getDEdxError()); - sdtTrack.setRadiusOfInnermostHit(siTrack.getRadiusOfInnermostHit()); - debug()<<"siTrack trackerHits_size="<<siTrack.trackerHits_size()<<endmsg; - for(unsigned int iSiTackerHit=0;iSiTackerHit<siTrack.trackerHits_size(); - iSiTackerHit++){ - sdtTrack.addToTrackerHits(siTrack.getTrackerHits(iSiTackerHit)); - } - //TODO tracks - for(auto digiDC:*digiDCHitsCol){ - //if(Sim->MCParti!=current) continue;//TODO - sdtTrack.addToTrackerHits(digiDC); + if(m_siSubsetTrackCol.exist()) siTrackCol=m_siSubsetTrackCol.get(); + if(nullptr!=siTrackCol) { + ///New SDT track + for(auto siTrack:*siTrackCol){ + edm4hep::Track sdtTrack=sdtTrackCol->create(); + edm4hep::TrackState sdtTrackState; + edm4hep::TrackState siTrackStat=siTrack.getTrackStates(0);//FIXME? + sdtTrackState.location=siTrackStat.location; + sdtTrackState.D0=siTrackStat.D0; + sdtTrackState.phi=siTrackStat.phi; + sdtTrackState.omega=siTrackStat.omega; + sdtTrackState.Z0=siTrackStat.Z0; + sdtTrackState.tanLambda=siTrackStat.tanLambda; + sdtTrackState.referencePoint=siTrackStat.referencePoint; + for(int k=0;k<15;k++){ + sdtTrackState.covMatrix[k]=siTrackStat.covMatrix[k]; + } + sdtTrack.addToTrackStates(sdtTrackState); + sdtTrack.setType(siTrack.getType()); + sdtTrack.setChi2(siTrack.getChi2()); + sdtTrack.setNdf(siTrack.getNdf()); + sdtTrack.setDEdx(siTrack.getDEdx()); + sdtTrack.setDEdxError(siTrack.getDEdxError()); + sdtTrack.setRadiusOfInnermostHit(siTrack.getRadiusOfInnermostHit()); + debug()<<"siTrack trackerHits_size="<<siTrack.trackerHits_size()<<endmsg; + for(unsigned int iSiTackerHit=0;iSiTackerHit<siTrack.trackerHits_size(); + iSiTackerHit++){ + sdtTrack.addToTrackerHits(siTrack.getTrackerHits(iSiTackerHit)); + } + //TODO tracks + for(auto digiDC:*digiDCHitsCol){ + //if(Sim->MCParti!=current) continue;//TODO + sdtTrack.addToTrackerHits(digiDC); + } } } @@ -211,7 +212,7 @@ StatusCode TruthTrackerAlg::execute() //dcTrack.setDEdx();//TODO //set hits double radiusOfInnermostHit=1e9; - debug()<<digiDCHitsCol->size()<<endmsg; + debug()<<"digiDCHitsCol size"<<digiDCHitsCol->size()<<endmsg; for(auto digiDC : *digiDCHitsCol){ //if(Sim->MCParti!=current) continue;//TODO edm4hep::Vector3d digiPos=digiDC.getPosition(); -- GitLab