From aad34dca827f8eca8b81512d8afc4c2d1201e81d Mon Sep 17 00:00:00 2001
From: Zhang Yao <zhangyao@ihep.ac.cn>
Date: Tue, 5 Jan 2021 11:15:50 +0800
Subject: [PATCH] Update TruthTrackerAlg for CRD. Add job options and compact
file for CRD fitting.
---
.../compact/CRD_o1_v01/CRD_o1_v01_noECAL.xml | 51 ++++
.../src/TruthTracker/TruthTrackerAlg.cpp | 233 ++++++++++-----
.../src/TruthTracker/TruthTrackerAlg.h | 38 ++-
sim_fit_CRD.py | 279 ++++++++++++++++++
4 files changed, 514 insertions(+), 87 deletions(-)
create mode 100644 Detector/DetCRD/compact/CRD_o1_v01/CRD_o1_v01_noECAL.xml
create mode 100644 sim_fit_CRD.py
diff --git a/Detector/DetCRD/compact/CRD_o1_v01/CRD_o1_v01_noECAL.xml b/Detector/DetCRD/compact/CRD_o1_v01/CRD_o1_v01_noECAL.xml
new file mode 100644
index 00000000..256dcf00
--- /dev/null
+++ b/Detector/DetCRD/compact/CRD_o1_v01/CRD_o1_v01_noECAL.xml
@@ -0,0 +1,51 @@
+<lccdd xmlns:compact="http://www.lcsim.org/schemas/compact/1.0"
+ xmlns:xs="http://www.w3.org/2001/XMLSchema"
+ xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/compact/1.0/compact.xsd">
+ <info name="CRD_i1_v01"
+ title="CepC reference detctor with coil inside Hcal"
+ author="C.D.Fu, "
+ url="http://cepc.ihep.ac.cn"
+ status="developing"
+ version="v01">
+ <comment>CepC reference detector simulation models used for detector study </comment>
+ </info>
+
+ <includes>
+ <gdmlFile ref="${DD4hepINSTALL}/DDDetectors/compact/elements.xml"/>
+ <gdmlFile ref="../CRD_common_v01/materials.xml"/>
+ </includes>
+
+ <define>
+ <constant name="world_size" value="2.6*m"/>
+ <constant name="world_x" value="world_size"/>
+ <constant name="world_y" value="world_size"/>
+ <constant name="world_z" value="world_size"/>
+
+ <include ref="${DD4hepINSTALL}/DDDetectors/compact/detector_types.xml"/>
+ </define>
+
+ <include ref="./CRD_Dimensions_v01_01.xml"/>
+
+ <include ref="../CRD_common_v01/Beampipe_v01_01.xml"/>
+ <include ref="../CRD_common_v01/VXD_v01_01.xml"/>
+ <include ref="../CRD_common_v01/FTD_SimpleStaggered_v01_01.xml"/>
+ <include ref="../CRD_common_v01/SIT_SimplePlanar_v01_01.xml"/>
+ <include ref="../CRD_common_v01/DC_Simple_v01_01.xml"/>
+ <include ref="../CRD_common_v01/SET_SimplePlanar_v01_01.xml"/>
+ <!-- <include ref="../CRD_common_v01/Ecal_Crystal_Barrel_v01_01.xml"/> -->
+ <!--include ref="../CRD_common_v01/Ecal_Crystal_Endcap_v01_01.xml"/-->
+ <!--include ref="../CRD_common_v01/Coil_v01_01.xml"/-->
+ <!--include ref="../CRD_common_v01/Hcal_v01_01.xml"/-->
+ <!--include ref="../CRD_common_v01/Yoke_v01_01.xml"/-->
+ <!--include ref="../CRD_common_v01/Lcal_v01_01.xml"/-->
+
+ <fields>
+ <field name="GlobalSolenoid" type="solenoid"
+ inner_field="Field_nominal_value"
+ outer_field="Field_outer_nominal_value"
+ zmax="SolenoidCoil_half_length"
+ outer_radius="SolenoidCoil_center_radius">
+ </field>
+ </fields>
+
+</lccdd>
diff --git a/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp b/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp
index 6a9842e0..9e24a601 100644
--- a/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp
+++ b/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.cpp
@@ -1,5 +1,15 @@
#include "TruthTrackerAlg.h"
#include "DataHelper/HelixClass.h"
+#include "DDSegmentation/Segmentation.h"
+#include "DetSegmentation/GridDriftChamber.h"
+#include "DetInterface/IGeomSvc.h"
+#include "DD4hep/Detector.h"
+#include "DD4hep/IDDescriptor.h"
+#include "DD4hep/Plugins.h"
+#include "DD4hep/DD4hepUnits.h"
+
+//external
+#include "CLHEP/Random/RandGauss.h"
#include "edm4hep/EventHeaderCollection.h"
#include "edm4hep/MCParticleCollection.h"
#include "edm4hep/SimTrackerHitCollection.h"
@@ -8,37 +18,30 @@
#include "edm4hep/MCRecoTrackerAssociationCollection.h"
#include "edm4hep/MCRecoParticleAssociationCollection.h"
#include "edm4hep/ReconstructedParticleCollection.h"
-#include "DDSegmentation/Segmentation.h"
-#include "DetSegmentation/GridDriftChamber.h"
-#include "DetInterface/IGeomSvc.h"
-//#include "edm4hep/SimCalorimeterHitCollection.h"
-//#include "edm4hep/CaloHitContributionCollection.h"
-
-#include "DD4hep/Detector.h"
-#include "DD4hep/IDDescriptor.h"
-#include "DD4hep/Plugins.h"
-#include "DD4hep/DD4hepUnits.h"
DECLARE_COMPONENT(TruthTrackerAlg)
TruthTrackerAlg::TruthTrackerAlg(const std::string& name, ISvcLocator* svcLoc)
-: GaudiAlgorithm(name,svcLoc),m_dd4hep(nullptr),m_gridDriftChamber(nullptr),
+: GaudiAlgorithm(name, svcLoc),m_dd4hep(nullptr),m_gridDriftChamber(nullptr),
m_decoder(nullptr)
{
declareProperty("MCParticle", m_mcParticleCol,
"Handle of the input MCParticle collection");
- declareProperty("DigiDCHitCollection", m_digiDCHitsCol,
- "Handle of digi DCHit collection");
- declareProperty("DCHitAssociationCollection", m_dcHitAssociationCol,
- "Handle of association collection");
- declareProperty("DCTrackCollection", m_dcTrackCol,
- "Handle of association collection");
- declareProperty("DCRecParticleCollection", m_dcRecParticleCol,
+ 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,
+ "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,
+ m_DCRecParticleAssociationCol,
"Handle of drift chamber reconstructed particle collection");
- declareProperty("debug", m_debug=false);
}
StatusCode TruthTrackerAlg::initialize()
@@ -87,17 +90,22 @@ StatusCode TruthTrackerAlg::execute()
}
///Retrieve DC digi
const edm4hep::TrackerHitCollection* digiDCHitsCol=nullptr;
- digiDCHitsCol=m_digiDCHitsCol.get();//FIXME DEBUG
+ digiDCHitsCol=m_DCDigiCol.get();//FIXME DEBUG
if(nullptr==digiDCHitsCol){
debug()<<"TrackerHitCollection not found"<<endmsg;
- //return StatusCode::SUCCESS;
+ //return StatusCode::SUCCESS;//FIXME return when no hits in DC + silicon
}
if((int) digiDCHitsCol->size()>m_maxDCDigiCut) return StatusCode::SUCCESS;
///Output Track collection
- edm4hep::TrackCollection* dcTrackCol=m_dcTrackCol.createAndPut();
+ edm4hep::TrackCollection* dcTrackCol=m_DCTrackCol.createAndPut();
+ edm4hep::TrackCollection* sdtTrackCol=m_SDTTrackCol.createAndPut();
+ edm4hep::ReconstructedParticleCollection* dcRecParticleCol(nullptr);
+ if(m_writeRecParticle){
+ dcRecParticleCol=m_DCRecParticleCol.createAndPut();
+ }
////TODO
- /////Output MCRecoTrackerAssociationCollection collection
+ //Output MCRecoTrackerAssociationCollection collection
//const edm4hep::MCRecoTrackerAssociationCollection*
// mcRecoTrackerAssociationCol=nullptr;
//if(nullptr==mcRecoTrackerAssociationCol){
@@ -107,29 +115,85 @@ StatusCode TruthTrackerAlg::execute()
//}
//mcRecoTrackerAssociationCol=m_mcRecoParticleAssociation.get();
- ///Convert MCParticle to Track and ReconstructedParticle
- if(m_debug){
- debug()<<"MCParticleCol size="<<mcParticleCol->size()<<endmsg;
+ ///Retrieve silicon Track
+ const edm4hep::TrackCollection* siTrackCol=nullptr;
+ siTrackCol=m_siSubsetTrackCol.get();
+
+ ///New SDT track
+ for(auto siTrack:*siTrackCol){
+ edm4hep::Track sdtTrack=sdtTrackCol->create();
+ edm4hep::TrackState sdtTrackState;
+ edm4hep::TrackState siTrackStat=siTrack.getTrackStates(0);//FIXME?
+ sdtTrackState.location=siTrackStat.location;
+ sdtTrackState.D0=siTrackStat.D0;
+ sdtTrackState.phi=siTrackStat.phi;
+ sdtTrackState.omega=siTrackStat.omega;
+ sdtTrackState.Z0=siTrackStat.Z0;
+ sdtTrackState.tanLambda=siTrackStat.tanLambda;
+ sdtTrackState.referencePoint=siTrackStat.referencePoint;
+ for(int k=0;k<15;k++){
+ sdtTrackState.covMatrix[k]=siTrackStat.covMatrix[k];
+ }
+ sdtTrack.addToTrackStates(sdtTrackState);
+ sdtTrack.setType(siTrack.getType());
+ sdtTrack.setChi2(siTrack.getChi2());
+ sdtTrack.setNdf(siTrack.getNdf());
+ sdtTrack.setDEdx(siTrack.getDEdx());
+ sdtTrack.setDEdxError(siTrack.getDEdxError());
+ sdtTrack.setRadiusOfInnermostHit(siTrack.getRadiusOfInnermostHit());
+ debug()<<"siTrack trackerHits_size="<<siTrack.trackerHits_size()<<endmsg;
+ for(unsigned int iSiTackerHit=0;iSiTackerHit<siTrack.trackerHits_size();
+ iSiTackerHit++){
+ sdtTrack.addToTrackerHits(siTrack.getTrackerHits(iSiTackerHit));
+ }
+ //TODO tracks
+ for(auto digiDC:*digiDCHitsCol){
+ //if(Sim->MCParti!=current) continue;//TODO
+ sdtTrack.addToTrackerHits(digiDC);
+ }
}
+
+ ///Convert MCParticle to Track and ReconstructedParticle
+ debug()<<"MCParticleCol size="<<mcParticleCol->size()<<endmsg;
for(auto mcParticle : *mcParticleCol){
- //if(fabs(mcParticle.getCharge()<1e-6) continue;//Skip neutral particles
- edm4hep::Vector3d posV= mcParticle.getVertex();
- edm4hep::Vector3f momV= mcParticle.getMomentum();//GeV
- float pos[3]={(float) posV.x, (float) posV.y, (float) posV.z};
- float mom[3]={momV.x, momV.y, momV.z};
- //FIXME
- //pivotToFirstLayer(mcPocaPos,mcPocaMom,firstLayerPos,firstLayerMom);
+ /// skip mcParticleVertex do not have enough associated hits TODO
+
+ ///Vertex
+ const edm4hep::Vector3d mcParticleVertex=mcParticle.getVertex();//mm
+ edm4hep::Vector3f mcParticleVertexSmeared;//mm
+ mcParticleVertexSmeared.x=
+ CLHEP::RandGauss::shoot(mcParticleVertex.x,m_resVertexX);
+ mcParticleVertexSmeared.y=
+ CLHEP::RandGauss::shoot(mcParticleVertex.y,m_resVertexY);
+ mcParticleVertexSmeared.z=
+ CLHEP::RandGauss::shoot(mcParticleVertex.z,m_resVertexZ);
+ ///Momentum
+ edm4hep::Vector3f mcParticleMom=mcParticle.getMomentum();//GeV
+ double mcParticlePt=sqrt(mcParticleMom.x*mcParticleMom.x+
+ mcParticleMom.y*mcParticleMom.y);
+ //double mcParticlePtSmeared=
+ // CLHEP::RandGauss::shoot(mcParticlePt,m_resPT);
+ double mcParticleMomPhi=atan2(mcParticleMom.y,mcParticleMom.x);
+ double mcParticleMomPhiSmeared=
+ CLHEP::RandGauss::shoot(mcParticleMomPhi,m_resMomPhi);
+ edm4hep::Vector3f mcParticleMomSmeared;
+ mcParticleMomSmeared.x=mcParticlePt*cos(mcParticleMomPhiSmeared);
+ mcParticleMomSmeared.y=mcParticlePt*sin(mcParticleMomPhiSmeared);
+ mcParticleMomSmeared.z=
+ CLHEP::RandGauss::shoot(mcParticleMom.z,m_resPz);
+
+ ///Converted to Helix
double B[3]={1e9,1e9,1e9};
m_dd4hepField.magneticField({0.,0.,0.},B);
- if(m_debug)std::cout<<" PDG "<<mcParticle.getPDG()<<" charge "
- << mcParticle.getCharge()
- <<" pos "<<pos[0]<<" "<<pos[1]<<" "<<pos[2]
- <<" mom "<<mom[0]<<" "<<mom[1]<<" "<<mom[2]
- <<" Bxyz "<<B[0]/dd4hep::tesla<<" "<<B[1]/dd4hep::tesla
- <<" "<<B[2]/dd4hep::tesla<<" tesla"<<std::endl;
-
HelixClass helix;
- helix.Initialize_VP(pos,mom,(int) mcParticle.getCharge(),B[2]/dd4hep::tesla);
+ float pos[3]={mcParticleVertexSmeared.x,
+ mcParticleVertexSmeared.y,mcParticleVertexSmeared.z};
+ float mom[3]={mcParticleMomSmeared.x,mcParticleMomSmeared.y,
+ mcParticleMomSmeared.z};
+ helix.Initialize_VP(pos,mom,mcParticle.getCharge(),B[2]/dd4hep::tesla);
+
+ ///new Track
+ edm4hep::Track dcTrack=dcTrackCol->create();
edm4hep::TrackState trackState;
trackState.D0=helix.getD0();
trackState.phi=helix.getPhi0();
@@ -140,40 +204,61 @@ StatusCode TruthTrackerAlg::execute()
std::array<float,15> covMatrix;
for(int i=0;i<15;i++){covMatrix[i]=999.;}//FIXME
trackState.covMatrix=covMatrix;
-
- if(m_debug){
- std::cout<<" helix d0 "<<trackState.D0<<
- " phi "<<trackState.phi<<
- " omega "<<trackState.omega<<
- " z0 "<<trackState.Z0<<
- " tanLambda "<<trackState.tanLambda<<
- " referencePoint "<<trackState.referencePoint<< std::endl;
- }
- //new Track
- edm4hep::Track track = dcTrackCol->create();
- track.addToTrackStates(trackState);
- if(nullptr!=digiDCHitsCol) {
- //track.setType();//TODO
- //track.setChi2(trackState->getChi2());
- track.setNdf(digiDCHitsCol->size()-5);
- //track.setDEdx();//TODO
- //track.setRadiusOfInnermostHit();//TODO
- for(unsigned int i=0; i<digiDCHitsCol->size(); i++ ){
- edm4hep::TrackerHit digiDC=digiDCHitsCol->at(i);
- //if(Sim->MCParti!=current) continue;//TODO
- track.addToTrackerHits(digiDC);
- }
+ 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()<<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);
}
- //TODO
- //new ReconstructedParticle
- //recParticle->setType();
- //dcRecParticle->setEnergy();
- //edm4hep::Vector3f momVec3(helix.getMomentum()[0],
- // helix.getMomentum()[1],helix.getMomentum()[2]);
- //recParticle->setMomentum(momVec3);
- }
+ dcTrack.setRadiusOfInnermostHit(radiusOfInnermostHit);//TODO
+
+ edm4hep::ReconstructedParticle 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
+
+ debug()<<"mcParticle "<<mcParticle
+ <<" momPhi "<<mcParticleMomPhi
+ <<" mcParticleVertex("<<mcParticleVertex<<")mm "
+ <<" mcParticleVertexSmeared("<<mcParticleVertexSmeared<<")mm "
+ <<" mcParticleMom("<<mcParticleMom<<")GeV "
+ <<" 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
- if(m_debug) std::cout<<"Output DCTrack size="<<dcTrackCol->size()<<std::endl;
+ debug()<<"Output DCTrack size="<<dcTrackCol->size()<<endmsg;
return StatusCode::SUCCESS;
}
diff --git a/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.h b/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.h
index 36a87d65..30efd0a8 100644
--- a/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.h
+++ b/Reconstruction/Tracking/src/TruthTracker/TruthTrackerAlg.h
@@ -27,9 +27,9 @@ class TruthTrackerAlg: public GaudiAlgorithm
public:
TruthTrackerAlg(const std::string& name, ISvcLocator* svcLoc);
- virtual StatusCode initialize();
- virtual StatusCode execute();
- virtual StatusCode finalize();
+ virtual StatusCode initialize() override;
+ virtual StatusCode execute() override;
+ virtual StatusCode finalize() override;
private:
SmartIF<IGeomSvc> m_geomSvc;
@@ -41,26 +41,38 @@ class TruthTrackerAlg: public GaudiAlgorithm
//reader
DataHandle<edm4hep::MCParticleCollection> m_mcParticleCol{
"MCParticle", Gaudi::DataHandle::Reader, this};
- DataHandle<edm4hep::TrackerHitCollection> m_digiDCHitsCol{
- "DigiDCHitsCollection", Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::TrackerHitCollection> m_DCDigiCol{
+ "DigiDCHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<edm4hep::MCRecoTrackerAssociationCollection>
- m_dcHitAssociationCol{ "DCHitAssociationCollection",
+ m_DCHitAssociationCol{ "DCHitAssociationCollection",
+ Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::TrackCollection>
+ m_siSubsetTrackCol{ "SiSubsetTrackCollection",
Gaudi::DataHandle::Reader, this};
//writer
- DataHandle<edm4hep::TrackCollection> m_dcTrackCol{"DCTrackCollection",
- Gaudi::DataHandle::Writer, this};
- DataHandle<edm4hep::ReconstructedParticleCollection> m_dcRecParticleCol{
+ 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",
+ m_DCRecParticleAssociationCol{
+ "DCRecMCRecoParticleAssociationCollection",
Gaudi::DataHandle::Writer, this};
//readout for getting segmentation
Gaudi::Property<std::string> m_readout_name{this, "readout",
"DriftChamberHitsCollection"};
-
- Gaudi::Property<int> m_debug{this, "debug", false};
- Gaudi::Property<int> m_maxDCDigiCut{this, "maxDCDigiCut",150};
+ Gaudi::Property<bool> m_writeRecParticle{this,"writeRecParticle",false};
+ Gaudi::Property<float> m_resPT{this,"resPT",0};//ratio
+ Gaudi::Property<float> m_resPz{this,"resPz",0};//ratio
+ 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};
};
#endif
diff --git a/sim_fit_CRD.py b/sim_fit_CRD.py
new file mode 100644
index 00000000..93c8615f
--- /dev/null
+++ b/sim_fit_CRD.py
@@ -0,0 +1,279 @@
+#!/usr/bin/env python
+from Gaudi.Configuration import *
+
+##############################################################################
+# Event service
+##############################################################################
+from Configurables import k4DataSvc
+dsvc = k4DataSvc("EventDataSvc")
+
+##############################################################################
+# Random seed
+##############################################################################
+from Configurables import RndmGenSvc, HepRndm__Engine_CLHEP__RanluxEngine_
+# rndmengine = HepRndm__Engine_CLHEP__RanluxEngine_() # The default engine in Gaudi
+rndmengine = HepRndm__Engine_CLHEP__HepJamesRandom_() # The default engine in Geant4
+rndmengine.SetSingleton = True
+rndmengine.Seeds = [10]
+
+from Configurables import MarlinEvtSeeder
+evtseeder = MarlinEvtSeeder("EventSeeder")
+
+##############################################################################
+# Detector geometry
+##############################################################################
+geometry_option = "CRD_o1_v01/CRD_o1_v01_noECAL.xml"
+
+if not os.getenv("DETCRDROOT"):
+ print("Can't find the geometry. Please setup envvar DETCRDROOT." )
+ sys.exit(-1)
+
+geometry_path = os.path.join(os.getenv("DETCRDROOT"), "compact", geometry_option)
+if not os.path.exists(geometry_path):
+ print("Can't find the compact geometry file: %s"%geometry_path)
+ sys.exit(-1)
+
+from Configurables import GeomSvc
+geosvc = GeomSvc("GeomSvc")
+geosvc.compact = geometry_path
+
+from Configurables import GearSvc
+gearsvc = GearSvc("GearSvc")
+
+##############################################################################
+# Physics Generator
+##############################################################################
+from Configurables import GenAlgo
+from Configurables import GtGunTool
+from Configurables import StdHepRdr
+from Configurables import SLCIORdr
+from Configurables import HepMCRdr
+from Configurables import GenPrinter
+gun = GtGunTool("GtGunTool")
+gun.Particles = ["mu-"]
+gun.EnergyMins = [100.] # GeV
+gun.EnergyMaxs = [100.] # GeV
+gun.ThetaMins = [0] # deg
+gun.ThetaMaxs = [180] # deg
+gun.PhiMins = [0] # deg
+gun.PhiMaxs = [360] # deg
+# stdheprdr = StdHepRdr("StdHepRdr")
+# stdheprdr.Input = "/cefs/data/stdhep/CEPC250/2fermions/E250.Pbhabha.e0.p0.whizard195/bhabha.e0.p0.00001.stdhep"
+# lciordr = SLCIORdr("SLCIORdr")
+# lciordr.Input = "/cefs/data/stdhep/lcio250/signal/Higgs/E250.Pbbh.whizard195/E250.Pbbh_X.e0.p0.whizard195/Pbbh_X.e0.p0.00001.slcio"
+# hepmcrdr = HepMCRdr("HepMCRdr")
+# hepmcrdr.Input = "example_UsingIterators.txt"
+
+genprinter = GenPrinter("GenPrinter")
+
+genalg = GenAlgo("GenAlgo")
+genalg.GenTools = ["GtGunTool"]
+#genalg.GenTools = ["StdHepRdr"]
+# genalg.GenTools = ["StdHepRdr", "GenPrinter"]
+# genalg.GenTools = ["SLCIORdr", "GenPrinter"]
+# genalg.GenTools = ["HepMCRdr", "GenPrinter"]
+
+##############################################################################
+# Detector Simulation
+##############################################################################
+from Configurables import DetSimSvc
+detsimsvc = DetSimSvc("DetSimSvc")
+
+from Configurables import DetSimAlg
+detsimalg = DetSimAlg("DetSimAlg")
+#detsimalg.VisMacs = ["vis.mac"]
+detsimalg.RunCmds = [
+# "/tracking/verbose 1",
+]
+detsimalg.AnaElems = [
+ # example_anatool.name()
+ # "ExampleAnaElemTool"
+ "Edm4hepWriterAnaElemTool"
+]
+detsimalg.RootDetElem = "WorldDetElemTool"
+
+from Configurables import AnExampleDetElemTool
+example_dettool = AnExampleDetElemTool("AnExampleDetElemTool")
+
+from Configurables import CalorimeterSensDetTool
+from Configurables import DriftChamberSensDetTool
+
+calo_sensdettool = CalorimeterSensDetTool("CalorimeterSensDetTool")
+driftchamber_sensdettool = DriftChamberSensDetTool("DriftChamberSensDetTool")
+
+##############################################################################
+# DedxSimTool
+##############################################################################
+# dedxoption = "DummyDedxSimTool"
+dedxoption = "BetheBlochEquationDedxSimTool"
+
+driftchamber_sensdettool.DedxSimTool = dedxoption
+
+from Configurables import DummyDedxSimTool
+from Configurables import BetheBlochEquationDedxSimTool
+
+if dedxoption == "DummyDedxSimTool":
+ dedx_simtool = DummyDedxSimTool("DummyDedxSimTool")
+elif dedxoption == "BetheBlochEquationDedxSimTool":
+ dedx_simtool = BetheBlochEquationDedxSimTool("BetheBlochEquationDedxSimTool")
+ dedx_simtool.material_Z = 2
+ dedx_simtool.material_A = 4
+ dedx_simtool.scale = 10
+ dedx_simtool.resolution = 0.0001
+
+##############################################################################
+# DCHDigiAlg
+##############################################################################
+from Configurables import DCHDigiAlg
+dCHDigiAlg = DCHDigiAlg("DCHDigiAlg")
+dCHDigiAlg.readout = "DriftChamberHitsCollection"
+dCHDigiAlg.drift_velocity = 40#um/ns
+dCHDigiAlg.mom_threshold = 0 #GeV
+dCHDigiAlg.WriteAna = True
+
+##############################################################################
+# PODIO
+##############################################################################
+from Configurables import PodioOutput
+out = PodioOutput("outputalg")
+out.filename = "CRD-o1-v01-sim00.root"
+out.outputCommands = ["keep *"]
+
+
+##############################################################################
+# Silicon digitization
+##############################################################################
+vxdhitname = "VXDTrackerHits"
+sithitname = "SITTrackerHits"
+sitspname = "SITSpacePoints"
+tpchitname = "TPCTrackerHits"
+sethitname = "SETTrackerHits"
+setspname = "SETSpacePoints"
+ftdspname = "FTDSpacePoints"
+ftdhitname = "FTDTrackerHits"
+from Configurables import PlanarDigiAlg
+digiVXD = PlanarDigiAlg("VXDDigi")
+digiVXD.SimTrackHitCollection = "VXDCollection"
+digiVXD.TrackerHitCollection = vxdhitname
+digiVXD.ResolutionU = [0.0028, 0.006, 0.004, 0.004, 0.004, 0.004]
+digiVXD.ResolutionV = [0.0028, 0.006, 0.004, 0.004, 0.004, 0.004]
+
+digiSIT = PlanarDigiAlg("SITDigi")
+digiSIT.IsStrip = 1
+digiSIT.SimTrackHitCollection = "SITCollection"
+digiSIT.TrackerHitCollection = sithitname
+digiSIT.TrackerHitAssociationCollection = "SITTrackerHitAssociation"
+digiSIT.ResolutionU = [0.007]
+digiSIT.ResolutionV = [0.000]
+
+digiSET = PlanarDigiAlg("SETDigi")
+digiSET.IsStrip = 1
+digiSET.SimTrackHitCollection = "SETCollection"
+digiSET.TrackerHitCollection = sethitname
+digiSET.TrackerHitAssociationCollection = "SETTrackerHitAssociation"
+digiSET.ResolutionU = [0.007]
+digiSET.ResolutionV = [0.000]
+
+digiFTD = PlanarDigiAlg("FTDDigi")
+digiFTD.SimTrackHitCollection = "FTDCollection"
+digiFTD.TrackerHitCollection = ftdhitname
+digiFTD.TrackerHitAssociationCollection = "FTDTrackerHitAssociation"
+digiFTD.ResolutionU = [0.003, 0.003, 0.007, 0.007, 0.007, 0.007, 0.007, 0.007]
+digiFTD.ResolutionV = [0.003, 0.003, 0, 0, 0, 0, 0, 0 ]
+#digiFTD.OutputLevel = DEBUG
+
+from Configurables import SpacePointBuilderAlg
+spSIT = SpacePointBuilderAlg("SITBuilder")
+spSIT.TrackerHitCollection = sithitname
+spSIT.TrackerHitAssociationCollection = "SITTrackerHitAssociation"
+spSIT.SpacePointCollection = sitspname
+spSIT.SpacePointAssociationCollection = "SITSpacePointAssociation"
+#spSIT.OutputLevel = DEBUG
+
+spFTD = SpacePointBuilderAlg("FTDBuilder")
+spFTD.TrackerHitCollection = ftdhitname
+spFTD.TrackerHitAssociationCollection = "FTDTrackerHitAssociation"
+spFTD.SpacePointCollection = ftdspname
+spFTD.SpacePointAssociationCollection = "FTDSpacePointAssociation"
+#spFTD.OutputLevel = DEBUG
+
+
+##############################################################################
+# Silicon reconstruction
+##############################################################################
+from Configurables import TrackSystemSvc
+tracksystemsvc = TrackSystemSvc("TrackSystemSvc")
+
+from Configurables import SiliconTrackingAlg
+tracking = SiliconTrackingAlg("SiliconTracking")
+tracking.HeaderCol = "EventHeader"
+tracking.VTXHitCollection = vxdhitname
+tracking.SITHitCollection = sitspname
+tracking.FTDPixelHitCollection = ftdhitname
+tracking.FTDSpacePointCollection = ftdspname
+tracking.SITRawHitCollection = sithitname
+tracking.FTDRawHitCollection = ftdhitname
+tracking.UseSIT = 1
+tracking.SmoothOn = 0
+#tracking.OutputLevel = DEBUG
+
+from Configurables import ForwardTrackingAlg
+forward = ForwardTrackingAlg("ForwardTracking")
+forward.FTDPixelHitCollection = ftdhitname
+forward.FTDSpacePointCollection = ftdspname
+forward.FTDRawHitCollection = ftdhitname
+forward.Chi2ProbCut = 0.0
+forward.HitsPerTrackMin = 3
+forward.BestSubsetFinder = "SubsetSimple"
+forward.Criteria = ["Crit2_DeltaPhi","Crit2_StraightTrackRatio","Crit3_3DAngle","Crit3_ChangeRZRatio","Crit3_IPCircleDist","Crit4_3DAngleChange","Crit4_DistToExtrapolation",
+ "Crit2_DeltaRho","Crit2_RZRatio","Crit3_PT"]
+forward.CriteriaMin = [0, 0.9, 0, 0.995, 0, 0.8, 0, 20, 1.002, 0.1, 0, 0.99, 0, 0.999, 0, 0.99, 0]
+forward.CriteriaMax = [30, 1.02, 10, 1.015, 20, 1.3, 1.0, 150, 1.08, 99999999, 0.8, 1.01, 0.35, 1.001, 1.5, 1.01, 0.05]
+#forward.OutputLevel = DEBUG
+
+from Configurables import TrackSubsetAlg
+subset = TrackSubsetAlg("TrackSubset")
+subset.TrackInputCollections = ["ForwardTracks", "SiTracks"]
+subset.RawTrackerHitCollections = [vxdhitname, sithitname, ftdhitname, sitspname, ftdspname]
+subset.TrackSubsetCollection = "SubsetTracks"
+#subset.OutputLevel = DEBUG
+
+##############################################################################
+# TruthTrackerAlg
+##############################################################################
+from Configurables import TruthTrackerAlg
+truthTrackerAlg = TruthTrackerAlg("TruthTrackerAlg")
+truthTrackerAlg.SiSubsetTrackCollection = "SubsetTracks"
+
+##############################################################################
+# RecGenfitAlgSDT
+##############################################################################
+#from Configurables import RecGenfitAlgSDT
+#recGenfitAlgSDT = RecGenfitAlgSDT("RecGenfitAlgSDT")
+#recGenfitAlgSDT.debug=10
+
+##############################################################################
+# NTupleSvc
+##############################################################################
+from Configurables import NTupleSvc
+ntsvc = NTupleSvc("NTupleSvc")
+ntsvc.Output = [
+#"MyTuples DATAFILE='DCH_digi_ana.root' OPT='NEW' TYP='ROOT'",
+ "RecGenfitAlgSDT DATAFILE='fit_SDT.root' OPT='NEW' TYP='ROOT'"]
+
+
+##############################################################################
+# ApplicationMgr
+##############################################################################
+from Configurables import ApplicationMgr
+ApplicationMgr(
+ TopAlg = [genalg, detsimalg, digiVXD, digiSIT, digiSET, digiFTD, spSIT,
+ spFTD, tracking, forward, subset, dCHDigiAlg, truthTrackerAlg,
+ #recGenfitAlgSDT,
+ out],
+ EvtSel = 'NONE',
+ EvtMax = 1,
+ ExtSvc = [rndmengine, dsvc, ntsvc, evtseeder, geosvc, gearsvc, tracksystemsvc],
+ HistogramPersistency = "ROOT",
+ OutputLevel=DEBUG
+)
--
GitLab