|
| 1 | +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. |
| 2 | +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. |
| 3 | +// All rights not expressly granted are reserved. |
| 4 | +// |
| 5 | +// This software is distributed under the terms of the GNU General Public |
| 6 | +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". |
| 7 | +// |
| 8 | +// In applying this license CERN does not waive the privileges and immunities |
| 9 | +// granted to it by virtue of its status as an Intergovernmental Organization |
| 10 | +// or submit itself to any jurisdiction. |
| 11 | + |
| 12 | +/// \file checkKeepStepCylinders.macro |
| 13 | +/// \brief Verify that a generated hull contains every recorded hit |
| 14 | +/// |
| 15 | +/// A hull that excludes a position where a hit was produced will lose that hit |
| 16 | +/// in a real run. This reads the hull straight out of a generated |
| 17 | +/// KeepStepCylinders.macro and tests it against the hit positions of an |
| 18 | +/// o2-sim output, per detector. It is a necessary condition, cheap to run, and |
| 19 | +/// independent of the replay machinery. |
| 20 | +/// |
| 21 | +/// root -l -b -q 'checkKeepStepCylinders.macro("KeepStepCylinders.macro", "o2sim")' |
| 22 | + |
| 23 | +#include <TFile.h> |
| 24 | +#include <TTree.h> |
| 25 | +#include <TSystem.h> |
| 26 | +#include <TSystemDirectory.h> |
| 27 | +#include <TString.h> |
| 28 | +#include <TH1.h> |
| 29 | +#include <TObjArray.h> |
| 30 | +#include <TBranch.h> |
| 31 | + |
| 32 | +#include <cstdio> |
| 33 | +#include <fstream> |
| 34 | +#include <regex> |
| 35 | +#include <string> |
| 36 | +#include <vector> |
| 37 | + |
| 38 | +namespace |
| 39 | +{ |
| 40 | +/// pull the rSq[] and edges[] initialisers out of a generated macro |
| 41 | +bool readHull(const char* macro, std::vector<double>& rad, std::vector<double>& edges) |
| 42 | +{ |
| 43 | + std::ifstream in(macro); |
| 44 | + if (!in) { |
| 45 | + printf("[check] cannot open %s\n", macro); |
| 46 | + return false; |
| 47 | + } |
| 48 | + std::string all((std::istreambuf_iterator<char>(in)), std::istreambuf_iterator<char>()); |
| 49 | + auto grab = [&all](const char* name, std::vector<double>& out, bool sqrtIt) { |
| 50 | + std::smatch m; |
| 51 | + std::regex e(std::string(name) + R"(\[\]\s*=\s*\{([^}]*)\})"); |
| 52 | + if (!std::regex_search(all, m, e)) { |
| 53 | + return false; |
| 54 | + } |
| 55 | + std::string body = m[1].str(); |
| 56 | + std::regex num(R"(-?[0-9]+(\.[0-9]+)?)"); |
| 57 | + for (auto it = std::sregex_iterator(body.begin(), body.end(), num); it != std::sregex_iterator(); ++it) { |
| 58 | + out.push_back(std::stod(it->str())); |
| 59 | + } |
| 60 | + if (sqrtIt) { // rSq is written as "R.0f * R.0f", so every value appears twice |
| 61 | + std::vector<double> half; |
| 62 | + for (size_t i = 0; i < out.size(); i += 2) { |
| 63 | + half.push_back(out[i]); |
| 64 | + } |
| 65 | + out.swap(half); |
| 66 | + } |
| 67 | + return true; |
| 68 | + }; |
| 69 | + return grab("rSq", rad, true) && grab("edges", edges, false); |
| 70 | +} |
| 71 | +} // namespace |
| 72 | + |
| 73 | +void checkKeepStepCylinders(const char* macro = "KeepStepCylinders.macro", |
| 74 | + const char* prefix = "o2sim") |
| 75 | +{ |
| 76 | + std::vector<double> rad, edges; |
| 77 | + if (!readHull(macro, rad, edges) || edges.size() != rad.size() + 1) { |
| 78 | + printf("[check] could not parse a hull out of %s\n", macro); |
| 79 | + return; |
| 80 | + } |
| 81 | + printf("[check] hull from %s:\n", macro); |
| 82 | + for (size_t k = 0; k < rad.size(); ++k) { |
| 83 | + printf(" z in [%8.1f, %8.1f) R < %6.1f\n", edges[k], edges[k + 1], rad[k]); |
| 84 | + } |
| 85 | + |
| 86 | + auto inside = [&](double x, double y, double z) { |
| 87 | + if (z < edges.front() || z >= edges.back()) { |
| 88 | + return true; |
| 89 | + } |
| 90 | + for (size_t k = 0; k < rad.size(); ++k) { |
| 91 | + if (z >= edges[k] && z < edges[k + 1]) { |
| 92 | + return (x * x + y * y) < rad[k] * rad[k]; |
| 93 | + } |
| 94 | + } |
| 95 | + return true; |
| 96 | + }; |
| 97 | + |
| 98 | + TString dir = gSystem->DirName(prefix); |
| 99 | + TString base = gSystem->BaseName(prefix); |
| 100 | + TSystemDirectory sdir(dir, dir); |
| 101 | + TList* files = sdir.GetListOfFiles(); |
| 102 | + if (!files) { |
| 103 | + printf("[check] no such directory %s\n", dir.Data()); |
| 104 | + return; |
| 105 | + } |
| 106 | + long long nTot = 0, nOut = 0; |
| 107 | + TIter next(files); |
| 108 | + while (auto* f = (TSystemFile*)next()) { |
| 109 | + TString name = f->GetName(); |
| 110 | + if (!name.BeginsWith(base + "_Hits") || !name.EndsWith(".root")) { |
| 111 | + continue; |
| 112 | + } |
| 113 | + TString det = name; |
| 114 | + det.ReplaceAll(base + "_Hits", ""); |
| 115 | + det.ReplaceAll(".root", ""); |
| 116 | + TFile* file = TFile::Open(dir + "/" + name); |
| 117 | + TTree* t = file ? (TTree*)file->Get("o2sim") : nullptr; |
| 118 | + if (!t) { |
| 119 | + if (file) { |
| 120 | + file->Close(); |
| 121 | + } |
| 122 | + continue; |
| 123 | + } |
| 124 | + // BasicXYZxHit and friends store the position as a ROOT::Math::PositionVector3D |
| 125 | + TString br; |
| 126 | + for (int i = 0; i < t->GetListOfBranches()->GetEntries(); ++i) { |
| 127 | + TString b = t->GetListOfBranches()->At(i)->GetName(); |
| 128 | + if (b.Contains("Hit")) { |
| 129 | + br = b; |
| 130 | + break; |
| 131 | + } |
| 132 | + } |
| 133 | + long long n = 0, bad = 0; |
| 134 | + const TString pos = br + ".mPos.fCoordinates.f"; |
| 135 | + if (!br.IsNull() && t->GetLeaf(pos + "X")) { |
| 136 | + const TString expr = Form("%sX:%sY:%sZ", pos.Data(), pos.Data(), pos.Data()); |
| 137 | + // one hit is one value, not one entry: the value buffer has to be sized |
| 138 | + // to the number of selected values or GetV1..3 return truncated arrays |
| 139 | + Long64_t sel = t->Draw(expr, "", "goff"); |
| 140 | + if (sel > t->GetEstimate()) { |
| 141 | + t->SetEstimate(sel + 1); |
| 142 | + sel = t->Draw(expr, "", "goff"); |
| 143 | + } |
| 144 | + const double *vx = t->GetV1(), *vy = t->GetV2(), *vz = t->GetV3(); |
| 145 | + for (Long64_t i = 0; i < sel; ++i) { |
| 146 | + ++n; |
| 147 | + if (!inside(vx[i], vy[i], vz[i])) { |
| 148 | + ++bad; |
| 149 | + } |
| 150 | + } |
| 151 | + } |
| 152 | + if (n > 0) { |
| 153 | + printf("[check] %-5s hits %10lld outside hull %8lld (%.3f%%)%s\n", |
| 154 | + det.Data(), n, bad, 100. * bad / n, bad ? " <-- HULL TOO SMALL" : ""); |
| 155 | + nTot += n; |
| 156 | + nOut += bad; |
| 157 | + } else { |
| 158 | + printf("[check] %-5s no position leaf (%s), skipped\n", det.Data(), br.Data()); |
| 159 | + } |
| 160 | + file->Close(); |
| 161 | + } |
| 162 | + printf("[check] TOTAL hits %lld, outside hull %lld (%.4f%%)\n", nTot, nOut, |
| 163 | + nTot ? 100. * nOut / nTot : 0.); |
| 164 | +} |
0 commit comments