From 77b223381ca9f23c930aae6ddcb92dcbb7da8fad Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Wed, 11 Feb 2026 08:20:27 +0100 Subject: [PATCH 1/5] ITS3: define alignable volumes Signed-off-by: Felix Schlepper --- .../ITS/base/include/ITSBase/GeometryTGeo.h | 2 +- .../ITSMFT/ITS/base/src/GeometryTGeo.cxx | 8 +--- .../ITSMFT/ITS/simulation/src/Detector.cxx | 16 ++----- .../DescriptorInnerBarrelITS3.h | 8 ++-- .../src/DescriptorInnerBarrelITS3.cxx | 44 +++++++++++++++++++ 5 files changed, 54 insertions(+), 24 deletions(-) diff --git a/Detectors/ITSMFT/ITS/base/include/ITSBase/GeometryTGeo.h b/Detectors/ITSMFT/ITS/base/include/ITSBase/GeometryTGeo.h index e236c898851f5..c8ef445e273d3 100644 --- a/Detectors/ITSMFT/ITS/base/include/ITSBase/GeometryTGeo.h +++ b/Detectors/ITSMFT/ITS/base/include/ITSBase/GeometryTGeo.h @@ -314,7 +314,7 @@ class GeometryTGeo : public o2::itsmft::GeometryTGeo static const char* getITS3PixelArrayPattern(int layer) { return Form("%s%d", getITS3PixelArrayPatternRaw(), layer); }; /// sym name of the layer - static const char* composeSymNameITS(bool isITS3 = false); + static const char* composeSymNameITS(); /// sym name of the layer static const char* composeSymNameLayer(int lr, bool isITS3 = false); diff --git a/Detectors/ITSMFT/ITS/base/src/GeometryTGeo.cxx b/Detectors/ITSMFT/ITS/base/src/GeometryTGeo.cxx index 60570b2f204c5..5dc499d05037e 100644 --- a/Detectors/ITSMFT/ITS/base/src/GeometryTGeo.cxx +++ b/Detectors/ITSMFT/ITS/base/src/GeometryTGeo.cxx @@ -290,14 +290,8 @@ bool GeometryTGeo::getChipId(int index, int& lay, int& hba, int& sta, int& hsta, } //__________________________________________________________________________ -const char* GeometryTGeo::composeSymNameITS(bool isITS3) +const char* GeometryTGeo::composeSymNameITS() { - if (isITS3) { -#ifdef ENABLE_UPGRADES - return o2::detectors::DetID(o2::detectors::DetID::IT3).getName(); -#endif - } - return o2::detectors::DetID(o2::detectors::DetID::ITS).getName(); } diff --git a/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx b/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx index 8cfe13097d581..f98bf8ef19cf5 100644 --- a/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx +++ b/Detectors/ITSMFT/ITS/simulation/src/Detector.cxx @@ -1106,7 +1106,7 @@ void Detector::addAlignableVolumes() const TString detName = GetName(); TString path = Form("/cave_1/barrel_1/%s_2", GeometryTGeo::getITSVolPattern()); - TString sname = GeometryTGeo::composeSymNameITS((detName == "IT3")); + TString sname = GeometryTGeo::composeSymNameITS(); LOG(debug) << sname << " <-> " << path; @@ -1119,13 +1119,13 @@ void Detector::addAlignableVolumes() const if (lr < mNumberInnerLayers) { if (detName == "ITS") { ((DescriptorInnerBarrelITS2*)mDescriptorIB.get())->addAlignableVolumesLayer(lr, mWrapperLayerId[lr], path, lastUID); + } else { + ((DescriptorInnerBarrelITS3*)mDescriptorIB.get())->addAlignableVolumesLayer(lr, mWrapperLayerId[lr], path, lastUID); } } else { addAlignableVolumesLayer(lr, path, lastUID); } } - - return; } void Detector::addAlignableVolumesLayer(int lr, TString& parent, Int_t& lastUID) const @@ -1148,8 +1148,6 @@ void Detector::addAlignableVolumesLayer(int lr, TString& parent, Int_t& lastUID) for (Int_t hb = start; hb < nhbarrel; hb++) { addAlignableVolumesHalfBarrel(lr, hb, path, lastUID); } - - return; } void Detector::addAlignableVolumesHalfBarrel(Int_t lr, Int_t hb, TString& parent, Int_t& lastUID) const @@ -1177,8 +1175,6 @@ void Detector::addAlignableVolumesHalfBarrel(Int_t lr, Int_t hb, TString& parent for (int st = 0; st < nstaves; st++) { addAlignableVolumesStave(lr, hb, st, path, lastUID); } - - return; } void Detector::addAlignableVolumesStave(Int_t lr, Int_t hb, Int_t st, TString& parent, Int_t& lastUID) const @@ -1205,8 +1201,6 @@ void Detector::addAlignableVolumesStave(Int_t lr, Int_t hb, Int_t st, TString& p for (Int_t sst = start; sst < nhstave; sst++) { addAlignableVolumesHalfStave(lr, hb, st, sst, path, lastUID); } - - return; } void Detector::addAlignableVolumesHalfStave(Int_t lr, Int_t hb, Int_t st, Int_t hst, TString& parent, Int_t& lastUID) const @@ -1236,8 +1230,6 @@ void Detector::addAlignableVolumesHalfStave(Int_t lr, Int_t hb, Int_t st, Int_t for (Int_t md = start; md < nmodules; md++) { addAlignableVolumesModule(lr, hb, st, hst, md, path, lastUID); } - - return; } void Detector::addAlignableVolumesModule(Int_t lr, Int_t hb, Int_t st, Int_t hst, Int_t md, TString& parent, Int_t& lastUID) const @@ -1266,8 +1258,6 @@ void Detector::addAlignableVolumesModule(Int_t lr, Int_t hb, Int_t st, Int_t hst for (Int_t ic = 0; ic < nchips; ic++) { addAlignableVolumesChip(lr, hb, st, hst, md, ic, path, lastUID); } - - return; } void Detector::addAlignableVolumesChip(Int_t lr, Int_t hb, Int_t st, Int_t hst, Int_t md, Int_t ch, TString& parent, diff --git a/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/DescriptorInnerBarrelITS3.h b/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/DescriptorInnerBarrelITS3.h index 80565df55d154..3e230cee474bd 100644 --- a/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/DescriptorInnerBarrelITS3.h +++ b/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/DescriptorInnerBarrelITS3.h @@ -40,17 +40,19 @@ class DescriptorInnerBarrelITS3 : public o2::its::DescriptorInnerBarrel void createLayer(int idLayer, TGeoVolume* dest); void createServices(TGeoVolume* dest); void configure() {} + void addAlignableVolumesLayer(int idLayer, int wrapperLayerId, TString& parentPath, int& lastUID) const; protected: - int mNumLayers{constants::nLayers}; - // wrapper volume properties static constexpr double mTolerance{1e-3}; static constexpr double mWrapperMinRadiusITS3{constants::radiiInner[0] - mTolerance}; static constexpr double mWrapperMaxRadiusITS3{constants::services::radiusOuter + mTolerance}; - static constexpr double mWrapperZSpanITS3{constants::services::length * 2 + mTolerance}; // z length is divided in half + static constexpr double mWrapperZSpanITS3{(constants::services::length * 2) + mTolerance}; // z length is divided in half private: + void addAlignableVolumesHalfBarrel(int idLayer, int iHB, TString& parentPath, int& lastUID) const; + void addAlignableVolumesChips(int idLayer, int iHalfBarrel, TString& parentPath, int& lastUID) const; + std::array, constants::nLayers> mIBLayers; std::unique_ptr mServices; diff --git a/Detectors/Upgrades/ITS3/simulation/src/DescriptorInnerBarrelITS3.cxx b/Detectors/Upgrades/ITS3/simulation/src/DescriptorInnerBarrelITS3.cxx index 04f244284d5b6..42644fbfe0c38 100644 --- a/Detectors/Upgrades/ITS3/simulation/src/DescriptorInnerBarrelITS3.cxx +++ b/Detectors/Upgrades/ITS3/simulation/src/DescriptorInnerBarrelITS3.cxx @@ -10,6 +10,8 @@ // or submit itself to any jurisdiction. #include "ITS3Simulation/DescriptorInnerBarrelITS3.h" +#include "ITSBase/GeometryTGeo.h" +#include "Framework/Logger.h" using namespace o2::its3; @@ -26,3 +28,45 @@ void DescriptorInnerBarrelITS3::createServices(TGeoVolume* dest) mServices = std::make_unique(); mServices->createCYSSAssembly(dest); } + +void DescriptorInnerBarrelITS3::addAlignableVolumesLayer(int idLayer, int wrapperLayerId, TString& parentPath, int& lastUID) const +{ + TString wrpV = wrapperLayerId != -1 ? Form("%s%d_1", its::GeometryTGeo::getITSWrapVolPattern(), wrapperLayerId) : ""; + TString path = Form("%s/%s/%s%d_0", parentPath.Data(), wrpV.Data(), its::GeometryTGeo::getITS3LayerPattern(), idLayer); + TString sname = its::GeometryTGeo::composeSymNameLayer(idLayer, true); + + for (int iHalfBarrel{0}; iHalfBarrel < 2; ++iHalfBarrel) { + addAlignableVolumesHalfBarrel(idLayer, iHalfBarrel, path, lastUID); + } +} + +void DescriptorInnerBarrelITS3::addAlignableVolumesHalfBarrel(int idLayer, int iHB, TString& parentPath, int& lastUID) const +{ + // for ITS3 smallest alignable volume is the half-barrel (e.g., the carbon-form composite structure with the sensors) + TString path = Form("%s/%s%d_%d", parentPath.Data(), its::GeometryTGeo::getITS3HalfBarrelPattern(), idLayer, iHB); + TString sname = its::GeometryTGeo::composeSymNameHalfBarrel(idLayer, iHB, true); + if (!gGeoManager->SetAlignableEntry(sname.Data(), path.Data())) { + LOG(fatal) << "Unable to set alignable entry ! " << sname << " : " << path; + } + addAlignableVolumesChips(idLayer, iHB, path, lastUID); +} + +void DescriptorInnerBarrelITS3::addAlignableVolumesChips(int idLayer, int iHB, TString& parentPath, int& lastUID) const +{ + for (int seg{0}; seg < constants::nSegments[idLayer]; ++seg) { + for (int rsu{0}; rsu < constants::segment::nRSUs; ++rsu) { + for (int tile{0}; tile < constants::rsu::nTiles; ++tile) { + TString path = parentPath; + path += Form("/%s_0/", its::GeometryTGeo::getITS3ChipPattern(idLayer)); + path += Form("%s_%d/", its::GeometryTGeo::getITS3SegmentPattern(idLayer), seg); + path += Form("%s_%d/", its::GeometryTGeo::getITS3RSUPattern(idLayer), rsu); + path += Form("%s_%d/", its::GeometryTGeo::getITS3TilePattern(idLayer), tile); + TString sname = its::GeometryTGeo::composeSymNameChip(idLayer, iHB, 0, seg, rsu, tile, true); + if (!gGeoManager->SetAlignableEntry(sname.Data(), path.Data())) { + LOG(fatal) << "Unable to set alignable entry ! " << sname << " : " << path; + } + ++lastUID; + } + } + } +} From f8e23e0e2b42a382084453f8fe9c41560b60050b Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Wed, 11 Feb 2026 08:21:02 +0100 Subject: [PATCH 2/5] ITS3: load chip response functions from ccdb Signed-off-by: Felix Schlepper --- .../ITS3/macros/test/CheckChipResponseFile.C | 23 ++++++++----------- 1 file changed, 10 insertions(+), 13 deletions(-) diff --git a/Detectors/Upgrades/ITS3/macros/test/CheckChipResponseFile.C b/Detectors/Upgrades/ITS3/macros/test/CheckChipResponseFile.C index 32d5bad87ce21..5bc053c516079 100644 --- a/Detectors/Upgrades/ITS3/macros/test/CheckChipResponseFile.C +++ b/Detectors/Upgrades/ITS3/macros/test/CheckChipResponseFile.C @@ -22,6 +22,7 @@ #include #include +#include "CCDB/BasicCCDBManager.h" #define ENABLE_UPGRADES #include "ITSMFTSimulation/AlpideSimResponse.h" #include "ITS3Simulation/ChipSimResponse.h" @@ -37,16 +38,12 @@ double cm2um(double cm) { return cm * 1e+4; } std::unique_ptr mAlpSimResp0, mAlpSimResp1, mAptSimResp1; -std::unique_ptr loadResponse(const std::string& fileName, const std::string& respName) +std::unique_ptr loadResponse(const std::string& path) { - TFile* f = TFile::Open(fileName.data()); - if (!f) { - std::cerr << fileName << " not found" << std::endl; - return nullptr; - } - auto base = f->Get(respName.c_str()); + auto& cdb = o2::ccdb::BasicCCDBManager::instance(); + o2::itsmft::AlpideSimResponse* base = cdb.get(path); if (!base) { - std::cerr << respName << " not found in " << fileName << std::endl; + std::cerr << path << " not found in " << '\n'; return nullptr; } return std::make_unique(base); @@ -54,24 +51,24 @@ std::unique_ptr loadResponse(const std::string& fileN void LoadRespFunc() { - std::string AptsFile = "$(O2_ROOT)/share/Detectors/Upgrades/ITS3/data/ITS3ChipResponseData/APTSResponseData.root"; - std::string AlpideFile = "$(O2_ROOT)/share/Detectors/ITSMFT/data/AlpideResponseData/AlpideResponseData.root"; + auto& cdb = o2::ccdb::BasicCCDBManager::instance(); + cdb.setURL("https://alice-ccdb.cern.ch/"); std::cout << "=====================\n"; LOGP(info, "ALPIDE Vbb=0V response"); - mAlpSimResp0 = loadResponse(AlpideFile, "response0"); // Vbb=0V + mAlpSimResp0 = loadResponse("ITSMFT/Calib/ALPIDEResponseVbb0"); // Vbb=0V mAlpSimResp0->computeCentreFromData(); mAlpSimResp0->print(); LOGP(info, "Response Centre {}", mAlpSimResp0->getRespCentreDep()); std::cout << "=====================\n"; LOGP(info, "ALPIDE Vbb=-3V response"); - mAlpSimResp1 = loadResponse(AlpideFile, "response1"); // Vbb=-3V + mAlpSimResp1 = loadResponse("ITSMFT/Calib/ALPIDEResponseVbbM3"); // Vbb=-3V mAlpSimResp1->computeCentreFromData(); mAlpSimResp1->print(); LOGP(info, "Response Centre {}", mAlpSimResp1->getRespCentreDep()); std::cout << "=====================\n"; LOGP(info, "APTS response"); - mAptSimResp1 = loadResponse(AptsFile, "response1"); // APTS + mAptSimResp1 = loadResponse("IT3/Calib/APTSResponse"); // APTS mAptSimResp1->computeCentreFromData(); mAptSimResp1->print(); LOGP(info, "Response Centre {}", mAptSimResp1->getRespCentreDep()); From 22759e0919af4b9cadf58a9340e4c6f3c2049a9a Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Tue, 3 Feb 2026 12:39:10 +0100 Subject: [PATCH 3/5] ITS3: improve stepping speed & trk frame Signed-off-by: Felix Schlepper --- .../ITS3/base/include/ITS3Base/SpecsV2.h | 2 +- .../ITS3/reconstruction/src/IOUtils.cxx | 37 ++++++---- .../include/ITS3Simulation/ITS3Layer.h | 5 +- .../ITS3/simulation/src/ITS3Layer.cxx | 68 +++++++++---------- 4 files changed, 59 insertions(+), 53 deletions(-) diff --git a/Detectors/Upgrades/ITS3/base/include/ITS3Base/SpecsV2.h b/Detectors/Upgrades/ITS3/base/include/ITS3Base/SpecsV2.h index a7422c55e72b8..937fa8d2e982c 100644 --- a/Detectors/Upgrades/ITS3/base/include/ITS3Base/SpecsV2.h +++ b/Detectors/Upgrades/ITS3/base/include/ITS3Base/SpecsV2.h @@ -104,7 +104,7 @@ namespace carbonfoam // TODO: Waiting for the further information from WP5(Corrado) constexpr double HringLength{6.0 * mm}; // from blueprint constexpr double longeronsWidth{2.0 * mm}; // what is the height of the longerons? -constexpr double longeronsLength{segment::length - 2 * HringLength}; // 263mm from blueprint; overrriden to be consitent +constexpr double longeronsLength{segment::length - (2 * HringLength)}; // 263mm from blueprint; overrriden to be consitent constexpr double edgeBetwChipAndFoam{1.0 * mm}; // from blueprint but not used cause forms are already overlapping constexpr double gapBetwHringsLongerons{0.05 * mm}; // from blueprint constexpr std::array nHoles{11, 11, 11}; // how many holes for each layer? diff --git a/Detectors/Upgrades/ITS3/reconstruction/src/IOUtils.cxx b/Detectors/Upgrades/ITS3/reconstruction/src/IOUtils.cxx index d7ba4d48dbce4..e898631837169 100644 --- a/Detectors/Upgrades/ITS3/reconstruction/src/IOUtils.cxx +++ b/Detectors/Upgrades/ITS3/reconstruction/src/IOUtils.cxx @@ -10,17 +10,19 @@ // or submit itself to any jurisdiction. #include "ITS3Reconstruction/IOUtils.h" -#include "ITStracking/IOUtils.h" #include "ITStracking/TimeFrame.h" #include "ITStracking/BoundedAllocator.h" #include "DataFormatsITSMFT/CompCluster.h" #include "DataFormatsITSMFT/ROFRecord.h" #include "ITS3Reconstruction/TopologyDictionary.h" #include "ITSBase/GeometryTGeo.h" -#include "ITS3Base/SpecsV2.h" #include "ITStracking/TrackingConfigParam.h" #include "Framework/Logger.h" +#include + +#include + #include namespace o2::its3::ioutils @@ -45,7 +47,7 @@ void convertCompactClusters(gsl::span clusters, } for (auto& c : clusters) { - float sigmaY2, sigmaZ2, sigmaYZ = 0; + float sigmaY2 = NAN, sigmaZ2 = NAN; auto locXYZ = extractClusterData(c, pattIt, dict, sigmaY2, sigmaZ2); const auto detID = c.getSensorID(); auto& cl3d = output.emplace_back(detID, geom->getMatrixT2L(detID) ^ locXYZ); // local --> tracking @@ -54,7 +56,7 @@ void convertCompactClusters(gsl::span clusters, sigmaY2 += conf.sysErrY2[lrID]; sigmaZ2 += conf.sysErrZ2[lrID]; } - cl3d.setErrors(sigmaY2, sigmaZ2, sigmaYZ); + cl3d.setErrors(sigmaY2, sigmaZ2, 0.f); } } @@ -76,26 +78,31 @@ int loadROFrameDataITS3(its::TimeFrame<7>* tf, for (size_t iRof{0}; iRof < rofs.size(); ++iRof) { const auto& rof = rofs[iRof]; for (int clusterId{rof.getFirstEntry()}; clusterId < rof.getFirstEntry() + rof.getNEntries(); ++clusterId) { - auto& c = clusters[clusterId]; - auto sensorID = c.getSensorID(); - auto layer = geom->getLayer(sensorID); + const auto& c = clusters[clusterId]; + const auto sensorID = c.getSensorID(); + const auto layer = geom->getLayer(sensorID); float sigmaY2{0}, sigmaZ2{0}, sigmaYZ{0}; uint8_t clusterSize{0}; - auto locXYZ = extractClusterData(c, pattIt, dict, sigmaY2, sigmaZ2, clusterSize); + const auto locXYZ = extractClusterData(c, pattIt, dict, sigmaY2, sigmaZ2, clusterSize); clusterSizeVec.push_back(clusterSize); // Transformation to the local --> global - auto gloXYZ = geom->getMatrixL2G(sensorID) * locXYZ; + const auto gloXYZ = geom->getMatrixL2G(sensorID) * locXYZ; // Inverse transformation to the local --> tracking - o2::math_utils::Point3D trkXYZ = geom->getMatrixT2L(sensorID) ^ locXYZ; + const o2::math_utils::Point3D trkXYZ = geom->getMatrixT2L(sensorID) ^ locXYZ; // Tracking alpha angle - float alpha = geom->getSensorRefAlpha(sensorID); - - tf->addTrackingFrameInfoToLayer(layer, gloXYZ.x(), gloXYZ.y(), gloXYZ.z(), trkXYZ.x(), alpha, - std::array{trkXYZ.y(), trkXYZ.z()}, + // We want that each cluster rotates its tracking frame to the clusters phi + // that way the track linearization around the measurement is less biases to the arc + // this means automatically that the measurement on the arc is at 0 + // const float alpha = geom->getSensorRefAlpha(sensorID); + const float radius = std::hypot(gloXYZ.x(), gloXYZ.y()); + const float alpha = std::atan2(gloXYZ.y(), gloXYZ.x()); + + tf->addTrackingFrameInfoToLayer(layer, gloXYZ.x(), gloXYZ.y(), gloXYZ.z(), radius, alpha, + std::array{0, trkXYZ.z()}, std::array{sigmaY2, sigmaYZ, sigmaZ2}); /// Rotate to the global frame @@ -103,7 +110,7 @@ int loadROFrameDataITS3(its::TimeFrame<7>* tf, tf->addClusterExternalIndexToLayer(layer, clusterId); } for (unsigned int iL{0}; iL < tf->getUnsortedClusters().size(); ++iL) { - tf->mROFramesClusters[iL][iRof + 1] = tf->getUnsortedClusters()[iL].size(); + tf->mROFramesClusters[iL][iRof + 1] = (int)tf->getUnsortedClusters()[iL].size(); } } diff --git a/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/ITS3Layer.h b/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/ITS3Layer.h index fd9195f9ee228..f45a4469ae2b8 100644 --- a/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/ITS3Layer.h +++ b/Detectors/Upgrades/ITS3/simulation/include/ITS3Simulation/ITS3Layer.h @@ -26,7 +26,7 @@ namespace o2::its3 { /// This class defines the geometry for the ITS3 IB layers. -class ITS3Layer +class ITS3Layer final { // The hierarchy will be the following: // ITS2 -> ITS3 @@ -76,7 +76,6 @@ class ITS3Layer void buildPartial(TGeoVolume* motherVolume, TGeoMatrix* mat = nullptr, BuildLevel level = BuildLevel::kAll, bool createMaterials = false); private: - bool mBuilt{false}; TGeoMedium* mSilicon{nullptr}; TGeoMedium* mAir{nullptr}; TGeoMedium* mCarbon{nullptr}; @@ -91,7 +90,7 @@ class ITS3Layer void createSegment(); void createChip(); void createCarbonForm(); - TGeoCompositeShape* getHringShape(TGeoTubeSeg* Hring); + TGeoCompositeShape* getHringShape(TGeoTubeSeg* Hring) const; void createLayerImpl(); uint8_t mNLayer{0}; // Layer number diff --git a/Detectors/Upgrades/ITS3/simulation/src/ITS3Layer.cxx b/Detectors/Upgrades/ITS3/simulation/src/ITS3Layer.cxx index e0be011096450..c0f8fdc19d03b 100644 --- a/Detectors/Upgrades/ITS3/simulation/src/ITS3Layer.cxx +++ b/Detectors/Upgrades/ITS3/simulation/src/ITS3Layer.cxx @@ -67,11 +67,11 @@ void ITS3Layer::createLayer(TGeoVolume* motherVolume) // Create one layer of ITS3 and attach it to the motherVolume. getMaterials(); createLayerImpl(); - mBuilt = true; if (motherVolume == nullptr) { return; } + // Add it to motherVolume auto* trans = new TGeoTranslation(0, 0, -constants::segment::lengthSensitive / 2.); motherVolume->AddNode(mLayer, 0, trans); @@ -122,8 +122,8 @@ void ITS3Layer::createTile() mTile->AddNode(mPixelArray, 0, phiRotPixelArray); // Biasing - double biasPhi1 = constants::pixelarray::width / mR * o2m::Rad2Deg + readoutPhi2; - double biasPhi2 = biasing::width / mR * o2m::Rad2Deg + biasPhi1; + double biasPhi1 = (constants::pixelarray::width / mR * o2m::Rad2Deg) + readoutPhi2; + double biasPhi2 = (biasing::width / mR * o2m::Rad2Deg) + biasPhi1; auto biasing = new TGeoTubeSeg(mRmin, mRmax, biasing::length / 2, biasPhi1, biasPhi2); auto biasingVol = new TGeoVolume(Form("biasing%d", mNLayer), biasing, mSilicon); biasingVol->SetLineColor(biasing::color); @@ -131,9 +131,9 @@ void ITS3Layer::createTile() mTile->AddNode(biasingVol, 0); // Power Switches are on the side right side of the pixel array and biasing. - auto zMovePowerSwitches = new TGeoTranslation(0, 0, +powerswitches::length / 2. + constants::pixelarray::length / 2.); + auto zMovePowerSwitches = new TGeoTranslation(0, 0, (+powerswitches::length / 2.) + (constants::pixelarray::length / 2.)); double powerPhi1 = readoutPhi2; - double powerPhi2 = powerswitches::width / mR * o2m::Rad2Deg + powerPhi1; + double powerPhi2 = (powerswitches::width / mR * o2m::Rad2Deg) + powerPhi1; auto powerSwitches = new TGeoTubeSeg(mRmin, mRmax, powerswitches::length / 2, powerPhi1, powerPhi2); auto powerSwitchesVol = new TGeoVolume(Form("powerswitches%d", mNLayer), powerSwitches, mSilicon); powerSwitchesVol->SetLineColor(powerswitches::color); @@ -166,7 +166,7 @@ void ITS3Layer::createRSU() // Lower Left auto zMoveLL1 = new TGeoTranslation(0, 0, constants::tile::length); auto zMoveLL2 = new TGeoTranslation(0, 0, constants::tile::length * 2.); - auto zMoveLLDB = new TGeoTranslation(0, 0, -databackbone::length / 2. - constants::pixelarray::length / 2.); + auto zMoveLLDB = new TGeoTranslation(0, 0, (-databackbone::length / 2.) - (constants::pixelarray::length / 2.)); // Lets attach the tiles to the QS. mRSU->AddNode(mTile, nCopyRSU++, nullptr); mRSU->AddNode(mTile, nCopyRSU++, zMoveLL1); @@ -175,9 +175,9 @@ void ITS3Layer::createRSU() // Lower Right auto zMoveLR0 = new TGeoTranslation(0, 0, +length / 2.); - auto zMoveLR1 = new TGeoTranslation(0, 0, constants::tile::length + length / 2.); - auto zMoveLR2 = new TGeoTranslation(0, 0, constants::tile::length * 2. + length / 2.); - auto zMoveLRDB = new TGeoTranslation(0, 0, -databackbone::length / 2. + length / 2. - constants::pixelarray::length / 2.); + auto zMoveLR1 = new TGeoTranslation(0, 0, constants::tile::length + (length / 2.)); + auto zMoveLR2 = new TGeoTranslation(0, 0, (constants::tile::length * 2.) + (length / 2.)); + auto zMoveLRDB = new TGeoTranslation(0, 0, (-databackbone::length / 2.) + (length / 2.) - (constants::pixelarray::length / 2.)); // Lets attach the tiles to the QS. mRSU->AddNode(mTile, nCopyRSU++, zMoveLR0); mRSU->AddNode(mTile, nCopyRSU++, zMoveLR1); @@ -192,7 +192,7 @@ void ITS3Layer::createRSU() // Upper Left auto zMoveUL1 = new TGeoCombiTrans(0, 0, constants::tile::length, rot); auto zMoveUL2 = new TGeoCombiTrans(0, 0, constants::tile::length * 2., rot); - auto zMoveULDB = new TGeoCombiTrans(0, 0, -databackbone::length / 2. - constants::pixelarray::length / 2., rot); + auto zMoveULDB = new TGeoCombiTrans(0, 0, (-databackbone::length / 2.) - (constants::pixelarray::length / 2.), rot); // Lets attach the tiles to the QS. mRSU->AddNode(mTile, nCopyRSU++, rot); mRSU->AddNode(mTile, nCopyRSU++, zMoveUL1); @@ -201,9 +201,9 @@ void ITS3Layer::createRSU() // Upper Right auto zMoveUR0 = new TGeoCombiTrans(0, 0, +length / 2., rot); - auto zMoveUR1 = new TGeoCombiTrans(0, 0, constants::tile::length + length / 2., rot); - auto zMoveUR2 = new TGeoCombiTrans(0, 0, constants::tile::length * 2. + length / 2., rot); - auto zMoveURDB = new TGeoCombiTrans(0, 0, -databackbone::length / 2. + length / 2. - constants::pixelarray::length / 2., rot); + auto zMoveUR1 = new TGeoCombiTrans(0, 0, constants::tile::length + (length / 2.), rot); + auto zMoveUR2 = new TGeoCombiTrans(0, 0, (constants::tile::length * 2.) + (length / 2.), rot); + auto zMoveURDB = new TGeoCombiTrans(0, 0, (-databackbone::length / 2.) + (length / 2.) - (constants::pixelarray::length / 2.), rot); // Lets attach the tiles to the QS. mRSU->AddNode(mTile, nCopyRSU++, zMoveUR0); mRSU->AddNode(mTile, nCopyRSU++, zMoveUR1); @@ -225,9 +225,9 @@ void ITS3Layer::createSegment() mSegment = new TGeoVolumeAssembly(its3TGeo::getITS3SegmentPattern(mNLayer)); mSegment->VisibleDaughters(); - for (size_t i{0}; i < nRSUs; ++i) { - auto zMove = new TGeoTranslation(0, 0, +i * constants::rsu::length + constants::rsu::databackbone::length + constants::pixelarray::length / 2.); - mSegment->AddNode(mRSU, i, zMove); + for (unsigned int i{0}; i < nRSUs; ++i) { + auto zMove = new TGeoTranslation(0, 0, (i * constants::rsu::length) + constants::rsu::databackbone::length + (constants::pixelarray::length / 2.)); + mSegment->AddNode(mRSU, (int)i, zMove); } // LEC @@ -242,7 +242,7 @@ void ITS3Layer::createSegment() mSegment->AddNode(lecVol, 0, zMoveLEC); // REC; reuses lecPhi1,2 - auto zMoveREC = new TGeoTranslation(0, 0, nRSUs * constants::rsu::length + rec::length / 2.); + auto zMoveREC = new TGeoTranslation(0, 0, (nRSUs * constants::rsu::length) + (rec::length / 2.)); auto rec = new TGeoTubeSeg(mRmin, mRmax, rec::length / 2., lecPhi1, lecPhi2); auto recVol = new TGeoVolume(Form("rec%d", mNLayer), rec, mSilicon); @@ -266,11 +266,11 @@ void ITS3Layer::createChip() auto phiOffset = constants::segment::width / mR * o2m::Rad2Deg; for (unsigned int i{0}; i < constants::nSegments[mNLayer]; ++i) { auto rot = new TGeoRotation(Form("its3PhiSegmentOffset_%d_%d", mNLayer, i), 0, 0, phiOffset * i); - mChip->AddNode(mSegment, i, rot); + mChip->AddNode(mSegment, (int)i, rot); } // Add metal stack positioned radially outward - auto zMoveMetal = new TGeoTranslation(0, 0, constants::metalstack::length / 2. - constants::segment::lec::length); + auto zMoveMetal = new TGeoTranslation(0, 0, (constants::metalstack::length / 2.) - constants::segment::lec::length); auto metal = new TGeoTubeSeg(mRmax, mRmax + constants::metalstack::thickness, constants::metalstack::length / 2., 0, constants::nSegments[mNLayer] * phiOffset); auto metalVol = new TGeoVolume(Form("metal%d", mNLayer), metal, mCopper); metalVol->SetLineColor(constants::metalstack::color); @@ -296,7 +296,7 @@ void ITS3Layer::createCarbonForm() dRadius = constants::carbonfoam::thicknessOuterFoam; // TODO: lack of carbon foam radius for layer 2, use 0.7 cm as a temporary value } double phiSta = edgeBetwChipAndFoam / (0.5 * constants::radii[mNLayer + 1] + constants::radii[mNLayer]) * o2m::Rad2Deg; - double phiEnd = (constants::nSegments[mNLayer] * constants::segment::width) / constants::radii[mNLayer] * o2m::Rad2Deg - phiSta; + double phiEnd = ((constants::nSegments[mNLayer] * constants::segment::width) / constants::radii[mNLayer] * o2m::Rad2Deg) - phiSta; double phiLongeronsCover = longeronsWidth / (0.5 * constants::radii[mNLayer + 1] + constants::radii[mNLayer]) * o2m::Rad2Deg; // H-rings foam @@ -308,35 +308,37 @@ void ITS3Layer::createCarbonForm() HringCVol->SetLineColor(color); auto HringAVol = new TGeoVolume(Form("hringA%d", mNLayer), HringAWithHoles, mCarbon); HringAVol->SetLineColor(color); - auto zMoveHringC = new TGeoTranslation(0, 0, -constants::segment::lec::length + HringLength / 2.); - auto zMoveHringA = new TGeoTranslation(0, 0, -constants::segment::lec::length + HringLength / 2. + constants::segment::length - HringLength); + auto zMoveHringC = new TGeoTranslation(0, 0, -constants::segment::lec::length + (HringLength / 2.)); + auto zMoveHringA = new TGeoTranslation(0, 0, -constants::segment::lec::length + (HringLength / 2.) + constants::segment::length - HringLength); // Longerons are made by same material + // added separately to make navigation faster [[maybe_unused]] auto longeronR = new TGeoTubeSeg(Form("longeronR%d", mNLayer), mRmax, mRmax + dRadius, longeronsLength / 2., phiSta, phiSta + phiLongeronsCover); [[maybe_unused]] auto longeronL = new TGeoTubeSeg(Form("longeronL%d", mNLayer), mRmax, mRmax + dRadius, longeronsLength / 2., phiEnd - phiLongeronsCover, phiEnd); - TString nameLongerons = Form("longeronR%d + longeronL%d", mNLayer, mNLayer); - auto longerons = new TGeoCompositeShape(nameLongerons); - auto longeronsVol = new TGeoVolume(Form("longerons%d", mNLayer), longerons, mCarbon); - longeronsVol->SetLineColor(color); - auto zMoveLongerons = new TGeoTranslation(0, 0, -constants::segment::lec::length + constants::segment::length / 2.); + auto longeronRVol = new TGeoVolume(Form("longeronR%d", mNLayer), longeronR, mCarbon); + longeronRVol->SetLineColor(color); + auto longeronLVol = new TGeoVolume(Form("longeronL%d", mNLayer), longeronL, mCarbon); + longeronLVol->SetLineColor(color); + auto zMoveLongerons = new TGeoTranslation(0, 0, -constants::segment::lec::length + (constants::segment::length / 2.)); mCarbonForm->AddNode(HringCVol, 0, zMoveHringC); mCarbonForm->AddNode(HringAVol, 0, zMoveHringA); - mCarbonForm->AddNode(longeronsVol, 0, zMoveLongerons); + mCarbonForm->AddNode(longeronRVol, 0, zMoveLongerons); + mCarbonForm->AddNode(longeronLVol, 0, zMoveLongerons); mCarbonForm->AddNode(mChip, 0); } -TGeoCompositeShape* ITS3Layer::getHringShape(TGeoTubeSeg* Hring) +TGeoCompositeShape* ITS3Layer::getHringShape(TGeoTubeSeg* Hring) const { // Function to dig holes in H-rings using namespace constants::carbonfoam; double stepPhiHoles = (Hring->GetPhi2() - Hring->GetPhi1()) / (nHoles[mNLayer]); - double phiHolesSta = Hring->GetPhi1() + stepPhiHoles / 2.; + double phiHolesSta = Hring->GetPhi1() + (stepPhiHoles / 2.); double radiusHring = 0.5 * (Hring->GetRmin() + Hring->GetRmax()); TGeoCompositeShape* HringWithHoles = nullptr; TString nameAllHoles = ""; for (int iHoles = 0; iHoles < nHoles[mNLayer]; iHoles++) { - double phiHole = phiHolesSta + stepPhiHoles * iHoles; + double phiHole = phiHolesSta + (stepPhiHoles * iHoles); TString nameHole = Form("hole_%d_%d", iHoles, mNLayer); [[maybe_unused]] auto hole = new TGeoTube(nameHole, 0, radiusHoles[mNLayer], 3 * Hring->GetDz()); // move hole to the hring radius @@ -376,9 +378,7 @@ void ITS3Layer::createLayerImpl() void ITS3Layer::buildPartial(TGeoVolume* motherVolume, TGeoMatrix* mat, BuildLevel level, bool createMaterials) { - if (!mBuilt) { - getMaterials(createMaterials); - } + getMaterials(createMaterials); switch (level) { case BuildLevel::kPixelArray: createPixelArray(); From 910fc5e1e009ae49b4af358c6d909c77201c68ac Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Thu, 5 Feb 2026 17:45:16 +0100 Subject: [PATCH 4/5] ITS3: alignment Signed-off-by: Felix Schlepper --- .../Upgrades/ITS3/alignment/CMakeLists.txt | 30 +- .../include/ITS3Align/AlignmentHierarchy.h | 258 ++++++ .../include/ITS3Align/AlignmentParams.h | 54 ++ .../include/ITS3Align/AlignmentSensors.h | 41 + .../include/ITS3Align/AlignmentSpec.h | 23 + .../include/ITS3Align/AlignmentTypes.h | 67 ++ .../ITS3/alignment/src/AlignmentHierarchy.cxx | 201 +++++ .../ITS3/alignment/src/AlignmentParams.cxx | 13 + .../ITS3/alignment/src/AlignmentSensors.cxx | 203 +++++ .../ITS3/alignment/src/AlignmentSpec.cxx | 810 ++++++++++++++++++ .../ITS3/alignment/src/AlignmentTypes.cxx | 24 + .../ITS3/alignment/src/ITS3AlignLinkDef.h | 8 + .../ITS3/alignment/src/alignment-workflow.cxx | 64 ++ .../ITS3/reconstruction/src/IOUtils.cxx | 27 +- dependencies/O2Dependencies.cmake | 18 + 15 files changed, 1823 insertions(+), 18 deletions(-) create mode 100644 Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentHierarchy.h create mode 100644 Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentParams.h create mode 100644 Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentSensors.h create mode 100644 Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentSpec.h create mode 100644 Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentTypes.h create mode 100644 Detectors/Upgrades/ITS3/alignment/src/AlignmentHierarchy.cxx create mode 100644 Detectors/Upgrades/ITS3/alignment/src/AlignmentParams.cxx create mode 100644 Detectors/Upgrades/ITS3/alignment/src/AlignmentSensors.cxx create mode 100644 Detectors/Upgrades/ITS3/alignment/src/AlignmentSpec.cxx create mode 100644 Detectors/Upgrades/ITS3/alignment/src/AlignmentTypes.cxx create mode 100644 Detectors/Upgrades/ITS3/alignment/src/alignment-workflow.cxx diff --git a/Detectors/Upgrades/ITS3/alignment/CMakeLists.txt b/Detectors/Upgrades/ITS3/alignment/CMakeLists.txt index f89ad821c65e7..5a3f3dd0d2d07 100644 --- a/Detectors/Upgrades/ITS3/alignment/CMakeLists.txt +++ b/Detectors/Upgrades/ITS3/alignment/CMakeLists.txt @@ -9,18 +9,42 @@ # granted to it by virtue of its status as an Intergovernmental Organization # or submit itself to any jurisdiction. +# add_compile_options(-O0 -g -fPIC -fno-omit-frame-pointer) o2_add_library(ITS3Align + TARGETVARNAME targetName SOURCES src/MisalignmentParameters.cxx src/MisalignmentHits.cxx src/MisalignmentManager.cxx src/Deformations.cxx + src/AlignmentHierarchy.cxx + src/AlignmentSensors.cxx + src/AlignmentParams.cxx + src/AlignmentTypes.cxx + src/AlignmentSpec.cxx PUBLIC_LINK_LIBRARIES O2::MathUtils O2::Steer O2::ITSBase - O2::ITSMFTSimulation) + O2::ITSMFTSimulation + O2::Framework + O2::GlobalTrackingWorkflowReaders + O2::GlobalTrackingWorkflowHelpers + O2::DataFormatsGlobalTracking + O2::DetectorsVertexing + GBL::GBL) +if (OpenMP_CXX_FOUND) + target_compile_definitions(${targetName} PRIVATE WITH_OPENMP) + target_link_libraries(${targetName} PRIVATE OpenMP::OpenMP_CXX) +endif() o2_target_root_dictionary(ITS3Align HEADERS include/ITS3Align/MisalignmentParameters.h include/ITS3Align/MisalignmentHits.h - include/ITS3Align/MisalignmentHits.h - include/ITS3Align/Deformations.h) + include/ITS3Align/Deformations.h + include/ITS3Align/AlignmentParams.h + include/ITS3Align/AlignmentTypes.h) + + +o2_add_executable(alignment-workflow + SOURCES src/alignment-workflow.cxx + COMPONENT_NAME its3 + PUBLIC_LINK_LIBRARIES O2::ITS3Align) diff --git a/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentHierarchy.h b/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentHierarchy.h new file mode 100644 index 0000000000000..7dea295f0b264 --- /dev/null +++ b/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentHierarchy.h @@ -0,0 +1,258 @@ +// Copyright 2019-2026 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. + +#ifndef O2_ITS3_ALIGNMENT_HIERARCHY_H +#define O2_ITS3_ALIGNMENT_HIERARCHY_H + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +namespace o2::its3::align +{ +using Matrix36 = Eigen::Matrix; +using Matrix66 = Eigen::Matrix; + +using RigidBodyDOFMask = uint8_t; +template +constexpr RigidBodyDOFMask DOF_BIT(Args... bits) +{ + return (RigidBodyDOFMask{0} | ... | (RigidBodyDOFMask{1} << bits)); +} +// DOFs are defined in LOC +enum RigidBodyDOF : RigidBodyDOFMask { + TX = 0, // Translation along local X + TY, // Translation along local Y + TZ, // Translation along local Z + RX, // Rotation around local X + RY, // Rotation around local Y + RZ, // Rotation around local Z + NDOF, +}; +constexpr RigidBodyDOFMask RigidBodyDOFTra = DOF_BIT(TX, TY, TZ); +constexpr RigidBodyDOFMask RigidBodyDOFRot = DOF_BIT(RX, RY, RZ); +constexpr RigidBodyDOFMask RigidBodyDOFAll = RigidBodyDOFTra | RigidBodyDOFRot; +constexpr RigidBodyDOFMask RigidBodyDOFNone = 0; +constexpr RigidBodyDOFMask RigidBodyDOFPseudo = std::numeric_limits::max(); // special value representing that the volume itself does not have any dofs but should not curtail the parent's ones +static constexpr const char* RigidBodyDOFNames[RigidBodyDOF::NDOF] = {"TX", "TY", "TZ", "RX", "RY", "RZ"}; +inline bool hasRigidBodyDOF(RigidBodyDOFMask m, RigidBodyDOFMask d) +{ + return (m != RigidBodyDOFPseudo) && (m & DOF_BIT(d)); +} +inline void enableRigidBodyDOF(RigidBodyDOFMask& m, RigidBodyDOFMask d) +{ + m |= DOF_BIT(d); +} +inline void disableRigidBodyDOF(RigidBodyDOFMask& m, RigidBodyDOFMask d) +{ + m &= ~DOF_BIT(d); +} +// return the rigid body derivatives +// trk has be at in the measurment frame +auto getRigidBodyDerivatives(const auto& trk) +{ + // calculate slopes + const double tgl = trk.getTgl(), snp = trk.getSnp(); + const double csp = 1. / sqrt(1. + (tgl * tgl)); + const double u = trk.getY(), v = trk.getZ(); + const double uP = snp * csp, vP = tgl * csp; + Matrix36 der; + der.setZero(); + // columns: Tt, Tu, Tv, Rt, Ru, Rv + // (X) (Y) (Z) (RX) (RY) (RZ) + der << uP, -1., 0., v, v * uP, -u * uP, + vP, 0., -1., -u, v * vP, -u * vP; + return der; +} + +class GlobalLabel +{ + // Millepede label is any positive integer [1....) + public: + using T = uint32_t; + static constexpr int DOF_BITS = 5; // 0...4 + static constexpr int ID_BITS = 23; // 5...27 + static constexpr int SENS_BITS = 1; // 28 + static constexpr int TOTAL_BITS = sizeof(T) * 8; + static constexpr int DET_BITS = TOTAL_BITS - (DOF_BITS + ID_BITS + SENS_BITS) - 1; // one less bit since GBL uses int! + static constexpr T bitMask(int b) noexcept + { + return (T(1) << b) - T(1); + } + static constexpr int DOF_SHIFT = 0; + static constexpr T DOF_MAX = (T(1) << DOF_BITS) - T(1); + static constexpr T DOF_MASK = DOF_MAX << DOF_SHIFT; + static constexpr int ID_SHIFT = DOF_BITS; + static constexpr T ID_MAX = (T(1) << ID_BITS) - T(1); + static constexpr T ID_MASK = ID_MAX << ID_SHIFT; + static constexpr int SENS_SHIFT = ID_BITS + DOF_BITS; + static constexpr T SENS_MAX = (T(1) << SENS_BITS) - T(1); + static constexpr T SENS_MASK = SENS_MAX << SENS_SHIFT; + static constexpr int DET_SHIFT = DOF_BITS + ID_BITS + SENS_BITS; + static constexpr T DET_MAX = (T(1) << DET_BITS) - T(1); + static constexpr T DET_MASK = DET_MAX << DET_SHIFT; + + GlobalLabel(T det, T id, bool sens) : mID(((((id + 1) & ID_MAX) << ID_SHIFT) | ((det & DET_MAX) << DET_SHIFT) | ((sens & SENS_MAX) << SENS_SHIFT))) + { + } + + constexpr T raw(T dof) const noexcept { return (mID & ~DOF_MASK) | ((dof & DOF_MAX) << DOF_SHIFT); } + constexpr int rawGBL(T dof) const noexcept { return static_cast(raw(dof)); } + constexpr T id() const noexcept + { + return ((mID - 1) & ID_MASK) >> ID_SHIFT; + } + constexpr T det() const noexcept + { + return (mID & DET_MASK) >> DET_SHIFT; + } + constexpr bool sens() const noexcept + { + return (mID & SENS_MASK) >> SENS_SHIFT; + } + + std::string asString() const + { + return std::format("Det:{} Id:{} Sens:{}", det(), id(), sens()); + } + + constexpr auto operator<=>(const GlobalLabel&) const noexcept = default; + + private: + T mID{0}; +}; + +// Rigid body constraints for the parents +class HierarchyConstraint +{ + public: + HierarchyConstraint(std::string name, double value) : mName(std::move(name)), mValue(value) {} + void add(uint32_t lab, double coeff) + { + mLabels.push_back(lab); + mCoeff.push_back(coeff); + } + void write(std::ostream& os) const; + auto getSize() const noexcept { return mLabels.size(); } + + private: + std::string mName; // name of the constraint + double mValue{0.0}; // constraint value + std::vector mLabels; // parameter labels + std::vector mCoeff; // their coefficients +}; + +class AlignableVolume +{ + public: + using Ptr = std::unique_ptr; + using SensorMapping = std::map; + + AlignableVolume(const AlignableVolume&) = default; + AlignableVolume(AlignableVolume&&) = delete; + AlignableVolume& operator=(const AlignableVolume&) = default; + AlignableVolume& operator=(AlignableVolume&&) = delete; + AlignableVolume(const char* symName, uint32_t label, uint32_t det, bool sens, RigidBodyDOFMask dof); + AlignableVolume(const char* symName, GlobalLabel label, RigidBodyDOFMask dof); + virtual ~AlignableVolume() = default; + + void finalise(uint8_t level = 0); + + // steering file output + void writeRigidBodyConstraints(std::ostream& os) const; + void writeParameters(std::ostream& os) const; + void writeTree(std::ostream& os, int indent = 0) const; + + // tree-like + auto getLevel() const noexcept { return mLevel; } + bool isRoot() const noexcept { return mParent == nullptr; } + bool isLeaf() const noexcept { return mChildren.empty(); } + template + requires std::derived_from + AlignableVolume* addChild(const char* symName, uint32_t label, uint32_t det, bool sens, RigidBodyDOFMask dof) + { + auto c = std::make_unique(symName, label, det, sens, dof); + return setParent(std::move(c)); + } + template + requires std::derived_from + AlignableVolume* addChild(const char* symName, GlobalLabel lbl, RigidBodyDOFMask dof) + { + auto c = std::make_unique(symName, lbl, dof); + return setParent(std::move(c)); + } + + // bfs traversal + void traverse(const std::function& visitor) + { + visitor(this); + for (auto& c : mChildren) { + c->traverse(visitor); + } + } + + std::string getSymName() const noexcept { return mSymName; } + GlobalLabel getLabel() const noexcept { return mLabel; } + AlignableVolume* getParent() const { return mParent; } + size_t getNChildren() const noexcept { return mChildren.size(); } + void setRigidBodyDOF(RigidBodyDOFMask m) noexcept { mDOF = m; } + RigidBodyDOFMask getRigidBodyDOF() const noexcept { return mDOF; } + + // transformation matrices + virtual void defineMatrixL2G() {} + virtual void defineMatrixT2L() {} + virtual void computeJacobianL2T(const double* pos, Matrix66& jac) const {}; + const TGeoHMatrix& getL2P() const { return mL2P; } + const TGeoHMatrix& getT2L() const { return mT2L; } + const Matrix66& getJL2P() const { return mJL2P; } + const Matrix66& getJP2L() const { return mJP2L; } + + protected: + /// matrices + AlignableVolume* mParent{nullptr}; // parent + TGeoPNEntry* mPNE{nullptr}; // physical entry + TGeoPhysicalNode* mPN{nullptr}; // physical node + TGeoHMatrix mL2G; // (LOC) -> (GLO) + TGeoHMatrix mL2P; // (LOC) -> (PAR) + Matrix66 mJL2P; // jac (LOC) -> (PAR) + Matrix66 mJP2L; // jac (PAR) -> (LOC) + TGeoHMatrix mT2L; // (TRK) -> (LOC) + + private: + double mPreSigma{0.0}; // asigned pre-sigma + std::string mSymName; // unique symbolic name + GlobalLabel mLabel; // internal global idetenifier + uint8_t mLevel{0}; // depth-in tree + RigidBodyDOFMask mDOF{RigidBodyDOFAll}; // allowed dofs + + AlignableVolume* setParent(Ptr c) + { + c->mParent = this; + mChildren.push_back(std::move(c)); + return mChildren.back().get(); + } + std::vector mChildren; // children + + void init(); +}; + +} // namespace o2::its3::align + +#endif diff --git a/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentParams.h b/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentParams.h new file mode 100644 index 0000000000000..f267a5fe0da7f --- /dev/null +++ b/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentParams.h @@ -0,0 +1,54 @@ +// 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. +#ifndef ALICEO2_ITS3_ALIGNMENTPARAMS_H_ +#define ALICEO2_ITS3_ALIGNMENTPARAMS_H_ + +#include "CommonUtils/ConfigurableParam.h" +#include "CommonUtils/ConfigurableParamHelper.h" +#include "DetectorsBase/Propagator.h" + +namespace o2::its3::align +{ +struct AlignmentParams : public o2::conf::ConfigurableParamHelper { + // Track selection + float minPt = 1.f; // minimum pt required + int minITSCls = 7; // minimum number of ITS clusters + float maxITSChi2Ndf = 1.2; // maximum ITS track chi2 + + // propagation opt + double maxSnp = 0.85; + double maxStep = 2.0; + // o2::base::PropagatorD::MatCorrType matCorrType = o2::base::PropagatorD::MatCorrType::USEMatCorrTGeo; + o2::base::PropagatorD::MatCorrType matCorrType = o2::base::PropagatorD::MatCorrType::USEMatCorrLUT; + + bool useStableRefit = true; // use input tracks as linearization point + float minMS = 1e-6f; // minimum scattering to account for + float maxChi2Ndf = 10; // maximum Chi2/Ndf allowed for GBL fit + + // Ridder options + int ridderMaxExtrap = 10; + double ridderRelIniStep[5] = {0.01, 0.01, 0.02, 0.02, 0.02}; + double ridderMaxIniStep[5] = {0.1, 0.1, 0.05, 0.05, 0.05}; + double ridderShrinkFac = 2.0; + double ridderEps = 1e-16; + double ridderMaxJacDiagTol = 0.1; // max tolerance of diagonal elements away from 1 + + // MillePede output + std::string milleBinFile = "mp2data.bin"; + std::string milleConFile = "mp2con.txt"; + std::string milleParamFile = "mp2param.txt"; + std::string milleTreeFile = "mp2tree.txt"; + + O2ParamDef(AlignmentParams, "ITS3AlignmentParams"); +}; +} // namespace o2::its3::align + +#endif diff --git a/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentSensors.h b/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentSensors.h new file mode 100644 index 0000000000000..535f67156a16c --- /dev/null +++ b/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentSensors.h @@ -0,0 +1,41 @@ +// Copyright 2019-2026 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. + +#ifndef O2_ITS3_ALIGNMENT_SENSORS_H +#define O2_ITS3_ALIGNMENT_SENSORS_H + +#include "ITS3Align/AlignmentHierarchy.h" + +namespace o2::its3::align +{ + +AlignableVolume::Ptr buildHierarchyITS(AlignableVolume::SensorMapping& sensorMap); +AlignableVolume::Ptr buildHierarchyIT3(AlignableVolume::SensorMapping& sensorMap); + +class AlignableSensorITS final : public AlignableVolume +{ + using AlignableVolume::AlignableVolume; + void defineMatrixL2G() final; + void defineMatrixT2L() final; + void computeJacobianL2T(const double* pos, Matrix66& jac) const final; +}; + +class AlignableSensorIT3 final : public AlignableVolume +{ + using AlignableVolume::AlignableVolume; + void defineMatrixL2G() final; + void defineMatrixT2L() final; + void computeJacobianL2T(const double* pos, Matrix66& jac) const final; +}; + +} // namespace o2::its3::align + +#endif diff --git a/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentSpec.h b/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentSpec.h new file mode 100644 index 0000000000000..84422f774b202 --- /dev/null +++ b/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentSpec.h @@ -0,0 +1,23 @@ +// 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. + +#ifndef O2_ITS3_ALIGNMENT_H +#define O2_ITS3_ALIGNMENT_H + +#include "ReconstructionDataFormats/GlobalTrackID.h" +#include "Framework/DataProcessorSpec.h" + +namespace o2::its3::align +{ +o2::framework::DataProcessorSpec getAlignmentSpec(o2::dataformats::GlobalTrackID::mask_t srcTracks, o2::dataformats::GlobalTrackID::mask_t srcClus, bool useMC, bool withPV, bool withITS3); +} // namespace o2::its3::align + +#endif diff --git a/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentTypes.h b/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentTypes.h new file mode 100644 index 0000000000000..605b50f350a24 --- /dev/null +++ b/Detectors/Upgrades/ITS3/alignment/include/ITS3Align/AlignmentTypes.h @@ -0,0 +1,67 @@ +// 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. + +#ifndef O2_ITS3_ALIGNMENT_TYPES_H +#define O2_ITS3_ALIGNMENT_TYPES_H + +#include +#include "ReconstructionDataFormats/Track.h" +#include "DataFormatsITS/TrackITS.h" + +namespace o2::its3::align +{ + +struct Measurement final { + double dy = 0.f; + double dz = 0.f; + double sig2y = 0.f; + double sig2z = 0.f; + double phi = 0.f; + double z = 0.f; + ClassDefNV(Measurement, 1) +}; + +struct FrameInfoExt final { + int16_t sens = -1; + int8_t lr = -1; // -1 = vtx + o2::math_utils::Point3D trk; + o2::math_utils::Point3D loc; + o2::math_utils::Point3D glo; + float x{-999.f}; + float alpha{-999.f}; + std::array positionTrackingFrame = {999., 999.}; + std::array covarianceTrackingFrame = {999., 999., 999.}; + + std::string asString() const; + + ClassDefNV(FrameInfoExt, 1) +}; + +struct FitInfo final { + float chi2Ndf{-1}; // Chi2/Ndf of track refit + float chi2{-1}; // Chi2 + int ndf; // ndf + ClassDefNV(FitInfo, 1) +}; + +struct Track { + o2::its::TrackITS its; // original ITS track + o2::track::TrackParCovD track; // track at innermost update point, refitted from outwards seed + FitInfo kfFit; // kf fit information + FitInfo gblFit; // gbl fit information + std::vector points; // measurment point + std::vector info; // frame info + ClassDefNV(Track, 1) +}; + +} // namespace o2::its3::align + +#endif diff --git a/Detectors/Upgrades/ITS3/alignment/src/AlignmentHierarchy.cxx b/Detectors/Upgrades/ITS3/alignment/src/AlignmentHierarchy.cxx new file mode 100644 index 0000000000000..c41ac5c90d1ab --- /dev/null +++ b/Detectors/Upgrades/ITS3/alignment/src/AlignmentHierarchy.cxx @@ -0,0 +1,201 @@ +// Copyright 2019-2026 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. + +#include + +#include +#include + +#include "ITS3Align/AlignmentHierarchy.h" +#include "ITSBase/GeometryTGeo.h" +#include "Framework/Logger.h" +#include "MathUtils/Utils.h" + +using namespace o2::its3::align; + +void HierarchyConstraint::write(std::ostream& os) const +{ + os << "!!! " << mName << '\n'; + os << "Constraint " << mValue << '\n'; + for (size_t i{0}; i < mLabels.size(); ++i) { + os << mLabels[i] << " " << mCoeff[i] << '\n'; + } + os << '\n'; +} + +AlignableVolume::AlignableVolume(const char* symName, uint32_t label, uint32_t det, bool sens, RigidBodyDOFMask dof) : mSymName(symName), mLabel(det, label, sens), mDOF(dof) +{ + init(); +} + +AlignableVolume::AlignableVolume(const char* symName, GlobalLabel label, RigidBodyDOFMask dof) : mSymName(symName), mLabel(label), mDOF(dof) +{ + init(); +} + +void AlignableVolume::init() +{ + // check if this sym volume actually exists + mPNE = gGeoManager->GetAlignableEntry(mSymName.c_str()); + if (mPNE == nullptr) { + LOGP(fatal, "Symbolic volume '{}' has no corresponding alignable entry!", mSymName); + } + mPN = mPNE->GetPhysicalNode(); + if (mPN == nullptr) { + LOGP(debug, "Adding physical node to {}", mSymName); + mPN = gGeoManager->MakePhysicalNode(mPNE->GetTitle()); + if (mPN == nullptr) { + LOGP(fatal, "Failed to make physical node for {}", mSymName); + } + } +} + +void AlignableVolume::finalise(uint8_t level) +{ + if (level == 0 && !isRoot()) { + LOGP(fatal, "Finalise should be called only from the root node!"); + } + mLevel = level; + if (!isLeaf()) { + // depth first + for (const auto& c : mChildren) { + c->finalise(level + 1); + } + // check if any children have any dofs otherwise disable the dof for this parent + for (uint8_t iDOF{0}; iDOF < RigidBodyDOF::NDOF; ++iDOF) { + if (!(mDOF & DOF_BIT(iDOF))) { + continue; + } + int nChildrenDOF{0}; // count of children with any dof active + for (const auto& c : mChildren) { + if (c->getRigidBodyDOF()) { + ++nChildrenDOF; + } + } + if (!nChildrenDOF) { + LOGP(warn, "Auto-disabling DOF {} for {} since no active children", RigidBodyDOFNames[iDOF], mSymName); + disableRigidBodyDOF(mDOF, iDOF); + } + } + } else { + // for sensors we need also to define the transformation from the measurment (TRK) to the local frame (LOC) + // need to it with including possible pre-alignment to allow for iterative convergence + // (TRK) is defined wrt global z-axis + defineMatrixL2G(); + defineMatrixT2L(); + } + if (!isRoot()) { + // prepare the transformation matrices, e.g. from child frame to parent frame + // this is not necessarily just one level transformation + TGeoHMatrix mat = *mPN->GetMatrix(); // global matrix (including possible pre-alignment) from this volume to the global frame + if (isLeaf()) { + mat = mL2G; // for sensor volumes they might have redefined the L2G definition + } + auto inv = mParent->mPN->GetMatrix()->Inverse(); // global (including possible pre-alignment) from this volume to the global frame + mat.MultiplyLeft(inv); // left mult. effectively subtracts the parent transformation which is included in the the childs + mL2P = mat; // now this is directly the child to the parent transformation (LOC) (including possible pre-alignment) + + // prepare jacobian from child to parent frame + Eigen::Map> rotL2P(mL2P.GetRotationMatrix()); + Eigen::Matrix3d rotInv = rotL2P.transpose(); // parent-to-child rotation + const double* t = mL2P.GetTranslation(); // child origin in parent frame + Eigen::Matrix3d skewT; + skewT << 0, -t[2], t[1], t[2], 0, -t[0], -t[1], t[0], 0; + mJL2P.setZero(); + mJL2P.topLeftCorner<3, 3>() = rotInv; + mJL2P.topRightCorner<3, 3>() = -rotInv * skewT; + mJL2P.bottomRightCorner<3, 3>() = rotInv; + mJP2L = mJL2P.inverse(); + } +} + +void AlignableVolume::writeRigidBodyConstraints(std::ostream& os) const +{ + // generate hierarchy constraints: parent = averge of children + // (1/N) * sum_k sum_j (J_hk^{-1})_{i,j} * child_k_j = 0 + + if (isLeaf()) { + return; + } + + for (uint8_t iDOF{0}; iDOF < RigidBodyDOF::NDOF; ++iDOF) { + if (!(mDOF & DOF_BIT(iDOF))) { + continue; + } + double nChildrenDOF{0.}; // count of children with this dof active + for (const auto& c : mChildren) { + if (c->getRigidBodyDOF()) { + ++nChildrenDOF; + } + } + if (nChildrenDOF == 0.) { + LOGP(fatal, "{} has dof {} active but no children with this active dof!", mSymName, RigidBodyDOFNames[iDOF]); + } + const double invN = 1.0 / nChildrenDOF; + HierarchyConstraint con(std::format("DOF {} for {}", RigidBodyDOFNames[iDOF], mSymName), 0.0); + for (const auto& c : mChildren) { + if (!c->getRigidBodyDOF()) { + continue; + } + for (uint8_t jDOF{0}; jDOF < RigidBodyDOF::NDOF; ++jDOF) { + if (!hasRigidBodyDOF(c->getRigidBodyDOF(), jDOF)) { + continue; + } + double coeff = invN * c->getJP2L()(iDOF, jDOF); + if (std::abs(coeff) > 1e-16f) { + con.add(c->getLabel().raw(jDOF), coeff); + } + } + } + + if (con.getSize() > 1) { + con.write(os); + } + } + for (const auto& c : mChildren) { + c->writeRigidBodyConstraints(os); + } +} + +void AlignableVolume::writeParameters(std::ostream& os) const +{ + if (isRoot()) { + os << "Parameter\n"; + } + for (uint8_t iDOF{0}; iDOF < RigidBodyDOF::NDOF; ++iDOF) { + os << std::format("{:<10} {:>+15g} {:>+15g} ! {} {} ", mLabel.raw(iDOF), mPreSigma, (hasRigidBodyDOF(mDOF, iDOF) ? 0.0 : -1.0), (hasRigidBodyDOF(mDOF, iDOF) ? 'V' : 'F'), RigidBodyDOFNames[iDOF]) << mSymName << '\n'; + } + for (const auto& c : mChildren) { + c->writeParameters(os); + } +} + +void AlignableVolume::writeTree(std::ostream& os, int indent) const +{ + os << std::string(static_cast(indent * 2), ' ') << mSymName << (mLabel.sens() ? " (sens)" : " (pasv)"); + if (mDOF == RigidBodyDOFPseudo) { + os << " is a pseudo-volume"; + } else if (mDOF == RigidBodyDOFNone) { + os << " no DOFs " << mLabel.raw(0); + } else { + os << " with DOFs ["; + for (uint8_t iDOF{0}; iDOF < RigidBodyDOF::NDOF; ++iDOF) { + if (hasRigidBodyDOF(mDOF, iDOF)) { + os << " " << RigidBodyDOFNames[iDOF] << " (" << mLabel.raw(iDOF) << ") "; + } + } + os << " ]"; + } + os << '\n'; + for (const auto& c : mChildren) { + c->writeTree(os, indent + 2); + } +} diff --git a/Detectors/Upgrades/ITS3/alignment/src/AlignmentParams.cxx b/Detectors/Upgrades/ITS3/alignment/src/AlignmentParams.cxx new file mode 100644 index 0000000000000..0d89cb4d4cffd --- /dev/null +++ b/Detectors/Upgrades/ITS3/alignment/src/AlignmentParams.cxx @@ -0,0 +1,13 @@ +// 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. + +#include "ITS3Align/AlignmentParams.h" +O2ParamImpl(o2::its3::align::AlignmentParams); diff --git a/Detectors/Upgrades/ITS3/alignment/src/AlignmentSensors.cxx b/Detectors/Upgrades/ITS3/alignment/src/AlignmentSensors.cxx new file mode 100644 index 0000000000000..b436ffaf9a6d5 --- /dev/null +++ b/Detectors/Upgrades/ITS3/alignment/src/AlignmentSensors.cxx @@ -0,0 +1,203 @@ +// Copyright 2019-2026 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. + +#include + +#include +#include + +#include "ITSMFTBase/SegmentationAlpide.h" +#include "ITS3Align/AlignmentSensors.h" +#include "ITSBase/GeometryTGeo.h" + +namespace o2::its3::align +{ + +AlignableVolume::Ptr buildHierarchyITS(AlignableVolume::SensorMapping& sensorMap) +{ + uint32_t gLbl{0}, det{0}; + auto geom = o2::its::GeometryTGeo::Instance(); + AlignableVolume *volHB{nullptr}, *volSt{nullptr}, *volHSt{nullptr}, *volMod{nullptr}; + std::unordered_map sym2vol; + + auto root = std::make_unique(geom->composeSymNameITS(), gLbl++, det, false, RigidBodyDOFNone); + sym2vol[root->getSymName()] = root.get(); + for (int ilr = 0; ilr < geom->getNumberOfLayers(); ilr++) { + for (int ihb = 0; ihb < geom->getNumberOfHalfBarrels(); ihb++) { + volHB = root->addChild(geom->composeSymNameHalfBarrel(ilr, ihb), gLbl++, det, false, RigidBodyDOFNone); + sym2vol[volHB->getSymName()] = volHB; + int nstavesHB = geom->getNumberOfStaves(ilr) / 2; + for (int ist = 0; ist < nstavesHB; ist++) { + volSt = volHB->addChild(geom->composeSymNameStave(ilr, ihb, ist), gLbl++, det, false, RigidBodyDOFNone); + sym2vol[volSt->getSymName()] = volSt; + for (int ihst = 0; ihst < geom->getNumberOfHalfStaves(ilr); ihst++) { + volHSt = volSt->addChild(geom->composeSymNameHalfStave(ilr, ihb, ist, ihst), gLbl++, det, false, RigidBodyDOFNone); + sym2vol[volHSt->getSymName()] = volHSt; + for (int imd = 0; imd < geom->getNumberOfModules(ilr); imd++) { + volMod = volHSt->addChild(geom->composeSymNameModule(ilr, ihb, ist, ihst, imd), gLbl++, det, false, RigidBodyDOFNone); + sym2vol[volMod->getSymName()] = volMod; + } + } + } + } + } + + // NOTE: for ITS sensors the local x and y are swapped + int lay = 0, hba = 0, sta = 0, ssta = 0, modd = 0, chip = 0; + for (int ich = 0; ich < geom->getNumberOfChips(); ich++) { + geom->getChipId(ich, lay, hba, sta, ssta, modd, chip); + GlobalLabel lbl(det, ich, true); + AlignableVolume* parVol = sym2vol[modd < 0 ? geom->composeSymNameStave(lay, hba, sta) : geom->composeSymNameModule(lay, hba, sta, ssta, modd)]; + if (!parVol) { + throw std::runtime_error(fmt::format("did not find parent for chip {}", ich)); + } + int nch = modd < 0 ? geom->getNumberOfChipsPerStave(lay) : geom->getNumberOfChipsPerModule(lay); + int jch = ich % nch; + if (ich == 148 || ich == 147) { + RigidBodyDOFMask m = DOF_BIT(RigidBodyDOF::TZ); + parVol->setRigidBodyDOF(m); + sensorMap[lbl] = parVol->addChild(geom->composeSymNameChip(lay, hba, sta, ssta, modd, jch), lbl, m); + } else { + sensorMap[lbl] = parVol->addChild(geom->composeSymNameChip(lay, hba, sta, ssta, modd, jch), lbl, RigidBodyDOFNone); + } + } + return root; +} + +AlignableVolume::Ptr buildHierarchyIT3(AlignableVolume::SensorMapping& sensorMap) +{ + uint32_t gLbl{0}, det{0}; + auto geom = o2::its::GeometryTGeo::Instance(); + AlignableVolume *volHB{nullptr}, *volSt{nullptr}, *volHSt{nullptr}, *volMod{nullptr}; + std::unordered_map sym2vol; + + auto root = std::make_unique(geom->composeSymNameITS(), gLbl++, det, false, RigidBodyDOFNone); + sym2vol[root->getSymName()] = root.get(); + for (int ilr = 0; ilr < geom->getNumberOfLayers(); ilr++) { + const bool isLayITS3 = (ilr < 3); + for (int ihb = 0; ihb < geom->getNumberOfHalfBarrels(); ihb++) { + volHB = root->addChild(geom->composeSymNameHalfBarrel(ilr, ihb, isLayITS3), gLbl++, det, false, RigidBodyDOFNone); + sym2vol[volHB->getSymName()] = volHB; + if (isLayITS3) { // for its3 there are no more alignable volumes after this! + continue; + } + int nstavesHB = geom->getNumberOfStaves(ilr) / 2; + for (int ist = 0; ist < nstavesHB; ist++) { + volSt = volHB->addChild(geom->composeSymNameStave(ilr, ihb, ist), gLbl++, det, false, RigidBodyDOFNone); + sym2vol[volSt->getSymName()] = volSt; + for (int ihst = 0; ihst < geom->getNumberOfHalfStaves(ilr); ihst++) { + volHSt = volSt->addChild(geom->composeSymNameHalfStave(ilr, ihb, ist, ihst), gLbl++, det, false, RigidBodyDOFNone); + sym2vol[volHSt->getSymName()] = volHSt; + for (int imd = 0; imd < geom->getNumberOfModules(ilr); imd++) { + volMod = volHSt->addChild(geom->composeSymNameModule(ilr, ihb, ist, ihst, imd), gLbl++, det, false, RigidBodyDOFNone); + sym2vol[volMod->getSymName()] = volMod; + } + } + } + } + } + + int lay = 0, hba = 0, sta = 0, ssta = 0, modd = 0, chip = 0; + for (int ich = 0; ich < geom->getNumberOfChips(); ich++) { + geom->getChipId(ich, lay, hba, sta, ssta, modd, chip); + const bool isLayITS3 = (lay < 3); + GlobalLabel lbl(det, ich, true); + if (isLayITS3) { + // ITS3 chips by construction do not have any DOFs still add them to have the measurment to alignable layer relation + AlignableVolume* parVol = sym2vol[geom->composeSymNameHalfBarrel(lay, hba, true)]; + if (!parVol) { + throw std::runtime_error(fmt::format("did not find parent for chip {}", ich)); + } + sensorMap[lbl] = parVol->addChild(geom->composeSymNameChip(lay, hba, sta, ssta, modd, chip, true), lbl, RigidBodyDOFPseudo); + if (o2::its3::constants::detID::getSensorID(ich) == 2) { + parVol->setRigidBodyDOF(RigidBodyDOFAll); + } + } else { + AlignableVolume* parVol = sym2vol[modd < 0 ? geom->composeSymNameStave(lay, hba, sta) : geom->composeSymNameModule(lay, hba, sta, ssta, modd)]; + if (!parVol) { + throw std::runtime_error(fmt::format("did not find parent for chip {}", ich)); + } + int nch = modd < 0 ? geom->getNumberOfChipsPerStave(lay) : geom->getNumberOfChipsPerModule(lay); + int jch = ich % nch; + sensorMap[lbl] = parVol->addChild(geom->composeSymNameChip(lay, hba, sta, ssta, modd, jch), lbl, RigidBodyDOFNone); + } + } + return root; +} + +void AlignableSensorITS::defineMatrixL2G() +{ + // the chip volume is not the measurment plane, need to correct for the epitaxial layer + const auto* chipL2G = mPN->GetMatrix(); + mL2G = *chipL2G; + double delta = itsmft::SegmentationAlpide::SensorLayerThickness - itsmft::SegmentationAlpide::SensorLayerThicknessEff; + TGeoTranslation tra(0., 0.5 * delta, 0.); + mL2G *= tra; +} + +void AlignableSensorITS::defineMatrixT2L() +{ + double locA[3] = {-100., 0., 0.}, locB[3] = {100., 0., 0.}, gloA[3], gloB[3]; + mL2G.LocalToMaster(locA, gloA); + mL2G.LocalToMaster(locB, gloB); + double dx = gloB[0] - gloA[0], dy = gloB[1] - gloA[1]; + double t = (gloB[0] * dx + gloB[1] * dy) / (dx * dx + dy * dy); + double xp = gloB[0] - (dx * t), yp = gloB[1] - (dy * t); + double alp = std::atan2(yp, xp); + o2::math_utils::bringTo02Pid(alp); + mT2L.RotateZ(alp * TMath::RadToDeg()); // mT2L before is identity and afterwards rotated + const TGeoHMatrix l2gI = mL2G.Inverse(); + mT2L.MultiplyLeft(l2gI); +} + +void AlignableSensorITS::computeJacobianL2T(const double* posLoc, Matrix66& jac) const +{ + jac.setZero(); + Eigen::Map> rotT2L(mT2L.GetRotationMatrix()); + Eigen::Matrix3d skew, rotL2T = rotT2L.transpose(); + skew << 0, -posLoc[2], posLoc[1], posLoc[2], 0, -posLoc[0], -posLoc[1], posLoc[0], 0; + jac.topLeftCorner<3, 3>() = rotL2T; + jac.topRightCorner<3, 3>() = -rotL2T * skew; + jac.bottomRightCorner<3, 3>() = rotL2T; +} + +void AlignableSensorIT3::defineMatrixL2G() +{ + mL2G = *mPN->GetMatrix(); +} + +void AlignableSensorIT3::defineMatrixT2L() +{ + double locA[3] = {-100., 0., 0.}, locB[3] = {100., 0., 0.}, gloA[3], gloB[3]; + mL2G.LocalToMaster(locA, gloA); + mL2G.LocalToMaster(locB, gloB); + double dx = gloB[0] - gloA[0], dy = gloB[1] - gloA[1]; + double t = (gloB[0] * dx + gloB[1] * dy) / (dx * dx + dy * dy); + double xp = gloB[0] - (dx * t), yp = gloB[1] - (dy * t); + double alp = std::atan2(yp, xp); + o2::math_utils::bringTo02Pid(alp); + mT2L.RotateZ(alp * TMath::RadToDeg()); + const TGeoHMatrix l2gI = mL2G.Inverse(); + mT2L.MultiplyLeft(l2gI); +} + +void AlignableSensorIT3::computeJacobianL2T(const double* posLoc, Matrix66& jac) const +{ + jac.setZero(); + Eigen::Map> rotT2L(mT2L.GetRotationMatrix()); + Eigen::Matrix3d skew, rotL2T = rotT2L.transpose(); + skew << 0, -posLoc[2], posLoc[1], posLoc[2], 0, -posLoc[0], -posLoc[1], posLoc[0], 0; + jac.topLeftCorner<3, 3>() = rotL2T; + jac.topRightCorner<3, 3>() = -rotL2T * skew; + jac.bottomRightCorner<3, 3>() = rotL2T; +} + +} // namespace o2::its3::align diff --git a/Detectors/Upgrades/ITS3/alignment/src/AlignmentSpec.cxx b/Detectors/Upgrades/ITS3/alignment/src/AlignmentSpec.cxx new file mode 100644 index 0000000000000..5469766ba6aed --- /dev/null +++ b/Detectors/Upgrades/ITS3/alignment/src/AlignmentSpec.cxx @@ -0,0 +1,810 @@ +// Copyright 2019-2026 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. + +#include +#include +#include + +#ifdef WITH_OPENMP +#include +#endif + +#include +#include +#include +#include +#include +#include + +#include "Framework/ConfigParamRegistry.h" +#include "Framework/DataProcessorSpec.h" +#include "Framework/Task.h" +#include "ITSBase/GeometryTGeo.h" +#include "CommonUtils/EnumFlags.h" +#include "DataFormatsGlobalTracking/RecoContainer.h" +#include "DetectorsCommonDataFormats/DetID.h" +#include "DetectorsBase/Propagator.h" +#include "DetectorsBase/GRPGeomHelper.h" +#include "ReconstructionDataFormats/PrimaryVertex.h" +#include "Steer/MCKinematicsReader.h" +#include "CommonUtils/TreeStreamRedirector.h" +#include "ReconstructionDataFormats/VtxTrackRef.h" +#include "ITS3Reconstruction/TopologyDictionary.h" +#include "DataFormatsITSMFT/TopologyDictionary.h" +#include "ITStracking/MathUtils.h" +#include "ITS3Reconstruction/IOUtils.h" +#include "ITS3Align/AlignmentParams.h" +#include "ITS3Align/AlignmentTypes.h" +#include "ITS3Align/AlignmentHierarchy.h" +#include "ITS3Align/AlignmentSensors.h" + +namespace o2::its3::align +{ +using namespace o2::framework; +using DetID = o2::detectors::DetID; +using DataRequest = o2::globaltracking::DataRequest; +using PVertex = o2::dataformats::PrimaryVertex; +using V2TRef = o2::dataformats::VtxTrackRef; +using VTIndex = o2::dataformats::VtxTrackIndex; +using GTrackID = o2::dataformats::GlobalTrackID; +using TrackD = o2::track::TrackParCovD; + +class AlignmentSpec final : public Task +{ + public: + ~AlignmentSpec() final = default; + AlignmentSpec(const AlignmentSpec&) = delete; + AlignmentSpec(AlignmentSpec&&) = delete; + AlignmentSpec& operator=(const AlignmentSpec&) = delete; + AlignmentSpec& operator=(AlignmentSpec&&) = delete; + AlignmentSpec(std::shared_ptr dr, std::shared_ptr gr, GTrackID::mask_t src, bool useMC, bool withPV) + : mDataRequest(dr), mGGCCDBRequest(gr), mTracksSrc(src), mUseMC(useMC), mWithPV(withPV) + { + } + + void init(InitContext& ic) final; + void run(ProcessingContext& pc) final; + void endOfStream(EndOfStreamContext& ec) final; + void finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) final; + void process(); + + private: + void updateTimeDependentParams(ProcessingContext& pc); + void buildHierarchy(); + + // calculate the transport jacobian for points FROM and TO numerically via ridder's method + // this assumes the track is already at point FROM and will be extrapolated to TO's x (xTo) + // method does not modify the original track + bool getTransportJacobian(const TrackD& track, double xTo, double alphaTo, gbl::Matrix5d& jac, gbl::Matrix5d& err); + + // refit ITS track with inward/outward fit (opt. impose pv as additional constraint) + // after this we have the refitted track at the innermost update point + bool prepareITSTrack(int iTrk, const o2::its::TrackITS& itsTrack, Track& resTrack); + + // prepare ITS measuremnt points + void prepareMeasurments(std::span clusters, std::span pattIt); + + // build track to vertex association + void buildT2V(); + + // create a track in double precision from single precision + TrackD convertTrack(const track::TrackParCov& trk) const noexcept; + + // Steering debug + enum class OutputOpt : uint8_t { + VerboseGBL = 0, + MilleData, + MilleSteer, + Debug, + }; + utils::EnumFlags mOutOpt; + std::unique_ptr mDBGOut; + std::vector mPVs; + std::vector mT2PV; + const o2::itsmft::TopologyDictionary* mITSDict{nullptr}; + const o2::its3::TopologyDictionary* mIT3Dict{nullptr}; + o2::globaltracking::RecoContainer* mRecoData = nullptr; + std::unique_ptr mcReader; + std::vector mITSTrackingInfo; + std::shared_ptr mDataRequest; + std::shared_ptr mGGCCDBRequest; + std::unique_ptr mHierarchy; // tree-hiearchy + AlignableVolume::SensorMapping mChip2Hiearchy; // global label mapping to leaves in the tree + bool mUseMC{false}; + bool mWithPV{false}; + GTrackID::mask_t mTracksSrc; + int mNThreads{1}; +}; + +void AlignmentSpec::init(InitContext& ic) +{ + o2::base::GRPGeomHelper::instance().setRequest(mGGCCDBRequest); + mNThreads = ic.options().get("nthreads"); + mOutOpt.set(ic.options().get("output")); + if (mOutOpt) { + LOG(info) << mOutOpt.pstring(); + mDBGOut = std::make_unique("its3_debug_alg.root", "recreate"); + } + if (mUseMC) { + mcReader = std::make_unique("collisioncontext.root"); + } +} + +void AlignmentSpec::run(ProcessingContext& pc) +{ + o2::globaltracking::RecoContainer recoData; + mRecoData = &recoData; + mRecoData->collectData(pc, *mDataRequest); + updateTimeDependentParams(pc); + process(); + mRecoData = nullptr; +} + +void AlignmentSpec::process() +{ + if (!mITSDict && !mIT3Dict) { + LOGP(fatal, "ITS data is not loaded"); + } + const auto& param = AlignmentParams::Instance(); + auto prop = o2::base::PropagatorD::Instance(); + const auto bz = prop->getNominalBz(); + const auto itsTracks = mRecoData->getITSTracks(); + const auto itsClRefs = mRecoData->getITSTracksClusterRefs(); + const auto clusITS = mRecoData->getITSClusters(); + const auto patterns = mRecoData->getITSClustersPatterns(); + std::span mcLbls; + if (mUseMC) { + mcLbls = mRecoData->getITSTracksMCLabels(); + } + prepareMeasurments(clusITS, patterns); + + if (mWithPV) { + buildT2V(); + } + + LOGP(info, "Starting fits with {} threads", mNThreads); + + // Data + std::vector> gblTrajSlots(mNThreads); + std::vector> resTrackSlots(mNThreads); + + auto timeStart = std::chrono::high_resolution_clock::now(); + int cFailedRefit{0}, cFailedProp{0}, cSelected{0}, cGBLFit{0}, cGBLFitFail{0}, cGBLChi2Rej{0}, cGBLConstruct{0}; + double chi2Sum{0}, lostWeightSum{0}; + int ndfSum{0}; +#ifdef WITH_OPENMP +#pragma omp parallel num_threads(mNThreads) \ + reduction(+ : cFailedRefit) \ + reduction(+ : cFailedProp) \ + reduction(+ : cSelected) \ + reduction(+ : cGBLFit) \ + reduction(+ : cGBLFitFail) \ + reduction(+ : cGBLChi2Rej) \ + reduction(+ : cGBLConstruct) \ + reduction(+ : chi2Sum) \ + reduction(+ : lostWeightSum) \ + reduction(+ : ndfSum) +#endif + { +#ifdef WITH_OPENMP + const int tid = omp_get_thread_num(); +#else + const int tid = 0; +#endif + auto& gblTrajSlot = gblTrajSlots[tid]; + auto& resTrackSlot = resTrackSlots[tid]; + +#ifdef WITH_OPENMP +#pragma omp for schedule(dynamic) +#endif + for (size_t iTrk = 0; iTrk < (int)itsTracks.size(); ++iTrk) { + const auto& trk = itsTracks[iTrk]; + if (trk.getNClusters() < param.minITSCls || + (trk.getChi2() / ((float)trk.getNClusters() * 2 - 5)) >= param.maxITSChi2Ndf || + trk.getPt() < param.minPt || + (mUseMC && (!mcLbls[iTrk].isValid() || !mcLbls[iTrk].isCorrect()))) { + continue; + } + ++cSelected; + Track& resTrack = resTrackSlot.emplace_back(); + if (!prepareITSTrack((int)iTrk, trk, resTrack)) { + ++cFailedRefit; + resTrackSlot.pop_back(); + continue; + } + + o2::track::TrackParD* refLin = nullptr; + if (param.useStableRefit) { + refLin = &resTrack.track; + } + + // outward stepping from track IU + auto wTrk = resTrack.track; + const bool hasPV = resTrack.info[0].lr == -1; + std::vector points; + bool first = true, failed = false; + const int np = (int)resTrack.points.size(); + track::TrackLTIntegral lt; + lt.setTimeNotNeeded(); + constexpr int perm[5] = {4, 2, 3, 0, 1}; // ALICE->GBL: Q/Pt,Snp,Tgl,Y,Z + + // calcualte seed precision from IU covariance + // auto covALICE = resTrack.track.getCov(); + // gbl::Matrix5d seedCov, seedPrec; + // for (int i = 0; i < 5; i++) { + // for (int j = 0; j < 5; j++) { + // // getCov returns lower triangle in packed form + // int iA = perm[i], jA = perm[j]; + // if (iA < jA) { + // std::swap(iA, jA); + // } + // seedCov(i, j) = covALICE[(iA * (iA + 1) / 2) + jA]; + // } + // } + // seedPrec = seedCov.inverse(); + + for (int ip{0}; ip < np; ++ip) { + const auto& frame = resTrack.info[ip]; + gbl::Matrix5d err = gbl::Matrix5d::Identity(), jacALICE = gbl::Matrix5d::Identity(), jacGBL; + float msErr = 0.f; + if (!first) { + // numerically calculates the transport jacobian from prev. point to this point + // then we actually do the step to the point and accumulate the material + if (!getTransportJacobian(wTrk, frame.x, frame.alpha, jacALICE, err) || + !prop->propagateToAlphaX(wTrk, refLin, frame.alpha, frame.x, false, param.maxSnp, param.maxStep, 1, param.matCorrType, <)) { + ++cFailedProp; + failed = true; + break; + } + msErr = its::math_utils::MSangle(trk.getPID().getMass(), trk.getP(), lt.getX2X0()); + } else { + first = false; + } + + // after computing jac, reorder to GBL convention + for (int i = 0; i < 5; i++) { + for (int j = 0; j < 5; j++) { + jacGBL(i, j) = jacALICE(perm[i], perm[j]); + } + } + + // wTrk is now in the measurment frame + gbl::GblPoint point(jacGBL); + // measurement + Eigen::Vector2d res, prec; + res << frame.positionTrackingFrame[0] - wTrk.getY(), frame.positionTrackingFrame[1] - wTrk.getZ(); + prec << 1. / resTrack.points[ip].sig2y, 1. / resTrack.points[ip].sig2z; + // the projection matrix is in the tracking frame the idendity so no need to diagonalize it + point.addMeasurement(res, prec); + if (msErr > param.minMS && ip < np - 1) { + Eigen::Vector2d scat(0., 0.), scatPrec = Eigen::Vector2d::Constant(1. / (msErr * msErr)); + point.addScatterer(scat, scatPrec); + lt.clearFast(); // only clear if accounted + } + + if (frame.lr >= 0) { + GlobalLabel lbl(0, frame.sens, true); + if (mChip2Hiearchy.find(lbl) == mChip2Hiearchy.end()) { + LOGP(fatal, "Cannot find global label: {}", lbl.asString()); + } + + // derivatives for all sensitive volumes and their parents + // this is the derivative in TRK but we want to align in LOC + // so dr/da_(LOC) = dr/da_(TRK) * da_(TRK)/da_(LOC) + const auto* sensorVol = mChip2Hiearchy.at(lbl); + Matrix36 der = getRigidBodyDerivatives(wTrk); + + // count columns: only volumes with real DOFs (not DOFPseudo) + int nCol{0}; + for (const auto* v = sensorVol; v && !v->isRoot(); v = v->getParent()) { + if (v->getRigidBodyDOF() != RigidBodyDOFPseudo) { + nCol += RigidBodyDOF::NDOF; + } + } + std::vector gLabels; + gLabels.reserve(nCol); + Eigen::MatrixXd gDer(3, nCol); + Eigen::Index curCol{0}; + + // 1) sensor: TRK -> LOC via precomputed T2L and J_L2T + const double posTrk[3] = {frame.x, 0., 0.}; + double posLoc[3]; + sensorVol->getT2L().LocalToMaster(posTrk, posLoc); + Matrix66 jacL2T; + sensorVol->computeJacobianL2T(posLoc, jacL2T); + der *= jacL2T; + if (sensorVol->getRigidBodyDOF() != RigidBodyDOFPseudo) { + for (uint8_t iDOF{0}; iDOF < RigidBodyDOF::NDOF; ++iDOF) { + gLabels.push_back(sensorVol->getLabel().rawGBL(iDOF)); + } + gDer.middleCols(curCol++ * RigidBodyDOF::NDOF, RigidBodyDOF::NDOF) = der; + } + + // 2) chain through parents: child's J_L2P + for (const auto* child = sensorVol; child->getParent() && !child->getParent()->isRoot(); child = child->getParent()) { + der *= child->getJL2P(); + const auto* parent = child->getParent(); + if (parent->getRigidBodyDOF() != RigidBodyDOFPseudo) { + for (uint8_t iDOF{0}; iDOF < RigidBodyDOF::NDOF; ++iDOF) { + gLabels.push_back(parent->getLabel().rawGBL(iDOF)); + } + gDer.middleCols(curCol++ * RigidBodyDOF::NDOF, RigidBodyDOF::NDOF) = der; + } + } + point.addGlobals(gLabels, gDer); + } + + if (mOutOpt[OutputOpt::VerboseGBL]) { + static Eigen::IOFormat fmt(4, 0, ", ", "\n", "[", "]"); + LOGP(info, "WORKING-POINT {}", ip); + LOGP(info, "Track: {}", wTrk.asString()); + LOGP(info, "FrameInfo: {}", frame.asString()); + std::cout << "jacALICE:\n" + << jacALICE.format(fmt) << '\n'; + std::cout << "jacGBL:\n" + << jacGBL.format(fmt) << '\n'; + LOGP(info, "Point {}: GBL res=({}, {}), KF stored res=({}, {})", + ip, res[0], res[1], resTrack.points[ip].dy, resTrack.points[ip].dz); + LOGP(info, "residual: dy={} dz={}", res[0], res[1]); + LOGP(info, "precision: precY={} precZ={}", prec[0], prec[1]); + point.printPoint(5); + } + points.push_back(point); + } + if (!failed) { + // TODO: check why external seed worsens fit? + // gbl::GblTrajectory traj(points, 1, seedPrec, std::abs(bz) > 0.001); + gbl::GblTrajectory traj(points, std::abs(bz) > 0.001); + if (traj.isValid()) { + double chi2 = NAN, lostWeight = NAN; + int ndf = 0; + if (auto ierr = traj.fit(chi2, ndf, lostWeight); !ierr) { + if (mOutOpt[OutputOpt::VerboseGBL]) { + LOGP(info, "GBL FIT chi2 {} ndf {}", chi2, ndf); + traj.printTrajectory(5); + } + if (chi2 / ndf > param.maxChi2Ndf) { + LOGP(error, "GBL fit exceeded red chi2 {}", chi2 / ndf); + ++cGBLChi2Rej; + if (std::abs(resTrack.kfFit.chi2Ndf - 1) < 0.02) { + LOGP(error, "\tGBL is far away from KF fit!!!!"); + continue; + } + } else { + ++cGBLFit; + chi2Sum += chi2; + lostWeightSum += lostWeight; + ndfSum += ndf; + if (mOutOpt[OutputOpt::MilleData]) { + gblTrajSlot.push_back(traj); + } + FitInfo fit{.ndf = ndf, .chi2 = (float)chi2, .chi2Ndf = (float)chi2 / (float)ndf}; + resTrack.gblFit = fit; + } + } else { + ++cGBLFitFail; + } + } else { + ++cGBLConstruct; + } + } + } + } + auto timeEnd = std::chrono::high_resolution_clock::now(); + auto duration = std::chrono::duration_cast(timeEnd - timeStart); + LOGP(info, "Fitted {} tracks out of {} (selected {}) in {} sec", cGBLFit, itsTracks.size(), cSelected, duration.count() / 1e3); + LOGP(info, "\tRefit failed for {} tracks; Failed prop for {} tracks", cFailedRefit, cFailedProp); + LOGP(info, "\tGBL SUMMARY:"); + LOGP(info, "\t\tGBL construction failed {}", cGBLConstruct); + LOGP(info, "\t\tGBL fit failed {}", cGBLFitFail); + LOGP(info, "\t\tGBL chi2Ndf rejected {}", cGBLChi2Rej); + if (!ndfSum) { + LOGP(info, "\t\tGBL Chi2/Ndf = NDF IS 0"); + } else { + LOGP(info, "\t\tGBL Chi2/Ndf = {}", chi2Sum / ndfSum); + } + LOGP(info, "\t\tGBL LostWeight = {}", lostWeightSum); + LOGP(info, "Streaming results to output"); + if (mOutOpt[OutputOpt::MilleData]) { + gbl::MilleBinary mille(param.milleBinFile, true); + for (auto& slot : gblTrajSlots) { + for (auto& traj : slot) { + traj.milleOut(mille); + } + } + } + if (mOutOpt[OutputOpt::Debug]) { + for (auto& slot : resTrackSlots) { + for (auto& res : slot) { + (*mDBGOut) << "res" + << "trk=" << res + << "\n"; + } + } + } +} + +void AlignmentSpec::updateTimeDependentParams(ProcessingContext& pc) +{ + o2::base::GRPGeomHelper::instance().checkUpdates(pc); + if (static bool initOnce{false}; !initOnce) { + initOnce = true; + auto geom = o2::its::GeometryTGeo::Instance(); + o2::its::GeometryTGeo::Instance()->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::T2L, o2::math_utils::TransformType::L2G, o2::math_utils::TransformType::T2G)); + AlignmentParams::Instance().printKeyValues(true, true); + + buildHierarchy(); + } +} + +void AlignmentSpec::buildHierarchy() +{ + if (mIT3Dict) { + mHierarchy = buildHierarchyIT3(mChip2Hiearchy); + } else { + mHierarchy = buildHierarchyITS(mChip2Hiearchy); + } + + const auto& param = AlignmentParams::Instance(); + mHierarchy->finalise(); + if (mOutOpt[OutputOpt::MilleSteer]) { + std::ofstream tree(param.milleTreeFile); + mHierarchy->writeTree(tree); + std::ofstream cons(param.milleConFile); + mHierarchy->writeRigidBodyConstraints(cons); + std::ofstream par(param.milleParamFile); + mHierarchy->writeParameters(par); + } +} + +bool AlignmentSpec::getTransportJacobian(const TrackD& track, double xTo, double alphaTo, gbl::Matrix5d& jac, gbl::Matrix5d& err) +{ + auto prop = o2::base::PropagatorD::Instance(); + const auto bz = prop->getNominalBz(); + const auto& param = AlignmentParams::Instance(); + const auto minStep = std::sqrt(std::numeric_limits::epsilon()); + const gbl::Vector5d x0(track.getParams()); + auto trackC = track; + o2::track::TrackParD* refLin{nullptr}; + if (param.useStableRefit) { + refLin = &trackC; + } + + auto propagate = [&](gbl::Vector5d& p) -> bool { + TrackD tmp(track); + for (int i{0}; i < track::kNParams; ++i) { + tmp.setParam(p[i], i); + } + if (!prop->propagateToAlphaX(tmp, refLin, alphaTo, xTo, false, param.maxSnp, param.maxStep, 1, param.matCorrType)) { + return false; + } + p = gbl::Vector5d(tmp.getParams()); + return true; + }; + + for (int iPar{0}; iPar < track::kNParams; ++iPar) { + // step size + double h = std::min(param.ridderMaxIniStep[iPar], std::max(minStep, std::abs(track.getParam(iPar)) * param.ridderRelIniStep[iPar]) * std::pow(param.ridderShrinkFac, param.ridderMaxExtrap / 2)); + ; + // romberg tableu + Eigen::MatrixXd cur(track::kNParams, param.ridderMaxExtrap); + Eigen::MatrixXd pre(track::kNParams, param.ridderMaxExtrap); + double normErr = std::numeric_limits::max(); + gbl::Vector5d bestDeriv = gbl::Vector5d::Constant(std::numeric_limits::max()); + for (int iExt{0}; iExt < param.ridderMaxExtrap; ++iExt) { + gbl::Vector5d xPlus = x0, xMinus = x0; + xPlus(iPar) += h; + xMinus(iPar) -= h; + if (!propagate(xPlus) || !propagate(xMinus)) { + return false; + } + cur.col(0) = (xPlus - xMinus) / (2.0 * h); + if (!iExt) { + bestDeriv = cur.col(0); + } + // shrink step in next iteration + h /= param.ridderShrinkFac; + // richardson extrapolation + double fac = param.ridderShrinkFac * param.ridderShrinkFac; + for (int k{1}; k <= iExt; ++k) { + cur.col(k) = (fac * cur.col(k - 1) - pre.col(k - 1)) / (fac - 1.0); + fac *= param.ridderShrinkFac * param.ridderShrinkFac; + double e = std::max((cur.col(k) - cur.col(k - 1)).norm(), (cur.col(k) - pre.col(k - 1)).norm()); + if (e <= normErr) { + normErr = e; + bestDeriv = cur.col(k); + if (normErr < param.ridderEps) { + break; + } + } + } + if (normErr < param.ridderEps) { + break; + } + // check stability + if (iExt > 0) { + double tableauErr = (cur.col(iExt) - pre.col(iExt - 1)).norm(); + if (tableauErr >= 2.0 * normErr) { + break; + } + } + std::swap(cur, pre); + } + if (bestDeriv.isApproxToConstant(std::numeric_limits::max())) { + return false; + } + jac.col(iPar) = bestDeriv; + err.col(iPar) = gbl::Vector5d::Constant(normErr); + } + + if (jac.isIdentity(1e-8)) { + LOGP(error, "Near jacobian idendity for taking track from {} to {}", track.getX(), xTo); + return false; + } + // only check for elements where there is no change expected q/pt, tgl + if (std::abs(1. - jac(4, 4)) > param.ridderMaxJacDiagTol || std::abs(1. - jac(3, 3)) > param.ridderMaxJacDiagTol) { + LOGP(error, "Diagonal element not expected to change in propagation jacobian is to far away from 1"); + LOG(info) << jac; + return false; + } + + return true; +} + +bool AlignmentSpec::prepareITSTrack(int iTrk, const o2::its::TrackITS& itsTrack, align::Track& resTrack) +{ + const auto& param = AlignmentParams::Instance(); + const auto itsClRefs = mRecoData->getITSTracksClusterRefs(); + auto trFit = convertTrack(itsTrack.getParamOut()); // take outer track fit as start of refit + auto prop = o2::base::PropagatorD::Instance(); + auto geom = o2::its::GeometryTGeo::Instance(); + const auto bz = prop->getNominalBz(); + std::array frameArr{}; + o2::track::TrackParD trkOut = trFit, *refLin = nullptr; + if (param.useStableRefit) { + refLin = &trkOut; + } + + auto accountCluster = [&](int i, TrackD& tr, float& chi2, Measurement& meas, o2::track::TrackParD* refLin) { + if (frameArr[i]) { // update with cluster + if (!prop->propagateToAlphaX(tr, refLin, frameArr[i]->alpha, frameArr[i]->x, false, param.maxSnp, param.maxStep, 1, param.matCorrType)) { + return 2; + } + meas.dy = frameArr[i]->positionTrackingFrame[0] - tr.getY(); + meas.dz = frameArr[i]->positionTrackingFrame[1] - tr.getZ(); + meas.sig2y = frameArr[i]->covarianceTrackingFrame[0]; + meas.sig2z = frameArr[i]->covarianceTrackingFrame[2]; + meas.z = tr.getZ(); + meas.phi = tr.getPhi(); + o2::math_utils::bringTo02Pid(meas.phi); + chi2 += (float)tr.getPredictedChi2Quiet(frameArr[i]->positionTrackingFrame, frameArr[i]->covarianceTrackingFrame); + if (!tr.update(frameArr[i]->positionTrackingFrame, frameArr[i]->covarianceTrackingFrame)) { + return 2; + } + if (refLin) { // displace the reference to the last updated cluster + refLin->setY(frameArr[i]->positionTrackingFrame[0]); + refLin->setZ(frameArr[i]->positionTrackingFrame[1]); + } + return 0; + } + return 1; + }; + + // add PV as constraint + FrameInfoExt* pvInfo{nullptr}; + if (mWithPV) { + const int iPV = mT2PV[iTrk]; + if (iPV < 0) { + return false; + } + const auto& pv = mPVs[iPV]; + auto tmp = convertTrack(itsTrack.getParamIn()); + if (!prop->propagateToDCA(pv, tmp, bz)) { + return false; + } + pvInfo = new FrameInfoExt; + pvInfo->alpha = (float)tmp.getAlpha(); + float ca{0}, sa{0}; + o2::math_utils::bringTo02Pi(pvInfo->alpha); + o2::math_utils::sincos(pvInfo->alpha, sa, ca); + pvInfo->x = tmp.getX(); + pvInfo->positionTrackingFrame[0] = -pv.getX() * sa + pv.getY() * ca; + pvInfo->positionTrackingFrame[1] = pv.getZ(); + pvInfo->covarianceTrackingFrame[0] = 0.5 * (pv.getSigmaX2() + pv.getSigmaY2()); + pvInfo->covarianceTrackingFrame[2] = pv.getSigmaY2(); + pvInfo->sens = -1; + pvInfo->lr = -1; + } + frameArr[0] = pvInfo; + + // collect all track clusters to array, placing them to layer+1 slot + int nCl = itsTrack.getNClusters(); + for (int i = 0; i < nCl; i++) { // clusters are ordered from the outermost to the innermost + const auto& curInfo = mITSTrackingInfo[itsClRefs[itsTrack.getClusterEntry(i)]]; + frameArr[1 + curInfo.lr] = &curInfo; + } + + // start refit + resTrack.points.clear(); + resTrack.info.clear(); + trFit.resetCovariance(); + trFit.setCov(trFit.getQ2Pt() * trFit.getQ2Pt() * trFit.getCov()[14], 14); + float chi2{0}; + for (int i{7}; i >= 0; --i) { + Measurement point; + int res = accountCluster(i, trFit, chi2, point, refLin); + if (res == 2) { + return false; + } else if (res == 0) { + resTrack.points.push_back(point); + resTrack.info.push_back(*frameArr[i]); + resTrack.track = trFit; // put track to whatever the IU is + } + } + // reverse inserted points so they are in the same order as the track + std::reverse(resTrack.info.begin(), resTrack.info.end()); + std::reverse(resTrack.points.begin(), resTrack.points.end()); + resTrack.kfFit.chi2 = chi2; + resTrack.kfFit.ndf = (int)resTrack.info.size() * 2 - 5; + resTrack.kfFit.chi2Ndf = chi2 / (float)resTrack.kfFit.ndf; + + delete pvInfo; + + return true; +} + +void AlignmentSpec::prepareMeasurments(std::span clusters, std::span patterns) +{ + LOGP(info, "Preparing {} measurments", clusters.size()); + auto geom = its::GeometryTGeo::Instance(); + geom->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::T2L, o2::math_utils::TransformType::L2G)); + mITSTrackingInfo.clear(); + mITSTrackingInfo.reserve(clusters.size()); + auto pattIt = patterns.begin(); + for (const auto& cls : clusters) { + const auto sens = cls.getSensorID(); + const auto lay = geom->getLayer(sens); + float sigmaY2{0}, sigmaZ2{0}; + math_utils::Point3D locXYZ; + if (mIT3Dict) { + locXYZ = o2::its3::ioutils::extractClusterData(cls, pattIt, mIT3Dict, sigmaY2, sigmaZ2); + } else { + locXYZ = o2::its::ioutils::extractClusterData(cls, pattIt, mITSDict, sigmaY2, sigmaZ2); + } + // Transformation to the local --> global + const auto gloXYZ = geom->getMatrixL2G(sens) * locXYZ; + // Inverse transformation to the local --> tracking + const o2::math_utils::Point3D trkXYZ = geom->getMatrixT2L(sens) ^ locXYZ; + // Tracking alpha angle + // We want that each cluster rotates its tracking frame to the clusters phi + // that way the track linearization around the measurement is less biases to the arc + // this means automatically that the measurement on the arc is at 0 for the curved layers + float alpha = geom->getSensorRefAlpha(sens); + float x = trkXYZ.x(), y = trkXYZ.y(); + if (mIT3Dict && lay < 3) { + y = 0.f; + // alpha&x always have to be defined wrt to the global Z axis! + x = std::hypot(gloXYZ.x(), gloXYZ.y()); + alpha = std::atan2(gloXYZ.y(), gloXYZ.x()); + } + mITSTrackingInfo.emplace_back(sens, lay, trkXYZ, locXYZ, gloXYZ, x, alpha, + std::array{y, trkXYZ.z()}, + std::array{sigmaY2, 0., sigmaZ2}); + } +} + +void AlignmentSpec::buildT2V() +{ + const auto& itsTracks = mRecoData->getITSTracks(); + mT2PV.clear(); + mT2PV.resize(itsTracks.size(), -1); + if (mUseMC) { + mPVs.reserve(mcReader->getNEvents(0)); + for (int iEve{0}; iEve < mcReader->getNEvents(0); ++iEve) { + const auto& eve = mcReader->getMCEventHeader(0, iEve); + dataformats::VertexBase vtx; + constexpr float err{30e-7f}; + vtx.setX((float)eve.GetX() + (float)gRandom->Gaus(0.f, err)); + vtx.setY((float)eve.GetY() + (float)gRandom->Gaus(0.f, err)); + vtx.setZ((float)eve.GetZ() + (float)gRandom->Gaus(0.f, err)); + vtx.setCov(err, 0.f, 0.f, err, 0.f, err); + mPVs.push_back(vtx); + } + const auto& mcLbls = mRecoData->getITSTracksMCLabels(); + for (size_t iTrk{0}; iTrk < mcLbls.size(); ++iTrk) { + const auto& lbl = mcLbls[iTrk]; + if (!lbl.isValid() || !lbl.isCorrect()) { + continue; + } + const auto& mcTrk = mcReader->getTrack(lbl); + if (mcTrk->isPrimary()) { + mT2PV[iTrk] = lbl.getEventID(); + } + } + } else { + LOGP(fatal, "Data PV to track TODO"); + } +} + +TrackD AlignmentSpec::convertTrack(const track::TrackParCov& trk) const noexcept +{ + TrackD dst; + dst.setX(trk.getX()); + dst.setAlpha(trk.getAlpha()); + for (int iPar{0}; iPar < track::kNParams; ++iPar) { + dst.setParam(trk.getParam(iPar), iPar); + } + dst.setAbsCharge(trk.getAbsCharge()); + dst.setPID(trk.getPID()); + dst.setUserField(trk.getUserField()); + for (int iCov{0}; iCov < track::kCovMatSize; ++iCov) { + dst.setCov(trk.getCov()[iCov], iCov); + } + return dst; +} + +void AlignmentSpec::endOfStream(EndOfStreamContext& /*ec*/) +{ + mDBGOut->Close(); + mDBGOut.reset(); +} + +void AlignmentSpec::finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) +{ + if (o2::base::GRPGeomHelper::instance().finaliseCCDB(matcher, obj)) { + return; + } + if (matcher == ConcreteDataMatcher("ITS", "CLUSDICT", 0)) { + LOG(info) << "its cluster dictionary updated"; + mITSDict = (const o2::itsmft::TopologyDictionary*)obj; + return; + } + if (matcher == ConcreteDataMatcher("IT3", "CLUSDICT", 0)) { + LOG(info) << "it3 cluster dictionary updated"; + mIT3Dict = (const o2::its3::TopologyDictionary*)obj; + return; + } +} + +DataProcessorSpec getAlignmentSpec(GTrackID::mask_t srcTracks, GTrackID::mask_t srcClusters, bool useMC, bool withPV, bool withITS) +{ + auto dataRequest = std::make_shared(); + dataRequest->requestTracks(srcTracks, useMC); + if (!withITS) { + dataRequest->requestIT3Clusters(useMC); + } else { + dataRequest->requestClusters(srcClusters, useMC); + } + if (withPV && !useMC) { + dataRequest->requestPrimaryVertices(useMC); + } + auto ggRequest = std::make_shared(false, // orbitResetTime + true, // GRPECS=true + true, // GRPLHCIF + true, // GRPMagField + true, // askMatLUT + o2::base::GRPGeomRequest::Aligned, // geometry + dataRequest->inputs, // inputs + true, // askOnce + true); // propagatorD + Options opts{ + {"nthreads", VariantType::Int, 1, {"number of threads"}}, + {"output", VariantType::String, "", {"output flags"}}, + }; + + return DataProcessorSpec{ + .name = "its3-alignment", + .inputs = dataRequest->inputs, + .outputs = {}, + .algorithm = AlgorithmSpec{adaptFromTask(dataRequest, ggRequest, srcTracks, useMC, withPV)}, + .options = opts}; +} +} // namespace o2::its3::align diff --git a/Detectors/Upgrades/ITS3/alignment/src/AlignmentTypes.cxx b/Detectors/Upgrades/ITS3/alignment/src/AlignmentTypes.cxx new file mode 100644 index 0000000000000..2231a715cdd12 --- /dev/null +++ b/Detectors/Upgrades/ITS3/alignment/src/AlignmentTypes.cxx @@ -0,0 +1,24 @@ +// 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. + +#include +#include + +#include "ITS3Align/AlignmentTypes.h" +ClassImp(o2::its3::align::Point); +ClassImp(o2::its3::align::FrameInfoExt); +ClassImp(o2::its3::align::FitInfo); +ClassImp(o2::its3::align::Track); + +std::string o2::its3::align::FrameInfoExt::asString() const +{ + return std::format("Sensor={} Layer={} TrkX={} Alpha={}\n\tTRK: y={} z={}\n\tMEAS: y={} z={}", sens, lr, x, alpha, trk.y(), trk.z(), positionTrackingFrame[0], positionTrackingFrame[1]); +} diff --git a/Detectors/Upgrades/ITS3/alignment/src/ITS3AlignLinkDef.h b/Detectors/Upgrades/ITS3/alignment/src/ITS3AlignLinkDef.h index ef526284f3a58..527eeebeb82d5 100644 --- a/Detectors/Upgrades/ITS3/alignment/src/ITS3AlignLinkDef.h +++ b/Detectors/Upgrades/ITS3/alignment/src/ITS3AlignLinkDef.h @@ -17,4 +17,12 @@ #pragma link C++ class o2::its3::align::MisalignmentParameters + ; +#pragma link C++ struct o2::its3::align::AlignmentParams + ; +#pragma link C++ class o2::conf::ConfigurableParamHelper < o2::its3::align::AlignmentParams> + ; + +#pragma link C++ struct o2::its3::align::Measurement + ; +#pragma link C++ struct o2::its3::align::FrameInfoExt + ; +#pragma link C++ struct o2::its3::align::FitInfo + ; +#pragma link C++ struct o2::its3::align::Track + ; + #endif diff --git a/Detectors/Upgrades/ITS3/alignment/src/alignment-workflow.cxx b/Detectors/Upgrades/ITS3/alignment/src/alignment-workflow.cxx new file mode 100644 index 0000000000000..29c0f380f319b --- /dev/null +++ b/Detectors/Upgrades/ITS3/alignment/src/alignment-workflow.cxx @@ -0,0 +1,64 @@ +// 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. + +#include "CommonUtils/ConfigurableParam.h" +#include "Framework/ConfigParamSpec.h" +#include "Framework/CallbacksPolicy.h" +#include "GlobalTrackingWorkflowHelpers/InputHelper.h" +#include "DetectorsRaw/HBFUtilsInitializer.h" +#include "ITS3Align/AlignmentSpec.h" + +using namespace o2::framework; +using GID = o2::dataformats::GlobalTrackID; +using DetID = o2::detectors::DetID; + +void customize(std::vector& policies) +{ + o2::raw::HBFUtilsInitializer::addNewTimeSliceCallback(policies); +} +void customize(std::vector& workflowOptions) +{ + std::vector options{ + {"disable-mc", o2::framework::VariantType::Bool, false, {"enable MC propagation"}}, + {"track-sources", VariantType::String, std::string{GID::ALL}, {"comma-separated list of track sources to use"}}, + {"cluster-sources", VariantType::String, "ITS", {"comma-separated list of cluster sources to use"}}, + {"with-its", VariantType::Bool, false, {"ITS alignment mode"}}, + {"without-pv", VariantType::Bool, false, {"Do not use in track refit the PV as an additional constraint"}}, + {"disable-root-input", VariantType::Bool, false, {"disable root-files input reader"}}, + {"configKeyValues", VariantType::String, "", {"Semicolon separated key=value strings ..."}}}; + o2::raw::HBFUtilsInitializer::addConfigOption(options); + std::swap(workflowOptions, options); +} +#include "Framework/runDataProcessing.h" + +WorkflowSpec defineDataProcessing(ConfigContext const& cfg) +{ + o2::conf::ConfigurableParam::updateFromString(cfg.options().get("configKeyValues")); + + const GID::mask_t allowedSourcesTrc = GID::getSourcesMask("ITS,TPC,ITS-TPC,ITS-TPC-TRD,ITS-TPC-TOF,ITS-TPC-TRD-TOF"); + const GID::mask_t allowedSourcesClus = GID::getSourcesMask("ITS"); + const GID::mask_t srcTrc = allowedSourcesTrc & GID::getSourcesMask(cfg.options().get("track-sources")); + const GID::mask_t srcCls = allowedSourcesClus & GID::getSourcesMask(cfg.options().get("cluster-sources")); + const auto useMC = !cfg.options().get("disable-mc"); + const auto withPV = !cfg.options().get("without-pv"); + const auto withITS = cfg.options().get("with-its"); + + WorkflowSpec specs; + o2::globaltracking::InputHelper::addInputSpecs(cfg, specs, srcCls, srcTrc, srcTrc, useMC); + if (withPV && !useMC) { + o2::globaltracking::InputHelper::addInputSpecsPVertex(cfg, specs, useMC); + } + + specs.emplace_back(o2::its3::align::getAlignmentSpec(srcTrc, srcCls, useMC, withPV, withITS)); + + o2::raw::HBFUtilsInitializer hbfIni(cfg, specs); + return std::move(specs); +} diff --git a/Detectors/Upgrades/ITS3/reconstruction/src/IOUtils.cxx b/Detectors/Upgrades/ITS3/reconstruction/src/IOUtils.cxx index e898631837169..ff306d83cee44 100644 --- a/Detectors/Upgrades/ITS3/reconstruction/src/IOUtils.cxx +++ b/Detectors/Upgrades/ITS3/reconstruction/src/IOUtils.cxx @@ -17,13 +17,6 @@ #include "ITS3Reconstruction/TopologyDictionary.h" #include "ITSBase/GeometryTGeo.h" #include "ITStracking/TrackingConfigParam.h" -#include "Framework/Logger.h" - -#include - -#include - -#include namespace o2::its3::ioutils { @@ -52,7 +45,7 @@ void convertCompactClusters(gsl::span clusters, const auto detID = c.getSensorID(); auto& cl3d = output.emplace_back(detID, geom->getMatrixT2L(detID) ^ locXYZ); // local --> tracking if (applyMisalignment) { - auto lrID = geom->getLayer(detID); + const auto lrID = geom->getLayer(detID); sigmaY2 += conf.sysErrY2[lrID]; sigmaZ2 += conf.sysErrZ2[lrID]; } @@ -96,13 +89,17 @@ int loadROFrameDataITS3(its::TimeFrame<7>* tf, // Tracking alpha angle // We want that each cluster rotates its tracking frame to the clusters phi // that way the track linearization around the measurement is less biases to the arc - // this means automatically that the measurement on the arc is at 0 - // const float alpha = geom->getSensorRefAlpha(sensorID); - const float radius = std::hypot(gloXYZ.x(), gloXYZ.y()); - const float alpha = std::atan2(gloXYZ.y(), gloXYZ.x()); - - tf->addTrackingFrameInfoToLayer(layer, gloXYZ.x(), gloXYZ.y(), gloXYZ.z(), radius, alpha, - std::array{0, trkXYZ.z()}, + // this means automatically that the measurement on the arc is at 0 for the curved layers + float alpha = geom->getSensorRefAlpha(sensorID); + float x = trkXYZ.x(), y = trkXYZ.y(); + if (layer < 3) { + y = 0.f; + x = std::hypot(gloXYZ.x(), gloXYZ.y()); + alpha = std::atan2(gloXYZ.y(), gloXYZ.x()); + } + + tf->addTrackingFrameInfoToLayer(layer, gloXYZ.x(), gloXYZ.y(), gloXYZ.z(), x, alpha, + std::array{y, trkXYZ.z()}, std::array{sigmaY2, sigmaYZ, sigmaZ2}); /// Rotate to the global frame diff --git a/dependencies/O2Dependencies.cmake b/dependencies/O2Dependencies.cmake index 26f381a4ef79f..c38044f6db935 100644 --- a/dependencies/O2Dependencies.cmake +++ b/dependencies/O2Dependencies.cmake @@ -113,6 +113,24 @@ set_package_properties(fmt PROPERTIES TYPE REQUIRED) find_package(nlohmann_json) set_package_properties(nlohmann_json PROPERTIES TYPE REQUIRED) +# Our Eigen3 install only provides the header files +# 'mock' the cmake target +add_library(Eigen3::Eigen INTERFACE IMPORTED) +set_target_properties(Eigen3::Eigen PROPERTIES + INTERFACE_INCLUDE_DIRECTORIES "${EIGEN3_ROOT}/include/eigen3" +) + +find_package(GBL) +set_package_properties(GBL PROPERTIES TYPE REQUIRED) +if(GBL_FOUND AND NOT TARGET GBL::GBL) + add_library(GBL::GBL INTERFACE IMPORTED) + target_include_directories(GBL::GBL INTERFACE ${GBL_INCLUDE_DIR}) + target_link_libraries(GBL::GBL INTERFACE + ${GBL_LIBRARIES} + Eigen3::Eigen + ) +endif() + find_package(Boost 1.70 COMPONENTS container thread From c7d5036190561d7094388f2049321bbb9394426c Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Thu, 19 Feb 2026 15:27:26 +0100 Subject: [PATCH 5/5] Study: add v0 mc Signed-off-by: Felix Schlepper --- .../study/CMakeLists.txt | 6 + .../include/GlobalTrackingStudy/SVMCStudy.h | 26 + .../study/include/GlobalTrackingStudy/V0Ext.h | 65 ++ .../study/macros/CheckSVMCStudy.C | 604 ++++++++++++ .../study/src/GlobalTrackingStudyLinkDef.h | 3 + .../study/src/SVMCStudy.cxx | 921 ++++++++++++++++++ .../study/src/sv-mc-study-workflow.cxx | 69 ++ 7 files changed, 1694 insertions(+) create mode 100644 Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/SVMCStudy.h create mode 100644 Detectors/GlobalTrackingWorkflow/study/macros/CheckSVMCStudy.C create mode 100644 Detectors/GlobalTrackingWorkflow/study/src/SVMCStudy.cxx create mode 100644 Detectors/GlobalTrackingWorkflow/study/src/sv-mc-study-workflow.cxx diff --git a/Detectors/GlobalTrackingWorkflow/study/CMakeLists.txt b/Detectors/GlobalTrackingWorkflow/study/CMakeLists.txt index df42af503db46..b2be54a4dc167 100644 --- a/Detectors/GlobalTrackingWorkflow/study/CMakeLists.txt +++ b/Detectors/GlobalTrackingWorkflow/study/CMakeLists.txt @@ -16,6 +16,7 @@ o2_add_library(GlobalTrackingStudy SOURCES src/TPCTrackStudy.cxx src/TrackingStudy.cxx src/SVStudy.cxx + src/SVMCStudy.cxx src/TrackMCStudy.cxx src/TPCDataFilter.cxx src/ITSOffsStudy.cxx @@ -50,6 +51,11 @@ o2_add_executable(study-workflow SOURCES src/sv-study-workflow.cxx PUBLIC_LINK_LIBRARIES O2::GlobalTrackingStudy) +o2_add_executable(study-workflow + COMPONENT_NAME sv-mc + SOURCES src/sv-mc-study-workflow.cxx + PUBLIC_LINK_LIBRARIES O2::GlobalTrackingStudy) + o2_add_executable(study-workflow COMPONENT_NAME tpc-track SOURCES src/tpc-track-study-workflow.cxx diff --git a/Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/SVMCStudy.h b/Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/SVMCStudy.h new file mode 100644 index 0000000000000..b7efd27d75012 --- /dev/null +++ b/Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/SVMCStudy.h @@ -0,0 +1,26 @@ +// 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. + +#ifndef O2_SV_GAMMA_STUDY_H +#define O2_SV_GAMMA_STUDY_H + +#include "ReconstructionDataFormats/GlobalTrackID.h" +#include "Framework/DataProcessorSpec.h" +#include "TPCCalibration/CorrectionMapsLoader.h" + +namespace o2::svstudy +{ + +o2::framework::DataProcessorSpec getSVMCStudySpec(o2::dataformats::GlobalTrackID::mask_t srcTracks, o2::dataformats::GlobalTrackID::mask_t srcCls, const o2::tpc::CorrectionMapsLoaderGloOpts& sclOpts); + +} // namespace o2::svstudy + +#endif diff --git a/Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/V0Ext.h b/Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/V0Ext.h index b1a9f6923f04d..c3f3692ead7e2 100644 --- a/Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/V0Ext.h +++ b/Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/V0Ext.h @@ -15,7 +15,9 @@ #define ALICEO2_V0EXT_H #include "ReconstructionDataFormats/V0.h" +#include "ReconstructionDataFormats/PrimaryVertex.h" #include "SimulationDataFormat/MCCompLabel.h" +#include "SimulationDataFormat/MCTrack.h" namespace o2::dataformats { @@ -44,6 +46,69 @@ struct V0Ext { ClassDefNV(V0Ext, 2); }; +// Ordered enum: higher value = got further in the selection. +// Use max-assign semantics so the deepest-reached stage is always recorded. +enum PoolInfo : uint8_t { + Unset = 0, // not yet processed + Unassigned, // seen only in unassigned-track bucket (no collision match) + Considered, // seen in vertex-associated loop + RejTPCMaxX, // processTPCTrack: X position cut + RejTPCCorr, // processTPCTrack: drift correction failed + RejTPCPhotonCls, // processTPCTrack: photon-tune min TPC clusters + RejTPCPhotonDPV, // processTPCTrack: photon-tune z-to-beam cut + RejTPCPhotonRD2, // processTPCTrack: photon-tune transverse radius cut + RejPVContrib, // acceptTrack: PV contributor cut + RejDCAToPV, // acceptTrack: DCA to PV too small + RejShortITS, // buildT2V: short ITS-only track + Accepted, // entered the track pool +}; +struct ProngMCInfoExt { + detectors::DetID::mask_t recoMask; // detector mask + o2::track::TrackPar trk; // most global reco'd track + MCTrack mcTrk; + PoolInfo poolStatus = PoolInfo::Unset; + ClassDefNV(ProngMCInfoExt, 1); +}; + +struct V0MCExt { + V0 v0; // reconstructed v0 + MCTrack mcTrk; + std::array prInfo{}; + const ProngMCInfoExt& getPrInfo(int i) const { return prInfo[i]; } + bool hasBothProngs() const noexcept { return prInfo[0].recoMask.any() && prInfo[1].recoMask.any(); } + enum PairStatus : uint8_t { + Unpaired = 0, // never reached pair loop + NoOverlap, // bracket overlap check failed + MaxPV, // PV contributor cut + // checkV0 rejections, ordered by depth in the function: + RejTgl, // |tgl_pos - tgl_neg| too large (TPC photon tune) + RejD2R, // circle center distance vs radii sum (TPC photon tune) + RejDR, // estimated conversion radius too large (TPC photon tune) + RejDCAFitter, // DCAFitter found no candidate + RejMinR, // V0 radius too close to beam line + RejCausality, // causality check (V0 radius vs track inner radius) + RejPropagation, // propagation to PCA failed + RejPt, // V0 pT too low + RejTgLambda, // V0 tgLambda too large + RejHypothesis, // no mass hypothesis matched + RejDCAToMV, // DCA or cosPAXY to mean vertex failed + RejCosPA, // no PV passed cosPA cut + Found, // reconstructed V0 added to output + }; + PairStatus pairStatus = Unpaired; + + ClassDefNV(V0MCExt, 1); +}; + +struct V0MCFull { + V0 v0; + std::array hypStatus{}; + int pdg = -1; + o2::MCCompLabel label; + + ClassDefNV(V0MCFull, 1); +}; + } // namespace o2::dataformats #endif diff --git a/Detectors/GlobalTrackingWorkflow/study/macros/CheckSVMCStudy.C b/Detectors/GlobalTrackingWorkflow/study/macros/CheckSVMCStudy.C new file mode 100644 index 0000000000000..f4ba800f01cad --- /dev/null +++ b/Detectors/GlobalTrackingWorkflow/study/macros/CheckSVMCStudy.C @@ -0,0 +1,604 @@ +// 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. + +// Analysis macro for sv_debug.root output from SVMCStudy. +// Prints per-species waterfall of signal losses at prong and V0 level, +// and produces diagnostic plots for TPC photon-tune cut optimisation. +// +// Usage: root -l -b -q 'macros/CheckSVMCStudy.C("sv_debug.root")' + +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "GlobalTrackingStudy/V0Ext.h" +#endif + +// ── label arrays (must match enum order in V0Ext.h) ───────────────────────── + +const char* kPoolLabels[] = { + "Unset", + "Unassigned", + "Considered", + "RejTPCMaxX", + "RejTPCCorr", + "RejTPCPhotonCls", + "RejTPCPhotonDPV", + "RejTPCPhotonRD2", + "RejPVContrib", + "RejDCAToPV", + "RejShortITS", + "Accepted", +}; +constexpr int kNPool = sizeof(kPoolLabels) / sizeof(kPoolLabels[0]); + +const char* kPairLabels[] = { + "Unpaired", + "NoOverlap", + "MaxPV", + "RejTgl", + "RejD2R", + "RejDR", + "RejDCAFitter", + "RejMinR", + "RejCausality", + "RejPropagation", + "RejPt", + "RejTgLambda", + "RejHypothesis", + "RejDCAToMV", + "RejCosPA", + "Found", +}; +constexpr int kNPair = sizeof(kPairLabels) / sizeof(kPairLabels[0]); + +// ── helpers ───────────────────────────────────────────────────────────────── + +void printBar(const char* label, int count, int total, int barWidth = 40) +{ + float frac = total > 0 ? (float)count / total : 0.f; + int filled = (int)(frac * barWidth); + printf(" %-18s %6d (%5.1f%%) |", label, count, 100.f * frac); + for (int i = 0; i < barWidth; ++i) { + putchar(i < filled ? '#' : ' '); + } + printf("|\n"); +} + +/// Draw a vertical cut line + label on the current pad +void drawCut(float x, float ymax, const char* label, Color_t col = kRed) +{ + auto* line = new TLine(x, 0, x, ymax); + line->SetLineColor(col); + line->SetLineWidth(2); + line->SetLineStyle(2); + line->Draw("same"); + auto* txt = new TLatex(x, ymax * 0.92, label); + txt->SetTextColor(col); + txt->SetTextSize(0.035); + txt->Draw("same"); +} + +// ── species waterfall ─────────────────────────────────────────────────────── + +void analyseSpecies(TTree* tree, const char* name) +{ + o2::dataformats::V0MCExt* v0{nullptr}; + tree->SetBranchAddress("v0", &v0); + + const int total = tree->GetEntries(); + int poolCount[2][kNPool] = {}; + int pairCount[kNPair] = {}; + int bothAccepted = 0; + // Break down Unpaired by prong pool status + int unpBothAcc = 0, unpOneRej = 0, unpBothRej = 0; + // For Unpaired with both accepted, track why they didn't pair + int unpBothAccUnset = 0; // at least one prong has Unset (not even seen) + + for (int i = 0; i < total; ++i) { + tree->GetEntry(i); + int ps0 = static_cast(v0->prInfo[0].poolStatus); + int ps1 = static_cast(v0->prInfo[1].poolStatus); + if (ps0 >= 0 && ps0 < kNPool) { + poolCount[0][ps0]++; + } + if (ps1 >= 0 && ps1 < kNPool) { + poolCount[1][ps1]++; + } + int pair = static_cast(v0->pairStatus); + if (pair >= 0 && pair < kNPair) { + pairCount[pair]++; + } + bool acc0 = (ps0 == o2::dataformats::PoolInfo::Accepted); + bool acc1 = (ps1 == o2::dataformats::PoolInfo::Accepted); + if (acc0 && acc1) { + bothAccepted++; + } + if (pair == 0) { // Unpaired + if (acc0 && acc1) { + unpBothAcc++; + } else if (!acc0 && !acc1) { + unpBothRej++; + } else { + unpOneRej++; + } + if (ps0 == 0 || ps1 == 0) { // Unset + unpBothAccUnset++; + } + } + } + + printf("\n"); + printf("========================================================\n"); + printf(" %s (both prongs reconstructed: %d)\n", name, total); + printf("========================================================\n"); + + const char* prongName[2] = {"Pos (+)", "Neg (-)"}; + for (int pn = 0; pn < 2; ++pn) { + printf("\n--- Prong %s poolStatus ---\n", prongName[pn]); + for (int j = 0; j < kNPool; ++j) { + if (poolCount[pn][j] > 0) { + printBar(kPoolLabels[j], poolCount[pn][j], total); + } + } + } + + printf("\n--- Both prongs in pool ---\n"); + printBar("BothAccepted", bothAccepted, total); + + printf("\n--- V0 pairStatus (where signal is lost in checkV0) ---\n"); + for (int j = 0; j < kNPair; ++j) { + if (pairCount[j] > 0) { + printBar(kPairLabels[j], pairCount[j], total); + } + } + + if (pairCount[0] > 0) { + int nUnp = pairCount[0]; + printf("\n--- Unpaired breakdown (%d entries) ---\n", nUnp); + printBar("Both rejected", unpBothRej, nUnp); + printBar("One rejected", unpOneRej, nUnp); + printBar("Both accepted", unpBothAcc, nUnp); + printBar("Has Unset prong", unpBothAccUnset, nUnp); + } + + printf("\n--- Cumulative survival (pair-level) ---\n"); + int remaining = total; + printf(" %-18s %6d\n", "Total (both reco)", remaining); + for (int j = 0; j < kNPair - 1; ++j) { + remaining -= pairCount[j]; + if (pairCount[j] > 0 || remaining != total) { + printf(" after %-12s %6d (%5.1f%% of total)\n", kPairLabels[j], remaining, 100.f * remaining / total); + } + } + printf(" %-18s %6d (%5.1f%%)\n", "Found", pairCount[kNPair - 1], 100.f * pairCount[kNPair - 1] / total); +} + +// ── TPC photon-tune diagnostics ───────────────────────────────────────────── + +void analyseTpcPhotonTune(TTree* tree) +{ + // Branches + Char_t dCls{false}, dDPV{false}, dRD2{false}, ok{false}; + float pt{0}, zPVExtrapDiff{0}, cR{0}, drd2{0}; + tree->SetBranchAddress("dCls", &dCls); + tree->SetBranchAddress("pt", &pt); + tree->SetBranchAddress("zPVExtrapDiff", &zPVExtrapDiff); + tree->SetBranchAddress("dDPV", &dDPV); + tree->SetBranchAddress("cR", &cR); + tree->SetBranchAddress("drd2", &drd2); + tree->SetBranchAddress("dRD2", &dRD2); + tree->SetBranchAddress("ok", &ok); + + const int n = tree->GetEntries(); + if (n == 0) { + printf("tpcPhotonTune tree is empty\n"); + return; + } + + // ── First pass: collect values into vectors for binning-free eff/rej ── + std::vector zSig, zBkg, rSig, rBkg; + zSig.reserve(n); + zBkg.reserve(n); + rSig.reserve(n); + rBkg.reserve(n); + + // Determine axis ranges from data (99.5th percentile), filtering NaN/Inf + std::vector allZ, allR; + allZ.reserve(n); + allR.reserve(n); + int nNanR = 0; + + for (int i = 0; i < n; ++i) { + tree->GetEntry(i); + float absZ = std::abs(zPVExtrapDiff); + float r = drd2; + bool rOk = std::isfinite(r); + if (!rOk) { + ++nNanR; + r = -1.f; // placeholder, will be skipped in histo fills + } + allZ.push_back(absZ); + if (rOk) { + allR.push_back(r); + } + if (ok) { + zSig.push_back(absZ); + if (rOk) { + rSig.push_back(r); + } + } else { + zBkg.push_back(absZ); + if (rOk) { + rBkg.push_back(r); + } + } + } + std::sort(allZ.begin(), allZ.end()); + std::sort(allR.begin(), allR.end()); + float zMax = allZ.size() > 0 ? allZ[(int)(allZ.size() * 0.995)] : 100.f; + float rMax = allR.size() > 0 ? allR[(int)(allR.size() * 0.995)] : 100.f; + // round up to a nice number + zMax = std::ceil(zMax / 10.f) * 10.f; + rMax = std::ceil(rMax / 10.f) * 10.f; + if (nNanR > 0) { + printf(" WARNING: %d entries have NaN/Inf in drd2 (cR < rC), skipped for d_r plots\n", nNanR); + } + + printf("\n--- TPC photon-tune diagnostics ---\n"); + printf(" Total tracks: %d (correct PV: %d, wrong PV: %d)\n", n, (int)zSig.size(), (int)zBkg.size()); + printf(" Auto-range: |z| up to %.0f cm, d_r up to %.0f cm\n", zMax, rMax); + + // ── Histograms with data-driven ranges ── + auto* hZSig = new TH1F("hZSig", ";|z_{extrap} #minus z_{PV}| (cm);Normalised", 100, 0, zMax); + auto* hZBkg = new TH1F("hZBkg", ";|z_{extrap} #minus z_{PV}| (cm);Normalised", 100, 0, zMax); + auto* hRSig = new TH1F("hRSig", ";d_{r} = #sqrt{c_{R}^{2} #minus r_{C}^{2}} (cm);Normalised", 100, 0, rMax); + auto* hRBkg = new TH1F("hRBkg", ";d_{r} = #sqrt{c_{R}^{2} #minus r_{C}^{2}} (cm);Normalised", 100, 0, rMax); + auto* hPtSig = new TH1F("hPtSig", ";p_{T} (GeV/#it{c});Normalised", 100, 0, 5); + auto* hPtBkg = new TH1F("hPtBkg", ";p_{T} (GeV/#it{c});Normalised", 100, 0, 5); + auto* h2Sig = new TH2F("h2Sig", "Correct PV match;|z_{extrap} #minus z_{PV}| (cm);d_{r} (cm)", 80, 0, zMax, 80, 0, rMax); + auto* h2Bkg = new TH2F("h2Bkg", "Wrong PV match;|z_{extrap} #minus z_{PV}| (cm);d_{r} (cm)", 80, 0, zMax, 80, 0, rMax); + // zoomed version + float zMaxZoom = 30.f; + float rMaxZoom = std::min(rMax, 120.f); + auto* h2SigZ = new TH2F("h2SigZ", "Correct PV match (zoom);|z_{extrap} #minus z_{PV}| (cm);d_{r} (cm)", 80, 0, zMaxZoom, 80, 0, rMaxZoom); + auto* h2BkgZ = new TH2F("h2BkgZ", "Wrong PV match (zoom);|z_{extrap} #minus z_{PV}| (cm);d_{r} (cm)", 80, 0, zMaxZoom, 80, 0, rMaxZoom); + + // Fill histograms (skip NaN/Inf in drd2) + for (int i = 0; i < n; ++i) { + tree->GetEntry(i); + float absZ = std::abs(zPVExtrapDiff); + bool rOk = std::isfinite(drd2); + if (ok) { + hZSig->Fill(absZ); + hPtSig->Fill(pt); + if (rOk) { + hRSig->Fill(drd2); + h2Sig->Fill(absZ, drd2); + h2SigZ->Fill(absZ, drd2); + } + } else { + hZBkg->Fill(absZ); + hPtBkg->Fill(pt); + if (rOk) { + hRBkg->Fill(drd2); + h2Bkg->Fill(absZ, drd2); + h2BkgZ->Fill(absZ, drd2); + } + } + } + + // Style + gStyle->SetOptStat(0); + auto styleSig = [](TH1* h) { + h->SetLineColor(kBlue); + h->SetLineWidth(2); + h->SetFillColorAlpha(kBlue, 0.15); + }; + auto styleBkg = [](TH1* h) { + h->SetLineColor(kRed); + h->SetLineWidth(2); + h->SetFillColorAlpha(kRed, 0.15); + }; + styleSig(hZSig); + styleBkg(hZBkg); + styleSig(hRSig); + styleBkg(hRBkg); + styleSig(hPtSig); + styleBkg(hPtBkg); + + // Normalise to unit area + if (hZSig->Integral() > 0) { + hZSig->Scale(1. / hZSig->Integral()); + } + if (hZBkg->Integral() > 0) { + hZBkg->Scale(1. / hZBkg->Integral()); + } + if (hRSig->Integral() > 0) { + hRSig->Scale(1. / hRSig->Integral()); + } + if (hRBkg->Integral() > 0) { + hRBkg->Scale(1. / hRBkg->Integral()); + } + if (hPtSig->Integral() > 0) { + hPtSig->Scale(1. / hPtSig->Integral()); + } + if (hPtBkg->Integral() > 0) { + hPtBkg->Scale(1. / hPtBkg->Integral()); + } + + auto makeLeg = [](TH1* hS, TH1* hB) { + auto* leg = new TLegend(0.55, 0.75, 0.88, 0.88); + leg->AddEntry(hS, "Correct PV", "f"); + leg->AddEntry(hB, "Wrong PV", "f"); + leg->SetBorderSize(0); + return leg; + }; + + // Helper to draw a pair with proper axis propagation + auto drawPair = [&](TH1* hS, TH1* hB) { + float ymax = std::max(hS->GetMaximum(), hB->GetMaximum()) * 5; + hB->SetMaximum(ymax); + hB->Draw("hist"); + hS->Draw("hist same"); + makeLeg(hS, hB)->Draw(); + }; + + // ── Canvas 1: 1D distributions ── + auto* c1 = new TCanvas("cTpcPhoton1D", "TPC photon-tune 1D", 1500, 500); + c1->Divide(3, 1); + + c1->cd(1)->SetLogy(); + drawPair(hZSig, hZBkg); + + c1->cd(2)->SetLogy(); + drawPair(hRSig, hRBkg); + + c1->cd(3)->SetLogy(); + drawPair(hPtSig, hPtBkg); + + c1->SaveAs("tpcPhotonTune_1D.png"); + + // ── Canvas 2: 2D signal vs background (full range) ── + auto* c2 = new TCanvas("cTpcPhoton2D", "TPC photon-tune 2D", 1200, 500); + c2->Divide(2, 1); + c2->cd(1); + h2Sig->Draw("colz"); + c2->cd(2); + h2Bkg->Draw("colz"); + c2->SaveAs("tpcPhotonTune_2D.png"); + + // ── Canvas 2b: 2D zoomed to signal region ── + auto* c2b = new TCanvas("cTpcPhoton2Dzoom", "TPC photon-tune 2D zoom", 1200, 500); + c2b->Divide(2, 1); + c2b->cd(1); + h2SigZ->Draw("colz"); + c2b->cd(2); + h2BkgZ->Draw("colz"); + c2b->SaveAs("tpcPhotonTune_2D_zoom.png"); + + // ── Canvas 3: cut scan — signal efficiency & background rejection vs cut ── + // For each variable: "accept if value < cut" + // Signal efficiency = fraction of signal with value < cut + // Background rejection = fraction of background with value >= cut + std::sort(zSig.begin(), zSig.end()); + std::sort(zBkg.begin(), zBkg.end()); + std::sort(rSig.begin(), rSig.end()); + std::sort(rBkg.begin(), rBkg.end()); + + const int nScan = 200; + + auto buildEffRej = [&](const std::vector& sig, const std::vector& bkg, + int nBins, float xMax, + std::vector& effOut, std::vector& rejOut, + std::vector& cutOut) { + float nS = sig.size(), nB = bkg.size(); + effOut.resize(nBins); + rejOut.resize(nBins); + cutOut.resize(nBins); + for (int ib = 0; ib < nBins; ++ib) { + float cut = xMax * (ib + 0.5f) / nBins; + cutOut[ib] = cut; + float sigKept = std::upper_bound(sig.begin(), sig.end(), cut) - sig.begin(); + float bkgKept = std::upper_bound(bkg.begin(), bkg.end(), cut) - bkg.begin(); + effOut[ib] = nS > 0 ? sigKept / nS : 0; + rejOut[ib] = nB > 0 ? 1.f - bkgKept / nB : 0; + } + }; + + std::vector effZ, rejZ, cutZ, effR, rejR, cutR; + buildEffRej(zSig, zBkg, nScan, zMax, effZ, rejZ, cutZ); + buildEffRej(rSig, rBkg, nScan, rMax, effR, rejR, cutR); + + // Build TGraphs for cut-scan plots + auto* gEffZ = new TGraph(nScan, cutZ.data(), effZ.data()); + auto* gRejZ = new TGraph(nScan, cutZ.data(), rejZ.data()); + auto* gEffR = new TGraph(nScan, cutR.data(), effR.data()); + auto* gRejR = new TGraph(nScan, cutR.data(), rejR.data()); + + gEffZ->SetLineColor(kBlue); + gEffZ->SetLineWidth(2); + gRejZ->SetLineColor(kRed); + gRejZ->SetLineWidth(2); + gEffR->SetLineColor(kBlue); + gEffR->SetLineWidth(2); + gRejR->SetLineColor(kRed); + gRejR->SetLineWidth(2); + + // Build ROC curve (signal eff vs bkg eff = 1-rejection) + auto* gROCz = new TGraph(nScan); + auto* gROCr = new TGraph(nScan); + for (int i = 0; i < nScan; ++i) { + gROCz->SetPoint(i, effZ[i], rejZ[i]); + gROCr->SetPoint(i, effR[i], rejR[i]); + } + gROCz->SetLineColor(kBlue); + gROCz->SetLineWidth(2); + gROCz->SetMarkerStyle(1); + gROCr->SetLineColor(kRed); + gROCr->SetLineWidth(2); + gROCr->SetMarkerStyle(1); + + // Helper: find cut value where eff*rej is maximised (best working point) + auto bestWP = [](const std::vector& eff, const std::vector& rej, + const std::vector& cut) -> std::tuple { + float best = -1; + int idx = 0; + for (int i = 0; i < (int)eff.size(); ++i) { + float score = eff[i] * rej[i]; + if (score > best) { + best = score; + idx = i; + } + } + return {cut[idx], eff[idx], rej[idx]}; + }; + auto [wpCutZ, wpEffZ, wpRejZ] = bestWP(effZ, rejZ, cutZ); + auto [wpCutR, wpEffR, wpRejR] = bestWP(effR, rejR, cutR); + + printf(" Best working points (max eff*rej):\n"); + printf(" |z| cut = %.1f cm → eff = %.1f%%, rej = %.1f%%\n", wpCutZ, 100 * wpEffZ, 100 * wpRejZ); + printf(" d_r cut = %.1f cm → eff = %.1f%%, rej = %.1f%%\n", wpCutR, 100 * wpEffR, 100 * wpRejR); + + auto* c3 = new TCanvas("cTpcPhotonEff", "TPC photon-tune cut optimisation", 1500, 900); + c3->Divide(3, 2); + + // Top row: cut scans + c3->cd(1); + auto* frZ = gPad->DrawFrame(0, 0, zMax, 1.05); + frZ->SetTitle("Cut scan: |z_{extrap} #minus z_{PV}|;Cut value (cm);Fraction"); + gEffZ->Draw("L"); + gRejZ->Draw("L"); + drawCut(wpCutZ, 1.0, Form("%.0f cm", wpCutZ), kGreen + 2); + auto* legZ = new TLegend(0.45, 0.15, 0.88, 0.35); + legZ->AddEntry(gEffZ, "Signal eff (keep if < cut)", "l"); + legZ->AddEntry(gRejZ, "Bkg rejection (remove if #geq cut)", "l"); + legZ->SetBorderSize(0); + legZ->Draw(); + + c3->cd(2); + auto* frR = gPad->DrawFrame(0, 0, rMax, 1.05); + frR->SetTitle("Cut scan: d_{r} = #sqrt{c_{R}^{2} #minus r_{C}^{2}};Cut value (cm);Fraction"); + gEffR->Draw("L"); + gRejR->Draw("L"); + drawCut(wpCutR, 1.0, Form("%.0f cm", wpCutR), kGreen + 2); + auto* legR = new TLegend(0.45, 0.15, 0.88, 0.35); + legR->AddEntry(gEffR, "Signal eff (keep if < cut)", "l"); + legR->AddEntry(gRejR, "Bkg rejection (remove if #geq cut)", "l"); + legR->SetBorderSize(0); + legR->Draw(); + + // Top-right: ROC + c3->cd(3); + auto* frROC = gPad->DrawFrame(0, 0, 1.05, 1.05); + frROC->SetTitle("ROC: signal efficiency vs bkg rejection;Signal efficiency;Background rejection"); + gROCz->Draw("L"); + gROCr->Draw("L"); + // mark working points + auto* mZ = new TMarker(wpEffZ, wpRejZ, 20); + mZ->SetMarkerColor(kBlue); + mZ->SetMarkerSize(1.5); + mZ->Draw(); + auto* mR = new TMarker(wpEffR, wpRejR, 20); + mR->SetMarkerColor(kRed); + mR->SetMarkerSize(1.5); + mR->Draw(); + auto* legROC = new TLegend(0.15, 0.15, 0.55, 0.35); + legROC->AddEntry(gROCz, "|z| cut", "l"); + legROC->AddEntry(gROCr, "d_{r} cut", "l"); + legROC->SetBorderSize(0); + legROC->Draw(); + + // Bottom row: zoomed cut scans (near signal region) + float zZoom = std::min(zMax, 50.f); + float rZoom = std::min(rMax, 50.f); + + c3->cd(4); + auto* frZz = gPad->DrawFrame(0, 0, zZoom, 1.05); + frZz->SetTitle("Cut scan zoom: |z_{extrap} #minus z_{PV}|;Cut value (cm);Fraction"); + gEffZ->Draw("L"); + gRejZ->Draw("L"); + drawCut(wpCutZ, 1.0, Form("%.0f cm", wpCutZ), kGreen + 2); + + c3->cd(5); + auto* frRz = gPad->DrawFrame(0, 0, rZoom, 1.05); + frRz->SetTitle("Cut scan zoom: d_{r};Cut value (cm);Fraction"); + gEffR->Draw("L"); + gRejR->Draw("L"); + drawCut(wpCutR, 1.0, Form("%.0f cm", wpCutR), kGreen + 2); + + // Bottom-right: significance-like metric (eff * rej) vs cut + auto* gProdZ = new TGraph(nScan); + auto* gProdR = new TGraph(nScan); + for (int i = 0; i < nScan; ++i) { + gProdZ->SetPoint(i, cutZ[i], effZ[i] * rejZ[i]); + gProdR->SetPoint(i, cutR[i], effR[i] * rejR[i]); + } + gProdZ->SetLineColor(kBlue); + gProdZ->SetLineWidth(2); + gProdR->SetLineColor(kRed); + gProdR->SetLineWidth(2); + c3->cd(6); + auto* frP = gPad->DrawFrame(0, 0, std::max(zMax, rMax), 1.05); + frP->SetTitle("eff #times rej vs cut (higher = better);Cut value (cm);eff #times rej"); + gProdZ->Draw("L"); + gProdR->Draw("L"); + auto* legP = new TLegend(0.55, 0.75, 0.88, 0.88); + legP->AddEntry(gProdZ, "|z| cut", "l"); + legP->AddEntry(gProdR, "d_{r} cut", "l"); + legP->SetBorderSize(0); + legP->Draw(); + + c3->SaveAs("tpcPhotonTune_effrej.png"); +} + +// ── main ──────────────────────────────────────────────────────────────────── + +void CheckSVMCStudy(const char* fname = "sv_debug.root") +{ + auto* f = TFile::Open(fname); + if (!f || f->IsZombie()) { + printf("Cannot open %s\n", fname); + return; + } + + printf("NOTE: Trees only contain V0s where both prongs were reconstructed.\n"); + printf(" V0s with missing prongs are not included (check log output).\n"); + + // const char* species[] = {"Photon", "K0"}; + const char* species[] = {"Photon"}; + for (auto* name : species) { + auto* tree = dynamic_cast(f->Get(name)); + if (!tree) { + printf("Tree '%s' not found, skipping\n", name); + continue; + } + analyseSpecies(tree, name); + } + + // TPC photon-tune cut diagnostics + auto* tpcTree = dynamic_cast(f->Get("tpcPhotonTune")); + if (tpcTree) { + analyseTpcPhotonTune(tpcTree); + } else { + printf("Tree 'tpcPhotonTune' not found, skipping diagnostics\n"); + } + + f->Close(); +} diff --git a/Detectors/GlobalTrackingWorkflow/study/src/GlobalTrackingStudyLinkDef.h b/Detectors/GlobalTrackingWorkflow/study/src/GlobalTrackingStudyLinkDef.h index 416820fc9aebb..42abcce781239 100644 --- a/Detectors/GlobalTrackingWorkflow/study/src/GlobalTrackingStudyLinkDef.h +++ b/Detectors/GlobalTrackingWorkflow/study/src/GlobalTrackingStudyLinkDef.h @@ -16,7 +16,10 @@ #pragma link off all functions; #pragma link C++ class o2::dataformats::ProngInfoExt + ; +#pragma link C++ class o2::dataformats::ProngMCInfoExt + ; #pragma link C++ class o2::dataformats::V0Ext + ; +#pragma link C++ class o2::dataformats::V0MCExt + ; +#pragma link C++ class o2::dataformats::V0MCFull + ; #pragma link C++ class o2::dataformats::TrackInfoExt + ; #pragma link C++ class std::vector < o2::dataformats::TrackInfoExt> + ; #pragma link C++ class std::vector < o2::dataformats::ProngInfoExt> + ; diff --git a/Detectors/GlobalTrackingWorkflow/study/src/SVMCStudy.cxx b/Detectors/GlobalTrackingWorkflow/study/src/SVMCStudy.cxx new file mode 100644 index 0000000000000..624aa5ac3d4ff --- /dev/null +++ b/Detectors/GlobalTrackingWorkflow/study/src/SVMCStudy.cxx @@ -0,0 +1,921 @@ +// 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. + +#include +#include +#include + +#include + +#include "GlobalTrackingStudy/SVMCStudy.h" +#include "CommonConstants/GeomConstants.h" +#include "GlobalTrackingStudy/V0Ext.h" +#include "ReconstructionDataFormats/V0.h" +#include "DataFormatsGlobalTracking/RecoContainer.h" +#include "DataFormatsGlobalTracking/RecoContainerCreateTracksVariadic.h" +#include "ReconstructionDataFormats/GlobalTrackID.h" +#include "DetectorsCommonDataFormats/DetID.h" +#include "DetectorsBase/GRPGeomHelper.h" +#include "ReconstructionDataFormats/PrimaryVertex.h" +#include "ReconstructionDataFormats/VtxTrackRef.h" +#include "Framework/CCDBParamSpec.h" +#include "Framework/Task.h" +#include "DetectorsVertexing/SVertexerParams.h" +#include "DetectorsVertexing/SVertexer.h" +#include "DetectorsVertexing/SVertexHypothesis.h" +#include "CommonUtils/TreeStreamRedirector.h" +#include "SimulationDataFormat/MCEventLabel.h" +#include "Steer/MCKinematicsReader.h" +#include "SimulationDataFormat/MCUtils.h" +#include "TPCCalibration/VDriftHelper.h" +#include "TPCCalibration/CorrectionMapsLoader.h" +#include "DataFormatsITSMFT/TrkClusRef.h" +#include "DCAFitter/DCAFitterN.h" + +namespace o2::svstudy +{ + +using DetID = o2::detectors::DetID; +using DataRequest = o2::globaltracking::DataRequest; + +using PVertex = o2::dataformats::PrimaryVertex; +using V2TRef = o2::dataformats::VtxTrackRef; +using VTIndex = o2::dataformats::VtxTrackIndex; +using GTrackID = o2::dataformats::VtxTrackIndex; +using TBracket = o2::math_utils::Bracketf_t; +using V0ID = o2::dataformats::V0Index; +using TrackCand = o2::vertexing::SVertexer::TrackCand; +using HypV0 = o2::vertexing::SVertexer::HypV0; + +class SVMCStudySpec final : public o2::framework::Task +{ + public: + SVMCStudySpec(std::shared_ptr dr, std::shared_ptr gr, const o2::tpc::CorrectionMapsLoaderGloOpts& sclOpts, GTrackID::mask_t src) : mDataRequest(dr), mGGCCDBRequest(gr), mSrc(src) + { + mTPCCorrMapsLoader.setLumiScaleType(sclOpts.lumiType); + mTPCCorrMapsLoader.setLumiScaleMode(sclOpts.lumiMode); + mTPCCorrMapsLoader.setCheckCTPIDCConsistency(sclOpts.checkCTPIDCconsistency); + } + void init(o2::framework::InitContext& ic) final; + void run(o2::framework::ProcessingContext& pc) final; + void endOfStream(o2::framework::EndOfStreamContext& ec) final; + void finaliseCCDB(o2::framework::ConcreteDataMatcher& matcher, void* obj) final; + + private: + void updateTimeDependentParams(o2::framework::ProcessingContext& pc); + void buildV0MC(); + void buildT2V0(); + void buildT2V(); + bool acceptTrack(const GTrackID gid, const o2::track::TrackParCov& trc, dataformats::ProngMCInfoExt* prInfo = nullptr) const; + bool processTPCTrack(const o2::tpc::TrackTPC& trTPC, GTrackID gid, int vtxid, dataformats::ProngMCInfoExt* prInfo); + float correctTPCTrack(TrackCand& trc, const o2::tpc::TrackTPC& tTPC, float tmus, float tmusErr) const; + bool checkV0(const TrackCand& seedP, const TrackCand& seedN); + + std::array, 2> mTracksPool{}; // pools of positive and negative seeds sorted in min VtxID + std::array, 2> mVtxFirstTrack{}; // 1st pos. and neg. track of the pools for each vertex + const vertexing::SVertexerParams* mSVParams{nullptr}; + + std::shared_ptr mDataRequest; + std::shared_ptr mGGCCDBRequest; + float mMUS2TPCBin{1.f / (8 * o2::constants::lhc::LHCBunchSpacingMUS)}; + float mTPCBin2Z{0}; + float mTPCVDrift{0}; + float mTPCVDriftCorrFact{1.}; ///< TPC nominal correction factort (wrt ref) + float mTPCVDriftRef{0}; + std::span mTPCTracksArray; ///< input TPC tracks span + const o2::tpc::ClusterNativeAccess* mTPCClusterIdxStruct{nullptr}; ///< struct holding the TPC cluster indices + std::span mTPCTrackClusIdx; ///< input TPC track cluster indices span + std::span mPVertices; + o2::tpc::VDriftHelper mTPCVDriftHelper; + o2::tpc::CorrectionMapsLoader mTPCCorrMapsLoader; + o2::globaltracking::RecoContainer mRecoData; + float mBz{-999}; + const o2::dataformats::MeanVertexObject* mMeanVertex{nullptr}; + std::vector mV0s; + std::unordered_map mV0MC; + std::unordered_map mDau2V0; + GTrackID::mask_t mSrc; + vertexing::DCAFitter2 mFitterV0; + std::array mV0Hyps; + TStopwatch mTimer; + + float mMinR2ToMeanVertex{-999.f}; + float mMaxDCAXY2ToMeanVertex{-999.f}; + float mMinPt2V0{-999.f}; + float mMaxTgl2V0{-999.f}; + + const std::map mPDGNames{{22, "Photon"}, {310, "K0"}}; + const std::set mPDGV0Acc{22, 310}; + std::unique_ptr mDBGTree{nullptr}; +}; + +void SVMCStudySpec::init(o2::framework::InitContext& ic) +{ + mTimer.Stop(); + mTimer.Reset(); + o2::base::GRPGeomHelper::instance().setRequest(mGGCCDBRequest); + if (mSrc[GTrackID::TPC]) { + mTPCCorrMapsLoader.init(ic); + } + mDBGTree = std::make_unique("sv_debug.root"); +} + +void SVMCStudySpec::run(o2::framework::ProcessingContext& pc) +{ + mTimer.Start(false); + mRecoData.collectData(pc, *mDataRequest); + updateTimeDependentParams(pc); + mPVertices = mRecoData.getPrimaryVertices(); + buildV0MC(); + buildT2V0(); + buildT2V(); + + // pair loop + int ntrP = mTracksPool[vertexing::SVertexer::POS].size(), ntrN = mTracksPool[vertexing::SVertexer::NEG].size(); + for (int itp = 0; itp < ntrP; itp++) { + auto& seedP = mTracksPool[vertexing::SVertexer::POS][itp]; + const int firstN = mVtxFirstTrack[vertexing::SVertexer::NEG][seedP.vBracket.getMin()]; + if (firstN < 0) { + continue; + } + const auto lblP = mRecoData.getTrackMCLabel(seedP.gid); + const auto dauPM = mDau2V0.find(lblP); + for (int itn = firstN; itn < ntrN; itn++) { // start from the 1st negative track of lowest-ID vertex of positive + auto& seedN = mTracksPool[vertexing::SVertexer::NEG][itn]; + dataformats::V0MCExt* moth{nullptr}; + const auto lblN = mRecoData.getTrackMCLabel(seedN.gid); + const auto dauNM = mDau2V0.find(lblN); + if (dauNM == dauPM && dauNM != mDau2V0.end()) { + moth = &mV0MC[mDau2V0[lblN]]; + } + + if (seedN.vBracket > seedP.vBracket) { // all vertices compatible with seedN are in future wrt that of seedP + LOG(debug) << "Brackets do not match"; + if (moth && dataformats::V0MCExt::NoOverlap > moth->pairStatus) { + moth->pairStatus = dataformats::V0MCExt::NoOverlap; + } + break; + } + if (mSVParams->maxPVContributors < 2 && seedP.gid.isPVContributor() + seedN.gid.isPVContributor() > mSVParams->maxPVContributors) { + if (moth && dataformats::V0MCExt::MaxPV > moth->pairStatus) { + moth->pairStatus = dataformats::V0MCExt::MaxPV; + } + continue; + } + checkV0(seedP, seedN); + } + } + std::map good; + std::map> dup; + for (const auto& v0 : mV0s) { + if (v0.pdg >= 0) { + ++good[v0.pdg]; + if (dup.find(v0.pdg) == dup.end()) { + dup[v0.pdg] = std::set(); + } + dup[v0.pdg].insert(v0.label); + } + } + LOGP(info, "Found {} v0s", mV0s.size()); + for (const auto [pdg, c] : good) { + LOGP(info, "\t{} has {} good ones with {} unique ones", mPDGNames.at(pdg), c, dup[pdg].size()); + } + + LOGP(info, "Creating output"); + for (const auto& [_, v0] : mV0MC) { + if (v0.hasBothProngs()) { // only interested in V0 where both prongs were reconstructed + (*mDBGTree) << mPDGNames.at(v0.mcTrk.GetPdgCode()) + << "v0=" << v0 + << "\n"; + } + } + for (const auto& v0 : mV0s) { + (*mDBGTree) << "out" + << "v0=" << v0 + << "\n"; + } + + mTimer.Stop(); +} + +void SVMCStudySpec::endOfStream(o2::framework::EndOfStreamContext& /*ec*/) +{ + mDBGTree->Close(); + mDBGTree.reset(); + LOGP(info, "SVMCStudy total timing: Cpu: {:.3f} Real: {:.3f} s", mTimer.CpuTime(), mTimer.RealTime()); +}; + +void SVMCStudySpec::finaliseCCDB(o2::framework::ConcreteDataMatcher& matcher, void* obj) +{ + if (o2::base::GRPGeomHelper::instance().finaliseCCDB(matcher, obj)) { + return; + } + if (mTPCVDriftHelper.accountCCDBInputs(matcher, obj)) { + return; + } + if (mTPCCorrMapsLoader.accountCCDBInputs(matcher, obj)) { + return; + } + if (matcher == o2::framework::ConcreteDataMatcher("GLO", "MEANVERTEX", 0)) { + LOG(info) << "Imposing new MeanVertex: " << ((const o2::dataformats::MeanVertexObject*)obj)->asString(); + mMeanVertex = (const o2::dataformats::MeanVertexObject*)obj; + return; + } + if (matcher == o2::framework::ConcreteDataMatcher("GLO", "SVPARAM", 0)) { + LOG(info) << "SVertexer Params updated from ccdb"; + return; + } +}; + +void SVMCStudySpec::updateTimeDependentParams(o2::framework::ProcessingContext& pc) +{ + o2::base::GRPGeomHelper::instance().checkUpdates(pc); + if (mSrc[GTrackID::TPC]) { + mTPCVDriftHelper.extractCCDBInputs(pc); + mTPCCorrMapsLoader.extractCCDBInputs(pc); + } + if (static bool initOnceDone{false}; !initOnceDone) { // these params need to be queried only once + pc.inputs().get("SVParam"); + o2::vertexing::SVertexerParams::Instance().printKeyValues(true, true); + mSVParams = &vertexing::SVertexerParams::Instance(); + } + if (mSrc[GTrackID::TPC]) { + bool updateMaps = false; + if (mTPCCorrMapsLoader.isUpdated()) { + mTPCCorrMapsLoader.acknowledgeUpdate(); + updateMaps = true; + } + if (mTPCVDriftHelper.isUpdated()) { + LOGP(info, "Updating TPC fast transform map with new VDrift factor of {} wrt reference {} and DriftTimeOffset correction {} wrt {} from source {}", + mTPCVDriftHelper.getVDriftObject().corrFact, mTPCVDriftHelper.getVDriftObject().refVDrift, + mTPCVDriftHelper.getVDriftObject().timeOffsetCorr, mTPCVDriftHelper.getVDriftObject().refTimeOffset, + mTPCVDriftHelper.getSourceName()); + + const auto& v = mTPCVDriftHelper.getVDriftObject(); + mTPCVDrift = v.refVDrift * v.corrFact; + mTPCVDriftCorrFact = v.corrFact; + mTPCVDriftRef = v.refVDrift; + mTPCBin2Z = mTPCVDrift / mMUS2TPCBin; + + mTPCVDriftHelper.acknowledgeUpdate(); + updateMaps = true; + } + if (updateMaps) { + mTPCCorrMapsLoader.updateVDrift(mTPCVDriftHelper.getVDriftObject().corrFact, mTPCVDriftHelper.getVDriftObject().refVDrift, mTPCVDriftHelper.getVDriftObject().getTimeOffset()); + } + + // Fitter setup + mBz = o2::base::Propagator::Instance()->getNominalBz(); + mFitterV0.setFitterID(0); + mFitterV0.setBz(mBz); + mFitterV0.setUseAbsDCA(mSVParams->useAbsDCA); + mFitterV0.setPropagateToPCA(false); + mFitterV0.setMaxR(mSVParams->maxRIni); + mFitterV0.setMinParamChange(mSVParams->minParamChange); + mFitterV0.setMinRelChi2Change(mSVParams->minRelChi2Change); + mFitterV0.setMaxDZIni(mSVParams->maxDZIni); + mFitterV0.setMaxDXYIni(mSVParams->maxDXYIni); + mFitterV0.setMaxChi2(mSVParams->maxChi2); + mFitterV0.setMatCorrType(o2::base::Propagator::MatCorrType(mSVParams->matCorr)); + mFitterV0.setUsePropagator(mSVParams->usePropagator); + mFitterV0.setRefitWithMatCorr(mSVParams->refitWithMatCorr); + mFitterV0.setMaxStep(mSVParams->maxStep); + mFitterV0.setMaxSnp(mSVParams->maxSnp); + mFitterV0.setMinXSeed(mSVParams->minXSeed); + mFitterV0.setBz(mBz); + + // Hypothesis + using PID = o2::track::PID; + mV0Hyps[HypV0::Photon].set(PID::Photon, PID::Electron, PID::Electron, mSVParams->pidCutsPhoton, mBz); + mV0Hyps[HypV0::K0].set(PID::K0, PID::Pion, PID::Pion, mSVParams->pidCutsK0, mBz); + mV0Hyps[HypV0::Lambda].set(PID::Lambda, PID::Proton, PID::Pion, mSVParams->pidCutsLambda, mBz); + mV0Hyps[HypV0::AntiLambda].set(PID::Lambda, PID::Pion, PID::Proton, mSVParams->pidCutsLambda, mBz); + + // precalculated selection cuts + mMinR2ToMeanVertex = mSVParams->minRToMeanVertex * mSVParams->minRToMeanVertex; + mMaxDCAXY2ToMeanVertex = mSVParams->maxDCAXYToMeanVertex * mSVParams->maxDCAXYToMeanVertex; + mMinPt2V0 = mSVParams->minPtV0 * mSVParams->minPtV0; + mMaxTgl2V0 = mSVParams->maxTglV0 * mSVParams->maxTglV0; + } + pc.inputs().get("meanvtx"); +} + +void SVMCStudySpec::buildV0MC() +{ + LOGP(info, "Building mc particle association"); + std::map mCounts; + steer::MCKinematicsReader reader("collisioncontext.root"); + for (int iSrc{0}; iSrc < reader.getNSources(); ++iSrc) { + for (int iEve{0}; iEve < reader.getNEvents(iSrc); ++iEve) { + const auto& mcTrks = reader.getTracks(iSrc, iEve); + const int nTrks = (int)mcTrks.size(); + for (int iTrk{0}; iTrk < nTrks; ++iTrk) { + const auto trk = reader.getTrack(iSrc, iEve, iTrk); + if (mPDGV0Acc.contains(trk->GetPdgCode())) { + if ((trk->getLastDaughterTrackId() - trk->getFirstDaughterTrackId()) != 1) { + continue; + } + MCTrack const *dau0{nullptr}, *dau1{nullptr}; + if (!(dau0 = mcutils::MCTrackNavigator::getDaughter0(*trk, mcTrks)) || + !(dau1 = mcutils::MCTrackNavigator::getDaughter1(*trk, mcTrks))) { + continue; + } + dataformats::V0MCExt v0; + v0.mcTrk = *trk; + v0.prInfo[0].mcTrk = *dau0; + v0.prInfo[1].mcTrk = *dau1; + o2::MCCompLabel lblM{iTrk, iEve, iSrc}, lblDau0{trk->getFirstDaughterTrackId(), iEve, iSrc}, lblDau1{trk->getLastDaughterTrackId(), iEve, iSrc}; + mV0MC[lblM] = v0; + mDau2V0[lblDau0] = lblM; + mDau2V0[lblDau1] = lblM; + ++mCounts[trk->GetPdgCode()]; + } + } + } + } + LOGP(info, "Found {} MC V0 total candidates", mV0MC.size()); + for (const auto [pdg, c] : mCounts) { + LOGP(info, "\t{} has {}", mPDGNames.at(pdg), c); + } +} + +void SVMCStudySpec::buildT2V0() +{ + LOGP(info, "Building track to v0 association"); + // find the most global true representation of the reconstructed prong + std::map mCounts, mCountsBoth; + int totalTracks{0}; + auto creator = [&](auto& trk, GTrackID gid, float /*time0*/, float /*terr*/) { + if constexpr (isBarrelTrack()) { + const auto lbl = mRecoData.getTrackMCLabel(gid); + if (mDau2V0.find(lbl) != mDau2V0.end()) { + auto& v0 = mV0MC[mDau2V0[lbl]]; + int posneg = trk.getSign() < 0 ? 1 : 0; + v0.prInfo[posneg].recoMask = gid.getSourceDetectorsMask(); + v0.prInfo[posneg].trk = trk; + if (v0.prInfo[((posneg) ? 0 : 1)].recoMask.any()) { + ++mCountsBoth[v0.mcTrk.GetPdgCode()]; + } + ++mCounts[v0.mcTrk.GetPdgCode()]; + ++totalTracks; + + return true; + } + } + return false; + }; + mRecoData.createTracksVariadic(creator); + LOGP(info, "Selected {} reconstructed tracks", totalTracks); + for (const auto [pdg, c] : mCounts) { + LOGP(info, "\t{} has {} prongs (single leg)", mPDGNames.at(pdg), c); + } + LOGP(info, "Findable v0 where both prongs were reconstructed", totalTracks); + for (const auto [pdg, c] : mCountsBoth) { + LOGP(info, "\t{} has {}", mPDGNames.at(pdg), c); + } +} + +void SVMCStudySpec::buildT2V() +{ + LOGP(info, "Building T2V"); + auto trackIndex = mRecoData.getPrimaryVertexMatchedTracks(); // Global ID's for associated tracks + auto vtxRefs = mRecoData.getPrimaryVertexMatchedTrackRefs(); // references from vertex to these track IDs + bool isTPCloaded = mRecoData.isTrackSourceLoaded(GTrackID::TPC); + bool isITSloaded = mRecoData.isTrackSourceLoaded(GTrackID::ITS); + bool isITSTPCloaded = mRecoData.isTrackSourceLoaded(GTrackID::ITSTPC); + if (isTPCloaded && !mSVParams->mExcludeTPCtracks) { + mTPCTracksArray = mRecoData.getTPCTracks(); + mTPCTrackClusIdx = mRecoData.getTPCTracksClusterRefs(); + mTPCClusterIdxStruct = &mRecoData.inputsTPCclusters->clusterIndex; + } + std::unordered_map> tmap; + std::unordered_map rejmap; + size_t nv = vtxRefs.size() - 1; // The last entry is for unassigned tracks, ignore them + for (int i = 0; i < 2; i++) { + mTracksPool[i].clear(); + mVtxFirstTrack[i].clear(); + mVtxFirstTrack[i].resize(nv, -1); + } + for (int iv = 0; iv < nv; iv++) { + const auto& vtref = vtxRefs[iv]; + int it = vtref.getFirstEntry(), itLim = it + vtref.getEntries(); + for (; it < itLim; it++) { + auto tvid = trackIndex[it]; + if (!mRecoData.isTrackSourceLoaded(tvid.getSource())) { + continue; + } + + const auto& trc = mRecoData.getTrackParam(tvid); + int posneg = trc.getSign() < 0 ? 1 : 0; + o2::MCCompLabel lbl = mRecoData.getTrackMCLabel(tvid); + dataformats::ProngMCInfoExt* prInfo{nullptr}; + if (mDau2V0.find(lbl) != mDau2V0.end()) { + prInfo = &mV0MC[mDau2V0[lbl]].prInfo[posneg]; + if (dataformats::PoolInfo::Considered > prInfo->poolStatus) { + prInfo->poolStatus = dataformats::PoolInfo::Considered; + } + } + + if (tvid.getSource() == GTrackID::TPC) { + if (mSVParams->mExcludeTPCtracks) { + continue; + } + if (processTPCTrack(mTPCTracksArray[tvid], tvid, iv, prInfo)) { + continue; + } + } + + if (tvid.isAmbiguous()) { // was this track already processed? + auto tref = tmap.find(tvid); + if (tref != tmap.end()) { + mTracksPool[tref->second.second][tref->second.first].vBracket.setMax(iv); // this track was already processed with other vertex, account the latter + continue; + } + // was it already rejected? + if (rejmap.find(tvid) != rejmap.end()) { + continue; + } + } + + bool hasTPC = false; + bool compatibleWithProton = false; + auto tpcGID = mRecoData.getTPCContributorGID(tvid); + if (tpcGID.isIndexSet() && isTPCloaded) { + hasTPC = true; + } + + // get Nclusters in the ITS if available + int nITSclu = -1; + bool shortOBITSOnlyTrack = false; + auto itsGID = mRecoData.getITSContributorGID(tvid); + if (itsGID.getSource() == GTrackID::ITS) { + if (isITSloaded) { + auto& itsTrack = mRecoData.getITSTrack(itsGID); + nITSclu = itsTrack.getNumberOfClusters(); + if (itsTrack.hasHitOnLayer(6) && itsTrack.hasHitOnLayer(5) && itsTrack.hasHitOnLayer(4) && itsTrack.hasHitOnLayer(3)) { + shortOBITSOnlyTrack = true; + } + } + } else if (itsGID.getSource() == GTrackID::ITSAB) { + if (isITSTPCloaded) { + auto& itsABTracklet = mRecoData.getITSABRef(itsGID); + nITSclu = itsABTracklet.getNClusters(); + } + } + if (!acceptTrack(tvid, trc, prInfo)) { + if (tvid.isAmbiguous()) { + rejmap[tvid] = true; + } + continue; + } + if ((isTPCloaded && !hasTPC) && (isITSloaded && (nITSclu < mSVParams->mITSSAminNclu && (!shortOBITSOnlyTrack || mSVParams->mRejectITSonlyOBtrack)))) { + if (prInfo && dataformats::PoolInfo::RejShortITS > prInfo->poolStatus) { + prInfo->poolStatus = dataformats::PoolInfo::RejShortITS; + } + continue; // reject short ITS-only + } + + float r = std::hypot(trc.getX(), trc.getY()); + mTracksPool[posneg].emplace_back(TrackCand{trc, tvid, {iv, iv}, r, hasTPC, (int8_t)nITSclu, compatibleWithProton}); + if (prInfo && dataformats::PoolInfo::Accepted > prInfo->poolStatus) { + prInfo->poolStatus = dataformats::PoolInfo::Accepted; + } + if (tvid.getSource() == GTrackID::TPC) { // constrained TPC track? + correctTPCTrack(mTracksPool[posneg].back(), mTPCTracksArray[tvid], -1, -1); + } + if (tvid.isAmbiguous()) { // track attached to >1 vertex, remember that it was already processed + tmap[tvid] = {mTracksPool[posneg].size() - 1, posneg}; + } + } + } + + // register 1st track of each charge for each vertex + for (int pn = 0; pn < 2; pn++) { + auto& vtxFirstT = mVtxFirstTrack[pn]; + const auto& tracksPool = mTracksPool[pn]; + for (unsigned i = 0; i < tracksPool.size(); i++) { + const auto& t = tracksPool[i]; + for (int j{t.vBracket.getMin()}; j <= t.vBracket.getMax(); ++j) { + if (vtxFirstT[j] == -1) { + vtxFirstT[j] = i; + } + } + } + } + + { + // check the unassigned tracks + std::map mCounts; + const auto& vtref = vtxRefs[nv]; + int it = vtref.getFirstEntry(), itLim = it + vtref.getEntries(); + for (; it < itLim; it++) { + auto tvid = trackIndex[it]; + if (!mRecoData.isTrackSourceLoaded(tvid.getSource())) { + continue; + } + const auto& trc = mRecoData.getTrackParam(tvid); + int posneg = trc.getSign() < 0 ? 1 : 0; + o2::MCCompLabel lbl = mRecoData.getTrackMCLabel(tvid); + if (mDau2V0.find(lbl) == mDau2V0.end()) { + continue; + } + auto& prInfo = mV0MC[mDau2V0[lbl]].prInfo[posneg]; + if (dataformats::PoolInfo::Unassigned > prInfo.poolStatus) { + prInfo.poolStatus = dataformats::PoolInfo::Unassigned; + } + ++mCounts[mV0MC[mDau2V0[lbl]].mcTrk.GetPdgCode()]; + } + LOGP(info, "Unassigned to any collision:"); + for (const auto [pdg, c] : mCounts) { + LOGP(info, "\t{} has {}", mPDGNames.at(pdg), c); + } + } + + LOG(info) << "Collected " << mTracksPool[vertexing::SVertexer::POS].size() << " positive and " << mTracksPool[vertexing::SVertexer::NEG].size() << " negative seeds"; +} + +bool SVMCStudySpec::acceptTrack(const GTrackID gid, const o2::track::TrackParCov& trc, dataformats::ProngMCInfoExt* prInfo) const +{ + auto updatePrInfo = [&](dataformats::PoolInfo s) { + if (prInfo && s > prInfo->poolStatus) { + prInfo->poolStatus = s; + } + }; + + if (gid.isPVContributor() && mSVParams->maxPVContributors < 1) { + updatePrInfo(dataformats::PoolInfo::RejPVContrib); + return false; + } + + // DCA to mean vertex + if (mSVParams->minDCAToPV > 0.f) { + o2::track::TrackPar trp(trc); + std::array dca{}; + auto* prop = o2::base::Propagator::Instance(); + if (mSVParams->usePropagator) { + if (trp.getX() > mSVParams->minRFor3DField && !prop->PropagateToXBxByBz(trp, mSVParams->minRFor3DField, mSVParams->maxSnp, mSVParams->maxStep, o2::base::Propagator::MatCorrType(mSVParams->matCorr))) { + return true; // we don't need actually to propagate to the beam-line + } + if (!prop->propagateToDCA(mMeanVertex->getXYZ(), trp, prop->getNominalBz(), mSVParams->maxStep, o2::base::Propagator::MatCorrType(mSVParams->matCorr), &dca)) { + return true; + } + } else { + if (!trp.propagateParamToDCA(mMeanVertex->getXYZ(), prop->getNominalBz(), &dca)) { + return true; + } + } + if (std::abs(dca[0]) < mSVParams->minDCAToPV) { + updatePrInfo(dataformats::PoolInfo::RejDCAToPV); + return false; + } + } + return true; +} + +bool SVMCStudySpec::processTPCTrack(const o2::tpc::TrackTPC& trTPC, GTrackID gid, int vtxid, dataformats::ProngMCInfoExt* prInfo) +{ + int posneg = trTPC.getSign() < 0 ? 1 : 0; + auto updatePrInfo = [&](dataformats::PoolInfo s) { + if (prInfo && s > prInfo->poolStatus) { + prInfo->poolStatus = s; + } + }; + + if (mSVParams->mTPCTrackMaxX > 0. && trTPC.getX() > mSVParams->mTPCTrackMaxX) { + updatePrInfo(dataformats::PoolInfo::RejTPCMaxX); + return true; + } + // if TPC track is unconstrained, try to create in the tracks pool a clone constrained to vtxid vertex time. + if (trTPC.hasBothSidesClusters()) { // this is effectively constrained track + return false; // let it be processed as such + } + const auto& vtx = mPVertices[vtxid]; + auto twe = vtx.getTimeStamp(); + + bool compatibleWithProton = false; + auto& trLoc = mTracksPool[posneg].emplace_back(TrackCand{trTPC, gid, {vtxid, vtxid}, 0., true, -1, compatibleWithProton}); + auto err = correctTPCTrack(trLoc, trTPC, twe.getTimeStamp(), twe.getTimeStampError()); + if (err < 0) { + mTracksPool[posneg].pop_back(); // discard + updatePrInfo(dataformats::PoolInfo::RejTPCCorr); + return true; + } + + if (mSVParams->mTPCTrackPhotonTune) { + // require minimum of tpc clusters + bool dCls = trTPC.getNClusters() < mSVParams->mTPCTrackMinNClusters; + // check track z cuts + float zPVExtrapDiff = (trLoc.getX() * trLoc.getTgl()) - trLoc.getZ() + vtx.getZ(); + bool dDPV = std::abs(zPVExtrapDiff) > mSVParams->mTPCTrack2Beam; + // check track transveres cuts + float sna{0}, csa{0}; + o2::math_utils::CircleXYf_t trkCircle; + trLoc.getCircleParams(mBz, trkCircle, sna, csa); + float cR = std::hypot(trkCircle.xC, trkCircle.yC); + float drd2 = std::sqrt((cR * cR) - (trkCircle.rC * trkCircle.rC)); + bool dRD2 = drd2 > mSVParams->mTPCTrackXY2Radius; + + { + bool ok{false}; + const auto& pvLbl = mRecoData.getPrimaryVertexMCLabel(vtxid); + if (pvLbl.getEventID() == mRecoData.getTrackMCLabel(gid).getEventID()) { + ok = true; + } + (*mDBGTree) << "tpcPhotonTune" + << "dCls=" << dCls + << "pt=" << trTPC.getPt() + << "zPVExtrapDiff=" << zPVExtrapDiff + << "dDPV=" << dDPV + << "trkCircle=" << trkCircle + << "cR=" << cR + << "drd2=" << drd2 + << "dRD2=" << dRD2 + << "ok=" << ok + << "\n"; + } + + if (dCls) { + mTracksPool[posneg].pop_back(); + updatePrInfo(dataformats::PoolInfo::RejTPCPhotonCls); + return true; + } + if (dDPV) { + mTracksPool[posneg].pop_back(); + updatePrInfo(dataformats::PoolInfo::RejTPCPhotonDPV); + return true; + } + if (dRD2) { + mTracksPool[posneg].pop_back(); + updatePrInfo(dataformats::PoolInfo::RejTPCPhotonRD2); + return true; + } + } + + updatePrInfo(dataformats::PoolInfo::Accepted); + return true; +} + +float SVMCStudySpec::correctTPCTrack(TrackCand& trc, const o2::tpc::TrackTPC& tTPC, float tmus, float tmusErr) const +{ + // Correct the track copy trc of the TPC track for the assumed interaction time + // return extra uncertainty in Z due to the interaction time uncertainty + float tTB = NAN, tTBErr = NAN; + if (tmusErr < 0) { // use track data + tTB = tTPC.getTime0(); + tTBErr = 0.5f * (tTPC.getDeltaTBwd() + tTPC.getDeltaTFwd()); + } else { + tTB = tmus * mMUS2TPCBin; + tTBErr = tmusErr * mMUS2TPCBin; + } + float dDrift = (tTB - tTPC.getTime0()) * mTPCBin2Z; + float driftErr = tTBErr * mTPCBin2Z; + if (driftErr < 0.) { // early return will be discarded anyway + return driftErr; + } + trc.setZ(tTPC.getZ() + (tTPC.hasASideClustersOnly() ? dDrift : -dDrift)); + trc.setCov(trc.getSigmaZ2() + (driftErr * driftErr), o2::track::kSigZ2); + uint8_t sector = 0, row = 0; + auto cl = &tTPC.getCluster(mTPCTrackClusIdx, tTPC.getNClusters() - 1, *mTPCClusterIdxStruct, sector, row); + float x = 0, y = 0, z = 0; + mTPCCorrMapsLoader.Transform(sector, row, cl->getPad(), cl->getTime(), x, y, z, tTB); + x = std::max(x, o2::constants::geom::XTPCInnerRef); + trc.minR = std::hypot(x, y); + return driftErr; +} + +bool SVMCStudySpec::checkV0(const TrackCand& seedP, const TrackCand& seedN) +{ + auto lblP = mRecoData.getTrackMCLabel(seedP.gid); + auto lblN = mRecoData.getTrackMCLabel(seedN.gid); + dataformats::V0MCExt* moth{nullptr}; + if (mDau2V0.find(lblP) != mDau2V0.end() && mDau2V0.find(lblN) != mDau2V0.end() && mDau2V0[lblP] == mDau2V0[lblN]) { + moth = &mV0MC[mDau2V0[lblN]]; + } + // Update pairStatus to the deepest rejection reached across all calls for this physical pair. + // Since TPC-only tracks have one pool entry per vertex, checkV0 may be called multiple times + // for the same physical pair; we keep the maximum (furthest-along) status. + auto updateMoth = [&](dataformats::V0MCExt::PairStatus s) { + if (moth && s > moth->pairStatus) { + moth->pairStatus = s; + } + }; + + // Fast rough cuts on pairs before feeding to DCAFitter, tracks are not in the same Frame or at same X + bool isTPConly = (seedP.gid.getSource() == GTrackID::TPC || seedN.gid.getSource() == GTrackID::TPC); + if (mSVParams->mTPCTrackPhotonTune && isTPConly) { + // Check if Tgl is close enough + if (std::abs(seedP.getTgl() - seedN.getTgl()) > mSVParams->maxV0TglAbsDiff) { + LOG(debug) << "RejTgl"; + updateMoth(dataformats::V0MCExt::RejTgl); + return false; + } + // Check in transverse plane + float sna = NAN, csa = NAN; + o2::math_utils::CircleXYf_t trkPosCircle; + seedP.getCircleParams(mBz, trkPosCircle, sna, csa); + o2::math_utils::CircleXYf_t trkEleCircle; + seedN.getCircleParams(mBz, trkEleCircle, sna, csa); + // Does the radius of both tracks compare to their absolute circle center distance + float c2c = std::hypot(trkPosCircle.xC - trkEleCircle.xC, + trkPosCircle.yC - trkEleCircle.yC); + float r2r = trkPosCircle.rC + trkEleCircle.rC; + float dcr = c2c - r2r; + if (std::abs(dcr) > mSVParams->mTPCTrackD2R) { + LOG(debug) << "RejD2R " << c2c << " " << r2r << " " << dcr; + updateMoth(dataformats::V0MCExt::RejD2R); + return false; + } + // Will the conversion point look somewhat reasonable + float r1_r = trkPosCircle.rC / r2r; + float r2_r = trkEleCircle.rC / r2r; + float dR = std::hypot((r2_r * trkPosCircle.xC) + (r1_r * trkEleCircle.xC), (r2_r * trkPosCircle.yC) + (r1_r * trkEleCircle.yC)); + if (dR > mSVParams->mTPCTrackDR) { + LOG(debug) << "RejDR" << dR; + updateMoth(dataformats::V0MCExt::RejDR); + return false; + } + + // Setup looser cuts for the DCAFitter + mFitterV0.setMaxDZIni(mSVParams->mTPCTrackMaxDZIni); + mFitterV0.setMaxDXYIni(mSVParams->mTPCTrackMaxDXYIni); + mFitterV0.setMaxChi2(mSVParams->mTPCTrackMaxChi2); + mFitterV0.setCollinear(true); + } + + // feed DCAFitter + int nCand = mFitterV0.process(seedP, seedN); + if (mSVParams->mTPCTrackPhotonTune && isTPConly) { + // Reset immediately to the defaults + mFitterV0.setMaxDZIni(mSVParams->maxDZIni); + mFitterV0.setMaxDXYIni(mSVParams->maxDXYIni); + mFitterV0.setMaxChi2(mSVParams->maxChi2); + mFitterV0.setCollinear(false); + } + if (nCand == 0) { // discard this pair + LOG(debug) << "RejDCAFitter"; + updateMoth(dataformats::V0MCExt::RejDCAFitter); + return false; + } + + const auto& v0XYZ = mFitterV0.getPCACandidate(); + // validate V0 radial position + // check closeness to the beam-line + float dxv0 = v0XYZ[0] - mMeanVertex->getX(), dyv0 = v0XYZ[1] - mMeanVertex->getY(), r2v0 = (dxv0 * dxv0) + (dyv0 * dyv0); + if (r2v0 < mMinR2ToMeanVertex) { + LOG(debug) << "RejMinR2ToMeanVertex"; + updateMoth(dataformats::V0MCExt::RejMinR); + return false; + } + + float rv0 = std::sqrt(r2v0), drv0P = rv0 - seedP.minR, drv0N = rv0 - seedN.minR; + if (drv0P > mSVParams->causalityRTolerance || drv0P < -mSVParams->maxV0ToProngsRDiff || + drv0N > mSVParams->causalityRTolerance || drv0N < -mSVParams->maxV0ToProngsRDiff) { + LOG(debug) << "RejCausality " << drv0P << " " << drv0N; + updateMoth(dataformats::V0MCExt::RejCausality); + return false; + } + + const int cand = 0; + if (!mFitterV0.isPropagateTracksToVertexDone(cand) && !mFitterV0.propagateTracksToVertex(cand)) { + LOG(debug) << "RejProp failed"; + updateMoth(dataformats::V0MCExt::RejPropagation); + return false; + } + + const auto& trPProp = mFitterV0.getTrack(0, cand); + const auto& trNProp = mFitterV0.getTrack(1, cand); + std::array pP{}, pN{}; + trPProp.getPxPyPzGlo(pP); + trNProp.getPxPyPzGlo(pN); + // estimate DCA of neutral V0 track to beamline: straight line with parametric equation + // x = X0 + pV0[0]*t, y = Y0 + pV0[1]*t reaches DCA to beamline (Xv, Yv) at + // t = -[ (x0-Xv)*pV0[0] + (y0-Yv)*pV0[1]) ] / ( pT(pV0)^2 ) + // Similar equation for 3D distance involving pV0[2] + std::array pV0 = {pP[0] + pN[0], pP[1] + pN[1], pP[2] + pN[2]}; + float pt2V0 = (pV0[0] * pV0[0]) + (pV0[1] * pV0[1]), prodXYv0 = (dxv0 * pV0[0]) + (dyv0 * pV0[1]), tDCAXY = prodXYv0 / pt2V0; + if (pt2V0 < mMinPt2V0) { // pt cut + LOG(debug) << "RejPt2 " << pt2V0; + updateMoth(dataformats::V0MCExt::RejPt); + return false; + } + if (pV0[2] * pV0[2] / pt2V0 > mMaxTgl2V0) { // tgLambda cut + LOG(debug) << "RejTgL " << pV0[2] * pV0[2] / pt2V0; + updateMoth(dataformats::V0MCExt::RejTgLambda); + return false; + } + float p2V0 = pt2V0 + (pV0[2] * pV0[2]), ptV0 = std::sqrt(pt2V0); + // apply mass selections + float p2Pos = (pP[0] * pP[0]) + (pP[1] * pP[1]) + (pP[2] * pP[2]), p2Neg = (pN[0] * pN[0]) + (pN[1] * pN[1]) + (pN[2] * pN[2]); + + bool goodHyp = false, photonOnly = mSVParams->mTPCTrackPhotonTune && isTPConly; + std::array hypCheckStatus{}; + int nPID = photonOnly ? (HypV0::Photon + 1) : HypV0::NHypV0; + for (int ipid = 0; (ipid < nPID) && mSVParams->checkV0Hypothesis; ipid++) { + if (mV0Hyps[ipid].check(p2Pos, p2Neg, p2V0, ptV0)) { + goodHyp = hypCheckStatus[ipid] = true; + } + } + + // check tight lambda mass only + if (!goodHyp && mSVParams->checkV0Hypothesis) { + updateMoth(dataformats::V0MCExt::RejHypothesis); + return false; + } + + float dcaX = dxv0 - (pV0[0] * tDCAXY), dcaY = dyv0 - (pV0[1] * tDCAXY), dca2 = (dcaX * dcaX) + (dcaY * dcaY); + float cosPAXY = prodXYv0 / std::sqrt(r2v0 * pt2V0); + + if (dca2 > mMaxDCAXY2ToMeanVertex || cosPAXY < mSVParams->minCosPAXYMeanVertex) { + if (mSVParams->mTPCTrackPhotonTune && isTPConly) { + // Check for looser cut for tpc-only photons only + if (dca2 > mSVParams->mTPCTrackMaxDCAXY2ToMeanVertex) { + updateMoth(dataformats::V0MCExt::RejDCAToMV); + return false; + } + } else { + updateMoth(dataformats::V0MCExt::RejDCAToMV); + return false; + } + } + + auto vlist = seedP.vBracket.getOverlap(seedN.vBracket); // indices of vertices shared by both seeds + bool candFound = false; + auto bestCosPA = mSVParams->minCosPA; + dataformats::V0 v0new; + + for (int iv = vlist.getMin(); iv <= vlist.getMax(); iv++) { + const auto& pv = mPVertices[iv]; + const auto v0XYZ = mFitterV0.getPCACandidatePos(cand); + // check cos of pointing angle + float dx = v0XYZ[0] - pv.getX(), dy = v0XYZ[1] - pv.getY(), dz = v0XYZ[2] - pv.getZ(), prodXYZv0 = (dx * pV0[0]) + (dy * pV0[1]) + (dz * pV0[2]); + float cosPA = prodXYZv0 / std::sqrt((dx * dx + dy * dy + dz * dz) * p2V0); + if (cosPA < bestCosPA) { + LOG(debug) << "Rej. cosPA: " << cosPA; + continue; + } + if (!candFound) { + new (&v0new) dataformats::V0(v0XYZ, pV0, mFitterV0.calcPCACovMatrixFlat(cand), trPProp, trNProp); + v0new.setDCA(mFitterV0.getChi2AtPCACandidate(cand)); + candFound = true; + } + v0new.setCosPA(cosPA); + bestCosPA = cosPA; + } + if (!candFound) { + updateMoth(dataformats::V0MCExt::RejCosPA); + return false; + } + + dataformats::V0MCFull v0; + v0.hypStatus = hypCheckStatus; + v0.v0 = v0new; + if (moth) { + v0.label = mDau2V0[lblP]; + v0.pdg = moth->mcTrk.GetPdgCode(); + } + updateMoth(dataformats::V0MCExt::Found); + mV0s.push_back(v0); + + return true; +} + +o2::framework::DataProcessorSpec getSVMCStudySpec(GTrackID::mask_t srcTracks, GTrackID::mask_t srcCls, const o2::tpc::CorrectionMapsLoaderGloOpts& sclOpts) +{ + o2::framework::Options opts; + std::vector outputs; + auto dataRequest = std::make_shared(); + + dataRequest->requestTracks(srcTracks, true); + dataRequest->requestClusters(srcCls, true); + dataRequest->requestPrimaryVertices(true); + dataRequest->inputs.emplace_back("meanvtx", "GLO", "MEANVERTEX", 0, o2::framework::Lifetime::Condition, o2::framework::ccdbParamSpec("GLO/Calib/MeanVertex", {}, 1)); + dataRequest->inputs.emplace_back("SVParam", "GLO", "SVPARAM", 0, o2::framework::Lifetime::Condition, o2::framework::ccdbParamSpec("GLO/Config/SVertexerParam")); + auto ggRequest = std::make_shared(true, // orbitResetTime + true, // GRPECS=true + false, // GRPLHCIF + true, // GRPMagField + true, // askMatLUT + o2::base::GRPGeomRequest::Aligned, // geometry + dataRequest->inputs, + true); + if (srcTracks[GTrackID::TPC]) { + o2::tpc::VDriftHelper::requestCCDBInputs(dataRequest->inputs); + o2::tpc::CorrectionMapsLoader::requestCCDBInputs(dataRequest->inputs, opts, sclOpts); + } + return o2::framework::DataProcessorSpec{ + .name = "sv-mc-study", + .inputs = dataRequest->inputs, + .outputs = outputs, + .algorithm = o2::framework::AlgorithmSpec{o2::framework::adaptFromTask(dataRequest, ggRequest, sclOpts, srcTracks)}, + .options = opts}; +} + +} // namespace o2::svstudy diff --git a/Detectors/GlobalTrackingWorkflow/study/src/sv-mc-study-workflow.cxx b/Detectors/GlobalTrackingWorkflow/study/src/sv-mc-study-workflow.cxx new file mode 100644 index 0000000000000..9657ec3d71e36 --- /dev/null +++ b/Detectors/GlobalTrackingWorkflow/study/src/sv-mc-study-workflow.cxx @@ -0,0 +1,69 @@ +// 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. + +// Simplified version of the SVertexer code with dumpable mc output +// to study inefficiencies and cut tuning. + +#include "GlobalTrackingStudy/SVMCStudy.h" +#include "ReconstructionDataFormats/GlobalTrackID.h" +#include "CommonUtils/ConfigurableParam.h" +#include "Framework/ConfigParamSpec.h" +#include "Framework/CallbacksPolicy.h" +#include "GlobalTrackingWorkflowHelpers/InputHelper.h" +#include "DetectorsRaw/HBFUtilsInitializer.h" +#include "TPCWorkflow/TPCScalerSpec.h" +#include "TPCCalibration/CorrectionMapsLoader.h" + +using namespace o2::framework; +using GID = o2::dataformats::GlobalTrackID; + +void customize(std::vector& policies) +{ + o2::raw::HBFUtilsInitializer::addNewTimeSliceCallback(policies); +} + +void customize(std::vector& workflowOptions) +{ + // option allowing to set parameters + std::vector options{ + {"track-sources", VariantType::String, std::string{GID::ALL}, {"comma-separated list of track sources to use"}}, + {"disable-root-input", VariantType::Bool, false, {"disable root-files input reader"}}, + {"disable-mc", VariantType::Bool, false, {"disable mc ?"}}, + {"configKeyValues", VariantType::String, "", {"Semicolon separated key=value strings ..."}}}; + o2::tpc::CorrectionMapsLoader::addGlobalOptions(options); + o2::raw::HBFUtilsInitializer::addConfigOption(options); + std::swap(workflowOptions, options); +} + +#include "Framework/runDataProcessing.h" + +WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) +{ + auto sclOpt = o2::tpc::CorrectionMapsLoader::parseGlobalOptions(configcontext.options()); + GID::mask_t allowedSourcesTrc = GID::getSourcesMask("ITS,TPC,ITS-TPC,TPC-TOF,TPC-TRD,ITS-TPC-TRD,TPC-TRD-TOF,ITS-TPC-TOF,ITS-TPC-TRD-TOF"); + GID::mask_t allowedClsTrc = GID::getSourcesMask("ITS,TPC"); + o2::conf::ConfigurableParam::updateFromString(configcontext.options().get("configKeyValues")); + GID::mask_t srcTrc = allowedSourcesTrc & GID::getSourcesMask(configcontext.options().get("track-sources")); + GID::mask_t srcCls = allowedClsTrc & srcTrc; + if (sclOpt.requestCTPLumi) { + srcTrc = srcTrc | GID::getSourcesMask("CTP"); + } + + WorkflowSpec specs; + o2::globaltracking::InputHelper::addInputSpecs(configcontext, specs, srcCls, srcTrc, srcTrc, true); + o2::globaltracking::InputHelper::addInputSpecsPVertex(configcontext, specs, true); // P-vertex is always needed + if (sclOpt.needTPCScalersWorkflow() && !configcontext.options().get("disable-root-input")) { + specs.emplace_back(o2::tpc::getTPCScalerSpec(sclOpt.lumiType == 2, sclOpt.enableMShapeCorrection)); + } + specs.emplace_back(o2::svstudy::getSVMCStudySpec(srcTrc, srcCls, sclOpt)); + o2::raw::HBFUtilsInitializer hbfIni(configcontext, specs); + return std::move(specs); +}