diff --git a/Detector/DetCRD/compact/CRD_common_v01/VXD_StaggeredLadder_v01_01.xml b/Detector/DetCRD/compact/CRD_common_v01/VXD_StaggeredLadder_v01_01.xml
index dde9f4808f438affd76e9fb752438b987492e6e5..2c3c63d33dfebfe092b71deddd759918a9a39228 100644
--- a/Detector/DetCRD/compact/CRD_common_v01/VXD_StaggeredLadder_v01_01.xml
+++ b/Detector/DetCRD/compact/CRD_common_v01/VXD_StaggeredLadder_v01_01.xml
@@ -19,9 +19,18 @@
<detectors>
<detector id="DetID_VXD" name="VXD" type="SiTrackerStaggeredLadder_v01" vis="VXDVis" readout="VXDCollection" insideTrackingVolume="true">
- <envelope>
- <shape type="Tube" rmin="VXD_inner_radius" rmax="VXD_outer_radius" dz="VXD_half_length" material="Air"/>
- </envelope>
+ <envelope>
+ <shape type="BooleanShape" operation="Subtraction" material="Air" >
+ <shape type="BooleanShape" operation="Subtraction" material="Air" >
+ <shape type="Tube" rmin="VXD_inner_radius+1*mm" rmax="VXD_outer_radius" dz="VXD_half_length" />
+ <shape type="Cone" rmin1="0" rmax1="BeamPipe_VertexRegion_rmax" rmin2="0" rmax2="Vertex_Side_rmin" z="(VXD_half_length-BeamPipe_CentralAl_zmax)/2." />
+ <position x="0" y="0" z="VXD_half_length-(VXD_half_length-BeamPipe_CentralAl_zmax)/2."/>
+ </shape>
+ <shape type="Cone" rmin1="0" rmax1="BeamPipe_VertexRegion_rmax" rmin2="0" rmax2="Vertex_Side_rmin" z="(VXD_half_length-BeamPipe_CentralAl_zmax)/2." />
+ <position x="0" y="0" z="-(VXD_half_length-(VXD_half_length-BeamPipe_CentralAl_zmax)/2.)"/>
+ <rotation x="0" y="180.*deg" z="0" />
+ </shape>
+ </envelope>
<type_flags type="DetType_TRACKER + DetType_BARREL + DetType_PIXEL "/>
@@ -92,7 +101,7 @@
<readouts>
<readout name="VXDCollection">
- <id>system:5,side:-2,layer:9,module:8,sensor:8,barrelside:-2</id>
+ <id>system:5,side:-2,layer:9,module:8,active:8,sensor:8</id>
</readout>
</readouts>
</lccdd>
diff --git a/Detector/DetCRD/src/Tracker/SiTrackerStaggeredLadder_v01_geo.cpp b/Detector/DetCRD/src/Tracker/SiTrackerStaggeredLadder_v01_geo.cpp
index 120bbab25b83cce80d98b0345faaa566220910c9..1496aff196289e0d5700f48255054a38bcc5f3f9 100644
--- a/Detector/DetCRD/src/Tracker/SiTrackerStaggeredLadder_v01_geo.cpp
+++ b/Detector/DetCRD/src/Tracker/SiTrackerStaggeredLadder_v01_geo.cpp
@@ -274,11 +274,11 @@ static dd4hep::Ref_t create_element(dd4hep::Detector& theDetector, xml_h e, dd4h
double zpos = -sensor_total_z/2.0 + sensor_active_len/2.0 + isensor*(sensor_active_len + dead_gap);
pv = SensorTopEnvelopeLogical.placeVolume(SensorLogical, Position(xpos,ypos_active,zpos));
//pv.addPhysVolID("topsensor", isensor ) ;
- pv.addPhysVolID("sensor", isensor ).addPhysVolID("barrelside", 1) ;
+ pv.addPhysVolID("layer", layer_id*2+1).addPhysVolID("active", 0).addPhysVolID("sensor", isensor) ;
TopSensor_pv.push_back(pv);
pv = SensorBottomEnvelopeLogical.placeVolume(SensorLogical, Position(xpos,ypos_active,zpos));
//pv.addPhysVolID("bottomsensor", isensor ) ;
- pv.addPhysVolID("sensor", isensor ).addPhysVolID("barrelside", -1) ;
+ pv.addPhysVolID("layer", layer_id*2 ).addPhysVolID("active", 0).addPhysVolID("sensor", isensor) ;
BottomSensor_pv.push_back(pv);
pv = SensorTopEnvelopeLogical.placeVolume(SensorDeadLogical, Position(xpos,ypos_dead,zpos));
pv = SensorBottomEnvelopeLogical.placeVolume(SensorDeadLogical, Position(xpos,ypos_dead,zpos));
@@ -348,10 +348,10 @@ static dd4hep::Ref_t create_element(dd4hep::Detector& theDetector, xml_h e, dd4h
}
Transform3D tr (RotationZYX(ladder_dphi*i,0.,0.),Position(ladder_radius*cos(ladder_phi0+ladder_dphi*i), ladder_radius*sin(ladder_phi0+ladder_dphi*i), 0.));
pv = layer_assembly.placeVolume(LadderLogical,tr);
- pv.addPhysVolID("layer", layer_id ).addPhysVolID("module", i ) ;
+ pv.addPhysVolID("module", i ) ;
ladderDE.setPlacement(pv);
std::cout << ladder_enum.str() << " done." << endl;
-
+ if(i==0) std::cout << "xy=" << ladder_radius*cos(ladder_phi0) << " " << ladder_radius*sin(ladder_phi0) << std::endl;
}
// package the reconstruction data
@@ -364,12 +364,12 @@ static dd4hep::Ref_t create_element(dd4hep::Detector& theDetector, xml_h e, dd4h
topLayer.lengthSensor = sensor_active_len;
topLayer.distanceSupport = sensitive_radius;
topLayer.thicknessSupport = support_thickness / 2.0;
- topLayer.offsetSupport = ladder_offset;
+ topLayer.offsetSupport = -ladder_offset;
topLayer.widthSupport = support_width;
topLayer.zHalfSupport = support_length / 2.0;
- topLayer.distanceSensitive = sensitive_radius + support_thickness / 2.0;
+ topLayer.distanceSensitive = sensitive_radius + support_height / 2.0 + flex_thickness;
topLayer.thicknessSensitive = sensor_thickness;
- topLayer.offsetSensitive = ladder_offset;
+ topLayer.offsetSensitive = -ladder_offset + (support_width/2.0 - sensor_active_width/2.0);
topLayer.widthSensitive = sensor_active_width;
topLayer.zHalfSensitive = (n_sensors_per_side*(sensor_active_len + dead_gap) - dead_gap) / 2.0;
@@ -377,19 +377,19 @@ static dd4hep::Ref_t create_element(dd4hep::Detector& theDetector, xml_h e, dd4h
bottomLayer.phi0 = 0.;
bottomLayer.sensorsPerLadder = n_sensors_per_side;
bottomLayer.lengthSensor = sensor_active_len;
- bottomLayer.distanceSupport = sensitive_radius - support_thickness / 2.0;
+ bottomLayer.distanceSupport = sensitive_radius - support_height / 2.0 - flex_thickness;
bottomLayer.thicknessSupport = support_thickness / 2.0;
- bottomLayer.offsetSupport = ladder_offset;
+ bottomLayer.offsetSupport = -ladder_offset;
bottomLayer.widthSupport = support_width;
bottomLayer.zHalfSupport = support_length / 2.0;
- bottomLayer.distanceSensitive = sensitive_radius - support_thickness / 2.0 - sensor_thickness;
+ bottomLayer.distanceSensitive = sensitive_radius - support_height / 2.0 - sensor_thickness - flex_thickness;
bottomLayer.thicknessSensitive = sensor_thickness;
- bottomLayer.offsetSensitive = ladder_offset;
+ bottomLayer.offsetSensitive = -ladder_offset + (support_width/2.0 - sensor_active_width/2.0);
bottomLayer.widthSensitive = sensor_active_width;
bottomLayer.zHalfSensitive = (n_sensors_per_side*(sensor_active_len + dead_gap) - dead_gap) / 2.0;
- zPlanarData->layers.push_back(topLayer);
zPlanarData->layers.push_back(bottomLayer);
+ zPlanarData->layers.push_back(topLayer);
}
std::cout << (*zPlanarData) << endl;
vxd.addExtension< ZPlanarData >(zPlanarData);
diff --git a/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.cpp b/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.cpp
index 7e4dda5ce5b68c3bd6f9b03605b6b980e98b1654..df95acc6cbe9cfa3fe53902b659792f367b5c558 100644
--- a/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.cpp
+++ b/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.cpp
@@ -32,6 +32,8 @@
//#include "Tools/KiTrackMarlinCEDTools.h"
#include "Tools/FTDHelixFitter.h"
+#include <TStopwatch.h>
+
using namespace MarlinTrk ;
// Used to fedine the quality of the track output collection
@@ -64,6 +66,23 @@ StatusCode ForwardTrackingAlg::initialize(){
// can be printed. As this is mainly used for debugging it is not a steerable parameter.
//if( _useCED )MarlinCED::init(this) ; //CED
+ if(m_dumpTime){
+ NTuplePtr nt1(ntupleSvc(), "MyTuples/Time"+name());
+ if ( !nt1 ) {
+ m_tuple = ntupleSvc()->book("MyTuples/Time"+name(),CLID_ColumnWiseTuple,"Tracking time");
+ if ( 0 != m_tuple ) {
+ m_tuple->addItem ("timeTotal", m_timeTotal ).ignore();
+ }
+ else {
+ fatal() << "Cannot book MyTuples/Time"+name() <<endmsg;
+ return StatusCode::FAILURE;
+ }
+ }
+ else{
+ m_tuple = nt1;
+ }
+ }
+
// Now set min and max values for all the criteria
for( unsigned i=0; i < _criteriaNames.size(); i++ ){
std::vector< float > emptyVec;
@@ -181,6 +200,8 @@ StatusCode ForwardTrackingAlg::initialize(){
StatusCode ForwardTrackingAlg::execute(){
debug() << " processing event number " << _nEvt << endmsg;
+ auto stopwatch = TStopwatch();
+
auto trkCol = _outColHdl.createAndPut();
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// //
@@ -683,6 +704,12 @@ StatusCode ForwardTrackingAlg::execute(){
//if( _useCED ) MarlinCED::draw(this);
_nEvt ++ ;
+
+ if(m_dumpTime&&m_tuple){
+ m_timeTotal = stopwatch.RealTime()*1000;
+ m_tuple->write();
+ }
+
return StatusCode::SUCCESS;
}
diff --git a/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.h b/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.h
index 56de6af9f2bb3ea94b24554b3139a10b658b0600..bd7a5173be02214e734579e3d4c025f8006ec1a2 100644
--- a/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.h
+++ b/Reconstruction/SiliconTracking/src/ForwardTrackingAlg.h
@@ -21,6 +21,8 @@
#include "Criteria/Criteria.h"
#include "ILDImpl/SectorSystemFTD.h"
+#include "GaudiKernel/NTuple.h"
+
using namespace KiTrack;
using namespace KiTrackMarlin;
@@ -234,7 +236,10 @@ class ForwardTrackingAlg : public GaudiAlgorithm {
int _nRun ;
int _nEvt ;
-
+
+ NTuple::Tuple* m_tuple;
+ NTuple::Item<float> m_timeTotal;
+
/** B field in z direction */
double _Bz;
@@ -255,6 +260,7 @@ class ForwardTrackingAlg : public GaudiAlgorithm {
Gaudi::Property<std::vector<std::string> > _criteriaNames{this, "Criteria", Criteria::getAllCriteriaNamesVec()};
Gaudi::Property<std::vector<float> > _critMinimaInit{this, "CriteriaMin", {} };
Gaudi::Property<std::vector<float> > _critMaximaInit{this, "CriteriaMax", {} };
+ Gaudi::Property<bool> m_dumpTime{this, "DumpTime", true};
std::map<std::string, std::vector<float> > _critMinima;
std::map<std::string, std::vector<float> > _critMaxima;
diff --git a/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp b/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp
index 7fbb98e9fe962a79542bbbd31b789f3a057f3c32..f9b6a5466aba7cf08527d77d2f3b6c5d429ba3d8 100644
--- a/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp
+++ b/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.cpp
@@ -34,6 +34,7 @@
#include "TrackSystemSvc/HelixFit.h"
#include "TrackSystemSvc/IMarlinTrack.h"
+#include <TStopwatch.h>
//#include "TrackSystemSvc/MarlinTrkDiagnostics.h"
//#ifdef MARLINTRK_DIAGNOSTICS_ON
//#include "TrackSystemSvc/DiagnosticsController.h"
@@ -104,7 +105,23 @@ StatusCode SiliconTrackingAlg::initialize() {
_nRun = -1 ;
_nEvt = 0 ;
//printParameters() ;
-
+ if(m_dumpTime){
+ NTuplePtr nt1(ntupleSvc(), "MyTuples/Time"+name());
+ if ( !nt1 ) {
+ m_tuple = ntupleSvc()->book("MyTuples/Time"+name(),CLID_ColumnWiseTuple,"Tracking time");
+ if ( 0 != m_tuple ) {
+ m_tuple->addItem ("timeTotal", m_timeTotal ).ignore();
+ m_tuple->addItem ("timeKalman", m_timeKalman ).ignore();
+ }
+ else {
+ fatal() << "Cannot book MyTuples/Time"+name() <<endmsg;
+ return StatusCode::FAILURE;
+ }
+ }
+ else{
+ m_tuple = nt1;
+ }
+ }
// set up the geometery needed by KalTest
//FIXME: for now do KalTest only - make this a steering parameter to use other fitters
auto _trackSystemSvc = service<ITrackSystemSvc>("TrackSystemSvc");
@@ -157,6 +174,7 @@ StatusCode SiliconTrackingAlg::initialize() {
}
StatusCode SiliconTrackingAlg::execute(){
+ auto stopwatch = TStopwatch();
Navigation::Instance()->Initialize();
//_current_event = evt;
//_allHits.reserve(1000);
@@ -355,6 +373,10 @@ StatusCode SiliconTrackingAlg::execute(){
CleanUp();
debug() << "Event is done " << endmsg;
_nEvt++;
+ if(m_dumpTime&&m_tuple){
+ m_timeTotal = stopwatch.RealTime()*1000;
+ m_tuple->write();
+ }
return StatusCode::SUCCESS;
}
@@ -907,7 +929,7 @@ StatusCode SiliconTrackingAlg::finalize(){
//delete _trksystem ; _trksystem = 0;
//delete _histos ; _histos = 0;
info() << "Processed " << _nEvt << " events " << endmsg;
- info() << lcio::ILDCellID0::encoder_string << " " << UTIL::ILDCellID0::encoder_string << endmsg;
+
return GaudiAlgorithm::finalize();
}
@@ -2545,6 +2567,7 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) {
float pyTot = 0.;
float pzTot = 0.;
debug() << "Total " << nTracks << " candidate tracks will be dealed" << endmsg;
+ auto stopwatch = TStopwatch();
for (int iTrk=0;iTrk<nTracks;++iTrk) {
TrackExtended * trackAR = _trackImplVec[iTrk];
@@ -2886,7 +2909,8 @@ void SiliconTrackingAlg::FinalRefit(edm4hep::TrackCollection* trk_col) {
}
}
}
-
+ if(m_dumpTime&&m_tuple) m_timeKalman = stopwatch.RealTime()*1000;
+
debug() << " -> run " << _nRun << " event " << _nEvt << endmsg;
debug() << "Number of Si tracks = " << nSiSegments << endmsg;
debug() << "Total 4-momentum of Si Tracks : E = " << std::setprecision(7) << eTot
diff --git a/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.h b/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.h
index 7f584b264d3fe274bd7a793f3ac6ea0081249e28..e8e915b56530c595f1ed01040c2d5d6421af4bfb 100644
--- a/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.h
+++ b/Reconstruction/SiliconTracking/src/SiliconTrackingAlg.h
@@ -27,6 +27,7 @@
#include <UTIL/BitField64.h>
#include <UTIL/ILDConf.h>
+#include "GaudiKernel/NTuple.h"
//using namespace edm4hep ;
namespace gear{
@@ -252,7 +253,7 @@ class SiliconTrackingAlg : public GaudiAlgorithm {
Gaudi::Property<bool> _ElossOn{this, "EnergyLossOn", true};
Gaudi::Property<bool> _SmoothOn{this, "SmoothOn", true};
Gaudi::Property<float> _helix_max_r{this, "HelixMaxR", 2000.};
-
+ Gaudi::Property<bool> m_dumpTime{this, "DumpTime", true};
//std::vector<int> _colours;
/** helper function to get collection using try catch block */
@@ -300,7 +301,11 @@ class SiliconTrackingAlg : public GaudiAlgorithm {
std::vector<TrackerHitExtendedVec> _sectors;
std::vector<TrackerHitExtendedVec> _sectorsFTD;
-
+
+ NTuple::Tuple* m_tuple;
+ NTuple::Item<float> m_timeTotal;
+ NTuple::Item<float> m_timeKalman;
+
/**
* A helper class to allow good code readability by accessing tracks with N hits.
* As the smalest valid track contains three hits, but the first index in a vector is 0,
diff --git a/Reconstruction/SiliconTracking/src/TrackSubsetAlg.cpp b/Reconstruction/SiliconTracking/src/TrackSubsetAlg.cpp
index 792685092237ca9248c16bc63dd7fca8a8ce9f1f..d4741f2eb165767d51fdaf7afee0612dacf942b7 100644
--- a/Reconstruction/SiliconTracking/src/TrackSubsetAlg.cpp
+++ b/Reconstruction/SiliconTracking/src/TrackSubsetAlg.cpp
@@ -15,6 +15,8 @@
#include "TrackSystemSvc/MarlinTrkUtils.h"
+#include <TStopwatch.h>
+
using namespace KiTrack;
DECLARE_COMPONENT(TrackSubsetAlg)
@@ -39,6 +41,23 @@ StatusCode TrackSubsetAlg::initialize() {
_nRun = 0 ;
_nEvt = 0 ;
+ if(m_dumpTime){
+ NTuplePtr nt1(ntupleSvc(), "MyTuples/Time"+name());
+ if ( !nt1 ) {
+ m_tuple = ntupleSvc()->book("MyTuples/Time"+name(),CLID_ColumnWiseTuple,"Tracking time");
+ if ( 0 != m_tuple ) {
+ m_tuple->addItem ("timeTotal", m_timeTotal ).ignore();
+ }
+ else {
+ fatal() << "Cannot book MyTuples/Time"+name() <<endmsg;
+ return StatusCode::FAILURE;
+ }
+ }
+ else{
+ m_tuple = nt1;
+ }
+ }
+
for(unsigned i=0; i<_trackInputColNames.size(); i++){
_inTrackColHdls.push_back(new DataHandle<edm4hep::TrackCollection> (_trackInputColNames[i], Gaudi::DataHandle::Reader, this));
}
@@ -96,6 +115,8 @@ StatusCode TrackSubsetAlg::finalize(){
}
StatusCode TrackSubsetAlg::execute(){
+ auto stopwatch = TStopwatch();
+
std::vector<edm4hep::Track> tracks;
auto trkCol = _outColHdl.createAndPut();
@@ -331,6 +352,12 @@ StatusCode TrackSubsetAlg::execute(){
Navigation::Instance()->Initialize();
_nEvt ++ ;
+
+ if(m_dumpTime&&m_tuple){
+ m_timeTotal = stopwatch.RealTime()*1000;
+ m_tuple->write();
+ }
+
return StatusCode::SUCCESS;
}
diff --git a/Reconstruction/SiliconTracking/src/TrackSubsetAlg.h b/Reconstruction/SiliconTracking/src/TrackSubsetAlg.h
index 952d08fa2607dee095444f3dcf0563901c28dc63..1023eba976e861461174d93a7158294328ef4854 100644
--- a/Reconstruction/SiliconTracking/src/TrackSubsetAlg.h
+++ b/Reconstruction/SiliconTracking/src/TrackSubsetAlg.h
@@ -11,6 +11,7 @@
#include "Math/ProbFunc.h"
+#include "GaudiKernel/NTuple.h"
/** Processor that takes tracks from multiple sources and outputs them (or modified versions, or a subset of them)
* as one track collection.
*
@@ -74,11 +75,15 @@ class TrackSubsetAlg : public GaudiAlgorithm {
Gaudi::Property<float> _initialTrackError_tanL{this, "InitialTrackErrorTanL",1e2};
Gaudi::Property<double> _maxChi2PerHit{this, "MaxChi2PerHit", 1e2};
Gaudi::Property<double> _omega{this, "Omega", 0.75};
+ Gaudi::Property<bool> m_dumpTime{this, "DumpTime", true};
float _bField;
int _nRun ;
int _nEvt ;
+
+ NTuple::Tuple* m_tuple;
+ NTuple::Item<float> m_timeTotal;
};
/** A functor to return whether two tracks are compatible: The criterion is if the share a TrackerHit or more */
diff --git a/Reconstruction/Tracking/src/Clupatra/clupatra_new.cpp b/Reconstruction/Tracking/src/Clupatra/clupatra_new.cpp
index c4ed1fed79fad18e61d46ae29e7afcb338f79abe..9ecee56418ae4387e73572d73b42b757a92bd2b2 100644
--- a/Reconstruction/Tracking/src/Clupatra/clupatra_new.cpp
+++ b/Reconstruction/Tracking/src/Clupatra/clupatra_new.cpp
@@ -1238,31 +1238,35 @@ start:
bool reverse_order = ( std::abs( hf->first->pos.z() ) > std::abs( hb->first->pos.z()) + 3. ) ;
unsigned nHit = 0 ;
-
+ int code = 0;
if( reverse_order ){
//std::cout << "It is true order" << std::endl;
for( CluTrack::reverse_iterator it=clu->rbegin() ; it != clu->rend() ; ++it){
edm4hep::TrackerHit ph = (*it)->first->edm4hepHit;
trk->addHit(ph) ;
++nHit ;
- //std::cout << " hit added " << (*it)->first->edm4hepHit << std::endl ;
+ //std::cout << " hit added " << (*it)->first->edm4hepHit.id() << std::endl ;
}
- trk->initialise( MarlinTrk::IMarlinTrack::forward ) ;
+ code = trk->initialise( MarlinTrk::IMarlinTrack::forward ) ;
} else {
//std::cout << "It is reverse order" << std::endl;
for( CluTrack::iterator it=clu->begin() ; it != clu->end() ; ++it){
edm4hep::TrackerHit ph = (*it)->first->edm4hepHit;
- trk->addHit(ph) ;
+ if( trk->addHit(ph) == MarlinTrk::IMarlinTrack::success ){
+ //std::cout << " hit added " << (*it)->first->edm4hepHit.id() << std::endl;
+ }
+ else{
+ //std::cout << " hit not added " << (*it)->first->edm4hepHit.id() << std::endl;
+ }
++nHit ;
- //std::cout << " hit added "<< (*it)->first->edm4hepHit << std::endl ;
}
- trk->initialise( MarlinTrk::IMarlinTrack::backward ) ;
+ code = trk->initialise( MarlinTrk::IMarlinTrack::backward ) ;
}
- int code = trk->fit( maxChi2 ) ;
+ if( code != MarlinTrk::IMarlinTrack::error ) code = trk->fit( maxChi2 ) ;
// for (auto hit : hitsInFit) std::cout << hit.first << std::endl;
if( code != MarlinTrk::IMarlinTrack::success ){
@@ -1270,7 +1274,8 @@ start:
std::cout << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> IMarlinTrkFitter : problem fitting track "
<< " error code : " << MarlinTrk::errorCode( code )
<< std::endl ;
-
+ delete trk;
+ return 0;
}
diff --git a/Reconstruction/Tracking/src/Clupatra/clupatra_new.h b/Reconstruction/Tracking/src/Clupatra/clupatra_new.h
index 46247d0b48c36b6975da0ba3411f27f9363e99d8..6ace630a7dcf6bdc552426f519f62c51a992da1a 100644
--- a/Reconstruction/Tracking/src/Clupatra/clupatra_new.h
+++ b/Reconstruction/Tracking/src/Clupatra/clupatra_new.h
@@ -467,11 +467,12 @@ namespace clupatra_new{
// FIXME Mingrui debug
// streamlog_out( DEBUG3 ) << " -- extrapolate TrackState : " << lcshort( ts ) << std::endl ;
-
- edm4hep::TrackerHit ht = trk.getTrackerHits(0);
- //need to add a dummy hit to the track
- mTrk->addHit( ht ) ; // is this the right hit ??????????
-
+ // fucd: getTrackerHits(0) is possible to miss ILDVTrackHit
+ for(int ih=0;ih<nHit;ih++){
+ edm4hep::TrackerHit ht = trk.getTrackerHits(ih);
+ //need to add a dummy hit to the track
+ if(mTrk->addHit( ht ) == MarlinTrk::IMarlinTrack::success) break; // is this the right hit ??????????
+ }
mTrk->initialise( ts , _b , MarlinTrk::IMarlinTrack::backward ) ;
double deltaChi ;
diff --git a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp
index b00f1d051256fc450499a21e30d90ed74e982243..08e1fb022fdfc240dd127957085fb954ca2edf54 100755
--- a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp
+++ b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp
@@ -53,6 +53,8 @@
#include <vector>
#include <bitset>
+#include <TStopwatch.h>
+
typedef std::vector<edm4hep::TrackerHit> TrackerHitVec;
using namespace edm4hep ;
@@ -143,7 +145,24 @@ StatusCode FullLDCTrackingAlg::initialize() {
PI = acos(-1.);
PIOVER2 = 0.5*PI;
TWOPI = 2*PI;
-
+
+ if(m_dumpTime){
+ NTuplePtr nt1(ntupleSvc(), "MyTuples/Time"+name());
+ if ( !nt1 ) {
+ m_tuple = ntupleSvc()->book("MyTuples/Time"+name(),CLID_ColumnWiseTuple,"Tracking time");
+ if ( 0 != m_tuple ) {
+ m_tuple->addItem ("timeTotal", m_timeTotal ).ignore();
+ m_tuple->addItem ("timeKalman", m_timeKalman ).ignore();
+ }
+ else {
+ fatal() << "Cannot book MyTuples/Time"+name() <<endmsg;
+ return StatusCode::FAILURE;
+ }
+ }
+ else{
+ m_tuple = nt1;
+ }
+ }
// set up the geometery needed by KalTest
//FIXME: for now do KalTest only - make this a steering parameter to use other fitters
auto _trackSystemSvc = service<ITrackSystemSvc>("TrackSystemSvc");
@@ -181,6 +200,8 @@ StatusCode FullLDCTrackingAlg::execute() {
// debug() << endmsg;
info() << "FullLDCTrackingAlg -> run = " << 0/*evt->getRunNumber()*/ << " event = " << _nEvt << endmsg;
+
+ auto stopwatch = TStopwatch();
// debug() << endmsg;
auto outCol = _OutputTrackColHdl.createAndPut();
@@ -222,6 +243,11 @@ StatusCode FullLDCTrackingAlg::execute() {
debug() << "Cleanup is done." << endmsg;
_nEvt++;
// getchar();
+ if(m_tuple){
+ m_timeTotal = stopwatch.RealTime()*1000;
+ m_tuple->write();
+ }
+
debug() << "FullLDCTrackingAlg execute() finished" << endmsg;
return StatusCode::SUCCESS;
@@ -248,6 +274,7 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr
//SJA:FIXME: So here we are going to do one final refit. This can certainly be optimised, but rather than worry about the mememory management right now lets make it work, and optimise it later ...
debug() << "Total " << nTrkCand << " trkCand to deal" << endmsg;
+ auto stopwatch = TStopwatch();
for (int iTRK=0;iTRK<nTrkCand;++iTRK) {
TrackExtended * trkCand = trkVec[iTRK];
@@ -385,9 +412,7 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr
int error_code = 0;
try {
-
error_code = MarlinTrk::createFinalisedLCIOTrack(marlinTrk, trkHits, &track, fit_backwards, &ts_initial, _bField, _maxChi2PerHit);
-
} catch (...) {
// delete track;
@@ -554,7 +579,7 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr
debug() << " Add Track to final Collection: ID = " << track.id() << " for trkCand "<< trkCand << endmsg;
}
}
-
+ if(m_tuple) m_timeKalman = stopwatch.RealTime()*1000;
// streamlog_out(DEBUG5) << endmsg;
debug() << "Number of accepted " << _OutputTrackColHdl.fullKey() << " = " << nTotTracks << endmsg;
debug() << "Total 4-momentum of " << _OutputTrackColHdl.fullKey() << " : E = " << eTot
@@ -1118,7 +1143,6 @@ void FullLDCTrackingAlg::prepareVectors() {
auto Cov = getCovMatrix(tpcTrack);
int NC = int(Cov.size());
for (int ic=0;ic<NC;ic++) {
- if (ic>=(sizeof(cov)/sizeof(float))) break;
cov[ic] = Cov[ic];
}
@@ -1183,7 +1207,6 @@ void FullLDCTrackingAlg::prepareVectors() {
auto Cov = getCovMatrix(siTrack);
int NC = int(Cov.size());
for (int ic=0;ic<NC;ic++) {
- if (ic>=(sizeof(cov)/sizeof(float))) break;
cov[ic] = Cov[ic];
}
trackExt->setCovMatrix(cov);
@@ -3413,7 +3436,7 @@ void FullLDCTrackingAlg::AddNotAssignedHits() {
nonAssignedTPCHits.push_back(trkHitExt);
}
}
- debug() << "AddNotAssignedHits : Number of Non Assigned TPC hits = " << nonAssignedTPCHits.size() << endmsg;
+ debug() << "AddNotAssignedHits : Number of Non Assigned TPC hits = " << nonAssignedTPCHits.size() << " distance cut = " << _distCutForTPCHits << endmsg;
AssignTPCHitsToTracks(nonAssignedTPCHits, _distCutForTPCHits);
}
@@ -3856,6 +3879,7 @@ void FullLDCTrackingAlg::AssignTPCHitsToTracks(TrackerHitExtendedVec hitVec,
debug() << "AssignTPCHitsToTracks: Starting loop " << nTrk << " tracks and " << nHits << " hits" << endmsg;
for (int iT=0;iT<nTrk;++iT) { // loop over all tracks
+ debug() << " track " << iT << endmsg;
TrackExtended * foundTrack = _trkImplVec[iT];
int tanlambdaSign = std::signbit(foundTrack->getTanLambda());//we only care about positive or negative
GroupTracks * group = foundTrack->getGroupTracks();
@@ -3874,7 +3898,7 @@ void FullLDCTrackingAlg::AssignTPCHitsToTracks(TrackerHitExtendedVec hitVec,
HelixClass helix;
helix.Initialize_Canonical(phi0,d0,z0,omega,tanLambda,_bField);
float OnePFivehalfPeriodZ = 1.5*fabs(acos(-1.)*tanLambda/omega);
-
+ debug() << " OnePFivehalfPeriodZ = " << OnePFivehalfPeriodZ << endmsg;
for (int iH=0;iH<nHits;++iH) { // loop over leftover TPC hits
//check if the hit and the track or on the same side
@@ -3886,9 +3910,10 @@ void FullLDCTrackingAlg::AssignTPCHitsToTracks(TrackerHitExtendedVec hitVec,
bool consider = DeltaStart <= OnePFivehalfPeriodZ;
consider = consider || (DeltaEnd <= OnePFivehalfPeriodZ);
consider = consider || ( (HitPositions[iH][2]>=startPointZ) && (HitPositions[iH][2]<=endPointZ) );
-
+ debug() << " hit " << iH << " " << consider << " DeltaStart = " << DeltaStart << " DeltaEnd = " << DeltaEnd << endmsg;
if(consider){
float distance = helix.getDistanceToPoint(HitPositions[iH], minDistances[iH]);
+ debug() << " distance = " << distance << " cut = " << minDistances[iH] << endmsg;
if (distance < minDistances[iH]) {
minDistances[iH] = distance;
tracksToAttach[iH] = foundTrack;
@@ -3903,7 +3928,9 @@ void FullLDCTrackingAlg::AssignTPCHitsToTracks(TrackerHitExtendedVec hitVec,
if (tracksToAttach[iH]!=NULL) {
tracksToAttach[iH]->addTrackerHitExtended(trkHitExt);
trkHitExt->setTrackExtended( tracksToAttach[iH] );
- trkHitExt->setUsedInFit( false );
+ //by fucd
+ //trkHitExt->setUsedInFit( false );
+ trkHitExt->setUsedInFit( true );
}
}
diff --git a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h
index 6822d60f4dd05007b757256eb2e57efe13435ecc..cfc804879d44a7e453d069e79abe261640494120 100755
--- a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h
+++ b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h
@@ -27,6 +27,7 @@
#include <UTIL/BitField64.h>
#include <UTIL/ILDConf.h>
+#include "GaudiKernel/NTuple.h"
//using namespace edm4hep ;
//class gear::GearMgr ;
@@ -394,7 +395,7 @@ protected:
Gaudi::Property<float> _vetoMergeMomentumCut{this, "VetoMergeMomentumCut", 2.5};
Gaudi::Property<float> _maxAllowedPercentageOfOutliersForTrackCombination{this, "MaxAllowedPercentageOfOutliersForTrackCombination", 0.3};
Gaudi::Property<int> _maxAllowedSiHitRejectionsForTrackCombination{this, "MaxAllowedSiHitRejectionsForTrackCombination", 2};
-
+ Gaudi::Property<bool> m_dumpTime{this, "DumpTime", true};
//float _dPCutForForcedMerging;
bool _runMarlinTrkDiagnostics = false;
@@ -514,6 +515,10 @@ protected:
DataHandle<edm4hep::TrackCollection> _SiTrackColHdl{"SiTracks", Gaudi::DataHandle::Reader, this};
DataHandle<edm4hep::TrackCollection> _OutputTrackColHdl{"MarlinTrkTracks", Gaudi::DataHandle::Writer, this};
+
+ NTuple::Tuple* m_tuple;
+ NTuple::Item<float> m_timeTotal;
+ NTuple::Item<float> m_timeKalman;
};
#endif
diff --git a/Service/GearSvc/src/GearSvc.cpp b/Service/GearSvc/src/GearSvc.cpp
index 81c3611f63eafe0e0eb537b127b4cb44c0045a74..e96ae681694de49f30b430dedf1ec1ac9568f0e3 100644
--- a/Service/GearSvc/src/GearSvc.cpp
+++ b/Service/GearSvc/src/GearSvc.cpp
@@ -12,11 +12,13 @@
#include "gearimpl/FixedPadSizeDiskLayout.h"
#include "gearimpl/CalorimeterParametersImpl.h"
#include "gearimpl/SimpleMaterialImpl.h"
+#include "gear/VXDLayerLayout.h"
#include "gearxml/tinyxml.h"
#include "DD4hep/Detector.h"
#include "DD4hep/DetElement.h"
#include "DDRec/DetectorData.h"
+#include "DDRec/MaterialManager.h"
#include "DD4hep/DD4hepUnits.h"
#include "CLHEP/Units/SystemOfUnits.h"
@@ -194,6 +196,49 @@ StatusCode GearSvc::convertVXD(dd4hep::DetElement& vxd){
extensionDataValid = false;
info() << e.what() << " " << vxdData << endmsg;
}
+ if(vxdData){
+ int vxdType = gear::ZPlanarParameters::CMOS;
+ gear::ZPlanarParametersImpl* gearVXD = new gear::ZPlanarParametersImpl( vxdType, vxdData->rInnerShell/dd4hep::mm, vxdData->rOuterShell/dd4hep::mm,
+ vxdData->zHalfShell/dd4hep::mm , vxdData->gapShell/dd4hep::mm , 0. ) ;
+ for(unsigned i=0,n=vxdData->layers.size() ; i<n; ++i){
+ const dd4hep::rec::ZPlanarData::LayerLayout& l = vxdData->layers[i];
+ // FIXME set rad lengths to 0 -> need to get from dd4hep ....
+ gearVXD->addLayer(l.ladderNumber, l.phi0,
+ l.distanceSupport/dd4hep::mm, l.offsetSupport/dd4hep::mm, l.thicknessSupport/dd4hep::mm, l.zHalfSupport/dd4hep::mm, l.widthSupport/dd4hep::mm, 0.,
+ l.distanceSensitive/dd4hep::mm, l.offsetSensitive/dd4hep::mm, l.thicknessSensitive/dd4hep::mm, l.zHalfSensitive/dd4hep::mm, l.widthSensitive/dd4hep::mm, 0.);
+ }
+ m_gearMgr->setVXDParameters(gearVXD);
+
+ dd4hep::rec::MaterialManager matMgr( dd4hep::Detector::getInstance().world().volume() ) ;
+ const dd4hep::rec::ZPlanarData::LayerLayout& l = vxdData->layers[0] ;
+ dd4hep::rec::Vector3D a( l.distanceSensitive + l.thicknessSensitive, l.phi0 , 0. , dd4hep::rec::Vector3D::cylindrical ) ;
+ dd4hep::rec::Vector3D b( l.distanceSupport + l.thicknessSupport, l.phi0 , 0. , dd4hep::rec::Vector3D::cylindrical ) ;
+ const dd4hep::rec::MaterialVec& materials = matMgr.materialsBetween( a , b ) ;
+ dd4hep::rec::MaterialData mat = ( materials.size() > 1 ? matMgr.createAveragedMaterial( materials ) : materials[0].first ) ;
+
+ std::cout << " ####### found materials between points : " << a << " and " << b << " : " ;
+ for( unsigned i=0,n=materials.size();i<n;++i){
+ std::cout << materials[i].first.name() << "[" << materials[i].second << "], " ;
+ }
+ std::cout << std::endl ;
+ std::cout << " averaged material : " << mat << std::endl ;
+ gear::SimpleMaterialImpl* VXDSupportMaterial = new gear::SimpleMaterialImpl("VXDSupportMaterial", mat.A(), mat.Z(),
+ mat.density()/(dd4hep::kg/(dd4hep::g*dd4hep::m3)),
+ mat.radiationLength()/dd4hep::mm,
+ mat.interactionLength()/dd4hep::mm);
+ m_gearMgr->registerSimpleMaterial(VXDSupportMaterial);
+
+ info() << vxdData->rInnerShell << " " << vxdData->rOuterShell << " " << vxdData->zHalfShell << " " << vxdData->gapShell << endmsg;
+ for(int i=0,n=vxdData->layers.size(); i<n; i++){
+ const dd4hep::rec::ZPlanarData::LayerLayout& thisLayer = vxdData->layers[i];
+ info() << i << ": " << thisLayer.ladderNumber << "," << thisLayer.phi0 << ","
+ << thisLayer.distanceSupport/dd4hep::mm << "," << thisLayer.offsetSupport/dd4hep::mm << "," << thisLayer.thicknessSupport/dd4hep::mm << ","
+ << thisLayer.zHalfSupport/dd4hep::mm << "," << thisLayer.widthSupport/dd4hep::mm << "," << "NULL,"
+ << thisLayer.distanceSensitive/dd4hep::mm << "," << thisLayer.offsetSensitive/dd4hep::mm << "," << thisLayer.thicknessSensitive/dd4hep::mm << ","
+ << thisLayer.zHalfSensitive/dd4hep::mm << "," << thisLayer.widthSensitive/dd4hep::mm << ",NULL" << endmsg;
+ }
+ return sc;
+ }
std::vector<helpLayer> helpSensitives;
std::vector<helpLayer> helpLadders;
@@ -575,30 +620,34 @@ StatusCode GearSvc::convertVXD(dd4hep::DetElement& vxd){
gearParameters->setDoubleVal("CryostatAlHalfZ", dzSty+drSty);
m_gearMgr->setGearParameters("VXDInfra", gearParameters);
}
+
+ dd4hep::rec::MaterialManager matMgr( dd4hep::Detector::getInstance().world().volume() ) ;
+ const gear::VXDLayerLayout& l = m_gearMgr->getVXDParameters().getVXDLayerLayout();
+ dd4hep::rec::Vector3D a( l.getSensitiveDistance(0)+l.getSensitiveThickness(0), l.getPhi0(0), 0. , dd4hep::rec::Vector3D::cylindrical ) ;
+ dd4hep::rec::Vector3D b( l.getLadderDistance(0)+l.getLadderThickness(0), l.getPhi0(0), 0. , dd4hep::rec::Vector3D::cylindrical ) ;
+ const dd4hep::rec::MaterialVec& materials = matMgr.materialsBetween( a , b ) ;
+ dd4hep::rec::MaterialData mat = ( materials.size() > 1 ? matMgr.createAveragedMaterial( materials ) : materials[0].first ) ;
+
+ std::cout << " ####### found materials between points : " << a << " and " << b << " : " ;
+ for( unsigned i=0,n=materials.size();i<n;++i){
+ std::cout << materials[i].first.name() << "[" << materials[i].second << "], " ;
+ }
+ std::cout << std::endl ;
+ std::cout << " averaged material : " << mat << std::endl ;
+ gear::SimpleMaterialImpl* VXDSupportMaterial = new gear::SimpleMaterialImpl("VXDSupportMaterial", mat.A(), mat.Z(),
+ mat.density()/(dd4hep::kg/(dd4hep::g*dd4hep::m3)),
+ mat.radiationLength()/dd4hep::mm,
+ mat.interactionLength()/dd4hep::mm);
+
//effective A different with what in Mokka, fix them as Mokka's
gear::SimpleMaterialImpl* VXDFoamShellMaterial_old = new gear::SimpleMaterialImpl("VXDFoamShellMaterial_old", 1.043890843e+01, 5.612886646e+00, 2.500000000e+01, 1.751650267e+04, 0);
m_gearMgr->registerSimpleMaterial(VXDFoamShellMaterial_old);
gear::SimpleMaterialImpl* VXDFoamShellMaterial = new gear::SimpleMaterialImpl("VXDFoamShellMaterial", styAeff, styZeff, styDensity*(CLHEP::g/CLHEP::cm3)/(CLHEP::kg/CLHEP::m3),
styRadLen, styIntLen);
m_gearMgr->registerSimpleMaterial(VXDFoamShellMaterial);
- gear::SimpleMaterialImpl* VXDSupportMaterial_old = new gear::SimpleMaterialImpl("VXDSupportMaterial_old", 2.075865162e+01, 1.039383117e+01, 2.765900000e+02, 1.014262421e+03, 3.341388059e+03);
- m_gearMgr->registerSimpleMaterial(VXDSupportMaterial_old);
- gear::SimpleMaterialImpl* VXDSupportMaterial = new gear::SimpleMaterialImpl("VXDSupportMaterial", VXDSupportAeff, VXDSupportZeff, VXDSupportDensity, VXDSupportRadLen, VXDSupportIntLen);
m_gearMgr->registerSimpleMaterial(VXDSupportMaterial);
- info() << "=====================from ZPlanarData==============================" << endmsg;
- if(vxdData){
- info() << vxdData->rInnerShell << " " << vxdData->rOuterShell << " " << vxdData->zHalfShell << " " << vxdData->gapShell << endmsg;
- const std::vector<dd4hep::rec::ZPlanarData::LayerLayout>& layers = vxdData->layers;
- for(int i=0;i<layers.size();i++){
- const dd4hep::rec::ZPlanarData::LayerLayout& thisLayer = layers[i];
- info() << i << ": " << thisLayer.ladderNumber << "," << thisLayer.phi0 << "," << thisLayer.distanceSupport << "," << thisLayer.offsetSupport << ","
- << thisLayer.thicknessSupport << "," << thisLayer.zHalfSupport << "," << thisLayer.widthSupport << "," << "NULL,"
- << thisLayer.distanceSensitive << "," << thisLayer.offsetSensitive << "," << thisLayer.thicknessSensitive << "," << thisLayer.zHalfSensitive << ","
- << thisLayer.widthSensitive << ",NULL" << endmsg;
- }
- }
- info() << rAlu << " " << drAlu << " " << rInner << " " << aluEndcapZ << " " << aluHalfZ << endmsg;
- //info() << m_materials["VXDSupportMaterial"] << endmsg;
+ info() << "cryostat parameters: " << rAlu << " " << drAlu << " " << rInner << " " << aluEndcapZ << " " << aluHalfZ << endmsg;
+
return sc;
}
@@ -853,7 +902,7 @@ StatusCode GearSvc::convertDC(dd4hep::DetElement& dc){
}
info() << (*dcData) << endmsg;
}
-
+ debug() << dcData->maxRow << ": " << dcData->rMinReadout*CLHEP::cm << " " << dcData->rMaxReadout*CLHEP::cm << endmsg;
// regard as TPCParameters, TODO: drift chamber parameters
gear::TPCParametersImpl *tpcParameters = new gear::TPCParametersImpl();
gear::PadRowLayout2D *padLayout = new gear::FixedPadSizeDiskLayout(dcData->rMinReadout*CLHEP::cm, dcData->rMaxReadout*CLHEP::cm,