diff --git a/Digitisers/DCHDigi/src/DCHDigiAlg.cpp b/Digitisers/DCHDigi/src/DCHDigiAlg.cpp
index 41f70a4b9708a5d34dd89eb76037d8b076ed482a..d516b50f43fdd115bb155da0be277849c1793fd7 100644
--- a/Digitisers/DCHDigi/src/DCHDigiAlg.cpp
+++ b/Digitisers/DCHDigi/src/DCHDigiAlg.cpp
@@ -1,3 +1,4 @@
+/* -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*- */
#include "DCHDigiAlg.h"
@@ -52,6 +53,8 @@ StatusCode DCHDigiAlg::initialize()
return StatusCode::FAILURE;
}
+ fRandom.SetSeed(105105);//FIXME: set by users
+
if(m_WriteAna){
NTuplePtr nt( ntupleSvc(), "MyTuples/DCH_digi_evt" );
@@ -66,6 +69,10 @@ StatusCode DCHDigiAlg::initialize()
m_tuple->addItem( "simhit_x", m_n_sim, m_simhit_x).ignore();
m_tuple->addItem( "simhit_y", m_n_sim, m_simhit_y).ignore();
m_tuple->addItem( "simhit_z", m_n_sim, m_simhit_z).ignore();
+ m_tuple->addItem( "Simdca" , m_n_sim, m_Simdca ).ignore();
+ m_tuple->addItem( "simhitT", m_n_sim, m_simhitT ).ignore();
+ m_tuple->addItem( "simhitmom", m_n_sim, m_simhitmom ).ignore();
+ m_tuple->addItem( "simPDG", m_n_sim, m_simPDG ).ignore();
m_tuple->addItem( "chamber" , m_n_digi, m_chamber ).ignore();
m_tuple->addItem( "layer" , m_n_digi, m_layer ).ignore();
m_tuple->addItem( "cell" , m_n_digi, m_cell ).ignore();
@@ -83,6 +90,7 @@ StatusCode DCHDigiAlg::initialize()
m_tuple->addItem( "poca_y" , m_n_digi, m_poca_y ).ignore();
m_tuple->addItem( "hit_dE" , m_n_digi,m_hit_dE ).ignore();
m_tuple->addItem( "hit_dE_dx", m_n_digi,m_hit_dE_dx ).ignore();
+ m_tuple->addItem( "truthlength", m_n_digi,m_truthlength).ignore();
} else { // did not manage to book the N tuple....
info() << " Cannot book N-tuple:" << long( m_tuple ) << endmsg;
}
@@ -110,6 +118,7 @@ StatusCode DCHDigiAlg::execute()
auto SimHit0 = SimHitCol->at(0);
std::map<unsigned long long, std::vector<decltype(SimHit0)>> id_hits_map;
+
for( int i = 0; i < SimHitCol->size(); i++ )
{
auto SimHit = SimHitCol->at(i);
@@ -119,6 +128,10 @@ StatusCode DCHDigiAlg::execute()
if(sim_hit_mom > m_mom_threshold_high) continue;
if(SimHit.getEDep() <= m_edep_threshold) continue;
+ //Wire efficiency
+ double hitProb = fRandom.Uniform(1.);
+ if(hitProb > m_wireEff) continue;
+
if ( id_hits_map.find(id) != id_hits_map.end()) id_hits_map[id].push_back(SimHit);
else
{
@@ -133,116 +146,119 @@ StatusCode DCHDigiAlg::execute()
}
for(auto iter = id_hits_map.begin(); iter != id_hits_map.end(); iter++)
{
- unsigned long long wcellid = iter->first;
- auto trkHit = Vec->create();
- trkHit.setCellID(wcellid);
- double tot_edep = 0 ;
- double tot_length = 0 ;
- double pos_x = 0 ;
- double pos_y = 0 ;
- double pos_z = 0 ;
- double momx,momy = 0;
- int simhit_size = iter->second.size();
- for(unsigned int i=0; i< simhit_size; i++)
- {
- tot_edep += iter->second.at(i).getEDep();//GeV
- }
- int chamber = m_decoder->get(wcellid, "chamber");
- int layer = m_decoder->get(wcellid, "layer" );
- int cellID = m_decoder->get(wcellid, "cellID" );
- TVector3 Wstart(0,0,0);
- TVector3 Wend (0,0,0);
- m_segmentation->cellposition(wcellid, Wstart, Wend);
- float dd4hep_mm = dd4hep::mm;
- //std::cout<<"dd4hep_mm="<<dd4hep_mm<<std::endl;
- // Wstart =(1/dd4hep_mm)* Wstart;// from DD4HEP cm to mm
- // Wend =(1/dd4hep_mm)* Wend ;
- if(m_debug) std::cout<<"DCHDigi wcellid ="<<wcellid<< ",chamber="<<chamber<<",layer="<<layer<<",cellID="<<cellID<<",s_x="<<Wstart.x()<<",s_y="<<Wstart.y()<<",s_z="<<Wstart.z()<<",E_x="<<Wend.x()<<",E_y="<<Wend.y()<<",E_z="<<Wend.z()<<std::endl;
-
- TVector3 denominator = (Wend-Wstart) ;
- float min_distance = 999 ;
- float min_SM_distance = 999 ;
- float min_D_distance = 999 ;
- float tmp_distance =0;
- float SMdca = 0;
- float distance =0;
- TVector3 hitPosition;
- TVector3 PCA;
- for(unsigned int i=0; i< simhit_size; i++)
- {
- float sim_hit_mom = sqrt( iter->second.at(i).getMomentum()[0]*iter->second.at(i).getMomentum()[0] + iter->second.at(i).getMomentum()[1]*iter->second.at(i).getMomentum()[1] + iter->second.at(i).getMomentum()[2]*iter->second.at(i).getMomentum()[2] );//GeV
- float sim_hit_pt = sqrt( iter->second.at(i).getMomentum()[0]*iter->second.at(i).getMomentum()[0] + iter->second.at(i).getMomentum()[1]*iter->second.at(i).getMomentum()[1] );//GeV
- TVector3 pos(iter->second.at(i).getPosition()[0]*dd4hep_mm, iter->second.at(i).getPosition()[1]*dd4hep_mm, iter->second.at(i).getPosition()[2]*dd4hep_mm);
-
- TVector3 sim_mon(iter->second.at(i).getMomentum()[0],iter->second.at(i).getMomentum()[1],iter->second.at(i).getMomentum()[2]);
- float Steplength = iter->second.at(i).getPathLength();
- TVector3 pos_start = pos - 0.5 * Steplength * sim_mon.Unit();
- TVector3 pos_end = pos + 0.5 * Steplength * sim_mon.Unit();
- tmp_distance = m_segmentation->Distance(wcellid,pos_start,pos_end,hitPosition,PCA);
- tmp_distance = tmp_distance/dd4hep_mm; //mm
-
- // std::cout << " Steplength= " << Steplength << std::endl;
-
- if(tmp_distance < min_distance){
- min_distance = tmp_distance;
- pos_x = hitPosition.x(); //pos.x();
- pos_y = hitPosition.y(); //pos.y();
- pos_z = pos.z();
- momx = iter->second.at(i).getMomentum()[0];
- momy = iter->second.at(i).getMomentum()[1];
- }
- tot_length += iter->second.at(i).getPathLength();//mm
- auto asso = AssoVec->create();
- asso.setRec(trkHit);
- asso.setSim(iter->second.at(i));
- asso.setWeight(iter->second.at(i).getEDep()/tot_edep);
- //std::cout<<" asso setRec setSim "<<trkHit<<" "<<iter->second.at(i)<<std::endl;
-
- if(m_WriteAna && (nullptr!=m_tuple)) { // && min_distance <0.3){
- m_simhit_x[m_n_sim] = pos.x();
- m_simhit_y[m_n_sim] = pos.y();
- m_simhit_z[m_n_sim] = pos.z();
- m_n_sim ++ ;
- }
- }
-
- trkHit.setTime(min_distance*1e3/m_velocity);//m_velocity is um/ns, drift time in ns
- trkHit.setEDep(tot_edep);// GeV
- trkHit.setEdx (tot_edep/tot_length); // GeV/mm
- trkHit.setPosition (edm4hep::Vector3d(pos_x, pos_y, pos_z));//position of closest sim hit
- trkHit.setCovMatrix(std::array<float, 6>{m_res_x, 0, m_res_y, 0, 0, m_res_z});//cov(x,x) , cov(y,x) , cov(y,y) , cov(z,x) , cov(z,y) , cov(z,z) in mm
-
- if(m_WriteAna && (nullptr!=m_tuple)) { // && min_distance <0.3){
- m_chamber [m_n_digi] = chamber;
- m_layer [m_n_digi] = layer ;
- m_cell [m_n_digi] = cellID;
- m_cell_x [m_n_digi] = Wstart.x();
- m_cell_y [m_n_digi] = Wstart.y();
- m_cell1_x [m_n_digi] = Wend.x();
- m_cell1_y [m_n_digi] = Wend.y();
- m_hit_x [m_n_digi] = pos_x;
- m_hit_y [m_n_digi] = pos_y;
- m_hit_z [m_n_digi] = pos_z;
- m_mom_x [m_n_digi] = momx ;
- m_mom_y [m_n_digi] = momy ;
- m_dca [m_n_digi] = min_distance;
- m_poca_x [m_n_digi] = PCA.x();
- m_poca_y [m_n_digi] = PCA.y();
- m_hit_dE [m_n_digi] = trkHit.getEDep();
- m_hit_dE_dx[m_n_digi] = trkHit.getEdx() ;
- m_n_digi ++ ;
- }
- }
+ unsigned long long wcellid = iter->first;
+ auto trkHit = Vec->create();
+ trkHit.setCellID(wcellid);
+ double tot_edep = 0 ;
+ double tot_length = 0 ;
+ double pos_x = 0 ;
+ double pos_y = 0 ;
+ double pos_z = 0 ;
+ double momx,momy = 0;
+ int simhit_size = iter->second.size();
+ for(unsigned int i=0; i< simhit_size; i++)
+ {
+ tot_edep += iter->second.at(i).getEDep();//GeV
+ }
+ int chamber = m_decoder->get(wcellid, "chamber");
+ int layer = m_decoder->get(wcellid, "layer" );
+ int cellID = m_decoder->get(wcellid, "cellID" );
+ TVector3 Wstart(0,0,0);
+ TVector3 Wend (0,0,0);
+ m_segmentation->cellposition(wcellid, Wstart, Wend);
+ float dd4hep_mm = dd4hep::mm;
+ //std::cout<<"dd4hep_mm="<<dd4hep_mm<<std::endl;
+ // Wstart =(1/dd4hep_mm)* Wstart;// from DD4HEP cm to mm
+ // Wend =(1/dd4hep_mm)* Wend ;
+ if(m_debug) std::cout<<"DCHDigi wcellid ="<<wcellid<< ",chamber="<<chamber<<",layer="<<layer<<",cellID="<<cellID<<",s_x="<<Wstart.x()<<",s_y="<<Wstart.y()<<",s_z="<<Wstart.z()<<",E_x="<<Wend.x()<<",E_y="<<Wend.y()<<",E_z="<<Wend.z()<<std::endl;
+
+ TVector3 denominator = (Wend-Wstart) ;
+ float min_distance = 999 ;
+ float sim_distance = 999 ;
+ float tmp_distance =0;
+ float SMdca = 0;
+ float distance =0;
+ TVector3 hitPosition;
+ TVector3 PCA;
+ for(unsigned int i=0; i< simhit_size; i++)
+ {
+ float sim_hit_mom = sqrt( iter->second.at(i).getMomentum()[0]*iter->second.at(i).getMomentum()[0] + iter->second.at(i).getMomentum()[1]*iter->second.at(i).getMomentum()[1] + iter->second.at(i).getMomentum()[2]*iter->second.at(i).getMomentum()[2] );//GeV
+ float sim_hit_pt = sqrt( iter->second.at(i).getMomentum()[0]*iter->second.at(i).getMomentum()[0] + iter->second.at(i).getMomentum()[1]*iter->second.at(i).getMomentum()[1] );//GeV
+ TVector3 pos(iter->second.at(i).getPosition()[0]*dd4hep_mm, iter->second.at(i).getPosition()[1]*dd4hep_mm, iter->second.at(i).getPosition()[2]*dd4hep_mm);
+
+ TVector3 sim_mon(iter->second.at(i).getMomentum()[0],iter->second.at(i).getMomentum()[1],iter->second.at(i).getMomentum()[2]);
+ float Steplength = iter->second.at(i).getPathLength();
+ TVector3 pos_start = pos - 0.5 * Steplength * sim_mon.Unit();
+ TVector3 pos_end = pos + 0.5 * Steplength * sim_mon.Unit();
+ tmp_distance = m_segmentation->Distance(wcellid,pos_start,pos_end,hitPosition,PCA);
+ tmp_distance = tmp_distance/dd4hep_mm; //mm
+ sim_distance = tmp_distance;
+ if(tmp_distance < min_distance){
+ min_distance = tmp_distance;
+ pos_x = hitPosition.x(); //pos.x();
+ pos_y = hitPosition.y(); //pos.y();
+ pos_z = pos.z();
+ momx = iter->second.at(i).getMomentum()[0];
+ momy = iter->second.at(i).getMomentum()[1];
+ }
+ tot_length += iter->second.at(i).getPathLength();//mm
+ auto asso = AssoVec->create();
+ asso.setRec(trkHit);
+ asso.setSim(iter->second.at(i));
+ asso.setWeight(iter->second.at(i).getEDep()/tot_edep);
+ //std::cout<<" asso setRec setSim "<<trkHit<<" "<<iter->second.at(i)<<std::endl;
+
+ if(m_WriteAna && (nullptr!=m_tuple)) {
+ m_simhit_x[m_n_sim] = pos.x();
+ m_simhit_y[m_n_sim] = pos.y();
+ m_simhit_z[m_n_sim] = pos.z();
+ m_Simdca[m_n_sim] = sim_distance;
+ m_simhitT[m_n_sim] = iter->second.at(i).getTime();
+ m_simhitmom[m_n_sim] = sim_hit_mom;
+ m_simPDG[m_n_sim] = iter->second.at(i).getMCParticle().getPDG();
+ m_n_sim ++ ;
+ }
+ }
+
+ trkHit.setTime(min_distance*1e3/m_velocity);//m_velocity is um/ns, drift time in ns
+ trkHit.setEDep(tot_edep);// GeV
+ trkHit.setEdx (tot_edep/tot_length); // GeV/mm
+ trkHit.setPosition (edm4hep::Vector3d(pos_x, pos_y, pos_z));//position of closest sim hit
+ trkHit.setCovMatrix(std::array<float, 6>{m_res_x, 0, m_res_y, 0, 0, m_res_z});//cov(x,x) , cov(y,x) , cov(y,y) , cov(z,x) , cov(z,y) , cov(z,z) in mm
+
+ if(m_WriteAna && (nullptr!=m_tuple)) {
+ m_chamber [m_n_digi] = chamber;
+ m_layer [m_n_digi] = layer ;
+ m_cell [m_n_digi] = cellID;
+ m_cell_x [m_n_digi] = Wstart.x();
+ m_cell_y [m_n_digi] = Wstart.y();
+ m_cell1_x [m_n_digi] = Wend.x();
+ m_cell1_y [m_n_digi] = Wend.y();
+ m_hit_x [m_n_digi] = pos_x;
+ m_hit_y [m_n_digi] = pos_y;
+ m_hit_z [m_n_digi] = pos_z;
+ m_mom_x [m_n_digi] = momx ;
+ m_mom_y [m_n_digi] = momy ;
+ m_dca [m_n_digi] = min_distance;
+ m_poca_x [m_n_digi] = PCA.x();
+ m_poca_y [m_n_digi] = PCA.y();
+ m_hit_dE [m_n_digi] = trkHit.getEDep();
+ m_hit_dE_dx[m_n_digi] = trkHit.getEdx() ;
+ m_truthlength[m_n_digi] = tot_length ;
+ m_n_digi ++ ;
+ }
+
+ }
debug()<<"output digi DCHhit size="<< Vec->size() <<endmsg;
_nEvt ++ ;
if(m_WriteAna && (nullptr!=m_tuple)){
StatusCode status = m_tuple->write();
if ( status.isFailure() ) {
- error() << " Cannot fill N-tuple:" << long( m_tuple ) << endmsg;
- return StatusCode::FAILURE;
+ error() << " Cannot fill N-tuple:" << long( m_tuple ) << endmsg;
+ return StatusCode::FAILURE;
}
}
m_end = clock();
@@ -255,6 +271,6 @@ StatusCode DCHDigiAlg::execute()
StatusCode DCHDigiAlg::finalize()
{
- info() << "Processed " << _nEvt << " events " << endmsg;
- return GaudiAlgorithm::finalize();
+ info() << "Processed " << _nEvt << " events " << endmsg;
+ return GaudiAlgorithm::finalize();
}
diff --git a/Digitisers/DCHDigi/src/DCHDigiAlg.h b/Digitisers/DCHDigi/src/DCHDigiAlg.h
index d674dc7753c0bc308d99b3e0ef8e9196f26a4633..4dcf76944852db0580e99a1e9558ffb131fdee5b 100644
--- a/Digitisers/DCHDigi/src/DCHDigiAlg.h
+++ b/Digitisers/DCHDigi/src/DCHDigiAlg.h
@@ -7,6 +7,7 @@
#include "edm4hep/SimTrackerHitCollection.h"
#include "edm4hep/TrackerHitCollection.h"
#include "edm4hep/MCRecoTrackerAssociationCollection.h"
+#include "edm4hep/MCParticleConst.h"
#include <DDRec/DetectorData.h>
#include <DDRec/CellIDPositionConverter.h>
@@ -15,8 +16,7 @@
#include "DetSegmentation/GridDriftChamber.h"
#include "TVector3.h"
-
-
+#include "TRandom3.h"
class DCHDigiAlg : public GaudiAlgorithm
{
@@ -44,6 +44,8 @@ protected:
typedef std::vector<float> FloatVec;
int _nEvt ;
+ TRandom3 fRandom;
+
NTuple::Tuple* m_tuple = nullptr ;
NTuple::Item<int> m_evt;
NTuple::Item<long> m_n_sim;
@@ -59,16 +61,21 @@ protected:
NTuple::Array<float> m_simhit_x ;
NTuple::Array<float> m_simhit_y ;
NTuple::Array<float> m_simhit_z ;
+ NTuple::Array<float> m_simhitT ;
+ NTuple::Array<float> m_simhitmom ;
+ NTuple::Array<int> m_simPDG ;
NTuple::Array<float> m_hit_x ;
NTuple::Array<float> m_hit_y ;
NTuple::Array<float> m_hit_z ;
NTuple::Array<float> m_mom_x ;
NTuple::Array<float> m_mom_y ;
NTuple::Array<float> m_dca ;
+ NTuple::Array<float> m_Simdca ;
NTuple::Array<float> m_poca_x ;
NTuple::Array<float> m_poca_y ;
NTuple::Array<float> m_hit_dE ;
NTuple::Array<float> m_hit_dE_dx ;
+ NTuple::Array<double> m_truthlength ;
clock_t m_start,m_end;
@@ -87,6 +94,7 @@ protected:
Gaudi::Property<float> m_edep_threshold{ this, "edep_threshold", 0};// GeV
Gaudi::Property<bool> m_WriteAna { this, "WriteAna", false};
Gaudi::Property<bool> m_debug{ this, "debug", false};
+ Gaudi::Property<double> m_wireEff{ this, "wireEff", 1.0};
// Input collections
diff --git a/Reconstruction/RecGenfitAlg/CMakeLists.txt b/Reconstruction/RecGenfitAlg/CMakeLists.txt
index 566801d154c2fb36a98a609303ade9b288e6cd01..6136171a0361d29af7c4541b92b1702bebb1b399 100644
--- a/Reconstruction/RecGenfitAlg/CMakeLists.txt
+++ b/Reconstruction/RecGenfitAlg/CMakeLists.txt
@@ -2,11 +2,15 @@
if (GenFit_FOUND)
gaudi_add_module(RecGenfitAlg
- SOURCES src/RecGenfitAlgDC.cpp
+ SOURCES src/RecGenfitAlgSDT.cpp
src/GenfitTrack.cpp
+ src/GenfitHit.cpp
src/GenfitField.cpp
src/GenfitFitter.cpp
src/GenfitMaterialInterface.cpp
+ src/LSFitting.cpp
+ src/WireMeasurementDC.cpp
+ src/PlanarMeasurementSDT.cpp
LINK GearSvc
Gaudi::GaudiAlgLib
Gaudi::GaudiKernel
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitField.cpp b/Reconstruction/RecGenfitAlg/src/GenfitField.cpp
index 744084254143d4480d2cb096970944aa2b6f0a58..e0b8e3548f6be8f7969805790f2201e6aef8246f 100644
--- a/Reconstruction/RecGenfitAlg/src/GenfitField.cpp
+++ b/Reconstruction/RecGenfitAlg/src/GenfitField.cpp
@@ -1,4 +1,5 @@
#include "GenfitField.h"
+#include "GenfitUnit.h"
//External
#include "DD4hep/DD4hepUnits.h"
@@ -39,16 +40,21 @@ GenfitField::get(const double& posX, const double& posY, const double& posZ,
{
/// get field from dd4hepField
const dd4hep::Direction& field=m_dd4hepField.magneticField(
- dd4hep::Position(posX/dd4hep::cm,posY/dd4hep::cm,posZ/dd4hep::cm));
- //m_dd4hepField.magneticField(pos,B);
-
- Bx=field.X()/dd4hep::kilogauss;
- By=field.Y()/dd4hep::kilogauss;
- Bz=field.Z()/dd4hep::kilogauss;
-
-// std::cout<<"GenfitField "
-// <<Form("xyz(%f,%f,%f)cm B(%f,%f,%f)kilogauss",posX,posY,posZ,Bx,By,Bz)
-// <<std::endl;
+ dd4hep::Position(
+ posX/GenfitUnit::cm*dd4hep::cm,
+ posY/GenfitUnit::cm*dd4hep::cm,
+ posZ/GenfitUnit::cm*dd4hep::cm));
+
+ Bx=field.X()/dd4hep::kilogauss*GenfitUnit::kilogauss;
+ By=field.Y()/dd4hep::kilogauss*GenfitUnit::kilogauss;
+ Bz=field.Z()/dd4hep::kilogauss*GenfitUnit::kilogauss;
+
+ //std::cout<<"GenfitField "
+ // <<Form("xyz(%f,%f,%f)cm B(%f,%f,%f)kilogauss",
+ // posX/GenfitUnit::cm*dd4hep::cm,
+ // posY/GenfitUnit::cm*dd4hep::cm,
+ // posZ/GenfitUnit::cm*dd4hep::cm,Bx,By,Bz)
+ // <<std::endl;
}
double GenfitField::getBz(const TVector3& pos) const
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitFitter.cpp b/Reconstruction/RecGenfitAlg/src/GenfitFitter.cpp
index 4850cf4d360e4ed71c705a32d5a18e4cdc8f319c..501b5af337f9c5d113f7eee329c26b6dd1c5c3f0 100644
--- a/Reconstruction/RecGenfitAlg/src/GenfitFitter.cpp
+++ b/Reconstruction/RecGenfitAlg/src/GenfitFitter.cpp
@@ -1,8 +1,15 @@
+
+#undef GENFIT_MY_DEBUG
+//#define GENFIT_MY_DEBUG 1
#include "GenfitFitter.h"
#include "GenfitTrack.h"
#include "GenfitField.h"
#include "GenfitMaterialInterface.h"
+#ifdef GENFIT_MY_DEBUG
+#include "GenfitHist.h"
+#endif
+
//Gaudi
#include "GaudiKernel/StatusCode.h"
#include "GaudiKernel/ISvcLocator.h"
@@ -33,13 +40,15 @@
//STL
#include <iostream>
#include <string>
+#include <string.h>
+//#define GENFIT_MY_DEBUG 1
GenfitFitter::~GenfitFitter(){
delete m_absKalman;
}
-GenfitFitter::GenfitFitter(const char* type, const char* name):
+GenfitFitter::GenfitFitter(const char* type,int debug,const char* name):
m_absKalman(nullptr)
,m_genfitField(nullptr)
,m_geoMaterial(nullptr)
@@ -49,7 +58,7 @@ GenfitFitter::GenfitFitter(const char* type, const char* name):
,m_maxIterations(10)
,m_deltaPval(1e-3)
,m_relChi2Change(0.2)
- ,m_blowUpFactor(1e3)
+ ,m_blowUpFactor(500)
,m_resetOffDiagonals(true)
,m_blowUpMaxVal(1.e6)
,m_multipleMeasurementHandling(genfit::unweightedClosestToPredictionWire)
@@ -66,9 +75,10 @@ GenfitFitter::GenfitFitter(const char* type, const char* name):
,m_noiseBrems(false)
,m_ignoreBoundariesBetweenEqualMaterials(true)
,m_mscModelName("GEANE")
- ,m_debug(0)
+ //,m_debug(0)
,m_hist(0)
{
+ m_debug=debug;
/// Initialize genfit fitter
init();
}
@@ -78,15 +88,20 @@ void GenfitFitter::setField(const GenfitField* field)
if(nullptr==m_genfitField) m_genfitField=field;
}
+
/// Set geometry for material, use geometry from IOADatabase
-void GenfitFitter::setGeoMaterial(const dd4hep::Detector* dd4hepGeo)
+void GenfitFitter::setGeoMaterial(const dd4hep::Detector* dd4hepGeo,
+ double extDistCut, bool skipWireMaterial)
{
if(nullptr==m_geoMaterial){
m_geoMaterial=GenfitMaterialInterface::getInstance(dd4hepGeo);
}
+ m_geoMaterial->setMinSafetyDistanceCut(extDistCut);
+ m_geoMaterial->setSkipWireMaterial(skipWireMaterial);
+ m_geoMaterial->setDebugLvl(m_debug);
}
-/// initialize genfit fitter
+/// initialize genfit fitter, old fitter will be deleted
int GenfitFitter::init(bool deleteOldFitter)
{
if(deleteOldFitter && m_absKalman) delete m_absKalman;
@@ -95,17 +110,21 @@ int GenfitFitter::init(bool deleteOldFitter)
<<m_fitterType<<std::endl;
if (m_fitterType=="DAFRef") {
+ if(m_debug>=2)std::cout<<" m_fitterType==DAFRef "<<std::endl;
m_absKalman = new genfit::DAF(true,getDeltaPval(),
getConvergenceDeltaWeight());
}
else if (m_fitterType=="DAF") {
+ if(m_debug>=2)std::cout<<" m_fitterType==DAF"<<std::endl;
m_absKalman = new genfit::DAF(false,getDeltaPval(),
getConvergenceDeltaWeight());
}
else if (m_fitterType=="KalmanFitter") {
+ if(m_debug>=2)std::cout<<" m_fitterType==KalmanFitter"<<std::endl;
m_absKalman = new genfit::KalmanFitter(getMaxIterations());
}
else if (m_fitterType=="KalmanFitterRefTrack") {
+ if(m_debug>=2)std::cout<<" m_fitterType==KalmanFitterRefTrack"<<std::endl;
m_absKalman = new genfit::KalmanFitterRefTrack(getMaxIterations());
}
else {
@@ -116,6 +135,9 @@ int GenfitFitter::init(bool deleteOldFitter)
}
if(m_debug>=2)std::cout<<"Fitter type is "<<m_fitterType<<std::endl;
m_absKalman->setDebugLvl(m_debug);
+#ifdef GENFIT_MY_DEBUG
+ //m_absKalman->setDebugLvlLocal(m_debugLocal);
+#endif
return 0;
}
@@ -125,6 +147,10 @@ int GenfitFitter::processTrackWithRep(GenfitTrack* track,int repID,bool resort)
{
if(m_debug>=2)std::cout<< "In ProcessTrackWithRep rep "<<repID<<std::endl;
if(getDebug()>2) print("");
+ if(track->getNumPoints()<=0){
+ if(m_debug>=2)std::cout<<"skip track w.o. hit"<<std::endl;
+ return false;
+ }
if(getDebug()>0){
if(m_debug>=2)std::cout<<"Print track seed "<<std::endl;
@@ -147,6 +173,10 @@ int GenfitFitter::processTrackWithRep(GenfitTrack* track,int repID,bool resort)
int GenfitFitter::processTrack(GenfitTrack* track, bool resort)
{
if(m_debug>=2)std::cout<<"In ProcessTrack"<<std::endl;
+ if(track->getNumPoints()<=0){
+ if(m_debug>=2)std::cout<<"skip track w.o. hit"<<std::endl;
+ return false;
+ }
if(getDebug()>2) print("");
/// Do the fitting
@@ -170,142 +200,6 @@ void GenfitFitter::setFitterType(const char* val)
init(oldFitterType==val);
}
-/// Extrapolate track to the cloest point of approach(POCA) to the wire of hit,
-/// return StateOnPlane of this POCA
-/// inputs
-/// pos,mom ... position & momentum at starting point, unit= [mm]&[GeV/c]
-/// (converted to [cm]&[GeV/c] inside this function)
-/// hit ... destination
-/// outputs poca [mm] (converted from [cm] in this function) ,global
-double GenfitFitter::extrapolateToHit( TVector3& poca, TVector3& pocaDir,
- TVector3& pocaOnWire, double& doca, const GenfitTrack* track,
- TVector3 pos, TVector3 mom,
- TVector3 end0,//one end of the hit wire
- TVector3 end1,//the orhter end of the hit wire
- int debug,
- int repID,
- bool stopAtBoundary,
- bool calcJacobianNoise)
-{
-
- return track->extrapolateToHit(poca,pocaDir,pocaOnWire,doca,pos,mom,end0,
- end1,debug,repID,stopAtBoundary,calcJacobianNoise);
-}//end of ExtrapolateToHit
-
-
-/// Extrapolate the track to the cyliner at fixed raidus
-/// position & momentum as starting point
-/// position and momentum at global coordinate in dd4hepUnit
-/// return trackLength in dd4hepUnit
- double
-GenfitFitter::extrapolateToCylinder(TVector3& pos, TVector3& mom,
- GenfitTrack* track, double radius, const TVector3 linePoint,
- const TVector3 lineDirection, int hitID, int repID,
- bool stopAtBoundary, bool calcJacobianNoise)
-{
- double trackLength(1e9*dd4hep::cm);
- if(!track->getFitStatus(repID)->isFitted()) return trackLength;
- try{
- // get track rep
- genfit::AbsTrackRep* rep = track->getRep(repID);
- if(nullptr == rep) {
- if(m_debug>=2)std::cout<<"In ExtrapolateToCylinder rep is null"
- <<std::endl;
- return trackLength*dd4hep::cm;
- }
-
- // get track point with fitter info
- genfit::TrackPoint* tp =
- track->getTrack()->getPointWithFitterInfo(hitID,rep);
- if(nullptr == tp) {
- if(m_debug>=2)std::cout<<
- "In ExtrapolateToCylinder TrackPoint is null"<<std::endl;
- return trackLength*dd4hep::cm;
- }
-
- // get fitted state on plane of this track point
- genfit::KalmanFittedStateOnPlane* state =
- static_cast<genfit::KalmanFitterInfo*>(
- tp->getFitterInfo(rep))->getBackwardUpdate();
-
- if(nullptr == state){
- if(m_debug>=2)std::cout<<"In ExtrapolateToCylinder "
- <<"no KalmanFittedStateOnPlane in backwardUpdate"<<std::endl;
- return trackLength*dd4hep::cm;
- }
- rep->setPosMom(*state, pos*(1/dd4hep::cm), mom*(1/dd4hep::GeV));
-
- /// extrapolate
- trackLength = rep->extrapolateToCylinder(*state,
- radius/dd4hep::cm, linePoint*(1/dd4hep::cm), lineDirection,
- stopAtBoundary, calcJacobianNoise);
- // get pos&mom at extrapolated point on the cylinder
- rep->getPosMom(*state,pos,mom);//FIXME exception exist
- pos = pos*dd4hep::cm;
- mom = mom*dd4hep::GeV;
- } catch(genfit::Exception& e){
- if(m_debug>=3)std::cout
- <<"Exception in GenfitFitter::extrapolateToCylinder "
- << e.what()<<std::endl;
- trackLength = 1e9*dd4hep::cm;
- }
- return trackLength*dd4hep::cm;
-}
-
-double GenfitFitter::extrapolateToPoint(TVector3& pos, TVector3& mom,
- const GenfitTrack* track,
- const TVector3& point,
- int repID,// same with pidType
- bool stopAtBoundary,
- bool calcJacobianNoise) const
-{
- double trackLength(1e9*dd4hep::cm);
- if(!track->getFitStatus(repID)->isFitted()) return trackLength;
- try{
- // get track rep
- genfit::AbsTrackRep* rep = track->getRep(repID);
- if(nullptr == rep) {
- if(m_debug>=2)std::cout<<"In ExtrapolateToPoint rep "
- <<repID<<" not exist!"<<std::endl;
- return trackLength*dd4hep::cm;
- }
-
- /// extrapolate to point
- //genfit::StateOnPlane state(*(&(track->getTrack()->getFittedState(0,rep))));
-
- // get track point with fitter info
- genfit::TrackPoint* tp =
- track->getTrack()->getPointWithFitterInfo(0,rep);
- if(nullptr == tp) {
- if(m_debug>=2)std::cout<<
- "In ExtrapolateToPoint TrackPoint is null"<<std::endl;
- return trackLength*dd4hep::cm;
- }
-
- // get fitted state on plane of this track point
- genfit::KalmanFittedStateOnPlane* state =
- static_cast<genfit::KalmanFitterInfo*>(
- tp->getFitterInfo(rep))->getBackwardUpdate();
-
- if(nullptr == state) {
- if(m_debug>=2)std::cout<<
- "In ExtrapolateToPoint KalmanFittedStateOnPlane is null"<<std::endl;
- return trackLength*dd4hep::cm;
- }
- trackLength = rep->extrapolateToPoint(*state,
- point*(1/dd4hep::cm),stopAtBoundary, calcJacobianNoise);
- rep->getPosMom(*state,pos,mom);//FIXME exception exist
- pos = pos*dd4hep::cm;
- mom = mom*dd4hep::GeV;
- } catch(genfit::Exception& e){
- if(m_debug>=3)std::cout
- <<"Exception in GenfitFitter::extrapolateToPoint"
- << e.what()<<std::endl;
- trackLength = 1e9*dd4hep::cm;
- }
- return trackLength*dd4hep::cm;
-}//end of ExtrapolateToPoint
-
GenfitFitter& GenfitFitter::operator=(
const GenfitFitter& r)
{
@@ -370,7 +264,7 @@ void GenfitFitter::print(const char* name)
<<" m_noiseBrems = " << m_noiseBrems<<std::endl;
if(m_debug>=2)std::cout<<name
<<" m_ignoreBoundariesBetweenEqualMaterials= "
- << m_ignoreBoundariesBetweenEqualMaterials<<std::endl;
+ << m_ignoreBoundariesBetweenEqualMaterials<<std::endl;
if(m_debug>=2)std::cout<<name
<<" m_mscModelName = " << m_mscModelName<<std::endl;
if(m_debug>=2)std::cout<<name
@@ -458,18 +352,27 @@ void GenfitFitter::setBlowUpFactor(double val)
{
m_absKalman->setBlowUpFactor(val);
m_blowUpFactor = val;
+ if (m_fitterType=="DAFRef" || m_fitterType=="DAF") {
+ getDAF()->getKalman()->setBlowUpFactor(m_blowUpFactor);
+ }
}
void GenfitFitter::setResetOffDiagonals(bool val)
{
m_absKalman->setResetOffDiagonals(val);
m_resetOffDiagonals = val;
+ if (m_fitterType=="DAFRef" || m_fitterType=="DAF") {
+ getDAF()->getKalman()->setResetOffDiagonals(m_resetOffDiagonals);
+ }
}
void GenfitFitter::setBlowUpMaxVal(double val)
{
m_absKalman->setBlowUpMaxVal(val);
m_blowUpMaxVal = val;
+ if (m_fitterType=="DAFRef" || m_fitterType=="DAF") {
+ getDAF()->getKalman()->setBlowUpMaxVal(m_blowUpMaxVal);
+ }
}
void GenfitFitter::setMultipleMeasurementHandling(
@@ -558,16 +461,33 @@ void GenfitFitter::setMaterialDebugLvl(unsigned int val)
genfit::MaterialEffects::getInstance()->setDebugLvl(val);
}
-///Set GenfitFitter parameters
void GenfitFitter::setDebug(unsigned int val)
{
- if(m_debug>=2)std::cout<<"set fitter debugLvl "<<val<<std::endl;
- m_absKalman->setDebugLvl(val);
- if(nullptr!=getDAF()) getDAF()->setDebugLvl(val);
- if(val>10) genfit::MaterialEffects::getInstance()->setDebugLvl(val);
m_debug = val;
}
+void GenfitFitter::setDebugGenfit(unsigned int val)
+{
+ //std::cout<<"set fitter debugGenfit "<<val<<std::endl;
+ m_debugGenfit=val;
+ m_absKalman->setDebugLvl(val);
+}
+
+void GenfitFitter::setDebugLocal(unsigned int val)
+{
+ //std::cout<<"set fitter debugLvlLocal "<<val<<std::endl;
+ m_debugLocal=val;
+#ifdef GENFIT_MY_DEBUG
+ ////std::cout<<" GenfitFitter::setDebugLvlLocal "<<val<<std::endl;
+ //m_absKalman->setDebugLvlLocal(val);
+ //if(0==strncmp(m_fitterType.c_str(),"DAF",3)){
+ //std::cout<<" GenfitFitter::setDebugLvlLocal DAF "<<val<<std::endl;
+ // getDAF()->setDebugLvlLocal(val+1);
+ //}
+ //getDAF()->setDebugLvlLocal();
+#endif
+}
+
void GenfitFitter::setHist(unsigned int val) {m_hist = val;}
genfit::DAF* GenfitFitter::getDAF()
@@ -591,21 +511,33 @@ genfit::KalmanFitterRefTrack* GenfitFitter::getKalRef()
}catch(...){
if(m_debug>=3)std::cout
<< "dynamic_cast m_rom AbsFitter to KalmanFitterRefTrack m_ailed!"
- <<std::endl;
+ <<std::endl;
}
return ref;
}
void GenfitFitter::initHist(std::string name)
{
- if(m_debug>=2)std::cout<<"GenfitFitter::initHist "<<name<<std::endl;
- //getDAF()->initHist(name);
+ if(m_debug)std::cout<<"GenfitFitter::initHist "<<name<<std::endl;
+#ifdef GENFIT_MY_DEBUG
+ genfit::GenfitHist::instance()->initHist(name);
+#endif
}
void GenfitFitter::writeHist()
{
- if(m_debug>=2)std::cout<<"GenfitFitter::writeHist "<<std::endl;
- //getDAF()->writeHist();
+ if(m_debug)std::cout<<"GenfitFitter::writeHist "<<std::endl;
+#ifdef GENFIT_MY_DEBUG
+ if(genfit::GenfitHist::instance()->initialized()){
+ genfit::GenfitHist::instance()->writeHist();
+ }
+#endif
}
+void GenfitFitter::SetRunEvent(int event){
+#ifdef GENFIT_MY_DEBUG
+ m_absKalman->SetRunEvent(event);
+#endif
+}
+
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitFitter.h b/Reconstruction/RecGenfitAlg/src/GenfitFitter.h
index 9ca4664456dab2f9484f1b1e9ad2c8d33fc5474a..7b3363e6a27f22030d9acd953fc812521fe851bd 100644
--- a/Reconstruction/RecGenfitAlg/src/GenfitFitter.h
+++ b/Reconstruction/RecGenfitAlg/src/GenfitFitter.h
@@ -18,6 +18,7 @@
#include "AbsKalmanFitter.h"
#include <string>
+#include "MaterialEffects.h"
class GenfitField;
class GenfitMaterialInterface;
@@ -29,6 +30,7 @@ namespace genfit{
class AbsKalmanFitter;
class KalmanFitterRefTrack;
class DAF;
+ class MaterialEffects;
}
namespace dd4hep{
class OverlayedField;
@@ -40,49 +42,22 @@ namespace dd4hep{
class GenfitFitter{
public:
/// Type of fitters are :DAFRef,DAF,KalmanFitter,KalmanFitterRefTrack
- GenfitFitter(const char* type="DAFRef",const char* name="GenfitFitter");
+ GenfitFitter(const char* type="DAFRef",int debug=0,const char* name="GenfitFitter");
virtual ~GenfitFitter();
/// Magnetic field and geometry for material effect in genfit
/// please SET before use !!!!
void setField(const GenfitField* field);
/// please SET before use !!!!
- void setGeoMaterial(const dd4hep::Detector* dd4hepGeo);
+ void setGeoMaterial(const dd4hep::Detector* dd4hepGeo,
+ double extDistCut=1e-4, bool skipWireMaterial=false);
/// Main fitting function
- int processTrack(GenfitTrack* track, bool resort=false);
+ int processTrack(GenfitTrack* track, bool resort=true);
/// fitting with rep
int processTrackWithRep(GenfitTrack* track,int repID=0,
- bool resort=false);
-
- /// Extrapolate the track to the CDC hit
- /// Output: poca pos and dir and poca distance to the hit wire
- /// Input: genfit track, pos and mom, two ends of a wire
- /// pos, and mom are position & momentum at starting point
- double extrapolateToHit(TVector3& poca, TVector3& pocaDir,
- TVector3& pocaOnWire, double& doca, const GenfitTrack* track,
- TVector3 pos, TVector3 mom, TVector3 end0, TVector3 end1,
- int debug=0, int repID=0, bool stopAtBoundary=false,
- bool calcJacobianNoise=false);
-
- /// Extrapolate the track to the cyliner at fixed raidus
- /// Output: pos and mom at fixed radius
- /// Input: genfitTrack, radius of cylinder at center of the origin,
- /// repID, stopAtBoundary and calcAverageState
- double extrapolateToCylinder(TVector3& pos, TVector3& mom,
- GenfitTrack* track, double radius, const TVector3 linePoint,
- const TVector3 lineDirection, int hitID =0, int repID=0,
- bool stopAtBoundary=false, bool calcJacobianNoise=false);
-
- /// Extrapolate the track to the point
- /// Output: pos and mom of POCA point to point
- /// Input: genfitTrack,point,repID,stopAtBoundary and calcAverageState
- /// repID same with pidType
- double extrapolateToPoint(TVector3& pos, TVector3& mom,
- const GenfitTrack* genfitTrack, const TVector3& point,
- int repID=0, bool stopAtBoundary = false,
- bool calcJacobianNoise = false) const;
+ bool resort=true);
/// setters of fitter properties
void setFitterType(const char* val);
@@ -110,10 +85,16 @@ class GenfitFitter{
void setMscModelName(std::string val);
void setMaterialDebugLvl(unsigned int val);
void setDebug(unsigned int val);
+ void setDebugGenfit(unsigned int val);
+ void setDebugLocal(unsigned int val);
void setHist(unsigned int val);
/// getters of fitter properties
std::string getFitterType() const {return m_fitterType;}
+ GenfitMaterialInterface* getGeoMaterial() const {return m_geoMaterial;}
+ genfit::MaterialEffects* getMaterialEffects() const {
+ return genfit::MaterialEffects::getInstance();
+ }
unsigned int getMinIterations() const { return m_minIterations; }
unsigned int getMaxIterations() const { return m_maxIterations; }
double getDeltaPval() const { return m_deltaPval; }
@@ -137,7 +118,10 @@ class GenfitFitter{
{return m_ignoreBoundariesBetweenEqualMaterials;}
std::string getMscModelName(){return m_mscModelName;}
int getDebug() const {return m_debug;}
+ int getDebugGenfit() const {return m_debugGenfit;}
+ int getDebugLocal() const {return m_debugLocal;}
int getHist() const {return m_hist;}
+ void SetRunEvent(int event);
/// Printer
void print(const char* name="");
@@ -194,7 +178,9 @@ class GenfitFitter{
std::string m_mscModelName;
int m_debug; /// debug
- int m_hist; /// hist
+ int m_debugGenfit; /// debug
+ int m_debugLocal; /// debug
+ int m_hist; /// hist
};
#endif
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitMaterialInterface.cpp b/Reconstruction/RecGenfitAlg/src/GenfitMaterialInterface.cpp
index 73a59316e0dbbe52c1b1743eaa6bb6bcc4b4b6c9..031cd7bfa06edc0a805f1ee97f02eb63c1199692 100644
--- a/Reconstruction/RecGenfitAlg/src/GenfitMaterialInterface.cpp
+++ b/Reconstruction/RecGenfitAlg/src/GenfitMaterialInterface.cpp
@@ -1,8 +1,8 @@
///
/// Authors: ZHANG Yao (zhangyao@ihep.ac.cn)
///
-
#include "GenfitMaterialInterface.h"
+#include "GenfitUnit.h"
//Gaudi
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/SmartIF.h"
@@ -21,7 +21,7 @@
#include "TGeoMaterial.h"
#include "TGeoManager.h"
#include "MaterialEffects.h"
-#include "Material.h"
+//#include "Material.h"
//STL
#include <assert.h>
#include <math.h>
@@ -39,6 +39,8 @@ GenfitMaterialInterface::GenfitMaterialInterface(
assert(nullptr!=dd4hepGeo);
m_geoManager=&(dd4hepGeo->manager());
assert(nullptr!=m_geoManager);
+ m_skipWireMaterial=false;
+ debugLvl_=0;
}
GenfitMaterialInterface* GenfitMaterialInterface::getInstance(
@@ -64,7 +66,6 @@ TGeoManager* GenfitMaterialInterface::getGeoManager()
GenfitMaterialInterface::initTrack(double posX, double posY,
double posZ, double dirX, double dirY, double dirZ)
{
- //debugLvl_ = 1;
#ifdef DEBUG_GENFITGEO
std::cout << "GenfitMaterialInterface::initTrack. \n";
std::cout << "Pos "; TVector3(posX, posY, posZ).Print();
@@ -77,11 +78,13 @@ GenfitMaterialInterface::initTrack(double posX, double posY,
// Set the intended direction.
setCurrentDirection(dirX, dirY, dirZ);
+#ifdef DEBUG_GENFITGEO
if (debugLvl_ > 0) {
std::cout << " GenfitMaterialInterface::initTrack at \n";
std::cout << " position: "; TVector3(posX, posY, posZ).Print();
std::cout << " direction: "; TVector3(dirX, dirY, dirZ).Print();
}
+#endif
return result;
}
@@ -91,6 +94,9 @@ genfit::Material GenfitMaterialInterface::getMaterialParameters()
{
TGeoMaterial* mat =
getGeoManager()->GetCurrentVolume()->GetMedium()->GetMaterial();
+#ifdef DEBUG_GENFITGEO
+ getGeoManager()->GetCurrentVolume()->Print();
+#endif
//Scalar density; /// Density in g / cm^3
//Scalar Z; /// Atomic number
//Scalar A; /// Mass number in g / mol
@@ -98,6 +104,17 @@ genfit::Material GenfitMaterialInterface::getMaterialParameters()
//Scalar mEE; /// Mean excitaiton energy in eV
//Material from CEPCSW is NOT follow the dd4hep?FIXME
+ //getGeoManager()->GetCurrentVolume()->Print();
+ if(m_skipWireMaterial){//FIXME
+ if((strncmp(getGeoManager()->GetCurrentVolume()->GetName(),"FieldWire",9) == 0) ||
+ (strncmp(getGeoManager()->GetCurrentVolume()->GetName(),"SenseWire",9) == 0)){
+ //Air: den 0.0012 radlen 30528.8 Z 7.366 A 14.7844
+ std::cout<<" skipWireMateria "<<std::endl;
+ std::cout<<"CurrentVolume "<<getGeoManager()->GetCurrentVolume()->GetName()<<std::endl;
+ return genfit::Material(0.0012,7.366,14.8744,30528.8,MeanExcEnergy_get(mat));
+ }
+ }
+
//std::cout<<__FILE__<<" "<<__LINE__<<" yzhang debug material "<<std::endl;
//mat->Print();
return genfit::Material(mat->GetDensity(), mat->GetZ(), mat->GetA(),
@@ -110,7 +127,6 @@ GenfitMaterialInterface::findNextBoundary(const genfit::RKTrackRep* rep,
double sMax, // signed
bool varField)
{
- //debugLvl_ = 1;
// cm, distance limit beneath which straight-line steps are taken.
const double delta(1.E-2);
const double epsilon(1.E-1); // cm, allowed upper bound on arch
@@ -131,6 +147,8 @@ GenfitMaterialInterface::findNextBoundary(const genfit::RKTrackRep* rep,
findNextBoundary(fabs(sMax) - s);
double safety = getSafeDistance(); // >= 0
double slDist = getStep();
+ if (debugLvl_ > 0)
+ std::cout << " slDist=getStep()= " << slDist<< "; safety=getSafeDistance()=" << safety<< "\n";
// this should not happen, but it happens sometimes.
// The reason are probably overlaps in the geometry.
@@ -175,6 +193,11 @@ GenfitMaterialInterface::findNextBoundary(const genfit::RKTrackRep* rep,
state7 = stateOrig;
rep->RKPropagate(state7, nullptr, SA, stepSign*(s + step), varField);
+ if(debugLvl_){
+ std::cout<<" RKTrackRep at state "<<state7[0]<<" "<<state7[1]
+ <<" "<<state7[2] <<" "<<state7[3] <<" "<<state7[4]
+ <<" "<<state7[5] <<" "<<state7[6]<<std::endl;
+ }
// Straight line distance虏 between extrapolation finish and
// the end of the previously determined safe segment.
double dist2 = (pow(state7[0] - oldState7[0], 2)
@@ -321,12 +344,12 @@ MeanExcEnergy_get(TGeoMaterial* mat) {
double GenfitMaterialInterface::getSafeDistance()
{
- return getGeoManager()->GetSafeDistance();
+ return getGeoManager()->GetSafeDistance();//yzhang debug FIXME
}
double GenfitMaterialInterface::getStep()
{
- return getGeoManager()->GetSafeDistance();
+ return getGeoManager()->GetStep();//yzhang debug FIXME
}
TGeoNode* GenfitMaterialInterface::findNextBoundary(double stepmax,
@@ -340,8 +363,16 @@ bool GenfitMaterialInterface::isSameLocation(double posX, double posY,
{
//std::cout<<" MatInt "<<__LINE__<<" posXYZ*dd4hep::cm "<<posX*dd4hep::cm
//<<" posY "<<posY*dd4hep::cm<<" posZ "<<posZ*dd4hep::cm<<std::endl;
- return getGeoManager()->IsSameLocation(posX,posY,
- posZ,change);
+#ifdef DEBUG_GENFITGEO
+ std::cout<<" MatInt "<<" posXYZ "<<posX<<" "<<posY<<" "<<posZ<<std::endl;
+#endif
+
+ double rootUnitCM=1.;//FIXME
+ //std::cout<<" MatInt "<<" posXYZ cm "<<posX/rootUnitCM<<" "<<posY/rootUnitCM<<" "<<posZ/rootUnitCM<<std::endl;
+ return getGeoManager()->IsSameLocation(
+ posX/GenfitUnit::cm*rootUnitCM,
+ posY/GenfitUnit::cm*rootUnitCM,
+ posZ/GenfitUnit::cm*rootUnitCM,change);
}
void GenfitMaterialInterface::setCurrentDirection(double nx, double ny,
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitMaterialInterface.h b/Reconstruction/RecGenfitAlg/src/GenfitMaterialInterface.h
index 83da9fc1e4cbb372e7c7d7ee35efde3bd5d22646..dbf60261099ef1259a0825aebc16828d41a8ee24 100644
--- a/Reconstruction/RecGenfitAlg/src/GenfitMaterialInterface.h
+++ b/Reconstruction/RecGenfitAlg/src/GenfitMaterialInterface.h
@@ -11,6 +11,7 @@
#define RECGENFITALG_GENFITMATERIALINTERFACE_H
#include "AbsMaterialInterface.h"
+//#include "MaterialProperties.h"
class RKTrackRep;
class TGeoManager;
@@ -33,6 +34,9 @@ class GenfitMaterialInterface : public genfit::AbsMaterialInterface{
//Set Detector pointer of dd4hep
void setDetector(dd4hep::Detector*);
+ //void getMaterialParameters(double& density,double& Z,double& A,
+ // double& radiationLength, double& mEE) {return;}
+ //void getMaterialParameters(genfit::MaterialProperties& parameters) {return;}
/** @brief Initialize the navigator at given position and with given
direction. Returns true if the volume changed.
*/
@@ -54,6 +58,11 @@ class GenfitMaterialInterface : public genfit::AbsMaterialInterface{
bool varField = true) override;
// ClassDef(GenfitMaterialInterface, 1);
+ void setMinSafetyDistanceCut(double safeDistCut=1e-7)
+ {m_safeDistCut=safeDistCut;}
+ void setSkipWireMaterial(bool skipWireMateria)
+ {m_skipWireMaterial=skipWireMateria;}
+ virtual void setDebugLvl(unsigned int lvl = 1) {debugLvl_ = lvl;}
private:
static GenfitMaterialInterface* m_instance;
@@ -66,6 +75,8 @@ class GenfitMaterialInterface : public genfit::AbsMaterialInterface{
bool isSameLocation(double posX, double posY, double posZ,
bool change=false);
void setCurrentDirection(double nx, double ny, double nz);
+ double m_safeDistCut;
+ bool m_skipWireMaterial;
};
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp b/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp
index 7209847b8864166328e6f454662ccd28aaaf559c..5f2c9475fa3e55115137aa1ce26c3bf95092c2cc 100644
--- a/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp
+++ b/Reconstruction/RecGenfitAlg/src/GenfitTrack.cpp
@@ -1,120 +1,142 @@
#include "GenfitTrack.h"
+#include "GenfitHit.h"
#include "GenfitField.h"
//CEPCSW
#include "DataHelper/HelixClass.h"
+#include "DataHelper/TrackHelper.h"
+#include "DataHelper/TrackerHitHelper.h"
+#include "DetInterface/IGeomSvc.h"
+#include "GenfitUnit.h"
#include "UTIL/ILDConf.h"
+#include "WireMeasurementDC.h"
//Externals
#include "DD4hep/DD4hepUnits.h"
-#include "edm4hep/MCParticle.h"
+#include "DD4hep/Detector.h"
+#include "DD4hep/DetElement.h"
+#include "DD4hep/Segmentations.h"
+#include "DDRec/ISurface.h"
+#include "DDRec/SurfaceManager.h"
+#include "DDRec/Vector3D.h"
+#include "DetSegmentation/GridDriftChamber.h"
+#include "edm4hep/MutableReconstructedParticle.h"
+#include "edm4hep/SimTrackerHit.h"
#include "edm4hep/Track.h"
#include "edm4hep/MutableTrack.h"
#include "edm4hep/TrackerHit.h"
-#include "edm4hep/SimTrackerHit.h"
-#include "edm4hep/ReconstructedParticle.h"
-#include "edm4hep/MutableReconstructedParticle.h"
#include "edm4hep/TrackerHitCollection.h"
+#include "edm4hep/MCParticle.h"
#include "edm4hep/MCRecoTrackerAssociationCollection.h"
#include "edm4hep/Vector3d.h"
-#include "DetSegmentation/GridDriftChamber.h"
+#include "GaudiKernel/SmartIF.h"
//genfit
-#include "Track.h"
-#include "MeasuredStateOnPlane.h"
-#include "RKTrackRep.h"
-#include "TrackPoint.h"
-#include "StateOnPlane.h"
-#include "KalmanFitterInfo.h"
-#include "KalmanFittedStateOnPlane.h"
+#include "AbsMeasurement.h"
#include "AbsTrackRep.h"
#include "FitStatus.h"
+#include "KalmanFitterInfo.h"
+#include "KalmanFittedStateOnPlane.h"
+#include "RKTrackRep.h"
+//#include "PlanarMeasurement.h"
+#include "PlanarMeasurementSDT.h"
#include "SpacepointMeasurement.h"
+#include "StateOnPlane.h"
+#include "RectangularFinitePlane.h"
+#include "Track.h"
+#include "TrackPoint.h"
+#include "MeasuredStateOnPlane.h"
#include "WireMeasurementNew.h"
+#include "AbsMeasurement.h"
+#include "TrackPoint.h"
//ROOT
-#include "TRandom.h"
-#include "TVector3.h"
#include "TLorentzVector.h"
#include "TMatrixDSym.h"
+#include "TRandom.h"
+#include "TVector3.h"
//cpp
#include <cfloat>
+#undef GENFIT_MY_DEBUG
+//#define GENFIT_MY_DEBUG 1
+
const int GenfitTrack::s_PDG[2][5]
={{-11,-13,211,321,2212},{11,13,-211,-321,-2212}};
bool
-sortDCHit(edm4hep::SimTrackerHit hit1,edm4hep::SimTrackerHit hit2)
+sortDCSimHit(edm4hep::SimTrackerHit hit1,edm4hep::SimTrackerHit hit2)
{
//std::cout<<"hit1"<<hit1<<std::endl;
//std::cout<<"hit2"<<hit2<<std::endl;
bool isEarly=hit1.getTime()<hit2.getTime();
return isEarly;
}
-
- GenfitTrack::GenfitTrack(const GenfitField* genfitField, const dd4hep::DDSegmentation::GridDriftChamber* seg, const char* name)
-:m_name(name),m_track(nullptr),m_reps(),m_debug(0),
-m_genfitField(genfitField),m_gridDriftChamber(seg)
+ bool
+sortDCDigi(std::pair<double,edm4hep::TrackerHit*> hitPair1,std::pair<double,edm4hep::TrackerHit*> hitPair2)
+{
+ bool isEarly=hitPair1.first<hitPair2.first;
+ return isEarly;
+}
+//bool sortDCDigiLayer(std::pair<int,edm4hep::TrackerHit> hitPair1,std::pair<int,edm4hep::TrackerHit> hitPair2)
+//{
+// bool isEarly=hitPair1.first<hitPair2.first;
+// return isEarly;
+//}
+
+
+GenfitTrack::GenfitTrack(const GenfitField* genfitField,
+ const dd4hep::DDSegmentation::GridDriftChamber* seg,
+ SmartIF<IGeomSvc> geom, const char* name)
+:m_name(name),m_track(nullptr),m_debug(0),m_debugLocal(0),m_geomSvc(geom),
+ m_genfitField(genfitField),m_gridDriftChamber(seg),
+ m_decoderDC(geom->getDecoder("DriftChamberHitsCollection"))
{
-
}
GenfitTrack::~GenfitTrack()
{
+ clearGenfitHitVec();
///Note: track reps and points will be deleted in the destructor of track
///implemented in genfit::Track::Clear()
delete m_track;
}
-/// create a Genfit track from track state, without trackRep
-/// Initialize track with seed states
-/// NO unit conversion here
+void GenfitTrack::clearGenfitHitVec(){
+ for(auto h:m_genfitHitVec){delete h;}
+ m_genfitHitVec.clear();
+ std::vector<GenfitHit*>().swap(m_genfitHitVec);
+}
+
+/// create a Genfit track from track state
+/// Initialize track with seed state and cov
+/// unit conversion here!!!
bool GenfitTrack::createGenfitTrack(int pdgType,int charge,
- TLorentzVector posInit, TVector3 momInit, TMatrixDSym covMInit)
+ TVectorD trackParam, TMatrixDSym covMInit_6)
{
- TVectorD seedState(6);
- TMatrixDSym seedCov(6);
-
- //for(int i = 0; i < 6; ++i) {
- // for(int j = 0; j < 6; ++j) {
- // seedCov(i,j)=covMInit(i,j);
- // }
- //}
- //yzhang FIXME
- //seed position
- for(int i = 0; i < 3; ++i) {
- seedState(i)=posInit[i];
- //yzhang TODO from covMInit to seedCov
- double resolution = 0.1;//*dd4hep::mm/dd4hep::cm;
- seedCov(i,i)=resolution*resolution;
- if(i==2) seedCov(i,i)=0.5*0.5;
- }
- //seed momentum
- for(int i = 3; i < 6; ++i){
- //seedState(i)=momInit[i-3]*(dd4hep::GeV);
- seedState(i)=momInit[i-3];
- //yzhang TODO from covMInit to seedCov
- seedCov(i,i)=0.01;//pow(resolution / sqrt(3),2);
- }
-
- if(nullptr==m_track) m_track=new genfit::Track();
- m_track->setStateSeed(seedState);
- m_track->setCovSeed(seedCov);
-
- /// new a track representation and add to the track
- int chargeId=0;
- charge>0 ? chargeId=0 : chargeId=1;
+ if(m_debug){
+ std::cout<<"createGenfitTrack pdgType "<<pdgType<<" charge "<<charge
+ <<" trackParam "<<std::endl; trackParam.Print();
+ std::cout<<"Cov"<<std::endl; covMInit_6.Print();
+ }
+ ///new a track and set seed state and cov
+ if(nullptr!=m_track) {
+ delete m_track;
+ clearGenfitHitVec();
+ }
+ m_track=new genfit::Track();
+ m_track->setStateSeed(trackParam);
+ m_track->setCovSeed(covMInit_6);
- if(m_debug>=2)std::cout<<m_name<<" CreateGenfitTrack seed pos("
- <<seedState[0]<<" "<<seedState[1]<<" "<<seedState[2]<<")cm ("
- <<seedState[3]<<" "<<seedState[4]<<" "<<seedState[5]<<")GeV charge "
- <<charge<<" pdg "<<s_PDG[chargeId][pdgType]<<std::endl;
- if(m_debug>=2)std::cout<<"seedCov "<<std::endl;
+#ifdef GENFIT_MY_DEBUG
+ //m_track->setDebugLvlLocal(m_debugLocal);
+#endif
- addTrackRep(s_PDG[chargeId][pdgType]);
+ ///new a track representation and add to the track
+ addTrackRep(pdgType,charge);
- if(m_debug>0) seedCov.Print();
+ if(m_debug>0) printSeed();
return true;
}
@@ -122,227 +144,292 @@ bool GenfitTrack::createGenfitTrack(int pdgType,int charge,
bool GenfitTrack::createGenfitTrackFromMCParticle(int pidType,
const edm4hep::MCParticle& mcParticle, double eventStartTime)
{
- ///get track parameters from McParticle
- edm4hep::Vector3d mcPocaPos = mcParticle.getVertex();//mm
- edm4hep::Vector3f mcPocaMom = mcParticle.getMomentum();//GeV
- if(m_debug>=2)std::cout<<"seedPos poca "<< mcPocaPos.x
- <<" "<<mcPocaPos.y<<" "<<mcPocaPos.z<<" mm "<<std::endl;
- if(m_debug>=2)std::cout<<"seedMom poca "<< mcPocaMom.x
- <<" "<<mcPocaMom.y<<" "<<mcPocaMom.z<<" GeV "<<std::endl;
-
- ///Pivot to first layer to avoid correction of beam pipe
- edm4hep::Vector3d firstLayerPos(1e9,1e9,1e9);
- edm4hep::Vector3f firstLayerMom(1e9,1e9,1e9);
- pivotToFirstLayer(mcPocaPos,mcPocaMom,firstLayerPos,firstLayerMom);
-
- //TODO convert unit
- ///Get seed position and momentum
- TLorentzVector seedPos(firstLayerPos.x,firstLayerPos.y,firstLayerPos.z,
- eventStartTime);
- TVector3 seedMom(firstLayerMom.x,firstLayerMom.y,firstLayerMom.z);
- if(m_debug>=2)std::cout<<"seedPos "<< firstLayerPos.x
- <<" "<<firstLayerPos.y<<" "<<firstLayerPos.z<<std::endl;
- if(m_debug>=2)std::cout<<"seedMom "<< firstLayerMom.x
- <<" "<<firstLayerMom.y<<" "<<firstLayerMom.z<<std::endl;
-
///Get error matrix of seed track
- TMatrixDSym covM(5);//FIXME, TODO
-
- ///Create a genfit track with seed
- if(m_debug>=2)std::cout<<"createGenfitTrack " ;
- if(!GenfitTrack::createGenfitTrack(pidType,mcParticle.getCharge(),
- seedPos,seedMom,covM)){
- if(m_debug>=2)std::cout<<"GenfitTrack"
- <<" Error in createGenfitTrack" <<std::endl;
- return false;
- }else{
- if(m_debug>=2)std::cout<<"GenfitTrack "
- <<"createGenfitTrackFromMCParticle track created" <<std::endl;
- }
- return true;
+ TMatrixDSym covMInit_6(6);
+ getSeedCov(covMInit_6);
+ TVectorD seedState(6);
+ getTrackFromMCPartile(mcParticle,seedState,covMInit_6);
+
+ return createGenfitTrack(pidType,mcParticle.getCharge(),seedState,covMInit_6);
}//end of createGenfitTrackFromMCParticle
-///Create a Genfit track with MCParticle, unit conversion here
+///Create a Genfit track with edm4hep Track
bool GenfitTrack::createGenfitTrackFromEDM4HepTrack(int pidType,
- edm4hep::Track track, double eventStartTime)
+ const edm4hep::Track& track, double eventStartTime, bool isUseCovTrack)
{
- //std::cout<<__FILE__<<" "<<__LINE__<<" bz kilogauss "<<m_genfitField->getBz({0.,0.,0.})/dd4hep::kilogauss<<std::endl;
- //std::cout<<__FILE__<<" "<<__LINE__<<" bz tesla "<<m_genfitField->getBz({0.,0.,0.})/dd4hep::tesla<<std::endl;
- //std::cout<<__FILE__<<" "<<__LINE__<<" bz "<<m_genfitField->getBz({0.,0.,0.})<< dd4hep::kilogauss <<" "<<dd4hep::tesla<<std::endl;
- //TODO
- //pivotToFirstLayer(mcPocaPos,mcPocaMom,firstLayerPos,firstLayerMom);
- //Get track parameters
- edm4hep::TrackState trackStat=track.getTrackStates(0);//FIXME?
- HelixClass helixClass;
- helixClass.Initialize_Canonical(trackStat.phi,trackStat.D0,
- trackStat.Z0,trackStat.omega,trackStat.tanLambda,
- m_genfitField->getBz({0.,0.,0.})*dd4hep::kilogauss/dd4hep::tesla);
- TLorentzVector posInit(helixClass.getReferencePoint()[0],
- helixClass.getReferencePoint()[1],
- helixClass.getReferencePoint()[2],eventStartTime);
- posInit.SetX(posInit.X()*dd4hep::mm);
- posInit.SetY(posInit.Y()*dd4hep::mm);
- posInit.SetZ(posInit.Z()*dd4hep::mm);
- TVector3 momInit(helixClass.getMomentum()[0],
- helixClass.getMomentum()[1],helixClass.getMomentum()[2]);
- momInit.SetX(momInit.x()*dd4hep::GeV);
- momInit.SetY(momInit.y()*dd4hep::GeV);
- momInit.SetZ(momInit.z()*dd4hep::GeV);
- TMatrixDSym covMInit;
- if(!createGenfitTrack(pidType,
- int(trackStat.omega/fabs(trackStat.omega)),//charge,FIXME?
- posInit,momInit,covMInit)){
+
+ ///Skip track w.o. hit
+ if(track.trackerHits_size()<=0) {
+ std::cout << " track.trackerHits_size = " << track.trackerHits_size() << std::endl;
+ if(m_debug){
+ std::cout<<"createGenfitTrackFromEDM4HepTrack skip track w/o hit"<<std::endl;
+ }
return false;
}
- return true;
+
+ //pivotToFirstLayer(mcPocaPos,mcPocaMom,firstLayerPos,firstLayerMom);
+ ///Get track parameters
+ double charge=0;
+ TVectorD seedState(6);
+ TMatrixDSym covMInit_6(6);
+
+ getTrackFromEDMTrack(track,charge,seedState,covMInit_6);///unit conversion
+ if(!isUseCovTrack) getSeedCov(covMInit_6);
+
+ bool status=createGenfitTrack(pidType,charge,seedState,covMInit_6);
+ return status;
}
-/// Add a 3d SpacepointMeasurement on TrackerHit
-bool GenfitTrack::addSpacePointTrakerHit(edm4hep::TrackerHit hit,
- int hitID)
+/// Add a 3d SpacepointMeasurement
+bool GenfitTrack::addSpacePointMeasurement(const TVector3& pos,
+ std::vector<float> sigmaUVec,std::vector<float> sigmaVVec,int cellID,int hitID)
{
- edm4hep::Vector3d pos=hit.getPosition();
- TVectorD p(3);
- p[0]=pos.x*dd4hep::mm;
- p[1]=pos.y*dd4hep::mm;
- p[2]=pos.z*dd4hep::mm;
+ TVectorD pos_smeared(3);
+ for(int i=0;i<3;i++) pos_smeared[i]=pos(i)*GenfitUnit::mm;
- if(m_debug>=2)std::cout<<m_name<<" addSpacePointTrakerHit"<<hitID
- <<"pos "<<p[0]<<" "<<p[1]<<" "<<p[2]<<" cm"<<std::endl;
/// New a SpacepointMeasurement
- double cov[6];
- for(int i=0;i<6;i++) {
- cov[i]=hit.getCovMatrix(i);
- if(m_debug>=2)std::cout<<"cov "<<cov[i]<<std::endl;
- }
- 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,
- nullptr);
- genfit::TrackPoint* trackPoint = new genfit::TrackPoint(sMeas,m_track);
- m_track->insertPoint(trackPoint);
+ /// space point error matrix and smear, unit cm
+ TMatrixDSym hitCov(3);
+ hitCov.Zero();
+
+ //get sigmas
+ float sigmaU,sigmaV;
+ getSigmas(cellID,sigmaUVec,sigmaVVec,sigmaU,sigmaV);
+ float sigmaX=sigmaU;//sigmaU*cos(atan2(pos_smeared[1],pos_smeared[0]));
+ float sigmaY=sigmaU;//*sin(atan2(pos_smeared[1],pos_smeared[0]));
+ float sigmaZ=sigmaV;
+
+ //smear 3d track position
+ bool smear=(sigmaUVec[0]>0);
+ if(smear){
+ pos_smeared[0]+=gRandom->Gaus(0,sigmaX);
+ pos_smeared[1]+=gRandom->Gaus(0,sigmaX);
+ pos_smeared[2]+=gRandom->Gaus(0,sigmaX);
+ }
+ hitCov(0,0)=sigmaX*sigmaX;//FIXME
+ hitCov(1,1)=sigmaX*sigmaX;
+ hitCov(2,2)=sigmaX*sigmaX;
+
+ if(m_debug>=2){
+ std::cout<<" addSpacePointMeasurement pos "<<std::endl;
+ pos.Print();
+ std::cout<<m_name<<" hit "<<hitID
+ <<" pos_smeared "<<pos_smeared[0]<<" "<<pos_smeared[1]
+ <<" "<<pos_smeared[2]<<" "<<" hitCov (0,0) "<<hitCov(0,0)<<" cm"<<std::endl;
+ hitCov.Print();
+ }
- if(m_debug>=2)std::cout<<"end of addSpacePointTrakerHit"<<std::endl;
+ /// new space point and TrackPoint, unit in cm
+ genfit::TrackPoint* trackPoint = new genfit::TrackPoint(
+ new genfit::SpacepointMeasurement(pos_smeared,hitCov,
+ cellID ,hitID, nullptr), m_track);
+ if(m_debug>=2){
+ std::cout<<m_name<<" add TrackPoint \n";
+ trackPoint->Print();
+ }
+ m_track->insertPoint(trackPoint);
return true;
+}//end of addSpacePointMeasurement
+
+/// Return isurface of a silicon hit
+const dd4hep::rec::ISurface*
+GenfitTrack::getISurface(edm4hep::TrackerHit* hit){
+ dd4hep::rec::SurfaceManager surfaceManager(*m_geomSvc->lcdd());
+
+ std::string detectorName;
+ unsigned long long cellID=hit->getCellID();
+ int detTypeID=getDetTypeID(hit->getCellID());
+ if(detTypeID==lcio::ILDDetID::VXD){
+ detectorName="VXD";
+ if(hit->getPosition().z>0){
+ cellID+=0x100000000;
+ }else{
+ cellID+=0x300000000;
+ }
+ }else if(detTypeID==lcio::ILDDetID::SIT){
+ detectorName="SIT";
+ }else if(detTypeID==lcio::ILDDetID::SET){
+ detectorName="SET";
+ }else if(detTypeID==lcio::ILDDetID::FTD){
+ detectorName="FTD";
+ }else{
+ std::cout << "ERROR:getISurface iSurface = NULL!" << std::endl;
+ return nullptr;
+ }
+ if(m_debug>2) std::cout<<m_name<<" detectorName "<<detectorName
+ <<" hit position mm "<<hit->getPosition()<<" hit.cellID "<<hit->getCellID()
+ <<" cellId "<<cellID<<" detTypeID "<<detTypeID<<std::endl;
+ const dd4hep::rec::SurfaceMap* surfaceMap=surfaceManager.map(detectorName);
+ auto iter=surfaceMap->find(cellID);
+ dd4hep::rec::ISurface* iSurface=nullptr;
+ if(iter!=surfaceMap->end()){iSurface=(*iter).second;}
+
+ //std::multimap< unsigned long, dd4hep::rec::ISurface*>::const_iterator it,itend;
+ //it=surfaceMap->begin();
+ //itend= surfaceMap->end();
+ //std::cout<<" print map "<<detectorName<<std::endl;
+ //for(; it!=itend; it++){
+ // dd4hep::rec::ISurface* surf = it->second;
+ // dd4hep::rec::Vector3D origin = surf->origin();
+ // std::cout <<"surf id "<< surf->id() << " origin xyz " << origin.x()
+ // << " " << origin.y() << " " << origin.z() << std::endl;
+ //}
+ return iSurface;
}
-/// Add a 3d SpacepointMeasurement with MC truth position smeared by sigma
-bool GenfitTrack::addSpacePointMeasurement(const TVectorD& pos,
- double sigma, int detID, int hitID, bool smear)
+/// Add a 1d strip or 2d pixel smeared by sigma
+ bool
+GenfitTrack::addSiliconMeasurement(edm4hep::TrackerHit* hit,
+ float sigmaU,float sigmaV,int cellID,int hitID)
{
- double sigma_t=sigma*dd4hep::mm;
- /// Convert from CEPCSW unit to genfit unit, cm
- TVectorD pos_t(3);
- pos_t(0)=pos(0)*dd4hep::mm;
- pos_t(1)=pos(1)*dd4hep::mm;
- pos_t(2)=pos(2)*dd4hep::mm;
+ if(m_debug>0) std::cout<<"addSiliconMeasurement "<<*hit<<std::endl;
- /// smear hit position with same weight
- TVectorD pos_smeared(3);
- for (int i=0;i<3;i++){
- pos_smeared[i]=pos_t(i);
- if(smear) pos_smeared[i]+=gRandom->Gaus(0,sigma_t/TMath::Sqrt(3.));
+ ///get surface by cellID
+ const dd4hep::rec::ISurface* iSurface = getISurface(hit);
+ if(nullptr==iSurface){
+ std::cout<<m_name<<" addSiliconMeasurement get surface ERROR!"<<std::endl;
+ return false;
}
+ TVector3 o,u,v;
+ double lengthU,lengthV;
+ getISurfaceOUV(iSurface,o,u,v,lengthU,lengthV);
- /// New a SpacepointMeasurement
- TMatrixDSym hitCov(3);
- hitCov(0,0)=sigma_t*sigma_t;
- hitCov(1,1)=sigma_t*sigma_t;
- hitCov(2,2)=sigma_t*sigma_t;
-
- if(m_debug>=2)std::cout<<m_name<<" addSpacePointMeasurement detID "
- <<detID<<" hitId "<<hitID<<" " <<pos_t[0]<<" "<<pos_t[1]<<" "<<pos_t[2]
- <<" cm smeared "<<pos_smeared[0]<<" "<<pos_smeared[1]<<" "
- <<pos_smeared[2]<<" sigma_t "<<sigma_t<<" cm"<<std::endl;
-
- genfit::SpacepointMeasurement* sMeas =
- new genfit::SpacepointMeasurement(pos_smeared,hitCov,detID,hitID,nullptr);
- genfit::TrackPoint* trackPoint = new genfit::TrackPoint(sMeas,m_track);
- m_track->insertPoint(trackPoint);
+ ///Get detector plane parameter
+ //VXD
+ //SIT
+ //SET
+
+ ///Get measurement and cov
+ TVectorD hitCoords(2);
+ TVector3 p;
+ TMatrixDSym hitCov(2);
+ getMeasurementAndCov(hit,p,hitCov);
+ hitCoords(0)=(p-o).Dot(u);
+ hitCoords(1)=(p-o).Dot(v);
+ if(m_debug) std::cout<<"yzhang debug hitCoords cm "<<hitCoords(0)<<" "<<hitCoords(1)<<std::endl;
+ hitCov.Zero();
+ hitCov(0,0)=sigmaU*sigmaU;
+ hitCov(1,1)=sigmaV*sigmaV;
+
+ ///Create planer finite detector plane, measurement and TrackPoint
+ genfit::RectangularFinitePlane* pixelOrStripPlane=
+ new genfit::RectangularFinitePlane(
+ -lengthU/2.,lengthU/2.,-lengthV/2.,lengthV/2.);
+ genfit::SharedPlanePtr plane(new genfit::DetPlane(o,u,v));
+ plane->setFinitePlane(pixelOrStripPlane);
+ genfit::PlanarMeasurementSDT* planarMeasurement=new genfit::PlanarMeasurementSDT(
+ hitCoords,hitCov,cellID,hitID,nullptr);
+ planarMeasurement->setTrackerHit(hit);
+ //genfit::PlanarMeasurement* planarMeasurement=new genfit::PlanarMeasurement(
+ // hitCoords,hitCov,cellID,hitID,nullptr);
+ planarMeasurement->setPlane(plane);
+ m_track->insertPoint(new genfit::TrackPoint(planarMeasurement,m_track));
+
+ if(m_debug>2){
+ std::cout<<"hitID "<<hitID<<" unit cm"<<std::endl;
+ std::cout<<"sigmaU "<<sigmaU<<" sigmaV "<<sigmaV<<std::endl;
+ std::cout<<"lengthU "<<lengthU<<" lengthV "<<lengthV<<std::endl;
+ std::cout<<"u "<<u.x()<<" "<<u.y()<<" "<<u.z()<<std::endl;
+ std::cout<<"v "<<v.x()<<" "<<v.y()<<" "<<v.z()<<std::endl;
+ std::cout<<"o "<<o.x()<<" "<<o.y()<<" "<<o.z()<<std::endl;
+ std::cout<<"p "<<p.x()<<" "<<p.y()<<" "<<p.z()<<std::endl;
+ std::cout<<"hitCoords "<<hitCoords(0)<<" "<<hitCoords(1)<<std::endl;
+ std::cout<<"hitCov "<<hitCov(0,0)<<" "<<hitCov(1,1)<<std::endl;
+ }
return true;
}
+int GenfitTrack::addSiliconMeasurements(edm4hep::Track& track,
+ std::vector<float> sigmaUVec,std::vector<float> sigmaVVec)
+{
+ ///Get TrackerHit on Track
+ int nHitAdd=0;
+ for(unsigned int iHit=0;iHit<track.trackerHits_size();iHit++){
+ edm4hep::TrackerHit* trackerHit=
+ new edm4hep::TrackerHit(track.getTrackerHits(iHit));
+ unsigned long long cellID=trackerHit->getCellID();
+ float sigmaU,sigmaV;
+ int sigmaUID=getSigmas(cellID,sigmaUVec,sigmaVVec,sigmaU,sigmaV);
+ if(0==sigmaUID) continue;//skip DC track
+ addSiliconMeasurement(trackerHit,sigmaU,sigmaV,cellID,nHitAdd++);
+ GenfitHit* genfitHit=
+ new GenfitHit(trackerHit,nullptr,nullptr,
+ nullptr,0.,0.);
+ if(nullptr==genfitHit) continue;
+ m_genfitHitVec.push_back(genfitHit);
+ }
+ return nHitAdd;
+}
-/// Add a WireMeasurement, no Unit conversion here
-void GenfitTrack::addWireMeasurement(double driftDistance,
- double driftDistanceError, const TVector3& endPoint1,
- const TVector3& endPoint2, int lrAmbig, int detID, int hitID)
+//Add wire measurement on wire, unit conversion here
+//int GenfitTrack::addWireMeasurementsFromList(podio::RelationRange<edm4hep::TrackerHit> hits,float sigma,
+int GenfitTrack::addWireMeasurementsFromList(std::vector<edm4hep::TrackerHit*>& hits,float sigma,
+ const edm4hep::MCRecoTrackerAssociationCollection* assoHits,
+ int sortMethod, bool truthAmbig,float skipCorner,float skipNear)
{
- try{
- /// New a WireMeasurement
- genfit::WireMeasurementNew* wireMeas = new genfit::WireMeasurementNew(
- driftDistance, driftDistanceError, endPoint1, endPoint2, detID,
- hitID, nullptr);
- wireMeas->setMaxDistance(0.5);//0.5 cm FIXME
- wireMeas->setLeftRightResolution(lrAmbig);
-
- if(m_debug>=2)std::cout<<m_name<<" Add wire measurement(cm) "<<hitID
- <<" ep1("<<endPoint1[0]<<" "<<endPoint1[1]<<" "<<endPoint1[2]
- <<") ep2("<<endPoint2[0]<<" "<<endPoint2[1]<<" "<<endPoint2[2]
- <<") drift "<<driftDistance<<" driftErr "<<driftDistanceError
- <<" lr "<<lrAmbig<<" detId "<<detID << " hitId "<< hitID
- <<std::endl;
+ if(m_debug>0){ std::cout<<"addWireMeasurementsFromList"<<std::endl; }
+ //podio::RelationRange<edm4hep::TrackerHit> hits_t=track.getTrackerHits();
+ std::vector<edm4hep::TrackerHit*> sortedTrackerHits;
+// sortedTrackerHits.reserve(100);
+ getSortedTrackerHits(hits,assoHits,sortedTrackerHits,sortMethod);
+
+ if(m_debug>0){
+ std::cout<<"n sortedTrackerHits "<<sortedTrackerHits.size()<<std::endl;
+ }
+ int nHitAdd=0;
+ for(auto trackerHit : sortedTrackerHits){
+ edm4hep::SimTrackerHit simTrackerHitAsso;
+ getAssoSimTrackerHit(assoHits,trackerHit,simTrackerHitAsso);
+ //FIXME driftVelocity
+ GenfitHit* genfitHit=
+ makeAGenfitHit(trackerHit,&simTrackerHitAsso,sigma,truthAmbig,
+ skipCorner,skipNear);
+ if(nullptr==genfitHit) continue;
+ m_genfitHitVec.push_back(genfitHit);
///New a TrackPoint,create connection between meas. and trackPoint
- genfit::TrackPoint* trackPoint=new genfit::TrackPoint(wireMeas,m_track);
- wireMeas->setTrackPoint(trackPoint);
-
- m_track->insertPoint(trackPoint);
-
- }catch(genfit::Exception& e){
- if(m_debug>=2)std::cout<<m_name
- <<"Add wire measurement exception"<<std::endl;
- e.what();
- }
-}//end of addWireMeasurementOnTrack
+ WireMeasurementDC* wireMeasurementDC=
+ new WireMeasurementDC(genfitHit,nHitAdd++);
+ m_track->insertPoint(new genfit::TrackPoint(wireMeasurementDC,m_track));
+ if(m_debug>=2){ std::cout<<nHitAdd-1; wireMeasurementDC->Print(); }
+ }//end of loop over sotred hits
+ return nHitAdd;
+}//end of addWireMeasurementsFromList
//Add wire measurement on wire, unit conversion here
-bool GenfitTrack::addWireMeasurementOnTrack(edm4hep::Track track,double sigma)
+int GenfitTrack::addWireMeasurementsOnTrack(edm4hep::Track& track,float sigma,
+ const edm4hep::MCRecoTrackerAssociationCollection* assoHits,
+ int sortMethod, bool truthAmbig,float skipCorner,float skipNear)
{
- for(unsigned int iHit=0;iHit<track.trackerHits_size();iHit++){
- edm4hep::TrackerHit hit=track.getTrackerHits(iHit);
+ if(m_debug>0){ std::cout<<"addWireMeasurementsOnTrack"<<std::endl; }
- double driftVelocity=40;//FIXME, TODO, um/ns
- double driftDistance=hit.getTime()*driftVelocity*dd4hep::um;
- TVector3 endPointStart(0,0,0);
- TVector3 endPointEnd(0,0,0);
- m_gridDriftChamber->cellposition(hit.getCellID(),endPointStart,
- endPointEnd);
- int lrAmbig=0;
- if(m_debug>=2)std::cout<<m_name<<" time "<<hit.getTime()
- <<" driftVelocity " <<driftVelocity<<std::endl;
- if(m_debug>=2)std::cout<<m_name<<" wire pos " <<endPointStart.X()
- <<" "<<endPointStart.Y()<<" " <<endPointStart.Z()<<" "
- <<endPointEnd.X()<<" " <<endPointEnd.Y()<<" "
- <<endPointEnd.Z()<<" " <<std::endl;
- endPointStart.SetX(endPointStart.x()*dd4hep::cm);
- endPointStart.SetY(endPointStart.y()*dd4hep::cm);
- endPointStart.SetZ(endPointStart.z()*dd4hep::cm);
- endPointEnd.SetX(endPointEnd.x()*dd4hep::cm);
- endPointEnd.SetY(endPointEnd.y()*dd4hep::cm);
- endPointEnd.SetZ(endPointEnd.z()*dd4hep::cm);
- addWireMeasurement(driftDistance,sigma*dd4hep::cm,endPointStart,
- endPointEnd,lrAmbig,hit.getCellID(),iHit);
+ if(m_debug>0){
+ std::cout<<"n trackerHits size"<<track.trackerHits_size()<<std::endl;
}
- return true;
-}//end of addWireMeasurementOnTrack of Track
+ int nHitAdd=0;
+ if(0==track.trackerHits_size()){
+ if(m_debug>0) std::cout<<"No hit on track"<<std::endl;
+ return nHitAdd;
+ }
+
+ podio::RelationRange<edm4hep::TrackerHit> hits_t=track.getTrackerHits();
+ std::vector<edm4hep::TrackerHit*> hits;
+ //hits.reserve(1000);
+ for(auto &h:hits_t){
+ edm4hep::TrackerHit* hit = const_cast<edm4hep::TrackerHit*>(&h);
+ hits.push_back(hit);
+ }
+ nHitAdd=addWireMeasurementsFromList(hits,sigma,assoHits,sortMethod,truthAmbig,skipCorner,skipNear);
+
+ return nHitAdd;
+}//end of addWireMeasurements
/// Get MOP
-bool GenfitTrack::getMOP(int hitID,
- genfit::MeasuredStateOnPlane& mop, genfit::AbsTrackRep* trackRep) const
+bool GenfitTrack::getMOP(int hitID,genfit::MeasuredStateOnPlane& mop,
+ genfit::AbsTrackRep* trackRep) const
{
+
if(nullptr == trackRep) trackRep = getRep();
try{
mop = m_track->getFittedState(hitID,trackRep);
@@ -354,30 +441,22 @@ bool GenfitTrack::getMOP(int hitID,
}
/// New and add a track representation to track
-genfit::RKTrackRep* GenfitTrack::addTrackRep(int pdg)
+genfit::RKTrackRep* GenfitTrack::addTrackRep(int pdgType,int charge)
{
- /// create a new track representation
- genfit::RKTrackRep* rep = new genfit::RKTrackRep(pdg);
- m_reps.push_back(rep);
+ int chargeId=0;
+ charge>0 ? chargeId=0 : chargeId=1;//s_PDG[0]: positive particle
+
+ genfit::RKTrackRep* rep=new genfit::RKTrackRep(s_PDG[chargeId][pdgType]);
m_track->addTrackRep(rep);
- //try{
- // genfit::MeasuredStateOnPlane stateInit(rep);
- // rep->setPosMomCov(stateInit, pos, mom, covM);
- //}catch(genfit::Exception e){
- // if(m_debug>=2)std::cout<<m_name<<" Exception in set track status"
- // <<std::endl ;
- // std::cout<<e.what()<<std::endl;
- // return false;
- //}
return rep;
}
/// Get the position from genfit::Track::getStateSeed
const TLorentzVector GenfitTrack::getSeedStatePos()const
{
- TVectorD seedStat(6);
- seedStat = m_track->getStateSeed();
- TVector3 p(seedStat[0],seedStat[1],seedStat[2]);
+ TVectorD seedState(6);
+ seedState = m_track->getStateSeed();
+ TVector3 p(seedState[0],seedState[1],seedState[2]);
p = p*dd4hep::cm;
TLorentzVector pos(p[0],p[1],p[2],9999);//FIXME
return pos;
@@ -386,18 +465,18 @@ const TLorentzVector GenfitTrack::getSeedStatePos()const
/// Get the momentum from genfit::Track::getStateSeed
const TVector3 GenfitTrack::getSeedStateMom() const
{
- TVectorD seedStat(6); seedStat = m_track->getStateSeed();
- TVector3 mom(seedStat[3],seedStat[4],seedStat[5]);
+ TVectorD seedState(6); seedState = m_track->getStateSeed();
+ TVector3 mom(seedState[3],seedState[4],seedState[5]);
return mom*dd4hep::GeV;
}
/// Get the seed states of momentum and position
void GenfitTrack::getSeedStateMom(TLorentzVector& pos, TVector3& mom) const
{
- TVectorD seedStat(6); seedStat = m_track->getStateSeed();
- mom = TVector3(seedStat[3],seedStat[4],seedStat[5])*dd4hep::GeV;
- seedStat = m_track->getStateSeed();
- TVector3 p = TVector3(seedStat[0],seedStat[1],seedStat[2])*dd4hep::cm;
+ TVectorD seedState(6); seedState = m_track->getStateSeed();
+ mom = TVector3(seedState[3],seedState[4],seedState[5])*dd4hep::GeV;
+ seedState = m_track->getStateSeed();
+ TVector3 p = TVector3(seedState[0],seedState[1],seedState[2])*dd4hep::cm;
pos.SetXYZT(p[0],p[1],p[2],9999);//FIXME time
}
@@ -406,6 +485,23 @@ unsigned int GenfitTrack::getNumPoints() const
return m_track->getNumPoints();
}
+GenfitHit* GenfitTrack::GetHit(long unsigned int i) const {
+ if(i>=m_genfitHitVec.size())return nullptr;
+ return m_genfitHitVec[i];
+}
+
+unsigned int GenfitTrack::getNumPointsDet(int detTypeID) const
+{
+ unsigned int nHit=0;
+ const std::vector<genfit::TrackPoint*> tps=m_track->getPoints();
+ for(auto tp:tps){
+ const genfit::AbsMeasurement* m=nullptr;
+ if(tp->hasRawMeasurements()) m=tp->getRawMeasurement();
+ if(nullptr!=m && detTypeID==getDetTypeID(m->getDetId())) nHit++;
+ }
+ return nHit;
+}
+
/// Test the fit result FIXME
bool GenfitTrack::fitSuccess(int repID) const
{
@@ -417,8 +513,8 @@ bool GenfitTrack::fitSuccess(int repID) const
||fitStatus->isFitConvergedFully()) {
if(m_debug>=2)std::cout<<m_name<< "Fitting is failed... isFitted"
<<fitStatus->isFitted()<<" , isFitConverged "
- <<fitStatus->isFitConverged()<<", isFitConvergedFully "
- <<fitStatus->isFitConvergedFully()<<std::endl;
+ <<fitStatus->isFitConverged()<<", isFitConvergedFully "
+ <<fitStatus->isFitConvergedFully()<<std::endl;
return false;
}
@@ -436,18 +532,40 @@ bool GenfitTrack::fitSuccess(int repID) const
return true;
}
-void GenfitTrack::setDebug(int debug)
+void GenfitTrack::setDebugGenfit(int debug)
{
- m_debug = debug;
- for(unsigned int i=0;i<m_reps.size();i++){ m_reps[i]->setDebugLvl(debug); }
+ if(m_track){
+ for(unsigned int i=0;i<m_track->getNumReps();i++){
+ m_track->getTrackRep(i)->setDebugLvl(debug);
+ }
+ }
+}
+void GenfitTrack::setDebugLocal(int debug){
+#ifdef GENFIT_MY_DEBUG
+ if(m_track){
+ for(unsigned int i=0;i<m_track->getNumReps();i++){
+ //m_track->getTrackRep(i)->setDebugLvlLocal(debug);
+ }
+ }
+ m_debugLocal=debug;
+#endif
+ if(m_debug)std::cout<<"GenfitTrack:setDebugLvlLocal "<<debug<<std::endl;
}
void GenfitTrack::printSeed() const
{
TLorentzVector pos = getSeedStatePos();
TVector3 mom = getSeedStateMom();
+ std::cout << " Seed pos " << std::endl;
+ pos.Print();
+ std::cout << " Seed Mom = " << std::endl;
+ mom.Print();
print(pos,mom);
- if(m_debug>=2)std::cout<<m_name<<" NumPoints "<<getNumPoints()<<std::endl;
+ TMatrixDSym covSeed=m_track->getCovSeed();
+ //if(m_debug>1)
+ std::cout << " covSeed = " << std::endl;
+ covSeed.Print();
+ //std::cout<<" pdg "<<0<<getRep(0)<<std::endl;//FIXME
}
void GenfitTrack::printFitted(int repID) const
@@ -479,7 +597,9 @@ void GenfitTrack::print( TLorentzVector pos, TVector3 mom,
if(m_debug>=2)std::cout<<m_name<<" "<<comment<<std::endl;
if(m_debug>1){
- for(unsigned int i=0;i<m_reps.size();i++){ m_reps[i]->Print(); }
+ for(unsigned int i=0;i<m_track->getNumReps();i++){
+ m_track->getTrackRep(i)->Print();
+ }
}
//for(unsigned int i=0; i<m_track->getNumPoints(); i++){
// m_track->getPoint(i)->print();
@@ -488,7 +608,7 @@ void GenfitTrack::print( TLorentzVector pos, TVector3 mom,
/// Get position, momentum, cov on plane of hitID-th hit
bool GenfitTrack::getPosMomCovMOP(int hitID, TLorentzVector& pos,
- TVector3& mom, TMatrixDSym& cov, int repID) const
+ TVector3& mom, TMatrixDSym& cov,int repID) const
{
TVector3 p;
genfit::MeasuredStateOnPlane mop;
@@ -505,7 +625,9 @@ bool GenfitTrack::getPosMomCovMOP(int hitID, TLorentzVector& pos,
int GenfitTrack::getNumPointsWithFittedInfo(int repID) const
{
int nHitWithFittedInfo = 0;
+ //check number of hit with fitter info
int nHit = m_track->getNumPointsWithMeasurement();
+ //check number of hit with fitter info
for(int i=0; i<nHit; i++){
if(nullptr != m_track->getPointWithFitterInfo(i,getRep(repID))){
nHitWithFittedInfo++;
@@ -515,7 +637,7 @@ int GenfitTrack::getNumPointsWithFittedInfo(int repID) const
}
int GenfitTrack::getFittedState(TLorentzVector& pos, TVector3& mom,
- TMatrixDSym& cov, int repID, bool biased) const
+ TMatrixDSym& cov, int trackPointId, int repID, bool biased) const
{
//check number of hit with fitter info
if(getNumPointsWithFittedInfo(repID)<=2) return 1;
@@ -527,14 +649,13 @@ int GenfitTrack::getFittedState(TLorentzVector& pos, TVector3& mom,
//get first or last measured state on plane
genfit::MeasuredStateOnPlane mop;
try{
- mop = m_track->getFittedState(biased);
+ mop = m_track->getFittedState(trackPointId,rep,biased);
}catch(genfit::Exception& e){
- std::cout<<" getNumPointsWithFittedInfo "
+ std::cout<<" getNumPointsWithFittedInfo="
<<getNumPointsWithFittedInfo(repID)
<<" no TrackPoint with fitted info "<<std::endl;
if(m_debug>=2)std::cout<<m_name
<<"Exception in getFittedState"<<std::endl;
- std::cout<<e.what()<<std::endl;
return 3;
}
@@ -557,12 +678,12 @@ int GenfitTrack::getDetIDWithFitterInfo(int hitID, int idRaw) const
int GenfitTrack::getPDG(int id) const
{
- return m_reps[id]->getPDG();
+ return m_track->getTrackRep(id)->getPDG();
}
int GenfitTrack::getPDGCharge(int id) const
{
- return m_reps[id]->getPDGCharge();
+ return m_track->getTrackRep(id)->getPDGCharge();
}
const genfit::FitStatus*
@@ -572,264 +693,701 @@ GenfitTrack::getFitStatus(int repID) const
}
/// Extrapolate track to the cloest point of approach(POCA) to the wire of hit,
-/// return StateOnPlane of this POCA
-/// inputs
-/// pos,mom ... position & momentum at starting point, unit= [mm]&[GeV/c]
-/// (converted to [cm]&[GeV/c] inside this function)
-/// hit ... destination
-/// outputs poca [mm] (converted from [cm] in this function) ,global
-double GenfitTrack::extrapolateToHit( TVector3& poca, TVector3& pocaDir,
- TVector3& pocaOnWire, double& doca, TVector3 pos, TVector3 mom,
- TVector3 end0,//one end of the hit wire
- TVector3 end1,//the orhter end of the hit wire
- int debug,
- int repID,
- bool stopAtBoundary,
- bool calcJacobianNoise)const
+/// return doca, poca on wire and track
+/// outputs poca [cm]
+/// unit conversion here
+double GenfitTrack::extrapolateToHit(TVector3& poca, TVector3& pocaDir,
+ TVector3& pocaOnWire, double& doca,edm4hep::MCParticle mcParticle,
+ int cellID, int repID, bool stopAtBoundary, bool calcJacobianNoise)const
{
-
- //genfit::MeasuredStateOnPlane state = getMOP(iPoint); // better?
- //genfit::MeasuredStateOnPlane state = getMOP(0); // better?
- //To do the extrapolation without IHitSelection,above 2 lines are not used.
- pos = pos*dd4hep::cm;//FIXME
- mom = mom*dd4hep::GeV;
-
- //std::cout<<__LINE__<<" extrapolate to Hit pos and mom"<<std::endl;
- //pos.Print();
- //mom.Print();
-
- genfit::AbsTrackRep* rep = new genfit::RKTrackRep(getRep(repID)->getPDG());
- rep->setDebugLvl(debug);
- genfit::MeasuredStateOnPlane state(rep);
+ ///Get MCParticle and convert unit
+ TVector3 pos;
+ TVector3 mom;
+ getPosMomFromMCPartile(mcParticle,pos,mom);
+
+ TVector3 end0(0,0,0);
+ TVector3 end1(0,0,0);
+ getEndPointsOfWire(cellID,end0,end1);
+
+ int chargeId;
+ mcParticle.getCharge() >0 ? chargeId=0 : chargeId=1;//s_PDG[0]: positive particle
+ genfit::RKTrackRep* rep = new genfit::RKTrackRep(s_PDG[chargeId][repID]);
+ //genfit::MeasuredStateOnPlane state(rep);
+ genfit::StateOnPlane state(rep);
rep->setPosMom(state, pos, mom);
-
- //genfit::MeasuredStateOnPlane state(m_track->getRep(repID));
- //m_track->getRep(repID)->setPosMom(state, pos, mom);
-
- //m_track->PrintSeed();
- double extrapoLen(0);
- //std::cout<<" wire1 "<<std::endl;
- //end0.Print();
- //std::cout<<" wire0 "<<std::endl;
- //end1.Print();
- //std::cout<<" state "<<std::endl;
- //state.Print();
+ if(m_debug){
+ std::cout<<"GenfitterTrack::extrapolateToHit before pos and mom "
+ <<" charge "<<(int)mcParticle.getCharge()<<" repID "<<repID
+ <<" pdg "<<s_PDG[chargeId][repID]<<std::endl;
+ pos.Print();
+ mom.Print();
+ std::cout<<" end point"<<std::endl;
+ end0.Print();
+ end1.Print();
+ std::cout<<" state "<<std::endl;
+ state.Print();
+ std::cout<<" stopAtBoundary "<<stopAtBoundary<<std::endl;
+ std::cout<<" calcJacobianNoise "<<calcJacobianNoise<<std::endl;
+ }
- // forth
+ double extrapoLen=0.;
try {
- //genfit::RKTrackRep* rkRep =
- // dynamic_cast<genfit::RKTrackRep*> (m_track->getRep(repID));
- //extrapoLen = rkRep->extrapolateToLine(state, end0, end1, poca,
- // pocaDir, pocaOnWire, stopAtBoundary, calcJacobianNoise);
- extrapoLen = rep->extrapolateToLine(state, end0, end1, poca,
- pocaDir, pocaOnWire, stopAtBoundary, calcJacobianNoise);
+ extrapoLen = rep->extrapolateToLine(state,end0,end0-end1,
+ poca,pocaDir,pocaOnWire,stopAtBoundary,calcJacobianNoise);
}catch(genfit::Exception& e) {
if(m_debug>=3)std::cout<<
- "Exception in GenfigFitter::ExtrapolateToHit"<<e.what()<<std::endl;
- return extrapoLen;
+ "Exception in GenfitterTrack::extrapolateToHit"<<e.what()<<std::endl;
+ return extrapoLen/GenfitUnit::mm*dd4hep::mm;
}
- //poca = poca*(dd4hep::cm);
- //pocaOnWire = pocaOnWire*(dd4hep::cm);
- pocaOnWire = pocaOnWire;
doca = (pocaOnWire-poca).Mag();
- //TODO: debug pocaOnWire
- if(m_debug>=2)std::cout<< " poca "<<poca.x()<<","<<poca.y()
- <<" "<<poca.z()<<" doca "<<doca<<std::endl;
- if(m_debug>=2)std::cout<< " pocaOnWire "<<pocaOnWire.x()
- <<","<<pocaOnWire.y()<<" "<<pocaOnWire.z()<<" doca "<<doca<<std::endl;
+ if(m_debug>=2){
+ std::cout<< " poca "<<poca.x()<<","<<poca.y()<<" "<<poca.z()
+ <<" pocaDir "<<pocaDir.x()<<","<<pocaDir.y()<<" "<<pocaDir.z()
+ <<" pocaOnWire "<<pocaOnWire.x()<<","<<pocaOnWire.y()<<" "<<pocaOnWire.z()
+ <<" doca "<<doca<<" cm"<<std::endl;
+ }
- return extrapoLen*(dd4hep::cm);
-}//end of ExtrapolateToHit
+ delete rep;
+ return extrapoLen/GenfitUnit::mm*dd4hep::mm;
+}//end of extrapolateToHit
+
+/////Add space point measurement from edm4hep::Track to genfit track
+//int GenfitTrack::addHitsOnEdm4HepTrack(const edm4hep::Track& track,
+// const edm4hep::MCRecoTrackerAssociationCollection* assoHits,
+// std::vector<float> sigma,bool smear,
+// bool measurementTypeSi, bool measurementTypeDC){
+// ///Get TrackerHit on Track
+// int hitID=0;
+// for(unsigned int iHit=0;iHit<track.trackerHits_size();iHit++){
+// edm4hep::TrackerHit hit=track.getTrackerHits(iHit);
+// ///Get hit type
+// int detTypeID=getDetTypeID(hit.getCellID());
+// if(m_debug>=2)std::cout<<"addHitsOnEdm4HepTrack "<<iHit<<" hit "<<hit
+// <<" detTypeID "<<detTypeID<<" type "<<hit.getType()<<std::endl;
+//
+// bool hitIsSpapcePoint=UTIL::BitSet32(hit.getType())[
+// UTIL::ILDTrkHitTypeBit::COMPOSITE_SPACEPOINT];
+// bool hitIsPlanar=UTIL::BitSet32(hit.getType())[
+// UTIL::ILDTrkHitTypeBit::ONE_DIMENSIONAL];
+// if(m_debug>2){
+// std::cout<<detTypeID<<" COMPOSITE_SPACEPOINT "<<hitIsSpapcePoint
+// <<std::endl;
+// std::cout<<detTypeID<<" ONE_DIMENSIONAL "<<hitIsPlanar<<std::endl;
+// }
+//
+// }
+//
+// return 1;
+//}
+
+///Add space point measurement of silicon from edm4hep::Track to genfit track
+int GenfitTrack::addSpacePointsSi(const edm4hep::Track& track,
+ std::vector<float> sigmaU,std::vector<float> sigmaV)
+{
+ ///Get TrackerHit on Track
+ int nHitAdd=0;
+ for(unsigned int iHit=0;iHit<track.trackerHits_size();iHit++){
+ int detTypeID=getDetTypeID(track.getTrackerHits(iHit).getCellID());
+ if(m_geomSvc->lcdd()->constant<int>("DetID_DC")==detTypeID) continue;
+ edm4hep::TrackerHit hit=track.getTrackerHits(iHit);
+ edm4hep::Vector3d pos=hit.getPosition();
+
+ //edm4hep::Vector3d pos=simTrackerHitAsso.getPosition();
+
+ TVector3 p(pos.x,pos.y,pos.z);
+ p.Print();
+ unsigned long long cellID = hit.getCellID();
+ if(addSpacePointMeasurement(p,sigmaU,sigmaV,cellID,nHitAdd)){
+ if(m_debug>=2)std::cout<<"add silicon space point"<<std::endl;
+ nHitAdd++;
+ }else{
+ if(m_debug>=2)std::cout<<"addSpacePointMeasurement"
+ <<cellID<<" faieled" <<std::endl;
+ }
+ }
+ if(m_debug>=2) std::cout<<m_name<<" Si addHitAdd="<<nHitAdd<<std::endl;
+ return nHitAdd;
+}//end of addSpacePointsSi
-///Add space point measurement from edm4hep::Track to genfit track
-int GenfitTrack::addSimTrackerHits(edm4hep::Track track,
+///Add drift chamber space point from edm4hep::track
+int GenfitTrack::addSpacePointsDC(const edm4hep::Track& track,
const edm4hep::MCRecoTrackerAssociationCollection* assoHits,
- float sigma,bool smear){
- //A TrakerHit collection
- std::vector<edm4hep::SimTrackerHit> sortedDCTrackHitCol;
+ std::vector<float> sigmaU,std::vector<float> sigmaV)
+{
+ if(m_debug>=2){
+ std::cout<<m_name<<" 5. addSpacePointsDC nTrackerHit="
+ <<track.trackerHits_size()<<std::endl;
+ }
- if(m_debug>=2)std::cout<<m_name<<" VXD "
- <<lcio::ILDDetID::VXD<<" SIT "
- <<lcio::ILDDetID::SIT<<" SET "
- <<lcio::ILDDetID::SET<<" FTD "
- <<lcio::ILDDetID::FTD<<" "<<std::endl;
- ///Get TrackerHit on Track
- int hitID=0;
+ ///Get TrackerHit with min. time in Cell
+ std::vector<edm4hep::SimTrackerHit> sortedDCTrackHitCol;
for(unsigned int iHit=0;iHit<track.trackerHits_size();iHit++){
edm4hep::TrackerHit hit=track.getTrackerHits(iHit);
-
- UTIL::BitField64 encoder(lcio::ILDCellID0::encoder_string);
- encoder.setValue(hit.getCellID());
- int detID=encoder[lcio::ILDCellID0::subdet];
- if(m_debug>=2)std::cout<<m_name<<" "<<iHit<<" hit "<<hit
- <<" detID "<<detID<<std::endl;
-
- if(detID==lcio::ILDDetID::VXD || detID==lcio::ILDDetID::SIT ||
- detID==lcio::ILDDetID::SET || detID==lcio::ILDDetID::FTD){
- if(addSpacePointTrakerHit(hit,hitID)){
- if(m_debug>=2)std::cout<<"add slicon space point"<<std::endl;
- hitID++;
- }else{
- if(m_debug>=2)std::cout<<"silicon addSpacePointTrakerHit"
- <<detID<<" "<<hit.getCellID()<<" faieled"<<std::endl;
- }
- }else if(detID==7){
- //if(addSpacePointMeasurement(p,sigma,hit.getCellID(),hitID)){
- // if(m_debug>=2)std::cout<<"add DC space point"<<std::endl;
- // hitID++;
- //}else{
- // if(m_debug>=2)std::cout<<"addSpacePointMeasurement"
- // <<detID<<" faieled" <<std::endl;
- //}
- float minTime=FLT_MAX;
- edm4hep::SimTrackerHit 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();
- minTime=assoHits->at(iSimHit).getSim().getTime();
- }
+ int detTypeID=getDetTypeID(track.getTrackerHits(iHit).getCellID());
+ if(m_geomSvc->lcdd()->constant<int>("DetID_DC")!=detTypeID) continue;
+
+ ///Select the SimTrakerHit with least time
+ float minTime=FLT_MAX;
+ edm4hep::SimTrackerHit minTimeSimHit;
+ 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();
+ minTime=assoHits->at(iSimHit).getSim().getTime();
}
- //std::cout<<"minTimeSimHit "<<minTimeSimHit<<std::endl;
+ }
+ if(!minTimeSimHit.isProducedBySecondary()){
sortedDCTrackHitCol.push_back(minTimeSimHit);
- }else{
- if(m_debug>=2)std::cout<<" Skip add this hit!"<<std::endl;
}
- }//end loop over hit on track
+ }//end loop over TrackerHit
- if(m_debug>=2)std::cout<<" addSimTrakerHits trackerHits_size="
- <<track.trackerHits_size()<<std::endl;
+ ///Sort DC SimTrackerHit by time
+ std::sort(sortedDCTrackHitCol.begin(),sortedDCTrackHitCol.end(),sortDCSimHit);
///Add DC hits to track
- //Sort sim DC hits by time
- //std::sort(sortedDCTrackHitCol.begin(),sortedDCTrackHitCol.end(),sortDCHit);
+ int nHitAdd=0;
for(auto dCTrackerHit: sortedDCTrackHitCol){
edm4hep::Vector3d pos=dCTrackerHit.getPosition();
- TVectorD p(3);
- p[0]=pos.x;
- p[1]=pos.y;
- p[2]=pos.z;
- unsigned long long detID = dCTrackerHit.getCellID();
- if(addSpacePointMeasurement(p,sigma,detID,hitID,smear)){
+ TVector3 p(pos.x,pos.y,pos.z);
+
+ if(m_debug) std::cout<<"6. addSpacePointsDC hit "<<dCTrackerHit<<std::endl;
+ unsigned long long cellID=dCTrackerHit.getCellID();
+ if(addSpacePointMeasurement(p,sigmaU,sigmaV,dCTrackerHit.getCellID()
+ ,nHitAdd)){
if(m_debug>=2)std::cout<<"add DC space point"<<std::endl;
- hitID++;
+ nHitAdd++;
}else{
if(m_debug>=2)std::cout<<"addSpacePointMeasurement"
- <<detID<<" faieled" <<std::endl;
+ <<cellID<<" faieled" <<std::endl;
+ }
+ }
+
+ if(m_debug>=2) std::cout<<m_name<<" DC addHitAdd="<<nHitAdd<<std::endl;
+ return nHitAdd;
+}//end of addSpacePointsDC
+
+//double GenfitTrack::extrapolateToPoint(TVector3& pos, TVector3& mom,
+// TMatrixDSym& cov,
+// const TVector3& point,
+// int repID,// same with pidType
+// bool stopAtBoundary,
+// bool calcJacobianNoise) const
+//{
+// return extrapolateToPoint(pos,mom,cov,point,repID,stopAtBoundary,
+// calcJacobianNoise);
+//
+//}//end of extrapolateToPoint
+
+double GenfitTrack::extrapolateToPoint(TVector3& pos, TVector3& mom,
+ TMatrixDSym& cov, const TVector3& point,
+ int repID,// same with pidType
+ bool stopAtBoundary,
+ bool calcJacobianNoise) const
+{
+ double trackLength(1e9*dd4hep::cm);
+ if(!getFitStatus(repID)->isFitted()) return trackLength;
+ try{
+ // get track rep
+ genfit::AbsTrackRep* rep = getRep(repID);
+ if(nullptr == rep) {
+ if(m_debug>=2)std::cout<<"In extrapolateToPoint rep "
+ <<repID<<" not exist!"<<std::endl;
+ return trackLength*dd4hep::cm;
+ }
+
+ /// extrapolate to point
+ //genfit::StateOnPlane state(*(&(track->getTrack()->getFittedState(0,rep))));
+
+ // get track point with fitter info
+ genfit::TrackPoint* tp = getTrack()->getPointWithFitterInfo(0,rep);
+ if(nullptr == tp) {
+ if(m_debug>=2)std::cout<<
+ "In extrapolateToPoint TrackPoint is null"<<std::endl;
+ return trackLength*dd4hep::cm;//FIXME unit
+ }
+
+ // get fitted state on plane of this track point
+ genfit::KalmanFittedStateOnPlane* state =
+ static_cast<genfit::KalmanFitterInfo*>(
+ tp->getFitterInfo(rep))->getBackwardUpdate();
+ genfit::StateOnPlane orignalState(*state);
+ if(m_debug>3){
+ tp->Print();
+ std::cout<<" original state before extrapolate "<<std::endl;
+ state->Print();
+ }
+
+ if(nullptr == state) {
+ if(m_debug>=2)std::cout<<
+ "In extrapolateToPoint KalmanFittedStateOnPlane is null"<<std::endl;
+ return trackLength*dd4hep::cm;//FIXME unit
+ }
+ //rep->setDebugLvl(10);
+ trackLength = rep->extrapolateToPoint(*state,
+ point*(1/dd4hep::cm),stopAtBoundary, calcJacobianNoise);
+
+ rep->getPosMomCov(*state,pos,mom,cov);//FIXME exception exist
+
+ pos = pos*dd4hep::cm;//FIXME unit
+ mom = mom*dd4hep::GeV;//FIXME unit
+ if(m_debug>3){
+ std::cout<<" original state before extrapolate "<<std::endl;
+ orignalState.Print();
+ std::cout<<" extrapolated state "<<std::endl;
+ state->Print();
}
+
+ } catch(genfit::Exception& e){
+ if(m_debug>=3)std::cout
+ <<"Exception in GenfitTrack::extrapolateToPoint"
+ << e.what()<<std::endl;
+ trackLength = 1e9*dd4hep::cm;
}
+ return trackLength*dd4hep::cm;
+}//end of extrapolateToPoint
+
+/// Extrapolate the track to the cyliner at fixed raidus
+/// position & momentum as starting point
+/// position and momentum at global coordinate in dd4hepUnit
+/// return trackLength in dd4hepUnit
+ double
+GenfitTrack::extrapolateToCylinder(TVector3& pos, TVector3& mom,
+ double radius, const TVector3 linePoint,
+ const TVector3 lineDirection, int hitID, int repID,
+ bool stopAtBoundary, bool calcJacobianNoise)
+{
+ double trackLength(1e9*dd4hep::cm);//FIXME unit
+ if(!getFitStatus(repID)->isFitted()) return trackLength;
+ try{
+ // get track rep
+ genfit::AbsTrackRep* rep = getRep(repID);
+ if(nullptr == rep) {
+ if(m_debug>=2)std::cout<<"In extrapolateToCylinder rep is null"
+ <<std::endl;
+ return trackLength*dd4hep::cm;//FIXME unit
+ }
- if(m_debug>=2)std::cout<<"GenfitTrack nHitAdded "<<hitID<<std::endl;
- return hitID;
+ // get track point with fitter info
+ genfit::TrackPoint* tp =
+ getTrack()->getPointWithFitterInfo(hitID,rep);
+ if(nullptr == tp) {
+ if(m_debug>=2)std::cout<<
+ "In extrapolateToCylinder TrackPoint is null"<<std::endl;
+ return trackLength*dd4hep::cm;//FIXME unit
+ }
+
+ // get fitted state on plane of this track point
+ genfit::KalmanFittedStateOnPlane* state =
+ static_cast<genfit::KalmanFitterInfo*>(
+ tp->getFitterInfo(rep))->getBackwardUpdate();
+
+ if(nullptr == state){
+ if(m_debug>=2)std::cout<<"In extrapolateToCylinder "
+ <<"no KalmanFittedStateOnPlane in backwardUpdate"<<std::endl;
+ return trackLength*dd4hep::cm;
+ }
+ rep->setPosMom(*state, pos*(1/dd4hep::cm), mom*(1/dd4hep::GeV));//FIXME unit
+
+ /// extrapolate
+ trackLength = rep->extrapolateToCylinder(*state,
+ radius/dd4hep::cm, linePoint*(1/dd4hep::cm), lineDirection,
+ stopAtBoundary, calcJacobianNoise);
+ // get pos&mom at extrapolated point on the cylinder
+ rep->getPosMom(*state,pos,mom);//FIXME exception exist
+ pos = pos*dd4hep::cm;
+ mom = mom*dd4hep::GeV;
+ } catch(genfit::Exception& e){
+ if(m_debug>=3)std::cout
+ <<"Exception in GenfitTrack::extrapolateToCylinder "
+ << e.what()<<std::endl;
+ trackLength = 1e9*dd4hep::cm;
+ }
+ return trackLength*dd4hep::cm;
+}
+
+bool GenfitTrack::debugDistance(const edm4hep::TrackerHitCollection* dCDigiCol,
+ int& nFittedDC,int& nFittedSDT,int& ngenfitHit,
+ std::vector<double>& smearDistance,
+ std::vector<double>& truthDistance,double driftVelocity)
+{
+
+ int DCHit = 0;
+ int SDTHit = 0;
+ unsigned int nPoints = m_track->getNumPoints();
+ for(unsigned int i = 0; i<nPoints; i++){
+ genfit::TrackPoint* point = m_track->getPoint(i);
+ genfit::AbsMeasurement* absMea = point->getRawMeasurement();
+ genfit::PlanarMeasurementSDT* sdtMea =
+ dynamic_cast<genfit::PlanarMeasurementSDT*>(absMea);
+ if(sdtMea){
+ const edm4hep::TrackerHit* TrackerHit_ = sdtMea->getTrackerHit();
+ SDTHit++;
+ }else{
+ WireMeasurementDC* dcMea =
+ dynamic_cast<WireMeasurementDC*>(absMea);
+ if(dcMea){
+ const edm4hep::TrackerHit* TrackerHit_ = dcMea->getTrackerHit();
+ smearDistance.push_back(1e-3*driftVelocity*TrackerHit_->getTime());
+ DCHit++;
+ for(auto dcDigi: *dCDigiCol){
+ if(dcDigi.getCellID() == TrackerHit_->getCellID())
+ {
+ truthDistance.push_back(1e-3*driftVelocity*(dcDigi.getTime()));
+ }
+ }
+ }
+ }
+ }
+
+ DCHit = nFittedDC;
+ SDTHit = nFittedSDT;
+ ngenfitHit = nFittedDC+nFittedDC;
+
+ return true;
+}
+
+/// Get doca on plane of hitID-th hit
+bool GenfitTrack::GetDocaRes(int hitID, double& DriftDis,double& fittedDoca,
+ double& res, int repID, bool biased) const
+{
+ genfit::TrackPoint* trackPoint = m_track->getPointWithFitterInfo(hitID, getRep(repID));
+ if(!trackPoint)return false;
+ genfit::AbsMeasurement* thismea = trackPoint->getRawMeasurement(0);
+ const TVectorD measereddoca = thismea->getRawHitCoords();
+ DriftDis = abs(measereddoca[0]);
+ genfit::MeasuredStateOnPlane mop;
+ if(!getMOP(hitID,mop,getRep(repID))) return false;
+ TVectorD state = mop.getState();
+ if(state.GetNrows() != 5)return false;
+ fittedDoca = abs(state[3]);
+ res = 10*(fittedDoca - DriftDis);
+ return true;
}
bool GenfitTrack::storeTrack(edm4hep::MutableReconstructedParticle& recParticle,
- int pidType, int ndfCut, double chi2Cut)
+ edm4hep::MutableTrack& track,TVector3& pocaToOrigin_pos,
+ TVector3& pocaToOrigin_mom,TMatrixDSym& pocaToOrigin_cov,
+ // edm4hep::TrackState& trackState,
+ int pidType, int ndfCut, double chi2Cut,
+ int& nFittedDC, int& nFittedSDT, int& ngenfitHit,
+ std::vector<double>& trackL, std::vector<double>& hitMom,
+ std::vector<float>& truthMomEdep,
+ const edm4hep::MCRecoTrackerAssociationCollection* assoHits,
+ std::vector<double>& driftDis,
+ std::vector<double>& FittedDoca,
+ std::vector<double>& Res)
{
+ int id = 0;
+ int fitid = 0;
+ int DCHit = 0;
+ int dcFit = 0;
+ int sdtFit = 0;
+ int SDTHit = 0;
+
+ // getNumRawMeasurements
+ unsigned int nPoints = m_track->getNumPoints();
+ std::cout << " nPoints = " << nPoints << std::endl;
+
+ std::vector<double> hitMomMag;
+
+ while(1){
+ genfit::TrackPoint* point = m_track->getPointWithFitterInfo(id);
+ if(!point)break;
+ id++;
+ genfit::AbsMeasurement* absMea = point->getRawMeasurement();
+
+ // getPoint FitterInfo
+ genfit::AbsFitterInfo* FitterInfo = point->getFitterInfo();
+ genfit::KalmanFitterInfo *KalmanfitterInfo =
+ dynamic_cast<genfit::KalmanFitterInfo*>(FitterInfo);
+
+ unsigned int nNumMea = KalmanfitterInfo->getNumMeasurements();
+
+ bool flag = false;
+ for(int i=0;i<nNumMea;i++)
+ {
+ genfit::MeasurementOnPlane* MeaOnPlane =
+ KalmanfitterInfo->getMeasurementOnPlane(i);
+ double weight = MeaOnPlane->getWeight();
+ if(weight>0.8) flag = true;
+ }
+ if(flag) fitid++;
+
+ genfit::MeasurementOnPlane * Mea =
+ KalmanfitterInfo->getMeasurementOnPlane();
+
+ genfit::PlanarMeasurementSDT* sdtMea =
+ dynamic_cast<genfit::PlanarMeasurementSDT*>(absMea);
+ if(sdtMea){
+ SDTHit++;
+ if(flag) sdtFit++;
+ const edm4hep::TrackerHit* TrackerHit_ = sdtMea->getTrackerHit();
+ track.addToTrackerHits(*TrackerHit_);
+ }else{
+ WireMeasurementDC* dcMea =
+ dynamic_cast<WireMeasurementDC*>(absMea);
+ if(dcMea){
+ const edm4hep::TrackerHit* TrackerHit_ = dcMea->getTrackerHit();
+ track.addToTrackerHits(*TrackerHit_);
+
+
+ DCHit++;
+ if(flag)
+ {
+
+ truthMomEdep.push_back(TrackerHit_->getEDep());
+ double DriftDis,fittedDoca,res = 0;
+ GetDocaRes(dcFit,DriftDis,fittedDoca,res);
+ driftDis.push_back(DriftDis);
+ FittedDoca.push_back(fittedDoca);
+ Res.push_back(res);
+
+ // TMatrixDSym fittedHitCov(6);//cm, GeV
+ // TLorentzVector fittedHitPos;
+ // TVector3 fittedHitMom;
+ // int fittedState=getFittedState(fittedHitPos,fittedHitMom,fittedHitCov,dcFit);
+ // hitMomMag.push_back(fittedHitMom.Mag());
+ //
+ // for(int iSim=0;iSim<assoHits->size();iSim++)
+ // {
+ // //if(*TrackerHit_ == assoHits->at(iSim).getRec())
+ // if(TrackerHit_->getCellID() == assoHits->at(iSim).getRec().getCellID())
+ // {
+ // // std::cout << " if TrackerHit_ = " << TrackerHit_->getCellID() << std::endl;
+ // // std::cout << " if assoHits->at(iSim).getRec() = " << assoHits->at(iSim).getRec().getCellID() << std::endl;
+ // // std::cout << " Sim Mom = " << assoHits->at(iSim).getSim().getMomentum() << std::endl;
+ // // std::cout << " Sim MomMag = " << sqrt(assoHits->at(iSim).getSim().getMomentum()[0]*assoHits->at(iSim).getSim().getMomentum()[0]+assoHits->at(iSim).getSim().getMomentum()[1]*assoHits->at(iSim).getSim().getMomentum()[1]+assoHits->at(iSim).getSim().getMomentum()[2]*assoHits->at(iSim).getSim().getMomentum()[2]) << std::endl;
+ //
+ // break;
+ // }
+ // }
+
+ //std::cout << " i = " << dcFit << "fittedMom = " << fittedHitMom.X() << " " << fittedHitMom.Y() << " " << fittedHitMom.Z() << std::endl;
+ //std::cout << " i = " << dcFit << "fittedMomMag = " << fittedHitMom.Mag() << std::endl;
+ dcFit++;
+ }
+ }
+ }
+ }
+
+ std::cout << " id = " << id << std::endl;
+ std::cout << " fitid = " << fitid << std::endl;
+
+ nFittedDC = dcFit;
+ nFittedSDT = SDTHit;
+
+ std::cout<<"nDC: "<<nFittedDC<<", nSDT: "<<nFittedSDT<<std::endl;
+ std::cout<<"nFittedDC: "<<dcFit<<", nFittedSDT: "<<sdtFit<<std::endl;
+ if(m_debug>0)std::cout<<m_name<<" store track ndfCut "<<ndfCut<<" chi2Cut "
+ <<chi2Cut<<std::endl;
+
/// Get fit status
- const genfit::FitStatus* fitState = m_track->getFitStatus();
+ const genfit::FitStatus* fitState = getFitStatus();
double ndf = fitState->getNdf();
- if(ndf>ndfCut) return false;
+ if(ndf>ndfCut){
+ if(m_debug>0){
+ std::cout<<m_name<<" cut by ndf="<<ndf<<">"<<ndfCut<<std::endl;
+ }
+ return false;
+ }
double chi2 = fitState->getChi2();
- if(chi2>chi2Cut) return false;
+ if(chi2>chi2Cut){
+ if(m_debug>0){
+ std::cout<<m_name<<" cut by chi2="<<chi2<<">"<<chi2Cut<<std::endl;
+ }
+ return false;
+ }
double charge = fitState->getCharge();
int isFitted = fitState->isFitted();
int isConverged = fitState->isFitConverged();
int isConvergedFully = fitState->isFitConvergedFully();
- TMatrixDSym fittedCov;
+
+ TMatrixDSym fittedCov(6);//cm, GeV
TLorentzVector fittedPos;
TVector3 fittedMom;
int fittedState=getFittedState(fittedPos,fittedMom,fittedCov);
- if(m_debug>=2)std::cout<<"fit result: get status OK? pidType "
- <<pidType<<" fittedState==0 " <<(0==fittedState)<<" isFitted "<<isFitted
- <<" isConverged "<<isConverged<<" ndf "<<ndf<<std::endl;
- if((0!=fittedState) || (!isFitted) || !isConvergedFully || (ndf<ndfCut)){
- if(m_debug>=2)std::cout<<" fitting failed"<<std::endl;
+
+ //get dx
+ TVector3 pos;
+ TVector3 mom;
+ pos.SetXYZ(fittedPos.X(),fittedPos.Y(),fittedPos.Z());
+ mom.SetXYZ(fittedMom.X(),fittedMom.Y(),fittedMom.Z());
+ //std::cout << "fitted momx = " << mom.X() << "momy = " << mom.Y() << "momz = " << mom.Z() << " mom.Mag = " << mom.Mag() << std::endl;
+ //std::cout << "fitted posx = " << pos.X() << "posy = " << pos.Y() << "posz = " << pos.Z() << std::endl;
+ double radius = 80.022; // cm
+ const TVector3 linePoint(0,0,0);
+ const TVector3 lineDirection(0,0,1);
+
+ for(int i = 0;i<fitid;i++){
+ int repID=1;
+ radius+=1.741321;
+ double tracklength =
+ extrapolateToCylinder(pos,mom,radius,linePoint,lineDirection,repID);
+ hitMomMag.push_back(mom.Mag());
+ //std::cout << "momx = " << mom.X() << "momy = " << mom.Y() << "momz = " << mom.Z() << " mom.Mag = " << mom.Mag() << std::endl;
+ //std::cout << "posx = " << pos.X() << "posy = " << pos.Y() << "posz = " << pos.Z() << std::endl;
+ trackL.push_back(tracklength);
+ }
+
+ for(int j=0;j<hitMomMag.size()-1;j++)
+ {
+ hitMom.push_back(hitMomMag[j]-hitMomMag[j+1]);
+ //std::cout << " Error Dep = " << hitMomMag[j]-hitMomMag[j+1] << std::endl;
+ }
+
+ if(m_debug>0)std::cout<<m_name<<" fit result: get status OK? pidType "
+ <<pidType<<" fittedState"<<fittedState<<"==0? "<<(0==fittedState)
+ <<" isFitted "<<isFitted
+ <<" isConverged "<<isConverged<<" ndf "<<ndf<<std::endl;
+ if((0!=fittedState)||(!isFitted)||(!isConvergedFully)||(ndf>ndfCut)){
+ if(m_debug>0)std::cout<<m_name<<" fitting FAILED!=========="<<std::endl;
}else{
- if(m_debug>=2)std::cout<<"fit result: Pos("<<
- fittedPos.X()<<" "<<
- fittedPos.Y()<<" "<<
- fittedPos.Z()<<") mom("<<
- fittedMom.X()<<" "<<
- fittedMom.Y()<<" "<<
- fittedMom.Z()<<") p_tot "<<
- fittedMom.Mag()<<" pt "<<
- fittedMom.Perp()<<
- " chi2 "<<chi2<<
- " ndf "<<ndf
- <<std::endl;
- }
- float pos[3]={float(fittedPos.X()),float(fittedPos.Y()),float(fittedPos.Z())};
- float mom[3]={float(fittedMom.X()),float(fittedMom.Y()),float(fittedMom.Z())};
- HelixClass helix;
- helix.Initialize_VP(pos,mom,charge,m_genfitField->getBz(fittedPos.Vect()));
-
-
- /////track status at POCA to origin
- //TVector3 origin(pivot);
- //TVector3 pocaToOrigin_pos(1e9*dd4hep::cm,1e9*dd4hep::cm,1e9*dd4hep::cm);
- //TVector3 pocaToOrigin_mom(1e9*dd4hep::GeV,1e9*dd4hep::GeV,1e9*dd4hep::GeV);
-
- //if(ExtrapolateToPoint(pocaToOrigin_pos,pocaToOrigin_mom,
- // m_track,origin) > 1e6*dd4hep::cm){
- // log<<"extrapolate to origin failed";
- // return false;
- //}
+ if(m_debug>0){
+ const TLorentzVector seedPos=getSeedStatePos();
+ const TVector3 seedMom=getSeedStateMom();
+ std::cout<<m_name<<"===fit-seed:delta pos("
+ <<fittedPos.X()-seedPos.X()<<","
+ <<fittedPos.Y()-seedPos.Y()<<","
+ <<fittedPos.Z()-seedPos.Z()<<") delta mom("
+ <<fittedMom.X()-seedMom.X()<<","
+ <<fittedMom.Y()-seedMom.Y()<<","
+ <<fittedMom.Z()-seedMom.Z()<<")"<<" delta mag "<<fittedMom.Mag()-seedMom.Mag()<<std::endl;
+ std::cout<<m_name<<" seed pos("
+ <<seedPos.X()<<","
+ <<seedPos.Y()<<","
+ <<seedPos.Z()<<") seed mom("
+ <<seedMom.X()<<","
+ <<seedMom.Y()<<","
+ <<seedMom.Z()<<std::endl;
+ std::cout<<m_name<<" fit result: Pos("<<
+ fittedPos.X()<<" "<<
+ fittedPos.Y()<<" "<<
+ fittedPos.Z()<<") mom("<<
+ fittedMom.X()<<" "<<
+ fittedMom.Y()<<" "<<
+ fittedMom.Z()<<") p_tot "<<
+ fittedMom.Mag()<<" pt "<<
+ fittedMom.Perp()<<
+ " chi2 "<<chi2<<
+ " ndf "<<ndf
+ <<std::endl;
+ std::cout<<"fittedCov "<<std::endl;
+ fittedCov.Print();
+ }
+ }
- //new TrackState
- edm4hep::TrackState* trackState = new edm4hep::TrackState();
- trackState->location=0;//edm4hep::AtIP;
- trackState->D0=helix.getD0();
- trackState->phi=helix.getPhi0();
- trackState->omega=helix.getOmega();
- trackState->Z0=helix.getZ0();
- trackState->tanLambda=helix.getTanLambda();
- trackState->referencePoint=helix.getReferencePoint();
-
- // std::array<float,15> covMatrix;
- // int k=0;
- // for(int i=0;i<5;i++){
- // for(int j=0;j<5;j++){
- // if(i<=j) covMatrix[k]=;//FIXME
- // }
- // }
- // trackState.covMatrix=
+ //TVectorD seedState5(5);
+ //getSeedState5(seedState5);
+
+ ///track status at POCA to referencePoint: origin
+ const TVector3 referencePoint(0,0,0);
+ //TVector3
+ pocaToOrigin_pos.SetXYZ(1e9*dd4hep::cm,1e9*dd4hep::cm,1e9*dd4hep::cm);
+ //TVector3
+ pocaToOrigin_mom.SetXYZ(1e9*dd4hep::GeV,1e9*dd4hep::GeV,1e9*dd4hep::GeV);
+ //TMatrixDSym pocaToOrigin_cov;
+ if(extrapolateToPoint(pocaToOrigin_pos,pocaToOrigin_mom,pocaToOrigin_cov,
+ referencePoint) > 1e6*dd4hep::cm){
+ if(m_debug>0)std::cout<<m_name<<" extrapolate to origin failed"<<std::endl;
+ return false;
+ }
- //new Track
- edm4hep::MutableTrack* track = new edm4hep::MutableTrack();
- //track->setType();
- track->setChi2(fitState->getChi2());
- track->setNdf(fitState->getNdf());
- //track->setDEdx();
- //track->setRadiusOfInnermostHit();//FIXME
- //track->addToTrackerHits();
+ //Unit conversion of position and momentum cm->mm
+ pocaToOrigin_pos.SetXYZ(pocaToOrigin_pos.X()/dd4hep::mm,
+ pocaToOrigin_pos.Y()/dd4hep::mm,
+ pocaToOrigin_pos.Z()/dd4hep::mm);
+ pocaToOrigin_mom.SetXYZ(pocaToOrigin_mom.X(),
+ pocaToOrigin_mom.Y(),
+ pocaToOrigin_mom.Z());
+
+ //unit conversion of error matrix
+ TMatrixDSym covMatrix_6=fittedCov;
+ for(int i=0;i<5;i++){
+ covMatrix_6[0][i]=fittedCov[0][i]/dd4hep::mm;//d0 column
+ covMatrix_6[2][i]=fittedCov[1][i]/dd4hep::mm;//omega column
+ covMatrix_6[3][i]=fittedCov[2][i]/dd4hep::mm;//z0 column
+ covMatrix_6[i][0]=fittedCov[i][0]/dd4hep::mm;//d0 row
+ covMatrix_6[i][2]=fittedCov[i][1]/dd4hep::mm;//omega row
+ covMatrix_6[i][3]=fittedCov[i][2]/dd4hep::mm;//z0 row
+ //covMatrix_6[0][i]=fittedCov[0][i]/GenfitUnit::cm*dd4hep::cm;//d0 column
+ //covMatrix_6[2][i]=fittedCov[2][i]*GenfitUnit::cm/dd4hep::cm;//omega column
+ //covMatrix_6[3][i]=fittedCov[3][i]/GenfitUnit::cm*dd4hep::cm;//z0 column
+ //covMatrix_6[i][0]=fittedCov[i][0]/GenfitUnit::cm*dd4hep::cm;//d0 row
+ //covMatrix_6[i][2]=fittedCov[i][2]*GenfitUnit::cm/dd4hep::cm;//omega row
+ //covMatrix_6[i][3]=fittedCov[i][3]/GenfitUnit::cm*dd4hep::cm;//z0 row
+ }
+
+ double Bz=m_genfitField->getBz(referencePoint)/GenfitUnit::tesla;
+
+ if(m_debug>0){
+ std::cout<<m_name<<" fit result poca: Pos"<<std::endl;
+ pocaToOrigin_pos.Print();
+ pocaToOrigin_mom.Print();
+ std::cout<<" chi2 "<<chi2<< " ndf "<<ndf <<std::endl;
+ std::cout<<"fittedCov "<<std::endl;
+ fittedCov.Print();
+ std::cout<<"covMatrix_6 "<<std::endl;
+ covMatrix_6.Print();
+ std::cout<<__LINE__<<" debug Bz tesla"<<Bz<<std::endl;
+ }
+
+ edm4hep::TrackState trackState;
+ CEPC::getTrackStateFromPosMom(trackState,Bz,pocaToOrigin_pos,
+ pocaToOrigin_mom,charge,covMatrix_6);
+
+ trackState.location=0;//edm4hep::AtIP;//FIXME
+ if(m_debug>2){
+ std::cout<<m_name<<" trackState "<<trackState<<std::endl;
+ }
+ track.addToTrackStates(trackState);
+
+ //ngenfitHit = m_genfitHitVec.size();
+ ngenfitHit = nPoints;
+ std::cout << " m_genfitHitVec size = " << m_genfitHitVec.size() << std::endl;
+ // for(long unsigned int i=0; i<m_genfitHitVec.size();i++)
+ // {
+ // GenfitHit * genfitHit = GetHit(i);
+ // edm4hep::TrackerHit* trackHit =
+ // const_cast<edm4hep::TrackerHit*>(genfitHit->getTrackerHit());
+ //
+ // track.addToTrackerHits(*trackHit);
+ // }
+ //track.setType();
+ track.setChi2(chi2);
+ track.setNdf(ndf);
+ //track.setDEdx();
+ //track.setRadiusOfInnermostHit();//FIXME
+ //track.addToTrackerHits();
//new ReconstructedParticle
//recParticle->setType();
//dcRecParticle->setEnergy();
- edm4hep::Vector3f momVec3(helix.getMomentum()[0],
- helix.getMomentum()[1],helix.getMomentum()[2]);
- recParticle.setMomentum(momVec3);
- //recParticle->setReferencePoint(referencePoint);
- recParticle.setCharge(helix.getCharge());
- // recParticle->setMass();
- // recParticle->setCovMatrix();
- // rcRecParticle->setStartVertex();
- //recParticle->addToTracks(track);
+ recParticle.setMomentum(edm4hep::Vector3f(pocaToOrigin_mom.X(),
+ pocaToOrigin_mom.Y(),pocaToOrigin_mom.Z()));
+ recParticle.setReferencePoint(edm4hep::Vector3f(referencePoint.X(),
+ referencePoint.Y(),referencePoint.Z()));
+ recParticle.setCharge(charge);
+ //recParticle->setMass();
+ //recParticle.setCovMatrix(trackState->covMatrix);
+ //recParticle->setStartVertex();
+ recParticle.addToTracks(track);
+ if(m_debug>2){
+ std::cout<<m_name<<" storeTrack trackState "<<trackState<<std::endl;
+ std::cout<<m_name<<" storeTrack track "<<track<<std::endl;
+ }
return true;
}
-void GenfitTrack::pivotToFirstLayer(edm4hep::Vector3d& pos,
- edm4hep::Vector3f& mom, edm4hep::Vector3d& firstPos,
+void GenfitTrack::pivotToFirstLayer(const edm4hep::Vector3d& pos,
+ const edm4hep::Vector3f& mom, edm4hep::Vector3d& firstPos,
edm4hep::Vector3f& firstMom)
{
//FIXME, TODO
@@ -837,3 +1395,246 @@ void GenfitTrack::pivotToFirstLayer(edm4hep::Vector3d& pos,
firstMom=mom;
}
+int GenfitTrack::getDetTypeID(unsigned long long cellID) const
+{
+ UTIL::BitField64 encoder(lcio::ILDCellID0::encoder_string);
+ encoder.setValue(cellID);
+ return encoder[lcio::ILDCellID0::subdet];
+}
+
+void GenfitTrack::getAssoSimTrackerHit(
+ const edm4hep::MCRecoTrackerAssociationCollection*& assoHits,
+ edm4hep::TrackerHit* trackerHit,
+ edm4hep::SimTrackerHit& simTrackerHit) const
+{
+ for(auto assoHit: *assoHits){
+ if(assoHit.getRec()==*trackerHit)
+ {
+ simTrackerHit=assoHit.getSim();
+
+ }
+ }
+}
+
+genfit::AbsTrackRep* GenfitTrack::getRep(int id) const
+{
+ return m_track->getTrackRep(id);
+}
+
+void GenfitTrack::getSeedCov(TMatrixDSym& cov){
+ cov.Zero();
+ for(int i=0;i<3;i++) {
+ double posResolusion=1.;
+ cov(i,i)=posResolusion*posResolusion; //seed position
+ double momResolusion=5.;
+ cov(i+3,i+3)=momResolusion*momResolusion; //seed momentum
+ }
+}
+
+//unit conversion here
+void GenfitTrack::getEndPointsOfWire(int cellID,TVector3& end0,TVector3& end1) const{
+ m_gridDriftChamber->cellposition(cellID,end0,end1);//dd4hep unit
+ end0*=(1./dd4hep::cm*GenfitUnit::cm);
+ end1*=(1./dd4hep::cm*GenfitUnit::cm);
+}
+//unit conversion here
+void GenfitTrack::getTrackFromMCPartile(const edm4hep::MCParticle mcParticle,
+ TVectorD& trackParam, TMatrixDSym& cov) const{
+ TVector3 pos;
+ TVector3 mom;
+ getPosMomFromMCPartile(mcParticle,pos,mom);
+ trackParam[0]=pos[0];
+ trackParam[1]=pos[1];
+ trackParam[2]=pos[2];
+ trackParam[3]=mom[0];
+ trackParam[4]=mom[1];
+ trackParam[5]=mom[2];
+ //cov TODO
+}
+//unit conversion here
+void GenfitTrack::getPosMomFromMCPartile(const edm4hep::MCParticle mcParticle,
+ TVector3& pos,TVector3& mom) const{
+ const edm4hep::Vector3d mcParticleVertex=mcParticle.getVertex();//mm
+ const edm4hep::Vector3f mcParticleMom=mcParticle.getMomentum();//GeV
+ pos[0]=mcParticleVertex.x*GenfitUnit::mm;
+ pos[1]=mcParticleVertex.y*GenfitUnit::mm;
+ pos[2]=mcParticleVertex.z*GenfitUnit::mm;
+ mom[0]=mcParticleMom.x*GenfitUnit::GeV;
+ mom[1]=mcParticleMom.y*GenfitUnit::GeV;
+ mom[2]=mcParticleMom.z*GenfitUnit::GeV;
+ if(m_debug>2){
+ std::cout<<"getPosMomFromTrackState pos("<<pos.X()<<" "<<pos.Y()
+ <<" "<<pos.Z();
+ std::cout<<") "<<" ("<<mom.X()<<" "<<mom.Y()<<" "<<mom.Z()<<")"
+ <<" mom "<<mom.Mag()<<" pt "<<mom.Perp()<<std::endl;
+ }
+}
+
+//unit conversion here
+void GenfitTrack::getTrackFromEDMTrack(const edm4hep::Track& edm4HepTrack,
+ double& charge, TVectorD& trackParam, TMatrixDSym& cov) const{
+ TVector3 pos;
+ TVector3 mom;
+ double Bz=m_genfitField->getBz(TVector3{0.,0.,0.})/GenfitUnit::tesla;
+ double charge_double;
+ CEPC::getPosMomFromTrackState(edm4HepTrack.getTrackStates(0),Bz,pos,mom,charge_double,cov);
+
+ //std::cout<<__LINE__<<" Bz "<<Bz<<" charge "<<charge_double<<std::endl;
+ //pos.Print();
+ //mom.Print();
+ charge=(int) charge_double;
+ trackParam[0]=pos[0]*GenfitUnit::mm;
+ trackParam[1]=pos[1]*GenfitUnit::mm;
+ trackParam[2]=pos[2]*GenfitUnit::mm;
+ trackParam[3]=mom[0]*GenfitUnit::GeV;
+ trackParam[4]=mom[1]*GenfitUnit::GeV;
+ trackParam[5]=mom[2]*GenfitUnit::GeV;
+
+ //cov unit conversion
+ for(int i=0;i<6;i++){
+ for(int j=0;j<6;j++){
+ cov(i,j) = cov(i,j)*GenfitUnit::mm;
+ }
+ }
+}
+
+//to genfit unit
+void GenfitTrack::getISurfaceOUV(const dd4hep::rec::ISurface* iSurface,TVector3& o,
+ TVector3& u,TVector3& v, double& lengthU,double& lengthV){
+ o.SetXYZ(iSurface->origin().x()/dd4hep::mm*GenfitUnit::mm,
+ iSurface->origin().y()/dd4hep::mm*GenfitUnit::mm,
+ iSurface->origin().z()/dd4hep::mm*GenfitUnit::mm);
+ u.SetXYZ(iSurface->u().x()/dd4hep::mm*GenfitUnit::mm,
+ iSurface->u().y()/dd4hep::mm*GenfitUnit::mm,
+ iSurface->u().z()/dd4hep::mm*GenfitUnit::mm);
+ v.SetXYZ(iSurface->v().x()/dd4hep::mm*GenfitUnit::mm,
+ iSurface->v().y()/dd4hep::mm*GenfitUnit::mm,
+ iSurface->v().z()/dd4hep::mm*GenfitUnit::mm);
+
+ lengthU=iSurface->length_along_u();
+ lengthV=iSurface->length_along_v();
+}
+
+void GenfitTrack::getMeasurementAndCov(edm4hep::TrackerHit* hit,TVector3& pos,TMatrixDSym& cov){
+
+ pos.SetXYZ(hit->getPosition().x*GenfitUnit::mm,
+ hit->getPosition().y*GenfitUnit::mm,
+ hit->getPosition().z*GenfitUnit::mm);
+}
+
+
+int GenfitTrack::getSigmas(int cellID,std::vector<float> sigmaUVec,
+ std::vector<float> sigmaVVec,float& sigmaU,float& sigmaV)const{
+ int detTypeID=getDetTypeID(cellID);
+ int sigmaUID=0;
+ int sigmaVID=0;
+ int layer=0;
+ dd4hep::DDSegmentation::BitFieldCoder* decoder;
+ if(m_geomSvc->lcdd()->constant<int>("DetID_DC")==detTypeID){
+ sigmaUID=0;
+ sigmaVID=0;
+ }else if(detTypeID==lcio::ILDDetID::VXD){
+ decoder=m_geomSvc->getDecoder("VXDCollection");//FIXME
+ layer=decoder->get(cellID,"layer");//FIXME
+ sigmaUID=layer+1;
+ sigmaVID=layer+1;
+ }else if(detTypeID==lcio::ILDDetID::SIT){
+ sigmaUID=7;
+ sigmaVID=7;
+ }else if(detTypeID==lcio::ILDDetID::SET){
+ sigmaUID=8;
+ sigmaVID=8;
+ }else if(detTypeID==lcio::ILDDetID::FTD){
+ decoder=m_geomSvc->getDecoder("FTDCollection");//FIXME
+ layer=decoder->get(cellID,"layer");//FIXME
+ sigmaUID=layer+9;
+ sigmaVID=layer+9;
+ }else{
+ if(m_debug>=0) std::cout<<m_name<<" Error: no detType!"<<std::endl;
+ }
+ if(m_debug){
+ std::cout<<"sigmaUID "<<sigmaUID<<" sigmaVID "<<sigmaVID<<std::endl;
+ //std::cout<<"pos "<<pos_smeared[0]<<" "<<pos_smeared[1]<<" "<<pos_smeared[2]<<std::endl;
+ //std::cout<<"angle "<<atan2(pos_smeared[1],pos_smeared[0])<<std::endl;
+ std::cout<<"sigmaU "<<sigmaUVec[sigmaUID]*GenfitUnit::mm
+ <<" sigmaV "<<sigmaVVec[sigmaVID]*GenfitUnit::mm<<std::endl;
+ }
+ sigmaU=sigmaUVec[sigmaUID]*GenfitUnit::mm;
+ sigmaV=sigmaVVec[sigmaVID]*GenfitUnit::mm;
+ return sigmaUID;
+}
+
+bool GenfitTrack::isCDCHit(edm4hep::TrackerHit* hit){
+ return m_geomSvc->lcdd()->constant<int>("DetID_DC")==
+ getDetTypeID(hit->getCellID());
+}
+
+void GenfitTrack::getSortedTrackerHits(
+ std::vector<edm4hep::TrackerHit*>& trackerHits,
+ const edm4hep::MCRecoTrackerAssociationCollection* assoHits,
+ std::vector<edm4hep::TrackerHit*>& sortedDCTrackerHits,
+ int sortMethod){
+
+ std::vector<std::pair<double,edm4hep::TrackerHit*> > sortedDCTrackerHitPair;
+ for(auto trackerHit : trackerHits){
+ edm4hep::TrackerHit* thisHit = trackerHit;
+ if(!isCDCHit(thisHit))continue;//skip non-DC trackerHit
+
+ edm4hep::SimTrackerHit simTrackerHitAsso;
+ getAssoSimTrackerHit(assoHits,trackerHit,simTrackerHitAsso);
+ double time=simTrackerHitAsso.getTime();
+ if(0==sortMethod){
+ //by time
+ sortedDCTrackerHitPair.push_back(std::make_pair(time,trackerHit));
+ if(m_debug>0){ std::cout<<"sorted DC digi by time"<<std::endl;}
+ }else if(1==sortMethod){
+ //by layer
+ sortedDCTrackerHitPair.push_back(std::make_pair(
+ m_decoderDC->get(trackerHit->getCellID(),"layer"),trackerHit));
+ if(m_debug>0){ std::cout<<"sorted DC digi by layer"<<std::endl;}
+ }else{
+ sortedDCTrackerHits.push_back(trackerHit);
+ }
+ }
+ if(0==sortMethod || 1==sortMethod){
+ std::sort(sortedDCTrackerHitPair.begin(),sortedDCTrackerHitPair.end(),
+ sortDCDigi);
+ for(auto trackerHit:sortedDCTrackerHitPair){
+ sortedDCTrackerHits.push_back(trackerHit.second);
+ }
+ }
+ if(m_debug>0){
+ std::cout<<"trackerHits on track after sort\n";
+ for(auto trackerHit:sortedDCTrackerHits){
+ std::cout<<trackerHit->getCellID()<<"("<<std::setw(2)
+ <<m_decoderDC->get(trackerHit->getCellID(),"layer")
+ <<","<<std::setw(3)
+ <<m_decoderDC->get(trackerHit->getCellID(),"cellID")<<")\n";
+ }
+ std::cout<<"\n"; std::cout.unsetf(std::ios_base::floatfield);
+ }
+ return;
+}//end of getSortedTrackerHits
+
+//make a genfit hit and do hits selection
+GenfitHit* GenfitTrack::makeAGenfitHit(edm4hep::TrackerHit* trackerHit,
+ edm4hep::SimTrackerHit* simTrackerHitAsso,
+ double sigma,bool truthAmbig,double skipCorner,double skipNear){
+
+ //TODO truthAmbig
+ double driftVelocity=40.;//FIXME
+ GenfitHit* genfitHit=new GenfitHit(trackerHit,simTrackerHitAsso,m_decoderDC,
+ m_gridDriftChamber,driftVelocity,sigma*GenfitUnit::mm);
+ //skip corner hit
+ if(fabs(genfitHit->getDriftDistance())>skipCorner){
+ if(m_debug) std::cout<<" skip hit dd > skipCorner"<<std::endl;
+ delete genfitHit;
+ return nullptr;
+ }
+ if(genfitHit->getDriftDistance()<skipNear){
+ if(m_debug) std::cout<<" skip hit dd < skipCorner"<<std::endl;
+ delete genfitHit;
+ return nullptr;
+ }
+ return genfitHit;
+}
diff --git a/Reconstruction/RecGenfitAlg/src/GenfitTrack.h b/Reconstruction/RecGenfitAlg/src/GenfitTrack.h
index bcff947944fbca3966818333c6496a6217f1d42c..bfcb7a99524e826b3f219a12f1ade4e9a6da7e80 100644
--- a/Reconstruction/RecGenfitAlg/src/GenfitTrack.h
+++ b/Reconstruction/RecGenfitAlg/src/GenfitTrack.h
@@ -11,8 +11,6 @@
///
/// Authors:
/// Zhang Yao (zhangyao@ihep.ac.cn)
-/// Y.Fujii (yfujii@ihep.ac.cn)
-/// Yohei Nakatsugawa (yohei@ihep.ac.cn)
///
//////////////////////////////////////////////////////////////////
@@ -20,15 +18,23 @@
#define RECGENFITALG_GENFITTRACK_H
#include "GenfitFitter.h"
+#include "GenfitHit.h"
//ROOT
+#include "TVector3.h"
#include "TVectorD.h"
#include "TMatrixDSym.h"
+#include "TLorentzVector.h"
+
+//Gaudi
+#include "GaudiKernel/SmartIF.h"
//STL
#include <vector>
class TLorentzVector;
+class IGeomSvc;
+class WireMeasurementDC;
namespace genfit{
class Track;
@@ -44,101 +50,132 @@ namespace edm4hep{
class MutableReconstructedParticle;
class MCRecoTrackerAssociationCollection;
class Track;
- class Track;
+ class MutableTrack;
+ class TrackCollection;
class TrackerHit;
+ class SimTrackerHit;
class Vector3d;
class Vector3f;
+ class TrackerHitCollection;
}
namespace dd4hep {
namespace DDSegmentation{
class GridDriftChamber;
+ class BitFieldCoder;
+ }
+ namespace rec{
+ class ISurface;
}
}
class GenfitTrack {
friend int GenfitFitter::processTrack(
- GenfitTrack* track, bool resort);
+ GenfitTrack* track, bool resort=true);
friend int GenfitFitter::processTrackWithRep(
- GenfitTrack* track, int repID, bool resort);
-
- friend double GenfitFitter::extrapolateToHit(TVector3& poca,
- TVector3& pocaDir,
- TVector3& pocaOnWire, double& doca, const GenfitTrack* track,
- TVector3 pos, TVector3 mom, TVector3 end0, TVector3 end1, int debug,
- int repID, bool stopAtBoundary, bool calcJacobianNoise);
-
- friend double GenfitFitter::extrapolateToCylinder(TVector3& pos,
- TVector3& mom,
- GenfitTrack* track, double radius, const TVector3 linePoint,
- const TVector3 lineDirection, int hitID, int repID,
- bool stopAtBoundary, bool calcJacobianNoise);
-
- friend double GenfitFitter::extrapolateToPoint(TVector3& pos, TVector3& mom,
- const GenfitTrack* genfitTrack, const TVector3& point, int repID,
- bool stopAtBoundary, bool calcJacobianNoise) const;
+ GenfitTrack* track, int repID, bool resort=true);
public:
GenfitTrack(const GenfitField* field,
const dd4hep::DDSegmentation::GridDriftChamber* seg,
+ SmartIF<IGeomSvc> geom,
const char* name="GenfitTrack");
virtual ~GenfitTrack();
- /// Add a Genfit track
- virtual bool createGenfitTrack(int pdgType,int charge,TLorentzVector pos, TVector3 mom,
- TMatrixDSym covM);
- //virtual bool createGenfitTrack(TLorentzVector posInit,TVector3 momInit,
- //TMatrixDSym covMInit);
-
///Create genfit track from MCParticle
bool createGenfitTrackFromMCParticle(int pidTyep,const edm4hep::MCParticle&
mcParticle, double eventStartTime=0.);
- bool createGenfitTrackFromEDM4HepTrack(int pidType, edm4hep::Track track,
- double eventStartTime);
+ bool createGenfitTrackFromEDM4HepTrack(int pidType,const edm4hep::Track& track,
+ double eventStartTime,bool isUseCovTrack);
- // /// Prepare a hit list, return number of hits on track
- // int PrepareHits();//TODO
+ /// ---------Add measurements---------
+ ///Add one space point measurement, return number of hits on track
+ virtual bool addSpacePointMeasurement(const TVector3&,std::vector<float>
+ sigmaU,std::vector<float> sigmaV,int cellID,int hitID);
- /// Add a space point measurement, return number of hits on track
- bool addSpacePointTrakerHit(edm4hep::TrackerHit hit, int hitID);
+ ///Add silicon space points from edm4hep::track
+ int addSpacePointsSi(const edm4hep::Track& track,
+ std::vector<float> sigmaU,std::vector<float> sigmaV);
- /// Add a space point measurement, return number of hits on track
- virtual bool addSpacePointMeasurement(const TVectorD&, double,
- int detID=-1, int hitID=-1, bool smear=false);
+ ///Add drift chamber space points from edm4hep::track
+ int addSpacePointsDC(const edm4hep::Track& track,
+ const edm4hep::MCRecoTrackerAssociationCollection* assoHits,
+ std::vector<float> sigmaU,std::vector<float> sigmaV);
- /// Add a WireMeasurement with MC truth position smeared by sigma
- virtual void addWireMeasurement(double driftDistance,
- double driftDistanceError, const TVector3& endPoint1,
- const TVector3& endPoint2, int lrAmbig, int detID, int hitID);
- /// Add a WireMeasurement with DC digi
- virtual bool addWireMeasurementOnTrack(edm4hep::Track track, double sigma);
+ ///Add WireMeasurements of hits on track
+ virtual int addWireMeasurementsOnTrack(edm4hep::Track& track,float sigma,
+ const edm4hep::MCRecoTrackerAssociationCollection* assoHits,
+ int sortMethod,bool truthAmbig,float skipCorner, float skipNear);
- ///Add space point from truth to track
- int addSimTrackerHits( edm4hep::Track track,
- const edm4hep::MCRecoTrackerAssociationCollection* assoHits,
- float sigma,bool smear=false);// float nSigmaSelection
+ ///Add WireMeasurements of hits on track from hit selection
+ virtual int addWireMeasurementsFromList(std::vector<edm4hep::TrackerHit*>& hits,float sigma,
+ const edm4hep::MCRecoTrackerAssociationCollection* assoHits,
+ int sortMethod,bool truthAmbig,float skipCorner, float skipNear);
+
+ ///Add one silicon hits
+ bool addSiliconMeasurement(edm4hep::TrackerHit* hit,
+ float sigmaU,float sigmaV,int cellID,int hitID);
+
+ ///Add silicon measurements, return number of hits on track
+ int addSiliconMeasurements(edm4hep::Track& track,
+ std::vector<float> sigmaU,std::vector<float> sigmaV);
+
+ bool debugDistance(const edm4hep::TrackerHitCollection* dCDigiCol,
+ int& nFittedDC,int& nFittedSDT,int& ngenfitHit,
+ std::vector<double>& smearDistance,
+ std::vector<double>& truthDistance,double driftVelocity);
+ bool GetDocaRes(int hitID, double& DriftDis,double& fittedDoca,
+ double& res,int repID=0, bool biased=true) const;
///Store track to ReconstructedParticle
- bool storeTrack(edm4hep::MutableReconstructedParticle& dcRecParticle,int pidType,
- int ndfCut=1e9, double chi2Cut=1.e9);
+ bool storeTrack(edm4hep::MutableReconstructedParticle& recParticle,
+ edm4hep::MutableTrack& track,
+ TVector3& pocaToOrigin_pos,
+ TVector3& pocaToOrigin_mom,
+ TMatrixDSym& pocaToOrigin_cov,
+ int pidType, int ndfCut, double chi2Cut,
+ int& nFittedDC, int& nFittedSDT,int& ngenfitHit,
+ std::vector<double>& trackL, std::vector<double>& hitMom,
+ std::vector<float>& truthMomEdep,
+ const edm4hep::MCRecoTrackerAssociationCollection* assoHits,
+ std::vector<double>& driftDis,
+ std::vector<double>& FittedDoca,
+ std::vector<double>& Res);
///A tool to convert track to the first layer of DC
- void pivotToFirstLayer(edm4hep::Vector3d& pos,edm4hep::Vector3f& mom,
- edm4hep::Vector3d& firstPos, edm4hep::Vector3f& firstMom);
+ void pivotToFirstLayer(const edm4hep::Vector3d& pos,
+ const edm4hep::Vector3f& mom, edm4hep::Vector3d& firstPos,
+ edm4hep::Vector3f& firstMom);
/// Copy a track to event
//void CopyATrack()const;
- ///Extrapolate to Hit
- double extrapolateToHit( TVector3& poca, TVector3& pocaDir,
- TVector3& pocaOnWire, double& doca, TVector3 pos, TVector3 mom,
- TVector3 end0,//one end of the hit wire
- TVector3 end1,//the orhter end of the hit wire
- int debug,
- int repID,
- bool stopAtBoundary,
- bool calcJacobianNoise) const;
+ //return dd4hep unit
+ double extrapolateToHit(TVector3& poca, TVector3& pocaDir,
+ TVector3& pocaOnWire, double& doca,edm4hep::MCParticle mcParticle,
+ int cellID, int repID, bool stopAtBoundary, bool calcJacobianNoise)const;
+ /// Extrapolate the track to the point
+ /// Output: pos and mom of POCA point to point
+ /// Input: genfitTrack,point,repID,stopAtBoundary and calcAverageState
+ /// repID same with pidType
+// double extrapolateToPoint(TVector3& pos, TVector3& mom,TMatrixDSym& cov,
+// const TVector3& point, int repID=0, bool stopAtBoundary = false,
+// bool calcJacobianNoise = true) const;
+
+ double extrapolateToPoint(TVector3& pos, TVector3& mom, TMatrixDSym& cov,
+ const TVector3& point, int repID=0,
+ bool stopAtBoundary = false, bool calcJacobianNoise = true) const;
+
+ /// Extrapolate the track to the cyliner at fixed raidus
+ /// Output: pos and mom at fixed radius
+ /// Input: genfitTrack, radius of cylinder at center of the origin,
+ /// repID, stopAtBoundary and calcAverageState
+ double extrapolateToCylinder(TVector3& pos, TVector3& mom,
+ double radius, const TVector3 linePoint,
+ const TVector3 lineDirection, int hitID =0, int repID=0,
+ bool stopAtBoundary=false, bool calcJacobianNoise=true);
+
bool getPosMomCovMOP(int hitID, TLorentzVector& pos, TVector3& mom,
TMatrixDSym& cov, int repID=0) const;
@@ -148,11 +185,12 @@ class GenfitTrack {
const TVector3 getSeedStateMom() const;
void getSeedStateMom(TLorentzVector& pos, TVector3& mom) const;
unsigned int getNumPoints() const;
+ unsigned int getNumPointsDet(int cellID) const;
/// get the fitted track status
const genfit::FitStatus* getFitStatus(int repID=0) const;
int getFittedState(TLorentzVector& pos, TVector3& mom, TMatrixDSym& cov,
- int repID=0, bool biased=true) const;
+ int trackPointId=0, int repID=0, bool biased=true) const;
int getNumPointsWithFittedInfo(int repID=0) const;
bool getFirstPointWithFittedInfo(int repID=0) const;
bool fitSuccess(int repID) const;
@@ -169,7 +207,9 @@ class GenfitTrack {
void print(TLorentzVector pos, TVector3 mom, const char* comment="") const;
/// set and get debug level
- void setDebug(int debug);
+ void setDebug(int debug){m_debug=debug;}
+ void setDebugGenfit(int debug);
+ void setDebugLocal(int debug);
int getDebug(void) const { return m_debug;}
/// get name of this object
@@ -177,30 +217,67 @@ class GenfitTrack {
genfit::Track* getTrack() const{return m_track;}
/// Add a track representation
- genfit::RKTrackRep* addTrackRep(int pdg);
+ genfit::RKTrackRep* addTrackRep(int pdgType,int charge);
+ /// Get a hit according to index
+ GenfitHit* GetHit(long unsigned int i) const;
protected:
//genfit::Track* getTrack() {return m_track;}
private:
+ /// ---------Add a Genfit track-------
+ bool createGenfitTrack(int pdgType,int charge,
+ TVectorD trackParam, TMatrixDSym covMInit_6);
+
+ int getDetTypeID(unsigned long long cellID) const;
const char* m_name;
///Note! private functions are using genfit unit, cm and MeV
- genfit::AbsTrackRep* getRep(int id=0) const {return m_reps[id];}
+ genfit::AbsTrackRep* getRep(int id=0) const;
bool getMOP(int hitID, genfit::MeasuredStateOnPlane& mop,
genfit::AbsTrackRep* trackRep=nullptr) const;
+ const dd4hep::rec::ISurface* getISurface(edm4hep::TrackerHit* hit);
+ void getSeedCov(TMatrixDSym& cov);
+ void getAssoSimTrackerHit(
+ const edm4hep::MCRecoTrackerAssociationCollection*& assoHits,
+ edm4hep::TrackerHit* trackerHit,
+ edm4hep::SimTrackerHit& simTrackerHit) const;
+ void getEndPointsOfWire(int cellID,TVector3& end0,TVector3& end1)const;
+ void getTrackFromEDMTrack(const edm4hep::Track& edm4HepTrack,
+ double& charge, TVectorD& trackParam, TMatrixDSym& cov) const;
+ void getTrackFromMCPartile(const edm4hep::MCParticle mcParticle,
+ TVectorD& trackParam, TMatrixDSym& cov) const;
+ void getPosMomFromMCPartile(const edm4hep::MCParticle mcParticle,
+ TVector3& pos,TVector3& mom) const;
+ void clearGenfitHitVec();
+ void getISurfaceOUV(const dd4hep::rec::ISurface* iSurface,TVector3& o,
+ TVector3& u,TVector3& v,double& lengthU,double& lengthV);
+ void getMeasurementAndCov(edm4hep::TrackerHit* hit,TVector3& pos,TMatrixDSym& cov);
+ int getSigmas(int cellID,std::vector<float> sigmaUVec,
+ std::vector<float> sigmaVVec,float& sigmaU,float& sigmaV)const;
+ bool isCDCHit(edm4hep::TrackerHit* hit);
+ GenfitHit* makeAGenfitHit(edm4hep::TrackerHit* trackerHit,
+ edm4hep::SimTrackerHit* simTrackerHitAsso,
+ double sigma,bool truthAmbig,double skipCorner,double skipNear);
+ void getSortedTrackerHits(std::vector<edm4hep::TrackerHit*>& hits,
+ const edm4hep::MCRecoTrackerAssociationCollection* assoHits,
+ std::vector<edm4hep::TrackerHit*>& sortedDCTrackerHits,
+ int sortMethod);
genfit::Track* m_track;/// track
- std::vector<genfit::AbsTrackRep*> m_reps;/// track repesentations
int m_debug;/// debug level
+ int m_debugLocal;/// debug level local
+ SmartIF<IGeomSvc> m_geomSvc;
const GenfitField* m_genfitField;//pointer to genfit field
const dd4hep::DDSegmentation::GridDriftChamber* m_gridDriftChamber;
+ const dd4hep::DDSegmentation::BitFieldCoder* m_decoderDC;
static const int s_PDG[2][5];
-};
+ std::vector<GenfitHit*> m_genfitHitVec;
+};
#endif
diff --git a/Reconstruction/Tracking/include/Tracking/TrackingHelper.h b/Reconstruction/Tracking/include/Tracking/TrackingHelper.h
index beed38e6231ae56a7dd6e45c13f9df4eb98bc0a4..c12f55174033068d259fa520eed26bd16fd47eef 100644
--- a/Reconstruction/Tracking/include/Tracking/TrackingHelper.h
+++ b/Reconstruction/Tracking/include/Tracking/TrackingHelper.h
@@ -62,5 +62,4 @@ inline int getLayer(const edm4hep::TrackerHit hit) {
return layer;
}
-
#endif
diff --git a/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp b/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp
index 29e714e3b4ec2ddc5281c8f08700a117c2439102..c87ebc39a212dfaa2ac67d6780f580aa6a17728e 100644
--- a/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp
+++ b/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp
@@ -7,6 +7,7 @@
#include "DD4hep/IDDescriptor.h"
#include "DD4hep/Plugins.h"
#include "DD4hep/DD4hepUnits.h"
+#include "UTIL/ILDConf.h"
//external
#include "CLHEP/Random/RandGauss.h"
@@ -25,37 +26,71 @@ TruthTrackerAlg::TruthTrackerAlg(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc),m_dd4hep(nullptr),m_gridDriftChamber(nullptr),
m_decoder(nullptr)
{
+ declareProperty("NoiseSimHitsCollection", w_SimNoiseHCol,
+ "Handle of the output simulation DC Noise hits collection");
+ declareProperty("NoiseDCHitsCollection", w_NoiseHitCol,
+ "Handle of the input DC Noise hits collection");
+ declareProperty("NoiseDCHitAssociationCollection", w_NoiseAssociationCol,
+ "Handle of the simNoiseDCHit and noiseDC Noise hits collection");
declareProperty("MCParticle", m_mcParticleCol,
"Handle of the input MCParticle collection");
+ declareProperty("DriftChamberHitsCollection", m_DCSimTrackerHitCol,
+ "Handle of DC SimTrackerHit collection");
declareProperty("DigiDCHitCollection", m_DCDigiCol,
"Handle of DC digi(TrackerHit) collection");
declareProperty("DCHitAssociationCollection", m_DCHitAssociationCol,
"Handle of association collection");
declareProperty("DCTrackCollection", m_DCTrackCol,
"Handle of DC track collection");
- declareProperty("SiSubsetTrackCollection", m_siSubsetTrackCol,
+ declareProperty("SubsetTracks", m_siSubsetTrackCol,
"Handle of silicon subset track collection");
declareProperty("SDTTrackCollection", m_SDTTrackCol,
"Handle of SDT track collection");
- declareProperty("DCRecParticleCollection", m_DCRecParticleCol,
- "Handle of drift chamber reconstructed particle collection");
- declareProperty("DCRecParticleAssociationCollection",
- m_DCRecParticleAssociationCol,
- "Handle of drift chamber reconstructed particle collection");
+ declareProperty("VXDTrackerHits", m_VXDTrackerHits,
+ "Handle of input VXD tracker hit collection");
+ declareProperty("SITTrackerHits", m_SITTrackerHits,
+ "Handle of input SIT tracker hit collection");
+ declareProperty("SETTrackerHits", m_SETTrackerHits,
+ "Handle of input SET tracker hit collection");
+ declareProperty("FTDTrackerHits", m_FTDTrackerHits,
+ "Handle of input FTD tracker hit collection");
+ declareProperty("SITSpacePoints", m_SITSpacePointCol,
+ "Handle of input SIT hit collection");
+// declareProperty("SETSpacePoints", m_SETSpacePointCol,
+// "Handle of input SET hit collection");
+ declareProperty("FTDSpacePoints", m_FTDSpacePointCol,
+ "Handle of input FTD hit collection");
+ declareProperty("VXDCollection", m_VXDCollection,
+ "Handle of input VXD hit collection");
+ declareProperty("SITCollection", m_SITCollection,
+ "Handle of input SIT hit collection");
+ declareProperty("SETCollection", m_SETCollection,
+ "Handle of input SET hit collection");
+ declareProperty("FTDCollection", m_FTDCollection,
+ "Handle of input FTD hit collection");
+ declareProperty("TruthTrackerHitCollection", m_truthTrackerHitCol,
+ "Handle of output truth TrackerHit collection");
+ declareProperty("SmearTrackerHitCollection", m_SmeartruthTrackerHitCol,
+ "Handle of output smear truth TrackerHit collection");
+ declareProperty("SmearSimHitsCollection", w_SimSmearHCol,
+ "Handle of output smear SimTrackerHit collection");
+ declareProperty("SmearDCHitAssociationCollection", w_SmearAssociationCol,
+ "Handle of output smear simulationhit and TrackerHit collection");
}
+
StatusCode TruthTrackerAlg::initialize()
{
///Get geometry
m_geomSvc=service<IGeomSvc>("GeomSvc");
if (!m_geomSvc) {
- error() << "Failed to get GeomSvc." << endmsg;
+ error()<<"Failed to get GeomSvc."<<endmsg;
return StatusCode::FAILURE;
}
///Get Detector
m_dd4hep=m_geomSvc->lcdd();
if (nullptr==m_dd4hep) {
- error() << "Failed to get dd4hep::Detector." << endmsg;
+ error()<<"Failed to get dd4hep::Detector."<<endmsg;
return StatusCode::FAILURE;
}
//Get Field
@@ -66,44 +101,125 @@ StatusCode TruthTrackerAlg::initialize()
m_gridDriftChamber=dynamic_cast<dd4hep::DDSegmentation::GridDriftChamber*>
(readout.segmentation().segmentation());
if(nullptr==m_gridDriftChamber){
- error() << "Failed to get the GridDriftChamber" << endmsg;
+ error()<<"Failed to get the GridDriftChamber"<<endmsg;
return StatusCode::FAILURE;
}
///Get Decoder
m_decoder = m_geomSvc->getDecoder(m_readout_name);
if (nullptr==m_decoder) {
- error() << "Failed to get the decoder" << endmsg;
+ error()<<"Failed to get the decoder"<<endmsg;
return StatusCode::FAILURE;
}
+ m_tuple=nullptr;
+ ///book tuple
+ if(m_hist){
+ NTuplePtr nt(ntupleSvc(), "TruthTrackerAlg/truthTrackerAlg");
+ if(nt){
+ m_tuple=nt;
+ }else{
+ m_tuple=ntupleSvc()->book("TruthTrackerAlg/truthTrackerAlg",
+ CLID_ColumnWiseTuple,"TruthTrackerAlg");
+ if(m_tuple){
+ StatusCode sc;
+ sc=m_tuple->addItem("run",m_run);
+ sc=m_tuple->addItem("evt",m_evt);
+ sc=m_tuple->addItem("siMom",3,m_siMom);
+ sc=m_tuple->addItem("siPos",3,m_siPos);
+ sc=m_tuple->addItem("mcMom",3,m_mcMom);
+
+ sc=m_tuple->addItem("mcPos",3,m_mcPos);
+ sc=m_tuple->addItem("nDCTrackHit",m_nDCTrackHit,0,100000);
+ sc=m_tuple->addItem("DriftDistance",m_nDCTrackHit,m_DriftDistance);
+ sc=m_tuple->addItem("nSmearDCTrackHit",m_nSmearDCTrackHit,0,100000);
+ sc=m_tuple->addItem("SmearDriftDistance",m_nSmearDCTrackHit,m_SmearDriftDistance);
+ sc=m_tuple->addItem("nNoiseDCTrackHit",m_nNoiseDCTrackHit,0,100000);
+ sc=m_tuple->addItem("NoiseDriftDistance",m_nNoiseDCTrackHit,m_NoiseDriftDistance);
+
+ sc=m_tuple->addItem("nSimTrackerHitVXD",m_nSimTrackerHitVXD);
+ sc=m_tuple->addItem("nSimTrackerHitSIT",m_nSimTrackerHitSIT);
+ sc=m_tuple->addItem("nSimTrackerHitSET",m_nSimTrackerHitSET);
+ sc=m_tuple->addItem("nSimTrackerHitFTD",m_nSimTrackerHitFTD);
+ sc=m_tuple->addItem("nSimTrackerHitDC",m_nSimTrackerHitDC);
+ sc=m_tuple->addItem("nTrackerHitVXD",m_nTrackerHitVXD);
+ sc=m_tuple->addItem("nTrackerHitSIT",m_nTrackerHitSIT);
+ sc=m_tuple->addItem("nTrackerHitSET",m_nTrackerHitSET);
+ sc=m_tuple->addItem("nTrackerHitFTD",m_nTrackerHitFTD);
+ sc=m_tuple->addItem("nTrackerHitDC",m_nTrackerHitDC);
+ sc=m_tuple->addItem("nSpacePointSIT",m_nSpacePointSIT);
+ sc=m_tuple->addItem("nSpacePointSET",m_nSpacePointSET);
+ sc=m_tuple->addItem("nSpacePointFTD",m_nSpacePointFTD);
+ sc=m_tuple->addItem("nHitOnSiTkXVD",m_nHitOnSiTkVXD);
+ sc=m_tuple->addItem("nHitOnSiTkSIT",m_nHitOnSiTkSIT);
+ sc=m_tuple->addItem("nHitOnSiTkSET",m_nHitOnSiTkSET);
+ sc=m_tuple->addItem("nHitOnSiTkFTD",m_nHitOnSiTkFTD);
+ sc=m_tuple->addItem("nHitOnSdtTkVXD",m_nHitOnSdtTkVXD);
+ sc=m_tuple->addItem("nHitOnSdtTkSIT",m_nHitOnSdtTkSIT);
+ sc=m_tuple->addItem("nHitOnSdtTkSET",m_nHitOnSdtTkSET);
+ sc=m_tuple->addItem("nHitOnSdtTkFTD",m_nHitOnSdtTkFTD);
+ sc=m_tuple->addItem("nHitOnSdtTkDC",m_nHitOnSdtTkDC);
+ sc=m_tuple->addItem("nHitSdt",m_nHitOnSdtTk);
+ sc=m_tuple->addItem("nNoiseOnSdt",m_nNoiseOnSdtTk);
+ }
+ }
+ }
return GaudiAlgorithm::initialize();
}
StatusCode TruthTrackerAlg::execute()
{
+ info()<<"In execute()"<<endmsg;
+
+ ///Output DC Track collection
+ edm4hep::TrackCollection* dcTrackCol=m_DCTrackCol.createAndPut();
+
+ ///Output SDT Track collection
+ edm4hep::TrackCollection* sdtTkCol=m_SDTTrackCol.createAndPut();
+
+ ///Output Hit collection
+ auto truthTrackerHitCol=m_truthTrackerHitCol.createAndPut();
+
///Retrieve MC particle(s)
const edm4hep::MCParticleCollection* mcParticleCol=nullptr;
mcParticleCol=m_mcParticleCol.get();
+ //if(m_mcParticleCol.exist()){mcParticleCol=m_mcParticleCol.get();}
if(nullptr==mcParticleCol){
debug()<<"MCParticleCollection not found"<<endmsg;
return StatusCode::SUCCESS;
}
///Retrieve DC digi
const edm4hep::TrackerHitCollection* digiDCHitsCol=nullptr;
- digiDCHitsCol=m_DCDigiCol.get();//FIXME DEBUG
- if(nullptr==digiDCHitsCol){
- debug()<<"TrackerHitCollection not found"<<endmsg;
- //return StatusCode::SUCCESS;//FIXME return when no hits in DC + silicon
+ const edm4hep::SimTrackerHitCollection* dcSimHitCol
+ =m_DCSimTrackerHitCol.get();
+ const edm4hep::MCRecoTrackerAssociationCollection* assoHits
+ = m_DCHitAssociationCol.get();
+ if(m_useDC){
+ if(nullptr==dcSimHitCol){
+ debug()<<"DC SimTrackerHitCollection not found"<<endmsg;
+ }else{
+ debug()<<"DriftChamberHitsCollection size "
+ <<dcSimHitCol->size()<<endmsg;
+ if(m_tuple)m_nSimTrackerHitDC=dcSimHitCol->size();
+ }
+ digiDCHitsCol=m_DCDigiCol.get();
+ m_nDCTrackHit = digiDCHitsCol->size();
+ if(nullptr==digiDCHitsCol){
+ debug()<<"TrackerHitCollection not found"<<endmsg;
+ }else{
+ debug()<<"digiDCHitsCol size "<<digiDCHitsCol->size()<<endmsg;
+ if((int) digiDCHitsCol->size()>m_maxDCDigiCut){
+ debug()<<"Track cut by m_maxDCDigiCut "<<m_maxDCDigiCut<<endmsg;
+ return StatusCode::SUCCESS;
+ }
+ }
}
- if((int) digiDCHitsCol->size()>m_maxDCDigiCut) return StatusCode::SUCCESS;
- ///Output Track collection
- edm4hep::TrackCollection* dcTrackCol=m_DCTrackCol.createAndPut();
- edm4hep::TrackCollection* sdtTrackCol=m_SDTTrackCol.createAndPut();
- edm4hep::ReconstructedParticleCollection* dcRecParticleCol(nullptr);
- if(m_writeRecParticle){
- dcRecParticleCol=m_DCRecParticleCol.createAndPut();
- }
+ // make noise
+ edm4hep::SimTrackerHitCollection* SimVec = w_SimNoiseHCol.createAndPut();
+ edm4hep::MCRecoTrackerAssociationCollection* AssoVec =
+ w_NoiseAssociationCol.createAndPut();
+ edm4hep::TrackerHitCollection* Vec = w_NoiseHitCol.createAndPut();
+
////TODO
//Output MCRecoTrackerAssociationCollection collection
//const edm4hep::MCRecoTrackerAssociationCollection*
@@ -115,50 +231,184 @@ StatusCode TruthTrackerAlg::execute()
//}
//mcRecoTrackerAssociationCol=m_mcRecoParticleAssociation.get();
- ///Retrieve silicon Track
- const edm4hep::TrackCollection* siTrackCol=nullptr;
- if(m_siSubsetTrackCol.exist()) siTrackCol=m_siSubsetTrackCol.get();
- if(nullptr!=siTrackCol) {
- ///New SDT track
- for(auto siTrack:*siTrackCol){
- auto 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));
+ ///New SDT track
+ edm4hep::MutableTrack sdtTk=sdtTkCol->create();
+
+ int nVXDHit=0;
+ int nSITHit=0;
+ int nSETHit=0;
+ int nFTDHit=0;
+ int nDCHitDCTk=0;
+ int nDCHitSDTTk=0;
+
+ ///Create track with mcParticle
+ edm4hep::TrackState trackStateMc;
+ getTrackStateFromMcParticle(mcParticleCol,trackStateMc);
+ if(m_useTruthTrack.value()||!m_useSi){sdtTk.addToTrackStates(trackStateMc);}
+
+ if(m_useSi){
+ ///Retrieve silicon Track
+ const edm4hep::TrackCollection* siTrackCol=nullptr;
+ //if(m_siSubsetTrackCol.exist()){
+ siTrackCol=m_siSubsetTrackCol.get();
+ if(nullptr==siTrackCol){
+ debug()<<"SDTTrackCollection is empty"<<endmsg;
+ if(!m_useTruthTrack.value()) return StatusCode::SUCCESS;
+ }else{
+ debug()<<"SiSubsetTrackCol size "<<siTrackCol->size()<<endmsg;
+ if(!m_useTruthTrack.value()&&0==siTrackCol->size()){
+ return StatusCode::SUCCESS;
}
- //TODO tracks
- for(auto digiDC:*digiDCHitsCol){
- //if(Sim->MCParti!=current) continue;//TODO
- sdtTrack.addToTrackerHits(digiDC);
+ }
+ //}
+
+ for(auto siTk:*siTrackCol){
+ if(!m_useTruthTrack.value()){
+ debug()<<"siTk: "<<siTk<<endmsg;
+ edm4hep::TrackState siTrackStat=siTk.getTrackStates(0);//FIXME?
+ sdtTk.addToTrackStates(siTrackStat);
+ sdtTk.setType(siTk.getType());
+ sdtTk.setChi2(siTk.getChi2());
+ sdtTk.setNdf(siTk.getNdf());
+ sdtTk.setDEdx(siTk.getDEdx());
+ sdtTk.setDEdxError(siTk.getDEdxError());
+ sdtTk.setRadiusOfInnermostHit(
+ siTk.getRadiusOfInnermostHit());
}
+ if(!m_useSiTruthHit){
+ debug()<<"use Si hit on track"<<endmsg;
+ nVXDHit=addHotsToTk(siTk,sdtTk,lcio::ILDDetID::VXD,"VXD",nVXDHit);
+ nSITHit=addHotsToTk(siTk,sdtTk,lcio::ILDDetID::SIT,"SIT",nSITHit);
+ nSETHit=addHotsToTk(siTk,sdtTk,lcio::ILDDetID::SET,"SET",nSETHit);
+ nFTDHit=addHotsToTk(siTk,sdtTk,lcio::ILDDetID::FTD,"FTD",nFTDHit);
+ }//end of loop over hits on siTk
+ }//end of loop over siTk
+
+ if(m_useSiTruthHit){
+ ///Add silicon SimTrackerHit
+ debug()<<"Add silicon SimTrackerHit"<<endmsg;
+ nVXDHit=addSimHitsToTk(m_VXDCollection,truthTrackerHitCol,sdtTk,"VXD",nVXDHit);
+ nSITHit=addSimHitsToTk(m_SITCollection,truthTrackerHitCol,sdtTk,"SIT",nSITHit);
+ nSETHit=addSimHitsToTk(m_SETCollection,truthTrackerHitCol,sdtTk,"SET",nSETHit);
+ nFTDHit=addSimHitsToTk(m_FTDCollection,truthTrackerHitCol,sdtTk,"FTD",nFTDHit);
+ }else{
+ ///Add reconstructed hit or digi
+ debug()<<"Add VXD TrackerHit"<<endmsg;
+ nVXDHit=addHitsToTk(m_VXDTrackerHits,sdtTk,"VXD digi",nVXDHit);
+ nSITHit=addHitsToTk(m_SITTrackerHits,sdtTk,"SIT digi",nSITHit);
+ if(m_useSiSpacePoint.value()){
+ ///Add silicon SpacePoint
+ debug()<<"Add silicon SpacePoint"<<endmsg;
+ if(m_useSiSpacePoint){
+ nSITHit=addHitsToTk(m_SITSpacePointCol,sdtTk,"SIT sp",nSITHit);
+ }
+ // nSETHit=addHitsToTk(m_SETSpacePointCol,sdtTk,"SET sp",nSETHit);
+ nFTDHit=addHitsToTk(m_FTDSpacePointCol,sdtTk,"FTD sp",nFTDHit);
+ }else{
+ ///Add silicon TrackerHit
+ debug()<<"Add silicon TrackerHit"<<endmsg;
+ nSITHit=addHitsToTk(m_SITTrackerHits,sdtTk,"SIT digi",nSITHit);
+ nSETHit=addHitsToTk(m_SETTrackerHits,sdtTk,"SET digi",nSETHit);
+ nFTDHit=addHitsToTk(m_FTDTrackerHits,sdtTk,"FTD digi",nFTDHit);
+ }//end of use space point
+ }
+ }//end of use silicon
+
+ if(m_useDC){
+ ///Create DC Track
+ edm4hep::MutableTrack dcTrack=dcTrackCol->create();
+
+ //Create TrackState
+ edm4hep::TrackState trackStateFirstDCHit;
+ float charge=trackStateMc.omega/fabs(trackStateMc.omega);
+ if(m_useFirstHitForDC&&getTrackStateFirstHit(m_DCSimTrackerHitCol,
+ charge,trackStateFirstDCHit)){
+ dcTrack.addToTrackStates(trackStateFirstDCHit);
+ dcTrack.addToTrackStates(trackStateMc);
+ }else{
+ dcTrack.addToTrackStates(trackStateMc);
+ dcTrack.addToTrackStates(trackStateFirstDCHit);
+ }
+ //sdtTk.addToTrackStates(trackStateFirstDCHit);
+
+ ///Add other track properties
+ dcTrack.setNdf(dcTrack.trackerHits_size()-5);
+
+ ///Add DC hits to tracks after track state set
+ if(m_useIdealHit){
+ nDCHitDCTk=addIdealHitsToTk(m_DCDigiCol,truthTrackerHitCol,dcTrack,
+ "DC digi",nDCHitDCTk);
+ }else{
+ nDCHitDCTk=addHitsToTk(m_DCDigiCol,dcTrack,"DC digi",nDCHitDCTk);
+ }
+// if(m_useSi) nDCHitSDTTk=addHitsToTk(m_DCDigiCol,sdtTk,"DC digi",nDCHitSDTTk);
+ // if(m_useSi)
+ // {
+ if(m_smearHits){
+ m_nSmearDCTrackHit = smearDCTkhit(m_DCDigiCol,
+ m_SmeartruthTrackerHitCol,m_DCSimTrackerHitCol,
+ w_SimSmearHCol,m_DCHitAssociationCol,
+ w_SmearAssociationCol,m_resX,m_resY,m_resZ);
}
+
+ if(!m_useNoiseHits)
+ {
+ //nDCHitSDTTk=addHitsToTk(m_DCDigiCol,sdtTk,"DC digi",nDCHitSDTTk);
+ nDCHitSDTTk=addHitsToTk(m_SmeartruthTrackerHitCol,sdtTk,"DC digi",nDCHitSDTTk);
+ } else {
+ m_nNoiseOnSdtTk = makeNoiseHit(SimVec,Vec,AssoVec,
+ m_SmeartruthTrackerHitCol.get(),w_SmearAssociationCol.get());
+ if(debugNoiseHitsCol(Vec))std::cout << " Make noise hits successfully!" << std::endl;
+ nDCHitSDTTk=addHitsToTk(w_NoiseHitCol,sdtTk,
+ "DC digi",nDCHitSDTTk);
+ }
+ //}
+
+ //track.setType();//TODO
+ //track.setChi2(gauss(digiDCHitsCol->size-5(),1));//FIXME
+ //track.setDEdx();//TODO
+
+ debug()<<"dcTrack nHit "<<dcTrack.trackerHits_size()<<dcTrack<<endmsg;
+ }
+
+ ///Set other track parameters
+ //sdtTk.setNdf(sdtTk.trackerHits_size()-5);
+ //double radiusOfInnermostHit=1e9;
+ //edm4hep::Vector3d digiPos=digiDC.getPosition();
+ //double r=sqrt(digiPos.x*digiPos.x+digiPos.y*digiPos.y);
+ //if(r<radiusOfInnermostHit) radiusOfInnermostHit=r;
+
+ debug()<<"sdtTk nHit "<<sdtTk.trackerHits_size()<<sdtTk<<endmsg;
+ debug()<<"nVXDHit "<<nVXDHit<<" nSITHit "<<nSITHit<<" nSETHit "<<nSETHit
+ <<" nFTDHit "<<nFTDHit<<" nDCHitSDTTk "<<nDCHitSDTTk<<endmsg;
+
+ if(m_tuple){
+ m_nHitOnSdtTkVXD=nVXDHit;
+ m_nHitOnSdtTkSIT=nSITHit;
+ m_nHitOnSdtTkSET=nSETHit;
+ m_nHitOnSdtTkFTD=nFTDHit;
+ m_nHitOnSdtTkDC=nDCHitSDTTk;
+ //m_nHitOnSdtTk=sdtTk.trackerHits_size();
+ debugEvent();
+ StatusCode sc=m_tuple->write();
}
- ///Convert MCParticle to Track and ReconstructedParticle
+ return StatusCode::SUCCESS;
+}
+
+StatusCode TruthTrackerAlg::finalize()
+{
+ return GaudiAlgorithm::finalize();
+}
+
+void TruthTrackerAlg::getTrackStateFromMcParticle(
+ const edm4hep::MCParticleCollection* mcParticleCol,
+ edm4hep::TrackState& trackState)
+{
+ ///Convert MCParticle to DC Track and ReconstructedParticle
debug()<<"MCParticleCol size="<<mcParticleCol->size()<<endmsg;
for(auto mcParticle : *mcParticleCol){
/// skip mcParticleVertex do not have enough associated hits TODO
-
///Vertex
const edm4hep::Vector3d mcParticleVertex=mcParticle.getVertex();//mm
edm4hep::Vector3f mcParticleVertexSmeared;//mm
@@ -169,7 +419,7 @@ StatusCode TruthTrackerAlg::execute()
mcParticleVertexSmeared.z=
CLHEP::RandGauss::shoot(mcParticleVertex.z,m_resVertexZ);
///Momentum
- edm4hep::Vector3f mcParticleMom=mcParticle.getMomentum();//GeV
+ const edm4hep::Vector3f mcParticleMom=mcParticle.getMomentum();//GeV
double mcParticlePt=sqrt(mcParticleMom.x*mcParticleMom.x+
mcParticleMom.y*mcParticleMom.y);
//double mcParticlePtSmeared=
@@ -180,8 +430,7 @@ StatusCode TruthTrackerAlg::execute()
edm4hep::Vector3f mcParticleMomSmeared;
mcParticleMomSmeared.x=mcParticlePt*cos(mcParticleMomPhiSmeared);
mcParticleMomSmeared.y=mcParticlePt*sin(mcParticleMomPhiSmeared);
- mcParticleMomSmeared.z=
- CLHEP::RandGauss::shoot(mcParticleMom.z,m_resPz);
+ mcParticleMomSmeared.z=CLHEP::RandGauss::shoot(mcParticleMom.z,m_resPz);
///Converted to Helix
double B[3]={1e9,1e9,1e9};
@@ -192,15 +441,17 @@ StatusCode TruthTrackerAlg::execute()
//float mom[3]={mcParticleMomSmeared.x,mcParticleMomSmeared.y,
// mcParticleMomSmeared.z};
////FIXME DEBUG
- float pos[3]={(float)mcParticleVertex.x,
- (float)mcParticleVertex.y,(float)mcParticleVertex.z};
- float mom[3]={(float)mcParticleMom.x,(float)mcParticleMom.y,
- (float)mcParticleMom.z};
+ double pos[3]={mcParticleVertex.x,mcParticleVertex.y,mcParticleVertex.z};//mm
+ double mom[3]={mcParticleMom.x,mcParticleMom.y,mcParticleMom.z};//mm
helix.Initialize_VP(pos,mom,mcParticle.getCharge(),B[2]/dd4hep::tesla);
+ if(m_tuple) {
+ for(int ii=0;ii<3;ii++) {
+ m_mcMom[ii]=mom[ii];//GeV
+ m_mcPos[ii]=pos[ii];//mm
+ }
+ }
///new Track
- auto dcTrack=dcTrackCol->create();
- edm4hep::TrackState trackState;
trackState.D0=helix.getD0();
trackState.phi=helix.getPhi0();
trackState.omega=helix.getOmega();
@@ -208,49 +459,19 @@ StatusCode TruthTrackerAlg::execute()
trackState.tanLambda=helix.getTanLambda();
trackState.referencePoint=helix.getReferencePoint();
std::array<float,15> covMatrix;
- for(int i=0;i<15;i++){covMatrix[i]=999.;}//FIXME
+ for(int i=0;i<15;i++){covMatrix[i]=1.;}//FIXME
trackState.covMatrix=covMatrix;
- dcTrack.addToTrackStates(trackState);
- //dcTrack.setType();//TODO
- //dcTrack.setChi2(gauss(digiDCHitsCol->size-5(),1));//FIXME
- dcTrack.setNdf(digiDCHitsCol->size()-5);
- //dcTrack.setDEdx();//TODO
- //set hits
- double radiusOfInnermostHit=1e9;
- debug()<<"digiDCHitsCol size"<<digiDCHitsCol->size()<<endmsg;
- for(auto digiDC : *digiDCHitsCol){
- //if(Sim->MCParti!=current) continue;//TODO
- edm4hep::Vector3d digiPos=digiDC.getPosition();
- double r=sqrt(digiPos.x*digiPos.x+digiPos.y*digiPos.y);
- if(r<radiusOfInnermostHit) radiusOfInnermostHit=r;
- dcTrack.addToTrackerHits(digiDC);
- }
- dcTrack.setRadiusOfInnermostHit(radiusOfInnermostHit);//TODO
-
- edm4hep::MutableReconstructedParticle dcRecParticle;
- if(m_writeRecParticle){
- dcRecParticle=dcRecParticleCol->create();
- ///new ReconstructedParticle
- //dcRecParticle.setType();//TODO
- double mass=mcParticle.getMass();
- double p=sqrt(mcParticleMomSmeared.x*mcParticleMomSmeared.x+
- mcParticleMomSmeared.y*mcParticleMomSmeared.y+
- mcParticleMomSmeared.z*mcParticleMomSmeared.z);
- dcRecParticle.setEnergy(sqrt(mass*mass+p*p));
- dcRecParticle.setMomentum(mcParticleMomSmeared);
- dcRecParticle.setReferencePoint(mcParticleVertexSmeared);
- dcRecParticle.setCharge(mcParticle.getCharge());
- dcRecParticle.setMass(mass);
- //dcRecParticle.setGoodnessOfPID();//TODO
- //std::array<float>,10> covMat=?;//TODO
- //dcRecParticle.setCovMatrix(covMat);//TODO
- //dcRecParticle.setStartVertex();//TODO
- edm4hep::ParticleID particleID(0,mcParticle.getPDG(),0,1);//FIXME
- dcRecParticle.setParticleIDUsed(particleID);
- dcRecParticle.addToTracks(dcTrack);
- }//end of write RecParticle
+ getCircleFromPosMom(pos,mom,B[2]/dd4hep::tesla,mcParticle.getCharge(),m_helixRadius,m_helixXC,m_helixYC);
+
+ debug()<<"dd4hep::mm "<<dd4hep::mm<<" dd4hep::cm "<<dd4hep::cm<<endmsg;
debug()<<"mcParticle "<<mcParticle
+ <<" helix radius "<<helix.getRadius()<<" "<<helix.getXC()<<" "
+ <<helix.getYC()<<" mm "
+ <<" myhelix radius "<<m_helixRadius<<" "<<m_helixXC<<" "
+ <<m_helixYC<<" mm "
+ <<" momMC "<<mom[0]<<" "<<mom[1]<<" "<<mom[2]<<"GeV"
+ <<" posMC "<<pos[0]<<" "<<pos[1]<<" "<<pos[2]<<"mm"
<<" momPhi "<<mcParticleMomPhi
<<" mcParticleVertex("<<mcParticleVertex<<")mm "
<<" mcParticleVertexSmeared("<<mcParticleVertexSmeared<<")mm "
@@ -258,17 +479,443 @@ StatusCode TruthTrackerAlg::execute()
<<" mcParticleMomSmeared("<<mcParticleMomSmeared<<")GeV "
<<" Bxyz "<<B[0]/dd4hep::tesla<<" "<<B[1]/dd4hep::tesla
<<" "<<B[2]/dd4hep::tesla<<" tesla"<<endmsg;
- debug()<<"trackState:location,D0,phi,omega,Z0,tanLambda"
- <<",referencePoint,cov"<<std::endl<<trackState<<std::endl;
- debug()<<"dcTrack"<<dcTrack<<endmsg;
- if(m_writeRecParticle) debug()<<"dcRecParticle"<<dcRecParticle<<endmsg;
}//end loop over MCParticleCol
+}//end of getTrackStateFromMcParticle
- debug()<<"Output DCTrack size="<<dcTrackCol->size()<<endmsg;
- return StatusCode::SUCCESS;
+bool TruthTrackerAlg::getTrackStateFirstHit(
+ DataHandle<edm4hep::SimTrackerHitCollection>& dcSimTrackerHitCol,
+ float charge,edm4hep::TrackState& trackState)
+{
+
+ const edm4hep::SimTrackerHitCollection* col=nullptr;
+ col=dcSimTrackerHitCol.get();
+ debug()<<"TruthTrackerAlg::getTrackStateFirstHit"<<endmsg;
+ debug()<<"simTrackerHitCol size "<<col->size()<<endmsg;
+ float minHitTime=1e9;
+ if(nullptr!=col||0==col->size()){
+ edm4hep::SimTrackerHit firstHit;
+ for(auto dcSimTrackerHit:*col){
+ const edm4hep::Vector3f mom=dcSimTrackerHit.getMomentum();
+ if(abs(sqrt(mom[0]*mom[0]+mom[1]*mom[1]))<m_momentumCut)continue;//yzhang TEMP skip hits with momentum <0.5GeV/c
+ if(dcSimTrackerHit.getTime()<minHitTime) {
+ minHitTime=dcSimTrackerHit.getTime();
+ firstHit=dcSimTrackerHit;
+ }
+ //debug()<<"simTrackerHit time "<<dcSimTrackerHit.getTime()
+ // <<" pos "<<dcSimTrackerHit.getPosition()
+ // <<" mom "<<dcSimTrackerHit.getMomentum()<<endmsg;
+ }
+ const edm4hep::Vector3d pos=firstHit.getPosition();
+ const edm4hep::Vector3f mom=firstHit.getMomentum();
+ debug()<<"first Hit pos "<<pos<<" mom "<<mom<<" time "<<minHitTime<<endmsg;
+ float pos_t[3]={(float)pos[0],(float)pos[1],(float)pos[2]};
+ float mom_t[3]={(float)mom[0],(float)mom[1],(float)mom[2]};
+ ///Converted to Helix
+ double B[3]={1e9,1e9,1e9};
+ m_dd4hepField.magneticField({0.,0.,0.},B);
+ HelixClass helix;
+ helix.Initialize_VP(pos_t,mom_t,charge,B[2]/dd4hep::tesla);
+ m_helixRadiusFirst=helix.getRadius();
+ m_helixXCFirst=helix.getXC();
+ m_helixYCFirst=helix.getYC();
+
+ ///new Track
+ trackState.D0=helix.getD0();
+ trackState.phi=helix.getPhi0();
+ trackState.omega=helix.getOmega();
+ trackState.Z0=helix.getZ0();
+ trackState.tanLambda=helix.getTanLambda();
+ trackState.referencePoint=helix.getReferencePoint();
+ std::array<float,15> covMatrix;
+ for(int i=0;i<15;i++){covMatrix[i]=100.;}//FIXME
+ trackState.covMatrix=covMatrix;
+ debug()<<"first hit trackState "<<trackState<<endmsg;
+ return true;
+ }
+ return false;
+}//end of getTrackStateFirstHit
+
+void TruthTrackerAlg::debugEvent()
+{
+ if(m_useSi){
+ ///Retrieve silicon Track
+ const edm4hep::TrackCollection* siTrackCol=nullptr;
+ siTrackCol=m_siSubsetTrackCol.get();
+ if(nullptr!=siTrackCol){
+ for(auto siTk:*siTrackCol){
+ debug()<<"siTk: "<<siTk<<endmsg;
+ edm4hep::TrackState trackStat=siTk.getTrackStates(0);//FIXME?
+ double B[3]={1e9,1e9,1e9};
+ m_dd4hepField.magneticField({0.,0.,0.},B);
+ HelixClass helix;
+ helix.Initialize_Canonical(trackStat.phi, trackStat.D0, trackStat.Z0,
+ trackStat.omega, trackStat.tanLambda, B[2]/dd4hep::tesla);
+
+ if(m_tuple){
+ m_siMom[0]=helix.getMomentum()[0];
+ m_siMom[1]=helix.getMomentum()[1];
+ m_siMom[2]=helix.getMomentum()[2];
+ m_siPos[0]=helix.getReferencePoint()[0];
+ m_siPos[1]=helix.getReferencePoint()[1];
+ m_siPos[2]=helix.getReferencePoint()[2];
+ m_nHitOnSiTkVXD=nHotsOnTrack(siTk,lcio::ILDDetID::VXD);
+ m_nHitOnSiTkSIT=nHotsOnTrack(siTk,lcio::ILDDetID::SIT);
+ m_nHitOnSiTkSET=nHotsOnTrack(siTk,lcio::ILDDetID::SET);
+ m_nHitOnSiTkFTD=nHotsOnTrack(siTk,lcio::ILDDetID::FTD);
+ }
+ }//end of loop over siTk
+ }
+ if(m_tuple){
+ //SimTrackerHits
+ m_nSimTrackerHitVXD=simTrackerHitColSize(m_VXDCollection);
+ m_nSimTrackerHitSIT=simTrackerHitColSize(m_SITCollection);
+ m_nSimTrackerHitSET=simTrackerHitColSize(m_SETCollection);
+ m_nSimTrackerHitFTD=simTrackerHitColSize(m_FTDCollection);
+
+ //TrackerHits
+ m_nTrackerHitVXD=trackerHitColSize(m_VXDTrackerHits);
+ m_nTrackerHitSIT=trackerHitColSize(m_SITTrackerHits);
+ m_nTrackerHitSET=trackerHitColSize(m_SETTrackerHits);
+ m_nTrackerHitFTD=trackerHitColSize(m_FTDTrackerHits);
+ m_nTrackerHitDC=trackerHitColSize(m_DCDigiCol);
+
+ //SpacePoints
+ if(m_useSiSpacePoint){
+ m_nSpacePointSIT=trackerHitColSize(m_SITSpacePointCol);
+ }
+ //m_nSpacePointSET=trackerHitColSize(m_SETSpacePointCol);
+ m_nSpacePointFTD=trackerHitColSize(m_FTDSpacePointCol);
+ }
+ }
}
-StatusCode TruthTrackerAlg::finalize()
+int TruthTrackerAlg::addIdealHitsToTk(DataHandle<edm4hep::TrackerHitCollection>&
+ colHandle, edm4hep::TrackerHitCollection*& truthTrackerHitCol,
+ edm4hep::MutableTrack& track, const char* msg,int nHitAdded)
{
- return GaudiAlgorithm::finalize();
+ if(nHitAdded>0) return nHitAdded;
+ int nHit=0;
+ const edm4hep::TrackerHitCollection* col=colHandle.get();
+ debug()<<"add "<<msg<<" "<<col->size()<<" trackerHit"<<endmsg;
+ debug()<<track<<endmsg;
+ for(auto hit:*col){
+ //get end point of this wire
+ TVector3 endPointStart(0,0,0);
+ TVector3 endPointEnd(0,0,0);
+ m_gridDriftChamber->cellposition(hit.getCellID(),endPointStart,
+ endPointEnd);//cm
+
+ //calc. doca of helix to wire
+ TVector3 wire(endPointStart.X()/dd4hep::mm,endPointStart.Y()/dd4hep::mm,0);//to mm
+ TVector3 center(m_helixXC,m_helixYC,0);//mm
+ double docaIdeal=(center-wire).Mag()-m_helixRadius;//mm
+ TVector3 centerFirst(m_helixXCFirst,m_helixYCFirst,0);//mm
+ double docaIdealFirst=(centerFirst-wire).Mag()-m_helixRadiusFirst;//mm
+
+ //add modified hit
+ auto tmpHit = truthTrackerHitCol->create();
+ tmpHit=hit.clone();
+ tmpHit.setTime(fabs(docaIdeal)*1e3/40.);//40#um/ns, drift time in ns
+ track.addToTrackerHits(tmpHit);
+
+ long long int detID=hit.getCellID();
+ debug()<<" addIdealHitsToTk "<<m_helixRadius<<" center "<<m_helixXC
+ <<" "<<m_helixYC<<" mm wire("<<m_decoder->get(detID,"layer")<<","
+ <<m_decoder->get(detID,"cellID")<<") "<<wire.X()<<" "<<wire.Y()
+ <<"mm docaIdeal "<<docaIdeal<<" docaIdealFirst "<<docaIdealFirst<<"mm "
+ <<"hit.Time orignal "<<hit.getTime()<<" new Time "
+ <<fabs(docaIdeal)*1e3/40.<<endmsg;
+ ++nHit;
+ }
+ return nHit;
+}
+
+bool TruthTrackerAlg::debugNoiseHitsCol(edm4hep::TrackerHitCollection* Vec)
+{
+ m_nNoiseDCTrackHit = Vec->size();
+std::cout << " m_nNoiseDCTrackHit = " << m_nNoiseDCTrackHit << std::endl;
+ int ihit=0;
+ for(auto Vechit: *Vec)
+ {
+ m_NoiseDriftDistance[ihit] = (Vechit.getTime())*m_driftVelocity*1e-3;
+ ihit++;
+ }
+
+ return true;
+}
+
+int TruthTrackerAlg::makeNoiseHit(edm4hep::SimTrackerHitCollection* SimVec,
+ edm4hep::TrackerHitCollection* Vec,
+ edm4hep::MCRecoTrackerAssociationCollection* AssoVec,
+ const edm4hep::TrackerHitCollection* digiDCHitsCol,
+ const edm4hep::MCRecoTrackerAssociationCollection* assoHits)
+{
+ int nNoise = 0;
+ //loop all hits
+ for(unsigned int i = 0; i<(digiDCHitsCol->size()); i++)
+ {
+ edm4hep::TrackerHit mcHit = digiDCHitsCol->at(i);
+ unsigned long long wcellid = mcHit.getCellID();
+
+ float pocaTime = mcHit.getTime();
+
+ // store trackHit
+ auto trkHit = Vec->create();
+
+ trkHit.setCellID(wcellid);
+ trkHit.setTime(pocaTime);
+ trkHit.setEDep(mcHit.getEDep());
+ trkHit.setEdx(mcHit.getEdx());
+ trkHit.setPosition(mcHit.getPosition());
+ trkHit.setCovMatrix(mcHit.getCovMatrix());
+ for(int iAsso=0;iAsso<(int) assoHits->size();iAsso++)
+ {
+
+ if(assoHits->at(iAsso).getRec().getCellID()==wcellid )
+ {
+ auto SimtrkHit = SimVec->create();
+
+ SimtrkHit.setCellID(assoHits->at(iAsso).getSim().getCellID());
+ SimtrkHit.setEDep(assoHits->at(iAsso).getSim().getEDep());
+ SimtrkHit.setTime(assoHits->at(iAsso).getSim().getTime());
+ SimtrkHit.setPathLength(assoHits->at(iAsso).getSim().getPathLength());
+ SimtrkHit.setQuality(assoHits->at(iAsso).getSim().getQuality());
+ SimtrkHit.setPosition(assoHits->at(iAsso).getSim().getPosition());
+ SimtrkHit.setMomentum(assoHits->at(iAsso).getSim().getMomentum());
+ SimtrkHit.setMCParticle(assoHits->at(iAsso).getSim().getMCParticle());
+
+ auto asso = AssoVec->create();
+ asso.setRec(trkHit);
+ asso.setSim(SimtrkHit);
+ asso.setWeight(SimtrkHit.getEDep()/trkHit.getEDep());
+ }
+ }
+
+ double prob1 = fRandom.Uniform(1.);
+ if(prob1 > m_fHitPurity) continue;
+
+ float noiseTime = m_pocaTime*fRandom.Uniform(1.);
+ if(noiseTime>pocaTime) continue;
+ nNoise++;
+ trkHit.setTime(noiseTime);
+
+ for(int iAsso=0;iAsso<(int) assoHits->size();iAsso++)
+ {
+
+ if(assoHits->at(iAsso).getRec().getCellID()==wcellid )
+ {
+
+ AssoVec->at(iAsso).setRec(trkHit);
+
+ }
+
+ }
+
+ }
+
+std::cout << " Truth Noise number = " << nNoise << std::endl;
+ return nNoise;
+}
+
+int TruthTrackerAlg::smearDCTkhit(DataHandle<edm4hep::TrackerHitCollection>&
+ colHandle,DataHandle<edm4hep::TrackerHitCollection>& smearCol,
+ DataHandle<edm4hep::SimTrackerHitCollection>& SimDCHitCol,
+ DataHandle<edm4hep::SimTrackerHitCollection>& SimSmearDCHitCol,
+ DataHandle<edm4hep::MCRecoTrackerAssociationCollection>& AssoDCHitCol,
+ DataHandle<edm4hep::MCRecoTrackerAssociationCollection>& AssoSmearDCHitCol,
+ double resX, double resY, double resZ)
+{
+ int nHit=0;
+ const edm4hep::TrackerHitCollection* col=colHandle.get();
+ auto smearTrackerHitCol = smearCol.createAndPut();
+
+ auto simDCCol = SimDCHitCol.get();
+ auto SimVec = SimSmearDCHitCol.createAndPut();
+
+ auto assoHits = AssoDCHitCol.get();
+ auto AssoVec = AssoSmearDCHitCol.createAndPut();
+
+ //sort,FIXME
+ for(auto hit:*col){
+
+ auto smearHit=smearTrackerHitCol->create();
+
+ float truthD = (hit.getTime())*m_driftVelocity*1e-3;
+ m_DriftDistance[nHit] = truthD;
+
+ truthD+=gRandom->Gaus(0,fabs(m_resX));
+ m_SmearDriftDistance[nHit] = truthD;
+ float smearT = (truthD*1e3)/m_driftVelocity;
+
+ smearHit.setTime(smearT);
+ smearHit.setCellID(hit.getCellID());
+ smearHit.setType(hit.getCellID());
+ smearHit.setQuality(hit.getQuality());
+ smearHit.setEDep(hit.getEDep());
+ smearHit.setEDepError(hit.getEDepError());
+ smearHit.setEdx(hit.getEdx());
+ smearHit.setPosition(hit.getPosition());
+ smearHit.setCovMatrix(hit.getCovMatrix());
+ smearHit.addToRawHits(hit.getObjectID());
+
+ ++nHit;
+ for(int iAsso=0;iAsso<(int) assoHits->size();iAsso++)
+ {
+
+ if(assoHits->at(iAsso).getRec().getCellID()== smearHit.getCellID())
+ {
+ auto SimtrkHit = SimVec->create();
+
+ SimtrkHit.setCellID(assoHits->at(iAsso).getSim().getCellID());
+ SimtrkHit.setEDep(assoHits->at(iAsso).getSim().getEDep());
+ SimtrkHit.setTime(assoHits->at(iAsso).getSim().getTime());
+ SimtrkHit.setPathLength(assoHits->at(iAsso).getSim().getPathLength());
+ SimtrkHit.setQuality(assoHits->at(iAsso).getSim().getQuality());
+ SimtrkHit.setPosition(assoHits->at(iAsso).getSim().getPosition());
+ SimtrkHit.setMomentum(assoHits->at(iAsso).getSim().getMomentum());
+ SimtrkHit.setMCParticle(assoHits->at(iAsso).getSim().getMCParticle());
+
+ auto asso = AssoVec->create();
+ asso.setRec(smearHit);
+ asso.setSim(SimtrkHit);
+ asso.setWeight(SimtrkHit.getEDep()/smearHit.getEDep());
+ }
+ }
+ }
+ std::cout << " SimSmear size = " << SimVec->size() << " Asso size = " << AssoVec->size() << std::endl;
+ std::cout << " Sim size = " << simDCCol->size() << " Asso size = " << assoHits->size() << std::endl;
+
+ return nHit;
+}
+
+int TruthTrackerAlg::addHitsToTk(DataHandle<edm4hep::TrackerHitCollection>&
+ colHandle, edm4hep::MutableTrack& track, const char* msg,int nHitAdded)
+{
+ if(nHitAdded>0) return nHitAdded;
+ int nHit=0;
+ const edm4hep::TrackerHitCollection* col=colHandle.get();
+ debug()<<"add "<<msg<<" "<<col->size()<<" trackerHit"<<endmsg;
+ //sort,FIXME
+ for(auto hit:*col){
+ track.addToTrackerHits(hit);
+ ++nHit;
+ }
+ return nHit;
+}
+
+int TruthTrackerAlg::addSimHitsToTk(
+ DataHandle<edm4hep::SimTrackerHitCollection>& colHandle,
+ edm4hep::TrackerHitCollection*& truthTrackerHitCol,
+ edm4hep::MutableTrack& track, const char* msg,int nHitAdded)
+{
+ if(nHitAdded>0) return nHitAdded;
+ int nHit=0;
+ const edm4hep::SimTrackerHitCollection* col=colHandle.get();
+ for(auto simTrackerHit:*col){
+ auto trackerHit=truthTrackerHitCol->create();
+ if(m_skipSecondaryHit&&simTrackerHit.isProducedBySecondary()) {
+ debug()<<"skip secondary simTrackerHit "<<msg<<endmsg;
+ continue;
+ }
+ auto& pos = simTrackerHit.getPosition();
+ debug()<<" addSimHitsToTk "<<msg<<" "<<sqrt(pos.x*pos.x+pos.y*pos.y)<<endmsg;
+ UTIL::BitField64 encoder(lcio::ILDCellID0::encoder_string) ;
+ int detID=encoder[lcio::ILDCellID0::subdet] ;
+ double resolution[3];
+ if(lcio::ILDDetID::VXD==detID){
+ for(int i=0;i<3;i++)resolution[i]=m_resVXD[i];
+ }else if(lcio::ILDDetID::SIT==detID){
+ for(int i=0;i<3;i++)resolution[i]=m_resSIT[i];
+ }else if(lcio::ILDDetID::SET==detID){
+ for(int i=0;i<3;i++)resolution[i]=m_resSET[i];
+ }else if(lcio::ILDDetID::FTD==detID){
+ for(int i=0;i<3;i++)resolution[i]=m_resFTDPixel[i];//FIXME
+ }else{
+ for(int i=0;i<3;i++)resolution[i]=0.003;
+ }
+ edm4hep::Vector3d posSmeared;//mm
+ posSmeared.x=CLHEP::RandGauss::shoot(pos.x,resolution[0]);
+ posSmeared.y=CLHEP::RandGauss::shoot(pos.y,resolution[1]);
+ posSmeared.z=CLHEP::RandGauss::shoot(pos.z,resolution[2]);
+ trackerHit.setPosition(posSmeared) ;
+ encoder.setValue(simTrackerHit.getCellID()) ;
+ trackerHit.setCellID(encoder.lowWord());//?FIXME
+ std::array<float, 6> cov;
+ cov[0]=resolution[0]*resolution[0];
+ cov[1]=0.;
+ cov[2]=resolution[1]*resolution[1];
+ cov[3]=0.;
+ cov[4]=0.;
+ cov[5]=resolution[2]*resolution[2];
+ trackerHit.setCovMatrix(cov);
+ debug()<<"add simTrackerHit "<<msg<<" trackerHit "<<trackerHit<<endmsg;
+ ///Add hit to track
+ track.addToTrackerHits(trackerHit);
+ trackerHit.setEDep(simTrackerHit.getEDep());
+ trackerHit.addToRawHits(simTrackerHit.getObjectID());
+ trackerHit.setType(-8);//FIXME?
+ ++nHit;
+ }
+ debug()<<"add simTrackerHit "<<msg<<" "<<nHit<<endmsg;
+ return nHit;
+}
+
+int TruthTrackerAlg::addHotsToTk(edm4hep::Track& sourceTrack,
+ edm4hep::MutableTrack& targetTrack, int hitType,const char* msg,int nHitAdded)
+{
+ if(nHitAdded>0) return nHitAdded;
+ int nHit=0;
+ for(unsigned int iHit=0;iHit<sourceTrack.trackerHits_size();iHit++){
+ edm4hep::TrackerHit hit=sourceTrack.getTrackerHits(iHit);
+ UTIL::BitField64 encoder(lcio::ILDCellID0::encoder_string);
+ encoder.setValue(hit.getCellID());
+ if(encoder[lcio::ILDCellID0::subdet]==hitType){
+ targetTrack.addToTrackerHits(hit);
+ debug()<<endmsg<<" add siHit "<<msg<<" "<<iHit<<" "<<hit
+ <<" pos "<<hit.getPosition().x<<" "<<hit.getPosition().y<<" "
+ <<hit.getPosition().z<<" " <<endmsg;
+ ++nHit;
+ }
+ }
+ debug()<<endmsg<<" "<<nHit<<" "<<msg<<" hits add on track"<<endmsg;
+ return nHit;
+}
+
+int TruthTrackerAlg::nHotsOnTrack(edm4hep::Track& track, int hitType)
+{
+ int nHit=0;
+ for(unsigned int iHit=0;iHit<track.trackerHits_size();iHit++){
+ edm4hep::TrackerHit hit=track.getTrackerHits(iHit);
+ UTIL::BitField64 encoder(lcio::ILDCellID0::encoder_string);
+ encoder.setValue(hit.getCellID());
+ if(encoder[lcio::ILDCellID0::subdet]==hitType){
+ ++nHit;
+ }
+ }
+ return nHit;
+}
+
+int TruthTrackerAlg::trackerHitColSize(DataHandle<edm4hep::TrackerHitCollection>& col)
+{
+ const edm4hep::TrackerHitCollection* c=col.get();
+ if(nullptr!=c) return c->size();
+ return 0;
+}
+
+int TruthTrackerAlg::simTrackerHitColSize(DataHandle<edm4hep::SimTrackerHitCollection>& col)
+{
+ const edm4hep::SimTrackerHitCollection* c=col.get();
+ if(nullptr!=c) return c->size();
+ return 0;
+}
+//unit length is mm
+void TruthTrackerAlg::getCircleFromPosMom(double pos[3],double mom[3],
+ double Bz,double q,double& helixRadius,double& helixXC,double& helixYC)
+{
+ double FCT = 2.99792458E-4;//mm
+ double pxy = sqrt(mom[0]*mom[0]+mom[1]*mom[1]);
+ helixRadius = pxy / (FCT*Bz);
+ double phiMomRefPoint = atan2(mom[1],mom[0]);
+ helixXC= pos[0] + helixRadius*cos(phiMomRefPoint-M_PI*0.5*q);
+ helixYC= pos[1] + helixRadius*sin(phiMomRefPoint-M_PI*0.5*q);
}
diff --git a/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.h b/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.h
index 30efd0a8a71f7d4ec4211fcb6e074e781ecc9fa0..563d98543e54e94128d4846940d170a080b3a21a 100644
--- a/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.h
+++ b/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.h
@@ -4,6 +4,9 @@
#include "GaudiAlg/GaudiAlgorithm.h"
#include "k4FWCore/DataHandle.h"
#include "DD4hep/Fields.h"
+#include "GaudiKernel/NTuple.h"
+
+#include "TRandom3.h"
class IGeomSvc;
namespace dd4hep {
@@ -15,10 +18,14 @@ namespace dd4hep {
}
namespace edm4hep {
class MCParticleCollection;
+ class SimTrackerHitCollection;
class TrackerHitCollection;
class TrackCollection;
- class MCRecoTrackerAssociationCollection;
+ class Track;
+ class MutableTrack;
+ class TrackState;
class ReconstructedParticleCollection;
+ class MCRecoTrackerAssociationCollection;
class MCRecoParticleAssociationCollection;
}
@@ -32,47 +39,199 @@ class TruthTrackerAlg: public GaudiAlgorithm
virtual StatusCode finalize() override;
private:
+ void getTrackStateFromMcParticle(const edm4hep::MCParticleCollection*
+ mcParticleCol, edm4hep::TrackState& stat);
+ int addSimHitsToTk(DataHandle<edm4hep::SimTrackerHitCollection>&
+ colHandle, edm4hep::TrackerHitCollection*& truthTrackerHitCol,
+ edm4hep::MutableTrack& track, const char* msg,int nHitAdded);
+ int smearDCTkhit(DataHandle<edm4hep::TrackerHitCollection>&
+ colHandle,DataHandle<edm4hep::TrackerHitCollection>& smearCol,
+ DataHandle<edm4hep::SimTrackerHitCollection>& SimDCHitCol,
+ DataHandle<edm4hep::SimTrackerHitCollection>& SimSmearDCHitCol,
+ DataHandle<edm4hep::MCRecoTrackerAssociationCollection>& AssoDCHitCol,
+ DataHandle<edm4hep::MCRecoTrackerAssociationCollection>& AssoSmearDCHitCol,
+ double resX, double resY, double resZ);
+ int addHitsToTk(DataHandle<edm4hep::TrackerHitCollection>&
+ colHandle, edm4hep::MutableTrack& track, const char* msg,int nHitAdded);
+ int addIdealHitsToTk(DataHandle<edm4hep::TrackerHitCollection>&
+ colHandle, edm4hep::TrackerHitCollection*& truthTrackerHitCol,
+ edm4hep::MutableTrack& track, const char* msg,int nHitAdded);
+
+ int addHotsToTk(edm4hep::Track& sourceTrack,edm4hep::MutableTrack&
+ targetTrack, int hitType,const char* msg,int nHitAdded);
+ int nHotsOnTrack(edm4hep::Track& track, int hitType);
+ int trackerHitColSize(DataHandle<edm4hep::TrackerHitCollection>& hitCol);
+ int simTrackerHitColSize(DataHandle<edm4hep::SimTrackerHitCollection>& hitCol);
+ bool getTrackStateFirstHit(
+ DataHandle<edm4hep::SimTrackerHitCollection>& dcSimTrackerHitCol,
+ float charge,edm4hep::TrackState& trackState);
SmartIF<IGeomSvc> m_geomSvc;
dd4hep::Detector* m_dd4hep;
dd4hep::OverlayedField m_dd4hepField;
dd4hep::DDSegmentation::GridDriftChamber* m_gridDriftChamber;
dd4hep::DDSegmentation::BitFieldCoder* m_decoder;
+ void debugEvent();
+ //unit length is mm
+ void getCircleFromPosMom(double pos[3],double mom[3],
+ double Bz,double q,double& helixRadius,double& helixXC,double& helixYC);
+ int makeNoiseHit(edm4hep::SimTrackerHitCollection* SimVec,
+ edm4hep::TrackerHitCollection* Vec,
+ edm4hep::MCRecoTrackerAssociationCollection* AssoVec,
+ const edm4hep::TrackerHitCollection* digiDCHitsCol,
+ const edm4hep::MCRecoTrackerAssociationCollection* assoHits);
+ bool debugNoiseHitsCol(edm4hep::TrackerHitCollection* Vec);
- //reader
+ //reader
+ // DataHandle<edm4hep::TrackerHitCollection> m_NoiseHitCol{
+ // "NoiseDCHitsCollection", Gaudi::DataHandle::Reader, this};
DataHandle<edm4hep::MCParticleCollection> m_mcParticleCol{
"MCParticle", Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::SimTrackerHitCollection> m_DCSimTrackerHitCol{
+ "DriftChamberHitsCollection", Gaudi::DataHandle::Reader, this};
DataHandle<edm4hep::TrackerHitCollection> m_DCDigiCol{
"DigiDCHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<edm4hep::MCRecoTrackerAssociationCollection>
m_DCHitAssociationCol{ "DCHitAssociationCollection",
Gaudi::DataHandle::Reader, this};
DataHandle<edm4hep::TrackCollection>
- m_siSubsetTrackCol{ "SiSubsetTrackCollection",
+ m_siSubsetTrackCol{ "SubsetTracks",
Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::TrackerHitCollection> m_SITSpacePointCol{
+ "SITSpacePoints" , Gaudi::DataHandle::Reader, this};
+ // DataHandle<edm4hep::TrackerHitCollection> m_SETSpacePointCol{
+ // "SETSpacePoints" , Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::TrackerHitCollection> m_FTDSpacePointCol{
+ "FTDSpacePoints" , Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::TrackerHitCollection> m_VXDTrackerHits{
+ "VXDTrackerHits" , Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::TrackerHitCollection> m_SETTrackerHits{
+ "SETTrackerHits" , Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::TrackerHitCollection> m_SITTrackerHits{
+ "SITTrackerHits" , Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::TrackerHitCollection> m_FTDTrackerHits{
+ "FTDTrackerHits" , Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::SimTrackerHitCollection> m_VXDCollection{
+ "VXDCollection" , Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::SimTrackerHitCollection> m_SETCollection{
+ "SETCollection" , Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::SimTrackerHitCollection> m_SITCollection{
+ "SITCollection" , Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::SimTrackerHitCollection> m_FTDCollection{
+ "FTDCollection" , Gaudi::DataHandle::Reader, this};
//writer
DataHandle<edm4hep::TrackCollection> m_DCTrackCol{
"DCTrackCollection", Gaudi::DataHandle::Writer, this};
DataHandle<edm4hep::TrackCollection> m_SDTTrackCol{
"SDTTrackCollection", Gaudi::DataHandle::Writer, this};
- DataHandle<edm4hep::ReconstructedParticleCollection> m_DCRecParticleCol{
- "DCRecParticleCollection", Gaudi::DataHandle::Writer, this};
- DataHandle<edm4hep::MCRecoParticleAssociationCollection>
- m_DCRecParticleAssociationCol{
- "DCRecMCRecoParticleAssociationCollection",
- Gaudi::DataHandle::Writer, this};
+ DataHandle<edm4hep::TrackerHitCollection> m_truthTrackerHitCol{
+ "TruthTrackerHitCollection", Gaudi::DataHandle::Writer, this};
+
+ // Smear hit
+ DataHandle<edm4hep::SimTrackerHitCollection> w_SimSmearHCol{
+ "SmearSimHitsCollection", Gaudi::DataHandle::Writer, this};
+ DataHandle<edm4hep::TrackerHitCollection> m_SmeartruthTrackerHitCol{
+ "SmearTrackerHitCollection", Gaudi::DataHandle::Writer, this};
+ DataHandle<edm4hep::MCRecoTrackerAssociationCollection> w_SmearAssociationCol{
+ "SmearDCHitAssociationCollection", Gaudi::DataHandle::Writer, this};
+ // make noise hit
+ DataHandle<edm4hep::SimTrackerHitCollection> w_SimNoiseHCol{
+ "NoiseSimHitsCollection", Gaudi::DataHandle::Writer, this};
+ DataHandle<edm4hep::TrackerHitCollection> w_NoiseHitCol{
+ "NoiseDCHitsCollection", Gaudi::DataHandle::Writer, this};
+ DataHandle<edm4hep::MCRecoTrackerAssociationCollection> w_NoiseAssociationCol{
+ "NoiseDCHitAssociationCollection", Gaudi::DataHandle::Writer, this};
+
//readout for getting segmentation
Gaudi::Property<std::string> m_readout_name{this, "readout",
"DriftChamberHitsCollection"};
- Gaudi::Property<bool> m_writeRecParticle{this,"writeRecParticle",false};
+ Gaudi::Property<bool> m_hist{this,"hist",false};
+ Gaudi::Property<bool> m_useDC{this,"useDC",true};
+ Gaudi::Property<bool> m_useSi{this,"useSi",true};
+ Gaudi::Property<bool> m_useTruthTrack{this,"useTruthTrack",false};
+ Gaudi::Property<bool> m_useSiTruthHit{this,"useSiTruthHit",false};
+ Gaudi::Property<bool> m_skipSecondaryHit{this,"skipSecondaryHit",true};
+ Gaudi::Property<bool> m_useFirstHitForDC{this,"useFirstHitForDC",false};
+ Gaudi::Property<bool> m_useSiSpacePoint{this,"useSiSpacePoint",false};
+ Gaudi::Property<bool> m_useIdealHit{this,"useIdealHit",false};
+
+ Gaudi::Property<bool> m_useNoiseHits{this,"useNoiseHits",false};
+ Gaudi::Property<bool> m_smearHits{this,"smearHits",false};
+
+ Gaudi::Property<float> m_momentumCut{this,"momentumCut",0.1};//momentum cut for the first hit
Gaudi::Property<float> m_resPT{this,"resPT",0};//ratio
Gaudi::Property<float> m_resPz{this,"resPz",0};//ratio
+ Gaudi::Property<float> m_resX{this,"resX",0.11};//mm
+ Gaudi::Property<float> m_resY{this,"resY",0.11};//mm
+ Gaudi::Property<float> m_resZ{this,"resZ",0.11};//mm
+ Gaudi::Property<float> m_driftVelocity{this,"driftVelocity",40};// um/us
Gaudi::Property<float> m_resMomPhi{this,"resMomPhi",0};//radian
Gaudi::Property<float> m_resMomTheta{this,"resMomTheta",0};//radian
Gaudi::Property<float> m_resVertexX{this,"resVertexX",0.003};//3um
Gaudi::Property<float> m_resVertexY{this,"resVertexY",0.003};//3um
Gaudi::Property<float> m_resVertexZ{this,"resVertexZ",0.003};//3um
Gaudi::Property<int> m_maxDCDigiCut{this,"maxDigiCut",1e6};
+ Gaudi::Property<std::vector<float> > m_resVXD{this,"resVXD",{0.003,0.003,0.003}};//mm
+ Gaudi::Property<std::vector<float> > m_resSIT{this,"resSIT",{0.003,0.003,0.003}};//mm
+ Gaudi::Property<std::vector<float> > m_resSET{this,"resSET",{0.003,0.003,0.003}};//mm
+ Gaudi::Property<std::vector<float> > m_resFTDPixel{this,"resFTDPixel",{0.003,0.003,0.003}};//mm
+ Gaudi::Property<std::vector<float> > m_resFTDStrip{this,"resFTDStrip",{0.003,0.003,0.003}};//mm
+
+ Gaudi::Property<double> m_fHitPurity{this,"fHitPurity",0.1};
+ Gaudi::Property<float> m_pocaTime { this, "pocaTime", 225};// ns
+
+ double m_helixRadius,m_helixXC,m_helixYC;
+ double m_helixRadiusFirst,m_helixXCFirst,m_helixYCFirst;
+
+ NTuple::Tuple* m_tuple;
+ NTuple::Item<int> m_run;
+ NTuple::Item<int> m_evt;
+ NTuple::Array<double> m_siMom;
+ NTuple::Array<double> m_siPos;
+ NTuple::Array<double> m_mcMom;
+ NTuple::Array<double> m_mcPos;
+
+ NTuple::Item<int> m_nDCTrackHit;
+ NTuple::Item<int> m_nSmearDCTrackHit;
+ NTuple::Item<int> m_nNoiseDCTrackHit;
+ NTuple::Array<float> m_DriftDistance;
+ NTuple::Array<float> m_SmearDriftDistance;
+ NTuple::Array<float> m_NoiseDriftDistance;
+
+ NTuple::Item<int> m_nSimTrackerHitVXD;
+ NTuple::Item<int> m_nSimTrackerHitSIT;
+ NTuple::Item<int> m_nSimTrackerHitSET;
+ NTuple::Item<int> m_nSimTrackerHitFTD;
+ NTuple::Item<int> m_nSimTrackerHitDC;
+ NTuple::Item<int> m_nTrackerHitVXD;
+ NTuple::Item<int> m_nTrackerHitSIT;
+ NTuple::Item<int> m_nTrackerHitSET;
+ NTuple::Item<int> m_nTrackerHitFTD;
+ NTuple::Item<int> m_nTrackerHitDC;
+ NTuple::Item<int> m_nTrackerHitErrVXD;
+ NTuple::Item<int> m_nTrackerHitErrSIT;
+ NTuple::Item<int> m_nTrackerHitErrSET;
+ NTuple::Item<int> m_nTrackerHitErrFTD;
+ NTuple::Item<int> m_nSpacePointSIT;
+ NTuple::Item<int> m_nSpacePointSET;
+ NTuple::Item<int> m_nSpacePointFTD;
+ NTuple::Item<int> m_nSpacePointErrVXD;
+ NTuple::Item<int> m_nSpacePointErrSIT;
+ NTuple::Item<int> m_nSpacePointErrSET;
+ NTuple::Item<int> m_nSpacePointErrFTD;
+ NTuple::Item<int> m_nHitOnSiTkVXD;
+ NTuple::Item<int> m_nHitOnSiTkSIT;
+ NTuple::Item<int> m_nHitOnSiTkSET;
+ NTuple::Item<int> m_nHitOnSiTkFTD;
+ NTuple::Item<int> m_nHitOnSdtTkVXD;
+ NTuple::Item<int> m_nHitOnSdtTkSIT;
+ NTuple::Item<int> m_nHitOnSdtTkSET;
+ NTuple::Item<int> m_nHitOnSdtTkFTD;
+ NTuple::Item<int> m_nHitOnSdtTkDC;
+ NTuple::Item<int> m_nHitOnSdtTk;
+ NTuple::Item<int> m_nNoiseOnSdtTk;
+
+ TRandom3 fRandom;
};
#endif
diff --git a/Utilities/DataHelper/CMakeLists.txt b/Utilities/DataHelper/CMakeLists.txt
index 8398a6557a29db9dd46c8cb5b8518ae6ff7ecc3e..74ec620fb6ddf31c1795ac9f0b99f7b386f1b375 100644
--- a/Utilities/DataHelper/CMakeLists.txt
+++ b/Utilities/DataHelper/CMakeLists.txt
@@ -15,9 +15,11 @@ gaudi_add_library(DataHelperLib
src/TrackerHitExtended.cc
src/TrackExtended.cc
src/TrackHitPair.cc
- src/TrackerHitHelper.cpp
+ src/TrackerHitHelper.cpp
+ src/TrackHelper.cc
LINK EDM4HEP::edm4hep
EDM4HEP::edm4hepDict
+ DetSegmentation
GSL::gsl
${CLHEP_LIBRARIES}
Identifier
diff --git a/Utilities/DataHelper/include/DataHelper/HelixClass.h b/Utilities/DataHelper/include/DataHelper/HelixClass.h
index 794f8158652a637180cb99de38e35a249c2142d9..f573ee9d519e759130c86ca935606f57eeb74371 100644
--- a/Utilities/DataHelper/include/DataHelper/HelixClass.h
+++ b/Utilities/DataHelper/include/DataHelper/HelixClass.h
@@ -69,6 +69,7 @@ class HelixClass {
* - particle charge : q;<br>
* - magnetic field (in Tesla) : B<br>
*/
+ void Initialize_VP(double * pos, double * mom, double q, double B);
void Initialize_VP(float * pos, float * mom, float q, float B);
/**
@@ -266,35 +267,35 @@ class HelixClass {
float getCharge();
private:
- float _momentum[3]; // momentum @ ref point
+ float _momentum[3]; // momentum @ ref point
float _referencePoint[3]; // coordinates @ ref point
- float _phi0; // phi0 in canonical parameterization
- float _d0; // d0 in canonical parameterisation
- float _z0; // z0 in canonical parameterisation
- float _omega; // signed curvuture in canonical parameterisation
- float _tanLambda; // TanLambda
- float _pxy; // Transverse momentum
- float _charge; // Particle Charge
- float _bField; // Magnetic field (assumed to point to Z>0)
- float _radius; // radius of circle in XY plane
- float _xCentre; // X of circle centre
- float _yCentre; // Y of circle centre
- float _phiRefPoint; // Phi w.r.t. (X0,Y0) of circle @ ref point
- float _phiAtPCA; // Phi w.r.t. (X0,Y0) of circle @ PCA
- float _xAtPCA; // X @ PCA
- float _yAtPCA; // Y @ PCA
- float _pxAtPCA; // PX @ PCA
- float _pyAtPCA; // PY @ PCA
- float _phiMomRefPoint; // Phi of Momentum vector @ ref point
- float _const_pi; // PI
- float _const_2pi; // 2*PI
- float _const_pi2; // PI/2
- float _FCT; // 2.99792458E-4
+ double _phi0; // phi0 in canonical parameterization
+ double _d0; // d0 in canonical parameterisation
+ double _z0; // z0 in canonical parameterisation
+ double _omega; // signed curvuture in canonical parameterisation
+ double _tanLambda; // TanLambda
+ double _pxy; // Transverse momentum
+ double _charge; // Particle Charge
+ double _bField; // Magnetic field (assumed to point to Z>0)
+ double _radius; // radius of circle in XY plane
+ double _xCentre; // X of circle centre
+ double _yCentre; // Y of circle centre
+ double _phiRefPoint; // Phi w.r.t. (X0,Y0) of circle @ ref point
+ double _phiAtPCA; // Phi w.r.t. (X0,Y0) of circle @ PCA
+ double _xAtPCA; // X @ PCA
+ double _yAtPCA; // Y @ PCA
+ double _pxAtPCA; // PX @ PCA
+ double _pyAtPCA; // PY @ PCA
+ double _phiMomRefPoint; // Phi of Momentum vector @ ref point
+ double _const_pi; // PI
+ double _const_2pi; // 2*PI
+ double _const_pi2; // PI/2
+ double _FCT; // 2.99792458E-4
float _xStart[3]; // Starting point of track segment
float _xEnd[3]; // Ending point of track segment
- float _bZ;
- float _phiZ;
+ double _bZ;
+ double _phiZ;
};
diff --git a/Utilities/DataHelper/src/HelixClass.cc b/Utilities/DataHelper/src/HelixClass.cc
index 591acc5148ce71b84eb79421cd5c841f423a5f2f..a9ef6ca7092f6231d9914f423666e743e6dd3ad3 100644
--- a/Utilities/DataHelper/src/HelixClass.cc
+++ b/Utilities/DataHelper/src/HelixClass.cc
@@ -12,6 +12,75 @@ HelixClass::HelixClass() {
HelixClass::~HelixClass() {}
+void HelixClass::Initialize_VP(double * pos, double * mom, double q, double B) {
+ _referencePoint[0] = (float) pos[0];
+ _referencePoint[1] = (float) pos[1];
+ _referencePoint[2] = (float) pos[2];
+ _momentum[0] = (float) mom[0];
+ _momentum[1] = (float) mom[1];
+ _momentum[2] = (float) mom[2];
+ _charge = q;
+ _bField = B;
+ _pxy = sqrt(mom[0]*mom[0]+mom[1]*mom[1]);
+ _radius = _pxy / (_FCT*B);
+ _omega = q/_radius;
+ _tanLambda = mom[2]/_pxy;
+ _phiMomRefPoint = atan2(mom[1],mom[0]);
+ _xCentre = pos[0] + _radius*cos(_phiMomRefPoint-_const_pi2*q);
+ _yCentre = pos[1] + _radius*sin(_phiMomRefPoint-_const_pi2*q);
+ _phiRefPoint = atan2(pos[1]-_yCentre,pos[0]-_xCentre);
+ _phiAtPCA = atan2(-_yCentre,-_xCentre);
+ _phi0 = -_const_pi2*q + _phiAtPCA;
+ while (_phi0<0) _phi0+=_const_2pi;
+ while (_phi0>=_const_2pi) _phi0-=_const_2pi;
+ _xAtPCA = _xCentre + _radius*cos(_phiAtPCA);
+ _yAtPCA = _yCentre + _radius*sin(_phiAtPCA);
+ // _d0 = -_xAtPCA*sin(_phi0) + _yAtPCA*cos(_phi0);
+
+ if (q>0) {
+ _d0 = double(q)*_radius - double(sqrt(_xCentre*_xCentre+_yCentre*_yCentre));
+ }
+ else {
+ _d0 = double(q)*_radius + double(sqrt(_xCentre*_xCentre+_yCentre*_yCentre));
+ }
+
+// if (fabs(_d0)>0.001 ) {
+// std::cout << "New helix : " << std::endl;
+// std::cout << " Position : " << pos[0]
+// << " " << pos[1]
+// << " " << pos[2] << std::endl;
+// std::cout << " Radius = " << _radius << std::endl;
+// std::cout << " RC = " << sqrt(_xCentre*_xCentre+_yCentre*_yCentre) << std::endl;
+// std::cout << " D0 = " << _d0 << std::endl;
+// }
+
+ _pxAtPCA = _pxy*cos(_phi0);
+ _pyAtPCA = _pxy*sin(_phi0);
+ float deltaPhi = _phiRefPoint - _phiAtPCA;
+ float xCircles = -pos[2]*q/(_radius*_tanLambda) - deltaPhi;
+ xCircles = xCircles/_const_2pi;
+ int nCircles;
+ int n1,n2;
+
+ if (xCircles >= 0.) {
+ n1 = int(xCircles);
+ n2 = n1 + 1;
+ }
+ else {
+ n1 = int(xCircles) - 1;
+ n2 = n1 + 1;
+ }
+
+ if (fabs(n1-xCircles) < fabs(n2-xCircles)) {
+ nCircles = n1;
+ }
+ else {
+ nCircles = n2;
+ }
+ _z0 = pos[2] + _radius*_tanLambda*q*(deltaPhi + _const_2pi*nCircles);
+
+}
+
void HelixClass::Initialize_VP(float * pos, float * mom, float q, float B) {
_referencePoint[0] = pos[0];
_referencePoint[1] = pos[1];
diff --git a/Utilities/DataHelper/src/TrackerHitHelper.cpp b/Utilities/DataHelper/src/TrackerHitHelper.cpp
index e12eb886c5aa7863f320b1157d278948cf21c880..f12b16873bd1d553023c684510ca56d21a7e0246 100644
--- a/Utilities/DataHelper/src/TrackerHitHelper.cpp
+++ b/Utilities/DataHelper/src/TrackerHitHelper.cpp
@@ -7,6 +7,8 @@
#include "CLHEP/Vector/ThreeVector.h"
#include "CLHEP/Vector/Rotation.h"
#include <bitset>
+#include "TVector3.h"
+#include "DD4hep/DD4hepUnits.h"
std::array<float,6> CEPC::GetCovMatrix(edm4hep::TrackerHit& hit, bool useSpacePointBuilderMethod){
if(hit.isAvailable()){
@@ -115,3 +117,44 @@ std::array<float, 6> CEPC::ConvertToCovXYZ(float dU, float thetaU, float phiU, f
}
return cov;
}
+
+const edm4hep::SimTrackerHit CEPC::getAssoClosestSimTrackerHit(
+ const edm4hep::MCRecoTrackerAssociationCollection* assoHits,
+ const edm4hep::TrackerHit trackerHit,
+ const dd4hep::DDSegmentation::GridDriftChamber* segmentation,
+ int docaMehtod)
+{
+ std::vector<edm4hep::SimTrackerHit> hits;
+ for(auto assoHit: *assoHits){
+ if(assoHit.getRec()==trackerHit){
+ hits.push_back(assoHit.getSim());
+ }
+ }
+ edm4hep::SimTrackerHit minSimTrackerHit;
+ double min_distance = 999 ;
+ double tmp_distance =0.;
+ for(auto hit:hits){
+ unsigned long long wcellid=hit.getCellID();
+ TVector3 pos(hit.getPosition()[0]*dd4hep::mm, hit.getPosition()[1]*dd4hep::mm, hit.getPosition()[2]*dd4hep::mm);
+
+ TVector3 sim_mon(hit.getMomentum()[0],hit.getMomentum()[1],hit.getMomentum()[2]);
+ float Steplength = hit.getPathLength();
+ TVector3 pos_start=pos - 0.5 * Steplength * sim_mon.Unit();
+ TVector3 pos_end=pos + 0.5 * Steplength * sim_mon.Unit();
+ TVector3 hitPosition;
+ TVector3 PCA;
+ tmp_distance = segmentation->Distance(wcellid,pos_start,pos_end,hitPosition,PCA);
+ tmp_distance = tmp_distance/dd4hep::mm; //mm
+
+ if(tmp_distance < min_distance){
+ min_distance = tmp_distance;
+ //pos_x = hitPosition.x(); //pos.x();
+ //pos_y = hitPosition.y(); //pos.y();
+ //pos_z = pos.z();
+ //momx = hit.getMomentum()[0];
+ //momy = hit.getMomentum()[1];
+ minSimTrackerHit=hit;
+ }
+ }
+ return minSimTrackerHit;
+}