diff --git a/RecFCCeeCalorimeter/src/components/CaloTowerToolFCCee.cpp b/RecFCCeeCalorimeter/src/components/CaloTowerToolFCCee.cpp index ecc5825..10197a4 100644 --- a/RecFCCeeCalorimeter/src/components/CaloTowerToolFCCee.cpp +++ b/RecFCCeeCalorimeter/src/components/CaloTowerToolFCCee.cpp @@ -11,6 +11,12 @@ // DD4hep #include "DD4hep/Detector.h" #include "DD4hep/Readout.h" +#include "DDSegmentation/MultiSegmentation.h" + +// k4geo +#include "detectorSegmentations/FCCSWGridModuleThetaMerged_k4geo.h" +#include "detectorSegmentations/FCCSWGridPhiTheta_k4geo.h" +#include "detectorSegmentations/FCCSWEndcapTurbine_k4geo.h" DECLARE_COMPONENT(CaloTowerToolFCCee) @@ -37,69 +43,56 @@ StatusCode CaloTowerToolFCCee::initialize() { return StatusCode::FAILURE; } - // check if readouts exist & retrieve Module-Theta segmentations - // if readout does not exist, reconstruction without this calorimeter part will be performed - std::pair tmpPair; + std::vector listPhiMax; + std::vector listThetaMax; + listPhiMax.reserve(7); + listThetaMax.reserve(7); - info() << "Retrieving Ecal barrel segmentation" << endmsg; - tmpPair = retrieveSegmentation(m_ecalBarrelReadoutName); - m_ecalBarrelSegmentation = tmpPair.first; - m_ecalBarrelSegmentationType = tmpPair.second; - if (tmpPair.first != nullptr && tmpPair.second == SegmentationType::kWrong) { - error() << "Wrong type of segmentation" << endmsg; - return StatusCode::FAILURE; - } - if (m_useHalfTower) { - m_decoder = m_geoSvc->getDetector()->readout(m_ecalBarrelReadoutName).idSpec().decoder(); + std::pair tmpPair; + if (retrievePhiThetaExtrema(m_ecalBarrelReadoutName, tmpPair) == StatusCode::SUCCESS) { + listPhiMax.push_back(tmpPair.first); + listThetaMax.push_back(tmpPair.second); } - info() << "Retrieving Ecal endcap segmentation" << endmsg; - tmpPair = retrieveSegmentation(m_ecalEndcapReadoutName); - m_ecalEndcapSegmentation = tmpPair.first; - m_ecalEndcapSegmentationType = tmpPair.second; - if (tmpPair.first != nullptr && tmpPair.second == SegmentationType::kWrong) { - error() << "Wrong type of segmentation" << endmsg; - return StatusCode::FAILURE; + if (retrievePhiThetaExtrema(m_ecalEndcapReadoutName, tmpPair) == StatusCode::SUCCESS) { + listPhiMax.push_back(tmpPair.first); + listThetaMax.push_back(tmpPair.second); } - info() << "Retrieving Ecal forward segmentation" << endmsg; - tmpPair = retrieveSegmentation(m_ecalFwdReadoutName); - m_ecalFwdSegmentation = tmpPair.first; - m_ecalFwdSegmentationType = tmpPair.second; - if (tmpPair.first != nullptr && tmpPair.second == SegmentationType::kWrong) { - error() << "Wrong type of segmentation" << endmsg; - return StatusCode::FAILURE; + if (retrievePhiThetaExtrema(m_ecalFwdReadoutName, tmpPair) == StatusCode::SUCCESS) { + listPhiMax.push_back(tmpPair.first); + listThetaMax.push_back(tmpPair.second); } - info() << "Retrieving Hcal barrel segmentation" << endmsg; - tmpPair = retrieveSegmentation(m_hcalBarrelReadoutName); - m_hcalBarrelSegmentation = tmpPair.first; - m_hcalBarrelSegmentationType = tmpPair.second; - if (tmpPair.first != nullptr && tmpPair.second == SegmentationType::kWrong) { - error() << "Wrong type of segmentation" << endmsg; - return StatusCode::FAILURE; + if (retrievePhiThetaExtrema(m_hcalBarrelReadoutName, tmpPair) == StatusCode::SUCCESS) { + listPhiMax.push_back(tmpPair.first); + listThetaMax.push_back(tmpPair.second); } - info() << "Retrieving Hcal extended barrel segmentation" << endmsg; - tmpPair = retrieveSegmentation(m_hcalExtBarrelReadoutName); - m_hcalExtBarrelSegmentation = tmpPair.first; - m_hcalExtBarrelSegmentationType = tmpPair.second; - if (tmpPair.first != nullptr && tmpPair.second == SegmentationType::kWrong) { - error() << "Wrong type of segmentation" << endmsg; - return StatusCode::FAILURE; + if (retrievePhiThetaExtrema(m_hcalExtBarrelReadoutName, tmpPair) == StatusCode::SUCCESS) { + listPhiMax.push_back(tmpPair.first); + listThetaMax.push_back(tmpPair.second); } - info() << "Retrieving Hcal endcap segmentation" << endmsg; - tmpPair = retrieveSegmentation(m_hcalEndcapReadoutName); - m_hcalEndcapSegmentation = tmpPair.first; - m_hcalEndcapSegmentationType = tmpPair.second; - if (tmpPair.first != nullptr && tmpPair.second == SegmentationType::kWrong) { - error() << "Wrong type of segmentation" << endmsg; - return StatusCode::FAILURE; + if (retrievePhiThetaExtrema(m_hcalEndcapReadoutName, tmpPair) == StatusCode::SUCCESS) { + listPhiMax.push_back(tmpPair.first); + listThetaMax.push_back(tmpPair.second); } - info() << "Retrieving Hcal forward segmentation" << endmsg; - tmpPair = retrieveSegmentation(m_hcalFwdReadoutName); - m_hcalFwdSegmentation = tmpPair.first; - m_hcalFwdSegmentationType = tmpPair.second; - if (tmpPair.first != nullptr && tmpPair.second == SegmentationType::kWrong) { - error() << "Wrong type of segmentation" << endmsg; - return StatusCode::FAILURE; + if (retrievePhiThetaExtrema(m_hcalFwdReadoutName, tmpPair) == StatusCode::SUCCESS) { + listPhiMax.push_back(tmpPair.first); + listThetaMax.push_back(tmpPair.second); } + // Maximum theta & phi of the calorimeter system + m_phiMax = *std::max_element(listPhiMax.begin(), listPhiMax.end()); + m_thetaMax = *std::max_element(listThetaMax.begin(), listThetaMax.end()); + debug() << "Detector limits: phiMax " << m_phiMax << " thetaMax " << m_thetaMax << endmsg; + + // very small number (epsilon) substructed from the edges to ensure correct division + float epsilon = 0.0001; + // number of phi bins + m_nPhiTower = ceil(2 * (m_phiMax - epsilon) / m_deltaPhiTower); + // number of theta bins + m_nThetaTower = ceil(2 * (fabs(m_thetaMax - M_PI/2.) - epsilon) / m_deltaThetaTower); + debug() << "Towers: thetaMax " << m_thetaMax << ", deltaThetaTower " << m_deltaThetaTower << ", nThetaTower " << m_nThetaTower + << endmsg; + debug() << "Towers: phiMax " << m_phiMax << ", deltaPhiTower " << m_deltaPhiTower << ", nPhiTower " << m_nPhiTower + << endmsg; + return StatusCode::SUCCESS; } @@ -110,24 +103,40 @@ StatusCode CaloTowerToolFCCee::finalize() { return AlgTool::finalize(); } -std::pair CaloTowerToolFCCee::retrievePhiThetaExtrema(dd4hep::DDSegmentation::Segmentation* aSegmentation, SegmentationType aType) { +StatusCode CaloTowerToolFCCee::retrievePhiThetaExtrema(std::string aReadoutName, std::pair &phiThetaPair) { + double phiMax = -1; double thetaMax = -1; - dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged_k4geo* segmentation = nullptr; - dd4hep::DDSegmentation::FCCSWGridPhiTheta_k4geo* segmentationHCal = nullptr; + + // check if readout exists & retrieve Module-Theta segmentation + // if readout does not exist, reconstruction without this calorimeter part will be performed + + std::pair tmpPair; + tmpPair = retrieveSegmentation(aReadoutName); + dd4hep::DDSegmentation::Segmentation* aSegmentation = tmpPair.first; + SegmentationType aType= tmpPair.second; + if (aSegmentation != nullptr && aType == SegmentationType::kWrong) { + error() << "Wrong type of segmentation" << endmsg; + return StatusCode::FAILURE; + } + if (m_useHalfTower) { + m_decoder = m_geoSvc->getDetector()->readout(aReadoutName).idSpec().decoder(); + } + if (aSegmentation != nullptr) { + switch (aType) { case SegmentationType::kModuleTheta: { - info() << "== Retrieving segmentation " << aSegmentation->name() << endmsg; - segmentation = dynamic_cast(aSegmentation); + + dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged_k4geo* segmentation = dynamic_cast(aSegmentation); phiMax = M_PI - M_PI/segmentation->nModules(); thetaMax = M_PI - fabs(segmentation->offsetTheta()) + segmentation->gridSizeTheta() * .5; break; } - case SegmentationType::kPhiTheta: { + case SegmentationType::kPhiTheta: { info() << "== Retrieving segmentation " << aSegmentation->name() << endmsg; - segmentationHCal = dynamic_cast(aSegmentation); + dd4hep::DDSegmentation::FCCSWGridPhiTheta_k4geo* segmentationHCal = dynamic_cast(aSegmentation); phiMax = M_PI - M_PI/segmentationHCal->phiBins(); thetaMax = M_PI - fabs(segmentationHCal->offsetTheta()) + segmentationHCal->gridSizeTheta() * .5; break; @@ -137,69 +146,39 @@ std::pair CaloTowerToolFCCee::retrievePhiThetaExtrema(dd4hep::DD double theta = -1; info() << "== Retrieving multi segmentation " << aSegmentation->name() << endmsg; for (const auto& subSegm: dynamic_cast(aSegmentation)->subSegmentations()) { - segmentation = dynamic_cast(subSegm.segmentation); + dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged_k4geo* segmentation = dynamic_cast(subSegm.segmentation); phi = M_PI - M_PI/segmentation->nModules(); - theta = M_PI - fabs(segmentation->offsetTheta()) + segmentation->gridSizeTheta() * .5; + theta = M_PI - fabs(segmentation->offsetTheta()) + segmentation->gridSizeTheta() * .5; if (theta > thetaMax) { thetaMax = theta;} if (phi > phiMax) { phiMax = phi;} } break; } + case SegmentationType::kEndcapTurbine: { + dd4hep::DDSegmentation::FCCSWEndcapTurbine_k4geo* ECTsegmentation = dynamic_cast(aSegmentation); + thetaMax = M_PI - fabs(ECTsegmentation->offsetTheta()); + phiMax = M_PI; + break; + } case SegmentationType::kWrong: { info() << "== Retrieving WRONG segmentation" << endmsg; phiMax = -1; thetaMax = -1; - break; + return StatusCode::FAILURE; + } + default: { + error() << " Unsupported segmentation" << endmsg; + phiMax = -1; + thetaMax = -1; + return StatusCode::FAILURE; } } } - return std::make_pair(phiMax, thetaMax); + phiThetaPair = std::make_pair(phiMax, thetaMax); + return StatusCode::SUCCESS; } void CaloTowerToolFCCee::towersNumber(int& nTheta, int& nPhi) { - std::vector listPhiMax; - std::vector listThetaMax; - listPhiMax.reserve(7); - listThetaMax.reserve(7); - - std::pair tmpPair; - tmpPair = retrievePhiThetaExtrema(m_ecalBarrelSegmentation, m_ecalBarrelSegmentationType); - listPhiMax.push_back(tmpPair.first); - listThetaMax.push_back(tmpPair.second); - tmpPair = retrievePhiThetaExtrema(m_ecalEndcapSegmentation, m_ecalEndcapSegmentationType); - listPhiMax.push_back(tmpPair.first); - listThetaMax.push_back(tmpPair.second); - tmpPair = retrievePhiThetaExtrema(m_ecalFwdSegmentation, m_ecalFwdSegmentationType); - listPhiMax.push_back(tmpPair.first); - listThetaMax.push_back(tmpPair.second); - tmpPair = retrievePhiThetaExtrema(m_hcalBarrelSegmentation, m_hcalBarrelSegmentationType); - listPhiMax.push_back(tmpPair.first); - listThetaMax.push_back(tmpPair.second); - tmpPair = retrievePhiThetaExtrema(m_hcalExtBarrelSegmentation, m_hcalExtBarrelSegmentationType); - listPhiMax.push_back(tmpPair.first); - listThetaMax.push_back(tmpPair.second); - tmpPair = retrievePhiThetaExtrema(m_hcalEndcapSegmentation, m_hcalEndcapSegmentationType); - listPhiMax.push_back(tmpPair.first); - listThetaMax.push_back(tmpPair.second); - tmpPair = retrievePhiThetaExtrema(m_hcalFwdSegmentation, m_hcalFwdSegmentationType); - listPhiMax.push_back(tmpPair.first); - listThetaMax.push_back(tmpPair.second); - - // Maximum theta & phi of the calorimeter system - m_phiMax = *std::max_element(listPhiMax.begin(), listPhiMax.end()); - m_thetaMax = *std::max_element(listThetaMax.begin(), listThetaMax.end()); - debug() << "Detector limits: phiMax " << m_phiMax << " thetaMax " << m_thetaMax << endmsg; - - // very small number (epsilon) substructed from the edges to ensure correct division - float epsilon = 0.0001; - // number of phi bins - m_nPhiTower = ceil(2 * (m_phiMax - epsilon) / m_deltaPhiTower); - // number of theta bins - m_nThetaTower = ceil(2 * (fabs(m_thetaMax - M_PI/2.) - epsilon) / m_deltaThetaTower); - debug() << "Towers: thetaMax " << m_thetaMax << ", deltaThetaTower " << m_deltaThetaTower << ", nThetaTower " << m_nThetaTower - << endmsg; - debug() << "Towers: phiMax " << m_phiMax << ", deltaPhiTower " << m_deltaPhiTower << ", nPhiTower " << m_nPhiTower - << endmsg; nTheta = m_nThetaTower; nPhi = m_nPhiTower; @@ -215,45 +194,44 @@ uint CaloTowerToolFCCee::buildTowers(std::vector>& aTowers, b const edm4hep::CalorimeterHitCollection* ecalBarrelCells = m_ecalBarrelCells.get(); debug() << "Input Ecal barrel cell collection size: " << ecalBarrelCells->size() << endmsg; // Loop over a collection of calorimeter cells and build calo towers - if (m_ecalBarrelSegmentation != nullptr) { - CellsIntoTowers(aTowers, ecalBarrelCells, m_ecalBarrelSegmentation, m_ecalBarrelSegmentationType, fillTowersCells); + if (ecalBarrelCells->size() >0) { + CellsIntoTowers(aTowers, ecalBarrelCells, fillTowersCells); totalNumberOfCells += ecalBarrelCells->size(); } -/* // 2. ECAL endcap calorimeter const edm4hep::CalorimeterHitCollection* ecalEndcapCells = m_ecalEndcapCells.get(); debug() << "Input Ecal endcap cell collection size: " << ecalEndcapCells->size() << endmsg; // Loop over a collection of calorimeter cells and build calo towers - if (m_ecalEndcapSegmentation != nullptr) { - CellsIntoTowers(aTowers, ecalEndcapCells, m_ecalEndcapSegmentation, m_ecalEndcapSegmentationType, fillTowersCells); + if (ecalEndcapCells->size() > 0) { + CellsIntoTowers(aTowers, ecalEndcapCells, fillTowersCells); totalNumberOfCells += ecalEndcapCells->size(); } - + /* // 3. ECAL forward calorimeter const edm4hep::CalorimeterHitCollection* ecalFwdCells = m_ecalFwdCells.get(); debug() << "Input Ecal forward cell collection size: " << ecalFwdCells->size() << endmsg; // Loop over a collection of calorimeter cells and build calo towers - if (m_ecalFwdSegmentation != nullptr) { - CellsIntoTowers(aTowers, ecalFwdCells, m_ecalFwdSegmentation, m_ecalFwdSegmentationType, fillTowersCells); + if (ecalFwdCells->size() > 0) { + CellsIntoTowers(aTowers, ecalFwdCells, fillTowersCells); totalNumberOfCells += ecalFwdCells->size(); } -*/ + */ // 4. HCAL barrel const edm4hep::CalorimeterHitCollection* hcalBarrelCells = m_hcalBarrelCells.get(); debug() << "Input hadronic barrel cell collection size: " << hcalBarrelCells->size() << endmsg; // Loop over a collection of calorimeter cells and build calo towers - if (m_hcalBarrelSegmentation != nullptr) { - CellsIntoTowers(aTowers, hcalBarrelCells, m_hcalBarrelSegmentation, m_hcalBarrelSegmentationType, fillTowersCells); + if (hcalBarrelCells->size()>0) { + CellsIntoTowers(aTowers, hcalBarrelCells, fillTowersCells); totalNumberOfCells += hcalBarrelCells->size(); } -/* + /* // 5. HCAL extended barrel const edm4hep::CalorimeterHitCollection* hcalExtBarrelCells = m_hcalExtBarrelCells.get(); debug() << "Input hadronic extended barrel cell collection size: " << hcalExtBarrelCells->size() << endmsg; // Loop over a collection of calorimeter cells and build calo towers - if (m_hcalExtBarrelSegmentation != nullptr) { - CellsIntoTowers(aTowers, hcalExtBarrelCells, m_hcalExtBarrelSegmentation, m_hcalExtBarrelSegmentationType, fillTowersCells); + if (hcalExtBarrelCells->size() >0) { + CellsIntoTowers(aTowers, hcalExtBarrelCells, fillTowersCells); totalNumberOfCells += hcalExtBarrelCells->size(); } @@ -261,8 +239,8 @@ uint CaloTowerToolFCCee::buildTowers(std::vector>& aTowers, b const edm4hep::CalorimeterHitCollection* hcalEndcapCells = m_hcalEndcapCells.get(); debug() << "Input Hcal endcap cell collection size: " << hcalEndcapCells->size() << endmsg; // Loop over a collection of calorimeter cells and build calo towers - if (m_hcalEndcapSegmentation != nullptr) { - CellsIntoTowers(aTowers, hcalEndcapCells, m_hcalEndcapSegmentation, m_hcalEndcapSegmentationType, fillTowersCells); + if (hcalEndcapCells->size() > 0) { + CellsIntoTowers(aTowers, hcalEndcapCells, fillTowersCells); totalNumberOfCells += hcalEndcapCells->size(); } @@ -270,8 +248,8 @@ uint CaloTowerToolFCCee::buildTowers(std::vector>& aTowers, b const edm4hep::CalorimeterHitCollection* hcalFwdCells = m_hcalFwdCells.get(); debug() << "Input Hcal forward cell collection size: " << hcalFwdCells->size() << endmsg; // Loop over a collection of calorimeter cells and build calo towers - if (m_hcalFwdSegmentation != nullptr) { - CellsIntoTowers(aTowers, hcalFwdCells, m_hcalFwdSegmentation, m_hcalFwdSegmentationType, fillTowersCells); + if (hcalFwdCells->size()>0) { + CellsIntoTowers(aTowers, hcalFwdCells, fillTowersCells); totalNumberOfCells += hcalFwdCells->size(); } */ @@ -318,31 +296,15 @@ std::map, std::vector> CaloTowerT // to fill the cell infomation into towers void CaloTowerToolFCCee::CellsIntoTowers(std::vector>& aTowers, const edm4hep::CalorimeterHitCollection* aCells, - dd4hep::DDSegmentation::Segmentation* aSegmentation, SegmentationType aType, bool fillTowersCells) { // Loop over a collection of calorimeter cells and build calo towers // tower index of the borders of the cell int iTheta = 0; int iPhi = 0; - - const dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged_k4geo* segmentation = nullptr; - const dd4hep::DDSegmentation::MultiSegmentation* multisegmentation = nullptr; - const dd4hep::DDSegmentation::FCCSWGridPhiTheta_k4geo* segmentationHCal = nullptr; - - if( aType == SegmentationType::kModuleTheta) { - segmentation = dynamic_cast(aSegmentation); - } else if( aType == SegmentationType::kMulti) { - multisegmentation = dynamic_cast(aSegmentation); - } else if ( aType == SegmentationType::kPhiTheta) { - segmentationHCal = dynamic_cast(aSegmentation); - } + bool pass = true; for (const auto& cell : *aCells) { pass = true; - // if multisegmentation is used - first find out which segmentation to use - if( aType == SegmentationType::kMulti) { - segmentation = dynamic_cast(&multisegmentation->subsegmentation(cell.getCellID())); - } if (m_useHalfTower) { uint layerId = m_decoder->get(cell.getCellID(), "layer"); if ( layerId > m_max_layer ) { @@ -357,15 +319,15 @@ void CaloTowerToolFCCee::CellsIntoTowers(std::vector>& aTower float cellPhi = atan2(cellY, cellX); iTheta = idTheta(cellTheta); iPhi = idPhi(cellPhi); - if (aType == SegmentationType::kPhiTheta) { - aTowers[iTheta][phiNeighbour(iPhi)] += cell.getEnergy() * sin(segmentationHCal->theta(cell.getCellID())); - } else { - aTowers[iTheta][phiNeighbour(iPhi)] += cell.getEnergy() * sin(segmentation->theta(cell.getCellID())); - } + aTowers[iTheta][phiNeighbour(iPhi)] += + cell.getEnergy() * sin(cellTheta); + if (fillTowersCells) { m_cellsInTowers[std::make_pair(iTheta, phiNeighbour(iPhi))].push_back(cell.clone()); - if ( m_cellsInTowers[std::make_pair(iTheta, phiNeighbour(iPhi))].size() > 5 ) - verbose() << "NUM CELLs IN TOWER : " << m_cellsInTowers[std::make_pair(iTheta, phiNeighbour(iPhi))].size() << endmsg; + int ncells = m_cellsInTowers[std::make_pair(iTheta, phiNeighbour(iPhi))].size(); + + if ( ncells > 5) + verbose() << "NUM CELLs IN TOWER : " << ncells << endmsg; } } } @@ -383,23 +345,27 @@ std::pair( m_geoSvc->getDetector()->readout(aReadoutName).segmentation().segmentation()); if (segmentation == nullptr) { - segmentation = dynamic_cast( - m_geoSvc->getDetector()->readout(aReadoutName).segmentation().segmentation()); - if (segmentation == nullptr) { - warning() << "There is no module-theta, phi-theta or multi- segmentation for the readout " << aReadoutName << " defined." << endmsg; - return std::make_pair(nullptr, SegmentationType::kWrong); - } else { - // check if multisegmentation contains only module-theta sub-segmentations - dd4hep::DDSegmentation::Segmentation* subsegmentation = nullptr; - for (const auto& subSegm: dynamic_cast(segmentation)->subSegmentations()) { - subsegmentation = dynamic_cast(subSegm.segmentation); - if (subsegmentation == nullptr) { - warning() << "At least one of the sub-segmentations in MultiSegmentation named " << aReadoutName << " is not a module-theta grid." << endmsg; - return std::make_pair(nullptr, SegmentationType::kWrong); - } - } - return std::make_pair(segmentation, SegmentationType::kMulti); - } + segmentation = dynamic_cast(m_geoSvc->getDetector()->readout(aReadoutName).segmentation().segmentation()); + if (segmentation == nullptr) { + segmentation = dynamic_cast( m_geoSvc->getDetector()->readout(aReadoutName).segmentation().segmentation()); + if (segmentation == nullptr) { + warning() << "There is no module-theta, phi-theta, endcap turbine, or multi- segmentation for the readout " << aReadoutName << " defined." << endmsg; + return std::make_pair(nullptr, SegmentationType::kWrong); + } else { + // check if multisegmentation contains only module-theta sub-segmentations + dd4hep::DDSegmentation::Segmentation* subsegmentation = nullptr; + for (const auto& subSegm: dynamic_cast(segmentation)->subSegmentations()) { + subsegmentation = dynamic_cast(subSegm.segmentation); + if (subsegmentation == nullptr) { + warning() << "At least one of the sub-segmentations in MultiSegmentation named " << aReadoutName << " is not a module-theta grid." << endmsg; + return std::make_pair(nullptr, SegmentationType::kWrong); + } + } + return std::make_pair(segmentation, SegmentationType::kMulti); + } + } else { + return std::make_pair(segmentation, SegmentationType::kEndcapTurbine); + } } else { return std::make_pair(segmentation, SegmentationType::kPhiTheta); } diff --git a/RecFCCeeCalorimeter/src/components/CaloTowerToolFCCee.h b/RecFCCeeCalorimeter/src/components/CaloTowerToolFCCee.h index 0809859..758544d 100644 --- a/RecFCCeeCalorimeter/src/components/CaloTowerToolFCCee.h +++ b/RecFCCeeCalorimeter/src/components/CaloTowerToolFCCee.h @@ -4,17 +4,12 @@ // from Gaudi #include "GaudiKernel/AlgTool.h" -// k4geo -#include "detectorSegmentations/FCCSWGridModuleThetaMerged_k4geo.h" -#include "detectorSegmentations/FCCSWGridPhiTheta_k4geo.h" // k4FWCore #include "k4FWCore/DataHandle.h" #include "k4Interface/ITowerToolThetaModule.h" class IGeoSvc; -// dd4hep -#include "DDSegmentation/MultiSegmentation.h" namespace dd4hep { namespace DDSegmentation { @@ -109,7 +104,7 @@ class CaloTowerToolFCCee : public AlgTool, virtual public ITowerToolThetaModule private: /// Type of the segmentation - enum class SegmentationType {kWrong, kModuleTheta, kMulti, kPhiTheta}; + enum class SegmentationType {kWrong, kModuleTheta, kMulti, kPhiTheta, kEndcapTurbine}; /** Correct way to access the neighbour of the phi tower, taking into account * the full coverage in phi. * Full coverage means that first tower in phi, with ID = 0 is a direct @@ -121,15 +116,18 @@ class CaloTowerToolFCCee : public AlgTool, virtual public ITowerToolThetaModule /** This is where the cell info is filled into towers * @param[in] aTowers Calorimeter towers. * @param[in] aCells Calorimeter cells collection. - * @param[in] aSegmentation Segmentation of the calorimeter + * @param[in] fillTowerCells If true, make a list of the cells in each tower */ void CellsIntoTowers(std::vector>& aTowers, const edm4hep::CalorimeterHitCollection* aCells, - dd4hep::DDSegmentation::Segmentation* aSegmentation, SegmentationType aType, bool fillTowersCells); - /** Check if the readout name exists. If so, it returns the segmentation. + /** Find the maximum phi, theta covered by a readout + * @param[in] aReadoutName Readout name to be checked for maximum phi, theta + * @param[out] phiThetaPair Values of the maximum phi and theta + */ + StatusCode retrievePhiThetaExtrema(std::string aReadoutName, std::pair &phiThetaPair); + /** Check if the readout name exists. If so, it returns the segmentation. * @param[in] aReadoutName Readout name to be retrieved */ - std::pair retrievePhiThetaExtrema(dd4hep::DDSegmentation::Segmentation* aSegmentation, SegmentationType aType); std::pair retrieveSegmentation(std::string aReadoutName); /// Handle for electromagnetic barrel cells (input collection) mutable DataHandle m_ecalBarrelCells{"ecalBarrelCells", Gaudi::DataHandle::Reader, this}; @@ -168,20 +166,6 @@ class CaloTowerToolFCCee : public AlgTool, virtual public ITowerToolThetaModule /// Name of the hcal forward calorimeter readout Gaudi::Property m_hcalFwdReadoutName{this, "hcalFwdReadoutName", "", "name of the hcal fwd readout"}; - /// ModuleTheta segmentation of the electromagnetic barrel (owned by DD4hep) - dd4hep::DDSegmentation::Segmentation* m_ecalBarrelSegmentation; - /// ModuleTheta?? segmentation of the ecal endcap calorimeter (owned by DD4hep) - dd4hep::DDSegmentation::Segmentation* m_ecalEndcapSegmentation; - /// ModuleTheta?? segmentation of the ecal forward calorimeter (owned by DD4hep) - dd4hep::DDSegmentation::Segmentation* m_ecalFwdSegmentation; - /// ModuleTheta?? segmentation of the hadronic barrel (owned by DD4hep) - dd4hep::DDSegmentation::Segmentation* m_hcalBarrelSegmentation; - /// ModuleTheta?? segmentation of the hadronic extended barrel (owned by DD4hep) - dd4hep::DDSegmentation::Segmentation* m_hcalExtBarrelSegmentation; - /// ModuleTheta?? segmentation of the hcal endcap calorimeter (owned by DD4hep) - dd4hep::DDSegmentation::Segmentation* m_hcalEndcapSegmentation; - /// ModuleTheta?? segmentation of the hcal forward calorimeter (owned by DD4hep) - dd4hep::DDSegmentation::Segmentation* m_hcalFwdSegmentation; /// Type of segmentation of the electromagnetic barrel SegmentationType m_ecalBarrelSegmentationType; /// Type of segmentation of the ecal endcap calorimeter @@ -198,6 +182,15 @@ class CaloTowerToolFCCee : public AlgTool, virtual public ITowerToolThetaModule SegmentationType m_hcalFwdSegmentationType; /// decoder: only for barrel dd4hep::DDSegmentation::BitFieldCoder* m_decoder; + // Set of bools to record whether a proper segmentation was found + bool m_ecalBarrelSegmentationOK; + bool m_ecalEndcapSegmentationOK; + bool m_ecalFwdSegmentationOK; + bool m_hcalBarrelSegmentationOK; + bool m_hcalExtBarrelSegmentationOK; + bool m_hcalEndcapSegmentationOK; + bool m_hcalFwdSegmentationOK; + /// Maximum theta of detector float m_thetaMax; /// Maximum phi of the detector diff --git a/RecFCCeeCalorimeter/src/components/CellPositionsECalEndcapTurbineSegTool.cpp b/RecFCCeeCalorimeter/src/components/CellPositionsECalEndcapTurbineSegTool.cpp new file mode 100644 index 0000000..17b99cd --- /dev/null +++ b/RecFCCeeCalorimeter/src/components/CellPositionsECalEndcapTurbineSegTool.cpp @@ -0,0 +1,128 @@ +#include "CellPositionsECalEndcapTurbineSegTool.h" + +// EDM +#include "edm4hep/CalorimeterHitCollection.h" + +#include + +DECLARE_COMPONENT(CellPositionsECalEndcapTurbineSegTool) + +CellPositionsECalEndcapTurbineSegTool::CellPositionsECalEndcapTurbineSegTool(const std::string& type, const std::string& name, + const IInterface* parent) + : AlgTool(type, name, parent) { + declareInterface(this); +} + +StatusCode CellPositionsECalEndcapTurbineSegTool::initialize() { + StatusCode sc = AlgTool::initialize(); + if (sc.isFailure()) return sc; + m_geoSvc = service("GeoSvc"); + if (!m_geoSvc) { + error() << "Unable to locate Geometry service." << endmsg; + return StatusCode::FAILURE; + } + + // get segmentation + m_segmentation = dynamic_cast(m_geoSvc->getDetector()->readout(m_readoutName).segmentation().segmentation()); + if (m_segmentation == nullptr) { + error() << "There is no endcap turbine segmentation!!!!" << endmsg; + return StatusCode::FAILURE; + } + debug() << "Found endcap turbine segmentation" << endmsg; + + // Take readout bitfield decoder from GeoSvc + m_decoder = m_geoSvc->getDetector()->readout(m_readoutName).idSpec().decoder(); + m_volman = m_geoSvc->getDetector()->volumeManager(); + + // check if decoder contains "layer" + std::vector fields; + for (uint itField = 0; itField < m_decoder->size(); itField++) { + fields.push_back((*m_decoder)[itField].name()); + debug() << "In positioning tool, field is " << (*m_decoder)[itField].name() << endmsg; + } + auto iter = std::find(fields.begin(), fields.end(), "layer"); + if (iter == fields.end()) { + error() << "Readout does not contain field: 'layer'" << endmsg; + } + + return sc; +} + +void CellPositionsECalEndcapTurbineSegTool::getPositions(const edm4hep::CalorimeterHitCollection& aCells, + edm4hep::CalorimeterHitCollection& outputColl) { + + debug() << "Input collection size : " << aCells.size() << endmsg; + + // Loop through input cell collection, call xyzPosition method for each cell + // and assign position to cloned hit to be saved in outputColl + for (const auto& cell : aCells) { + auto outSeg = CellPositionsECalEndcapTurbineSegTool::xyzPosition(cell.getCellID()); + auto edmPos = edm4hep::Vector3f(); + edmPos.x = outSeg.x() / dd4hep::mm; + edmPos.y = outSeg.y() / dd4hep::mm; + edmPos.z = outSeg.z() / dd4hep::mm; + + auto positionedHit = cell.clone(); + positionedHit.setPosition(edmPos); + outputColl.push_back(positionedHit); + + // Debug information about cell position + debug() << "Cell energy (GeV) : " << positionedHit.getEnergy() << "\tcellID " << positionedHit.getCellID() << endmsg; + debug() << "Position of cell (mm) : \t" << outSeg.x() / dd4hep::mm << "\t" << outSeg.y() / dd4hep::mm << "\t" + << outSeg.z() / dd4hep::mm << "\n" + << endmsg; + } + debug() << "Output positions collection size: " << outputColl.size() << endmsg; +} + +dd4hep::Position CellPositionsECalEndcapTurbineSegTool::xyzPosition(const uint64_t& aCellId) const { + + dd4hep::Position outSeg; + double radius; + + // find position of volume corresponding to first of group of merged cells + debug() << "cellID: " << aCellId << endmsg; + dd4hep::DDSegmentation::CellID volumeId = aCellId; + debug() << "volumeId: " << volumeId << endmsg; + m_decoder->set(volumeId, "rho", 0); + debug() << "volumeId: " << volumeId << endmsg; + m_decoder->set(volumeId, "z", 0); + debug() << "volumeId: " << volumeId << endmsg; + auto detelement = m_volman.lookupDetElement(volumeId); + const auto& transformMatrix = detelement.nominal().worldTransformation(); + double outGlobal[3]; + double inLocal[] = {0, 0, 0}; + transformMatrix.LocalToMaster(inLocal, outGlobal); + debug() << "Position of volume (mm) : \t" + << outGlobal[0] / dd4hep::mm << "\t" + << outGlobal[1] / dd4hep::mm << "\t" + << outGlobal[2] / dd4hep::mm << endmsg; + + // get R, phi of volume + radius = std::sqrt(std::pow(outGlobal[0], 2) + std::pow(outGlobal[1], 2)); + double phi = std::atan2(outGlobal[1], outGlobal[0]); + debug() << "R (mm), phi of volume : \t" << radius / dd4hep::mm << " , " << phi << endmsg; + + // now get offset in theta and in phi from cell ID (due to theta grid + merging in theta/modules) + // the local position is normalised to r_xy=1 so theta is atan(1/z) + dd4hep::DDSegmentation::Vector3D inSeg = m_segmentation->position(aCellId); + debug() << "Local position of cell (mm) : \t" + << inSeg.x() / dd4hep::mm << "\t" + << inSeg.y() / dd4hep::mm << "\t" + << inSeg.z() / dd4hep::mm << endmsg; + outSeg = inSeg; + debug() << "Position of cell (mm) : \t" << outSeg.x() / dd4hep::mm << "\t" << outSeg.y() / dd4hep::mm << "\t" + << outSeg.z() / dd4hep::mm << "\n" + << endmsg; + + return outSeg; +} + +int CellPositionsECalEndcapTurbineSegTool::layerId(const uint64_t& aCellId) { + int layer; + dd4hep::DDSegmentation::CellID cID = aCellId; + layer = m_decoder->get(cID, "layer"); + return layer; +} + +StatusCode CellPositionsECalEndcapTurbineSegTool::finalize() { return AlgTool::finalize(); } diff --git a/RecFCCeeCalorimeter/src/components/CellPositionsECalEndcapTurbineSegTool.h b/RecFCCeeCalorimeter/src/components/CellPositionsECalEndcapTurbineSegTool.h new file mode 100644 index 0000000..6da21e1 --- /dev/null +++ b/RecFCCeeCalorimeter/src/components/CellPositionsECalEndcapTurbineSegTool.h @@ -0,0 +1,67 @@ +#ifndef RECCALORIMETER_CELLPOSITIONSECALENDCAPTURBINESEGTOOL_H +#define RECCALORIMETER_CELLPOSITIONSECALENDCAPTURBINESEGTOOL_H + +// GAUDI +#include "GaudiKernel/AlgTool.h" +#include "GaudiKernel/ServiceHandle.h" + +// k4FWCore +#include "k4FWCore/DataHandle.h" +#include "k4Interface/IGeoSvc.h" +#include "k4Interface/ICellPositionsTool.h" + +// k4geo +#include "detectorCommon/DetUtils_k4geo.h" +#include "detectorSegmentations/FCCSWEndcapTurbine_k4geo.h" + +// DD4hep +#include "DD4hep/Readout.h" +#include "DD4hep/Volumes.h" +#include "DDSegmentation/Segmentation.h" +#include "TGeoManager.h" + +class IGeoSvc; +namespace DD4hep { +namespace DDSegmentation { +class Segmentation; +} +} + +/** @class CellPositionsECalEndcapTurbineSegTool Reconstruction/RecFCCeeCalorimeter/src/components/CellPositionsECalEndcapTurbineSegTool.h + * CellPositionsECalEndcapTurbineSegTool.h + * + * Tool to determine each Calorimeter cell position. + * + * For the FCCee Endcap ECAL, determined from the placed volumes and the FCCSW endcap turbine segmentation. + * + * @author Erich Varnes + */ + +class CellPositionsECalEndcapTurbineSegTool : public AlgTool, virtual public ICellPositionsTool { +public: + CellPositionsECalEndcapTurbineSegTool(const std::string& type, const std::string& name, const IInterface* parent); + ~CellPositionsECalEndcapTurbineSegTool() = default; + + virtual StatusCode initialize() final; + + virtual StatusCode finalize() final; + + virtual void getPositions(const edm4hep::CalorimeterHitCollection& aCells, edm4hep::CalorimeterHitCollection& outputColl) final; + + virtual dd4hep::Position xyzPosition(const uint64_t& aCellId) const final; + + virtual int layerId(const uint64_t& aCellId) final; + +private: + /// Pointer to the geometry service + SmartIF m_geoSvc; + /// Name of the electromagnetic calorimeter readout + Gaudi::Property m_readoutName{this, "readoutName", "ECalEndcapTurbine"}; + /// Merged module-theta segmentation + dd4hep::DDSegmentation::FCCSWEndcapTurbine_k4geo* m_segmentation; + /// Cellid decoder + dd4hep::DDSegmentation::BitFieldCoder* m_decoder; + /// Volume manager + dd4hep::VolumeManager m_volman; +}; +#endif /* RECCALORIMETER_CELLPOSITIONSECALENDCAPTURBINESEGTOOL_H */ diff --git a/RecFCCeeCalorimeter/tests/options/runEcalBarrel_ReconstructionTopoClusters_crosstalk.py b/RecFCCeeCalorimeter/tests/options/runEcalBarrel_ReconstructionTopoClusters_crosstalk.py new file mode 100644 index 0000000..80056b9 --- /dev/null +++ b/RecFCCeeCalorimeter/tests/options/runEcalBarrel_ReconstructionTopoClusters_crosstalk.py @@ -0,0 +1,445 @@ +# +# IMPORTS +# +from Configurables import ApplicationMgr +#from Configurables import EventCounter +from Configurables import AuditorSvc, ChronoAuditor +# Input/output +from Configurables import k4DataSvc, PodioInput +from Configurables import PodioOutput +# Geometry +from Configurables import GeoSvc +# Create cells +from Configurables import CreateCaloCells +from Configurables import CreateEmptyCaloCellsCollection +# Cell positioning tools +from Configurables import CreateCaloCellPositionsFCCee +from Configurables import CellPositionsECalBarrelModuleThetaSegTool +from Configurables import CellPositionsECalEndcapTurbineSegTool +# Redo segmentation for ECAL +from Configurables import RedoSegmentation +# Change HCAL segmentation +from Configurables import RewriteBitfield +# Apply sampling fraction corrections +from Configurables import CalibrateCaloHitsTool +from Configurables import CalibrateInLayersTool +# Up/down stream correction +from Configurables import CorrectCaloClusters +# SW clustering +from Configurables import CaloTowerToolFCCee +from Configurables import CreateCaloClustersSlidingWindowFCCee +# Topo clustering +from Configurables import CaloTopoClusterInputTool +from Configurables import TopoCaloNeighbours +from Configurables import TopoCaloNoisyCells +from Configurables import CaloTopoClusterFCCee +# Decorate clusters with shower shape parameters +from Configurables import AugmentClustersFCCee +# Read crosstalk map +from Configurables import ReadCaloCrosstalkMap +# Logger +from Gaudi.Configuration import INFO, VERBOSE, DEBUG +# units and physical constants +from GaudiKernel.SystemOfUnits import GeV, tesla, mm +from GaudiKernel.PhysicalConstants import pi, halfpi, twopi +# python libraries +import os +from math import cos, sin, tan + +# +# SETTINGS +# + +# - general settings +# +inputfile = "https://fccsw.web.cern.ch/fccsw/filesForSimDigiReco/ALLEGRO/ALLEGRO_o1_v03/ALLEGRO_sim.root" +# note - this file probably contains the old ecal endcap segmentation so we disable the endcap digitisation later +Nevts = 50 # -1 means all events +dumpGDML = False + +# - what to save in output file +# +# for big productions, save significant space removing hits and cells +# however, hits and cluster cells might be wanted for small productions for detailed event displays +# also, cluster cells are needed for the MVA training +# saveHits = False +# saveCells = False +saveHits = True +saveCells = True +saveClusterCells = True +doCrosstalk = True # switch on/off the crosstalk + +# ECAL barrel parameters for digitisation +samplingFraction=[0.37586625991994105] * 1 + [0.13459486704309379] * 1 + [0.142660085165352] * 1 + [0.14768106642302886] * 1 + [0.15205230356024715] * 1 + [0.15593671843591686] * 1 + [0.15969313426201745] * 1 + [0.16334257010426537] * 1 + [0.16546584993953908] * 1 + [0.16930439771304764] * 1 + [0.1725913708958098] * 1 +upstreamParameters = [[0.025582045561310333, -0.9524128168665387, -53.10089405478649, 1.283851527438571, -295.30650178662637, -284.8945817377308]] +downstreamParameters = [[0.0018280333929494054, 0.004932212590963076, 0.8409676097173655, -1.2676690014715288, 0.005347798049886769, 4.161741293789687]] + +ecalBarrelLayers = len(samplingFraction) +resegmentECalBarrel = False + +# - parameters for clustering +# +doSWClustering = True +doTopoClustering = True + +# calculate cluster energy and barycenter per layer and save it as extra parameters +addShapeParameters = True +ecalBarrelThetaWeights = [-1, 3.0, 3.0, 3.0, 4.25, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0] # to be recalculated for V03, separately for topo and calo clusters... + +# +# ALGORITHMS AND SERVICES SETUP +# + +# Input: load the output of the SIM step +evtsvc = k4DataSvc('EventDataSvc') +evtsvc.input = inputfile +input_reader = PodioInput('InputReader') +podioevent = k4DataSvc("EventDataSvc") + +# Detector geometry +# prefix all xmls with path_to_detector +# if K4GEO is empty, this should use relative path to working directory +geoservice = GeoSvc("GeoSvc") +path_to_detector = os.environ.get("K4GEO", "") +detectors_to_use = [ + 'FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ALLEGRO_o1_v03.xml' +] +geoservice.detectors = [ + os.path.join(path_to_detector, _det) for _det in detectors_to_use +] +geoservice.OutputLevel = INFO + +# GDML dump of detector model +if dumpGDML: + from Configurables import GeoToGdmlDumpSvc + gdmldumpservice = GeoToGdmlDumpSvc("GeoToGdmlDumpSvc") + +# Digitisation (merging hits into cells, EM scale calibration via sampling fractions) + +# - ECAL readouts +ecalBarrelReadoutName = "ECalBarrelModuleThetaMerged" +ecalBarrelReadoutName2 = "ECalBarrelModuleThetaMerged2" +ecalEndcapReadoutName = "ECalEndcapTurbine" + +hcalBarrelReadoutName = "" +hcalBarrelReadoutName2 = "" +hcalEndcapReadoutName = "" + +# - EM scale calibration (sampling fraction) +# * ECAL barrel +calibEcalBarrel = CalibrateInLayersTool("CalibrateECalBarrel", + samplingFraction=samplingFraction, + readoutName=ecalBarrelReadoutName, + layerFieldName="layer") +# * ECAL endcap +calibEcalEndcap = CalibrateCaloHitsTool( + "CalibrateECalEndcap", invSamplingFraction="4.27") # FIXME: to be updated for ddsim + +# Create cells in ECal barrel (needed if one wants to apply cell calibration, +# which is not performed by ddsim) + +# read the crosstalk map +readCrosstalkMap = ReadCaloCrosstalkMap("ReadCrosstalkMap", + fileName="https://fccsw.web.cern.ch/fccsw/filesForSimDigiReco/ALLEGRO/ALLEGRO_o1_v03/xtalk_neighbours_map_ecalB_thetamodulemerged.root", + OutputLevel=INFO) + +# - merge hits into cells according to initial segmentation +ecalBarrelCellsName = "ECalBarrelCells" +createEcalBarrelCells = CreateCaloCells("CreateECalBarrelCells", + doCellCalibration=True, + calibTool=calibEcalBarrel, + crosstalksTool=readCrosstalkMap, + addCrosstalk=doCrosstalk, + addCellNoise=False, + filterCellNoise=False, + addPosition=True, + OutputLevel=INFO, + hits=ecalBarrelReadoutName, + cells=ecalBarrelCellsName) + +# - add to Ecal barrel cells the position information +# (good for physics, all coordinates set properly) +cellPositionEcalBarrelTool = CellPositionsECalBarrelModuleThetaSegTool( + "CellPositionsECalBarrel", + readoutName=ecalBarrelReadoutName, + OutputLevel=INFO +) +ecalBarrelPositionedCellsName = ecalBarrelReadoutName + "Positioned" +createEcalBarrelPositionedCells = CreateCaloCellPositionsFCCee( + "CreateECalBarrelPositionedCells", + OutputLevel=INFO +) +createEcalBarrelPositionedCells.positionsTool = cellPositionEcalBarrelTool +createEcalBarrelPositionedCells.hits.Path = ecalBarrelCellsName +createEcalBarrelPositionedCells.positionedHits.Path = ecalBarrelPositionedCellsName + +# Create cells in ECal endcap (needed if one wants to apply cell calibration, +# which is not performed by ddsim) +#ecalEndcapCellsName = "ECalEndcapCells" +#createEcalEndcapCells = CreateCaloCells("CreateEcalEndcapCaloCells", +# doCellCalibration=True, +# calibTool=calibEcalEndcap, +# addCellNoise=False, +# filterCellNoise=False, +# OutputLevel=INFO, +# hits=ecalEndcapReadoutName, +# cells=ecalEndcapCellsName) + +# Add to Ecal endcap cells the position information +# (good for physics, all coordinates set properly) +#cellPositionEcalEndcapTool = CellPositionsECalEndcapTurbineSegTool( +# "CellPositionsECalEndcap", +# readoutName=ecalEndcapReadoutName, +# OutputLevel=INFO +#) +#ecalEndcapPositionedCellsName = "ECalEndcapPositionedCells" +#createEcalEndcapPositionedCells = CreateCaloCellPositionsFCCee( +# "CreateECalEndcapPositionedCells", +# OutputLevel=INFO +#) +#createEcalEndcapPositionedCells.positionsTool = cellPositionEcalEndcapTool +#createEcalEndcapPositionedCells.hits.Path = ecalEndcapCellsName +#createEcalEndcapPositionedCells.positionedHits.Path = ecalEndcapPositionedCellsName + +hcalBarrelCellsName = "emptyCaloCells" +hcalBarrelPositionedCellsName = "emptyCaloCells" +hcalBarrelCellsName2 = "emptyCaloCells" +hcalBarrelPositionedCellsName2 = "emptyCaloCells" +cellPositionHcalBarrelTool = None +cellPositionHcalBarrelTool2 = None + +# Empty cells for parts of calorimeter not implemented yet +createemptycells = CreateEmptyCaloCellsCollection("CreateEmptyCaloCells") +createemptycells.cells.Path = "emptyCaloCells" + +# Produce sliding window clusters +if doSWClustering: + towers = CaloTowerToolFCCee("towers", + deltaThetaTower=4 * 0.009817477 / 4, deltaPhiTower=2 * 2 * pi / 1536., + ecalBarrelReadoutName=ecalBarrelReadoutName, + #ecalEndcapReadoutName=ecalEndcapReadoutName, + ecalEndcapReadoutName="", + ecalFwdReadoutName="", + hcalBarrelReadoutName=hcalBarrelReadoutName2, + hcalExtBarrelReadoutName="", + hcalEndcapReadoutName="", + hcalFwdReadoutName="", + OutputLevel=INFO) + towers.ecalBarrelCells.Path = ecalBarrelPositionedCellsName + #towers.ecalEndcapCells.Path = ecalEndcapPositionedCellsName + towers.ecalEndcapCells.Path = "emptyCaloCells" + towers.ecalFwdCells.Path = "emptyCaloCells" + towers.hcalBarrelCells.Path = hcalBarrelPositionedCellsName2 + towers.hcalExtBarrelCells.Path = "emptyCaloCells" + towers.hcalEndcapCells.Path = "emptyCaloCells" + towers.hcalFwdCells.Path = "emptyCaloCells" + + # Cluster variables + windT = 9 + windP = 17 + posT = 5 + posP = 11 + dupT = 7 + dupP = 13 + finT = 9 + finP = 17 + # Minimal energy to create a cluster in GeV (FCC-ee detectors have to reconstruct low energy particles) + threshold = 0.040 + + createClusters = CreateCaloClustersSlidingWindowFCCee("CreateClusters", + towerTool=towers, + nThetaWindow=windT, nPhiWindow=windP, + nThetaPosition=posT, nPhiPosition=posP, + nThetaDuplicates=dupT, nPhiDuplicates=dupP, + nThetaFinal=finT, nPhiFinal=finP, + energyThreshold=threshold, + energySharingCorrection=False, + attachCells=True, + OutputLevel=INFO + ) + createClusters.clusters.Path = "CaloClusters" + createClusters.clusterCells.Path = "CaloClusterCells" + + if addShapeParameters: + augmentCaloClusters = AugmentClustersFCCee("augmentCaloClusters", + inClusters=createClusters.clusters.Path, + outClusters="Augmented" + createClusters.clusters.Path, + systemIDs=[4], + systemNames=["EMB"], + numLayers=[ecalBarrelLayers], + readoutNames=[ecalBarrelReadoutName], + layerFieldNames=["layer"], + thetaRecalcWeights=[ecalBarrelThetaWeights], + OutputLevel=INFO + ) + +if doTopoClustering: + # Produce topoclusters (ECAL only or ECAL+HCAL) + createTopoInput = CaloTopoClusterInputTool("CreateTopoInput", + ecalBarrelReadoutName=ecalBarrelReadoutName, + ecalEndcapReadoutName="", + ecalFwdReadoutName="", + hcalBarrelReadoutName=hcalBarrelReadoutName2, + hcalExtBarrelReadoutName="", + hcalEndcapReadoutName="", + hcalFwdReadoutName="", + OutputLevel=INFO) + + createTopoInput.ecalBarrelCells.Path = ecalBarrelPositionedCellsName + createTopoInput.ecalEndcapCells.Path = "emptyCaloCells" + createTopoInput.ecalFwdCells.Path = "emptyCaloCells" + createTopoInput.hcalBarrelCells.Path = hcalBarrelPositionedCellsName2 + createTopoInput.hcalExtBarrelCells.Path = "emptyCaloCells" + createTopoInput.hcalEndcapCells.Path = "emptyCaloCells" + createTopoInput.hcalFwdCells.Path = "emptyCaloCells" + cellPositionHcalBarrelNoSegTool = None + cellPositionHcalExtBarrelTool = None + + neighboursMap = "https://fccsw.web.cern.ch/fccsw/filesForSimDigiReco/ALLEGRO/ALLEGRO_o1_v03/neighbours_map_ecalB_thetamodulemerged.root" + noiseMap = "https://fccsw.web.cern.ch/fccsw/filesForSimDigiReco/ALLEGRO/ALLEGRO_o1_v03/cellNoise_map_electronicsNoiseLevel_ecalB_thetamodulemerged.root" + + readNeighboursMap = TopoCaloNeighbours("ReadNeighboursMap", + fileName=neighboursMap, + OutputLevel=INFO) + + # Noise levels per cell + readNoisyCellsMap = TopoCaloNoisyCells("ReadNoisyCellsMap", + fileName=noiseMap, + OutputLevel=INFO) + + createTopoClusters = CaloTopoClusterFCCee("CreateTopoClusters", + TopoClusterInput=createTopoInput, + # expects neighbours map from cellid->vec < neighbourIds > + neigboursTool=readNeighboursMap, + # tool to get noise level per cellid + noiseTool=readNoisyCellsMap, + # cell positions tools for all sub - systems + positionsECalBarrelTool=cellPositionEcalBarrelTool, + positionsHCalBarrelTool=cellPositionHcalBarrelTool2, + # positionsHCalBarrelNoSegTool=cellPositionHcalBarrelNoSegTool, + # positionsHCalExtBarrelTool=cellPositionHcalExtBarrelTool, + # positionsHCalExtBarrelTool = HCalExtBcells, + # positionsEMECTool = EMECcells, + # positionsHECTool = HECcells, + # positionsEMFwdTool = ECalFwdcells, + # positionsHFwdTool = HCalFwdcells, + noSegmentationHCal=False, + seedSigma=4, + neighbourSigma=2, + lastNeighbourSigma=0, + OutputLevel=INFO) + createTopoClusters.clusters.Path = "CaloTopoClusters" + createTopoClusters.clusterCells.Path = "CaloTopoClusterCells" + + + # Correction below is for EM-only clusters + # Need something different for EM+HCAL + if addShapeParameters: + augmentCaloTopoClusters = AugmentClustersFCCee("augmentCaloTopoClusters", + inClusters=createTopoClusters.clusters.Path, + outClusters="Augmented" + createTopoClusters.clusters.Path, + systemIDs=[4], + systemNames=["EMB"], + numLayers=[ecalBarrelLayers], + readoutNames=[ecalBarrelReadoutName], + layerFieldNames=["layer"], + thetaRecalcWeights=[ecalBarrelThetaWeights], + OutputLevel=INFO) +# Output +out = PodioOutput("out", + OutputLevel=INFO) +out.filename = "ALLEGRO_sim_digi_reco.root" + +# drop the unpositioned ECal barrel cells +out.outputCommands = ["keep *", "drop HCal*", "drop emptyCaloCells", "drop ECalBarrelCells*"] +out.outputCommands.append("drop %s" % ecalBarrelReadoutName) +out.outputCommands.append("drop %s" % ecalBarrelReadoutName2) +out.outputCommands.append("drop ECalBarrelCellsMerged") + +# drop lumi, vertex, DCH, Muons (unless want to keep for event display) +out.outputCommands.append("drop Lumi*") +# out.outputCommands.append("drop Vertex*") +# out.outputCommands.append("drop DriftChamber_simHits*") +out.outputCommands.append("drop MuonTagger*") + +# drop hits/positioned cells/cluster cells if desired +if not saveHits: + out.outputCommands.append("drop %s_contributions" % ecalBarrelReadoutName) + out.outputCommands.append("drop %s_contributions" % ecalBarrelReadoutName2) +if not saveCells: + out.outputCommands.append("drop %s" % ecalBarrelPositionedCellsName) + out.outputCommands.append("drop %s" % ecalBarrelPositionedCellsName2) + +if not saveClusterCells: + out.outputCommands.append("drop Calo*ClusterCells*") + +# CPU information +chra = ChronoAuditor() +audsvc = AuditorSvc() +audsvc.Auditors = [chra] +out.AuditExecute = True + +# Event counter +#event_counter = EventCounter('event_counter') +#event_counter.Frequency = 10 + +# Configure list of external services +ExtSvc = [geoservice, podioevent, audsvc] +if dumpGDML: + ExtSvc += [gdmldumpservice] + +# Setup alg sequence +TopAlg = [ +# event_counter, + input_reader, + createEcalBarrelCells, + createEcalBarrelPositionedCells, +# createEcalEndcapCells, +# createEcalEndcapPositionedCells +] +createEcalBarrelCells.AuditExecute = True +createEcalBarrelPositionedCells.AuditExecute = True +#createEcalEndcapCells.AuditExecute = True + +if resegmentECalBarrel: + TopAlg += [ + resegmentEcalBarrelTool, + createEcalBarrelCells2, + createEcalBarrelPositionedCells2, + ] + resegmentEcalBarrelTool.AuditExecute = True + createEcalBarrelCells2.AuditExecute = True + createEcalBarrelPositionedCells2.AuditExecute = True + +if doSWClustering or doTopoClustering: + TopAlg += [createemptycells] + createemptycells.AuditExecute = True + + if doSWClustering: + TopAlg += [createClusters] + createClusters.AuditExecute = True + + if addShapeParameters: + TopAlg += [augmentCaloClusters] + augmentCaloClusters.AuditExecute = True + + if doTopoClustering: + TopAlg += [createTopoClusters] + createTopoClusters.AuditExecute = True + + if addShapeParameters: + TopAlg += [augmentCaloTopoClusters] + augmentCaloTopoClusters.AuditExecute = True + +TopAlg += [ + out +] + +ApplicationMgr( + TopAlg=TopAlg, + EvtSel='NONE', + EvtMax=Nevts, + ExtSvc=ExtSvc, + StopOnSignal=True, +) + diff --git a/RecFCCeeCalorimeter/tests/options/run_thetamodulemerged.py b/RecFCCeeCalorimeter/tests/options/run_thetamodulemerged.py index fd46261..077272c 100644 --- a/RecFCCeeCalorimeter/tests/options/run_thetamodulemerged.py +++ b/RecFCCeeCalorimeter/tests/options/run_thetamodulemerged.py @@ -1,11 +1,9 @@ from Configurables import ApplicationMgr from Configurables import AuditorSvc, ChronoAuditor from Configurables import PodioOutput -from Configurables import CaloTowerToolFCCee -from Configurables import CreateCaloClustersSlidingWindowFCCee from Configurables import CorrectCaloClusters -from Configurables import CalibrateCaloClusters -from Configurables import AugmentClustersFCCee +from Configurables import CreateCaloClustersSlidingWindow +from Configurables import CaloTowerTool from Configurables import CreateEmptyCaloCellsCollection from Configurables import CreateCaloCellPositionsFCCee from Configurables import CellPositionsECalBarrelModuleThetaSegTool @@ -24,36 +22,20 @@ from Configurables import HepMCToEDMConverter from Configurables import GenAlg from Configurables import FCCDataSvc -from Configurables import RewriteBitfield -from Configurables import ReadCaloCrosstalkMap -from Gaudi.Configuration import INFO +from Configurables import CaloTopoClusterInputTool +from Configurables import TopoCaloNeighbours +from Configurables import TopoCaloNoisyCells +from Configurables import CaloTopoClusterFCCee +from Gaudi.Configuration import * import os -from GaudiKernel.SystemOfUnits import GeV, tesla, mm -from GaudiKernel.PhysicalConstants import pi, halfpi, twopi -from math import cos, sin, tan +from GaudiKernel.SystemOfUnits import GeV, tesla use_pythia = False addNoise = False dumpGDML = False runHCal = False -# for big productions, save significant space removing hits and cells -# however, hits and cluster cells might be wanted for small productions for detailed event displays -# also, cluster cells are needed for the MVA training -saveHits = False -saveCells = False -saveClusterCells = False - -# clustering -doClustering = True -# NOTE: since topoclustering requires root files with noise and neighbours that are -# not in the release, topoclusters are disabled -# cluster energy corrections -# simple parametrisations of up/downstream losses -applyUpDownstreamCorrections = True -# calculate cluster energy and barycenter per layer and save it as extra parameters -addShapeParameters = True # Input for simulations (momentum is expected in GeV!) # Parameters for the particle gun simulations, dummy if use_pythia is set @@ -63,30 +45,23 @@ # (in strips: 0.5625/4=0.14) # Nevts = 20000 -Nevts = 2 +Nevts = 100 # Nevts = 1 # Nevts=1000 - -# particle momentum and direction -# momentum = 100 # in GeV -momentum = 10 # in GeV -# momentum = 10 # in GeV -thetaMin = 45 # degrees -thetaMax = 135 # degrees +# momentum = 100 # in GeV +# momentum = 50 # in GeV +momentum = 10 # in GeV +_pi = 3.14159 +thetaMin = 40 # degrees +thetaMax = 140 # degrees # thetaMin = 89 # thetaMax = 91 -# thetaMin = 90 # degrees -# thetaMax = 90 # degrees -# phiMin = halfpi -# phiMax = halfpi +# thetaMin = 90 # degrees +# thetaMax = 90 # degrees +# phiMin = _pi/2. +# phiMax = _pi/2. phiMin = 0 -phiMax = twopi - -# particle origin -# origR = 1000.0*mm -origR = 0.0 * mm -origTheta = halfpi -origPhi = 0.0 +phiMax = 2 * _pi # particle type: 11 electron, 13 muon, 22 photon, 111 pi0, 211 pi+ pdgCode = 11 @@ -119,6 +94,14 @@ pythia8gentool.printPythiaStatistics = False pythia8gentool.pythiaExtraSettings = [""] genAlg.SignalProvider = pythia8gentool + # to smear the primary vertex position: + # from Configurables import GaussSmearVertex + # smeartool = GaussSmearVertex() + # smeartool.xVertexSigma = 0.5*units.mm + # smeartool.yVertexSigma = 0.5*units.mm + # smeartool.zVertexSigma = 40.0*units.mm + # smeartool.tVertexSigma = 180.0*units.picosecond + # genAlg.VertexSmearingTool = smeartool else: from Configurables import MomentumRangeParticleGun pgun = MomentumRangeParticleGun("ParticleGun") @@ -127,30 +110,12 @@ pgun.MomentumMax = momentum * GeV pgun.PhiMin = phiMin pgun.PhiMax = phiMax - pgun.ThetaMin = thetaMin * pi / 180. - pgun.ThetaMax = thetaMax * pi / 180. + pgun.ThetaMin = thetaMin * _pi / 180. + pgun.ThetaMax = thetaMax * _pi / 180. genAlg.SignalProvider = pgun genAlg.hepmc.Path = "hepmc" -# smear/shift vertex -if origR > 0.0: - origX = origR * cos(origPhi) - origY = origR * sin(origPhi) - origZ = origR / tan(origTheta) - print("Particle gun will be moved to %f %f %f" % (origX, origY, origZ)) - from Configurables import GaussSmearVertex - vertexSmearAndShiftTool = GaussSmearVertex() - vertexSmearAndShiftTool.xVertexSigma = 0. - vertexSmearAndShiftTool.yVertexSigma = 0. - vertexSmearAndShiftTool.tVertexSigma = 0. - vertexSmearAndShiftTool.xVertexMean = origX - vertexSmearAndShiftTool.yVertexMean = origY - vertexSmearAndShiftTool.zVertexMean = origZ - vertexSmearAndShiftTool.tVertexMean = 0. - genAlg.VertexSmearingTool = vertexSmearAndShiftTool - -# hepMC -> EDM converter hepmc_converter = HepMCToEDMConverter() hepmc_converter.hepmc.Path = "hepmc" genParticlesOutputName = "genParticles" @@ -165,7 +130,8 @@ path_to_detector = os.environ.get("K4GEO", "") print(path_to_detector) detectors_to_use = [ - 'FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ALLEGRO_o1_v03.xml' + # 'Detector/DetFCCeeIDEA-LAr/compact/FCCee_DectMaster_thetamodulemerged.xml', + 'FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/ALLEGRO_o1_v02.xml' ] # prefix all xmls with path_to_detector geoservice.detectors = [ @@ -208,7 +174,7 @@ if magneticField == 1: field = SimG4ConstantMagneticFieldTool( "SimG4ConstantMagneticFieldTool", - FieldComponentZ=-2 * tesla, + FieldComponentZ=-2*tesla, FieldOn=True, IntegratorStepper="ClassicalRK4" ) @@ -226,15 +192,13 @@ # ECAL ecalBarrelReadoutName = "ECalBarrelModuleThetaMerged" ecalBarrelReadoutName2 = "ECalBarrelModuleThetaMerged2" -ecalEndcapReadoutName = "ECalEndcapTurbine" +ecalEndcapReadoutName = "ECalEndcapPhiEta" # HCAL if runHCal: hcalBarrelReadoutName = "HCalBarrelReadout" - hcalBarrelReadoutName2 = "BarHCal_Readout_phitheta" hcalEndcapReadoutName = "HCalEndcapReadout" else: hcalBarrelReadoutName = "" - hcalBarrelReadoutName2 = "" hcalEndcapReadoutName = "" # Configure saving of calorimeter hits @@ -291,7 +255,8 @@ # Digitization (Merging hits into cells, EM scale calibration) # EM scale calibration (sampling fraction) calibEcalBarrel = CalibrateInLayersTool("CalibrateECalBarrel", - samplingFraction=[0.3775596654349802] * 1 + [0.13400227700041234] * 1 + [0.14390509963164044] * 1 + [0.14998482026270935] * 1 + [0.15457673722531148] * 1 + [0.15928098152159675] * 1 + [0.1635367867767212] * 1 + [0.16801070646031507] * 1 + [0.1713409944779989] * 1 + [0.17580195406064622] * 1 + [0.17966699467772812] * 1, + samplingFraction=[0.36599110182660616] * 1 + [0.1366222373338866] * 1 + [0.1452035173747207] * 1 + [0.1504319190969367] * 1 + [0.15512713637727382] * 1 + [0.1592916726494782] * 1 + [ + 0.16363478857307595] * 1 + [0.1674697333180323] * 1 + [0.16998205747422343] * 1 + [0.1739146363733975] * 1 + [0.17624609543603845] * 1 + [0.1768613530850488] * 1, readoutName=ecalBarrelReadoutName, layerFieldName="layer") @@ -310,18 +275,11 @@ # (merging several modules and severla theta readout cells). # Add noise at this step if you derived the noise already assuming merged cells -# read the crosstalk map -readCrosstalkMap = ReadCaloCrosstalkMap("ReadCrosstalkMap", - fileName="https://fccsw.web.cern.ch/fccsw/filesForSimDigiReco/ALLEGRO/ALLEGRO_o1_v03/xtalk_neighbours_map_ecalB_thetamodulemerged.root", - OutputLevel=INFO) - # Step 1: merge hits into cells according to initial segmentation ecalBarrelCellsName = "ECalBarrelCells" createEcalBarrelCells = CreateCaloCells("CreateECalBarrelCells", doCellCalibration=True, calibTool=calibEcalBarrel, - crosstalksTool=readCrosstalkMap, - addCrosstalk=False, addCellNoise=False, filterCellNoise=False, addPosition=True, @@ -398,7 +356,7 @@ if runHCal: # Create cells in HCal - # 1 - merge hits into cells with the default readout + # 1. step - merge hits into cells with the default readout hcalBarrelCellsName = "HCalBarrelCells" createHcalBarrelCells = CreateCaloCells("CreateHCalBarrelCells", doCellCalibration=True, @@ -410,14 +368,22 @@ cells=hcalBarrelCellsName, OutputLevel=INFO) - # 2 - attach positions to the cells + # createHcalEndcapCells = CreateCaloCells("CreateHcalEndcapCaloCells", + # doCellCalibration=True, + # calibTool=calibHcalEndcap, + # addCellNoise=False, + # filterCellNoise=False, + # OutputLevel=INFO) + # createHcalEndcapCells.hits.Path="HCalEndcapHits" + # createHcalEndcapCells.cells.Path="HCalEndcapCells" + from Configurables import CellPositionsHCalBarrelPhiThetaSegTool cellPositionHcalBarrelTool = CellPositionsHCalBarrelPhiThetaSegTool( "CellPositionsHCalBarrel", readoutName=hcalBarrelReadoutName, OutputLevel=INFO ) - hcalBarrelPositionedCellsName = "HCalBarrelPositionedCells" + hcalBarrelPositionedCellsName = "HCalBarrelPositionedCells" createHcalBarrelPositionedCells = CreateCaloCellPositionsFCCee( "CreateHcalBarrelPositionedCells", OutputLevel=INFO @@ -425,158 +391,196 @@ createHcalBarrelPositionedCells.positionsTool = cellPositionHcalBarrelTool createHcalBarrelPositionedCells.hits.Path = hcalBarrelCellsName createHcalBarrelPositionedCells.positionedHits.Path = hcalBarrelPositionedCellsName - - # 3 - compute new cellID of cells based on new readout - removing row information - hcalBarrelCellsName2 = "HCalBarrelCells2" - rewriteHCalBarrel = RewriteBitfield("RewriteHCalBarrel", - # old bitfield (readout) - oldReadoutName=hcalBarrelReadoutName, - # specify which fields are going to be deleted - removeIds=["row"], - # new bitfield (readout), with new segmentation - newReadoutName=hcalBarrelReadoutName2, - debugPrint=10, - OutputLevel=INFO) - # clusters are needed, with deposit position and cellID in bits - rewriteHCalBarrel.inhits.Path = hcalBarrelCellsName - rewriteHCalBarrel.outhits.Path = hcalBarrelCellsName2 - - # 4 - attach positions to the new cells - from Configurables import CellPositionsHCalBarrelPhiThetaSegTool - hcalBarrelPositionedCellsName2 = "HCalBarrelPositionedCells2" - cellPositionHcalBarrelTool2 = CellPositionsHCalBarrelPhiThetaSegTool( - "CellPositionsHCalBarrel2", - readoutName=hcalBarrelReadoutName2, - OutputLevel=INFO - ) - createHcalBarrelPositionedCells2 = CreateCaloCellPositionsFCCee( - "CreateHCalBarrelPositionedCells2", - OutputLevel=INFO - ) - createHcalBarrelPositionedCells2.positionsTool = cellPositionHcalBarrelTool2 - createHcalBarrelPositionedCells2.hits.Path = hcalBarrelCellsName2 - createHcalBarrelPositionedCells2.positionedHits.Path = hcalBarrelPositionedCellsName2 - - # createHcalEndcapCells = CreateCaloCells("CreateHcalEndcapCaloCells", - # doCellCalibration=True, - # calibTool=calibHcalEndcap, - # addCellNoise=False, - # filterCellNoise=False, - # OutputLevel=INFO) - # createHcalEndcapCells.hits.Path="HCalEndcapHits" - # createHcalEndcapCells.cells.Path="HCalEndcapCells" - else: hcalBarrelCellsName = "emptyCaloCells" hcalBarrelPositionedCellsName = "emptyCaloCells" - hcalBarrelCellsName2 = "emptyCaloCells" - hcalBarrelPositionedCellsName2 = "emptyCaloCells" cellPositionHcalBarrelTool = None - cellPositionHcalBarrelTool2 = None - + # Empty cells for parts of calorimeter not implemented yet createemptycells = CreateEmptyCaloCellsCollection("CreateEmptyCaloCells") createemptycells.cells.Path = "emptyCaloCells" -# Produce sliding window clusters (ECAL only) -towers = CaloTowerToolFCCee("towers", - deltaThetaTower=4 * 0.009817477/4, deltaPhiTower=2 * 2 * pi / 1536., - ecalBarrelReadoutName=ecalBarrelReadoutName, - ecalEndcapReadoutName=ecalEndcapReadoutName, - ecalFwdReadoutName="", - hcalBarrelReadoutName="", - hcalExtBarrelReadoutName="", - hcalEndcapReadoutName="", - hcalFwdReadoutName="", - OutputLevel=INFO) -towers.ecalBarrelCells.Path = ecalBarrelPositionedCellsName +# Produce sliding window clusters +towers = CaloTowerTool("towers", + deltaEtaTower=0.01, deltaPhiTower=2*_pi/768, + radiusForPosition=2160 + 40 / 2.0, + # ecalBarrelReadoutName = ecalBarrelReadoutNamePhiEta, + ecalBarrelReadoutName=ecalBarrelReadoutName, + ecalEndcapReadoutName=ecalEndcapReadoutName, + ecalFwdReadoutName="", + hcalBarrelReadoutName=hcalBarrelReadoutName, + hcalExtBarrelReadoutName="", + hcalEndcapReadoutName="", + hcalFwdReadoutName="", + OutputLevel=INFO) +towers.ecalBarrelCells.Path = ecalBarrelCellsName towers.ecalEndcapCells.Path = "ECalEndcapCells" towers.ecalFwdCells.Path = "emptyCaloCells" -towers.hcalBarrelCells.Path = "emptyCaloCells" +towers.hcalBarrelCells.Path = hcalBarrelCellsName towers.hcalExtBarrelCells.Path = "emptyCaloCells" towers.hcalEndcapCells.Path = "emptyCaloCells" towers.hcalFwdCells.Path = "emptyCaloCells" # Cluster variables -windT = 9 +windE = 9 windP = 17 -posT = 5 +posE = 5 posP = 11 -dupT = 7 +dupE = 7 dupP = 13 -finT = 9 +finE = 9 finP = 17 # Minimal energy to create a cluster in GeV (FCC-ee detectors have to reconstruct low energy particles) threshold = 0.040 -createClusters = CreateCaloClustersSlidingWindowFCCee("CreateClusters", - towerTool=towers, - nThetaWindow=windT, nPhiWindow=windP, - nThetaPosition=posT, nPhiPosition=posP, - nThetaDuplicates=dupT, nPhiDuplicates=dupP, - nThetaFinal=finT, nPhiFinal=finP, - energyThreshold=threshold, - energySharingCorrection=False, - attachCells=True, - OutputLevel=INFO - ) +createClusters = CreateCaloClustersSlidingWindow("CreateClusters", + towerTool=towers, + nEtaWindow=windE, nPhiWindow=windP, + nEtaPosition=posE, nPhiPosition=posP, + nEtaDuplicates=dupE, nPhiDuplicates=dupP, + nEtaFinal=finE, nPhiFinal=finP, + energyThreshold=threshold, + energySharingCorrection=False, + attachCells=True, + OutputLevel=INFO + ) createClusters.clusters.Path = "CaloClusters" createClusters.clusterCells.Path = "CaloClusterCells" +createEcalBarrelPositionedCaloClusterCells = CreateCaloCellPositionsFCCee( + "ECalBarrelPositionedCaloClusterCells", + OutputLevel=INFO +) +createEcalBarrelPositionedCaloClusterCells.positionsTool = cellPositionEcalBarrelTool +createEcalBarrelPositionedCaloClusterCells.hits.Path = "CaloClusterCells" +createEcalBarrelPositionedCaloClusterCells.positionedHits.Path = "PositionedCaloClusterCells" + correctCaloClusters = CorrectCaloClusters("correctCaloClusters", inClusters=createClusters.clusters.Path, - outClusters="Corrected" + createClusters.clusters.Path, - systemIDs=[4], - numLayers=[11], + outClusters="Corrected"+createClusters.clusters.Path, + numLayers=[0], firstLayerIDs=[0], - lastLayerIDs=[10], + lastLayerIDs=[-1], + # readoutNames = [ecalBarrelReadoutNamePhiEta], readoutNames=[ecalBarrelReadoutName], - # do not split the following line or it will break scripts that update the values of the corrections - upstreamParameters = [[0.025582045561310333, -0.9524128168665387, -53.10089405478649, 1.283851527438571, -295.30650178662637, -284.8945817377308]], + # upstreamParameters = [[0.02729094887360858, -1.378665489864182, -68.40424543618059, 3.6930827214130053, -5528.714729126099, -1630.7911298009794]], + upstreamParameters=[ + [0.02729094887360858, -1.378665489864182, -68.40424543618059, 3.6930827214130053, -5528.714729126099, -1630.7911298009794]], upstreamFormulas=[ ['[0]+[1]/(x-[2])', '[0]+[1]/(x-[2])']], - # do not split the following line or it will break scripts that update the values of the corrections - downstreamParameters = [[0.0018280333929494054, 0.004932212590963076, 0.8409676097173655, -1.2676690014715288, 0.005347798049886769, 4.161741293789687]], + # downstreamParameters = [[-0.0032351643028483354, 0.006597484738888312, 0.8972024981692965, -1.0207168610322181, 0.017878133854084398, 9.108099243443101]], + downstreamParameters=[ + [-0.0032351643028483354, 0.006597484738888312, 0.8972024981692965, -1.0207168610322181, 0.017878133854084398, 9.108099243443101]], downstreamFormulas=[ ['[0]+[1]*x', '[0]+[1]/sqrt(x)', '[0]+[1]/x']], OutputLevel=INFO ) -augmentCaloClusters = AugmentClustersFCCee("augmentCaloClusters", - inClusters=createClusters.clusters.Path, - outClusters="Augmented" + createClusters.clusters.Path, - systemIDs=[4], - systemNames=["EMB"], - numLayers=[11], - readoutNames=[ecalBarrelReadoutName], - layerFieldNames=["layer"], - thetaRecalcWeights=[ - # [-1, 3.0, 3.0, 3.0, 4.25, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0] # -1 : use linear weights - [-1, 3.0, 3.0, 3.0, 4.25, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0] # -1 : use linear weights - ], - OutputLevel=INFO - ) +# TOPO CLUSTERS PRODUCTION +createTopoInput = CaloTopoClusterInputTool("CreateTopoInput", + ecalBarrelReadoutName=ecalBarrelReadoutName, + ecalEndcapReadoutName="", + ecalFwdReadoutName="", + hcalBarrelReadoutName=hcalBarrelReadoutName, + hcalExtBarrelReadoutName="", + hcalEndcapReadoutName="", + hcalFwdReadoutName="", + OutputLevel=INFO) + +createTopoInput.ecalBarrelCells.Path = ecalBarrelPositionedCellsName +# createTopoInput.ecalBarrelCells.Path = "ECalBarrelPositionedCells2" +createTopoInput.ecalEndcapCells.Path = "emptyCaloCells" +createTopoInput.ecalFwdCells.Path = "emptyCaloCells" +createTopoInput.hcalBarrelCells.Path = hcalBarrelPositionedCellsName +createTopoInput.hcalExtBarrelCells.Path = "emptyCaloCells" +createTopoInput.hcalEndcapCells.Path = "emptyCaloCells" +createTopoInput.hcalFwdCells.Path = "emptyCaloCells" +cellPositionHcalBarrelNoSegTool = None +cellPositionHcalExtBarrelTool = None + +readNeighboursMap = TopoCaloNeighbours("ReadNeighboursMap", + fileName=os.environ['FCCBASEDIR'] + + "/LAr_scripts/data/neighbours_map_barrel_thetamodulemerged.root", + #"/LAr_scripts/data/neighbours_map_HCalBarrel.root", + OutputLevel=INFO) +# Noise levels per cell +readNoisyCellsMap = TopoCaloNoisyCells("ReadNoisyCellsMap", + fileName=os.environ['FCCBASEDIR'] + + "/LAr_scripts/data/cellNoise_map_electronicsNoiseLevel_thetamodulemerged.root", + OutputLevel=INFO) + +createTopoClusters = CaloTopoClusterFCCee("CreateTopoClusters", + TopoClusterInput=createTopoInput, + # expects neighbours map from cellid->vec < neighbourIds > + neigboursTool=readNeighboursMap, + # tool to get noise level per cellid + noiseTool=readNoisyCellsMap, + # cell positions tools for all sub - systems + positionsECalBarrelTool=cellPositionEcalBarrelTool, + positionsHCalBarrelTool=cellPositionHcalBarrelTool, + positionsHCalBarrelNoSegTool=cellPositionHcalBarrelNoSegTool, + positionsHCalExtBarrelTool=cellPositionHcalExtBarrelTool, + # positionsHCalExtBarrelTool = HCalExtBcells, + # positionsEMECTool = EMECcells, + # positionsHECTool = HECcells, + # positionsEMFwdTool = ECalFwdcells, + # positionsHFwdTool = HCalFwdcells, + noSegmentationHCal=False, + seedSigma=4, + neighbourSigma=2, + lastNeighbourSigma=0, + OutputLevel=INFO) +createTopoClusters.clusters.Path = "CaloTopoClusters" +createTopoClusters.clusterCells.Path = "CaloTopoClusterCells" + +createEcalBarrelPositionedCaloTopoClusterCells = CreateCaloCellPositionsFCCee( + "ECalBarrelPositionedCaloTopoClusterCells", + OutputLevel=INFO +) +# createEcalBarrelPositionedCaloTopoClusterCells.positionsTool = cellPositionEcalBarrelTool2 +createEcalBarrelPositionedCaloTopoClusterCells.positionsTool = cellPositionEcalBarrelTool +createEcalBarrelPositionedCaloTopoClusterCells.hits.Path = "CaloTopoClusterCells" +createEcalBarrelPositionedCaloTopoClusterCells.positionedHits.Path = "PositionedCaloTopoClusterCells" + +correctCaloTopoClusters = CorrectCaloClusters( + "correctCaloTopoClusters", + inClusters=createTopoClusters.clusters.Path, + outClusters="Corrected"+createTopoClusters.clusters.Path, + numLayers=[0], + firstLayerIDs=[0], + lastLayerIDs=[-1], + # readoutNames = [ecalBarrelReadoutNamePhiEta], + readoutNames=[ecalBarrelReadoutName], + # upstreamParameters = [[0.02729094887360858, -1.378665489864182, -68.40424543618059, 3.6930827214130053, -5528.714729126099, -1630.7911298009794]], + upstreamParameters=[[0.02729094887360858, -1.378665489864182, -68.40424543618059, + 3.6930827214130053, -5528.714729126099, -1630.7911298009794]], + upstreamFormulas=[['[0]+[1]/(x-[2])', '[0]+[1]/(x-[2])']], + # downstreamParameters = [[-0.0032351643028483354, 0.006597484738888312, 0.8972024981692965, -1.0207168610322181, 0.017878133854084398, 9.108099243443101]], + downstreamParameters=[[-0.0032351643028483354, 0.006597484738888312, + 0.8972024981692965, -1.0207168610322181, 0.017878133854084398, 9.108099243443101]], + downstreamFormulas=[['[0]+[1]*x', '[0]+[1]/sqrt(x)', '[0]+[1]/x']], + OutputLevel=INFO +) # Output out = PodioOutput("out", OutputLevel=INFO) +# out.outputCommands = ["keep *"] +# out.outputCommands = ["keep *", "drop ECalBarrelHits", "drop HCal*", "drop ECalBarrelCellsStep*", "drop ECalBarrelPositionedHits", "drop emptyCaloCells", "drop CaloClusterCells"] +# out.outputCommands = ["keep *", "drop ECalBarrelHits", "drop HCal*", "drop ECalBarrelCellsStep*", "drop ECalBarrelPositionedHits", "drop emptyCaloCells", "drop CaloClusterCells", "drop %s"%ecalBarrelCellsName, "drop %s"%createEcalBarrelPositionedCells.positionedHits.Path] +# out.outputCommands = ["keep *", "drop ECalBarrelHits", "drop HCal*", "drop *ells*", "drop ECalBarrelPositionedHits", "drop emptyCaloCells"] +# out.outputCommands = ["keep *", "drop HCal*", "drop ECalBarrel*", "drop emptyCaloCells"] if runHCal: out.outputCommands = ["keep *", "drop emptyCaloCells"] else: out.outputCommands = ["keep *", "drop HCal*", "drop emptyCaloCells"] -if not saveCells: - out.outputCommands.append("drop ECal*Cells*") -if not saveClusterCells: - out.outputCommands.append("drop *ClusterCells*") -if not saveHits: - out.outputCommands.append("drop ECal*Hits*") - -out.filename = "./output_evts_" + str(Nevts) + "_pdg_" + str(pdgCode) + "_" + str(momentum) + "_GeV" + "_ThetaMinMax_" + str(thetaMin) + "_" + str( - thetaMax) + "_PhiMinMax_" + str(phiMin) + "_" + str(phiMax) + "_MagneticField_" + str(magneticField) + "_Noise" + str(addNoise) + ".root" +# out.filename = "root/output_fullCalo_SimAndDigi_withTopoCluster_MagneticField_"+str(magneticField)+"_pMin_"+str(momentum*1000)+"_MeV"+"_ThetaMinMax_"+str(thetaMin)+"_"+str(thetaMax)+"_pdgId_"+str(pdgCode)+"_pythia"+str(use_pythia)+"_Noise"+str(addNoise)+".root" +out.filename = "./root_merge/output_evts_"+str(Nevts)+"_pdg_"+str(pdgCode)+"_"+str(momentum)+"_GeV"+"_ThetaMinMax_"+str(thetaMin)+"_"+str( + thetaMax)+"_PhiMinMax_"+str(phiMin)+"_"+str(phiMax)+"_MagneticField_"+str(magneticField)+"_Noise"+str(addNoise)+".root" # CPU information chra = ChronoAuditor() @@ -589,6 +593,7 @@ createEcalBarrelPositionedCells.AuditExecute = True if runHCal: createHcalBarrelCells.AuditExecute = True +createTopoClusters.AuditExecute = True out.AuditExecute = True ExtSvc = [geoservice, podioevent, geantservice, audsvc] @@ -606,33 +611,20 @@ createEcalBarrelPositionedCells2, createEcalEndcapCells ] - if runHCal: TopAlg += [ createHcalBarrelCells, createHcalBarrelPositionedCells, - rewriteHCalBarrel, - createHcalBarrelPositionedCells2, # createHcalEndcapCells ] - -if doClustering: - TopAlg += [ - createemptycells, - createClusters, - ] - - if applyUpDownstreamCorrections: - TopAlg += [ - correctCaloClusters, - ] - - if addShapeParameters: - TopAlg += [ - augmentCaloClusters, - ] - TopAlg += [ + createemptycells, + # createClusters, + # createEcalBarrelPositionedCaloClusterCells, + # correctCaloClusters, + createTopoClusters, + createEcalBarrelPositionedCaloTopoClusterCells, + correctCaloTopoClusters, out ] diff --git a/RecFCCeeCalorimeter/tests/options/run_thetamodulemerged_SW_benchmarkCorr.py b/RecFCCeeCalorimeter/tests/options/run_thetamodulemerged_SW_benchmarkCorr.py index 979ab1f..9119563 100644 --- a/RecFCCeeCalorimeter/tests/options/run_thetamodulemerged_SW_benchmarkCorr.py +++ b/RecFCCeeCalorimeter/tests/options/run_thetamodulemerged_SW_benchmarkCorr.py @@ -8,6 +8,7 @@ from Configurables import CreateEmptyCaloCellsCollection from Configurables import CreateCaloCellPositionsFCCee from Configurables import CellPositionsECalBarrelModuleThetaSegTool +from Configurables import CellPositionsECalEndcapTurbineSegTool from Configurables import RedoSegmentation from Configurables import CreateCaloCells from Configurables import CalibrateCaloHitsTool @@ -288,8 +289,14 @@ readoutName=ecalBarrelReadoutName, layerFieldName="layer") -calibEcalEndcap = CalibrateCaloHitsTool( - "CalibrateECalEndcap", invSamplingFraction="4.27") +#calibEcalEndcap = CalibrateCaloHitsTool( +# "CalibrateECalEndcap", invSamplingFraction="4.27") +calibEcalEndcap = CalibrateInLayersTool("CalibrateECalEndcap", + samplingFraction = [0.16419] * 1 + [0.192898] * 1 + [0.18783] * 1 + [0.193203] * 1 + [0.193928] * 1 + [0.192286] * 1 + [0.199959] * 1 + [0.200153] * 1 + [0.212635] * 1 + [0.180345] * 1 + [0.18488] * 1 + [0.194762] * 1 + [0.197775] * 1 + [0.200504] * 1 + [0.205555] * 1 + [0.203601] * 1 + [0.210877] * 1 + [0.208376] * 1 + [0.216345] * 1 + [0.201452] * 1 + [0.202134] * 1 + [0.207566] * 1 + [0.208152] * 1 + [0.209889] * 1 + [0.211743] * 1 + [0.213188] * 1 + [0.215864] * 1 + [0.22972] * 1 + [0.192515] * 1 + [0.0103233] * 1, + readoutName=ecalEndcapReadoutName, + layerFieldName="layer") + + if runHCal: calibHcells = CalibrateCaloHitsTool( "CalibrateHCal", invSamplingFraction="30.4") @@ -389,6 +396,23 @@ createEcalEndcapCells.hits.Path = "ECalEndcapHits" createEcalEndcapCells.cells.Path = "ECalEndcapCells" +# Add to Ecal endcap cells the position information +# (good for physics, all coordinates set properly) +# not yet merged!! + cellPositionEcalEndcapTool = CellPositionsECalEndcapTurbineSegTool( + "CellPositionsECalEndcap", + readoutName=ecalEndcapReadoutName, + OutputLevel=INFO +) +ecalEndcapPositionedCellsName = "ECalEndcapPositionedCells" +createEcalEndcapPositionedCells = CreateCaloCellPositionsFCCee( + "CreateECalEndcapPositionedCells", + OutputLevel=INFO +) +createEcalEndcapPositionedCells.positionsTool = cellPositionEcalEndcapTool +createEcalEndcapPositionedCells.hits.Path = ecalEndcapCellsName +createEcalEndcapPositionedCells.positionedHits.Path = ecalEndcapPositionedCellsName + if runHCal: # Create cells in HCal # 1 - merge hits into cells with the default readout @@ -493,7 +517,7 @@ deltaPhiTower=2 * pi / 256., max_layer=25, ecalBarrelReadoutName=ecalBarrelReadoutName, - ecalEndcapReadoutName="", + ecalEndcapReadoutName=ecalEndcapReadoutName, ecalFwdReadoutName="", hcalBarrelReadoutName=hcalBarrelReadoutName2, hcalExtBarrelReadoutName="", @@ -501,7 +525,7 @@ hcalFwdReadoutName="", OutputLevel=INFO) towers.ecalBarrelCells.Path = ecalBarrelPositionedCellsName -towers.ecalEndcapCells.Path = "emptyCaloCells" +towers.ecalEndcapCells.Path = ecalEndcapPositionedCellsName towers.ecalFwdCells.Path = "emptyCaloCells" towers.hcalBarrelCells.Path = hcalBarrelPositionedCellsName2 @@ -610,7 +634,8 @@ resegmentEcalBarrel, createEcalBarrelCells2, createEcalBarrelPositionedCells2, - createEcalEndcapCells + createEcalEndcapCells, + createEcalEndcapPositionedCells, ] if runHCal: