From f55870002f26c835c80ab4f58d3564cb4667e10a Mon Sep 17 00:00:00 2001 From: Ke Li <ke.li@cern.ch> Date: Tue, 21 Jun 2022 20:10:01 +0200 Subject: [PATCH] update the trackinggeometry and materialmapping --- .../TrackerSPFit/src/TrackerSPFit.cxx | 408 ++++--- .../TrackerSPFit/src/TrackerSPFit.h | 5 + .../TrackerSpacePoint/TrackerSeed.h | 10 +- .../TrackerSpacePoint/src/TrackerSeed.cxx | 14 +- .../Acts/FaserActsGeometry/CMakeLists.txt | 138 +-- .../FaserActsJsonGeometryConverter.h | 261 ---- .../FaserActsGeometry/FaserActsLayerBuilder.h | 14 +- .../python/ActsGeometryConfig.py | 26 +- .../python/FaserActsExtrapolationConfig.py | 15 +- .../FaserActsMaterialMapping_jobOptions.py | 2 +- .../python/FaserActsMaterialValidationAlg.py | 66 + .../FaserActsWriteTrackingGeometryConfig.py | 23 +- .../src/CuboidVolumeBuilder.cxx | 23 +- .../src/FaserActsDetectorElement.cxx | 1 - .../src/FaserActsExtrapolationAlg.cxx | 38 +- .../src/FaserActsExtrapolationAlg.h | 8 +- .../src/FaserActsExtrapolationTool.cxx | 2 +- .../src/FaserActsJsonGeometryConverter.cxx | 1071 ----------------- .../src/FaserActsLayerBuilder.cxx | 318 +++-- .../src/FaserActsMaterialJsonWriterTool.cxx | 11 +- .../src/FaserActsMaterialJsonWriterTool.h | 4 +- .../src/FaserActsMaterialMapping.cxx | 108 +- .../src/FaserActsMaterialMapping.h | 4 +- .../src/FaserActsMaterialTrackWriterSvc.cxx | 354 ++++++ .../src/FaserActsMaterialTrackWriterSvc.h | 109 ++ .../src/FaserActsPropStepRootWriterSvc.h | 4 - .../src/FaserActsSurfaceMappingTool.cxx | 2 +- .../src/FaserActsTrackingGeometrySvc.cxx | 51 +- .../src/FaserActsTrackingGeometrySvc.h | 7 +- .../src/FaserActsWriteTrackingGeometry.cxx | 5 + .../src/FaserActsWriteTrackingGeometry.h | 2 + .../components/FaserActsGeometry_entries.cxx | 17 +- .../test/FaserActsWriteTrackingGeometry.py | 15 +- 33 files changed, 1290 insertions(+), 1846 deletions(-) delete mode 100644 Tracking/Acts/FaserActsGeometry/FaserActsGeometry/FaserActsJsonGeometryConverter.h create mode 100644 Tracking/Acts/FaserActsGeometry/python/FaserActsMaterialValidationAlg.py delete mode 100644 Tracking/Acts/FaserActsGeometry/src/FaserActsJsonGeometryConverter.cxx create mode 100644 Tracking/Acts/FaserActsGeometry/src/FaserActsMaterialTrackWriterSvc.cxx create mode 100644 Tracking/Acts/FaserActsGeometry/src/FaserActsMaterialTrackWriterSvc.h diff --git a/Tracker/TrackerRecAlgs/TrackerSPFit/src/TrackerSPFit.cxx b/Tracker/TrackerRecAlgs/TrackerSPFit/src/TrackerSPFit.cxx index e78f7a8c..8daed1f9 100755 --- a/Tracker/TrackerRecAlgs/TrackerSPFit/src/TrackerSPFit.cxx +++ b/Tracker/TrackerRecAlgs/TrackerSPFit/src/TrackerSPFit.cxx @@ -1,6 +1,6 @@ /* Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration -*/ + */ /*************************************************************************** ------------------- @@ -111,6 +111,7 @@ namespace Tracker declareProperty("UseBiasedResidual", m_bias=true); declareProperty("SaveAllClusterInfo", m_saveallcluster=true); declareProperty("SaveAllSPInfo", m_saveallsp=true); + declareProperty("MakeDoublets", m_doublets=true); } //----------------------------------------------------------------------- @@ -135,6 +136,12 @@ namespace Tracker m_tree= new TTree("spfit","tree"); TrackerSPFit::initializeTree(); + ATH_CHECK( m_trackerSeedContainerKey.initialize() ); + if (m_trackerSeedContainerKey.key().empty()){ + ATH_MSG_FATAL( "No name set for output track seeds"); + return StatusCode::FAILURE; + } + // create containers (requires the Identifier Helpers) ATH_CHECK(detStore()->retrieve(m_idHelper,"FaserSCT_ID")); @@ -422,6 +429,12 @@ namespace Tracker return StatusCode::SUCCESS; } + + SG::WriteHandle<TrackerSeedCollection> seedContainer(m_trackerSeedContainerKey, ctx); + ATH_CHECK(seedContainer.record( std::make_unique<TrackerSeedCollection>() ) ); + ATH_MSG_DEBUG("Created track seed container " << m_trackerSeedContainerKey.key()); + + //read cluster SG::ReadHandle<Tracker::FaserSCT_ClusterContainer> sct_clcontainer( m_Sct_clcontainerKey, ctx ); if ((!sct_clcontainer.isValid())&&m_saveallcluster){ @@ -438,6 +451,7 @@ namespace Tracker const Tracker::FaserSCT_ClusterCollection *colNext=&(**clusit); FaserSCT_ClusterCollection::const_iterator clusters {colNext->begin()} ; FaserSCT_ClusterCollection::const_iterator clustersEnd {colNext->end() }; + ATH_MSG_DEBUG("loop over all clusters " ); for (; clusters != clustersEnd; ++clusters) { ++m_numberOfCluster; @@ -459,6 +473,7 @@ namespace Tracker m_cluster_all_side.push_back(m_idHelper->side(elementID)); m_cluster_all_strip.push_back(m_idHelper->strip(elementID)); m_cluster_all_size.push_back(clus->rdoList().size()); + ATH_MSG_DEBUG("cluster found in station/layer/phi/etamodule "<<m_idHelper->station(elementID)<<"/"<<m_idHelper->layer(elementID)<<"/"<<ietamodule<<"/"<<iphimodule ); } } @@ -466,13 +481,14 @@ namespace Tracker // retrieve SCT spacepoint container - SG::ReadHandle<FaserSCT_SpacePointContainer> sct_spcontainer( m_Sct_spcontainerKey, ctx ); + SG::ReadHandle<FaserSCT_SpacePointContainer> sct_spcontainer{ m_Sct_spcontainerKey, ctx }; if (!sct_spcontainer.isValid()){ msg(MSG:: FATAL) << "Could not find the data object "<< sct_spcontainer.name() << " !" << endmsg; return StatusCode::RECOVERABLE; } // save all the sp information + int nsp_sta[4] = {0}; if(m_saveallsp){ ATH_MSG_DEBUG( "SCT spacepoint container found: " << sct_spcontainer->size() << " collections" ); FaserSCT_SpacePointContainer::const_iterator it = sct_spcontainer->begin(); @@ -480,12 +496,15 @@ namespace Tracker for (; it != itend; ++it){ ++m_numberOfSPCol; + ATH_MSG_DEBUG( "loop over spacepoint collection " << m_numberOfSPCol ); const FaserSCT_SpacePointCollection *colNext=&(**it); FaserSCT_SpacePointCollection::const_iterator sp_begin= colNext->begin(); FaserSCT_SpacePointCollection::const_iterator sp_end= colNext->end(); Identifier elementID = colNext->identify(); + ATH_MSG_DEBUG("loop over all spacepoint in collection with size = "<<colNext->size()); for (; sp_begin != sp_end; ++sp_begin){ ++m_numberOfSP; + ATH_MSG_DEBUG("reading "<<m_numberOfSP<<" Spacepoint at station "<<m_idHelper->station(elementID)); const Tracker::FaserSCT_SpacePoint* sp=&(**sp_begin); m_sp_all_x_local.push_back(sp->localParameters().x()); m_sp_all_y_local.push_back(sp->localParameters().y()); @@ -493,6 +512,7 @@ namespace Tracker m_sp_all_y.push_back(sp->globalPosition().y()); m_sp_all_z.push_back(sp->globalPosition().z()); m_sp_all_station.push_back(m_idHelper->station(elementID)); + nsp_sta[m_idHelper->station(elementID)] += 1; m_sp_all_layer.push_back(m_idHelper->layer(elementID)); m_sp_all_identify0.push_back(sp->clusterList().first->identify().get_compact()); m_sp_all_identify1.push_back(sp->clusterList().second->identify().get_compact()); @@ -519,12 +539,13 @@ namespace Tracker FaserSCT_SpacePointContainer::const_iterator it = sct_spcontainer->begin(); FaserSCT_SpacePointContainer::const_iterator itend = sct_spcontainer->end(); - std::vector<SP_Seed> sp_sta0; - std::vector<SP_Seed> sp_sta1; - std::vector<SP_Seed> sp_sta2; - sp_sta0.clear(); - sp_sta1.clear(); - sp_sta2.clear(); + std::vector<std::vector<SP_Seed>> sp_sta; + sp_sta.resize(4); + sp_sta[0].resize(nsp_sta[0]); + sp_sta[1].resize(nsp_sta[1]); + sp_sta[2].resize(nsp_sta[2]); + sp_sta[3].resize(nsp_sta[3]); + int index_sp_sta[4] = {0}; for (; it != itend; ++it){ const FaserSCT_SpacePointCollection *colNext=&(**it); // int nReceivedSPSCT = colNext->size(); @@ -535,6 +556,7 @@ namespace Tracker ATH_MSG_VERBOSE("reading SpacePoints collection "<<elementID.get_compact()<<" , and hash = "<<colNext->identifyHash() ); int istation = m_idHelper->station(elementID); +// if(istation !=0) continue; int ilayer = m_idHelper->layer(elementID); int ietamodule = m_idHelper->eta_module(elementID); int iphimodule = m_idHelper->phi_module(elementID); @@ -556,217 +578,192 @@ namespace Tracker m_hist_eta->Fill(sp->eta()); m_hist_phi->Fill(sp->phi()); Amg::Vector3D gloPos=sp->globalPosition(); - Amg::Vector3D newgloPos; + if(gloPos.x()>150||gloPos.x()<-150||gloPos.y()>150||gloPos.y()<-150||gloPos.z()==0)continue; + if(sqrt(gloPos.x()*gloPos.x()+gloPos.y()*gloPos.y()+gloPos.z()*gloPos.z())<0.1)continue; m_hist_x->Fill(gloPos.x()); m_hist_y->Fill(gloPos.y()); m_hist_z->Fill(gloPos.z()); struct SP_Seed tmp{gloPos,sp->globCovariance(),ilayer,imodule}; - ATH_MSG_DEBUG( "TrackerSPFit algorithm found no space points" ); - if(istation==1) - sp_sta0.push_back(tmp); - else if(istation==2) - sp_sta1.push_back(tmp); - else - sp_sta2.push_back(tmp); - - if ( ((istation-1)*3+ilayer) == 0 ) m_hist_x_y_plane0->Fill(gloPos.x(),gloPos.y()); - if ( ((istation-1)*3+ilayer) == 1 ) m_hist_x_y_plane1->Fill(gloPos.x(),gloPos.y()); - if ( ((istation-1)*3+ilayer) == 2 ) m_hist_x_y_plane2->Fill(gloPos.x(),gloPos.y()); - if ( ((istation-1)*3+ilayer) == 3 ) m_hist_x_y_plane3->Fill(gloPos.x(),gloPos.y()); - if ( ((istation-1)*3+ilayer) == 4 ) m_hist_x_y_plane4->Fill(gloPos.x(),gloPos.y()); - if ( ((istation-1)*3+ilayer) == 5 ) m_hist_x_y_plane5->Fill(gloPos.x(),gloPos.y()); - if ( ((istation-1)*3+ilayer) == 6 ) m_hist_x_y_plane6->Fill(gloPos.x(),gloPos.y()); - if ( ((istation-1)*3+ilayer) == 7 ) m_hist_x_y_plane7->Fill(gloPos.x(),gloPos.y()); - if ( ((istation-1)*3+ilayer) == 8 ) m_hist_x_y_plane8->Fill(gloPos.x(),gloPos.y()); + sp_sta[istation][index_sp_sta[istation]]=tmp; + index_sp_sta[istation] += 1; + + if ( ((istation)*3+ilayer) == 0 ) m_hist_x_y_plane0->Fill(gloPos.x(),gloPos.y()); + if ( ((istation)*3+ilayer) == 1 ) m_hist_x_y_plane1->Fill(gloPos.x(),gloPos.y()); + if ( ((istation)*3+ilayer) == 2 ) m_hist_x_y_plane2->Fill(gloPos.x(),gloPos.y()); + if ( ((istation)*3+ilayer) == 3 ) m_hist_x_y_plane3->Fill(gloPos.x(),gloPos.y()); + if ( ((istation)*3+ilayer) == 4 ) m_hist_x_y_plane4->Fill(gloPos.x(),gloPos.y()); + if ( ((istation)*3+ilayer) == 5 ) m_hist_x_y_plane5->Fill(gloPos.x(),gloPos.y()); + if ( ((istation)*3+ilayer) == 6 ) m_hist_x_y_plane6->Fill(gloPos.x(),gloPos.y()); + if ( ((istation)*3+ilayer) == 7 ) m_hist_x_y_plane7->Fill(gloPos.x(),gloPos.y()); + if ( ((istation)*3+ilayer) == 8 ) m_hist_x_y_plane8->Fill(gloPos.x(),gloPos.y()); } ATH_MSG_VERBOSE( size << " SpacePoints successfully added to Container !" ); } } - ATH_MSG_INFO( "TrackerSPFit number of spacepoints in each stations: "<<sp_sta0.size()<<" "<<sp_sta1.size()<<" "<<sp_sta2.size() ); + ATH_MSG_INFO( "TrackerSPFit number of spacepoints in each stations: "<<sp_sta[0].size()<<" "<<sp_sta[1].size()<<" "<<sp_sta[2].size() ); bool ifTrack=false; - if(sp_sta0.size()>2&&sp_sta0.size()<100) { - auto tracks_sta0=makeTrackSeg(sp_sta0,m_maxchi2); - ATH_MSG_INFO( "TrackerSPFit found "<<tracks_sta0.size()<<" track segments in station 0 " ); - if(tracks_sta0.size()>0){ - ifTrack=true; - ATH_MSG_DEBUG( "TrackerSPFit fill the track info0 " ); - for(auto track:tracks_sta0){ - auto pre=track.predicted; - auto pos=track.pos; - ++m_numberOfTrack; - m_trackId.push_back(int(m_numberOfTrack/3)); - m_sp_track_x.push_back(pos.x()); - m_sp_track_y.push_back(pos.y()); - m_sp_track_z.push_back(pos.z()); - m_sp_track_x_err.push_back(track.err.x()); - m_sp_track_y_err.push_back(track.err.y()); - m_sp_track_z_err.push_back(track.err.z()); - m_sp_track_x_predicted.push_back(pre.x()); - m_sp_track_y_predicted.push_back(pre.y()); - m_sp_track_z_predicted.push_back(pre.z()); - m_sp_track_x_residual.push_back(pre.x()-pos.x()); - m_sp_track_y_residual.push_back(pre.y()-pos.y()); - m_sp_track_z_residual.push_back(pre.z()-pos.z()); - ATH_MSG_DEBUG(" fill the result "<<pre.x()<<" "<<pre.y()<<" "<<pre.z()<<" "<<pos.x()<<" "<<pos.y()<<" "<<pos.z()); - ATH_MSG_DEBUG(" Station/layer/module =1"<<"/"<<track.layer<<"/"<<track.module<<" , and track ID "<<m_numberOfTrack/3); - m_sp_track_station.push_back(0); - m_sp_track_layer.push_back(track.layer); - m_sp_track_module.push_back(track.module); - + for(int ista=0;ista<4;ista++){ + if(sp_sta[ista].size()>1&&sp_sta[ista].size()<100) { + auto tracks_sta0=makeTrackSeg(sp_sta[ista],m_maxchi2); + ATH_MSG_INFO( "TrackerSPFit found "<<tracks_sta0.size()<<" track segments in station "<<ista ); + if(tracks_sta0.size()>0){ + ifTrack=true; + ATH_MSG_DEBUG( "TrackerSPFit fill the track info0 " ); + + //make TrackerSeed if(m_bias){ - m_track_chi2.push_back(track.chi2); - m_track_edm.push_back(track.edm); - m_track_ndf.push_back(track.ndf); - m_track_p0.push_back(track.p0); - m_track_p1.push_back(track.p1); - m_track_p2.push_back(track.p2); - m_track_p3.push_back(track.p3); - } - if(track.layer==0){ m_residual_x_plane0->Fill(pre.x()-pos.x());m_residual_y_plane0->Fill(pre.y()-pos.y());} - if(track.layer==1){ m_residual_x_plane1->Fill(pre.x()-pos.x());m_residual_y_plane1->Fill(pre.y()-pos.y());} - if(track.layer==2){ m_residual_x_plane2->Fill(pre.x()-pos.x());m_residual_y_plane2->Fill(pre.y()-pos.y());} - if(track.layer==1){ - ATH_MSG_DEBUG(" only fill the layer 1, module = "<<track.module); - ATH_MSG_DEBUG(" measured = "<<pos.x()<<" , "<<pos.y()<<" , "<<pos.z()); - ATH_MSG_DEBUG(" predicted = "<<pre.x()<<" , "<<pre.y()<<" , "<<pre.z()); - if(track.module==0){ - m_preYhistX_s1l1m0->Fill(pos.x(),pos.y()-pre.y()); - m_preYhistY_s1l1m0->Fill(pos.y(),pos.y()-pre.y()); - m_preXhistX_s1l1m0->Fill(pos.x(),pos.x()-pre.x()); - m_preXhistY_s1l1m0->Fill(pos.y(),pos.x()-pre.x()); - } - if(track.module==1){ - m_preYhistX_s1l1m1->Fill(pos.x(),pos.y()-pre.y()); - m_preYhistY_s1l1m1->Fill(pos.y(),pos.y()-pre.y()); - m_preXhistX_s1l1m1->Fill(pos.x(),pos.x()-pre.x()); - m_preXhistY_s1l1m1->Fill(pos.y(),pos.x()-pre.x()); - } - if(track.module==2){ - m_preYhistX_s1l1m2->Fill(pos.x(),pos.y()-pre.y()); - m_preYhistY_s1l1m2->Fill(pos.y(),pos.y()-pre.y()); - m_preXhistX_s1l1m2->Fill(pos.x(),pos.x()-pre.x()); - m_preXhistY_s1l1m2->Fill(pos.y(),pos.x()-pre.x()); - } - if(track.module==3){ - m_preYhistX_s1l1m3->Fill(pos.x(),pos.y()-pre.y()); - m_preYhistY_s1l1m3->Fill(pos.y(),pos.y()-pre.y()); - m_preXhistX_s1l1m3->Fill(pos.x(),pos.x()-pre.x()); - m_preXhistY_s1l1m3->Fill(pos.y(),pos.x()-pre.x()); + for(size_t iseed=0;iseed<tracks_sta0.size();iseed+=3){ + Tracker::TrackerSeed* trackerSeed = new Tracker::TrackerSeed(); + trackerSeed->setStation(ista); + trackerSeed->set_id(TrackerSeed::TRIPLET_SP); + double ndf=tracks_sta0[iseed].ndf; + double* param=new double[6]; + param[0] = tracks_sta0[iseed].p0; + param[1] = tracks_sta0[iseed].p1; + param[2] = tracks_sta0[iseed].p2; + param[3] = tracks_sta0[iseed].p3; + param[4] = tracks_sta0[iseed].chi2; + param[5] = ndf; + // double* param[] = {tracks_sta0[iseed].p0, tracks_sta0[iseed].p1, tracks_sta0[iseed].p2, tracks_sta0[iseed].p3, tracks_sta0[iseed].chi2, ndf}; + trackerSeed->setParameters(param); + std::vector<const Tracker::FaserSCT_SpacePoint*> spseed; + spseed.clear(); + spseed.emplace_back(findSpacePoint(&*sct_spcontainer, tracks_sta0[iseed].pos.x(), tracks_sta0[iseed].pos.y(), tracks_sta0[iseed].pos.z())); + spseed.emplace_back(findSpacePoint(&*sct_spcontainer, tracks_sta0[iseed+1].pos.x(), tracks_sta0[iseed+1].pos.y(), tracks_sta0[iseed+1].pos.z())); + spseed.emplace_back(findSpacePoint(&*sct_spcontainer, tracks_sta0[iseed+2].pos.x(), tracks_sta0[iseed+2].pos.y(), tracks_sta0[iseed+2].pos.z())); + trackerSeed->add(spseed); + seedContainer->push_back(trackerSeed); } - if(track.module==4){ - m_preYhistX_s1l1m4->Fill(pos.x(),pos.y()-pre.y()); - m_preYhistY_s1l1m4->Fill(pos.y(),pos.y()-pre.y()); - m_preXhistX_s1l1m4->Fill(pos.x(),pos.x()-pre.x()); - m_preXhistY_s1l1m4->Fill(pos.y(),pos.x()-pre.x()); - } - if(track.module==5){ - m_preYhistX_s1l1m5->Fill(pos.x(),pos.y()-pre.y()); - m_preYhistY_s1l1m5->Fill(pos.y(),pos.y()-pre.y()); - m_preXhistX_s1l1m5->Fill(pos.x(),pos.x()-pre.x()); - m_preXhistY_s1l1m5->Fill(pos.y(),pos.x()-pre.x()); - } - if(track.module==6){ - m_preYhistX_s1l1m6->Fill(pos.x(),pos.y()-pre.y()); - m_preYhistY_s1l1m6->Fill(pos.y(),pos.y()-pre.y()); - m_preXhistX_s1l1m6->Fill(pos.x(),pos.x()-pre.x()); - m_preXhistY_s1l1m6->Fill(pos.y(),pos.x()-pre.x()); + + } + + for(auto track:tracks_sta0){ + auto pre=track.predicted; + auto pos=track.pos; + ++m_numberOfTrack; + m_trackId.push_back(int(m_numberOfTrack/3)); + m_sp_track_x.push_back(pos.x()); + m_sp_track_y.push_back(pos.y()); + m_sp_track_z.push_back(pos.z()); + m_sp_track_x_err.push_back(track.err.x()); + m_sp_track_y_err.push_back(track.err.y()); + m_sp_track_z_err.push_back(track.err.z()); + m_sp_track_x_predicted.push_back(pre.x()); + m_sp_track_y_predicted.push_back(pre.y()); + m_sp_track_z_predicted.push_back(pre.z()); + m_sp_track_x_residual.push_back(pre.x()-pos.x()); + m_sp_track_y_residual.push_back(pre.y()-pos.y()); + m_sp_track_z_residual.push_back(pre.z()-pos.z()); + ATH_MSG_DEBUG(" fill the result "<<pre.x()<<" "<<pre.y()<<" "<<pre.z()<<" "<<pos.x()<<" "<<pos.y()<<" "<<pos.z()); + ATH_MSG_DEBUG(" Station/layer/module = "<<ista<<"/"<<track.layer<<"/"<<track.module<<" , and track ID "<<m_numberOfTrack/3); + m_sp_track_station.push_back(0); + m_sp_track_layer.push_back(track.layer); + m_sp_track_module.push_back(track.module); + + if(m_bias){ + m_track_chi2.push_back(track.chi2); + m_track_edm.push_back(track.edm); + m_track_ndf.push_back(track.ndf); + m_track_p0.push_back(track.p0); + m_track_p1.push_back(track.p1); + m_track_p2.push_back(track.p2); + m_track_p3.push_back(track.p3); + } - if(track.module==7){ - m_preYhistX_s1l1m7->Fill(pos.x(),pos.y()-pre.y()); - m_preYhistY_s1l1m7->Fill(pos.y(),pos.y()-pre.y()); - m_preXhistX_s1l1m7->Fill(pos.x(),pos.x()-pre.x()); - m_preXhistY_s1l1m7->Fill(pos.y(),pos.x()-pre.x()); + if(track.layer==0){ m_residual_x_plane0->Fill(pre.x()-pos.x());m_residual_y_plane0->Fill(pre.y()-pos.y());} + if(track.layer==1){ m_residual_x_plane1->Fill(pre.x()-pos.x());m_residual_y_plane1->Fill(pre.y()-pos.y());} + if(track.layer==2){ m_residual_x_plane2->Fill(pre.x()-pos.x());m_residual_y_plane2->Fill(pre.y()-pos.y());} + if(track.layer==1){ + ATH_MSG_DEBUG(" only fill the layer 1, module = "<<track.module); + ATH_MSG_DEBUG(" measured = "<<pos.x()<<" , "<<pos.y()<<" , "<<pos.z()); + ATH_MSG_DEBUG(" predicted = "<<pre.x()<<" , "<<pre.y()<<" , "<<pre.z()); + if(track.module==0){ + m_preYhistX_s1l1m0->Fill(pos.x(),pos.y()-pre.y()); + m_preYhistY_s1l1m0->Fill(pos.y(),pos.y()-pre.y()); + m_preXhistX_s1l1m0->Fill(pos.x(),pos.x()-pre.x()); + m_preXhistY_s1l1m0->Fill(pos.y(),pos.x()-pre.x()); + } + if(track.module==1){ + m_preYhistX_s1l1m1->Fill(pos.x(),pos.y()-pre.y()); + m_preYhistY_s1l1m1->Fill(pos.y(),pos.y()-pre.y()); + m_preXhistX_s1l1m1->Fill(pos.x(),pos.x()-pre.x()); + m_preXhistY_s1l1m1->Fill(pos.y(),pos.x()-pre.x()); + } + if(track.module==2){ + m_preYhistX_s1l1m2->Fill(pos.x(),pos.y()-pre.y()); + m_preYhistY_s1l1m2->Fill(pos.y(),pos.y()-pre.y()); + m_preXhistX_s1l1m2->Fill(pos.x(),pos.x()-pre.x()); + m_preXhistY_s1l1m2->Fill(pos.y(),pos.x()-pre.x()); + } + if(track.module==3){ + m_preYhistX_s1l1m3->Fill(pos.x(),pos.y()-pre.y()); + m_preYhistY_s1l1m3->Fill(pos.y(),pos.y()-pre.y()); + m_preXhistX_s1l1m3->Fill(pos.x(),pos.x()-pre.x()); + m_preXhistY_s1l1m3->Fill(pos.y(),pos.x()-pre.x()); + } + if(track.module==4){ + m_preYhistX_s1l1m4->Fill(pos.x(),pos.y()-pre.y()); + m_preYhistY_s1l1m4->Fill(pos.y(),pos.y()-pre.y()); + m_preXhistX_s1l1m4->Fill(pos.x(),pos.x()-pre.x()); + m_preXhistY_s1l1m4->Fill(pos.y(),pos.x()-pre.x()); + } + if(track.module==5){ + m_preYhistX_s1l1m5->Fill(pos.x(),pos.y()-pre.y()); + m_preYhistY_s1l1m5->Fill(pos.y(),pos.y()-pre.y()); + m_preXhistX_s1l1m5->Fill(pos.x(),pos.x()-pre.x()); + m_preXhistY_s1l1m5->Fill(pos.y(),pos.x()-pre.x()); + } + if(track.module==6){ + m_preYhistX_s1l1m6->Fill(pos.x(),pos.y()-pre.y()); + m_preYhistY_s1l1m6->Fill(pos.y(),pos.y()-pre.y()); + m_preXhistX_s1l1m6->Fill(pos.x(),pos.x()-pre.x()); + m_preXhistY_s1l1m6->Fill(pos.y(),pos.x()-pre.x()); + } + if(track.module==7){ + m_preYhistX_s1l1m7->Fill(pos.x(),pos.y()-pre.y()); + m_preYhistY_s1l1m7->Fill(pos.y(),pos.y()-pre.y()); + m_preXhistX_s1l1m7->Fill(pos.x(),pos.x()-pre.x()); + m_preXhistY_s1l1m7->Fill(pos.y(),pos.x()-pre.x()); + } } } } - } - } - if(sp_sta1.size()>2&&sp_sta1.size()<100) { - auto tracks_sta1=makeTrackSeg(sp_sta1,m_maxchi2); - ATH_MSG_INFO( "TrackerSPFit found "<<tracks_sta1.size()<<" track segments in station 1 " ); - if(tracks_sta1.size()>0){ - ifTrack=true; - for(auto track:tracks_sta1){ - auto pre=track.predicted; - auto pos=track.pos; - ++m_numberOfTrack; - m_trackId.push_back(int(m_numberOfTrack/3)); - m_sp_track_x.push_back(pos.x()); - m_sp_track_y.push_back(pos.y()); - m_sp_track_z.push_back(pos.z()); - m_sp_track_x_err.push_back(track.err.x()); - m_sp_track_y_err.push_back(track.err.y()); - m_sp_track_z_err.push_back(track.err.z()); - m_sp_track_x_predicted.push_back(pre.x()); - m_sp_track_y_predicted.push_back(pre.y()); - m_sp_track_z_predicted.push_back(pre.z()); - m_sp_track_x_residual.push_back(pre.x()-pos.x()); - m_sp_track_y_residual.push_back(pre.y()-pos.y()); - m_sp_track_z_residual.push_back(pre.z()-pos.z()); - m_sp_track_station.push_back(1); - m_sp_track_layer.push_back(track.layer); - m_sp_track_module.push_back(track.module); - if(m_bias){ - m_track_chi2.push_back(track.chi2); - m_track_edm.push_back(track.edm); - m_track_ndf.push_back(track.ndf); - m_track_p0.push_back(track.p0); - m_track_p1.push_back(track.p1); - m_track_p2.push_back(track.p2); - m_track_p3.push_back(track.p3); - } - if(track.layer==0){ m_residual_x_plane3->Fill(pre.x()-pos.x());m_residual_y_plane3->Fill(pre.y()-pos.y());} - if(track.layer==1){ m_residual_x_plane4->Fill(pre.x()-pos.x());m_residual_y_plane4->Fill(pre.y()-pos.y());} - if(track.layer==2){ m_residual_x_plane5->Fill(pre.x()-pos.x());m_residual_y_plane5->Fill(pre.y()-pos.y());} - } - } - } - if(sp_sta2.size()>2&&sp_sta2.size()<100) { - auto tracks_sta2=makeTrackSeg(sp_sta2,m_maxchi2); - ATH_MSG_INFO( "TrackerSPFit found "<<tracks_sta2.size()<<" track segments in station 2 " ); - if(tracks_sta2.size()>0){ - ifTrack=true; - for(auto track:tracks_sta2){ - auto pre=track.predicted; - auto pos=track.pos; - ++m_numberOfTrack; - m_trackId.push_back(int(m_numberOfTrack/3)); - m_sp_track_x.push_back(pos.x()); - m_sp_track_y.push_back(pos.y()); - m_sp_track_z.push_back(pos.z()); - m_sp_track_x_err.push_back(track.err.x()); - m_sp_track_y_err.push_back(track.err.y()); - m_sp_track_z_err.push_back(track.err.z()); - m_sp_track_x_predicted.push_back(pre.x()); - m_sp_track_y_predicted.push_back(pre.y()); - m_sp_track_z_predicted.push_back(pre.z()); - m_sp_track_x_residual.push_back(pre.x()-pos.x()); - m_sp_track_y_residual.push_back(pre.y()-pos.y()); - m_sp_track_z_residual.push_back(pre.z()-pos.z()); - m_sp_track_station.push_back(2); - m_sp_track_layer.push_back(track.layer); - m_sp_track_module.push_back(track.module); - if(m_bias){ - m_track_chi2.push_back(track.chi2); - m_track_edm.push_back(track.edm); - m_track_ndf.push_back(track.ndf); - m_track_p0.push_back(track.p0); - m_track_p1.push_back(track.p1); - m_track_p2.push_back(track.p2); - m_track_p3.push_back(track.p3); + else if(m_doublets){ + for(unsigned int isp=0;isp<sp_sta[ista].size()-1;isp++){ + for(unsigned int jsp=isp+1;jsp<sp_sta[ista].size();jsp++){ + if(sp_sta[ista][isp].layer != sp_sta[ista][jsp].layer){ + Tracker::TrackerSeed* trackerSeed = new Tracker::TrackerSeed(); + trackerSeed->setStation(ista); + trackerSeed->set_id(TrackerSeed::DOUBLET_SP); + double* param = new double[6]; + double tmp_x_slope = (sp_sta[ista][jsp].pos.x()-sp_sta[ista][isp].pos.x())/(sp_sta[ista][jsp].pos.z()-sp_sta[ista][isp].pos.z()); + double tmp_x_offset = sp_sta[ista][isp].pos.x()-sp_sta[ista][isp].pos.x()*tmp_x_slope; + double tmp_y_slope = (sp_sta[ista][jsp].pos.y()-sp_sta[ista][isp].pos.y())/(sp_sta[ista][jsp].pos.z()-sp_sta[ista][isp].pos.z()); + double tmp_y_offset = sp_sta[ista][isp].pos.y()-sp_sta[ista][isp].pos.z()*tmp_y_slope; + param[0] = tmp_x_offset; + param[1] = tmp_x_slope; + param[2] = tmp_y_offset; + param[3] = tmp_y_slope; + param[4] = -999; + param[5] = -999; + trackerSeed->setParameters(param); + std::vector<const Tracker::FaserSCT_SpacePoint*> spseed; + spseed.clear(); + spseed.emplace_back(findSpacePoint(&*sct_spcontainer, sp_sta[ista][isp].pos.x(), sp_sta[ista][isp].pos.y(), sp_sta[ista][isp].pos.z())); + spseed.emplace_back(findSpacePoint(&*sct_spcontainer, sp_sta[ista][jsp].pos.x(), sp_sta[ista][jsp].pos.y(), sp_sta[ista][jsp].pos.z())); + trackerSeed->add(spseed); + seedContainer->push_back(trackerSeed); + } + } } - if(track.layer==0){ m_residual_x_plane6->Fill(pre.x()-pos.x());m_residual_y_plane6->Fill(pre.y()-pos.y());} - if(track.layer==1){ m_residual_x_plane7->Fill(pre.x()-pos.x());m_residual_y_plane7->Fill(pre.y()-pos.y());} - if(track.layer==2){ m_residual_x_plane8->Fill(pre.x()-pos.x());m_residual_y_plane8->Fill(pre.y()-pos.y());} } } } - if(ifTrack) - m_tree->Fill(); - else if(m_saveallcluster||m_saveallsp) + if(ifTrack||m_saveallcluster||m_saveallsp) m_tree->Fill(); + return StatusCode::SUCCESS; } @@ -907,7 +904,6 @@ namespace Tracker m_ndf->Fill(ndf); if(chi2<maxchi2){ - // result.Print(std::cout); const double *fitParam=result.GetParams(); ATH_MSG_DEBUG(" fit status: ook "<<chi2<<" "<<edm<<" "<<ndf<<" "<<fitParam[0]<<" "<<fitParam[1]<<" "<<fitParam[2]<<" "<<fitParam[3]); Amg::Vector3D err0(sqrt((sp0.cov)(0,0)),sqrt((sp0.cov)(1,1)),sqrt((sp0.cov)(2,2))); @@ -929,6 +925,22 @@ namespace Tracker return Amg::Vector3D(p[0]+p[1]*z,p[2]+p[3]*z,z); } + const FaserSCT_SpacePoint* TrackerSPFit::findSpacePoint(const FaserSCT_SpacePointContainer* spcontainer, double x, double y, double z) const{ + FaserSCT_SpacePointContainer::const_iterator it = spcontainer->begin(); + FaserSCT_SpacePointContainer::const_iterator itend = spcontainer->end(); + + for (; it != itend; ++it){ + const FaserSCT_SpacePointCollection *colNext=&(**it); + FaserSCT_SpacePointCollection::const_iterator sp_begin= colNext->begin(); + FaserSCT_SpacePointCollection::const_iterator sp_end= colNext->end(); + Identifier elementID = colNext->identify(); + for (; sp_begin != sp_end; ++sp_begin){ + const FaserSCT_SpacePoint* sp=&(**sp_begin); + if(pow(sp->globalPosition().x()-x,2)+pow(sp->globalPosition().y()-y,2)+pow(sp->globalPosition().z()-z,2)<0.001)return sp; + } + } + return nullptr; + } /* * need to be updated diff --git a/Tracker/TrackerRecAlgs/TrackerSPFit/src/TrackerSPFit.h b/Tracker/TrackerRecAlgs/TrackerSPFit/src/TrackerSPFit.h index 3db42689..40c06ec4 100755 --- a/Tracker/TrackerRecAlgs/TrackerSPFit/src/TrackerSPFit.h +++ b/Tracker/TrackerRecAlgs/TrackerSPFit/src/TrackerSPFit.h @@ -18,6 +18,8 @@ #include "TrackerSpacePoint/FaserSCT_SpacePoint.h" #include "TrackerSpacePoint/FaserSCT_SpacePointContainer.h" #include "TrackerSpacePoint/FaserSCT_SpacePointOverlapCollection.h" +#include "TrackerSpacePoint/TrackerSeedCollection.h" +#include "TrackerSpacePoint/TrackerSeed.h" #include "TrkEventPrimitives/LocalParameters.h" #include "GaudiKernel/ServiceHandle.h" @@ -112,6 +114,8 @@ namespace Tracker SG::ReadHandleKey<Tracker::FaserSCT_ClusterContainer> m_Sct_clcontainerKey{this, "SCT_ClustersName", "SCT clContainer"}; SG::ReadCondHandleKey<TrackerDD::SiDetectorElementCollection> m_SCTDetEleCollKey{this, "SCTDetEleCollKey", "SCT_DetectorElementCollection", "Key of SiDetectorElementCollection for SCT"}; + SG::WriteHandleKey<Tracker::TrackerSeedCollection> m_trackerSeedContainerKey{this, "FaserTrackerSeedName", "FaserTrackerSeedCollection" }; + const Tracker::FaserSCT_SpacePoint* findSpacePoint(const FaserSCT_SpacePointContainer* spcontainer, double x, double y, double z) const; const FaserSCT_ID* m_idHelper{nullptr}; mutable std::atomic<int> m_numberOfEvents{0}; @@ -247,6 +251,7 @@ namespace Tracker mutable std::vector<int> m_sp_all_module; mutable std::vector<long long int> m_sp_all_identify0; mutable std::vector<long long int> m_sp_all_identify1; + bool m_doublets; }; diff --git a/Tracker/TrackerRecEvent/TrackerSpacePoint/TrackerSpacePoint/TrackerSeed.h b/Tracker/TrackerRecEvent/TrackerSpacePoint/TrackerSpacePoint/TrackerSeed.h index 80983149..dea1aa95 100755 --- a/Tracker/TrackerRecEvent/TrackerSpacePoint/TrackerSpacePoint/TrackerSeed.h +++ b/Tracker/TrackerRecEvent/TrackerSpacePoint/TrackerSpacePoint/TrackerSeed.h @@ -15,13 +15,13 @@ namespace Tracker { public: - enum StrategyId{NULLID=0, TRIPLET_SP_FIRSTSTATION=1}; + enum StrategyId{NULLID=0, TRIPLET_SP_FIRSTSTATION=1, TRIPLET_SP=2, DOUBLET_SP=3};//, DOUBLET=2}; TrackerSeed(); TrackerSeed(const StrategyId, const TrackerSeed &); virtual ~TrackerSeed(); - TrackerSeed(const StrategyId, vector<const FaserSCT_SpacePoint*> seed); + TrackerSeed(const StrategyId, vector<const FaserSCT_SpacePoint*> seed, int stationId= 1 , double* param=0); void set_id(const StrategyId id) { m_strategyId = id; } StrategyId id() const { return m_strategyId; } @@ -30,6 +30,10 @@ namespace Tracker { vector<const FaserSCT_SpacePoint*> getSpacePoints() const {return m_seed;}; int size() const; + double* getParameters() const{return m_param;}; + int getStation() const{return m_stationId;}; + void setParameters(double* param); + void setStation(int station); TrackerSeed &operator=(const TrackerSeed &); @@ -40,6 +44,8 @@ namespace Tracker { StrategyId m_strategyId; vector<const FaserSCT_SpacePoint*> m_seed; + double* m_param; + int m_stationId; }; diff --git a/Tracker/TrackerRecEvent/TrackerSpacePoint/src/TrackerSeed.cxx b/Tracker/TrackerRecEvent/TrackerSpacePoint/src/TrackerSeed.cxx index 2e0f2551..5d9f9e3b 100755 --- a/Tracker/TrackerRecEvent/TrackerSpacePoint/src/TrackerSeed.cxx +++ b/Tracker/TrackerRecEvent/TrackerSpacePoint/src/TrackerSeed.cxx @@ -4,32 +4,38 @@ namespace Tracker { TrackerSeed::TrackerSeed() : m_strategyId(NULLID) {} - TrackerSeed::TrackerSeed(const StrategyId id, const TrackerSeed& trackerSeed) : m_strategyId(id), m_seed(trackerSeed.m_seed) {} + TrackerSeed::TrackerSeed(const StrategyId id, const TrackerSeed& trackerSeed) : m_strategyId(id), m_seed(trackerSeed.m_seed), m_param(trackerSeed.m_param), m_stationId(trackerSeed.m_stationId) {} TrackerSeed::~TrackerSeed() {} - TrackerSeed::TrackerSeed(const StrategyId id, vector<const FaserSCT_SpacePoint*> seed) { m_strategyId = id; m_seed = seed; } + TrackerSeed::TrackerSeed(const StrategyId id, vector<const FaserSCT_SpacePoint*> seed, int stationId, double* param) { m_strategyId = id; m_seed = seed; m_param = param; m_stationId = stationId; } void TrackerSeed::add(vector<const FaserSCT_SpacePoint*> seed) { m_seed = seed; } int TrackerSeed::size() const { return m_seed.size(); } + void TrackerSeed::setParameters(double* param) { m_param = param; } + + void TrackerSeed::setStation(int station) { m_stationId= station; } + TrackerSeed& TrackerSeed::operator=(const TrackerSeed& trackSeed){ if(&trackSeed != this) { TrackerSeed::operator=(trackSeed); m_seed = trackSeed.m_seed; + m_param = trackSeed.m_param; + m_stationId = trackSeed.m_stationId; } return *this; } MsgStream& TrackerSeed::dump(MsgStream& stream) const { - stream << "TrackerSeed object" << endl; + stream << "TrackerSeed object at station " <<m_stationId<< " with strategy = "<<m_strategyId<< endl; this->TrackerSeed::dump(stream); return stream; } ostream& TrackerSeed::dump(ostream& stream) const { - stream << "TrackerSeed object" << endl; + stream << "TrackerSeed object at station " <<m_stationId<< " with strategy = "<<m_strategyId<< endl; this->TrackerSeed::dump(stream); return stream; } diff --git a/Tracking/Acts/FaserActsGeometry/CMakeLists.txt b/Tracking/Acts/FaserActsGeometry/CMakeLists.txt index 48c9abd9..e15adcec 100755 --- a/Tracking/Acts/FaserActsGeometry/CMakeLists.txt +++ b/Tracking/Acts/FaserActsGeometry/CMakeLists.txt @@ -8,92 +8,76 @@ find_package( Eigen ) find_package( Boost ) find_package( nlohmann_json ) -#set( Acts_DIR /home/tboeckh/Documents/acts/run ) -#find_package( Acts REQUIRED COMPONENTS Core PATHS /home/tboeckh/Documents/acts/run NO_DEFAULT_PATH ) find_package( Acts COMPONENTS Core ) atlas_add_library( FaserActsGeometryLib - src/FaserActsSurfaceMappingTool.cxx - src/FaserActsMaterialMapping.cxx - src/FaserActsAlignmentStore.cxx - src/FaserActsDetectorElement.cxx - src/FaserActsLayerBuilder.cxx - src/CuboidVolumeBuilder.cxx -# src/FaserActsJsonGeometryConverter.cxx - src/util/*.cxx - PUBLIC_HEADERS FaserActsGeometry - INCLUDE_DIRS ${CLHEP_INCLUDE_DIRS} ${EIGEN_INCLUDE_DIRS} ${BOOST_INCLUDE_DIRS} - LINK_LIBRARIES ${CLHEP_LIBRARIES} ${EIGEN_LIBRARIES} ${BOOST_LIBRARIES} nlohmann_json::nlohmann_json - TrackerIdentifier - TrackerReadoutGeometry - ActsInteropLib - FaserActsGeometryInterfacesLib - AthenaKernel - ActsCore - GeoModelFaserUtilities - ActsGeometryLib - ActsGeometryInterfacesLib - MagFieldInterfaces - MagFieldElements - MagFieldConditions - TrackerRawData - TrackerSimData - GeneratorObjects - TrackerSimEvent - TrackerSpacePoint - TrackerIdentifier - TrackerPrepRawData - TrkSpacePoint - TrkGeometry + src/FaserActsSurfaceMappingTool.cxx + src/FaserActsMaterialJsonWriterTool.cxx + src/FaserActsMaterialMapping.cxx + src/FaserActsAlignmentStore.cxx + src/FaserActsDetectorElement.cxx + src/FaserActsLayerBuilder.cxx + src/CuboidVolumeBuilder.cxx + src/util/*.cxx + PUBLIC_HEADERS FaserActsGeometry + INCLUDE_DIRS ${CLHEP_INCLUDE_DIRS} ${EIGEN_INCLUDE_DIRS} ${BOOST_INCLUDE_DIRS} + LINK_LIBRARIES ${CLHEP_LIBRARIES} ${EIGEN_LIBRARIES} ${BOOST_LIBRARIES} nlohmann_json::nlohmann_json + TrackerIdentifier + TrackerReadoutGeometry + ActsInteropLib + FaserActsGeometryInterfacesLib + AthenaKernel + ActsCore + GeoModelFaserUtilities + ActsGeometryLib + ActsGeometryInterfacesLib + MagFieldInterfaces + MagFieldElements + MagFieldConditions + GeneratorObjects + TrackerIdentifier + TrkGeometry ) # Component(s) in the package: atlas_add_component( FaserActsGeometry -##src/*.cxx - src/FaserActsSurfaceMappingTool.cxx - src/FaserActsMaterialMapping.cxx -# src/FaserActsMaterialJsonWriterTool.cxx -# src/FaserActsJsonGeometryConverter.cxx - src/FaserActsTrackingGeometrySvc.cxx - src/FaserActsTrackingGeometryTool.cxx - src/FaserActsWriteTrackingGeometry.cxx - src/FaserActsObjWriterTool.cxx - src/FaserActsExtrapolationAlg.cxx - src/FaserActsExtrapolationTool.cxx - src/FaserActsPropStepRootWriterSvc.cxx - src/FaserActsAlignmentCondAlg.cxx - src/NominalAlignmentCondAlg.cxx -#src/FaserActsKalmanFilterAlg.cxx - src/FaserActsVolumeMappingTool.cxx - src/components/*.cxx - PUBLIC_HEADERS FaserActsGeometry - INCLUDE_DIRS ${CLHEP_INCLUDE_DIRS} ${EIGEN_INCLUDE_DIRS} ${BOOST_INCLUDE_DIRS} - LINK_LIBRARIES ${CLHEP_LIBRARIES} ${EIGEN_LIBRARIES} ${BOOST_LIBRARIES} - EventInfo - AthenaBaseComps - GaudiKernel - FaserActsGeometryLib - ActsInteropLib - FaserActsGeometryInterfacesLib - ActsCore - ActsGeometryLib - ActsGeometryInterfacesLib - GeoModelFaserUtilities - MagFieldInterfaces - MagFieldElements - MagFieldConditions - TrackerRawData - TrackerSimData - GeneratorObjects - TrackerSimEvent - TrackerSpacePoint - TrackerIdentifier - TrackerPrepRawData - TrkSpacePoint - TrkGeometry + src/FaserActsSurfaceMappingTool.cxx + src/FaserActsMaterialMapping.cxx + src/FaserActsMaterialJsonWriterTool.cxx + src/FaserActsTrackingGeometrySvc.cxx + src/FaserActsTrackingGeometryTool.cxx + src/FaserActsWriteTrackingGeometry.cxx + src/FaserActsObjWriterTool.cxx + src/FaserActsExtrapolationAlg.cxx + src/FaserActsExtrapolationTool.cxx + src/FaserActsPropStepRootWriterSvc.cxx + src/FaserActsMaterialTrackWriterSvc.cxx + src/FaserActsAlignmentCondAlg.cxx + src/NominalAlignmentCondAlg.cxx + src/FaserActsVolumeMappingTool.cxx + src/components/*.cxx + PUBLIC_HEADERS FaserActsGeometry + INCLUDE_DIRS ${CLHEP_INCLUDE_DIRS} ${EIGEN_INCLUDE_DIRS} ${BOOST_INCLUDE_DIRS} + LINK_LIBRARIES ${CLHEP_LIBRARIES} ${EIGEN_LIBRARIES} ${BOOST_LIBRARIES} + EventInfo + AthenaBaseComps + GaudiKernel + FaserActsGeometryLib + ActsInteropLib + FaserActsGeometryInterfacesLib + ActsCore + ActsGeometryLib + ActsGeometryInterfacesLib + GeoModelFaserUtilities + MagFieldInterfaces + MagFieldElements + MagFieldConditions + GeneratorObjects + TrackerIdentifier + TrkGeometry - ) +) # Install files from the package: atlas_install_headers( FaserActsGeometry ) diff --git a/Tracking/Acts/FaserActsGeometry/FaserActsGeometry/FaserActsJsonGeometryConverter.h b/Tracking/Acts/FaserActsGeometry/FaserActsGeometry/FaserActsJsonGeometryConverter.h deleted file mode 100644 index f665a575..00000000 --- a/Tracking/Acts/FaserActsGeometry/FaserActsGeometry/FaserActsJsonGeometryConverter.h +++ /dev/null @@ -1,261 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2017-2019 CERN for the benefit of the Acts project -// -// This Source Code Form is subject to the terms of the Mozilla Public -// License, v. 2.0. If a copy of the MPL was not distributed with this -// file, You can obtain one at http://mozilla.org/MPL/2.0/. - -#pragma once - -#include "Acts/Geometry/TrackingGeometry.hpp" -#include "Acts/Material/ISurfaceMaterial.hpp" -#include "Acts/Material/IVolumeMaterial.hpp" -#include "Acts/Material/MaterialSlab.hpp" -#include "Acts/Utilities/BinUtility.hpp" -#include "Acts/Definitions/Common.hpp" -#include "Acts/Utilities/Logger.hpp" -#include <Acts/Geometry/TrackingVolume.hpp> -#include <Acts/Surfaces/Surface.hpp> - -#include <map> - -#include <nlohmann/json.hpp> - -namespace Acts { - -/// @class FaserActsJsonGeometryConverter -/// -/// @brief read the material from Json -class FaserActsJsonGeometryConverter { - public: - using SurfaceMaterialMap = - std::map<GeometryIdentifier, std::shared_ptr<const ISurfaceMaterial>>; - - using VolumeMaterialMap = - std::map<GeometryIdentifier, std::shared_ptr<const IVolumeMaterial>>; - - using DetectorMaterialMaps = std::pair<SurfaceMaterialMap, VolumeMaterialMap>; - - using geo_id_value = uint64_t; - - using SurfaceMaterialRep = std::map<geo_id_value, const ISurfaceMaterial*>; - using SurfaceRep = std::map<geo_id_value, const Surface*>; - using VolumeMaterialRep = std::map<geo_id_value, const IVolumeMaterial*>; - - /// @brief Layer representation for Json writing - struct LayerRep { - // the layer id - GeometryIdentifier layerID; - - SurfaceMaterialRep sensitives; - SurfaceRep sensitiveSurfaces; - SurfaceMaterialRep approaches; - SurfaceRep approacheSurfaces; - const ISurfaceMaterial* representing = nullptr; - const Surface* representingSurface = nullptr; - - /// The LayerRep is actually worth it to write out - operator bool() const { - return (!sensitives.empty() or !approaches.empty() or - representing != nullptr); - } - }; - - /// @brief Volume representation for Json writing - struct VolumeRep { - // The geometry id - GeometryIdentifier volumeID; - - /// The namne - std::string volumeName; - - std::map<geo_id_value, LayerRep> layers; - SurfaceMaterialRep boundaries; - SurfaceRep boundarySurfaces; - const IVolumeMaterial* material = nullptr; - - /// The VolumeRep is actually worth it to write out - operator bool() const { - return (!layers.empty() or !boundaries.empty() or material != nullptr); - } - }; - - /// @brief Detector representation for Json writing - struct DetectorRep { - std::map<geo_id_value, VolumeRep> volumes; - }; - - /// @class Config - /// Configuration of the Reader - class Config { - public: - /// The geometry version - std::string geoversion = "undefined"; - /// The detector tag - std::string detkey = "detector"; - /// The volume identification string - std::string volkey = "volumes"; - /// The name identification - std::string namekey = "Name"; - /// The boundary surface string - std::string boukey = "boundaries"; - /// The layer identification string - std::string laykey = "layers"; - /// The volume material string - std::string matkey = "material"; - /// The approach identification string - std::string appkey = "approach"; - /// The sensitive identification string - std::string senkey = "sensitive"; - /// The representing idntification string - std::string repkey = "representing"; - /// The bin keys - std::string bin0key = "bin0"; - /// The bin1 key - std::string bin1key = "bin1"; - /// The bin2 key - std::string bin2key = "bin2"; - /// The local to global tranfo key - std::string transfokeys = "tranformation"; - /// The type key -> proto, else - std::string typekey = "type"; - /// The data key - std::string datakey = "data"; - /// The geoid key - std::string geometryidkey = "Geoid"; - /// The surface geoid key - std::string surfacegeometryidkey = "SGeoid"; - /// The mapping key, add surface to map if true - std::string mapkey = "mapMaterial"; - /// The surface type key - std::string surfacetypekey = "stype"; - /// The surface position key - std::string surfacepositionkey = "sposition"; - /// The surface range key - std::string surfacerangekey = "srange"; - /// The default logger - std::shared_ptr<const Logger> logger; - /// The name of the writer - std::string name = ""; - - /// Steering to handle sensitive data - bool processSensitives = true; - /// Steering to handle approach data - bool processApproaches = true; - /// Steering to handle representing data - bool processRepresenting = true; - /// Steering to handle boundary data - bool processBoundaries = true; - /// Steering to handle volume data - bool processVolumes = true; - /// Steering to handle volume data - bool processDenseVolumes = false; - /// Add proto material to all surfaces - bool processnonmaterial = false; - /// Write out data - bool writeData = true; - - /// Constructor - /// - /// @param lname Name of the writer tool - /// @param lvl The output logging level - Config(const std::string& lname = "FaserActsJsonGeometryConverter", - Logging::Level lvl = Logging::INFO) - : logger(getDefaultLogger(lname, lvl)), name(lname) {} - }; - - /// Constructor - /// - /// @param cfg configuration struct for the reader - FaserActsJsonGeometryConverter(const Config& cfg); - - /// Destructor - ~FaserActsJsonGeometryConverter() = default; - - /// Convert method - /// - /// @param surfaceMaterialMap The indexed material map collection - std::pair< - std::map<GeometryIdentifier, std::shared_ptr<const ISurfaceMaterial>>, - std::map<GeometryIdentifier, std::shared_ptr<const IVolumeMaterial>>> - jsonToMaterialMaps(const nlohmann::json& materialmaps); - - /// Convert method - /// - /// @param surfaceMaterialMap The indexed material map collection - nlohmann::json materialMapsToJson(const DetectorMaterialMaps& maps); - - /// Write method - /// - /// @param tGeometry is the tracking geometry which contains the material - nlohmann::json trackingGeometryToJson(const TrackingGeometry& tGeometry); - - private: - /// Convert to internal representation method, recursive call - /// - /// @param tGeometry is the tracking geometry which contains the material - void convertToRep(DetectorRep& detRep, const TrackingVolume& tVolume); - - /// Convert to internal representation method - /// - /// @param tGeometry is the tracking geometry which contains the material - LayerRep convertToRep(const Layer& tLayer); - - /// Create the Surface Material from Json - /// - factory method, ownership given - /// @param material is the json part representing a material object - const ISurfaceMaterial* jsonToSurfaceMaterial(const nlohmann::json& material); - - /// Create the Volume Material from Json - /// - factory method, ownership given - /// @param material is the json part representing a material object - const IVolumeMaterial* jsonToVolumeMaterial(const nlohmann::json& material); - - /// Create the Material Matrix from Json - /// - /// @param data is the json part representing a material data array - MaterialSlabMatrix jsonToMaterialMatrix(const nlohmann::json& data); - - /// Create the BinUtility for from Json - BinUtility jsonToBinUtility(const nlohmann::json& bin); - - /// Create the local to global transform for from Json - Transform3 jsonToTransform(const nlohmann::json& transfo); - - /// Create Json from a detector represenation - nlohmann::json detectorRepToJson(const DetectorRep& detRep); - - /// SurfaceMaterial to Json - /// - /// @param the SurfaceMaterial - nlohmann::json surfaceMaterialToJson(const ISurfaceMaterial& sMaterial); - - /// VolumeMaterial to Json - /// - /// @param the VolumeMaterial - nlohmann::json volumeMaterialToJson(const IVolumeMaterial& vMaterial); - - /// Add surface information to json surface - /// - /// @param The json surface The surface - void addSurfaceToJson(nlohmann::json& sjson, const Surface* surface); - - /// Default BinUtility to create proto material - /// - /// @param the Surface - Acts::BinUtility DefaultBin(const Acts::Surface& surface); - - /// Default BinUtility to create proto material - /// - /// @param the Volume - Acts::BinUtility DefaultBin(const Acts::TrackingVolume& volume); - - /// The config class - Config m_cfg; - - /// Private access to the logging instance - const Logger& logger() const { return *m_cfg.logger; } -}; - -} // namespace Acts diff --git a/Tracking/Acts/FaserActsGeometry/FaserActsGeometry/FaserActsLayerBuilder.h b/Tracking/Acts/FaserActsGeometry/FaserActsGeometry/FaserActsLayerBuilder.h index fbb91545..9ab2ae89 100644 --- a/Tracking/Acts/FaserActsGeometry/FaserActsGeometry/FaserActsLayerBuilder.h +++ b/Tracking/Acts/FaserActsGeometry/FaserActsGeometry/FaserActsLayerBuilder.h @@ -10,6 +10,7 @@ // ATHENA #include "TrackerReadoutGeometry/SiDetectorElement.h" +#include "GeoModelFaserUtilities/GeoModelExperiment.h" // ACTS #include "Acts/Geometry/GeometryContext.hpp" @@ -51,9 +52,12 @@ public: FaserActsDetectorElement::Subdetector subdetector = FaserActsDetectorElement::Subdetector::SCT; const TrackerDD::SCT_DetectorManager* mng; + const GeoVDetectorManager* vetomng; + const GeoVDetectorManager* triggermng; std::shared_ptr<const Acts::LayerCreator> layerCreator = nullptr; std::shared_ptr<std::vector<std::shared_ptr<const FaserActsDetectorElement>>> elementStore; std::shared_ptr<std::map<Identifier, Acts::GeometryIdentifier>> identifierMap; + std::pair<size_t, size_t> MaterialBins = {30, 30}; int station; int plane; }; @@ -66,7 +70,6 @@ public: = Acts::getDefaultLogger("GMLayBldr", Acts::Logging::INFO)) : m_logger(std::move(logger)) { - // std::cout << "GMLB construct" << std::endl; m_cfg = cfg; } @@ -105,13 +108,16 @@ public: buildLayers(const Acts::GeometryContext&, std::vector<std::shared_ptr<const Acts::Surface> >&, FaserActs::CuboidVolumeBuilder::LayerConfig&, FaserActs::CuboidVolumeBuilder::SurfaceConfig&); private: + + FaserActs::CuboidVolumeBuilder::VolumeConfig buildScintVolume(double , double , std::string ); + double m_ModuleWidth; double m_ModuleLength; /// configruation object Config m_cfg; - Acts::Vector3 m_worldDimensions = { 400.0_mm, 400.0_mm, 6000.0_mm }; - Acts::Vector3 m_worldCenter = {0.0, 0.0, 1276.0_mm}; - Acts::Vector3 m_trackerDimensions = { 400.0_mm, 400.0_mm, 1200.0_mm }; + Acts::Vector3 m_worldDimensions = { 400.0_mm, 400.0_mm, 8000.0_mm }; + Acts::Vector3 m_worldCenter = {0.0, 0.0, 0.0}; + Acts::Vector3 m_trackerDimensions = { 400.0_mm, 400.0_mm, 50.0_mm }; /// Private access to the logger const Acts::Logger& diff --git a/Tracking/Acts/FaserActsGeometry/python/ActsGeometryConfig.py b/Tracking/Acts/FaserActsGeometry/python/ActsGeometryConfig.py index a1c3c9f1..d81e6fd9 100644 --- a/Tracking/Acts/FaserActsGeometry/python/ActsGeometryConfig.py +++ b/Tracking/Acts/FaserActsGeometry/python/ActsGeometryConfig.py @@ -1,7 +1,7 @@ # Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator from AthenaConfiguration.ComponentFactory import CompFactory -FaserActsWriteTrackingGeometry,FaserActsTrackingGeometrySvc,FaserActsTrackingGeometryTool,FaserActsSurfaceMappingTool,FaserActsVolumeMappingTool,FaserActsObjWriterTool,FaserActsAlignmentCondAlg = CompFactory.getComps("FaserActsWriteTrackingGeometry","FaserActsTrackingGeometrySvc","FaserActsTrackingGeometryTool","FaserActsSurfaceMappingTool","FaserActsVolumeMappingTool","FaserActsObjWriterTool","FaserActsAlignmentCondAlg") +FaserActsExtrapolationTool,FaserActsMaterialTrackWriterSvc,FaserActsMaterialJsonWriterTool,FaserActsWriteTrackingGeometry,FaserActsTrackingGeometrySvc,FaserActsTrackingGeometryTool,FaserActsSurfaceMappingTool,FaserActsVolumeMappingTool,FaserActsObjWriterTool,FaserActsAlignmentCondAlg = CompFactory.getComps("FaserActsExtrapolationTool","FaserActsMaterialTrackWriterSvc","FaserActsMaterialJsonWriterTool","FaserActsWriteTrackingGeometry","FaserActsTrackingGeometrySvc","FaserActsTrackingGeometryTool","FaserActsSurfaceMappingTool","FaserActsVolumeMappingTool","FaserActsObjWriterTool","FaserActsAlignmentCondAlg") @@ -35,11 +35,10 @@ def ActsPropStepRootWriterSvcCfg(configFlags, def ActsTrackingGeometryToolCfg(configFlags, name = "ActsTrackingGeometryTool" ) : result = ComponentAccumulator() - from FaserActsGeometry.FaserActsWriteTrackingGeometryConfig import FaserActsTrackingGeometrySvcCfg - acc = FaserActsTrackingGeometrySvcCfg(configFlags) + acc = ActsTrackingGeometrySvcCfg(configFlags) result.merge(acc) -# Acts_ActsTrackingGeometryTool = CompFactory.ActsTrackingGeometryTool - actsTrackingGeometryTool = FaserActsTrackingGeometryTool("TrackingGeometryTool") + Acts_ActsTrackingGeometryTool = FaserActsTrackingGeometryTool + actsTrackingGeometryTool = Acts_ActsTrackingGeometryTool("TrackingGeometryTool") result.addPublicTool(actsTrackingGeometryTool) return result, actsTrackingGeometryTool @@ -69,7 +68,7 @@ def ActsAlignmentCondAlgCfg(configFlags, name = "ActsAlignmentCondAlg", **kwargs return result from MagFieldServices.MagFieldServicesConfig import MagneticFieldSvcCfg -def ActsExtrapolationToolCfg(configFlags, name="ActsExtrapolationTool", **kwargs) : +def ActsExtrapolationToolCfg(configFlags, name="FaserActsExtrapolationTool", **kwargs) : result=ComponentAccumulator() acc = MagneticFieldSvcCfg(configFlags) @@ -78,18 +77,18 @@ def ActsExtrapolationToolCfg(configFlags, name="ActsExtrapolationTool", **kwargs acc, actsTrackingGeometryTool = ActsTrackingGeometryToolCfg(configFlags) result.merge(acc) - Acts_ActsExtrapolationTool = CompFactory.ActsExtrapolationTool + Acts_ActsExtrapolationTool = CompFactory.FaserActsExtrapolationTool actsExtrapolationTool = Acts_ActsExtrapolationTool(name, **kwargs) result.addPublicTool(actsExtrapolationTool, primary=True) return result -def ActsMaterialTrackWriterSvcCfg(name="ActsMaterialTrackWriterSvc", +def ActsMaterialTrackWriterSvcCfg(name="FaserActsMaterialTrackWriterSvc", FilePath="MaterialTracks_mapping.root", TreeName="material-tracks") : result = ComponentAccumulator() - Acts_ActsMaterialTrackWriterSvc = CompFactory.ActsMaterialTrackWriterSvc + Acts_ActsMaterialTrackWriterSvc = CompFactory.FaserActsMaterialTrackWriterSvc ActsMaterialTrackWriterSvc = Acts_ActsMaterialTrackWriterSvc(name, FilePath=FilePath, TreeName=TreeName) @@ -115,10 +114,10 @@ def ActsSurfaceMappingToolCfg(configFlags, name = "FaserActsSurfaceMappingTool" result=ComponentAccumulator() acc, actsTrackingGeometryTool = ActsTrackingGeometryToolCfg(configFlags) - result.merge(acc) - + result.merge(acc) Acts_ActsSurfaceMappingTool = CompFactory.FaserActsSurfaceMappingTool ActsSurfaceMappingTool = Acts_ActsSurfaceMappingTool(name) + ActsSurfaceMappingTool.TrackingGeometryTool = actsTrackingGeometryTool from AthenaCommon.Constants import INFO ActsSurfaceMappingTool.OutputLevel = INFO @@ -134,6 +133,7 @@ def ActsVolumeMappingToolCfg(configFlags, name = "ActsVolumeMappingTool" ) : Acts_ActsVolumeMappingTool = CompFactory.FaserActsVolumeMappingTool FaserActsVolumeMappingTool = Acts_ActsVolumeMappingTool(name) + FaserActsVolumeMappingTool.TrackingGeometryTool = actsTrackingGeometryTool from AthenaCommon.Constants import INFO FaserActsVolumeMappingTool.OutputLevel = INFO @@ -144,7 +144,7 @@ def ActsVolumeMappingToolCfg(configFlags, name = "ActsVolumeMappingTool" ) : def ActsMaterialJsonWriterToolCfg(name= "ActsMaterialJsonWriterTool", **kwargs) : result=ComponentAccumulator() - Acts_ActsMaterialJsonWriterTool = CompFactory.ActsMaterialJsonWriterTool + Acts_ActsMaterialJsonWriterTool = CompFactory.FaserActsMaterialJsonWriterTool ActsMaterialJsonWriterTool = Acts_ActsMaterialJsonWriterTool(name, **kwargs) from AthenaCommon.Constants import INFO @@ -156,7 +156,7 @@ def ActsMaterialJsonWriterToolCfg(name= "ActsMaterialJsonWriterTool", **kwargs) def ActsObjWriterToolCfg(name= "ActsObjWriterTool", **kwargs) : result=ComponentAccumulator() - Acts_ActsObjWriterTool = CompFactory.ActsObjWriterTool + Acts_ActsObjWriterTool = FaserActsObjWriterTool ActsObjWriterTool = Acts_ActsObjWriterTool(name, **kwargs) from AthenaCommon.Constants import INFO diff --git a/Tracking/Acts/FaserActsGeometry/python/FaserActsExtrapolationConfig.py b/Tracking/Acts/FaserActsGeometry/python/FaserActsExtrapolationConfig.py index 8720049d..b16eb1df 100644 --- a/Tracking/Acts/FaserActsGeometry/python/FaserActsExtrapolationConfig.py +++ b/Tracking/Acts/FaserActsGeometry/python/FaserActsExtrapolationConfig.py @@ -10,10 +10,10 @@ from MagFieldServices.MagFieldServicesConfig import MagneticFieldSvcCfg FaserActsExtrapolationAlg,FaserActsExtrapolationTool,FaserActsTrackingGeometrySvc,FaserActsTrackingGeometryTool,FaserActsAlignmentCondAlg = CompFactory.getComps("FaserActsExtrapolationAlg","FaserActsExtrapolationTool","FaserActsTrackingGeometrySvc","FaserActsTrackingGeometryTool","FaserActsAlignmentCondAlg") -def FaserActsTrackingGeometrySvcCfg(flags, **kwargs): - acc = ComponentAccumulator() - acc.addService(FaserActsTrackingGeometrySvc(name = "FaserActsTrackingGeometrySvc", **kwargs)) - return acc +#def FaserActsTrackingGeometrySvcCfg(flags, **kwargs): +# acc = ComponentAccumulator() +# acc.addService(FaserActsTrackingGeometrySvc(name = "FaserActsTrackingGeometrySvc", **kwargs)) +# return acc def FaserActsAlignmentCondAlgCfg(flags, **kwargs): @@ -29,11 +29,14 @@ def FaserActsExtrapolationAlgBasicCfg(flags, **kwargs): actsExtrapolationTool=FaserActsExtrapolationTool("FaserActsExtrapolationTool") actsExtrapolationTool.FieldMode="FASER" #actsExtrapolationTool.FieldMode="Constant" - actsExtrapolationTool.TrackingGeometryTool=FaserActsTrackingGeometryTool("TrackingGeometryTool") + from FaserActsGeometry.ActsGeometryConfig import ActsTrackingGeometryToolCfg + result, actsTrackingGeometryTool = ActsTrackingGeometryToolCfg(flags) + actsExtrapolationTool.TrackingGeometryTool=actsTrackingGeometryTool actsExtrapolationAlg.ExtrapolationTool=actsExtrapolationTool actsExtrapolationAlg.EtaRange=[1, 10] actsExtrapolationAlg.PtRange=[0.1, 100] actsExtrapolationAlg.NParticlesPerEvent=int(1e4) + acc.merge(result) acc.addEventAlgo(actsExtrapolationAlg) return acc @@ -51,7 +54,7 @@ def FaserActsExtrapolationAlgCfg(flags, **kwargs): # Initialize field service acc.merge(MagneticFieldSvcCfg(flags)) - acc.merge(FaserActsTrackingGeometrySvcCfg(flags, **kwargs)) + #acc.merge(FaserActsTrackingGeometrySvcCfg(flags, **kwargs)) #acc.merge(FaserActsAlignmentCondAlgCfg(flags)) acc.merge(FaserActsExtrapolationAlgBasicCfg(flags, **kwargs)) acc.merge(FaserActsExtrapolationAlg_OutputCfg(flags)) diff --git a/Tracking/Acts/FaserActsGeometry/python/FaserActsMaterialMapping_jobOptions.py b/Tracking/Acts/FaserActsGeometry/python/FaserActsMaterialMapping_jobOptions.py index 095edc76..99e81e8c 100644 --- a/Tracking/Acts/FaserActsGeometry/python/FaserActsMaterialMapping_jobOptions.py +++ b/Tracking/Acts/FaserActsGeometry/python/FaserActsMaterialMapping_jobOptions.py @@ -76,7 +76,7 @@ if "__main__" == __name__: cfg = MainServicesCfg(ConfigFlags) cfg.merge(FaserGeometryCfg(ConfigFlags)) - cfg.merge(ActsMaterialTrackWriterSvcCfg("ActsMaterialTrackWriterSvc", + cfg.merge(ActsMaterialTrackWriterSvcCfg("FaserActsMaterialTrackWriterSvc", "MaterialTracks_mapping.root")) cfg.merge(PoolReadCfg(ConfigFlags)) diff --git a/Tracking/Acts/FaserActsGeometry/python/FaserActsMaterialValidationAlg.py b/Tracking/Acts/FaserActsGeometry/python/FaserActsMaterialValidationAlg.py new file mode 100644 index 00000000..640d9ee0 --- /dev/null +++ b/Tracking/Acts/FaserActsGeometry/python/FaserActsMaterialValidationAlg.py @@ -0,0 +1,66 @@ +""" +This job options file will run an example extrapolation using the +Acts tracking geometry and the Acts extrapolation toolchain. + +Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration +""" + +# start from scratch with component accumulator + + +from FaserActsGeometry.ActsGeometryConfig import FaserActsExtrapolationToolCfg +from FaserActsGeometry.FaserActsExtrapolationConfig import FaserActsExtrapolationAlgCfg + + +if "__main__" == __name__: + from AthenaCommon.Configurable import Configurable + from AthenaCommon.Logging import log + from AthenaCommon.Constants import INFO + from CalypsoConfiguration.AllConfigFlags import ConfigFlags + from CalypsoConfiguration.MainServicesConfig import MainServicesCfg + from FaserActsGeometry.ActsGeometryConfig import FaserActsMaterialTrackWriterSvcCfg + + Configurable.configurableRun3Behavior = True + + ## Just enable ID for the moment. + ConfigFlags.Input.isMC = True + ConfigFlags.GeoModel.FaserVersion = "FASER-01" + ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-01" + ConfigFlags.TrackingGeometry.MaterialSource = "material-maps.json" + #ConfigFlags.Concurrency.NumThreads = 10 + #ConfigFlags.Concurrency.NumConcurrentEvents = 10 + + ConfigFlags.lock() + ConfigFlags.dump() + + cfg = MainServicesCfg(ConfigFlags) + +# alignCondAlgCfg = ActsAlignmentCondAlgCfg(ConfigFlags) +# cfg.merge(alignCondAlgCfg) + cfg.merge(FaserActsMaterialTrackWriterSvcCfg(ConfigFlags, "FaserActsMaterialTrackWriterSvc", "MaterialTracks_mapped.root")) + + print('DEF WRITER : ') + extrapol = FaserActsExtrapolationToolCfg(ConfigFlags, + InteractionMultiScatering = True, + InteractionEloss = True, + InteractionRecord = True) + cfg.merge(extrapol) + + alg = FaserActsExtrapolationAlgCfg(ConfigFlags, + OutputLevel=INFO, + NParticlesPerEvent=int(1e3), + EtaRange=[4.8, 100], + PtRange=[999, 1001], + XRange=[-150, 150], + YRange=[-150, 150], + WriteMaterialTracks = True, + ExtrapolationTool=extrapol.getPrimary()) + + cfg.merge(alg) + + #tgSvc = cfg.getService("FaserActsTrackingGeometrySvc") + #cfg.printConfig() + + log.info("CONFIG DONE") + + cfg.run(100) diff --git a/Tracking/Acts/FaserActsGeometry/python/FaserActsWriteTrackingGeometryConfig.py b/Tracking/Acts/FaserActsGeometry/python/FaserActsWriteTrackingGeometryConfig.py index 2aef1dac..c760be88 100644 --- a/Tracking/Acts/FaserActsGeometry/python/FaserActsWriteTrackingGeometryConfig.py +++ b/Tracking/Acts/FaserActsGeometry/python/FaserActsWriteTrackingGeometryConfig.py @@ -6,12 +6,12 @@ from AthenaConfiguration.ComponentFactory import CompFactory #from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg from FaserSCT_GeoModel.FaserSCT_GeoModelConfig import FaserSCT_GeometryCfg -FaserActsWriteTrackingGeometry,FaserActsTrackingGeometrySvc,FaserActsTrackingGeometryTool,FaserActsObjWriterTool,FaserActsAlignmentCondAlg = CompFactory.getComps("FaserActsWriteTrackingGeometry","FaserActsTrackingGeometrySvc","FaserActsTrackingGeometryTool","FaserActsObjWriterTool","FaserActsAlignmentCondAlg") +FaserActsMaterialJsonWriterTool,FaserActsWriteTrackingGeometry,FaserActsTrackingGeometrySvc,FaserActsTrackingGeometryTool,FaserActsObjWriterTool,FaserActsAlignmentCondAlg = CompFactory.getComps("FaserActsMaterialJsonWriterTool","FaserActsWriteTrackingGeometry","FaserActsTrackingGeometrySvc","FaserActsTrackingGeometryTool","FaserActsObjWriterTool","FaserActsAlignmentCondAlg") -def FaserActsTrackingGeometrySvcCfg(flags, **kwargs): - acc = ComponentAccumulator() - acc.addService(FaserActsTrackingGeometrySvc(name = "FaserActsTrackingGeometrySvc", **kwargs)) - return acc +#def FaserActsTrackingGeometrySvcCfg(flags, **kwargs): +# acc = ComponentAccumulator() +# acc.addService(FaserActsTrackingGeometrySvc(name = "FaserActsTrackingGeometrySvc", **kwargs)) +# return acc def FaserActsAlignmentCondAlgCfg(flags, **kwargs): @@ -24,8 +24,15 @@ def FaserActsAlignmentCondAlgCfg(flags, **kwargs): def FaserActsWriteTrackingGeometryBasicCfg(flags, **kwargs): acc = ComponentAccumulator() actsWriteTrackingGeometry=FaserActsWriteTrackingGeometry(**kwargs) - actsWriteTrackingGeometry.TrackingGeometryTool=FaserActsTrackingGeometryTool("TrackingGeometryTool") - actsWriteTrackingGeometry.ObjWriterTool=FaserActsObjWriterTool("FaserActsObjWriterTool",OutputDirectory="./", SubDetectors=["Station_1","Station_2","Station_3"]) + from FaserActsGeometry.ActsGeometryConfig import ActsTrackingGeometryToolCfg + result, actsTrackingGeometryTool = ActsTrackingGeometryToolCfg(flags) + acc.merge(result) + actsWriteTrackingGeometry.TrackingGeometryTool=actsTrackingGeometryTool + actsWriteTrackingGeometry.ObjWriterTool=FaserActsObjWriterTool("FaserActsObjWriterTool",OutputDirectory="./", SubDetectors=["Station_0","Station_1","Station_2","Station_3"]) + actsWriteTrackingGeometry.MaterialJsonWriterTool= FaserActsMaterialJsonWriterTool(OutputFile = "geometry-maps.json", + processSensitives = False, + processnonmaterial = True) + acc.addEventAlgo(actsWriteTrackingGeometry) return acc @@ -38,7 +45,7 @@ def FaserActsWriteTrackingGeometry_OutputCfg(flags): def FaserActsWriteTrackingGeometryCfg(flags, **kwargs): #acc = FaserGeometryCfg(flags) acc = FaserSCT_GeometryCfg(flags) - acc.merge(FaserActsTrackingGeometrySvcCfg(flags, **kwargs)) + #acc.merge(FaserActsTrackingGeometrySvcCfg(flags, **kwargs)) #acc.merge(FaserActsAlignmentCondAlgCfg(flags)) acc.merge(FaserActsWriteTrackingGeometryBasicCfg(flags, **kwargs)) acc.merge(FaserActsWriteTrackingGeometry_OutputCfg(flags)) diff --git a/Tracking/Acts/FaserActsGeometry/src/CuboidVolumeBuilder.cxx b/Tracking/Acts/FaserActsGeometry/src/CuboidVolumeBuilder.cxx index a8a8bd2d..df96a311 100644 --- a/Tracking/Acts/FaserActsGeometry/src/CuboidVolumeBuilder.cxx +++ b/Tracking/Acts/FaserActsGeometry/src/CuboidVolumeBuilder.cxx @@ -75,12 +75,11 @@ std::shared_ptr<const Acts::Layer> CuboidVolumeBuilder::buildLayer( if ( !cfg.surfaces.empty() ) { return layerCreator.planeLayer(gctx, cfg.surfaces, cfg.binsX, cfg.binsY, - Acts::BinningValue::binZ, std::nullopt, trafo); -// std::move(ap)); + Acts::BinningValue::binZ, std::nullopt, trafo, + std::move(ap)); } else { return layerCreator.planeLayer(gctx, {cfg.surface}, cfg.binsX, cfg.binsY, - Acts::BinningValue::binZ, std::nullopt, trafo); -// std::move(ap)); + Acts::BinningValue::binZ, std::nullopt, trafo, std::move(ap)); } } @@ -124,21 +123,23 @@ std::shared_ptr<Acts::TrackingVolume> CuboidVolumeBuilder::buildVolume( SurfaceConfig sCfg; sCfg.position = cfg.position; // Rotation of the surfaces: +pi/2 around axis y - Acts::Vector3 xPos(0., 0., 1.); - Acts::Vector3 yPos(0., 1., 0.); - Acts::Vector3 zPos(-1., 0., 0.); - sCfg.rotation.col(0) = xPos; - sCfg.rotation.col(1) = yPos; - sCfg.rotation.col(2) = zPos; + //Acts::Vector3 xPos(0., 0., 1.); + //Acts::Vector3 yPos(0., 1., 0.); + //Acts::Vector3 zPos(-1., 0., 0.); + //sCfg.rotation.col(0) = xPos; + //sCfg.rotation.col(1) = yPos; + //sCfg.rotation.col(2) = zPos; // Bounds sCfg.rBounds = std::make_shared<const Acts::RectangleBounds>( - Acts::RectangleBounds(cfg.length.y() * 0.5, cfg.length.z() * 0.5)); + Acts::RectangleBounds(cfg.length.x() * 0.5, cfg.length.y() * 0.5)); LayerConfig lCfg; lCfg.surfaceCfg = sCfg; cfg.layerCfg.push_back(lCfg); } + else{ + } // Gather the layers Acts::LayerVector layVec; diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsDetectorElement.cxx b/Tracking/Acts/FaserActsGeometry/src/FaserActsDetectorElement.cxx index 58ad59dc..55453348 100644 --- a/Tracking/Acts/FaserActsGeometry/src/FaserActsDetectorElement.cxx +++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsDetectorElement.cxx @@ -60,7 +60,6 @@ FaserActsDetectorElement::FaserActsDetectorElement(const TrackerDD::SiDetectorEl } else if (boundsType == Trk::SurfaceBounds::Trapezoid) { std::string shapeString = detElem->getMaterialGeom()->getLogVol()->getShape()->type(); - //std::cout << __FUNCTION__ << "tapezoid, GeoLogVol -> shape says: " << shapeString << std::endl; const TrackerDD::SiDetectorDesign &design = detElem->design(); diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsExtrapolationAlg.cxx b/Tracking/Acts/FaserActsGeometry/src/FaserActsExtrapolationAlg.cxx index 429b5cb7..6fd900f6 100755 --- a/Tracking/Acts/FaserActsGeometry/src/FaserActsExtrapolationAlg.cxx +++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsExtrapolationAlg.cxx @@ -18,14 +18,14 @@ #include "Acts/Definitions/Algebra.hpp" #include "Acts/Utilities/Logger.hpp" #include "Acts/Surfaces/PerigeeSurface.hpp" +#include "ActsInterop/Logger.h" // PACKAGE +#include "ActsGeometryInterfaces/IActsMaterialTrackWriterSvc.h" +#include "FaserActsGeometry/FaserActsGeometryContext.h" #include "FaserActsGeometryInterfaces/IFaserActsExtrapolationTool.h" #include "FaserActsGeometryInterfaces/IFaserActsTrackingGeometryTool.h" -#include "ActsInterop/Logger.h" -#include "FaserActsGeometry/FaserActsGeometryContext.h" #include "FaserActsGeometryInterfaces/IFaserActsPropStepRootWriterSvc.h" -//#include "FaserActsGeometry/IFaserActsMaterialTrackWriterSvc.h" // OTHER #include "CLHEP/Random/RandomEngine.h" @@ -40,8 +40,8 @@ FaserActsExtrapolationAlg::FaserActsExtrapolationAlg(const std::string& name, ISvcLocator* pSvcLocator) : AthReentrantAlgorithm(name, pSvcLocator), m_propStepWriterSvc("FaserActsPropStepRootWriterSvc", name), - m_rndmGenSvc("AthRNGSvc", name)//, - //m_materialTrackWriterSvc("FaserActsMaterialTrackWriterSvc", name) + m_rndmGenSvc("AthRNGSvc", name), + m_materialTrackWriterSvc("FaserActsMaterialTrackWriterSvc", name) { } @@ -53,9 +53,9 @@ StatusCode FaserActsExtrapolationAlg::initialize() { ATH_CHECK( m_extrapolationTool.retrieve() ); ATH_CHECK( m_propStepWriterSvc.retrieve() ); - //if (m_writeMaterialTracks) { - //ATH_CHECK( m_materialTrackWriterSvc.retrieve() ); - //} + if (m_writeMaterialTracks) { + ATH_CHECK( m_materialTrackWriterSvc.retrieve() ); + } m_objOut = std::make_unique<std::ofstream>("steps.obj"); @@ -75,14 +75,23 @@ StatusCode FaserActsExtrapolationAlg::execute(const EventContext& ctx) const for (size_t i = 0; i < m_nParticlePerEvent; i++) { - double d0 = 0; - double z0 = 0; +// double d0 = 0; +// double z0 = 0; double phi = rngEngine->flat() * 2*M_PI - M_PI; std::vector<double> etaRange = m_etaRange; + std::vector<double> xRange = m_xRange; + std::vector<double> yRange = m_yRange; double etaMin = etaRange.at(0); double etaMax = etaRange.at(1); double eta = rngEngine->flat() * std::abs(etaMax - etaMin) + etaMin; + double xMin = xRange.at(0); + double xMax = xRange.at(1); + double yMin = yRange.at(0); + double yMax = yRange.at(1); + double x0 = rngEngine->flat() * std::abs(xMax - xMin) + xMin; + double y0 = rngEngine->flat() * std::abs(yMax - yMin) + yMin; + std::vector<double> ptRange = m_ptRange; double ptMin = ptRange.at(0) * 1_GeV; double ptMax = ptRange.at(1) * 1_GeV; @@ -104,7 +113,7 @@ StatusCode FaserActsExtrapolationAlg::execute(const EventContext& ctx) const double t = 0; Acts::BoundVector pars; - pars << d0, z0, phi, theta, qop, t; + pars << x0, y0, phi, theta, qop, t; std::optional<Acts::BoundSymMatrix> cov = std::nullopt; if(charge != 0.) { @@ -116,6 +125,13 @@ StatusCode FaserActsExtrapolationAlg::execute(const EventContext& ctx) const auto output = m_extrapolationTool->propagationSteps(ctx, startParameters); m_propStepWriterSvc->write(output.first); writeStepsObj(output.first); + if(m_writeMaterialTracks){ + Acts::RecordedMaterialTrack track; + track.first.first = Acts::Vector3(0,0,0); + track.first.second = momentum; + track.second = std::move(output.second); + m_materialTrackWriterSvc->write(track); + } } diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsExtrapolationAlg.h b/Tracking/Acts/FaserActsGeometry/src/FaserActsExtrapolationAlg.h index 96bbda65..75887540 100755 --- a/Tracking/Acts/FaserActsGeometry/src/FaserActsExtrapolationAlg.h +++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsExtrapolationAlg.h @@ -35,7 +35,7 @@ namespace Acts { } } -//class IFaserActsMaterialTrackWriterSvc; +class IActsMaterialTrackWriterSvc; class EventContext; class IAthRNGSvc; @@ -56,13 +56,15 @@ private: ToolHandle<IFaserActsExtrapolationTool> m_extrapolationTool{this, "ExtrapolationTool", "FaserActsExtrapolationTool"}; // poor-mans Particle Gun is included here right now + Gaudi::Property<std::vector<double>> m_xRange{this, "XRange", {-150, 150}, ""}; + Gaudi::Property<std::vector<double>> m_yRange{this, "YRange", {-150, 150}, ""}; Gaudi::Property<std::vector<double>> m_etaRange{this, "EtaRange", {4, 10}, ""}; Gaudi::Property<std::vector<double>> m_ptRange{this, "PtRange", {0.1, 1000}, ""}; Gaudi::Property<size_t> m_nParticlePerEvent{this, "NParticlesPerEvent", 1, "The number of particles per event"}; // this does not work right now - //Gaudi::Property<bool> m_writeMaterialTracks{this, "WriteMaterialTracks", false, ""}; - //ServiceHandle<IFaserActsMaterialTrackWriterSvc> m_materialTrackWriterSvc; + Gaudi::Property<bool> m_writeMaterialTracks{this, "WriteMaterialTracks", false, ""}; + ServiceHandle<IActsMaterialTrackWriterSvc> m_materialTrackWriterSvc; mutable std::mutex m_writeMutex{}; mutable std::unique_ptr<std::ofstream> m_objOut; diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsExtrapolationTool.cxx b/Tracking/Acts/FaserActsGeometry/src/FaserActsExtrapolationTool.cxx index 77485eb7..0e92508b 100644 --- a/Tracking/Acts/FaserActsGeometry/src/FaserActsExtrapolationTool.cxx +++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsExtrapolationTool.cxx @@ -228,7 +228,7 @@ FaserActsExtrapolationTool::propagate(const EventContext& ctx, auto parameters = boost::apply_visitor([&](const auto& propagator) -> std::unique_ptr<const Acts::CurvilinearTrackParameters> { auto result = propagator.propagate(startParameters, options); if (!result.ok()) { - ATH_MSG_ERROR("Got error during propagation:" << result.error() + ATH_MSG_WARNING("Got error during propagation:" << result.error() << ". Returning empty parameters."); return nullptr; } diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsJsonGeometryConverter.cxx b/Tracking/Acts/FaserActsGeometry/src/FaserActsJsonGeometryConverter.cxx deleted file mode 100644 index c707a004..00000000 --- a/Tracking/Acts/FaserActsGeometry/src/FaserActsJsonGeometryConverter.cxx +++ /dev/null @@ -1,1071 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2017-2019 CERN for the benefit of the Acts project -// -// This Source Code Form is subject to the terms of the Mozilla Public -// License, v. 2.0. If a copy of the MPL was not distributed with this -// file, You can obtain one at http://mozilla.org/MPL/2.0/. - -#include "FaserActsGeometry/FaserActsJsonGeometryConverter.h" - -#include "Acts/Geometry/ApproachDescriptor.hpp" -#include "Acts/Geometry/CuboidVolumeBounds.hpp" -#include "Acts/Geometry/CutoutCylinderVolumeBounds.hpp" -#include "Acts/Geometry/CylinderVolumeBounds.hpp" -#include "Acts/Geometry/GeometryIdentifier.hpp" -#include "Acts/Geometry/TrackingVolume.hpp" -#include "Acts/Material/BinnedSurfaceMaterial.hpp" -#include "Acts/Material/HomogeneousSurfaceMaterial.hpp" -#include "Acts/Material/HomogeneousVolumeMaterial.hpp" -#include "Acts/Material/InterpolatedMaterialMap.hpp" -#include "Acts/Material/MaterialGridHelper.hpp" -#include "Acts/Material/ProtoSurfaceMaterial.hpp" -#include "Acts/Material/ProtoVolumeMaterial.hpp" -#include "Acts/Surfaces/RectangleBounds.hpp" -#include "Acts/Surfaces/SurfaceArray.hpp" -#include "Acts/Utilities/BinUtility.hpp" -#include "Acts/Utilities/BinningType.hpp" -#include <Acts/Surfaces/AnnulusBounds.hpp> -#include <Acts/Surfaces/CylinderBounds.hpp> -#include <Acts/Surfaces/RadialBounds.hpp> -#include <Acts/Surfaces/SurfaceBounds.hpp> -#include <Acts/Definitions/Algebra.hpp> - - -#include <cstdio> -#include <fstream> -#include <iostream> -#include <map> -#include <sstream> -#include <stdexcept> -#include <string> - -#include <boost/algorithm/string.hpp> -#include <boost/algorithm/string/finder.hpp> -#include <boost/algorithm/string/iter_find.hpp> - -namespace { - -using json = nlohmann::json; - -// helper functions to encode/decode indefinite material -// -// encoded either as `null` for vacuum or to an array of material parameters - -json encodeMaterial(const Acts::Material& material) { - if (!material) { - return nullptr; - } - json encoded = json::array(); - for (unsigned i = 0; i < material.parameters().size(); ++i) { - encoded.push_back(material.parameters()[i]); - } - return encoded; -} - -Acts::Material decodeMaterial(const json& encoded) { - if (encoded.is_null()) { - return {}; - } - Acts::Material::ParametersVector params = - Acts::Material::ParametersVector::Zero(); - for (auto i = params.size(); 0 < i--;) { - // .at(...) ensures bound checks - params[i] = encoded.at(i); - } - return Acts::Material(params); -} - -// helper functions to encode/decode concrete material slabs -// -// encoded as an object w/ two entries: `material` and `thickness` - -json encodeMaterialSlab(const Acts::MaterialSlab& slab) { - return { - {"material", encodeMaterial(slab.material())}, - {"thickness", slab.thickness()}, - }; -} - -Acts::MaterialSlab decodeMaterialSlab(const json& encoded) { - return Acts::MaterialSlab(decodeMaterial(encoded.at("material")), - encoded.at("thickness").get<float>()); -} - -} // namespace - -Acts::FaserActsJsonGeometryConverter::FaserActsJsonGeometryConverter( - const Acts::FaserActsJsonGeometryConverter::Config& cfg) - : m_cfg(std::move(cfg)) { - // Validate the configuration - if (!m_cfg.logger) { - throw std::invalid_argument("Missing logger"); - } -} - -std::pair<std::map<Acts::GeometryIdentifier, - std::shared_ptr<const Acts::ISurfaceMaterial>>, - std::map<Acts::GeometryIdentifier, - std::shared_ptr<const Acts::IVolumeMaterial>>> -Acts::FaserActsJsonGeometryConverter::jsonToMaterialMaps(const json& materialmaps) { - auto& j = materialmaps; - // The return maps - std::pair<SurfaceMaterialMap, VolumeMaterialMap> maps; - ACTS_VERBOSE("j2a: Reading material maps from json file."); - ACTS_VERBOSE("j2a: Found entries for " << j.count(m_cfg.volkey) - << " volume(s)."); - // Structured binding - for (auto& [key, value] : j.items()) { - // Check if this the volume key - if (key == m_cfg.volkey) { - // Get the volume json - auto volj = value; - for (auto& [vkey, vvalue] : volj.items()) { - // Create the volume id - int vid = std::stoi(vkey); - Acts::GeometryIdentifier volumeID; - volumeID.setVolume(vid); - ACTS_VERBOSE("j2a: -> Found Volume " << vid); - // Loop through the information in the volume - for (auto& [vckey, vcvalue] : vvalue.items()) { - if (vckey == m_cfg.boukey and m_cfg.processBoundaries and - not vcvalue.empty()) { - ACTS_VERBOSE("j2a: --> BoundarySurface(s) to be parsed"); - for (auto& [bkey, bvalue] : vcvalue.items()) { - // Create the boundary id - int bid = std::stoi(bkey); - Acts::GeometryIdentifier boundaryID(volumeID); - boundaryID.setBoundary(bid); - ACTS_VERBOSE("j2a: ---> Found boundary surface " << bid); - if (bvalue[m_cfg.mapkey] == true) { - auto boumat = jsonToSurfaceMaterial(bvalue); - maps.first[boundaryID] = - std::shared_ptr<const ISurfaceMaterial>(boumat); - } - } - } else if (vckey == m_cfg.laykey) { - ACTS_VERBOSE("j2a: --> Layer(s) to be parsed"); - // Loop over layers and repeat - auto layj = vcvalue; - for (auto& [lkey, lvalue] : layj.items()) { - // Create the layer id - int lid = std::stoi(lkey); - Acts::GeometryIdentifier layerID(volumeID); - layerID.setLayer(lid); - ACTS_VERBOSE("j2a: ---> Found Layer " << lid); - // Finally loop over layer components - for (auto& [lckey, lcvalue] : lvalue.items()) { - if (lckey == m_cfg.repkey and m_cfg.processRepresenting and - not lcvalue.empty()) { - ACTS_VERBOSE("j2a: ----> Found representing surface"); - if (lcvalue[m_cfg.mapkey] == true) { - auto repmat = jsonToSurfaceMaterial(lcvalue); - maps.first[layerID] = - std::shared_ptr<const Acts::ISurfaceMaterial>(repmat); - } - } else if (lckey == m_cfg.appkey and m_cfg.processApproaches and - not lcvalue.empty()) { - ACTS_VERBOSE("j2a: ----> Found approach surface(s)"); - // Loop over approach surfaces - for (auto& [askey, asvalue] : lcvalue.items()) { - // Create the layer id, todo set to max value - int aid = (askey == "*") ? 0 : std::stoi(askey); - Acts::GeometryIdentifier approachID(layerID); - approachID.setApproach(aid); - ACTS_VERBOSE("j2a: -----> Approach surface " << askey); - if (asvalue[m_cfg.mapkey] == true) { - auto appmat = jsonToSurfaceMaterial(asvalue); - maps.first[approachID] = - std::shared_ptr<const Acts::ISurfaceMaterial>(appmat); - } - } - } else if (lckey == m_cfg.senkey and m_cfg.processSensitives and - not lcvalue.empty()) { - ACTS_VERBOSE("j2a: ----> Found sensitive surface(s)"); - // Loop over sensitive surfaces - for (auto& [sskey, ssvalue] : lcvalue.items()) { - // Create the layer id, todo set to max value - int sid = (sskey == "*") ? 0 : std::stoi(sskey); - Acts::GeometryIdentifier senisitiveID(layerID); - senisitiveID.setSensitive(sid); - ACTS_VERBOSE("j2a: -----> Sensitive surface " << sskey); - if (ssvalue[m_cfg.mapkey] == true) { - auto senmat = jsonToSurfaceMaterial(ssvalue); - maps.first[senisitiveID] = - std::shared_ptr<const Acts::ISurfaceMaterial>(senmat); - } - } - } - } - } - - } else if (m_cfg.processVolumes and vckey == m_cfg.matkey and - not vcvalue.empty()) { - ACTS_VERBOSE("--> VolumeMaterial to be parsed"); - if (vcvalue[m_cfg.mapkey] == true) { - auto intermat = jsonToVolumeMaterial(vcvalue); - maps.second[volumeID] = - std::shared_ptr<const Acts::IVolumeMaterial>(intermat); - } - } - } - } - } else if (key == m_cfg.geoversion) { - ACTS_VERBOSE("Detector version: " << m_cfg.geoversion); - } - } - - // Return the filled maps - return maps; -} - -/// Convert method -/// -json Acts::FaserActsJsonGeometryConverter::materialMapsToJson( - const DetectorMaterialMaps& maps) { - DetectorRep detRep; - // Collect all GeometryIdentifiers per VolumeID for the formatted output - for (auto& [key, value] : maps.first) { - geo_id_value vid = key.volume(); - auto volRep = detRep.volumes.find(vid); - if (volRep == detRep.volumes.end()) { - detRep.volumes.insert({vid, VolumeRep()}); - volRep = detRep.volumes.find(vid); - volRep->second.volumeID = key; - } - geo_id_value lid = key.layer(); - if (lid != 0) { - // we are on a layer, get the layer rep - auto layRep = volRep->second.layers.find(lid); - if (layRep == volRep->second.layers.end()) { - volRep->second.layers.insert({lid, LayerRep()}); - layRep = volRep->second.layers.find(lid); - layRep->second.layerID = key; - } - // now insert appropriately - geo_id_value sid = key.sensitive(); - geo_id_value aid = key.approach(); - if (sid != 0) { - layRep->second.sensitives.insert({sid, value.get()}); - } else if (aid != 0) { - layRep->second.approaches.insert({aid, value.get()}); - } else { - layRep->second.representing = value.get(); - } - - } else { - // not on a layer can only be a boundary surface - geo_id_value bid = key.boundary(); - volRep->second.boundaries.insert({bid, value.get()}); - } - } - for (auto& [key, value] : maps.second) { - // find the volume representation - geo_id_value vid = key.volume(); - auto volRep = detRep.volumes.find(vid); - if (volRep == detRep.volumes.end()) { - detRep.volumes.insert({vid, VolumeRep()}); - volRep = detRep.volumes.find(vid); - volRep->second.volumeID = key; - } - volRep->second.material = value.get(); - } - // convert the detector representation to json format - return detectorRepToJson(detRep); -} - -/// Create Json from a detector represenation -json Acts::FaserActsJsonGeometryConverter::detectorRepToJson(const DetectorRep& detRep) { - json detectorj; - ACTS_VERBOSE("a2j: Writing json from detector representation"); - ACTS_VERBOSE("a2j: Found entries for " << detRep.volumes.size() - << " volume(s)."); - - json volumesj; - for (auto& [key, value] : detRep.volumes) { - json volj; - ACTS_VERBOSE("a2j: -> Writing Volume: " << key); - volj[m_cfg.namekey] = value.volumeName; - std::ostringstream svolumeID; - svolumeID << value.volumeID; - volj[m_cfg.geometryidkey] = svolumeID.str(); - if (m_cfg.processVolumes && value.material) { - volj[m_cfg.matkey] = volumeMaterialToJson(*value.material); - } - // Write the layers - if (not value.layers.empty()) { - ACTS_VERBOSE("a2j: ---> Found " << value.layers.size() << " layer(s) "); - json layersj; - for (auto& [lkey, lvalue] : value.layers) { - ACTS_VERBOSE("a2j: ----> Convert layer " << lkey); - json layj; - std::ostringstream slayerID; - slayerID << lvalue.layerID; - layj[m_cfg.geometryidkey] = slayerID.str(); - // First check for approaches - if (not lvalue.approaches.empty() and m_cfg.processApproaches) { - ACTS_VERBOSE("a2j: -----> Found " << lvalue.approaches.size() - << " approach surface(s)"); - json approachesj; - for (auto& [akey, avalue] : lvalue.approaches) { - ACTS_VERBOSE("a2j: ------> Convert approach surface " << akey); - approachesj[std::to_string(akey)] = surfaceMaterialToJson(*avalue); - if (lvalue.approacheSurfaces.find(akey) != - lvalue.approacheSurfaces.end()) - addSurfaceToJson(approachesj[std::to_string(akey)], - lvalue.approacheSurfaces.at(akey)); - } - // Add to the layer json - layj[m_cfg.appkey] = approachesj; - } - // Then check for sensitive - if (not lvalue.sensitives.empty() and m_cfg.processSensitives) { - ACTS_VERBOSE("a2j: -----> Found " << lvalue.sensitives.size() - << " sensitive surface(s)"); - json sensitivesj; - for (auto& [skey, svalue] : lvalue.sensitives) { - ACTS_VERBOSE("a2j: ------> Convert sensitive surface " << skey); - sensitivesj[std::to_string(skey)] = surfaceMaterialToJson(*svalue); - if (lvalue.sensitiveSurfaces.find(skey) != - lvalue.sensitiveSurfaces.end()) - addSurfaceToJson(sensitivesj[std::to_string(skey)], - lvalue.sensitiveSurfaces.at(skey)); - } - // Add to the layer json - layj[m_cfg.senkey] = sensitivesj; - } - // Finally check for representing - if (lvalue.representing != nullptr and m_cfg.processRepresenting) { - ACTS_VERBOSE("a2j: ------> Convert representing surface "); - layj[m_cfg.repkey] = surfaceMaterialToJson(*lvalue.representing); - if (lvalue.representingSurface != nullptr) - addSurfaceToJson(layj[m_cfg.repkey], lvalue.representingSurface); - } - layersj[std::to_string(lkey)] = layj; - } - volj[m_cfg.laykey] = layersj; - } - // Write the boundary surfaces - if (not value.boundaries.empty()) { - ACTS_VERBOSE("a2j: ---> Found " << value.boundaries.size() - << " boundary/ies "); - json boundariesj; - for (auto& [bkey, bvalue] : value.boundaries) { - ACTS_VERBOSE("a2j: ----> Convert boundary " << bkey); - boundariesj[std::to_string(bkey)] = surfaceMaterialToJson(*bvalue); - if (value.boundarySurfaces.find(bkey) != value.boundarySurfaces.end()) - addSurfaceToJson(boundariesj[std::to_string(bkey)], - value.boundarySurfaces.at(bkey)); - } - volj[m_cfg.boukey] = boundariesj; - } - - volumesj[std::to_string(key)] = volj; - } - // Assign the volume json to the detector json - detectorj[m_cfg.volkey] = volumesj; - - return detectorj; -} - -/// Create the Surface Material -const Acts::ISurfaceMaterial* -Acts::FaserActsJsonGeometryConverter::jsonToSurfaceMaterial(const json& material) { - Acts::ISurfaceMaterial* sMaterial = nullptr; - // The bin utility for deescribing the data - Acts::BinUtility bUtility; - for (auto& [key, value] : material.items()) { - if (key == m_cfg.transfokeys and not value.empty()) { - bUtility = Acts::BinUtility(jsonToTransform(value)); - break; - } - } - // Convert the material - Acts::MaterialSlabMatrix mpMatrix; - // Structured binding - for (auto& [key, value] : material.items()) { - // Check json keys - if (key == m_cfg.bin0key and not value.empty()) { - bUtility += jsonToBinUtility(value); - } else if (key == m_cfg.bin1key and not value.empty()) { - bUtility += jsonToBinUtility(value); - } - if (key == m_cfg.datakey and not value.empty()) { - mpMatrix = jsonToMaterialMatrix(value); - } - } - - // We have protoMaterial - if (mpMatrix.empty()) { - sMaterial = new Acts::ProtoSurfaceMaterial(bUtility); - } else if (bUtility.bins() == 1) { - sMaterial = new Acts::HomogeneousSurfaceMaterial(mpMatrix[0][0]); - } else { - sMaterial = new Acts::BinnedSurfaceMaterial(bUtility, mpMatrix); - } - // return what you have - return sMaterial; -} - -/// Create the Volume Material -const Acts::IVolumeMaterial* Acts::FaserActsJsonGeometryConverter::jsonToVolumeMaterial( - const json& material) { - Acts::IVolumeMaterial* vMaterial = nullptr; - // The bin utility for deescribing the data - Acts::BinUtility bUtility; - for (auto& [key, value] : material.items()) { - if (key == m_cfg.transfokeys and not value.empty()) { - bUtility = Acts::BinUtility(jsonToTransform(value)); - break; - } - } - // Convert the material - std::vector<Material> mmat; - // Structured binding - for (auto& [key, value] : material.items()) { - // Check json keys - if (key == m_cfg.bin0key and not value.empty()) { - bUtility += jsonToBinUtility(value); - } else if (key == m_cfg.bin1key and not value.empty()) { - bUtility += jsonToBinUtility(value); - } else if (key == m_cfg.bin2key and not value.empty()) { - bUtility += jsonToBinUtility(value); - } - if (key == m_cfg.datakey and not value.empty()) { - for (const auto& bin : value) { - mmat.push_back(decodeMaterial(bin)); - } - } - } - - // We have protoMaterial - if (mmat.empty()) { - vMaterial = new Acts::ProtoVolumeMaterial(bUtility); - } else if (mmat.size() == 1) { - vMaterial = new Acts::HomogeneousVolumeMaterial(mmat[0]); - } else { - if (bUtility.dimensions() == 2) { - std::function<Acts::Vector2(Acts::Vector3)> transfoGlobalToLocal; - Acts::Grid2D grid = createGrid2D(bUtility, transfoGlobalToLocal); - - Acts::Grid2D::point_t min = grid.minPosition(); - Acts::Grid2D::point_t max = grid.maxPosition(); - Acts::Grid2D::index_t nBins = grid.numLocalBins(); - - Acts::EAxis axis1(min[0], max[0], nBins[0]); - Acts::EAxis axis2(min[1], max[1], nBins[1]); - - // Build the grid and fill it with data - MaterialGrid2D mGrid(std::make_tuple(axis1, axis2)); - - for (size_t bin = 0; bin < mmat.size(); bin++) { - mGrid.at(bin) = mmat[bin].parameters(); - } - MaterialMapper<MaterialGrid2D> matMap(transfoGlobalToLocal, mGrid); - vMaterial = - new Acts::InterpolatedMaterialMap<MaterialMapper<MaterialGrid2D>>( - std::move(matMap), bUtility); - } else if (bUtility.dimensions() == 3) { - std::function<Acts::Vector3(Acts::Vector3)> transfoGlobalToLocal; - Acts::Grid3D grid = createGrid3D(bUtility, transfoGlobalToLocal); - - Acts::Grid3D::point_t min = grid.minPosition(); - Acts::Grid3D::point_t max = grid.maxPosition(); - Acts::Grid3D::index_t nBins = grid.numLocalBins(); - - Acts::EAxis axis1(min[0], max[0], nBins[0]); - Acts::EAxis axis2(min[1], max[1], nBins[1]); - Acts::EAxis axis3(min[2], max[2], nBins[2]); - - // Build the grid and fill it with data - MaterialGrid3D mGrid(std::make_tuple(axis1, axis2, axis3)); - - for (size_t bin = 0; bin < mmat.size(); bin++) { - mGrid.at(bin) = mmat[bin].parameters(); - } - MaterialMapper<MaterialGrid3D> matMap(transfoGlobalToLocal, mGrid); - vMaterial = - new Acts::InterpolatedMaterialMap<MaterialMapper<MaterialGrid3D>>( - std::move(matMap), bUtility); - } - } - // return what you have - return vMaterial; -} - -json Acts::FaserActsJsonGeometryConverter::trackingGeometryToJson( - const Acts::TrackingGeometry& tGeometry) { - DetectorRep detRep; - convertToRep(detRep, *tGeometry.highestTrackingVolume()); - return detectorRepToJson(detRep); -} - -void Acts::FaserActsJsonGeometryConverter::convertToRep( - DetectorRep& detRep, const Acts::TrackingVolume& tVolume) { - // The writer reader volume representation - VolumeRep volRep; - volRep.volumeName = tVolume.volumeName(); - // there are confined volumes - if (tVolume.confinedVolumes() != nullptr) { - // get through the volumes - auto& volumes = tVolume.confinedVolumes()->arrayObjects(); - // loop over the volumes - for (auto& vol : volumes) { - // recursive call - convertToRep(detRep, *vol); - } - } - // there are dense volumes - if (m_cfg.processDenseVolumes && !tVolume.denseVolumes().empty()) { - // loop over the volumes - for (auto& vol : tVolume.denseVolumes()) { - // recursive call - convertToRep(detRep, *vol); - } - } - // Get the volume Id - Acts::GeometryIdentifier volumeID = tVolume.geometryId(); - geo_id_value vid = volumeID.volume(); - - // Write the material if there's one - if (tVolume.volumeMaterial() != nullptr) { - volRep.material = tVolume.volumeMaterial(); - } else if (m_cfg.processnonmaterial == true) { - Acts::BinUtility bUtility = DefaultBin(tVolume); - Acts::IVolumeMaterial* bMaterial = new Acts::ProtoVolumeMaterial(bUtility); - volRep.material = bMaterial; - } - // there are confied layers - if (tVolume.confinedLayers() != nullptr) { - // get the layers - auto& layers = tVolume.confinedLayers()->arrayObjects(); - // loop of the volumes - for (auto& lay : layers) { - auto layRep = convertToRep(*lay); - if (layRep) { - // it's a valid representation so let's go with it - Acts::GeometryIdentifier layerID = lay->geometryId(); - geo_id_value lid = layerID.layer(); - volRep.layers.insert({lid, std::move(layRep)}); - } - } - } - // Let's finally check the boundaries - for (auto& bsurf : tVolume.boundarySurfaces()) { - // the surface representation - auto& bssfRep = bsurf->surfaceRepresentation(); - if (bssfRep.surfaceMaterial() != nullptr) { - Acts::GeometryIdentifier boundaryID = bssfRep.geometryId(); - geo_id_value bid = boundaryID.boundary(); - // Ignore if the volumeID is not correct (i.e. shared boundary) - // if (boundaryID.value(Acts::GeometryIdentifier::volume_mask) == vid){ - volRep.boundaries[bid] = bssfRep.surfaceMaterial(); - volRep.boundarySurfaces[bid] = &bssfRep; - // } - } else if (m_cfg.processnonmaterial == true) { - // if no material suface exist add a default one for the mapping - // configuration - Acts::GeometryIdentifier boundaryID = bssfRep.geometryId(); - geo_id_value bid = boundaryID.boundary(); - Acts::BinUtility bUtility = DefaultBin(bssfRep); - Acts::ISurfaceMaterial* bMaterial = - new Acts::ProtoSurfaceMaterial(bUtility); - volRep.boundaries[bid] = bMaterial; - volRep.boundarySurfaces[bid] = &bssfRep; - } - } - // Write if it's good - if (volRep) { - volRep.volumeName = tVolume.volumeName(); - volRep.volumeID = volumeID; - detRep.volumes.insert({vid, std::move(volRep)}); - } - return; -} - -Acts::FaserActsJsonGeometryConverter::LayerRep Acts::FaserActsJsonGeometryConverter::convertToRep( - const Acts::Layer& tLayer) { - LayerRep layRep; - // fill layer ID information - layRep.layerID = tLayer.geometryId(); - if (m_cfg.processSensitives and tLayer.surfaceArray() != nullptr) { - for (auto& ssf : tLayer.surfaceArray()->surfaces()) { - if (ssf != nullptr && ssf->surfaceMaterial() != nullptr) { - Acts::GeometryIdentifier sensitiveID = ssf->geometryId(); - geo_id_value sid = sensitiveID.sensitive(); - layRep.sensitives.insert({sid, ssf->surfaceMaterial()}); - layRep.sensitiveSurfaces.insert({sid, ssf}); - } else if (m_cfg.processnonmaterial == true) { - // if no material suface exist add a default one for the mapping - // configuration - Acts::GeometryIdentifier sensitiveID = ssf->geometryId(); - geo_id_value sid = sensitiveID.sensitive(); - Acts::BinUtility sUtility = DefaultBin(*ssf); - Acts::ISurfaceMaterial* sMaterial = - new Acts::ProtoSurfaceMaterial(sUtility); - layRep.sensitives.insert({sid, sMaterial}); - layRep.sensitiveSurfaces.insert({sid, ssf}); - } - } - } - // the representing - if (!(tLayer.surfaceRepresentation().geometryId() == GeometryIdentifier())) { - if (tLayer.surfaceRepresentation().surfaceMaterial() != nullptr) { - layRep.representing = tLayer.surfaceRepresentation().surfaceMaterial(); - layRep.representingSurface = &tLayer.surfaceRepresentation(); - } else if (m_cfg.processnonmaterial == true) { - // if no material suface exist add a default one for the mapping - // configuration - Acts::BinUtility rUtility = DefaultBin(tLayer.surfaceRepresentation()); - Acts::ISurfaceMaterial* rMaterial = - new Acts::ProtoSurfaceMaterial(rUtility); - layRep.representing = rMaterial; - layRep.representingSurface = &tLayer.surfaceRepresentation(); - } - } - // the approach - if (tLayer.approachDescriptor() != nullptr) { - for (auto& asf : tLayer.approachDescriptor()->containedSurfaces()) { - // get the surface and check for material - if (asf->surfaceMaterial() != nullptr) { - Acts::GeometryIdentifier approachID = asf->geometryId(); - geo_id_value aid = approachID.approach(); - layRep.approaches.insert({aid, asf->surfaceMaterial()}); - layRep.approacheSurfaces.insert({aid, asf}); - } else if (m_cfg.processnonmaterial == true) { - // if no material suface exist add a default one for the mapping - // configuration - Acts::GeometryIdentifier approachID = asf->geometryId(); - geo_id_value aid = approachID.approach(); - Acts::BinUtility aUtility = DefaultBin(*asf); - Acts::ISurfaceMaterial* aMaterial = - new Acts::ProtoSurfaceMaterial(aUtility); - layRep.approaches.insert({aid, aMaterial}); - layRep.approacheSurfaces.insert({aid, asf}); - } - } - } - // return the layer representation - return layRep; -} - -json Acts::FaserActsJsonGeometryConverter::surfaceMaterialToJson( - const Acts::ISurfaceMaterial& sMaterial) { - json smj; - // A bin utility needs to be written - const Acts::BinUtility* bUtility = nullptr; - // Check if we have a proto material - auto psMaterial = dynamic_cast<const Acts::ProtoSurfaceMaterial*>(&sMaterial); - if (psMaterial != nullptr) { - // Type is proto material - smj[m_cfg.typekey] = "proto"; - // by default the protoMaterial is not used for mapping - smj[m_cfg.mapkey] = false; - bUtility = &(psMaterial->binUtility()); - } else { - // Now check if we have a homogeneous material - auto hsMaterial = - dynamic_cast<const Acts::HomogeneousSurfaceMaterial*>(&sMaterial); - if (hsMaterial != nullptr) { - // type is homogeneous - smj[m_cfg.typekey] = "homogeneous"; - smj[m_cfg.mapkey] = true; - if (m_cfg.writeData) { - smj[m_cfg.datakey] = json::array({ - json::array({ - encodeMaterialSlab(hsMaterial->materialSlab(0, 0)), - }), - }); - } - } else { - // Only option remaining: BinnedSurface material - auto bsMaterial = - dynamic_cast<const Acts::BinnedSurfaceMaterial*>(&sMaterial); - if (bsMaterial != nullptr) { - // type is binned - smj[m_cfg.typekey] = "binned"; - smj[m_cfg.mapkey] = true; - bUtility = &(bsMaterial->binUtility()); - // convert the data - // get the material matrix - if (m_cfg.writeData) { - json mmat = json::array(); - for (const auto& mpVector : bsMaterial->fullMaterial()) { - json mvec = json::array(); - for (const auto& mp : mpVector) { - mvec.push_back(encodeMaterialSlab(mp)); - } - mmat.push_back(std::move(mvec)); - } - smj[m_cfg.datakey] = std::move(mmat); - } - } - } - } - // add the bin utility - if (bUtility != nullptr && !bUtility->binningData().empty()) { - std::vector<std::string> binkeys = {m_cfg.bin0key, m_cfg.bin1key}; - // loop over dimensions and write - auto& binningData = bUtility->binningData(); - // loop over the dimensions - for (size_t ibin = 0; ibin < binningData.size(); ++ibin) { - json binj; - auto cbData = binningData[ibin]; - binj.push_back(Acts::binningValueNames[cbData.binvalue]); - if (cbData.option == Acts::closed) { - binj.push_back("closed"); - } else { - binj.push_back("open"); - } - binj.push_back(cbData.bins()); - // If protoMaterial has a non uniform binning (non default) then it is - // used by default in the mapping - if (smj[m_cfg.typekey] == "proto" && cbData.bins() > 1) - smj[m_cfg.mapkey] = true; - // If it's not a proto map, write min / max - if (smj[m_cfg.typekey] != "proto") { - std::pair<double, double> minMax = {cbData.min, cbData.max}; - binj.push_back(minMax); - } - smj[binkeys[ibin]] = binj; - } - std::vector<double> transfo; - Acts::Transform3 transfo_matrix = bUtility->transform(); - if (not transfo_matrix.isApprox(Acts::Transform3::Identity())) { - for (int i = 0; i < 4; i++) { - for (int j = 0; j < 4; j++) { - transfo.push_back(transfo_matrix(j, i)); - } - } - smj[m_cfg.transfokeys] = transfo; - } - } - return smj; -} - -json Acts::FaserActsJsonGeometryConverter::volumeMaterialToJson( - const Acts::IVolumeMaterial& vMaterial) { - json smj; - // A bin utility needs to be written - const Acts::BinUtility* bUtility = nullptr; - // Check if we have a proto material - auto pvMaterial = dynamic_cast<const Acts::ProtoVolumeMaterial*>(&vMaterial); - if (pvMaterial != nullptr) { - // Type is proto material - smj[m_cfg.typekey] = "proto"; - // by default the protoMaterial is not used for mapping - smj[m_cfg.mapkey] = false; - bUtility = &(pvMaterial->binUtility()); - } else { - // Now check if we have a homogeneous material - auto hvMaterial = - dynamic_cast<const Acts::HomogeneousVolumeMaterial*>(&vMaterial); - if (hvMaterial != nullptr) { - // type is homogeneous - smj[m_cfg.typekey] = "homogeneous"; - smj[m_cfg.mapkey] = true; - if (m_cfg.writeData) { - // array of encoded materials w/ one entry - smj[m_cfg.datakey] = json::array({ - encodeMaterial(hvMaterial->material({0, 0, 0})), - }); - } - } else { - // Only option remaining: material map - auto bvMaterial2D = dynamic_cast< - const Acts::InterpolatedMaterialMap<MaterialMapper<MaterialGrid2D>>*>( - &vMaterial); - // Now check if we have a 2D map - if (bvMaterial2D != nullptr) { - // type is binned - smj[m_cfg.typekey] = "interpolated2D"; - smj[m_cfg.mapkey] = true; - bUtility = &(bvMaterial2D->binUtility()); - // convert the data - if (m_cfg.writeData) { - json mmat = json::array(); - MaterialGrid2D grid = bvMaterial2D->getMapper().getGrid(); - for (size_t bin = 0; bin < grid.size(); bin++) { - mmat.push_back(encodeMaterial(grid.at(bin))); - } - smj[m_cfg.datakey] = std::move(mmat); - } - } else { - // Only option remaining: material map - auto bvMaterial3D = dynamic_cast<const Acts::InterpolatedMaterialMap< - MaterialMapper<MaterialGrid3D>>*>(&vMaterial); - // Now check if we have a 3D map - if (bvMaterial3D != nullptr) { - // type is binned - smj[m_cfg.typekey] = "interpolated3D"; - smj[m_cfg.mapkey] = true; - bUtility = &(bvMaterial3D->binUtility()); - // convert the data - if (m_cfg.writeData) { - json mmat = json::array(); - MaterialGrid3D grid = bvMaterial3D->getMapper().getGrid(); - for (size_t bin = 0; bin < grid.size(); bin++) { - mmat.push_back(encodeMaterial(grid.at(bin))); - } - smj[m_cfg.datakey] = std::move(mmat); - } - } - } - } - } - // add the bin utility - if (bUtility != nullptr && !bUtility->binningData().empty()) { - std::vector<std::string> binkeys = {m_cfg.bin0key, m_cfg.bin1key, - m_cfg.bin2key}; - // loop over dimensions and write - auto& binningData = bUtility->binningData(); - // loop over the dimensions - for (size_t ibin = 0; ibin < binningData.size(); ++ibin) { - json binj; - auto cbData = binningData[ibin]; - binj.push_back(Acts::binningValueNames[cbData.binvalue]); - if (cbData.option == Acts::closed) { - binj.push_back("closed"); - } else { - binj.push_back("open"); - } - binj.push_back(cbData.bins()); - // If protoMaterial has a non uniform binning (non default) then it is - // used by default in the mapping - if (smj[m_cfg.typekey] == "proto" && cbData.bins() > 1) - smj[m_cfg.mapkey] = true; - // If it's not a proto map, write min / max - if (smj[m_cfg.typekey] != "proto") { - std::pair<double, double> minMax = {cbData.min, cbData.max}; - binj.push_back(minMax); - } - smj[binkeys[ibin]] = binj; - } - std::vector<double> transfo; - Acts::Transform3 transfo_matrix = bUtility->transform(); - for (int i = 0; i < 4; i++) { - for (int j = 0; j < 4; j++) { - transfo.push_back(transfo_matrix(j, i)); - } - } - smj[m_cfg.transfokeys] = transfo; - } - return smj; -} - -void Acts::FaserActsJsonGeometryConverter::addSurfaceToJson(json& sjson, - const Surface* surface) { - // Get the ID of the surface (redundant but help readability) - std::ostringstream SurfaceID; - SurfaceID << surface->geometryId(); - sjson[m_cfg.surfacegeometryidkey] = SurfaceID.str(); - - // Cast the surface bound to both disk and cylinder - const Acts::SurfaceBounds& surfaceBounds = surface->bounds(); - auto sTransform = surface->transform(GeometryContext()); - - const Acts::RadialBounds* radialBounds = - dynamic_cast<const Acts::RadialBounds*>(&surfaceBounds); - const Acts::CylinderBounds* cylinderBounds = - dynamic_cast<const Acts::CylinderBounds*>(&surfaceBounds); - const Acts::AnnulusBounds* annulusBounds = - dynamic_cast<const Acts::AnnulusBounds*>(&surfaceBounds); - - if (radialBounds != nullptr) { - sjson[m_cfg.surfacetypekey] = "Disk"; - sjson[m_cfg.surfacepositionkey] = sTransform.translation().z(); - sjson[m_cfg.surfacerangekey] = {radialBounds->rMin(), radialBounds->rMax()}; - } - if (cylinderBounds != nullptr) { - sjson[m_cfg.surfacetypekey] = "Cylinder"; - sjson[m_cfg.surfacepositionkey] = cylinderBounds->get(CylinderBounds::eR); - sjson[m_cfg.surfacerangekey] = { - -1 * cylinderBounds->get(CylinderBounds::eHalfLengthZ), - cylinderBounds->get(CylinderBounds::eHalfLengthZ)}; - } - if (annulusBounds != nullptr) { - sjson[m_cfg.surfacetypekey] = "Annulus"; - sjson[m_cfg.surfacepositionkey] = sTransform.translation().z(); - sjson[m_cfg.surfacerangekey] = { - {annulusBounds->rMin(), annulusBounds->rMax()}, - {annulusBounds->phiMin(), annulusBounds->phiMax()}}; - } -} - -/// Create the Material Matrix -Acts::MaterialSlabMatrix Acts::FaserActsJsonGeometryConverter::jsonToMaterialMatrix( - const json& data) { - Acts::MaterialSlabMatrix mpMatrix; - // the input data must be array[array[object]] - for (auto& outer : data) { - Acts::MaterialSlabVector mpVector; - for (auto& inner : outer) { - mpVector.emplace_back(decodeMaterialSlab(inner)); - } - mpMatrix.push_back(std::move(mpVector)); - } - return mpMatrix; -} - -/// Create the BinUtility for this -Acts::BinUtility Acts::FaserActsJsonGeometryConverter::jsonToBinUtility( - const json& bin) { - if (bin.size() >= 3) { - // finding the iterator position to determine the binning value - auto bit = std::find(Acts::binningValueNames.begin(), - Acts::binningValueNames.end(), bin[0]); - size_t indx = std::distance(Acts::binningValueNames.begin(), bit); - Acts::BinningValue bval = Acts::BinningValue(indx); - Acts::BinningOption bopt = bin[1] == "open" ? Acts::open : Acts::closed; - unsigned int bins = bin[2]; - float min = 0; - float max = 0; - if (bin.size() >= 4 && bin[3].size() == 2) { - min = bin[3][0]; - max = bin[3][1]; - } - return Acts::BinUtility(bins, min, max, bopt, bval); - } - return Acts::BinUtility(); -} - -/// Create the local to global transform -Acts::Transform3 Acts::FaserActsJsonGeometryConverter::jsonToTransform( - const json& transfo) { - Transform3 transform; - int i = 0; - int j = 0; - for (auto& element : transfo) { - transform(j, i) = element; - j++; - if (j == 4) { - i++; - j = 0; - } - } - return transform; -} - -Acts::BinUtility Acts::FaserActsJsonGeometryConverter::DefaultBin( - const Acts::Surface& surface) { - Acts::BinUtility bUtility; - - const Acts::SurfaceBounds& surfaceBounds = surface.bounds(); - const Acts::RadialBounds* radialBounds = - dynamic_cast<const Acts::RadialBounds*>(&surfaceBounds); - const Acts::CylinderBounds* cylinderBounds = - dynamic_cast<const Acts::CylinderBounds*>(&surfaceBounds); - const Acts::AnnulusBounds* annulusBounds = - dynamic_cast<const Acts::AnnulusBounds*>(&surfaceBounds); - const Acts::RectangleBounds* rectangleBounds = - dynamic_cast<const Acts::RectangleBounds*>(&surfaceBounds); - - if (radialBounds != nullptr) { - bUtility += BinUtility( - 1, - radialBounds->get(RadialBounds::eAveragePhi) - - radialBounds->get(RadialBounds::eHalfPhiSector), - radialBounds->get(RadialBounds::eAveragePhi) + - radialBounds->get(RadialBounds::eHalfPhiSector), - (radialBounds->get(RadialBounds::eHalfPhiSector) - M_PI) < s_epsilon - ? Acts::closed - : Acts::open, - Acts::binPhi); - bUtility += BinUtility(1, radialBounds->rMin(), radialBounds->rMax(), - Acts::open, Acts::binR); - return bUtility; - } - if (cylinderBounds != nullptr) { - bUtility += BinUtility( - 1, - cylinderBounds->get(CylinderBounds::eAveragePhi) - - cylinderBounds->get(CylinderBounds::eHalfPhiSector), - cylinderBounds->get(CylinderBounds::eAveragePhi) + - cylinderBounds->get(CylinderBounds::eHalfPhiSector), - (cylinderBounds->get(CylinderBounds::eHalfPhiSector) - M_PI) < s_epsilon - ? Acts::closed - : Acts::open, - Acts::binPhi); - bUtility += - BinUtility(1, -1 * cylinderBounds->get(CylinderBounds::eHalfLengthZ), - cylinderBounds->get(CylinderBounds::eHalfLengthZ), - Acts::open, Acts::binZ); - return bUtility; - } - if (annulusBounds != nullptr) { - bUtility += BinUtility(1, annulusBounds->get(AnnulusBounds::eMinPhiRel), - annulusBounds->get(AnnulusBounds::eMaxPhiRel), - Acts::open, Acts::binPhi); - bUtility += BinUtility(1, annulusBounds->rMin(), annulusBounds->rMax(), - Acts::open, Acts::binR); - return bUtility; - } - if (rectangleBounds != nullptr) { - bUtility += BinUtility(1, rectangleBounds->get(RectangleBounds::eMinX), - rectangleBounds->get(RectangleBounds::eMaxX), - Acts::open, Acts::binX); - bUtility += BinUtility(1, rectangleBounds->get(RectangleBounds::eMinY), - rectangleBounds->get(RectangleBounds::eMaxY), - Acts::open, Acts::binY); - return bUtility; - } - ACTS_INFO( - "No corresponding bound found for the surface : " << surface.name()); - return bUtility; -} - -Acts::BinUtility Acts::FaserActsJsonGeometryConverter::DefaultBin( - const Acts::TrackingVolume& volume) { - Acts::BinUtility bUtility; - - auto cyBounds = - dynamic_cast<const CylinderVolumeBounds*>(&(volume.volumeBounds())); - auto cutcylBounds = - dynamic_cast<const CutoutCylinderVolumeBounds*>(&(volume.volumeBounds())); - auto cuBounds = - dynamic_cast<const CuboidVolumeBounds*>(&(volume.volumeBounds())); - - if (cyBounds != nullptr) { - bUtility += BinUtility(1, cyBounds->get(CylinderVolumeBounds::eMinR), - cyBounds->get(CylinderVolumeBounds::eMaxR), - Acts::open, Acts::binR); - bUtility += BinUtility( - 1, -cyBounds->get(CylinderVolumeBounds::eHalfPhiSector), - cyBounds->get(CylinderVolumeBounds::eHalfPhiSector), - (cyBounds->get(CylinderVolumeBounds::eHalfPhiSector) - M_PI) < s_epsilon - ? Acts::closed - : Acts::open, - Acts::binPhi); - bUtility += - BinUtility(1, -cyBounds->get(CylinderVolumeBounds::eHalfLengthZ), - cyBounds->get(CylinderVolumeBounds::eHalfLengthZ), - Acts::open, Acts::binZ); - return bUtility; - } - if (cutcylBounds != nullptr) { - bUtility += - BinUtility(1, cutcylBounds->get(CutoutCylinderVolumeBounds::eMinR), - cutcylBounds->get(CutoutCylinderVolumeBounds::eMaxR), - Acts::open, Acts::binR); - bUtility += BinUtility(1, -M_PI, M_PI, Acts::closed, Acts::binPhi); - bUtility += BinUtility( - 1, -cutcylBounds->get(CutoutCylinderVolumeBounds::eHalfLengthZ), - cutcylBounds->get(CutoutCylinderVolumeBounds::eHalfLengthZ), Acts::open, - Acts::binZ); - return bUtility; - } else if (cuBounds != nullptr) { - bUtility += BinUtility(1, -cuBounds->get(CuboidVolumeBounds::eHalfLengthX), - cuBounds->get(CuboidVolumeBounds::eHalfLengthX), - Acts::open, Acts::binX); - bUtility += BinUtility(1, -cuBounds->get(CuboidVolumeBounds::eHalfLengthY), - cuBounds->get(CuboidVolumeBounds::eHalfLengthY), - Acts::open, Acts::binY); - bUtility += BinUtility(1, -cuBounds->get(CuboidVolumeBounds::eHalfLengthZ), - cuBounds->get(CuboidVolumeBounds::eHalfLengthZ), - Acts::open, Acts::binZ); - return bUtility; - } - ACTS_INFO( - "No corresponding bound found for the volume : " << volume.volumeName()); - return bUtility; -} \ No newline at end of file diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsLayerBuilder.cxx b/Tracking/Acts/FaserActsGeometry/src/FaserActsLayerBuilder.cxx index f8803a27..8ed2d789 100644 --- a/Tracking/Acts/FaserActsGeometry/src/FaserActsLayerBuilder.cxx +++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsLayerBuilder.cxx @@ -1,6 +1,6 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration -*/ + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + */ // ATHENA #include "TrackerReadoutGeometry/SiDetectorElement.h" #include "TrackerReadoutGeometry/SiDetectorElementCollection.h" @@ -11,6 +11,7 @@ #include "FaserActsGeometry/FaserActsDetectorElement.h" #include "FaserActsGeometry/CuboidVolumeBuilder.h" #include "ActsInterop/IdentityHelper.h" +#include "GeoModelKernel/GeoBox.h" // ACTS #include "Acts/Material/ProtoSurfaceMaterial.hpp" @@ -23,6 +24,7 @@ #include "Acts/Geometry/LayerCreator.hpp" #include "Acts/Geometry/GeometryContext.hpp" #include "Acts/Definitions/Common.hpp" +#include "Acts/Geometry/PassiveLayerBuilder.hpp" #include "Acts/Definitions/Algebra.hpp" #include "Acts/Definitions/Units.hpp" #include "Acts/Utilities/BinningType.hpp" @@ -35,69 +37,110 @@ using namespace Acts::UnitLiterals; FaserActs::CuboidVolumeBuilder::Config FaserActsLayerBuilder::buildVolume(const Acts::GeometryContext& gctx) { - //Acts::CuboidVolumeBuilder::Config cvbConfig; FaserActs::CuboidVolumeBuilder::Config cvbConfig; - //std::vector<Acts::CuboidVolumeBuilder::VolumeConfig> volumeConfigs = {}; std::vector<FaserActs::CuboidVolumeBuilder::VolumeConfig> volumeConfigs = {}; + //build the tracker stations // iStation starts from 1 for the main detector (Faser-01 geometry) and from 0 if we include the IFT (Faser-02 geometry) auto siDetMng = static_cast<const TrackerDD::SCT_DetectorManager*>(m_cfg.mng); int nStations = siDetMng->numerology().numStations(); + for (int iStation=4-nStations; iStation<4; iStation++) { - //Acts::CuboidVolumeBuilder::VolumeConfig volumeConfig; - FaserActs::CuboidVolumeBuilder::VolumeConfig volumeConfig; + FaserActs::CuboidVolumeBuilder::VolumeConfig volumeConfig; + + std::vector<FaserActs::CuboidVolumeBuilder::LayerConfig> layerConfigs; + layerConfigs.clear(); - //std::vector<Acts::CuboidVolumeBuilder::LayerConfig> layerConfigs; - std::vector<FaserActs::CuboidVolumeBuilder::LayerConfig> layerConfigs; - layerConfigs.clear(); + for (int iPlane=0; iPlane<3; iPlane++) { - for (int iPlane=0; iPlane<3; iPlane++) { + m_cfg.station = iStation; + m_cfg.plane = iPlane; + std::vector<std::shared_ptr<const Acts::Surface>> surfaces; + surfaces.clear(); - m_cfg.station = iStation; - m_cfg.plane = iPlane; - std::vector<std::shared_ptr<const Acts::Surface>> surfaces; - surfaces.clear(); + FaserActs::CuboidVolumeBuilder::SurfaceConfig surfacecfg; - //Acts::CuboidVolumeBuilder::SurfaceConfig surfacecfg; - FaserActs::CuboidVolumeBuilder::SurfaceConfig surfacecfg; - - //Acts::CuboidVolumeBuilder::LayerConfig layercfg; - FaserActs::CuboidVolumeBuilder::LayerConfig layercfg; - layercfg.binsX = 2; - layercfg.binsY = 4; + FaserActs::CuboidVolumeBuilder::LayerConfig layercfg; + layercfg.binsX = 2; + layercfg.binsY = 4; - buildLayers(gctx, surfaces, layercfg, surfacecfg); + buildLayers(gctx, surfaces, layercfg, surfacecfg); - layercfg.surfaceCfg = surfacecfg; - layercfg.active = true; + layercfg.surfaceCfg = surfacecfg; + layercfg.active = true; - layercfg.surfaces = surfaces; - layerConfigs.push_back(layercfg); + layercfg.surfaces = surfaces; + layerConfigs.push_back(layercfg); + + if (iPlane == 1) { + volumeConfig.position = Acts::Vector3(0, 0, surfacecfg.position.z()); + } + } + volumeConfig.length = m_trackerDimensions; + volumeConfig.layerCfg = layerConfigs; + volumeConfig.name = "Station_" + std::to_string(iStation); + volumeConfigs.emplace_back(volumeConfig); - if (iPlane == 1) { - volumeConfig.position = Acts::Vector3(0, 0, surfacecfg.position.z()); - } - if (iStation == 0 && iPlane == 1) { - cvbConfig.position = Acts::Vector3(0, 0, surfacecfg.position.z()); - } + } + //build the veto/trigger stations + //first build veto station, then trigger station + //veto station + //Veto A, VetoRadiator, Veto B + auto vetoManager = static_cast<const GeoVDetectorManager*>(m_cfg.vetomng); + for(unsigned int i=0; i<vetoManager->getNumTreeTops(); i++){ + auto vol0 = vetoManager->getTreeTop(i)->getLogVol(); + //get the shape params and translation + const GeoBox* shape = dynamic_cast<const GeoBox*>(vol0->getShape()); + auto trans = vetoManager->getTreeTop(i)->getX(); + if(vetoManager->getTreeTop(i)->getNChildVols()==0){ + volumeConfigs.emplace_back(buildScintVolume(trans.translation().z(), shape->getZHalfLength(),vol0->getName())); + } + else{ + for(size_t j =0 ;j< vetoManager->getTreeTop(i)->getNChildVols();j++){ + auto childtrans=vetoManager->getTreeTop(i)->getXToChildVol(j); + const GeoBox* childshape = dynamic_cast<const GeoBox*>(vetoManager->getTreeTop(i)->getChildVol(j)->getLogVol()->getShape()); + volumeConfigs.emplace_back(buildScintVolume(trans.translation().z() + childtrans.translation().z(), childshape->getZHalfLength(),vol0->getName()+std::string("_")+std::to_string(j))); + } + } + } + //build trigger stations + auto triggerManager = static_cast<const GeoVDetectorManager*>(m_cfg.triggermng); + for(unsigned int i=0; i<triggerManager->getNumTreeTops(); i++){ + auto vol0 = triggerManager->getTreeTop(i)->getLogVol(); + //get the shape params and translation + const GeoBox* shape = dynamic_cast<const GeoBox*>(vol0->getShape()); + auto trans = triggerManager->getTreeTop(i)->getX(); + if(triggerManager->getTreeTop(i)->getNChildVols()==0){ + volumeConfigs.emplace_back(buildScintVolume(trans.translation().z(), shape->getZHalfLength(),vol0->getName())); + } + else{ + for(size_t j =0 ;j< triggerManager->getTreeTop(i)->getNChildVols();j++){ + auto childtrans=triggerManager->getTreeTop(i)->getXToChildVol(j); + const GeoBox* childshape = dynamic_cast<const GeoBox*>(triggerManager->getTreeTop(i)->getChildVol(j)->getLogVol()->getShape()); + volumeConfigs.emplace_back(buildScintVolume(trans.translation().z() + childtrans.translation().z(), childshape->getZHalfLength(),vol0->getName()+std::string("_")+std::to_string(j))); + } } + } + - volumeConfig.length = m_trackerDimensions; - volumeConfig.layerCfg = layerConfigs; - volumeConfig.name = "Station_" + std::to_string(iStation); - volumeConfigs.push_back(volumeConfig); + //sort the geometry by the Z position, neccessary to have the correct boundary + auto sortFunc = [](FaserActs::CuboidVolumeBuilder::VolumeConfig& v1, FaserActs::CuboidVolumeBuilder::VolumeConfig& v2){return v1.position.z()<v2.position.z();}; + std::sort(volumeConfigs.begin(),volumeConfigs.end(), sortFunc); + //check volume positions + if (logger().doPrint(Acts::Logging::VERBOSE)) { + for(auto iter: volumeConfigs) ACTS_VERBOSE(" build volume centerred at Z = " << iter.position.z()); } - - //cvbConfig.position = m_worldCenter; + + cvbConfig.position = m_worldCenter; cvbConfig.length = m_worldDimensions; cvbConfig.volumeCfg = volumeConfigs; return cvbConfig; } -void + void FaserActsLayerBuilder::buildLayers(const Acts::GeometryContext& gctx, std::vector<std::shared_ptr<const Surface>>& surfaces, FaserActs::CuboidVolumeBuilder::LayerConfig& layercfg, @@ -105,36 +148,35 @@ FaserActsLayerBuilder::buildLayers(const Acts::GeometryContext& gctx, { auto siDetMng = static_cast<const TrackerDD::SCT_DetectorManager*>(m_cfg.mng); - + for (int iRow = 0; iRow < 4; iRow++) { - for (int iModule = -1; iModule < 2; iModule++) { - for (int iSensor = 0; iSensor < 2; iSensor++) { - //for (int iSensor = 0; iSensor < 1; iSensor++) { // only use the first sensor to construct the surface + for (int iModule = -1; iModule < 2; iModule++) { + for (int iSensor = 0; iSensor < 2; iSensor++) { - if (iModule == 0) continue; - const TrackerDD::SiDetectorElement* siDetElement = siDetMng->getDetectorElement(m_cfg.station, m_cfg.plane, iRow, iModule, iSensor) ; + if (iModule == 0) continue; + const TrackerDD::SiDetectorElement* siDetElement = siDetMng->getDetectorElement(m_cfg.station, m_cfg.plane, iRow, iModule, iSensor) ; - if (logger().doPrint(Acts::Logging::VERBOSE)) { - ACTS_VERBOSE("Found SCT sensor (" << m_cfg.station << "/" << m_cfg.plane << "/" << iRow << "/" << iModule << "/" << iSensor << ") with global center at (" << siDetElement->center().x() << ", " << siDetElement->center().y() << ", " << siDetElement->center().z() << ")." ); - } + if (logger().doPrint(Acts::Logging::VERBOSE)) { + ACTS_VERBOSE("Found SCT sensor (" << m_cfg.station << "/" << m_cfg.plane << "/" << iRow << "/" << iModule << "/" << iSensor << ") with global center at (" << siDetElement->center().x() << ", " << siDetElement->center().y() << ", " << siDetElement->center().z() << ")." ); + } - auto element = std::make_shared<const FaserActsDetectorElement>(siDetElement); + auto element = std::make_shared<const FaserActsDetectorElement>(siDetElement); - surfaces.push_back(element->surface().getSharedPtr()); + surfaces.push_back(element->surface().getSharedPtr()); - m_cfg.elementStore->push_back(element); - - m_ModuleWidth = siDetElement->width(); - m_ModuleLength = siDetElement->length(); - } - } + m_cfg.elementStore->push_back(element); + + m_ModuleWidth = siDetElement->width(); + m_ModuleLength = siDetElement->length(); + } + } } Acts::ProtoLayer pl(gctx, surfaces); if (logger().doPrint(Acts::Logging::VERBOSE)) { - ACTS_VERBOSE(" Plane's zMin / zMax: " << pl.min(Acts::binZ) << " / " << pl.max(Acts::binZ)); + ACTS_VERBOSE(" Plane's zMin / zMax: " << pl.min(Acts::binZ) << " / " << pl.max(Acts::binZ)); } std::shared_ptr<const Acts::ProtoSurfaceMaterial> materialProxy = nullptr; @@ -148,54 +190,118 @@ FaserActsLayerBuilder::buildLayers(const Acts::GeometryContext& gctx, if (std::abs(layerZInner) > std::abs(layerZOuter)) std::swap(layerZInner, layerZOuter); - auto rBounds = std::make_shared<const Acts::RectangleBounds>( 0.5*layercfg.binsY*m_ModuleWidth, 0.5*layercfg.binsX*m_ModuleLength ) ; - - Transform3 transformNominal(Translation3(0., 0., layerZ)); - - Transform3 transformInner(Translation3(0., 0., layerZInner)); - - Transform3 transformOuter(Translation3(0., 0., layerZOuter)); - - std::shared_ptr<Acts::PlaneSurface> innerBoundary - = Acts::Surface::makeShared<Acts::PlaneSurface>(transformInner, rBounds); - - std::shared_ptr<Acts::PlaneSurface> nominalSurface - = Acts::Surface::makeShared<Acts::PlaneSurface>(transformNominal, rBounds); - - std::shared_ptr<Acts::PlaneSurface> outerBoundary - = Acts::Surface::makeShared<Acts::PlaneSurface>(transformOuter, rBounds); - - size_t matBinsX = layercfg.binsX; - size_t matBinsY = layercfg.binsY; - - Acts::BinUtility materialBinUtil( - matBinsX, -0.5*layercfg.binsY*m_ModuleWidth, 0.5*layercfg.binsY*m_ModuleWidth, Acts::open, Acts::binX); - materialBinUtil += Acts::BinUtility( - matBinsY, -0.5*layercfg.binsX*m_ModuleLength, 0.5*layercfg.binsX*m_ModuleLength, Acts::open, Acts::binY, transformInner); - - materialProxy - = std::make_shared<const Acts::ProtoSurfaceMaterial>(materialBinUtil); - - ACTS_VERBOSE("[L] Layer is marked to carry support material on Surface ( " - "inner=0 / center=1 / outer=2 ) : " << "inner"); - ACTS_VERBOSE("with binning: [" << matBinsX << ", " << matBinsY << "]"); - - ACTS_VERBOSE("Created ApproachSurfaces for layer at:"); - ACTS_VERBOSE(" - inner: Z=" << layerZInner); - ACTS_VERBOSE(" - central: Z=" << layerZ); - ACTS_VERBOSE(" - outer: Z=" << layerZOuter); - - - // set material on inner - // @TODO: make this configurable somehow - innerBoundary->assignSurfaceMaterial(materialProxy); - - std::vector<std::shared_ptr<const Acts::Surface>> aSurfaces; - aSurfaces.push_back(std::move(innerBoundary)); - aSurfaces.push_back(std::move(nominalSurface)); - aSurfaces.push_back(std::move(outerBoundary)); - - layercfg.approachDescriptor = new Acts::GenericApproachDescriptor(aSurfaces); - + auto rBounds = std::make_shared<const Acts::RectangleBounds>( 0.5*layercfg.binsY*m_ModuleWidth, 0.5*layercfg.binsX*m_ModuleLength ) ; + + Transform3 transformNominal(Translation3(0., 0., layerZ)); + + Transform3 transformInner(Translation3(0., 0., layerZInner)); + + Transform3 transformOuter(Translation3(0., 0., layerZOuter)); + + std::shared_ptr<Acts::PlaneSurface> innerBoundary + = Acts::Surface::makeShared<Acts::PlaneSurface>(transformInner, rBounds); + + std::shared_ptr<Acts::PlaneSurface> nominalSurface + = Acts::Surface::makeShared<Acts::PlaneSurface>(transformNominal, rBounds); + + std::shared_ptr<Acts::PlaneSurface> outerBoundary + = Acts::Surface::makeShared<Acts::PlaneSurface>(transformOuter, rBounds); + + size_t matBinsX = layercfg.binsX*m_cfg.MaterialBins.first; + size_t matBinsY = layercfg.binsY*m_cfg.MaterialBins.second; + + Acts::BinUtility materialBinUtil( + matBinsX, -0.5*layercfg.binsY*m_ModuleWidth, 0.5*layercfg.binsY*m_ModuleWidth, Acts::open, Acts::binX); + materialBinUtil += Acts::BinUtility( + matBinsY, -0.5*layercfg.binsX*m_ModuleLength, 0.5*layercfg.binsX*m_ModuleLength, Acts::open, Acts::binY, transformInner); + + materialProxy + = std::make_shared<const Acts::ProtoSurfaceMaterial>(materialBinUtil); + + ACTS_VERBOSE("[L] Layer is marked to carry support material on Surface ( " + "inner=0 / center=1 / outer=2 ) : " << "inner"); + ACTS_VERBOSE("with binning: [" << matBinsX << ", " << matBinsY << "]"); + + ACTS_VERBOSE("Created ApproachSurfaces for layer at:"); + ACTS_VERBOSE(" - inner: Z=" << layerZInner); + ACTS_VERBOSE(" - central: Z=" << layerZ); + ACTS_VERBOSE(" - outer: Z=" << layerZOuter); + + + // set material on inner + // @TODO: make this configurable somehow + innerBoundary->assignSurfaceMaterial(materialProxy); + + + std::vector<std::shared_ptr<const Acts::Surface>> aSurfaces; + aSurfaces.push_back(std::move(innerBoundary)); + aSurfaces.push_back(std::move(nominalSurface)); + aSurfaces.push_back(std::move(outerBoundary)); + + layercfg.approachDescriptor = new Acts::GenericApproachDescriptor(aSurfaces); + } +FaserActs::CuboidVolumeBuilder::VolumeConfig FaserActsLayerBuilder::buildScintVolume(double zpos, double zhalflength, std::string name){ + + FaserActs::CuboidVolumeBuilder::VolumeConfig volumeConfig; + + std::vector<FaserActs::CuboidVolumeBuilder::LayerConfig> layerConfigs; + layerConfigs.clear(); + + FaserActs::CuboidVolumeBuilder::LayerConfig layercfg; + FaserActs::CuboidVolumeBuilder::SurfaceConfig surfacecfg; + surfacecfg.position = Acts::Vector3(0, 0, zpos); + surfacecfg.thickness=zhalflength; + layercfg.active = true; + + std::vector<std::shared_ptr<const Acts::Surface>> surfaces; + surfaces.clear(); + + layercfg.binsX = 1; + layercfg.binsY = 1; + + layercfg.active = true; + + //bounds are hard coded + auto rBounds = std::make_shared<const Acts::RectangleBounds>( 130.*1_mm, 130.*1_mm) ; + surfacecfg.rBounds=rBounds; + layercfg.surfaceCfg = surfacecfg; + + + Transform3 transformInner(Translation3(0., 0., (zpos)*1_mm)); + + std::shared_ptr<Acts::PlaneSurface> innerBoundary + = Acts::Surface::makeShared<Acts::PlaneSurface>(transformInner, rBounds); + + + size_t matBinsX = m_cfg.MaterialBins.first; + size_t matBinsY = m_cfg.MaterialBins.second; + + Acts::BinUtility materialBinUtil( + matBinsX, -130*1_mm, 130.*1_mm, Acts::open, Acts::binX); + materialBinUtil += Acts::BinUtility( + matBinsY, -130*1_mm, 130.*1_mm, Acts::open, Acts::binY, transformInner); + + std::shared_ptr<const Acts::ProtoSurfaceMaterial> materialProxy + = std::make_shared<const Acts::ProtoSurfaceMaterial>(materialBinUtil); + innerBoundary->assignSurfaceMaterial(materialProxy); + + std::vector<std::shared_ptr<const Acts::Surface>> aSurfaces; + aSurfaces.push_back(std::move(innerBoundary)); + + layercfg.surfaces = aSurfaces; + layerConfigs.push_back(layercfg); + + volumeConfig.position = Acts::Vector3(0, 0, zpos); + Acts::Vector3 length= { 300.0*1_mm, 300.0*1_mm, zhalflength*2.*1_mm }; + volumeConfig.length = length; + //volumeConfig.layerCfg = {}; + volumeConfig.layerCfg = layerConfigs; + volumeConfig.name = name; + + return volumeConfig; + +} + + diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsMaterialJsonWriterTool.cxx b/Tracking/Acts/FaserActsGeometry/src/FaserActsMaterialJsonWriterTool.cxx index 86c97a2c..b138bc60 100644 --- a/Tracking/Acts/FaserActsGeometry/src/FaserActsMaterialJsonWriterTool.cxx +++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsMaterialJsonWriterTool.cxx @@ -6,6 +6,7 @@ #include "ActsInterop/Logger.h" + #include <fstream> #include <ios> #include <iostream> @@ -26,8 +27,8 @@ FaserActsMaterialJsonWriterTool::initialize() { ATH_MSG_INFO("Starting Material writer"); - m_cfg.name = "FaserActsJsonGeometryConverter"; - m_cfg.logger = makeActsAthenaLogger(this, "FaserActsJsonGeometryConverter"); +// m_cfg.name = "FaserActsJsonGeometryConverter"; +// m_cfg.logger = makeActsAthenaLogger(this, "FaserActsJsonGeometryConverter"); m_cfg.processSensitives = m_processSensitives; m_cfg.processApproaches = m_processApproaches; m_cfg.processRepresenting = m_processRepresenting; @@ -42,8 +43,9 @@ FaserActsMaterialJsonWriterTool::initialize() void FaserActsMaterialJsonWriterTool::write(const Acts::MaterialMapJsonConverter::DetectorMaterialMaps& detMaterial) const { + ATH_MSG_INFO("Start writing material"); // Evoke the converter - Acts::MaterialMapJsonConverter jmConverter(m_cfg); + Acts::MaterialMapJsonConverter jmConverter(m_cfg, Acts::Logging::INFO); auto jout = jmConverter.materialMapsToJson(detMaterial); // And write the file std::ofstream ofj(m_filePath); @@ -53,8 +55,9 @@ FaserActsMaterialJsonWriterTool::write(const Acts::MaterialMapJsonConverter::Det void FaserActsMaterialJsonWriterTool::write(const Acts::TrackingGeometry& tGeometry) const { + ATH_MSG_INFO("Start writing geometry"); // Evoke the converter - Acts::MaterialMapJsonConverter jmConverter(m_cfg); + Acts::MaterialMapJsonConverter jmConverter(m_cfg, Acts::Logging::INFO); auto jout = jmConverter.trackingGeometryToJson(tGeometry); // And write the file std::ofstream ofj(m_filePath); diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsMaterialJsonWriterTool.h b/Tracking/Acts/FaserActsGeometry/src/FaserActsMaterialJsonWriterTool.h index c29a42c5..2e3b95fa 100644 --- a/Tracking/Acts/FaserActsGeometry/src/FaserActsMaterialJsonWriterTool.h +++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsMaterialJsonWriterTool.h @@ -2,8 +2,8 @@ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration */ -#ifndef ACTSGEOMETRY_FASERACTSMATERIALJSONWRITERTOOL_H -#define ACTSGEOMETRY_FASERACTSMATERIALJSONWRITERTOOL_H +#ifndef FASERACTSGEOMETRY_FASERACTSMATERIALJSONWRITERTOOL_H +#define FASERACTSGEOMETRY_FASERACTSMATERIALJSONWRITERTOOL_H // ATHENA #include "AthenaBaseComps/AthAlgTool.h" diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsMaterialMapping.cxx b/Tracking/Acts/FaserActsGeometry/src/FaserActsMaterialMapping.cxx index fd3e6e15..9edaa72f 100644 --- a/Tracking/Acts/FaserActsGeometry/src/FaserActsMaterialMapping.cxx +++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsMaterialMapping.cxx @@ -1,6 +1,6 @@ /* - Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration -*/ + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration + */ #include "FaserActsMaterialMapping.h" @@ -33,12 +33,12 @@ #include "Acts/Propagator/StraightLineStepper.hpp" FaserActsMaterialMapping::FaserActsMaterialMapping(const std::string &name, - ISvcLocator *pSvcLocator) - : AthReentrantAlgorithm(name, pSvcLocator), - m_materialTrackWriterSvc("ActsMaterialTrackWriterSvc", name), - m_inputMaterialStepCollection("MaterialStepRecords"), - m_mappingState(m_gctx,m_mctx), - m_mappingStateVol(m_gctx,m_mctx) + ISvcLocator *pSvcLocator) + : AthReentrantAlgorithm(name, pSvcLocator), + m_materialTrackWriterSvc("FaserActsMaterialTrackWriterSvc", name), + m_inputMaterialStepCollection("MaterialStepRecords"), + m_mappingState(m_gctx,m_mctx), + m_mappingStateVol(m_gctx,m_mctx) {} StatusCode FaserActsMaterialMapping::initialize() { @@ -49,7 +49,6 @@ StatusCode FaserActsMaterialMapping::initialize() { return StatusCode::FAILURE; } - ATH_CHECK(m_materialStepConverterTool.retrieve() ); ATH_CHECK(m_materialTrackWriterSvc.retrieve() ); if(m_mapSurfaces){ ATH_CHECK(m_surfaceMappingTool.retrieve() ); @@ -60,7 +59,7 @@ StatusCode FaserActsMaterialMapping::initialize() { m_mappingStateVol = m_volumeMappingTool->mappingState(); } ATH_CHECK(m_materialJsonWriterTool.retrieve() ); - ATH_CHECK( m_inputMaterialStepCollection.initialize() ); + ATH_CHECK(m_inputMaterialStepCollection.initialize() ); return StatusCode::SUCCESS; } @@ -69,16 +68,87 @@ StatusCode FaserActsMaterialMapping::execute(const EventContext &ctx) const { Acts::RecordedMaterialTrack mTrack; SG::ReadHandle<Trk::MaterialStepCollection> materialStepCollection(m_inputMaterialStepCollection, ctx); - mTrack = m_materialStepConverterTool->convertToMaterialTrack(*materialStepCollection); + + //make material track collection + Trk::MaterialStepCollection::const_iterator iter = materialStepCollection->begin(); + Trk::MaterialStepCollection::const_iterator iter_end = materialStepCollection->end(); + Trk::MaterialStepCollection* newmaterialStepCollection= new Trk::MaterialStepCollection; + //remove the material at Z>2500mm + //need further study + for(; iter != iter_end; ++iter){ + const Trk::MaterialStep* mat = *iter; + Trk::MaterialStep* newmat = const_cast<Trk::MaterialStep*>(mat); + //remove the magnets, need to be updated + if((newmat->hitZ()<-1600&&newmat->hitZ()>-1920)||(newmat->hitZ()>-200&&newmat->hitZ()<100)||(newmat->hitZ()>1000&&newmat->hitZ()<1300)||(newmat->hitZ()>2000&&newmat->hitZ()<2500)) + newmaterialStepCollection->push_back(newmat); + } + if(newmaterialStepCollection->size()<1) + return StatusCode::SUCCESS; + + /* test */ + std::vector<Acts::MaterialInteraction> nStep; + Acts::RecordedMaterial recorded; + double sum_X0 = 0; + double sum_L0 = 0; + + double x_length = newmaterialStepCollection->back()->hitX() - newmaterialStepCollection->front()->hitX(); + double y_length = newmaterialStepCollection->back()->hitY() - newmaterialStepCollection->front()->hitY(); + double z_length = newmaterialStepCollection->back()->hitZ() - newmaterialStepCollection->front()->hitZ(); + + double norm = 1/(std::sqrt(x_length*x_length + + y_length*y_length + + z_length*z_length)); + //set the initial z position to -4000 + double z_ini=-4000.; + Acts::Vector3 v_pos{newmaterialStepCollection->front()->hitX()-(newmaterialStepCollection->front()->hitZ()-z_ini)*(newmaterialStepCollection->front()->hitX()-newmaterialStepCollection->back()->hitX())/(newmaterialStepCollection->front()->hitZ()-newmaterialStepCollection->back()->hitZ()),newmaterialStepCollection->front()->hitY()-(newmaterialStepCollection->front()->hitZ()-z_ini)*(newmaterialStepCollection->front()->hitY()-newmaterialStepCollection->back()->hitY())/(newmaterialStepCollection->front()->hitZ()-newmaterialStepCollection->back()->hitZ()),z_ini}; + Acts::Vector3 v_imp{x_length*norm, y_length*norm, z_length*norm}; + Acts::Vector3 prev_pos = v_pos; + for(auto const step: *newmaterialStepCollection) { + + Acts::MaterialInteraction interaction; + + Acts::Vector3 pos{step->hitX(), step->hitY(), step->hitZ()}; + Acts::MaterialSlab matProp(Acts::Material::fromMassDensity(step->x0(), step->l0(), step->A(), step->Z(), (step->rho() * Acts::UnitConstants::g) ),step->steplength()); + interaction.position = pos; + + double x_dir = pos.x() - prev_pos.x(); + double y_dir = pos.y() - prev_pos.y(); + double z_dir = pos.z() - prev_pos.z(); + double norm_dir = 1/(std::sqrt(x_dir*x_dir + + y_dir*y_dir + + z_dir*z_dir)); + + Acts::Vector3 dir{x_dir * norm_dir, y_dir * norm_dir, z_dir * norm_dir}; + interaction.direction = dir; + prev_pos = pos; + interaction.materialSlab = matProp; + sum_X0 += step->steplengthInX0(); + sum_L0 += step->steplengthInL0(); + nStep.push_back(interaction); + } + recorded.materialInX0 = sum_X0; + recorded.materialInL0 = sum_L0; + recorded.materialInteractions = nStep; + + mTrack = std::make_pair(std::make_pair(v_pos, v_imp), recorded); if(m_mapSurfaces){ + ATH_MSG_INFO(name() << " start surface mapping"); auto mappingState - = const_cast<Acts::SurfaceMaterialMapper::State *>(&m_mappingState); + = const_cast<Acts::SurfaceMaterialMapper::State *>(&m_mappingState); + auto context = m_surfaceMappingTool->trackingGeometryTool()->getNominalGeometryContext().context(); + std::reference_wrapper<const Acts::GeometryContext> geoContext(context); + mappingState->geoContext = geoContext; + m_surfaceMappingTool->mapper()->mapMaterialTrack(*mappingState, mTrack); } if(m_mapVolumes){ + ATH_MSG_INFO(name() << " start volume mapping"); auto mappingStateVol - = const_cast<Acts::VolumeMaterialMapper::State *>(&m_mappingStateVol); + = const_cast<Acts::VolumeMaterialMapper::State *>(&m_mappingStateVol); + auto context = m_volumeMappingTool->trackingGeometryTool()->getNominalGeometryContext().context(); + std::reference_wrapper<const Acts::GeometryContext> geoContext(context); + mappingStateVol->geoContext = geoContext; m_volumeMappingTool->mapper()->mapMaterialTrack(*mappingStateVol, mTrack); } m_materialTrackWriterSvc->write(mTrack); @@ -108,28 +178,28 @@ StatusCode FaserActsMaterialMapping::finalize() { m_surfaceMappingTool->mapper()->finalizeMaps(m_mappingState); // Loop over the state, and collect the maps for surfaces for (auto& [key, value] : m_mappingState.surfaceMaterial) { - detectorMaterial.first.insert({key, std::move(value)}); + detectorMaterial.first.insert({key, std::move(value)}); } // Loop over the state, and collect the maps for volumes for (auto& [key, value] : m_mappingState.volumeMaterial) { - detectorMaterial.second.insert({key, std::move(value)}); + detectorMaterial.second.insert({key, std::move(value)}); } } if(m_mapVolumes){ m_volumeMappingTool->mapper()->finalizeMaps(m_mappingStateVol); // Loop over the state, and collect the maps for surfaces for (auto& [key, value] : m_mappingStateVol.surfaceMaterial) { - detectorMaterial.first.insert({key, std::move(value)}); + detectorMaterial.first.insert({key, std::move(value)}); } // Loop over the state, and collect the maps for volumes for (auto& [key, value] : m_mappingStateVol.volumeMaterial) { - detectorMaterial.second.insert({key, std::move(value)}); + detectorMaterial.second.insert({key, std::move(value)}); } } } - + m_materialJsonWriterTool->write(detectorMaterial); return StatusCode::SUCCESS; -} \ No newline at end of file +} diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsMaterialMapping.h b/Tracking/Acts/FaserActsGeometry/src/FaserActsMaterialMapping.h index bdeb9a7b..8ddd2a38 100644 --- a/Tracking/Acts/FaserActsGeometry/src/FaserActsMaterialMapping.h +++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsMaterialMapping.h @@ -22,7 +22,6 @@ #include "Acts/Material/VolumeMaterialMapper.hpp" // PACKAGE #include "ActsGeometryInterfaces/IActsMaterialTrackWriterSvc.h" -#include "ActsGeometryInterfaces/IActsMaterialStepConverterTool.h" #include "FaserActsGeometryInterfaces/IFaserActsSurfaceMappingTool.h" #include "FaserActsGeometryInterfaces/IFaserActsVolumeMappingTool.h" #include "FaserActsGeometryInterfaces/IFaserActsMaterialJsonWriterTool.h" @@ -59,7 +58,6 @@ private: ServiceHandle<IActsMaterialTrackWriterSvc> m_materialTrackWriterSvc; Gaudi::Property<bool> m_mapSurfaces{this, "mapSurfaces", true, "Map the material onto surfaces"}; Gaudi::Property<bool> m_mapVolumes{this, "mapVolumes", true, "Map the material onto volumes"}; - ToolHandle<IActsMaterialStepConverterTool> m_materialStepConverterTool{this, "MaterialStepConverterTool", "ActsMaterialStepConverterTool"}; SG::ReadHandleKey<Trk::MaterialStepCollection> m_inputMaterialStepCollection; ToolHandle<IFaserActsSurfaceMappingTool> m_surfaceMappingTool{this, "SurfaceMappingTool", "FaserActsSurfaceMappingTool"}; ToolHandle<IFaserActsVolumeMappingTool> m_volumeMappingTool{this, "VolumeMappingTool", "FaserActsVolumeMappingTool"}; @@ -71,4 +69,4 @@ private: Acts::VolumeMaterialMapper::State m_mappingStateVol; }; -#endif // ActsGeometry_ActsExtrapolation_h \ No newline at end of file +#endif // ActsGeometry_ActsExtrapolation_h diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsMaterialTrackWriterSvc.cxx b/Tracking/Acts/FaserActsGeometry/src/FaserActsMaterialTrackWriterSvc.cxx new file mode 100644 index 00000000..4be42568 --- /dev/null +++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsMaterialTrackWriterSvc.cxx @@ -0,0 +1,354 @@ +/* + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration + */ + +#include "FaserActsMaterialTrackWriterSvc.h" +#include "FaserActsGeometry/FaserActsGeometryContext.h" +#include "GaudiKernel/IInterface.h" + +#include "TTree.h" +#include "TFile.h" + +#include "Acts/Material/MaterialSlab.hpp" +#include <Acts/Geometry/GeometryIdentifier.hpp> +#include <Acts/Surfaces/CylinderBounds.hpp> +#include <Acts/Surfaces/RadialBounds.hpp> +#include <Acts/Utilities/Helpers.hpp> + +using namespace Acts::VectorHelpers; + +#include <vector> +#include <deque> +#include <mutex> +#include <thread> + +FaserActsMaterialTrackWriterSvc::FaserActsMaterialTrackWriterSvc( const std::string& name, ISvcLocator* svc ) + : base_class(name, svc), + m_trackingGeometrySvc("FaserActsTrackingGeometrySvc", name) { + } + + StatusCode +FaserActsMaterialTrackWriterSvc::initialize() +{ + + std::string filePath = m_filePath; + std::string treeName = m_treeName; + p_tFile = TFile::Open(filePath.c_str(), "RECREATE"); + p_tFile->cd(); + p_tree = new TTree(treeName.c_str(), treeName.c_str()); + + p_tree->Branch("v_x", &m_v_x); + p_tree->Branch("v_y", &m_v_y); + p_tree->Branch("v_z", &m_v_z); + p_tree->Branch("v_px", &m_v_px); + p_tree->Branch("v_py", &m_v_py); + p_tree->Branch("v_pz", &m_v_pz); + p_tree->Branch("v_phi", &m_v_phi); + p_tree->Branch("v_eta", &m_v_eta); + + p_tree->Branch("t_X0", &m_tX0); + p_tree->Branch("t_L0", &m_tL0); + + p_tree->Branch("mat_x", &m_step_x); + p_tree->Branch("mat_y", &m_step_y); + p_tree->Branch("mat_z", &m_step_z); + p_tree->Branch("mat_dx", &m_step_dx); + p_tree->Branch("mat_dy", &m_step_dy); + p_tree->Branch("mat_dz", &m_step_dz); + p_tree->Branch("mat_step_length", &m_step_length); + p_tree->Branch("mat_X0", &m_step_X0); + p_tree->Branch("mat_L0", &m_step_L0); + p_tree->Branch("mat_A", &m_step_A); + p_tree->Branch("mat_Z", &m_step_Z); + p_tree->Branch("mat_rho", &m_step_rho); + + if (m_storeSurface) { + p_tree->Branch("sur_id", &m_sur_id); + p_tree->Branch("sur_type", &m_sur_type); + p_tree->Branch("sur_x", &m_sur_x); + p_tree->Branch("sur_y", &m_sur_y); + p_tree->Branch("sur_z", &m_sur_z); + p_tree->Branch("sur_range_min", &m_sur_range_min); + p_tree->Branch("sur_range_max", &m_sur_range_max); + } + if (m_storeVolume) { + p_tree->Branch("vol_id", &m_vol_id); + } + + ATH_MSG_INFO("Starting writer thread"); + ATH_MSG_DEBUG("Maximum queue size is set to:" << m_maxQueueSize); + m_doEnd = false; + m_writeThread = std::thread(&FaserActsMaterialTrackWriterSvc::writerThread, this); + + return StatusCode::SUCCESS; +} + + StatusCode +FaserActsMaterialTrackWriterSvc::finalize() +{ + + ATH_MSG_INFO("Waiting for writer thread to finish."); + m_doEnd = true; + m_writeThread.join(); + ATH_MSG_INFO("Writer thread has terminated."); + + ATH_MSG_INFO("Closing TFile"); + p_tFile->cd(); + p_tree->FlushBaskets(); + p_tree->AutoSave(); + p_tree->Write(); + p_tFile->Write(); + p_tFile->Close(); + + return StatusCode::SUCCESS; +} + + void +FaserActsMaterialTrackWriterSvc::write(const Acts::RecordedMaterialTrack& mTrack) +{ + std::lock_guard<std::mutex> lock(m_writeMutex); + + ATH_MSG_VERBOSE("Appending material track to write queue"); + m_mTracks.push_back(mTrack); +} + + void +FaserActsMaterialTrackWriterSvc::writerThread() +{ + using namespace std::chrono_literals; + // wait until we have events + while(m_mTracks.size() == 0) { + std::this_thread::sleep_for(2s); + if (m_doEnd) return; + } + + ATH_MSG_DEBUG("Begin regular write loop"); + while(true) { + ATH_MSG_VERBOSE("Obtaining write lock"); + std::unique_lock<std::mutex> lock(m_writeMutex); + + if (m_mTracks.empty()) { + lock.unlock(); + if (!m_doEnd) { + ATH_MSG_VERBOSE("Queue was empty, delay next execution"); + std::this_thread::sleep_for(0.1s); + continue; + } else { + ATH_MSG_INFO("Writer thread caught termination signal. Shutting down."); + + return; + } + } + + + if(m_mTracks.size() < m_maxQueueSize) { + // just pop one + ATH_MSG_VERBOSE("Queue at " << m_mTracks.size() << "/" << m_maxQueueSize + << ": Pop entry and write"); + Acts::RecordedMaterialTrack mTrack = std::move(m_mTracks.front()); + m_mTracks.pop_front(); + // writing can now happen without lock + lock.unlock(); + doWrite(std::move(mTrack)); + } + else { + ATH_MSG_DEBUG("Queue at " << m_mTracks.size() << "/" << m_maxQueueSize + << ": Lock and write until empty"); + while(!m_mTracks.empty()) { + ATH_MSG_VERBOSE("Pop entry and write"); + // keep the lock! + Acts::RecordedMaterialTrack mTrack = std::move(m_mTracks.front()); + m_mTracks.pop_front(); + doWrite(std::move(mTrack)); + } + ATH_MSG_DEBUG("Queue is empty, continue"); + + } + } +} + + void +FaserActsMaterialTrackWriterSvc::doWrite(const Acts::RecordedMaterialTrack& mTrack) +{ + ATH_MSG_VERBOSE("Write to tree"); + size_t mints = mTrack.second.materialInteractions.size(); + + // Clearing the vector first + m_step_sx.clear(); + m_step_sy.clear(); + m_step_sz.clear(); + m_step_x.clear(); + m_step_y.clear(); + m_step_z.clear(); + m_step_ex.clear(); + m_step_ey.clear(); + m_step_ez.clear(); + m_step_dx.clear(); + m_step_dy.clear(); + m_step_dz.clear(); + m_step_length.clear(); + m_step_X0.clear(); + m_step_L0.clear(); + m_step_A.clear(); + m_step_Z.clear(); + m_step_rho.clear(); + + if(m_storeSurface){ + m_sur_id.clear(); + m_sur_type.clear(); + m_sur_x.clear(); + m_sur_y.clear(); + m_sur_z.clear(); + m_sur_range_min.clear(); + m_sur_range_max.clear(); + } + + if (m_storeVolume) { + m_vol_id.clear(); + } + // Reserve the vector then + m_step_sx.reserve(mints); + m_step_sy.reserve(mints); + m_step_sz.reserve(mints); + m_step_x.reserve(mints); + m_step_y.reserve(mints); + m_step_ez.reserve(mints); + m_step_ex.reserve(mints); + m_step_ey.reserve(mints); + m_step_ez.reserve(mints); + m_step_dx.reserve(mints); + m_step_dy.reserve(mints); + m_step_dz.reserve(mints); + m_step_length.reserve(mints); + m_step_X0.reserve(mints); + m_step_L0.reserve(mints); + m_step_A.reserve(mints); + m_step_Z.reserve(mints); + m_step_rho.reserve(mints); + + if(m_storeSurface){ + m_sur_id.reserve(mints); + m_sur_type.reserve(mints); + m_sur_x.reserve(mints); + m_sur_y.reserve(mints); + m_sur_z.reserve(mints); + m_sur_range_min.reserve(mints); + m_sur_range_max.reserve(mints); + } + if (m_storeVolume) { + m_vol_id.reserve(mints); + } + + // reset the global counter + m_tX0 = mTrack.second.materialInX0; + m_tL0 = mTrack.second.materialInL0; + + // set the track information at vertex + m_v_x = mTrack.first.first.x(); + m_v_y = mTrack.first.first.y(); + m_v_z = mTrack.first.first.z(); + m_v_px = mTrack.first.second.x(); + m_v_py = mTrack.first.second.y(); + m_v_pz = mTrack.first.second.z(); + m_v_phi = phi(mTrack.first.second); + m_v_eta = eta(mTrack.first.second); + + // an now loop over the material + for (auto& mint : mTrack.second.materialInteractions) { + // The material step position information + m_step_x.push_back(mint.position.x()); + m_step_y.push_back(mint.position.y()); + m_step_z.push_back(mint.position.z()); + m_step_dx.push_back(mint.direction.x()); + m_step_dy.push_back(mint.direction.y()); + m_step_dz.push_back(mint.direction.z()); + + Acts::Vector3 prePos + = mint.position - 0.5 * mint.pathCorrection * mint.direction; + Acts::Vector3 posPos + = mint.position + 0.5 * mint.pathCorrection * mint.direction; + m_step_sx.push_back(prePos.x()); + m_step_sy.push_back(prePos.y()); + m_step_sz.push_back(prePos.z()); + m_step_ex.push_back(posPos.x()); + m_step_ey.push_back(posPos.y()); + m_step_ez.push_back(posPos.z()); + + if (m_storeSurface) { + const Acts::Surface* surface = mint.surface; + Acts::GeometryIdentifier layerID; + if (surface) { + auto gctx = std::make_unique<FaserActsGeometryContext>(); + gctx->alignmentStore = m_trackingGeometrySvc->getNominalAlignmentStore(); + auto sfIntersection = surface->intersect( + gctx->context(), mint.position, mint.direction, true); + layerID = surface->geometryId(); + m_sur_id.push_back(layerID.value()); + m_sur_type.push_back(surface->type()); + m_sur_x.push_back(sfIntersection.intersection.position.x()); + m_sur_y.push_back(sfIntersection.intersection.position.y()); + m_sur_z.push_back(sfIntersection.intersection.position.z()); + + const Acts::SurfaceBounds& surfaceBounds = surface->bounds(); + const Acts::RadialBounds* radialBounds = + dynamic_cast<const Acts::RadialBounds*>(&surfaceBounds); + const Acts::CylinderBounds* cylinderBounds = + dynamic_cast<const Acts::CylinderBounds*>(&surfaceBounds); + + if (radialBounds) { + m_sur_range_min.push_back(radialBounds->rMin()); + m_sur_range_max.push_back(radialBounds->rMax()); + } else if (cylinderBounds) { + m_sur_range_min.push_back( + -cylinderBounds->get(Acts::CylinderBounds::eHalfLengthZ)); + m_sur_range_max.push_back( + cylinderBounds->get(Acts::CylinderBounds::eHalfLengthZ)); + } else { + m_sur_range_min.push_back(0); + m_sur_range_max.push_back(0); + } + } else { + layerID.setVolume(0); + layerID.setBoundary(0); + layerID.setLayer(0); + layerID.setApproach(0); + layerID.setSensitive(0); + m_sur_id.push_back(layerID.value()); + m_sur_type.push_back(-1); + + m_sur_x.push_back(0); + m_sur_y.push_back(0); + m_sur_z.push_back(0); + m_sur_range_min.push_back(0); + m_sur_range_max.push_back(0); + } + } + // store volume information + if (m_storeVolume) { + const Acts::Volume* volume = mint.volume; + Acts::GeometryIdentifier vlayerID; + if (volume) { + vlayerID = volume->geometryId(); + m_vol_id.push_back(vlayerID.value()); + } else { + vlayerID.setVolume(0); + vlayerID.setBoundary(0); + vlayerID.setLayer(0); + vlayerID.setApproach(0); + vlayerID.setSensitive(0); + m_vol_id.push_back(vlayerID.value()); + } + } + // the material information + const auto& mprops = mint.materialSlab; + m_step_length.push_back(mprops.thickness()); + m_step_X0.push_back(mprops.material().X0()); + m_step_L0.push_back(mprops.material().L0()); + m_step_A.push_back(mprops.material().Ar()); + m_step_Z.push_back(mprops.material().Z()); + m_step_rho.push_back(mprops.material().massDensity()); + + } + p_tree->Fill(); + ATH_MSG_VERBOSE("Write complete"); +} diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsMaterialTrackWriterSvc.h b/Tracking/Acts/FaserActsGeometry/src/FaserActsMaterialTrackWriterSvc.h new file mode 100644 index 00000000..962ffbb9 --- /dev/null +++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsMaterialTrackWriterSvc.h @@ -0,0 +1,109 @@ +/* + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration + */ + +#ifndef FASERACTSGEOMETRY_FASERACTSMATERIALTRACKWRITERSVC_H +#define FASERACTSGEOMETRY_FASERACTSMATERIALTRACKWRITERSVC_H + +#include "ActsGeometryInterfaces/IActsMaterialTrackWriterSvc.h" +#include "FaserActsGeometryInterfaces/IFaserActsTrackingGeometrySvc.h" + +#include "AthenaBaseComps/AthService.h" +#include "StoreGate/StoreGateSvc.h" +//#include "GaudiKernel/IInterface.h" +//#include "Gaudi/Property.h" /*no forward decl: typedef*/ + +#include <vector> +#include <deque> +#include <mutex> +#include <thread> +#include <atomic> + +#include "TTree.h" +#include "TFile.h" + +class FaserActsMaterialTrackWriterSvc : public extends<AthService, IActsMaterialTrackWriterSvc> { + public: + + virtual StatusCode initialize() override; + virtual StatusCode finalize() override; + + FaserActsMaterialTrackWriterSvc( const std::string& name, ISvcLocator* svc ); + + virtual void + write(const Acts::RecordedMaterialTrack& mTrack) override; + + private: + + std::deque<Acts::RecordedMaterialTrack> m_mTracks; + std::mutex m_writeMutex; + std::thread m_writeThread; + std::atomic<bool> m_doEnd; + TFile* p_tFile; + TTree* p_tree; + + float m_v_x; ///< start global x + float m_v_y; ///< start global y + float m_v_z; ///< start global z + float m_v_px; ///< start global momentum x + float m_v_py; ///< start global momentum y + float m_v_pz; ///< start global momentum z + float m_v_phi; ///< start phi direction + float m_v_eta; ///< start eta direction + float m_tX0; ///< thickness in X0/L0 + float m_tL0; ///< thickness in X0/L0 + + std::vector<float> m_step_sx; ///< step x (start) position (optional) + std::vector<float> m_step_sy; ///< step y (start) position (optional) + std::vector<float> m_step_sz; ///< step z (start) position (optional) + std::vector<float> m_step_x; ///< step x position + std::vector<float> m_step_y; ///< step y position + std::vector<float> m_step_z; ///< step z position + std::vector<float> m_step_ex; ///< step x (end) position (optional) + std::vector<float> m_step_ey; ///< step y (end) position (optional) + std::vector<float> m_step_ez; ///< step z (end) position (optional) + std::vector<float> m_step_dx; ///< step x direction + std::vector<float> m_step_dy; ///< step y direction + std::vector<float> m_step_dz; ///< step z direction + std::vector<float> m_step_length; ///< step length + std::vector<float> m_step_X0; ///< step material x0 + std::vector<float> m_step_L0; ///< step material l0 + std::vector<float> m_step_A; ///< step material A + std::vector<float> m_step_Z; ///< step material Z + std::vector<float> m_step_rho; ///< step material rho + + std::vector<std::uint64_t> + m_sur_id; ///< ID of the suface associated with the step + std::vector<int32_t> + m_sur_type; ///< Type of the suface associated with the step + std::vector<float> m_sur_x; ///< x position of the center of the suface + ///< associated with the step + std::vector<float> m_sur_y; ///< y position of the center of the suface + ///< associated with the step + std::vector<float> m_sur_z; ///< z position of the center of the suface + ///< associated with the step + + std::vector<float> + m_sur_range_min; ///< Min range of the suface associated with the step + std::vector<float> + m_sur_range_max; ///< Max range of the suface associated with the step + + std::vector<std::uint64_t> + m_vol_id; ///< ID of the volume associated with the step + + void writerThread(); + void doWrite(const Acts::RecordedMaterialTrack &mTrack); + + ServiceHandle<IFaserActsTrackingGeometrySvc> m_trackingGeometrySvc; + + // jobOptions properties + Gaudi::Property<std::string> m_filePath{this, "FilePath", "MaterialTracks.root", "Output root file for charged particle"}; + Gaudi::Property<std::string> m_treeName{this, "TreeName", "material-tracks", ""}; + Gaudi::Property<bool> m_storeSurface{this, "StoreSurface", true, "Store the surface info in the root file"}; + Gaudi::Property<bool> m_storeVolume{this, "StoreVolume", true, "Store the volume info in the root file"}; + Gaudi::Property<size_t> m_maxQueueSize{this, "MaxQueueSize", 5000, "Limit the write queue to this size"}; + +}; + + +#endif diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsPropStepRootWriterSvc.h b/Tracking/Acts/FaserActsGeometry/src/FaserActsPropStepRootWriterSvc.h index 42a3b526..2f840479 100644 --- a/Tracking/Acts/FaserActsGeometry/src/FaserActsPropStepRootWriterSvc.h +++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsPropStepRootWriterSvc.h @@ -62,10 +62,6 @@ private: // jobOptions properties Gaudi::Property<std::string> m_filePath{this, "FilePath", "propsteps.root", "Output root file for charged particle"}; Gaudi::Property<std::string> m_treeName{this, "TreeName", "propsteps", ""}; - //Gaudi::Property<bool> m_writeBoundary{this, "WriteBoundary", true, ""}; - //Gaudi::Property<bool> m_writeMaterial{this, "WriteMaterial", true, ""}; - //Gaudi::Property<bool> m_writeSensitive{this, "WriteSensitive", true, ""}; - //Gaudi::Property<bool> m_writePassive{this, "WritePassive", true, ""}; // root branch storage TFile* m_outputFile; ///< the output file diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsSurfaceMappingTool.cxx b/Tracking/Acts/FaserActsGeometry/src/FaserActsSurfaceMappingTool.cxx index ac5494ad..eff504bb 100644 --- a/Tracking/Acts/FaserActsGeometry/src/FaserActsSurfaceMappingTool.cxx +++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsSurfaceMappingTool.cxx @@ -64,4 +64,4 @@ FaserActsSurfaceMappingTool::mappingState() const m_geoContext, m_magFieldContext, *m_trackingGeometry); return mappingState; -} \ No newline at end of file +} diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsTrackingGeometrySvc.cxx b/Tracking/Acts/FaserActsGeometry/src/FaserActsTrackingGeometrySvc.cxx index a8091afa..c6f476e7 100644 --- a/Tracking/Acts/FaserActsGeometry/src/FaserActsTrackingGeometrySvc.cxx +++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsTrackingGeometrySvc.cxx @@ -9,6 +9,10 @@ #include "StoreGate/StoreGateSvc.h" #include "GaudiKernel/EventContext.h" #include "TrackerIdentifier/FaserSCT_ID.h" +//#include "DipoleGeoModel/DipoleManager.h" +#include "GeoModelKernel/GeoTube.h" +#include "GeoModelKernel/GeoBox.h" +#include "GeoModelFaserUtilities/GeoModelExperiment.h" // ACTS #include "Acts/Geometry/GeometryContext.hpp" @@ -57,6 +61,12 @@ FaserActsTrackingGeometrySvc::initialize() ATH_CHECK ( m_detStore->retrieve(p_SCTManager, "SCT") ); + //get all the subdetectors + const GeoModelExperiment* theExpt = nullptr; + ATH_CHECK( m_detStore->retrieve( theExpt )); + const GeoVDetectorManager *vetoManager = theExpt->getManager("Veto"); + const GeoVDetectorManager *triggerManager = theExpt->getManager("Trigger"); + Acts::LayerArrayCreator::Config lacCfg; auto layerArrayCreator = std::make_shared<const Acts::LayerArrayCreator>( lacCfg, makeActsAthenaLogger(this, "LayArrCrtr", "ActsTGSvc")); @@ -93,11 +103,10 @@ FaserActsTrackingGeometrySvc::initialize() try { // SCT tgbConfig.trackingVolumeBuilders.push_back([&]( - const auto& gctx, const auto& /*inner*/, const auto&) { - auto tv = makeVolumeBuilder(gctx, p_SCTManager); - return tv->trackingVolume(gctx); - }); - + const auto& gctx, const auto& /*inner*/, const auto&) { + auto tv = makeVolumeBuilder(gctx, p_SCTManager, vetoManager, triggerManager); + return tv->trackingVolume(gctx); + }); } catch (const std::invalid_argument& e) { ATH_MSG_ERROR(e.what()); @@ -106,8 +115,8 @@ FaserActsTrackingGeometrySvc::initialize() auto trackingGeometryBuilder - = std::make_shared<const Acts::TrackingGeometryBuilder>(tgbConfig, - makeActsAthenaLogger(this, "TrkGeomBldr", "ActsTGSvc")); + = std::make_shared<const Acts::TrackingGeometryBuilder>(tgbConfig, + makeActsAthenaLogger(this, "TrkGeomBldr", "ActsTGSvc")); // default geometry context, this is nominal @@ -119,10 +128,13 @@ FaserActsTrackingGeometrySvc::initialize() ATH_MSG_VERBOSE("Building nominal alignment store"); FaserActsAlignmentStore* nominalAlignmentStore = new FaserActsAlignmentStore(); + ATH_MSG_VERBOSE("finish Building nominal alignment store"); populateAlignmentStore(nominalAlignmentStore); + ATH_MSG_VERBOSE("finish alignment population"); populateIdentifierMap(m_identifierMap.get()); + ATH_MSG_VERBOSE("get alignment store"); // manage ownership m_nominalAlignmentStore = std::unique_ptr<const FaserActsAlignmentStore>(nominalAlignmentStore); @@ -139,7 +151,7 @@ FaserActsTrackingGeometrySvc::trackingGeometry() { } std::shared_ptr<const Acts::ITrackingVolumeBuilder> -FaserActsTrackingGeometrySvc::makeVolumeBuilder(const Acts::GeometryContext& gctx, const TrackerDD::SCT_DetectorManager* manager) +FaserActsTrackingGeometrySvc::makeVolumeBuilder(const Acts::GeometryContext& gctx, const TrackerDD::SCT_DetectorManager* manager, const GeoVDetectorManager* vetoManager, const GeoVDetectorManager* triggerManager) { std::string managerName = manager->getName(); @@ -148,8 +160,13 @@ FaserActsTrackingGeometrySvc::makeVolumeBuilder(const Acts::GeometryContext& gct FaserActsLayerBuilder::Config cfg; cfg.subdetector = FaserActsDetectorElement::Subdetector::SCT; cfg.mng = static_cast<const TrackerDD::SCT_DetectorManager*>(manager); + //get veto and trigger manager + cfg.vetomng = static_cast<const GeoVDetectorManager*>(vetoManager); + cfg.triggermng = static_cast<const GeoVDetectorManager*>(triggerManager); cfg.elementStore = m_elementStore; cfg.identifierMap = m_identifierMap; + std::vector<size_t> matBins(m_MaterialBins); + cfg.MaterialBins = { matBins.at(0), matBins.at(1) }; gmLayerBuilder = std::make_shared<FaserActsLayerBuilder>(cfg, makeActsAthenaLogger(this, "GMLayBldr", "ActsTGSvc")); auto cvbConfig = gmLayerBuilder->buildVolume(gctx); @@ -163,11 +180,15 @@ FaserActsTrackingGeometrySvc::populateAlignmentStore(FaserActsAlignmentStore *st { // populate the alignment store with all detector elements m_trackingGeometry->visitSurfaces( - [store](const Acts::Surface* srf) { - const Acts::DetectorElementBase* detElem = srf->associatedDetectorElement(); - const auto* gmde = dynamic_cast<const FaserActsDetectorElement*>(detElem); - gmde->storeTransform(store); - }); + [store](const Acts::Surface* srf) { + FaserActsGeometryContext constructionContext; + constructionContext.construction = true; + const Acts::DetectorElementBase* detElem = srf->associatedDetectorElement(); + if(detElem){ + const auto* gmde = dynamic_cast<const FaserActsDetectorElement*>(detElem); + gmde->storeTransform(store); + } + }); } const FaserActsAlignmentStore* @@ -182,8 +203,10 @@ FaserActsTrackingGeometrySvc::populateIdentifierMap(IdentifierMap *map) const m_trackingGeometry->visitSurfaces( [map](const Acts::Surface* srf) { const Acts::DetectorElementBase* detElem = srf->associatedDetectorElement(); + if(detElem){ const auto* faserDetElem = dynamic_cast<const FaserActsDetectorElement*>(detElem); map->insert(std::make_pair(faserDetElem->identify(), srf->geometryId())); + } }); } @@ -191,4 +214,4 @@ std::shared_ptr<IdentifierMap> FaserActsTrackingGeometrySvc::getIdentifierMap() const { return m_identifierMap; -} \ No newline at end of file +} diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsTrackingGeometrySvc.h b/Tracking/Acts/FaserActsGeometry/src/FaserActsTrackingGeometrySvc.h index 65e21a13..83ef0efa 100644 --- a/Tracking/Acts/FaserActsGeometry/src/FaserActsTrackingGeometrySvc.h +++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsTrackingGeometrySvc.h @@ -9,6 +9,7 @@ #include "AthenaBaseComps/AthService.h" #include "StoreGate/StoreGateSvc.h" #include "GaudiKernel/EventContext.h" +#include "GeoModelFaserUtilities/GeoModelExperiment.h" // ACTS #include "Acts/Geometry/GeometryContext.hpp" @@ -25,6 +26,7 @@ using namespace Acts::UnitLiterals; namespace TrackerDD { class SCT_DetectorManager; + class DipoleManager; } class FaserActsAlignmentStore; @@ -35,6 +37,7 @@ namespace Acts { class TrackingGeometry; class CylinderVolumeHelper; +class CylinderVolumeBuilder; class ITrackingVolumeBuilder; class GeometryID; @@ -68,7 +71,7 @@ public: private: std::shared_ptr<const Acts::ITrackingVolumeBuilder> - makeVolumeBuilder(const Acts::GeometryContext& gctx, const TrackerDD::SCT_DetectorManager* manager); + makeVolumeBuilder(const Acts::GeometryContext& gctx, const TrackerDD::SCT_DetectorManager* manager, const GeoVDetectorManager* vetoManager, const GeoVDetectorManager* triggerManager); ServiceHandle<StoreGateSvc> m_detStore; const TrackerDD::SCT_DetectorManager* p_SCTManager; @@ -81,7 +84,7 @@ private: Gaudi::Property<bool> m_useMaterialMap{this, "UseMaterialMap", false, ""}; Gaudi::Property<std::string> m_materialMapInputFile{this, "MaterialMapInputFile", "", ""}; - Gaudi::Property<std::vector<size_t>> m_MaterialBins{this, "BarrelMaterialBins", {10, 10}}; + Gaudi::Property<std::vector<size_t>> m_MaterialBins{this, "MaterialBins", {30, 30}}; }; diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsWriteTrackingGeometry.cxx b/Tracking/Acts/FaserActsGeometry/src/FaserActsWriteTrackingGeometry.cxx index 2ac697d2..0fe031df 100755 --- a/Tracking/Acts/FaserActsGeometry/src/FaserActsWriteTrackingGeometry.cxx +++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsWriteTrackingGeometry.cxx @@ -11,6 +11,7 @@ // PACKAGE #include "FaserActsGeometryInterfaces/IFaserActsTrackingGeometrySvc.h" +#include "FaserActsGeometryInterfaces/IFaserActsMaterialJsonWriterTool.h" #include "FaserActsGeometry/FaserActsGeometryContext.h" #include "ActsInterop/Logger.h" @@ -28,6 +29,7 @@ StatusCode FaserActsWriteTrackingGeometry::initialize() { ATH_CHECK(m_objWriterTool.retrieve()); ATH_CHECK(m_trackingGeometryTool.retrieve()); + ATH_CHECK(m_materialJsonWriterTool.retrieve() ); return StatusCode::SUCCESS; } @@ -40,7 +42,10 @@ StatusCode FaserActsWriteTrackingGeometry::execute(const EventContext& /*ctx*/ ) //const FaserActsGeometryContext& gctx = m_trackingGeometryTool->getGeometryContext(ctx); const FaserActsGeometryContext& gctx = m_trackingGeometryTool->getNominalGeometryContext(); + ATH_MSG_INFO("start to write object"); m_objWriterTool->write(gctx, *trackingGeometry); + ATH_MSG_INFO("start to write geometry"); + m_materialJsonWriterTool->write(*trackingGeometry); return StatusCode::SUCCESS; } diff --git a/Tracking/Acts/FaserActsGeometry/src/FaserActsWriteTrackingGeometry.h b/Tracking/Acts/FaserActsGeometry/src/FaserActsWriteTrackingGeometry.h index 4b877274..c6dc6fb7 100755 --- a/Tracking/Acts/FaserActsGeometry/src/FaserActsWriteTrackingGeometry.h +++ b/Tracking/Acts/FaserActsGeometry/src/FaserActsWriteTrackingGeometry.h @@ -26,6 +26,7 @@ namespace Acts { } class IFaserActsTrackingGeometrySvc; +class IFaserActsMaterialJsonWriterTool; class FaserActsWriteTrackingGeometry : public AthReentrantAlgorithm { public: @@ -39,6 +40,7 @@ private: ToolHandle<IFaserActsTrackingGeometryTool> m_trackingGeometryTool{this, "TrackingGeometryTool", "FaserActsTrackingGeometryTool"}; ToolHandle<FaserActsObjWriterTool> m_objWriterTool{this, "ObjWriterTool", "FaserActsObjWriterTool"}; + ToolHandle<IFaserActsMaterialJsonWriterTool> m_materialJsonWriterTool{this, "MaterialJsonWriterTool", "FaserActsMaterialJsonWriterTool"}; }; diff --git a/Tracking/Acts/FaserActsGeometry/src/components/FaserActsGeometry_entries.cxx b/Tracking/Acts/FaserActsGeometry/src/components/FaserActsGeometry_entries.cxx index 1adfac54..e7b7fb2f 100755 --- a/Tracking/Acts/FaserActsGeometry/src/components/FaserActsGeometry_entries.cxx +++ b/Tracking/Acts/FaserActsGeometry/src/components/FaserActsGeometry_entries.cxx @@ -11,15 +11,11 @@ #include "../FaserActsPropStepRootWriterSvc.h" #include "../FaserActsAlignmentCondAlg.h" #include "../FaserNominalAlignmentCondAlg.h" -//#include "FaserActsGeometry/FaserActsKalmanFilterAlg.h" -//#include "FaserActsGeometry/FaserActsExCellWriterSvc.h" -//#include "FaserActsGeometry/FaserActsMaterialTrackWriterSvc.h" -//#include "FaserActsGeometry/GeomShiftCondAlg.h" -//#include "FaserActsGeometry/FaserActsWriteTrackingGeometryTransforms.h" #include "../FaserActsVolumeMappingTool.h" -//#include "FaserActsGeometry/FaserActsMaterialJsonWriterTool.h" +#include "../FaserActsMaterialJsonWriterTool.h" #include "../FaserActsMaterialMapping.h" #include "../FaserActsSurfaceMappingTool.h" +#include "../FaserActsMaterialTrackWriterSvc.h" DECLARE_COMPONENT( FaserActsTrackingGeometrySvc ) DECLARE_COMPONENT( FaserActsTrackingGeometryTool ) @@ -29,13 +25,8 @@ DECLARE_COMPONENT( FaserActsPropStepRootWriterSvc ) DECLARE_COMPONENT( FaserActsObjWriterTool ) DECLARE_COMPONENT( FaserActsWriteTrackingGeometry ) DECLARE_COMPONENT( FaserActsAlignmentCondAlg ) -//DECLARE_COMPONENT( FaserActsKalmanFilterAlg ) -//DECLARE_COMPONENT( FaserNominalAlignmentCondAlg ) -//DECLARE_COMPONENT( FaserActsExCellWriterSvc ) -//DECLARE_COMPONENT( FaserActsMaterialTrackWriterSvc ) -//DECLARE_COMPONENT( FaserActsWriteTrackingGeometryTransforms ) -//DECLARE_COMPONENT( FaserGeomShiftCondAlg ) DECLARE_COMPONENT( FaserActsVolumeMappingTool ) -//DECLARE_COMPONENT( FaserActsMaterialJsonWriterTool ) +DECLARE_COMPONENT( FaserActsMaterialJsonWriterTool ) +DECLARE_COMPONENT( FaserActsMaterialTrackWriterSvc ) DECLARE_COMPONENT( FaserActsMaterialMapping ) DECLARE_COMPONENT( FaserActsSurfaceMappingTool ) diff --git a/Tracking/Acts/FaserActsGeometry/test/FaserActsWriteTrackingGeometry.py b/Tracking/Acts/FaserActsGeometry/test/FaserActsWriteTrackingGeometry.py index 61f1e624..21dfba7e 100644 --- a/Tracking/Acts/FaserActsGeometry/test/FaserActsWriteTrackingGeometry.py +++ b/Tracking/Acts/FaserActsGeometry/test/FaserActsWriteTrackingGeometry.py @@ -17,16 +17,17 @@ log.setLevel(DEBUG) Configurable.configurableRun3Behavior = True # Configure -ConfigFlags.Input.Files = ["/cvmfs/atlas-nightlies.cern.ch/repo/data/data-art/esd/100evts10lumiblocks.ESD.root"] +ConfigFlags.Input.Files = ["my03.HITS.pool.root"] #ConfigFlags.Output.RDOFileName = "myRDO_sp.pool.root" -ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-01" # Always needed; must match FaserVersion -ConfigFlags.GeoModel.FaserVersion = "FASER-01" # Always needed +ConfigFlags.IOVDb.GlobalTag = "OFLCOND-FASER-02" # Always needed; must match FaserVersion +ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # Always needed +#ConfigFlags.GeoModel.FaserVersion = "FASERNU-03" # Always needed # Workaround for bug/missing flag; unimportant otherwise ConfigFlags.addFlag("Input.InitialTimeStamp", 0) # Workaround to avoid problematic ISF code ConfigFlags.GeoModel.Layout = "Development" ConfigFlags.Detector.GeometryFaserSCT = True -ConfigFlags.GeoModel.AtlasVersion = "ATLAS-R2-2016-01-00-01" # Always needed to fool autoconfig; value ignored +#ConfigFlags.GeoModel.AtlasVersion = "ATLAS-R2-2016-01-00-01" # Always needed to fool autoconfig; value ignored ConfigFlags.GeoModel.Align.Dynamic = False #ConfigFlags.Concurrency.NumThreads = 1 #ConfigFlags.Beam.NumberOfCollisions = 0. @@ -35,7 +36,9 @@ ConfigFlags.lock() # Core components acc = MainServicesCfg(ConfigFlags) acc.merge(PoolReadCfg(ConfigFlags)) -acc.merge(PoolWriteCfg(ConfigFlags)) +#acc.merge(PoolWriteCfg(ConfigFlags)) +from FaserGeoModel.FaserGeoModelConfig import FaserGeometryCfg +acc.merge(FaserGeometryCfg(ConfigFlags)) # Inner Detector acc.merge(FaserActsWriteTrackingGeometryCfg(ConfigFlags)) @@ -61,6 +64,6 @@ acc.getService("ConditionStore").Dump = True acc.printConfig(withDetails=True) ConfigFlags.dump() # Execute and finish -sc = acc.run(maxEvents=-1) +sc = acc.run(maxEvents=1) # Success should be 0 sys.exit(not sc.isSuccess()) -- GitLab