diff --git a/Detectors/TPC/calibration/SpacePoints/CMakeLists.txt b/Detectors/TPC/calibration/SpacePoints/CMakeLists.txt index 47bb9c09a9951..ac33a0b632dab 100644 --- a/Detectors/TPC/calibration/SpacePoints/CMakeLists.txt +++ b/Detectors/TPC/calibration/SpacePoints/CMakeLists.txt @@ -43,9 +43,40 @@ o2_add_test_root_macro(macro/staticMapCreator.C PUBLIC_LINK_LIBRARIES O2::SpacePoints LABELS tpc COMPILE_ONLY) +o2_add_test_root_macro(macro/staticMapCreatorCPM.C + PUBLIC_LINK_LIBRARIES O2::SpacePoints + O2::CCDB + O2::Algorithm + O2::Framework + O2::CommonConstants + O2::DataFormatsParameters + O2::ReconstructionDataFormats + O2::DetectorsBase + LABELS tpc COMPILE_ONLY) + +o2_add_test_root_macro(macro/SmoothingExtrapolate.C + PUBLIC_LINK_LIBRARIES O2::SpacePoints + O2::Algorithm + O2::TPCCalibration + O2::CCDB + O2::TPCBaseRecSim + LABELS tpc COMPILE_ONLY) + +o2_add_test_root_macro(macro/voxResQA.C + PUBLIC_LINK_LIBRARIES O2::SpacePoints + O2::TPCBaseRecSim + O2::CommonUtils + O2::Framework + LABELS tpc COMPILE_ONLY) + install(FILES macro/staticMapCreator.C DESTINATION share/macro/) +install(FILES macro/staticMapCreatorCPM.C + macro/SmoothingExtrapolate.C + macro/voxResQA.C + DESTINATION share/macro/) + o2_add_test(TrackResiduals COMPONENT_NAME calibration PUBLIC_LINK_LIBRARIES O2::SpacePoints diff --git a/Detectors/TPC/calibration/SpacePoints/macro/SmoothingExtrapolate.C b/Detectors/TPC/calibration/SpacePoints/macro/SmoothingExtrapolate.C new file mode 100644 index 0000000000000..2af42e68898ea --- /dev/null +++ b/Detectors/TPC/calibration/SpacePoints/macro/SmoothingExtrapolate.C @@ -0,0 +1,1002 @@ +// 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. + +#if !defined(__CLING__) || defined(__ROOTCLING__) + +#include +#include +#include +#include "TFile.h" +#include "TSystem.h" +#include "TTree.h" +#include "TF1.h" +#include "TGraph.h" +#include "TProfile.h" +#include "Math/MinimizerOptions.h" + +#include "Algorithm/RangeTokenizer.h" +#include "SpacePoints/TrackResiduals.h" +#include "TPCCalibration/TPCFastSpaceChargeCorrectionHelper.h" +#include "CCDB/BasicCCDBManager.h" +#include "TPCCalibration/TPCScaler.h" +#if __has_include("TPCBaseRecSim/CDBTypes.h") +#include "TPCBaseRecSim/CDBTypes.h" +#else +#include "TPCBase/CDBTypes.h" +#endif + +#endif + +using namespace o2::tpc; +using namespace o2::gpu; + +//------------------------------------------------------------------------------------------------------------ +static const Float_t RowX[153] = + { + 85.225, 85.975, 86.725, 87.475, 88.225, 88.975, 89.725, 90.475, 91.225, 91.975, 92.725, 93.475, 94.225, 94.975, 95.725, 96.475, + 97.225, 97.975, 98.725, 99.475, 100.225, 100.975, 101.725, 102.475, 103.225, 103.975, 104.725, 105.475, 106.225, 106.975, 107.725, + 108.475, 109.225, 109.975, 110.725, 111.475, 112.225, 112.975, 113.725, 114.475, 115.225, 115.975, 116.725, 117.475, 118.225, 118.975, + 119.725, 120.475, 121.225, 121.975, 122.725, 123.475, 124.225, 124.975, 125.725, 126.475, 127.225, 127.975, 128.725, 129.475, 130.225, + 130.975, 131.725, 135.200, 136.200, 137.200, 138.200, 139.200, 140.200, 141.200, 142.200, 143.200, 144.200, 145.200, 146.200, 147.200, + 148.200, 149.200, 150.200, 151.200, 152.200, 153.200, 154.200, 155.200, 156.200, 157.200, 158.200, 159.200, 160.200, 161.200, 162.200, + 163.200, 164.200, 165.200, 166.200, 167.200, 168.200, 171.400, 172.600, 173.800, 175.000, 176.200, 177.400, 178.600, 179.800, 181.000, + 182.200, 183.400, 184.600, 185.800, 187.000, 188.200, 189.400, 190.600, 191.800, 193.000, 194.200, 195.400, 196.600, 197.800, 199.000, + 200.200, 201.400, 202.600, 203.800, 205.000, 206.200, 209.650, 211.150, 212.650, 214.150, 215.650, 217.150, 218.650, 220.150, 221.650, + 223.150, 224.650, 226.150, 227.650, 229.150, 230.650, 232.150, 233.650, 235.150, 236.650, 238.150, 239.650, 241.150, 242.650, 244.150, + 245.650, 246.650}; // last value added +//------------------------------------------------------------------------------------------------------------ + +//---------------------------------------------------------------------------------------- +Double_t PolyFitFunc(Double_t* x_val, Double_t* par) +{ + Double_t x, y, par0, par1, par2, par3, par4, par5; + par0 = par[0]; + par1 = par[1]; + par2 = par[2]; + par3 = par[3]; + par4 = par[4]; + par5 = par[5]; + x = x_val[0]; + y = par0 + par1 * x + par2 * x * x + par3 * x * x * x + par4 * x * x * x * x + par5 * x * x * x * x * x; + return y; +} +//---------------------------------------------------------------------------------------- + +// Function to create Gaussian filter +void vec_FilterCreation(std::vector>>& vec_GKernel, Int_t Delta_X, Int_t Delta_Y, Int_t Delta_Z, Double_t sigma) +{ + // initialising standard deviation to 1.0 + // double sigma = 1.0; + double r, s = 2.0 * sigma * sigma; + + // sum is for normalization + double sum = 0.0; + + // generating 5x5 kernel + for (int x = -Delta_X; x <= Delta_X; x++) { + for (int y = -Delta_Y; y <= Delta_Y; y++) { + for (int z = -Delta_Z; z <= Delta_Z; z++) { + r = sqrt(x * x + y * y + z * z); + vec_GKernel[x + Delta_X][y + Delta_Y][z + Delta_Z] = (exp(-(r * r) / s)) / (M_PI * s); + sum += vec_GKernel[x + Delta_X][y + Delta_Y][z + Delta_Z]; + } + } + } + + // normalising the Kernel + for (int i = 0; i < (Delta_X * 2 + 1); ++i) { + for (int j = 0; j < (Delta_Y * 2 + 1); ++j) { + for (int k = 0; k < (Delta_Z * 2 + 1); ++k) { + vec_GKernel[i][j][k] /= sum; + } + } + } +} + +// Mean TPC scaler (IDC) values for one timestamp, from the standard CCDB CalScaler/CalScalerWeights +// objects. Used for the offline IDC join below -- see the "Offline IDC join" block in +// SmoothingExtrapolate() for why this is done here rather than during map creation. +bool getScalerValues(o2::ccdb::BasicCCDBManager& ccdbmgr, long tfTimeInMS, float& scA, float& scC) +{ + auto* scalerTree = ccdbmgr.getForTimeStamp(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalScaler), long(std::ceil(tfTimeInMS))); + auto* scalerWeights = ccdbmgr.getForTimeStamp(o2::tpc::CDBTypeMap.at(o2::tpc::CDBType::CalScalerWeights), long(std::ceil(tfTimeInMS))); + + if (!scalerTree) { + LOGP(error, "Could not get 'TPC/Calib/Scaler' for time stamp {}", tfTimeInMS); + return false; + } + // The caller sets setFatalWhenNull(false), so a missing object comes back as nullptr rather than + // aborting -- the weights have to be checked before being dereferenced below. + if (!scalerWeights) { + LOGP(error, "Could not get 'TPC/Calib/ScalerWeights' for time stamp {}", tfTimeInMS); + return false; + } + + o2::tpc::TPCScaler scaler; + scaler.setFromTree(*(scalerTree)); + scaler.setScalerWeights(*scalerWeights); + scaler.useWeights(true); + scaler.setIonDriftTimeMS(500); + + static bool defaultScalerReported = false; + static bool badScalerValueReported = false; + if (scaler.getRun() == 0) { + if (!defaultScalerReported) { + LOGP(error, "Retrieved default scaler entry 'TPC/Calib/Scaler' for time stamp {}", tfTimeInMS); + defaultScalerReported = true; + } + return false; + } + + scA = scaler.getMeanScaler(tfTimeInMS, o2::tpc::Side::A); + scC = scaler.getMeanScaler(tfTimeInMS, o2::tpc::Side::C); + + if ((scA <= 0) || (scC <= 0)) { + if (!badScalerValueReported) { + LOGP(error, "Bad scaler value, first seen for time stamp {}, scA: {}, scC: {}", tfTimeInMS, scA, scC); + badScalerValueReported = true; + } + scA = 0; + scC = 0; + return false; + } + + return true; +} + +void SmoothingExtrapolate(const char* fileName = "debugVoxRes.root", TString fileOutName = "SmoothVoxRes.root", + int do_smoothing = 1, int do_extrapolation = 2, + int A11maxZ2X = -1, bool maskIA11 = false, + Int_t N_bins_X_GF = 1, Int_t N_bins_Y_GF = 2, + Int_t N_bins_Z_GF = 0, Float_t sigma_GF = 1.2) +{ + + // do_extrapolation: + // 0 -> no extrapolation to low radii done + // 1 -> only smoothed values are extrapolated + // 2 -> smoothed and raw values are extrapolated + // 3 -> only raw values are extrapolated + + // Example, smoothing and extrapolating both smoothed and raw values: + // SmoothingExtrapolate("voxRes.__.it0.root", "voxRes._smooth.root", 1, 2, -1, false, 1, 2, 0, 1.2) + + const float maxDeltaCut = 25.0; // maximum value for any Delta to be accepted + const float min_statistics = 20; + const float max_extrapolation_value = 20.0; // maximum value for extrapolation in DX, DY, DZ + //---------------------------------------------------------------- + // input + if (gSystem->AccessPathName(fileName)) { + LOGP(error, "input file {} does not exist", fileName); + return; + } + + auto file = std::unique_ptr(TFile::Open(fileName, "READ")); + if (!file || !file->IsOpen()) { + LOGP(error, "input file {} does not exist", fileName); + return; + } + + TTree* voxResTree = nullptr; + file->cd(); + gDirectory->GetObject("voxResTree", voxResTree); + if (!voxResTree) { + LOGP(error, "tree voxResTree does not exist in {}", fileName); + return; + } + + o2::tpc::TrackResiduals::VoxRes* voxRes_map = nullptr; + Long64_t entries_input_map = voxResTree->GetEntries(); + LOGP(info, "entries_input_map: {}", entries_input_map); + voxResTree->SetBranchAddress("voxRes", &voxRes_map); + + // required for the binning that was used + auto userInfo = voxResTree->GetUserInfo(); + if (!userInfo->FindObject("y2xBinning") || !userInfo->FindObject("z2xBinning")) { + LOGP(error, "'y2xBinning' or 'z2xBinning' not found in UserInfo, but required to get the correct binning"); + return; + } + + // Obtain configuration + const SpacePointsCalibConfParam& params = SpacePointsCalibConfParam::Instance(); + if (std::filesystem::exists("scdconfig.ini")) { + params.updateFromFile("scdconfig.ini"); + } + // TrackResiduals::setZ2XBinning() (called below) reads scdcalib.maxZ2X directly and uses it to scale + // the physical z/x bin boundaries -- it is baked into what each z2x voxel index in the input tree + // actually means, not a cosmetic knob. There is no scdconfig.ini on the GRID, so without this the + // code default (1.0) would silently apply instead of whatever stage 1 actually used (production 1.4), + // misaligning this macro's re-derived binning against the tree's real geometry. Must happen BEFORE + // setZ2XBinning() below. See staticMapCreatorCPM.C's UserInfo::Add("maxZ2X", ...) for where this comes + // from. + if (auto* maxZ2XObj = userInfo->FindObject("maxZ2X")) { + const std::string maxZ2XStr = maxZ2XObj->GetTitle(); + o2::conf::ConfigurableParam::setValue("scdcalib.maxZ2X", maxZ2XStr); + LOGP(info, "Set scdcalib.maxZ2X = {} from input UserInfo (matches stage 1)", maxZ2XStr); + } else { + LOGP(warning, + "'maxZ2X' not found in input UserInfo (older input file?) -- using scdcalib.maxZ2X = {} " + "(scdconfig.ini/code default), which may NOT match the value stage 1 actually used to " + "build this tree's z2x binning!", + params.maxZ2X); + } + + LOGP(info, "----- Dumping configuration values START -----"); + params.printKeyValues(); + LOGP(info, "----- Dumping configuration values END -----"); + + LOGP(info, "Get binning from userInfo"); + o2::tpc::TrackResiduals trackResiduals; + auto y2xBins = o2::RangeTokenizer::tokenize(userInfo->FindObject("y2xBinning")->GetTitle()); + auto z2xBins = o2::RangeTokenizer::tokenize(userInfo->FindObject("z2xBinning")->GetTitle()); + trackResiduals.setY2XBinning(y2xBins); + trackResiduals.setZ2XBinning(z2xBins); + trackResiduals.init(); + LOGP(info, "trackResiduals initialized"); + //---------------------------------------------------------------- + + //---------------------------------------------------------------- + // Offline IDC join: staticMapCreatorCPM.C on the GRID intentionally does not access IDCs + // live (useCTPLumi=2 -- CCDB access from inside the hot per-TF loop was slow/unreliable), so + // meanIDC/medianIDC in the input's UserInfo are placeholder 0. Per-TF timestamps are recorded + // regardless (the "OrbitLumiInfo" tree's timeMSsel branch) specifically so this can be joined back + // offline, with a properly-cached CCDB client. Timestamps are sorted before querying: + // TPCScaler/TPCScalerWeights CCDB objects are valid over a time range, and BasicCCDBManager's cache + // (setCaching(true)) only hits when consecutive queries land in the same validity window -- + // unsorted access would bounce between windows and force a real CCDB fetch almost every call. + float meanIDCReal = 0.f; + float medianIDCReal = 0.f; + { + TTree* orbitLumiTree = nullptr; + file->GetObject("OrbitLumiInfo", orbitLumiTree); + if (!orbitLumiTree || orbitLumiTree->GetEntries() == 0) { + LOGP(warning, "Offline IDC join: no 'OrbitLumiInfo' tree (or it's empty) in the input file -- meanIDC/medianIDC stay 0"); + } else { + std::vector* timeMSselPtr = nullptr; + orbitLumiTree->SetBranchAddress("timeMSsel", &timeMSselPtr); + orbitLumiTree->GetEntry(0); + if (!timeMSselPtr || timeMSselPtr->empty()) { + LOGP(warning, "Offline IDC join: 'timeMSsel' branch missing or empty -- meanIDC/medianIDC stay 0"); + } else { + std::vector sortedTimes(*timeMSselPtr); + std::sort(sortedTimes.begin(), sortedTimes.end()); + + auto& ccdbmgr = o2::ccdb::BasicCCDBManager::instance(); + ccdbmgr.setCaching(true); + ccdbmgr.setFatalWhenNull(false); + ccdbmgr.setURL("http://alice-ccdb.cern.ch"); + + std::vector averageIDCs; + averageIDCs.reserve(sortedTimes.size()); + for (long tfTimeInMS : sortedTimes) { + float scA = 0.f, scC = 0.f; + if (getScalerValues(ccdbmgr, tfTimeInMS, scA, scC)) { + averageIDCs.emplace_back((scA + scC) / 2.f); + } + } + if (averageIDCs.empty()) { + LOGP(warning, "Offline IDC join: no valid scaler values found for any of {} TFs -- meanIDC/medianIDC stay 0", sortedTimes.size()); + } else { + double sum = 0.0; + for (float v : averageIDCs) { + sum += v; + } + meanIDCReal = static_cast(sum / averageIDCs.size()); + medianIDCReal = static_cast(TMath::Median(static_cast(averageIDCs.size()), averageIDCs.data())); + LOGP(info, "Offline IDC join: {} of {} TFs gave a valid scaler value, meanIDC={}, medianIDC={}", + averageIDCs.size(), sortedTimes.size(), meanIDCReal, medianIDCReal); + } + } + } + } + //---------------------------------------------------------------- + + //---------------------------------------------------------------- + // Get voxel map binning from map + + const int nXBins = trackResiduals.getNXBins(); + const int nY2XBins = trackResiduals.getNY2XBins(); + const int nZ2XBins = trackResiduals.getNZ2XBins(); + LOGP(info, "binning X,Y2X,Z2X: {}, {}, {}", nXBins, nY2XBins, nZ2XBins); + //---------------------------------------------------------------- + + //---------------------------------------------------------------- + // output + // Create a new file + a clone of old tree in new file + TFile* outputfile = new TFile(fileOutName.Data(), "RECREATE"); + LOGP(info, "Output file: {} created", fileOutName.Data()); + + o2::tpc::TrackResiduals::VoxRes mVoxelResultsOut{}; ///< the results from mVoxelResults are copied in here to be able to stream them + o2::tpc::TrackResiduals::VoxRes* mVoxelResultsOutPtr{&mVoxelResultsOut}; ///< pointer to set the branch address to for the output + std::unique_ptr mTreeOut; + //---------------------------------------------------------------- + + //---------------------------------------------------------------- + ROOT::Math::MinimizerOptions::SetDefaultMinimizer("GSLSimAn"); + + TF1* func_PolyFitFunc = new TF1("func_PolyFitFunc", PolyFitFunc, 0, 150, 6); + TF1* func_PolyFitFunc_raw = new TF1("func_PolyFitFunc_raw", PolyFitFunc, 0, 150, 6); + TProfile* tp_DX_vs_X_raw = new TProfile("tp_DX_vs_X_raw", "tp_DX_vs_X_raw", 250, 0, 250); + TProfile* tp_Stat_vs_X = new TProfile("tp_Stat_vs_X", "tp_Stat_vs_X;row;", 250, 0, 250); + TProfile* tp_Stat_vs_X_single = new TProfile("tp_Stat_vs_X_single", "tp_Stat_vs_X;row;", 250, 0, 250); + TProfile* tp_DX_vs_X_single = new TProfile("tp_DX_vs_X_single", "tp_DX_vs_X;row;", 250, 0, 250); + TGraph* tg_Stat_vs_X_slice = new TGraph(); + TProfile* tp_DX_vs_Row_raw = new TProfile("tp_DX_vs_Row_raw", "tp_DX_vs_Row_raw", 250, 0, 250); + TProfile* tp_DX_vs_X_smooth = new TProfile("tp_DX_vs_X_smooth", "tp_DX_vs_X_smooth", 250, 0, 250); + TProfile* tp_DX_vs_X_smooth_extr = new TProfile("tp_DX_vs_X_smooth_extr", "tp_DX_vs_X_smooth_extr", 250, 0, 250); + TProfile* tp_DY_vs_X_smooth_extr = new TProfile("tp_DY_vs_X_smooth_extr", "tp_DY_vs_X_smooth_extr", 250, 0, 250); + TProfile* tp_DZ_vs_X_smooth_extr_A = new TProfile("tp_DZ_vs_X_smooth_extr_A", "tp_DZ_vs_X_smooth_extr_A", 250, 0, 250); + TProfile* tp_DZ_vs_X_smooth_extr_C = new TProfile("tp_DZ_vs_X_smooth_extr_C", "tp_DZ_vs_X_smooth_extr_C", 250, 0, 250); + TProfile* tp_DZ_vs_X_smooth_A = new TProfile("tp_DZ_vs_X_smooth_A", "tp_DZ_vs_X_smooth_A", 250, 0, 250); + TProfile* tp_DZ_vs_X_smooth_C = new TProfile("tp_DZ_vs_X_smooth_C", "tp_DZ_vs_X_smooth_C", 250, 0, 250); + TProfile* tp_DZ_vs_X_raw_A = new TProfile("tp_DZ_vs_X_raw_A", "tp_DZ_vs_X_raw_A", 250, 0, 250); + TProfile* tp_DZ_vs_X_raw_C = new TProfile("tp_DZ_vs_X_raw_C", "tp_DZ_vs_X_raw_C", 250, 0, 250); + TProfile* tp_Stat_vs_row = new TProfile("tp_Stat_vs_row", "tp_Stat_vs_row;row;", 152, 0, 152); + int n_bins_z_phi_sector = 36 * nY2XBins * nZ2XBins; + TH1D* h_x_start_fit_vs_z_phi_sector = new TH1D("h_x_start_fit_vs_z_phi_sector", "h_x_start_fit_vs_z_phi_sector", n_bins_z_phi_sector, 0, n_bins_z_phi_sector); + //---------------------------------------------------------------- + + //-------------------------------------------------------------------------- + // Prepare output data + std::vector>>>> vec_DXYZ_vox; + std::vector>>>> vec_DXYZ_vox_GF; + vec_DXYZ_vox.resize(6); + vec_DXYZ_vox_GF.resize(6); + for (Int_t i_xyz = 0; i_xyz < 6; i_xyz++) { + vec_DXYZ_vox[i_xyz].resize(36); + vec_DXYZ_vox_GF[i_xyz].resize(36); + for (Int_t i_sector = 0; i_sector < 36; i_sector++) { + vec_DXYZ_vox[i_xyz][i_sector].resize(152); + vec_DXYZ_vox_GF[i_xyz][i_sector].resize(152); + for (Int_t voxX = 0; voxX < 152; voxX++) { + vec_DXYZ_vox[i_xyz][i_sector][voxX].resize((nY2XBins)); + vec_DXYZ_vox_GF[i_xyz][i_sector][voxX].resize((nY2XBins)); + for (Int_t voxY = 0; voxY < (nY2XBins); voxY++) { + vec_DXYZ_vox[i_xyz][i_sector][voxX][voxY].resize((nZ2XBins)); + vec_DXYZ_vox_GF[i_xyz][i_sector][voxX][voxY].resize((nZ2XBins)); + for (Int_t voxZ = 0; voxZ < (nZ2XBins); voxZ++) { + vec_DXYZ_vox[i_xyz][i_sector][voxX][voxY][voxZ] = 0.0; + vec_DXYZ_vox_GF[i_xyz][i_sector][voxX][voxY][voxZ] = 0.0; + } + } + } + } + } + std::vector>> vec_max_X_fit; + vec_max_X_fit.resize(36); + for (Int_t i_sector = 0; i_sector < 36; i_sector++) { + vec_max_X_fit[i_sector].resize(nY2XBins); + for (Int_t voxY = 0; voxY < (nY2XBins); voxY++) { + vec_max_X_fit[i_sector][voxY].resize(nZ2XBins); + for (Int_t voxZ = 0; voxZ < (nZ2XBins); voxZ++) { + vec_max_X_fit[i_sector][voxY][voxZ] = 0.0; + } + } + } + + for (Long64_t jentry = 0; jentry < entries_input_map; jentry++) { + voxResTree->GetEntry(jentry); + + const auto bvox_X = voxRes_map->bvox[o2::tpc::TrackResiduals::VoxX]; // bin number in x (= pad row) + const auto bvox_F = voxRes_map->bvox[o2::tpc::TrackResiduals::VoxF]; // bin number in y/x 0..14 + const auto bvox_Z = voxRes_map->bvox[o2::tpc::TrackResiduals::VoxZ]; // bin number in z/x 0..4 + const int sector = (int)voxRes_map->bsec; + const float xAV = voxRes_map->stat[o2::tpc::TrackResiduals::VoxX]; + const float z2xAV = voxRes_map->stat[o2::tpc::TrackResiduals::VoxZ]; + const float zAV = z2xAV * xAV; + + vec_DXYZ_vox[0][sector][bvox_X][bvox_F][bvox_Z] = voxRes_map->D[0]; // dX + vec_DXYZ_vox[1][sector][bvox_X][bvox_F][bvox_Z] = voxRes_map->D[1]; // dY + vec_DXYZ_vox[2][sector][bvox_X][bvox_F][bvox_Z] = voxRes_map->D[2]; // dZ + vec_DXYZ_vox[3][sector][bvox_X][bvox_F][bvox_Z] = voxRes_map->stat[3]; // #entries + vec_DXYZ_vox[4][sector][bvox_X][bvox_F][bvox_Z] = xAV; // xAV + vec_DXYZ_vox[5][sector][bvox_X][bvox_F][bvox_Z] = zAV; // zAV + + vec_DXYZ_vox_GF[0][sector][bvox_X][bvox_F][bvox_Z] = voxRes_map->DS[0]; // dXS + vec_DXYZ_vox_GF[1][sector][bvox_X][bvox_F][bvox_Z] = voxRes_map->DS[1]; // dYS + vec_DXYZ_vox_GF[2][sector][bvox_X][bvox_F][bvox_Z] = voxRes_map->DS[2]; // dZS + vec_DXYZ_vox_GF[3][sector][bvox_X][bvox_F][bvox_Z] = voxRes_map->stat[3]; // #entries + vec_DXYZ_vox_GF[4][sector][bvox_X][bvox_F][bvox_Z] = xAV; // xAV + vec_DXYZ_vox_GF[5][sector][bvox_X][bvox_F][bvox_Z] = zAV; // zAV + + tp_Stat_vs_row->Fill(bvox_X, voxRes_map->stat[3]); + float voxX_pos = RowX[bvox_X]; + if (fabs(bvox_F - (int)(nY2XBins / 2)) <= 1) { + tp_Stat_vs_X->Fill(voxX_pos, voxRes_map->stat[3]); + } + if (bvox_Z == 0) { + // if(voxX_pos < (85.0+32.0)) + { + tp_DX_vs_X_raw->Fill(voxX_pos, voxRes_map->D[0]); + tp_DX_vs_Row_raw->Fill(bvox_X, voxRes_map->D[0]); + if (sector < 18) { + tp_DZ_vs_X_raw_A->Fill(voxX_pos, voxRes_map->D[2]); + } else { + tp_DZ_vs_X_raw_C->Fill(voxX_pos, voxRes_map->D[2]); + } + } + } + } + + //-------------------------------------------------------------------------------- + tp_Stat_vs_row->GetXaxis()->SetRangeUser(20, 62); + const float meanEntries = tp_Stat_vs_row->GetMean(2); + tp_Stat_vs_row->GetXaxis()->SetRangeUser(2, 1); + const int startRowGoodEntries = tp_Stat_vs_row->FindFirstBinAbove(meanEntries * 0.7) - 1; + + // don't trust bins with too low statistics + tp_DX_vs_X_raw->GetXaxis()->SetRangeUser(RowX[startRowGoodEntries], RowX[63]); + // const float max_DX = tp_DX_vs_X_raw ->GetBinContent(tp_DX_vs_X_raw->GetMaximumBin()); + // const float max_X_DX = tp_DX_vs_X_raw ->GetBinCenter(tp_DX_vs_X_raw->GetMaximumBin()); + tp_DX_vs_Row_raw->GetXaxis()->SetRangeUser(startRowGoodEntries, 63); + // const int max_Row_DX = tp_DX_vs_Row_raw ->GetBinCenter(tp_DX_vs_Row_raw->GetMaximumBin()); + + float max_X_stat = 0.0; + for (int ibin = (tp_Stat_vs_X->GetNbinsX() - 3); ibin >= 0; ibin--) { + float X_val = tp_Stat_vs_X->GetBinCenter(ibin); + float stat = tp_Stat_vs_X->GetBinContent(ibin); + float DX = tp_DX_vs_X_raw->GetBinContent(ibin); + if (X_val < (85.0 + 32.0)) { + Double_t stat_previous[3] = {tp_Stat_vs_X->GetBinContent(ibin + 1), tp_Stat_vs_X->GetBinContent(ibin + 2), tp_Stat_vs_X->GetBinContent(ibin + 3)}; + Double_t Xpos_previous[3] = {tp_Stat_vs_X->GetBinCenter(ibin + 1), tp_Stat_vs_X->GetBinCenter(ibin + 2), tp_Stat_vs_X->GetBinCenter(ibin + 3)}; + if ((stat - stat_previous[0]) < 0.0 && (stat - stat_previous[1]) < 0.0 && (stat - stat_previous[2]) < 0.0 && stat > 0.0) { + Double_t ratio_stat[3] = {stat / stat_previous[0], stat / stat_previous[1], stat / stat_previous[2]}; + if (ratio_stat[0] < 0.8 && ratio_stat[1] < 0.5 && ratio_stat[2] < 0.5) { + max_X_stat = Xpos_previous[2]; + break; + } + } + if (fabs(stat < 0.1)) { + max_X_stat = Xpos_previous[2]; + break; + } + } + } + + const float max_X_DX = max_X_stat; + int max_Row_DX = 0; + + // find corresponding row + for (int i = 0; i < 50; ++i) { + if (RowX[i] > max_X_DX) { + break; + } + max_Row_DX = i; + } + const float max_DX = tp_DX_vs_X_raw->GetBinContent(tp_DX_vs_X_raw->FindBin(max_X_DX)); + + LOGP(info, "max DX value: {:.3f}, X-position of max DX value: {:.3f} (row: {}), first row checked: {}, mean entries: {:.2f}, max_X_stat: {:.3f}", max_DX, max_X_DX, max_Row_DX, startRowGoodEntries, meanEntries, max_X_stat); + //-------------------------------------------------------------------------------- + + //-------------------------------------------------------------------------------- + LOGP(info, "Calculating extrapolation fit start values for every phi slice"); + for (Int_t i_sector = 0; i_sector < 36; i_sector++) { + for (Int_t voxY = 0; voxY < (nY2XBins); voxY++) { + for (Int_t voxZ = 0; voxZ < (nZ2XBins); voxZ++) { + // fill the TGraph used for fitting + int ipoint = 0; + for (Int_t voxX = 0; voxX <= 151; voxX++) { + float voxX_pos = RowX[voxX]; + float statistics = vec_DXYZ_vox_GF[3][i_sector][voxX][voxY][voxZ]; + float DX = vec_DXYZ_vox_GF[0][i_sector][voxX][voxY][voxZ]; + tg_Stat_vs_X_slice->SetPoint(ipoint, voxX_pos, statistics); + ipoint++; + + if (i_sector == 11 && voxY == 10 && voxZ == 27) { + tp_Stat_vs_X_single->Fill(voxX_pos, statistics); + tp_DX_vs_X_single->Fill(voxX_pos, DX); + } + } + + float max_X_stat = 0.0; + for (int ibin = (tg_Stat_vs_X_slice->GetN() - 3); ibin >= 0; ibin--) { + double X_val = 0.0; + double stat = 0.0; + tg_Stat_vs_X_slice->GetPoint(ibin, X_val, stat); + if (stat < min_statistics) + continue; + if (X_val < (85.0 + 32.0)) { + Double_t stat_previous[3] = {0.0, 0.0, 0.0}; + Double_t Xpos_previous[3] = {0.0, 0.0, 0.0}; + + tg_Stat_vs_X_slice->GetPoint(ibin + 1, Xpos_previous[0], stat_previous[0]); + tg_Stat_vs_X_slice->GetPoint(ibin + 2, Xpos_previous[1], stat_previous[1]); + tg_Stat_vs_X_slice->GetPoint(ibin + 3, Xpos_previous[2], stat_previous[2]); + + if ((stat - stat_previous[0]) < 0.0 && (stat - stat_previous[1]) < 0.0 && (stat - stat_previous[2]) < 0.0 && stat > 0.0) { + Double_t ratio_stat[3] = {stat / stat_previous[0], stat / stat_previous[1], stat / stat_previous[2]}; + // if(i_sector == 9 && voxY == 19 && voxZ == 27) + //{ + // } + if (ratio_stat[0] < 0.8 && ratio_stat[1] < 0.5 && ratio_stat[2] < 0.5) { + // first check if there aren't any bins at lower radii with statistics + int flag_low_bin = 0; + for (int ibinB = ibin - 1; ibinB >= 0; ibinB--) { + double X_valB = 0.0; + double statB = 0.0; + tg_Stat_vs_X_slice->GetPoint(ibinB, X_valB, statB); + if (statB > 0.0 && statB / stat_previous[0] > 0.8) { + ibin = ibinB; + flag_low_bin = 1; + break; + } + } + if (!flag_low_bin) { + max_X_stat = Xpos_previous[2]; + break; + } + } + } + if (fabs(stat < 0.1)) { + // first check if there aren't any bins at lower radii with statistics + int flag_low_bin = 0; + for (int ibinB = ibin - 1; ibinB >= 0; ibinB--) { + double X_valB = 0.0; + double statB = 0.0; + tg_Stat_vs_X_slice->GetPoint(ibinB, X_valB, statB); + if (statB > 0.0) { + ibin = ibinB; + flag_low_bin = 1; + break; + } + } + if (!flag_low_bin) { + max_X_stat = Xpos_previous[2]; + break; + } + } + } + } + + if (max_X_stat < 85.0) { + max_X_stat = max_X_DX; // set to average value + } + int i_bin_z_phi_sector = voxZ * 36 * nY2XBins + i_sector * nY2XBins + voxY; + h_x_start_fit_vs_z_phi_sector->SetBinContent(i_bin_z_phi_sector, max_X_stat); + tg_Stat_vs_X_slice->Set(0); + + vec_max_X_fit[i_sector][voxY][voxZ] = max_X_stat; + } // end of Z loop + } // end of Y loop + } // end of sector loop + LOGP(info, "Done calculating extrapolation fit start values for every phi slice"); + //-------------------------------------------------------------------------------- + + // ========================================================================= + // treat acceptance edge in z-direction + // replace values very close to or beyond the pad plane by a value a bit further away + const float maxZ = 242.f; + for (Int_t i_sector = 0; i_sector < 36; i_sector++) { + for (Int_t voxX = 0; voxX < nXBins; voxX++) { + for (Int_t voxY = 0; voxY < nY2XBins; voxY++) { + float lastValue[4] = {0.f, 0.f, 0.f, 0.f}; + float lastValueS[4] = {0.f, 0.f, 0.f, 0.f}; + float lastValueA11[4] = {0.f, 0.f, 0.f, 0.f}; + float lastValueSA11[4] = {0.f, 0.f, 0.f, 0.f}; + for (Int_t voxZ = 0; voxZ < nZ2XBins; voxZ++) { + const float absZ = std::abs(vec_DXYZ_vox[5][i_sector][voxX][voxY][voxZ]); + // if (i_sector==0&&voxX>149&&voxY==5) { + //} + if (absZ < maxZ) { + for (int i = 0; i < 4; ++i) { + lastValue[i] = vec_DXYZ_vox[i][i_sector][voxX][voxY][voxZ]; + lastValueS[i] = vec_DXYZ_vox_GF[i][i_sector][voxX][voxY][voxZ]; + } + } else { + for (int i = 0; i < 4; ++i) { + vec_DXYZ_vox[i][i_sector][voxX][voxY][voxZ] = lastValue[i]; + vec_DXYZ_vox_GF[i][i_sector][voxX][voxY][voxZ] = lastValueS[i]; + } + } + + if (i_sector == 11 && A11maxZ2X > -1) { + if (voxZ <= A11maxZ2X) { + for (int i = 0; i < 4; ++i) { + lastValueA11[i] = vec_DXYZ_vox[i][i_sector][voxX][voxY][voxZ]; + lastValueSA11[i] = + vec_DXYZ_vox_GF[i][i_sector][voxX][voxY][voxZ]; + } + } else { + for (int i = 0; i < 4; ++i) { + vec_DXYZ_vox[i][i_sector][voxX][voxY][voxZ] = lastValueA11[i]; + vec_DXYZ_vox_GF[i][i_sector][voxX][voxY][voxZ] = + lastValueSA11[i]; + } + } + } + } + } + } + } + + //---------------------------------------------------------------- + // Gaussian filtering + + if (do_smoothing) { + LOGP(info, "Gaussian filtering started"); + std::vector>> vec_GKernel; + vec_GKernel.resize(N_bins_X_GF * 2 + 1); + for (Int_t i_X = 0; i_X < (Int_t)vec_GKernel.size(); i_X++) { + vec_GKernel[i_X].resize(N_bins_Y_GF * 2 + 1); + for (Int_t i_Y = 0; i_Y < (Int_t)vec_GKernel[i_X].size(); i_Y++) { + vec_GKernel[i_X][i_Y].resize(N_bins_Z_GF * 2 + 1); + } + } + vec_FilterCreation(vec_GKernel, N_bins_X_GF, N_bins_Y_GF, N_bins_Z_GF, sigma_GF); + + for (Int_t i_sector = 0; i_sector < 36; i_sector++) { + std::vector>>> arr_values; + std::vector>>> arr_values_used; + arr_values.resize(3); + arr_values_used.resize(3); + + for (Int_t i_xyz = 0; i_xyz < 3; i_xyz++) { + arr_values[i_xyz].resize(N_bins_X_GF * 2 + 1); + arr_values_used[i_xyz].resize(N_bins_X_GF * 2 + 1); + for (Int_t i_X = 0; i_X < (Int_t)arr_values[i_xyz].size(); i_X++) { + arr_values[i_xyz][i_X].resize(N_bins_Y_GF * 2 + 1); + arr_values_used[i_xyz][i_X].resize(N_bins_Y_GF * 2 + 1); + for (Int_t i_Y = 0; i_Y < (Int_t)arr_values[i_xyz][i_X].size(); i_Y++) { + arr_values[i_xyz][i_X][i_Y].resize(N_bins_Z_GF * 2 + 1); + arr_values_used[i_xyz][i_X][i_Y].resize(N_bins_Z_GF * 2 + 1); + for (Int_t i_Z = 0; i_Z < (Int_t)arr_values[i_xyz][i_X][i_Y].size(); i_Z++) { + arr_values[i_xyz][i_X][i_Y][i_Z] = 0.0; + arr_values_used[i_xyz][i_X][i_Y][i_Z] = 0.0; + } + } + } + } + + // CRU 0 : 000 - 016 (IROC) + // CRU 1 : 017 - 031 (IROC) + // CRU 2 : 032 - 047 (IROC) + // CRU 3 : 048 - 062 (IROC) + + // CRU 4 : 063 - 080 (OROC 1) + // CRU 5 : 081 - 096 (OROC 1) + + // CRU 6 : 097 - 112 (OROC 2) + // CRU 7 : 113 - 126 (OROC 2) + + // CRU 8 : 127 - 139 (OROC 3) + // CRU 9 : 140 - 151 (OROC 3) + + std::vector> vec_ROC_row; + vec_ROC_row.resize(4); + for (int iRoc = 0; iRoc < 4; iRoc++) { + vec_ROC_row[iRoc].resize(2); + } + vec_ROC_row[0][0] = 0; + vec_ROC_row[0][1] = 62; + vec_ROC_row[1][0] = 63; + vec_ROC_row[1][1] = 96; + vec_ROC_row[2][0] = 97; + vec_ROC_row[2][1] = 126; + vec_ROC_row[3][0] = 127; + vec_ROC_row[3][1] = 151; + + for (int iRoc = 0; iRoc < 4; iRoc++) { + int minRow = vec_ROC_row[iRoc][0]; + if (do_extrapolation && (iRoc == 0)) { + minRow = max_Row_DX + 2; // don't smooth over low radii large distortions + } + for (Int_t voxX = minRow; voxX <= vec_ROC_row[iRoc][1]; voxX++) { + for (Int_t voxY = 0; voxY < (nY2XBins); voxY++) { + for (Int_t voxZ = 0; voxZ < (nZ2XBins); voxZ++) { + Float_t sum_weight[3] = {0.0}; + Float_t sum_values[3] = {0.0}; + for (Int_t index_voxXB = -N_bins_X_GF; index_voxXB <= N_bins_X_GF; index_voxXB++) { + Int_t voxXB = voxX + index_voxXB; + if (voxXB < vec_ROC_row[iRoc][0]) + continue; + if (voxXB > vec_ROC_row[iRoc][1]) + continue; + for (Int_t index_voxYB = -N_bins_Y_GF; index_voxYB <= N_bins_Y_GF; index_voxYB++) { + Int_t voxYB = voxY + index_voxYB; + if (voxYB < 0) + continue; + if (voxYB >= nY2XBins) + continue; + for (Int_t index_voxZB = -N_bins_Z_GF; index_voxZB <= N_bins_Z_GF; index_voxZB++) { + Int_t voxZB = voxZ + index_voxZB; + if (voxZB < 0) + continue; + if (voxZB >= (nZ2XBins)) + continue; + float statistics = vec_DXYZ_vox[3][i_sector][voxXB][voxYB][voxZB]; + if ((int)statistics == 0) + continue; + if (TMath::IsNaN(vec_DXYZ_vox[0][i_sector][voxXB][voxYB][voxZB])) + continue; // NaN check + if (TMath::IsNaN(vec_DXYZ_vox[1][i_sector][voxXB][voxYB][voxZB])) + continue; // NaN check + if (TMath::IsNaN(vec_DXYZ_vox[2][i_sector][voxXB][voxYB][voxZB])) + continue; // NaN check + if (fabs(vec_DXYZ_vox[0][i_sector][voxXB][voxYB][voxZB]) > maxDeltaCut) + continue; + if (fabs(vec_DXYZ_vox[1][i_sector][voxXB][voxYB][voxZB]) > maxDeltaCut) + continue; + if (fabs(vec_DXYZ_vox[2][i_sector][voxXB][voxYB][voxZB]) > maxDeltaCut) + continue; + for (Int_t i_xyz = 0; i_xyz < 3; i_xyz++) { + arr_values_used[i_xyz][index_voxXB + N_bins_X_GF][index_voxYB + N_bins_Y_GF][index_voxZB + N_bins_Z_GF] = 1.0; + arr_values[i_xyz][index_voxXB + N_bins_X_GF][index_voxYB + N_bins_Y_GF][index_voxZB + N_bins_Z_GF] = vec_GKernel[index_voxXB + N_bins_X_GF][index_voxYB + N_bins_Y_GF][index_voxZB + N_bins_Z_GF] * vec_DXYZ_vox[i_xyz][i_sector][voxXB][voxYB][voxZB]; + sum_weight[i_xyz] += vec_GKernel[index_voxXB + N_bins_X_GF][index_voxYB + N_bins_Y_GF][index_voxZB + N_bins_Z_GF]; + sum_values[i_xyz] += arr_values[i_xyz][index_voxXB + N_bins_X_GF][index_voxYB + N_bins_Y_GF][index_voxZB + N_bins_Z_GF]; + if (TMath::IsNaN(vec_DXYZ_vox[i_xyz][i_sector][voxXB][voxYB][voxZB])) { + LOGP(error, "NaN vec_DXYZ_vox detected in xyz {}, sec {}, voxX {}, voxY {}, voxZ {}", i_xyz, i_sector, voxXB, voxYB, voxZB); + } + if (TMath::IsNaN(sum_values[i_xyz])) { + LOGP(error, "NaN sum_values detected in xyz {}, sec {}, voxX {}, voxY {}, voxZ {}", i_xyz, i_sector, voxXB, voxYB, voxZB); + } + } + } + } + } + + for (Int_t i_xyz = 0; i_xyz < 3; i_xyz++) { + if (sum_weight[i_xyz] > 0.0) { + sum_values[i_xyz] /= sum_weight[i_xyz]; + vec_DXYZ_vox_GF[i_xyz][i_sector][voxX][voxY][voxZ] = sum_values[i_xyz]; + if (TMath::IsNaN(vec_DXYZ_vox_GF[i_xyz][i_sector][voxX][voxY][voxZ])) { + LOGP(error, "NaN detected during smoothing process in xyz {}, sec {}, voxX {}, voxY {}, voxZ {}", i_xyz, i_sector, voxX, voxY, voxZ); + } + + float voxX_pos = RowX[voxX]; + if (voxZ == 0) { + if (i_xyz == 0) + tp_DX_vs_X_smooth->Fill(voxX_pos, sum_values[i_xyz]); + if (i_sector < 18) { + if (i_xyz == 2) + tp_DZ_vs_X_smooth_A->Fill(voxX_pos, sum_values[i_xyz]); + } else { + if (i_xyz == 2) + tp_DZ_vs_X_smooth_C->Fill(voxX_pos, sum_values[i_xyz]); + } + } + } + } + } + } + } + } + } // end loop over sectors + } // do_smoothing + //---------------------------------------------------------------- + + //---------------------------------------------------------------- + // Do extrapolation to low radii + if (do_extrapolation > 0) { + LOGP(info, "Starting extrapolation to small radii"); + // const float start_fit = max_X_DX + 0.0; + // const float stop_fit = max_X_DX + 10.0; + TGraph* tg_data_for_fit = new TGraph(); + TGraph* tg_data_for_fit_raw = new TGraph(); + for (Int_t i_sector = 0; i_sector < 36; i_sector++) { + for (Int_t voxY = 0; voxY < (nY2XBins); voxY++) { + for (Int_t voxZ = 0; voxZ < (nZ2XBins); voxZ++) { + const float start_fit = vec_max_X_fit[i_sector][voxY][voxZ] + 0.0; + const float stop_fit = vec_max_X_fit[i_sector][voxY][voxZ] + 10.0; + // find maximum dX + // float maxXdX = 0; + // float maxdX = 0; + // for(Int_t voxX = startRowGoodEntries; voxX < 63; voxX++) // enough to search in IROC + //{ + // const float dX = vec_DXYZ_vox[0][i_sector][voxX][voxY][voxZ]; + // if (dX > maxdX) { + // maxdX = dX; + // maxXdX = RowX[voxX]; + //} + //} + // const float start_fit = maxXdX + 1.0; + // const float stop_fit = maxXdX + 10.0; + + for (Int_t i_xyz = 0; i_xyz < 3; i_xyz++) { + // fill the TGraph used for fitting + int ipoint = 0; + int meanStat = 1; // statistics to use in extrapolation region + for (Int_t voxX = 0; voxX <= 151; voxX++) { + // float voxX_pos = vec_DXYZ_vox_GF[4][i_sector][voxX][voxY][voxZ]; + float voxX_pos = RowX[voxX]; + if (voxX_pos < start_fit) + continue; + if (voxX_pos > stop_fit) + break; + float statistics = vec_DXYZ_vox_GF[3][i_sector][voxX][voxY][voxZ]; + if (statistics < min_statistics) + continue; + float DXYZval = vec_DXYZ_vox_GF[i_xyz][i_sector][voxX][voxY][voxZ]; + tg_data_for_fit->SetPoint(ipoint, voxX_pos, DXYZval); + float DXYZval_raw = vec_DXYZ_vox[i_xyz][i_sector][voxX][voxY][voxZ]; + tg_data_for_fit_raw->SetPoint(ipoint, voxX_pos, DXYZval_raw); + ipoint++; + } + if (ipoint < 5) + continue; + + // fit the TGraph + for (Int_t i = 0; i < 6; i++) { + func_PolyFitFunc->SetParameter(i, 0.0); + func_PolyFitFunc->SetParError(i, 0.0); + func_PolyFitFunc_raw->SetParameter(i, 0.0); + func_PolyFitFunc_raw->SetParError(i, 0.0); + if (i > 2) { + func_PolyFitFunc->FixParameter(i, 0.0); + func_PolyFitFunc_raw->FixParameter(i, 0.0); + } + } + func_PolyFitFunc->SetParameter(0, 0.2); + func_PolyFitFunc->SetParameter(1, 0.3); + func_PolyFitFunc->SetParameter(2, 0.4); + func_PolyFitFunc->SetRange(start_fit, stop_fit); + tg_data_for_fit->Fit("func_PolyFitFunc", "QWMN", "", start_fit, stop_fit); + + func_PolyFitFunc_raw->SetParameter(0, 0.2); + func_PolyFitFunc_raw->SetParameter(1, 0.3); + func_PolyFitFunc_raw->SetParameter(2, 0.4); + func_PolyFitFunc_raw->SetRange(start_fit, stop_fit); + tg_data_for_fit_raw->Fit("func_PolyFitFunc_raw", "QWMN", "", start_fit, stop_fit); + + // if (ipoint>0) { + // meanStat /= ipoint; + // } + + // do the low radii extrapolation + for (Int_t voxX = 0; voxX <= 151; voxX++) { + // float voxX_pos = vec_DXYZ_vox_GF[4][i_sector][voxX][voxY][voxZ]; + float voxX_pos = RowX[voxX]; + if (voxX_pos > start_fit) { + break; + } + double extrapolation_value = func_PolyFitFunc->Eval(voxX_pos); + if (fabs(extrapolation_value) > max_extrapolation_value) + extrapolation_value = TMath::Sign(1, extrapolation_value) * max_extrapolation_value; + vec_DXYZ_vox_GF[i_xyz][i_sector][voxX][voxY][voxZ] = extrapolation_value; + vec_DXYZ_vox_GF[3][i_sector][voxX][voxY][voxZ] = meanStat; // set statistics -> important for spline creation + + double extrapolation_value_raw = func_PolyFitFunc_raw->Eval(voxX_pos); + if (fabs(extrapolation_value_raw) > max_extrapolation_value) + extrapolation_value_raw = TMath::Sign(1, extrapolation_value_raw) * max_extrapolation_value; + vec_DXYZ_vox[i_xyz][i_sector][voxX][voxY][voxZ] = extrapolation_value_raw; + vec_DXYZ_vox[3][i_sector][voxX][voxY][voxZ] = meanStat; // set statistics -> important for spline creation + } + + tg_data_for_fit->Set(0); + tg_data_for_fit_raw->Set(0); + } // end of xyz + + for (Int_t voxX = 0; voxX <= 151; voxX++) { + if (voxZ == 0) // for QA + { + float voxX_pos = RowX[voxX]; + tp_DX_vs_X_smooth_extr->Fill(voxX_pos, vec_DXYZ_vox_GF[0][i_sector][voxX][voxY][voxZ]); + if (i_sector < 18) { + tp_DZ_vs_X_smooth_extr_A->Fill(voxX_pos, vec_DXYZ_vox_GF[2][i_sector][voxX][voxY][voxZ]); + } else { + tp_DZ_vs_X_smooth_extr_C->Fill(voxX_pos, vec_DXYZ_vox_GF[2][i_sector][voxX][voxY][voxZ]); + } + } + } + } + } + } + } + //---------------------------------------------------------------- + + //---------------------------------------------------------------- + // IROC A11 maskign + if (maskIA11) { + int i_sector = 11; + for (Int_t voxY = 0; voxY < (nY2XBins); voxY++) { + for (Int_t voxZ = 0; voxZ < (nZ2XBins); voxZ++) { + for (Int_t voxX = 0; voxX <= 62; voxX++) { + for (Int_t i_xyz = 0; i_xyz < 4; i_xyz++) { + vec_DXYZ_vox_GF[i_xyz][i_sector][voxX][voxY][voxZ] = 0; + vec_DXYZ_vox[i_xyz][i_sector][voxX][voxY][voxZ] = 0; + } + } + } + } + } + + //---------------------------------------------------------------- + mTreeOut = std::make_unique("voxResTree", "Voxel results and statistics"); + mTreeOut->Branch("voxRes", &mVoxelResultsOutPtr); + // copy user info + auto userInfoOut = mTreeOut->GetUserInfo(); + for (auto o : *userInfo) { + userInfoOut->Add(o->Clone()); + } + // Overwrite the placeholder meanIDC/medianIDC just cloned above with the real offline-joined + // values computed in the "Offline IDC join" block earlier in this function. + if (auto* stale = userInfoOut->FindObject("meanIDC")) { + userInfoOut->Remove(stale); + delete stale; + } + if (auto* stale = userInfoOut->FindObject("medianIDC")) { + userInfoOut->Remove(stale); + delete stale; + } + userInfoOut->Add(new TNamed("meanIDC", std::to_string(meanIDCReal).data())); + userInfoOut->Add(new TNamed("medianIDC", std::to_string(medianIDCReal).data())); + userInfoOut->Add(new TNamed("startRowGoodEntries", std::to_string(startRowGoodEntries).data())); + userInfoOut->Add(new TNamed("maxDX", std::to_string(max_DX).data())); + userInfoOut->Add(new TNamed("maxDX_lx", std::to_string(max_X_DX).data())); + userInfoOut->Add(new TNamed("maxDX_row", std::to_string(max_Row_DX).data())); + + // copy aliases + if (voxResTree->GetListOfAliases()) { + for (auto o : *voxResTree->GetListOfAliases()) { + mTreeOut->SetAlias(o->GetName(), o->GetTitle()); + } + } + + for (Long64_t jentry = 0; jentry < entries_input_map; jentry++) { + voxResTree->GetEntry(jentry); + + auto bvox_X = voxRes_map->bvox[o2::tpc::TrackResiduals::VoxX]; // bin number in x (= pad row) + auto bvox_F = voxRes_map->bvox[o2::tpc::TrackResiduals::VoxF]; // bin number in y/x 0..14 + auto bvox_Z = voxRes_map->bvox[o2::tpc::TrackResiduals::VoxZ]; // bin number in z/x 0..4 + int sector = (int)voxRes_map->bsec; + // Int_t index_map = bvox_X+152*bvox_F+152*(nY2XBins)*bvox_Z+152*(nY2XBins)*(nZ2XBins)*sector; + + mVoxelResultsOut = *voxRes_map; // copy the entry from the input map + + for (int ixyz = 0; ixyz < 3; ixyz++) { + if (do_extrapolation == 1 || do_extrapolation == 2) + mVoxelResultsOut.DS[ixyz] = (float)vec_DXYZ_vox_GF[ixyz][sector][bvox_X][bvox_F][bvox_Z]; // overwrite the smoothed values + if (do_extrapolation == 2 || do_extrapolation == 3) + mVoxelResultsOut.D[ixyz] = (float)vec_DXYZ_vox[ixyz][sector][bvox_X][bvox_F][bvox_Z]; // overwrite the raw values + } + if (do_extrapolation == 1 || do_extrapolation == 2) + mVoxelResultsOut.stat[3] = (float)vec_DXYZ_vox_GF[3][sector][bvox_X][bvox_F][bvox_Z]; + if (do_extrapolation == 2 || do_extrapolation == 3) + mVoxelResultsOut.stat[3] = (float)vec_DXYZ_vox[3][sector][bvox_X][bvox_F][bvox_Z]; + for (int ixyz = 0; ixyz < 3; ixyz++) { + if (TMath::IsNaN(mVoxelResultsOut.DS[ixyz])) // NaN + { + LOGP(error, "NaN detected in smoothed value xyz {}, sec {}, voxX {}, voxY {}, voxZ {}", ixyz, sector, bvox_X, bvox_F, bvox_Z); + mVoxelResultsOut.DS[ixyz] = 0.0; + mVoxelResultsOut.stat[3] = 0; + mVoxelResultsOut.flags |= TrackResiduals::Masked; + } else { + mVoxelResultsOut.flags |= TrackResiduals::SmoothDone; + } + if (TMath::IsNaN(mVoxelResultsOut.D[ixyz])) // NaN + { + LOGP(error, "NaN detected in raw value xyz {}, sec {}, voxX {}, voxY {}, voxZ {}", ixyz, sector, bvox_X, bvox_F, bvox_Z); + mVoxelResultsOut.D[ixyz] = 0.0; + } + } + mTreeOut->Fill(); + } + //---------------------------------------------------------------- + + //---------------------------------------------------------------- + outputfile->cd(); + mTreeOut->Write(); + mTreeOut.release(); + tp_DX_vs_X_raw->Write(); + tp_DX_vs_X_smooth->Write(); + tp_DX_vs_X_smooth_extr->Write(); + tp_DZ_vs_X_raw_A->Write(); + tp_DZ_vs_X_raw_C->Write(); + tp_DZ_vs_X_smooth_A->Write(); + tp_DZ_vs_X_smooth_C->Write(); + tp_DZ_vs_X_smooth_extr_A->Write(); + tp_DZ_vs_X_smooth_extr_C->Write(); + tp_Stat_vs_X->Write(); + tp_Stat_vs_row->Write(); + h_x_start_fit_vs_z_phi_sector->Write(); + tp_Stat_vs_X_single->Write(); + tp_DX_vs_X_single->Write(); + outputfile->Close(); + //---------------------------------------------------------------- +} diff --git a/Detectors/TPC/calibration/SpacePoints/macro/staticMapCreatorCPM.C b/Detectors/TPC/calibration/SpacePoints/macro/staticMapCreatorCPM.C new file mode 100644 index 0000000000000..dea2b2f98c2ef --- /dev/null +++ b/Detectors/TPC/calibration/SpacePoints/macro/staticMapCreatorCPM.C @@ -0,0 +1,2564 @@ +// 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. + +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include "CCDB/CcdbApi.h" +#include "CCDB/BasicCCDBManager.h" + +#include "TSystem.h" +#include "Algorithm/RangeTokenizer.h" +#include "Framework/Logger.h" +#include "CommonConstants/LHCConstants.h" +#include "SpacePoints/SpacePointsCalibConfParam.h" +#include "SpacePoints/TrackResiduals.h" +#include "SpacePoints/TrackInterpolation.h" +#include "DataFormatsParameters/GRPMagField.h" +#include "DataFormatsParameters/GRPLHCIFData.h" +#include "DataFormatsTPC/Defs.h" +#include "ReconstructionDataFormats/GlobalTrackID.h" +#include "DetectorsBase/MatLayerCylSet.h" +#include "DetectorsBase/Propagator.h" +#include "TPCBase/Mapper.h" +#include "ReconstructionDataFormats/TrackUtils.h" + +#include +#include +#include +#include +#include "TTreePerfStats.h" +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +// For multiple threads +#include +#include +#include +#include +#include +#include "TROOT.h" + +// Compile-time detection of TrackData::filterFlag / UnbinnedResid::rejected, for back-compat with O2 +// builds that predate them. +template +struct HasFilterFlagMember : std::false_type { +}; +template +struct HasFilterFlagMember().filterFlag)>> : std::true_type { +}; + +template +struct HasRejectedMember : std::false_type { +}; +template +struct HasRejectedMember().rejected)>> : std::true_type { +}; + +template +bool hasPositiveFilterFlag(const T& t) +{ + if constexpr (HasFilterFlagMember::value) { + return t.filterFlag > 0; + } else { + return false; + } +} + +template +bool isRejectedResidual(const T& t) +{ + if constexpr (HasRejectedMember::value) { + return t.rejected; + } else { + return false; + } +} + +#else + +#error This macro must run in compiled mode + +#endif + +using namespace o2::tpc; +using GID = o2::dataformats::GlobalTrackID; +namespace fs = std::filesystem; + +constexpr int NSectors = SECTORSPERSIDE * SIDES; +constexpr int NRows = Mapper::PADROWS; + +// Portable peak-RSS report via getrusage(), which works on both Linux and macOS (unlike parsing +// /proc/self/status, which is Linux-only procfs). ru_maxrss's UNIT differs by platform though (bytes on +// macOS, KB on Linux) -- handled below. Only reports peak resident-set size (one number), not the +// separate VmRSS/VmPeak/VmSize that /proc/self/status exposes on Linux. +void printMemoryUsage(const std::string& label = "") +{ + struct rusage ru; + if (getrusage(RUSAGE_SELF, &ru) != 0) { + return; + } +#ifdef __APPLE__ + const double maxRssGB = ru.ru_maxrss / 1024.0 / 1024.0 / 1024.0; // macOS: ru_maxrss in bytes +#else + const double maxRssGB = ru.ru_maxrss / 1024.0 / 1024.0; // Linux: ru_maxrss in KB +#endif + if (label.empty()) { + LOGP(info, "VmPeakRSS {:.2f} GB", maxRssGB); + } else { + LOGP(info, "[{}] VmPeakRSS {:.2f} GB", label, maxRssGB); + } +} + +template +double calculateMean(const std::vector& vec) +{ + if (vec.empty()) { + LOGP(error, "vector is empty"); + return 0.; + } + + return std::accumulate(vec.begin(), vec.end(), 0.0) / vec.size(); +} + +template +double calculateMedian(std::vector vec) +{ + if (vec.empty()) { + LOGP(error, "vector is empty"); + return 0.; + } + + const size_t size = vec.size(); + const size_t midIndex = size / 2; + + std::nth_element(vec.begin(), vec.begin() + midIndex, vec.end()); + + if (size % 2 != 0) { + return vec[midIndex]; + } + + return (vec[midIndex - 1] + vec[midIndex]) / 2.0; +} + +// Lightweight double-precision 3-vector for the hot per-residual/per-voxel-flush loop, ~10x faster +// than TVector3 here (0.550s -> 0.054s over 20M calls) for numerically identical math -- matters +// because getIntCircles() runs while the per-voxel mutex is held. Double (not float) because this is +// live geometric computation (sqrt- and division-heavy) rather than storage -- contrast Vec3f below, +// which is storage-only and fine in single precision. +struct Vec3d { + double x = 0.0, y = 0.0, z = 0.0; + double X() const { return x; } + double Y() const { return y; } + double Z() const { return z; } + void SetXYZ(double xx, double yy, double zz) + { + x = xx; + y = yy; + z = zz; + } + double Perp() const { return std::sqrt(x * x + y * y); } + double Mag() const { return std::sqrt(x * x + y * y + z * z); } + void RotateZ(double angle) + { + const double c = std::cos(angle), s = std::sin(angle); + const double xn = x * c - y * s, yn = x * s + y * c; + x = xn; + y = yn; + } + Vec3d& operator-=(const Vec3d& o) + { + x -= o.x; + y -= o.y; + z -= o.z; + return *this; + } + Vec3d& operator+=(const Vec3d& o) + { + x += o.x; + y += o.y; + z += o.z; + return *this; + } + Vec3d& operator*=(double a) + { + x *= a; + y *= a; + z *= a; + return *this; + } +}; +inline Vec3d operator-(const Vec3d& a, const Vec3d& b) { return {a.x - b.x, a.y - b.y, a.z - b.z}; } +inline Vec3d operator*(const Vec3d& a, double s) { return {a.x * s, a.y * s, a.z * s}; } + +//------------------------------------------------------------------------------------------------------------ +// Crossing point of two circles in the transverse plane. +// +// Two sentinel returns, both of which callers reject via their own transverse-radius band cut: +// (9999, 9999, 9999) -- the circles do not intersect at all (Perp() far above any accepted radius) +// (0, 0, 0) -- they intersect, but neither solution falls in the valid TPC radial band below +// A third, implicit outcome is NaN: when d lands within a few ulps of a tangency bound (d == r1+r2 or +// d == |r1-r2|) the guards below still pass, but rounding can drive the sqrt argument marginally +// negative, since a == r1 exactly at both bounds in exact arithmetic. NaN then fails the callers' +// band cut too (every comparison against it is false), so such a degenerate pair is simply dropped -- +// which is the wanted behaviour, hence no clamping here. +Vec3d getIntCircles(double r1, double r2, Vec3d circleCenter1, Vec3d circleCenter2, double voxX, double voxY) +{ + Vec3d vecSp; // default-constructed to (0, 0, 0), which is the "no solution accepted" sentinel + + const Vec3d dVec = (circleCenter1 - circleCenter2); + const double d = dVec.Perp(); // dist. circle center 1 & 2 + + if (d > r1 + r2) { + // no solutions, the circles are separate + vecSp.SetXYZ(9999, 9999, 9999); + return vecSp; + } + if (d < std::fabs(r1 - r2)) { + // no solutions, one circle is contained in the other + vecSp.SetXYZ(9999, 9999, 9999); + return vecSp; + } + if (d < 0.0001) { + // no solutions, same circle center + vecSp.SetXYZ(9999, 9999, 9999); + return vecSp; + } + + const double dx = circleCenter2.X() - circleCenter1.X(); + const double dy = circleCenter2.Y() - circleCenter1.Y(); + + const double a = (r1 * r1 - r2 * r2 + d * d) / (2 * d); + const double h = std::sqrt(r1 * r1 - a * a); + + // Midpoint of the chord joining the two intersection points. Deliberately written as a reciprocal + // multiply rather than `a * dx / d`: floating-point reciprocal-then-multiply and multiply-then-divide + // are not guaranteed to give the same last-bit result. Keep this exact form. + const double Mx = circleCenter1.X() + (1 / d) * a * dx; + const double My = circleCenter1.Y() + (1 / d) * a * dy; + + const double sp1x = Mx + h * dy / d; + const double sp2x = Mx - h * dy / d; + + const double sp1y = My - h * dx / d; + const double sp2y = My + h * dx / d; + + // The two solutions take the z of the circle they came from, not a common z. + const double sp1z = circleCenter1.Z(); + const double sp2z = circleCenter2.Z(); + + const double sp1Perp = std::sqrt(sp1x * sp1x + sp1y * sp1y); + const double sp2Perp = std::sqrt(sp2x * sp2x + sp2y * sp2y); + + const bool sp1In = (sp1Perp > 60.0 && sp1Perp < 280.0); + const bool sp2In = (sp2Perp > 60.0 && sp2Perp < 280.0); + + if (sp1In && sp2In) { + // Both solutions are physically plausible -- pick whichever is closer to the voxel center rather + // than an arbitrary, data-independent choice that could just as easily discard the better of two + // valid solutions. + const double d1 = std::sqrt((sp1x - voxX) * (sp1x - voxX) + (sp1y - voxY) * (sp1y - voxY)); + const double d2 = std::sqrt((sp2x - voxX) * (sp2x - voxX) + (sp2y - voxY) * (sp2y - voxY)); + if (d1 <= d2) { + vecSp.SetXYZ(sp1x, sp1y, sp1z); + } else { + vecSp.SetXYZ(sp2x, sp2y, sp2z); + } + } else if (sp1In) { + vecSp.SetXYZ(sp1x, sp1y, sp1z); + } else if (sp2In) { + vecSp.SetXYZ(sp2x, sp2y, sp2z); + } + // else: neither in band, vecSp stays at its default (0,0,0) sentinel + + return vecSp; +} +//------------------------------------------------------------------------------------------------------------ + +struct range { + long from{-1}; + long to{-1}; + + bool operator<(const range& other) + { + return from < other.from; + } + + void sort() + { + if (from > to) { + std::swap(from, to); + } + } +}; + +// ---- Fastest-replica selection for alien:// residual files (used by getInputFileList below) ---- +// Different Storage Elements serving the same LFN can have very different real-world throughput +// depending on where this job actually lands on the network (verified: ALICE::FZK::SE faster than +// ALICE::CERN::EOS from one machine, the reverse from another) -- a name-based heuristic can't predict +// that, so probe with a real timed read instead. Residual files are O(10GB), so getting this wrong once +// costs far more than the probe itself. +struct SEReplica { + std::string se; + std::string pfn; +}; + +std::vector getAlienReplicas(const std::string& plainLFN) +{ + std::vector result; + std::string cmd = "alien.py whereis -r " + plainLFN + " 2>/dev/null"; + std::unique_ptr pipe(popen(cmd.c_str(), "r"), pclose); + if (!pipe) { + return result; + } + std::array buffer; + std::string output; + while (fgets(buffer.data(), buffer.size(), pipe.get()) != nullptr) { + output += buffer.data(); + } + + auto trim = [](std::string& s) { + s.erase(0, s.find_first_not_of(" \t")); + s.erase(s.find_last_not_of(" \t\r\n") + 1); + }; + std::istringstream iss(output); + std::string line; + while (std::getline(iss, line)) { + auto sePos = line.find("SE =>"); + auto pfnPos = line.find("pfn =>"); + if (sePos == std::string::npos || pfnPos == std::string::npos) { + continue; + } + std::string se = line.substr(sePos + 5, pfnPos - sePos - 5); + std::string pfn = line.substr(pfnPos + 6); + trim(se); + trim(pfn); + result.push_back({se, pfn}); + } + return result; +} + +// Generates a tiny standalone probe macro on disk (ROOT requires the file stem to match the top-level +// function name) that opens one alien:// replica and reads a few real entries from the 'unbinnedResid' +// tree (the same tree/branch doFileProcessing itself reads, with the same two dominant unused +// sub-branches disabled -- see there), timing the real transfer via TFile::GetBytesRead(). Prints +// "PROBE_OK " or "PROBE_FAIL" to stdout -- this is invoked as a subprocess by +// probeOneReplicaSubprocess below, never called in-process. +std::string writeProbeChildMacro(const std::string& funcName) +{ + static const char* templateSrc = + "#include \n" + "#include \n" + "#include \n" + "#include \n" + "#include \"SpacePoints/TrackInterpolation.h\"\n" + "#include \n" + "#include \n" + "#include \n" + "#include \n" + "#include \n" + "\n" + "void @FUNC@(const char* plainLFN, const char* se, Long64_t probeEntries = 3)\n" + "{\n" + " // This runs in a fresh, separate ROOT process (see probeOneReplicaSubprocess) -- unlike the parent\n" + " // process, gGrid is never already set up here, so it must be connected explicitly.\n" + " if (!gGrid && !TGrid::Connect(\"alien://\")) {\n" + " printf(\"PROBE_FAIL\\n\");\n" + " return;\n" + " }\n" + " std::string url = std::string(\"alien://\") + plainLFN + \"?se=\" + se;\n" + " std::unique_ptr f(TFile::Open(url.c_str(), \"READ\"));\n" + " if (!f || f->IsZombie()) { printf(\"PROBE_FAIL\\n\"); return; }\n" + " TTreeReader reader(\"unbinnedResid\", f.get());\n" + " TTree* residTree = reader.GetTree();\n" + " if (!residTree) { printf(\"PROBE_FAIL\\n\"); return; }\n" + " TTreeReaderValue> res(reader, \"res\");\n" + " residTree->SetBranchStatus(\"res.tgSlp\", 0);\n" + " residTree->SetBranchStatus(\"res.channel\", 0);\n" + "\n" + " Long64_t bytesBefore = f->GetBytesRead();\n" + " auto t0 = std::chrono::steady_clock::now();\n" + " Long64_t nRead = 0;\n" + " for (Long64_t i = 0; i < probeEntries && reader.Next(); ++i) {\n" + " if (static_cast(res.GetSetupStatus()) < 0) break;\n" + " (void)res->size();\n" + " ++nRead;\n" + " }\n" + " auto t1 = std::chrono::steady_clock::now();\n" + " Long64_t bytesRead = f->GetBytesRead() - bytesBefore;\n" + " if (nRead == 0 || bytesRead <= 0) { printf(\"PROBE_FAIL\\n\"); return; }\n" + " double sec = std::chrono::duration(t1 - t0).count();\n" + " printf(\"PROBE_OK %lld %f\\n\", (long long)bytesRead, sec);\n" + "}\n"; + + std::string src(templateSrc); + const std::string placeholder = "@FUNC@"; + auto pos = src.find(placeholder); + if (pos != std::string::npos) { + src.replace(pos, placeholder.size(), funcName); + } + + const std::string path = "/tmp/" + funcName + ".C"; + std::ofstream out(path); + out << src; + out.close(); + return path; +} + +// Probes one candidate replica in a SEPARATE OS PROCESS, bounded by the `timeout` command, so a genuine +// hang (not just a slow-but-working transfer) can be killed without touching this process's own +// TGrid/JAlien connection. An in-process std::async-based timeout was tried and rejected: +// std::future's destructor from std::launch::async BLOCKS until the task finishes regardless of what +// wait_for() returned, so it doesn't actually bound wall-clock time -- and leaking the future to dodge +// that reopens a real concurrent-JAlien-access crash. A real OS process boundary is the only mechanism +// that is both a genuine timeout AND safe against that crash. +// +// Confirmed necessary by a real GRID failure: a job hung completely (~1% CPU for 15 minutes, no +// progress) inside an in-process probe's first candidate until AliEn's idle-CPU watchdog killed the +// whole job. A bounded probe lets it fall through to the next candidate instead of losing the slot. +bool probeOneReplicaSubprocess(const std::string& plainLFN, const std::string& se, double& bytesPerSec, + int timeoutSec = 30, Long64_t probeEntries = 3) +{ + bytesPerSec = -1.0; + const std::string funcName = fmt::format("probeSEChild{}", static_cast(getpid())); + const std::string macroPath = writeProbeChildMacro(funcName); + + const std::string cmd = fmt::format( + "timeout {}s root.exe -b -q -l -x '{}(\"{}\", \"{}\", {})' 2>&1", + timeoutSec, macroPath, plainLFN, se, probeEntries); + + std::string output; + { + std::unique_ptr pipe(popen(cmd.c_str(), "r"), pclose); + if (pipe) { + std::array buffer; + while (fgets(buffer.data(), buffer.size(), pipe.get()) != nullptr) { + output += buffer.data(); + } + } + } + std::remove(macroPath.c_str()); + + // Parse line-by-line rather than a single sscanf over the whole output -- ROOT's own startup banner + // and "Info in " lines precede the actual PROBE_OK/PROBE_FAIL line. + std::istringstream iss(output); + std::string line; + while (std::getline(iss, line)) { + long long bytesRead = 0; + double sec = 0.0; + if (sscanf(line.c_str(), "PROBE_OK %lld %lf", &bytesRead, &sec) == 2) { + if (sec > 0 && bytesRead > 0) { + bytesPerSec = static_cast(bytesRead) / sec; + return true; + } + return false; + } + if (line.rfind("PROBE_FAIL", 0) == 0) { + return false; + } + } + return false; // empty/no matching line -- timed out (killed by `timeout`) or crashed +} + +// Probes each candidate replica in turn (via probeOneReplicaSubprocess above) and returns the SE with +// the highest measured throughput, or "" if every candidate failed or timed out. +std::string probeFastestSE(const std::string& plainLFN, const std::vector& replicas, + Long64_t probeEntries = 3, int timeoutSec = 30) +{ + std::string bestSE; + double bestBytesPerSec = -1.0; + for (const auto& r : replicas) { + double bps = -1.0; + if (!probeOneReplicaSubprocess(plainLFN, r.se, bps, timeoutSec, probeEntries)) { + LOGP(warning, "SE probe: {} failed or timed out after {}s", r.se, timeoutSec); + continue; + } + LOGP(info, "SE probe: {} -> {:.1f} MB/s", r.se, bps / (1024.0 * 1024.0)); + if (bps > bestBytesPerSec) { + bestBytesPerSec = bps; + bestSE = r.se; + } + } + return bestSE; +} + +std::vector loadRunTimeSpans(const std::string& flname, int onlyRun, const std::string& selection); +std::vector getInputFileList(const std::string& fileInput) +{ + std::vector fileList; + std::vector fileListVerified; + // check if only one input file (a txt file contaning a list of files is provided) + if (fileInput.length() > 3 && fileInput.substr(fileInput.length() - 3, 3) == "txt") { + LOGP(info, "Reading files from input file list {}", fileInput); + std::ifstream is(fileInput); + std::istream_iterator start(is); + std::istream_iterator end; + fileList.insert(fileList.begin(), start, end); + } else { + fileList.push_back(fileInput); + } + + // fastestSE is probed once, for the first alien:// file, and reused for the rest of this slot's files + // (same production run -> same SE set, in practice), avoiding a ~10GB-file probe cost per file. If a + // later file isn't actually hosted on the cached SE, doFileProcessing's own open has a fallback to + // unforced resolution -- so a stale cache can't silently drop a file, only cost a little speed for + // that one file. + std::string fastestSE; + for (auto file : fileList) { + if ((file.find("alien://") == 0) && !gGrid && !TGrid::Connect("alien://")) { + LOGP(fatal, "Failed to open alien connection"); + } + if (gSystem->Getenv("FORCESE") && !TString(file.data()).EndsWith(gSystem->Getenv("FORCESE"))) { + file += "?se="; + file += gSystem->Getenv("FORCESE"); + } else if (!gSystem->Getenv("FORCESE") && file.rfind("alien://", 0) == 0) { + if (fastestSE.empty()) { + std::string plainLFN = file.substr(std::string("alien://").size()); + auto slash = plainLFN.find_first_not_of('/'); + plainLFN = (slash == std::string::npos) ? "/" : "/" + plainLFN.substr(slash); + auto replicas = getAlienReplicas(plainLFN); + if (!replicas.empty()) { + fastestSE = (replicas.size() == 1) ? replicas.front().se : probeFastestSE(plainLFN, replicas); + if (!fastestSE.empty()) { + LOGP(info, "Auto-selected fastest SE {} for this slot's alien:// files", fastestSE); + } + } + } + if (!fastestSE.empty()) { + file += "?se="; + file += fastestSE; + } + } + fileListVerified.push_back(file); + } + + if (fileListVerified.size() == 0) { + LOGP(error, "No input files to process"); + } + return fileListVerified; +} + +bool revalidateTrack(const TrackData& trk, const SpacePointsCalibConfParam& params) +{ + + if (hasPositiveFilterFlag(trk)) { + return false; + } + + if (fabs(trk.par.getTgl()) > params.maxZ2X) { + return false; + } + if (trk.nClsITS < params.minITSNCls) { + return false; + } + if (trk.nClsTPC < params.minTPCNCls) { + return false; + } + // No TRD-based cuts here (neither on the tracklet count nor on chi2TRD): this macro does not use TRD. + // track quality cuts + if (trk.chi2ITS / trk.nClsITS > params.maxITSChi2) { + return false; + } + if (trk.chi2TPC / trk.nClsTPC > params.maxTPCChi2) { + return false; + } + + if (params.cutOnDCA) { + auto propagator = o2::base::Propagator::Instance(); + // o2::track::TrackPar trkPar(trk.x, trk.alpha, trk.p); // use this line, in case ClassDef version of TrackData < 4 + o2::track::TrackPar trkPar = trk.par; + if (!propagator->propagateToX(trkPar, 0, propagator->getNominalBz())) { + return false; + } + if (trkPar.getX() * trkPar.getX() + trkPar.getY() * trkPar.getY() > params.maxDCA * params.maxDCA) { + return false; + } + } + return true; +} + +// Check that a TTreeReaderValue was actually bound to a branch of the expected type. +// Unlike SetBranchAddress (which silently tolerated a missing/mismatching branch), dereferencing a +// TTreeReaderValue that failed to set up returns a null proxy and segfaults, so this has to be +// checked once per file before the data is used. The setup status is only final after the first +// entry has been loaded, so call this after SetEntry(). +template +bool checkReaderValue(const TTreeReaderValue& value, const int iThread, const std::string& fileName) +{ + // all failure codes of ESetupStatus are negative (kSetupMatch is 0, the other success codes are positive) + const auto status = value.GetSetupStatus(); + if (static_cast(status) < 0) { + LOGP(warning, "[Thread{}] Branch '{}' could not be set up (setup status {}) in file {}", iThread, value.GetBranchName(), static_cast(status), fileName); + return false; + } + return true; +} + +// Pool size per voxel/charge. Compile-time rather than a runtime parameter so that circleCenters and +// circleRadii below can be fixed-size std::array instead of heap-allocated std::vector: with ~3M +// voxels in a realistic bin configuration, per-voxel heap allocations add up to roughly 4.2 GB of +// resident memory and 12M+ small mallocs for these two members alone. Changing the pool size means +// editing this line and recompiling. +constexpr int NPool = 15; + +// Warm-up threshold for the no-input-map DZ rolling-average correction: minimum cumulative NP/PP/NN +// sample count (see VoxelData below) required before a voxel's residualsAll[0] (dX) is trusted enough +// to shift the Z propagation. See the call site in doFileProcessing for the full reasoning. +constexpr int MinDxSamplesForZCorr = 20; + +// Storage for circleCenters: pure scratch coordinates (never serialized/drawn), populated from +// float-precision inputs (xycircle.xC/yC, track/voxel positions) in the first place, so float is +// enough -- consistent with circleRadii already being float. Converts implicitly to Vec3d at the point +// of use (getIntCircles computes in double internally regardless of the float input). +struct Vec3f { + float x = 0.f, y = 0.f, z = 0.f; + operator Vec3d() const { return Vec3d{x, y, z}; } +}; + +// shared circle pool for one voxel; one instance per voxel (not per thread), guarded by its own mutex +struct VoxelData { + std::mutex mtx; //! per-voxel mutex, shared across threads + std::array, 2> circleCenters; // [charge] was vec_TV3_circle_center_thread + std::array, 2> circleRadii; // [charge] was vec_TV3_circle_radius_thread + std::array poolCounter{}; // [charge] was vec_counter_thread + + std::array residualsAll{}; // was vec_residualsAll_thread; [2] (Z) is a running sum, see counterZAll + std::array residualsNP{}; // was vec_residualsNP_thread + std::array residualsPP{}; // was vec_residualsPP_thread + std::array residualsNN{}; // was vec_residualsNN_thread + int counterNP = 0; // was vec_residuals_counterNP_thread + int counterPP = 0; // was vec_residuals_counterPP_thread + int counterNN = 0; // was vec_residuals_counterNN_thread + int counterZAll = 0; // was vec_residuals_counterZAll_thread +}; + +// One TimeFrame's data, fully OWNED (copied out of the TTreeReaderValues rather than referencing them). +// TTreeReaderValue::operator* reuses the same underlying storage on every SetEntry -- a background +// producer thread reading TF N+1 into that storage while the consumer is still processing TF N's tracks +// would race and corrupt data, the same class of hazard already found once in this file (a lazy- +// deserialization race across concurrent track-worker threads). Copying the three vectors out per TF +// avoids that; the copy itself is small next to the per-TF track-processing time. +struct TFPackage { + int iEntry = 0; + std::vector trackRefsVec; + std::vector trackDataVec; + std::vector unbinnedResidualsVec; +}; + +// Bounded single-producer/single-consumer queue of TFPackages. Lets one background thread stay a few +// TimeFrames ahead of the (I/O-free -- see doFileProcessing's own track-loop-parallelism notes) track +// processing, so a TTreeCache refill on some later TF can overlap with the current TF's track-worker +// compute instead of blocking it. This is meant to hide the periodic TTreeCache-refill spikes this +// pipeline's I/O shows in practice, and only pays off paired with a moderate TTreeCache size -- too +// large a cache makes individual refills bigger than any reasonable queue depth can absorb. +class BoundedTFQueue +{ + public: + explicit BoundedTFQueue(size_t maxDepth) : mMaxDepth(maxDepth) {} + + // Producer side. Blocks while the queue is already at capacity. + void push(std::unique_ptr pkg) + { + std::unique_lock lock(mMtx); + mNotFull.wait(lock, [this] { return mQueue.size() < mMaxDepth; }); + mQueue.push_back(std::move(pkg)); + lock.unlock(); + mNotEmpty.notify_one(); + } + + // Producer side, called once (file fully read, or the maxTracks quota was reached). + void setDone() + { + { + std::lock_guard lock(mMtx); + mDone = true; + } + mNotEmpty.notify_one(); + } + + // Consumer side. Returns nullptr once the producer is done AND the queue has been fully drained. + std::unique_ptr pop() + { + std::unique_lock lock(mMtx); + mNotEmpty.wait(lock, [this] { return !mQueue.empty() || mDone; }); + if (mQueue.empty()) { + return nullptr; + } + auto pkg = std::move(mQueue.front()); + mQueue.pop_front(); + lock.unlock(); + mNotFull.notify_one(); + return pkg; + } + + // Diagnostic only: how many packages are sitting ready right now. Lets the consumer log whether the + // producer is comfortably ahead (queue usually near mMaxDepth) or struggling to keep up (queue usually + // near empty) -- distinguishes "the buffer is the wrong depth" from "the producer itself is too slow + // to ever fill it, no matter how deep it is". + size_t size() const + { + std::lock_guard lock(mMtx); + return mQueue.size(); + } + + private: + const size_t mMaxDepth; + mutable std::mutex mMtx; + std::condition_variable mNotEmpty; + std::condition_variable mNotFull; + std::deque> mQueue; + bool mDone = false; +}; + +// How many TimeFrames the background producer may stay ahead of track-processing. Absorbs some of the +// periodic TTreeCache-refill spikes this pipeline's I/O shows in practice, though 10 (the default +// below) isn't enough to fully hide the biggest ones -- a much larger depth would be needed for that, +// at a real memory cost (a single TF can carry tens of thousands of tracks). Overridable via +// SCDCALIB_TF_QUEUE_DEPTH (no recompile) to tune against a real workload's spike size. +constexpr size_t TFQueueDepthDefault = 10; + +void doFileProcessing(const int iThread, + const int nFileThreads, + const int maxTrackWorkers, + const long firstTFTime, + const long lastTFTime, + const bool invertBadRange, + const float maxdEdx, + const float maxdEdxExp, + const float maxDevdEdxOverExp, + const float skipEdgePads, + std::vector& nEdgeClustersSkipped_thread, + std::vector& nTracksSkippedByBadRangeList_thread, + std::vector& nTFs_thread, + std::vector& nTFsSkippedByBadRangeList_thread, + std::vector& nTFsSkippedByTimeWindow_thread, + const std::string voxMapInput, + const GID::mask_t sources, + const int64_t orbitResetTimeMS, + const float magfieldvalue, + const std::vector fileList, + const int maxTracksPerSlice, + const Long64_t maxTracks, + std::atomic& nTracksProcessed, + const std::array, NSectors>& voxelResults, // read-only input correction map, shared across threads (no per-thread copy) + std::vector>& badRanges_thread, + const TrackResiduals& trackResiduals, // read-only: findVoxelBin/getVoxelCoordinates/getGlbVoxBin are all const, safe to share across threads + const float maxDistIntCls, + const int nY2XBins, + const int nZ2XBins, + std::vector& voxels, // flat [sec*152*nY2XBins*nZ2XBins + ix*nY2XBins*nZ2XBins + iy*nZ2XBins + iz] -- shared across threads, one instance per voxel + std::vector& totalBytesReadPerf_thread, + std::vector& lumiEntriesCTP_thread, + std::vector& lumiSumCTP_thread, + std::vector>& orbitsSel_thread, + std::vector>& ctpLumiSel_thread, + std::vector>& timeMSsel_thread) +{ + // Get Mapper + const Mapper& mapper = Mapper::instance(); + + // yMaxCentrePadByRow only depends on the pad row (152 values) -- precomputed once per file-thread + // here, gated on skipEdgePads since that's its only use, rather than via a Mapper lookup per residual. + std::array yMaxCentrePadByRow{}; + if (skipEdgePads) { + for (int irow = 0; irow < NRows; ++irow) { + yMaxCentrePadByRow[irow] = mapper.getPadCentre(o2::tpc::PadPos(irow, 0)).Y() - mapper.getPadRegionInfo(o2::tpc::Mapper::REGION[irow]).getPadWidth() / 2; + } + } + + // Obtain configuration per thread + const SpacePointsCalibConfParam& params_thread = SpacePointsCalibConfParam::Instance(); + + // Per-thread input handles and I/O-monitoring state. These are plain locals: each file-thread only ever + // touches its own, so there is nothing to share with the other threads or hand back to the caller. + // + // DECLARATION ORDER IS LOAD-BEARING. Locals are destroyed in reverse order of declaration, and these + // objects reference each other: a TTreeReaderValue refers to its TTreeReader, which refers to a TTree + // owned by the TFile, and TTreePerfStats refers to that tree too. Declaring the file first and the + // reader values last therefore tears them down in the only safe order -- values, then perf stats, then + // readers, then the file. Do not reorder these. + std::unique_ptr inputFile; + std::unique_ptr treeUnbinnedResiduals; + std::unique_ptr treeTrackData; + std::unique_ptr treeRecords; + std::unique_ptr perfStats; + std::unique_ptr>> unbinnedResiduals; // unbinned residuals input + std::unique_ptr>> trackRefs; // track references for the unbinned residuals + std::unique_ptr>> trackData; // additional track info (chi2, nClusters, track parameters) + std::unique_ptr>> orbits; // first orbit of each TF in the input data + std::unique_ptr> lumiTF; // lumi info + + // Previous I/O sample, for the instantaneous-rate log inside the TF loop. + Long64_t perfLastBytes = 0; + std::chrono::steady_clock::time_point perfLastSample; + + int trackCounter_local{0}; + + // Track-loop parallelism, WITHIN this one file-thread only -- active when nFileThreads==1 (GRID/ + // alien:// mode, forced by getInputFileList() for TGrid safety) or when there's only one input file + // (otherwise every other core would sit idle). Safe because the track/cluster loop body touches no + // TGrid/CCDB/file I/O (SetEntry() already happened before this point), only in-memory TF data plus the + // same per-voxel mutex (vox.mtx) multi-file-thread mode already relies on. Each worker gets its own + // copy of every piece of mutated state (RNG, scratch coordinates, counters) -- sharing any of it + // across workers would be a silent data race. Otherwise nTrackWorkers is forced to 1 (serial, via + // worker index 0) to avoid oversubscribing on top of the file-threads. + int nTrackWorkers = 1; + if (nFileThreads == 1 || fileList.size() == 1) { + if (maxTrackWorkers > 0) { + // Explicit override (e.g. the number of cores actually allocated to a batch/GRID job). + // hardware_concurrency() reports the machine's core count, not the job's allocation -- an 8-core + // GRID job auto-detects 32 and oversubscribes 4x -- so when the caller knows, trust it instead. + nTrackWorkers = maxTrackWorkers; + } else { + unsigned hc = std::thread::hardware_concurrency(); + nTrackWorkers = (hc == 0) ? 8 : static_cast(hc); + if (nTrackWorkers > 32) { + nTrackWorkers = 32; + } + } + } + LOGP(info, "[Thread_{}] Using {} track-worker thread(s) for the track loop", iThread, nTrackWorkers); + std::vector clsPosWorker(nTrackWorkers); + // Per-worker scratch: this track's position at the current row. + std::vector trackPosAtRow_worker(nTrackWorkers); + std::vector nEdgeClustersSkipped_worker(nTrackWorkers, 0); + std::vector trackCounter_local_worker(nTrackWorkers, 0); + + ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + // ===| FILE LOOP |=================================================================================================================== + ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + // This thread's share of the input files, as an explicit work list: this is what lets a file that + // looks unhealthy be pushed to the BACK and retried later instead of being processed now or dropped -- + // see the health checks below for why deferring beats both. Appending while iterating by index is + // safe: the bound is re-read every iteration and no iterators are held. + std::vector workList; + for (int i = iThread; i < int(fileList.size()); i += nFileThreads) { + workList.push_back(i); + } + // How many times each file has already been deferred, so a persistently sick one cannot cycle forever. + std::vector deferCount(fileList.size(), 0); + int maxFileDeferrals = 1; + if (const char* envMaxDeferrals = gSystem->Getenv("SCDCALIB_MAX_FILE_DEFERRALS")) { + maxFileDeferrals = std::atoi(envMaxDeferrals); + } + // Health-probe timings (seconds) of the files accepted so far, used as this job's own baseline. An + // absolute threshold cannot work here: measured healthy timings differ by ~6x between GRID regimes, + // so what matters is how a file compares to the others in ITS OWN slot, not to a constant. + std::vector probeTimes; + double probeSlowFactor = 5.0; + if (const char* envSlowFactor = gSystem->Getenv("SCDCALIB_PROBE_SLOW_FACTOR")) { + probeSlowFactor = std::atof(envSlowFactor); + } + // Fallback used until enough samples exist for a median to mean anything (and if the very first file + // of a slot is the sick one, this is the only thing standing between us and the stall). + double probeAbsMaxSec = 15.0; + if (const char* envProbeAbsMax = gSystem->Getenv("SCDCALIB_PROBE_ABS_MAX_SEC")) { + probeAbsMaxSec = std::atof(envProbeAbsMax); + } + // Floor under the median-relative threshold once >=3 samples exist. Without it, a fast-regime median + // (tens of ms) makes ordinary network jitter look "5x slower than usual" and trip the probe -- and a + // file unlucky enough to trip it twice is dropped for good (see deferFile below), losing real data to + // noise. The probe's actual target is stalls an order of magnitude bigger (~20s against a ~1.5s + // baseline, in the real case this was built around), so a small floor can't mask that while still + // absorbing sub-second jitter in a fast regime. + double probeMinThresholdSec = 2.0; + if (const char* envProbeMinThreshold = gSystem->Getenv("SCDCALIB_PROBE_MIN_THRESHOLD_SEC")) { + probeMinThresholdSec = std::atof(envProbeMinThreshold); + } + + for (size_t iWork = 0; iWork < workList.size(); ++iWork) { + const int iFile = workList[iWork]; + // Get filename from fileList + auto fileName = fileList[iFile]; + + // Check if enough tracks are processed + if ((maxTracks > 0) && (nTracksProcessed.load(std::memory_order_relaxed) > maxTracks)) { + LOGP(info, "[Thread_{}] Maximum number of requested tracks processed {} > {} ({}), will not process further files", iThread, nTracksProcessed.load(std::memory_order_relaxed), maxTracks, maxTracksPerSlice); + break; + } + + if (gSystem->Getenv("FORCESE") && !TString(fileName.data()).EndsWith(gSystem->Getenv("FORCESE"))) { + fileName += "?se="; + fileName += gSystem->Getenv("FORCESE"); + } + + // Open tree and set branches + LOGP(info, "[Thread{}] Processing input file {}", iThread, fileName); + unbinnedResiduals.reset(nullptr); + trackRefs.reset(nullptr); + lumiTF.reset(nullptr); + trackData.reset(nullptr); + orbits.reset(nullptr); + // Must be destroyed here, before the readers/file below: TTreePerfStats registers itself as the + // process-wide gPerfStats global and keeps raw TFile*/TTree* pointers to the file it watches. + // Leaving it alive past this point means gPerfStats still points at this (about to be destroyed) + // file while the NEXT file's setup does real reads -- TFile::ReadBuffer() unconditionally calls + // gPerfStats->FileReadEvent(thatFile, ...), which compares thatFile against the dangling fFile + // pointer; if the allocator hands the new TFile the same address the old one was just freed from + // (a real, common allocator pattern for same-sized immediate reuse), the comparison spuriously + // matches and it dereferences the also-dangling fTree -- a real, reproduced segfault on the 2nd+ + // file of a multi-file run, verified against the installed ROOT's TTreePerfStats.cxx/TFile.cxx + // source. perfStats's own destructor is safe to call here (only clears gPerfStats if it's still the + // registered one; never dereferences fTree/fFile), so this alone fixes it. + perfStats.reset(nullptr); + treeUnbinnedResiduals.reset(nullptr); + treeTrackData.reset(nullptr); + treeRecords.reset(nullptr); + + const auto openStart = std::chrono::steady_clock::now(); + inputFile.reset(TFile::Open(fileName.c_str())); + if ((!inputFile || inputFile->IsZombie())) { + // A forced-SE URL (FORCESE or the auto-picked/cached fastest SE, see getInputFileList) can fail if + // this particular file isn't actually hosted there -- fall back to unforced alien:// resolution + // rather than skipping the file outright, since a stale SE cache would otherwise silently drop data. + auto sePos = fileName.find("?se="); + if (sePos != std::string::npos) { + std::string fallbackName = fileName.substr(0, sePos); + LOGP(warning, "[Thread_{}] Forced-SE open failed for {}, retrying without SE override: {}", iThread, fileName, fallbackName); + inputFile.reset(TFile::Open(fallbackName.c_str())); + } + } + // Gates the "[prefetch diag]" per-TF/per-file I/O diagnostics further below: real signal for + // diagnosing GRID stalls/CPU-idle-watchdog kills (the reason this machinery exists at all), but + // pure noise on a local/Lustre run with many file-threads where none of that risk applies -- a + // 36-file-thread local run was producing an unreadable flood of per-TF consumer lines otherwise. + const bool isAlienFile = fileName.find("alien://") != std::string::npos; + if (isAlienFile) { + if (inputFile && !inputFile->IsZombie()) { + inputFile->SetBufferSize(4000000); + LOGP(info, "[Thread_{}] Set buffer size to {}", iThread, inputFile->GetBufferSize()); + } else { + LOGP(info, "[Thread_{}] TFile {} is empty", iThread, fileName); + } + } + + if (!inputFile || inputFile->IsZombie()) { + LOGP(warning, "[Thread{}] Skipping file {}", iThread, fileName); + continue; + } + + // Deferring an unhealthy-looking file, rather than skipping it. Observed for real: a file that was + // reproducibly slow in one job read at full speed a short time later -- the slowness is transient + // storage-server state, not a property of the file. So dropping it outright throws away data that + // would very likely have been fine. Pushing it to the back of this thread's work list costs exactly + // what skipping costs right now, but gives the server time to recover, and if maxTracks stops the + // job first the file is never touched again at all. The caller always moves on to the next file. + auto deferFile = [&](const char* reason, double measured, double threshold) { + if (deferCount[iFile] < maxFileDeferrals) { + ++deferCount[iFile]; + workList.push_back(iFile); + LOGP(warning, "[Thread{}] {} for {} ({:.1f} s vs {:.1f} s threshold) -- deferring it to the end of this thread's file list (attempt {} of {}); storage slowness has been seen to be transient, so it may read fine later, and maxTracks may stop the job before we return to it", + iThread, reason, fileName, measured, threshold, deferCount[iFile], maxFileDeferrals); + return; + } + LOGP(warning, "[Thread{}] {} for {} ({:.1f} s vs {:.1f} s threshold) and it has already been deferred {} time(s) -- skipping this file for good", + iThread, reason, fileName, measured, threshold, deferCount[iFile]); + }; + + // --- Unhealthy-replica gate ------------------------------------------------------------------- + // An abnormally slow TFile::Open is the earliest sign a replica's storage server is struggling, + // and the last point we can walk away cheaply: the warm-up read right after this pulls a whole + // TTreeCache-sized chunk in one uninterruptible call, so if the server stalls there the job hangs + // until AliEn's ~15-minute idle-CPU watchdog kills it, discarding all work already done. Skipping + // one file here costs a fraction of a slot's statistics, so the trade is heavily one-sided. 15 s + // threshold, calibrated against real GRID opens (absolute, not relative to this job's own timings + // -- may need revisiting on very different links). Tune via SCDCALIB_MAX_FILE_OPEN_SEC; <= 0 + // disables it. + double maxFileOpenSec = 15.0; + if (const char* envMaxOpenSec = gSystem->Getenv("SCDCALIB_MAX_FILE_OPEN_SEC")) { + maxFileOpenSec = std::atof(envMaxOpenSec); + } + const double openSec = std::chrono::duration(std::chrono::steady_clock::now() - openStart).count(); + if (maxFileOpenSec > 0 && openSec > maxFileOpenSec) { + deferFile("Slow file open (SCDCALIB_MAX_FILE_OPEN_SEC)", openSec, maxFileOpenSec); + continue; + } + + treeUnbinnedResiduals = std::make_unique("unbinnedResid", inputFile.get()); + if (!treeUnbinnedResiduals->GetTree()) { + LOGP(warning, "[Thread{}] Could not get tree 'unbinnedResid' from file {}. Skipping file!", iThread, fileName); + continue; + } + + // GetEntries() only reads TTree header metadata, not branch data -- cheap even on alien://, unlike + // the TTreeCache/branch setup further below. Fetched here (still before that setup) so the fail-fast + // time-window check below can run before anything expensive touches the network. + const auto nTFEntries = treeUnbinnedResiduals->GetEntries(); + if (nTFEntries <= 0) { + LOGP(warning, "[Thread{}] Tree 'unbinnedResid' in file {} has no entries. Skipping file!", iThread, fileName); + continue; + } + + treeRecords = std::make_unique("records", inputFile.get()); + if (!treeRecords->GetTree()) { + LOGP(warning, "[Thread{}] Could not get tree 'records' from file {}. Skipping file!", iThread, fileName); + continue; + } + orbits = std::make_unique>>(*treeRecords, "firstTForbit"); + // Real data has exactly one entry in 'records', but MC input can have several (e.g. one per + // simulation chunk merged into this file) -- each entry's own 'firstTForbit' only covers that + // chunk's TFs, so reading just entry 0 silently truncated the orbit list on MC, tripping the + // length check below and skipping the whole file. Concatenate every entry's vector instead, in + // entry order, to rebuild the same flat, TF-index-ordered list this file's real-data path already + // produced from its single entry (verified for real: MC input observed with 5 'records' entries). + std::vector combinedOrbits; + { + const Long64_t nRecordsEntries = treeRecords->GetEntries(); + bool recordsOk = true; + for (Long64_t ie = 0; ie < nRecordsEntries; ++ie) { + if (treeRecords->SetEntry(ie) != TTreeReader::kEntryValid || + !checkReaderValue(*orbits, iThread, fileName)) { + recordsOk = false; + break; + } + const auto& thisEntryOrbits = **orbits; + combinedOrbits.insert(combinedOrbits.end(), thisEntryOrbits.begin(), thisEntryOrbits.end()); + } + if (!recordsOk) { + LOGP(warning, "[Thread{}] Could not load the orbits from tree 'records' in file {}. Skipping file!", iThread, fileName); + continue; + } + } + // the orbit list is indexed with the TF index of the unbinnedResid tree below, so it has to be at least as long + if (static_cast(combinedOrbits.size()) < nTFEntries) { + LOGP(error, "[Thread{}] 'firstTForbit' has fewer entries than the residual tree has TFs ({} vs {}) in file {}. Skipping file!", iThread, + combinedOrbits.size(), nTFEntries, fileName); + continue; + } + + // Set timeStamp for processing, get this file's [min,max] orbit-derived time range, and find the + // first TF entry actually inside the requested window -- all in one pass over the (already + // downloaded) 'firstTForbit' array. Bounded to the first nTFEntries entries: orbits can have extra + // trailing entries beyond what the residual tree actually has (see the size check above) that don't + // correspond to any real TF and must not feed either the fail-fast check below or the warm-up entry. + uint32_t minFirstOrbit = -1; + uint32_t maxFirstOrbit = 0; + Long64_t warmupEntry = 0; + bool foundWarmupEntry = (firstTFTime <= 0); // no time filter set: entry 0 is always fine to warm up on + { + const auto& orbitsVec = combinedOrbits; + for (Long64_t i = 0; i < nTFEntries; ++i) { + const uint32_t orbit = orbitsVec[i]; + if (orbit < minFirstOrbit) { + minFirstOrbit = orbit; + } + if (orbit > maxFirstOrbit) { + maxFirstOrbit = orbit; + } + if (!foundWarmupEntry) { + const int64_t t = orbitResetTimeMS + orbit * o2::constants::lhc::LHCOrbitMUS * 1.e-3; + if (t >= firstTFTime && t <= lastTFTime) { + warmupEntry = i; + foundWarmupEntry = true; + } + } + } + } + // ---| Fail fast on files with zero overlap with the requested time window |--- + if (firstTFTime > 0) { + const int64_t fileMinTimeMS = orbitResetTimeMS + minFirstOrbit * o2::constants::lhc::LHCOrbitMUS * 1.e-3; + const int64_t fileMaxTimeMS = orbitResetTimeMS + maxFirstOrbit * o2::constants::lhc::LHCOrbitMUS * 1.e-3; + if (fileMaxTimeMS < firstTFTime || fileMinTimeMS > lastTFTime) { + LOGP(warning, "[Thread{}] File {} has no TF inside the requested time window [{}, {}] ms (file spans [{}, {}] ms) -- skipping the whole file without downloading its residual tree", + iThread, fileName, firstTFTime, lastTFTime, fileMinTimeMS, fileMaxTimeMS); + // Counted as if the per-TF loop below had visited and skipped each entry, so nTFs stays a + // consistent denominator for the skip-fraction summary at the end. + nTFs_thread[iThread] += nTFEntries; + nTFsSkippedByTimeWindow_thread[iThread] += nTFEntries; + continue; + } + } + + // --- Read-health probe, BEFORE the TTreeCache is enabled -------------------------------------- + // A slow open alone can miss a replica that opens fine but stalls on the first real read, so this + // probes a small read directly, before SetCacheSize/AddBranchToCache below turn the first read into + // a whole cache-sized fetch that can hang with nothing able to interrupt it. Probes 'trackData' (a + // small member-split branch) rather than the already-read 'records' tree, which sits next to the + // file header and reads fast regardless of whether the replica then stalls on real data. Uses + // TBranch::GetEntry directly rather than a TTreeReader, to avoid creating a second reader on a tree + // that gets its real one further below. A missing branch skips the probe rather than failing it -- + // this is a health check, not a validity check. + { + TTree* probeTree = dynamic_cast(inputFile->Get("trackData")); + TBranch* probeBranch = probeTree ? probeTree->GetBranch("trk.nClsTPC") : nullptr; + if (probeBranch && probeTree->GetEntries() > 0) { + const Long64_t probeEntry = std::min(warmupEntry, probeTree->GetEntries() - 1); + const auto probeStart = std::chrono::steady_clock::now(); + const Int_t probeBytes = probeBranch->GetEntry(probeEntry); + const double probeSec = std::chrono::duration(std::chrono::steady_clock::now() - probeStart).count(); + + // If the read returned nothing, the measurement says nothing about the replica's health -- fall + // through to the normal path rather than judging the file on it. Deliberately fail-open: a probe + // that cannot measure must not be able to reject files, or an unexpected branch layout would + // quietly defer every file in the slot. + if (probeBytes <= 0) { + LOGP(info, "[Thread{}] Read-health probe returned no data for {} (entry {}) -- skipping the health check for this file", iThread, fileName, probeEntry); + } else { + // Baseline: the median probe time of files already accepted in this slot. Below 3 samples a + // median is meaningless, so fall back to the absolute guard. + double probeThreshold = probeAbsMaxSec; + if (probeTimes.size() >= 3) { + std::vector sorted(probeTimes); + std::nth_element(sorted.begin(), sorted.begin() + sorted.size() / 2, sorted.end()); + const double median = sorted[sorted.size() / 2]; + probeThreshold = std::max(probeSlowFactor * median, probeMinThresholdSec); + } + if (probeThreshold > 0 && probeSec > probeThreshold) { + deferFile("Slow read probe (SCDCALIB_PROBE_SLOW_FACTOR/SCDCALIB_PROBE_ABS_MAX_SEC)", probeSec, probeThreshold); + continue; + } + // Only healthy files feed the baseline, so one sick file cannot raise the bar for the next. + probeTimes.push_back(probeSec); + } + } + } + + // I/O throughput monitor -- large cache so async prefetch has room to read many baskets ahead; only + // branches actually accessed get cached/prefetched. TTreePerfStats records raw bytes read from the + // remote file, so its rate is the actual download speed. + // + // 256 MB, not 512 MB: a real GRID job hit an 11+ minute stall refilling a single 512MB chunk (~0.7 + // MB/s vs. 45-100 MB/s for every other chunk on the same SE), then a second stall that never + // recovered, losing the whole job to AliEn's idle-CPU watchdog. A smaller chunk halves the + // worst-case single-chunk wait and lets a persistently slow file be abandoned sooner -- a + // reliability trade-off against 512MB's better throughput (~17% vs ~13% wall-clock reduction) when + // nothing is stalling. Overridable via SCDCALIB_CACHE_SIZE_MB. + { + TTree* residTree = treeUnbinnedResiduals->GetTree(); + int dbgCacheSizeMB = 256; + if (const char* envCacheSizeMB = gSystem->Getenv("SCDCALIB_CACHE_SIZE_MB")) { + dbgCacheSizeMB = std::atoi(envCacheSizeMB); + } + if (isAlienFile) { + LOGP(info, "[Thread{}] [prefetch diag] TTreeCache size = {} MB (SCDCALIB_CACHE_SIZE_MB, default 256)", iThread, dbgCacheSizeMB); + } + residTree->SetCacheSize(static_cast(dbgCacheSizeMB) * 1024 * 1024); + residTree->AddBranchToCache("*", true); + perfStats = std::make_unique(fmt::format("ioperf_{}_{}", iThread, iFile).data(), residTree); + } + perfLastBytes = 0; + perfLastSample = std::chrono::steady_clock::now(); + + unbinnedResiduals = std::make_unique>>(*treeUnbinnedResiduals, "res"); + trackRefs = std::make_unique>>(*treeUnbinnedResiduals, "trackInfo"); + lumiTF = std::make_unique>(*treeUnbinnedResiduals, "CTPLumi"); + + // Skip reading UnbinnedResid/TrackDataCompact members that are never used below. 'res' and + // 'trackInfo' are member-split branches, one sub-branch per struct field, and the unused ones + // dominate the file (res.tgSlp alone can be >1 GB in a single input file), so disabling them means + // they are never transferred at all -- measured ~17% fewer bytes read. Must come AFTER the + // TTreeReaderValues above so that it is the final word on these sub-branches' status. + // Do not disable res.dy/dz/y/z/row/sec/rejected or trackInfo.idxFirstResidual/nResiduals/sourceId -- + // all of those are read below. Note trackInfo.filterFlag (on TrackDataCompact) is unused, while the + // separate trk.filterFlag (on TrackData, read further down) is used; they are different fields. + { + TTree* residTree = treeUnbinnedResiduals->GetTree(); + for (const char* br : {"res.tgSlp", "res.channel", "trackInfo.multStack*", + "trackInfo.nExtDetResid", "trackInfo.filterFlag"}) { + residTree->SetBranchStatus(br, 0); + } + } + + // Load one entry once so the reader values get bound, then verify all branches are there before + // anything below dereferences them. Warms up on warmupEntry (computed above, the first entry + // actually inside the requested time window) rather than always entry 0 -- for a file only partially + // overlapping the window, entry 0 is often outside it, and fetching its residual data would be + // exactly the kind of wasted network read the whole-file fail-fast check above targets, just at the + // scale of one TF. + if (treeUnbinnedResiduals->SetEntry(warmupEntry) != TTreeReader::kEntryValid) { + LOGP(warning, "[Thread{}] Could not load entry {} of 'unbinnedResid' from file {}. Skipping file!", iThread, warmupEntry, fileName); + continue; + } + if (!checkReaderValue(*unbinnedResiduals, iThread, fileName) || + !checkReaderValue(*trackRefs, iThread, fileName) || + !checkReaderValue(*lumiTF, iThread, fileName)) { + LOGP(warning, "[Thread{}] Skipping file {}", iThread, fileName); + continue; + } + + // Re-prime the reader values after the pre-scan (which leaves the reader past the last entry), + // so a stale-read (e.g. the bad-range-skip stats below) at warmupEntry dereferences valid data. + if (treeUnbinnedResiduals->SetEntry(warmupEntry) != TTreeReader::kEntryValid) { + LOGP(warning, "[Thread{}] Could not re-load entry {} of 'unbinnedResid' from file {}. Skipping file!", iThread, warmupEntry, fileName); + continue; + } + + { + treeTrackData = std::make_unique("trackData", inputFile.get()); + if (!treeTrackData->GetTree()) { + LOGP(warning, "[Thread{}] Could not get tree 'trackData' from file {}. Skipping file!", iThread, fileName); + continue; + } + { + TTree* trackDataTree = treeTrackData->GetTree(); + for (const char* br : {"trk.gid*", "trk.chi2TRD", "trk.deltaTOF", "trk.nTrkltsTRD", + "trk.clAvailTOF", "trk.TRDTrkltSlope*", "trk.nExtDetResid", + "trk.clIdx.*", "trk.multStack*"}) { + trackDataTree->SetBranchStatus(br, 0); + } + } + trackData = std::make_unique>>(*treeTrackData, "trk"); + if (treeTrackData->GetEntries() != nTFEntries) { + LOGP(error, "[Thread{}] The input trees with unbinned residuals and track information have a different number of entries ({} vs {}). Skipping file!", iThread, + nTFEntries, treeTrackData->GetEntries()); + continue; + } + // Same TF indexing as 'unbinnedResid' -- warm up on the same entry for the same reason (see above). + if (treeTrackData->SetEntry(warmupEntry) != TTreeReader::kEntryValid || + !checkReaderValue(*trackData, iThread, fileName)) { + LOGP(warning, "[Thread{}] Skipping file {}", iThread, fileName); + continue; + } + } + + ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + // ===| TIME FRAME LOOP |============================================================================================================= + // Split into a background PRODUCER thread (skip-checks, SetEntry/warm-up, sync check, lumi/orbit + // bookkeeping, and reading each TF's data) pushing TFPackages into a bounded queue, and this thread + // as CONSUMER, popping a package and dispatching the track workers on it -- see BoundedTFQueue/ + // TFPackage above for why. The producer is the only thread that ever touches + // treeUnbinnedResiduals/treeTrackData/the TTreeReaderValues (exactly one thread doing TGrid/file I/O + // at a time); the consumer never touches them at all, only the packages it pops. + ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + size_t tfQueueDepth = TFQueueDepthDefault; + if (const char* envQueueDepth = gSystem->Getenv("SCDCALIB_TF_QUEUE_DEPTH")) { + tfQueueDepth = static_cast(std::atoi(envQueueDepth)); + } + if (isAlienFile) { + LOGP(info, "[Thread{}] [prefetch diag] TFQueueDepth = {} (SCDCALIB_TF_QUEUE_DEPTH, default {})", iThread, tfQueueDepth, TFQueueDepthDefault); + } + BoundedTFQueue tfQueue(tfQueueDepth); + + // Set by the consumer below when a single tfQueue.pop() waited unreasonably long, checked by the + // producer between TF reads (never inside one -- see the long comment at the consumer's check site + // for why). Assumes the fetch has become persistently slow rather than truly dead: a producer stuck + // forever never reaches this check at all and would need a separate, not-yet-built process-level + // watchdog -- this only shortens the "slow but eventually returns" case, not a genuine hang. + std::atomic abandonFile{false}; + double popWaitAbandonMs = 60000.0; // 60s -- matches the scale of the real stalls that motivated this + if (const char* envAbandonMs = gSystem->Getenv("SCDCALIB_POP_WAIT_ABANDON_MS")) { + popWaitAbandonMs = std::atof(envAbandonMs); + } + + std::thread producerThread([&]() { + // A failed SetEntry() (checked below) only means the tree could not be repositioned -- it says + // nothing about whether a given lazily-read branch's basket actually arrived. TTreeReaderValue + // fetches each branch on its own first dereference for the current entry, and a mid-stream storage + // hiccup there doesn't throw: it prints a ROOT error and leaves the underlying object in an + // unspecified state, which then segfaults wherever it's first used with no indication of why + // (observed for real: an xrootd "Operation expired" mid basket-read, immediately followed by a + // SIGSEGV with only the ROOT error line as a clue). Must be called right after the value's first + // dereference for this entry -- GetReadStatus() reports the status of that specific read. + auto checkReadOk = [&](ROOT::Internal::TTreeReaderValueBase& val, const char* branchName, int iEntry) { + if (val.GetReadStatus() != ROOT::Internal::TTreeReaderValueBase::kReadSuccess) { + LOGP(warning, "[Thread{}] Storage read error on branch '{}' at entry {} of file {} (GetReadStatus={}) -- skipping TF!", + iThread, branchName, iEntry, fileName, static_cast(val.GetReadStatus())); + return false; + } + return true; + }; + for (int iEntry = 0; iEntry < nTFEntries; ++iEntry) { + // Checked here, between TF reads, never inside one -- this point is only ever reached right + // after the previous push() succeeded, so the producer is never mid-call when this fires. See + // the consumer's check site for the full reasoning. + if (abandonFile.load(std::memory_order_relaxed)) { + LOGP(warning, "[Thread{}] Abandoning the rest of file {} ({}/{} TFs read) after a persistently slow fetch", iThread, fileName, iEntry, nTFEntries); + break; + } + ++nTFs_thread[iThread]; + + // Periodic I/O throughput progress log. Now reports the producer's own progress through the + // file, which can run ahead of what the consumer has actually finished processing. + if (iEntry % 50 == 0) { + double instMBps = 0.0, totMB = 0.0; + if (perfStats) { + const auto nowSample = std::chrono::steady_clock::now(); + const double dt = std::chrono::duration(nowSample - perfLastSample).count(); + const Long64_t br = perfStats->GetBytesRead(); + instMBps = (dt > 0) ? (br - perfLastBytes) / (1024.0 * 1024.0) / dt : 0.0; + totMB = br / (1024.0 * 1024.0); + perfLastBytes = br; + perfLastSample = nowSample; + } + const Long64_t nProcessedNow = nTracksProcessed.load(std::memory_order_relaxed); + if (maxTracks > 0) { + LOGP(info, "[Thread{}] TF entry {}/{} | read {:.1f} MB, inst. {:.1f} MB/s | tracks processed {}/{} ({:.1f}%)", + iThread, iEntry, nTFEntries, totMB, instMBps, nProcessedNow, maxTracks, 100.0 * nProcessedNow / maxTracks); + } else { + LOGP(info, "[Thread{}] TF entry {}/{} | read {:.1f} MB, inst. {:.1f} MB/s | tracks processed {} (no limit set)", + iThread, iEntry, nTFEntries, totMB, instMBps, nProcessedNow); + } + } + + // Check if enough tracks are processed + if ((maxTracks > 0) && (nTracksProcessed.load(std::memory_order_relaxed) > maxTracks)) { + LOGP(info, "[Thread{}] Maximum number of requested tracks processed {} > {} ({}), will not process further TFs", iThread, nTracksProcessed.load(std::memory_order_relaxed), maxTracks, maxTracksPerSlice); + break; + } + + // ---| check for TF time acceptance |--- + const int64_t tfTimeInMS = orbitResetTimeMS + combinedOrbits[iEntry] * o2::constants::lhc::LHCOrbitMUS * 1.e-3; + if ((firstTFTime > 0) && (tfTimeInMS < firstTFTime || tfTimeInMS > lastTFTime)) { + if (nTFsSkippedByTimeWindow_thread[iThread] == 0) { + // Log once per thread rather than per TF: this can legitimately fire for every TF of every + // file, e.g. when the requested [firstTFTime,lastTFTime] window contains no data at all + // because a time-slice boundary landed past the end of the run. A per-TF log would then + // produce one line per TF for the entire job. + LOGP(warning, "[Thread{}] TF at index {} (time {} ms, orbit {}) outside requested window [{}, {}] ms -- skipping (will keep happening silently for further TFs outside the window, see final summary for the total count)", + iThread, iEntry, tfTimeInMS, combinedOrbits[iEntry], firstTFTime, lastTFTime); + } + ++nTFsSkippedByTimeWindow_thread[iThread]; + continue; + } + // ---| check for time exclusion list |--- + if (badRanges_thread[iThread].size() > 0) { + bool skip = false; + for (const auto& range : badRanges_thread[iThread]) { + if ((combinedOrbits[iEntry] >= range.from) && (combinedOrbits[iEntry] <= range.to)) { + skip = true; + break; + } + } + if (invertBadRange) { + skip = !skip; + } + if (skip) { + nTracksSkippedByBadRangeList_thread[iThread] += (*trackRefs)->size(); + ++nTFsSkippedByBadRangeList_thread[iThread]; + continue; + } + } + if (params_thread.timeFilter) { + if (tfTimeInMS < params_thread.startTimeMS || tfTimeInMS > params_thread.endTimeMS) { + continue; + } + } + + // --- [prefetch diag]: producer-side read/deserialize timing, kept to confirm the periodic-spike + // I/O pattern still looks the same underneath the producer/consumer split -- expected to be + // mostly hidden from the consumer by TFQueueDepth, not eliminated. --- + const auto dbgTIoStart = std::chrono::steady_clock::now(); + + // ---| Read entries |--- + if (treeUnbinnedResiduals->SetEntry(iEntry) != TTreeReader::kEntryValid) { + LOGP(warning, "[Thread{}] Could not load entry {} of 'unbinnedResid' from file {}. Skipping TF!", iThread, iEntry, fileName); + continue; + } + if (treeTrackData->SetEntry(iEntry) != TTreeReader::kEntryValid) { + LOGP(warning, "[Thread{}] Could not load entry {} of 'trackData' from file {}. Skipping TF!", iThread, iEntry, fileName); + continue; + } + + // First dereference of 'trackInfo' for this entry -- triggers the actual branch read; must be + // checked before .size() (or anything else) trusts the result, see checkReadOk above. + (void)(**trackRefs); + if (!checkReadOk(*trackRefs, "trackInfo", iEntry)) { + continue; + } + const auto nTracks = (*trackRefs)->size(); + + // Materialize this TF's 'res' branch here, on the producer thread, before it's copied out below. + // TTreeReaderValue::operator* deserializes lazily on the first dereference per entry. + (void)(**unbinnedResiduals).size(); + if (!checkReadOk(*unbinnedResiduals, "res", iEntry)) { + continue; + } + + const double dbgIoMs = std::chrono::duration(std::chrono::steady_clock::now() - dbgTIoStart).count(); + { + thread_local double dbgSumIoMs = 0.0; + thread_local uint64_t dbgNProduced = 0; + dbgSumIoMs += dbgIoMs; + ++dbgNProduced; + if (isAlienFile && dbgNProduced % 10 == 0) { + LOGP(info, "[Thread{}] [prefetch diag][producer] TF {} nTracks={} ioMs={:.1f} | running: {} TF(s) read, sumIoMs={:.0f}", + iThread, iEntry, nTracks, dbgIoMs, dbgNProduced, dbgSumIoMs); + } + } + + // First dereference of 'trackData' for this entry -- same lazy-read/checkReadOk requirement as + // 'trackInfo'/'res' above. + (void)(**trackData); + if (!checkReadOk(*trackData, "trackData", iEntry)) { + continue; + } + + // the track loop below indexes trackData with the trackRefs index, so both have to be in sync + if ((**trackData).size() < nTracks) { + LOGP(warning, "[Thread{}] TF {} of file {} has fewer track data entries than track references ({} vs {}). Skipping TF!", iThread, iEntry, fileName, + (**trackData).size(), nTracks); + continue; + } + + lumiSumCTP_thread[iThread] += (*lumiTF)->getLumi(); + ++lumiEntriesCTP_thread[iThread]; + + timeMSsel_thread[iThread].emplace_back(tfTimeInMS); + orbitsSel_thread[iThread].emplace_back((*lumiTF)->orbit); + ctpLumiSel_thread[iThread].emplace_back((*lumiTF)->getLumi()); + + // Copy the three vectors out (see TFPackage) and hand the package to the consumer. push() blocks + // here if the queue is already at TFQueueDepth, which is exactly the back-pressure that keeps the + // producer from running arbitrarily far ahead. + auto pkg = std::make_unique(); + pkg->iEntry = iEntry; + pkg->trackRefsVec = **trackRefs; + pkg->trackDataVec = **trackData; + pkg->unbinnedResidualsVec = **unbinnedResiduals; + tfQueue.push(std::move(pkg)); + } + tfQueue.setDone(); + }); + + while (true) { + const auto dbgTPopStart = std::chrono::steady_clock::now(); + std::unique_ptr pkg = tfQueue.pop(); + const double dbgPopWaitMs = std::chrono::duration(std::chrono::steady_clock::now() - dbgTPopStart).count(); + const size_t dbgQueueSizeAfterPop = tfQueue.size(); // how many the producer still has ready right now + + // A pop() this slow only ever returns *after* the producer's own push() for this package already + // succeeded -- i.e. the producer is guaranteed to be between TF reads right now, not blocked + // inside one, so signalling it here can never race a live TFile/TGrid call. Real motivation: a + // real GRID job had one file's chunk refill alone take 11+ minutes (~0.7 MB/s vs. 45-100 MB/s for + // every other chunk on the same SE), then a second chunk on the same file never returned and the + // whole job was killed by AliEn's idle-CPU watchdog. This won't catch that second, truly-dead case + // (the producer never reaches this check then -- needs a separate process-level watchdog, not yet + // built), but it does mean a merely very slow file gets abandoned after one bad chunk instead of + // risking a permanent stall. + if (pkg && dbgPopWaitMs > popWaitAbandonMs && !abandonFile.load(std::memory_order_relaxed)) { + LOGP(warning, "[Thread{}] popWait {:.0f} ms exceeds abandon threshold {:.0f} ms (SCDCALIB_POP_WAIT_ABANDON_MS) -- signalling the producer to give up on the rest of this file, assuming a persistently slow fetch rather than a dead one", + iThread, dbgPopWaitMs, popWaitAbandonMs); + abandonFile.store(true, std::memory_order_relaxed); + } + if (!pkg) { + break; // producer is done and the queue is drained + } + const int iEntry = pkg->iEntry; + const auto nTracks = pkg->trackRefsVec.size(); + + ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + // ===| TRACK LOOP |================================================================================================================== + // Body extracted into a per-track lambda so it can be dispatched across nTrackWorkers compute + // threads (see the setup + rationale above the FILE LOOP). A `return` inside this lambda skips to + // the next track -- the lambda body IS one track's worth of work. The cluster loop's own + // `continue`s still mean what they always do: that for-loop lives inside the lambda unchanged. + ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + auto processTrack = [&](size_t iTrack, int iWorker) { + const auto& trkInfo = pkg->trackRefsVec[iTrack]; + if (!GID::includesSource(trkInfo.sourceId, sources)) { + return; + } + const auto& trk = pkg->trackDataVec[iTrack]; + if (!revalidateTrack(trk, params_thread)) { + return; + } + + auto propagator = o2::base::Propagator::Instance(); + o2::track::TrackPar trkPar = trk.par; + int sign = trkPar.getSign(); + + int charge = 0; + if (sign > 0) { + charge = 1; + } + + // dE/dx cut + if (maxdEdx > 0 && trk.dEdxTPC > maxdEdx) { + return; + } + + if (maxdEdxExp > 0 || maxDevdEdxOverExp > 0) { + // propagate to the beginning of the inner containment vessel, to use the momentum for dE/dx expected + if (!propagator->PropagateToXBxByBz(trkPar, 63.2, 0.99, 2., o2::base::Propagator::MatCorrType::USEMatCorrLUT)) { // USEMatCorrTGeo, USEMatCorrLUT, USEMatCorrNONE + return; + } + + const auto dEdxExp = o2::track::BetheBlochSolidOpt(trk.par.getP() / trk.par.getPID().getMass()) * 3e4; + + if (maxdEdxExp > 0 && dEdxExp > maxdEdxExp) { + return; + } + + if (maxDevdEdxOverExp > 0 && std::abs(trk.dEdxTPC / dEdxExp - 1) > maxDevdEdxOverExp) { + return; + } + } + if (!propagator->PropagateToXBxByBz(trkPar, 85., 0.99, 2., o2::base::Propagator::MatCorrType::USEMatCorrLUT)) { // USEMatCorrTGeo, USEMatCorrLUT, USEMatCorrNONE + return; + } + + // INCREASE LOCAL TRACK COUNTER + ++trackCounter_local_worker[iWorker]; + + ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + // ===| CLUSTER & RESIDUAL LOOP |===================================================================================================== + ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + for (unsigned int i = trkInfo.idxFirstResidual; i < trkInfo.idxFirstResidual + trkInfo.nResiduals; ++i) { + const auto& residIn = pkg->unbinnedResidualsVec[i]; + int sec = residIn.sec; + if (residIn.row >= NRows || sec >= NSectors || sec < 0) { // non TPC residuals have row>=160, though to see them one should loop until i<(trc.clIdx.getFirstEntry() + trc.clIdx.getEntries() + trc.nExtDetResid) + continue; + } + + if (isRejectedResidual(residIn)) { + continue; + } + + float angleSec = TMath::DegToRad() * (10.0 + 20.0 * sec); + std::array bvox; + // cluster position + float xPos = param::RowX[residIn.row]; + float yPos = residIn.y * param::MaxY / 0x7fff + residIn.dy * param::MaxResid / 0x7fff; + float zPos = residIn.z * param::MaxZ / 0x7fff + residIn.dz * param::MaxResid / 0x7fff; + // exclude the edge pads as they are biased! + // get max y-position of edge pad: pad centre last pad - pad width/2 + if (skipEdgePads && std::abs(yPos) > yMaxCentrePadByRow[residIn.row]) { + ++nEdgeClustersSkipped_worker[iWorker]; + continue; + } + + clsPosWorker[iWorker].SetXYZ(xPos, yPos, zPos); + if (!trackResiduals.findVoxelBin(sec, xPos, yPos, zPos, bvox)) { + // we are not inside any voxel + continue; + } + + // circle pool for this voxel is shared across threads -> lock for the rest of this iteration + auto& vox = voxels[((sec * NRows + bvox[2]) * nY2XBins + bvox[1]) * nZ2XBins + bvox[0]]; + std::unique_lock voxLock(vox.mtx); + + // XALEX + if (vox.poolCounter[charge] == NPool) { + continue; + } + + //--------------------------------------------------------- + // Main part of the new code + float xposvox, yoverxpos, zoverxpos; + trackResiduals.getVoxelCoordinates(sec, bvox[2], bvox[1], bvox[0], xposvox, yoverxpos, zoverxpos); + if (fabs(xposvox) < 5.0) { + continue; + } + float yposvox = yoverxpos * xposvox; + float zposvox = zoverxpos * xposvox; + + float dxclsvoxel = clsPosWorker[iWorker].X() - xposvox; + float dyclsvoxel = clsPosWorker[iWorker].Y() - yposvox; + float dzclsvoxel = clsPosWorker[iWorker].Z() - zposvox; + + Vec3d deltaClsVoxel; + deltaClsVoxel.SetXYZ(dxclsvoxel, dyclsvoxel, dzclsvoxel); + + //----------------------------------------- + trkPar.rotate(o2::math_utils::sector2Angle(sec)); + + propagator->PropagateToXBxByBz(trkPar, xPos, 0.99, 2., o2::base::Propagator::MatCorrType::USEMatCorrLUT); // USEMatCorrTGeo, USEMatCorrLUT, USEMatCorrNONE + trackPosAtRow_worker[iWorker].SetXYZ(trkPar.getX(), trkPar.getY(), trkPar.getZ()); + + float sna, csa; + o2::math_utils::CircleXY::value_t> xycircle; + trkPar.getCircleParams(magfieldvalue, xycircle, sna, csa); // in global coordinates + + Vec3d circleCenterEstimate; + circleCenterEstimate.SetXYZ(xycircle.xC, xycircle.yC, 0.0); + circleCenterEstimate.RotateZ(-angleSec); + float radius_estimate = xycircle.rC; + + if ((trackPosAtRow_worker[iWorker] - clsPosWorker[iWorker]).Perp() > maxDistIntCls) { + continue; + } + + //----------------------------------------- + // DeltaZ corrections, two methods + if (voxMapInput.size()) // with input map from first itteration, should be more precise than second method + { + // we already have a correction map available + const auto& voxRes = voxelResults[sec][trackResiduals.getGlbVoxBin(bvox)]; // bvox: z,y,x + float DX_input_map = voxRes.D[TrackResiduals::ResX]; + + propagator->PropagateToXBxByBz(trkPar, xPos - DX_input_map, 0.99, 2., o2::base::Propagator::MatCorrType::USEMatCorrLUT); // USEMatCorrTGeo, USEMatCorrLUT, USEMatCorrNONE + float DZdist = (zPos - trkPar.getZ()); // distortion + // accumulate Z sum (shared across threads for this voxel, guarded by vox.mtx) and increment counter + vox.residualsAll[2] += DZdist; + vox.counterZAll++; + } else { + // without input map, use rolling average -- gate on the cumulative NP/PP/NN counters (never + // reset across flushes, unlike poolCounter) so this only fires once residualsAll[0] has + // actually been computed from a decent number of samples, not merely whenever the in-flight + // pool for this charge happens to be non-empty (poolCounter resets to 0 on every flush, so + // checking it here both under- and over-fires relative to residualsAll[0]'s real validity). + if (vox.counterNP > MinDxSamplesForZCorr || vox.counterPP > MinDxSamplesForZCorr || vox.counterNN > MinDxSamplesForZCorr) { + float DXrollingaverage = vox.residualsAll[0]; + propagator->PropagateToXBxByBz(trkPar, xPos - DXrollingaverage, 0.99, 2., o2::base::Propagator::MatCorrType::USEMatCorrLUT); // USEMatCorrTGeo, USEMatCorrLUT, USEMatCorrNONE + float DZdist = (zPos - trkPar.getZ()); // distortion + vox.residualsAll[2] += DZdist; + vox.counterZAll++; + } + } + //----------------------------------------- + + circleCenterEstimate -= deltaClsVoxel; // shift track to voxel center to avoid smearing within voxel + + int counter = vox.poolCounter[charge]; + vox.circleCenters[charge][counter] = Vec3f{static_cast(circleCenterEstimate.X()), static_cast(circleCenterEstimate.Y()), static_cast(circleCenterEstimate.Z())}; + vox.circleRadii[charge][counter] = radius_estimate; + vox.poolCounter[charge]++; + + ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + // ===| PROCESSING VOXEL |============================================================================================================ + ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + if ((vox.poolCounter[0] >= NPool || vox.poolCounter[1] >= NPool)) { + //////////////////////////////// + // Same polarity combinations // + //////////////////////////////// + float weightCurvatureInt = 0.0; + Vec3d averageInt; + averageInt.SetXYZ(0.0, 0.0, 0.0); + + for (int counterPos = 0; counterPos < (vox.poolCounter[1] - 1); counterPos++) { + for (int counterPosB = (counterPos + 1); counterPosB < vox.poolCounter[1]; counterPosB++) { + float radiusA = vox.circleRadii[1][counterPos]; + float radiusB = vox.circleRadii[1][counterPosB]; + + if (fabs((radiusA - radiusB) / (radiusA + radiusB)) < 0.2) { + continue; // to similar bending radii + } + + Vec3d intCircles = getIntCircles(radiusA, radiusB, vox.circleCenters[1][counterPos], vox.circleCenters[1][counterPosB], xposvox, yposvox); + if (intCircles.Perp() < 280.0 && intCircles.Perp() > 60.0) { + float distIntCls = (intCircles - clsPosWorker[iWorker]).Perp(); // why? + if (distIntCls > maxDistIntCls) { + continue; + } + + float weight = fabs((radiusA - radiusB) / (radiusA + radiusB)); + float weight_curvature = weight; + + weightCurvatureInt += weight_curvature; + averageInt += intCircles * weight_curvature; + } + } + } + + if (weightCurvatureInt > 0) { + averageInt *= 1.0 / weightCurvatureInt; + + // distortion + // accumulate PP sums (x,y) and increment counter + vox.residualsPP[0] += (xposvox - averageInt.X()); + vox.residualsPP[1] += (yposvox - averageInt.Y()); + vox.counterPP++; + } + + weightCurvatureInt = 0.0; + averageInt.SetXYZ(0.0, 0.0, 0.0); + + for (int counterNeg = 0; counterNeg < (vox.poolCounter[0] - 1); counterNeg++) { + for (int counterNegB = (counterNeg + 1); counterNegB < vox.poolCounter[0]; counterNegB++) { + float radiusA = vox.circleRadii[0][counterNeg]; + float radiusB = vox.circleRadii[0][counterNegB]; + + if (fabs((radiusA - radiusB) / (radiusA + radiusB)) < 0.2) { + continue; // to similar bending radii + } + + Vec3d intCircles = getIntCircles(radiusA, radiusB, vox.circleCenters[0][counterNeg], vox.circleCenters[0][counterNegB], xposvox, yposvox); + if (intCircles.Perp() < 280.0 && intCircles.Perp() > 60.0) { + float distIntCls = (intCircles - clsPosWorker[iWorker]).Perp(); // why? + if (distIntCls > maxDistIntCls) { + continue; + } + + float weight = fabs((radiusA - radiusB) / (radiusA + radiusB)); + float weight_curvature = weight; + + weightCurvatureInt += weight_curvature; + averageInt += intCircles * weight_curvature; + } + } + } + + if (weightCurvatureInt > 0) { + averageInt *= 1.0 / weightCurvatureInt; + + // distortion + // accumulate NN sums (x,y) and increment counter + vox.residualsNN[0] += (xposvox - averageInt.X()); + vox.residualsNN[1] += (yposvox - averageInt.Y()); + vox.counterNN++; + } + + //////////////////////////////////// + // Opposite polarity combinations // + //////////////////////////////////// + weightCurvatureInt = 0.0; + averageInt.SetXYZ(0.0, 0.0, 0.0); + + for (int counterPos = 0; counterPos < vox.poolCounter[1]; counterPos++) { + for (int counterNeg = 0; counterNeg < vox.poolCounter[0]; counterNeg++) { + float radiusA = vox.circleRadii[1][counterPos]; + float radiusB = vox.circleRadii[0][counterNeg]; + Vec3d intCircles = getIntCircles(radiusA, radiusB, vox.circleCenters[1][counterPos], vox.circleCenters[0][counterNeg], xposvox, yposvox); + if (intCircles.Perp() < 280.0 && intCircles.Perp() > 60.0) { + float distIntCls = (intCircles - clsPosWorker[iWorker]).Perp(); // why? + if (distIntCls > maxDistIntCls) { + continue; + } + + // float weight_curvature = (1.0/radiusA)*(1.0/radiusB); // the larger the curvature the more precise the intersection can be calculated + // float weight_curvature = (radiusA)*(radiusB); // the larger the curvature the more precise the intersection can be calculated + + float weight_curvature = 1.0; // (1.0/radiusA)*(1.0/radiusB); + + weightCurvatureInt += weight_curvature; + averageInt += intCircles * weight_curvature; + } + } + } + + if (weightCurvatureInt > 0) { + averageInt *= 1.0 / weightCurvatureInt; + + // distortion + // accumulate NP sums (x,y) and increment counter + vox.residualsNP[0] += (xposvox - averageInt.X()); + vox.residualsNP[1] += (yposvox - averageInt.Y()); + vox.counterNP++; + } + + float weightNP = vox.counterNP * 1.0; + float weightPP = vox.counterPP * 0.01; + float weightNN = vox.counterNN * 0.01; + + float sum_weight = weightNP + weightPP + weightNN; + + if (sum_weight > 0.0f) { + // convert sums to means before weighting + float meanNPx = (vox.counterNP > 0) ? (vox.residualsNP[0] / static_cast(vox.counterNP)) : 0.0f; + float meanNPy = (vox.counterNP > 0) ? (vox.residualsNP[1] / static_cast(vox.counterNP)) : 0.0f; + float meanPPx = (vox.counterPP > 0) ? (vox.residualsPP[0] / static_cast(vox.counterPP)) : 0.0f; + float meanPPy = (vox.counterPP > 0) ? (vox.residualsPP[1] / static_cast(vox.counterPP)) : 0.0f; + float meanNNx = (vox.counterNN > 0) ? (vox.residualsNN[0] / static_cast(vox.counterNN)) : 0.0f; + float meanNNy = (vox.counterNN > 0) ? (vox.residualsNN[1] / static_cast(vox.counterNN)) : 0.0f; + + vox.residualsAll[0] = (weightNP * meanNPx + weightPP * meanPPx + weightNN * meanNNx) / sum_weight; + vox.residualsAll[1] = (weightNP * meanNPy + weightPP * meanPPy + weightNN * meanNNy) / sum_weight; + } + + // reset the pools: only the counter is reset. circleCenters/circleRadii are fixed-size + // std::array now (not std::vector) -- there's no clear()/resize() to even call; + // stale entries beyond the new poolCounter are simply overwritten as the pool refills. + for (int ic = 0; ic < 2; ic++) { + vox.poolCounter[ic] = 0; + } + } + //--------------------------------------------------------- + + // // update COG for voxel bvox (update for X only needed in case binning is not per pad row) + + } // end of cluster loop + }; // end of processTrack lambda + + const auto dbgTCpuStart = std::chrono::steady_clock::now(); + if (nTrackWorkers <= 1) { + for (size_t iTrack = 0; iTrack < nTracks; ++iTrack) { + processTrack(iTrack, 0); + } + } else { + std::vector trackThreads; + trackThreads.reserve(nTrackWorkers); + for (int iWorker = 0; iWorker < nTrackWorkers; ++iWorker) { + trackThreads.emplace_back([&, iWorker]() { + for (size_t iTrack = iWorker; iTrack < nTracks; iTrack += nTrackWorkers) { + processTrack(iTrack, iWorker); + } + }); + } + for (auto& th : trackThreads) { + th.join(); + } + } + const double dbgCpuMs = std::chrono::duration(std::chrono::steady_clock::now() - dbgTCpuStart).count(); + { + // [prefetch diag][consumer]: the number that actually matters. dbgPopWaitMs is how long this + // thread blocked in tfQueue.pop() waiting for the producer -- with the pipeline working, this + // should be near zero most of the time (the producer stays ahead), spiking only when it can't + // keep up (e.g. a cache-refill spike bigger than TFQueueDepth's cushion). + thread_local double dbgSumPopWaitMs = 0.0; + thread_local double dbgSumCpuMs = 0.0; + thread_local uint64_t dbgNTFs = 0; + dbgSumPopWaitMs += dbgPopWaitMs; + dbgSumCpuMs += dbgCpuMs; + ++dbgNTFs; + // Unlike the producer log above, this one had no sampling gate at all -- printed every single + // TF, on every thread, which is the dominant source of "[prefetch diag]" log volume. Matched to + // the producer's every-10th cadence (the cumulative sum* fields lose nothing from sampling) and + // gated on isAlienFile for the same reason as the producer log. + if (isAlienFile && dbgNTFs % 10 == 0) { + LOGP(info, "[Thread{}] [prefetch diag][consumer] TF {} nTracks={} popWaitMs={:.1f} cpuMs={:.1f} queueSizeAfterPop={} | running totals over {} TF(s): sumPopWaitMs={:.0f} sumCpuMs={:.0f}", + iThread, iEntry, nTracks, dbgPopWaitMs, dbgCpuMs, dbgQueueSizeAfterPop, dbgNTFs, dbgSumPopWaitMs, dbgSumCpuMs); + } + } + // Merge this TF's per-worker counters into the file-thread-level counters the rest of + // doFileProcessing (and the caller's merge loop, for nEdgeClustersSkipped_thread) expect. + for (int iWorker = 0; iWorker < nTrackWorkers; ++iWorker) { + trackCounter_local += trackCounter_local_worker[iWorker]; + trackCounter_local_worker[iWorker] = 0; + nEdgeClustersSkipped_thread[iThread] += nEdgeClustersSkipped_worker[iWorker]; + nEdgeClustersSkipped_worker[iWorker] = 0; + } + } // end of TF loop (consumer) + + // producerThread only ever reaches here via tfQueue.setDone() (end of file or maxTracks reached), + // which the consumer's pop() == nullptr check above already waited for -- this join is therefore + // just reclaiming the thread, not a real wait. Must happen before perfStats->Finish()/Print() below, + // since perfStats is only ever touched by the producer thread. + producerThread.join(); + + if (perfStats) { + perfStats->Finish(); + if (isAlienFile) { + perfStats->Print(); // full TTreePerfStats I/O report (incl. "Disk IO = ... MBytes/s") -- GRID-diagnostic only, see isAlienFile above + } + totalBytesReadPerf_thread[iThread] += perfStats->GetBytesRead(); + } + + // Occasionally flush local count to global + if (trackCounter_local >= 1000) { + nTracksProcessed.fetch_add(trackCounter_local, std::memory_order_relaxed); + trackCounter_local = 0; + } + + } // end of file loop + + if (trackCounter_local > 0) { + nTracksProcessed.fetch_add(trackCounter_local, std::memory_order_relaxed); + } +} + +void staticMapCreatorCPM(std::string fileInput = "files.txt", + int runNumber = 527976, + std::string fileOutput = "voxRes.root", + std::string voxMapInput = "", + std::string trackSources = static_cast(GID::ALL), + std::string z2xBinning = "", // empty: default binning; single number: uniform binning; otherwise bin boundaries, e.g. "0.,0.02, 0.04, 0.06, 1" + std::string y2xBinning = "", // empty: default binning; single number: uniform binning; otherwise bin boundaries, e.g. "-1,-0.998, -0.996, ... 0.996, 0.998, 1" + bool useSmoothed = true, // use smoothed residuals as input + bool createSpline = true, // create the splines + int maxTracksPerSlice = -1, // limit the number of total tracks processed to maxTracksPerSlice * nBinsZ2X * nBinxY2X * 36 + int minTracksPerSlice = -1, // request a minimum number of tracks per slice. Otherwise the calibration is not created + std::string badRangeList = "", // list of bad time ranges to be excluded in the calibration + long firstTFTime = -1, // First TF time to accept + long lastTFTime = -1, // Last TF time to accept + float maxdEdx = -1, // dE/dx cut above which tracks will be rejected + float maxdEdxExp = -1, // dE/dx expected cut above which tracks will be rejected + float maxDevdEdxOverExp = -1, // maximum deviation of dE/dx / expected value, above the track is rejected + bool skipEdgePads = 1, // skip edge pads in the calibration, by default on + std::string badRangeSelection = "ALL", // use bad time ranges only for specific comment e.g. C1 + float maxZ2XCut = 1.f, // overrides scdcalib.maxZ2X (the track tgl cut applied in revalidateTrack). + // Defaults to the SpacePointsCalibConfParam compiled-in value of 1.0 + int maxTrackWorkers = -1, // number of threads used for the track loop within one file. <=0 (default) + // auto-detects via hardware_concurrency(), which reports the machine's core + // count rather than the cores actually allocated to a batch job -- pass the + // real allocation explicitly when running on a batch system or the GRID + int nThreads = 8) // number of file-level threads, i.e. how many input files are processed + // concurrently. Forced to 1 for alien:// input regardless, since TGrid is + // not thread-safe (see below) +{ + LOGP(info, "TrackData::filterFlag {}, UnbinnedResid::rejected {}", + HasFilterFlagMember::value ? "available" : "NOT available (old O2, cut skipped)", + HasRejectedMember::value ? "available" : "NOT available (old O2, cut skipped)"); + + // Enable multiple threads + ROOT::EnableThreadSafety(); + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + // ==== TIME =============================================================================================================== + auto t_start = std::chrono::high_resolution_clock::now(); + // ========================================================================================================================= + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + fair::Logger::SetVerbosity(fair::Verbosity::medium); + fair::Logger::SetConsoleSeverity(fair::Severity::info); + fair::Logger::SetFileSeverity(fair::Severity::info); + + const Mapper& mapper = Mapper::instance(); + + // Obtain configuration + const SpacePointsCalibConfParam& params = SpacePointsCalibConfParam::Instance(); + if (!std::filesystem::exists("scdconfig.ini")) { + LOGP(warning, "Did not find configuration file. Using default parameters and storing them in scdconfig.ini"); + params.writeINI("scdconfig.ini", "scdcalib"); // to write default parameters to a file + } else { + params.updateFromFile("scdconfig.ini"); + } + // Explicit override from the macro's own argument, applied AFTER any scdconfig.ini load so it wins + // regardless of whether one happens to exist in CWD -- see maxZ2XCut's doc comment above. Using the + // string-based setValue overload (implemented out-of-line in ConfigurableParam.cxx), not the + // templated one -- that one needs boost::property_tree fully instantiated at the call site, which + // isn't included here. + o2::conf::ConfigurableParam::setValue("scdcalib.maxZ2X", std::to_string(maxZ2XCut)); + LOGP(info, "----- Dumping configuration values START -----"); + params.printKeyValues(); + LOGP(info, "----- Dumping configuration values END -----"); + + GID::mask_t allowedSources = GID::getSourcesMask("ITS-TPC,ITS-TPC-TRD,ITS-TPC-TOF,ITS-TPC-TRD-TOF"); + GID::mask_t sources = allowedSources & GID::getSourcesMask(trackSources); + + // Get CCDB objects + auto& ccdbmgr = o2::ccdb::BasicCCDBManager::instance(); + ccdbmgr.setCaching(true); + ccdbmgr.setFatalWhenNull(false); + ccdbmgr.setURL("http://alice-ccdb.cern.ch"); + auto runDuration = ccdbmgr.getRunDuration(runNumber); + auto tRun = runDuration.first + (runDuration.second - runDuration.first) / 2; // time stamp for the middle of the run duration + ccdbmgr.setTimestamp(tRun); + + // CTP orbit reset time, does not change during the run + const auto orbitResetTimeNS = ccdbmgr.get>("CTP/Calib/OrbitReset"); + const int64_t orbitResetTimeMS = (*orbitResetTimeNS)[0] * 1e-3; + LOGP(info, "Orbit reset time in MS is {}", orbitResetTimeMS); + + // Geometry, material budget and B-field + auto geoAligned = ccdbmgr.get("GLO/Config/GeometryAligned"); + auto magField = ccdbmgr.get("GLO/Config/GRPMagField"); + const o2::base::MatLayerCylSet* matLut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdbmgr.get("GLO/Param/MatLUT")); + o2::base::Propagator::initFieldFromGRP(magField); + auto prop = o2::base::Propagator::Instance(); + prop->setMatLUT(matLut); + float magfieldvalue = static_cast(magField->getNominalL3Field()); + LOGP(info, "Nominal L3 field: {:.3f}", magfieldvalue); + + // GRP LHC and beam type + auto grplhc = ccdbmgr.get("GLO/Config/GRPLHCIF"); + const auto beamA = grplhc->getBeamZ(o2::constants::lhc::BeamA); + const auto beamC = grplhc->getBeamZ(o2::constants::lhc::BeamC); + const auto eCM = grplhc->getSqrtS(); + bool isPbPb = (beamA == 82 && beamC == 82); + LOGP(info, "BeamA: {}, BeamC: {}, isPbPb: {}, Ecm: {}", beamA, beamC, isPbPb, eCM); + + // Input + auto fileList = getInputFileList(fileInput); + + // Single-threaded when reading from AliEn: TGrid/xrootd access is not thread-safe. This overrides the + // nThreads argument. Detected from the actual input file entries rather than from fileInput itself, + // which for a GRID job is a local list of alien:// entries, so that local/Lustre input still gets the + // requested parallelism. + if (!fileList.empty() && fileList[0].rfind("alien://", 0) == 0) { + LOGP(info, "AliEn input detected (alien:// prefix) -- TGrid access is not thread-safe, forcing nThreads=1"); + nThreads = 1; + } + int nFileThreads = nThreads; + LOGP(info, "Using {} threads for processing", nFileThreads); + + // Set up binning + const auto z2xBins = o2::RangeTokenizer::tokenize(z2xBinning); + const auto y2xBins = o2::RangeTokenizer::tokenize(y2xBinning); + + TrackResiduals trackResiduals; + + trackResiduals.setZ2XBinning(z2xBins); + trackResiduals.setY2XBinning(y2xBins); + trackResiduals.init(); + + const int nY2XBins = trackResiduals.getNY2XBins(); + const int nZ2XBins = trackResiduals.getNZ2XBins(); + LOGP(info, "nY2XBins: {}, nZ2XBins: {}", nY2XBins, nZ2XBins); + + const float maxDistIntCls = 25.0; + + std::vector>>>> vec_residualsAll; // sector, voxX, voxF, voxZ, xyz + std::vector>>> vec_residuals_counterAll; // sector, voxX, voxF, voxZ + + vec_residualsAll.resize(NSectors); + vec_residuals_counterAll.resize(NSectors); + for (int isec = 0; isec < NSectors; isec++) { + vec_residualsAll[isec].resize(NRows); + vec_residuals_counterAll[isec].resize(NRows); + for (int ix = 0; ix < NRows; ix++) { + vec_residualsAll[isec][ix].resize(nY2XBins); + vec_residuals_counterAll[isec][ix].resize(nY2XBins); + for (int iy = 0; iy < nY2XBins; iy++) { + vec_residualsAll[isec][ix][iy].resize(nZ2XBins); + vec_residuals_counterAll[isec][ix][iy].resize(nZ2XBins); + for (int iz = 0; iz < nZ2XBins; iz++) { + vec_residualsAll[isec][ix][iy][iz].resize(3); + for (int ixyz = 0; ixyz < 3; ixyz++) { + vec_residualsAll[isec][ix][iy][iz][ixyz] = 0.0; + } + vec_residuals_counterAll[isec][ix][iy][iz] = 0; + } + } + } + } + + const int nSlices = NSectors * trackResiduals.getNY2XBins() * trackResiduals.getNZ2XBins(); + const Long64_t maxTracks = maxTracksPerSlice * nSlices; + const int minTracks = minTracksPerSlice * nSlices; + std::atomic nTracksProcessed{0}; + Long64_t nTracksProcessed_final{0}; + + // Do we have a correction map available that we should apply to the clusters before the map extraction? + std::array, NSectors> voxelResults{}; + if (voxMapInput.size()) { + LOGP(info, "[InputCorrMap]: A correction map has been provided. Will apply the corrections to the cluster residuals"); + LOGP(info, "[InputCorrMap]: Resizing voxelResults to number of voxels"); + for (int iSec = 0; iSec < NSectors; ++iSec) { + voxelResults[iSec].resize(trackResiduals.getNVoxelsPerSector()); + } + TrackResiduals::VoxRes* voxResPtr = nullptr; + std::unique_ptr fIn = std::make_unique(voxMapInput.c_str()); + if (!fIn->IsOpen() || fIn->IsZombie()) { + LOGP(fatal, "[InputCorrMap]: Could not open input file {}", voxMapInput); + } + LOGP(info, "[InputCorrMap]: Getting TTree of voxels"); + std::unique_ptr treeIn; + treeIn.reset((TTree*)fIn->Get("voxResTree")); + if (!treeIn) { + LOGP(fatal, "[InputCorrMap]: Could not extract voxResTree from input file {}", voxMapInput); + } + treeIn->SetBranchAddress("voxRes", &voxResPtr); + LOGP(info, "[InputCorrMap]: Getting voxel results and filling voxelResults"); + for (int iEntry = 0; iEntry < treeIn->GetEntries(); ++iEntry) { + treeIn->GetEntry(iEntry); + auto& voxRes = *voxResPtr; + voxelResults[voxRes.bsec][trackResiduals.getGlbVoxBin(voxRes.bvox)] = voxRes; + } + } + // voxelResults is only read inside doFileProcessing and is not touched between thread start and join, + // so it is shared by const reference instead of copied once per thread (copying would cost + // nFileThreads x 36 x nVoxPerSector VoxRes objects). + LOGP(info, "Sharing the provided input Map with all threads"); + + // Check for bad ranges + std::vector badRanges; + std::vector> badRanges_thread(nFileThreads); + bool invertBadRange = false; + if (badRangeList.length() > 0) { + if (badRangeList[0] == '-') { + LOGP(info, "Inverting badRange list!"); + invertBadRange = true; + badRangeList.erase(0, 1); + } + badRanges = loadRunTimeSpans(badRangeList, runNumber, badRangeSelection); + } + for (int iThread = 0; iThread < nFileThreads; ++iThread) { + badRanges_thread[iThread] = badRanges; + } + + // ---| Lumi estimators |--- + size_t lumiEntriesCTP = 0; + double lumiSumCTP = 0; + + // vector of selected values + std::vector orbitsSel; + // This macro does not read IDC scalers from CCDB (see the "IDC values" note above the OrbitLumiInfo + // tree below). These two stay empty and are written out only to keep the output format stable for + // readers that expect the branches; the values are filled in offline from timeMSsel. + std::vector idcScalerASel; + std::vector idcScalerCSel; + std::vector ctpLumiSel; + std::vector timeMSsel; // time in ms collected over all selected TFs + + //---------------------------- + const std::filesystem::path pFileOutput(fileOutput); + std::string outPath(pFileOutput.parent_path().c_str()); + if (outPath.empty()) { + outPath = "."; + } + + fair::Logger::SetConsoleSeverity(fair::Severity::error); + /////////////////////////// + // Create thread vectors // + /////////////////////////// + // trackResiduals (declared earlier, already configured) is shared read-only across threads -- no per-thread copy needed + + long int nEdgeClustersSkipped{0}; + std::vector nEdgeClustersSkipped_thread(nFileThreads, 0); // Does NEED to be merged + long int nTracksSkippedByBadRangeList{0}; + std::vector nTracksSkippedByBadRangeList_thread(nFileThreads, 0); // Does NEED to be merged + long int nTFs{0}; + long int nTFsSkippedByBadRangeList{0}; + long int nTFsSkippedByTimeWindow{0}; + std::vector nTFs_thread(nFileThreads, 0); // Does NEED to be merged + std::vector nTFsSkippedByBadRangeList_thread(nFileThreads, 0); // Does NEED to be merged + std::vector nTFsSkippedByTimeWindow_thread(nFileThreads, 0); // Does NEED to be merged + + std::vector voxels(NSectors * NRows * nY2XBins * nZ2XBins); // flat [sec][ix][iy][iz] -- shared circle pool, ONE instance for all threads (guarded by per-voxel mutex); sized directly since VoxelData (holds a mutex) cannot be resize()'d after construction // Does NOT need to be merged + + // residualsAll/NP/PP/NN + their counters now live directly in VoxelData (shared, per-voxel mutex), so no per-thread copies or later merge pass are needed for them. + + // The per-file input handles (TFile, TTreeReaders, TTreeReaderValues) and the I/O-rate sampling state + // are plain locals inside doFileProcessing: each thread only ever used its own slot, so they never + // needed to be shared vectors here. Only totalBytesReadPerf is still per-thread, because it is summed + // across threads after the join below. + std::vector totalBytesReadPerf_thread(nFileThreads, 0); // Does NEED to be merged + + std::vector lumiEntriesCTP_thread(nFileThreads, 0); // Does NEED to be merged + std::vector lumiSumCTP_thread(nFileThreads, 0); // Does NEED to be merged + + // vector of selected values + std::vector> orbitsSel_thread(nFileThreads); // Does NEED to be merged + std::vector> ctpLumiSel_thread(nFileThreads); // Does NEED to be merged + std::vector> timeMSsel_thread(nFileThreads); // time in ms collected over all selected TFs // Does NEED to be merged + + fair::Logger::SetConsoleSeverity(fair::Severity::info); + printMemoryUsage("Memory usage after init buffers"); + + // Start threads + std::vector threads(nFileThreads); + for (int i = 0; i < nFileThreads; i++) { + threads[i] = std::thread(doFileProcessing, + i, + nFileThreads, + maxTrackWorkers, + firstTFTime, + lastTFTime, + invertBadRange, + maxdEdx, + maxdEdxExp, + maxDevdEdxOverExp, + skipEdgePads, + std::ref(nEdgeClustersSkipped_thread), + std::ref(nTracksSkippedByBadRangeList_thread), + std::ref(nTFs_thread), + std::ref(nTFsSkippedByBadRangeList_thread), + std::ref(nTFsSkippedByTimeWindow_thread), + voxMapInput, + sources, + orbitResetTimeMS, + magfieldvalue, + fileList, + maxTracksPerSlice, + maxTracks, + std::ref(nTracksProcessed), + std::cref(voxelResults), + std::ref(badRanges_thread), + std::cref(trackResiduals), + maxDistIntCls, + nY2XBins, + nZ2XBins, + std::ref(voxels), + std::ref(totalBytesReadPerf_thread), + std::ref(lumiEntriesCTP_thread), + std::ref(lumiSumCTP_thread), + std::ref(orbitsSel_thread), + std::ref(ctpLumiSel_thread), + std::ref(timeMSsel_thread)); + } + + // Wait for the threads to finish + for (auto& th : threads) { + th.join(); + } + + ////////////////////////////////////////////////////////////////////// + // ===| CODE TO MERGE VECTORS |======================================= + ////////////////////////////////////////////////////////////////////// + // START OF MERGE + nTracksProcessed_final = nTracksProcessed.load(); + Long64_t totalBytesReadPerf = 0; // sum of TTreePerfStats bytes read across all threads/files (network volume) + for (int iThread = 0; iThread < nFileThreads; ++iThread) { + // Merge counters + totalBytesReadPerf += totalBytesReadPerf_thread[iThread]; + lumiEntriesCTP += lumiEntriesCTP_thread[iThread]; + lumiSumCTP += lumiSumCTP_thread[iThread]; + nEdgeClustersSkipped += nEdgeClustersSkipped_thread[iThread]; + nTracksSkippedByBadRangeList += nTracksSkippedByBadRangeList_thread[iThread]; + nTFs += nTFs_thread[iThread]; + nTFsSkippedByBadRangeList += nTFsSkippedByBadRangeList_thread[iThread]; + nTFsSkippedByTimeWindow += nTFsSkippedByTimeWindow_thread[iThread]; + + // Merge vector data + orbitsSel.insert(orbitsSel.end(), + orbitsSel_thread[iThread].begin(), + orbitsSel_thread[iThread].end()); + + ctpLumiSel.insert(ctpLumiSel.end(), + ctpLumiSel_thread[iThread].begin(), + ctpLumiSel_thread[iThread].end()); + + timeMSsel.insert(timeMSsel.end(), + timeMSsel_thread[iThread].begin(), + timeMSsel_thread[iThread].end()); + } + + { + const double totalMBReadPerf = totalBytesReadPerf / (1024.0 * 1024.0); + LOGP(info, "All threads done | read {:.1f} MB total (unbinnedResid, across {} thread(s))", totalMBReadPerf, nFileThreads); + } + + // Compute final binned residuals directly from the shared per-voxel accumulators in `voxels` + // (each thread already accumulated straight into the single shared VoxelData under vox.mtx, + // so there is no per-thread data left to sum over here). + for (int isec = 0; isec < NSectors; isec++) { + for (int iz = 0; iz < nZ2XBins; iz++) { + for (int iy = 0; iy < nY2XBins; iy++) { + for (int ix = 0; ix < NRows; ix++) { + const auto& vox = voxels[((isec * NRows + ix) * nY2XBins + iy) * nZ2XBins + iz]; + + // weights as in original code + double weightNP = vox.counterNP * 1.0; + double weightPP = vox.counterPP * 0.01; + double weightNN = vox.counterNN * 0.01; + double sum_weight = weightNP + weightPP + weightNN; + + // For ixyz == 0,1: weighted average as before + for (int ixyz = 0; ixyz < 2; ++ixyz) { + if (sum_weight > 0.0) { + double meanNP = (vox.counterNP > 0) ? (vox.residualsNP[ixyz] / vox.counterNP) : 0.0; + double meanPP = (vox.counterPP > 0) ? (vox.residualsPP[ixyz] / vox.counterPP) : 0.0; + double meanNN = (vox.counterNN > 0) ? (vox.residualsNN[ixyz] / vox.counterNN) : 0.0; + vec_residualsAll[isec][ix][iy][iz][ixyz] = static_cast((weightNP * meanNP + weightPP * meanPP + weightNN * meanNN) / sum_weight); + } else { + vec_residualsAll[isec][ix][iy][iz][ixyz] = 0.0f; + } + } + // For ixyz == 2: simple mean; vox.residualsAll[2] is the running sum of DZdist, vox.counterZAll the count + vec_residualsAll[isec][ix][iy][iz][2] = (vox.counterZAll > 0) ? static_cast(vox.residualsAll[2] / vox.counterZAll) : 0.0f; + vec_residuals_counterAll[isec][ix][iy][iz] = vox.counterNP + vox.counterPP + vox.counterNN; + } + } + } + } + // END OF MERGE + ////////////////////////////////////////////////////////////////////// + + bool isBadCalib = false; + if ((minTracksPerSlice > 0) && (nTracksProcessed_final < minTracks)) { + LOGP(error, "Processed tracks: {} ({}), max requested tracks: {} ({}), minimum number of tracks not reached {} ({}), calibration will be marked as bad, skipped {} edge clusters, skipped tracks by badRangeList {} ({}), skipped TFs outside requested time window {} of {} ({:.1f}%)", nTracksProcessed_final, nTracksProcessed_final / nSlices, maxTracks, maxTracksPerSlice, minTracks, minTracksPerSlice, nEdgeClustersSkipped, nTracksSkippedByBadRangeList, nTracksSkippedByBadRangeList / nSlices, nTFsSkippedByTimeWindow, nTFs, (nTFs > 0) ? (100.0 * nTFsSkippedByTimeWindow / nTFs) : 0.0); + LOGP(info, "Processed time: {}ms, skipped time by badRangeList {}ms ({})", nTFs * o2::constants::lhc::LHCOrbitMUS * 1.e-3 * 32, nTFsSkippedByBadRangeList * o2::constants::lhc::LHCOrbitMUS * 1.e-3 * 32, float(nTFsSkippedByBadRangeList) / float(nTFs)); + const std::string stem = fs::path(fileOutput.data()).stem().c_str(); + std::ofstream(fmt::format("badCalib.{}", stem)).close(); + isBadCalib = true; + } else { + LOGP(info, "Processed tracks: {} ({}), max requested tracks: {} ({}), skipped {} edge clusters, skipped tracks by badRangeList {} ({}), skipped TFs outside requested time window {} of {} ({:.1f}%)", nTracksProcessed_final, nTracksProcessed_final / nSlices, maxTracks, maxTracksPerSlice, nEdgeClustersSkipped, nTracksSkippedByBadRangeList, nTracksSkippedByBadRangeList / nSlices, nTFsSkippedByTimeWindow, nTFs, (nTFs > 0) ? (100.0 * nTFsSkippedByTimeWindow / nTFs) : 0.0); + LOGP(info, "Processed time: {}ms, skipped time by badRangeList {}ms ({})", nTFs * o2::constants::lhc::LHCOrbitMUS * 1.e-3 * 32, nTFsSkippedByBadRangeList * o2::constants::lhc::LHCOrbitMUS * 1.e-3 * 32, float(nTFsSkippedByBadRangeList) / float(nTFs)); + } + //---------------------------------------------------------------- + + // IDC values: this macro deliberately does not query the TPC scalers from CCDB. Doing so per TF from + // inside the processing loop is slow and unreliable (uncached fetch, object reload on every jump in + // time), so meanIDC/medianIDC are written as 0 placeholders here. The per-TF timestamps needed to + // recover them are written to the OrbitLumiInfo tree below (timeMSsel), so the real values can be + // joined in afterwards, offline, against a properly cached CCDB client. + float meanIDC = 0.f; // not const: written to a TTree branch below, which needs a mutable address + + double meanCTP = 0; + if (lumiEntriesCTP > 0) { + meanCTP = lumiSumCTP / lumiEntriesCTP * (isPbPb ? 2.414 : 1); // 2.414 for PbPb + } + + TFile* outputfile = new TFile(fileOutput.c_str(), "RECREATE"); + LOGP(info, "Output file: {} created", fileOutput); + + o2::tpc::TrackResiduals::VoxRes mVoxelResultsOut{}; ///< the results from mVoxelResults are copied in here to be able to stream them + o2::tpc::TrackResiduals::VoxRes* mVoxelResultsOutPtr{&mVoxelResultsOut}; ///< pointer to set the branch address to for the output + std::unique_ptr mTreeOut; + + // Same tree-alias set as the real TrackResiduals::createOutputFile() (SpacePoints/TrackResiduals.cxx), + // so this tree can be TTree::Draw()'n the same way downstream regardless of which of the two wrote it. + if (trackResiduals.getNVoxelsPerSector() == 0) { + LOGP(warning, "For the tree aliases to work you must initialize the binning before calling createOutputFile()"); + } + mTreeOut = std::make_unique("voxResTree", "voxRes map results and statistics"); + mTreeOut->SetAlias("z2xBin", "bvox[0]"); + mTreeOut->SetAlias("y2xBin", "bvox[1]"); + mTreeOut->SetAlias("xBin", "bvox[2]"); + mTreeOut->SetAlias("z2xAV", "stat[0]"); + mTreeOut->SetAlias("y2xAV", "stat[1]"); + mTreeOut->SetAlias("xAV", "stat[2]"); + mTreeOut->SetAlias("fsector", "bsec+0.5+9.*(y2xAV)/pi"); + mTreeOut->SetAlias("phi", "(bsec%18+0.5+9.*(stat[1])/pi)/9*pi"); + mTreeOut->SetAlias("r", "stat[2]"); + mTreeOut->SetAlias("z", "z2xAV*xAV"); + mTreeOut->SetAlias("dX", "D[0]"); + mTreeOut->SetAlias("dY", "D[1]"); + mTreeOut->SetAlias("dZ", "D[2]"); + mTreeOut->SetAlias("dXS", "DS[0]"); + mTreeOut->SetAlias("dYS", "DS[1]"); + mTreeOut->SetAlias("dZS", "DS[2]"); + mTreeOut->SetAlias("dXE", "E[0]"); + mTreeOut->SetAlias("dYE", "E[1]"); + mTreeOut->SetAlias("dZE", "E[2]"); + mTreeOut->SetAlias("voxelIndex", Form("xBin + %i * (y2xBin + %i * z2xBin) + %i * bsec", trackResiduals.getNXBins(), trackResiduals.getNY2XBins(), trackResiduals.getNVoxelsPerSector())); + mTreeOut->SetAlias("entries", "stat[3]"); + mTreeOut->SetAlias("fitOK", Form("(flags & %u) == %u", TrackResiduals::DistDone, TrackResiduals::DistDone)); + mTreeOut->SetAlias("dispOK", Form("(flags & %u) == %u", TrackResiduals::DispDone, TrackResiduals::DispDone)); + mTreeOut->SetAlias("smtOK", Form("(flags & %u) == %u", TrackResiduals::SmoothDone, TrackResiduals::SmoothDone)); + mTreeOut->SetAlias("masked", Form("(flags & %u) == %u", TrackResiduals::Masked, TrackResiduals::Masked)); + mTreeOut->Branch("voxRes", &mVoxelResultsOutPtr); + + // Placeholder, filled in offline together with meanIDC -- see the note above. + float medianIDC = 0.f; + float medianCTP = static_cast(calculateMedian(ctpLumiSel)); + long medianTimeMS = static_cast(calculateMedian(timeMSsel)); + long meanTimeMS = static_cast(calculateMean(timeMSsel)); + + auto userInfo = mTreeOut->GetUserInfo(); + userInfo->Add(new TNamed("meanIDC", std::to_string(meanIDC).data())); + userInfo->Add(new TNamed("meanCTP", std::to_string(meanCTP).data())); + userInfo->Add(new TNamed("meanTimeMS", std::to_string(meanTimeMS).data())); + userInfo->Add(new TNamed("medianIDC", std::to_string(medianIDC).data())); + userInfo->Add(new TNamed("medianCTP", std::to_string(medianCTP).data())); + userInfo->Add(new TNamed("medianTimeMS", std::to_string(medianTimeMS).data())); + userInfo->Add(new TNamed("y2xBinning", y2xBinning.data())); + userInfo->Add(new TNamed("z2xBinning", z2xBinning.data())); + // TrackResiduals::setZ2XBinning() scales the physical z/x bin boundaries by scdcalib.maxZ2X -- this + // value is baked into what each z2x voxel index means in THIS tree, not just a cosmetic config knob. + // Stage 2 has no scdconfig.ini on the GRID and would otherwise silently reconstruct the binning with + // the O2 code default (1.0) instead of maxZ2XCut (1.4 in production), misaligning its geometry against + // the one this tree's voxels were actually filled with. Stored here so stage 2 can apply the exact + // same value it was built with, not a separately-configured guess. + userInfo->Add(new TNamed("maxZ2X", std::to_string(maxZ2XCut).data())); + userInfo->Add(new TNamed("nSlicesPhiZ", std::to_string(nSlices))); + userInfo->Add(new TNamed("maxTracks", std::to_string(maxTracks))); + userInfo->Add(new TNamed("minTracks", std::to_string(minTracks))); + userInfo->Add(new TNamed("nTracksProcessed", std::to_string(nTracksProcessed_final))); + if (isBadCalib) { + userInfo->Add(new TNamed("badCalib", "1")); + } + + for (int isec = 0; isec < NSectors; isec++) { + for (int iz = 0; iz < nZ2XBins; iz++) { + for (int iy = 0; iy < nY2XBins; iy++) { + for (int ix = 0; ix < NRows; ix++) { + for (int ixyz = 0; ixyz < 3; ixyz++) { + mVoxelResultsOut.D[ixyz] = vec_residualsAll[isec][ix][iy][iz][ixyz]; + mVoxelResultsOut.DS[ixyz] = vec_residualsAll[isec][ix][iy][iz][ixyz]; + mVoxelResultsOut.DC[ixyz] = vec_residualsAll[isec][ix][iy][iz][ixyz]; + mVoxelResultsOut.E[ixyz] = 0.1; + } + + float xposvox, yoverxpos, zoverxpos; + trackResiduals.getVoxelCoordinates(isec, ix, iy, iz, xposvox, yoverxpos, zoverxpos); + + mVoxelResultsOut.stat[0] = static_cast(zoverxpos); // z/x, y/x, x, entries + mVoxelResultsOut.stat[1] = static_cast(yoverxpos); + mVoxelResultsOut.stat[2] = static_cast(xposvox); + mVoxelResultsOut.stat[3] = static_cast(vec_residuals_counterAll[isec][ix][iy][iz]); // number of entries used + + mVoxelResultsOut.EXYCorr = 1.0; + mVoxelResultsOut.dYSigMAD = 1.0; + mVoxelResultsOut.dZSigLTM = 1.0; + + mVoxelResultsOut.bvox[0] = iz; + mVoxelResultsOut.bvox[1] = iy; + mVoxelResultsOut.bvox[2] = ix; + mVoxelResultsOut.bsec = isec; + mVoxelResultsOut.flags = 7; + + mTreeOut->Fill(); + } + } + } + } + + // write orbit and lumi info + outputfile->cd(); + TTree tOrbitLumi("OrbitLumiInfo", "Orbit and Lumi Info"); + int64_t orbitResetMS = orbitResetTimeMS; + tOrbitLumi.Branch("orbitResetTimeMS", &orbitResetMS); + tOrbitLumi.Branch("timeMSsel", &timeMSsel); + tOrbitLumi.Branch("orbitsSel", &orbitsSel); + tOrbitLumi.Branch("idcScalerASel", &idcScalerASel); + tOrbitLumi.Branch("idcScalerCSel", &idcScalerCSel); + tOrbitLumi.Branch("ctpLumiSel", &ctpLumiSel); + tOrbitLumi.Branch("meanIDC", &meanIDC); + tOrbitLumi.Branch("meanCTP", &meanCTP); + tOrbitLumi.Branch("meanTimeMS", &meanTimeMS); + tOrbitLumi.Branch("medianIDC", &medianIDC); + tOrbitLumi.Branch("medianCTP", &medianCTP); + tOrbitLumi.Branch("medianTimeMS", &medianTimeMS); + tOrbitLumi.Fill(); + tOrbitLumi.Write(); + + { + TTree tMetaData("MetaData", "Meta data information"); + + int nSlicesP = nSlices; + int maxTracksP = maxTracks; + int minTracksP = minTracks; + + tMetaData.Branch("runNumber", &runNumber); + tMetaData.Branch("nSlicesPhiZ", &nSlicesP); + tMetaData.Branch("maxTracks", &maxTracksP); + tMetaData.Branch("minTracks", &minTracksP); + tMetaData.Branch("nTracksProcessed", &nTracksProcessed_final); + tMetaData.Branch("fileOutput", &fileOutput); + tMetaData.Branch("voxMapInput", &voxMapInput); + tMetaData.Branch("trackSources", &trackSources); + tMetaData.Branch("z2xBinning", &z2xBinning); + tMetaData.Branch("y2xBinning", &y2xBinning); + tMetaData.Branch("useSmoothed", &useSmoothed); + tMetaData.Branch("createSpline", &createSpline); + tMetaData.Branch("maxTracksPerSlice", &maxTracksPerSlice); + tMetaData.Branch("minTracksPerSlice", &minTracksPerSlice); + tMetaData.Branch("badRangeList", &badRangeList); + tMetaData.Branch("firstTFTime", &firstTFTime); + tMetaData.Branch("lastTFTime", &lastTFTime); + tMetaData.Fill(); + tMetaData.Write(); + } + + outputfile->Write(); + //---------------------------------------------------------------- + + const std::string fileOutputInfo = fmt::format("{}/{}.txt", outPath, pFileOutput.stem().c_str()); + std::ofstream fInfo(fileOutputInfo); + fInfo << "meanIDC: " << meanIDC << "\n"; + fInfo << "meanCTP: " << meanCTP << "\n"; + fInfo << "meanTimeMS: " << meanTimeMS << "\n"; + fInfo << "medianIDC: " << medianIDC << "\n"; + fInfo << "medianCTP: " << medianCTP << "\n"; + fInfo << "medianTimeMS: " << medianTimeMS << "\n"; + fInfo.close(); + LOGP(info, "Found meanIDC: {}", meanIDC); + LOGP(info, "Found meanCTP: {}", meanCTP); + LOGP(info, "Found meanTimeMS: {}", meanTimeMS); + LOGP(info, "Found medianIDC: {}", medianIDC); + LOGP(info, "Found medianCTP: {}", medianCTP); + LOGP(info, "Found medianTimeMS: {}", medianTimeMS); + + // This macro only produces the raw voxel-residual map. Turning that into a TPCFastTransform is a + // separate step, done afterwards by TPCFastTransformInitCPM.C on this output, so that the two can be + // rerun independently. The createSpline argument is kept because it is recorded in the output + // metadata and consumed by that later step. + + // The per-thread input handles are locals inside doFileProcessing and are already released, in + // dependency order, when each thread returns -- nothing to tear down here. + mTreeOut.reset(); + + const auto t_end_1 = std::chrono::high_resolution_clock::now(); + LOGP(info, "Wall-clock time for the whole application: {:.1f} s", + std::chrono::duration(t_end_1 - t_start).count()); + + LOGP(info, "Done processing"); + + printMemoryUsage("Memory usage at end"); +} + +std::vector loadRunTimeSpans(const std::string& flname, int onlyRun, const std::string& selection = "ALL") +{ + std::ifstream inputFile(flname); + if (!inputFile) { + LOGP(fatal, "Failed to open selected run/timespans file {}", flname); + } + LOGP(info, "Reading bad ranges from file {}, for run {}", flname, onlyRun); + auto& ccdbmgr = o2::ccdb::BasicCCDBManager::instance(); + ccdbmgr.setURL("http://alice-ccdb.cern.ch"); + + ccdbmgr.setCaching(true); + ccdbmgr.setFatalWhenNull(false); + + std::vector badRanges; + + std::string line; + size_t cntl = 0, cntr = 0; + int64_t orbitResetTimeMS = 0; + int lastRunOrbitReset = -1; + while (std::getline(inputFile, line)) { + cntl++; + for (char& ch : line) { // Replace semicolons and tabs with spaces for uniform processing + if (ch == ';' || ch == '\t' || ch == ',') { + ch = ' '; + } + } + o2::utils::Str::trim(line); + if (line.size() < 1 || line[0] == '#') { + continue; + } + auto tokens = o2::utils::Str::tokenize(line, ' '); + auto logError = [&cntl, &line]() { LOGP(error, "Expected format for selection is tripplet , failed on line#{}: {}", cntl, line); }; + if (tokens.size() >= 3) { + int run = 0; + long rmin, rmax; + try { + run = std::stoi(tokens[0]); + rmin = std::stol(tokens[1]); + rmax = std::stol(tokens[2]); + } catch (...) { + logError(); + continue; + } + + if (onlyRun != run) { + continue; + } + + if (selection != "ALL") { + bool isSelection = false; + for (int iToken = 3; iToken < int(tokens.size()); ++iToken) { + if (tokens[iToken] == selection) { + isSelection = true; + } + } + if (isSelection == false) { + continue; + } + } + + constexpr long ISTimeStamp = 1514761200000L; + int isTimeStampMin = rmin > ISTimeStamp ? 1 : 0, isTimeStampMax = rmax > ISTimeStamp ? 1 : 0; // values above ISTimeStamp are timestamps (need to be converted to orbits) + if (rmin > rmax) { + LOGP(fatal, "Provided range limits are not in increasing order, entry is {}", line); + } + if (isTimeStampMin != isTimeStampMax) { + LOGP(fatal, "Provided range limits should be both consistent either with orbit number or with unix timestamp in ms, entry is {}", line); + } + if (isTimeStampMin) { + if (lastRunOrbitReset != run) { + LOGP(info, "Input needs conversion from time stamps to orbit"); + const auto [sor, eor] = ccdbmgr.getRunDuration(run); + const long timeMeanRun = (sor + eor) / 2.; + const double lengthRun = (eor - sor); + const auto orbitResetTimeNS = ccdbmgr.getSpecific>("CTP/Calib/OrbitReset", timeMeanRun); + orbitResetTimeMS = (*orbitResetTimeNS)[0] * 1e-3; + + LOGP(info, "Run {}, sor {}, eor {}, duration {} (min)", run, sor, eor, lengthRun / 1000. / 60.); + LOGP(info, "Orbit reset time in MS is {}", orbitResetTimeMS); + lastRunOrbitReset = run; + } + const auto orbitToMS = o2::constants::lhc::LHCOrbitMUS * 1e-3; + const auto rMinIn = rmin, rMaxIn = rmax; + rmin = long((rmin - orbitResetTimeMS) / orbitToMS); + rmax = long(std::ceil((rmax - orbitResetTimeMS) / orbitToMS)); + LOGP(info, "Run {} input range [{} - {}] ms -> [{} - {}] orbits", run, rMinIn, rMaxIn, rmin, rmax); + } + + badRanges.emplace_back(rmin, rmax); + cntr++; + } else { + logError(); + } + } + return badRanges; +} diff --git a/Detectors/TPC/calibration/SpacePoints/macro/voxResQA.C b/Detectors/TPC/calibration/SpacePoints/macro/voxResQA.C new file mode 100644 index 0000000000000..dbc766b2571d4 --- /dev/null +++ b/Detectors/TPC/calibration/SpacePoints/macro/voxResQA.C @@ -0,0 +1,445 @@ +// 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 "TChain.h" +#include "TMath.h" +#include "TCanvas.h" +#include "THStack.h" +#include "TLegend.h" +#include "TProfile.h" +#include "TProfile2D.h" +#include "TObjArray.h" +#include "TStyle.h" +#include "TSystem.h" +#include "TLatex.h" +#include "TPaletteAxis.h" +#include "TPCBase/Utils.h" +#include "TPCBaseRecSim/Painter.h" +#include "CommonUtils/StringUtils.h" +#include "Framework/Logger.h" +#include "SpacePoints/TrackResiduals.h" + +std::string getString(TList* l, const std::string& name, std::string defaultVal = "0"); +int getNbins(const std::string& name, int defaultVal = 20, const char delim = ','); +void setStyle(TH1* hist, int idir); +void setSizes(TAxis* axis, float titleSize, float titleOffset, float labelSize); + +using namespace o2::tpc; +namespace fs = std::filesystem; + +void voxResQA(std::string inFilesCmd, std::string outFileNameAdd = "", bool drawErrors = false, bool useSmoothed = true, int z2xBinSel = -1) +{ + gStyle->SetOptStat(0); + + TObjArray arrCanvases1D; + arrCanvases1D.SetName("voxResQA"); + TObjArray arrCanvases2D; + arrCanvases2D.SetName("voxResQA"); + + auto hsdXvsRowA = new THStack; + auto hsdXvsRowC = new THStack; + auto hsdZvsRowA = new THStack; + auto hsdZvsRowC = new THStack; + + const std::string dXtitle = useSmoothed ? "dXS" : "dX"; + const std::string dYtitle = useSmoothed ? "dYS" : "dY"; + const std::string dZtitle = useSmoothed ? "dZS" : "dZ"; + + const std::string sFiles(gSystem->GetFromPipe(inFilesCmd.data())); + const auto arrFiles = o2::utils::Str::tokenize(sFiles, '\n'); + if (arrFiles.size() == 0) { + return; + } + int idir = 0; + for (const auto& inFile : arrFiles) { + const fs::path inPath(inFile); + const std::string fileTitle(std::string(inPath.stem()).substr(17)); + auto fileIDName = fileTitle; + std::replace(fileIDName.begin(), fileIDName.end(), '.', '_'); + std::replace(fileIDName.begin(), fileIDName.end(), '-', '_'); + TChain cVoxRes("voxResTree"); + cVoxRes.AddFile(inFile.data()); + + if (!cVoxRes.GetBranch("voxRes")) { + LOGP(error, "Could not find branch voxRes in file '{}'", inFile); + continue; + } + + o2::tpc::TrackResiduals::VoxRes* vox{nullptr}; + cVoxRes.SetBranchAddress("voxRes", &vox); + + // retrieve binning + // OBJ: TNamed meanIDC 0.295421 + // OBJ: TNamed meanCTP 118553.744497 + // OBJ: TNamed y2xBinning 20 + // OBJ: TNamed z2xBinning 20 + + cVoxRes.GetEntry(0); + const auto userInfo = cVoxRes.GetTree()->GetUserInfo(); + userInfo->Print(); + const auto y2xBinning = getString(userInfo, "y2xBinning"); + const auto z2xBinning = getString(userInfo, "z2xBinning"); + const auto meanIDC = std::stof(getString(userInfo, "meanIDC")); + const auto meanCTP = std::stof(getString(userInfo, "meanCTP")); + + const int nbinsY2X = getNbins(y2xBinning); + const int nBinsSector = nbinsY2X * 36; + const int nbinsZ2X = getNbins(z2xBinning); + + // ===| 1D histograms |===================================================== + TObjArray arrHists1D; + auto hEntriesDist = new TH1F(("hEntriesDist" + fileIDName).data(), (fileTitle + ";#entries").data(), 500, 0, 5000); + arrHists1D.Add(hEntriesDist); + auto hdXDist = new TH1F(("hdXDist" + fileIDName).data(), (fileTitle + ";" + dXtitle + " (cm)").data(), 100, -10, 20); + arrHists1D.Add(hdXDist); + auto hdYDist = new TH1F(("hdYDist" + fileIDName).data(), (fileTitle + ";" + dYtitle + " (cm)").data(), 100, -10, 10); + arrHists1D.Add(hdYDist); + auto hdZDist = new TH1F(("hdZDist" + fileIDName).data(), (fileTitle + ";" + dZtitle + " (cm)").data(), 100, -10, 10); + arrHists1D.Add(hdZDist); + auto hdYSigmaMAD = new TProfile(("hdYSigmaMAD" + fileIDName).data(), (fileTitle + ";sector;<#sigma_{MAD}(dY)> (cm)").data(), nBinsSector, 0, nBinsSector); + arrHists1D.Add(hdYSigmaMAD); + + auto hdXvsRowA = new TProfile(("hdXvsRowA" + fileIDName).data(), (fileTitle + " (A-Side);row;<" + dXtitle + "> (cm) z2xBin 0").data(), 152, 0, 152); + hsdXvsRowA->Add(hdXvsRowA); + hsdXvsRowA->SetTitle("(A-Side);row; (cm) z2xBin 0"); + setStyle(hdXvsRowA, idir); + + auto hdXvsRowC = new TProfile(("hdXvsRowC" + fileIDName).data(), (fileTitle + " (C-Side);row;<" + dXtitle + "> (cm) z2xBin 0").data(), 152, 0, 152); + hsdXvsRowC->Add(hdXvsRowC); + hsdXvsRowC->SetTitle("(C-Side);row; (cm) z2xBin 0"); + setStyle(hdXvsRowC, idir); + + auto hdZvsRowA = new TProfile(("hdZvsRowA" + fileIDName).data(), (fileTitle + " (A-Side);row;<" + dZtitle + "> (cm) z2xBin 0").data(), 152, 0, 152); + hsdZvsRowA->Add(hdZvsRowA); + hsdZvsRowA->SetTitle("(A-Side);row; (cm) z2xBin 0"); + setStyle(hdZvsRowA, idir); + + auto hdZvsRowC = new TProfile(("hdZvsRowC" + fileIDName).data(), (fileTitle + " (C-Side);row;<" + dZtitle + "> (cm) z2xBin 0").data(), 152, 0, 152); + hsdZvsRowC->Add(hdZvsRowC); + hsdZvsRowC->SetTitle(("(C-Side);row;" + dZtitle + " (cm) z2xBin 0").data()); + setStyle(hdZvsRowC, idir); + + auto hdZvsZ2X = new TProfile(("hdZvsZ2X" + fileIDName).data(), (fileTitle + ";z2x-bin;<" + dZtitle + "> (cm) row 80").data(), 2 * nbinsZ2X, -nbinsZ2X, nbinsZ2X); + arrHists1D.Add(hdZvsZ2X); + + // ===| 2D histograms |===================================================== + TObjArray arrHists2D; + auto hMeanEntries = new TProfile2D(("hMeanEntries_" + fileIDName).data(), (fileTitle + ";sector;row;").data(), nBinsSector, 0, nBinsSector, 152, 0, 152); + arrHists2D.Add(hMeanEntries); + auto hMeanEntrieszRow = new TProfile2D(("hMeanEntrieszRow" + fileIDName).data(), (fileTitle + ";z-bin;row; (cm)").data(), 2 * nbinsZ2X, -nbinsZ2X, nbinsZ2X, 152, 0, 152); + arrHists2D.Add(hMeanEntrieszRow); + auto hSigmaMADSecRow = new TProfile2D(("hSigmaMADSecRow" + fileIDName).data(), (fileTitle + ";sector;row;<#sigma_{MAD}(dY)> (cm)").data(), nBinsSector, 0, nBinsSector, 152, 0, 152); + arrHists2D.Add(hSigmaMADSecRow); + auto hdXSecRow = new TProfile2D(("hdXSecRow" + fileIDName).data(), (fileTitle + ";sector;row;<" + dXtitle + "> (cm)").data(), nBinsSector, 0, nBinsSector, 152, 0, 152); + arrHists2D.Add(hdXSecRow); + auto hdYSecRow = new TProfile2D(("hdYSecRow" + fileIDName).data(), (fileTitle + ";sector;row;<" + dYtitle + "> (cm)").data(), nBinsSector, 0, nBinsSector, 152, 0, 152); + arrHists2D.Add(hdYSecRow); + auto hdZzRow = new TProfile2D(("hdZzRow" + fileIDName).data(), (fileTitle + ";z2x-bin;row;<" + dZtitle + "> (cm)").data(), 2 * nbinsZ2X, -nbinsZ2X, nbinsZ2X, 152, 0, 152); + arrHists2D.Add(hdZzRow); + + TProfile2D* hdXESecRow = nullptr; + TProfile2D* hdYESecRow = nullptr; + TProfile2D* hdZEzRow = nullptr; + if (drawErrors) { + hdXESecRow = new TProfile2D(("hdXESecRow" + fileIDName).data(), (fileTitle + ";sector;row; (cm)").data(), nBinsSector, 0, nBinsSector, 152, 0, 152); + arrHists2D.Add(hdXESecRow); + hdYESecRow = new TProfile2D(("hdYESecRow" + fileIDName).data(), (fileTitle + ";sector;row; (cm)").data(), nBinsSector, 0, nBinsSector, 152, 0, 152); + arrHists2D.Add(hdYESecRow); + hdZEzRow = new TProfile2D(("hdZEzRow" + fileIDName).data(), (fileTitle + ";z-bin;row; (cm)").data(), 2 * nbinsZ2X, -nbinsZ2X, nbinsZ2X, 152, 0, 152); + arrHists2D.Add(hdZEzRow); + } + + /* + OBJ: TNamed z2xBin bvox[0] + OBJ: TNamed y2xBin bvox[1] + OBJ: TNamed xBin bvox[2] + OBJ: TNamed z2xAV stat[0] + OBJ: TNamed y2xAV stat[1] + OBJ: TNamed xAV stat[2] + OBJ: TNamed fsector bsec+0.5+9.*(y2xAV)/pi + OBJ: TNamed phi (bsec%18+0.5+9.*(stat[1])/pi)/9*pi + OBJ: TNamed r stat[2] + OBJ: TNamed z z2xAV*xAV + OBJ: TNamed dX D[0] + OBJ: TNamed dY D[1] + OBJ: TNamed dZ D[2] + OBJ: TNamed dXS DS[0] + OBJ: TNamed dYS DS[1] + OBJ: TNamed dZS DS[2] + OBJ: TNamed dXE E[0] + OBJ: TNamed dYE E[1] + OBJ: TNamed dZE E[2] + OBJ: TNamed voxelIndex xBin + 152 * (y2xBin + 20 * z2xBin) + 60800 * bsec + OBJ: TNamed entries stat[3] + OBJ: TNamed fitOK (flags & 1) == 1 + OBJ: TNamed dispOK (flags & 2) == 2 + OBJ: TNamed smtOK (flags & 4) == 4 + OBJ: TNamed masked (flags & 128) == 128 + */ + + float z2xAV = 0; + + int isNaNdX = 0; + int isNaNdY = 0; + int isNaNdZ = 0; + int isNaNdXS = 0; + int isNaNdYS = 0; + int isNaNdZS = 0; + + for (Long64_t iEntry = 0; iEntry < cVoxRes.GetEntries(); ++iEntry) { + cVoxRes.GetEntry(iEntry); + const auto y2xBin = vox->bvox[1]; + const auto z2xBin = vox->bvox[0]; + const auto xBin = vox->bvox[2]; + const auto bsec = vox->bsec; + const auto entries = vox->stat[3]; + const auto dX = vox->D[0]; + const auto dY = vox->D[1]; + const auto dZ = vox->D[2]; + const auto dXS = vox->DS[0]; + const auto dYS = vox->DS[1]; + const auto dZS = vox->DS[2]; + const auto dXE = vox->E[0]; + const auto dYE = vox->E[1]; + const auto dZE = vox->E[2]; + const auto sectorFine = y2xBin + bsec * nbinsY2X; + const auto z2xBinSides = (z2xBin + 0.5) * (1 - 2 * (bsec > 17)); + const auto z = vox->stat[0] * vox->stat[2]; + + const auto dXdraw = useSmoothed ? dXS : dX; + const auto dYdraw = useSmoothed ? dYS : dY; + const auto dZdraw = useSmoothed ? dZS : dZ; + + isNaNdX += TMath::IsNaN(dX); + isNaNdY += TMath::IsNaN(dY); + isNaNdZ += TMath::IsNaN(dZ); + isNaNdXS += TMath::IsNaN(dXS); + isNaNdYS += TMath::IsNaN(dYS); + isNaNdZS += TMath::IsNaN(dZS); + + // only fill values in the acceptance + if (std::abs(z) < 248) { + if (z2xBinSel < 0 || z2xBinSel == z2xBin) { + if (z2xAV == 0) { + z2xAV = vox->stat[0]; + } + hMeanEntries->Fill(sectorFine, xBin, entries); + hdXSecRow->Fill(sectorFine, xBin, dXdraw); + hdYSecRow->Fill(sectorFine, xBin, dYdraw); + hSigmaMADSecRow->Fill(sectorFine, xBin, vox->dYSigMAD); + hdYSigmaMAD->Fill(sectorFine, vox->dYSigMAD); + } + } + + hMeanEntrieszRow->Fill(z2xBinSides, xBin, entries); + hdZzRow->Fill(z2xBinSides, xBin, dZdraw); + if (drawErrors && (z2xBinSel < 0 || z2xBinSel == z2xBin)) { + hdXESecRow->Fill(sectorFine, xBin, dXE); + hdYESecRow->Fill(sectorFine, xBin, dYE); + hdZEzRow->Fill(z2xBinSides, xBin, dZE); + } + + hEntriesDist->Fill(entries); + hdXDist->Fill(dXdraw); + hdYDist->Fill(dYdraw); + hdZDist->Fill(dZdraw); + if (xBin >= 79 && xBin <= 81) { + hdZvsZ2X->Fill(z2xBinSides, dZdraw); + } + + if (z2xBin == 0) { + if (bsec < 18) { + hdXvsRowA->Fill(xBin, dXdraw); + hdZvsRowA->Fill(xBin, dZdraw); + } else { + hdXvsRowC->Fill(xBin, dXdraw); + hdZvsRowC->Fill(xBin, dZdraw); + } + } + } + + std::vector hXadjust{hMeanEntries, hdXSecRow, hdXESecRow, hdYSecRow, hdYESecRow}; + for (auto h : hXadjust) { + if (!h) { + continue; + } + h->GetXaxis()->SetLimits(0, 36); + // 36 exact divisions means a label for every single sector -- unreadable at this pad size, they + // overlap into an illegible block. 12 (label every 3 sectors) still shows the A-/C-side structure + // without the collision. + h->GetXaxis()->SetNdivisions(12, false); + } + + // hEntriesDist's fixed 0-5000 range is mostly empty for real data (per-voxel entry counts rarely get + // anywhere near 5000) -- zoom to just past the last populated bin instead of showing mostly blank + // axis. + if (hEntriesDist->GetEntries() > 0) { + const int lastBin = hEntriesDist->FindLastBinAbove(0); + if (lastBin > 0) { + hEntriesDist->GetXaxis()->SetRangeUser(0, hEntriesDist->GetXaxis()->GetBinUpEdge(lastBin) * 1.1); + } + } + + // ===| output canvases |=================================================== + // + auto c1D = new TCanvas(("c1D_" + fileIDName + outFileNameAdd).data(), fileTitle.data(), 1500, 900); + arrCanvases1D.Add(c1D); + + int ipad = 1; + c1D->DivideSquare(arrHists1D.GetEntries()); + + for (auto o : arrHists1D) { + c1D->cd(ipad++); + o->Draw(); + const std::string name(o->GetName()); + if (o->IsA() != TProfile::Class()) { + gPad->SetLogy(); + } + // Without this, saveCanvas()'s plain c.SaveAs() (Detectors/TPC/base/src/Utils.cxx) can export a + // pad before its log-scale range is actually recomputed in batch mode -- confirmed real: the + // stored TCanvas in voxResQA*_1D.root has all real data (reopening and redrawing it interactively + // shows every panel fine), but the direct PNG export came out blank. The c2D loop below already + // does this after its own Draw() calls; c1D's was missing it. + gPad->Modified(); + gPad->Update(); + } + + auto c2D = new TCanvas(("c2D_" + fileIDName + outFileNameAdd).data(), fileTitle.data(), 1500, 900); + arrCanvases2D.Add(c2D); + + ipad = 1; + c2D->DivideSquare(arrHists2D.GetEntries()); + + TLatex l; + l.SetTextFont(42); + // Without this, DrawLatex's (x,y) below are interpreted in the pad's USER (data-axis) coordinates, + // not normalized pad fractions -- 0.75/0.85 then lands almost at the frame's bottom-left corner + // (sector~0.75 of 36, row~0.85 of 152), overlapping the plotted content instead of sitting clear of + // it. SetNDC() makes the coordinates pad-fraction-based, and the y moved down (below the frame, + // rather than "0.85" which was never actually near the top). + l.SetNDC(); + + for (auto o : arrHists2D) { + auto h = static_cast(o); + c2D->cd(ipad++); + h->Draw("colz"); + const std::string name(h->GetName()); + if (name.find("hdZzRow") == 0) { + l.DrawLatex(0.4, 0.02, "y2xBin averaged"); + } else { + if (z2xBinSel >= 0) { + l.DrawLatex(0.4, 0.02, fmt::format("z2xBin = {} ({})", z2xBinSel, z2xAV).data()); + } else { + l.DrawLatex(0.4, 0.02, "z2xBin averaged"); + } + } + gPad->Modified(); + gPad->Update(); + + auto palette = (TPaletteAxis*)h->GetListOfFunctions()->FindObject("palette"); + if (palette) { + painter::adjustPalette(h, 0.92); + } + } + + ++idir; + int nNaN = isNaNdX + isNaNdY + isNaNdZ; + int nNaNS = isNaNdXS + isNaNdYS + isNaNdZS; + const auto sNaN = fmt::format("NaN: {} - {} {} {} {}", inFile, nNaN, isNaNdX, isNaNdY, isNaNdZ); + const auto sNaNS = fmt::format("NaNS: {} - {} {} {} {}", inFile, nNaNS, isNaNdXS, isNaNdYS, isNaNdZS); + if (nNaN > 0) { + LOGP(error, "{}", sNaN); + } else { + LOGP(info, "{}", sNaN); + } + if (nNaNS > 0) { + LOGP(error, "{}", sNaNS); + } else { + LOGP(info, "{}", sNaNS); + } + } + + // ===| dX/dX vs row |======================================================== + // + auto cdXZvsRow = new TCanvas(fmt::format("cdXZvsRow{}", outFileNameAdd).data(), "dX/dZ vs row", 1500, 900); + cdXZvsRow->SetRightMargin(0.01); + cdXZvsRow->SetBottomMargin(0.15); + cdXZvsRow->Divide(1, 4, -1, -1); + cdXZvsRow->cd(1); + gPad->SetGrid(); + hsdXvsRowA->Draw("nostack"); + setSizes(hsdXvsRowA->GetYaxis(), 0.1, 0.4, 0.08); + setSizes(hsdXvsRowA->GetXaxis(), 0.1, 0.4, 0.08); + auto leg = gPad->BuildLegend(0.5, 0.5, 0.9, 0.9); + leg->SetMargin(0.05); + cdXZvsRow->cd(2); + gPad->SetGrid(); + hsdXvsRowC->Draw("nostack"); + setSizes(hsdXvsRowC->GetYaxis(), 0.1, 0.4, 0.08); + setSizes(hsdXvsRowC->GetXaxis(), 0.1, 0.4, 0.08); + cdXZvsRow->cd(3); + gPad->SetGrid(); + hsdZvsRowA->Draw("nostack"); + setSizes(hsdZvsRowA->GetYaxis(), 0.1, 0.4, 0.08); + setSizes(hsdZvsRowA->GetXaxis(), 0.1, 0.4, 0.08); + cdXZvsRow->cd(4); + gPad->SetGrid(); + hsdZvsRowC->Draw("nostack"); + setSizes(hsdZvsRowC->GetYaxis(), 0.1, 0.4, 0.08); + setSizes(hsdZvsRowC->GetXaxis(), 0.1, 0.4, 0.08); + + arrCanvases1D.Add(cdXZvsRow); + + // ===| save canvases |======================================================= + // + o2::tpc::utils::saveCanvases(arrCanvases1D, "./", "png,png", fmt::format("voxResQA{}_1D.root", outFileNameAdd.data())); + o2::tpc::utils::saveCanvases(arrCanvases2D, "./", "png,png", fmt::format("voxResQA{}_2D.root", outFileNameAdd.data())); +} + +std::string getString(TList* l, const std::string& name, std::string defaultVal) +{ + if (!l || !l->FindObject(name.data())) { + return defaultVal; + } + return l->FindObject(name.data())->GetTitle(); +} + +int getNbins(const std::string& name, int defaultVal, const char delim) +{ + if (name.find(delim) != name.npos) { + return std::count(name.begin(), name.end(), delim) + 1; + } + return std::stoi(name); +} + +void setStyle(TH1* hist, int idir) +{ + const std::vector colors = {kRed + 2, kOrange + 1, kGreen + 2, kAzure + 10, kBlue + 2, kMagenta + 1}; + const std::vector markers{20, 24, 21, 25, 47, 46, 34, 28}; + const std::vector styles{kSolid, kDashed, kDotted, kDashDotted}; + + hist->SetMarkerSize(1); + hist->SetMarkerColor(colors[idir % colors.size()]); + hist->SetLineColor(colors[idir % colors.size()]); + hist->SetMarkerStyle(markers[idir % markers.size()]); + hist->SetLineStyle(styles[(idir / colors.size()) % styles.size()]); +} + +void setSizes(TAxis* axis, float titleSize, float titleOffset, float labelSize) +{ + axis->SetTitleSize(titleSize); + axis->SetTitleOffset(titleOffset); + axis->SetLabelSize(labelSize); +} diff --git a/GPU/TPCFastTransformation/CMakeLists.txt b/GPU/TPCFastTransformation/CMakeLists.txt index c4fb7c04796f2..d48f7660c45b9 100644 --- a/GPU/TPCFastTransformation/CMakeLists.txt +++ b/GPU/TPCFastTransformation/CMakeLists.txt @@ -88,6 +88,22 @@ if(${ALIGPU_BUILD_TYPE} STREQUAL "O2") ${CMAKE_BINARY_DIR}/stage/include LABELS gpu tpc COMPILE_ONLY) endforeach() + + o2_add_test_root_macro(macro/TPCFastTransformInitCPM.C + PUBLIC_LINK_LIBRARIES O2::TPCFastTransformation + O2::DataFormatsTPC + O2::TPCSimulation + O2::TPCReconstruction + O2::TPCCalibration + O2::SpacePoints + O2::CommonUtils + O2::MathUtils + O2::Algorithm + O2::Framework + PUBLIC_INCLUDE_DIRECTORIES + ${CMAKE_BINARY_DIR}/stage/include + LABELS gpu tpc COMPILE_ONLY) + foreach(m IrregularSpline1DTest.C IrregularSpline2D3DCalibratorTest.C @@ -107,6 +123,9 @@ if(${ALIGPU_BUILD_TYPE} STREQUAL "O2") install(FILES macro/TPCFastTransformInit.C DESTINATION share/macro/) + + install(FILES macro/TPCFastTransformInitCPM.C + DESTINATION share/macro/) endif() if(ALIGPU_BUILD_TYPE STREQUAL "Standalone") diff --git a/GPU/TPCFastTransformation/macro/TPCFastTransformInitCPM.C b/GPU/TPCFastTransformation/macro/TPCFastTransformInitCPM.C new file mode 100644 index 0000000000000..a83efef0cdd75 --- /dev/null +++ b/GPU/TPCFastTransformation/macro/TPCFastTransformInitCPM.C @@ -0,0 +1,740 @@ +// Copyright 2019-2023 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file TPCFastTransformInitCPM.C +/// \brief A macro for generating TPC fast transformation +/// out of set of space charge correction voxels +/// +/// \author Sergey Gorbunov +/// + +/// how to run the macro: +/// +/// root -l TPCFastTransformInitCPM.C'("debugVoxRes.root")' +/// + +#if !defined(__CLING__) || defined(__ROOTCLING__) + +#include +#include +#include +#include +#include +#include "TFile.h" +#include "TSystem.h" +#include "TTree.h" +#include "TNtuple.h" +#include "TMath.h" +#include "Riostream.h" + +#include "CommonUtils/TreeStreamRedirector.h" +#include "MathUtils/fit.h" +#include "Algorithm/RangeTokenizer.h" +#include "Framework/Logger.h" +#include "GPU/TPCFastTransform.h" +#include "SpacePoints/TrackResiduals.h" +#include "TPCReconstruction/TPCFastTransformHelperO2.h" +#include "TPCCalibration/TPCFastSpaceChargeCorrectionHelper.h" + +#endif + +struct sMean { + float mean; + float rms; + float meanErr; + float rmsErr; + + ClassDef(sMean, 1); +}; + +#pragma link C++ class sMean + ; +#pragma link C++ class std::vector < sMean> + ; + +using namespace o2::tpc; +using namespace o2::gpu; +namespace mu = o2::math_utils; + +void createFastTransform(std::string outFileName, TTree* voxResTree, o2::tpc::TrackResiduals& trackResiduals, bool useSmoothed, bool invertSigns, float meanIDC = 0.f, float meanCTP = 0.f, int debug = 0, int useCTPLumi = 0); + +void TPCFastTransformInitCPM(const char* fileName = "debugVoxRes.root", + const char* outFileName = "TPCFastTransform_VoxRes.root", + bool useSmoothed = false, + int useCTPLumi = 0, + bool invertSigns = true, + float meanIDC = 0.f, + float meanCTP = 0.f, + int debug = 0, + int nThreads = 8) +{ + + // Initialise TPCFastTransform object from "voxRes" tree of + // o2::tpc::TrackResiduals::VoxRes track residual voxels + // + + /* + To visiualise the results: + + root -l transformDebug.root + corr->Draw("cx:y:z","iRoc==0&&iRow==10","") + grid->Draw("cx:y:z","iRoc==0&&iRow==10","same") + vox->Draw("vx:y:z","iRoc==0&&iRow==10","same") + */ + + o2::tpc::TPCFastSpaceChargeCorrectionHelper::instance()->setNthreads(nThreads); + + if (gSystem->AccessPathName(fileName)) { + LOGP(error, "input file {} does not exist!", fileName); + return; + } + + auto file = std::unique_ptr(TFile::Open(fileName, "READ")); + if (!file || !file->IsOpen()) { + LOGP(error, "could not open input file {}", fileName); + return; + } + + TTree* voxResTree = nullptr; + file->cd(); + gDirectory->GetObject("voxResTree", voxResTree); + if (!voxResTree) { + LOGP(error, "tree voxResTree does not exist!"); + return; + } + auto userInfo = voxResTree->GetUserInfo(); + if (!userInfo->FindObject("y2xBinning") || !userInfo->FindObject("z2xBinning")) { + LOGP(error, "'y2xBinning' or 'z2xBinning' not found in UserInfo, but required to get the correct binning"); + return; + } + + // Obtain configuration + const SpacePointsCalibConfParam& params = SpacePointsCalibConfParam::Instance(); + if (!std::filesystem::exists("scdconfig.ini")) { + LOGP(warning, "Did not find configuration file. Using default parameters and storing them in scdconfig.ini"); + params.writeINI("scdconfig.ini", "scdcalib"); // to write default parameters to a file + } else { + params.updateFromFile("scdconfig.ini"); + } + // TrackResiduals::setZ2XBinning() (called below) reads scdcalib.maxZ2X directly and uses it to scale + // the physical z/x bin boundaries -- baked into what each z2x voxel index in the input tree actually + // means, not a cosmetic knob. There is no scdconfig.ini on the GRID, so without this the code default + // (1.0) would silently apply instead of whatever stage 1 actually used (production 1.4), misaligning + // this macro's re-derived binning against the tree's real geometry. Must happen BEFORE setZ2XBinning() + // below. See staticMapCreatorCPM.C's UserInfo::Add("maxZ2X", ...) for where this comes from -- + // SmoothingExtrapolate.C clones the whole input UserInfo onto its own output, so it survives into this + // macro's input unchanged. + if (auto* maxZ2XObj = userInfo->FindObject("maxZ2X")) { + const std::string maxZ2XStr = maxZ2XObj->GetTitle(); + o2::conf::ConfigurableParam::setValue("scdcalib.maxZ2X", maxZ2XStr); + LOGP(info, "Set scdcalib.maxZ2X = {} from input UserInfo (matches stage 1)", maxZ2XStr); + } else { + LOGP(warning, + "'maxZ2X' not found in input UserInfo (older input file?) -- using scdcalib.maxZ2X = {} " + "(scdconfig.ini/code default), which may NOT match the value stage 1 actually used to " + "build this tree's z2x binning!", + params.maxZ2X); + } + + LOGP(info, "----- Dumping configuration values START -----"); + params.printKeyValues(); + LOGP(info, "----- Dumping configuration values END -----"); + + // required for the binning that was used + o2::tpc::TrackResiduals trackResiduals; + auto y2xBins = o2::RangeTokenizer::tokenize(userInfo->FindObject("y2xBinning")->GetTitle()); + const std::string z2xStr = userInfo->FindObject("z2xBinning")->GetTitle(); + auto z2xBins = o2::RangeTokenizer::tokenize(z2xStr); + LOGP(info, "z2xBins: {}", z2xStr); + trackResiduals.setY2XBinning(y2xBins); + trackResiduals.setZ2XBinning(z2xBins); + trackResiduals.init(); + + auto getFromUserInfo = [userInfo](std::string value, float& valueF) { + if (valueF != 0) { + LOGP(info, "{} set to {} via command line, not reading it from userInfo", value, valueF); + return; + } + + if (!userInfo || !userInfo->FindObject(value.data())) { + LOGP(error, "Could not find value for {} in userInfo", value); + valueF = 0.f; + return; + } + valueF = std::atof(userInfo->FindObject(value.data())->GetTitle()); + LOGP(info, "Found {} = {} in userInfo", value, valueF); + }; + + getFromUserInfo("meanIDC", meanIDC); + getFromUserInfo("meanCTP", meanCTP); + + if ((useCTPLumi != 2 && meanIDC == 0) || meanCTP == 0) { + LOGP(fatal, "meanCTP ({}) or meanIDC ({}) not set!", meanCTP, meanIDC); + } + if (useCTPLumi == 2) { + LOGP(warning, "Explicitly disabled IDCs!"); + } + + createFastTransform(outFileName, voxResTree, trackResiduals, useSmoothed, invertSigns, meanIDC, meanCTP, debug, useCTPLumi); +} + +void createFastTransform(std::string outFileName, TTree* voxResTree, o2::tpc::TrackResiduals& trackResiduals, bool useSmoothed, bool invertSigns, float meanIDC, float meanCTP, int debug, int useCTPLumi) +{ + LOGP(info, "create fast transformation ... "); + std::regex reg(".*FT_voxRes\\.residuals\\.([0-9]{6})_([0-9]{13})_([0-9]{13})_([0-9]+)_([0-9]+).*\\.root"); + std::smatch base_match; + int run = -1; + long validFrom = -1; + long validUntil = -1; + int firstTF = -1; + int lastTF = -1; + if (std::regex_match(outFileName, base_match, reg)) { + run = std::stoi(base_match[1].str()); + validFrom = std::stol(base_match[2].str()); + validUntil = std::stol(base_match[3].str()); + firstTF = std::stol(base_match[4].str()); + lastTF = std::stol(base_match[5].str()); + LOGP(info, "Found run {}, validFrom {}, validUntil {}, firstTF {}, lastTF {}", run, validFrom, validUntil, firstTF, lastTF); + } + + auto* helper = o2::tpc::TPCFastTransformHelperO2::instance(); + + o2::tpc::TPCFastSpaceChargeCorrectionHelper* corrHelper = o2::tpc::TPCFastSpaceChargeCorrectionHelper::instance(); + +#if __has_include("TPCFastTransformPOD.h") + TTree* voxResTreeInverse = nullptr; + o2::gpu::TPCFastSpaceChargeCorrectionMap mapDirect(0, 0), mapInverse(0, 0); + auto corrPtr = corrHelper->createFromTrackResiduals(trackResiduals, voxResTree, voxResTreeInverse, useSmoothed, invertSigns, &mapDirect, &mapInverse); +#else + auto corrPtr = corrHelper->createFromTrackResiduals(trackResiduals, voxResTree, useSmoothed, invertSigns); +#endif + + std::unique_ptr fastTransform(helper->create(0, *corrPtr)); + fastTransform->setLumi(meanCTP); + fastTransform->setIDC(meanIDC); // for SW version with IDC in FastTransfrom + + o2::gpu::TPCFastSpaceChargeCorrection& corr = fastTransform->getCorrection(); + + LOGP(info, "... create fast transformation completed"); + + if (!outFileName.empty()) { + fastTransform->writeToFile(outFileName.data(), "ccdb_object"); + } + + LOGP(info, "verify the results ..."); + + // the difference + + double maxDiff[3] = {0., 0., 0.}; + int maxDiffRoc[3] = {0, 0, 0}; + int maxDiffRow[3] = {0, 0, 0}; + + double sumDiff[3] = {0., 0., 0.}; + long nDiff = 0; + + // a debug file with some NTuples + + TDirectory* currDir = gDirectory; + + const std::filesystem::path pFileOutput(outFileName); + std::string outPath(pFileOutput.parent_path().c_str()); + if (outPath.empty()) { + outPath = "."; + } + const std::string fileOutputDebug = fmt::format("{}/{}.debug.root", outPath, pFileOutput.stem().c_str()); + const std::string fileOutputSummary = fmt::format("{}/{}.summary.root", outPath, pFileOutput.stem().c_str()); + + o2::utils::TreeStreamRedirector summary(fileOutputSummary.data(), "recreate"); + + TFile* debugFile = new TFile(fileOutputDebug.data(), "RECREATE"); + debugFile->cd(); + + // ntuple with the input data: voxel corrections + debugFile->cd(); + TNtuple* debugVox = new TNtuple("vox", "vox", "iRoc:iRow:y2xbin:z2xbin:x:y:z:vx:vy:vz:cx:cy:cz"); + + debugVox->SetMarkerStyle(8); + debugVox->SetMarkerSize(0.8); + debugVox->SetMarkerColor(kBlue); + + currDir->cd(); + + // check the difference in voxels and fill corresp. debug ntuple + + LOGP(info, "verify the results ..."); + + const o2::gpu::TPCFastTransformGeo& geo = helper->getGeometry(); + + o2::tpc::TrackResiduals::VoxRes* v = nullptr; + TBranch* branch = voxResTree->GetBranch("voxRes"); + branch->SetAddress(&v); + branch->SetAutoDelete(kTRUE); + + int nNaNdXV = 0; + int nNaNdYV = 0; + int nNaNdZV = 0; + int nNaNdXC = 0; + int nNaNdYC = 0; + int nNaNdZC = 0; + + int lastSector = -1; + + std::vector statsPerSecMean(36); + std::vector statsPerSecStdDev(36); + std::vector statsPerSecMedian(36); + + std::vector maxDiffPerSec[3]; + + std::vector deviationPerSecLTM95[3]; + std::vector deviationPerSecMedian[3]; + + std::vector entriesStats; + std::vector deviations[3]; + entriesStats.reserve(152 * trackResiduals.getNY2XBins() * trackResiduals.getNZ2XBins()); + for (int i = 0; i < 3; ++i) { + deviations[i].reserve(152 * trackResiduals.getNY2XBins() * trackResiduals.getNY2XBins()); + maxDiffPerSec[i].resize(36); + deviationPerSecLTM95[i].resize(36); + deviationPerSecMedian[i].resize(36); + } + + // retrieve infos from UserInfo + auto getFromUserInfo = [](TList* u, const char* name, int defVal = -1) { + const auto o = u->FindObject(name); + return o ? std::atoi(o->GetTitle()) : defVal; + }; + + auto userInfo = voxResTree->GetUserInfo(); + const int nSlicesPhiZ = getFromUserInfo(userInfo, "nSlicesPhiZ"); + const int maxTracks = getFromUserInfo(userInfo, "maxTracks"); + const int minTracks = getFromUserInfo(userInfo, "minTracks"); + const int nTracksProcessed = getFromUserInfo(userInfo, "nTracksProcessed"); + const bool badCalib = static_cast(getFromUserInfo(userInfo, "badCalib", 0)); + + for (int iVox = 0; iVox < voxResTree->GetEntriesFast(); iVox++) { + + voxResTree->GetEntry(iVox); + + const float voxEntries = v->stat[o2::tpc::TrackResiduals::VoxV]; + const int xBin = v->bvox[o2::tpc::TrackResiduals::VoxX]; // bin number in x (= pad row) + const int y2xBin = v->bvox[o2::tpc::TrackResiduals::VoxF]; // bin number in y/x 0..14 + const int z2xBin = v->bvox[o2::tpc::TrackResiduals::VoxZ]; // bin number in z/x 0..4 + const int iRoc = (int)v->bsec; + const int iRow = (int)xBin; + + const float x = trackResiduals.getX(xBin); // radius of the pad row + const float y2x = trackResiduals.getY2X(xBin, y2xBin); // y/x coordinate of the bin ~-0.15 ... 0.15 + const float z2x = trackResiduals.getZ2X(z2xBin); // z/x coordinate of the bin 0.1 .. 0.9 + const float y = x * y2x; +#if __has_include("TPCFastTransformPOD.h") + const float z = x * z2x * ((iRoc >= geo.getNumberOfSectorsA()) ? -1.f : 1.f); +#else + const float z = x * z2x * ((iRoc >= geo.getNumberOfSlicesA()) ? -1.f : 1.f); +#endif + + float correctionX = useSmoothed ? v->DS[o2::tpc::TrackResiduals::ResX] : v->D[o2::tpc::TrackResiduals::ResX]; + float correctionY = useSmoothed ? v->DS[o2::tpc::TrackResiduals::ResY] : v->D[o2::tpc::TrackResiduals::ResY]; + float correctionZ = useSmoothed ? v->DS[o2::tpc::TrackResiduals::ResZ] : v->D[o2::tpc::TrackResiduals::ResZ]; + + if (invertSigns) { + correctionX *= -1.; + correctionY *= -1.; + correctionZ *= -1.; + } + + entriesStats.emplace_back(voxEntries); + statsPerSecMean[iRoc] += voxEntries; + statsPerSecStdDev[iRoc] += voxEntries * voxEntries; + + nNaNdXV += TMath::IsNaN(correctionX); + nNaNdYV += TMath::IsNaN(correctionY); + nNaNdZV += TMath::IsNaN(correctionZ); + +#if __has_include("TPCFastTransformPOD.h") + float cx, cy, cz; + corr.getCorrectionLocal(iRoc, iRow, y, z, cx, cy, cz); +#else + float u, v, cx, cu, cv, cy, cz; + geo.convLocalToUV(iRoc, y, z, u, v); + corr.getCorrection(iRoc, iRow, u, v, cx, cu, cv); + geo.convUVtoLocal(iRoc, u + cu, v + cv, cy, cz); + cy -= y; + cz -= z; +#endif + + nNaNdXC += TMath::IsNaN(cx); + nNaNdYC += TMath::IsNaN(cy); + nNaNdZC += TMath::IsNaN(cz); + + const float d[3] = {cx - correctionX, cy - correctionY, cz - correctionZ}; + for (int i = 0; i < 3; i++) { + const float dAbs = std::abs(d[i]); + maxDiffPerSec[i][iRoc] = std::max(maxDiffPerSec[i][iRoc], dAbs); + deviations[i].emplace_back(dAbs); + if (std::abs(maxDiff[i]) < dAbs) { + maxDiff[i] = d[i]; + maxDiffRoc[i] = iRoc; + maxDiffRow[i] = iRow; + LOGP(info, "roc {} row {} xyz {} diff {}", iRoc, iRow, i, d[i]); + } + sumDiff[i] += d[i] * d[i]; + } + nDiff++; + + debugVox->Fill(iRoc, iRow, y2xBin, z2xBin, x, y, z, correctionX, correctionY, correctionZ, cx, cy, cz); + + if (lastSector > -1 && lastSector != iRoc) { + if (entriesStats.size() > 0) { + statsPerSecMean[lastSector] /= entriesStats.size(); + statsPerSecStdDev[lastSector] /= entriesStats.size(); + statsPerSecStdDev[lastSector] = (std::sqrt(std::abs(statsPerSecStdDev[lastSector] - statsPerSecMean[lastSector] * statsPerSecMean[lastSector]))); + statsPerSecMedian[lastSector] = mu::median(entriesStats); + } + entriesStats.clear(); + + static std::vector indexDev; + indexDev.resize(deviations[0].size()); + for (int i = 0; i < 3; i++) { + std::array fitRes; + mu::LTMUnbinned(deviations[i], indexDev, fitRes, 0.95); + deviationPerSecLTM95[i][lastSector].mean = fitRes[1]; + deviationPerSecLTM95[i][lastSector].rms = fitRes[2]; + deviationPerSecLTM95[i][lastSector].meanErr = fitRes[3]; + deviationPerSecLTM95[i][lastSector].rmsErr = fitRes[4]; + deviationPerSecMedian[i][lastSector] = mu::median(deviations[i]); + deviations[i].clear(); + } + } + lastSector = iRoc; + } + // last sector + if (lastSector > -1 && entriesStats.size() > 0) { + statsPerSecMean[lastSector] /= entriesStats.size(); + statsPerSecStdDev[lastSector] /= entriesStats.size(); + statsPerSecStdDev[lastSector] /= std::sqrt(std::abs(statsPerSecStdDev[lastSector] - statsPerSecMean[lastSector] * statsPerSecMean[lastSector])); + statsPerSecMedian[lastSector] = mu::median(entriesStats); + entriesStats.clear(); + + std::vector indexDev; + indexDev.resize(deviations[0].size()); + for (int i = 0; i < 3; i++) { + std::array fitRes; + mu::LTMUnbinned(deviations[i], indexDev, fitRes, 0.95); + deviationPerSecLTM95[i][lastSector].mean = fitRes[1]; + deviationPerSecLTM95[i][lastSector].rms = fitRes[2]; + deviationPerSecLTM95[i][lastSector].meanErr = fitRes[3]; + deviationPerSecLTM95[i][lastSector].rmsErr = fitRes[4]; + deviationPerSecMedian[i][lastSector] = mu::median(deviations[i]); + deviations[i].clear(); + } + } + + const int nNaNV = nNaNdXV + nNaNdYV + nNaNdZV; + const int nNaNC = nNaNdXC + nNaNdYC + nNaNdZC; + const auto sNaNV = fmt::format("NaNV: {} {} {} {}", nNaNV, nNaNdXV, nNaNdYV, nNaNdZV); + const auto sNaNC = fmt::format("NaNC: {} {} {} {}", nNaNC, nNaNdXC, nNaNdYC, nNaNdZC); + + if (nNaNV > 0) { + LOGP(error, "{}", sNaNV); + } else { + LOGP(info, "{}", sNaNV); + } + if (nNaNC > 0) { + LOGP(error, "{}", sNaNC); + } else { + LOGP(info, "{}", sNaNC); + } + + summary << "summary" + << "file=" << outFileName + // + << "meanIDC=" << meanIDC + << "meanCTP=" << meanCTP + << "run=" << run + << "validFrom=" << validFrom + << "validUntil=" << validUntil + << "firstTF=" << firstTF + << "lastTF =" << lastTF + // + << "nNaNdXV=" << nNaNdXV + << "nNaNdYV=" << nNaNdYV + << "nNaNdZV=" << nNaNdZV + << "nNaNdXC=" << nNaNdXC + << "nNaNdYC=" << nNaNdYC + << "nNaNdZC=" << nNaNdZC + // + << "statsMean=" << statsPerSecMean + << "statsStdDev=" << statsPerSecStdDev + << "statsMedian=" << statsPerSecMedian + // + << "DdXLTM95=" << deviationPerSecLTM95[0] + << "DdYLTM95=" << deviationPerSecLTM95[1] + << "DdZLTM95=" << deviationPerSecLTM95[2] + << "DdXMedian=" << deviationPerSecMedian[0] + << "DdYMedian=" << deviationPerSecMedian[1] + << "DdZMedian=" << deviationPerSecMedian[2] + << "DdXMax=" << maxDiffPerSec[0] + << "DdYMax=" << maxDiffPerSec[1] + << "DdZMax=" << maxDiffPerSec[2] + // + << "nSlicesPhiZ=" << nSlicesPhiZ + << "maxTracks=" << maxTracks + << "minTracks=" << minTracks + << "nTracksProcessed=" << nTracksProcessed + << "badCalib=" << badCalib + << "\n"; + + summary.Close(); + +#if __has_include("TPCFastTransformPOD.h") + if (debug > 0) { + debugFile->cd(); + TNtuple* ntAll = new TNtuple("all", "all", "sec:row:x:y:z:cx:cy:cz:ix:iy:iz"); + ntAll->SetMarkerStyle(8); + ntAll->SetMarkerSize(0.1); + ntAll->SetMarkerColor(kBlack); + + debugFile->cd(); + TNtuple* ntGrid = new TNtuple("grid", "grid", "sec:row:x:y:z:cx:cy:cz:ix:iy:iz"); + ntGrid->SetMarkerStyle(8); + ntGrid->SetMarkerSize(1.2); + ntGrid->SetMarkerColor(kBlack); + + debugFile->cd(); + TNtuple* ntFitPoints = new TNtuple("fitpoints", "fit points", "sec:row:x:y:z:px:py:pz:cx:cy:cz"); + ntFitPoints->SetMarkerStyle(8); + ntFitPoints->SetMarkerSize(0.4); + ntFitPoints->SetMarkerColor(kRed); + + currDir->cd(); + + auto getInvCorrections = [&](int iSector, int iRow, float realY, float realZ, float& ix, float& iy, float& iz) { + ix = corr.getCorrectionXatRealYZ(iSector, iRow, realY, realZ); + corr.getCorrectionYZatRealYZ(iSector, iRow, realY, realZ, iy, iz); + }; + + auto getAllCorrections = [&](int iSector, int iRow, float y, float z, float& cx, float& cy, float& cz, float& ix, float& iy, float& iz) { + corr.getCorrectionLocal(iSector, iRow, y, z, cx, cy, cz); + getInvCorrections(iSector, iRow, y + cy, z + cz, ix, iy, iz); + }; + + LOGP(info, "create debug ntuples at spline grid points and high granular ..."); + + for (int32_t iSector = 0; iSector < geo.getNumberOfSectors(); iSector++) { + LOGP(info, "debug ntuples for sector {}", iSector); + + for (int32_t iRow = 0; iRow < geo.getNumberOfRows(); iRow++) { + + double x = geo.getRowInfo(iRow).x; + + const auto& gridY = corr.getSplineForRow(iRow).getGridX1(); + const auto& gridZ = corr.getSplineForRow(iRow).getGridX2(); + + { + std::vector points[2], knots[2]; + auto [yMin, yMax] = geo.getRowInfo(iRow).getYrange(); + auto [zMin, zMax] = geo.getZrange(iSector); + + for (int32_t iu = 0; iu < gridY.getNumberOfKnots(); iu++) { + float y, z; + corr.convGridToLocal(iSector, iRow, gridY.getKnot(iu).getU(), 0., y, z); + knots[0].push_back(y); + points[0].push_back(y); + } + for (int32_t iv = 0; iv < gridZ.getNumberOfKnots(); iv++) { + float y, z; + corr.convGridToLocal(iSector, iRow, 0., gridZ.getKnot(iv).getU(), y, z); + knots[1].push_back(z); + points[1].push_back(z); + } + + for (int32_t iyz = 0; iyz <= 1; iyz++) { + std::sort(knots[iyz].begin(), knots[iyz].end()); + std::sort(points[iyz].begin(), points[iyz].end()); + int32_t n = points[iyz].size(); + int nsteps = (iyz == 0) ? 10 : 5; + for (int32_t i = 0; i < n - 1; i++) { + double d = (points[iyz][i + 1] - points[iyz][i]) / nsteps; + for (int32_t ii = 1; ii < nsteps; ii++) { + points[iyz].push_back(points[iyz][i] + d * ii); + } + } + } + points[0].push_back(yMin); + points[0].push_back(yMax); + points[1].push_back(zMin); + points[1].push_back(zMax); + for (int32_t iyz = 0; iyz <= 1; iyz++) { + std::sort(points[iyz].begin(), points[iyz].end()); + } + + for (int32_t iter = 0; iter < 2; iter++) { + std::vector& py = ((iter == 0) ? knots[0] : points[0]); + std::vector& pz = ((iter == 0) ? knots[1] : points[1]); + for (uint32_t iu = 0; iu < py.size(); iu++) { + for (uint32_t iv = 0; iv < pz.size(); iv++) { + float y = py[iu]; + float z = pz[iv]; + float cx{0}, cy{0}, cz{0}, ix{0}, iy{0}, iz{0}; + getAllCorrections(iSector, iRow, y, z, cx, cy, cz, ix, iy, iz); + if (iter == 0) { + ntGrid->Fill(iSector, iRow, x, y, z, cx, cy, cz, ix, iy, iz); + } else { + ntAll->Fill(iSector, iRow, x, y, z, cx, cy, cz, ix, iy, iz); + } + } + } + } + } + + // the data points used in spline fit + auto& fitPoints = mapDirect.getPoints(iSector, iRow); + for (uint32_t ip = 0; ip < fitPoints.size(); ip++) { + auto point = fitPoints[ip]; + float y = point.mY; + float z = point.mZ; + float correctionX = point.mDx; + float correctionY = point.mDy; + float correctionZ = point.mDz; + float cx, cy, cz; + corr.getCorrectionLocal(iSector, iRow, y, z, cx, cy, cz); + ntFitPoints->Fill(iSector, iRow, x, y, z, correctionX, correctionY, correctionZ, cx, cy, cz); + } + } + } + + debugFile->cd(); + ntAll->Write(); + ntGrid->Write(); + ntFitPoints->Write(); + } +#else + if (debug > 0) { + // ntuple with spline grid points + debugFile->cd(); + // ntuple with created TPC corrections + TNtuple* debugCorr = new TNtuple("corr", "corr", "iRoc:iRow:x:y:z:cx:cy:cz"); + + debugCorr->SetMarkerStyle(8); + debugCorr->SetMarkerSize(0.1); + debugCorr->SetMarkerColor(kBlack); + + TNtuple* debugGrid = new TNtuple("grid", "grid", "iRoc:iRow:x:y:z:cx:cy:cz"); + + debugGrid->SetMarkerStyle(8); + debugGrid->SetMarkerSize(1.2); + debugGrid->SetMarkerColor(kBlack); + + // ntuple with data points created from voxels (with data smearing and + // extension to the edges) + TNtuple* debugPoints = new TNtuple("points", "points", "iRoc:iRow:x:y:z:px:py:pz:cx:cy:cz"); + + debugPoints->SetMarkerStyle(8); + debugPoints->SetMarkerSize(0.4); + debugPoints->SetMarkerColor(kRed); + + currDir->cd(); + + LOGP(info, "create debug ntuples at spline grid points and high granular ..."); + + for (int iRoc = 0; iRoc < geo.getNumberOfSlices(); iRoc++) { + LOGP(info, "debug ntuples for roc {}", iRoc); + for (int iRow = 0; iRow < geo.getNumberOfRows(); iRow++) { + + double x = geo.getRowInfo(iRow).x; + + // the correction + + for (double su = 0.; su <= 1.0001; su += 0.01) { + for (double sv = 0.; sv <= 1.0001; sv += 0.1) { + float u, v; + geo.convScaledUVtoUV(iRoc, iRow, su, sv, u, v); + float y, z; + geo.convUVtoLocal(iRoc, u, v, y, z); + float cx, cu, cv; + corr.getCorrection(iRoc, iRow, u, v, cx, cu, cv); + float cy, cz; + geo.convUVtoLocal(iRoc, u + cu, v + cv, cy, cz); + cy -= y; + cz -= z; + debugCorr->Fill(iRoc, iRow, x, y, z, cx, cy, cz); + } + } + + // the spline grid + + const auto& gridU = corr.getSpline(iRoc, iRow).getGridX1(); + const auto& gridV = corr.getSpline(iRoc, iRow).getGridX2(); + for (int iu = 0; iu < gridU.getNumberOfKnots(); iu++) { + // double su = gridU.convUtoX(gridU.getKnot(iu).getU()); + for (int iv = 0; iv < gridV.getNumberOfKnots(); iv++) { + // double sv = gridV.convUtoX(gridV.getKnot(iv).getU()); + float u, v; + corr.convGridToUV(iRoc, iRow, iu, iv, u, v); + float y, z; + geo.convUVtoLocal(iRoc, u, v, y, z); + float cx, cu, cv; + corr.getCorrection(iRoc, iRow, u, v, cx, cu, cv); + float cy, cz; + geo.convUVtoLocal(iRoc, u + cu, v + cv, cy, cz); + cy -= y; + cz -= z; + debugGrid->Fill(iRoc, iRow, x, y, z, cx, cy, cz); + } + } + + // the data points used in spline fit + // (they are kept in + // TPCFastTransformHelperO2::instance()->getCorrectionMap() ) + + o2::gpu::TPCFastSpaceChargeCorrectionMap& map = corrHelper->getCorrectionMap(); + auto& points = map.getPoints(iRoc, iRow); + + for (unsigned int ip = 0; ip < points.size(); ip++) { + auto point = points[ip]; + float y = point.mY; + float z = point.mZ; + float correctionX = point.mDx; + float correctionY = point.mDy; + float correctionZ = point.mDz; + + float u, v, cx, cu, cv, cy, cz; + geo.convLocalToUV(iRoc, y, z, u, v); + corr.getCorrection(iRoc, iRow, u, v, cx, cu, cv); + geo.convUVtoLocal(iRoc, u + cu, v + cv, cy, cz); + cy -= y; + cz -= z; + + debugPoints->Fill(iRoc, iRow, x, y, z, correctionX, correctionY, correctionZ, cx, cy, cz); + } + } + } + + debugFile->cd(); + debugCorr->Write(); + debugGrid->Write(); + debugPoints->Write(); + } +#endif + + for (int i = 0; i < 3; i++) { + sumDiff[i] = sqrt(sumDiff[i]) / nDiff; + } + + LOGP(info, "Max difference in x : {} at ROC {} row {}", maxDiff[0], maxDiffRoc[0], maxDiffRow[0]); + LOGP(info, "Max difference in y : {} at ROC {} row {}", maxDiff[1], maxDiffRoc[1], maxDiffRow[1]); + LOGP(info, "Max difference in z : {} at ROC {} row {}", maxDiff[2], maxDiffRoc[2], maxDiffRow[2]); + LOGP(info, "Mean difference in x,y,z : {} {} {}", sumDiff[0], sumDiff[1], sumDiff[2]); + + corr.testInverse(0); + + debugFile->cd(); + debugVox->Write(); + debugFile->Close(); +}