From 55f1a04961964c04d288790e57a101b3d26ae088 Mon Sep 17 00:00:00 2001 From: chizhang Date: Fri, 8 Dec 2023 17:16:45 +0100 Subject: [PATCH 1/6] [PWGDQ] A first version of MillePede2 record creation task for muon alignment --- PWGDQ/Tasks/CMakeLists.txt | 26 ++- PWGDQ/Tasks/mchAlignRecord.cxx | 373 +++++++++++++++++++++++++++++++++ 2 files changed, 391 insertions(+), 8 deletions(-) create mode 100644 PWGDQ/Tasks/mchAlignRecord.cxx diff --git a/PWGDQ/Tasks/CMakeLists.txt b/PWGDQ/Tasks/CMakeLists.txt index 0995bb549e1..f039bc02a13 100644 --- a/PWGDQ/Tasks/CMakeLists.txt +++ b/PWGDQ/Tasks/CMakeLists.txt @@ -14,11 +14,6 @@ o2physics_add_dpl_workflow(table-reader PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore COMPONENT_NAME Analysis) -o2physics_add_dpl_workflow(table-reader-with-assoc - SOURCES tableReader_withAssoc.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2Physics::AnalysisCCDB O2Physics::PWGDQCore - COMPONENT_NAME Analysis) - o2physics_add_dpl_workflow(efficiency SOURCES dqEfficiency.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore @@ -64,7 +59,22 @@ o2physics_add_dpl_workflow(task-j-psi-hf PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore COMPONENT_NAME Analysis) -o2physics_add_dpl_workflow(correlation - SOURCES dqCorrelation.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore O2Physics::GFWCore +o2physics_add_dpl_workflow(task-mch-align-record + SOURCES mchAlignRecord.cxx + PUBLIC_LINK_LIBRARIES + O2::Framework + O2Physics::AnalysisCore + O2Physics::PWGDQCore + O2::CommonUtils + O2::MCHClustering + O2::DPLUtils + O2::CCDB + O2::DataFormatsParameters + O2::MCHBase + O2::MCHTracking + O2::DataFormatsMCH + O2::DetectorsBase + O2::MCHGeometryTransformer + O2::MathUtils + O2::MCHAlign COMPONENT_NAME Analysis) diff --git a/PWGDQ/Tasks/mchAlignRecord.cxx b/PWGDQ/Tasks/mchAlignRecord.cxx new file mode 100644 index 00000000000..08e0fcf0915 --- /dev/null +++ b/PWGDQ/Tasks/mchAlignRecord.cxx @@ -0,0 +1,373 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file mchAlignRecord.cxx +/// \brief A useful task to get MillePede record for alignment and apply different geometries +/// +/// \author Chi ZHANG, CEA-Saclay, chi.zhang@cern.ch + + +#include "Framework/AnalysisTask.h" +#include "Framework/runDataProcessing.h" + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "CommonConstants/LHCConstants.h" +#include "CommonUtils/NameConf.h" +#include "Common/DataModel/EventSelection.h" +#include "PWGDQ/Core/VarManager.h" + +#include "DataFormatsParameters/GRPObject.h" +#include "DataFormatsParameters/GRPMagField.h" +#include "DetectorsBase/GeometryManager.h" +#include "DetectorsBase/GRPGeomHelper.h" +#include "DetectorsBase/Propagator.h" +#include "Framework/Logger.h" +#include "Framework/CallbackService.h" +#include "CCDB/BasicCCDBManager.h" + +#include "MCHGeometryTransformer/Transformations.h" +#include "DataFormatsMCH/Cluster.h" +#include "DataFormatsMCH/TrackMCH.h" +#include "MCHTracking/Track.h" +#include "MCHTracking/TrackExtrap.h" +#include "MCHTracking/TrackParam.h" +#include "MCHTracking/TrackFitter.h" +#include "ReconstructionDataFormats/TrackMCHMID.h" +#include "MCHAlign/Aligner.h" +#include "DetectorsCommonDataFormats/AlignParam.h" +#include "DetectorsCommonDataFormats/DetID.h" +#include "DetectorsCommonDataFormats/DetectorNameConf.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::aod; + +using namespace std; +using std::cout; +using std::endl; + +const int fgNCh = 10; +const int fgNDetElemCh[fgNCh] = {4, 4, 4, 4, 18, 18, 26, 26, 26, 26}; +const int fgSNDetElemCh[fgNCh + 1] = {0, 4, 8, 12, 16, 34, 52, 78, 104, 130, 156}; + +//using MyEvents = soa::Join; +using MyEvents = aod::Collisions; + +struct mchAlignRecordTask { + + Service fCCDB; + parameters::GRPMagField* grpmag; + TGeoManager* geo; + + mch::TrackFitter trackFitter; + mch::Aligner mAlign{}; + Double_t weightRecord{1.0}; + + double Reso_X = 0.0; + double Reso_Y = 0.0; + double ImproveCut = 6.0; + int fCurrentRun; + + map transformOld; + map transformNew; + mch::geo::TransformationCreator transformation; + + Configurable fConfigCcdbUrl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable geoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"}; + Configurable grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; + Configurable fConfigCollType{"collision-type", "pp", "Resolution specification for trackfitter"}; + Configurable fFixChamber{"fix-chamber", "", "Fixing chamber"}; + Configurable fDoNewGeo{"do-realign", false, "Transform to a given new geometry"}; + Configurable fDoEvaluation{"do-evaluation", false, "Enable storage of residuals"}; + Configurable fConfigNewGeoFile{"new-geo", "o2sim_geometry-aligned.root", "New geometry for transformation"}; + Configurable fOutputRecFileName{"outfile-data", "recDataFile.root", "Name of output data record file"}; + Configurable fOutputConsFileName{"outfile-constraint", "recConsFile.root", "Name of output constraint record file"}; + + void init(InitContext& ic) + { + + // Load field and geometry informations here + fCCDB->setURL(fConfigCcdbUrl); + fCCDB->setCaching(true); + fCCDB->setLocalObjectValidityChecking(); + + // Congiguration of resolution for track fitter + if (fConfigCollType.value == "pp") { + Reso_X = 0.4; + Reso_Y = 0.4; + ImproveCut = 6.0; + LOG(info) << "Using pp parameter set for TrackFitter"; + } else { + Reso_X = 0.2; + Reso_Y = 0.2; + ImproveCut = 4.0; + LOG(info) << "Using PbPb parameter set for TrackFitter"; + } + trackFitter.setChamberResolution(Reso_X, Reso_Y); + trackFitter.useChamberResolution(); + + // Configuration for alignment object + mAlign.SetDoEvaluation(fDoEvaluation.value); + mAlign.SetAllowedVariation(0, 2.0); + mAlign.SetAllowedVariation(1, 0.3); + mAlign.SetAllowedVariation(2, 0.002); + mAlign.SetAllowedVariation(3, 2.0); + + // Configuration for chamber fixing + auto chambers = fFixChamber.value; + for (int i = 0; i < chambers.length(); ++i) { + if (chambers[i] == ',') + continue; + int chamber = chambers[i] - '0'; + LOG(info) << Form("%s%d", "Fixing chamber: ", chamber); + mAlign.FixChamber(chamber); + } + + // Init for output saving + mAlign.init(fOutputRecFileName.value, fOutputConsFileName.value); + + ic.services().get().set([this]() { + LOG(info) << "Saving records into ROOT file"; + mAlign.terminate(); + }); + } + + //_________________________________________________________________________________________________ + Int_t GetDetElemId(Int_t iDetElemNumber) + { + // make sure detector number is valid + if (!(iDetElemNumber >= fgSNDetElemCh[0] && + iDetElemNumber < fgSNDetElemCh[fgNCh])) { + LOG(fatal) << "Invalid detector element number: " << iDetElemNumber; + } + /// get det element number from ID + // get chamber and element number in chamber + int iCh = 0; + int iDet = 0; + for (int i = 1; i <= fgNCh; i++) { + if (iDetElemNumber < fgSNDetElemCh[i]) { + iCh = i; + iDet = iDetElemNumber - fgSNDetElemCh[i - 1]; + break; + } + } + + // make sure detector index is valid + if (!(iCh > 0 && iCh <= fgNCh && iDet < fgNDetElemCh[iCh - 1])) { + LOG(fatal) << "Invalid detector element id: " << 100 * iCh + iDet; + } + + // add number of detectors up to this chamber + return 100 * iCh + iDet; + } + + //_________________________________________________________________________________________________ + bool RemoveTrack(mch::Track& track, double ImproveCut) + { + double maxChi2Cluster = 2 * ImproveCut * ImproveCut; + bool removeTrack = false; + + try { + trackFitter.fit(track, false); + } catch (exception const& e) { + removeTrack = true; + return removeTrack; + } + + auto itStartingParam = std::prev(track.rend()); + + while (true) { + try { + trackFitter.fit(track, true, false, (itStartingParam == track.rbegin()) ? nullptr : &itStartingParam); + } catch (exception const&) { + removeTrack = true; + break; + } + + double worstLocalChi2 = -1.0; + + track.tagRemovableClusters(0x1F, false); + auto itWorstParam = track.end(); + + for (auto itParam = track.begin(); itParam != track.end(); ++itParam) { + if (itParam->getLocalChi2() > worstLocalChi2) { + worstLocalChi2 = itParam->getLocalChi2(); + itWorstParam = itParam; + } + } + + if (worstLocalChi2 < maxChi2Cluster) + break; + + if (!itWorstParam->isRemovable()) { + removeTrack = true; + track.removable(); + break; + } + + auto itNextParam = track.removeParamAtCluster(itWorstParam); + auto itNextToNextParam = (itNextParam == track.end()) ? itNextParam : std::next(itNextParam); + itStartingParam = track.rbegin(); + + if (track.getNClusters() < 10) { + removeTrack = true; + break; + } else { + while (itNextToNextParam != track.end()) { + if (itNextToNextParam->getClusterPtr()->getChamberId() != itNextParam->getClusterPtr()->getChamberId()) { + itStartingParam = std::make_reverse_iterator(++itNextParam); + break; + } + ++itNextToNextParam; + } + } + } + + if (!removeTrack) { + for (auto& param : track) { + param.setParameters(param.getSmoothParameters()); + param.setCovariances(param.getSmoothCovariances()); + } + } + + return removeTrack; + } + + template + void runProcessTracks(TEvent const& collision, aod::BCsWithTimestamps const&, TTracks const& tracks, TClusters const& clusters) + { + + auto bc = collision.template bc_as(); + + if (fCurrentRun != bc.runNumber()) { + + grpmag = fCCDB->getForTimeStamp(grpmagPath, bc.timestamp()); + if (grpmag != nullptr) { + o2::base::Propagator::initFieldFromGRP(grpmag); + } + // Configuration for magnetic field and track extrapolation + VarManager::SetupMuonMagField(); + mch::TrackExtrap::useExtrapV2(); + mAlign.SetBFieldOn(mch::TrackExtrap::isFieldON()); + trackFitter.smoothTracks(true); + + // Load reference geometry + geo = fCCDB->getForTimeStamp(geoPath, bc.timestamp()); + transformation = mch::geo::transformationFromTGeoManager(*geo); + for (int i = 0; i < 156; i++) { + int iDEN = GetDetElemId(i); + transformOld[iDEN] = transformation(iDEN); + } + + if (fDoNewGeo.value) { + // Load new geometry with which we want to check + base::GeometryManager::loadGeometry(fConfigNewGeoFile.value); + transformation = mch::geo::transformationFromTGeoManager(*gGeoManager); + for (int i = 0; i < 156; i++) { + int iDEN = GetDetElemId(i); + transformNew[iDEN] = transformation(iDEN); + } + } + + fCurrentRun = bc.runNumber(); + } + + // Loop over forward tracks + for (auto const& track : tracks) { + + int clIndex = -1; + mch::Track convertedTrack; + + // Loop over attached clusters + for (auto const& cluster : clusters) { + + if (cluster.template fwdtrack_as() != track) { + continue; + } + + clIndex += 1; + + mch::Cluster* mch_cluster = new mch::Cluster(); + mch_cluster->x = cluster.x(); + mch_cluster->y = cluster.y(); + mch_cluster->z = cluster.z(); + + if (fDoNewGeo.value) { + math_utils::Point3D local; + math_utils::Point3D master; + + master.SetXYZ(cluster.x(), cluster.y(), cluster.z()); + + transformOld[cluster.deId()].MasterToLocal(master, local); + transformNew[cluster.deId()].LocalToMaster(local, master); + + mch_cluster->x = master.x(); + mch_cluster->y = master.y(); + mch_cluster->z = master.z(); + } + + uint32_t ClUId = mch::Cluster::buildUniqueId(int(cluster.deId() / 100) - 1, cluster.deId(), clIndex); + mch_cluster->uid = ClUId; + + mch_cluster->ex = cluster.isGoodX() ? 0.2 : 10.0; + mch_cluster->ey = cluster.isGoodY() ? 0.2 : 10.0; + + convertedTrack.createParamAtCluster(*mch_cluster); + } + + if (convertedTrack.getNClusters() > 9) { + // Erase removable track + if (!RemoveTrack(convertedTrack, ImproveCut)) { + mAlign.ProcessTrack(convertedTrack, transformation, true, weightRecord); + } + } + } + } + + void processTracks(MyEvents::iterator const& collision, aod::BCsWithTimestamps const& bcs, aod::FwdTracks const& tracks, aod::FwdTrkCls const& clusters) + { + runProcessTracks(collision, bcs, tracks, clusters); + } + + PROCESS_SWITCH(mchAlignRecordTask, processTracks, "Process tracks", true); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +} From 5c99890ef2675797c005e21c1ef7e644aa1ebee5 Mon Sep 17 00:00:00 2001 From: czhang Date: Sun, 10 Mar 2024 23:02:50 +0100 Subject: [PATCH 2/6] Formatting --- PWGDQ/Tasks/CMakeLists.txt | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/PWGDQ/Tasks/CMakeLists.txt b/PWGDQ/Tasks/CMakeLists.txt index f039bc02a13..b87c10fd1d8 100644 --- a/PWGDQ/Tasks/CMakeLists.txt +++ b/PWGDQ/Tasks/CMakeLists.txt @@ -61,10 +61,10 @@ o2physics_add_dpl_workflow(task-j-psi-hf o2physics_add_dpl_workflow(task-mch-align-record SOURCES mchAlignRecord.cxx - PUBLIC_LINK_LIBRARIES - O2::Framework - O2Physics::AnalysisCore - O2Physics::PWGDQCore + PUBLIC_LINK_LIBRARIES + O2::Framework + O2Physics::AnalysisCore + O2Physics::PWGDQCore O2::CommonUtils O2::MCHClustering O2::DPLUtils @@ -76,5 +76,5 @@ o2physics_add_dpl_workflow(task-mch-align-record O2::DetectorsBase O2::MCHGeometryTransformer O2::MathUtils - O2::MCHAlign + O2::MCHAlign COMPONENT_NAME Analysis) From 0347686fe0a0cfc383c30a493fd8eb4d606e0fc9 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Sun, 10 Mar 2024 22:03:23 +0000 Subject: [PATCH 3/6] Please consider the following formatting changes --- PWGDQ/Tasks/mchAlignRecord.cxx | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/PWGDQ/Tasks/mchAlignRecord.cxx b/PWGDQ/Tasks/mchAlignRecord.cxx index 08e0fcf0915..629b61f93b0 100644 --- a/PWGDQ/Tasks/mchAlignRecord.cxx +++ b/PWGDQ/Tasks/mchAlignRecord.cxx @@ -14,7 +14,6 @@ /// /// \author Chi ZHANG, CEA-Saclay, chi.zhang@cern.ch - #include "Framework/AnalysisTask.h" #include "Framework/runDataProcessing.h" @@ -85,7 +84,7 @@ const int fgNCh = 10; const int fgNDetElemCh[fgNCh] = {4, 4, 4, 4, 18, 18, 26, 26, 26, 26}; const int fgSNDetElemCh[fgNCh + 1] = {0, 4, 8, 12, 16, 34, 52, 78, 104, 130, 156}; -//using MyEvents = soa::Join; +// using MyEvents = soa::Join; using MyEvents = aod::Collisions; struct mchAlignRecordTask { From 9b13cf4547931fd4a6941d9b62c912c48f7693a7 Mon Sep 17 00:00:00 2001 From: czhang Date: Sun, 10 Mar 2024 23:15:37 +0100 Subject: [PATCH 4/6] Fixing errors from MegaLinter --- PWGDQ/Tasks/mchAlignRecord.cxx | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/PWGDQ/Tasks/mchAlignRecord.cxx b/PWGDQ/Tasks/mchAlignRecord.cxx index 08e0fcf0915..1180a3438ba 100644 --- a/PWGDQ/Tasks/mchAlignRecord.cxx +++ b/PWGDQ/Tasks/mchAlignRecord.cxx @@ -14,10 +14,6 @@ /// /// \author Chi ZHANG, CEA-Saclay, chi.zhang@cern.ch - -#include "Framework/AnalysisTask.h" -#include "Framework/runDataProcessing.h" - #include #include #include @@ -25,6 +21,9 @@ #include #include +#include "Framework/AnalysisTask.h" +#include "Framework/runDataProcessing.h" + #include #include #include @@ -85,7 +84,7 @@ const int fgNCh = 10; const int fgNDetElemCh[fgNCh] = {4, 4, 4, 4, 18, 18, 26, 26, 26, 26}; const int fgSNDetElemCh[fgNCh + 1] = {0, 4, 8, 12, 16, 34, 52, 78, 104, 130, 156}; -//using MyEvents = soa::Join; +// using MyEvents = soa::Join; using MyEvents = aod::Collisions; struct mchAlignRecordTask { @@ -341,7 +340,7 @@ struct mchAlignRecordTask { mch_cluster->z = master.z(); } - uint32_t ClUId = mch::Cluster::buildUniqueId(int(cluster.deId() / 100) - 1, cluster.deId(), clIndex); + uint32_t ClUId = mch::Cluster::buildUniqueId(static_cast(cluster.deId() / 100) - 1, cluster.deId(), clIndex); mch_cluster->uid = ClUId; mch_cluster->ex = cluster.isGoodX() ? 0.2 : 10.0; From 4e8e4a18f0e14fd5363e65b3adda35c0cadeb1c7 Mon Sep 17 00:00:00 2001 From: czhang Date: Sun, 10 Mar 2024 23:19:03 +0100 Subject: [PATCH 5/6] Fixing MegaLinter errors --- PWGDQ/Tasks/mchAlignRecord.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/PWGDQ/Tasks/mchAlignRecord.cxx b/PWGDQ/Tasks/mchAlignRecord.cxx index 629b61f93b0..1180a3438ba 100644 --- a/PWGDQ/Tasks/mchAlignRecord.cxx +++ b/PWGDQ/Tasks/mchAlignRecord.cxx @@ -14,9 +14,6 @@ /// /// \author Chi ZHANG, CEA-Saclay, chi.zhang@cern.ch -#include "Framework/AnalysisTask.h" -#include "Framework/runDataProcessing.h" - #include #include #include @@ -24,6 +21,9 @@ #include #include +#include "Framework/AnalysisTask.h" +#include "Framework/runDataProcessing.h" + #include #include #include @@ -340,7 +340,7 @@ struct mchAlignRecordTask { mch_cluster->z = master.z(); } - uint32_t ClUId = mch::Cluster::buildUniqueId(int(cluster.deId() / 100) - 1, cluster.deId(), clIndex); + uint32_t ClUId = mch::Cluster::buildUniqueId(static_cast(cluster.deId() / 100) - 1, cluster.deId(), clIndex); mch_cluster->uid = ClUId; mch_cluster->ex = cluster.isGoodX() ? 0.2 : 10.0; From 2cc21c0d685f06f6cafdf42fcfed7312769cb806 Mon Sep 17 00:00:00 2001 From: czhang Date: Sun, 10 Mar 2024 23:28:57 +0100 Subject: [PATCH 6/6] Fixing MegaLinter errors --- PWGDQ/Tasks/mchAlignRecord.cxx | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/PWGDQ/Tasks/mchAlignRecord.cxx b/PWGDQ/Tasks/mchAlignRecord.cxx index 1180a3438ba..cd7730fa86a 100644 --- a/PWGDQ/Tasks/mchAlignRecord.cxx +++ b/PWGDQ/Tasks/mchAlignRecord.cxx @@ -14,17 +14,16 @@ /// /// \author Chi ZHANG, CEA-Saclay, chi.zhang@cern.ch +#include #include #include #include #include #include -#include #include "Framework/AnalysisTask.h" #include "Framework/runDataProcessing.h" -#include #include #include #include