diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/include/Utility.h b/Reconstruction/PFA/Pandora/GaudiPandora/include/Utility.h index b29092bf53c2496149d1e304f2bae19377c113be..ed30bf2f653d59b6cbff735caafc451943bc293d 100644 --- a/Reconstruction/PFA/Pandora/GaudiPandora/include/Utility.h +++ b/Reconstruction/PFA/Pandora/GaudiPandora/include/Utility.h @@ -16,5 +16,7 @@ namespace PanUtil{ int getLayer(float x, float y, float z, std::vector<float>& layers); int getLayer_v1(float x, float y, float z, std::vector<float>& layers); int getStave(float x, float y); + double getFieldFromCompact(); + std::vector<double> getTrackingRegionExtent(); } #endif diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/src/GeometryCreator.cpp b/Reconstruction/PFA/Pandora/GaudiPandora/src/GeometryCreator.cpp index 890c0f8b2bb658aa7c506701385f1c2e901fa2bd..2e996bfcdf35bdd963a5c214e48f08e082bac6d4 100644 --- a/Reconstruction/PFA/Pandora/GaudiPandora/src/GeometryCreator.cpp +++ b/Reconstruction/PFA/Pandora/GaudiPandora/src/GeometryCreator.cpp @@ -100,12 +100,21 @@ void GeometryCreator::SetMandatorySubDetectorParameters(SubDetectorTypeMap &subD const gear::TPCParameters &tpcParameters(_GEAR->getTPCParameters()); trackerParameters.m_subDetectorName = "Tracker"; trackerParameters.m_subDetectorType = pandora::INNER_TRACKER; - trackerParameters.m_innerRCoordinate = tpcParameters.getPadLayout().getPlaneExtent()[0]; + if(m_settings.m_use_dd4hep_geo){ + std::vector<double> tmp_extent = PanUtil::getTrackingRegionExtent(); + trackerParameters.m_innerRCoordinate = tmp_extent[0]; + trackerParameters.m_outerRCoordinate = tmp_extent[1]; + trackerParameters.m_outerZCoordinate = tmp_extent[2]; + std::cout<<"DD m_innerRCoordinate="<<trackerParameters.m_innerRCoordinate.Get()<<",m_outerRCoordinate="<<trackerParameters.m_outerRCoordinate.Get()<<",m_outerZCoordinate="<<trackerParameters.m_outerZCoordinate.Get()<<std::endl; + } + else{ + trackerParameters.m_innerRCoordinate = tpcParameters.getPadLayout().getPlaneExtent()[0]; + trackerParameters.m_outerRCoordinate = tpcParameters.getPadLayout().getPlaneExtent()[1]; + trackerParameters.m_outerZCoordinate = tpcParameters.getMaxDriftLength(); + } trackerParameters.m_innerZCoordinate = 0.f; trackerParameters.m_innerPhiCoordinate = 0.f; trackerParameters.m_innerSymmetryOrder = 0; - trackerParameters.m_outerRCoordinate = tpcParameters.getPadLayout().getPlaneExtent()[1]; - trackerParameters.m_outerZCoordinate = tpcParameters.getMaxDriftLength(); trackerParameters.m_outerPhiCoordinate = 0.f; trackerParameters.m_outerSymmetryOrder = 0; trackerParameters.m_isMirroredInZ = true; diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/src/PandoraPFAlg.cpp b/Reconstruction/PFA/Pandora/GaudiPandora/src/PandoraPFAlg.cpp index 5c5ba411b4c6c3a44649c7d22d7a4f659ec5e5ba..92eff92534fe76494250e43867e63b28049f1deb 100644 --- a/Reconstruction/PFA/Pandora/GaudiPandora/src/PandoraPFAlg.cpp +++ b/Reconstruction/PFA/Pandora/GaudiPandora/src/PandoraPFAlg.cpp @@ -10,6 +10,7 @@ #include <algorithm> #include "gear/BField.h" #include <gear/GEAR.h> +#include "Utility.h" #include "LCContent.h" @@ -41,6 +42,7 @@ PandoraPFAlg::PandoraPFAlg(const std::string& name, ISvcLocator* svcLoc) declareProperty("WriteClusterCollection" , m_ClusterCollection_w, "Handle of the ClusterCollection output collection" ); declareProperty("WriteReconstructedParticleCollection", m_ReconstructedParticleCollection_w, "Handle of the ReconstructedParticleCollection output collection" ); declareProperty("WriteVertexCollection" , m_VertexCollection_w, "Handle of the VertexCollection output collection" ); + declareProperty("WriteMCRecoParticleAssociation" , m_MCRecoParticleAssociation_w, "Handle of the MCRecoParticleAssociation output collection" ); } @@ -70,10 +72,17 @@ void PandoraPFAlg::FinaliseSteeringParameters(ISvcLocator* svcloc) { throw "Failed to find GearSvc ..."; } - gear::GearMgr* _GEAR = iSvc->getGearMgr(); - m_settings.m_innerBField = _GEAR->getBField().at(gear::Vector3D(0., 0., 0.)).z(); - std::cout<<"m_innerBField="<<m_settings.m_innerBField<<std::endl; - m_mcParticleCreatorSettings.m_bField = m_settings.m_innerBField; + if(m_use_dd4hep_geo){ + std::cout<<"DD m_innerBField="<<PanUtil::getFieldFromCompact()<<std::endl; + m_settings.m_innerBField = PanUtil::getFieldFromCompact(); + m_mcParticleCreatorSettings.m_bField = PanUtil::getFieldFromCompact(); + } + else{ + gear::GearMgr* _GEAR = iSvc->getGearMgr(); + m_settings.m_innerBField = _GEAR->getBField().at(gear::Vector3D(0., 0., 0.)).z(); + std::cout<<"m_innerBField="<<m_settings.m_innerBField<<std::endl; + m_mcParticleCreatorSettings.m_bField = m_settings.m_innerBField; + } } diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/src/TrackCreator.cpp b/Reconstruction/PFA/Pandora/GaudiPandora/src/TrackCreator.cpp index 6bf803ec4f07c15445faa1c9d57ce58fa000a4f0..8e51d16c863803f9f1b97b5b7596fa7b189f4a5e 100644 --- a/Reconstruction/PFA/Pandora/GaudiPandora/src/TrackCreator.cpp +++ b/Reconstruction/PFA/Pandora/GaudiPandora/src/TrackCreator.cpp @@ -6,7 +6,6 @@ */ -//#include "UTIL/ILDConf.h" #include "edm4hep/Vertex.h" #include "edm4hep/ReconstructedParticle.h" @@ -41,22 +40,9 @@ TrackCreator::TrackCreator(const Settings &settings, const pandora::Pandora *con m_pPandora(pPandora) { - IGearSvc* iSvc = 0; - StatusCode sc = svcloc->service("GearSvc", iSvc, false); - if ( !sc ) - { - throw "Failed to find GearSvc ..."; - } - _GEAR = iSvc->getGearMgr(); - - m_bField = (_GEAR->getBField().at(gear::Vector3D(0., 0., 0.)).z()); - m_tpcInnerR = (_GEAR->getTPCParameters().getPadLayout().getPlaneExtent()[0]); - m_tpcOuterR = (_GEAR->getTPCParameters().getPadLayout().getPlaneExtent()[1]); - m_tpcMaxRow = (_GEAR->getTPCParameters().getPadLayout().getNRows()); - m_tpcZmax = (_GEAR->getTPCParameters().getMaxDriftLength()); - // fg: FTD description in GEAR has changed ... if(m_settings.m_use_dd4hep_geo){ + m_bField = PanUtil::getFieldFromCompact(); const dd4hep::rec::LayeredCalorimeterData * eCalBarrelExtension= PanUtil::getExtension( ( dd4hep::DetType::CALORIMETER | dd4hep::DetType::ELECTROMAGNETIC | dd4hep::DetType::BARREL), ( dd4hep::DetType::AUXILIARY | dd4hep::DetType::FORWARD ) ); const dd4hep::rec::LayeredCalorimeterData * eCalEndcapExtension= PanUtil::getExtension( ( dd4hep::DetType::CALORIMETER | dd4hep::DetType::ELECTROMAGNETIC | dd4hep::DetType::ENDCAP), @@ -65,40 +51,128 @@ TrackCreator::TrackCreator(const Settings &settings, const pandora::Pandora *con m_eCalBarrelInnerSymmetry = eCalBarrelExtension->inner_symmetry; m_eCalBarrelInnerR = eCalBarrelExtension->extent[0]/dd4hep::mm; m_eCalEndCapInnerZ = eCalEndcapExtension->extent[2]/dd4hep::mm; - } + + dd4hep::Detector & mainDetector = dd4hep::Detector::getInstance(); + try{ + dd4hep::rec::FixedPadSizeTPCData * theExtension = 0; + //Get the TPC, make sure not to get the vertex + const std::vector< dd4hep::DetElement>& tpcDets= dd4hep::DetectorSelector(mainDetector).detectors( ( dd4hep::DetType::TRACKER | dd4hep::DetType::BARREL | dd4hep::DetType::GASEOUS ), dd4hep::DetType::VERTEX) ; + + //There should only be one TPC + if(tpcDets.size()==1) { + theExtension = tpcDets[0].extension<dd4hep::rec::FixedPadSizeTPCData>(); + m_tpcInnerR = theExtension->rMin/dd4hep::mm ; + m_tpcOuterR = theExtension->rMax/dd4hep::mm ; + m_tpcMaxRow = theExtension->maxRow; + m_tpcZmax = theExtension->zHalf/dd4hep::mm ; + } + else{ + m_tpcInnerR = 100 ; + m_tpcOuterR = 1800; + m_tpcMaxRow = 100 ; + m_tpcZmax = 2500; + std::cout<<"TrackCreator WARNING:Does not find TPC parameter from dd4hep and set it to dummy value"<<std::endl; + } + } + catch (std::runtime_error &exception){ + std::cout<<"TrackCreator WARNING:exception during TPC parameter construction:"<<exception.what()<<std::endl; + } + //Instead of gear, loop over a provided list of forward (read: endcap) tracking detectors. + const std::vector< dd4hep::DetElement>& endcapDets = dd4hep::DetectorSelector(mainDetector).detectors( ( dd4hep::DetType::TRACKER | dd4hep::DetType::ENDCAP )) ; + if(endcapDets.size()==0){ + m_ftdInnerRadii.push_back(100); + m_ftdOuterRadii.push_back(200); + m_ftdZPositions.push_back(2000); + m_nFtdLayers = 1; + std::cout<<"TrackCreator WARNING:Does not find forward tracking parameter from dd4hep, so set it to dummy value"<<std::endl; + } + for (std::vector< dd4hep::DetElement>::const_iterator iter = endcapDets.begin(), iterEnd = endcapDets.end();iter != iterEnd; ++iter){ + try{ + dd4hep::rec::ZDiskPetalsData * theExtension = 0; + const dd4hep::DetElement& theDetector = *iter; + theExtension = theDetector.extension<dd4hep::rec::ZDiskPetalsData>(); + unsigned int N = theExtension->layers.size(); + std::cout << " Filling FTD-like parameters from DD4hep for "<< theDetector.name() << "- n layers: " << N<< std::endl; + for(unsigned int i = 0; i < N; ++i){ + dd4hep::rec::ZDiskPetalsData::LayerLayout thisLayer = theExtension->layers[i]; + // Create a disk to represent even number petals front side + //FIXME! VERIFY THAT TIS MAKES SENSE! + m_ftdInnerRadii.push_back(thisLayer.distanceSensitive/dd4hep::mm); + m_ftdOuterRadii.push_back(thisLayer.distanceSensitive/dd4hep::mm+thisLayer.lengthSensitive/dd4hep::mm); + // Take the mean z position of the staggered petals + const double zpos(thisLayer.zPosition/dd4hep::mm); + m_ftdZPositions.push_back(zpos); + std::cout << " layer " << i << " - mean z position = " << zpos << std::endl; + } + m_nFtdLayers = m_ftdZPositions.size() ; + } + catch (std::runtime_error &exception){ + std::cout<<"TrackCreator WARNING: exception during Forward Tracking Disk parameter construction for detector: "<<exception.what()<<std::endl; + } + } + // Calculate etd and set parameters + // fg: make SET and ETD optional - as they might not be in the model ... + //FIXME: THINK OF A UNIVERSAL WAY TO HANDLE EXISTENCE OF ADDITIONAL DETECTORS + std::cout << " ETDLayerZ or SETLayerRadius parameters Not being handled! both will be set to " << std::numeric_limits<float>::quiet_NaN() << std::endl; + m_minEtdZPosition = std::numeric_limits<float>::quiet_NaN(); + m_minSetRadius = std::numeric_limits<float>::quiet_NaN(); + + }// if m_use_dd4hep_geo else{ + IGearSvc* iSvc = 0; + StatusCode sc = svcloc->service("GearSvc", iSvc, false); + if ( !sc ) throw "Failed to find GearSvc ..."; + _GEAR = iSvc->getGearMgr(); + m_bField = (_GEAR->getBField().at(gear::Vector3D(0., 0., 0.)).z()); m_eCalBarrelInnerSymmetry = (_GEAR->getEcalBarrelParameters().getSymmetryOrder()); m_eCalBarrelInnerPhi0 = (_GEAR->getEcalBarrelParameters().getPhi0()); m_eCalBarrelInnerR = (_GEAR->getEcalBarrelParameters().getExtent()[0]); m_eCalEndCapInnerZ = (_GEAR->getEcalEndcapParameters().getExtent()[2]); - } - try - { - m_ftdInnerRadii = _GEAR->getGearParameters("FTD").getDoubleVals("FTDInnerRadius"); - m_ftdOuterRadii = _GEAR->getGearParameters("FTD").getDoubleVals("FTDOuterRadius"); - m_ftdZPositions = _GEAR->getGearParameters("FTD").getDoubleVals("FTDZCoordinate"); - m_nFtdLayers = m_ftdZPositions.size(); - } - catch (gear::UnknownParameterException &) - { - const gear::FTDLayerLayout &ftdLayerLayout(_GEAR->getFTDParameters().getFTDLayerLayout()); - std::cout << " Filling FTD parameters from gear::FTDParameters - n layers: " << ftdLayerLayout.getNLayers() << std::endl; - - for(unsigned int i = 0, N = ftdLayerLayout.getNLayers(); i < N; ++i) - { - // Create a disk to represent even number petals front side - m_ftdInnerRadii.push_back(ftdLayerLayout.getSensitiveRinner(i)); - m_ftdOuterRadii.push_back(ftdLayerLayout.getMaxRadius(i)); - - // Take the mean z position of the staggered petals - const double zpos(ftdLayerLayout.getZposition(i)); - m_ftdZPositions.push_back(zpos); - std::cout << " layer " << i << " - mean z position = " << zpos << std::endl; + m_tpcInnerR = (_GEAR->getTPCParameters().getPadLayout().getPlaneExtent()[0]); + m_tpcOuterR = (_GEAR->getTPCParameters().getPadLayout().getPlaneExtent()[1]); + m_tpcMaxRow = (_GEAR->getTPCParameters().getPadLayout().getNRows()); + m_tpcZmax = (_GEAR->getTPCParameters().getMaxDriftLength()); + try{ + m_ftdInnerRadii = _GEAR->getGearParameters("FTD").getDoubleVals("FTDInnerRadius"); + m_ftdOuterRadii = _GEAR->getGearParameters("FTD").getDoubleVals("FTDOuterRadius"); + m_ftdZPositions = _GEAR->getGearParameters("FTD").getDoubleVals("FTDZCoordinate"); + m_nFtdLayers = m_ftdZPositions.size(); + } + catch (gear::UnknownParameterException &){ + const gear::FTDLayerLayout &ftdLayerLayout(_GEAR->getFTDParameters().getFTDLayerLayout()); + std::cout << " Filling FTD parameters from gear::FTDParameters - n layers: " << ftdLayerLayout.getNLayers() << std::endl; + for(unsigned int i = 0, N = ftdLayerLayout.getNLayers(); i < N; ++i) + { + // Create a disk to represent even number petals front side + m_ftdInnerRadii.push_back(ftdLayerLayout.getSensitiveRinner(i)); + m_ftdOuterRadii.push_back(ftdLayerLayout.getMaxRadius(i)); + // Take the mean z position of the staggered petals + const double zpos(ftdLayerLayout.getZposition(i)); + m_ftdZPositions.push_back(zpos); + std::cout << " layer " << i << " - mean z position = " << zpos << std::endl; + } + m_nFtdLayers = m_ftdZPositions.size() ; + } + // Calculate etd and set parameters + // fg: make SET and ETD optional - as they might not be in the model ... + try{ + const DoubleVector &etdZPositions(_GEAR->getGearParameters("ETD").getDoubleVals("ETDLayerZ")); + const DoubleVector &setInnerRadii(_GEAR->getGearParameters("SET").getDoubleVals("SETLayerRadius")); + if (etdZPositions.empty() || setInnerRadii.empty()) + throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER); + m_minEtdZPosition = *(std::min_element(etdZPositions.begin(), etdZPositions.end())); + m_minSetRadius = *(std::min_element(setInnerRadii.begin(), setInnerRadii.end())); + } + catch(gear::UnknownParameterException &){ + std::cout << "Warnning, ETDLayerZ or SETLayerRadius parameters missing from GEAR parameters!" << std::endl + << " -> both will be set to " << std::numeric_limits<float>::quiet_NaN() << std::endl; + //fg: Set them to NAN, so that they cannot be used to set trackParameters.m_reachesCalorimeter = true; + m_minEtdZPosition = std::numeric_limits<float>::quiet_NaN(); + m_minSetRadius = std::numeric_limits<float>::quiet_NaN(); } - - m_nFtdLayers = m_ftdZPositions.size() ; } + // Check tpc parameters if ((std::fabs(m_tpcZmax) < std::numeric_limits<float>::epsilon()) || (std::fabs(m_tpcInnerR) < std::numeric_limits<float>::epsilon()) || (std::fabs(m_tpcOuterR - m_tpcInnerR) < std::numeric_limits<float>::epsilon())) @@ -125,28 +199,6 @@ TrackCreator::TrackCreator(const Settings &settings, const pandora::Pandora *con m_tanLambdaFtd = m_ftdZPositions[0] / m_ftdOuterRadii[0]; - // Calculate etd and set parameters - // fg: make SET and ETD optional - as they might not be in the model ... - try - { - const DoubleVector &etdZPositions(_GEAR->getGearParameters("ETD").getDoubleVals("ETDLayerZ")); - const DoubleVector &setInnerRadii(_GEAR->getGearParameters("SET").getDoubleVals("SETLayerRadius")); - - if (etdZPositions.empty() || setInnerRadii.empty()) - throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER); - - m_minEtdZPosition = *(std::min_element(etdZPositions.begin(), etdZPositions.end())); - m_minSetRadius = *(std::min_element(setInnerRadii.begin(), setInnerRadii.end())); - } - catch(gear::UnknownParameterException &) - { - std::cout << "Warnning, ETDLayerZ or SETLayerRadius parameters missing from GEAR parameters!" << std::endl - << " -> both will be set to " << std::numeric_limits<float>::quiet_NaN() << std::endl; - - //fg: Set them to NAN, so that they cannot be used to set trackParameters.m_reachesCalorimeter = true; - m_minEtdZPosition = std::numeric_limits<float>::quiet_NaN(); - m_minSetRadius = std::numeric_limits<float>::quiet_NaN(); - } } //------------------------------------------------------------------------------------------------------------------------------------------ @@ -893,7 +945,13 @@ int TrackCreator::GetNTpcHits(const edm4hep::Track *const pTrack) const // trk->subdetectorHitNumbers()[ 2 * ILDDetID::TPC - 2 ] = hitCount ; // ---- use hitsInFit : //return pTrack->getSubdetectorHitNumbers()[ 2 * lcio::ILDDetID::TPC - 1 ]; - return pTrack->getSubDetectorHitNumbers(2 * 4 - 1);// lcio::ILDDetID::TPC=4, still use LCIO code now + if(m_settings.m_use_dd4hep_geo){ + //According to FG: [ 2 * lcio::ILDDetID::TPC - 2 ] is the first number and it is supposed to + //be the number of hits in the fit and this is what should be used ! + // at least for DD4hep/DDSim + return pTrack->getSubDetectorHitNumbers(2 * 4 - 2 );//FIXME + } + else return pTrack->getSubDetectorHitNumbers(2 * 4 - 1);// lcio::ILDDetID::TPC=4, still use LCIO code now } //------------------------------------------------------------------------------------------------------------------------------------------ @@ -906,7 +964,10 @@ int TrackCreator::GetNFtdHits(const edm4hep::Track *const pTrack) const // trk->subdetectorHitNumbers()[ 2 * ILDDetID::TPC - 2 ] = hitCount ; // ---- use hitsInFit : //return pTrack->getSubdetectorHitNumbers()[ 2 * lcio::ILDDetID::FTD - 1 ]; - return pTrack->getSubDetectorHitNumbers( 2 * 3 - 1 );// lcio::ILDDetID::FTD=3 + if(m_settings.m_use_dd4hep_geo){ + return pTrack->getSubDetectorHitNumbers(2 * 3 - 1);//FIXME + } + else return pTrack->getSubDetectorHitNumbers( 2 * 3 - 1 );// lcio::ILDDetID::FTD=3 } //------------------------------------------------------------------------------------------------------------------------------------------ diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/src/Utility.cpp b/Reconstruction/PFA/Pandora/GaudiPandora/src/Utility.cpp index 179cd6e6a31e835ef5c5783f0bb1178346f1073a..02797f01c188f175867e7e1af78ac0fc1e239322 100644 --- a/Reconstruction/PFA/Pandora/GaudiPandora/src/Utility.cpp +++ b/Reconstruction/PFA/Pandora/GaudiPandora/src/Utility.cpp @@ -5,6 +5,36 @@ std::string PanUtil::Convert (float number){ buff<<number; return buff.str(); } + +double PanUtil::getFieldFromCompact(){ + dd4hep::Detector& mainDetector = dd4hep::Detector::getInstance(); + const double position[3]={0,0,0}; // position to calculate magnetic field at (the origin in this case) + double magneticFieldVector[3]={0,0,0}; // initialise object to hold magnetic field + mainDetector.field().magneticField(position,magneticFieldVector); // get the magnetic field vector from DD4hep + return magneticFieldVector[2]/dd4hep::tesla; // z component at (0,0,0) +} + +std::vector<double> PanUtil::getTrackingRegionExtent(){ + + ///Rmin, Rmax, Zmax + std::vector<double> extent; + + extent.reserve(3); + try{ + dd4hep::Detector & mainDetector = dd4hep::Detector::getInstance(); + extent[0]=0.1; ///FIXME! CLIC-specific: Inner radius was set to 0 for SiD-type detectors + extent[1]=mainDetector.constantAsDouble("tracker_region_rmax")/dd4hep::mm; + extent[2]=mainDetector.constantAsDouble("tracker_region_zmax")/dd4hep::mm; + } + catch(std::runtime_error &exception){ + std::cout<<"WARNING, does not find TrackingRegion info from dd4hep, set it to dummy value:"<<exception.what()<<std::endl; + extent[0]=0.1; + extent[1]=1000; + extent[2]=2000; + } + return extent; +} + dd4hep::rec::LayeredCalorimeterData * PanUtil::getExtension(unsigned int includeFlag, unsigned int excludeFlag=0) {